Collective behavior of squirmers in thin films

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.

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 [31,37,38], hydrodynamic interactions [38][39][40], the presence of confinement [32,41], and dimensionality [12,42].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 spheroidallyshaped 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 microswim-mers.The squirmer model has been used to study collective behavior of microswimmers in quasi two-dimensional (2D) systems (i.e., a monolayer of swimmers) [31,38,39,41,[46][47][48][49] and three-dimensional (3D) settings with periodic boundary conditions (BCs) [37,[50][51][52][53].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 ϕ of particles [38,40].Pullers show the onset of MIPS at lower ϕ values in comparison to pushers.For spheroidal squirmers in quasi-2D systems [31,38], 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 ϕ values.For the 3D systems with periodic BCs [37,[50][51][52][53], 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 λ 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 ϕ values in comparison with pushers and neutral swimmers.At moderate volume fractions, all squirmers with rotlet dipole and pushers with λ = 0 show swarming behavior.Pullers prefer a nearly perpendicular orientation to the wall, which leads to the formation of flower-like structures at low ϕ and provides the nucleation for larger clusters as ϕ is increased.Pullers rarely switch between the two walls, while pushers do so frequently, indicating that a collection of pushers would quickly develop a multilayered 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

A. Squirmer model
We consider a spheroidal squirmer model, whose surface is described by 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]] where e z , e ζ , and e φ are unit vectors in Cartesian (x, y, z) or spheroidal (ζ, τ, φ) coordinates, which are related to each other as Here, c = b 2 z − b 2 x , and the spheroidal coordinates have the ranges [45,60].The coefficient β defines different swimming modes of the squirmer, including a pusher (β < 0), a neutral swimmer (β = 0), and a puller (β > 0).The second term in Eq. ( 2) defines a rotlet dipole with the strength λ, where r s = (x s , y s , z s ), r s = |r s |, rs = 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].

B. 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 bounceback 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 at the fluid-solid interface.Here, ϵ LJ is the potential strength, σ 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 σ 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 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 where rij = (r i − r j )/r ij , r ij = |r i − r j |, u i sq is the slip velocity at the position of squirmer particle i, while j corresponds to an outer-fluid DPD particle.Furthermore, the friction coefficient γ between fluid and squirmer particles is properly adjusted [55] to ensure the imposition of u sq at the squirmer surface.

