# Collective behavior of squirmers in thin films

Bohan Wu-Zhang, Dmitry A. Fedosov, and Gerhard Gompper  Theoretical Physics of
Living Matter, Institute of Biological Information Processing and Institute
for Advanced Simulation, Forschungszentrum Jülich, 52425 Jülich, Germany
Email<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Bacteria in biofilms form complex structures and can collectively migrate
within mobile aggregates, which is referred to as swarming. This behavior is
influenced by a combination of various factors, including morphological
characteristics and propulsive forces of swimmers, their volume fraction
within a confined environment, and hydrodynamic and steric interactions
between them. In our study, we employ the squirmer model for microswimmers and
the dissipative particle dynamics method for fluid modeling to investigate the
collective motion of swimmers in thin films. The film thickness permits a free
orientation of non-spherical squirmers, but constraints them to form a two-
layered structure at maximum. Structural and dynamic properties of squirmer
suspensions confined within the slit are analyzed for different volume
fractions of swimmers, motility types (e.g., pusher, neutral squirmer,
puller), and the presence of a rotlet dipolar flow field, which mimics the
counter-rotating flow generated by flagellated bacteria. Different states are
characterized, including a gas-like phase, swarming, and motility-induced
phase separation, as a function of increasing volume fraction. Our study
highlights the importance of an anisotropic swimmer shape, hydrodynamic
interactions between squirmers, and their interaction with the walls for the
emergence of different collective behaviors. Interestingly, the formation of
collective structures may not be symmetric with respect to the two walls.
Furthermore, the presence of a rotlet dipole significantly mitigates
differences in the collective behavior between various swimmer types. These
results contribute to a better understanding of the formation of bacterial
biofilms and the emergence of collective states in confined active matter.

##  I Introduction

Collective motion of microswimmers is a popular topic in the field of active
matter due to its wide applicability in the context of both biological and
artificial systems as well as the richness of observed behaviors and physical
mechanisms [1, 2, 3]. A prominent example of the collective behavior of
biological microswimmers is biofilms, which represent complex dynamic
communities of microorganisms at surfaces [4, 5]. Biofilms are often
associated with various infectious diseases, such as dental plaque formation
on teeth [6] and chronic wounds that resist healing [7], driving the research
to better understand surface colonization by bacteria and biofilm development
[5, 8]. In the context of artificial microswimmers, there exist a variety of
active systems, including collectives of diffusiophoretic and thermophoretic
Janus particles [9, 10, 11], active droplets [12, 13], and Quincke rollers
[14]. Studies with artificial microswimmers primarily focus on understanding
the emergence of collective motion and the governing physical mechanisms [1,
2, 3].

One of the prominent examples of collective behavior is motility-induced phase
separation (MIPS), which can occur without any attractive or alignment
interactions between self-propelled particles [15, 16, 17, 18, 19]. MIPS are
characterized by the co-existence of low and high density phases of active
particles, where the latter is generally represented by nearly immobile large
clusters of particles [9, 10, 20]. In the initial stage of biofilm formation,
a collective motility of bacteria known as swarming is frequently observed [5,
21]. In contrast to MIPS clusters, bacterial swarms are stable and highly
mobile aggregates, which generally migrate along surfaces. The formation of
swarms often requires some type of alignment interactions between
microswimmers [22, 23]. Another collective phenomenon observed for bacterial
suspensions is called active turbulence [24, 25, 26], which is similar to the
traditional Kolmogorov-Kraichnan-type hydrodynamic turbulence [27], but occurs
at very low Reynolds numbers. The state of active turbulence is observed for a
sufficiently large density of microswimmers, which each generates a force-
dipolar flow field around it [28, 29, 30, 31]. Finally, the motion of
microswimmers in confinement, such as in channels and pores, can result in a
hydrodynamic instability and the formation of complex flow patterns [24, 32,
33].

Different properties of microswimmers and their environment affect the
emerging collective behavior. These include shape anisotropy [34, 35, 36], the
generated flow field around a swimmer [37, 38, 31], hydrodynamic interactions
[39, 38, 40], the presence of confinement [32, 41], and dimensionality [42,
12]. Biological microswimmers (e.g., bacteria) have complex geometries and
propulsion mechanisms, such that their detailed modeling is computationally
challenging, especially in systems with a large number of swimmers. One of the
popular models of simplified swimmers is the squirmer model [43, 44], which
represents a microswimmer by a spherically- or spheroidally-shaped active
particle with a prescribed slip velocity at its surface [45]. Despite its
simplicity, the squirmer model is flexible enough to capture various swimming
modes, including pushers (e.g., E. coli), neutral microswimmers (e.g.,
Paramecium), and pullers (e.g., Chlamydomonas). In particular, the squirmer
model mimics well both the far- and near-field flow of a variety of realistic
microswimmers.

The squirmer model has been used to study collective behavior of microswimmers
in quasi two-dimensional (2D) systems (i.e., a monolayer of swimmers) [38, 31,
41, 46, 47, 39, 48, 49] and three-dimensional (3D) settings with periodic
boundary conditions (BCs) [50, 51, 52, 53, 37]. In quasi-2D systems, the
positions of spherical squirmers are restricted to a monolayer, while their
orientation can still cover a large range of angles in 3D, depending on their
aspect ratio [38, 41, 49]. As a result, the in-plane velocity of squirmers is
not constant, but has a wide distribution. Nevertheless, these systems show
MIPS at large enough concentrations of squirmers, which is qualitatively
consistent with quasi-2D simulations where the squirmer orientation is
restricted to a plane [46, 47, 48]. Hydrodynamic interactions are found to
suppress MIPS, such that the phase separation for active Brownian particles
takes place at lower volume fractions $\phi$ of particles [38, 40]. Pullers
show the onset of MIPS at lower $\phi$ values in comparison to pushers. For
spheroidal squirmers in quasi-2D systems [38, 31], the collective behavior is
qualitatively similar to spherical squirmers, but a larger aspect ratio of the
spheroidal body favors cluster formation, leading to shape-induced jamming and
alignment. The alignment interaction due to the anisotropy of spheroidal shape
can lead to swarming behavior at moderate $\phi$ values. For the 3D systems
with periodic BCs [50, 51, 52, 53, 37], the MIPS phase also develops at large
volume fractions of squirmers.

In our study, we investigate the collective behavior of spheroidal squirmers
in thin films that are thick enough to allow a full freedom of squirmer
orientation. However, the squirmers are still restricted in their vertical
position, so that they can form only a two-layered structure. Thus, we address
some aspects of biofilm formation beyond the monolayer structure, when the
development of further bacterial layers takes place. In particular, we address
the following questions:

  * •

Do swimmers within a thin film primarily assume orientations parallel or
perpendicular to the wall?

  * •

Under which conditions can swimmers spontaneously leave the wall to form
multi-layered structures?

  * •

Does the collective behavior of squirmers in thin films differ qualitatively
from that of quasi-2D systems?

In our simulations, the volume fraction of squirmers, their swimming mode, and
the strength $\lambda$ of rotlet dipole, which mimics a counter-rotating flow
field of swimmers whose body and flagella bundle rotate in opposite
directions, are varied. In agreement with the quasi-2D systems, pullers
display MIPS phase at lower $\phi$ values in comparison with pushers and
neutral swimmers. At moderate volume fractions, all squirmers with rotlet
dipole and pushers with $\lambda=0$ show swarming behavior. Pullers prefer a
nearly perpendicular orientation to the wall, which leads to the formation of
flower-like structures at low $\phi$ and provides the nucleation for larger
clusters as $\phi$ is increased. Pullers rarely switch between the two walls,
while pushers do so frequently, indicating that a collection of pushers would
quickly develop a multi-layered structure. In all investigated cases, a two-
layered structure is dominant with the orientation of squirmers parallel to
the walls. Interestingly, the presence of rotlet dipole significantly reduces
differences in the collective behavior of squirmers with different swimming
modes. These results provide the first steps in bridging very confined
quasi-2D systems with much less confined situations of the collective behavior
of swimmers in 3D.

The paper is organized as follows. Section II contains all necessary details
about the employed methods and models, including parameters used in
simulations. In Section III.1, structural properties of squirmer suspensions
are analysed, including cluster size distribution, position and orientation of
squirmers within the slit, and the radial distribution function. Dynamic
properties are characterized in Section III.2, where effective rotational
diffusion, mean-squared displacement, and the average speed of squirmers are
presented. The main results are discussed in Section IV, with short
conclusions.

Figure 1: Schematic of the simulation setup. (a) Sketch of a spheroidal
squirmer model. The orientation vector $\mathbf{e}$ is aligned with the
squirmer’s major axis $z$, and $b_{x}$ and $b_{z}=2b_{x}$ denote minor and
major radii of the spheroidal shape. $\mathbf{e}_{\tau}$ and
$\mathbf{e}_{\zeta}$ represent the local normal and tangential unit vectors,
respectively. (b) Several squirmers within a slit of thickness $L_{y}$. Note
that $L_{y}>2b_{z}$, so that the squirmers can freely rotate within the slit
and form a two-layered structure. $L$ is the size of the simulation domain in
the periodic dimensions $x$ and $z$.

##  II Methods and models

###  II.1 Squirmer model

We consider a spheroidal squirmer model, whose surface is described by

$\left(\frac{x}{b_{x}}\right)^{2}+\left(\frac{y}{b_{y}}\right)^{2}+\left(\frac{z}{b_{z}}\right)^{2}=1,$
(1)

where $b_{z}$ and $b_{x}=b_{y}$ are the major and minor radii of the
spheroidal squirmer [see Fig. 1(a)]. The orientational vector e of the
squirmer is aligned with its major axis. The aspect ratio of the spheroidal
shape is set to $b_{z}/b_{x}=2$, which is similar to the aspect ratio of $2-3$
for the body of E. coli bacteria [54]. The squirmer surface is discretized by
$N_{p}$ particles connected by springs into a triangulated network. To
maintain the spheroidal shape, the squirmer consists of a membrane-like
surface with shear and curvature elasticity, and has constraints for its
surface area and enclosed volume [55, 56]. Details of the membrane model are
described in Appendix A.

Locomotion of a spheroidal squirmer is imposed through the prescribed surface
slip velocity given by [38, 57, 58, 59]

$\textbf{u}_{sq}=-B_{1}(\textbf{e}_{\zeta}\cdot\textbf{e}_{z})(1+\beta\zeta)\textbf{e}_{\zeta}+\frac{3\lambda
z_{s}\bar{r}_{s}}{r_{s}^{5}}\textbf{e}_{\varphi},$ (2)

where $\textbf{e}_{z}$, $\textbf{e}_{\zeta}$, and $\textbf{e}_{\varphi}$ are
unit vectors in Cartesian $(x,y,z)$ or spheroidal $(\zeta,\tau,\varphi)$
coordinates, which are related to each other as

$\displaystyle x$
$\displaystyle=c\sqrt{\tau^{2}-1}\sqrt{1-\zeta^{2}}\cos\varphi,$ (3)
$\displaystyle y$
$\displaystyle=c\sqrt{\tau^{2}-1}\sqrt{1-\zeta^{2}}\sin\varphi,$ (4)
$\displaystyle z$ $\displaystyle=c\tau\zeta.$ (5)

