Motility-Induced Clustering and Meso-Scale Turbulence in Active Polar Fluids

Meso-scale turbulence was originally observed experimentally in various suspensions of swimming bacteria, as well as in the collective motion of active colloids. The corresponding large-scale dynamical patterns were reproduced in a simple model of a polar fluid, assuming a constant density of active particles. Recent, more detailed experimental studies revealed additional interesting aspects, such as anomalous velocity statistics and clustering phenomena. Those phenomena cannot be explained by currently available models for active polar fluids. Herein, we extend the continuum model suggested by Dunkel et al. to include density variations and a feedback between the local density and self-propulsion speed of the active particles. If the velocity decreases strong enough with the density, a linear stability analysis of the resulting model shows that, in addition to the short-wavelength instability of the original model, a long-wavelength instability occurs. This is typically observed for high densities of polar active particles and is analogous to the well-known phenomenon of motility-induced phase separation (MIPS) in scalar active matter. We determine a simple phase diagram indicating the linear instabilities and perform systematic numerical simulations for the various regions in the corresponding parameter space. The interplay between the well understood short-range instability and the long-range instability leads to interesting dynamics and novel phenomena concerning nucleation and coarsening processes. Our simulation results display a rich variety of novel patterns, including phase separation into domains with dynamically changing irregularly shaped boundaries. Anomalous velocity statistics are observed in all phases where the system segregates into regions of high and low densities. This offers a simple explanation for their occurrence in recent experiments with bacterial suspensions.


Introduction
Exploring active matter has become a popular subject in contemporary physics leading to many new insights into a large variety of intriguing systems. Active systems are ubiquitous in nature, thereby drawing interest from different scientific communities (physics, chemistry, biology, material science, ecology, robotics) and offering a wealth of surprising dynamic phenomena [1,2,3,4,5,6]. Moreover, they provide many new challenges for our understanding of nonequilibrium systems. Realizations of active matter systems range from intracellular processes and bacterial suspensions [1,6,7,8], artificial Janus particles [9,10,11] to schools of fish and flocks of birds [12,13]. More recent reviews have focused on the prominent role of models with alignment interaction [14] and on anisotropic, self-propelled particles [15] as well as on the large variety of computational approaches to active matter [16] and on a roadmap outlining a multitude of promising directions for the field [17].
Among various collective states that characterize active matter, one phenomenon of general interest is meso-scale turbulence. Meso-scale turbulence is reported in various experimental studies, for example for suspensions of Bacillus subtilis [18,19], Escherichia coli [20,21] and Serratia marcescens [22]. The main feature (and difference to ordinary inertial turbulence) of meso-scale turbulence in bacterial suspensions is the appearance of a characteristic length scale [18,20,23,24,25,26,27]. A continuum model that agrees with experimental findings for wild-type Bacillus subtilis suspensions was presented in [23,28,29,30,31]. When solving this model numerically for a broad range of parameters, the observed velocity statistic is close to Gaussian [32] in agreement with experimental findings [23].
Recent work on discrete models with self-propelled rods revealed that effective polar alignment is often observed in self-propelled rods with steric repulsion [33,34]. It is important to note that the model for meso-scale turbulence analyzed in this paper is appropriate for the description of polar fluids and does not apply to active turbulence in so called "active nematics", see e. g. [35,36], wherein a different kind of turbulence without a characteristic length-scale, compare discussion in [15], is observed. Current experiments on engineering of vortex lattices in bacterial suspensions [37] indeed showed that a typical length-scale is controlling the behavior in meso-scale turbulence. As a result, the continuum model of meso-scale turbulence described in detail above allowed to reproduce the characteristics of these experiments [38].
However, in recent experiments on Bacillus subtilis suspensions anomalous velocity statistics have been observed. For example, swarms of very short or very long cells (compared to the wild-type) show anomalous velocity statistics [39]. Moreover, adding sublethal concentrations of antibiotics to wild-type swarms produces anomalous velocity statistics [40]. The deviations from normal statistics are quantified by measuring the kurtosis κ (scaled fourth moment), where κ = 3 indicates anomalous statistics. Such anomalous statistics are not reported for the theory presented in [23,29]. Hence, a theory accounting for anomalous statistics in meso-scale turbulent systems is still lacking.
We present and analyse a minimal model based on [23,29], exhibiting mesoscale turbulence and anomalous statistics. The main idea is summarized as follows: The model introduced in [23,29] assumes constant density and constant self-propulsion speed. We relax these assumptions by allowing for velocity variations mediated by density variations. This is a very natural assumption, as a dependency of the speed on the density is reported for Bacillus subtilis suspensions in several experimental studies [32,41]. More specifically, we borrow ideas from motility-induced phase separation (MIPS) [42,43,44] to model density variations.
Combining ideas from meso-scale turbulence and MIPS is intriguing from a general point of view as well. In bacterial suspensions steric interactions (through volume exclusion), alignment (through elongated shapes) and hydrodynamic interactions (through self-propulsion in the surrounding medium) are assumed to be present simultaneously. Considering steric interactions and alignment individually gives rise to MIPS [42,43,44] and global order in Vicsek-like models [45,46,47,48] respectively, while meso-scale turbulence results from the combination of alignment and hydrodynamics [23,29,30,31]. However, a theory aiming to comprehensively describe the dynamics of bacterial suspensions needs to incorporate all three of these effects and the interplay between them. Recently, several studies focusing on the interplay of steric interactions and alignment in active matter [49,33,34,50,51,52,53,54,55] report interesting and sometimes contradicting results. However, hydrodynamic interactions are commonly neglected. Our model features elements from all three of these prominent theories of active matter (MIPS, global order, meso-scale turbulence) and connects them in a minimalist fashion. Hence, our work contributes to the current discussion on how to connect different branches of active matter and provides an insight into the expected dynamics. We remark that in this study we mainly focus on the interplay between MIPS and the finite-wavelength instability arising due to hydrodynamic interactions, while global order will be of less relevance.
The structure of the paper is as follows: In the second section we propose a phenomenological model based on continuum models well established in the literature, which combines their central features. In the third section we present a linear stability analysis, hinting at the expected dynamics. In the fourth section we present numerical solutions of our model, sketch a phase portrait and discuss the observed anomalous velocity statistics.

