Uncovering novel phase transitions in dense dry polar active fluids using a lattice Boltzmann method

Symmetry serves a foundational role in all areas of physics today. By identifying the underlying symmetries of a classical many-body system, one can derive the hydrodynamic equations of motion (EOM) that govern the macroscopic dynamics of that system, and, crucially, any other system that respects the same symmetries and conservation laws [1]. In the case of simple (passive) fluids, temporal, rotational, translational, chiral and Galilean symmetries, and mass conservation lead to the Navier–Stokes equations in the hydrodynamic limit [2]. Removing the Galilean symmetry and breaking the conservation of momentum lead instead to the Toner–Tu EOM [3–5], which generically describe dry polar active fluids [6]—systems typically composed of particles that generate non-momentum conserving propulsion through interaction with a fixed frictional medium.

Analysis of a hydrodynamic theory can elucidate the universal behaviour exhibited by all generic systems respecting the prescribed set of symmetries; conversely, any particular many-body system defined by microscopic rules that respect the same set of symmetries can also be used to study the associated universal behaviour in the hydrodynamic limit. An example of the former is the use of dynamic renormalisation group methods to confirm the universal long-time tail phenomenon [7] in thermal fluids first discovered by simulations [8, 9]; and an example of the latter is the use of lattice gas cellular automata to study the Navier–Stokes equations [10]. In dry active fluids, similar lattice gas models have helped to clarify the nature of the order-disorder transition in polar active fluids [11], to elucidate the emergence of diverse dynamical structures [12], and to determine the critical behaviour of motility-induced phase separation (MIPS) [13]. Superseding the lattice gas cellular automata are the celebrated lattice Boltzmann methods (LBM) [14–18], which led to a drastic improvement in the computational efficiency of the simulation requirement.

With their suitability to the almost incompressible regime, LBM have unsurprisingly been employed to simulate the fluid part of suspensions of liquid crystals and active swimmers [19–35]. These examples all describe wet active matter with momentum conserving forces [6], however the development of the LBM for dry polar active fluids is still lacking, which is our focus here (figure 1). Specifically, we present a method to study two-dimensional dry active fluids in a dense regime using a simple modification of the LBM. We then apply the simulation platform to a model that exhibits MIPS and an active model with contact inhibition of locomotion (CIL). The latter is of direct relevance to migrating cells [36, 37]. We recover all known salient features of these two models. In addition, we uncover in the CIL model multiple novel phase transitions: two first order and one potentially critical. We further support our simulation results with an analytical treatment of the hydrodynamic equations obtained via a Chapman–Enskog coarse-graining procedure.

Zoom In Zoom Out Reset image size

Our system involves the evolution of a distribution function ${f}_{i}\left(t,\mathbf{r}\right)$, i ∈ {0, 1, ..., 6}, defined on a two-dimensional triangular lattice (D2Q7 in standard LBM notation [38]) with corresponding lattice vectors e_i given by e₀ = 0, and ${\mathbf{e}}_{i}=\mathrm{cos}\left[\left(i-1\right)\pi /6\right]\hat{\mathbf{x}}+\mathrm{sin}\left[\left(i-1\right)\pi /6\right]\hat{\mathbf{y}}$ for i > 0. The distribution function f_i is related to the local mass density and local velocity of the system as follows:

$$\begin{equation}\rho \left(t,\mathbf{r}\right)=\sum\limits _{i=0}^{6}{f}_{i}\left(t,\mathbf{r}\right),\hspace{25.0pt}\mathbf{u}\left(t,\mathbf{r}\right)=\frac{{\sum }_{i=0}^{6}{f}_{i}\left(t,\mathbf{r}\right)c{\mathbf{e}}_{i}}{\rho \left(t,\mathbf{r}\right)},\end{equation} \tag{ 1 }$$

where c is the ratio of the lattice spacing Δx to the time increment Δt. Typically, scales are chosen so that c = 1, and |u| ≪ c is required for stable solutions (see appendix A.1).

We then evolve the discretized Boltzmann equation using the Bhatnagar–Gross–Krook collision operator [39], albeit with fluctuations included, as follows:

$$\begin{equation}{f}_{i}\left(t+{\Delta}t,\mathbf{r}+c{\mathbf{e}}_{i}{\Delta}t\right)-{f}_{i}\left(t,\mathbf{r}\right)=-\frac{1}{\tau }\left[{f}_{i}\left(t,\mathbf{r}\right)-{f}_{i}^{\text{SS}}\left(t,\mathbf{r}\right)\right]+{\eta }_{i}\left(t,\mathbf{r}\right),\end{equation} \tag{ 2 }$$

where τ is a relaxation parameter, ${f}_{i}^{\text{SS}}$ is the steady-state distribution, and

$$\begin{equation}{\eta }_{i}\left(t,\mathbf{r}\right)={\tilde {\eta }}_{i}\left(t,\mathbf{r}\right)-\frac{1}{7}\sum\limits _{i=0}^{6}{\tilde {\eta }}_{i}\left(t,\mathbf{r}\right),\end{equation} \tag{ 3 }$$

where ${\tilde {\eta }}_{i}$ are temporally and spatially uncorrelated random variables uniformly distributed in the interval [−σ, σ]. This particular choice of noise alters the local velocity, but equation (3) ensures the density is conserved. Note that the fluctuations do not conserve momentum as we are considering dry active fluids. Equation (2) is evolved in two parts, streaming and redistribution, as illustrated in figure 1.

In standard LBM, ${f}_{i}^{\text{SS}}$ is constructed so as to preserve the moments ρ, u at each site, accomplished in a D2Q7 model by

$$\begin{equation}{f}_{i}^{\text{SS}}={w}_{i}\rho \left[1+4\frac{{\mathbf{e}}_{i}\cdot {\mathbf{u}}^{{\ast}}}{c}+8\frac{{\left({\mathbf{e}}_{i}\cdot {\mathbf{u}}^{{\ast}}\right)}^{2}}{{c}^{2}}-2\frac{\vert {\mathbf{u}}^{{\ast}}{\vert }^{2}}{{c}^{2}}\right],\end{equation} \tag{ 4 }$$

with u* = u, where w_i are lattice specific weights: w₀ = 1/2 and w_i≠0 = 1/12. However, in breaking the conservation of local momentum, the restriction on the form of the steady-state distribution ${f}_{i}^{\text{SS}}$ is drastically relaxed. We make a simple modification: maintaining the form of ${f}_{i}^{\text{SS}}$ in equation (4), but taking u* to be a function of ρ and u, that respects the symmetries of our system of interest. This update rule conserves mass density, but does not conserve momentum nor energy, thus rendering the system in the realm of the hydrodynamic theory of Toner and Tu, rather than Navier–Stokes (see appendix B). This method can produce a variety of behaviours associated with dry active matter, depending on the choice of u*.

To demonstrate the versatility of our method, we first employ it to produce MIPS. MIPS can occur when self-driven particles that maintain their direction of propulsion with some persistency, interact through short-range repulsive forces [13, 40–46]. The particles congregate in high density clusters, resembling a liquid–gas phase separation.

