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

TL;DR

This paper introduces a modified lattice Boltzmann method to simulate dense dry active fluids, revealing new phase transitions and providing analytical insights into their hydrodynamics.

Contribution

It presents a novel simulation approach for dense dry active matter and uncovers previously unknown phase transitions in these systems.

Findings

Discovered multiple new phase transitions, including two first-order and one potentially critical.

Validated simulation results with analytical hydrodynamic equations.

Applied method to models relevant to biological collective behaviors.

Abstract

The dynamics of dry active matter have implications for a diverse collection of biological phenomena spanning a range of length and time scales, such as animal flocking, cell tissue dynamics, and swarming of inserts and bacteria. Uniting these systems are a common set of symmetries and conservation laws, defining dry active fluids as a class of physical system. Many interesting behaviours have been observed at high densities, which remain difficult to simulate due to the computational demand. Here, we show how two-dimensional dry active fluids in a dense regime can be studied using a simple modification of the lattice Boltzmann method. We apply our method on a model that exhibits motility-induced phase separation, and an active model with contact inhibition of locomotion, which has relevance to collective cell migration. For the latter, we uncover 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.

Tables

Table 1. Table 1: Parameters used to obtain the results depicted in Fig.  2 . Values are non-dimensional.

Parameter	Value
$\tau$	2
$C$	50
${\rho }_{\ast }$	0.2
$\Delta x$	1
$U_{\max }$	0.3
$c$	1

Equations

$$\rho (t,\mathbf{r})=\sum _{i=0}^6{f}_i(t,\mathbf{r}),\mathbf{u}(t,\mathbf{r})=\frac{{\sum }_{i=0}^6{f}_i(t,\mathbf{r})c{\mathbf{e}}_i}{\rho (t,\mathbf{r})},$$

$$f_i(t+\Delta t,\mathbf{r}+c{\mathbf{e}}_i\Delta t)-f_i(t,\mathbf{r})=-\frac{1}{\tau }\left[f_i(t,\mathbf{r})-f_i^{\mathrm{SS}}(t,\mathbf{r})\right]+{\eta }_i(t,\mathbf{r}),$$

$${\eta }_i(t,\mathbf{r})={\tilde{\eta }}_i(t,\mathbf{r})-\frac{1}{7}\sum _{i=0}^6{\eta }_i(t,\mathbf{r}),$$

$$f_i^{\mathrm{SS}}=w_i\rho \left[1+4\frac{{\mathbf{e}}_i\cdot {\mathbf{u}}^{\ast }}{c}+8\frac{({\mathbf{e}}_i\cdot {\mathbf{u}}^{\ast }{)}^2}{c^2}-2\frac{|{\mathbf{u}}^{\ast }{|}^2}{c^2}\right],$$