Here, $c=\sqrt{b_{z}^{2}-b_{x}^{2}}$, and the spheroidal coordinates have the
ranges $-1\leq\zeta\leq 1$, $1\leq\tau<\infty$ and $0\leq\varphi<2\pi$. The
parameter $B_{1}$ determines self-propulsion speed of the squirmer
$U_{0}=B_{1}\tau_{0}(\tau_{0}-(\tau_{0}^{2}-1)\coth^{-1}\tau_{0})$ with
$\tau_{0}=b_{z}/c$ [45, 60]. The coefficient $\beta$ defines different
swimming modes of the squirmer, including a pusher ($\beta<0$), a neutral
swimmer ($\beta=0$), and a puller ($\beta>0$). The second term in Eq. (2)
defines a rotlet dipole with the strength $\lambda$, where
$\textbf{r}_{s}=(x_{s},y_{s},z_{s})$, $r_{s}=|\textbf{r}_{s}|$,
$\bar{r}_{s}=\sqrt{x^{2}_{s}+y^{2}_{s}}$ with the subscript $s$ denoting
points at the squirmer surface. The rotlet dipole mimics a counter-rotating
flow field, e.g. of E. coli whose body and flagella bundle rotate in opposite
directions [61].

###  II.2 Boundary conditions

The suspension of squirmers is confined between two walls in the $y$
direction, while periodic BCs are applied along the $x$ and $z$ directions.
Fluid flow is modeled by the dissipative particle dynamics (DPD) method [62,
63], which is a particle-based hydrodynamics simulation technique (see
Appendix B for details). Solid walls are represented by frozen DPD particles
within a layer of thickness $r_{c}$ (the cutoff radius for DPD interactions)
with the same number density as that for the DPD fluid. To prevent wall
penetration by fluid particles, a reflective surface is placed at the fluid-
solid interface, where bounce-back reflection of fluid particles is enforced.
No-slip BCs at the walls are imposed by the dissipative interaction between
fluid and frozen-wall particles, and through the bounce-back reflection at the
interface.

To restrict the motion of squirmers between the two walls, particles within
the squirmer discretization are subject to the repulsive Lennard-Jones (LJ)
potential

$U_{LJ}(r)=4\epsilon_{LJ}\left[\left(\frac{\sigma_{LJ}}{r}\right)^{12}-\left(\frac{\sigma_{LJ}}{r}\right)^{6}\right]$
(6)

at the fluid-solid interface. Here, $\epsilon_{LJ}$ is the potential strength,
$\sigma_{LJ}$ sets a characteristic repulsion length, and $r$ is the distance
to the wall. Note that only the repulsive part of the LJ potential is
considered by setting the cutoff distance to $2^{1/6}\sigma_{LJ}$. Excluded-
volume interactions between different squirmers are also implemented through
the LJ interactions, where $r$ becomes a distance between two surface
particles belonging to distinct squirmers.

Squirmers are submerged within a DPD fluid, and also filled by DPD fluid
particles due to their membrane-like representation. The membrane surfaces of
all suspended squirmers serve as a boundary separating DPD particles inside
and outside of the membranes. This is achieved through the reflection of fluid
particles at the membrane surfaces from inside and outside. Note that the
dissipative and random forces between the internal and external fluid
particles are deactivated, and only the conservative force is employed to
maintain uniform fluid pressure across the membranes. To enforce the slip
velocity $\textbf{u}_{sq}$ at the squirmer surface, the dissipative
interaction (see Appendix B) between the squirmer particles and those of the
surrounding fluid is altered as follows

$\textbf{F}^{D}(r_{ij})=-\gamma
W^{D}(r_{ij})(\hat{\textbf{r}}_{ij}\cdot\textbf{v}_{ij}^{\ast})\hat{\textbf{r}}_{ij},\quad\quad\textbf{v}_{ij}^{\ast}=\textbf{v}_{i}-\textbf{v}_{j}+\textbf{u}_{sq}^{i},$
(7)

where $\hat{\textbf{r}}_{ij}=(\textbf{r}_{i}-\textbf{r}_{j})/r_{ij}$,
$r_{ij}=|\textbf{r}_{i}-\textbf{r}_{j}|$, $\textbf{u}_{sq}^{i}$ is the slip
velocity at the position of squirmer particle $i$, while $j$ corresponds to an
outer-fluid DPD particle. Furthermore, the friction coefficient $\gamma$
between fluid and squirmer particles is properly adjusted [55] to ensure the
imposition of $\textbf{u}_{sq}$ at the squirmer surface.

###  II.3 Simulation setup and parameters

The simulation setup corresponds to a domain of dimensions $L\times
L_{y}\times L$, where $L=28b_{z}$ and $L_{y}=2.5b_{z}$ ($b_{z}=4$ in
simulations), see Fig. 1(b). The number density of fluid particles is
$n_{f}=320/b_{z}^{3}$, with a particle mass $m=1$. The energy unit is
$k_{B}T=1$. Parameters selected for the DPD fluid (see Appendix B) yield a
fluid dynamic viscosity of $\eta=403.2\sqrt{mk_{B}T}/b_{z}^{2}$.

Each squirmer in simulations consists of $N_{p}=1024$ particles. Excluded-
volume interactions between different squirmers and between squirmers and the
walls are implemented through the LJ potential with parameters
$\epsilon_{LJ}=k_{B}T$ and $\sigma_{LJ}=0.125b_{z}$. The computed three-
dimensional rotational diffusion coefficient around the major axis of a
spheroid without confinement is $D_{R}=3.28\times
10^{-4}\sqrt{k_{B}T/m}/b_{z}$, which is close to the theoretical prediction of
$D_{R}=3.52\times 10^{-4}\sqrt{k_{B}T/m}/b_{z}$ [45]. To present simulation
results, the half length $b_{z}$ of the semi-major axis is chosen as a length
scale, and the characteristic rotational time $\tau_{R}=1/D_{R}$ as a time
scale. In all simulations, the activity parameter $B_{1}$ is fixed at
$56.8b_{z}D_{R}$. We have verified that the swimming velocity of a single
squirmer is $U_{0}=45.5b_{z}D_{R}$, in agreement with the theoretical
prediction [45]. Using these values, we can define a dimensionless Péclet
number $Pe=U_{0}/(2b_{z}D_{R})\approx 28.4$. Furthermore, the Reynolds number
$Re=2b_{z}U_{0}mn_{f}/\eta\approx 0.03$ is small enough to eliminate possible
inertial effects.

Different volume fractions $\phi=4N_{sq}\pi
b_{x}b_{y}b_{z}/(3L^{2}L_{y})\in\\{0.18,0.35,0.44,0.56\\}$ are considered,
corresponding to different numbers of squirmers
$N_{sq}\in\\{228,456,576,720\\}$. The three values of active stress
$\beta\in\\{-5,0,5\\}$ represent simulation systems with pusher, neutral, and
puller swimmers, respectively. Finally, we also consider two values of the
non-dimensional rotlet dipole strength
$\tilde{\lambda}=\lambda/(b_{z}^{4}D_{R})\in\\{0,133.5\\}$. For
$\tilde{\lambda}=133.5$, a single squirmer moves in a circular trajectory at a
wall with a radius of approximately $2b_{z}$. For comparison, the rotlet
dipole strength of $\tilde{\lambda}\in\\{0,561.1\\}$ has been used in Ref.
[31], where squirmers were confined to a single layer. All simulations are
first run for a time of at least $1.5\tau_{R}$ to reach a steady state, and
afterwards various structural and dynamical characteristics of the system are
measured during the time $1.5\tau_{R}$.

###  II.4 Mapping to experimental systems

Squirmer parameters used in simulations can also be mapped to the properties
of real swimmers. In the far-field approximation, the hydrodynamic field of a
microswimmer is dominated by its force-dipole strength
$\chi=f_{d}l_{d}/(8\pi\eta)$ [64, 65, 66], where $f_{d}$ and $l_{d}$ are the
characteristic force and length of the dipole. The active stress $\beta$ of a
spheroidal squirmer can be expressed as [45]

$\beta=-\frac{\chi}{U_{0}(b_{z}^{2}-b_{x}^{2})}\frac{[3\tau_{0}+(1-3\tau_{0}^{2})\coth^{-1}\tau_{0}][\tau_{0}-(\tau_{0}^{2}-1)\coth^{-1}\tau_{0}]}{2/3-\tau_{0}^{2}+(\tau_{0}^{2}-1)\coth^{-1}\tau_{0}}.$
(8)

Average properties of E. coli bacteria swimming in water correspond to
$U_{0}=29$ $\mu$m/s, $f_{d}=0.42$ pN, $l_{d}=1.9$ $\mu$m, $b_{z}=1.5$ $\mu$m,
$b_{x}=0.5$ $\mu$m, and $\eta=10^{-3}$ Pa$\cdot$s [65, 61, 54], resulting in
$\beta\approx-3$. Note that a direct comparison is possible only in the far-
field limit, while the near-field flow of each swimmer depends on its
geometric and propulsion details.

##  III Results

In our investigation, we analyze the collective structural and dynamical
properties of swimmer ensembles in thin fluid films as a function of volume
fraction $\phi$ and squirmer characteristics of active stress $\beta$ and
rotlet dipole strength $\tilde{\lambda}$. Snapshots of the emergent structures
are displayed in Fig. 2. We characterize these phases through the analysis of
cluster sizes and radial distribution functions of the squirmers. Furthermore,
the distribution of squirmers along the $y$-direction (perpendicular to the
walls) and their orientation are considered to distinguish a two-layered
arrangement with an average orientation within the $x-z$ plane from the
stacked packing with the orientation perpendicular to the walls. Finally,
dynamical properties, such as velocity distribution, and effective rotational
and translational diffusivities, are computed to characterize squirmer
motility within collective states.

An important question for any two-layered structure is to which extent the two
layers interact with and affect each other. The snapshots in Fig. 2 nicely
demonstrate that the correlation between the two layers strongly depends on
volume fraction and strength of the rotlet dipole. Without rotlet dipole, i.e.
for $\tilde{\lambda}=0$, the layer correlation is very weak for small $\phi$,
but becomes very significant for large $\phi$, where the clusters in the top
and bottom layers are essentially in registry. With strong rotlet dipole, i.e.
for $\tilde{\lambda}=133.5$, the situation is very different, as there is
hardly any visible correlation between the two layers.

Previous studies of structure formation of artificial and biological
microswimmers in two- and three-dimensional systems [1, 31, 2, 67, 68] report
the existence of various phases, including

  * •

a gas of small clusters, where the distribution of swimmers is homogeneous and
dynamic clusters are formed by a few swimmers;

  * •

large clusters or motility-induced phase separation (MIPS), where the cluster
size is often comparable with the size of the entire system; such large
clusters are nearly immobile;

  * •

swarming and flocking, which is characterized by the collective locomotion of
dense swimmer clusters with swirling and streaming patterns.

A similar behavior is found in our system (see Fig. 2); however, the
possibility of the formation of two distinct layers at the walls gives rise to
novel structures and dynamics.