The phase separation arises from two related mechanisms. Firstly, it has been observed that active particles tend to congregate in regions where their speeds are slow [47]. This is well demonstrated using light sensitive strains of E. coli bacteria, with heterogeneous distributions of light leading to matching density profiles [48]. Secondly, in systems with repulsive forces, the converse is also true; in areas with high densities, the speeds slow down. The exact reason for the slowdown differs depending on the system of interest. Run-and-tumble bacteria respond to chemical secretions by their neighbours, which affects their own behaviour [49, 50]. Active Brownian particles on the other hand do not alter their behaviour in the presence of other particles. Rather, at high densities, they experience many collisions in a short time, disrupting their motion. If their trajectory is averaged over a short interval, they are seen to effectively slow down at higher densities [41, 42].

These two mechanisms combine to form a positive feedback loop, and an instability occurs. Particles in a region of high density slow down, resulting in more particles congregating in that region, increasing the density further. The cluster grows and merges with others until eventually a stable equilibrium in the system is achieved.

To generate MIPS using our method we take

$$\begin{equation}{\mathbf{u}}^{{\ast}}=\begin{cases}C\left(1-\frac{\rho }{{\rho }^{{\ast}}}\right)\nabla \rho \enspace ,\quad \hfill & \mathrm{i}\mathrm{f}\;\rho {< }{\rho }^{{\ast}},\hfill \\ 0\enspace ,\quad \hfill & \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e},\hfill \end{cases}\end{equation} \tag{ 5 }$$

where C is a positive constant, dictating the sensitivity of the system to density gradients. The parameter ρ* acts as an upper bound for the local density as permitted by the system; above this there is no net flow, regardless of density gradients, hence the density cannot rise further. Since ∇ρ could in principle become very large and any physical system will not have a divergent speed, we impose further an upper bound on |u*|. Explicitly, we have

$$\begin{equation}\vert {\mathbf{u}}^{{\ast}}\vert {\mapsto}\mathrm{min}\left(\vert {\mathbf{u}}^{{\ast}}\vert ,{U}_{\mathrm{max}}\right).\end{equation} \tag{ 6 }$$

Standard LBM do not typically feature gradients. Indeed, apart from the streaming step, all of the calculations and operations are strictly localised to each lattice site. To calculate a gradient we must include information from neighbouring nodes. We do so by incorporating a finite difference procedure into the routine [51]. It is important to use a finite difference stencil that is isotropic when projected to the lattice. A suitable choice is given by

$$\begin{equation}\nabla \rho =a\left[\left({\rho }_{1}-{\rho }_{4}\right){\mathbf{e}}_{1}+\left({\rho }_{2}-{\rho }_{5}\right){\mathbf{e}}_{2}+\left({\rho }_{3}-{\rho }_{6}\right){\mathbf{e}}_{3}\right],\end{equation} \tag{ 7 }$$

for a = 1/(3Δx), where in this notation ρ_i = ρ(t, r + Δx e_i ) [50]. The resulting stencils are

$$\begin{equation}\frac{\partial {\rho }_{0}}{\partial x}=\frac{2\left({\rho }_{1}-{\rho }_{4}\right)+\left({\rho }_{2}+{\rho }_{6}\right)-\left({\rho }_{3}+{\rho }_{5}\right)}{6{\Delta}x},\end{equation} \tag{ 8 }$$

$$\begin{equation}\frac{\partial {\rho }_{0}}{\partial y}=\frac{\left({\rho }_{2}+{\rho }_{3}\right)-\left({\rho }_{5}+{\rho }_{6}\right)}{4{\Delta}y}.\end{equation} \tag{ 9 }$$

We motivate the expression in equation (5) firstly by intuition, and secondly through comparison with other studies. For our system to experience MIPS the speed of particles must decrease in regions of high density. The speeds of individual particles do not feature in our system, as u* is a macroscopic variable. However, a consequence of all particles' speeds decreasing is that the net speed |u*| will also decrease. This is clearly achieved by the term $\left(1-\rho /{\rho }^{{\ast}}\right)$, with the net speed hitting zero altogether for ρ ⩾ ρ*. Physical systems of hard particles have an upper limit as to how tightly they can be packed, hence it is natural to have a maximum density ρ*. There is no alignment mechanism in a MIPS model and directions of particles are typically random and isotropic. This raises a question as to how we can have a non-zero net flow at all. The answer comes in the form of a statistical average. Suppose at some location in our system there is a non-zero density gradient. Nearby particles travelling in the same direction as the gradient will be entering a region of higher density, therefore they will experience more collisions and thus be impeded more than particles moving in other directions. As the particles maintain their current direction with some persistency, after a short time the impeded particles will still be in place, pointing towards the density gradient, while the non-impeded particles will have left the vicinity entirely. The average velocity of the particles will therefore align with the density gradient, and as a result there is a net flow towards regions of higher density. This is also consistent with the observation that particles will congregate in areas of lower speeds or, equivalently, higher densities; we expect there to be a flow towards those regions.

We can also arrive at the expression in equation (5) through comparison with previous studies in the field. In the notation of [40]: for a system of non-interacting active particles propagating with fixed speed v(r), and orientational relaxation time τ(r), where v and τ can vary in space, it has been shown that the macroscopic density can be approximated by the following

$$\begin{equation}\frac{\partial \rho }{\partial t}=-\nabla \cdot \mathbf{J},\end{equation} \tag{ 10 }$$

$$\begin{equation}\mathbf{J}=-D\nabla \rho +\mathbf{\text{V}}\rho +\sqrt{2D\rho }{\Lambda},\end{equation} \tag{ 11 }$$

where V is a drift velocity, D is a diffusivity, and Λ is a white noise [40]. With negligible translational diffusivity these can be expressed as

$$\begin{equation}D=\frac{{v}^{2}\tau }{d},\end{equation} \tag{ 12 }$$

$$\begin{equation}\frac{\mathbf{V}}{D}=-\nabla \enspace \mathrm{log}\enspace v,\end{equation} \tag{ 13 }$$

where d is the dimension. Generalizing this system to include interaction effects, allows v and τ to now depend on the local density ρ, as well as r. By extension so do D and V. For simplicity we assume τ to be constant and v to depend only on ρ. Numerical studies of active Brownian particles have shown that the effective speed has an approximately linear dependence on the density [41, 42, 45, 46]. Hence we adopt

$$\begin{equation}v\left(\rho \right)={v}_{0}\left(1-\frac{\rho }{{\rho }^{{\ast}}}\right).\end{equation} \tag{ 14 }$$

The drift velocity V becomes

$$\begin{equation}\mathbf{V}=\frac{\tau {v}_{0}^{2}}{{\rho }^{{\ast}}d}\left(1-\frac{\rho }{{\rho }^{{\ast}}}\right)\nabla \rho .\end{equation} \tag{ 15 }$$

Finally we observe that equation (10) is a continuity equation and equate J ≡ ρ u in our system. We take the drift velocity to be the steady state that u tends towards, that is V ≡ u*. The term D∇ρ features in our model as the 'pressure' term in our hydrodynamic equations (see appendix B), and we further incorporate fluctuations in the system. Our system is very similar to that presented in equations (10) and (11), albeit our velocity u is dynamic, while J has no time dependence. In that sense, equation (11) is the overdamped version of our model.

