Hydrodynamic instabilities in a two-dimensional sheet of microswimmers embedded in a three-dimensional fluid

Abstract

A collection of microswimmers immersed in an incompressible fluid is characterised by strong interactions due to the long-range nature of the hydrodynamic fields generated by individual organisms. As a result, suspensions of rear-actuated ‘pusher’ swimmers such as bacteria exhibit a collective motion state often referred to as ‘bacterial turbulence’, characterised by large-scale chaotic flows. The onset of collective motion in pusher suspensions is classically understood within the framework of mean-field kinetic theories for dipolar swimmers. In bulk two and three dimensions, the theory predicts that the instability leading to bacterial turbulence is due to mutual swimmer reorientation and sets in at the largest length scale available to the suspension. Here, we construct a similar kinetic theory for the case of a dipolar microswimmer suspension restricted to a two-dimensional plane embedded in a three-dimensional incompressible fluid. This setting qualitatively mimics the effect of swimming close to a two-dimensional interface. We show that the in-plane flow fields are effectively compressible in spite of the incompressibility of the three-dimensional bulk fluid, and that microswimmers on average act as sources (pushers) or sinks (pullers). We analyse the stability of the homogeneous and isotropic state, and find two types of instability that are qualitatively different from the bulk, three-dimensional case: first, we show that the analogue of the orientational pusher instability leading to bacterial turbulence in bulk systems instead occurs at the smallest length scale available to the system. Second, an instability associated with density variations arises in puller suspensions as a generic consequence of the effective in-plane compressibility. Given these qualitative differences with respect to the standard bulk setting, we conclude that confinement can have a crucial role in determining the collective behaviour of microswimmer suspensions.

1. Introduction

To move through viscous fluids, motile micro-organisms exert forces and torques on their environment. Due to their small size, the resulting motion is dominated by the viscous stresses in the fluid, and is well described by the Stokes equation (Lauga & Powers 2009). The velocity fields created by micro-organisms in this regime can thus be described by a combination of the singular solutions of the Stokes equation, which decay algebraically in space (Spagnolie & Lauga 2012). The leading singularity for force- and torque-free microswimmers suspended in an infinite viscous fluid is given by a force dipole, which exhibits slow spatial decay, $r^{-d+1}$, where $r$ is the distance from the microorganism and $d$ is the dimensionality of space. Such velocity fields are long ranged, and even rather dilute suspensions of microorganisms can therefore exhibit significant fluid motion. Such ‘self-stirring’ in the absence of an external forcing is a strongly non-equilibrium phenomenon characteristic of active matter (Marchetti et al. 2013) and results in a unique set of transport and mechanical properties of microswimmer suspensions, such as enhanced diffusivity of tracer particles (Wu & Libchaber 2000; Leptos et al. 2009; Miño et al. 2013; Jepson et al. 2013) and significant changes of the apparent suspension viscosity (Rafaï, Jibuti & Peyla 2010; López et al. 2015; Saintillan 2018; Martinez et al. 2020).

Depending on the symmetry of their hydrodynamic flow, microswimmers can be divided into two distinct classes: (i) ‘pushers’, which accurately describes most swimming bacteria such as E. coli (Drescher et al. 2011), and (ii) ‘pullers’, whose reversed hydrodynamic flow field is typically exemplified by the alga C. reinhardtii (Guasto, Johnson & Gollub 2010). While these two flow fields lead to equivalent statistical properties for very dilute microswimmer suspensions, where swimmers can be treated as effectively non-interacting, they lead to strikingly different forms of swimmer–swimmer correlations beyond this limit. Arguably, the most profound effect of such swimmer–swimmer correlations is the transition to large-scale collective motion in elongated pusher-like microswimmers such as bacteria, often referred to as bacterial turbulence (Wensink et al. 2012; Dunkel et al. 2013). The mechanism of the transition to bacterial turbulence is usually rationalised based on mean-field kinetic theories that take into account the presence of long-range hydrodynamic interactions between microswimmers (Saintillan & Shelley 2008a). For three-dimensional (3-D) suspensions of pusher bacteria, these theories predict the onset of large-scale flows above a well-defined critical microswimmer volume fraction, and identify mutual particle reorientation due to hydrodynamic interactions as the main mechanism behind the instability, while no such instability occurs in 3-D puller suspensions. Even far below the onset of collective motion, the long-ranged hydrodynamic interactions result in strong correlations between microswimmers, leading to an enhancement of the effective tracer diffusivity in pusher suspensions compared with puller suspensions at the same density (Stenhammar et al. 2017; Škultéty et al. 2020).

In practice, observations of 3-D bulk collective flows are complicated by the unavoidable presence of boundaries in experimental set-ups used to study motile microorganisms. Importantly, motile organisms accumulate in close vicinity to solid (Berke et al. 2008) and liquid (Vladescu et al. 2014) boundaries, significantly depleting the bulk in between. Therefore, from an experimental perspective, it is much more convenient to study the ensuing collective motion in a layer of microswimmers confined by either type of boundary (Zhang et al. 2010; Chen et al. 2012; Sokolov & Aranson 2012; Gachelin et al. 2014).

The theoretical description of confined microswimmer suspensions is significantly more involved than their bulk fluid counterparts discussed above, as it includes several different effects of the confinement. Firstly, the spatial decay of Stokesian singularities is modified by the presence of boundaries (Spagnolie & Lauga 2012), with the corresponding far-field fluid velocity being dominated by flow singularities other than the force dipole which is dominant in bulk fluids (Brotto et al. 2013; Mathijssen et al. 2016; Jeanneret, Pushkin & Polin 2019). Secondly, the interactions between microswimmers and the boundaries might depend on the precise shape of the microswimmers, their mechanism of propulsion, surface charge and roughness, etc. The corresponding theory would need to sufficiently resolve the near-field fluid velocity generated by the microswimmers and include appropriate steric interactions (Pessot, Löwen & Menzel 2018; Zantop & Stark 2020) and the above-mentioned effects. Such a theory would be very challenging and strongly dependent on the specific system under investigation.

An additional important physical effect inherent to confined systems, although rarely discussed explicitly, is the effective compressibility of the in-plane flow fields generated by a layer of microswimmers (Salbreux, Prost & Joanny 2009; Jülicher, Grill & Salbreux 2018; Maitra 2020; Huang et al. 2021; Maitra 2023). A similar effect was previously discussed in the context of hydrodynamically interacting colloidal particles under partial confinement (Bleibel et al. 2014; Bleibel, Domínguez & Oettel 2015, 2017). To illustrate the point, we consider a force dipole oriented parallel to a solid wall at a distance $h$ from it. The fluid velocity component $\boldsymbol {u}_{\parallel }$ parallel to the wall can be deduced from the image system developed by Blake (1971) for a point force next to a solid boundary with vanishing boundary condition $\boldsymbol {u}_{\parallel }(z=h) = 0$, and is given by (Spagnolie & Lauga 2012)