Figure 2: Snapshots of the emergent structures for different active stresses
$\beta$, volume fractions $\phi$, and rotlet dipole strengths
$\tilde{\lambda}$. (a) Simulated structures for $\tilde{\lambda}=0$,
$\phi=0.18$ and $0.35$, and $\beta\in\\{0,\pm 5\\}$. Squirmers in the upper
(bottom) half of the slit are colored in orange (blue). See also corresponding
Movies S1 - S6. (b) Squirmer structures within the upper and bottom layers of
the simulated system for $\phi=0.35$, $\beta=5$, and $\tilde{\lambda}=0$. See
also Movie S6. (c) Snapshots of the simulated systems for
$\tilde{\lambda}=133.5$, $\phi=0.35$ and $\beta\in\\{0,\pm 5\\}$. See also
Movies S7 - S9. (d) Squirmer structures within the upper and bottom layers of
the slit for $\phi=0.35$, $\beta=5$, and $\tilde{\lambda}=133.5$. Squirmers
belong to the upper (orange) or bottom (blue) layers, when the $y$-coordinate
of their center of mass is in the upper or lower half of the slit,
respectively. See also Movie S9. Figure 3: Cluster size distributions for
different volume fractions of squirmers. (a) $\phi=0.18$, (b) $\phi=0.35$, and
(c) $\phi=0.56$. The curves with different colors represent various swimming
modes: pusher with $\beta=-5$ (red), neutral with $\beta=0$ (blue), and puller
with $\beta=5$ (green). The blue solid lines are the power-law fits, while the
black-dotted lines are the fits using Eq. (11) with the parameters shown in
Table 1. (d) The average cluster size for various conditions as a function of
volume fraction. The average cluster size is calculated using Eq. (10). All
plots are for the simulations without rotlet dipole ($\tilde{\lambda}=0$),
while the insets show data from the corresponding simulations with
$\tilde{\lambda}=133.5$.

###  III.1 Structural properties

####  III.1.1 Cluster size distribution

The cluster-size distribution function $\mathcal{N}(n)$ is calculated as

$\mathcal{N}(n)=\frac{1}{N_{sq}}np(n),$ (9)

where $\mathcal{N}(n)$ represents the fraction of squirmers belonging to
clusters of size $n$, and $p(n)$ denotes the number of clusters of size $n$.
The distribution is normalized, such that $\sum_{n}\mathcal{N}(n)=1$.
Different squirmers belong to the same cluster when the nearest surface-to-
surface distance $d_{s}$ between them satisfies $d_{s}/b_{z}<0.25$. The
average cluster size $\langle n\rangle$ is then [68]:

$<n>=\sum_{n}n\mathcal{N}(n).$ (10)

The cluster-size distribution is used to distinguish the different collective
phases [68, 21] mentioned above. For the gas of small clusters,
$\mathcal{N}(n)$ exhibits an exponential decay, as shown in Fig. 3(a) for
$\phi=0.18$. For MIPS or large clusters, a bimodal distribution of
$\mathcal{N}(n)$ is observed, with the second peak at large $n$ signaling the
presence of large clusters, see Fig. 3(b,c) for $\phi=0.35$ and $\phi=0.56$.
Finally, swarming can be characterized by a power-law decay of
$\mathcal{N}(n)$ with an exponential cutoff and without the presence of a
distinct peak [47, 31]. This implies a slow decline of cluster sizes over a
broad range, compared to the fast decay for a gas of small clusters.

$\phi$ | $\beta$ | $\tilde{\lambda}$ | $A$ | $\alpha$ | $\gamma$ | phase
---|---|---|---|---|---|---
$0.18$ | $-5$ | $0$ | $0.33$ | $5.07$ | $0.19$ | gas
$0.18$ | $0$ | $0$ | $0.26$ | $13.1$ | $0.57$ | gas
$0.18$ | $5$ | $0$ | $0.09$ | $-$ | $0.49$ | MIPS
$0.18$ | $-5$ | $133.5$ | $0.32$ | $7.1$ | $0.4$ | gas
$0.18$ | $0$ | $133.5$ | $0.31$ | $7.61$ | $0.43$ | gas
$0.18$ | $5$ | $133.5$ | $0.34$ | $8.55$ | $0.54$ | gas
$0.35$ | $-5$ | $0$ | $0.07$ | $49.4$ | $0.4$ | swarming
$0.35$ | $0$ | $0$ | $0.1$ | $-$ | $1.1$ | MIPS
$0.35$ | $5$ | $0$ | $0.05$ | $-$ | $0.99$ | MIPS
$0.35$ | $-5$ | $133.5$ | $0.07$ | $45$ | $0.45$ | swarming
$0.35$ | $0$ | $133.5$ | $0.07$ | $72.5$ | $0.58$ | swarming
$0.35$ | $5$ | $133.5$ | $0.07$ | $90$ | $0.62$ | swarming
$0.44$ | $-5$ | $0$ | $0.04$ | $-$ | $0.89$ | MIPS
$0.44$ | $0$ | $0$ | $0.02$ | $-$ | $0.93$ | MIPS
$0.44$ | $5$ | $0$ | $0.03$ | $-$ | $1.28$ | MIPS
$0.44$ | $-5$ | $133.5$ | $0.03$ | $-$ | $0.62$ | MIPS
$0.44$ | $0$ | $133.5$ | $0.03$ | $-$ | $0.69$ | MIPS
$0.44$ | $5$ | $133.5$ | $0.03$ | $-$ | $0.83$ | MIPS
$0.56$ | $-5$ | $0$ | $0.01$ | $-$ | $1.02$ | MIPS
$0.56$ | $0$ | $0$ | $0.03$ | $-$ | $1.64$ | MIPS
$0.56$ | $5$ | $0$ | $0.01$ | $-$ | $0.81$ | MIPS
$0.56$ | $-5$ | $133.5$ | $0.01$ | $-$ | $0.9$ | MIPS
$0.56$ | $0$ | $133.5$ | $0.01$ | $-$ | $0.94$ | MIPS
$0.56$ | $5$ | $133.5$ | $0.01$ | $-$ | $1.04$ | MIPS
Table 1: Parameters of the cluster size distributions, obtained from fitting
the simulation result with Eq. (11), for various packing fractions $\phi$,
active stresses $\beta$, and dimensionless rotlet dipole strengths
$\tilde{\lambda}=\lambda/(b_{z}^{4}D_{R})$. A ”$-$” sign in the
$\alpha$-column indicates that the exponential term in Eq. (11) is omitted, so
that the fitting function becomes $f(x)=Ax^{-\gamma}$. The last column
provides the classification of different simulation cases into the defined
phases.

The described characteristics of different phases can be extracted by fitting
$\mathcal{N}(n)$ with the function [47]

$f(x)=Ax^{-\gamma}\exp^{-x/\alpha},$ (11)

where $A$, $\gamma$, and $\alpha$ are the fitting parameters. Table 1 presents
these parameters for different simulated conditions.

For $\phi=0.18$, the simulated systems exhibit a homogeneous gas-like phase,
except for the case of pullers ($\beta=5$) without rotlet dipole
($\tilde{\lambda}=0$), shown in Fig. 2(a). Here, the exponential term
dominates, as can be seen well in Fig. 3(a). At low packing fractions of
squirmers, there is a limited chance for squirmer collisions, leading to a
gas-like phase, compare also Fig. 2. For $\beta=5$ at $\phi=0.18$,
$\mathcal{N}(n)$ resembles a bimodal distribution, indicating the existence of
large clusters ($n\approx 60$), whose formation is primarily governed by the
attractive hydrodynamic field around pullers.

Swarming-like behavior is observed in several cases at $\phi=0.35$,
characterized by the relatively large values of $\alpha$ and moderate values
of $\gamma\in[0.4,0.7]$, indicating the dominance of the power-law term in Eq.
(11). Note that for the cases of neutral swimmers and pullers (without rotlet
dipole), the value of $\gamma$ is close to unity, and $\mathcal{N}(n)$
exhibits a bimodal distribution in Fig. 3(b), suggesting the presence of a
MIPS phase. For even larger volume fraction $\phi=0.56$, Fig. 3(c) shows a
prominent peak in $\mathcal{N}(n)$ at large $n$ for all cases, which is the
main characteristic for the MIPS phase. This is also confirmed in Table 1
through $\gamma$ values close to unity. At large packing fractions, steric
interactions between squirmers dominate due to crowding, so that the
differences in $\mathcal{N}(n)$ for different swimming modes nearly disappear.

For squirmers without rotlet dipole ($\tilde{\lambda}=0$), the swimming mode
affects their collective behavior when the volume fraction is small enough,
i.e. $\phi\lesssim 0.4$. In this case, the local hydrodynamic flow field
generated by the squirmers is relevant for cluster formation and dynamics, in
agreement with previous studies [38, 31]. However, when the rotlet dipole is
activated ($\tilde{\lambda}=133.5$), differences in $\mathcal{N}(n)$ for
various active stresses $\beta$ nearly disappear, as displayed in the insets
of Fig. 3. With rotlet dipole, we obtain the gas phase at $\phi=0.18$,
swarming at $\phi=0.35$, and MIPS at $\phi=0.56$, independently of active
stress – see also Fig. 2(c). Note that the systems with $\phi=0.56$ are very
crowded, and MIPS is more difficult to recognize for $\tilde{\lambda}=133.5$
in Fig. 2(c) than for $\tilde{\lambda}=0$ in Fig. 2(a). Additional simulations
(not shown) indicate that the suspension of squirmers with
$\tilde{\lambda}=133.5$ transits from swarming to the MIPS at $\phi\approx
0.4$.

The average cluster size $\langle n\rangle$ as a function of $\phi$ is shown
in Fig. 3(d). For $\phi\gtrsim 0.18$, $\langle n\rangle$ increases rapidly
with increasing volume fraction of squirmers for all $\beta$ and $\lambda$
values. The simulation systems with different active stresses result in a
similar range of average cluster sizes for both rotlet dipole strengths.

Figure 4: Distributions $p_{y}$ of the squirmers’ centre-of-mass position
along the slit height for volume fractions (a) $\phi=0.18$ and (b)
$\phi=0.56$. The walls are located at $y/b_{z}=\pm 1.25$. Colors indicate the
swimming modes: pusher with $\beta=-5$ (red), neutral with $\beta=0$ (blue),
and puller with $\beta=5$ (green). Main figures are for the simulations
without rotlet dipole $\tilde{\lambda}=0$, while the insets show the
corresponding simulations with rotlet dipole $\tilde{\lambda}=133.5$.

####  III.1.2 Squirmer distribution between the walls

Figure 4 shows distributions $p_{y}$ of the squirmers’ center-of-mass position
along the $y$ direction between the two confining walls. Notably, the
squirmers exhibit a distinct affinity for the walls, especially in the case of
pushers, as indicated by the pronounced peaks at $y/b_{z}\approx\pm 0.6$. This
phenomenon is primarily attributed to the well-known wall-trapping effect,
which is on the one hand due to motion, steric interactions, and slow
reorientation [69, 70, 71], on the other hand due to hydrodynamic attraction
[72, 66, 65].