We ran our simulation in a rectangular geometry with periodic boundary conditions. For simulation details, including parameter values, see appendix C. Figure 2 shows the results from our MIPS simulations. Figures 2(a)–(d) are still frames from our simulations displayed using the visualization scheme described in appendix E. In panels (b)–(d) MIPS has occurred. The multi-coloured nature of the dense region indicates there is no velocity alignment between neighbouring lattice sites. Therefore, collectively, the dense regions have effectively zero velocity.

Zoom In Zoom Out Reset image size

Figure 2(f) displays a phase diagram, where the orange squares are the mean gas densities, and the blue circles are the mean liquid densities measured in systems that exhibited phase separation. Error bars for the standard error of the mean are smaller than the marker size. The data points appear to trace out a convex hull, the uppermost limits of which we have estimated with the black dashed line. Any system whose mean density and noise strength, $\left(\bar{\rho },\sigma \right)$, fall within the yellow region will experience phase separation, with the density in their liquid and gas phases given by the values on the blue and orange curves for that σ, respectively.

These results are consistent with previous studies in MIPS [13]. As such, we have demonstrated that our LBM is adaptable, and could be used to study MIPS in more detail.

The remainder of this paper is concerned with a class of active systems that exhibit collective motion, for which we take

$$\begin{equation}{\mathbf{u}}^{{\ast}}=\frac{{U}_{0}}{\vert \mathbf{u}\vert }\mathbf{u},\end{equation} \tag{ 16 }$$

where U₀ is of the form

$$\begin{equation}{U}_{0}\left(\rho \right)=\begin{cases}-A{\rho }^{2}+B\quad \hfill & \mathrm{i}\mathrm{f}\;B{\geqslant}A{\rho }^{2},\hfill \\ 0\quad \hfill & \mathrm{o}\mathrm{t}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{w}\mathrm{i}\mathrm{s}\mathrm{e}.\hfill \end{cases}\end{equation} \tag{ 17 }$$

In the above, A, B are positive constants. We expect this to reproduce the dynamics of active polar matter because it mimics the key behaviour that local particles' movements tend to align. As ${f}_{i}^{\text{SS}}$ is defined at each lattice site, the alignment interactions are short-ranged. In addition, we assume here that collective motility can be impeded when the density gets too high or even stop altogether, which can be viewed as a form of CIL, as observed in cell monolayers [52–54] and traffic jams [55, 56]. By using a quadratic expression, we assume that contact inhibition effects are minimal for low densities, and only become relevant at higher values. We will show in section 7 that incorporating fluctuations will generally render the effective 'speed' function U₀ a much more complicated function of ρ.

As for our MIPS model, we simulate our active fluid model in a rectangular geometry with periodic boundary conditions (figure 3(d)). This ensures that any ensuing collective motion aligns itself along the x-axis, which helps to reduce noise in our data collection. We use the average speed of the whole system as an order parameter φ, given by

$$\begin{equation}\varphi =\left\vert \frac{{\sum }_{\mathbf{r}\in \mathrm{l}\mathrm{a}\mathrm{t}\mathrm{t}\mathrm{i}\mathrm{c}\mathrm{e}}\rho \left(\mathbf{r}\right)\mathbf{u}\left(\mathbf{r}\right)}{{\sum }_{\mathbf{r}\in \mathrm{l}\mathrm{a}\mathrm{t}\mathrm{t}\mathrm{i}\mathrm{c}\mathrm{e}}\rho \left(\mathbf{r}\right)}\right\vert .\end{equation} \tag{ 18 }$$

Necessarily 0 ⩽ φ ⩽ B. Figure 3(a) shows how φ behaves with varying average density $\bar{\rho }$. Two curves have been traced by slowly increasing (blue) and decreasing (red) $\bar{\rho }$. The two curves coincide with the exception of four regions, in which there exist hysteresis loops indicating the existence of meta-stable states. We will use these hysteresis loops as proxies for the locations of first-order phase transitions that separate distinct phases. At low densities noise dominates, and the system is in the homogeneous and disordered (HD) phase. As the density increases, alignment effects take over and dense locally aligned bands spontaneously develop, which propagate through a low density disordered region. We refer to this as the banding (B) regime, which corresponds to the coexistence of the HD phase in the dilute regions and a homogeneous, ordered (HO) phase in the dense regions. Increasing the density further causes the dilute HD phase to vanish, resulting in a single HO phase. The first HD-B transition is easily recognised from previous flocking studies [11, 57–61] which found it to be first order. The second B–HO transition is less studied, likely because higher density molecular dynamics simulations are more difficult to compute, but has recently been shown to also be a first order transition [61]. Our results are in perfect agreement with both of these as evidenced by the aforementioned hysteresis loops. This lends credence to the conclusion that the two remaining hysteresis loops at higher densities also indicate first-order phase transitions, in agreement with a recent experimental study concurrent to ours [62]. These transitions separate the HO region from another HD region, again through a coexistence regime similar to the B regime (figure 3(d)), but with the properties inverted—the high density regions are disordered, while the low density regions are aligned and moving. In the steady state, as material flows into a dense band it comes to a halt, extending the band on that side, while on the other side the band decays as the material flows away. The net result is that the dense band, although locally stationary, appears to propagate against the direction of flow, hence we call it reverse banding (RB). In other words, the dense band moves due to condensation on one side and evaporation on the other. In active matter, this phenomenon was first reported in a collective cell migration model [56], and likened to density waves observed in traffic flow [55]. Interestingly, the same mechanism also underlies the directed motion of condensed drops under the influence of a chemical gradient [63], as well as the formation of a drop lattice in a chemical reaction-controlled phase separated drop system [64].

Zoom In Zoom Out Reset image size

The RB regime also bears striking resemblance to the glider structure discovered in [12], which also moves backward due to a similar differential absorption and evaporation mechanism. However, we believe that the RB regime discussed is fundamentally distinct from the glider for the following two reasons: (i) the glider in [12] does not form a band but is triangular in shape; and (ii) a glider consists of two halves, each with a distinct mean orientation. That is, the mean orientation differs within the glider depending on which region one focuses on, while the reverse band discussed here is ordered uniquely and homogeneously within the band.

To further elucidate the phase behaviour of the system, we now vary the contact inhibition parameter A as well as $\bar{\rho }$ (figure 4(a)). When A = 0, we recover the expected phenomenology of a typical polar active model, such as the Vicsek model [65], in that the system transitions from the HD phase to the B regime followed by the HO phase as density increases. However, for A > 0, we see that the HO phase will transition back to the HD phase via the RB regime as density increases further. Most interestingly, as A increases, all four first-order phase transition lines converge to a single point, indicated by a star in figure 4(a), which is a putative critical transition. A snapshot of this critical behaviour is shown in figure 4(b). We will discuss this critical point further in our hydrodynamic analysis. In addition to the four regimes already discussed we recognise a fifth somewhat trivial phase of no motion (indicated by red asterisks), where |u| = 0 since $-A{\bar{\rho }}^{2}+B{\leqslant}0$.

Zoom In Zoom Out Reset image size

Finally, we note that although we focus on the phase separations within the ordered phase, the interesting interplays between MIPS and the polar phase have also been discussed in various recent works [66–73], which are complimentary to this study.