(1.1)$$\begin{gather} \boldsymbol{u}_{{\parallel}} (\boldsymbol{\mathcal{x}}) = \frac{\kappa}{8{\rm \pi}} \left( \frac{\boldsymbol{\mathcal{x}}}{|\boldsymbol{\mathcal{x}}|^3} \left[ 3 \frac{(\boldsymbol{\mathcal{x}}\boldsymbol{\cdot}\boldsymbol{\mathcal{p}})^2}{|\boldsymbol{\mathcal{x}}|^2} - 1 \right] + \frac{\boldsymbol{\mathcal{x}}}{R^3} \right.\nonumber\\ \left.- \frac{3\boldsymbol{\mathcal{x}} (\boldsymbol{\mathcal{x}}\boldsymbol{\cdot}\boldsymbol{\mathcal{p}})^2 + 6h^2 \left\{\boldsymbol{\mathcal{x}} + 2 \boldsymbol{\mathcal{p}} (\boldsymbol{\mathcal{x}}\boldsymbol{\cdot}\boldsymbol{\mathcal{p}})\right\}}{R^5} + \frac{30 h^2 (\boldsymbol{\mathcal{x}}\boldsymbol{\cdot}\boldsymbol{\mathcal{p}})^2 \boldsymbol{\mathcal{x}}}{R^7} \right). \end{gather}$$

Here, $\boldsymbol {\mathcal {p}}$ and $\boldsymbol {\mathcal {x}}$ are two-dimensional (2-D) vectors that lie in the plane parallel to the wall and denote respectively the dipole orientation and the point where the velocity is evaluated relative to the position of the swimmer. Furthermore, $R=\sqrt {|\boldsymbol {\mathcal {x}}|^2 + 4 h^2}$, and $\kappa$ is the strength of the dipole. Next, we calculate the total flux of the fluid through a circle of radius $X$ centred on the microswimmer and parallel to the wall, as shown in figure 1, yielding

(1.2)\begin{equation} \int_{|\boldsymbol{\mathcal{x}}|=X} \,\mathrm{d} \boldsymbol{\mathcal{x}}\boldsymbol{\cdot}\boldsymbol{u}_{{\parallel}} = \frac{\kappa}{8 X} \left[ 1 - \frac{X^4-10 X^2 h^2+64h^4}{\left(X^2+4h^2 \right)^{7/2}}X^3\right]. \end{equation}

Hence, the total flux through an arbitrary circle around a microswimmer is non-zero, indicating that, if we only consider the velocity components parallel to the wall, they represent an effectively compressible velocity field. Moreover, one can show that the prefactor multiplying $\kappa$ is strictly positive, and the sign of the flux is therefore determined by the sign of the dipolar strength: pushers with $\kappa >0$ on average correspond to hydrodynamic sources in the plane parallel to the wall, while pullers with $\kappa <0$ correspond to hydrodynamic sinks. Therefore, when averaged over all orientations, two pushers advect each other to maximise their mutual separation, while two pullers do the opposite. This argument illustrates that hydrodynamic interactions between dipolar microswimmers moving next to a boundary have a very different nature to their bulk counterparts, which has previously been shown to promote dynamic self-assembly and crystallisation in active-particle systems (Singh & Adhikari 2016; Thutupalli et al. 2018). We stress that this effective interaction does not correspond to extra forces or torques on the particles, which is a common assumption when describing ‘dry’ active particles moving on a frictional substrate (ten Hagen et al. 2015), as the particles immersed in a fluid are strictly force and torque free.

Figure 1. Schematic picture of the flow around a single pusher microswimmer restricted to a 2-D plane embedded in a 3-D fluid. Note that the net flow $\boldsymbol {u}_{\parallel }$ going through the dashed circle is non-zero; in the 2-D plane, pushers therefore, on average, act as fluid sources, while pullers act as effective sinks.

The aim of this work is to understand the effect of this in-plane compressibility on collective motion in a minimal setting. We consider a layer of microswimmers restricted to move in an infinitely thin, 2-D plane embedded in a 3-D bulk fluid in the absence of any boundaries, corresponding to taking $h\to \infty$ in the example above. We assume low number density of microswimmers and approximate their velocity fields by those of 3-D force dipoles. As can be seen from (1.2), even in this limit, the in-plane velocity field generated by the microswimmer is still effectively compressible. Therefore, although this specific set-up is difficult to realise experimentally, it allows us to single out the effect of in-plane compressibility without needing to deal with other effects due to the confinement discussed above.

Below, we demonstrate that a 2-D layer of microswimmers embedded in a 3-D fluid is unstable for both pushers and pullers. We show that pusher suspensions are prone to an orientational instability that sets in at the smallest length scale available to the system, while puller suspensions exhibit an instability associated with density variations as a generic consequence of the effective in-plane compressibility.

The paper is organised as follows: in § 2, we present a general kinetic theory describing a dilute suspension of microswimmers interacting through long-range velocity fields. In § 3 we review its predictions for the widely studied case of microswimmers suspended in a 3-D bulk fluid. In § 4 we report on linear stability of a 2-D layer of microswimmers embedded in a 3-D fluid, which constitutes the main novel results of our work. We conclude by summarising our main findings in § 5.

2. Mean-field kinetic theory for microswimmer suspensions

Mean-field kinetic theories of microswimmer suspensions have been extensively discussed in the literature (Saintillan & Shelley 2008a,b; Subramanian & Koch 2009; Hohenegger & Shelley 2010; Koch & Subramanian 2011; Saintillan & Shelley 2013; Krishnamurthy & Subramanian 2015), and here we give only a brief summary of their general set-up.

A suspension of $N$ microswimmers moving autonomously through a viscous fluid is described by the one-particle distribution function $\varPsi (\boldsymbol {x},\boldsymbol {p},t)$, which defines the instantaneous probability of finding a particle at a spatial position $\boldsymbol {x}$ with an orientation given by the unit vector $\boldsymbol {p}$. The distribution function is normalised such that

(2.1)\begin{equation} \int \,\mathrm{d}\kern0.06em \boldsymbol{x}\, \mathrm{d} \boldsymbol{p} \, \varPsi(\boldsymbol{x},\boldsymbol{p},t) = 1. \end{equation}

Its time evolution is assumed to be governed by the Smoluchowski equation

(2.2)\begin{equation} \partial_{t} \varPsi + \nabla^{\alpha} \{ \dot{x}^{\alpha} \varPsi \} + \partial^{\alpha} \{ \dot{p}^{\alpha} \varPsi \} ={-} \lambda \varPsi + \lambda \int \frac{ \mathrm{d} \boldsymbol{p} }{\varOmega_d} \, \varPsi, \end{equation}

where Greek superscripts denote Cartesian components of vectors, $\nabla ^{\alpha } = \partial /\partial x^{\alpha }$, $\partial ^{\alpha } = (\delta ^{\alpha \beta } - p^{\alpha }p^{\beta })\partial / \partial p^{\beta }$, and $\delta ^{\alpha \beta }$ denotes the Kronecker delta; $\varOmega _d$ is the surface area of a $d$-sphere of unit radius, where $d$ is the dimensionality of space, and $\lambda$ is the tumbling frequency, as further described below. The deterministic part of the single-swimmer dynamics is described by the following microscopic equations of motion:

(2.3)\begin{gather} \dot{x}^{\alpha} = v_{s} p^{\alpha} + \mathcal{U}^{\alpha}(\boldsymbol{x}), \end{gather}

(2.4)\begin{gather}\dot{p}^{\alpha} = (\delta^{\alpha\beta} - p^{\alpha}p^{\beta}) (\mathcal{W}^{\beta\gamma}(\boldsymbol{x}) + B \mathcal{E}^{\beta\gamma}(\boldsymbol{x})) p^{\gamma}, \end{gather}

where dots denote time derivatives. According to (2.3), each particle changes its spatial position due to its own swimming with a constant speed $v_s$ in the direction of its orientation $\boldsymbol {p}$ and due to advection by the fluid flow with velocity $\boldsymbol {\mathcal {U}}$ at its position $\boldsymbol {x}$, created by all other microswimmers. Microswimmer orientations change according to Jeffery's equation, (2.4), which describes particle rotation by the gradients of the fluid velocity at its position (Kim & Karrila 2005). Here, $\mathcal {W}^{\alpha \beta }(\boldsymbol {x}) = \tfrac {1}{2} (\nabla ^{\beta } \mathcal {U}^{\alpha }(\boldsymbol {x}) - \nabla ^{\alpha } \mathcal {U}^{\beta }(\boldsymbol {x})),$ and $\mathcal {E}^{\alpha \beta }(\boldsymbol {x}) = \tfrac {1}{2} (\nabla ^{\beta } \mathcal {U}^{\alpha }(\boldsymbol {x}) + \nabla ^{\alpha } \mathcal {U}^{\beta }(\boldsymbol {x}))$ are the vorticity and the rate-of-strain tensors. The Bretherton parameter $B$ encodes the shape of the particle (Kim & Karrila 2005), with $B=1$ and $B=0$ respectively corresponding to needle-like and spherical particles. Finally, each microswimmer randomly selects a new orientation (tumbles) with rate $\lambda$. This discrete stochastic process cannot be incorporated in the time evolution (2.3) and (2.4) in a straightforward manner, but is readily described within the kinetic theory (Subramanian & Koch 2009; Koch & Subramanian 2011), as given on the right-hand side of (2.2).

The final ingredient of the theory is provided by the relationship between the local velocity field and the one-particle distribution function, given by

(2.5)\begin{equation} \mathcal{U}^{\alpha}(\boldsymbol{x},t) = N \int \,\mathrm{d}\kern0.06em \boldsymbol{x}' \,\mathrm{d} \boldsymbol{p}' \, u^{\alpha}(\boldsymbol{x}-\boldsymbol{x}',\boldsymbol{p}') \varPsi(\boldsymbol{x}',\boldsymbol{p}',t). \end{equation}

Here, $\boldsymbol {u}(\boldsymbol {x}-\boldsymbol {x}',\boldsymbol {p}')$ is the microscopic velocity field created at $\boldsymbol {x}$ by a microswimmer located at $\boldsymbol {x}'$ with the orientation $\boldsymbol {p}'$. In the following, we approximate $\boldsymbol {u}$ by the dipolar field generated by two equal and opposite point forces applied to the fluid infinitesimally close to each other (Lauga & Powers 2009). While this is a good approximation for dilute suspensions of microswimmers in 3-D bulk systems (Lauga & Powers 2009; Spagnolie & Lauga 2012; Škultéty et al. 2020), its validity in the context of the present study is discussed in § 5.

2.1. Linear stability of the homogeneous and isotropic state

The Smoluchowski equation (2.2) and the normalisation condition (2.1) admit as a solution the homogeneous and isotropic state $\varPsi _{HI} = 1/( \varOmega _d V_d )$, where $V_d$ is the $d$-dimensional volume of the suspension. The existence of this solution relies on the condition

(2.6)\begin{equation} \int \,\mathrm{d}\kern0.06em \boldsymbol{x}' \,\mathrm{d} \boldsymbol{p}' \, u^{\alpha}(\boldsymbol{x}-\boldsymbol{x}',\boldsymbol{p}') = 0, \end{equation}

which is true in all cases considered in this work. Notice that, while integrals over $\boldsymbol {x}'$ and $\boldsymbol {p}'$ vanish independently in 2-D and 3-D bulk suspensions, only the positional integral vanishes for microswimmers restricted to a 2-D plane in a 3-D fluid.

Stability of the homogeneous and isotropic state is determined by time evolution of small perturbations around $\varPsi _{HI}$. Techniques to study such problems are well established and have been extensively applied to linear stability analysis of the homogeneous and isotropic state for dilute 3-D suspensions of microswimmers (Saintillan & Shelley 2008a,b; Subramanian & Koch 2009; Hohenegger & Shelley 2010). Here, we use a somewhat different methodology (Stenhammar et al. 2017; Martinez et al. 2020), which is better adapted for our purpose.

First, we introduce a perturbation $\delta \varPsi (\boldsymbol {x},\boldsymbol {p},t)$ of the one-particle distribution function around the homogeneous and isotropic state, i.e.

(2.7)\begin{equation} \varPsi(\boldsymbol{x},\boldsymbol{p},t) = \frac{1}{\varOmega_d V_d} + \delta \varPsi(\boldsymbol{x},\boldsymbol{p},t), \end{equation}

where we assume that $\delta \varPsi (\boldsymbol {x},\boldsymbol {p},t)$ is small and that its integral over $\boldsymbol {x}$ and $\boldsymbol {p}$ vanishes. In what follows, we employ the $d$-dimensional Fourier transform defined through

(2.8a,b)\begin{equation} \hat{f}(\boldsymbol{k}) = \int \,\mathrm{d}\kern0.06em \boldsymbol{x} \, f(\boldsymbol{x}) \mathrm{e}^{-{\rm i}\boldsymbol{k}\,\boldsymbol{\cdot}\,\boldsymbol{x}}, \quad f(\boldsymbol{x}) = \frac{1}{(2{\rm \pi})^{d}} \int \,\mathrm{d} \boldsymbol{k}\ \hat{f}(\boldsymbol{k}) \mathrm{e}^{{\rm i}\boldsymbol{k}\,\boldsymbol{\cdot}\,\boldsymbol{x}}, \end{equation}

where $f$ is an arbitrary function of $\boldsymbol {x}$, and $\hat f$ denotes its Fourier space representation. The linearised Smoluchowski (2.2) in Fourier space reads

(2.9)\begin{equation} \left[ \partial_{t} + \lambda + {\rm i} v_{s} \boldsymbol{p}\boldsymbol{\cdot}\boldsymbol{k} \right] \delta \hat{\varPsi} = \frac{\lambda}{\varOmega_d} \delta \hat{\rho} + \frac{ n}{\varOmega_d} \left[{\rm d} B p^{\alpha} p^{\beta} - (1+B) \delta^{\alpha\beta} \right] {\rm i}k^{\alpha} \delta \hat{\mathcal{U}}^{\beta}, \end{equation}

