Role of particle conservation in self-propelled particle systems

Let S(θ) and C(θ) refer to a particle moving in the direction of θ and being associated with the class of single particles or cluster particles, respectively. In the absence of interactions, single particles are assumed to perform a persistent random walk, which we model as a succession of ballistic straight flights, interspersed by self-diffusion ('tumble') events. These tumble events are assumed to occur at a constant rate λ and reorient the particle's orientation θ by a random amount ϑ₀:

$$\begin{equation} S(\theta) \overset{\lambda}{\rightarrow} S(\theta'= \theta + \vartheta_{0}). \end{equation} \tag{ 1 }$$

For simplicity we assume ϑ₀ to be Gaussian distributed,

$$\begin{equation} p_{0}(\vartheta_0) = \frac{1}{\sqrt{2 \pi \sigma_0^2}} \exp{(-\vartheta_0^2/2 \sigma_{0}^2)}, \end{equation} \tag{ 2 }$$

with σ₀ denoting the standard deviation. On time scales much larger than λ⁻¹, this tumbling behavior can be described as conventional particle diffusion, with the particles' diffusion constant being a function of λ and σ₀ [40].

When two single particles S(θ₁) and S(θ₂) collide, they are assumed to assemble a cluster, i.e. each of the two particles becomes a cluster particle (see figure 2(a)):

$$\begin{equation} S(\theta_1) + S(\theta_2) \rightarrow C(\bar \theta + \vartheta) + C(\bar \theta + \vartheta), \end{equation} \tag{ 3 }$$

where³

$$\begin{equation} \bar \theta(\theta_1, \theta_2) = \textstyle\frac{1}{2} (\theta_1 + \theta_2) \end{equation} \tag{ 4 }$$

denotes the average of both pre-collisional angles θ₁ and θ₂, and where ϑ is a random variable which we, again, assume to be Gaussian distributed:

$$\begin{equation} p(\vartheta) = \frac{1}{\sqrt{2 \pi \sigma^2}} \exp{\big(-\vartheta^2/2 \sigma^2\big)}. \end{equation} \tag{ 5 }$$

The rate of binary collisions, such as equation (3), is determined by a particle-shape-dependent differential scattering cross section, which will be discussed below; see section 2.2 and appendix A.

Zoom In Zoom Out Reset image size

Collisions involving cluster particles are distinct from single particle events. Due to the close spatial proximity of particles within each cluster, these collisions correspond to many-particle interactions. Needless to say, a detailed description of cluster formation and the ensuing particle dynamics represents a highly complex matter, requiring explicit consideration of such many-particle interactions. For simplicity, we will resort to the following simplified interaction picture. We assume that (binary) collisions between single particles and cluster particles lead to a condensation process during which the single particle aligns to the cluster particle without changing the direction of the cluster as a whole:

$$\begin{equation} S(\theta_1) + C(\theta_2) \rightarrow C(\theta_2) + C(\theta_2). \end{equation} \tag{ 6 }$$

Equation (6) thus captures the net effect of collisions between single particles and cluster particles, during which multiple collisions, involving neighboring particles belonging to the same cluster, stabilize the cluster's direction; cf figure 2(b) for an illustration.

Collisions among cluster particles is an even more intricate process, since they actually depend on the size and shape of both colliding clusters, and in general involve multi-particle interactions. In the framework of a Boltzmann-like description, correlations in the particle distribution are neglected and only binary interactions are considered. The frequency of interactions is determined by a geometrical construction called the 'Boltzmann cylinder', assuming that particle positions are homogeneously distributed on local scales. With regard to many-particle interactions during collisions among cluster particles, we thus have to resort to some kind of simplified, binary collision picture. Since our kinetic model lacks any direct notion of cluster size or shape, we will stick to the assumption that, on average, collisions between cluster particles are devoid of any directional bias, leading to the same type of collision rule as for single particles:

$$\begin{equation} C(\theta_1) + C(\theta_2) \rightarrow C(\bar \theta + \vartheta) + C(\bar \theta + \vartheta). \end{equation} \tag{ 7 }$$

Again, ϑ constitutes a Gaussian-distributed random variable given in (5). Moreover, due to external (e.g. thermal background) and internal (e.g. noisy propelling mechanism) noise, cluster particles evaporate to become single particles. In analogy to the self-diffusion of single particles, we thus introduce a rate⁴ characterizing the following evaporation process:

$$\begin{equation} C(\theta) \overset{\epsilon}{\rightarrow} S(\theta' =\theta + \vartheta_0). \end{equation} \tag{ 8 }$$

Also in this case, the strength of the angular changes is Gaussian distributed according to (2), and for simplicity we use the same standard deviation σ₀ as for the single particles' persistent random walk. As discussed above, cluster particles are strongly caged due to their close proximity to neighboring, collinearly moving particles. Reorientations of cluster particles due to noise are therefore strongly counteracted by realigning particle collisions, rendering cluster particles considerably less susceptible to random fluctuations than single particles. Hence, we assume

$$\begin{equation} \epsilon\ll\lambda, \end{equation} \tag{ 9 }$$

which is consistent with the observations in agent-based simulations slightly below the ordering transition [31], finding coherently moving clusters in an unpolarized background of randomly moving particles. In this regime individual particles exhibit superdiffusive behavior, performing quasi-ballistic 'flights' as long as they are part of a cluster, and conventional particle diffusion if they are not.

Building on the modeling framework defined above, we now set up a kinetic description for the canonical model. We denote by s(θ,x,t) and c(θ,x,t) the one-particle distribution functions within the class of single particles and cluster particles, respectively, i.e. s(θ,x,t) dθ d²x gives the number of single particles located in an infinitesimal region [x,x + dx] with orientations in the interval [θ,θ + dθ] (and likewise for c(θ,x,t) dθ d²x). Both one-particle distribution functions are subject to convection due to the propelling velocity v of each particle. Moreover, local fluctuations in the one-particle distribution functions due to self-diffusion and collision events are to be accounted for. We thus arrive at the following set of Boltzmann-like equations for the canonical model:

$$\begin{equation} \partial_t s(\theta,{\mathbf{x}},t)+\mathbf{v}\cdot\nabla s(\theta,{\mathbf{x}},t) = \dot s(\theta,{\mathbf{x}},t), \end{equation} \tag{ 10a }$$

$$\begin{equation} \partial_t c(\theta,{\mathbf{x}},t)+\mathbf{v}\cdot\nabla c(\theta,{\mathbf{x}},t) = \dot c(\theta,{\mathbf{x}},t), \end{equation} \tag{ 10b }$$

$\dot s(\theta ,{\mathbf {x}},t)$

$\dot c(\theta ,{\mathbf {x}},t)$

$$\begin{equation} \dot s = \lambda\left[\mathcal{D}^{(+)}_s(\theta)-\mathcal{D}^{(-)}_s(\theta)\right]+\epsilon\mathcal{D}^{(+)}_c(\theta) -\mathcal{C}^{(-)}_s(\theta)-\mathcal{A}[s,c;\theta], \end{equation} \tag{ 11a }$$

$$\begin{equation} \dot c = -\epsilon \mathcal{D}^{(-)}_c(\theta)+\mathcal{C}^{(+)}_s(\theta)+\mathcal{C}^{(+)}_c(\theta)+\mathcal{A}[c,s;\theta]-\mathcal{C}^{(-)}_{c}(\theta). \end{equation} \tag{ 11b }$$

• 1.  
  Self-diffusion and evaporation. In these cases the source terms are products of the corresponding rates and probability densities, with

  $$\begin{equation} \mathcal{D}^{(-)}_{f}(\theta) = f(\theta) \end{equation} \tag{ 12 }$$

  denoting the probability density for a particular species f to have a certain angle θ, and

  $$\begin{equation} \mathcal{D}^{(+)}_{f}(\theta) = \left \langle f(\theta-\vartheta_0) \right\rangle_0 \end{equation} \tag{ 13 }$$

  denoting the transition probability from θ' = θ − ϑ₀ to θ averaged over all ϑ₀ with respect to the Gaussian weight (2). Note that, here and in the following, the argument of f is understood modulo 2π.
• 2.  
  Collisions within the same class of particles. The collision integrals, representing the processes defined in equations (3) and (7), are given by standard expressions [16–18, 22]

  $$\begin{equation} \mathcal{C}^{(+)}_f(\theta) = \left\langle \int {\mathrm{d}}\mathcal{I}\,f(\theta')f(\theta'') \, \delta(\bar \theta(\theta', \theta'')+\vartheta-\theta) \right\rangle, \end{equation} \tag{ 14a }$$

  $$\begin{equation} \mathcal{C}^{(-)}_f(\theta) = \int {\mathrm{d}}\mathcal{I}\,f(\theta')f(\theta'')\delta(\theta'-\theta). \end{equation} \tag{ 14b }$$

  Here 〈...〉 denotes an average over ϑ∈(−∞,∞) with respect to the Gaussian weight (5) and the average angle $\bar \theta$ is given in equation (4). The integral measure

  $$\begin{equation} \int {\mathrm{d}}\mathcal{I}\, (\ldots) \equiv \int_{-\pi}^{\pi}{\mathrm{d}}\theta' \int_{\theta'-\pi}^{\theta'+\pi}{\mathrm{d}}\theta'' \, \Gamma(L,d,|\theta'-\theta''|)\,(\ldots) \end{equation} \tag{ 15 }$$

  contains the differential scattering cross section

  $$\begin{equation} \fl \Gamma(L,d,|\theta'-\theta''|)=4\,{\mathrm{d}}v\left|\sin\left(\frac{\theta'-\theta''}{2}\right)\right|\left[1+\frac{(L/d)-1}{2}\left|\sin(\theta'-\theta'')\right|\right] \end{equation} \tag{ 16 }$$

  characterizing the frequency of collisions (i.e. hard-core interactions) between rod-like particles. The scattering function Γ itself carries all information concerning the shape of the particles and is a function of the relative orientation of the colliding particles. Reminiscent of the Boltzmann scattering cylinder, Γ can be derived on the basis of purely geometric considerations assuming that all spatial coordinates within the cylinder are equally probable; for details see appendix A.
• 3.  
  Assembly events of a single particle joining a cluster. These events, represented by (6), occur through binary collisions between single particles and cluster particles and are thus represented by an analogous integral expression:

  $$\begin{equation} \mathcal{A}[f,g;\theta] = \int_{}^{}{\mathrm{d}}\mathcal{I}\,f(\theta')g(\theta'')\delta(\theta'-\theta). \end{equation} \tag{ 17 }$$

Equations (21a) and (21a) constitute an infinite set of coupled equations in Fourier space, which are fully equivalent to the Boltzmann-like equations (10a) and (10b). To derive a closed set of hydrodynamic equations, we need to consider some additional assumptions, allowing us to truncate this infinite Fourier space representation.

Here, our focus will be on virtually isotropic systems in the vicinity of an ordering transition breaking rotational symmetry. In this case, deviations of the one-particle distribution functions from the constant distribution ∼1/2π are small and contributions from large wavenumbers in the Fourier series (21a) and (21a) are negligible. We further consider sufficiently dilute systems, in which the number of (binary) particle collisions per unit time and area [ ∼ (ρ_c + ρ_s)² I_0,0] is much smaller than the corresponding number of single particle diffusion events [ ∼ λ ρ_s]. Together with  ≪ λ (equation (9)), stating that disassembly from a cluster is strongly hindered by particle caging, allows us to treat single particle diffusion as a fast process. The single particle phase thus acts as an isotropic sea of particles where particle orientations (but not necessarily particle densities) are equilibrated, and hence the net hydrodynamic flow vanishes (u_s = 0). Finally, from a dimensional analysis of equations (19a) and (19b), together with (20c) and (20d), one finds $c_k/\rho _c \sim \mathcal {O}(|{\mathbf {u}}_c|^k/v^k)$. Near the onset of order, where |u_c|/v ≪ 1, we only consider the density (c₀) and polarity (c₁) of cluster particles, and use the stationary equation for c₂ as a closure relation, neglecting all contributions from higher order coefficients.

In summary, we resort to the following truncation scheme, leading to a set of hydrodynamic equations, valid near the onset of the ordering transition:

$$\begin{equation} s_k = 0, \quad\forall |k|>0, \end{equation} \tag{ 23a }$$

$$\begin{equation} c_k = 0, \quad\forall |k|>2. \end{equation} \tag{ 23b }$$

With the above truncation scheme, (21a) and (21b) reduce to

$$\begin{equation} \fl \partial_ts_0 = \epsilon c_0-I_{0,0}\left(s_0^2+s_0c_0\right), \end{equation} \tag{ 24a }$$

$$\begin{equation}\fl \partial_tc_0 = -v\left[\partial_x \mathfrak{R}(c_1)+\partial_y\mathfrak{I}(c_1)\right]-\partial_ts_0, \end{equation} \noindent \tag{ 24b }$$

$$\begin{eqnarray} &&\fl \partial_tc_1 = -\frac{v}{2}[\partial_x(c_2+c_0)-\mathrm{i}\partial_y(c_2-c_0)]+ [(2\,\mathrm{e}^{-\sigma^2/2} I_{1,1} - I_{1,0} - I_{0,0}) c_0 -\epsilon + I_{0,0} s_0] c_1 \nonumber\\ &&+ [ 2\,\mathrm{e}^{-\sigma^2/2} I_{2,1} - I_{1,0} - I_{2,0}] c_1^* c_2, \end{eqnarray} \noindent \tag{ 24c }$$

$$\begin{eqnarray} &&\fl \partial_t c_2 = -\frac{v}{2}[\partial_x+\mathrm{i}\partial_y]c_1+ [ (2\,\mathrm{e}^{-2\sigma^2} I_{1,0} - I_{2,0} - I_{0,0} ) c_0 -\epsilon + I_{0,0} s_0 ] c_2 \nonumber\\ &&+ [ \mathrm{e}^{-2\sigma^2} I_{0,0} - I_{1,0} ] c_1 c_1, \end{eqnarray} \tag{ 24d }$$

$f(\theta )\in \mathbb R$

$\mathfrak {R}(a)$

$\mathfrak {I}(a)$

Table 1. Summary of relevant collision integrals I_n,k as a function of the aspect ratio ξ = L/d, where L and d denote particle length and diameter, and where v is the particle velocity. The quantities I_n,k/I_0,0 depend only weakly on the aspect ratio ξ. In particular, the signs of I_n,k/I_0,0 do not change with ξ, leaving all our present conclusions made on the basis of the kinetic coefficients qualitatively unchanged.

Integral	I0,0	I1,0/I0,0	I1,1/I0,0	I2,0/I0,0	I2,1/I0,0
Value	$\frac {8\,{\mathrm {d}}v(2+\xi )}{3\pi }$	$-\frac {4+\xi }{5(2+\xi )}$	$\frac {3}{16} \frac {8+\pi (\xi -1)}{2+\xi }$	$\frac {6-13\xi }{35(2+\xi )}$	$\frac {3}{16} \frac {\pi (1-\xi )-8}{2+\xi }$

For given particle densities, the time scales governing the dynamics of the polar and nematic order parameter fields, represented by c₁ and c₂, are given by the linear coefficients in the second line of (24c) and (24d), respectively. As will be detailed in section 4.2, the onset of collective motion is hallmarked by a change in sign of the linear coefficient in (24c), implying a diverging time scale for the dynamics of the polarity field. On the other hand, the time scale for c₂ is finite for all densities, which implies that the relaxation of the nematic order parameter field is fast compared to the polarity field. This allows us to set ∂_tc₂ ≈ 0 in (24d).

In the following, it will be convenient to write down equations in dimensionless form. To this end, we construct the following characteristic scales: time and space will be measured in units of the cluster evaporation time and length scale

$$\begin{equation} \hat \tau_{\rm e}=\epsilon^{-1} \quad\mbox{and}\quad \hat\ell_{\rm e}=v/\epsilon. \end{equation} \tag{ 25 }$$

From the cluster evaporation time scale $\hat {\tau }_{\mathrm { e}}$ and the total scattering cross section I_0,0, we can construct the characteristic density scale

$$\begin{equation} \hat\rho_b = \frac{1}{I_{0,0}\,\hat\tau_{\rm e}}. \end{equation} \tag{ 26 }$$

The single particle and cluster particle phases constantly exchange particles at rates that are determined by cluster evaporation () on the one hand (cluster particles → single particles) and cluster nucleation due to particle collisions on the other hand (single particles →cluster particles) which occur with a rate ∼ρ I_0,0. Therefore the characteristic density scale $\hat \rho _b$ marks the particle density, where both rates balance. In particular, $\rho /\hat {\rho }_b=(\rho _s+\rho _c)/ \hat {\rho }_b$ gives the rate of inter-particle collisions relative to cluster evaporation events. Thus, the numerical quantity $\rho / \hat {\rho }_b$ provides a direct measure expressing the competition between the randomizing effects of noise and the order creating effects of particle collisions, hallmarking the onset (and maintenance) of collective motion [29].

We thus arrive at the following rescaling scheme:

$$\begin{equation} t \rightarrow t\cdot\hat\tau_{\rm e}, \end{equation} \tag{ 27a }$$

$$\begin{equation} {\mathbf{x}} \rightarrow {\mathbf{x}}\cdot\hat\ell_{\rm e}, \end{equation} \tag{ 27b }$$

$$\begin{equation} \rho_{s/c} \rightarrow \rho_{s/c}\cdot\hat\rho_b, \end{equation} \tag{ 27c }$$

$$\begin{equation} {\mathbf{g}} \rightarrow {\mathbf{g}}\cdot \hat\rho_b\,\frac{\hat\ell_{\rm e}}{\hat\tau_{\rm e}}, \end{equation} \tag{ 27d }$$

$$\begin{equation} I_{n,k} \rightarrow I_{n,k}\cdot \frac{1}{\hat\rho_b\,\hat\tau_{\rm e}}, \end{equation} \tag{ 27e }$$

$\hat \rho _b$

Then, upon eliminating c₂ from (24c) as discussed above, and using the relations between Fourier modes and hydrodynamic fields (for details see appendix B), (20a)–(20d), equations (24a)–(24d) give rise to the hydrodynamic equations corresponding to the canonical model. In rescaled variables they read

$$\begin{equation}\fl \partial_t\rho_s = \rho_c-\left(\rho_s+\rho_c\right)\rho_s, \end{equation} \tag{ 28a }$$

$$\begin{equation}\fl \partial_t\rho_c = -\nabla \cdot {\mathbf{g}}-\rho_c+\left(\rho_s+\rho_c\right)\rho_s, \end{equation} \tag{ 28b }$$

$$\begin{eqnarray} &&\fl \partial_t{\mathbf{g}} = -\nu _1 {\mathbf{g}} -\frac{\mu \kappa }{\nu _2}{\mathbf{g}}^2 {\mathbf{g}}-\frac{1}{2}\nabla \rho _c+\frac{1}{4\nu _2}\nabla ^2{\mathbf{g}}+\frac{\zeta _+}{\nu _2}({\mathbf{g}}\cdot \nabla ){\mathbf{g}}+\frac{\zeta _-}{\nu _2}[(\nabla \cdot {\mathbf{g}}){\mathbf{g}}-\frac{1}{2}\nabla ({\mathbf{g}}^2)]\nonumber\\ &&+\frac{\mu }{\nu _2^2}[{\mathbf{g}}({\mathbf{g}}\cdot {\boldsymbol{\partial}} [\rho_c,\rho_s])-\frac{1}{2}{\mathbf{g}}^2{\boldsymbol{\partial}}[\rho_c,\rho_s]]+\frac{1}{4\nu _2^2}[(\nabla \cdot {\mathbf{g}}){\boldsymbol{\partial}}[\rho_c,\rho_s]\nonumber\\ &&- (\nabla {\mathbf{g}}+\nabla {\mathbf{g}}^t){\boldsymbol{\partial}}[\rho_c,\rho_s]], \end{eqnarray} \tag{ 28c }$$

$$\begin{equation} {\boldsymbol{\partial}}[f,g]=\left(\partial _{f}\nu _2\right)\nabla f+\left(\partial _{g}\nu _2\right)\nabla g, \end{equation} \tag{ 29 }$$

and where we have introduced the following abbreviations:

$$\begin{equation} \nu_1 = 1-(\rho _s-\rho _c)+(I_{1,0}-2e^{-\sigma ^2/2}I_{1,1})\rho _c, \end{equation} \tag{ 30a }$$

$$\begin{equation} \nu_2 = 1-(\rho _s-\rho _c)+(I_{2,0}-2e^{-2\sigma ^2}I_{1,0})\rho _c, \end{equation} \tag{ 30b }$$

$$\begin{equation} \mu = \mathrm{e}^{-2\sigma ^2}-I_{1,0}, \end{equation} \tag{ 30c }$$

$$\begin{equation} \kappa = I_{1,0}+I_{2,0}-2\,\mathrm{e}^{-\sigma ^2/2}I_{2,1}, \end{equation} \tag{ 30d }$$

$$\begin{equation} \zeta _{\pm }=-\mu \pm \frac{\kappa}{2}. \end{equation} \tag{ 30e }$$

Equations (28a)–(28c) capture the evolution of our canonical model system on a hydrodynamic level. More specifically, (28a) and (28b) describe the spatio-temporal evolution of the particle densities ρ_s and ρ_c. Since, by the assumptions underlying our model, no macroscopic flow of single particles can build up, only the density of cluster particles (ρ_c) is subject to convection. This implies that the genuine hydrodynamic momentum field g = g_c + g_s ≡ g_c is carried solely by the subset of cluster particles. Therefore we omit the subscript c in (28b) and (28c) and denote g ≡ g_c. The dynamics of both densities is, moreover, driven by source terms, as determined by the reactions discussed in section 2.1. The gain and loss parts in these source terms of ρ_c and ρ_s are exactly balanced, such that the total density ρ = ρ_c + ρ_s is conserved. As an aside we note that any distinction between single particles and cluster particles is a purely conceptual matter. Experimentally, only the total density ρ and the momentum field g are accessible.

Equation (28c), governing the evolution of the current density g, can be interpreted as a generalization of the Navier–Stokes equation to active systems. The terms on the right-hand side of (28c) can be given by the following interpretation. In the first line, the first two terms account for the local dynamics of g. They play a crucial role in establishing and maintaining a state of macroscopic flow, as will be detailed below. The Navier–Stokes equation itself, which conserves momentum, is devoid of these terms. In formal analogy to the Navier–Stokes equation, the density gradient in the first line together with the last term in the second line can be interpreted as a pressure gradient. This effective pressure is given by $\frac {1}{2}(\rho _c+\frac {\zeta _-}{\nu _2}\,{\mathbf {g}}^2)$, when neglecting the density dependence of ν₂. The last term in the first line is analogous to the shear stress term in the Navier–Stokes equation, with a kinematic viscosity ∼ν⁻¹₂. The second line in (28c) is a generalization of the convection term to systems not obeying Galilean invariance, where all combinations of ∇ and factors second order in g transforming as vectors are allowed [41]. Finally, the last two lines describe couplings of the current density g and gradients thereof to density gradients. Note that the density gradients in these coupling terms are all of the same generic structure (29).

As already noted, the canonical model equations (28a)–(28c) conserve the total number of particles. To make this explicit, we define

$$\begin{equation} \rho \equiv \rho_c + \rho_s, \end{equation} \tag{ 31a }$$

$$\begin{equation} \eta \equiv \rho_c - \rho_s, \end{equation} \tag{ 31b }$$

$$\begin{equation} \fl \partial_t\rho = -\nabla \cdot {\mathbf{g}}, \end{equation} \tag{ 32a }$$

$$\begin{equation} \fl \partial_t\eta = -\nabla \cdot {\mathbf{g}}+\rho^2-(\rho+1)\eta-\rho, \end{equation} \tag{ 32b }$$

$$\begin{eqnarray} &&\fl \partial_t{\mathbf{g}} = -\nu_1 {\mathbf{g}}-\frac{\mu \kappa }{\nu _2}{\mathbf{g}}^2{\mathbf{g}}-\frac{1}{4}\nabla (\rho+\eta)+\frac{1}{4\nu _2}\nabla ^2{\mathbf{g}}+\frac{\zeta _+}{\nu _2}({\mathbf{g}}\cdot \nabla ){\mathbf{g}}+\frac{\zeta _-}{\nu _2}\left[(\nabla \cdot {\mathbf{g}}){\mathbf{g}}-\frac{1}{2}\nabla \left({\mathbf{g}}^2\right)\right]\nonumber\\ &&+\frac{\mu }{\nu _2^2}\left[{\mathbf{g}}({\mathbf{g}}\cdot {\boldsymbol{\partial}} [\rho,\eta])-\frac{1}{2}{\mathbf{g}}^2{\boldsymbol{\partial}}[\rho,\eta]\right]+\frac{1}{4\nu _2^2}[(\nabla \cdot {\mathbf{g}}){\boldsymbol{\partial}}[\rho,\eta] - (\nabla {\mathbf{g}}+\nabla {\mathbf{g}}^t){\boldsymbol{\partial}}[\rho,\eta]].\nonumber\\ \end{eqnarray} \tag{ 32c }$$

Now we turn to the grand canonical model, where the single particle phase is coupled to a particle reservoir, resulting in a situation where single particles constitute an isotropic sea of particles that is maintained at a constant density ρ⁰_s. Particle number conservation is now violated, and the only non-trivial density dynamics takes place within the phase of cluster particles. The hydrodynamic equations corresponding to the grand canonical model can be obtained immediately by setting in (28a)–(28c) the density of single particles to a constant value, yielding

$$\begin{equation} \fl \rho_s = \rho_{s}^{0} = \mbox{const.,} \end{equation} \tag{ 33a }$$

$$\begin{equation} \fl \partial_t\rho_c = -\nabla \cdot {\mathbf{g}}-\rho_c+\left(\rho_{s}^{0}+\rho_c\right)\rho_{s}^{0}, \end{equation} \tag{ 33b }$$

$$\begin{eqnarray} &&\fl \partial_t{\mathbf{g}} = -\nu _1 {\mathbf{g}}-\frac{\mu \kappa }{\nu _2}{\mathbf{g}}^2 {\mathbf{g}}-\frac{1}{2}\nabla \rho _c+\frac{1}{4\nu _2}\nabla ^2{\mathbf{g}}+\frac{\zeta _+}{\nu _2}({\mathbf{g}}\cdot \nabla ){\mathbf{g}}+\frac{\zeta _-}{\nu _2}\left[(\nabla \cdot {\mathbf{g}}){\mathbf{g}}-\frac{1}{2}\nabla \left({\mathbf{g}}^2\right)\right]\nonumber\\ &&+\frac{\mu \partial_{\rho_c} \nu_2}{\nu _2^2}\left[{\mathbf{g}}({\mathbf{g}}\cdot \nabla \rho_c)-\frac{1}{2}{\mathbf{g}}^2\nabla \rho_c \right]+\frac{\partial_{\rho_c} \nu_2}{4\nu _2^2}[ (\nabla \cdot {\mathbf{g}})\nabla \rho_c - (\nabla {\mathbf{g}}+\nabla {\mathbf{g}}^t)\nabla \rho_c ]. \end{eqnarray} \tag{ 33c }$$

To assess the density difference between the cluster particle and the single particle phase η, we calculate the dynamical fixed point η* of (34b), attracting the dynamics of η(t) in the long time limit t → ∞:

$$\begin{equation} \eta^*(\rho) = \frac{\rho^2-\rho}{\rho+1}. \end{equation} \tag{ 36 }$$

The defining equations (31a)–(31b) can be used to determine the corresponding (stationary) fixed point densities of single particles (ρ*_s) and cluster particles (ρ*_c) as a function of the total density ρ. Figure 3 summarizes these findings. Upon increasing the total density ρ, the ratio η*/ρ continuously grows from η*/ρ = −1 at ρ = 0, asymptotically approaching η*/ρ = 1 as ρ → ∞. Based on the sign of η*, two density regimes can be distinguished. In the low-density regime (ρ ≪ 1, η* < 0), particle collisions, underlying the formation of clusters, occur at much smaller rates than cluster evaporation events. Only a small fraction of all particles organize themselves in clusters, leading to a relatively dense population of single particles and a correspondingly small density of cluster particles. In the high-density regime (ρ ≫ 1, η* > 0), the situation is reversed: large overall densities imply frequent particle collisions and, consequently, cluster formation and cluster growth dominate over cluster evaporation. In this regime, the number of cluster particles exceeds the number of single particles.

Zoom In Zoom Out Reset image size

The crossover between the single particle dominated low-density regime and the cluster particle dominated high-density regime occurs at the crossover density $\bar \rho =\hat \rho _b=1$, where both the single particle and the cluster particle populations are of equal size (i.e. $\eta ^*(\bar \rho )=0$). The relation between the crossover density to clustering, and the geometrical shape of the constituent particles has been addressed previously in [37], based on agent-based simulations and a mean-field-type analytical analysis. Using our definition of the crossover density $\bar \rho$, we can establish the corresponding relation simply by restoring original units (equation (26)). Using packing fraction $\bar p\simeq \bar \rho \,Ld$ instead of particle density, and assuming for the sake of simplicity L/d ≫ 1, which allows us to estimate the particle surface A₀ ≃ Ld, we find that

$$\begin{equation} \bar p \simeq \frac{\epsilon\,Ld}{I_{0,0}} = \frac{3\pi}{8v}\frac{\epsilon L}{2+L/d}, \end{equation} \tag{ 37 }$$

which correctly reproduces the findings of [37] (taking into account that the cluster evaporation rate is assumed to be proportional to the inverse particle length, ∝L⁻¹). For the sake of completeness, we note that the definition of the clustering crossover density in  [37] is based on the cluster size distribution, and thus does not necessarily coincide with our definition. We stress, however, that in our description the scaling structure in equation (37) is completely generic. It is an immediate consequence of the characteristic scales of our model and of the fact that the rescaled hydrodynamic model equations are (virtually) independent of particle shape. The structure of equation (37) is thus robust under an arbitrary redefinition of the (rescaled) crossover density $\bar \rho$.

Having examined the composition of the system in terms of single particle and cluster particle densities, we now turn to a discussion of the spatially homogeneous solutions for the momentum current density g. Due to rotational invariance of (34c), only the magnitude g = |g| of the momentum current density, but not its direction, evolves in time. We can thus concentrate on the scalar equation

$$\begin{equation} \partial_t g = -\nu _1\,g -\frac{\mu \kappa }{\nu _2}\,g^3, \end{equation} \tag{ 38 }$$

which leads to the following fixed points g* as the attractor of the dynamics of g in the limit of long times:

$$\begin{equation} g^*=\left\{\begin{array}{@{}ll@{}} 0 & {\rm for} \ \nu_1>0,\\ g_0=\sqrt{-\frac{\nu_1\nu_2}{\mu\kappa}} & {\rm for} \ \nu_1<0. \end{array}\right. \end{equation} \tag{ 39 }$$

It can be shown that the coefficient in front of the cubic term in (38) is indeed strictly positive for all control parameters of density ρ and noise σ consistent with ν₁ < 0, ensuring the existence of the non-trivial fixed point in the second line of (39).

Depending on the sign of the linear coefficient ν₁, two parameter regimes can thus be distinguished. Parameters leading to ν₁ > 0 render stable an overall homogeneous and isotropic state with vanishing macroscopic flow g = 0. Upon crossing the phase boundary

$$\begin{equation} \nu_1(\rho,\sigma) = 0 \end{equation} \tag{ 40 }$$

in parameter space, the isotropic solution gets unstable and a macroscopic current density of non-zero amplitude builds up. In equation (40) we used the fact that the density difference η, in the stationary limit, is a function of the total density ρ; cf equation (36). Hence, in the limit of long times, ν₁ is a function of the total density ρ and the noise parameter σ, only.

Using the definition of the coefficient ν₁, equation (30a), we can readily calculate the shape of the phase boundary in the σρ-plane:

$$\begin{equation} \sigma_c(\rho) = \sqrt{-2 \, \ln{\left(\frac{2}{3} + \rho^{-2}\right)}}\quad(\rho \geqslant \sqrt{3}), \end{equation} \tag{ 41 }$$

where we used $I_{1,0}=-\frac {1}{3}$ and $I_{1,1}=\frac {1}{2}$. The corresponding phase diagram is shown in figure 4.

Zoom In Zoom Out Reset image size

To conclude this section, we note that the analysis of spatially homogeneous systems corroborates the general physical picture of active systems, which was alluded to in the introduction (e.g. cf [31]). Even in the absence of noise, σ = 0, for which the threshold density ρ^(c) is lowest, the fully isotropic state g = 0 remains stable up to a critical density $\rho ^{(c)}(\sigma =0)=\sqrt {3}\,\bar \rho$, which lies well beyond the density $\bar \rho$ indicating the crossover to clustering. We thus extract the following physical picture; cf figure 4. For low densities, $\rho <\bar \rho$, cluster evaporation dominates over cluster assembly via particle collisions and clusters form only transiently. The system most closely resembles a structureless, isotropic 'sea of particles'. At intermediate densities, $\bar \rho <\rho <\rho ^{(c)}$, particle collisions are more frequent. The emergence of clusters is now a virtually persistent phenomenon, with cluster evaporation occurring at a lower rate than cluster formation and growth. Yet, the collision rates between clusters (i.e. collisions among cluster particles) are still too low to orchestrate macroscopic order, leading to an overall isotropic 'sea of clusters'. Finally, for large densities, ρ > ρ^(c), the frequency of collisions among clusters is high enough to establish collective motion even on macroscopic scales.

We start by linearizing the hydrodynamic equations for the canonical model. In the canonical model, the total number of particles is conserved, and the appropriate base state reads (cf section 4)

$$\begin{equation} \rho = \rho_{\mathrm{h}} = \mbox{const.}, \end{equation} \tag{ 42a }$$

$$\begin{equation} \eta = \eta^*(\rho_{\mathrm{h}}) = \frac{\rho_{\mathrm{h}}^2-\rho_{\mathrm{h}}}{\rho_{\mathrm{h}}+1}, \end{equation} \tag{ 42b }$$

$$\begin{equation} {\mathbf{g}} = {\mathbf{g}}_{\mathrm{h}} \in \left\{0,\;g_0=\sqrt{-\frac{\nu_1\nu_2}{\mu\kappa}}\right\}\hat{\mathbf{e}}_g, \end{equation} \tag{ 42c }$$

$\hat {\mathbf {e}}_g$

$$\begin{equation} \rho({{\mathbf{x}}},t) = \rho_{h} + \delta \rho({{\mathbf{x}}}, t), \end{equation} \tag{ 43a }$$

$$\begin{equation} \eta({{\mathbf{x}}},t) = \eta^{*} + \delta \eta({{\mathbf{x}}}, t), \end{equation} \tag{ 43b }$$

$$\begin{equation} {\mathbf{g}}({{\mathbf{x}}},t) = {\mathbf{g}}_{\mathrm{h}} + \delta {\mathbf{g}}({{\mathbf{x}}}, t), \end{equation} \tag{ 43c }$$

$$\begin{equation} \delta \rho ({\mathbf{x}},t) = \delta \rho_0\,\mathrm{e}^{{\rm st}+\mathrm{i} \mathbf{q}\cdot {\mathbf{x}}}, \end{equation} \tag{ 44a }$$

$$\begin{equation} \delta \eta ({\mathbf{x}},t) = \delta \eta_0 \,\mathrm{e}^{{\rm st}+\mathrm{i} \mathbf{q}\cdot {\mathbf{x}}}, \end{equation} \tag{ 44b }$$

$$\begin{equation} \delta {\mathbf{g}} ({\mathbf{x}},t) = \delta {\mathbf{g}}_0\,\mathrm{e}^{{\rm st}+\mathrm{i} \mathbf{q}\cdot {\mathbf{x}}}. \end{equation} \tag{ 44c }$$

$$\begin{equation}\fl s\,\delta \rho_0 = -\mathrm{i} {\mathbf{q}}\cdot \delta{\mathbf{g}}_0, \end{equation} \tag{ 45a }$$

$$\begin{equation}\fl s\,\delta\eta_0 = \left(2\rho_{\mathrm{h}}-\eta^*-1\right)\delta \rho_0 -\left(1+\rho _{\mathrm{h}}\right)\delta\eta_0-\mathrm{i} {\mathbf{q}}\cdot \delta{\mathbf{g}}_0, \end{equation} \tag{ 45b }$$

$$\begin{eqnarray} &&\fl s\,\delta{\mathbf{g}}_0 = \left[\frac{\partial_{\rho}\nu_2\,\mu }{\nu _2^2}\left(\kappa\,{\mathbf{g}}_{\mathrm{h}}^2\,{\mathbf{g}}_{\mathrm{h}}+\left({\mathbf{g}}_{\mathrm{h}}\cdot\mathrm{ i}\mathbf{q}\right){\mathbf{g}}_{\mathrm{h}}-\frac{\mathrm{i}\mathbf{q}}{2}\,{\mathbf{g}}_{\mathrm{h}}^2\right)-\partial_{\rho}\nu_1\,{\mathbf{g}}_{\mathrm{h}}-\frac{\mathrm{i}\mathbf{q}}{4}\right]\delta \rho\nonumber\\ &&+\left[\frac{\partial_{\eta}\nu_2\,\mu }{\nu _2^2}\left(\kappa\,{\mathbf{g}}_{\mathrm{h}}^2\,{\mathbf{g}}_{\mathrm{h}}+\left({\mathbf{g}}_{\mathrm{h}}\cdot \mathrm{i}\mathbf{q}\right){\mathbf{g}}_{\mathrm{h}}-\frac{\mathrm{i}\mathbf{q}}{2}\,{\mathbf{g}}_{\mathrm{h}}^2\right)-\partial_{\eta}\nu_1\,{\mathbf{g}}_{\mathrm{h}}-\frac{\mathrm{i}\mathbf{q}}{4}\right]\delta \eta\nonumber\\ &&+\left[\frac{\zeta _+}{\nu _2}\left(\mathrm{i}\mathbf{q}\cdot {\mathbf{g}}_{\mathrm{h}}\right)-\frac{\mathbf{q}^2}{4\nu _2}-\nu_1-\frac{\mu\kappa }{\nu _2}{\mathbf{g}}_{\mathrm{h}}^2\right]\delta{\mathbf{g}}_0-\frac{2\mu \kappa }{\nu _2}{\mathbf{g}}_{\mathrm{h}}\left({\mathbf{g}}_{\mathrm{h}}\cdot \delta{\mathbf{g}}_0\right)\nonumber\\ &&+\frac{\zeta _-}{\nu _2}[{\mathbf{g}}_{\mathrm{h}}(\mathrm{i}\mathbf{q}\cdot \delta{\mathbf{g}}_0)-\mathrm{i}\mathbf{q}\left({\mathbf{g}}_{\mathrm{h}}\cdot \delta{\mathbf{g}}_0\right)]. \end{eqnarray} \tag{ 45c }$$

Unlike the canonical model, the grand canonical model conserves the number of single particles, but not the total number of particles. The appropriate base state in this case reads

$$\begin{equation} \rho_s = \mbox{const.}, \end{equation} \tag{ 46a }$$

$$\begin{equation} \rho_c = \rho_c^*(\rho_s) = \frac{\rho_s^2}{1-\rho_s}, \end{equation} \tag{ 46b }$$

$$\begin{equation} {\mathbf{g}} = {\mathbf{g}}_{\mathrm{h}} \in \left\{0,\;g_0=\sqrt{-\frac{\nu_1\nu_2}{\mu\kappa}}\right\}\hat{\mathbf{e}}_g. \end{equation} \tag{ 46c }$$

$$\begin{equation} \rho_c({{\mathbf{x}}},t) = \rho_{c}^* + \delta \rho_c({{\mathbf{x}}}, t), \end{equation} \tag{ 47a }$$

$$\begin{equation} {\mathbf{g}}({{\mathbf{x}}},t) = {\mathbf{g}}_{h} + \delta {\mathbf{g}}({{\mathbf{x}}}, t), \end{equation} \tag{ 47b }$$

$$\begin{equation} \delta \rho_c({\mathbf{x}},t) = \delta \rho_{c}^{0}\,\mathrm{e}^{{\rm st}+\mathrm{i} \mathbf{q}\cdot {\mathbf{x}}}, \end{equation} \tag{ 48a }$$

$$\begin{equation} \delta {\mathbf{g}} ({\mathbf{x}},t)= \delta {\mathbf{g}}_0\,\mathrm{e}^{{\rm st}+\mathrm{i} \mathbf{q}\cdot {\mathbf{x}}}. \end{equation} \tag{ 48b }$$

$$\begin{equation} \fl s\,\delta\rho _c^0=\left(\rho _s-1\right)\delta\rho_c^0-\mathrm{i}\mathbf{q}\cdot \delta{\mathbf{g}}_0, \end{equation} \tag{ 49a }$$

$$\begin{eqnarray} &&\fl s\,\delta{\mathbf{g}}_0=\Bigg[\frac{\partial_{\rho_c}\nu_2\,\mu }{\nu _2^2}\left(\kappa\,{\mathbf{g}}_{\mathrm{h}}^2\,{\mathbf{g}}_{\mathrm{h}}+\left({\mathbf{g}}_{\mathrm{h}}\cdot \mathrm{i}\mathbf{q}\right){\mathbf{g}}_{\mathrm{h}}-\frac{\mathrm{i}\mathbf{q}}{2}{\mathbf{g}}_{\mathrm{h}}^2\mathbf{ }\right) -\partial_{\rho_c}\nu_1{\mathbf{g}}_{\mathrm{h}}-\frac{\mathrm{i}\mathbf{q}}{2}\Bigg]\delta\rho_c^0\nonumber\\ &&+\left[\frac{\zeta _+}{\nu _2}\left(\mathrm{i}\mathbf{q}\cdot {\mathbf{g}}_{\mathrm{h}}\right)-\frac{\mathbf{q}^2}{4\nu _2}-\nu_1-\frac{\mu\kappa }{\nu _2}{\mathbf{g}}_{\mathrm{h}}^2\right]\delta{\mathbf{g}}_0-\frac{2\mu \kappa }{\nu _2}{\mathbf{g}}_{\mathrm{h}}\left({\mathbf{g}}_{\mathrm{h}}\cdot \delta{\mathbf{g}}_0\right)\nonumber\\ &&+\frac{\zeta _-}{\nu _2}\Bigl[{\mathbf{g}}_{\mathrm{h}}(\mathrm{i}\mathbf{q}\cdot \delta{\mathbf{g}}_0)-\mathrm{i}\mathbf{q}\left({\mathbf{g}}_{\mathrm{h}}\cdot \delta{\mathbf{g}}_0\right)\Bigr]. \end{eqnarray} \tag{ 49b }$$

We start by considering the homogeneous and isotropic base states, which was shown to be stable against spatially homogeneous perturbations for ρ < ρ^(c)(σ); cf section 4.2. To assess the stability of this state with respect to perturbations of arbitrary (non-zero) wave vectors in the canonical model, we use the linearized hydrodynamic equations (45a)–(45c) with g_h = 0. The resulting eigenvalue problem is most conveniently expressed in matrix form:

$$\begin{equation} s\left(\begin{array}{@{}c@{}} \delta\rho_0\\ \delta\eta_0\\ \delta g_0 \end{array}\right) =\left(\begin{array}{@{}ccc@{}} 0 & 0 & -\mathrm{i}q\\ 2\rho_{\mathrm{h}}-\eta^*-1 & -(1+\rho_{\mathrm{h}}) & -\mathrm{i}q\\ -\mathrm{i}q/4 & -\mathrm{i}q/4 & -\nu_1-q^2/(4\nu_2)\end{array}\right) \left(\begin{array}{@{}c@{}} \delta\rho_0\\ \delta\eta_0\\ \delta g_0 \end{array}\right). \end{equation} \noindent \tag{ 50 }$$

The corresponding eigenvalue problem for the grand canonical model is found from equations (49a)–(49b), and attains the following form:

$$\begin{equation} s\left(\begin{array}{@{}c} \delta\rho_c^0\\ \delta g_0 \end{array}\right) = \left(\begin{array}{@{}cc@{}} (\rho_s-1) & -\mathrm{i}q\\ -\mathrm{i}q/2 & -\nu_1-q^2/(4\nu_2) \end{array}\right) \left(\begin{array}{@{}c@{}} \delta\rho_c^0\\ \delta g_0 \end{array}\right). \end{equation} \noindent \tag{ 51 }$$

For g_h = 0, (45c) or (49b), respectively, implies q∥δg₀, allowing us to replace the vectors q and δg₀ by their respective magnitudes q and δg₀. We solved both eigenvalue problems numerically for arbitrary wavenumbers q > 0 with the results shown in figure 5. Note that in both models the real parts of all eigenvalues are negative for all wavenumbers q > 0, provided the particle density ρ < ρ^(c). The spatially homogeneous, isotropic state is thus stable against small perturbations with arbitrary wave vectors.

Zoom In Zoom Out Reset image size

For densities ρ > ρ^(c), in contrast, a narrow band of positive eigenvalues emerges in both models, located at wavenumbers q ≪ 1. Equations (50) and (51) evaluated at q = 0 return nothing but the linearized versions of the homogeneous hydrodynamic equations, (34a)–(34c) and (35a)–(35b). To gain new insights, we will therefore examine the limit q → 0 and consider the eigenvalues of the above coefficient matrices to leading order in the wavenumber q.

In this limit of small wavenumbers, the grand canonical coefficient matrix, given in (51), approaches diagonal form and the dynamics of density fluctuations δρ⁰_c and momentum current density fluctuations δg₀ practically decouple. Since ρ_s < 1 (cf figure 3), the first eigenvalue $s_1^{({\scriptsize{\mathrm GC}})}=\rho _s-1+\mathcal {O}(q^2)$ is strictly negative and density fluctuations decay exponentially. The second eigenvalue, $s_2^{({\scriptsize{\mathrm GC}})}=-\nu _1+\mathcal {O}(q^2)$, is positive at small wavenumbers, leading to an instability in the momentum current density against long-wavelength fluctuations.

In the case of the canonical model (50), the coefficient matrix approaches block diagonal form in the limit of small wavenumbers. Again, the dynamics of momentum current density fluctuations δg₀ practically decouples from density fluctuations (δρ₀ and δη₀), with momentum current density fluctuations being amplified by virtue of a positive eigenvalue $s_3^{({\scriptsize{\mathrm C}})}=-\nu _1+\mathcal {O}(q^2)$ at small wavenumbers. In contrast to the grand canonical model, however, particle conservation entails a marginally stable mode s^(C)₁(q = 0) = 0, which turns positive for q ≳ 0: s^(C)₁∝q² (the remaining eigenvalue $s_2^{({\scriptsize{C}})}=-(1+\rho _{\mathrm {h}})+\mathcal {O}(q^2)$ is strictly negative).

To sum up, the study of the linear stability of the homogeneous, isotropic state against spatially inhomogeneous perturbations of arbitrary wave vectors strongly suggests that particle conservation plays a vital role in the context of pattern formation. Both models exhibit spontaneous symmetry breaking by establishing a state of macroscopic collective motion. In the canonical model, in addition, conservation of total particle number entails a marginally stable density mode at q = 0 which is absent in the grand canonical model. This mode, in turn, gives rise to a density instability at small, non-zero wavenumbers, accompanying the spontaneous symmetry breaking event for ρ > ρ^(c). We note, however, that, at this point of the discussions, the existence of a narrow band of unstable modes at small wavenumbers does not allow for any conclusions concerning the structure of the macroscopic density and momentum current density for ρ > ρ^(c). We will address this issue in greater detail in the following section.

Both the canonical and grand canonical models exhibit spontaneous symmetry breaking for overall densities ρ > ρ^(c)(σ). To shed light on the spatial structure of this broken symmetry state, we start from the most simple case of a spatially homogeneous state of collective motion, and examine its stability with respect to wave-like perturbations in the hydrodynamic particle and momentum current densities. Without loss of generality, we assume that the direction of the macroscopic momentum current density coincides with the x-direction and choose ${\mathbf {g}}_{\mathrm {h}} = g_0\,\hat {\mathbf {e}}_x$. The wave vector q of the perturbation fields is assumed to make an angle ψ with the macroscopic momentum current density g_h, yielding ψ = ∠(q,e_x) and q = q (cos(ψ),sin(ψ)) with q = |q|.

The linearized canonical model equations (45a)–(45c) then attain the following form:

$$\begin{equation} \fl s\,\delta \rho_0 = -\mathrm{i} q \cos(\psi)\delta g_{x,0}-\mathrm{i} q \sin(\psi)\delta g_{y,0}, \end{equation} \tag{ 52a }$$

$$\begin{equation} \fl s\,\delta\eta_0 = \left(2\rho_{\mathrm{h}}-\eta^*(\rho_{\mathrm{h}})-1\right)\delta \rho_0 -\left(1+\rho _0\right)\delta\eta_0-\mathrm{i} q \cos(\psi)\delta g_{x,0}-\mathrm{i} q \sin(\psi)\delta g_{y,0}, \end{equation} \tag{ 52b }$$

$$\begin{eqnarray} &&\fl s\,\delta g_{x,0} = \Biggl[\frac{1}{2}\mathrm{i} q \cos(\psi)\left(g_0^2\,\frac{\mu}{\nu _2^2}\partial _{\rho }\nu _2-\frac{1}{2}\right)-g_0\left(\partial _{\rho }\nu_1-\frac{\mu\kappa }{\nu_2^2}g_0^2\partial _{\rho}\nu_2\right)\Biggr]\delta \rho_0\nonumber\\ &&+\Biggl[\frac{1}{2}\mathrm{i} q \cos(\psi)\left( g_0^2\,\frac{\mu}{\nu _2^2}\partial _{\eta}\nu _2-\frac{1}{2}\right)-g_0\left(\partial_{\eta}\nu_1-\frac{\mu\kappa }{\nu_2^2}g_0^2\partial_{\eta}\nu _2\right)\Biggr]\delta\eta_0\nonumber\\ &&+\left(\mathrm{i} q \cos(\psi)\frac{\zeta _+}{\nu _2}g_0-\frac{q^2}{4\nu _2} -\frac{2\mu\kappa}{\nu_2}g_0^2\right)\delta g_{x,0}+\mathrm{i} q\frac{\zeta _-}{\nu _2}\sin(\psi)g_0\delta g_{y,0}, \end{eqnarray} \tag{ 52c }$$

$$\begin{eqnarray} &&\fl s\,\delta g_{y,0} =-\frac{1}{2}\mathrm{i} q \sin(\psi)\left( g_0^2\,\frac{\mu}{\nu _2^2}\partial _{\rho }\nu _2+\frac{1}{2}\right)\delta \rho_0-\frac{1}{2}\mathrm{i} q \sin(\psi)\left( g_0^2\,\frac{\mu}{\nu _2^2}\partial _{\eta}\nu _2+\frac{1}{2}\right)\delta\eta_0\nonumber\\ &&-\mathrm{i} q \sin(\psi)\frac{\zeta _-}{\nu _2}\,g_0\,\, \delta g_{x,0}+\left(\mathrm{i} q \cos(\psi)\frac{\zeta _+}{\nu _2}g_0-\frac{q^2}{4\nu _2}\right)\delta g_{y,0}. \end{eqnarray} \tag{ 52d }$$

$$\begin{equation} \fl s\, \delta\rho_{c}^{0} = \left(\rho_{s}-1\right)\delta\rho_{c}^{0} -\mathrm{i }q \cos(\psi)\delta g_{x,0}-\mathrm{i} q \sin(\psi)\delta g_{y,0}, \end{equation} \tag{ 53a }$$

$$\begin{eqnarray} &&\fl s\, \delta g_{x,0} = \Biggl[\frac{1}{2}\mathrm{i} q \cos(\psi)\left( g_0^2\,\frac{\mu}{\nu _2^2}\partial _{\rho_c}\nu _2-1\right) -g_0\left(\partial_{\rho_c}\nu_1-\frac{\mu\kappa }{\nu_2^2}g_0^2\partial_{\rho_c}\nu _2\right)\Biggr]\delta\rho_{c}^{0}\nonumber\\ &&+\left(\mathrm{i} q \cos(\psi)\frac{\zeta _+}{\nu _2}g_0-\frac{q^2}{4\nu _2} -\frac{2\mu\kappa}{\nu_2}g_0^2\right)\delta g_{x,0}+\mathrm{i} q\frac{\zeta _-}{\nu _2}\sin(\psi)g_0\, \, \delta g_{y,0}, \end{eqnarray} \tag{ 53b }$$

$$\begin{eqnarray} &&\fl s\, \delta g_{y,0} =-\frac{1}{2}\mathrm{i} q \sin(\psi)\left( g_0^2\,\frac{\mu}{\nu _2^2}\partial _{\rho_{c}}\nu _2+1\right)\delta\rho_{c}^{0}-\mathrm{i} q \sin(\psi)\frac{\zeta _-}{\nu _2}\,g_0\,\, \delta g_{x,0}\nonumber\\ &&+\left(\mathrm{i} q \cos(\psi)\frac{\zeta _+}{\nu _2}g_0-\frac{q^2}{4\nu _2}\right)\delta g_{y,0}. \end{eqnarray} \tag{ 53c }$$

$-\nu _1-\frac {\mu \kappa }{\nu _2}\,g_0^2=0$

In the case of the canonical model, we find that the most unstable mode occurs for longitudinal perturbations, i.e. perturbations with wave vectors parallel to the direction of macroscopic motion, q∥g₀ (ψ = 0). Figure 6(a) shows the corresponding eigenvalues as functions of the wavenumber q for a set of density values slightly beyond ρ = ρ^(c). Further inspection of the coupling coefficients in equations (52a)–(52d) reveals that this longitudinal instability only affects the amplitude of g leaving the direction unchanged. For ψ = 0, the dynamics of δg_y,0 decouples and momentum current density fluctuations perpendicular to the direction of macroscopic motion decay exponentially,

$$\begin{equation} \delta g_{y,0} = s_4^{({\scriptsize{\mathrm C}})}\,\delta g_{y,0}, \end{equation} \tag{ 54 }$$

with a rate

$$\begin{equation} \Re\left[s_4^{({\scriptsize{\mathrm C}})}\right] = -\frac{q^2}{4\nu _2} < 0, \end{equation} \tag{ 55 }$$

which approaches zero for q → 0, as expected for a broken symmetry variable. To assess the nature of the instability in greater detail, we calculated the eigenvector corresponding to the most unstable longitudinal mode (evaluated at the most unstable wavenumber). It turns out that this eigenvector has approximately equally large components along the remaining fluctuation amplitudes δg_x,0, δρ₀ and δη₀. This is consistent with our previous findings, indicating that the density mode, which was alluded to in section 5.2 and which turns unstable at ρ = ρ_c, renders the state of homogeneous collective motion unstable to fluctuations of the magnitude of the momentum current density. We further note that this picture is in agreement with previous numerical [31] and analytical [22] results (cf figure 1).

Zoom In Zoom Out Reset image size

The stability regions of the grand canonical model strongly deviate from the above picture. Setting ψ = 0 (longitudinal perturbations), we calculated the largest eigenvalue s^(GC)_max of the linear system of equations (53a)–(53c):

$$\begin{equation} \Re\left[s_{\scriptsize{\mathrm max}}^{({\scriptsize{\mathrm GC}})}\right] = -\frac{(1-\rho_s)\,q^2}{4\left[(1-\rho_s)^2+\rho_s^2\left(\frac{14}{15}+\frac{2}{3}\,\mathrm{e}^{-2\sigma^2}\right)\right]}, \end{equation} \tag{ 56 }$$

which is always negative since ρ_s < 1. In contrast to the canonical model, longitudinal perturbations thus always decay exponentially fast in the grand canonical model. For perturbations in transverse directions, in contrast, a positive eigenvalue can be found for sufficiently low noise levels σ < σ_r, with the fastest growing modes possessing wave vectors q ⊥ g₀. Figure 6(b) shows the eigenvalue of the most unstable modes, which occur for σ = 0. To assess the implications of this instability for the dynamics of the various fluctuation amplitudes, we numerically examined the eigenvector corresponding to the positive eigenvalue, evaluated at the most unstable wavenumber. For densities in the vicinity of the ordering transition, we find that this eigenvector has approximately equal components in both momentum current density fluctuation amplitudes, δg_x,0 and δg_y,0, but an essentially vanishing component along the direction of density fluctuations δρ⁰_c. The corresponding instability can thus be classified as a hybrid shear/splay instability, leaving the spatially homogeneously distributed particle density virtually unaffected.

Three remarks are in order. Firstly, for noise values σ_c(ρ → ∞) > σ > σ_r, the state of homogeneous collective motion becomes linearly stable with respect to arbitrary perturbations, including transverse perturbations. Secondly, even the most unstable eigenvalues 'restabilize' for densities that are in the vicinity of the ordering transition threshold, which is depicted in figure 7. Finally, a restabilization can also be 'observed' for the longitudinal instability in the canonical model. In this case, however, the restabilization occurs for relatively large densities and thus lies outside the range of validity of the linearized equations (52a)–(52d).

Zoom In Zoom Out Reset image size