In the hydrodynamic limit, our model is necessarily described by the Toner–Tu equations due to our imposed symmetries, and we confirm this analytically using the Chapman–Enskog expansion method in appendix B. We now describe our simulation findings analytically by performing a linear stability analysis on our derived version of the Toner–Tu equations:

$$\begin{equation}{\partial }_{t}\rho +\nabla \cdot \left(\rho \mathbf{v}\right)=0,\end{equation} \tag{ 19 }$$

$$\begin{equation}\rho \left({\partial }_{t}\mathbf{v}+\mathbf{v}\cdot \nabla \mathbf{v}\right)=-\nabla P+\nu \nabla \cdot \rho \left(\nabla \mathbf{v}+{\left(\nabla \mathbf{v}\right)}^{T}\right)+\rho \left[\alpha \left(\rho \right)-\beta \vert \mathbf{v}{\vert }^{2}\right]\mathbf{v},\end{equation} \tag{ 20 }$$

where $P={c}_{\mathrm{s}}^{2}\rho$.

Since the bands in both the B and RB regimes are perpendicular to the direction of collective motion in the dense and dilute regions, respectively, the most unstable wave vectors in the HO state are expected to be parallel to the direction of motion. Therefore, we need only consider a one-dimensional system, where the coordinate x is the direction of the collective motion. Expanding about a homogeneous collective motion state with ρ = ρ₀ + δρ exp(st − ikx), $v=\sqrt{{\alpha }_{0}/\beta }+\delta v\enspace \mathrm{exp}\left(st-\text{i}kx\right)$ and α = α₀ + α₁ δρ, we have to linear order:

$$\begin{equation}\left(s-\text{i}k\sqrt{\frac{{\alpha }_{0}}{\beta }}\right)\delta \rho =\text{i}k{\rho }_{0}\delta v,\end{equation} \tag{ 21 }$$

$$\begin{equation}{\rho }_{0}\left(s-\text{i}k\sqrt{\frac{{\alpha }_{0}}{\beta }}+2{\alpha }_{0}+2\nu {k}^{2}\right)\delta v=\left(\text{i}k{c}_{\mathrm{s}}^{2}+{\rho }_{0}{\alpha }_{1}\sqrt{\frac{{\alpha }_{0}}{\beta }}\right)\delta \rho .\end{equation} \tag{ 22 }$$

Note that for a stable system β > 0 necessarily, thus α₀ > 0 implies collective motion, with α₀ ⩽ 0 giving a stationary state.

Solving for s and focussing on the real part, which determines the stability of the system, we have

$$\begin{equation}\mathrm{R}\mathrm{e}\left[s\right]=\begin{cases}-2{\alpha }_{0}+\mathcal{O}\left({k}^{2}\right),\quad \hfill \\ \frac{{k}^{2}}{8{\alpha }_{0}^{2}\beta }\left[{\alpha }_{1}^{2}{\rho }_{0}^{2}-4{\alpha }_{0}{c}_{\mathrm{s}}^{2}\beta \right]+\mathcal{O}\left({k}^{3}\right),\quad \hfill \end{cases}\end{equation} \tag{ 23 }$$

where Re[s] > 0 implies instability. The first eigenvalue −2α₀ corresponds to the fast relaxation that occurs when the velocity deviates from the mean field value $\sqrt{{\alpha }_{0}/\beta }$ in the absence of spatial variations. The second eigenvalue quantifies when the instability sets in. In particular, if α₀ → 0 and α₁ ≠ 0, the system is always unstable. In the previously mentioned studies, [11, 58–61, 65], α₁ is typically assumed to be positive, so that alignment interactions increase with density, and as α goes from negative to positive the system will necessarily become unstable, resulting in the banding instability. This corresponds to the banding (B) regime in our system, since a higher density reduces the effects of the noise, and as such alignment increases with density. However, it is also clear that the instability condition does not depend on the sign of α₁. Hence, if α goes from positive to negative due to contact inhibition, the system will again be unstable. This is exactly what occurs in our simulation and corresponds to the reverse-banding (RB) regime in our phase diagram. In our simulation we found that all four first-order phase transition lines converge at a critical point. With our hydrodynamic argument, we can see that this critical point is reached by setting both α₀ and α₁ to zero, which is achieved by fine-tuning $\bar{\rho }$ and A in our simulation. Furthermore, as detailed in the following section, by capturing the effect of fluctuations at the one-loop level, we can obtain an expression for α(ρ) in our system. The blue dashed curve in figure 4(a) was obtained by fitting the curve ${\alpha }_{\mathrm{R}}\left(\bar{\rho }\right)=0$ to pass through the critical point. The resulting curve bisects both the B and RB regions, which is in good agreement with our stability analysis; close to this line α₀ ≈ 0 and α₁ ≠ 0, hence instability, with the exception of close to the critical point where α₁ ≈ 0.

The previous analysis did not incorporate noise into the picture, we will do so now and demonstrate how, when the noise is included, the renormalized coefficient α will naturally depend on the density ρ as well. Here, we aim mainly to connect our hydrodynamic EOM derived from our Chapman–Enskog expansion (see appendix B) to the results of the linear stability analysis, and so will not perform a full renormalisation procedure on the EOM. Specifically, we will show how the β term, together with the noise term with strength σ, renormalise α around the critical region when α ≈ 0.

Firstly, we must understand the macroscopic effect of the noise. Recall that

$$\begin{equation}\sum\limits _{i=0}^{7}{\eta }_{i}=0,\end{equation} \tag{ 24 }$$

hence there is no contribution to the continuity equation. Let us define

$$\begin{equation}\mathbf{h}=\sum\limits _{i=0}^{7}{\eta }_{i}c{\mathbf{e}}_{i}.\end{equation} \tag{ 25 }$$

From the definition of η_i it follows that

$$\begin{equation}\left\langle {h}_{x}\right\rangle =\left\langle {h}_{y}\right\rangle =\left\langle {h}_{x}{h}_{y}\right\rangle =0,\end{equation} \tag{ 26 }$$

$$\begin{equation}\left\langle {h}_{x}\left(\mathbf{x},t\right){h}_{x}\left({\mathbf{x}}^{\prime },{t}^{\prime }\right)\right\rangle =\left\langle {h}_{y}\left(\mathbf{x},t\right){h}_{y}\left({\mathbf{x}}^{\prime },{t}^{\prime }\right)\right\rangle ={c}^{2}{\sigma }^{2}\delta \left(\mathbf{x}-{\mathbf{x}}^{\prime }\right)\delta \left(t-{t}^{\prime }\right).\end{equation} \tag{ 27 }$$

Our derived equation (B.17) now becomes

$$\begin{equation}\rho \left({\partial }_{t}\mathbf{v}+\mathbf{v}\cdot \nabla \mathbf{v}\right)=-\nabla P+\nabla \cdot \rho \nu \left(\nabla \mathbf{v}+{\left(\nabla \mathbf{v}\right)}^{\mathrm{T}}\right)+{\mathbf{F}}_{\text{CIL}}+\mathbf{h}.\end{equation} \tag{ 28 }$$