where $n=N/V_d$ is the number density of particles, and we have defined the perturbations of the microswimmer density and fluid velocity as

(2.10)\begin{gather} \delta \hat{\rho}(\boldsymbol{k},t) = \int \,\mathrm{d} \boldsymbol{p}\ \delta \hat{\varPsi}(\boldsymbol{k},\boldsymbol{p},t), \end{gather}

(2.11)\begin{gather}\delta \hat{\mathcal{U}}^{\alpha}(\boldsymbol{k},t) = \int \,\mathrm{d} \boldsymbol{p}\ \hat{u}^{\alpha}(\boldsymbol{k},\boldsymbol{p}) \delta \hat{\varPsi}(\boldsymbol{k},\boldsymbol{p},t). \end{gather}

We proceed by assuming an exponential solution to (2.9) of the following form

(2.12)\begin{equation} \delta \hat{\varPsi}(\boldsymbol{k},\boldsymbol{p},t) = \delta \hat{\varPsi}(\boldsymbol{k},\boldsymbol{p}) \mathrm{e}^{\chi t}, \end{equation}

where the sign of the real part of the temporal eigenvalue $\chi$ determines the stability of the system against infinitesimal perturbations: for Re$[\chi ]<0$ the system is linearly stable, while for Re$[\chi ]>0$ the system is linearly unstable. It should be noted that a more general approach to the temporal linear stability analysis is to treat the problem as an initial value one and obtain the solution via the Laplace transformation, as was done in Stenhammar et al. (2017) and Škultéty et al. (2020). Working with the ansatz (2.12) is mathematically simpler but may not necessarily produce a dispersion law for all parameter values, as we discuss further below.

Inserting (2.12) into (2.9), we derive the following closed set of equations for the density and velocity perturbations:

(2.13)\begin{gather} \delta \hat{\rho} = {\rm i} n \delta \hat{\mathcal{U}}^{\alpha} \int \frac{\mathrm{d} \boldsymbol{p}}{\varOmega_d} \frac{ B d p^{\alpha} \boldsymbol{p}\boldsymbol{\cdot}\boldsymbol{k} - (1+B)k^{\alpha} }{\mathcal{L}(\chi,\boldsymbol{p}\boldsymbol{\cdot}\boldsymbol{k})} + \delta \hat{\rho} \int \frac{\mathrm{d} \boldsymbol{p}}{\varOmega_d} \frac{\lambda}{\mathcal{L}(\chi,\boldsymbol{p}\boldsymbol{\cdot}\boldsymbol{k})}, \end{gather}

(2.14)\begin{gather}\delta \hat{\mathcal{U}}^{\alpha} = {\rm i} n \delta \hat{\mathcal{U}}^{\beta} \int \frac{\mathrm{d} \boldsymbol{p}}{\varOmega_d} \frac{ \hat{u}^{\alpha}(\boldsymbol{k},\boldsymbol{p}) [B d p^{\beta} \boldsymbol{p}\boldsymbol{\cdot}\boldsymbol{k} - (1+B)k^{\beta}] }{\mathcal{L}(\chi,\boldsymbol{p}\boldsymbol{\cdot}\boldsymbol{k})} + \delta \hat{\rho} \int \frac{\mathrm{d} \boldsymbol{p}}{\varOmega_d} \frac{\lambda \hat{u}^{\alpha}(\boldsymbol{k},\boldsymbol{p})}{ \mathcal{L}(\chi,\boldsymbol{p}\boldsymbol{\cdot}\boldsymbol{k}) }, \end{gather}

where $\mathcal {L}(\chi,\boldsymbol {p}\boldsymbol {\cdot }\boldsymbol {k}) = \chi + \lambda + {\rm i} v_s \boldsymbol {p}\boldsymbol {\cdot }\boldsymbol {k}$. The solution to this eigenvalue problem depends on the precise form of the microscopic velocity field $\hat {\boldsymbol {u}}(\boldsymbol {k},\boldsymbol {p})$, which, in turn, depends on the dimensionality of space and the dimensionality of the vectors $\boldsymbol {p}$ and $\boldsymbol {k}$. In the following, we begin by revisiting the known results on linear stability of bulk 3-D suspensions (Saintillan & Shelley 2008a,b) before turning our attention to a 2-D layer of microswimmers embedded in a 3-D fluid.

3. Microswimmers in an unbounded bulk suspension

In this section, we consider the case of microswimmers that move freely in a 3-D space filled with an incompressible fluid. The dipolar velocity field for $d=3$ is given by (Lauga & Powers 2009)

(3.1)\begin{equation} u_{3d}^{\alpha}(\boldsymbol{x},\boldsymbol{p}) = \frac{\kappa}{8{\rm \pi}} \frac{x^{\alpha}}{|\boldsymbol{x}|^{3}} \left[ 3 \frac{(\boldsymbol{x}\boldsymbol{\cdot}\boldsymbol{p})^{2}}{|\boldsymbol{x}|^{2}} - 1 \right], \end{equation}

and, as shown in Appendix A, its Fourier transform is

(3.2)\begin{equation} \hat{u}_{3d}^{\alpha}(\boldsymbol{k},\boldsymbol{p}) ={-} {\rm i} \kappa \frac{\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{p}}{k^{2}} \mathbb{P}^{\alpha\beta} p^{\beta}. \end{equation}

Here, $\kappa$ is the dipolar strength and $\mathbb {P}^{\alpha \beta } = \delta ^{\alpha \beta } - k^{\alpha } k^{\beta }/k^{2}$ is the transverse projection operator that ensures that (3.2) satisfies the incompressibility condition $k^\alpha \hat {u}_{3d}^\alpha = 0$. Setting $d = 3$ in (2.13)–(2.14), using the fact that $\delta \hat {\mathcal {U}}^{\alpha } = \mathbb {P}^{\alpha \beta } \delta \hat {\mathcal {U}}^{\beta }$ and assuming that $\chi$ does not depend on the orientation $\boldsymbol {p}$ (see further discussion below), we obtain the following set of equations: (Martinez et al. 2020)

(3.3)\begin{gather} \delta \hat{\rho} = \frac{\lambda}{v_{s} k} \arctan(b) \delta \hat{\rho}, \end{gather}

(3.4)\begin{gather}\delta \hat{\mathcal{U}}^{\alpha} = \frac{B n \kappa}{v_{s}k} \frac{b(3+2b^{2}) - 3 (b^{2} + 1) \arctan(b)}{2b^{4}} \delta \hat{\mathcal{U}}^{\alpha}, \end{gather}

where we have introduced the dimensionless parameter

(3.5)\begin{equation} b = \frac{v_{s}k}{\chi + \lambda}. \end{equation}

We note that (3.3) and (3.4) are decoupled, and can be studied independently. The former equation involves only $\delta \hat {\rho }$, and we thus refer to it as the density eigenvalue problem. The latter equation is known to lead to a long-wavelength instability associated with orientational degrees of freedom only (Saintillan & Shelley 2008a,b; Subramanian & Koch 2009; Hohenegger & Shelley 2010; Stenhammar et al. 2017), and we thus refer to it as the orientational eigenvalue problem. We will now study these two problems separately, while the analogous calculation for the case of an infinite 2-D suspension is carried out in Appendix B with qualitatively identical conclusions.