Pushers ($\beta=-5$) exhibit a distinctive two-layered structure with the
orientation parallel to the walls, which is reminiscent of E. coli behavior in
in vitro observations [73, 74]. For pushers, the far-field flow re-orients the
squirmers parallel to the wall [72, 66, 65]. Furthermore, the elongated
spheroidal body of the squirmers favors their orientation parallel to the
walls [75, 76]. For pullers ($\beta=5$), Fig. 4(a) show two small peaks at
$y/b_{z}\approx\pm 0.2$ at small volume fraction $\phi=0.18$, which indicates
a nearly perpendicular orientation with respect to the walls [see also
corresponding conformations in Fig. 2(b)]. Note that the far-field flow
generated by pullers favors such a perpendicular-to-the-wall orientation [72,
66, 65], while the spheroidal shape promotes parallel-to-the-wall orientation.
Thus, pullers in Fig. 4(a) show their major peaks in $p_{y}$ at
$y/b_{z}\approx\pm 0.5$, such that their most probable orientation has about
$40$ degree angle with respect to the walls [see also Fig. 5(a)]. $p_{y}$ for
neutral swimmers at $\phi=0.18$ is close to that for pushers.

At the higher volume fraction of $\phi=0.56$, differences in $p_{y}$ for the
three swimming modes significantly diminish, and all cases essentially show a
two-layered structure, see Fig. 4(b). Note that the two peaks in $p_{y}$ at
$y/b_{z}\approx\pm 0.6$ become broader than those for $\phi=0.18$, which can
be attributed to an increasing importance of steric interactions at large
volume fractions, such that close packing introduces a larger deviation to the
two-layered structure. Insets in Fig. 4 demonstrate that the activation of
rotlet dipole mitigates the influence of active stress on $p_{y}$, in
agreement with the results for the cluster-size distribution in Section
III.1.1.

Figure 5: Orientational distribution $p_{\theta}$ of squirmers, where the
angle $\theta$ is between the swimmer orientation vector and the $y$ axis
(wall normal axis). $\theta=0$ and $\theta=180$ correspond to the orientation
perpendicular to the walls, while $\theta=90$ represent the orientation
parallel to the walls, see the insets. $p_{\theta}$ is presented for
$\phi=0.18$ with the rotlet dipole strength (a) $\tilde{\lambda}=0$ and (b)
$\tilde{\lambda}=133.5$. The three different $\beta$ values correspond to
pushers ($\beta=-5$), neutral swimmers ($\beta=0$), and pullers ($\beta=5$).
The shaded blue areas in (a) mark the angle ranges $[0,25]$ and $[155,180]$
used to compute the fraction of squirmers with a perpendicular orientation to
the walls.

####  III.1.3 Squirmer orientation in the slit

The orientational distribution function $p_{\theta}(\theta)$ of the angle of
the squirmer orientation vector e with the wall normal ($y$ axis) are
presented in Fig. 5. A perpendicular orientation to the walls corresponds to
$\theta=0$ and $\theta=180$, a parallel-to-the-wall orientation to
$\theta=90$. Figure 5(a) shows that pushers and neutral swimmers with
$\tilde{\lambda}=0$ at $\phi=0.18$ are primarily oriented parallel to the
wall, in agreement with the results in Fig. 4(a). In contrast, pullers show
several peaks at $\theta\approx 20^{o}\,\&\,160^{o}$ (small peaks) and
$\theta\approx 60^{o}\,\&\,120^{o}$ (large peaks), which correspond to a
nearly perpendicular-to-the-wall orientation and a tilted orientation with
about $40^{o}$ angle to the walls, respectively [compare also snapshots in
Fig. 2(b)]. These results are consistent with the position distributions
$p_{y}$ discussed in Section III.1.2.

Figure 5(b) shows $p_{\theta}$ for a non-zero rotlet dipole
($\tilde{\lambda}=133.5$) at $\phi=0.18$, which again demsontrates a
significant reducuction of differences between different swimming modes due to
the rotlet dipole. In particular, for all $\beta$, a single peak in
$p_{\theta}$ centered around $90^{o}$ is observed, confirming the preferred
squirmer orientation parallel to the walls, as illustrated in Fig. 2(d). As
the volume fraction of squirmers increases, the dominant role of steric
interactions also leads to a diminished effect of different swimming modes
(even for $\tilde{\lambda}=0$) with the formation of a two-layered structure.

Figure 6: Probability $P_{\perp}$ of squirmers to be aligned with the $y$ axis
(i.e., wall normal vector) for different volume fractions and swimming modes
at $\tilde{\lambda}=0$. Data for pushers with $\beta=-5$ (red), neutral
squirmers with $\beta=0$ (blue), and pullers with $\beta=5$ (green). The
probability is calculated by integrating the orientation distributions in Fig.
5(a) over the ranges $[0,\,25]$ and $[155,\,180]$ degrees.

The probability $P_{\perp}$ of squirmers to be oriented nearly perpendicular
to the walls is displayed in Fig. 6 for different values of $\beta$ as a
function of $\phi$. The probability is computed by integrating the orientation
distributions over the ranges $[0,\,25^{o}]$ and $[155^{o},\,180^{o}]$, as
indicated by the blue shaded areas in Fig. 5(a). At low $\phi$, the
probability of perpendicular orientation to the walls is zero for pushers and
neutral squirmers, while $P_{\perp}\approx 0.11$ for pullers. This simply
confirms the tendency of pushers and neutral swimmers to align with the walls,
due to their hydrodynamic interactions with the walls. Interestingly, the
majority of pullers does not have a perpendicular-to-the-wall orientation,
despite of such predictions for a single puller [72, 66]. This is due to
collective effects between pullers, which will be discussed in Section III.1.4
below.

As the volume fraction of squirmers increases, $P_{\perp}$ for pullers
decreases because steric interactions between squirmers become dominant,
forcing the formation of a two-layered structure, as discussed above. For
neural squirmers, $P_{\perp}$ first increases with increasing $\phi$, but then
closely follows $P_{\perp}$ for pullers when $\phi\gtrsim 0.4$. For pushers,
there is only a slight increase in $P_{\perp}$ with increasing $\phi$.
Nevertheless, we expect that for large volume fractions $\phi\gtrsim 0.6$,
orientational differences between squirmers with various swimming modes
essentially disappear.

Figure 7: (a) 2D radial distribution function $g(r)$ of squirmers with
$\phi=0.18$ and $\tilde{\lambda}=0$, for pushers (red), neutral squirmers
(blue), and pullers (green). The locations of peaks of $g(r)$ for pullers
correlate well with local structural arrangements illustrated next to the
peaks. The inset presents the corresponding $g(r)$ functions for
$\tilde{\lambda}=133.5$ and the same $\phi$. (b) Snapshot of the simulated
system for $\phi=0.18$, $\tilde{\lambda}=0$, and $\beta=5$, which illustrates
the most frequent squirmer structures, see also Movie S3. (c) $g(r)$ of
squirmers for $\phi=0.56$ and $\tilde{\lambda}=0$. The legend for the curves
is the same as in (a).

####  III.1.4 2D radial distribution function parallel to the walls

To characterize the internal structure of the suspension of squirmers along
the walls, we compute the 2D radial distribution function

$g(r)=\frac{M}{N_{sq}^{2}}\left\langle\sum_{i}\sum_{j\neq
i}\delta(\bm{r}-\bm{r}_{ij})\right\rangle,$ (12)

where $\bm{r}_{ij}$ is the vector between two centers of mass of squirmer
pairs, and $M=A/(2\pi r)$ is a normalization factor ($A$ is the area of the
considered 2D plane) such that $g(r)\to 1$ for $r\to\infty$. Note that $g(r)$
is calculated in 2D within the $x-z$ plane, where the measurements are
performed within the two layers (upper and lower halves of the slit)
separately. This is a reasonable simplification due to the prevalence of a
two-layered structure as shown in Fig. 4.

Figure 7(a) presents $g(r)$ for different swimming modes at low volume
fraction $\phi=0.18$ and without rotlet dipole, $\tilde{\lambda}=0$. Only the
suspension of pullers with $\beta=5$ exhibits pronounced peaks in $g(r)$,
which represent the most frequent structural elements illustrated in Fig.
7(b). The existence of structure for pullers, but not for pushers and neutral
squirmers, is primarily related to the fact that pullers form large clusters
already at $\phi=0.18$, while the other squirmer types yield a gas-like phase,
see Table 1. Interestingly, the x-y perspective plot in Fig. 2(b) shows that
pullers frequently form flower-like arrangements with one or two squirmers
oriented perpendicular to the walls and surrounded by several ”petal”
squirmers [also illustrated in Fig. 7(b)]. These structures form due to the
attractive flow field generated by pullers, which slide along the walls until
they collide and form the flower-like structures through the attractive
hydrodynamic interactions. When the rotlet dipole is turned on at $\phi=0.18$
[see the inset in Fig. 7(a)], the suspension of pullers becomes gas-like and
thus looses its internal structure.

Figure 7(c) shows $g(r)$ at the high volume fraction $\phi=0.56$, where
several peaks are observed for all squirmer types. At this high volume
fraction, all squirmer suspensions are in MIPS phase (see Table 1), so that
large clusters are present with an internal fluid structure. Furthermore, due
to the dominance of steric interactions, radial distribution functions are
again very similar for different $\beta$, indicating that the internal
structure is nearly independent of the swimming mode at large $\phi$.

Figure 8: Distribution $p_{\alpha}$ of relative angle between two neighboring
squirmers for different swimming modes at $\phi=0.18$. $\alpha$ is the angle
between orientation vectors of a pair of squirmers, whose separation is
smaller than $3b_{z}$. (a) $p_{\alpha}$ for the case of no rotlet dipole
($\tilde{\lambda}=0$), with the insets illustrating configurations of
squirmers at the two peaks with $\alpha\approx 63^{o}$ and $\alpha\approx
113^{o}$. (b) Relative angle distribution for $\tilde{\lambda}=133.5$.

####  III.1.5 Angle between two neighboring squirmers

Previous simulation studies [45, 46] for a monolayer of squirmers in a thin
fluid between two walls have analyzed the angles between pairs of squirmers
for different swimming modes. Theers et al. [45] report that pairs of pullers
tend to adopt a fixed relative angle after their initial encounter, which is
equal to approximately $45^{o}$. The stable angle formation does not occur for
neutral squirmers and pushers. Similarly, Kyoya et al. [46] suggest that
neutral squirmer pairs tend to align with each other, puller pairs assume a
stable angle between them in the range of $0^{o}$ to $90^{o}$, and pusher
pairs do not show any order and swim away from each other after a brief
encounter.

We calculate distributions $p_{\alpha}$ of the relative angle $\alpha$ between
orientation vectors of two neighboring squirmers within a defined cutoff
radius as

$p_{\alpha}=\frac{1}{N_{sq}C}\left\langle\sum_{i}\sum_{j\neq
i}^{|\bm{r}_{j}-\bm{r}_{j}|\leq
3b_{z}}\delta(\alpha-\alpha_{ij})\right\rangle,$ (13)

where and $C$ is the normalization factor so that the integral of $p_{\alpha}$
is equal to unity. Figure 8(a) shows the distribution $p_{\alpha}$ for various
$\beta$ at $\phi=0.18$. $p_{\alpha}$ for pullers displays two peaks located at
$\alpha\approx 63^{o}$ and $\alpha\approx 113^{o}$, where the former value is
not far from $45^{o}$ found for a monolayer of squirmers in Ref. [45].
Notably, the peak at $\alpha\approx 63^{o}$ is quite prominent, since it
represents the flower-like arrangements illustrated in Fig. 7(b). However, the
distribution for pullers spans a wide range of relative angles. Neutral
squirmers and pushers yield broad and nearly flat distributions, indicating
that no stable structures are present. The existence of two small peaks at
$\alpha\approx 25^{o}$ and $\alpha\approx 155^{o}$ suggests that pushers and
neutral squirmers tend to swim parallel to each other, which is consistent
with previous predictions for a squirmer pair [58, 77].