The calculation mirrors almost exactly the corresponding calculation in the study of incompressible active polar fluids at criticality [74], except that the projection operator is absent here since our system is not incompressible. Specifically, the shift in α is given by

$$\begin{equation}\delta \alpha =-\beta \int \frac{{\mathrm{d}}^{d}\mathbf{q}}{{\left(2\pi \right)}^{d}}\frac{\mathrm{d}{\Omega}}{2\pi }\frac{\left\langle \vert \mathbf{f}{\vert }^{2}\right\rangle }{{\left(\mu {q}^{2}+\alpha \right)}^{2}+{{\Omega}}^{2}}\end{equation} \tag{ 29 }$$

$$\begin{equation}=-\beta \frac{{c}^{2}{\sigma }^{2}}{{\rho }^{2}}\int \frac{{\mathrm{d}}^{d}\mathbf{q}}{{\left(2\pi \right)}^{d}}\frac{1}{\mu {q}^{2}+\alpha }\end{equation} \tag{ 30 }$$

$$\begin{equation}=-\frac{{c}^{2}{\sigma }^{2}\beta }{{\rho }^{2}\left(\mu {{\Lambda}}^{2}+\alpha \right)}\frac{{S}_{d}\mathrm{d}\ell }{{\left(2\pi \right)}^{d}},\end{equation} \tag{ 31 }$$

where $\mathbf{f}=\frac{\mathbf{h}}{\rho }$, S_d ≡ 2π^d/2/Γ(d/2) is the surface area of a unit sphere in d dimensions, Λ⁻¹ corresponds to the small length scale cutoff, and Λ dℓ corresponds to the width of the momentum shell being integrated over in the above coarse-graining process.

Adding this shift to α modifies it to, to $\mathcal{O}\left(\alpha \right)$:

$$\begin{equation}{\alpha }_{\mathrm{R}}=\alpha \left(1-\frac{\sigma }{{\rho }^{2}{U}_{0}^{2}\mu {{\Lambda}}^{2}}\frac{{S}_{d}\enspace \mathrm{d}\ell }{{\left(2\pi \right)}^{d}}\right)=\alpha \left(1-{C}_{1}\frac{1}{{\rho }^{2}{U}_{0}^{2}}\right),\end{equation} \tag{ 32 }$$

where we have used the expression for β in (B.22).

The above expression indicates that the critical transition occurs when U₀ and ρ are fine tuned such that α_R = 0 and ∂α_R/∂ρ = 0, which can be achieved by fine tuning $\bar{\rho }$ and A in the expression of U₀. These are given by

$$\begin{equation}\bar{\rho }=\frac{3\sqrt{{C}_{1}}}{2B},\hspace{25.0pt}A=\frac{4{B}^{3}}{27{C}_{1}}.\end{equation} \tag{ 33 }$$

In our simulations we found ${\bar{\rho }}_{\mathrm{c}}\approx 0.18$, A_c ≈ 2.75. These values give B ≈ 0.27, C₁ ≈ 0.001, which is fairly consistent with our actual value of B = 0.3. In fact we can go further. In figure 4(a) and (b) the blue dashed line depicts the full contour for α_R = 0 using these values. For parameter values below this curve α > 0 and we have collective motion, with α < 0 for parameter values above. By design this line cuts through our critical point (star), however it also bisects the banding and reverse banding regions, in good agreement with our stability analysis; in the vicinity of this curve α₀ ≈ 0 and α₁ ≠ 0 except at the critical point.

We are aware of critical phenomena analysed, either analytically or computationally, in the following active matter systems: (1) self-propelled particles with long-ranged metric-free alignment interactions [75, 76], (2) active Lévy matter [77, 78], (3) non-equilibrium diffusive particles with long-ranged dynamics [79], (4) incompressible active fluids [74], (5) critical MIPS [13, 80, 81], and (6) self-propelled particles with velocity reversals and alignment interactions [82, 83]. Our critical system is different from system (1) to (3) because ours does not have long-range interactions or dynamics, from (4) because our velocity field is not divergence free, from (5) because critical MIPS happens in the disordered phase while criticality in our system occurs at the onset of the ordered phase, from (6) because our system exhibits collective motion while system (6) exhibits only orientational order, but no collective motion. Therefore, we believe that the critical behaviour uncovered here is likely to be novel.

Through symmetry considerations, we have generalized LBM to simulate dry active fluid systems and uncovered novel phase transitions. Looking ahead, we believe that the drastic improvement on computational efficiency will help solve interesting open questions in dry active fluids, such as the dynamic scaling in the incompressible limit [84], and universal behaviour in chiral active fluids [85–87]. Furthermore, the fact that LBM are adaptable to complex geometries renders them suitable for investigating interfacial instabilities in a multiphase active fluid system, such as in the context of tissue regeneration [88–90]. Finally, we believe our method can be straightforwardly extended to higher dimensions.

The data that support the findings of this study are available upon reasonable request from the authors.

DN was supported by the EPSRC Centre for Doctoral Training in Fluid Dynamics Across Scales (Grant EP/L016230/1).

A.1. Conventional lattice Boltzmann methods

LBM refers to solving the Boltzmann equation on a lattice. The Boltzmann equation itself concerns the temporal evolution of a distribution function f(t, r, ξ ), which corresponds to the mass density of matter moving with velocity ξ at location r and time t.

The typical variables of interest are

$$\begin{equation}\rho \left(t,\mathbf{r}\right)=\int {\mathrm{d}}^{d}\boldsymbol{\xi }\enspace f\left(t,\mathbf{r},\boldsymbol{\xi }\right),\end{equation} \tag{ A.1 }$$

$$\begin{equation}\rho \left(t,\mathbf{r}\right)\mathbf{u}\left(t,\mathbf{r}\right)=\int {\mathrm{d}}^{d}\boldsymbol{\xi }\enspace f\left(t,\mathbf{r},\boldsymbol{\xi }\right)\boldsymbol{\xi },\end{equation} \tag{ A.2 }$$

where d is the spatial dimension, ρ is the mass density field, and u is the velocity field.

Discretising ξ into a finite number of velocity vectors, directed towards neighbouring spatial nodes, naturally constructs a lattice in the domain. There are many suitable choices of lattice, but they must contain enough degrees of freedom to preserve the moments of the system in question according to the desired level of accuracy [10, 91]. The DdQq notation is conventional for describing the lattice structure of the domain, where d is again the dimension, and q is the number of nodes each is connected to, including itself. Typically a 'zero' vector is included in the discrete directional space to help balance higher moments. Here, we will focus on a D2Q7 lattice, which is a triangular lattice in two dimensions with a discretized distribution function f_i (t, r), where i = 0, 1, ..., 6. The corresponding vectors are e₀ = 0, and ${\mathbf{e}}_{i}=\mathrm{cos}\left[\left(i-1\right)\pi /6\right]\hat{\mathbf{x}}+\mathrm{sin}\left[\left(i-1\right)\pi /6\right]\hat{\mathbf{y}}$ for i > 0. The local mass density and velocity are thus given by

$$\begin{equation}\rho =\sum\limits _{i=0}^{6}{f}_{i},\hspace{25.0pt}\mathbf{u}=\frac{1}{\rho }\sum\limits _{i=0}^{6}{f}_{i}c{\mathbf{e}}_{i},\end{equation} \tag{ A.3 }$$