Modeling Approach
We present our phenomenological model in the following steps: First, we revisit two continuum models. The first model describes MIPS, while the second model describes meso-scale turbulence. Based on these models, we propose a minimal phenomenological model that combines their main features.
When discussing continuum models of active matter, it is helpful to keep the microscopic picture in mind. We consider active particles that self-propel with some speed along their individual axes. Coarse-graining gives an averaged polarization, which, multiplied with the speed, coincides with the macroscopic velocity. Hence, the polarization plays a dual role of both order parameter and velocity field [1,4,56].

Revisiting established models
A minimal hydrodynamic model to describe MIPS was presented in [43,44]. It consists of coupled equations for the particle density ρ and polarization density field p where D and D r are diffusion coefficients. The polarization density p is given by the averaged orientation of the self-propelled particles. The coupling between the density ρ and the polarization density p is achieved through a density- modelling a decrease of the self-propulsion speed at very high density (jamming) [32,57]. A transition from a homogeneous density profile to phase separation is encountered for sufficiently high densities (or alternatively a strong enough damping constant ζ). The phase separation can be understood by assuming a slow variation of p in time and space. Setting the corresponding derivatives in Eq. (1b) to zero and substituting the resulting expression for p into Eq. (1a) gives Above a critical density, the sign of the effective diffusion coefficient D(ρ) changes due to v (ρ) < 0, triggering phase separation. For more details refer to [43,44,56]. Taking a different approach, in [23,29] the authors present a model which reproduces the statistical features of meso-scale turbulence as observed in dense suspensions of Bacillus subtilis. This model was proposed on a phenomenological basis and was later derived from a microscopic microswimmer model [30,31]. It can be regarded as the combination of a (simplified) Toner-Tu model [46,47] and a fourth-order term as in the Swift-Hohenberg equation [58]. The timedependent polarization density evolves according to The density is assumed to be constant, which leads to the incompressibility condition Eq. (4b) and introduces the Lagrange multiplier P enforcing this condition. Rewriting model (4) in potential form with shows that the dynamics of p is governed by pure relaxational dynamics derived from F and a convective part λ 0 (p·∇)p. For the bulk terms in F we distinguish between two regimes: For A > 0 and C > 0 these terms stabilize a disordered state. For A < 0 and C > 0 they represent a double-well potential, forcing a nonzero magnitude of the polarization density. Physically, the rotational symmetry is broken spontaneously and a globally ordered state with |p| = −A/C is stable in this case. Hence, for Γ 0 > 0 and Γ 2 = 0 the model (4) reduces to the Toner-Tu theory. However, if Γ 0 < 0 and Γ 2 > 0 a finite-wavelength instability is introduced, similar as in the Swift-Hohenberg equation [58]. This instability destabilizes both the disordered and globally ordered state simultaneously, leading to meso-scale turbulence. Setting Γ 0 < 0 is justified by physical arguments. Indeed, model (4) can be derived from microscopic considerations using a coarse-graining procedure [30,31], which indeed leads to Γ 0 < 0 due to activity and hydrodynamics. The competition between alignment and hydrodynamics sets the effective length scale of the evolving pattern. In this derivation it also becomes apparent that the coefficients in model (4) should, in general, also depend on density.