$$\mathbf{u}^{*}=\left\{\begin{array}[]{cc}C\left(1-\frac{\rho}{\rho^{*}}\right)\nabla\rho\ ,&{\rm if}\rho<\rho^{*},\\ 0\ ,&{\rm otherwise},\end{array}\right.$$

$$|{\mathbf{u}}^{\ast }|\mapsto \min \left(|{\mathbf{u}}^{\ast }|,U_{\max }\right).$$

$$\nabla \rho =a\left[({\rho }_1-{\rho }_4){\mathbf{e}}_1+({\rho }_2-{\rho }_5){\mathbf{e}}_2+({\rho }_3-{\rho }_6){\mathbf{e}}_3\right],$$

$$\frac{\partial {\rho }_0}{\partial x}$$

$$\frac{\partial {\rho }_0}{\partial y}$$

$$\frac{\partial \rho }{\partial t}$$

$$\mathbf{J}$$

$$D$$

$$\frac{\mathbf{V}}{D}$$

$$v(\rho )=v_0\left(1-\frac{\rho }{{\rho }^{\ast }}\right).$$

$$\mathbf{V}=\frac{\tau v_0^2}{{\rho }^{\ast }d}\left(1-\frac{\rho }{{\rho }^{\ast }}\right)\nabla \rho .$$

$${\mathbf{u}}^{\ast }=\frac{U_0}{|\mathbf{u}|}\mathbf{u},$$

$$U_{0}(\rho)=\left\{\begin{array}[]{cc}\hfill-A\rho^{2}+B&{\rm if}B\geq A\rho^{2}\ ,\\ \hfill 0&{\rm otherwise}.\end{array}\right.$$

$$\phi =\left|\frac{{\sum }_{\mathbf{r}\in \mathrm{lattice}}\rho (\mathbf{r})\mathbf{u}(\mathbf{r})}{{\sum }_{\mathbf{r}\in \mathrm{lattice}}\rho (\mathbf{r})}\right|.$$

$${\partial }_t\rho +\nabla \cdot \left(\rho \mathbf{v}\right)=0,$$

$$\rho \left({\partial }_t\mathbf{v}+\mathbf{v}\cdot \nabla \mathbf{v}\right)=-\nabla P+\nu \nabla \cdot \rho \left(\nabla \mathbf{v}+(\nabla \mathbf{v}{)}^T\right)+\rho \left[\alpha (\rho )-\beta |\mathbf{v}{|}^2\right]\mathbf{v},$$

$$\left(s-\mathrm{i}k\sqrt{\frac{{\alpha }_0}{\beta }}\right)\delta \rho =\mathrm{i}k{\rho }_0\delta v,$$

$${\rho }_0\left(s-\mathrm{i}k\sqrt{\frac{{\alpha }_0}{\beta }}+2{\alpha }_0+2\nu k^2\right)\delta v=\left(\mathrm{i}k{c}_s^2+{\rho }_0{\alpha }_1\sqrt{\frac{{\alpha }_0}{\beta }}\right)\delta \rho .$$

$${\rm Re}[s]=\left\{\begin{array}[]{l}-2\alpha_{0}+{\cal O}(k^{2})\ ,\\ \frac{k^{2}}{8\alpha_{0}^{2}\beta}\left[\alpha_{1}^{2}\rho_{0}^{2}-4\alpha_{0}c_{s}^{2}\beta\right]+{\cal O}(k^{3})\ ,\end{array}\right.$$

$$\sum _{i=0}^7{\eta }_i=0,$$

$$\mathbf{h}=\sum _{i=0}^7{\eta }_{i}c{\mathbf{e}}_i.$$

$$⟨h_x⟩=⟨h_y⟩=⟨h_x{h}_y⟩=0,$$

$$⟨h_x(\mathbf{x},t)h_x({\mathbf{x}}^{\prime },t^{\prime })⟩=⟨h_y(\mathbf{x},t)h_y({\mathbf{x}}^{\prime },t^{\prime })⟩=c^2{\sigma }^2\delta (\mathbf{x}-{\mathbf{x}}^{\prime })\delta (t-t^{\prime }).$$

$$\rho \left({\partial }_t\mathbf{v}+\mathbf{v}\cdot \nabla \mathbf{v}\right)=-\nabla P+\nabla \cdot \rho \nu \left(\nabla \mathbf{v}+(\nabla \mathbf{v}{)}^T\right)+{\mathbf{F}}_{\mathrm{CIL}}+\mathbf{h}.$$

$$\delta \alpha$$

$$=-\beta \frac{c^2{\sigma }^2}{{\rho }^2}\int \frac{{\mathrm{d}}^{\mathrm{D}}\mathbf{q}}{(2\pi {)}^{\mathrm{D}}}\frac{1}{\mu q^2+\alpha }$$

$$=-\frac{c^2{\sigma }^2\beta }{{\rho }^2\left(\mu {\Lambda }^2+\alpha \right)}\frac{S_{\mathrm{D}}\mathrm{d}\ell }{(2\pi {)}^{\mathrm{D}}},$$

$${\alpha }_R=\alpha \left(1-\frac{\sigma }{{\rho }^2{U}_0^2\mu {\Lambda }^2}\frac{S_{\mathrm{D}}\mathrm{d}\ell }{(2\pi {)}^{\mathrm{D}}}\right)=\alpha \left(1-C_1\frac{1}{{\rho }^2{U}_0^2}\right),$$

$$\bar{\rho }=\frac{3\sqrt{C_1}}{2B},A=\frac{4B^3}{27C_1}.$$

$$\rho (t,\mathbf{r})$$

Full text

††: New J. Phys.

1 Introduction

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, 4, 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 (DRG) 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 [13]. Superseding the lattice gas cellular automata are the celebrated lattice Boltzmann methods (LBM) [14, 15, 16, 17, 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, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 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. 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 motility-induced phase separation (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 the Chapman-Enskog coarse-graining procedure.

2 Dry polar active fluids

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

[TABLE]

where $c$ is the ratio of the lattice spacing $\Delta x$ to the time increment $\Delta t$. Typically, scales are chosen so that $c=1$, and $|\mathbf{u}|\ll c$ is required for stable solutions (see A.1).

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

[TABLE]

where $\tau$ is a relaxation parameter, $f_{i}^{\rm SS}$ is the steady-state distribution, and

[TABLE]

where $\tilde{\eta}_{i}$ are temporally and spatially uncorrelated random variables uniformly distributed in the interval $[-\sigma,\sigma]$. This particular choice of noise alters the local velocity, but Eq. 3 ensures the density is conserved. Note that the fluctuations do not conserve momentum as we are considering dry active fluids.

In standard LBM, $f_{i}^{\rm SS}$ is constructed so as to preserve the moments $\rho$, $\mathbf{u}$ at each site, accomplished in a D2Q7 model by

[TABLE]

with $\mathbf{u}^{*}=\mathbf{u}$, where $w_{i}$ are lattice specific weights: $w_{0}=1/2$ and $w_{i\neq 0}=1/12$. However, in breaking the conservation of local momentum, the restriction on the form of the steady-state distribution $f_{i}^{\rm SS}$ is drastically relaxed. We make a simple modification: maintaining the form of $f_{i}^{\rm SS}$ in Eq. 4, but taking $\mathbf{u}^{*}$ to be a function of $\rho$ and $\mathbf{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 B). This method can produce a variety of behaviours associated with dry active matter, depending on the choice of $\mathbf{u}^{*}$.

3 Motility-Induced Phase Separation

To demonstrate the versatility of our method, we first employ it to produce motility-induced phase separation (MIPS). MIPS can occur when self-driven particles that maintain their direction of propulsion with some persistency, interact through short-range repulsive forces [40, 41, 42, 43, 44, 45, 46, 13]. 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

[TABLE]

where $C$ is a positive constant, dictating the sensitivity of the system to density gradients. The parameter $\rho^{*}$ 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 $\nabla\rho$ could in principle become very large and any physical system will not have a divergent speed, we impose further an upper bound on $|\mathbf{u}^{*}|$. Explicitly, we have

[TABLE]

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

[TABLE]

for $a=1/(3\Delta x)$, where in this notation $\rho_{i}=\rho(t,\mathbf{r}+\Delta x\mathbf{e}_{i})$ [51]. The resulting stencils are

[TABLE]

We motivate the expression in Eq. 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 $\mathbf{u}^{*}$ is a macroscopic variable. However, a consequence of all particles’ speeds decreasing is that the net speed $|\mathbf{u}^{*}|$ will also decrease. This is clearly achieved by the term $\left(1-{\rho}/{\rho^{*}}\right)$, with the net speed hitting zero altogether for $\rho\geq\rho^{*}$. 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 $\rho^{*}$. 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, whilst 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 Eq. 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(\mathbf{r})$, and orientational relaxation time $\tau(\mathbf{r})$, where $v$ and $\tau$ can vary in space, it has been shown that the macroscopic density can be approximated by the following

[TABLE]

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

[TABLE]

where $d$ is the dimension. Generalizing this system to include interaction effects, allows $v$ and $\tau$ to now depend on the local density $\rho$, as well as $\mathbf{r}$. By extension so do $D$ and $\bf V$. For simplicity we assume $\tau$ to be constant and $v$ to depend only on $\rho$. 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

[TABLE]

The drift velocity $\bf V$ becomes

[TABLE]

Finally we observe that Eq. 10 is a continuity equation and equate $\mathbf{J}\equiv\rho\mathbf{u}$ in our system. We take the drift velocity to be the steady state that $\mathbf{u}$ tends towards, that is $\bf V\equiv\mathbf{u}^{*}$. The term $D\nabla\rho$ features in our model as the “pressure” term in our hydrodynamic equations (see B), and we further incorporate fluctuations in the system. Our system is very similar to that presented in Eq. 10, 11, albeit our velocity $\mathbf{u}$ is dynamic, whilst $\mathbf{J}$ has no time dependence. In that sense, Eq. 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 C. Fig. 2 shows the results from our MIPS simulations. Fig. 2(a)-(d) are still frames from our simulations displayed using the visualization scheme described in E. In panels (b)-(d) MIPS has occurred. The multi-coloured nature of the the dense region indicates there is no velocity alignment between neighbouring lattice sites. Therefore, collectively, the dense regions have effectively zero velocity.

Fig. 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, $(\bar{\rho},\sigma)$, 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 $\sigma$, respectively.

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

4 Contact Inhibition of Locomotion

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

[TABLE]

where $U_{0}$ is of the form

[TABLE]

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}^{\rm 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 contact inhibition of locomotion (CIL), as observed in cell monolayers [52, 53, 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_{0}$ a much more complicated function of $\rho$.

As for our MIPS model, we simulate our active fluid model in a rectangular geometry with periodic boundary conditions (Fig. 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 $\varphi$, given by

[TABLE]

Necessarily $0\leq\varphi\leq B$. Fig. 3(a) shows how $\varphi$ 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, 58, 59, 60, 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 (Fig. 3(d)), but with the properties inverted—the high density regions are disordered, whilst 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, whilst 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].

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.

5 Phase behaviour

To further elucidate the phase behaviour of the system, we now vary the contact inhibition parameter $A$ as well as $\bar{\rho}$ (Fig.  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 Fig.  4(a), which is a putative critical transition. A snapshot of this critical behaviour is shown in Fig. 4(b). We will discuss this critical point further in our hydrodynamic analysis. Finally, we note that in addition to the four regimes already discussed we recognise a fifth somewhat trivial phase of no motion (indicated by red asterisks), where $|\mathbf{u}|=0$ since $-A\bar{\rho}^{2}+B\leq 0$.

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, 67, 68, 69, 70, 71, 72, 73], which are complimentary to this study.

6 Hydrodynamic description of the transitions

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 B. We now describe our simulation findings analytically by performing a linear stability analysis on our derived version of the Toner-Tu equations:

[TABLE]

where $P=c_{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 $\rho=\rho_{0}+\delta\rho\exp(st-{\rm i}kx)$, $v=\sqrt{\alpha_{0}/\beta}+\delta v\exp(st-{\rm i}kx)$ and $\alpha=\alpha_{0}+\alpha_{1}\delta\rho$, we have to linear order:

[TABLE]

Note that for a stable system $\beta>0$ necessarily, thus $\alpha_{0}>0$ implies collective motion, with $\alpha_{0}\leq 0$ giving a stationary state.

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

[TABLE]

where ${\rm Re}[s]>0$ implies instability. The first eigenvalue $-2\alpha_{0}$ 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 $\alpha_{0}\to 0$ and $\alpha_{1}\neq 0$, the system is always unstable. In the previously mentioned studies, [11, 65, 58, 59, 60, 61], $\alpha_{1}$ is typically assumed to be positive, so that alignment interactions increase with density, and as $\alpha$ 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 $\alpha_{1}$. Hence, if $\alpha$ 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 $\alpha_{0}$ and $\alpha_{1}$ 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 $\alpha(\rho)$ in our system. The blue dashed curve in Fig. 4(a) was obtained by fitting the curve $\alpha_{R}(\bar{\rho})=0$ to pass through the critical point. The resulting curve bisects both the B and RB regions, which is in perfect agreement with our stability analysis; close to this line $\alpha_{0}\approx 0$ and $\alpha_{1}\neq 0$, hence instability, with the exception of close to the critical point where $\alpha_{1}\approx 0$.

7 Fluctuation-induced renormalization of $\alpha$

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 $\alpha$ will naturally depend on the density $\rho$ as well. Here, we aim mainly to connect our hydrodynamic EOM derived from our Chapman-Enskog expansion (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 $\beta$ term, together with the noise term with strength $\sigma$, renormalise $\alpha$ around the critical region when $\alpha\approx 0$.

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

[TABLE]

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

[TABLE]

From the definition of $\eta_{i}$ it follows that

[TABLE]

Our derived Eq. 57 now becomes

[TABLE]

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 $\alpha$ is given by

[TABLE]

where $\bf f=\frac{{\bf h}}{\rho}$, $S_{\rm D}\equiv 2\pi^{{\rm D}/2}/\Gamma({\rm D}/2)$ is the surface area of a unit sphere in D dimensions, $\Lambda^{-1}$ corresponds to the small length scale cutoff, and $\Lambda{\rm d}\ell$ corresponds to the width of the momentum shell being integrated over in the above coarse-graining process.

Adding this shift to $\alpha$ modifies it to, to ${\cal O}(\alpha)$:

[TABLE]

where we have used the expression for $\beta$ in 62.

The above expression indicates that the critical transition occurs when $U_{0}$ and $\rho$ are fine tuned such that $\alpha_{R}=0$ and $\partial\alpha_{R}/\partial\rho=0$, which can be achieved by fine tuning $\bar{\rho}$ and $A$ in the expression of $U_{0}$. These are given by

[TABLE]

In our simulations we found $\bar{\rho}_{c}\approx 0.18$, $A_{c}\approx 2.75$. These values give $B\approx 0.27,C_{1}\approx 0.001$, which is fairly consistent with our actual value of $B=0.3$. In fact we can go further. In Fig. 3(a) (and Fig. S3) the blue dashed line depicts the full contour for $\alpha_{R}=0$ using these values. For parameter values below this curve $\alpha>0$ and we have collective motion, with $\alpha<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 perfect agreement with our stability analysis; in the vicinity of this curve $\alpha_{0}\approx 0$ and $\alpha_{1}\neq 0$ except at the critical point.

8 Critical phenomena in active systems

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 motility-induced phase separation [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.

9 Conclusion & Outlook

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, 86, 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, 89, 90]. Finally, we believe our method can be straightforwardly extended to higher dimensions.

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

Appendix A Simulation Algorithm

A.1 Conventional Lattice Boltzmann Methods

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,\mathbf{r},\bm{\xi})$, which corresponds to the mass density of matter moving with velocity $\bm{\xi}$ at location $\mathbf{r}$ and time $t$.

The typical variables of interest are

[TABLE]

where D is the spatial dimension, $\rho$ is the mass density field, and $\mathbf{u}$ is the velocity field.

Discretising $\bm{\xi}$ 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 DQ 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 discretised distribution function $f_{i}(t,\mathbf{r})$, where $i=0,1,...,6$. The corresponding vectors are ${\mathbf{e}}_{0}=\mathbf{0}$, and ${{\mathbf{e}}_{i}}=\cos[(i-1)\pi/6]\hat{\mathbf{x}}+\sin[(i-1)\pi/6]\hat{{\bf y}}$ for $i>0$. The local mass density and velocity are thus given by

[TABLE]

where $c$ is the ratio of the lattice spacing $\Delta x$ to the time increment $\Delta t$. For accurate solutions velocities must obey $|\mathbf{u}|\ll 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]

[TABLE]

where $f_{i}^{eq}(t,\mathbf{r})$ is a lattice dependent steady-state (equilibrium) distribution and $\tau$ is a relaxation parameter. When simulating the Navier Stokes equations, $\tau$ 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. 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}^{eq}$ is obtained by discretising the Boltzmann distribution for velocities at thermal equilibrium [92] and in a D2Q7 framework is given by

[TABLE]

where $w_{i}$ are lattice specific weights: $w_{0}=1/2$ and $w_{i\neq 0}=1/12$. This discretisation is possible through the low Mach number assumption $|\mathbf{u}|\ll c$.

A.2 Summary of modified lattice Boltzmann model algorithm

A single time step is partitioned into the following steps

Streaming: $f_{i}^{*}(t+\Delta t,\mathbf{r}+c\mathbf{e}_{i}\Delta t)=f_{i}(t,\mathbf{r})$. 2. 2.

Calculate $\rho$, $\mathbf{u}$ from $f_{i}^{*}$. 3. 3.

Calculate $\mathbf{u}^{*}$ and hence $f_{i}^{\rm SS}$. 4. 4.

Redistribution (collision): $f_{i}=f_{i}^{*}-\frac{1}{\tau}\left[f_{i}^{*}-f_{i}^{\rm SS}\right]$. 5. 5.

Add noise: $f_{i}\mapsto f_{i}+\eta_{i}$. 6. 6.

Correct any negative values of $f_{i}$ (see Section A.3).

A.3 Negative Correction

In our modified LBM, we have included stochastic noise as detailed in the main text. It is possible for one or more of the $f_{i}$ values to become negative in this process. This is not permitted, hence we correct this using the following procedure at each lattice site.

If $f_{0}<0$, simply set $f_{0}=0$. 2. 2.

If $f_{i}<0$ for $i>0$, then take $j$, such that $\mathbf{e}_{j}$ is the reverse direction of $\mathbf{e}_{i}$ ( $1\longleftrightarrow 4,2\longleftrightarrow 5,3\longleftrightarrow 6$), and add $|f_{i}|$ to $f_{j}$, then set $f_{i}=0$. 3. 3.

By now $f_{i}\geq 0$ for all $i$, however it is possible that the density has changed. This is corrected by rescaling

[TABLE]

Also note that due to both the noise and the rescaling of $\mathbf{u}$ in $f_{i}^{\rm SS}$, it is possible that

[TABLE]

when $\tau<1$, even if $f_{i}>0$ for all $i$. Therefore to ensure that no values become negative, we only use $\tau\geq 1$. Different values of $\tau$ were tested, yielding results that were qualitatively similar to those described in the letter.

Appendix B Lattice Boltzmann to Toner-Tu equations via Chapman-Enskog expansion

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’s introduce the notation

[TABLE]

Hence in the conventional LBM $f_{i}^{eq}=\xi_{i}(\rho,\mathbf{u})$ and in our active model $f_{i}^{\rm SS}=\xi_{i}(\rho,\mathbf{u}^{*})$.

The evolution equation becomes

[TABLE]

where

[TABLE]

Note that

[TABLE]

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

[TABLE]

After expressing Eq. 42 in terms of the actual velocity,

[TABLE]

where

[TABLE]

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

[TABLE]

and the momentum equation

[TABLE]

where

[TABLE]

and

[TABLE]

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

[TABLE]

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 Eq. 47,

[TABLE]

Rearranging Eq. 45, we can express both $\mathbf{u}$ and $\mathbf{u}^{*}$ in terms of $\mathbf{v}$ by

[TABLE]

Thus

[TABLE]

Substituting all of these into $\Theta$, and factorizing gives

[TABLE]

Clearly $\nabla\cdot\Theta$ has parallels with the term $\nabla\cdot(\rho\mathbf{v}\mathbf{v})$. However, for our simulation to behave smoothly, we assume that $(\mathbf{u}^{*}-\mathbf{v})\sim O(\varepsilon)$ where $\varepsilon$ is small, thus we can safely ignore $\nabla\cdot\Theta$, as it is of order $O(\varepsilon^{2})$. Therefore the resulting momentum equation is

[TABLE]

This is valid for both the MIPS model and the CIL model; the differences stem from the choice of $\mathbf{u}^{*}$ in $\mathbf{F}$. Both are versions of the Toner-Tu equations [3, 4, 94, 5].

B.1 MIPS model

For our MIPS model we use $\mathbf{u}^{*}$ as given in Eq. 5, thus

[TABLE]

Since $\nabla P=c_{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 $\rho<\rho^{*}$ 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 $\alpha\rho\mathbf{v}$ with

[TABLE]

As $\alpha$ is negative, this term acts like friction, reducing the velocity in the system. As such $\mathbf{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 Eq. 11, thus our method is consistent with other models of MIPS.

B.2 CIL model

For our CIL model we use $\mathbf{u}^{*}$ as given in Eq. 16. Note that

[TABLE]

Thus

[TABLE]

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

[TABLE]

These are not well defined when $U_{0}=0$. In this case

[TABLE]

From this we can deduce nothing about $\beta$, but conclude that

[TABLE]

as in the MIPS model. These values of $\alpha$ motivate the idea that $\alpha\rho\mathbf{u}$ is a driving force in the system; $\alpha>0$ gives a non-zero flow $U_{0}$, whereas if $\alpha<0$ flow will decay to $\mathbf{u}=\mathbf{0}$.

Appendix C Simulation procedure for results

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 1. Note that our values are non-dimensional.

To obtain the still frames in Fig. 2(a)-(d), we used a lattice of $180\times 60$ nodes, and the parameters:

[TABLE]

To generate a phase diagram we reduced the system size to $90\times 30$ nodes, due to the long runtimes required to achieve the state in Fig. 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 $\sigma$. The simulation was initialized with $\sigma=1/70\approx 0.014$, in a state close to the theoretical equilibrium. It was allowed to evolve for a further $10^{6}$ time steps before measurements were taken. A sample was then taken every $10^{4}$ time steps, for $100$ samples. After taking the final sample, the noise strength was increased by $\delta\sigma=1/700$. The system was allowed to relax for another $10^{6}$ time steps before the measurements resumed at the same rate as before. After a further $100$ samples the noise strength was increased again by $\delta\sigma$ and this process was repeated until $\sigma=3/70$.

To systematically determine whether or not phase separation had occurred, we quantified the large scale density variation by a parameter $\widetilde{\Delta\rho}$ (see D). Fig. 2(e) shows the mean values of $\widetilde{\Delta\rho}$ for the $100$ samples for each value of $\sigma$. 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 $\widetilde{\Delta\rho}$ indicate phase separation, but given that the system is noisy, $\widetilde{\Delta\rho}$ is never exactly zero, hence we must impose a cutoff. The value of $\widetilde{\Delta\rho}$ falls with increasing noise. Furthermore it falls faster and faster, until it quickly levels out around $\sigma=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 $\widetilde{\Delta\rho}_{\rm threshold}=0.0430$, with $\widetilde{\Delta\rho}>\widetilde{\Delta\rho}_{\rm 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\Delta 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

We used a $180\times 60$ node triangular lattice with periodic boundary conditions. Increasing the size to $360\times 120$ nodes produced no obvious qualitative change in behaviour, hence we deduced that $180\times 60$ was large enough to negate finite-size effects. Every simulation was initialized by choosing 2 uniformly distributed random numbers for each lattice site—the first deciding the initial direction for $\mathbf{u}$ with each direction in $[0,2\pi]$ equally likely, and the other choosing the initial density such that it falls within $\pm 10$% of $\bar{\rho}$, the prescribed average density. The initial magnitude of $\mathbf{u}$ is determined by Eq. 17, and the initial distribution given by $f_{i}^{\rm SS}$ in Eq. 4. The system was then evolved 10000 steps before a “run” would formally begin. This marks the end of the initialization process.

As detailed below, during the course of each run the mean density would be deliberately changed by a small amount $\delta\bar{\rho}$. On each occurrence, this was achieved by adding $\delta\bar{\rho}/7$ to every $f_{i}$ across the lattice—the lattice was not reinitialized. For measurements involving averaging over multiple runs the initialization process was repeated between each run.

With the exception of $A$ and $\bar{\rho}$, the parameters were held constant with $\Delta x=1$, $\Delta t=1$, $c=1$, $B=0.3$, $\tau=1$, and $\sigma=1/70$. As with the MIPS model, these parameters are all non-dimensional. In conventional LBM $\tau$ is linearly proportional to the viscosity of the system. We expect a similar relationship here, but the effect of viscosity is not the main focus of this study.

C.3 Simulation procedure in Fig. 3

For Fig. 3(a) $A=2$, and only $\bar{\rho}$ is varied. The goal was to detect any hysteresis effect when adiabatically increasing or decreasing the density. For the increasing density run, we start with $\bar{\rho}=0.09$ and increase it by $|\delta\bar{\rho}|=0.0025$ every $t_{\rm sample}$ steps. A sample was taken just before each density change. For higher densities the system would reach equilibrium faster, hence in order to be able to detect the hysteresis effect between HO-RB and RB-HD, a lower value of $t_{\rm sample}=500$ was used for $\bar{\rho}>0.2$, with $t_{\rm sample}=1500$ otherwise. This is because in a finite system, running the system long enough would eliminate all hysteresis effects.

For the decreasing density run, we start with $\bar{\rho}=0.35$ and again decrease it by $|\delta\bar{\rho}|=0.0025$ every $t_{\rm sample}$ steps, which is again 500 for $\bar{\rho}>0.2$ and 1500 otherwise. The system was reinitialized between increasing and decreasing runs. As a result, they should be treated as independent.

The whole increasing and decreasing density cycle was repeated over 50 runs. The order parameter $\varphi$ was evaluated for each sample, and the values for samples at equal densities $\bar{\rho}$ taken in the same direction of the cycle were collated and averaged. The error bars are the standard deviation of $\bar{\varphi}$, that is, the estimate of the standard deviation of the samples $\sigma_{\varphi}$ scaled by the square root of the sample size, $\sigma_{\varphi}/\sqrt{50}$.

C.4 Simulation procedure in Fig. 4

For Fig. 4 $A$ was varied as well as $\bar{\rho}$. Results were obtained by fixing $A$ then slowly varying $\bar{\rho}$, much as for Fig. 3(a) except: there were no descending runs, $\Delta\rho=0.01$, and $t_{\rm sample}=10000$, which is typically sufficient for the system to reach its steady state. Each run was repeated 10 times resulting in 10 independent samples for each point $(\bar{\rho},A)$. For each of these samples, the parameters $\varphi,\widetilde{\Delta\rho},\zeta$ (see Section D) were calculated, then averaged for each point $(\bar{\rho},A)$. The resulting values were used to deduce the behaviour of the system for those parameters.

Appendix D Regime Categorisation in the phase diagram

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

The order parameter $\varphi$ measures the global velocity alignment. Since we have a finite system, even in the most disordered phase $\varphi$ 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 $\varphi_{\rm threshold}=0.024$, with $\varphi<\varphi_{\rm threshold}$ marked as HD. The MIPS model has no mechanism for collective motion, thus $\varphi\approx 0$ for all parameter values.

To determine whether or not phase separation had occurred, we quantified the large scale density variation by a parameter $\widetilde{\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 Fig. 2(d) and Fig. 3(d). In these cases, the system is largely invariant in the $y$-direction. Averaging over the density in that direction,

[TABLE]

thereby reduces the effects of the small scale stochastic fluctuations. We then took $\widetilde{\Delta\rho}$ to be the standard deviation of $\rho_{y}$.

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

The background colour in the phase diagram in Fig. 4(b) corresponds to the value of $\widetilde{\Delta\rho}$, and demonstrates that our chosen cutoff value is appropriate. The two strips of yellow (high $\widetilde{\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, $\varphi$ and $\widetilde{\Delta\rho}$, do not enable us to distinguish between the B and RB phases. Therefore we also measured the following parameter:

[TABLE]

where $\zeta$ is large if high density regions are stationary, as in RB. Hence we identified RB values with $\zeta>\zeta_{\rm threshold}=0.1$. In summary, for the CIL model, data points were classified using the following rules

[TABLE]

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 Fig. 4, anything above the red line $B-A\bar{\rho}^{2}=0$ was marked as motionless.

Appendix E Representation of static configurations

Each of the plots in Fig. 3(d) and Fig. 2(a)-(d) are scatter plots of $180\times 60$ circles, arranged in the triangular lattice configuration, as is Fig. 4(c) but for a lattice of $720\times 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 $\mathbf{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

[TABLE]

where $\rho_{\max},(\rho_{\min})$ is the maximum (minimum) value of $\rho$ over all of the lattice sites in that frame, $a_{0}$ is an appropriately chosen value for the size of the biggest circle, which depends on the spacing of the particular plot. $\epsilon$ is a small value to ensure no values are exactly zero.

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

[TABLE]

Fig. 4(c) is created with the values $A=2.7$, $\rho=0.175$ and a lattice size of $720\times 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 $100000$ steps after initialization.

Appendix F Movies

Four movies are provided to demonstrate the behaviour of the system as shown in Fig. 3(d) and Fig. 4(c), using the same visual representation method. The movies are available at https://www.dropbox.com/sh/0q16xf2t9tfi8h1/AACC8AYrVLwhkYegVimgBdega?dl=0.

For each of the movies the system was initialized with $A=2$ and $\bar{\rho}$ again given by

[TABLE]

respectively. A frame was captured every 5 time steps. Each movie is constructed of 200 frames displayed at 12 frames per second.