where c is the ratio of the lattice spacing Δx to the time increment Δt. For accurate solutions velocities must obey |u| ≪ c, which is known as the low Mach number assumption. Typically scales in the LBM scheme are chosen such that c = 1.

The discretized Boltzmann equation can then be evolved using the Bhatnagar–Gross–Krook collision operator [39]

$$\begin{equation}{f}_{i}\left(t+{\Delta}t,\mathbf{r}+c{\mathbf{e}}_{i}{\Delta}t\right)-{f}_{i}\left(t,\mathbf{r}\right)=-\frac{1}{\tau }\left[{f}_{i}\left(t,\mathbf{r}\right)-{f}_{i}^{\mathrm{e}\mathrm{q}}\left(t,\mathbf{r}\right)\right],\end{equation} \tag{ A.4 }$$

where ${f}_{i}^{\mathrm{e}\mathrm{q}}\left(t,\mathbf{r}\right)$ is a lattice dependent steady-state (equilibrium) distribution and τ is a relaxation parameter. When simulating the Navier Stokes equations, τ is linearly proportional to the viscosity. The evolution of this equation is typically split into two processes—the streaming step and the collision step, the latter of which is effectively a local redistribution (see figure 1). The design of the steady-state distribution is highly flexible and the only requirement is that it respects the symmetries and conservation laws imposed on the system. For the Navier–Stokes equation, ${f}_{i}^{\mathrm{e}\mathrm{q}}$ is obtained by discretising the Boltzmann distribution for velocities at thermal equilibrium [92] and in a D2Q7 framework is given by

$$\begin{equation}{f}_{i}^{\mathrm{e}\mathrm{q}}={w}_{i}\rho \left[1+4\frac{{\mathbf{e}}_{i}\cdot \mathbf{u}}{c}+8\frac{{\left({\mathbf{e}}_{i}\cdot \mathbf{u}\right)}^{2}}{{c}^{2}}-2\frac{\vert \mathbf{u}{\vert }^{2}}{{c}^{2}}\right],\end{equation} \tag{ A.5 }$$

where w_i are lattice specific weights: w₀ = 1/2 and w_i≠0 = 1/12. This discretisation is possible through the low Mach number assumption |u| ≪ c.

A.2. Summary of modified lattice Boltzmann model algorithm

A.3. Negative correction

In this section we will demonstrate that our LBM simulates a version of the Toner–Tu equations. By expressing our active LBM in the form of a conventional LBM with a body force term, we can use the results of [93] to derive the macroscopic equations.

Let us introduce the notation

$$\begin{equation}{\xi }_{i}\left(\rho ,\mathbf{u}\right)={w}_{i}\rho \left[1+4\frac{{\mathbf{e}}_{i}\cdot \mathbf{u}}{c}+8\frac{{\left({\mathbf{e}}_{i}\cdot \mathbf{u}\right)}^{2}}{{c}^{2}}-2\frac{\vert \mathbf{u}{\vert }^{2}}{{c}^{2}}\right].\end{equation} \tag{ B.1 }$$

Hence in the conventional LBM ${f}_{i}^{\mathrm{e}\mathrm{q}}={\xi }_{i}\left(\rho ,\mathbf{u}\right)$ and in our active model ${f}_{i}^{\text{SS}}={\xi }_{i}\left(\rho ,{\mathbf{u}}^{{\ast}}\right)$.

The evolution equation becomes

$$\begin{equation}{f}_{i}\left(t+{\Delta}t,\mathbf{r}+{\mathbf{e}}_{i}{\Delta}t\right)-{f}_{i}\left(t,\mathbf{r}\right)=-\frac{1}{\tau }\left[{f}_{i}-{\xi }_{i}\left(\rho ,\mathbf{u}\right)\right]+{F}_{i},\end{equation} \tag{ B.2 }$$

where

$$\begin{equation}{F}_{i}=\frac{1}{\tau }\left[{\xi }_{i}\left(\rho ,{\mathbf{u}}^{{\ast}}\right)-{\xi }_{i}\left(\rho ,\mathbf{u}\right)\right].\end{equation} \tag{ B.3 }$$

Note that

$$\begin{equation}\sum\limits _{i=0}\;{F}_{i}=0,\quad \sum\limits _{i=0}\;{F}_{i}c{\mathbf{e}}_{i}=\rho \frac{{\mathbf{u}}^{{\ast}}-\mathbf{u}}{\tau }={\Delta}t\enspace \mathbf{F}.\end{equation} \tag{ B.4 }$$

As in [93], we define the 'actual velocity' according to

$$\begin{equation}\mathbf{v}=\mathbf{u}+\frac{{\Delta}t\mathbf{F}}{2\rho }=\mathbf{u}-\frac{1}{2\tau }\left(\mathbf{u}-{\mathbf{u}}^{{\ast}}\right).\end{equation} \tag{ B.5 }$$

After expressing equation (B.2) in terms of the actual velocity,

$$\begin{equation}{f}_{i}\left(t+{\Delta}t,\mathbf{r}+{\mathbf{e}}_{i}{\Delta}t\right)-{f}_{i}\left(t,\mathbf{r}\right)=-\frac{1}{\tau }\left[{f}_{i}-{\xi }_{i}\left(\rho ,\mathbf{v}\right)\right]+{\tilde {F}}_{i},\end{equation} \tag{ B.6 }$$

where

$$\begin{equation}{\tilde {F}}_{i}={F}_{i}+\frac{1}{\tau }\left[{\xi }_{i}\left(\rho ,\mathbf{u}\right)-{\xi }_{i}\left(\rho ,\mathbf{v}\right)\right]=\frac{1}{\tau }\left[{\xi }_{i}\left(\rho ,{\mathbf{u}}^{{\ast}}\right)-{\xi }_{i}\left(\rho ,\mathbf{v}\right)\right],\end{equation} \tag{ B.7 }$$

we can now use the results of [93] to determine our macroscopic equations. They are the continuity equation

$$\begin{equation}{\partial }_{t}\rho +\nabla \cdot \left(\rho \mathbf{v}\right)=0,\end{equation} \tag{ B.8 }$$

and the momentum equation

$$\begin{equation}{\partial }_{t}\left(\rho \mathbf{v}\right)+\nabla \cdot \left(\rho \mathbf{v}\mathbf{v}\right)=-\nabla P+\nabla \cdot \tilde {{\Pi}}+\mathbf{F}+\nabla \cdot {\Theta},\end{equation} \tag{ B.9 }$$

where

$$\begin{equation}\tilde {{\Pi}}=\rho \nu \left(\nabla \mathbf{v}+{\left(\nabla \mathbf{v}\right)}^{\mathrm{T}}\right),\end{equation} \tag{ B.10 }$$

and

$$\begin{equation}{\Theta}=\left(\tau -\frac{1}{2}\right){\Delta}t\left(\mathbf{v}\mathbf{F}+\mathbf{F}\mathbf{v}\right)-\tau \sum\limits _{i}{\mathbf{e}}_{i}{\mathbf{e}}_{i}{\tilde {F}}_{i}.\end{equation} \tag{ B.11 }$$