The distribution of $p_{\alpha}$ for an active rotlet dipole with
$\tilde{\lambda}=133.5$ is shown in Fig. 8(b). For all $\beta$ values,
$p_{\alpha}$ is centered at $\alpha\approx 90^{o}$, indicating that
hydrodynamic interactions favor the perpendicular orientation between squirmer
pairs. Note that for pushers, the central region in $p_{\alpha}$ is flatter
than that for pullers and neutral squirmers.

Figure 9: Effective rotational diffusion $D_{R}$ as a function of $\phi$ for
(a) $\tilde{\lambda}=0$ and (b) $\tilde{\lambda}=133.5$. $D_{R}$ is obtained
by fitting auto-correlation functions of squirmer orientation using Eq. (14),
and normalized by the rotational diffusion $D_{R}^{0}$ of a spheroid. Data for
pushers (red), for neutral squirmers (blue), and for pullers (green).

###  III.2 Dynamical properties

####  III.2.1 Effective rotational diffusion coefficient

To characterize the rotational properties of squirmers, we compute their
effective rotational diffusion coefficient $D_{R}$, obtained by fitting the
auto-correlation function of squirmer orientation

$<\textbf{e}(t)\cdot\textbf{e}(t+\Delta t)>=e^{-D_{R}\Delta t},$ (14)

where $<\dots>$ denotes the average over all squirmers and all times in the
simulated system. Here, the exponential decay applies for small $\Delta
t\lesssim 1/D_{R}$. Figure 9 shows the dependence of $D_{R}$ on $\phi$,
$\beta$, and $\lambda$. In most cases, $D_{R}/D_{R}^{0}>1$, indicating that
the rotational dynamics of squirmers within a suspension is enhanced in
comparison with the rotational diffusion $D_{R}^{0}$ of a spheroid in an
unconfined system. This is due to hydrodynamic interactions between squirmers,
and their collisions during motion. Only for the case of $\beta=0$,
$\phi=0.56$, and $\tilde{\lambda}=0$, rotational diffusion is reduced,
$D_{R}<D_{R}^{0}$, which is likely due to the crowding of squirmers at large
volume fractions [compare Fig. 2(a)].

Figure 9(a) ($\tilde{\lambda}=0$) shows that the enhancement of rotational
diffusion coefficient for pushers is significantly larger than for pullers and
neutral squirmers. This occurs due to repulsive hydrodynamic interactions
between pushers [45, 46], which lead to the enhanced scattering between
squirmers. The dependence of $D_{R}$ on $\phi$ first exhibits an increase,
followed by a maximum and subsequent decrease. With increasing $\phi$, the
collision rate between squirmers increases, resulting in the enhancement of
$D_{R}$. However, at large $\phi$, squirmer clustering and the transition to
the MIPS phase take place [compare Fig. 2(a)], so that the effective
rotational diffusion is significantly slowed down.

Noteworthy is the difference in $D_{R}$ for $\tilde{\lambda}=0$ and
$\tilde{\lambda}=133.5$ in Fig. 9(a) and (b), respectively, where the
enhancement for squirmers with rotlet dipole is at least a factor three larger
than that without. Interestingly, squirmers with $\tilde{\lambda}=133.5$
switch between the two layers (i.e. migrate from one wall to the other) much
more frequently than squirmers with $\tilde{\lambda}=0$, which will be
discussed in Section III.2.4. The ability of interchange between the two
layers increases significantly the effective rotational diffusion. The $D_{R}$
curves for $\tilde{\lambda}>0$ are similar for different swimming modes
$\beta$, in agreement with the structural characteristics for suspensions of
squirmers with the rotlet dipole, as discussed in Sec. III.1 above.

Figure 10: Mean-squared displacement of squirmers for $\phi=0.18$ (dashed
lines) and $\phi=0.56$ (solid lines) (a) without the rotlet dipole
$\tilde{\lambda}=0$ and (b) with the rotlet dipole $\tilde{\lambda}=133.5$.
Curves for pushers are in red, for neutral squirmers in blue, and for pullers
in green. Black dotted lines indicate power-laws $\sim t^{2}$ in the ballistic
regime and $\sim t$ in the diffusive regime.

####  III.2.2 Mean-square displacement

The mean-square displacement of squirmers for various conditions are shown in
Fig. 10. The motion is ballistic at short times, with a quadratic time-
dependence, and diffusive at long times, with $<(r(t)-r(0))>^{2}\sim t$, as
expected. For $\phi=0.18$ and $\tilde{\lambda}=0$, pullers yield the lowest
effective diffusivity, followed by neutral squirmers and pushers, because the
suspension of pullers at these conditions is already in the MIPS phase (see
Table 1). The effective diffusivity of squirmers at $\phi=0.56$ is lower than
at $\phi=0.18$ due to crowding and the formation of large clusters.

The transition from the ballistic to the diffusive regime for
$\tilde{\lambda}=0$ in Fig. 10(a) occurs around $tD_{R}^{0}\simeq 0.2-0.5$.
Theoretically, this transition should take place around $tD_{R}\simeq 1$ [78],
which is consistent with the values of $D_{R}/D_{R}^{0}\simeq 2$ in Fig. 9(a).
Note that the ballistic-to-diffusive transition for squirmers with rotlet
dipole occurs at $tD_{R}^{0}\gtrsim 0.1$, which is also consistent with the
values of $D_{R}/D_{R}^{0}\simeq 10$ in Fig. 9(b). Furthermore, mean-square-
displacement curves for $\tilde{\lambda}=133.5$ in Fig. 10(b) are similar to
each other for various $\beta$ at a fixed volume fraction, confirming once
more that the rotlet dipole nearly offsets the effects of different swimming
modes.

Figure 11: Average squirmer speed $\bar{v}$ as a function of $\phi$ for (a)
$\tilde{\lambda}=0$ and (b) $\tilde{\lambda}=133.5$. The speed is normalized
by the swimming speed $U_{0}$ of a squirmer in an unconfined system. Data for
pushers (red), for neutral squirmers (blue), and for pullers (green).

####  III.2.3 Average squirmer speed

To further characterize the mobility of squirmers for different conditions, we
also study the average squirmer speed $\bar{v}$. For a single squirmer $i$,
the instantaneous swim speed is calculated as

$v_{i}(t)=|\mathbf{r}_{i}(t+\Delta t)-\mathbf{r}_{i}(t)|/\Delta t,$ (15)

with time interval $\Delta t=0.022/D_{R}$ (significantly smaller than
$1/D_{R}$), during which a free squirmer moves a distance $b_{z}$. $\bar{v}$
is then obtained by an ensemble and time average as a function of the volume
fraction $\phi$.

As shown in Fig. 11, pushers are the most mobile, displaying the largest
average speed $\bar{v}$ among all squirmer types. For the case without rotlet
dipole in Fig. 11(a), $\bar{v}$ decreases with increasing $\phi$ due to
crowding at large volume fractions and the transition to MIPS or clustering
regime. For all studied $\beta$ values, the average speed of squirmers is
below $U_{0}$, which is the speed of a single squirmer in an unconfined
system.

The average speeds of squirmers with rotlet dipole in Fig. 11(b) are larger
than those with $\tilde{\lambda}=0$. Furthermore, the decay in $\bar{v}$ for
$\tilde{\lambda}=133.5$ with increasing $\phi$ is slower than for
$\tilde{\lambda}=0$. Interestingly, for $\phi\lesssim 0.4$, the average speed
of squirmers with rotlet dipole is slightly larger than $U_{0}$, indicating
motility enhancement due to interactions with the walls. An increased mobility
near walls was reported previously for pushers without rotlet dipole [46, 79,
76]. The results in Fig. 11 confirm that the presence of rotlet dipole
enhances the mobility of squirmers and delays the transition to MIPS as a
function of $\phi$.

Figure 12: Fraction $P_{m}$ of squirmers that switch between the two layers at
least once during a simulation as a function of $\phi$ for
$\tilde{\lambda}=0$. Curves for pushers (red), for neutral squirmers (blue),
and for pullers (green). The black dotted line marks $P_{m}=1$, which is the
fraction of squirmers that migrate between the two layers for a non-zero
rotlet dipole with $\tilde{\lambda}=133.5$, independently of $\phi$ and
$\beta$.

####  III.2.4 Migration between the two layers

Due to the dominance of two-layered squirmer structures for nearly all
conditions (see Fig. 4), an interesting issue is the switching freqnuncy of
individual squirmers between the layers. Figure 12 presents the fraction
$P_{m}$ of squirmers that switch between the two layers at least once during
the simulation time of $1.5\tau_{R}$ for different $\phi$ without rotlet
dipole. For pushers, $P_{m}$ is close to unity at small $\phi$, and decreases
to $P_{m}\approx 0.6$ at $\phi=0.56$. Previous studies [80, 45, 79, 57] have
shown that the residence of pushers near a wall is unstable due to local
hydrodynamic interactions which orient them slightly away from the wall. The
drop in $P_{m}$ at large $\phi$ for pushers is due to significant crowding,
which limits the migration between the two layers. For neutral squirmers,
$P_{m}$ reduces from unity to $P_{m}\approx 0.1$ as the volume fraction is
increased, indicating that their residence at the walls is more long-lived
than for pushers.

For pullers, notably a relatively low fraction migrates between the two
layers. This is consistent with the preferred orientation toward a wall,
compare Sec. III.1.3, ass well as the formation of stable flower-like
structures, as shown in Sec. III.1.4. Thus, switching fo pullers between the
two layers by pullers is facilitated by squirmer collisions and does not occur
spontaneously.

Finally, for squirmers with rotlet dipole $\tilde{\lambda}=133.5$, $P_{m}=1$
for all simulated $\beta$ and $\phi$ values. Thus, a rotlet dipole
significantly destabilizes the residence of squirmers near the walls, and
thereby enhances squirmer motility, as shown in Figs. 9 and 11.

##  IV Discussion and conclusions

In this study, we have performed mesoscale hydrodynamic simulations of a
confined system of squirmers in a thin fluid film between two parallel walls,
in order to better understand the intricate interplay between hydrodynamic and
steric interactions and their impact on the collective behavior of squirmers.
In contrast, previous studies that considered hydrodynamic interactions
between swimmers have mainly focused on quasi-2D systems (i.e., a monolayer of
squirmers) [38, 31, 46, 47, 39, 48, 49] or 3D systems with periodic boundary
conditions [50, 51, 52, 53, 37]. The considered film thickness is large enough
to allow the formation of two-layered structures. This situation is of
interest when the formation of a biofilm, which starts from a monolayer of
bacteria, proceeds toward the development of further layers. In particular, an
interesting question is whether bacteria can spontaneously leave the surface
layer.