Extended model
Two major assumptions underlie model (4): That density is constant and that the particles propel with constant speed v 0 along their orientation. We relax these assumptions and replace them by expressions motivated by model (1). First, we assume that the velocity is density-dependent, i.e. we replace the constant speed v 0 with a density-dependent speed v(ρ). Such an assumption is very natural in realistic active matter systems as excluded volume as well as collective effects lead to a density-dependent speed [11,32,41,43]. Next, we have to replace the incompressibility condition by an evolution equation for the density. A natural choice is a continuity equation consisting of the divergence of the mass flux and a diffusive term Moreover, this equation agrees with the ones derived by [43,44] for MIPS (see Eq. (1a)) as well as the one derived for model (4) (in a suitable defined limit).
In model (4) the coupling to the (degenerate) density equation is accomplished via the Lagrange multiplier P acting as a pressure. As we replace the incompressibility condition Eq. (4b) with the continuity equation Eq. (7) an explicit coupling in terms of the density is needed. We choose Such a term appears naturally when coarse-graining microscopic models that incorporate self-propulsion. Indeed, this term is reported for all comparable systems we are aware of [43,44,33,34,52,59,60], see also Eq. (1b). While the details of the underlying microscopic model and coarse-graining procedure (especially the choice of an appropriate closure) might introduce additional terms to the dynamics of the polarization density p, a term as in Eq. (8) will always be present. Secondly, one can think of Eq. (8) as a low order approximation of the pressure, disregarding higher order coupling between ρ and p. We note that Eq. (8) can only formally be regarded as a pressure. Determining the pressure of active fluids is in general a complicated task [61,62,63], especially when accounting for hydrodynamic interactions. Alternatively, a virial expansion [64] or a treatment as in the Toner-Tu theory [46,47] is possible. As briefly mentioned in section 2.1, all coefficients of model (4) are, in principle, dependent on the density. An appropriate rescaling reduces the possible density-dependent parameters to λ 0 , Γ 0 and A. Experimental and numerical findings agree that the characteristic length scale set by Γ 0 in the turbulent regime does not depend on the (overall) density [65]. Hence, we drop that dependency. Furthermore, for simplicity we assume λ 0 to be independent of density as well. Note though that earlier studies suggest a non-monotone dependence of λ 0 on ρ [30]. This leaves A as the only density-dependent parameter. In fact, the transition from a dilute, disordered state to a dense, globally ordered state in the Vicsek model can be explained by a change of sign of A through an increased density. Altogether, our phenomenological model in its most general form is given by Overall, our model is essentially a minimal model that incorporates the main features of the models presented in the previous section. Moreover, the three instabilities can be tuned independently through the coupling terms and Γ 0 .

Stability Analysis
As a first insight into the dynamics expected from model (9) we perform a linear stability analysis for the steady states of the model. Clearly, a trivial steady state is given by (ρ, p) = (ρ 0 , p 0 ), where ρ 0 and p 0 are uniform in space. In the following we distinguish between the case p 0 = 0 (disorder) and |p 0 | > 0 (global polar order). Due to the special structure of the stability matrix (see Appendix A for details), the dispersion relations for the disordered state can be computed analytically. This is also possible for the polar state. However, finding the dispersion relation, i.e. solving for the eigenvalues of the stability matrix, produces lengthy expressions offering little insight. Hence, we only present numerical results for the polar state.

Disordered State
Linearizing Eq. (9) around the steady state (ρ, p) ≡ (ρ 0 , 0) and expanding perturbations into Fourier modes reveals the dispersion relations where k = |k| is the magnitude of the wavevector k = (k 1 , k 2 ). See Appendix A for details. Furthermore, we introduced where the constant γ, quantifying the coupling to the density-dependent propulsion speed, is given by The first eigenvalue σ 1 (k) does not contain any coupling terms nor any contributions from the density equation. Furthermore, the corresponding eigenvector is given by which is independent of density. Hence, σ 1 (k) solely affects the stability of the polarization density p. Moreover, the same dispersion relation is found in model (4), see [23,29]. Therefore, model (9) inherits the finite-wavelength instability of the polarization for sufficiently small Γ 0 < 0 from model (4). This instability is characterized by a band of unstable modes bounded away from zero as can be seen in figure 1A. It is straightforward to compute the critical parameter from Eq. (10a) as The qualitative behavior of the other eigenvalues σ 2,3 (k) cannot be read off directly from Eq. (10b). Instead, performing a small wavenumber expansion reveals a long-wavelength instability for D(ρ 0 ) < 0. From this expression and Eq. (12) we can calculate critical parameters by specifying v(ρ). The resulting long-wavelength instability is pictured in figure 1B. As expected from our modeling approach, this instability is similar to the one reported for MIPS, see section 2.1 and [43,60]. Note that D(ρ 0 ) coincides with D(ρ) introduced in Eq.
(3) for the choice ρ = ρ 0 . While both instabilities have been studied in great detail separately, to the authors knowledge, a situation as depicted in figure 1C, where both instabilities are present at the same time, has not been studied yet. As we will show in the next section, this leads to interesting dynamics.