Several O(u³) terms have been ignored, consistent with a small Mach number assumption. The pressure $P={c}_{\mathrm{s}}^{2}\rho$ and the viscosity

$$\begin{equation}\nu =\frac{2\tau -1}{2}{c}_{\mathrm{s}}^{2}{\Delta}t,\end{equation} \tag{ B.12 }$$

both feature the effective speed of sound c_s which depends on the choice of lattice. In the case of a D2Q7 lattice, c_s = c/2. From equation (B.7),

$$\begin{equation}\tau \sum\limits _{i}{\mathbf{e}}_{i}{\mathbf{e}}_{i}{\tilde {F}}_{i}=\rho {\mathbf{u}}^{{\ast}}{\mathbf{u}}^{{\ast}}-\rho \mathbf{v}\mathbf{v}.\end{equation} \tag{ B.13 }$$

Rearranging equation (B.5), we can express both u and u* in terms of v by

$$\begin{equation}\mathbf{u}=\frac{1}{2\tau -1}\left(2\tau \mathbf{v}-{\mathbf{u}}^{{\ast}}\right).\end{equation} \tag{ B.14 }$$

Thus

$$\begin{equation}\mathbf{F}=\frac{\rho }{\tau {\Delta}t}\left({\mathbf{u}}^{{\ast}}-\mathbf{u}\right)=\frac{\rho }{{\Delta}t}\frac{2}{2\tau -1}\left({\mathbf{u}}^{{\ast}}-\mathbf{v}\right).\end{equation} \tag{ B.15 }$$

Substituting all of these into Θ, and factorizing gives

$$\begin{equation}{\Theta}=-\rho \left({\mathbf{u}}^{{\ast}}-\mathbf{v}\right)\left({\mathbf{u}}^{{\ast}}-\mathbf{v}\right).\end{equation} \tag{ B.16 }$$

Clearly ∇ ⋅ Θ has parallels with the term ∇ ⋅ (ρ vv). However, for our simulation to behave smoothly, we assume that (u* − v) ∼ O(ɛ) where ɛ is small, thus we can safely ignore ∇ ⋅ Θ, as it is of order O(ɛ²). Therefore the resulting momentum equation is

$$\begin{equation}\rho \left({\partial }_{t}\mathbf{v}+\mathbf{v}\cdot \nabla \mathbf{v}\right)=-\nabla P+\nabla \cdot \rho \nu \left(\nabla \mathbf{v}+{\left(\nabla \mathbf{v}\right)}^{\mathrm{T}}\right)+\mathbf{F}.\end{equation} \tag{ B.17 }$$

This is valid for both the MIPS model and the CIL model; the differences stem from the choice of u* in F. Both are versions of the Toner–Tu equations [3–5, 94].

B.1. MIPS model

For our MIPS model we use u* as given in equation (5), thus

$$\begin{equation}{\mathbf{F}}_{\text{MIPS}}=\frac{2\rho }{{\Delta}t\left(2\tau -1\right)}\left(C\left(1-\frac{\rho }{{\rho }^{{\ast}}}\right)\nabla \rho -\mathbf{v}\right).\end{equation} \tag{ B.18 }$$

Since $\nabla P={c}_{\mathrm{s}}^{2}\nabla \rho$, the first term on the rhs can be viewed as an effective pressure in which the coefficient varies with density. Crucially, for ρ < ρ* this coefficient is positive, hence the effective pressure acts in the opposite manner from the conventional pressure, driving material towards areas of high density instead of away from them.

The second term can be expressed as αρ v with

$$\begin{equation}\alpha =-\frac{2}{2\tau -1}{< }0.\end{equation} \tag{ B.19 }$$

As α is negative, this term acts like friction, reducing the velocity in the system. As such v decays to a value determined by a balance between the effective pressure, the conventional pressure and the stochastic noise. This balance is also seen in equation (11), thus our method is consistent with other models of MIPS.

B.2. CIL model

For our CIL model we use u* as given in equation (16). Note that

$$\begin{equation}{\mathbf{u}}^{{\ast}}={U}_{0}\frac{\mathbf{u}}{\vert \mathbf{u}\vert }={U}_{0}\frac{\mathbf{v}}{\vert \mathbf{v}\vert }.\end{equation} \tag{ B.20 }$$

Thus

$$\begin{equation}{\mathbf{F}}_{\text{CIL}}=\frac{2\rho }{{\Delta}t\left(2\tau -1\right)}\left({U}_{0}\frac{\mathbf{v}}{\vert \mathbf{v}\vert }-\mathbf{v}\right).\end{equation} \tag{ B.21 }$$

F_CIL corresponds to a linearised version of $\rho \left[\alpha \left(\rho \right)-\beta \vert \mathbf{u}{\vert }^{2}\right]\mathbf{u}$ in the Toner–Tu equations, where the linearisation has occurred about the stationary points $\vert \mathbf{u}\vert =\sqrt{\alpha /\beta }$. Matching the terms gives

$$\begin{equation}\alpha =\frac{1}{{\Delta}t\left(2\tau -1\right)},\hspace{25.0pt}\beta =\frac{\alpha }{{U}_{0}^{2}}.\end{equation} \tag{ B.22 }$$

These are not well defined when U₀ = 0. In this case

$$\begin{equation}\frac{{\mathbf{F}}_{\text{CIL}}}{\rho }=-\frac{2}{{\Delta}t\left(2\tau -1\right)}\mathbf{v}.\end{equation} \tag{ B.23 }$$

From this we can deduce nothing about β, but conclude that

$$\begin{equation}\alpha =-\frac{2}{{\Delta}t\left(2\tau -1\right)},\end{equation} \tag{ B.24 }$$

as in the MIPS model. These values of α motivate the idea that αρ u is a driving force in the system; α > 0 gives a non-zero flow U₀, whereas if α < 0 flow will decay to u = 0.

C.1. MIPS model

We are primarily concerned with the effects of varying the mean density and noise strength; we held the remaining parameters constant with the values in table C1. Note that our values are non-dimensional. To obtain the still frames in figures 2(a)–(d), we used a lattice of 180 × 60 nodes, and the parameters:

	$\bar{\rho }$	σ
(a)	0.10	0.031
(b)	0.10	0.014
(c)	0.18	0.014
(d)	0.10	0.014

Table C1. Parameters used to obtain the results depicted in figure 2. Values are non-dimensional.

Parameter	Value
τ	2
C	50
ρ*	0.2
Δx	1
Umax	0.3
c	1

To generate a phase diagram we reduced the system size to 90 × 30 nodes, due to the long runtimes required to achieve the state in figure 2(d). Doing so had no obvious effect on the qualitative behaviour of the system, hence we concluded that this system size is large enough to account for any finite-size effects. We fixed $\bar{\rho }=0.12$, and slowly varied σ. The simulation was initialized with σ = 1/70 ≈ 0.014, in a state close to the theoretical equilibrium. It was allowed to evolve for a further 10⁶ time steps before measurements were taken. A sample was then taken every 10⁴ time steps, for 100 samples. After taking the final sample, the noise strength was increased by δσ = 1/700. The system was allowed to relax for another 10⁶ time steps before the measurements resumed at the same rate as before. After a further 100 samples the noise strength was increased again by δσ and this process was repeated until σ = 3/70.