By systematically varying relevant parameters, such as the volume fraction of
squirmers, their active stress, and their rotlet dipole strength, we show the
emergence of different structures and phases (Fig. 2). As expected, at low
volume fractions of squirmers, the system is in a gas-like phase. As $\phi$ is
increased, a swarming state is observed, with the formation of mobile clusters
with a wide range of sizes. The elongated shape of squirmers is essential for
the emergence of these swarming clusters, as it provides an alignment
interaction, which has been observed previously for self-propelled rods,
spheroids, and semi-flexible filaments [22, 38, 81]. Swarming clusters form at
intermediate values of $\phi$ for pushers and squirmers with the rotlet
dipole. The absence of the swarming state for pullers with $\tilde{\lambda}=0$
suggests that contractile hydrodynamic interactions between swimmers suppress
the formation of swarming clusters.

At large enough volume fractions, a single large cluster of squirmers emerges,
surrounded by a few mobile swimmers (see Fig. 2(a) at $\phi=0.44$ or $0.56$),
indicating formation of MIPS. Pullers exhibit MIPS already at low volume
fraction, $\phi\approx 0.18$, while neutral squirmers and pushers require
larger volume fractions, in agreement with previous studies [38, 46]. Our
simulations also suggest that hydrodynamic interactions suppress motility-
induced clustering and phase separation, in agreement with previous reports
[38, 40] for particles with aspect ratios not far from unity. Note that the
suspension of pushers at $\phi=0.56$ does not show a transition to the state
of active turbulence, because the strength of the induced force dipole is
likely too weak [29, 25].

An interesting aspect of our investigation is the interplay of simultaneous
hydrodynamic and steric interaction both between squirmers and of squirmers
with the confining walls, and its effect on the collective behavior. In
particular, pullers favor a nearly perpendicular orientation with the walls,
resulting in the formation of flower-like clusters [Fig. 2(b)] at low volume
frations, where pullers first slide along the surface, and upon collisions
form a structure with a single puller oriented nearly perpendicular to the
wall surrounded by a few petal pullers. This leads to the nucleation of
clusters for pullers already at low $\phi$, promoting MIPS at low $\phi$.
Pushers and neutral squirmers favor an orientation parallel to the walls,
which does not significantly hinder their mobility. Therefore, pushers can
frequently switch between the layers of a two-layer structure, substantially
delaying the formation of large clusters with increasing $\phi$. Thus, a
colony of pusher-like swimmers (e.g., E. coli) should be able to spontaneously
switch from an initial monolayer structure to a multi-layered configuration.
Recent experiments [74] suggest that E. coli bacteria also employ tumbling
(i.e., active turning due to the rotation reversal of one of its flagella) to
leave a surface, indicating that the aspect ratio of a swimmer is very
important for its interaction with the wall. In our simulations,
$b_{z}/b_{x}=2$, which is somewhat smaller than that for E. coli.

For all studied cases, there is a clear preference for the two-layered
structure. Only pullers at $\phi\lesssim 0.25$ show partial preference for the
perpendicular-to-the-wall orientation due to their hydrodynamic interactions
with the walls. Pushers and neutral squirmers display mostly parallel-to-wall
alignment.

Our comprehensive simulations reveal that the the effect of a rotlet dipole is
to offset the effect of active stress, characterized by $\beta$, to a large
extent. Even at relatively low volume fractions, the presence of a rotlet
dipole of dimensionless strength $\tilde{\lambda}\neq 0$ nearly eliminates
differences in the behavior of pullers, neutral squirmers, and pushers. With
the rotlet dipole, squirmers assume a parallel-to-the-wall orientation and
therefore, remain mobile up to moderate volume fractions, frequently switching
between the two layers. This enhancement of layer switching can be attributed
to the rotational flow field of rotors, which implies an induced motion around
their center of mass. Furthermore, the rotlet dipole leads to a circular
trajectories of squirmers near the walls [82, 83], which further contributes
to the mobility of squirmers within the confinement. In our study, a single
squirmer with $\tilde{\lambda}=133.5$ circles at the wall with a radius of
approximately $2b_{z}$. However, the circling motion near a wall is difficult
to detect at higher volume fractions due to frequent squimer-squirmer
scattering. As a result, the effective rotational diffusivity of squirmers
with rotlet dipole is approximately one order of magnitude larger than that
for squirmers with $\tilde{\lambda}=0$. Thus, swimmers with rotlet dipole are
expected to spontaneously initiate multi-layered structures within a biofilm.

Finally, we would like to discuss some limitations of our study. Despite the
fact that the squirmer model can produce various swimming modes, the near-
field flow of real microswimmers, like bacteria, is unavoidably species-
specific, arising from body shape, the geometry and dynamics of propelling
flagella, and the ability to navigate (e.g., E. coli tumbling). Furthermore,
steric interactions between microswimmers including their flagella are likely
different from interactions between idealized spheroidal shapes. However,
modeling of swimmers with explicit appendages [83] is computationally
expensive, significantly limiting the number of simulated swimmers in a study.
Therefore, the direct comparison of our simulation results with experiments
requires some calibration, like of the aspect ratio of squirmers, and the
strengths of active stress and rotlet dipole. Nevertheless, we expect that our
simulations provide a good qualitative characterization of swimmer suspensions
under confinement. Future work should also consider thicker films to bridge
these results with those from three-dimensional periodic systems [50, 51, 52,
53, 37].

## Appendix A: Vesicle-like squirmer model

Each squirmer surface is represented by $N_{p}=1024$ particles connected by
springs into a triangulated network. Shear elasticity of the squirmer membrane
is supplied by the spring potential [55]

$U_{bond}=\frac{k_{B}Tl_{m}}{4p}\frac{3x^{2}-2x^{3}}{1-x}+\frac{k_{p}}{l},$
(16)

where $x=l/l_{m}\in(0,1)$, $l$ is the spring length, $l_{m}$ is the maximum
spring extension, $p$ is the persistence length, $k_{p}$ is the force
coefficient, and $k_{B}T$ is the energy unit defined by the temperature $T$ in
the simulated system. The curvature elasticity that provides bending
resistance is implemented through the discretisation of the Helfrich bending
energy [84] as

$U_{bend}=\frac{\kappa}{2}\sum_{i}\sigma_{i}(H_{i}-H_{0}^{i})^{2},$ (17)

where $\kappa$ is the bending rigidity, $\sigma_{i}$ is the area corresponding
to vertex $i$ in the membrane triangulation, $H_{i}$ is the mean curvature at
vertex $i$, and $H_{0}^{i}$ is the spontaneous curvature at vertex $i$. The
mean curvature is discretized as $H_{i}={\bf
n}_{i}\cdot\sum_{j(i)}\sigma_{ij}{\bf r}_{ij}/(\sigma_{i}r_{ij})$, where ${\bf
n}_{i}$ is the unit normal at the vertex $i$,
$\sigma_{i}=\sum_{j(i)}\sigma_{ij}r_{ij}/4$, $j(i)$ spans all vertices linked
to vertex $i$, and $\sigma_{ij}=r_{ij}(\cot\theta_{1}+\cot\theta_{2})/2$ is
the bond length in the dual lattice with $\theta_{1}$ and $\theta_{2}$ being
the angles at the two vertices opposite to the edge $ij$ in the dihedral. The
spontaneous curvature $H_{0}^{i}$ is set locally after the triangulation of
the spheroidal surface of the squirmer.

The area and volume conservation constraints are represented by the potential
[55]

$U_{a+v}=\frac{k_{a}(A-A_{0}^{tot})^{2}}{2A_{0}^{tot}}+\sum_{m\in
1...N_{t}}\frac{k_{d}(A_{m}-A_{0}^{m})^{2}}{2A_{0}^{m}}+\frac{k_{v}(V-V_{0}^{tot})^{2}}{2V_{0}^{tot}},$
(18)

where $k_{a}$, $k_{d}$, and $k_{v}$ are the coefficients of global area, local
area and volume conservation constraints, respectively. $A$ and $V$ are the
instantaneous area and volume of the enclosed membrane, $A_{0}^{tot}$ and
$V_{0}^{tot}$ are the targeted global area and volume which are defined by the
spheroidal shape. $A_{m}$ is the area of the $m$-th triangle (or face) within
the triangulation, while $A_{m}^{0}$ is the targeted value. $N_{t}$ is the
number of triangles within the triangulated surface.

The membrane-like representation of squirmers has the shear modulus
$\mu_{0}=1.44\times 10^{5}k_{B}T/b_{z}^{2}$, the bending modulus
$\kappa=250k_{B}T$, the area-constraint coefficients $k_{d}=1.6\times
10^{3}k_{B}T/b_{z}^{2}$ and $k_{a}=8\times 10^{3}k_{B}T/b_{z}^{2}$, the
volume-constraint coefficient $k_{v}=3.2\times 10^{4}k_{B}T/b_{z}^{3}$, the
total area $A_{0}^{tot}=5.36b_{z}^{2}$, and the total volume
$V_{0}^{tot}=1.05b_{z}^{3}$. In all simulations, $k_{B}T=1$ and $b_{z}=4$.

## Appendix B: Modeling fluid flow

To model fluid flow, we employ the dissipative particle dynamics (DPD) method
[62, 63], which is a mesoscopic hydrodynamics simulation technique. DPD is a
particle-based Lagrangian method, where each particle represents a small fluid
volume. DPD particles $i$ and $j$ interact through three types (conservative,
dissipative, and random) of pairwise forces given by

$\displaystyle\textbf{F}^{C}(r_{ij})$
$\displaystyle=aW^{C}(r_{ij})\hat{\textbf{r}}_{ij},$ (19)
$\displaystyle\textbf{F}^{D}(r_{ij})$ $\displaystyle=-\gamma
W^{D}(r_{ij})(\hat{\textbf{r}}_{ij}\cdot\textbf{v}_{ij})\hat{\textbf{r}}_{ij},$
(20) $\displaystyle\textbf{F}^{R}(r_{ij})$ $\displaystyle=\sigma
W^{R}(r_{ij})\xi_{ij}\hat{\textbf{r}}_{ij}/\sqrt{\Delta t},$ (21)

where $a$, $\gamma$, and $\sigma$ are the force amplitudes,
$\textbf{r}_{ij}=\textbf{r}_{i}-\textbf{r}_{j}$ is the relative position
vector, $r_{ij}=|\textbf{r}_{ij}|$,
$\hat{\textbf{r}}_{ij}=\textbf{r}_{ij}/r_{ij}$, and
$\textbf{v}_{ij}=\textbf{v}_{i}-\textbf{v}_{j}$ is the velocity difference.
$\xi_{ij}=\xi_{ji}$ is a symmetric Gaussian random variable with zero mean and
unit variance, and $\Delta t$ is the time step. All forces act within a cutoff
radius $r_{c}$ and vanish beyond it. The conservative force controlls fluid
compressibility, while the dissipative and random forces form a thermostat, so
that the DPD fluid has an isotropic temperature $T$. Thus, $\textbf{F}^{D}$
and $\textbf{F}^{R}$ are related through the fluctuation-dissipation theorem
[63] as

$\sigma^{2}=2\gamma k_{B}T,\quad\quad W^{D}(r_{ij})=[W^{R}(r_{ij})]^{2}.$ (22)

The weight functions are defined as