Polar State
There is an additional long-wavelength instability of the disordered state for A(ρ 0 ) < 0. As discussed in section 2.1, a steady state with polar order, i.e. (ρ, p) ≡ (ρ 0 , p 0 ) and |p 0 | = −A(ρ 0 )/C, emerges in this situation. The stability analysis of the polar state can be carried out similar to the disordered case. However, the resulting dispersion relations are lengthy and intricate, hampering an intuitive interpretation. Numerical computation of the eigenvalues reveals instabilities similar to the disordered state. There is a finite-wavelength instability for Γ 0 ≤ Γ c 0 perpendicular to p 0 (see figure 1D) and a long-wavelength instability above a critical coupling constant, also perpendicular to p 0 (see figure 1E). Finally, both instabilities can be present as in the disordered case, see figure 1F. Interestingly, the numerically computed critical values for both instabilities are the same as in the disordered case. Hence, the disordered and polar states loose stability simultaneously, indicating the existence of a new dynamical attractor.

Numerical solution of the model equations
We now explore the dynamics produced by model (9) numerically. For this purpose, we have to specify the coupling terms v(ρ) and A(ρ). As the model is complex, we aim for simple coupling terms in order to ease the numerical burdens and reduce the amount of possible parameters. Hence, we choose v(ρ) = v 0 − ζρ, i.e. we study the linear case of Eq. (2). A monotone decrease with density is motivated by crowding effects, i.e. self-propulsion is counteracted by steric hindrance in dense areas. The linear model Eq. (2) with coefficient ζ was discussed in [43,44,59,60] and derived as a first-order approximation from microscopic coniderations. Other monotonically decreasing functions of density have been studied in the literature as well (see [42] and [52] for an exponential or hyperbolic tangent dependence respectively).
In addition, to further simplify the analysis, we choose A(ρ) ≡ A > 0. While such a choice might seem arbitrary at first glance, complex non-equilibrium dynamics can be observed for sufficiently small Γ 0 , even for A > 0, due to local shear stresses. A profound analysis on the influence of the potential terms with coefficients A and C on the dynamics of model (4) can be found in [29]. Therein, the authors conclude that the main difference is an absence of jets for A > 0. From a general point of view, the choice A(ρ) ≡ A > 0 disregards the polar state and its effects on the dynamics, which allows us to focus on the interplay between meso-scale turbulence and MIPS and to reduce the amount of parameters by setting C = 0.
First, we will present a numerical phase portrait. As the model includes several parameters, we have to restrict ourself to a low-dimensional cut in parameter space. We study the model in the space spanned by Γ 0 and ζ. These parameters can be used to control the instabilities independently. Furthermore, we can compare critical values computed numerically with the ones found from the stability analysis in section 3.1. The phenomenological parameters can be related to physical properties by examining microscopic models. While Γ 0 is determined by the characteristics of the surrounding fluid and the activity (see [30,31]), the parameter ζ depends on the details of the repulsive interactions (see [44]). We then give a qualitative description of the dynamical phases encountered, when the different instabilities are present. Finally, we discuss the anomalous velocity statistics observed in more detail.

Phase Portrait
To obtain a phase portrait, we numerically solve Eq. (9) for slightly perturbed homogeneous initial conditions. The exact simulation setup can be found in Appendix B, details on the numerical implementation are provided in Appendix C. To distinguish phases numerically, we introduce two quantifiers: the enstrophy as a measure for the presence of vortices and the modality of the density distribution to detect clustering. Details can be found in Appendix B. Using these indicators, the phase portrait figure 2 is numerically calculated, where the different colors correspond to the phases. Phase boundaries expected from the linear stability analysis are depicted as yellow lines. They can be obtained by calculating critical parameters, which we already determined for Γ c 0 in Eq. (14). Similarly, the critical damping parameter ζ c can be obtained by using Eq. (15) and Eq. (2), giving where the last approximation holds, if the product A(ρ 0 )D is small.

Qualitative phase descriptions
Model (9) produces a wealth of new dynamics, which cannot be covered entirely within this work. We therefore only give a qualitative description of the phases characterized by our coarse numerical measures.