3.1. Bulk orientational instability

Equation (3.4) yields the following eigenvalue problem:

(3.6)\begin{equation} \gamma = \frac{5}{2} \frac{b(3+2b^{2}) - 3(b^{2} + 1) \arctan(b)}{b^{4}}, \quad \gamma \equiv \frac{5v_{s}k}{Bn\kappa}. \end{equation}

We now solve (3.5) for $\chi$ and use the definition of $\gamma$ from (3.6), to arrive at

(3.7)\begin{equation} \chi ={-} \lambda + \frac{1}{5} F(\gamma) B \kappa n, \quad F(\gamma) = \frac{\gamma}{b(\gamma)} . \end{equation}

Equation (3.6) can be inverted numerically to obtain $b(\gamma )$ or, equivalently, $F(\gamma )$, while an analytical approximation to $F(\gamma )$ was developed in Škultéty et al. (2020) and Martinez et al. (2020). In figure 2 we plot the results of our numerical evaluation of $F(\gamma )$, showing that (i) the real part of $F$ is positive and decreases monotonically with increasing $\gamma$, and (ii) no solution is found for $\gamma >\gamma ^*$, a fact that is discussed further below. Together with (3.7), this implies that puller suspensions ($\kappa < 0$) are always stable while the instability of pusher suspensions ($\kappa > 0$) sets in at the largest possible length scale, corresponding to $\gamma \rightarrow 0$ and $F(\gamma ) \rightarrow 1$. It follows from (3.7), and from imposing $\chi >0$, that pusher suspensions are unstable for densities larger than

(3.8)\begin{equation} n_{c} = \frac{5\lambda}{B\kappa}. \end{equation}

This latter result is exact, and can be alternatively derived directly from (3.4) by setting $k \rightarrow 0$.

Figure 2. The real (a) and imaginary (b) parts of the function $F(\gamma )$ from (3.7).

The origin of this instability is the mutual reorientation of microswimmers (Saintillan & Shelley 2008a,b; Subramanian & Koch 2009; Hohenegger & Shelley 2010; Stenhammar et al. 2017). It is thus unsurprising to observe that no instability is present for spherical particles ($B=0$). Furthermore, the instability condition (3.8) does not involve the swimming speed $v_s$, implying that the orientational instability persists even for shakers – microswimmers that do not self-propel, yet are capable of generating dipolar fields (Stenhammar et al. 2017).

We finally comment on the fact that no solution for $F(\gamma )$ is found for $\gamma >\gamma ^*\approx 2.8$. This means that a solution to the linearised kinetic equation (2.9) obeying the ansatz (2.12) ceases to exist for $\gamma >\gamma ^*$. The reason for this is that, in the process of deriving (3.3)–(3.4), we have assumed that $\chi$ is independent of the swimmer orientation $\boldsymbol {p}$. Hohenegger & Shelley (2010) performed a detailed analysis of the eigenvalue problem in this parameter range avoiding such an assumption, and concluded that the system is linearly stable for $\gamma > \gamma ^{*}$. The most unstable eigenvalue is thus given by (3.7), leading to the instability threshold (3.8). The same conclusion is reached by solving (2.9) using Laplace transform techniques, as was done in Stenhammar et al. (2017) and Škultéty et al. (2020).

3.2. Density fluctuations in bulk suspensions

Both puller and pusher suspensions are stable against density fluctuations. For wavenumbers $k \leq {\rm \pi}\lambda /(2v_{s})$ this can be readily shown from (3.3), yielding

(3.9)\begin{equation} \chi ={-} \lambda + \frac{v_{s}k}{\tan(v_{s}k/\lambda)}. \end{equation}

from which we conclude that $\chi <0$. The corresponding eigenstates represent spatial modulation of the microswimmer density, while their orientations remain isotropically distributed regardless of their spatial positions. Similar to the orientational instability, (3.3) does not allow one to determine stability of the suspension when $k > {\rm \pi}\lambda /(2v_{s})$, since, in this regime, (3.3) has no solution. As discussed in § 3.1, this is a limitation of the ansatz (2.12). In Appendix C.1, we analyse this problem in some detail by solving (2.9) using Laplace transform techniques and show that the system is indeed stable with respect to infinitesimal density perturbations for all $k$.

4. Microswimmers restricted to a 2-D plane

We now turn our attention to the main problem of this study – the case of a 2-D layer of microswimmers embedded in a 3-D bulk fluid. As discussed in § 1, this spatial arrangement of microswimmers leads to the in-plane fluid velocity field being effectively compressible, with pushers acting on average as sources, while pullers act as sinks. Here, we study the consequences of this effective compressibility on the onset and type of collective motion expected in this arrangement of microswimmers.

The Fourier representation of the in-plane velocity field created by a microswimmer is given by (see Appendix A for details)

(4.1)\begin{equation} u_{plane}^{\alpha}(\boldsymbol{\mathcal{k}},\boldsymbol{\mathcal{p}}) ={-} \frac{{\rm i}\kappa}{2} \frac{\boldsymbol{\mathcal{k}}\boldsymbol{\cdot}\boldsymbol{\mathcal{p}}}{\mathcal{k}} \left[ \mathbb{P}^{\alpha\beta} + \frac{1}{2} \mathbb{Q}^{\alpha\beta} \right] \mathcal{p}^{\beta}, \end{equation}

where $\kappa$ is the dipolar strength which has the same dimensions as in the bulk 3-D case, and $\boldsymbol {\mathcal {k}} = ( \mathcal {k}_x, \mathcal {k}_y)$, $\boldsymbol {\mathcal {p}} = ( \mathcal {p}_x, \mathcal {p}_y)$ are respectively the in-plane wave and orientation vectors. In writing (4.1), we explicitly separated the term proportional to the longitudinal projection operator $\mathbb {Q}^{\alpha \beta } = \mathcal {k}^{\alpha } \mathcal {k}^{\beta } / \mathcal {k}^{2}$ to stress the compressible nature of $\boldsymbol {u}_{plane}$. We now define the transverse (incompressible) and longitudinal (compressible) velocity perturbations by

(4.2a,b)\begin{equation} \delta \hat{\mathcal{U}}_{{\perp}}^{\alpha} = \mathbb{P}^{\alpha \beta} \delta \hat{\mathcal{U}}^{\beta}, \quad \delta \hat{\mathcal{U}}_{{\parallel}}^{\alpha} = \mathbb{Q}^{\alpha \beta} \delta \hat{\mathcal{U}}^{\beta}, \end{equation}

where $\delta \hat {\mathcal {U}}^{\alpha } = \delta \hat {\mathcal {U}}_{\perp }^{\alpha } + \delta \hat {\mathcal {U}}_{\parallel }^{\alpha }$. Substituting (4.1) for $\boldsymbol {u}_{plane}$ into (2.13)–(2.14) results in the following compact set of equations determining linear stability of the system:

(4.3)\begin{equation} \begin{pmatrix} M_{11}-1 & 0 & 0 \\ 0 & M_{22}-1 & {\rm i} \mathcal{k}^{\alpha} M_{23} \\ 0 & {\rm i} \mathcal{k}^{\alpha} M_{32} & M_{33}-1 \end{pmatrix} {\cdot} \begin{pmatrix} \delta \hat{\mathcal{U}}_{{\perp}}^{\alpha} \\ \delta \hat{\mathcal{U}}_{{\parallel}}^{\alpha} \\ \delta \hat{\rho} \end{pmatrix} = 0, \end{equation}

where

(4.4)\begin{equation} \left. \begin{gathered} M_{11} = \frac{\kappa n}{v_{s}} \frac{B \left(b^2-2 \sqrt{b^2+1}+2\right) }{2 b^3}, \quad M_{22} = \frac{\kappa n}{v_{s}} \frac{B b^2 + (b^{2} + 2B)\left( 1 - \sqrt{b^{2} + 1} \right)}{4 b^3 \sqrt{b^2+1}}, \\ M_{23} = \frac{\kappa\lambda}{ v_{s} \mathcal{k}^{2} } \frac{1}{4 b} \left( \frac{1}{\sqrt{b^2+1}}-1 \right) , \quad M_{32} = \frac{n}{ v_{s} \mathcal{k}} \frac{2 B \left(\sqrt{b^2+1}-1\right) -b^2 (B+1)}{b \sqrt{b^2+1}} , \\ M_{33} = \frac{\lambda}{\mathcal{k} v_{s}} \frac{b}{\sqrt{1+b^{2}}}, \end{gathered} \right\} \end{equation}

and

(4.5)\begin{equation} b = \frac{v_s \mathcal{k}}{\chi + \lambda}. \end{equation}

As is apparent from (4.3), $\delta \hat {\mathcal {U}}_{\perp }^{\alpha }$ is decoupled from $\delta \hat {\mathcal {U}}_{\parallel }^{\alpha }$ and $\delta \hat {\rho }$. This allows us to analyse the two subsets of equations independently; as we show below, they correspond to different types of instabilities.

4.1. Orientational instability

We start by considering the transverse mode $\delta \mathcal {U}_{\perp }^{\alpha }$, which amounts to solving $M_{11} = 1$. This yields

(4.6)\begin{equation} \frac{b^2-2 \sqrt{b^2+1}+2 }{b^3} = \frac{2v_{s}}{B \kappa n}. \end{equation}

Assuming that $b \neq 0$, (4.6) can be transformed into a polynomial equation that can be solved analytically. Inserting the obtained roots into (4.6), we found that only two solutions to (4.6) exist, and only within the region $2\sqrt {2}v_{s} \le B \kappa n$. In analogy with (3.6) for the bulk suspension, the eigenvalue $\chi$ within this parameter range is found to have the following form:

(4.7)\begin{equation} \chi ={-} \lambda + G\left(\frac{2 v_{s}}{B\kappa n}\right) \frac{B \kappa n \mathcal{k}}{2}, \end{equation}

where

(4.8)\begin{gather} G(x) = 6x^{2} \left( 4 - \frac{\left(1\pm {\rm i} \sqrt{3}\right)}{ H(x)} - \left(1\mp {\rm i} \sqrt{3}\right) H(x) \right)^{{-}1}, \end{gather}

(4.9)\begin{gather}H(x) = \left[ 54 x^2 -1 + 6 \sqrt{3}x \sqrt{27 x^2 - 1}\right]^{1/3}. \end{gather}

Similar to the bulk suspension, in the region $2\sqrt {2} v_{s} > B \kappa n$, where the solution to (4.6) does not exist, stability analysis requires the use of an alternative approach. This is discussed in Appendix C.2, where we show that the system remains stable in this parameter range.

The function $G$ has the same qualitative features as the function $F$ in the case of bulk suspensions (compare figures 2 and 5 (see Appendix B)): its real part is positive and its largest branch is a monotonically decreasing function of its argument. Because (4.7) is an increasing function of $\mathcal {k}$, the orientational instability in a 2-D layer of microswimmers embedded in a 3-D fluid sets in at the largest possible value of $\mathcal {k}$, in contrast with bulk pusher suspensions. In the absence of any physical restrictions on the maximum value taken by $\mathcal {k}$, this instability occurs for

(4.10)\begin{equation} n > n_{c} = 2\sqrt{2} \frac{v_{s}}{B \kappa}. \end{equation}

Just as in bulk systems, the fact that $n_{c}$ diverges for spherical particles ($B=0$) shows that this is an instability driven by particle reorientations due to hydrodynamic interactions. However, unlike the 3-D bulk case (3.8), $n_c$ is now a function of the swimming speed $v_s$ rather than of the tumbling rate $\lambda$, meaning that increasing the swimming speed will stabilise a sheet of pushers, which become unstable for any density in the shaker limit $v_s \to 0$. While the stabilising role of swimming in the stability and pre-transitional correlations of finite microswimmer systems has previously been discussed for 3-D bulk suspensions (Saintillan & Shelley 2013; Liu, Zhang & Cheng 2019; Škultéty et al. 2020; Martinez et al. 2020; Albritton & Ohm 2023), its effect on the stability of a sheet of pusher microswimmers, (4.10), is significantly stronger. Finally, the independence of $n_{c}$ of the tumbling rate implies that even a suspension of straight swimmers has a non-vanishing critical density, in contrast to the case encountered in 3-D bulk suspensions.

The above reasoning for deriving the orientational instability assumes that the most unstable wave vector is $\mathcal {k}_c \to \infty$; this is clearly unphysical since it corresponds to microscopic length scales. In practice, the instability instead sets in at some finite length scale $l_c$, corresponding to a finite value of $\mathcal {k}_c$. As observed previously for 3-D bulk suspensions (Saintillan & Shelley 2008a; Subramanian & Koch 2009), one possible source of regularisation is provided by the spatial diffusivity of microswimmers. Repeating the above analysis in the presence of Brownian diffusion yields an additional term, $-D\mathcal {k}^2$, on the right-hand side of (4.7), where $D$ is the diffusion constant. Minimising the real part of the eigenvalue with respect to $k$ leads to the length scale $l_c \sim \sqrt {D/\lambda }$ being selected at the instability. Using the approximate values $D\sim 0.2- 0.4\ \mathrm {\mu }\textrm {m}^2\ \textrm {s}^{-1}$ (Jepson et al. 2013) and $\lambda \sim 1\ \textrm {s}^{-1}$, this length scale becomes comparable to the microswimmer size, $l_c \sim 1\ \mathrm {\mu }\textrm {m}$, which is an unphysically small length scale for the instability to occur. It is therefore unlikely that spatial diffusivity is the relevant mechanism of the length scale selection in (4.7).

The second relevant microscopic length scale in the problem is the (density-dependent) average distance between microswimmers, below which the system can no longer be viewed as a continuum. Thus, we assume that the length scale selected at the instability is now $l_c \sim 2/\sqrt {{\rm \pi} n_c}$, leading to a maximum wave vector

