# Hydrodynamic instabilities in a 2-D sheet of microswimmers embedded in a 3-D
fluid

Viktor Škultéty1  Dóra Bárdfalvy2  Joakim Stenhammar2
Cesare Nardini3  Alexander Morozov1<EMAIL_ADDRESS>1 SUPA, School
of Physics and Astronomy, The University of Edinburgh, James Clerk Maxwell
Building, Peter Guthrie Tait Road, Edinburgh, EH9 3FD, United Kingdom
2Division of Physical Chemistry, Lund University, Box 124, S-221 00 Lund,
Sweden 3Servicede Physique de l’État Condensé, CNRS UMR 3680, CEA-Saclay,
91191 Gif-sur-Yvette, France

###### 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 2-D and 3-D, 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 3-D bulk fluid, and that microswimmers on average act
as sources (pushers) or sinks (pullers). We analyse 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.

###### keywords:

Authors should not enter keywords on the manuscript, as these must be chosen
by the author during the online submission process and will then be added
during the typesetting process (see
http://journals.cambridge.org/data/relatedlink/jfm-keywords.pdf for the full
list)

##  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., 2013a) and
significant changes of the apparent suspension viscosity (Rafaï et al., 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 describe
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 et al., 2010). While these two flow fields
lead to equivalent statistical properties for very dilute microswimmer
suspensions, where swimmers can be treated as effectively noninteracting, 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 suspensions of pusher bacteria, these theories predict the onset
of large-scale flows above a well-defined critical microswimmer volume
fraction, and identifies mutual particle reorientation due to hydrodynamic
interactions as the main mechanism behind the instability, while no such
instability occurs in three-dimensional 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 to puller
suspensions at the same density (Stenhammar et al., 2017; Škultéty et al.,
2020).

In practice, observations of three-dimensional bulk collective flows are
complicated by the unavoidable presence of boundaries in experimental setups
used to study motile microorganisms. Importantly, motile organisms accumulate
in close vicinity of 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 (Mathijssen et al., 2016; Jeanneret et al.,
2019; Brotto et al., 2013). 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 et al., 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.

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