Disordered State
We start our analysis of the phase portrait around the disordered state (ρ, p 0 ) = (ρ 0 , 0). As expected from the stability analysis, the disordered state is stable for ζ < ζ c and Γ 0 > Γ c 0 , which encompasses region D of the phase diagram figure 2. We do not include snapshots of the dynamics in figure 3, since there are no notable dynamics or features to report.
Isotropic Turbulence For ζ ζ c and Γ 0 < Γ c 0 Isotropic Turbulence is observed (region IT in figure 2). This state is governed by vortex-like structures which split and merge but exhibit a characteristic length scale, see figure 3A. This length scale can be obtained by the dominant mode of Eq. (10a), which reveals Λ = 2π −2Γ 2 /Γ 0 . Numerically, the length scale manifests itself as a dip in the spatial (time averaged) velocity correlation function and as a peak in the power spectrum (not plotted here, see [23,29]) as these quantities are linked by the Wiener-Khinchin theorem [66]. The density stays almost constant throughout the simulation (narrow distribution around ρ 0 , see figure 3A and Appendix B). As briefly discussed in the introduction of this section, we label this state as Isotropic Turbulence (IT) to account for the absence of jets commonly found in bacterial turbulence, see [29].

Motility-Induced Clustering and Motility-Induced Phase Separation
Taking ζ ≥ ζ c and Γ 0 Γ c 0 leads to Motility-Induced Clustering (MIC) and Motility-Induced Phase Separation (MIPS), to be found in region MIC of the phase portrait. The dynamics are characterized by the emergence of dense clusters with v = 0 surrounded by a dilute phase, see figure 3B. We term the generic case MIC, but refer to MIPS when cluster coarsen over time, eventually reaching a completely phase-separated state, see figure 3D. In MIPS, clusters have an almost perfect spherical shape and their number decreases monotonically.
Isotropic Turbulence with Clustering In the lower right corner of the phase portrait (region ITC), a combination of the two states discussed previously is encountered. To be precise, we observe a phase separation into a dense phase with v = 0 and a dilute phase with v = 0. The dilute phase shows dynamics similar to isotropic turbulence. That is, we encounter vortices with a characteristic length scale in the dilute part of the simulation domain. Snapshots can be seen in 3C and E. Note that the interfaces between dilute and dense phases are highly irregular. Furthermore, we observe fluctuations in the number of clusters and their shape.
Most of the dynamics discussed previously can be expected from the linear stability analysis. However, there are two noteworthy exceptions which we label as regions ITC 2 and ITC 3 in figure 2. To understand the dynamics in these regimes, we have to study the nucleation and coarsening processes in more detail.

Nucleation through turbulence
Insights into region ITC 2 can be gained by studying the nucleation of clusters. Nucleation and coarsening in MIPS is well studied in the literature, see for ex-ample [59,60,67,68,69,70,71]. A central result is that these processes in MIPS are quite similar to passive gas-liquid phase separation. As this classical transition is known to be a first-order phase transition (with the critical point being a notable exception) the same applies for MIPS, see [72,73] for an extensive study. Accordingly, clustering and phase separation can occur via two different mechanisms: Either nucleation and growth (in the metastable region) or spinodal decomposition (in the spinodal region). In the latter case there is no energy barrier to form a new phase. Hence, small perturbations start growing almost instantly. Therefore, the boundaries of the spinodal region, also known as the spinodal, can be detected by a linear stability analysis. In our case this corresponds to vertical yellow line in figure 2. To the right of that line we observe spinodal decomposition, either with or without the presence of turbulence, see figure 4A. However, the spinodal region is accompanied by a metastable region. There, perturbations have to overcome an energy barrier to form nuclei. In classic gasliquid phase separation this energy barrier can be overcome over time by thermal fluctuations. As this is not possible in our model, we do not observe nucleation and growth in the absence of turbulence. In contrast, the turbulent motion for Γ 0 < Γ c 0 leads to local density inhomogeneities, which can overcome the energy barrier and trigger nucleation. This is what is found in region ITC 2 and illustrated in figure 4B. At t = 4, no clustering is observed, whereas for spinodal decomposition clusters appear throughout the entire simulation domain, see figure 4A. Turbulence sets in at approximately t ≈ 20, which leads to the formation of a few nuclei at random sites, see t = 24. We want to point out that this is a surprising result. Since inertial turbulence is associate with mixing, one would expect turbulence to inhibit nucleation.
This mechanism offers a possible explanation for the shape of the phase boundary between regions IT and ITC 2 . Clearly, nucleation and growth are only possible in the metastable regime. The extent of the metastable region can be numerically estimated by checking for hysteresis. We initiate simulation runs either with a homogeneous density profile or with a completely phaseseparated state consisting of a single droplet (in the non-turbulent regime). The lower boundary of the metastable region is reached when the droplet looses stability. From our simulations this point can be estimated to be at ζ b ≈ 0.68. The upper boundary is provided by the spinodal at ζ c = 1.25. While checking for hysteresis, we made another important observation. Close to the spinodal a droplet with slightly higher density than the average density ρ 0 is sufficient to trigger nucleation. At the binodal, only droplets with a density close enough to the maximal density are stable. Additionally, the maximal density increases when decreasing ζ due to Eq. (2). Hence, when moving away from the spinodal, larger density inhomogeneities have to be provided to trigger nucleation. Furthermore, from our simulations we observe that the amplitude of density variations, i.e. ρ max − ρ min , is proportional to Γ c 0 − Γ 0 , see Appendix E. Hence, for Γ 0 ≈ Γ c 0 , only small density inhomogenities are observed. These are enough to trigger nucleation close to the spinodal, but are not sufficient close to the binodal. Altogether, these observations explain the shape of the phase boundary in the metastable regime of the phase portrait.