(4.11)\begin{equation} \mathcal{k}_c = \frac{2{\rm \pi}}{l_c} = \sqrt{ {\rm \pi}^{3} n_c }. \end{equation}

The calculation of the corresponding instability threshold is outlined in Appendix D, leading to the following approximate expression for $n_c$:

(4.12)\begin{equation} n_{c} \approx \frac{4}{\rm \pi} \left( \frac{\lambda}{B\kappa} \right)^{2/3} + 2\sqrt{2} \frac{v_{s}}{B\kappa}. \end{equation}

We observe that the critical density (4.12) differs from its $\mathcal {k}_c \to \infty$ counterpart (4.10) by the presence of the first term, which dominates at low swimming speeds. Using the experimental values $v_s \approx 15\ \mathrm {\mu }\textrm {m}\ \textrm {s}^{-1}$, $\kappa \approx 800\ \mathrm {\mu }\textrm {m}^3\ \textrm {s}^{-1}$ measured in Drescher et al. (2011), together with $B \approx 1$ and $\lambda \approx 1\ \textrm {s}^{-1}$, (4.12) and (4.11) give the critical length scale $l_c\sim 4$ $\mathrm {\mu }$m, which is larger than, although comparable to, the value resulting from translational diffusion.

4.2. Density instability

The second subset of (4.4) couples $\delta \hat {\mathcal {U}}_{\parallel }^{\alpha }$ and $\delta \hat {\rho }$ and thus governs the appearance of density modulations within the microswimmer layer. After some algebra, this set of equations can be re-expressed as

(4.13)\begin{equation} \frac{\varPhi \tilde{\mathcal{k}} \left[ B \left(b^2 + 2- 2 \sqrt{b^2+1}\right) y - b^2 \left(\sqrt{b^2+1}-1\right) (y+1) \right]}{4 b^4 (y+1) \left[ \sqrt{b^2+1}(y + 1)-1 \right]} = 1, \end{equation}

where we have introduced the following dimensionless quantities:

(4.14a–d)\begin{equation} b = \frac{\tilde{\mathcal{k}}}{y + 1}, \quad y = \frac{\chi}{\lambda}, \quad \varPhi = \frac{\kappa}{v_{s}}n, \quad \tilde{\mathcal{k}} = \frac{v_{s}}{\lambda} \mathcal{k}. \end{equation}

A full analytical solution of (4.13) cannot be found, but can be achieved perturbatively for small $\tilde {\mathcal {k}}$, yielding

(4.15)\begin{equation} y ={-} \frac{1}{8} \varPhi \tilde{\mathcal{k}} - \left( \frac{1}{2} + \frac{B \varPhi^{2}}{128} \right) \tilde{\mathcal{k}}^{2} + \mathcal{O}(\tilde{\mathcal{k}}^{3}). \end{equation}

From this expression, we conclude that (i) pusher suspensions ($\varPhi >0$) are stable against density modulations at large spatial scales since $y < 0$ for small $\tilde {\mathcal {k}}$, and (ii) puller suspensions ($\varPhi < 0$) are linearly unstable at any density. Note that this instability is independent of $B$, and hence occurs even for spherical particles. Even though the puller density instability occurs at vanishingly small densities in an infinite system, for a finite system with linear size $H$ one can show that puller suspensions are unstable against large-scale perturbations only above a critical density $n_c$ given by

(4.16)\begin{equation} n_c ={-}\frac{8{\rm \pi}}{H} \frac{v_{s}^{2}}{\lambda\kappa}. \end{equation}

This result is obtained by assuming that the eigenvalue $y$ is real for small $\tilde {\mathcal {k}}$, setting $y = 0$ in (4.13), and assuming that the critical wave vector $\mathcal {k}_c$ is set by the system dimensions, i.e. $\mathcal {k}_c = 2{\rm \pi} /H$.

We now proceed by solving (4.13) numerically, which we have done for several values of $B$ and $\varPhi$, as shown in figure 3. This analysis indeed confirms both the existence of a finite-wavelength instability in puller suspensions at large spatial scales (small $\mathcal {k}$), and that a sheet of pushers is always stable against density perturbations. However, it also shows that, when $\varPhi$ is sufficiently negative, $\mathrm {Re}[y]$ increases in an unbounded fashion at large $\mathcal {k}$ for pullers. Thus, a second density instability, now at small spatial scales, emerges in puller suspensions for sufficiently large densities. To analyse this effect further, we numerically calculate the largest eigenvalue from (4.13). The resulting phase diagram is shown in figure 4, which displays the critical scale $\mathcal {k}_c$ as a function of the dimensionless system size $H \lambda /v_s$ and the reduced particle density $\varPhi$. Our findings confirm that the density instability in the confined suspension sets in either at the smallest or the largest available spatial scale, depending on the system parameters as encoded in $\varPhi$ and $B$. In drawing the phase diagram we have, in analogy to § 4.1, chosen a cutoff at small scales corresponding to the interparticle separation (i.e. $\mathcal {k} = \sqrt { {\rm \pi}^{3} n }$ ). In the dimensionless units used here, the maximal wavenumber available is set by $\tilde {\mathcal {k}}=\sqrt {{\rm \pi} ^3 \varPhi } \sqrt {v_s^3/\lambda ^2 \kappa }$, where the factor $\sqrt {v_s^3/\lambda ^2 \kappa } \sim \mathcal {O}(1)$ for free-swimming E. coli bacteria (Drescher et al. 2011).

Figure 3. The real (a) and imaginary (b) parts of the eigenvalue $y = \chi /\lambda$ corresponding to density perturbations, obtained by numerically solving (4.13). All plots correspond to needle-like particles ($B = 1$), but are qualitatively similar for all values of $B$. For pushers ($\varPhi > 0$), the real part of the eigenvalue is strictly negative. For pullers ($\varPhi < 0$), the real part of the eigenvalue becomes positive at small wave vectors. At larger densities, the global maximum of $\mathrm {Re} [y]$ moves from small $\tilde {\mathcal {k}}$ to $\tilde {\mathcal {k}} \to \infty$. No solution exists in the region Re$[y] < - 1$, whose stability is instead addressed in Appendix C.1.

Figure 4. Density instability in a 2-D puller suspension embedded in a 3-D fluid. The phase diagrams show regions of the two types of density instability obtained by numerical solution of (4.13). At high $\varPhi$ (grey regions), the instability sets in at the smallest physically relevant spatial scale, which we assume to be the particle–particle separation, i.e. $\tilde {\mathcal {k}}=\sqrt {{\rm \pi} ^3 \varPhi } \sqrt {v_s^3/\lambda ^2 \kappa }$, and set $\sqrt {v_s^3/\lambda ^2 \kappa }=1$. At moderate and low $\varPhi$ (red heat map), the instability sets in at larger spatial scales set by the maximum of the arc-like part of the dispersion law; see figure 3. If the system size $H$ is too small, the latter maximum in the dispersion law cannot be accessed, and the instability instead sets in at the scale of the system dimensions. (a) $B = 1$, (b) $B = 0$.