$W(r_{ij})=W^{R}(r_{ij})=\begin{cases}(1-r_{ij}/r_{c})^{s},&r_{ij}<r_{c},\\\
0,&r_{ij}\geq r_{c},\end{cases}$ (23)

with an exponent $s$. For the conservative force, $W^{C}(r_{ij})=W(r_{ij})$
with $s=1$.

Time evolution of each DPD particle follows the Newton’s second law

$\frac{d\textbf{r}_{i}}{dt}=\textbf{v}_{i},\quad\quad
m_{i}\frac{d\textbf{v}_{i}}{dt}=\sum_{j\neq
i}\left(\textbf{F}^{C}(r_{ij})+\textbf{F}^{D}(r_{ij})+\textbf{F}^{R}(r_{ij})\right),$
(24)

where $m_{i}$ is the mass of particle $i$. Time integration is performed using
the velocity-Verlet algorithm.

DPD parameters for the interactions between fluid particles and between fluid
and wall particles are $a=200k_{B}T/b_{z}$, $\gamma=80\sqrt{mk_{B}T}/b_{z}$
($m=1$ in all simulations), $s=0.15$, and $r_{c}=b_{z}/4$. Furthermore,
friction coupling of squirmers to fluid flow is performed using DPD parameters
$a=0$, $\gamma=100\sqrt{mk_{B}T}/b_{z}$, $s=0.1$, and $r_{c}=0.175b_{z}$. The
time step for integration is $1.25\times 10^{-3}b_{z}\sqrt{m/k_{B}T}$.

## Author Contributions

G.G. and D.A.F. designed the research project. B.W.-Z. performed the
simulations and analysed the obtained data. All authors participated in the
discussions and writing of the manuscript.

## Conflicts of interest

There are no conflicts to declare.

## Acknowledgements

We thank Roland G. Winkler for numerous discussions related to the study. We
acknowledge funding from the ETN-PHYMOT ”Physics of microbial motility” within
the European Union’s Horizon 2020 research and innovation programme under the
Marie Skłodowska-Curie grant agreement No 955910. The authors also gratefully
appreciate computing time on the supercomputer JURECA [85] at
Forschungszentrum Jülich under grant no. actsys.