Altered Ostwald Ripening
We will now focus on region ITC 3 of figure 2, i.e. we want to answer the question why we observe enstrophy above the critical value for the onset of turbulence Γ c 0 . To explain this, we focus on the coarsening kinetics. Coarsening in MIPS is similar to Ostwald ripening, see [42,67,68]. That is, clusters grow on the expense of smaller clusters, which dissolve and redeposit onto larger ones. What that process typically looks like for classical MIPS (model (1)) is shown in figure  5A. During the dissolution process, the mass flux streamlines (which indicate the direction of mass transport) are perpendicular to the dissolving surface. Almost immediately after the cluster disappears, the streamlines rearrange, leaving no trace of the dissolved cluster, see t = 19.
However, the presence of the convective term λ 0 (p · ∇)p in our simulations (see model (9)) seems to alter the dissolution dynamics. The snapshots in figure  5A are produced by setting λ 0 = 0, whereas the snapshots in figure 5B are obtained for λ 0 = 3 (all other parameters are unchanged). For the latter choice, the streamlines appear tangential to the cluster surface, leading to a vortex after the cluster disappears. As λ 0 = 3 is the parameter used to compute the phase diagram, this behaviour is observed in the entire region ITC. In region ITC 3 (as opposed to regions ITC 1 and ITC 2 ), these vortices are transient, i.e. they vanish some time ∆t after the cluster has dissolved. This is expected as the finite-wavelength instability is inactive in this regime, i.e. turbulence is not self-sustained. However, as other clusters dissolve, new vortices are formed constantly, leading to an increased enstrophy over a long time, see figure 10 in Appendix E.
While focusing on coarsening dynamics, we made another interesting observation concerning the phase portrait: For some simulation runs, there appears to be no coarsening (at least on the simulation time scale). That is, the number of clusters does not decrease nor does the size of the largest cluster increase, see Appendix F for details. The dynamics for these cases are shown in figure 3B and C without and with the presence of turbulence respectively.
Altogether, the linear stability analysis provides a viable intuition into the expected phases. However, it fails to cover metastable regimes and nonlinear effects. Nevertheless, these effects and the apparent arrest of coarsening pose an intriguing research opportunity, which we will pursue in the future.

Anomalous velocity statistics
The main motivation to study the combined model (9) was to observe and explain anomalous velocity statistics. To quantify the deviations from Gaussian statistics, we compute the kurtosis κ of the standardized velocitiesv where v and σ(v i ) are the mean and the standard deviation of the velocity component v i respectively. The kurtosis coincides with the fourth moment for a standardized distribution. In the lower left corner of figure 6 (almost) normal statistics are observed, i.e. a kurtosis κ ≈ 3 is reported. This is not surprising since this region coincides with the turbulent regime (region IT in figure 2), where we expect normal statistics as in the incompressible model (4). However, strongly anomalous statistics (κ ≥ 6) are reported for a large part of the phase portrait. The anomalous region matches with the clustering phases (ITC and MIC) of figure 2, indicating that the deviation from normal statistics stems from clustering. Figure 6: Kurtosis in the ζ-Γ 0 plane. For ζ < ζ c and Γ 0 > Γ c 0 no velocity statistics are computed as v ≡ 0 (up to numerical errors). Simulation runs with a kurtosis higher than 6 are yellow.
Indeed, that hypothesis is supported by investigating the velocity statistics in more detail. Compared to the (almost) normal statistics in the IT regime, a clear peak around zero is visible in the ITC phase in figure 7A. This peak can be explained by a subpopulation argument: Computing the velocity statistics in the dense clusters and dilute, turbulent regimes separately reveals a clear split, see figure 7B. Details on the thresholding can be found in Appendix F. As expected, the statistics in the dilute part are similar to the ones reported for IT, whereas in the dense phase v ≈ 0. Combining both distributions gives the blue curve in figure 7A. Note that the statistics in the dilute phase are not perfectly Gaussian. The offset could be due to the interface between phases and the interactions of turbulence and clustering. Moreover, the statistics in the IT phase already show a slight deviation from normal statistics, i.e. they are only approximately Gaussian. Furthermore, we want to point out that figure 7B was produced in the spinodal regime (region ITC 1 ), i.e. when clusters form due to spinodal decomposition. The subpopulation argument works worse in the metastable regime, i.e. when nucleation is triggered by turbulence (region ITC 2 ). We speculate that in this case the interactions between clustering and turbulence are pronounced, leading to stronger correlations. Also note that the transition from normal to anomalous statistics gets smeared out for lower values of ζ. This hints at a possible continuous transition at the nucleation point. Altogether, the statistical properties of the system in the metastable regime require a more fundamental analysis and deeper understanding of the interactions in the competing processes.