C. Simulation setup and parameters
The simulation setup corresponds to a domain of dimensions L × L y × L, where L = 28b z and L y = 2.5b z (b z = 4 in simulations), see Fig.  [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 τ 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.5bz D R , in agreement with the theoretical prediction [45].Using these values, we can define a dimensionless Péclet number P e = U 0 /(2b z D R ) ≈ 28.4.Furthermore, the Reynolds number Re = 2b z U 0 mn f /η ≈ 0.03 is small enough to eliminate possible inertial effects.Different volume fractions ϕ = 4N sq πb x b y b z /(3L 2 L y ) ∈ {0.18, 0.35, 0.44, 0.56} are considered, corresponding to different numbers of squirmers N sq ∈ {228, 456, 576, 720}.The three values of active stress β ∈ {−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 λ = λ/(b 4 z D R ) ∈ {0, 133.5}.For λ = 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 λ ∈ {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τ R to reach a steady state, and afterwards various structural and dynamical characteristics of the system are measured during the time 1.5τ R .

D. Mapping to experimental systems
Squirmer parameters used in simulations can also be mapped to the properties of real swimmers.In the farfield approximation, the hydrodynamic field of a microswimmer is dominated by its force-dipole strength χ = f d l d /(8πη) [64][65][66], where f d and l d are the characteristic force and length of the dipole.The active stress β of a spheroidal squirmer can be expressed as [45] [54,61,65], resulting in β ≈ −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 ϕ and squirmer characteristics of active stress β and rotlet dipole strength λ.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 λ = 0, the layer correlation is very weak for small ϕ, but becomes very significant for large ϕ, where the clusters in the top and bottom layers are essentially in registry.With strong rotlet dipole, i.e. for λ = 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,2,31,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.
A. Structural properties

Cluster size distribution
The cluster-size distribution function N (n) is calculated as where same cluster when the nearest surface-to-surface distance d s between them satisfies d s /b z < 0.25.The average cluster size ⟨n⟩ is then [68]: The cluster-size distribution is used to distinguish the different collective phases [21,68] mentioned above.For the gas of small clusters, N (n) exhibits an exponential decay, as shown in Fig. 3(a) for ϕ = 0.18.For MIPS or large clusters, a bimodal distribution of N (n) is observed, with the second peak at large n signaling the presence of large clusters, see Fig. 3(b,c) for ϕ = 0.35 and ϕ = 0.56.Finally, swarming can be characterized by a power-law decay of N (n) with an exponential cutoff and without the presence of a distinct peak [31,47].This implies a slow decline of cluster sizes over a broad range, compared to the fast decay for a gas of small clusters.The described characteristics of different phases can be extracted by fitting N (n) with the function [47] where A, γ, and α are the fitting parameters.Table I 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 I.(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 ( λ = 0), while the insets show data from the corresponding simulations with λ = 133.5.
presents these parameters for different simulated conditions.
For ϕ = 0.18, the simulated systems exhibit a homogeneous gas-like phase, except for the case of pullers (β = 5) without rotlet dipole ( λ = 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 gaslike phase, compare also Fig. 2. For β = 5 at ϕ = 0.18, N (n) resembles a bimodal distribution, indicating the existence of large clusters (n ≈ 60), whose formation is primarily governed by the attractive hydrodynamic field around pullers.
Swarming-like behavior is observed in several cases at ϕ = 0.35, characterized by the relatively large values of α and moderate values of γ ∈ [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 γ is close to unity, and N (n) exhibits a bimodal distribution in Fig. 3(b), suggesting the presence of a MIPS phase.For even larger volume fraction ϕ = 0.56, Fig. 3(c) shows a prominent peak in N (n) at large n for all cases, which is the main characteristic for the MIPS phase.This is also confirmed in Table I through γ values close to unity.At large packing fractions, steric interactions between squirmers dominate due to crowding, so that the differences in N (n) for different swimming modes nearly disappear.
For squirmers without rotlet dipole ( λ = 0), the swimming mode affects their collective behavior when the volume fraction is small enough, i.e. ϕ ≲ 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 [31,38].However, when the rotlet dipole is activated ( λ = 133.5),differences in N (n) for various active stresses β nearly disappear, as displayed in the insets of Fig. 3.With rotlet dipole, we obtain the gas phase at ϕ = 0.18, swarming at ϕ = 0.35, and MIPS at ϕ = 0.56, independently of active stress -see also Fig. 2(c).Note that the systems with   71], on the other hand due to hydrodynamic attraction [65,66,72].Pushers (β = −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 [65,66,72].Furthermore, the elongated spheroidal body of the squirmers favors their orientation parallel to the walls [75,76].For pullers (β = 5), Fig. 4(a) show two small peaks at y/b z ≈ ±0.2 at small volume fraction ϕ = 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 [65,66,72], 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 ≈ ±0.5, such that their most probable ori-entation has about 40 degree angle with respect to the walls [see also Fig. 5(a)].p y for neutral swimmers at ϕ = 0.18 is close to that for pushers.At the higher volume fraction of ϕ = 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 ≈ ±0.6 become broader than those for ϕ = 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 A 1.

Squirmer orientation in the slit
The orientational distribution function p θ (θ) 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 θ = 0 and θ = 180, a parallel-to-the-wall orientation to θ = 90.Figure 5(a) shows that pushers and neutral swimmers with λ = 0 at ϕ = 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 θ ≈ 20 o & 160 o (small peaks) and θ ≈ 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 A 2. Figure 5(b) shows p θ for a non-zero rotlet dipole ( λ = 133.5)at ϕ = 0.18, which again demsontrates a significant reducuction of differences between different swimming modes due to the rotlet dipole.In particular, for all β, a single peak in p θ 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 λ = 0) with the formation of a two-layered structure.The probability P ⊥ of squirmers to be oriented nearly perpendicular to the walls is displayed in Fig. 6 for different values of β as a function of ϕ.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 ϕ, the probability of perpendicular orientation to the walls is zero for pushers and neutral squirmers, while P ⊥ ≈ 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 perpendicularto-the-wall orientation, despite of such predictions for a single puller [66,72].This is due to collective effects between pullers, which will be discussed in Section III A 4 below.
As the volume fraction of squirmers increases, P ⊥ 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 ⊥ first increases with increasing ϕ, but then closely follows P ⊥ for pullers when ϕ ≳ 0.4.For pushers, there is only a slight increase in P ⊥ with increasing ϕ.Nevertheless, we expect that for large volume fractions ϕ ≳ 0.6, orientational differences between squirmers with various swimming modes essentially disappear.FIG.6: Probability P ⊥ of squirmers to be aligned with the y axis (i.e., wall normal vector) for different volume fractions and swimming modes at λ = 0. Data for pushers with β = −5 (red), neutral squirmers with β = 0 (blue), and pullers with β = 5 (green).The probability is calculated by integrating the orientation distributions in Fig. 5(a) over the ranges [0, 25] and [155, 180] degrees.

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 where r ij is the vector between two centers of mass of squirmer pairs, and M = A/(2πr) is a normalization factor (A is the area of the considered 2D plane) such that g(r) → 1 for r → ∞.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 ϕ = 0.18 and without rotlet dipole, λ = 0.Only the suspension of pullers with β = 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 ϕ = 0.18, while the other squirmer types yield a gas-like phase, see Table I.Interestingly, the x-y perspective plot in Fig. 2 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 ϕ = 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 ϕ = 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 I), 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 β, indicating that the internal structure is nearly independent of the swimming mode at large ϕ.

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 α of the relative angle α between orientation vectors of two neighboring squirmers within a defined cutoff radius as where and C is the normalization factor so that the integral of p α is equal to unity.Figure 8(a) shows the distribution p α for various β at ϕ = 0.18.p α for pullers displays two peaks located at α ≈ 63 o and α ≈ 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 α ≈ 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 α ≈ 25 o and α ≈ 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 α for an active rotlet dipole with λ = 133.5 is shown in Fig. 8(b).For all β values, p α is centered at α ≈ 90 o , indicating that hydrodynamic interactions favor the perpendicular orientation between squirmer pairs.Note that for pushers, the central region in p α is flatter than that for pullers and neutral squirmers.

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 where < • • • > denotes the average over all squirmers and all times in the simulated system.Here, the exponential decay applies for small ∆t ≲ 1/D R .Figure 9 shows the dependence of D R on ϕ, β, and λ.In most cases, D R /D 0 R > 1, indicating that the rotational dynamics of squirmers within a suspension is enhanced in comparison with the rotational diffusion D 0 R 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 β = 0, ϕ = 0.56, and λ = 0, rotational diffusion is reduced, D R < D 0 R , which is likely due to the crowding of squirmers at large volume fractions [compare Fig. 2(a)].Figure 9(a) ( λ = 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 ϕ first exhibits an increase, followed by a maximum and subsequent decrease.With increasing ϕ, the collision rate between squirmers increases, resulting in the enhance-ment of D R .However, at large ϕ, 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 λ = 0 and λ = 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 λ = 133.5 switch between the two layers (i.e.migrate from one wall to the other) much more frequently than squirmers with λ = 0, which will be discussed in Section III B 4. The ability of interchange between the two layers increases significantly the effective rotational diffusion.The D R curves for λ > 0 are similar for different swimming modes β, in agreement with the structural characteristics for suspensions of squirmers with the rotlet dipole, as discussed in Sec.III A above.FIG.10: Mean-squared displacement of squirmers for ϕ = 0.18 (dashed lines) and ϕ = 0.56 (solid lines) (a) without the rotlet dipole λ = 0 and (b) with the rotlet dipole λ = 133.5.Curves for pushers are in red, for neutral squirmers in blue, and for pullers in green.Black dotted lines indicate power-laws ∼ t 2 in the ballistic regime and ∼ t in the diffusive regime.

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 ∼ t, as expected.For ϕ = 0.18 and λ = 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 I).The effective diffusivity of squirmers at ϕ = 0.56 is lower than at ϕ = 0.18 due to crowding and the formation of large clusters.The transition from the ballistic to the diffusive regime for λ = 0 in Fig. 10(a) occurs around tD 0 R ≃ 0.2 − 0.5.Theoretically, this transition should take place around tD R ≃ 1 [78], which is consistent with the values of D R /D 0 R ≃ 2 in Fig. 9(a).Note that the ballistic-todiffusive transition for squirmers with rotlet dipole occurs at tD 0 R ≳ 0.1, which is also consistent with the values of D R /D 0 R ≃ 10 in Fig. 9(b).Furthermore, mean-squaredisplacement curves for λ = 133.5 in Fig. 10(b) are similar to each other for various β at a fixed volume fraction, confirming once more that the rotlet dipole nearly offsets the effects of different swimming modes.

Average squirmer speed
To further characterize the mobility of squirmers for different conditions, we also study the average squirmer speed v.For a single squirmer i, the instantaneous swim speed is calculated as with time interval ∆t = 0.022/D R (significantly smaller than 1/D R ), during which a free squirmer moves a distance b z .v is then obtained by an ensemble and time average as a function of the volume fraction ϕ.
As shown in Fig. 11, pushers are the most mobile, displaying the largest average speed v among all squirmer types.For the case without rotlet dipole in Fig. 11(a), v decreases with increasing ϕ due to crowding at large volume fractions and the transition to MIPS or clustering regime.For all studied β 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 λ = 0. Furthermore, the decay in v for λ = 133.5 with increasing ϕ is slower than for λ = 0. Interestingly, for ϕ ≲ 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,76,79].The results in Fig. 11  presence of rotlet dipole enhances the mobility of squirmers and delays the transition to MIPS as a function of ϕ.

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τ R for different ϕ without rotlet dipole.For pushers, P m is close to unity at small ϕ, and decreases to P m ≈ 0.6 at ϕ = 0.56.Previous studies [45,57,79,80] 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 ϕ 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 ≈ 0.1 as the volume fraction is increased, indicating that their residence at the walls is more longlived than for pushers.
FIG. 12: Fraction P m of squirmers that switch between the two layers at least once during a simulation as a function of ϕ for λ = 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 λ = 133.5,independently of ϕ and β.
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 A 3, ass well as the formation of stable flower-like structures, as shown in Sec.III A 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 λ = 133.5,P m = 1 for all simulated β and ϕ 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) [31,38,39,[46][47][48][49] or 3D systems with periodic boundary conditions [37,[50][51][52][53].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 ϕ 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] [38,46].Our simulations also suggest that hydrodynamic interactions suppress motilityinduced 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 ϕ = 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 [25,29].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 flowerlike 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 ϕ, promoting MIPS at low ϕ.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 ϕ.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 ϕ ≲ 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 β, to a large extent.Even at relatively low volume fractions, the presence of a rotlet dipole of dimensionless strength λ ̸ = 0 nearly eliminates differences in the behavior of pullers, neutral squirmers, and pushers.With the rotlet dipole, squirmers assume a parallelto-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 λ = 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 λ = 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 speciesspecific, 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 [37,[50][51][52][53].

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] where x = l/l m ∈ (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 where κ is the bending rigidity, σ i is the area corresponding to vertex i in the membrane triangulation, H i is the mean curvature at vertex i, and H i 0 is the spontaneous curvature at vertex i.The mean curvature is discretized as H i = n i • j(i) σ ij r ij /(σ i r ij ), where n i is the unit normal at the vertex i, σ i = j(i) σ ij r ij /4, j(i) spans all vertices linked to vertex i, and σ ij = r ij (cot θ 1 +cot θ 2 )/2 is the bond length in the dual lattice with θ 1 and θ 2 being the angles at the two vertices opposite to the edge ij in the dihedral.The spontaneous curvature H i 0 is set locally after the triangulation of the spheroidal surface of the squirmer.

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 particlebased 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 where a, γ, and σ are the force amplitudes, r ij = r i − r j is the relative position vector, r ij = |r ij |, rij = r ij /r ij , and v ij = v i − v j is the velocity difference.ξ ij = ξ ji is a symmetric Gaussian random variable with zero mean and unit variance, and ∆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, F D and F R are related through the fluctuation-dissipation theorem [63] as The weight functions are defined as 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 where m i is the mass of particle i.Time integration is performed using the velocity-Verlet algorithm.

CONFLICTS OF INTEREST
There are no conflicts to declare.

FIG. 1 :
FIG. 1: Schematic of the simulation setup.(a) Sketch of a spheroidal squirmer model.The orientation vector 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.e τ and e ζ 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.
1(b).The number density of fluid particles is n f = 320/b 3 z , 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 η = 403.2√ mk B T /b 2 z .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 ϵ LJ = k B T and σ 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 × 10 −4 k B T /m/b z , which is close to the theoretical prediction of D R = 3.52 × 10 −4 k B T /m/b z ) Average properties of E. coli bacteria swimming in water correspond to U 0 = 29 µm/s, f d = 0.42 pN, l d = 1.9 µm, b z = 1.5 µm, b x = 0.5 µm, and η = 10 −3 Pa•s FIG. 2: Snapshots of the emergent structures for different active stresses β, volume fractions ϕ, and rotlet dipole strengths λ.(a) Simulated structures for λ = 0, ϕ = 0.18 and 0.35, and β ∈ {0, ±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 ϕ = 0.35, β = 5, and λ = 0. See also Movie S6.(c) Snapshots of the simulated systems for λ = 133.5,ϕ = 0.35 and β ∈ {0, ±5}.See also Movies S7 -S9.(d) Squirmer structures within the upper and bottom layers of the slit for ϕ = 0.35, β = 5, and λ = 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.

FIG. 3 :
FIG. 3: Cluster size distributions for different volume fractions of squirmers.(a) ϕ = 0.18, (b) ϕ = 0.35, and (c) ϕ = 0.56.The curves with different colors represent various swimming modes: pusher with β = −5 (red), neutral with β = 0 (blue), and puller with β = 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 TableI.(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 ( λ = 0), while the insets show data from the corresponding simulations with λ = 133.5.

2 .
Squirmer distribution between the walls

Figure 4
Figure4shows distributions p y of the squirmers' centerof-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 ≈ ±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-

FIG. 5 :
FIG. 5: Orientational distribution p θ of squirmers, where the angle θ is between the swimmer orientation vector and the y axis (wall normal axis).θ = 0 and θ = 180 correspond to the orientation perpendicular to the walls, while θ = 90 represent the orientation parallel to the walls, see the insets.p θ is presented for ϕ = 0.18 with the rotlet dipole strength (a) λ = 0 and (b) λ = 133.5.The three different β values correspond to pushers (β = −5), neutral swimmers (β = 0), and pullers (β = 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.
FIG. 7: (a) 2D radial distribution function g(r) of squirmers with ϕ = 0.18 and λ = 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 λ = 133.5 and the same ϕ.(b) Snapshot of the simulated system for ϕ = 0.18, λ = 0, and β = 5, which illustrates the most frequent squirmer structures, see also Movie S3.(c) g(r) of squirmers for ϕ = 0.56 and λ = 0.The legend for the curves is the same as in (a).

FIG. 8 :
FIG. 8: Distribution p α of relative angle between two neighboring squirmers for different swimming modes at ϕ = 0.18.α is the angle between orientation vectors of a pair of squirmers, whose separation is smaller than 3b z .(a) p α for the case of no rotlet dipole ( λ = 0), with the insets illustrating configurations of squirmers at the two peaks with α ≈ 63 o and α ≈ 113 o .(b) Relative angle distribution for λ = 133.5.

FIG. 9 :
FIG. 9: Effective rotational diffusion D R as a function of ϕ for (a) λ = 0 and (b) λ = 133.5.D R is obtained by fitting auto-correlation functions of squirmer orientation using Eq.(14), and normalized by the rotational diffusion D 0 R of a spheroid.Data for pushers (red), for neutral squirmers (blue), and for pullers (green).
FIG. 11: Average squirmer speed v as a function of ϕ for (a) λ = 0 and (b) λ = 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).
The area and volume conservation constraints are represented by the potential [55] U a+v = k a (A − A tot 0 ) 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 tot 0 and V tot 0 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 0 m 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 µ 0 = 1.44×10 5 k B T /b 2 z , the bending modulus κ = 250k B T , the area-constraint coefficients k d = 1.6 × 10 3 k B T /b 2 z and k a = 8 × 10 3 k B T /b 2 z , the volumeconstraint coefficient k v = 3.2 × 10 4 k B T /b 3 z , the total area A tot 0 = 5.36b 2 z , and the total volume V tot 0 = 1.05b 3 z .In all simulations, k B T = 1 and b z = 4.
DPD parameters for the interactions between fluid particles and between fluid and wall particles are a = 200k B T /b z , γ = 80 √ 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, γ = 100 √ mk B T /b z , s = 0.1, and r c = 0.175b z .The time step for integration is 1.25 × 10 −3 b z 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.