To systematically determine whether or not phase separation had occurred, we quantified the large scale density variation by a parameter $\tilde {{\Delta}\rho }$ (see appendix D). Figure 2(e) shows the mean values of $\tilde {{\Delta}\rho }$ for the 100 samples for each value of σ. The error bars are the estimated standard deviation of the samples (not the sample mean), indicating the spread of the sample distribution. Large values of $\tilde {{\Delta}\rho }$ indicate phase separation, but given that the system is noisy, $\tilde {{\Delta}\rho }$ is never exactly zero, hence we must impose a cutoff. The value of $\tilde {{\Delta}\rho }$ falls with increasing noise. Furthermore it falls faster and faster, until it quickly levels out around σ = 0.0275, at which point the standard deviation of the samples drops as well. From inspection of the samples, this is the point at which the phase separation ceases to be obvious. Hence we use ${\tilde {{\Delta}\rho }}_{\text{threshold}}=0.0430$, with $\tilde {{\Delta}\rho }{ >}{\tilde {{\Delta}\rho }}_{\text{threshold}}$ taken to be phase separated.

Using this condition to identify cases in which phase separation had occurred, we measured the density in the liquid and gas phases. We made the assumption that the system had reached the theoretical steady state, so that the dense liquid phase is entirely contained in a single vertical strip. To obtain the typical density value, both within and outside of the strip, we averaged the density profile over a rectangle spanning the entire y-axis, but only 20Δx wide. This rectangle was shifted along the entire domain, with the highest and lowest values obtained taken as the density of the liquid and gas phases respectively.

C.2. CIL model: general notes

C.3. Simulation procedure in figure 3

C.4. Simulation procedure in figure 4

To systematically categorise the system's behaviour at each of our data points, we evaluate several parameters.

The order parameter φ measures the global velocity alignment. Since we have a finite system, even in the most disordered phase φ is never exactly zero, and as it is strictly positive it cannot average to zero either. Therefore we had to impose a cutoff to distinguish between a completely disordered state and one with an ordered phase (including the B and RB phases). For our CIL model we deduced the cutoff to be φ_threshold = 0.024, with φ < φ_threshold marked as HD. The MIPS model has no mechanism for collective motion, thus φ ≈ 0 for all parameter values.

To determine whether or not phase separation had occurred, we quantified the large scale density variation by a parameter $\tilde {{\Delta}\rho }$. In both the MIPS and CIL models, due to the aspect ratio of the geometry and the periodic boundary conditions, the steady state of the phase-separated regimes manifested as strips of high density, parallel to the y-axis, as demonstrated in figures 2(d) and 3(d). In these cases, the system is largely invariant in the y-direction. Averaging over the density in that direction,

$$\begin{equation}{\rho }_{y}\left(x\right)=\frac{1}{{L}_{y}}{\int }_{0}^{{L}_{y}}\mathrm{d}y\enspace \rho \left(x,y\right),\end{equation} \tag{ D.1 }$$

thereby reduces the effects of the small scale stochastic fluctuations. We then took $\tilde {{\Delta}\rho }$ to be the standard deviation of ρ_y .

As with φ we had to impose a cutoff to distinguish between phase separated regions and local noise. For our CIL model, we chose this to be ${\tilde {{\Delta}\rho }}_{\text{threshold}}=0.012$, with $\tilde {{\Delta}\rho }{ >}{\tilde {{\Delta}\rho }}_{\text{threshold}}$ considered to be either B or RB.

The background colour in the phase diagram in figure 4(b) corresponds to the value of $\tilde {{\Delta}\rho }$, and demonstrates that our chosen cutoff value is appropriate. The two strips of yellow (high $\tilde {{\Delta}\rho }$) clearly correspond to the B and RB regimes. These strips narrow as they approach each other, and vanish entirely at the critical point (star).

The two aforementioned quantities, φ and $\tilde {{\Delta}\rho }$, do not enable us to distinguish between the B and RB phases. Therefore we also measured the following parameter:

$$\begin{equation}\zeta =\frac{{\sum }_{\mathbf{r}\in \mathrm{l}\mathrm{a}\mathrm{t}\mathrm{t}\mathrm{i}\mathrm{c}\mathrm{e}}\rho \left(\mathbf{r}\right)\left[B-\vert \mathbf{u}\left(\mathbf{r}\right)\vert \right]}{{\sum }_{\mathbf{r}\in \mathrm{l}\mathrm{a}\mathrm{t}\mathrm{t}\mathrm{i}\mathrm{c}\mathrm{e}}\rho \left(\mathbf{r}\right)},\end{equation} \tag{ D.2 }$$

where ζ is large if high density regions are stationary, as in RB. Hence we identified RB values with ζ > ζ_threshold = 0.1. In summary, for the CIL model, data points were classified using the following rules

Regime	φ	$\tilde {{\Delta}\rho }$	ζ
HD	<0.024	—	—
HO	⩾0.024	⩽0.012	—
B	⩾0.024	>0.012	⩽0.1
RB	⩾0.024	>0.012	>0.1

Visual inspection of the samples at several distinct points in the phase diagram shows the categorization to be largely accurate.

For the phase diagrams in figure 4, anything above the red line $B-A{\bar{\rho }}^{2}=0$ was marked as motionless.

Each of the plots in figures 3(d) and 2(a)–(d) are scatter plots of 180 × 60 circles, arranged in the triangular lattice configuration, as is figure 4(c) but for a lattice of 720 × 240 nodes. They represent the density and velocity fields throughout our lattice for a single time step. The colour of each dot corresponds to the direction of the macroscopic velocity u measured at that site, according to the inset colour wheel. As our system is two dimensional a single angle (colour) is sufficient to determine the direction. The area of each dot is assigned by the following formulae

$$\begin{equation}a={a}_{0}{\left(\frac{\rho -{\rho }_{\mathrm{min}}}{{\rho }_{\mathrm{max}}-{\rho }_{\mathrm{min}}}\right)}^{3/2}+{\epsilon},\end{equation} \tag{ E.1 }$$

where ρ_max, (ρ_min) is the maximum (minimum) value of ρ over all of the lattice sites in that frame, a₀ is an appropriately chosen value for the size of the biggest circle, which depends on the spacing of the particular plot. is a small value to ensure no values are exactly zero.

In figure 3(d) the samples are obtained from one of the increasing runs used to create figure 3(a), specifically sampled with the mean densities $\bar{\rho }$:

Figure 4(c) is created with the values A = 2.7, ρ = 0.175 and a lattice size of 720 × 240, with all other parameters the same as before. We have chosen a higher resolution to better illustrate the critical behaviour. This particular frame is 100 000 steps after initialization.

Five movies are provided to demonstrate the behaviour of the system as shown in figures 3(d) and 4(c), using the same visual representation method. The movies are available as supplementary material.

For each of the movies that correspond to figure 3(d), the system was initialized with A = 2 and $\bar{\rho }$ again given by

respectively. A frame was captured every five time steps. Each movie is constructed of 200 frames displayed at 12 frames per second. The movie (M5) that corresponds to figure 4(c) is created with a lattice size of 180x60 around the critical point shown in figure 4(b) (pink star).