Conclusion, Discussion and Outlook
We have presented and studied a phenomenological model combining ideas from MIPS and a fourth-order theory proposed to describe active turbulence. The underlying theories are distinct in the type of main instability they describe. While MIPS is driven by a global (long-wavelength) instability, meso-scale turbulence is characterized by a specific length scale (short-wavelength instability). Our model inherits both instabilities, which results in rich dynamics and an interesting interplay between the two.
While we showed that turbulence can trigger nucleation, the effect of turbulence on the coarsening process is not yet clear. Important alterations concerning the shape of the phase boundaries, fluctuations and cluster fluidity are evident from our simulations. However, these aspects need to be investigated further to quantify their importance and physical origin. Furthermore, we observed situations where coarsening appears to be frozen on a certain length scale. As this happens also in the absence of turbulence, we speculate that such an arrested phase separation might be connected to the loss of pure relaxation dynamics by the convective term.
Spatial segregation into static clusters surrounded by swarming bacteria was reported for colonies of Bacillus subtilis when adding (sublethal concentrations of) antibiotics in [40]. Motile cells are diluted, resulting in a lower density of swarming bacteria. Our simulations show qualitatively similar results as clusters with v = 0 are surrounded by a dilute, turbulent phase. Moreover, we observe anomalous statistics with a kurtosis up to 6 or larger, which was also reported in [40]. Furthermore, a recent experimental study [74] indicates that the coexistence of swarming dynamics and MIPS might explain the stressinduced transition to biofilm.
Connecting phenomenological transport coefficients with experimentally traceable parameters is not straightforward. The derivation of model (4) presented in [30,31] can give a hint for most parameters. However, obtaining the exact form of v(ρ) for biological system from microscopic considerations is fairly complicated. Alternatively, the local v(ρ) can be fitted to (global) experimental data reported in [32,75] for example. This data suggests a non-monotone dependence of velocity on density, in contrast to the simple linear assumption we used in our numerical study. Nevertheless, preliminary simulations show that our general results still hold for more complicated functions v(ρ).
Our phenomenological approach allows us to understand anomalous statistics in meso-scale turbulent systems without specifying the microscopic details. Hence, our model might be applicable for different experimental setups. For example, anomalous velocity statistics were reported in [39] (different aspect ratio of mutated Bacillus subtilis) and [57] (monolayer swarming of Bacillus subtilis for low density). However, different mechanisms could be underlying these observations. For example, a chemical response could lead to velocity variations while leaving the density unchanged. Such a situation could be possibly modelled by allowing v(ρ) to depend on an external scalar field instead of the density.
Altogether, our results provide a simple explanation for how anomalous statistics can arise in meso-scale turbulence. While addressing that topic, new questions are raised concerning the role of turbulence for nucleation and coarsening. This poses a challenging research opportunity we will pursue in the future.

A Stability Analysis
The stability analysis of the system (9) can be done in quite general fashion, i.e. without specifying the exact from of the coupling terms v(ρ) and A(ρ) nor the steady state. Introducing perturbations p → p 0 + δp and ρ → ρ 0 + δρ and linearizing Eq. (9) close to the steady state (ρ, p) ≡ (ρ 0 , p 0 ) yields the equations of motion for the perturbations Expanding the perturbations in Fourier modes with k = (k 1 , k 2 ), i.e. setting leads to the coupled dispersion relations where we grouped terms belonging toρ andp. Hence, it is straightforward to read off the stability matrix from Eq. (19). Roughly speaking, a finitewavelength instability is expected for sufficiently negative values of Γ 0 and a long-wavelength instability for A(ρ 0 ) < 0 (disordered case). Another longwavelength instability can be found in a similar way as in section 2.1 for MIPS: Inserting the first coupling term from the equation forp into the equation for ρ results in an effective diffusion equation for the density similar to Eq. (3).
For v (ρ 0 ) < 0 sufficiently small, the sign of the diffusion coefficient changes, triggering a long-wavelength instability. As mentioned in the main text, solving for the eigenvalues of the stability matrix is possible, but yields complicated expressions which are not instructive. However, assuming a disordered state reduces the complexity of the calculation significantly. Rewriting Eq. (19) in matrix notation for the disordered state p 0 = 0 leads to the equation with stability matrix with coupling terms α, β given by Note that αβ = γ for γ introduced in Eq. (12). The stability matrix Eq. (21) contains a scalar multiple of the identity as a block matrix in the lower right corner. This makes the computation of the eigenvalues much handier as compared to the polar case and leads to the dispersion relations Eq. (10).