These results show that, while the density instability is present also for spherical particles, the particle shape strongly affects the emergence of small- and large-scale instabilities in that the small-scale instability dominates over a wider range of $\varPhi$ values compared with the elongated ($B=1$) particle case. Understanding the physical manifestations of these two instabilities requires detailed investigations using explicit numerical simulations and is left for future work.

5. Conclusion

In this work, we have analysed the linear stability of the homogeneous and isotropic distribution of microswimmers restricted to a 2-D plane embedded in a 3-D bulk viscous fluid. Our method, based on mean-field kinetic theory, is somewhat different from that used in the previous analysis of bulk suspensions (Saintillan & Shelley 2008a,b; Subramanian & Koch 2009; Hohenegger & Shelley 2010) and allows us to easily identify the origin of the instability as either an orientational mode, stemming from the mutual reorientation of microswimmers, or a density mode, stemming from local accumulation of microswimmers.

The instabilities we found for a sheet of microswimmers are qualitatively different from their bulk counterpart. In such infinite 2-D or 3-D systems, linear stability analysis shows the presence of an orientational instability for pushers above the onset density given respectively by (B5) (2-D) and (3.8) (3-D). This instability sets in at the largest scale available in the system, typically the dimension of the experimental set-up or of the simulation box. Furthermore, due to the incompressible nature of the embedding fluid, dilute 2-D and 3-D bulk suspensions of microswimmers are always stable towards fluctuations in the microswimmer density.

In contrast, we showed that a 2-D layer of microswimmers embedded in a 3-D fluid is unstable for both pushers and pullers. For the former swimmer class, an orientational instability similar to the pusher instability in a bulk fluid occurs for large enough densities. However, in contrast to the bulk case, the orientational instability in the 2-D layer sets in at the smallest spatial scale, which we interpret as the one below which our continuum description is not valid anymore, i.e. the average particle–particle separation. For parameter values roughly corresponding to E. coli bacteria, this leads to the instability setting in at a length scale of $\sim 5\ \mathrm {\mu }$m, comparable to the bacterial dimensions.

On the other hand, a layer of puller microswimmers is susceptible to a density instability characterised by spatial aggregation of microswimmers independent of their orientations. This instability is caused by the mutual advection of the microswimmers and is the direct consequence of the sink-like nature of their in-plane velocity fields, as discussed in § 1. In the thermodynamic limit, the puller suspension is unstable at any density, while a non-zero critical density emerges for finite system sizes, as shown in (4.16). The scale at which this puller instability sets in depends on the microswimmer density: at small and moderate densities, the instability occurs at long length scales, while at high enough densities it is replaced by a small-scale instability. While the cross-over between these two regimes depends on the particle shape through $B$, both instabilities are generically present regardless of particle shape. Whether this implies that the system reaches different steady states at long times is still an open question that we leave for future work. We also note that the density instability presented here is qualitatively similar to the instability reported by Baskaran & Marchetti (2009), as both are caused by the fluid's compressibility; in the case of Baskaran & Marchetti (2009), however, the compressibility arises as a result of an erroneous coarse-graining procedure, as discussed by Aranson (2022).

We would like to stress once more that the instability scenarios discussed here are not the consequence of the spatial dimensionality, as 2-D and 3-D bulk suspensions exhibit qualitatively a very similar behaviour. Instead, the qualitative differences to the corresponding bulk suspensions come from the fact that the in-plane fluid in which the microswimmers move is effectively compressible, an effect which is generically present in the case of motion restricted to a subset of the embedding space, including microswimmer layers in the vicinity of solid and liquid boundaries. Importantly, this effect is likely to be one of the main driving forces behind the boundary-induced clustering and phase separation observed in self-propelled suspension droplets (Thutupalli et al. 2018) and squirmer suspensions (Singh & Adhikari 2016). These authors highlight the importance of the presence of a boundary in inducing effective in-plane attractions among puller swimmers after accumulating at the surface or interface, although they rationalise their findings in terms of attractive hydrodynamic forces without explicitly mentioning the effective compressibility mechanism present in this geometry. Thus, even though our chosen geometry of a 2-D plane in an infinite 3-D fluid seems artificial, we argue that the generic mechanism demonstrated here is in fact an important driving force behind accumulation in active-particle suspensions near interfaces.

Furthermore, the phenomenon of effective in-plane compressibility is not limited to dipolar flow fields. As demonstrated by Spagnolie & Lauga (2012), all flow singularities relevant to self-propulsion of microorganisms can be obtained by repeatedly applying the $(\boldsymbol {p} \boldsymbol {\cdot } \boldsymbol {\nabla })$ operator to the stokeslet and the point source singularity. Therefore, these singularities will have alternating parities with respect to the angle defined by $\boldsymbol {x} \boldsymbol {\cdot } \boldsymbol {p}$. In the context of the calculation presented in § 1, this leads to a layer of force dipoles, force octupoles and so on having effectively compressible in-plane velocity fields, while the remaining flow singularities remain incompressible even when restricted to a lower-dimensional subset of the embedding space. The same reasoning holds for the source singularities, where sources, source quadrupoles and so on exhibit effectively compressible in-plane velocity fields. Our findings are therefore likely to be of more general importance for hydrodynamic interactions and collective motion of realistic microorganisms near boundaries.

One of the remaining open questions pertains to the stability of a sheet of microswimmers with respect to out-of-plane perturbations. At first glance, a sheet of pushers generates on average a fluid flow towards the sheet along the third direction and one thus expects out-of-plane perturbations to be suppressed by such a flow. By the same argument, the opposite is expected for pullers. Some support for this hypothesis can be drawn from Ishikawa & Pedley (2008), who simulated a collection of spherical squirmers confined to a single layer, as in our set-up. They observed no instability for pushers, while pullers exhibited formation of dense clusters visually consistent with our density instability. Interestingly, when bottom-heavy puller squirmers were allowed to move out of the original monolayer, the authors observed the formation of a dynamical steady state in a form of a microswimmer band, suggesting that pullers are indeed unstable with respect to out-of-plane perturbations. In contrast, the continuum simulations by Nejad & Yeomans (2022) showed coexisting regions of in-plane and out-of-plane director orientations in extensile active nematics, while no such instability was present in the corresponding contractile systems. This complex picture clearly has to be understood via a combination of linear stability analysis and fully nonlinear simulations, a question that we leave open to future studies. We stress, however, that we view the set-up considered in this work as a way to isolate the effect of in-plane compressibility on collective motion of experimental realisations of microswimmers accumulated close to a physical boundary. In such systems, there are always additional effects ensuring microswimmers’ presence at the boundary, e.g. gravity or surface tension, and the question of stability with respect to out-of-plane perturbations should be addressed in that context.

We conclude by observing that, while the linear stability analysis presented here identifies the origin, the onset conditions and the associated scales of the instabilities, their nonlinear evolution and ensuing non-equilibrium steady states must be assessed in numerical simulations of particle-based models (Bárdfalvy et al. 2019, 2020) or mean-field moment equations (Saintillan & Shelley 2013). We aim to further explore this and similar questions in future studies.