An additional important physical effect inherent to confined systems, though
rarely discussed explicitly, is the effective compressibility of the in-plane
flow fields generated by a layer of microswimmers (Salbreux et al., 2009;
Jülicher et al., 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, 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 ${\bm{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
${\bm{u}}_{\parallel}(z=h)=0$, and is given by Spagnolie & Lauga (2012):

$\displaystyle{\bm{u}}_{\parallel}({\bm{\mathcal{x}}})=$ $\displaystyle\
\frac{\kappa}{8\pi}\Bigg{(}\frac{{\bm{\mathcal{x}}}}{|{\bm{\mathcal{x}}}|^{3}}\left[3\frac{({\bm{\mathcal{x}}}\cdot{\bm{\mathcal{p}}})^{2}}{|{\bm{\mathcal{x}}}|^{2}}-1\right]+\frac{{\bm{\mathcal{x}}}}{R^{3}}$
$\displaystyle\hskip
56.9055pt-\frac{3{\bm{\mathcal{x}}}({\bm{\mathcal{x}}}\cdot{\bm{\mathcal{p}}})^{2}+6h^{2}\left\\{{\bm{\mathcal{x}}}+2{\bm{\mathcal{p}}}({\bm{\mathcal{x}}}\cdot{\bm{\mathcal{p}}})\right\\}}{R^{5}}+\frac{30h^{2}({\bm{\mathcal{x}}}\cdot{\bm{\mathcal{p}}})^{2}{\bm{\mathcal{x}}}}{R^{7}}\Bigg{)}.$
(1)

Here, ${\bm{\mathcal{p}}}$ and ${\bm{\mathcal{x}}}$ are two-dimensional
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{|{\bm{\mathcal{x}}}|^{2}+4h^{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
Fig. 1, yielding

$\displaystyle\int_{|{\bm{\mathcal{x}}}|=X}\mathrm{d}{\bm{\mathcal{x}}}\cdot{\bm{u}}_{\parallel}=\frac{\kappa}{8X}\Biggl{[}1-\frac{X^{4}-10X^{2}h^{2}+64h^{4}}{\left(X^{2}+4h^{2}\right)^{7/2}}X^{3}\Biggr{]}.$
(2)

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 than their bulk counterparts, which has
previously been shown to promote dynamic self-assembly and crystallisation in
active particle systems (Thutupalli et al., 2018; Singh & Adhikari, 2016). 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.

The aim of this work is to understand the effect of this in-plane
compressiblity on collective motion in a minimal setting. We consider a layer
of microswimmers restricted to move in an infinitely thin, two-dimensional
plane embedded in a three-dimensional 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 three-dimensional force dipoles. As can be seen from Eq.
(2), even in this limit, the in-plane velocity field generated by the
microswimmer is still effectively compressible. Therefore, although this
specific setup 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 two-dimensional layer of microswimmers embedded in
a three-dimensional 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 Section 2, we present a general kinetic
theory describing a dilute suspension of microswimmers interacting through
long-range velocity fields. In Section 3 we review its predictions for the the
widely studied case of microswimmers suspended in a three-dimensional bulk
fluid. In Section 4 we report on linear stability of a two-dimensional layer
of microswimmers embedded in a three-dimensional fluid, which constitutes the
main novel results of our work. We conclude by summarising our main findings
in Section 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 setup.

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

$\displaystyle\int\mathrm{d}{\bm{x}}\,\mathrm{d}{\bm{p}}\
\Psi({\bm{x}},{\bm{p}},t)=1.$ (3)

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

$\displaystyle\partial_{t}\Psi+\nabla^{\alpha}\\{\dot{x}^{\alpha}\Psi\\}+\partial^{\alpha}\\{\dot{p}^{\alpha}\Psi\\}=-\lambda\Psi+\lambda\int\frac{\mathrm{d}{\bm{p}}}{\Omega_{d}}\
\Psi,$ (4)

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;
$\Omega_{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:

$\displaystyle\dot{x}^{\alpha}$
$\displaystyle=v_{s}p^{\alpha}+\mathcal{U}^{\alpha}({\bm{x}}),$ (5)
$\displaystyle\dot{p}^{\alpha}$
$\displaystyle=(\delta^{\alpha\beta}-p^{\alpha}p^{\beta})({\mathcal{W}}^{\beta\gamma}({\bm{x}})+B{\mathcal{E}}^{\beta\gamma}({\bm{x}}))p^{\gamma},$
(6)

where dots denote time derivatives. According to Eq.(5), each particle changes
its spatial position due to its own swimming with a constant speed $v_{s}$ in
the direction of its orientation ${\bm{p}}$ and due to advection by the fluid
flow with velocity $\bm{\mathcal{U}}$ at its position ${\bm{x}}$, created by
all other microswimmers. Microswimmer orientations change according to
Jeffery’s equation, (6), which describes particle rotation by the gradients of
the fluid velocity at its position (Kim & Karrila, 2005). Here,
${\mathcal{W}}^{\alpha\beta}({\bm{x}})=\tfrac{1}{2}(\nabla^{\beta}\mathcal{U}^{\alpha}({\bm{x}})-\nabla^{\alpha}\mathcal{U}^{\beta}({\bm{x}})),$
and
${\mathcal{E}}^{\alpha\beta}({\bm{x}})=\tfrac{1}{2}(\nabla^{\beta}\mathcal{U}^{\alpha}({\bm{x}})+\nabla^{\alpha}\mathcal{U}^{\beta}({\bm{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 Eqs. (5) and (6) 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 (4).

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

$\displaystyle\mathcal{U}^{\alpha}({\bm{x}},t)=N\int\mathrm{d}{\bm{x}}^{\prime}\mathrm{d}{\bm{p}}^{\prime}\
u^{\alpha}({\bm{x}}-{\bm{x}}^{\prime},{\bm{p}}^{\prime})\Psi({\bm{x}}^{\prime},{\bm{p}}^{\prime},t).$
(7)

Here, ${\bm{u}}({\bm{x}}-{\bm{x}}^{\prime},{\bm{p}}^{\prime})$ is the
microscopic velocity field created at ${\bm{x}}$ by a microswimmer located at
${\bm{x}}^{\prime}$ with the orientation ${\bm{p}}^{\prime}$. In the
following, we approximate ${\bm{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 three-dimensional 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 Section 5.

###  2.1 Linear stability of the homogeneous and isotropic state

The Smoluchowski equation (4) and the normalisation condition (3) admit as a
solution the homogeneous and isotropic state
$\Psi_{\text{HI}}=1/(\Omega_{d}V_{d})$, where $V_{d}$ is the $d$-dimensional
volume of the suspension. The existence of this solution relies on the
condition

$\displaystyle\int\mathrm{d}{\bm{x}}^{\prime}\mathrm{d}{\bm{p}}^{\prime}\
u^{\alpha}({\bm{x}}-{\bm{x}}^{\prime},{\bm{p}}^{\prime})=0,$ (8)

which is true in all cases considered in this work. Notice that, while
integrals over ${\bm{x}}^{\prime}$ and ${\bm{p}}^{\prime}$ vanish
independently in two- and three-dimensional bulk suspensions, only the
positional integral vanishes for microswimmers restricted to a two-dimensional
plane in a three-dimensional fluid.

Stability of the homogeneous and isotropic state is determined by time
evolution of small perturbations around $\Psi_{\text{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,
three-dimensional 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\Psi({\bm{x}},{\bm{p}},t)$ of the
one-particle distribution function around the homogeneous and isotropic state,
_i.e._ ,

$\displaystyle\Psi({\bm{x}},{\bm{p}},t)=\frac{1}{\Omega_{d}V_{d}}+\delta\Psi({\bm{x}},{\bm{p}},t),$
(9)

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

$\displaystyle\hat{f}({\bm{k}})$ $\displaystyle=\int\mathrm{d}{\bm{x}}\
f({\bm{x}})\mathrm{e}^{-i{\bm{k}}\cdot{\bm{x}}},\quad
f({\bm{x}})=\frac{1}{(2\pi)^{d}}\int\mathrm{d}{\bm{k}}\
\hat{f}({\bm{k}})\mathrm{e}^{i{\bm{k}}\cdot{\bm{x}}},$ (10)

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

$\displaystyle\Bigl{[}\partial_{t}+\lambda+iv_{s}{\bm{p}}\cdot{\bm{k}}\Bigr{]}\delta\hat{\Psi}$
$\displaystyle=\frac{\lambda}{\Omega_{d}}\delta\hat{\rho}+\frac{n}{\Omega_{d}}\Bigl{[}dBp^{\alpha}p^{\beta}-(1+B)\delta^{\alpha\beta}\Bigr{]}ik^{\alpha}\delta\hat{\mathcal{U}}^{\beta},$
(11)

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

$\displaystyle\delta\hat{\rho}({\bm{k}},t)$
$\displaystyle=\int\mathrm{d}{\bm{p}}\ \delta\hat{\Psi}({\bm{k}},{\bm{p}},t),$
(12) $\displaystyle\delta\hat{\mathcal{U}}^{\alpha}({\bm{k}},t)$
$\displaystyle=\int\mathrm{d}{\bm{p}}\
\hat{u}^{\alpha}({\bm{k}},{\bm{p}})\delta\hat{\Psi}({\bm{k}},{\bm{p}},t).$
(13)

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

$\displaystyle\delta\hat{\Psi}({\bm{k}},{\bm{p}},t)=\delta\hat{\Psi}({\bm{k}},{\bm{p}})\mathrm{e}^{\chi
t},$ (14)

where the sign of the real part of the temporal eigenvalue $\chi$ determines
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; Škultéty et al., 2020). Working with the Ansatz (14) is
mathematically simpler but may not necessarily produce a dispersion law for
all parameter values, as we discuss further below.

Inserting (14) into (11), we derive the following closed set of equations for
the density and velocity perturbations

$\displaystyle\delta\hat{\rho}$
$\displaystyle=in\delta\hat{\mathcal{U}}^{\alpha}\int\frac{\mathrm{d}{\bm{p}}}{\Omega_{d}}\
\frac{Bdp^{\alpha}{\bm{p}}\cdot{\bm{k}}-(1+B)k^{\alpha}}{\mathcal{L}(\chi,{\bm{p}}\cdot{\bm{k}})}+\delta\hat{\rho}\int\frac{\mathrm{d}{\bm{p}}}{\Omega_{d}}\
\frac{\lambda}{\mathcal{L}(\chi,{\bm{p}}\cdot{\bm{k}})},$ (15)
$\displaystyle\delta\hat{\mathcal{U}}^{\alpha}$
$\displaystyle=in\delta\hat{\mathcal{U}}^{\beta}\int\frac{\mathrm{d}{\bm{p}}}{\Omega_{d}}\
\frac{\hat{u}^{\alpha}({\bm{k}},{\bm{p}})[Bdp^{\beta}{\bm{p}}\cdot{\bm{k}}-(1+B)k^{\beta}]}{\mathcal{L}(\chi,{\bm{p}}\cdot{\bm{k}})}+\delta\hat{\rho}\int\frac{\mathrm{d}{\bm{p}}}{\Omega_{d}}\
\frac{\lambda\hat{u}^{\alpha}({\bm{k}},{\bm{p}})}{\mathcal{L}(\chi,{\bm{p}}\cdot{\bm{k}})},$
(16)

where
$\mathcal{L}(\chi,{\bm{p}}\cdot{\bm{k}})=\chi+\lambda+iv_{s}{\bm{p}}\cdot{\bm{k}}$.
The solution to this eigenvalue problem depends on the precise form of the
microscopic velocity field $\hat{{\bm{u}}}({\bm{k}},{\bm{p}})$, which, in
turn, depends on the dimensionality of space and the dimensionality of the
vectors ${\bm{p}}$ and ${\bm{k}}$. In the following, we begin by revisiting
the known results on linear stability of bulk three-dimensional suspensions
(Saintillan & Shelley (2008a, b)) before turning our attention to a two-
dimensional layer of microswimmers embedded in a three-dimensional 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)

$\displaystyle
u_{3d}^{\alpha}({\bm{x}},{\bm{p}})=\frac{\kappa}{8\pi}\frac{x^{\alpha}}{|{\bm{x}}|^{3}}\left[3\frac{({\bm{x}}\cdot{\bm{p}})^{2}}{|{\bm{x}}|^{2}}-1\right],$
(17)

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

$\displaystyle\hat{u}_{3d}^{\alpha}({\bm{k}},{\bm{p}})=-i\kappa\frac{{\bm{k}}\cdot{\bm{p}}}{k^{2}}\mathbb{P}^{\alpha\beta}p^{\beta}.$
(18)

Here, $\kappa$ is the dipolar strength and
$\mathbb{P}^{\alpha\beta}=\delta^{\alpha\beta}-k^{\alpha}k^{\beta}/k^{2}$ is
the transversal projection operator that ensures that Eq.(18) satisfies the
incompressibility condition $k^{\alpha}\hat{u}_{3d}^{\alpha}=0$. Setting $d=3$
in (15)–(16), 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 ${\bm{p}}$ (see
further discussion below), we obtain the following set of equations (Martinez
et al., 2020)

$\displaystyle\delta\hat{\rho}$
$\displaystyle=\frac{\lambda}{v_{s}k}\arctan(b)\delta\hat{\rho},$ (19)
$\displaystyle\delta\hat{\mathcal{U}}^{\alpha}$
$\displaystyle=\frac{Bn\kappa}{v_{s}k}\frac{b(3+2b^{2})-3(b^{2}+1)\arctan(b)}{2b^{4}}\delta\hat{\mathcal{U}}^{\alpha},$
(20)

where we have introduced the dimensionless parameter

$\displaystyle b=\frac{v_{s}k}{\chi+\lambda}.$ (21)

We note that Eqs. (19) and (20) 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; Hohenegger & Shelley,
2010; Subramanian & Koch, 2009; 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

Figure 2: The real (left) and imaginary (right) part of the function
$F(\gamma)$ from (23).

Eq.(20) yields the following eigenvalue problem:

$\displaystyle\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}.$
(22)

We now solve Eq. (21) for $\chi$ and use the definition of $\gamma$ from Eq.
(22), to arrive at

$\displaystyle\chi=-\lambda+\frac{1}{5}F(\gamma)B\kappa n,\quad
F(\gamma)=\frac{\gamma}{b(\gamma)}.$ (23)

Eq. (22) 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 Fig. 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 Eq. (23), 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 Eq. (23), and from imposing $\chi>0$, that pusher suspensions are
unstable for densities larger than

$n_{c}=\frac{5\lambda}{B\kappa}.$ (24)

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

The origin of this instability is the mutual reorientation of microswimmers
(Saintillan & Shelley, 2008a, b; Hohenegger & Shelley, 2010; Subramanian &
Koch, 2009; Stenhammar et al., 2017). It is thus unsurprising to observe that
no instability is present for spherical particles ($B=0$). Furthermore, the
instability condition (24) 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 (11) obeying the Ansatz (14) ceases to exist for
$\gamma>\gamma^{*}$. The reason for this is that, in the process of deriving
(19)–(20), we have assumed that $\chi$ is independent of the swimmer
orientation ${\bm{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 (23),
leading to the instability threshold (24). The same conclusion is reached by
solving (11) using Laplace transform techniques as was done in Stenhammar et
al. (2017); Š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\pi\lambda/(2v_{s})$ this can be readily shown from Eq.
(19), yielding

$\displaystyle\chi=-\lambda+\frac{v_{s}k}{\tan(v_{s}k/\lambda)}.$ (25)

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, Eq. (19) does not allow one to
determine stability of the suspension when $k>\pi\lambda/(2v_{s})$, since in
this regime Eq. (19) has no solution. As discussed in Section 3.1, this is a
limitation of the Ansatz (14). In Appendix C.1, we analyse this problem in
some detail by solving Eq. (11) 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
Section 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)

$\displaystyle
u_{\text{plane}}^{\alpha}({\bm{\mathcal{k}}},{\bm{\mathcal{p}}})=-\frac{i\kappa}{2}\frac{{\bm{\mathcal{k}}}\cdot{\bm{\mathcal{p}}}}{\mathcal{k}}\left[\mathbb{P}^{\alpha\beta}+\frac{1}{2}\mathbb{Q}^{\alpha\beta}\right]\mathcal{p}^{\beta},$
(26)

where $\kappa$ is the dipolar strength which has the same dimensions as in the
bulk 3-D case, and
${\bm{\mathcal{k}}}=\left(\mathcal{k}_{x},\mathcal{k}_{y}\right)$,
${\bm{\mathcal{p}}}=\left(\mathcal{p}_{x},\mathcal{p}_{y}\right)$ are
respectively the in-plane wave and orientation vectors. In writing Eq. (26),
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 ${\bm{u}}_{\text{plane}}$. We now define
the transversal (incompressible) and longitudinal (compressible) velocity
perturbations by

$\displaystyle\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},$
(27)

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

$\displaystyle\begin{pmatrix}M_{11}-1&0&0\\\
0&M_{22}-1&i\mathcal{k}^{\alpha}M_{23}\\\
0&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,$ (28)

where

$\displaystyle M_{11}$ $\displaystyle=\frac{\kappa
n}{v_{s}}\frac{B\left(b^{2}-2\sqrt{b^{2}+1}+2\right)}{2b^{3}},\quad
M_{22}=\frac{\kappa
n}{v_{s}}\frac{Bb^{2}+(b^{2}+2B)\left(1-\sqrt{b^{2}+1}\right)}{4b^{3}\sqrt{b^{2}+1}},$
(29) $\displaystyle M_{23}$
$\displaystyle=\frac{\kappa\lambda}{v_{s}\mathcal{k}^{2}}\frac{1}{4b}\left(\frac{1}{\sqrt{b^{2}+1}}-1\right),\quad
M_{32}=\frac{n}{v_{s}\mathcal{k}}\frac{2B\left(\sqrt{b^{2}+1}-1\right)-b^{2}(B+1)}{b\sqrt{b^{2}+1}},$
$\displaystyle M_{33}$
$\displaystyle=\frac{\lambda}{\mathcal{k}v_{s}}\frac{b}{\sqrt{1+b^{2}}},$

and

$\displaystyle b=\frac{v_{s}\mathcal{k}}{\chi+\lambda}.$ (30)

As is apparent from Eq. (28), $\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 transversal mode
$\delta\mathcal{U}_{\perp}^{\alpha}$, which amounts to solving $M_{11}=1$.
This yields

$\displaystyle\frac{b^{2}-2\sqrt{b^{2}+1}+2}{b^{3}}=\frac{2v_{s}}{B\kappa n}.$
(31)

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

$\displaystyle\chi=-\lambda+G\left(\frac{2v_{s}}{B\kappa
n}\right)\frac{B\kappa n\mathcal{k}}{2},$ (32)

where

$\displaystyle G(x)$ $\displaystyle=6x^{2}\left(4-\frac{\left(1\pm
i\sqrt{3}\right)}{H(x)}-\left(1\mp i\sqrt{3}\right)H(x)\right)^{-1},$ (33)
$\displaystyle H(x)$
$\displaystyle=\left[54x^{2}-1+6\sqrt{3}x\sqrt{27x^{2}-1}\right]^{1/3}.$ (34)

Similar to the bulk suspension, in the region $2\sqrt{2}v_{s}>B\kappa n$,
where the solution to (31) 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 Figs. 2 and 5 (see Appendix B)): its real
part is positive and its largest branch is a monotonically decreasing function
of its argument. Because Eq. (32) is an increasing function of $\mathcal{k}$,
the orientational instability in a two-dimensional layer of microswimmers
embedded in a three-dimensional 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

$\displaystyle n>n_{c}=2\sqrt{2}\frac{v_{s}}{B\kappa}.$ (35)

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 (24), $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 et al., 2019; Škultéty et al., 2020; Martinez et al., 2020;
Albritton & Ohm, 2023), its effect on the stability of a sheet of pusher
microswimmers, Eq.(35), 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 three-dimensional 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 three-
dimensional 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 (32), where $D$ is the diffusion constant. Minimising the
real part of the eigenvalue with respect to $k$ leads to the lengthscale
$l_{c}\sim\sqrt{D/\lambda}$ being selected at the instability. Using the
approximate values $D\sim 0.2-0.4$ $\mu$m2/s (Jepson et al., 2013b) and
$\lambda\sim 1$ s-1, this lengthscale becomes comparable with the microswimmer
size, $l_{c}\sim 1$ $\mu$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 (32).

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{\pi n_{c}}$, leading to
a maximum wave vector

$\displaystyle\mathcal{k}_{c}=\frac{2\pi}{l_{c}}=\sqrt{\pi^{3}n_{c}}.$ (36)

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

$\displaystyle
n_{c}\approx\frac{4}{\pi}\left(\frac{\lambda}{B\kappa}\right)^{2/3}+2\sqrt{2}\frac{v_{s}}{B\kappa}.$
(37)

We observe that the critical density (37) differs from its
$\mathcal{k}_{c}\to\infty$ counterpart Eq. (35) by the presence of the first
term, which dominates at low swimming speeds. Using the experimental values
$v_{s}\approx 15$ $\mu$m/s, $\kappa\approx 800$ $\mu$m3/s measured in Drescher
et al. (2011), together with $B\approx 1$ and $\lambda\approx 1$ s-1, Eqs.
(37) and (36) give the critical length-scale $l_{c}\sim 4$ $\mu$m, which is
larger than, although comparable to, the value resulting from translational
diffusion.

###  4.2 Density instability

Figure 3:  The real (left) and imaginary (right) part of the eigenvalue
$y=\chi/\lambda$ corresponding to density perturbations, obtained by
numerically solving Eq. (38). All plots correspond to needle-like particles
($B=1$), but are qualitatively similar for all values of $B$. For pushers
($\Phi>0$), the real part of the eigenvalue is strictly negative. For pullers
($\Phi<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.

The second subset of Eqs. (29) 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

$\displaystyle\frac{\Phi\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]}{4b^{4}(y+1)\left[\sqrt{b^{2}+1}(y+1)-1\right]}=1,$
(38)

where we have introduced the following dimensionless quantities:

$\displaystyle b$ $\displaystyle=\frac{\tilde{\mathcal{k}}}{y+1},\quad
y=\frac{\chi}{\lambda},\quad\Phi=\frac{\kappa}{v_{s}}n,\quad\tilde{\mathcal{k}}=\frac{v_{s}}{\lambda}\mathcal{k}.$
(39)

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

$\displaystyle y$
$\displaystyle=-\frac{1}{8}\Phi\tilde{\mathcal{k}}-\left(\frac{1}{2}+\frac{B\Phi^{2}}{128}\right)\tilde{\mathcal{k}}^{2}+\mathcal{O}(\tilde{\mathcal{k}}^{3}).$
(40)

From this expression, we conclude that $(i)$ pusher suspensions ($\Phi>0$) are
stable against density modulations at large spatial scales since $y<0$ for
small $\tilde{\mathcal{k}}$, and $(ii)$ puller suspensions ($\Phi<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

$\displaystyle n_{c}=-\frac{8\pi}{H}\frac{v_{s}^{2}}{\lambda\kappa}.$ (41)

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

We now proceed by solving Eq. (38) numerically, which we have done for several
values of $B$ and $\Phi$, as shown in Fig. 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
$\Phi$ 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 Eq. (38). The resulting phase diagram is
shown in Fig. 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 $\Phi$. 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
$\Phi$ and $B$. In drawing the phase diagram we have, in analogy to Section
4.1, chosen a cutoff at small scales corresponding to the interparticle
separation (_i.e._ , $\mathcal{k}=\sqrt{\pi^{3}n}$ ). In the dimensionless
units used here, the maximal wave number available is set by
$\tilde{\mathcal{k}}=\sqrt{\pi^{3}\Phi}\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).

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 $\Phi$ values compared to 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.

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 Eq. (38). At high $\Phi$ (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{\pi^{3}\Phi}\sqrt{v_{s}^{3}/\lambda^{2}\kappa}$,
and set $\sqrt{v_{s}^{3}/\lambda^{2}\kappa}=1$. At moderate and low $\Phi$
(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 Fig. 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.

##  5 Conclusion

In this work, we have analysed linear stability of the homogeneous and
isotropic distribution of microswimmers restricted to a two-dimensional plane
embedded in a three-dimensional 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 two- or three-
dimensional systems, linear stability analysis shows the presence of an
orientational instability for pushers above the onset density given
respectively by Eqs. (54) (two-dimensional) and (24) (three-dimensional). This
instability sets in at the largest scale available in the system, typically
the dimension of the experimental setup or of the simulation box. Furthermore,
due to the incompressible nature of the embedding fluid, dilute two- and
three-dimensional bulk suspensions of microswimmers are always stable towards
fluctuations in the microswimmer density.

In contrast, we showed that a two-dimensional layer of microswimmers embedded
in a three-dimensional 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 two-
dimensional 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\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 Section 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 Eq.
(41). 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 crossover 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 two- and three-
dimensional 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 two-dimensional plane in
an infinite three-dimensional 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 $({\bm{p}}\cdot\nabla)$ operator to the
Stokeslet and the point source singularity. Therefore, these singularities
will have alternating parities with respect to the angle defined by
${\bm{x}}\cdot{\bm{p}}$. In the context of the calculation presented in
Section 1, this leads to a layer of force dipoles, force octupoles, _etc._ ,
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, _etc._ , 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 setup. 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 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 non-
linear simulations, a question that we leave open to future studies. We
stress, however, that we view the setup 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 non-linear 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.

## Acknowledgement

For the purpose of open access, the authors have applied a Creative Commons
Attribution (CC BY) licence to any Author Accepted Manuscript version arising
from this submission.

## Funding

JS acknowledges financial support from the Swedish Research Council (Project
No. 2019-03718). AM acknowledges financial support from EPSRC (grant number
EP/V048198/1). CN acknowledges funding from Institut de physique du CNRS.

## Competing Interestst

The authors report no competing interests.

##  Appendix A Dipolar fluid velocity

The far-field fluid velocity ${\bm{u}}$ generated by a swimming microorganism
at low Reynolds number is described by the hydrodynamic force dipole, (Lauga &
Powers, 2009; Spagnolie & Lauga, 2012). The latter satisfies Stokes’ equation
driven by two point-like forces $\pm f{\bm{p}}$ acting on the fluid at the
positions $\pm\tfrac{1}{2}l{\bm{p}}$,

$\displaystyle\mu\nabla^{2}u^{\alpha}({\bm{x}})-\nabla^{\alpha}P({\bm{x}})+fp^{\alpha}[\delta({\bm{x}}-\tfrac{1}{2}l{\bm{p}})-\delta({\bm{x}}+\tfrac{1}{2}l{\bm{p}})]=0,$
(42)
$\displaystyle\qquad\qquad\qquad\qquad\nabla^{\alpha}u^{\alpha}({\bm{x}})=0,$
(43)

Here, $\mu$ is the viscosity of the fluid, $P$ is the pressure, and
$\delta({\bm{x}})$ denotes the Dirac delta-function. The solution to this
system is readily found by performing the $d$-dimensional Fourier transform,
(10), which gives

$\displaystyle-\mu
k^{2}\hat{u}^{\alpha}-ik^{\alpha}\hat{P}+fp^{\alpha}[\mathrm{e}^{-\frac{1}{2}il{\bm{k}}\cdot{\bm{p}}}-\mathrm{e}^{\frac{1}{2}il{\bm{k}}\cdot{\bm{p}}}]=0,$
(44) $\displaystyle\qquad\qquad\qquad\qquad k^{\alpha}\hat{u}^{\alpha}=0.$
(45)

The point-dipole approximation, relevant at large scales where $kl\ll 1$, is
obtained to linear order in the dipolar length $l$, yielding

$\displaystyle\hat{u}^{\alpha}({\bm{k}},{\bm{p}})=-i\kappa\frac{{\bm{k}}\cdot{\bm{p}}}{k^{2}}\mathbb{P}^{\alpha\beta}p^{\beta}.$
(46)

Here,
$\mathbb{P}^{\alpha\beta}=\delta^{\alpha\beta}-k^{\alpha}k^{\beta}/k^{2}$ is
the transversal projection operator, $k=|{\bm{k}}|$, and $\kappa=fl/\mu$ is
the dipolar strength.

The derivation presented above is independent of the spatial dimensionality
$d$ and the Fourier transform of the dipolar field (46) has the same form for
bulk systems of all dimensions. Its real space representation, however,
strongly depends on $d$. In three dimensions, the inverse Fourier transform
(10) yields

$\displaystyle
u_{3d}^{\alpha}({\bm{r}},{\bm{p}})=\frac{\kappa}{8\pi}\frac{r^{\alpha}}{|{\bm{r}}|^{3}}\left[3\frac{({\bm{r}}\cdot{\bm{p}})^{2}}{|{\bm{r}}|^{2}}-1\right],\quad{\bm{r}}={\bm{x}}-{\bm{x}}_{0},$
(47)

where we have introduced an arbitrary position ${\bm{x}}_{0}$ of the force
dipole. In two dimensions, the real-space representation is

$\displaystyle
u_{2d}^{\alpha}({\bm{r}},{\bm{p}})=\frac{\kappa}{4\pi}\frac{r^{\alpha}}{|{\bm{r}}|^{2}}\left[2\frac{({\bm{r}}\cdot{\bm{p}})^{2}}{|{\bm{r}}|^{2}}-1\right],\quad{\bm{r}}={\bm{x}}-{\bm{x}}_{0},$
(48)

where we notice that the dipolar strength $\kappa=fl/\mu$ has different
physical dimensions in two- and three-dimensional bulk systems.

The velocity field of a hydrodynamic dipole restricted to a two-dimensional
plane embedded in an infinite three-dimensional fluid can be obtained from
(47). We select $z=0$ to be the plane of the microswimmers, and recall the
definitions of two-dimensional vectors ${\bm{\mathcal{x}}}=(x,y)$, and
${\bm{\mathcal{p}}}=(p_{x},p_{y})$, from the main text. Performing the
corresponding two-dimensional Fourier transform of (47) with respect to
${\bm{\mathcal{x}}}$, leads to

$\displaystyle\hat{u}_{\text{plane}}^{\alpha}({\bm{\mathcal{k}}},{\bm{\mathcal{p}}})=-\frac{i\kappa}{2}\frac{{\bm{\mathcal{k}}}\cdot{\bm{p}}}{\mathcal{k}}\left[\mathbb{P}^{\alpha\beta}+\frac{1}{2}\mathbb{Q}^{\alpha\beta}\right]p^{\beta},$
(49)

where ${\bm{\mathcal{k}}}=\left(\mathcal{k}_{x},\mathcal{k}_{y}\right)$ is the
wavevector associated to ${\bm{\mathcal{x}}}$, and we have introduced the
longitudinal projection operator
$\mathbb{Q}^{\alpha\beta}=\mathcal{k}^{\alpha}\mathcal{k}^{\beta}/\mathcal{k}^{2}$.

##  Appendix B Linear stability analysis of 2-D bulk suspensions

Here, we repeat the calculation presented in Section 3 for the case of
microswimmers suspended in two-dimensional bulk fluid. We set $d=2$ in
(15)–(16), and employ (46), since the form of the dipolar Fourier transform is
the same in all dimensions, as noted in Appendix A. We obtain the following
eigenvalue problems

$\displaystyle\delta\hat{\rho}$
$\displaystyle=\frac{\lambda}{v_{s}k}\frac{b}{\sqrt{1+b^{2}}}\delta\hat{\rho},$
(50) $\displaystyle\delta\hat{\mathcal{U}}^{\alpha}$
$\displaystyle=\frac{Bn\kappa}{v_{s}k}\frac{b^{2}-2\sqrt{1+b^{2}}+2}{b^{3}}\delta\hat{\mathcal{U}}^{\alpha}.$
(51)

which, similarly to the three-dimensional case, are decoupled. The
orientational eigenvalue problem can be re-written as

$\displaystyle\gamma_{2d}\equiv\frac{v_{s}k}{Bn\kappa}=\frac{b^{2}-2\sqrt{1+b^{2}}+2}{b^{3}},$
(52)

where, as before, $b=v_{s}k/(\chi+\lambda)$. In contrast to (22), the solution
to (52) can be found analytically, yielding

Figure 5: The real (left) and imaginary (right) part of the function
$G(\gamma)$ from (33).

$\displaystyle\chi=-\lambda+G(\gamma_{2d})B\kappa n,$ (53)

where $G(\gamma_{2d})$ is given in Eq. (33). The plot of $G$ is shown in Fig.
5, and it has the same qualitative features as $F(\gamma)$ shown in Fig. 2. As
in the 3-D case, the orientational instability only exists for pushers,
$\kappa>0$, and sets in at the largest possible scale, $k\rightarrow 0$. The
corresponding instability condition is given by

$\displaystyle n>n_{c}=\frac{4\lambda}{B\kappa},$ (54)

where the number density $n$ and the dipolar strength $\kappa$ have different
dimensions than their 3-D counterparts. The solution (53) ceases to exist for
$\gamma_{2d}>1/\sqrt{2}$, where one can instead prove stability using an
analysis completely analogous to the 3-D case (see Appendix C.1). The analysis
of density fluctuations is significantly simpler and the corresponding
eigenvalue is readily obtained from (50):

$\displaystyle\chi=-\lambda\pm\sqrt{\lambda^{2}-k^{2}v_{s}^{2}}.$ (55)

Clearly, $\chi$ is negative for all $k$, and thus the suspension is always
stable with respect to small fluctuations in microswimmer density.

##  Appendix C Linear stability analysis using Laplace transforms

###  C.1 Density fluctuations in bulk suspensions

As mentioned in Section 4.1, the linear stability analysis based on the Ansatz
(14) leads to a missing solution of the eigenvalue problem at small spatial
scales. This issue can be tackled by exactly solving the linear dynamics (11)
via Laplace transforms. We introduce the following notation for the forward
Laplace transform

$\displaystyle f(s)$
$\displaystyle=\int_{0}^{\infty}\mathrm{d}tf(t)\mathrm{e}^{-st}$ (56)

where $s$ has a positive real part. For the inverse Laplace transformation, we
use the definition

$\displaystyle f(t)=\frac{1}{2\pi
i}\lim\limits_{T\rightarrow\infty}\int_{\gamma-iT}^{\gamma+iT}\mathrm{d}s\
f(s)\mathrm{e}^{st},$ (57)

where $\gamma$ must be chosen such that the real part of every pole in $f(s)$
is smaller than $\gamma$.

We now consider the case of a 3-D bulk suspension, for which the Laplace
transform of (11) gives

$\displaystyle(s+\lambda+iv_{s}{\bm{p}}\cdot{\bm{k}})\delta\hat{\Psi}-\delta\hat{\Psi}_{0}$
$\displaystyle=\frac{\lambda}{4\pi}\delta\hat{\rho}+\frac{3Bn}{4\pi}ik^{\alpha}p^{\alpha}p^{\beta}\delta\hat{\mathcal{U}}^{\beta},$
(58)

where $\delta\hat{\Psi}_{0}=\delta\hat{\Psi}({\bm{k}},{\bm{p}},t=0)$
represents the initial condition. Dividing the last equation with the linear
operator on the left hand side and integrating over the orientation ${\bm{p}}$
leads to

$\displaystyle\delta\hat{\rho}(s)$
$\displaystyle=\frac{1}{1-\frac{\lambda}{v_{s}k}\arctan(\frac{v_{s}k}{s+\lambda})}\int\mathrm{d}{\bm{p}}\
\frac{\delta\hat{\Psi}_{0}}{s+\lambda+iv_{s}{\bm{p}}\cdot{\bm{k}}}.$ (59)

The expression (59) contains two separate poles

$\displaystyle s_{1}$ $\displaystyle=-\lambda-
iv_{s}{\bm{p}}\cdot{\bm{k}},\quad
s_{2}=-\lambda+\frac{v_{s}k}{\tan(\frac{v_{s}k}{\lambda})},$ (60)

where the second pole appears only in the limited region of parameter space
$v_{s}k/\lambda<\pi/2$. The result of the inverse Laplace transformation (57)
of (59) consists of two terms, one corresponding to each pole:

$\displaystyle\delta\hat{\rho}(t)=\delta\hat{\rho}_{1}(t)+\delta\hat{\rho}_{2}(t)$
(61)

which we now analyse separately. The first term in (61) exists for any
$v_{s}k/\lambda$ and reads

$\displaystyle\delta\hat{\rho}_{1}(t)$ $\displaystyle=\int\mathrm{d}{\bm{p}}\
\Bigg{[}\frac{1}{1-\frac{\lambda}{v_{s}k}\arctan(\frac{ik}{{\bm{p}}\cdot{\bm{k}}})}\mathrm{e}^{-\left(\lambda+iv_{s}{\bm{p}}\cdot{\bm{k}}\right)t}\Bigg{]}\delta\hat{\Psi}_{0}({\bm{k}},{\bm{p}}).$
(62)

The integral over ${\bm{p}}$ is conveniently expressed in spherical
coordinates with its polar axis being along ${\bm{k}}$, and reads

$\displaystyle\delta\hat{\rho}_{1}(t)$ $\displaystyle=\mathrm{e}^{-\lambda
t}\int_{-1}^{1}\mathrm{d}x\ F(x)\mathrm{e}^{-iv_{s}ktx},\quad
F(x)=\frac{\int_{0}^{2\pi}\mathrm{d}\phi\
\delta\hat{\Psi}_{0}({\bm{k}},\phi,\arccos(x))}{1-\frac{\lambda}{v_{s}k}\arctan(\frac{i}{x})}.$
(63)

We observe that $\int_{-1}^{1}\mathrm{d}x|F(x)|<\infty$, provided
$\max_{x}\Big{|}\int_{0}^{2\pi}\mathrm{d}\phi\
\delta\hat{\Psi}_{0}({\bm{k}},\phi,\arccos(x))\Big{|}<\infty$, as expected
from a physical initial condition $\delta\hat{\Psi}({\bm{k}},{\bm{p}})$.
Applying the Riemann-Lebesgue lemma (Bender & Orszag, 1999), we conclude that
$\delta\hat{\rho}_{1}(t)\to 0$, as $t\to\infty$ for $v_{s}k>0$.

The second term in (61) exists only for $v_{s}k/\lambda<\pi/2$ and has the
form

$\displaystyle\delta\hat{\rho}_{2}(t)$
$\displaystyle=\exp\left[-\left(\lambda-\frac{v_{s}k}{\tan(v_{s}k/\lambda)}\right)t\right]\int\mathrm{d}{\bm{p}}\
\Bigg{[}\frac{(\frac{v_{s}k}{\lambda})}{\tan(\frac{v_{s}k}{\lambda})}\frac{1+\tan^{2}(\frac{v_{s}k}{\lambda})}{1+i\tan(\frac{v_{s}k}{\lambda}){\bm{p}}\cdot{\bm{k}}}\Bigg{]}\delta\hat{\Psi}_{0}({\bm{k}},{\bm{p}}).$
(64)

The temporal evolution of (64) follows an exponential decay provided that

$\displaystyle\lambda-\frac{v_{s}k}{\tan(v_{s}k/\lambda)}>0,$ (65)

which always holds for $v_{s}k/\lambda<\pi/2$. $\delta\hat{\rho}_{2}$ thus
corresponds to the solution obtained when using the exponential Ansatz
$\delta\hat{\rho}(t)\sim\exp\\{st\\}$, as was analysed in Section 3.2.

We conclude that, although the use of the exponential Ansatz in (14) does not
allow for retrieving the sub-exponential contributions to the linear dynamics,
the missing part does not carry any new information relevant for the linear
stability analysis.

###  C.2 Orientational instability in a sheet of microswimmers

We now repeat the calculation for the case of transversal velocity
fluctuations in a two-dimensional suspension restricted to a plane in a three-
dimensional fluid. The starting point is the Smoluchowski equation (11) in
2-D:

$\displaystyle\Bigl{[}\partial_{t}+\lambda+iv_{s}{\bm{\mathcal{p}}}\cdot{\bm{\mathcal{k}}}\Bigr{]}\delta\hat{\Psi}$
$\displaystyle=\frac{\lambda}{2\pi}\delta\hat{\rho}+\frac{n}{2\pi}\Bigl{[}2B\mathcal{p}^{\alpha}\mathcal{p}^{\beta}-(1+B)\delta^{\alpha\beta}\Bigr{]}i\mathcal{k}^{\alpha}\delta\hat{\mathcal{U}}^{\beta},$
(66)

Eq. (66) is now transformed using (56), and after multiplying with
$\mathbb{P}^{\alpha\beta}\hat{u}^{\beta}_{\text{plane}}$ and integrating over
the orentation ${\bm{\mathcal{p}}}$, we arrive at

$\displaystyle\delta\hat{\mathcal{U}}_{\perp}^{\alpha}(s)=\frac{1}{1-M_{11}}\int_{0}^{2\pi}\frac{\mathrm{d}\theta}{2\pi}\frac{\mathbb{P}^{\alpha\beta}\hat{u}_{\text{plane}}^{\beta}\delta\hat{\Psi}_{0}({\bm{\mathcal{k}}},\theta)}{s+\lambda+iv_{s}\mathcal{k}\cos\theta},$
(67)

where $M_{11}$ is defined in (29). Similarly to before, the inverse Laplace
transform consists of two parts, each corresponding to one of the poles in Eq.
(67). The first one is the pole corresponding to $M_{11}=1$, which leads to
the dispersion law in Eq. (32). As discussed in the main text, the latter
solution ceases to exist for $2\sqrt{2}v_{s}>B\kappa n$. In this region, the
stability is instead determined by the second pole $s=-\lambda-
iv_{s}\mathcal{k}\cos\theta$, where the time representation of the
corresponding solution is given by

$\displaystyle\delta\hat{\mathcal{U}}_{\perp}^{\alpha}(t)=\mathrm{e}^{-\lambda
t}\int_{0}^{2\pi}\mathrm{d}\theta\
F(\theta)\mathrm{e}^{-iv_{s}\mathcal{k}t\cos\theta},\quad
F(\theta)=\frac{\mathbb{P}^{\alpha\beta}\hat{u}_{\text{plane}}^{\beta}\delta\hat{\Psi}_{0}({\bm{\mathcal{k}}},\theta)}{1-M_{11}|_{b\rightarrow\frac{i}{\cos\theta}}}.$
(68)

Again, it is possible to show that
$\int_{0}^{2\pi}\mathrm{d}\theta|F(\theta)|$ is bounded provided
$2\sqrt{2}v_{s}>B\kappa n$ and
$\delta\hat{\Psi}_{0}({\bm{\mathcal{k}}},\theta)$ satisfies the same physical
requirements as in Appendix C.1. Applying the Riemann-Lebesgue lemma, we
conclude that $\delta\hat{\mathcal{U}}_{\perp}^{\alpha}(t)$ vanishes as
$t\to\infty$, and the system is stable for $2\sqrt{2}v_{s}>B\kappa n$.

##  Appendix D Approximation of $n_{c}$ for a sheet of pushers

In this Appendix, we will derive the approximate expression (37) for the
critical density $n_{c}$ corresponding to the orientational instability
setting in at the scale of particle-particle separation in a 2-D layer of
pusher microswimmers. We set $\mathcal{k}_{c}=\sqrt{\pi^{3}n_{c}}$ and
$n=n_{c}$ in (32), giving

$\displaystyle\chi=-\lambda+G\left(\frac{2v_{s}}{B\kappa
n_{c}}\right)\frac{B\kappa(\pi n_{c})^{3/2}}{2}.$ (69)

We now separate the expression into two approximations valid asymptotically
for small and large values of $v_{s}$. We first set $v_{s}\rightarrow 0$ in
(69), yielding

$\displaystyle\chi(v_{s}\to 0)=-\lambda+\frac{1}{8}B\kappa(\pi n_{c})^{3/2},$
(70)

where $\text{Re}[\chi]=0$ now gives the ‘shaker’ approximation

$\displaystyle n_{c}(v_{s}\to
0)=\frac{4}{\pi}\left(\frac{\lambda}{B\kappa}\right)^{2/3}.$ (71)

For the fast-swimming limit, we use the result (35), _i.e._ , we assume that
the critical density is linearly proportional to the swimming speed,

$\displaystyle
n_{c}(v_{s}\to\infty)=2\sqrt{2}\frac{v_{s}}{B\kappa}+\mathcal{O}(1).$ (72)

Combining (71) and (72) gives the result (37). The latter shows a good
quantitative agreement with the numerical solution of $\text{Re}[\chi]=0$, as
shown in Fig. 6.

Figure 6: Comparison of the approximation (37) with the exact values obtained
by numerical solution of $\text{Re}[\chi]=0$, where $\chi$ is given by Eq.
(69). Our approximation shows good quantitative agreement with the exact
numerical values. The kink in the numerical data is not an artefact and
corresponds to a switch between two eigenvalue branches.

## References

  * Albritton & Ohm (2023) Albritton, Dallas & Ohm, Laurel 2023 On the stabilizing effect of swimming in an active suspension, arXiv: 2205.04922.
  * Aranson (2022) Aranson, Igor S 2022 Bacterial active matter. Rep. Prog. Phys. 85, 076601.
  * Bárdfalvy et al. (2020) Bárdfalvy, Dóra, Anjum, Shan, Nardini, Cesare, Morozov, Alexander & Stenhammar, Joakim 2020 Symmetric mixtures of pusher and puller microswimmers behave as noninteracting suspensions. Phys. Rev. Lett. 125, 018003.
  * Bárdfalvy et al. (2019) Bárdfalvy, Dóra, Nordanger, Henrik, Nardini, Cesare, Morozov, Alexander & Stenhammar, Joakim 2019 Particle-resolved lattice boltzmann simulations of 3-dimensional active turbulence. Soft Matter 15, 7747–7756.
  * Baskaran & Marchetti (2009) Baskaran, Aparna & Marchetti, M. Cristina 2009 Statistical mechanics and hydrodynamics of bacterial suspensions. Proc. Natl. Acad. Sci. USA 106 (37), 15567–15572.
  * Bender & Orszag (1999) Bender, C.M. & Orszag, S.A. 1999 Advanced Mathematical Methods for Scientists and Engineers I: Asymptotic Methods and Perturbation Theory. Springer.
  * Berke et al. (2008) Berke, Allison P., Turner, Linda, Berg, Howard C. & Lauga, Eric 2008 Hydrodynamic attraction of swimming microorganisms by surfaces. Phys. Rev. Lett. 101, 038102\.
  * Blake (1971) Blake, J. R. 1971 A note on the image system for a stokeslet in a no-slip boundary. Math. Proc. Camb. Philos. Soc. 70, 303–310.
  * Bleibel et al. (2017) Bleibel, J., Domínguez, Alvaro & Oettel, M. 2017 Onset of anomalous diffusion in colloids confined to quasimonolayers. Phys. Rev. E 95, 032604.
  * Bleibel et al. (2014) Bleibel, J., Domínguez, A., Günther, F., Harting, J. & Oettel, M. 2014 Hydrodynamic interactions induce anomalous diffusion under partial confinement. Soft Matter 10, 2945–2948.
  * Bleibel et al. (2015) Bleibel, Johannes, Domínguez, Alvaro & Oettel, Martin 2015 3d hydrodynamic interactions lead to divergences in 2d diffusion. J. Phys.: Condens. Matter 27, 194113.
  * Brotto et al. (2013) Brotto, Tommaso, Caussin, Jean-Baptiste, Lauga, Eric & Bartolo, Denis 2013 Hydrodynamics of confined active fluids. Phys. Rev. Lett. 110, 038101.
  * Chen et al. (2012) Chen, Xiao, Dong, Xu, Be’er, Avraham, Swinney, Harry L. & Zhang, H. P. 2012 Scale-invariant correlations in dynamic bacterial clusters. Phys. Rev. Lett. 108, 148101\.
  * Drescher et al. (2011) Drescher, K., Dunkel, J., Cisneros, L. H., Ganguly, S. & Goldstein, R. E. 2011 Fluid dynamics and noise in bacterial cell-cell and cell-surface scattering. Proc. Natl. Acad. Sci. USA 108, 10940.
  * Dunkel et al. (2013) Dunkel, J., Heidenreich, S., Dreschner, K., Wensink, H. H., Bär, M. & Goldstein, R. E. 2013 Fluid dynamics of bacterial turbulence. Phys. Rev. Lett. 110, 228102\.
  * Gachelin et al. (2014) Gachelin, J., Rousselet, A., Lindner, A. & Clement, E. 2014 Collective motion in an active suspension of Escherichia coli bacteria. New J. Phys. 16.
  * Guasto et al. (2010) Guasto, Jeffrey S., Johnson, Karl A. & Gollub, J. P. 2010 Oscillatory flows induced by microorganisms swimming in two dimensions. Phys. Rev. Lett. 105, 168102.
  * ten Hagen et al. (2015) ten Hagen, Borge, Wittkowski, Raphael, Takagi, Daisuke, Kümmel, Felix, Bechinger, Clemens & Löwen, Hartmut 2015 Can the self-propulsion of anisotropic microswimmers be described by using forces and torques? J. Condens. Matter Phys. 27, 194110.
  * Hohenegger & Shelley (2010) Hohenegger, C. & Shelley, M. J. 2010 Stability of active suspensions. Phys. Rev. E 81, 046311.
  * Huang et al. (2021) Huang, Mingji, Hu, Wensi, Yang, Siyuan, Liu, Quan-Xing & Zhang, H. P. 2021 Circular swimming motility and disordered hyperuniform state in an algae system. Proc. Natl. Acad. Sci. USA 118, e2100493118.
  * Ishikawa & Pedley (2008) Ishikawa, Takuji & Pedley, T. J. 2008 Coherent structures in monolayers of swimming particles. Phys. Rev. Lett. 100, 088103.
  * Jeanneret et al. (2019) Jeanneret, Raphaël, Pushkin, Dmitri O. & Polin, Marco 2019 Confinement enhances the diversity of microbial flow fields. Phys. Rev. Lett. 123, 248102.
  * Jepson et al. (2013a) Jepson, Alys, Martinez, Vincent A., Schwarz-Linek, Jana, Morozov, Alexander & Poon, Wilson C. K. 2013a Enhanced diffusion of nonswimmers in a three-dimensional bath of motile bacteria. Phys. Rev. E 88, 041002.
  * Jepson et al. (2013b) Jepson, Alys, Martinez, Vincent A., Schwarz-Linek, Jana, Morozov, Alexander & Poon, Wilson C. K. 2013b Enhanced diffusion of nonswimmers in a three-dimensional bath of motile bacteria. Phys. Rev. E 88, 041002.
  * Jülicher et al. (2018) Jülicher, Frank, Grill, Stephan W & Salbreux, Guillaume 2018 Hydrodynamic theory of active matter. Rep. Prog. Phys. 81 (7), 076601.
  * Kim & Karrila (2005) Kim, S. & Karrila, S. J. 2005 Microhydrodynamics: principles and selected applications. Dover Publications, Inc.
  * Koch & Subramanian (2011) Koch, Donald L. & Subramanian, Ganesh 2011 Collective Hydrodynamics of Swimming Microorganisms: Living Fluids. Annu. Rev. Fluid Mech. 43 (1), 637–659.
  * Krishnamurthy & Subramanian (2015) Krishnamurthy, Deepak & Subramanian, Ganesh 2015 Collective motion in a suspension of micro-swimmers that run-and-tumble and rotary diffuse. J. Fluid Mech. 781, 422–466.
  * Lauga & Powers (2009) Lauga, E. & Powers, T. R. 2009 The hydrodynamics of swimming microorganisms. Rep. Prog. Phys. 72, 096601.
  * Leptos et al. (2009) Leptos, Kyriacos C., Guasto, Jeffrey S., Gollub, J. P., Pesci, Adriana I. & Goldstein, Raymond E. 2009 Dynamics of enhanced tracer diffusion in suspensions of swimming eukaryotic microorganisms. Phys. Rev. Lett. 103, 198103.
  * Liu et al. (2019) Liu, Zhengyang, Zhang, Kechun & Cheng, Xiang 2019 Rheology of bacterial suspensions under confinement. Rheol. Acta 58, 439–451.
  * López et al. (2015) López, Héctor Matías, Gachelin, Jérémie, Douarche, Carine, Auradou, Harold & Clément, Eric 2015 Turning bacteria suspensions into superfluids. Phys. Rev. Lett. 115, 028301.
  * Maitra (2020) Maitra, Ananyo 2020 Stable long-range uniaxial order in active fluids at two-dimensional interfaces. arXiv:2010.15044 .
  * Maitra (2023) Maitra, A. 2023 Two-dimensional long-range uniaxial order in three-dimensional active fluids. Nat. Phys. 19, 733.
  * Marchetti et al. (2013) Marchetti, M. C., Joanny, J. F., Ramaswamy, S., Liverpool, T. B., Prost, J., Rao, Madan & Simha, R. Aditi 2013 Hydrodynamics of soft active matter. Rev. Mod. Phys. 85, 1143–1189.
  * Martinez et al. (2020) Martinez, Vincent A., Clément, Eric, Arlt, Jochen, Douarche, Carine, Dawson, Angela, Schwarz-Linek, Jana, Creppy, Adama K., Škultéty, Viktor, Morozov, Alexander N., Auradou, Harold & Poon, Wilson C. K. 2020 A combined rheometry and imaging study of viscosity reduction in bacterial suspensions. Proc. Natl. Acad. Sci. USA 117, 2326–2331.
  * Mathijssen et al. (2016) Mathijssen, A. J. T. M., Doostmohammadi, A., Yeomans, J. M. & Shendruk, T. N. 2016 Hydrodynamics of micro-swimmers in films. J. Fluid Mech. 806, 35–70.
  * Miño et al. (2013) Miño, G. L., Dunstan, J., Rousselet, A., Clément, E. & Soto, R. 2013 Induced diffusion of tracers in a bacterial suspension: Theory and experiments. J. Fluid Mech. 729, 423.
  * Nejad & Yeomans (2022) Nejad, Mehrana R. & Yeomans, Julia M. 2022 Active extensile stress promotes 3d director orientations and flows. Phys. Rev. Lett. 128, 048001.
  * Pessot et al. (2018) Pessot, Giorgio, Löwen, Hartmut & Menzel, Andreas M. 2018 Binary pusher–puller mixtures of active microswimmers and their collective behaviour. Mol. Phys. 116 (21-22), 3401–3408.
  * Rafaï et al. (2010) Rafaï, Salima, Jibuti, Levan & Peyla, Philippe 2010 Effective viscosity of microswimmer suspensions. Phys. Rev. Lett. 104, 098102.
  * Saintillan (2018) Saintillan, David 2018 Rheology of Active Fluids. Annu. Rev. Fluid Mech. 50 (1), 563–592.
  * Saintillan & Shelley (2008a) Saintillan, David & Shelley, Michael J. 2008a Instabilities and pattern formation in active particle suspensions: Kinetic theory and continuum simulations. Phys. Rev. Lett. 100, 178103.
  * Saintillan & Shelley (2008b) Saintillan, David & Shelley, Michael J. 2008b Instabilities, pattern formation, and mixing in active suspensions. Phys. Fluids 20 (12), 123304\.
  * Saintillan & Shelley (2013) Saintillan, David & Shelley, Michael J. 2013 Active suspensions and their nonlinear models. C. R. Phys. 14, 497 – 517.
  * Salbreux et al. (2009) Salbreux, G., Prost, J. & Joanny, J. F. 2009 Hydrodynamics of cellular cortical flows and the formation of contractile rings. Phys. Rev. Lett. 103, 058102.
  * Singh & Adhikari (2016) Singh, Rajesh & Adhikari, R. 2016 Universal hydrodynamic mechanisms for crystallization in active colloidal suspensions. Phys. Rev. Lett. 117, 228002.
  * Sokolov & Aranson (2012) Sokolov, Andrey & Aranson, Igor S. 2012 Physical properties of collective motion in suspensions of bacteria. Phys. Rev. Lett. 109, 248109.
  * Spagnolie & Lauga (2012) Spagnolie, Saverio E. & Lauga, Eric 2012 Hydrodynamics of self-propulsion near a boundary: predictions and accuracy of far-field approximations. J. Fluid Mech. 700, 105–147.
  * Stenhammar et al. (2017) Stenhammar, Joakim, Nardini, Cesare, Nash, Rupert W., Marenduzzo, Davide & Morozov, Alexander 2017 Role of Correlations in the Collective Behavior of Microswimmer Suspensions. Phys. Rev. Lett. 119 (2), 028005.
  * Subramanian & Koch (2009) Subramanian, Ganesh & Koch, Donald L. 2009 Critical bacterial concentration for the onset of collective swimming. J. Fluid Mech. 632, 359–400.
  * Thutupalli et al. (2018) Thutupalli, Shashi, Geyer, Delphine, Singh, Rajesh, Adhikari, Ronojoy & Stone, Howard A. 2018 Flow-induced phase separation of active particles is controlled by boundary conditions. Proc. Natl. Acad. Sci. USA 115 (21), 5403–5408.
  * Vladescu et al. (2014) Vladescu, I. D., Marsden, E. J., Schwarz-Linek, J., Martinez, V. A., Arlt, J., Morozov, A. N., Marenduzzo, D., Cates, M. E. & Poon, W. C. K. 2014 Filling an emulsion drop with motile bacteria. Phys. Rev. Lett. 113, 268101\.
  * Škultéty et al. (2020) Škultéty, Viktor, Nardini, Cesare, Stenhammar, Joakim, Marenduzzo, Davide & Morozov, Alexander 2020 Swimming suppresses correlations in dilute suspensions of pusher microorganisms. Phys. Rev. X 10, 031059\.
  * Wensink et al. (2012) Wensink, H. H., Dunkel, J., Heidenreich, S., Drescher, K., Goldstein, R. E., Löwen, H. & Yeomans, J. M. 2012 Meso-scale turbulence in living fluids. Proc. Natl. Acad. Sci. USA 109, 14308–14313.
  * Wu & Libchaber (2000) Wu, X. L. & Libchaber, A. 2000 Particle diffusion in a quasi-two-dimensional bacterial bath. Phys. Rev. Lett. 84, 3017–3020.
  * Zantop & Stark (2020) Zantop, Arne W. & Stark, Holger 2020 Squirmer rods as elongated microswimmers: flow fields and confinement. Soft Matter 16, 6400–6412.
  * Zhang et al. (2010) Zhang, H. P., Be’er, Avraham, Florin, E.-L. & Swinney, Harry L. 2010 Collective motion and density fluctuations in bacterial colonies. Proc. Natl. Acad. Sci. USA 107 (31), 13626–13630.