B Phase Identifiers
The phase diagram figure 2 is produced by solving Eq. (9) and assigning a phase based on two quantifiers: the enstrophy of the polarization density and the modality of the density distribution. Runs are initiated with small perturbations around the overall density ρ 0 and a perturbed smooth polarization field with low amplitude.
State Ω ρ (pdf) Disordered ≈ 0 unimodal Isotropic Turbulence 1 unimodal Motility-Induced Clustering ≈ 0 multimodal Isotropic Turbulence with Clustering 1 multimodal The enstrophy is defined as Ω = 1 2 |ω(x, t)| 2 , where ω is the vorticity field ω = ∂ x v y − ∂ y v x of some vector field v = (v x , v y ). Brackets and overbars denote spatial and temporal averages respectively. We compute the enstrophy of the polarization density p in the dilute phase as a measure for the presence of vortices. The density distribution is computed after an initial transient and checked for bimodality. The transient is needed to filter out the initial (almost homogeneous) density field. Results are shown in figure 9 for the five cases depicted in figure 3.

C Numerical Implementation
We solve the system (9) numerically using a pseudo-spectral method with operator splitting for time integration [76]. Nonlinearities are treated by applying a 2/3 dealias rule [77]. While pseudo-spectral methods have been proven to be powerful and reliable for incompressible flows, special care is needed when dealing with conservation laws [78]. The main culprit lies in sharp boundaries between phases. Those are ill-suited for a global method and result in numerical errors in the form of Gibbs phenomena [79,80]. This necessitates very fine grids. As we are studying coarsening kinetics (long time scale, large system size) finite computation resources and times restrict the resolution of the grid. A more feasible and numerically efficient approach is choosing the diffusion coefficient of the density equation big enough to sufficiently smooth out phase boundaries. This allows us to choose grid sizes ranging from 128x128 to 512x512 and time steps of ∆t = 10 −3 . . . 10 −1 . Alternatively, advanced techniques [81,82] could be applied to increase stability and accuracy.  Table 1: Simulation parameters for figure 1.   Table 2: Simulation parameters for figures 2, 3, 6, 7, 9 and 10.

D Simulation Parameters
The subplots in figure 3 are generated according to table 3. Figure 4 A and B show runs depicted in figure 3 D and C at earlier times respectively. Furthermore, the data for figure 7 is obtained from the simulation runs shown in figure 3 A and D for IT and ITC.   E Additional figures for region ITC 2 and ITC 3 In our simulations we observe the relationship ρ max − ρ min ∼ Γ c 0 − Γ 0 between the amplitude of density variations and the distance to the critical parameter for the onset of turbulence Γ c 0 . This can be deduced from figure 10B. The values ρ min and ρ max are computed for every simulation run over the whole time and space domain, i.e. ρ min = min ρ(x, t) and ρ max = max ρ(x, t). Furthermore, we plot the time-dependent enstrophy Ω(t) for runs in region ITC 1 and ITC 3 in figure 10A. This shows that in region ITC 3 meso-scale turbulence is not self-sustained as Ω(t) ≈ 0 for some t, while Ω(t) 0 in region ITC 1 .

F Coarsening
We track the number of dense clusters and the size of the largest one to examine coarsening. The detection of a cluster is done as follows: We define a density threshold as where the maximal density is given by ρ max = v 0 /ζ by Eq. (2). This threshold is sufficient if the interfaces are negligible. After each snapshot is labeled according to this threshold, we use methods from image recognition [83] to count and characterize clusters where special care is taken to account for the periodic boundary conditions. From this data linear trends are fitted for the number and mass fraction (mass of the cluster compared to the overall mass in the computation box). Since we are interested in the long time dynamics and not the onset of clustering, we start to fit after a transient of five time units since the emergence of the first cluster. We define a relatively generous criterion for coarsening by requiring non increasing numbers of clusters and an increase in mass of the largest cluster (applied to the linear fits).