## References

  * Elgeti  _et al._ [2015] J. Elgeti, R. G. Winkler, and G. Gompper, Physics of microswimmers \- single particle motion and collective behavior: a review, Rep. Prog. Phys. 78, 056601 (2015).
  * Bechinger  _et al._ [2016] C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe, and G. Volpe, Active particles in complex and crowded environments, Rev. Mod. Phys. 88, 045006 (2016).
  * Gompper  _et al._ [2020] G. Gompper, R. G. Winkler, T. Speck, A. Solon, C. Nardini, F. Peruani, H. Löwen, R. Golestanian, U. B. Kaupp, L. Alvarez, T. Kiorboe, E. Lauga, W. C. K. Poon, A. DeSimone, S. Muiños-Landin, A. Fischer, N. A. Söker, F. Cichos, R. Kapral, P. Gaspard, M. Ripoll, F. Sagues, A. Doostmohammadi, J. M. Yeomans, I. S. Aranson, C. Bechinger, H. Stark, C. K. Hemelrijk, F. J. Nedelec, T. Sarkar, T. Aryaksama, M. Lacroix, G. Duclos, V. Yashunsky, P. Silberzan, M. Arroyo, and S. Kale, The 2020 motile active matter roadmap, J. Phys.: Condens. Matter 32, 193001 (2020).
  * Hall-Stoodley  _et al._ [2004] L. Hall-Stoodley, J. W. Costerton, and P. Stoodley, Bacterial biofilms: from the natural environment to infectious diseases, Nat. Rev. Microbiol. 2, 95 (2004).
  * Verstraeten  _et al._ [2008] N. Verstraeten, K. Braeken, B. Debkumari, M. Fauvart, J. Fransaer, J. Vermant, and J. Michiels, Living on a surface: swarming and biofilm formation, Trends Microbiol. 16, 496 (2008).
  * Marsh [2006] P. D. Marsh, Dental plaque as a biofilm and a microbial community–implications for health and disease, BMC Oral Health 6, S14 (2006).
  * Attinger and Wolcott [2012] C. Attinger and R. Wolcott, Clinically addressing biofilm in chronic wounds, Adv. Wound Care 1, 127 (2012).
  * Mazza [2016] M. G. Mazza, The physics of biofilms-an introduction, J. Phys. D: Appl. Phys. 49, 203001 (2016).
  * Palacci  _et al._ [2013] J. Palacci, S. Sacanna, A. P. Steinberg, D. J. Pine, and P. M. Chaikin, Living crystals of light-activated colloidal surfers, Science 339, 936 (2013).
  * Buttinoni  _et al._ [2013] I. Buttinoni, J. Bialké, F. Kümmel, H. Löwen, C. Bechinger, and T. Speck, Dynamical clustering and phase separation in suspensions of self-propelled colloidal particles, Phys. Rev. Lett. 110, 238301 (2013).
  * Bäuerle  _et al._ [2020] T. Bäuerle, R. C. Löffler, and C. Bechinger, Formation of stable and responsive collective states in suspensions of active colloids, Nat. Comm. 11, 2547 (2020).
  * Krüger  _et al._ [2016] C. Krüger, C. Bahr, S. Herminghaus, and C. C. Maass, Dimensionality matters in the collective behaviour of active emulsions, Eur. Phys. J. E 39, 64 (2016).
  * Thutupalli  _et al._ [2018] S. Thutupalli, D. Geyer, R. Singh, R. Adhikari, and H. A. Stone, Flow-induced phase separation of active particles is controlled by boundary conditions, Proc. Natl. Acad. Sci. USA 115, 5403 (2018).
  * Bricard  _et al._ [2013] A. Bricard, J.-B. Caussin, N. Desreumaux, O. Dauchot, and D. Bartolo, Emergence of macroscopic directed motion in populations of motile colloids, Nature 503, 95 (2013).
  * Fily and Marchetti [2012] Y. Fily and M. C. Marchetti, Athermal phase separation of self-propelled particles with no alignment, Phys. Rev. Lett. 108, 235702 (2012).
  * Bialké  _et al._ [2013] J. Bialké, H. Löwen, and T. Speck, Microscopic theory for the phase separation of self-propelled repulsive disks, Europhys. Lett. 103, 30008 (2013).
  * Redner  _et al._ [2013] G. S. Redner, M. F. Hagan, and A. Baskaran, Structure and dynamics of a phase-separating active colloidal fluid, Phys. Rev. Lett. 110, 055701 (2013).
  * Cates and Tailleur [2015] M. E. Cates and J. Tailleur, Motility-induced phase separation, Annu. Rev. Condens. Matter Phys. 6, 219 (2015).
  * Digregorio  _et al._ [2018] P. Digregorio, D. Levis, A. Suma, L. F. Cugliandolo, G. Gonnella, and I. Pagonabarraga, Full phase diagram of active Brownian disks: from melting to motility-induced phase separation, Phys. Rev. Lett. 121, 098003 (2018).
  * Liu  _et al._ [2019] G. Liu, A. Patch, F. Bahar, D. Yllanes, R. D. Welch, M. C. Marchetti, S. Thutupalli, and J. W. Shaevitz, Self-driven phase transitions drive Myxococcus xanthus fruiting body formation, Phys. Rev. Lett. 122, 248102 (2019).
  * Be’er  _et al._ [2020] A. Be’er, B. Ilkanaiv, R. Gross, D. B. Kearns, S. Heidenreich, M. Bär, and G. Ariel, A phase diagram for bacterial swarming, Commun. Phys. 3, 66 (2020).
  * Peruani  _et al._ [2006] F. Peruani, A. Deutsch, and M. Bär, Nonequilibrium clustering of self-propelled rods, Phys. Rev. E 74, 030904 (2006).
  * Vicsek and Zafeiris [2012] T. Vicsek and A. Zafeiris, Collective motion, Phys. Rep. 517, 71 (2012).
  * Wioland  _et al._ [2013] H. Wioland, F. G. Woodhouse, J. Dunkel, J. O. Kessler, and R. E. Goldstein, Confinement stabilizes a bacterial suspension into a spiral vortex, Phys. Rev. Lett. 110, 268102 (2013).
  * Martinez  _et al._ [2020] V. A. Martinez, E. Clément, J. Arlt, C. Douarche, A. Dawson, J. Schwarz-Linek, A. K. Creppy, V. Skultéty, A. N. Morozov, H. Auradou, and W. C. K. Poon, A combined rheometry and imaging study of viscosity reduction in bacterial suspensions, Proc. Natl. Acad. Sci. USA 117, 2326 (2020).
  * Guo  _et al._ [2018] S. Guo, D. Samanta, Y. Peng, X. Xu, and X. Cheng, Symmetric shear banding and swarming vortices in bacterial superfluids, Proc. Natl. Acad. Sci. USA 115, 7212 (2018).
  * Kraichnan and Montgomery [1980] R. H. Kraichnan and D. Montgomery, Two-dimensional turbulence, Rep. Prog. Phys. 43, 547 (1980).
  * Wensink  _et al._ [2012] H. H. Wensink, J. Dunkel, S. Heidenreich, K. Drescher, R. E. Goldstein, H. Löwen, and J. M. Yeomans, Meso-scale turbulence in living fluids, Proc. Natl. Acad. Sci. USA 109, 14308 (2012).
  * Stenhammar  _et al._ [2017] J. Stenhammar, C. Nardini, R. W. Nash, D. Marenduzzo, and A. Morozov, Role of correlations in the collective behavior of microswimmer suspensions, Phys. Rev. Lett 119, 028005 (2017).
  * Bárdfalvy  _et al._ [2019] D. Bárdfalvy, H. Nordanger, C. Nardini, A. Morozov, and J. Stenhammar, Particle-resolved lattice Boltzmann simulations of 3-dimensional active turbulence, Soft Matter 15, 7747 (2019).
  * Qi  _et al._ [2022] K. Qi, E. Westphal, G. Gompper, and R. G. Winkler, Emergence of active turbulence in microswimmer suspensions due to active hydrodynamic stress and volume exclusion, Commun. Phys. 5, 49 (2022).
  * Lushi  _et al._ [2014] E. Lushi, H. Wioland, and R. E. Goldstein, Fluid flows created by swimming bacteria drive self-organization in confined suspensions, Proc. Natl. Acad. Sci. USA 111, 9733 (2014).
  * Theillard  _et al._ [2017] M. Theillard, R. Alonso-Matillaa, and D. Saintillan, Geometric control of active collective motion, Soft Matter 13, 363 (2017).
  * Suma  _et al._ [2014] A. Suma, G. Gonnella, D. Marenduzzo, and E. Orlandini, Motility-induced phase separation in an active dumbbell fluid, Europhys. Lett. 108, 56004 (2014).
  * Bär  _et al._ [2020] M. Bär, R. Großmann, S. Heidenreich, and F. Peruani, Self-propelled rods: insights and perspectives for active matter, Annu. Rev. Condens. Matter Phys. 11, 441 (2020).
  * Moran  _et al._ [2022] S. E. Moran, I. R. Bruss, P. W. A. Schönhöfer, and S. C. Glotzer, Particle anisotropy tunes emergent behavior in active colloidal systems, Soft Matter 18, 1044 (2022).
  * Alarcón and Pagonabarraga [2013] F. Alarcón and I. Pagonabarraga, Spontaneous aggregation and global polar ordering in squirmer suspensions, J. Mol. Liquids 185, 56 (2013).
  * Theers  _et al._ [2018] M. Theers, E. Westphal, K. Qi, R. G. Winkler, and G. Gompper, Clustering of microswimmers: interplay of shape and hydrodynamics, Soft Matter 14, 8590 (2018).
  * Zöttl and Stark [2014] A. Zöttl and H. Stark, Hydrodynamics determines collective motion and phase behavior of active colloids in quasi-two-dimensional confinement, Phys. Rev. Lett. 112, 118101 (2014).
  * Matas-Navarro  _et al._ [2014] R. Matas-Navarro, R. Golestanian, T. B. Liverpool, and S. M. Fielding, Hydrodynamic suppression of phase separation in active suspensions, Phys. Rev. E 90, 032304 (2014).
  * Kuhr  _et al._ [2019] J.-T. Kuhr, F. Rühle, and H. Stark, Collective dynamics in a monolayer of squirmers confined to a boundary by gravity, Soft Matter 15, 5685 (2019).
  * Stenhammar  _et al._ [2014] J. Stenhammar, D. Marenduzzo, R. J. Allen, and M. E. Cates, Phase behaviour of active Brownian particles: the role of dimensionality, Soft Matter 10, 1489 (2014).
  * Lighthill [1952] M. J. Lighthill, On the squirming motion of nearly spherical deformable bodies through liquids at very small Reynolds numbers, Commun. Pure Appl. Math. 5, 109 (1952).
  * Blake [1971] J. R. Blake, A spherical envelope approach to ciliary propulsion, J. Fluid Mech. 46, 199 (1971).
  * Theers  _et al._ [2016] M. Theers, E. Westphal, G. Gompper, and R. G. Winkler, Modeling a spheroidal microswimmer and cooperative swimming in a narrow slit, Soft Matter 12, 7372 (2016).
  * Kyoya  _et al._ [2015] K. Kyoya, D. Matsunaga, Y. Imai, T. Omori, and T. Ishikawa, Shape matters: near-field fluid mechanics dominate the collective motions of ellipsoidal squirmers, Phys. Rev. E 92, 063027 (2015).
  * Alarcón  _et al._ [2017] F. Alarcón, C. Valeriani, and I. Pagonabarraga, Morphology of clusters of attractive dry and wet self-propelled spherical particle suspensions, Soft Matter 13, 814 (2017).
  * Yoshinaga and Liverpool [2017] N. Yoshinaga and T. B. Liverpool, Hydrodynamic interactions in dense active suspensions: from polar order to dynamical clusters, Phys. Rev. E 96, 020603 (2017).
  * Blaschke  _et al._ [2016] J. Blaschke, M. Maurer, K. Menon, A. Zöttl, and H. Stark, Phase separation and coexistence of hydrodynamically interacting microswimmers, Soft Matter 12, 9821 (2016).
  * Ishikawa and Pedley [2007] T. Ishikawa and T. J. Pedley, Diffusion of swimming model micro-organisms in a semi-dilute suspension, J. Fluid Mech. 588, 437 (2007).
  * Delmotte  _et al._ [2015] B. Delmotte, E. E. Keaveny, F. Plouraboué, and E. Climent, Large-scale simulation of steady and time-dependent active suspensions with the force-coupling method, J. Comp. Phys. 302, 524 (2015).
  * Ishikawa  _et al._ [2008] T. Ishikawa, J. T. Locsei, and T. J. Pedley, Development of coherent structures in concentrated suspensions of swimming model micro-organisms, J. Fluid Mech. 615, 401 (2008).
  * Evans  _et al._ [2011] A. A. Evans, T. Ishikawa, T. Yamaguchi, and E. Lauga, Orientational order in concentrated suspensions of spherical microswimmers, Phys. Fluids 23, 111702 (2011).
  * Darnton  _et al._ [2007] N. C. Darnton, L. Turner, S. Rojevsky, and H. C. Berg, On torque and tumbling in swimming escherichia coli, J. Bacteriol. 189, 1756 (2007).
  * Fedosov  _et al._ [2010] D. A. Fedosov, B. Caswell, and G. E. Karniadakis, A multiscale red blood cell model with accurate mechanics, rheology, and dynamics, Biophys. J. 98, 2215 (2010).
  * Fedosov  _et al._ [2014] D. A. Fedosov, H. Noguchi, and G. Gompper, Multiscale modeling of blood flow: from single cells to blood rheology, Biomech. Model. Mechanobiol. 13, 239 (2014).
  * Zöttl and Stark [2018] A. Zöttl and H. Stark, Simulating squirmers with multiparticle collision dynamics, Eur. Phys. J. E 41, 61 (2018).
  * Ishikawa  _et al._ [2006] T. Ishikawa, M. P. Simmonds, and T. J. Pedley, Hydrodynamic interaction of two swimming model micro-organisms, J. Fluid Mech. 568, 119 (2006).
  * Pagonabarraga and Llopis [2013] I. Pagonabarraga and I. Llopis, The structure and rheology of sheared model swimmer suspensions, Soft Matter 9, 7174 (2013).
  * Qi  _et al._ [2020] K. Qi, H. Annepu, G. Gompper, and R. G. Winkler, Rheotaxis of spheroidal squirmers in microchannel flow: interplay of shape, hydrodynamics, active stress, and thermal fluctuations, Phys. Rev. Res. 2, 033275 (2020).
  * Hu  _et al._ [2015a] J. Hu, M. Yang, G. Gompper, and R. G. Winkler, Modelling the mechanics and hydrodynamics of swimming E. coli, Soft Matter 11, 7867 (2015a).
  * Hoogerbrugge and Koelman [1992] P. J. Hoogerbrugge and J. M. V. A. Koelman, Simulating microscopic hydrodynamic phenomena with dissipative particle dynamics, Europhys. Lett. 19, 155 (1992).
  * Español and Warren [1995] P. Español and P. Warren, Statistical mechanics of dissipative particle dynamics, Europhys. Lett. 30, 191 (1995).
  * Solon  _et al._ [2015] A. P. Solon, J. Stenhammar, R. Wittkowski, M. Kardar, Y. Kafri, M. E. Cates, and J. Tailleur, Pressure and phase equilibria in interacting active Brownian spheres, Phys. Rev. Lett. 114, 198301 (2015).
  * Drescher  _et al._ [2011] K. Drescher, J. Dunkel, L. H. Cisneros, S. Ganguly, and R. E. Goldstein, Fluid dynamics and noise in bacterial cell-cell and cell-surface scattering, Proc. Natl. Acad. Sci. USA 108, 10940 (2011).
  * Lauga and Powers [2009] E. Lauga and T. R. Powers, The hydrodynamics of swimming microorganisms, Rep. Prog. Phys. 72, 096601 (2009).
  * Bialké  _et al._ [2012] J. Bialké, T. Speck, and H. Löwen, Crystallization in a dense suspension of self-propelled particles, Phys. Rev. Lett. 108, 168301 (2012).
  * Levis and Berthier [2014] D. Levis and L. Berthier, Clustering and heterogeneous dynamics in a kinetic Monte Carlo model of self-propelled hard disks, Phys. Rev. E 89, 062301 (2014).
  * Li and Tang [2009] G. Li and J. X. Tang, Accumulation of microswimmers near a surface mediated by collision and rotational Brownian motion, Phys. Rev. Lett. 103, 078101 (2009).
  * Volpe  _et al._ [2011] G. Volpe, I. Buttinoni, D. Vogt, H.-J. Kümmerer, and C. Bechinger, Microswimmers in patterned environments, Soft Matter 7, 8810 (2011).
  * Elgeti and Gompper [2013] J. Elgeti and G. Gompper, Wall accumulation of self-propelled spheres, Europhys. Lett. 101, 48003 (2013).
  * Berke  _et al._ [2008] A. P. Berke, L. Turner, H. C. Berg, and E. Lauga, Hydrodynamic attraction of swimming microorganisms by surfaces, Phys. Rev. Lett. 101, 038102 (2008).
  * Frymier  _et al._ [1995] P. D. Frymier, R. M. Ford, H. C. Berg, and P. T. Cummings, Three-dimensional tracking of motile bacteria near a solid planar surface, Proc. Natl. Acad. Sci. USA 92, 6195 (1995).
  * Junot  _et al._ [2022] G. Junot, T. Darnige, A. Lindner, V. A. Martinez, J. Arlt, A. Dawson, W. C. K. Poon, H. Auradou, and E. Clément, Run-to-tumble variability controls the surface residence times of E. coli bacteria, Phys. Rev. Lett. 128, 248101 (2022).
  * Vigeant  _et al._ [2002] M. A.-S. Vigeant, R. M. Ford, M. Wagner, and L. K. Tamm, Reversible and irreversible adhesion of motile escherichia coli cells analyzed by total internal reflection aqueous fluorescence microscopy, Appl. Environ. Microbiol. 68, 2794 (2002).
  * Spagnolie and Lauga [2012] S. E. Spagnolie and E. Lauga, Hydrodynamics of self-propulsion near a boundary: predictions and accuracy of far-field approximations, J. Fluid Mech. 700, 105 (2012).
  * Götze and Gompper [2011] I. O. Götze and G. Gompper, Flow generation by rotating colloids in planar microchannels, Europhys. Lett. 92, 64003 (2011).
  * Romanczuk  _et al._ [2012] P. Romanczuk, M. Bär, W. Ebeling, B. Lindner, and L. Schimansky-Geier, Active brownian particles: From individual to collective stochastic dynamics, Eur. Phys. J. Special Topics 202, 1 (2012).
  * Lintuvuori  _et al._ [2016] J. S. Lintuvuori, A. T. Brown, K. Stratford, and D. Marenduzzo, Hydrodynamic oscillations and variable swimming speed in squirmers close to repulsive walls, Soft Matter 12, 7959 (2016).
  * Ishimoto and Gaffney [2013] K. Ishimoto and E. A. Gaffney, Squirmer dynamics near a boundary, Phys. Rev. E 88, 062702 (2013).
  * Duman  _et al._ [2018] Ö. Duman, R. E. Isele-Holder, J. Elgeti, and G. Gompper, Collective dynamics of self-propelled semiflexible filaments, Soft Matter 14, 4483 (2018).
  * Lauga  _et al._ [2006] E. Lauga, W. R. DiLuzio, G. M. Whitesides, and H. A. Stone, Swimming in circles: motion of bacteria near solid boundaries, Biophys. J. 90, 400 (2006).
  * Hu  _et al._ [2015b] J. Hu, A. Wysocki, R. G. Winkler, and G. Gompper, Physical sensing of surface properties by microswimmers – directing bacterial motion via wall slip, Sci. Rep. 5, 9586 (2015b).
  * Helfrich [1973] W. Helfrich, Elastic properties of lipid bilayers: theory and possible experiments, Z. Naturforsch. 28, 693 (1973).
  * Jülich Supercomputing Centre [2021] Jülich Supercomputing Centre, JURECA: Data Centric and Booster Modules implementing the Modular Supercomputing Architecture at Jülich Supercomputing Centre, J. Large-Scale Res. Facil. 7, A182 (2021).