Gyrotactic cluster formation of bottom-heavy squirmers

Squirmers that are bottom-heavy experience a torque that aligns them along the vertical so that they swim upwards. In a suspension of many squirmers, they also interact hydrodynamically via flow fields that are initiated by their swimming motion and by gravity. Swimming under the combined action of flow field vorticity and gravitational torque is called gyrotaxis. Using the method of multi-particle collision dynamics, we perform hydrodynamic simulations of a many-squirmer system floating above the bottom surface. Due to gyrotaxis they exhibit pronounced cluster formation with increasing gravitational torque. The clusters are more volatile at low values but compactify to smaller clusters at larger torques. The mean distance between clusters is mainly controlled by the gravitational torque and not the global density. Furthermore, we observe that neutral squirmers form clusters more easily, whereas pullers require larger gravitational torques due to their additional force-dipole flow fields. We do not observe clustering for pusher squirmers. Adding a rotlet dipole to the squirmer flow field induces swirling clusters. At high gravitational strengths, the hydrodynamic interactions with the no-slip boundary create an additional vertical alignment for neutral squirmers, which also supports cluster formation.

Gyrotaxis refers to directed locomotion resulting from a combination of gravitational and viscous torques in a flow [51][52][53]. In particular, it applies to the action of fluid vorticities on microswimmers [51,53]. As such, it is responsible for the dancing motion of Volvox algae [54], focussing in channel flow [55], the formation of layers [10] as well as patches of phytoplankton [56], and convective patterns in algae and bacteria [30,[57][58][59]. These examples show the importance of gyrotaxis for biological systems. They motivated us, in an earlier article, to perform hydrodynamic simulations using the method of multi-particle collision dynamics in order to systematically study the dynamics and gyrotaxis of bottom-heavy squirmers. The latter are generic model microswimmers [60]. Depending on the strength of gravity and bottom heaviness, we identified different dynamical states including inverted sedimentation, plumes and convective rolls, as well as spawning clusters.
In this article, we continue our work on gyrotaxis of bottom-heavy squirmers but now concentrate on moderate densities and strong gravity such that they float above but close to the bottom surface. Only with increasing gravitational torque gyrotactic cluster formation sets in as a subtle balance between hydrodynamic and gravitational torques. We thoroughly characterize the cluster formation of neutral squirmers for different areal densities and torques using different quantities such as the vertical density profile, mean cluster size, and radial distribution function. We also discuss the influence of squirmer flow fields from pusher and puller squirmers. While pusher squirmers do not exhibit noticeable cluster formation, larger torques are needed to observe it for puller squirmers. Furthermore, a hydrodynamic rotlet dipole realized, for example, by bacteria also weakens cluster formation.
The article is structured as follows. In Sect. 2 we present the squirmer model and describe in detail the hydrodynamics of this model swimmer under gravity. Furthermore, we introduce the method of multi-particle collision dynamics to simulate the flow field generated by the squirmers. In Sect. 3 we present our simulation results and we conclude in Sect. 4.

Methods
We simulate squirmer microswimmers under gravity using the mesoscale simulation technique of multiparticle collision dynamics. The squirmers move in a three-dimensional simulation box with no-slip walls at the bottom and at the top and periodic boundary conditions in lateral direction. Tables 1 and 2 sum-  Edge length of collision cell m0 Mass of fluid particle T0 Temperature of fluid Δt Duration of streaming step n fl Fluid particle number density η Dynamic fluid viscosity marize our physical and simulation parameters, respectively.

Squirmer model
We use the squirmer model introduced by Lighthill and Blake [61,62] in order to investigate gyrotactic cluster formation of microswimmers. The original idea of the model is to consider a spherical particle with a surface velocity field and expand it into appropriate base function to model the surface actuation of a ciliated microswimmer. A more general ansatz includes chiral surface flow and was presented in Ref. [63], which we follow here. We work with an axisymmetric tangential flow field. Then, the surface slip velocity for a squirmer with radius R that swims along the unit vector e can be expressed as Here P n refers to the first derivative of the nth Legendre polynomial P n . Often, one truncates this expansion after the second order. The mode B 1 is necessary for self-propulsion where v 0 = 2 3 B 1 is the swimming velocity. Additionally, we consider the modes B 2 and C 2 , which we express in terms of the squirmer parameter β = B 2 /B 1 and chirality parameter χ = C 2 /B 1 , see for example Ref. [64]. The parameter β switches between neutral squirmers (β = 0), pullers (β > 0) and pushers (β < 0).

Hydrodynamics of bottom-heavy squirmers
We collect all the flow and vorticity fields of a bottomheavy squirmer and thereby provide an overview on the hydrodynamics of our squirmer system. While we use periodic boundary conditions in our numerical scheme, we do not include them in the following presentation.
Flow and vorticity fields The flow field of the freely swimming squirmer resulting from the truncated surface field in eq. (1) is [63] where we have written the different contributions such that one immediately recognizes the source dipole (first term), force dipole (term with β) and rotlet dipole (term with χ).
Equation (2) shows that a squirmer always creates a source dipole field, which becomes the leading order for the neutral squirmer (∼ r −3 ). For a pusher or puller squirmer, a force dipole field is added (∼ r −2 ), whereas χ = 0 refers to the rotlet dipole. Pushers model propulsion originating from the back of a swimmer's body, such as rotating flagella, whereas pullers are a more suitable description for a breast-stroke type of motion. Pusher and puller fields are the most typical hydrodynamic modes considered in microswimmer systems [60,65]. Furthermore, the rotlet dipole field is an important aspect of many swimming organisms, for example, of bacteria with counter-rotating flagella and cell body. It has been studied in some previous works [64,66].
The complete hydrodynamic flow field generated and experienced by a collection of squirmers is strongly determined by their manifold hydrodynamic interactions and is therefore analytically untractable. However, we can still gain insights from the single squirmer fields and their different contributions, which can be ordered according to their radial decay. This can help to interpret numerical results of their collective dynamics. In our case, the self-created flow fields of the squirmers from Eq. (2) are not the only contributions. Also their gravitational force and bottom-heaviness generate flow fields. Furthermore, the no-slip boundary condition on the bottom surface affects the hydrodynamics by creating a wall-induced flow [67]. Fortunately, the Stokes equations are linear, which allows us to superimpose all occuring fields.
Contributions from gravity The first part comes from the gravitational force acting on the squirmer and is the same as the flow field of a passive sphere dragged through the fluid. Therefore, a squirmer under gravity always induces the flow fields of a stokeslet and a source dipole [68]. The latter adds to the source dipole of the squirmer [63]. Additionally, since the center of mass of the squirmer is shifted by a distance r 0 from its geometric center, it experiences bottom heaviness, which aligns the swimming direction e along the vertical e z . Due to the gravitational torque T bh = mgr 0 e × e z acting on the squirmer, it rotates with angular velocity Ω bh = T bh /(8πηR 3 ) and thereby induces a rotlet field. The velocity fields of both contributions are [60,63,67]: Here, we introduced the dimensionless parameters α = v 0 /v sed and r 0 /(Rα), with the bulk sedimentation velocity v sed = mg/(6πηR) and the center-ofmass shift r 0 . The parameter α measures the swimming velocity relative to the sedimentation velocity and r 0 /(Rα) is a unitless strength of the gravitational torque.
In these units the related angular velocity becomes [19]: The actual flow field of the squirmer together with the gravity-induced contributions constitute the hydrodynamic signature of our problem and give the relevant flow fields in leading order in 1/r. Importantly, the non-zero vorticities of these flow fields reorient nearby squirmers or the squirmer itself when it interacts hydrodynamically with a wall. Since the vorticity of a source dipole is zero, we only need to consider the vorticity from the gravity-induced stokeslet and rotlet as well as the vorticity from the force dipole of squirmers with β = 0. Using the definition of the vorticity, ω(r) = ∇ × u(r), the respective vorticities of the stokeslet, rotlet, and force dipole of a bottom-heavy squirmer are Note that the vorticity of the rotlet, ω R from eq. (7), vanishes for e → e z and therefore this contribution is typically small. Now, placing a second squirmer at position r 2 in the flow field of a sufficiently distant first squirmer, its angular and translational velocities are determined by the collected far field expressions. Using Faxén's theorem [68,69], the angular velocity up to terms proportional to 1/r 3 2 is given by 1 2 (ω S (r 2 ) + ω R (r 2 )+ω F d (r 2 ) +O(r −4 2 ) and the translational velocity becomes u sq (r 2 ) + u grav (r 2 ) +(R 2 /6) ∇ 2 u grav (r 2 ) + u bh (r 2 ) + O(r −4 2 ), where r 2 is the distance vector from the first squirmer to the second.
The collective dynamics in gyrotactic systems is determined by the competition and balancing of external torque and flow vorticity acting on the microswimmers. One can use this effect to focus swimmers in an external Poiseuille flow [55]. More recently, it has been shown that a prescribed periodic flow can lead to vertical trapping of gyrotactic motile cells [70]. In our case, however, the flow vorticity arises from the flow fields of the neighboring squirmers. Assuming two squirmers on the same height with distance r, they rotate each other away from the vertical by an angle ϑ due to the stokeslet vorticity of eq. (6), in leading order. Balancing this rotation by the angular velocity of eq. (5) due to the gravitational torque gives In the following we will identify stable clusters of squirmers in our simulations only when the gravitational torque is sufficiently large. We present a simple argument for this behavior. At the closest possible distance r = 2R, the angle where both angular velocities in eq. (9) balance each other becomes Since sin ϑ ≤ 1, such a stable balance requires the lower bound for the dimensionless torque. As we will show in Sect. 3.2.1 this gives a good estimate for the torque, where the clusters start to become compact.
Wall terms Close to a no-slip wall mirror multipoles have to be added to the hydrodynamic multipoles introduced above, in order to fulfill the no-slip boundary condition [67]. The leading far-field contribution of a squirmer under gravity is the stokeslet, which close to a no-slip wall turns into the Blake tensor [71,72]; the influence of its mirror multipoles decays with the inverse height. Furthermore, the squirming modes controlled by the parameters β and χ further affect the wall-induced flow field, which depends on the squirmer orientation. For example, vertically oriented pushers are repelled and pullers are attracted to the wall [21]. We are mainly interested in the vorticities of the wallinduced flow fields since they are directly related to gyrotaxis. We give the most relevant terms here following Refs. [21,67]. By symmetry the mirror image in the Blake tensor does not contribute to the wallinduced vorticity acting on a single squirmer. Therefore, the leading-order contribution for pushers and pullers is the reflected force dipole field, whereas for β = 0 and χ = 0 it is the reflected rotlet-dipole field. Second, a force quadrupole arises as a wall-induced mirror term in the three squirmer multipoles with strengths B 1 , B 2 and C 2 . In particular, this force quadrupole results from the reflected vorticity field of the source dipole [67] and, therefore, determines the orientation of the neutral squirmer in wall proximity [21]. We write the vorticity or angular velocity Ω wall as a function of the distance to the no-slip wall, Δz, and the angle ϑ of the squirmer orientation e with respect to the vertical [67]. Using cylindrical coordinates, we arrive at the following components For example, this angular velocity becomes zero for a vertically oriented neutral squirmer with β, χ = 0. The additional force-dipole flow field of pushers and pullers changes this stable orientation [21]. Furthermore, it is known that the rotlet dipole induces the wall component Ω wall z of eq. (14) that is responsible for circular swimming [66,73].

Multi-particle collision dynamics
Various successful strategies exist for treating the low-Reynolds-number fluid dynamics of suspended microswimmers. This includes direct discrete simulations, using finite-volume [74] or smoothed-profile [64] methods, the boundary-element method, where an integral equation is solved [66,75], or mesoscale techniques such as lattice-Boltzmann simulations [76,77]. In this article, we use the mesoscopic simulation technique of multi-particle collision dynamics [78][79][80], which is well established for hydrodynamic simulations of microswimmers, in particular, squirmers [34,41,42,[81][82][83][84][85]. In particular, this technique has successfully reproduced the correct squirmer flow fields [81,86] as well as their near-and far-field interactions [83,87]. It simulates hydrodynamic flow fields as solutions of the Navier-Stokes equations and also includes thermal noise [88,89]. The fluid consists of coarse-grained but point-like fluid particles of mass m 0 . The algorithm of multi-particle collision dynamics applies a sequence of alternating streaming and collision steps to these particles [80].
During the streaming step with duration Δt, the fluid particles move ballistically. During this step one also implements the no-slip boundary condition at bounding walls using the bounce-back rule [90]. It is also used to prescribe the tangential velocity field of Eq. (1) on the squirmer surface, which enables self-propulsion. Furthermore, the suspended squirmers are treated as rigid bodies and their positions and orientations are updated using velocity-Verlet integration [91], which accounts for the gravitational forces and torques acting on them. In order to model the steric repulsion between squirmers and between squirmers and bounding walls, we use a steep WCA potential. We always choose the Reynolds number as small as computationally feasible and arrive at Re = 0.17 [60]. This is still a good approximation of the Stokes flow regime. Thus, the acceleration of the swimmer due to the external force is balanced by a friction force, consistent with force-free swimming. For the collision step we introduce a length scale a 0 to indicate the range of fluid-particle collisions, which take place in cubic cells of side length a 0 . On average, we place n fl fluid particles into the cells. Both values of Δt, as well as n fl , determine the fluid properties such as viscosity. Further simulation details are described more exhaustively in our previous works [42,60,91].
The collisions of the fluid particles inside a collision cell are governed by a specific collision operator or collision rule, which redistributes the relative momentum between the particles keeping the total momentum fixed. Several different collision rules exist for multiparticle collision dynamics [80]. We use the MPC-AT+a algorithm which assures Galilean invariance and angular momentum conservation [80,92]. It is based on the Andersen thermostat [80,89], which keeps the thermal energy of the fluid at the fixed value k B T 0 . If a collision cell overlaps with a suspended squirmer or a bounding wall, this volume is filled with virtual particles [91,93] to ensure that the mean density in the collision step is always n fl . This approach also improves the accuracy of the flow fields and also prevents depletion forces [94]. All in all, the cell length a 0 , the mass m 0 of a fluid particle, and the thermal energy k B T 0 form the units of length, mass, and energy. From this we can derive other units, in particular, for time a 0 m 0 /k B T 0 and velocity k B T 0 /m 0 .

Implementation and parameters
We simulate the MPCD fluid inside a rectangular system with height H and a quadratic cross section with edge length L. We use periodic boundary conditions in lateral direction and two no-slip walls in vertical direction. Usually, we have H = L = 108a 0 , i.e., a cubic simulation box. For the fluid we use the duration of the streaming step Δt = 0.02a 0 m 0 /k B T 0 and the fluidparticle number density n fl = 10. This choice of parameters results in a viscosity of η = 16.05 [95]. Our squirmers have B 1 = 0.1 k B T 0 /m 0 and, unless stated otherwise, β = χ = 0. The resulting swimming velocity is v 0 = 0.067 k B T 0 /m 0 , which corresponds to an active Péclet number Pe = Rv0 DT = 330, comparable, for example, to bacteria [37,96]. Furthermore, the squirmers have a radius R = 4a 0 , and the relevant ballistic time scale v 0 /R amounts to around 3000Δt [60]. To set the parameter α, we choose a moderate strength of gravity that always creates a sedimentation velocity larger than v 0 , typically α = 0.8. This prevents squirmers from escaping from the bottom of the system to the top wall.
Due to the high computational cost of hydrodynamic simulations, we compile software for running on multiple processing units simultaneously, i.e., the program is significantly parallelized. Earlier simulations were performed on a high-performance computing cluster using approximately 150 CPU cores. For later simulations we switched to a CUDA implementation, which we ran on graphics processor units.

Initial condition
We tested two different initial conditions. Either, the squirmers were distributed randomly in the simulation box or we initialized them on a regular grid in the midplane z = H/2, all pointing upwards with cos ϑ = 1. This had little effect on the clustering behaviour that we focus on in the following. The only difference we observed was that for the random initialization more squirmers escaped to the top wall for strong bottom-heaviness. Otherwise, these squirm-ers have no influence on the rest of the simulation, as they stay at the top plate without interfering during the rest of the simulation. Choosing random orientations at the beginning of the simulation or initializing the squirmers in another plane with constant z also brought no qualitative change.

Results
We present our results on gyrotactic clusters and the associated parameter studies in the following. The decisive parameters are the velocity ratio α and the torque measure r 0 /Rα that we introduced above. In previous studies [34,60] we investigated large-scale convection and plumes that move over the entire system height. In contrast, we here look at the cluster formation of bottom-heavy squirmers that are constrained to the bottom wall due to strong gravity (α < 1). However, the squirmers are still capable of moving vertically, in contrast to monolayer formation for α below 0.1 [47]. Thus, we observe the emergence of dynamic clusters from a squirmer suspension that float above the bottom wall.

Phenomenology of gyrotactic clusters
At the outset we discuss the basic properties of the emergent structures that form under gravity and bottom-heaviness. We present snapshots of the squirmer system and discuss the vertical density profiles and dynamics of the clusters. However, we start by illustrating the coupled dynamics of two squirmers.

Gyrotaxis of a squirmer pair
Gyrotaxis is often discussed in the context of collective dynamics of microswimmers, but already a twoparticle system can illustrate the combined effect of a gravitational torque and fluid flow.  (6), just by considering the stokeslet flow fields and the gravitational torque, and thereby reached good agreement with their experimental data [54]. They found that the dynamical system of two hydrodynamically interacting microswimmers contains oscillatory trajectories when reorientation due to the neighbor's stokeslet vorticity and due to the own gravitational torque balance each other. In the following, we reproduce such oscillations in our squirmer simulations. We show the simulated oscillations for a pair of squirmers in Fig. 1 and outline the mechanism in the following.
The simulations are performed for α = 0.8, where single neutral squirmers float above the bottom wall at a height of several squirmer radii R [21]. Due to bottomheaviness, they experience a gravitational torque proportional to − sin ϑ, where ϑ is the polar angle between squirmer orientation and vertical (see last snapshot in  Fig. 1). Thus, squirmers in bulk tend towards the vertical where the gravitational torque is zero. Two squirmers with such a vertical orientation are shown in the first snapshot in Fig. 1. The neighbor can tilt the orientation vector towards itself due to the vorticity of the stokeslet, as shown in the second snapshot. Balancing reorientations due to vorticity and torque, the squirmers are tilted and thus move towards each other. They eventually pass each other (snapshots 3 and 4), whereupon the sense of the vorticity reverses and the motion repeats itself in the other direction, which is shown in snapshots 5-8. The external torque is important to stabilize the pair so that the distance of the squirmers does not increase too much after they have passed each other, which would otherwise destroy the oscillations [54]. During a passing trajectory of a squirmer, it sinks to a lower height due its tilted orientation and swims around and slightly underneath the other squirmer. The squirmers do not penetrate each other and the overlap seen in the snapshots is merely an effect of the projection. Interestingly, even for large torques the squirmer orientations are not completely frozen, but pair formation still occurs with a smaller oscillation amplitude. Since the gravitational torque approaches zero for a vertical orientation, a small tilt and thereby horizontal motion can still be achieved by the vorticity in the flow field.

Formation of clusters
Now, we give an overview over what happens to the oscillatory motion in systems with a larger number of neutral squirmers.
In Fig. 2 we show snapshots of squirmer configurations seen from the top for different torque strengths and squirmer numbers. We used a constant α = 0.8. The third row with r 0 /Rα = 0.32 is available in the supplement as videos M1-M3. These videos also show a side view of the system.
The snapshots show that squirmer clusters emerge for different squirmer numbers and torque values and spread over the whole plane with a typical cluster distance Δ as indicated in the central snapshot for They tilt more strongly and, therefore, exhibit stronger horizontal motion. As explained above, squirmers swim past each other and cannot form stable structures.
In Fig. 3 we show the density profiles of the squirmer system in the horizontal plane. In the top panel the torque increases from top left to bottom right while N = 110.
The density profiles were obtained by averaging over 10 5 time steps and over all squirmers in the bottom half of the system. At zero torque the density profile is unstructured and almost uniform. There is some patchiness which we attribute to the aligning effect of the bottom wall, leading to floating squirmers that stay in the same spot for longer times [21].
At r 0 /Rα = 0.16 we observe that squirmers form clusters which are stable during a longer time (Fig. 3, top). There is still considerable horizontal motion and squirmers frequently join and leave the clusters, therefore the space around the density peaks is still well explored. With increasing torque (e.g., last two rows in Fig. 2) the clusters become more static and compact. This is understandable since a stronger vertical alignment obstructs the escape of squirmers from clusters because the horizontal velocities are smaller. According to the two-dimensional density profile for r 0 /Rα = 0.32 in Fig. 3, top, the squirmers explore less space in the horizontal plane and the clusters indeed become more compact, which is most clearly seen for r 0 /Rα = 0.63. The same effect occurs for the oscillating squirmer pairs, i.e., the horizontal extent of the clusters decreases with increasing torque. In parallel, the total number of clusters increases as the right-most column of Fig. 2 shows for N = 110.
The bottom plot of Fig 3 shows how the clustering of squirmers evolves when their number increases at constant torque. An increasing N can also be expressed as an increasing area fraction NπR 2 /L 2 , where L is the edge length of the simulation box. We observe a stronger horizontal motion, which is visible from the disappearance of unvisited space and less pronounced peaks in density. This makes sense since each additional squirmer introduces a new stokeslet, which reorients neighbors away from the vertical.
Cluster dynamics Drescher et al. coined the term "minuet dance" for the oscillating two-squirmer motion [54]. For larger clusters the motion looks less regular and consists of squirmers switching their positions inside the cluster as videos M1-M3 show. We stress again that the clusters become more static with increasing acting torque. The cluster reorganization sometimes creates symmetrical structures, such as pentagons or hexagons, as video M2 illustrates. However, these configurations do not persist. A similar effect was reported in experiments with active emulsion droplets [97]. They form clusters that also rearrange frequently at high Péclet number.

Vertical density profile
Stacking The top-view snapshots in Fig. 2 show occasional stacking of squirmers, visible as overlaps, especially at the higher squirmer number in the right column. We stress that these overlaps are due to the projection of the squirmers on one plane. They only occur when the squirmers move on different heights and, of course, the squirmers do not penetrate each other. In order to investigate the vertical structure more closely, we show the vertical density profiles in Fig. 4 for different torque values. The stacking can clearly be seen as a double peak at the two higher torque values. More precisely, a bilayer is formed since the distance between the two peaks in Fig. 4 for r 0 /Rα = 0.32 and 0.63 is approximately one squirmer diameter.
Collective sinking Note that the vertical density profile at zero rescaled torque reaches to larger heights compared to the non-zero values. Thus, the squirmers are more often found at larger heights. This is the case, even though the orientational distribution gets more strongly peaked at the vertical orientation (cos θ = 1) with increasing r 0 /Rα (see inset of Fig. 4). The reason for this counter-intuitive behaviour is a reduced translational friction coefficient per squirmer when they form clusters [48,98]. At zero torque clusters do not exist and squirmers sink individually due to gravity. In contrast, for higher torques they sediment faster within clusters since hydrodynamic friction is reduced due to hydrodynamic interactions, as was studied in detail in Ref. [60]. This leads to stronger effective sinking, even though the alignment is more vertical.

Parameter study
We present a more detailed parameter study of squirmer clustering due to gyrotaxis. We characterize clusters by their radius and mean particle number, look at the radial distribution function, study the influence of the squirmer flow field beyond the neutral squirmer, and look closer at the influence of gravity by varying α.

Cluster radius and size
In order to demonstrate how the gravitational torque influences the sizes of the emerging clusters, we define measures for the cluster extension and size. We plot them versus the torque value in Fig. 5 for different squirmer numbers N . First, to quantify the cluster extension, we introduce a mean cluster radius by determining the average distance of a squirmer to the cluster center at r, |r − r| cl . Here, besides the time average, . . . cl also means averaging over all squirmers in a cluster and then over all clusters. Second, for the cluster size we calculate the mean number of squirmers in a cluster, N cl . We also consider the standard deviation ΔN cl = [ (N cl − N cl ) 2 ] 1/2 to have a measure how much the cluster size varies between the different clusters and also in time. To determine all these quantities, we define a cluster as a set of squirmers, where for each squirmer the distance to at least one other squirmer in this set is less than 2R + R/4. Finally, we only take the last 2.5 · 10 5 timesteps into account where the system has reached steady state.
Both, the cluster extension (a) and size [inset of (b)], show a pronounced maximum, especially for N = 110, where clusters are largest. The increase at small torque values leading up to this maximum is quite sharp. As already stated, in this regime clusters start to form due to gyrotaxis and are more short-lived, because the external torque is not strong enough to stabilize them. At and close to the maxima in extension and size, the clusters are rather loosely bound, which is visible in the snapshots for torque value r 0 /Rα = 0.32 in Fig. 2 and in the corresponding videos M1-M3. The clusters strongly vary in size and also in shape with ongoing time. This behavior is reflected by the peak in the normalized standard deviation ΔN cl / N cl of the mean cluster size as shown in the main plot of Fig. 5(b).
Beyond the maxima of cluster radius and size, both quantities decrease with increasing torque. The clusters become more compact and exhibit less variation also in time, as the decrease in the normalized standard deviation ΔN cl / N cl shows. This behavior is most clearly visible for N = 110. For further increasing torque all three quantities then become nearly constant. While the decrease in N cl is not very pronounced, it is stronger for the mean cluster radius, in particular, for N = 110. Thus the clusters become compactified, while their total number increases, as the last row in Fig. 2 nicely demonstrates. Interestingly, the compactification of the clusters agrees with the location of the dashed line. According to eq. (11) it shows the value of the torque that is just able to balance squirmer rotation due to the stokeslet vorticity of a neighboring touching squirmer. Since squirmers in a cluster also hinder each other from moving, the cluster is stabilized and forms compact objects.
The reason for the decrease of the cluster radius is similar to the two-squirmer case described in Sect. 3.1.1. As the torque increases, squirmers are more vertically oriented and thus their horizontal mobility is further reduced. This has the effect that the motion of squirmers within each cluster is more constrained just as the oscillation amplitude of the squirmer pair also decreases with the torque. In order to follow the vorticity acting on the squirmers, they eventually also evade in the negative z-direction and the cluster assumes a stacked arrangement, as we have seen in Fig. 4.

Radial distribution function
To investigate squirmer ordering within the clusters and also how the clusters position themselves relative to each other, we calculate the radial distribution function g(r). Here, r is the squirmer distance in a horizontal plane, after projecting all squirmer positions on this plane. Thus, we treat the squirmer system as a twodimensional system. The two-dimensional radial distribution function for N squirmers in an area A with positions r i is defined as [99]  where r = |r i − r j | and the ensemble average is performed over all recorded squirmer configurations in the last 2.5 · 10 5 time steps of the simulations. Note that a constant value of g(r) = 1 corresponds to the random ordering of an ideal gas. The radial distribution function gives meaningful results for distances up to ca. L/2, thus half the linear box size. Figure 2 shows that the cluster distance Δ can assume such values. In order to rule out finite-size effects, we have doubled the linear size of our simulation box and use L = 216a 0 . This increases the crosssectional area by a factor of 4. Since we simulate 256 squirmers, the area fraction in the horizontal plane is the same as when we use N = 64 squirmers for our regular box size.
In the main plot of Fig. 6 we show the radial distribution function g(r) for the larger system for several torque values. In all curves a peak at the nearest neighbor distance r = 2R is visible indicated by the dashed vertical line on the left. The blue curve for zero torque then assumes the constant value one at larger distances reflecting the randomly distributed squirmer positions. Already at r 0 /Rα = 0.16 a small maximum develops around the next-nearest neighbor distance 2 √ 3R = 3.46R (right dashed line), which becomes more pronounced with increasing torque. In particular, for the two highest torque values, where compact squirmer clusters occur, sharp peaks are visible. Note that both peaks are broadened towards smaller distances due to the partial overlap of the squirmers in the projection plane, when they assume different heights. The overlap is clearly visible in the snapshots of Fig. 2 as already discussed earlier.
Now we discuss the characteristic distances Δ between the clusters as revealed by the radial distribution function. Already at r 0 /Rα = 0.16 a weak minimum appears around r = 6R followed by a very shallow maximum located between r = 10R and r = 15R. This signals the appearance of squirmer clusters separated from each other. Increasing the torque to r 0 /Rα = 0.32 (green curve), this third peak becomes more pronounced and its position shifts to smaller values. Finally, for the two highest torque values (r 0 /Rα = 0.63 and 0.95), where clusters show compact packing, the peaks are identical. They are further shifted to smaller distances and their heights are larger. Interestingly, a further maximum for both torque values is visible at around r = 17R and 18R, respectively. It corresponds to the next-nearest neighbor shell of the cluster ordering. The maximum cannot be seen in the system with regular size, since all clusters are nearest neighbors due to the periodic boundary conditions. Now, we focus on the the average cluster-cluster distance Δ in the inset of Fig. 6. We plot it versus the dimensionless torque for different squirmer numbers N in our regular system size and for the larger system. For Δ we use the position of the first maximum in g(r) connected to cluster ordering. With increasing N the curves slightly shift downwards but overall we can state that the cluster distance does not depend on the area density of the squirmers. As we already realized, at small torques, when the clusters start to form, the cluster distance decreases with increasing torque but once compact clusters are formed at r 0 /Rα = 0.64 and beyond it stays constant.
The probable explanation for the decrease in Δ is that clusters get smaller at larger torques [see Fig. 5 (a)]. Thus, there are more clusters for the same number of squirmers, and the space between them has to shrink. This connects the decrease in cluster distance to the clustering mechanism and to gyrotaxis, which suggests that Δ is an instrinsic quantity which emerges from the system parameters. We confirmed this by letting the end configuration for torque r 0 /(Rα) = 0.28 evolve further at an increased torque r 0 /(Rα) = 0.63. Indeed, the typical cluster distance decreases to the same value as observed for the random initial condition.

Influence of squirmer flow fields
We consider two squirmer modes in addition to the B 1 mode that modify the squirmer type. First, the B 2 mode quantified by the parameter β = B 2 /B 1 creates a force dipole in the far field and second, the C 2 mode introduces a rotlet dipole in the far field and is tuned via the parameter χ = C 2 /B 1 , as we explain in Sect. 2.1.1.
We again probe the order of the system by calculating the radial distribution functions for pushers, pullers, and rotlet-dipole squirmers. They are shown in Figs. 7 and 8 together with density profiles. We also provide videos of the puller squirmers for β = 2 at torques r 0 /Rα = 0.32 (M5) and r 0 /Rα = 0.63 (M6) as well as the rotlet-dipole squirmers with χ = 1.0 (M7). B 2 mode Pusher and puller squirmers generate forcedipole flow fields through which they also interact with each other [65,67]. In particular, the additional flow vorticity presented in Eq. (8) competes with the external and other hydrodynamic torques and thus contributes to the gyrotactic mechanism. As a consequence, the squirmer type β has a profound influence on the radial distribution function g(r), as Figs. 7 a) and b) demonstrate. At the lower torque value the peak corresponding to the mean cluster distance for neutral squirmers has disappeared completely for β = ±2 and clustering does not exist. The vertical density profiles in Fig. 7 c) underline this: pushers and pullers reach larger heights while neutral squirmers are concentrated close to the botton since their clusters cause collective sinking. In accordance with these findings, video M5 for puller squirmers with β = 2 does not show clustering.
Doubling the external torque [Figs. 7 b) and d)] still does not lead to visible clustering for the pusher system. This agrees with previous findings that cluster formation in pusher systems is generally not favored [40,60,100]. However, the behavior of the puller system (β = 2) changes drastically for the higher torque. First, we now see a weak peak at a distance 8R in the radial distribution function in Fig. 7 b) indicating the existence of squirmer clusters. Indeed, we also observe small clusters in video M6. This is in agreement with the density profile in Fig. 7 d), which is now more concentrated close to the bottom wall due to the collective sinking in clusters in contrast to the smaller torque value. The density profile also has two pronounced peaks, one of which is directly at the bottom wall. They indicate the stacked configurations of clusters typically observed for neutral squirmers already at smaller torques. Interestingly, in g(r) an additional, pronounced peak at a distance below the nearest-neighbor peak at r = 2R has developed, which we discuss further below.
The puller clusters have a different structure compared to the neutral squirmer clusters. Video M6 shows pullers in a cluster typically pointing towards the cluster center. This agrees with the typical hydrodynamic puller-puller attraction along the swimmer axis and mutual reorientations towards each other due to the force-dipole vorticity [101]. The tilt of the puller squirmer towards the horizontal is also caused by its near field close to a no-slip wall [102]. The pronounced tilt of the puller axis away from the vertical is also visible in the distribution of vertical orientations in the inset of Fig. 7 d). Even a weak local maximum at cos ϑ < 1 appears.
Like neutral squirmers, puller clusters develop a stacked configuration. The second peak in the density profile of the pullers in Fig. 7 d) belongs to a second layer of squirmers on top of the bottom layer. Note that the inwards pointing squirmers in the bottom layer induce fluid flow towards the cluster center such that the squirmer from the second layer settles there. Thus, the cluster has a pyramidal three-dimensional structure. Assuming a pyramid with an equiliateral triangle as the base gives a horizontal distance of √ 3/2R between the nearest squirmers from the upper and lower layers and a vertical distance of 3/2R. We indicate these values as vertical dashed lines in Figs. 7 b) and d), respectively. They agree well with the peaks found from the simulations. C 2 mode In Fig. 8 a)-c) we respectively show radial distribution functions, vertical density profiles, and distributions of horizontal velocities for neutral squirmers (β = 0) and several rotlet-dipole parameters χ. Clusters form for all values except χ = 3.3, where no peak is visible in the radial distribution function beyond the nearest-neighbor distance. Furthermore, only the vertical distribution of squirmers with χ = 3.3 extends to larger heights, whereas it is concentrated close to the bottom for the other χ values, typical for squirmer clusters. Thus, strong rotlet-dipole flow prevents clus-ter formation. In addition, Fig. 8 c) shows that with increasing χ the horizontal velocity components increase, which also contributes to the disappearance of clusters at high χ. This is also indicated in Fig. 8 d), which plots N cl versus χ (green triangles). The disappearance of the clusters also depends on other parameters. At r 0 /(Rα) = 0.32 the clusters already disappear at χ = 1.0, whereas puller clusters at β = 2 were not observed for the simulated non-zero χ.
For the value r 0 /(Rα) = 0.63 depicted in Fig. 8, an inspection of video M7 reveals that already at χ = 1.0 squirmers stay more loosely bound to each other and clusters are less compact. Furthermore, we observe that squirmers in a cluster swirl around the cluster center in anti-clockwise direction. In Eq. (14) we show the angular velocity, which a rotlet dipole close to a no-slip wall experiences due to the image flow. It is known to induce circular trajectories [66], for example, of E.coli bacteria [73,103]. While E.coli swim in circles with their orientation vector parallel to the wall, our squirmers have a strong vertical bias. According to Eq. (14) this corresponds to an angular velocity with reversed sign, as long as cos ϑ exceeds 1/3 ≈ 0.58. This is indeed the case, as the inset in Fig. 8 (b) shows. The constant inplane reorientation combined with the squirmers' selfpropulsion results in circular swimming trajectories and in our simulations creates the swirling clusters.
We introduce the swirling parameter W for a given cluster with a cluster center at r by averaging the angular velocity of the circling squirmers in the cluster: We only take clusters with a cluster size N cl > 2 into account. In Fig. 8 d) we show this quantity averaged over all clusters and 2.5 · 10 5 timesteps, normalized by the inverse ballistic time scale v 0 /R. Clearly, at χ = 0, the squirmers do not swirl. For non-zero χ the normalized swirling parameter jumps to a value with magnitude around 0.15v 0 /R, and afterwards fluctuates in a narrow region between 0.15v 0 /R and 0.2v 0 /R. For |χ| = 3.3, squirmers group together for short times, even though no long-term stable clusters emerge. Therefore we can still measure a value of W . Since the chirality of the squirmer system changes with the sign of χ, we can induce clock-wise motion with χ < 0, corresponding to W < 0. We show this case in video M8, where χ = −1.

Influence of gravity
The parameter α = v 0 /v sed determines how strongly squirmers are confined to the bottom wall. On the one hand, gravity induces the stokeslet flow via which squirmers also interact with each other hydrodynamically. Thus, at smaller α one expects squirmers to be more tilted against the vertical due to the stronger flow vorticity. On the other hand, at stronger gravity squirmers are closer to the bottom wall and there- fore the hydrodynamic interaction with the no-slip wall also forces a stronger upright orientation in addition to the gravitational torque, see Eq. (13) [21,67]. Thus, the behavior of the system is determined by these two effects.
In Fig. 9 we plot the mean cluster size N cl versus α. It increases monotonically for N = 32 and N = 64. Thus, clusters become smaller at higher gravity and finally disintegrate into single squirmers. In order to find the reason for this behavior, we measured the orientational distributions and show the mean vertical orientation in the inset. The general trend is that for decreasing α (increasing gravity), the squirmer orientation is more aligned with the vertical. Thus, the stronger vertical alignment of neutral squirmers due to their hydrodynamic interactions with the wall is able to counterbalance the action of the stronger stokeslet vorticities. This results in more compact and smaller clusters, similar to the effect of larger gravitational torques in Sect. 3.2.1.
For the system with N = 110 squirmers the mean cluster size seems to increase again for the lowest value α = 0.3. However, due to the large error bars this increase is not secured. A possible reason could be that due to the stronger confinement to the bottom wall and the larger areal density, the squirmers are in close contact to each other which increases the measured cluster size.
We note that as the clusters become smaller, they also start to rotate as a whole. In particular, we observe rotating squirmer trimers at large gravity. We mention them since they resemble the spinner states found at very high gravity in lattice-Boltzmann simulations [104]. In these studies the vertical alignment of neutral squirmes is solely due to their hydrodynamic interactions with the wall [see Eq. (13)]. Furthermore, in the limit of very large gravity (α → 0), our squirmer system reaches the regime of Ref. [47], where at zero external torque a quasi-hexagonal monolayer (Wigner fluid) of squirmers with upright orientation is formed.

Conclusions
We have investigated the cluster formation of bottomheavy squirmers under strong gravity floating above the bottom surface. Already a squirmer pair mimicking the minuet dance, found in experiments with the Volvox algae, reveals the basic gyrotactic mechanism [54,105]. While the stokeslet vorticity from the squirmer neighbor tilts the squirmer towards the neighbor so that they swim towards each other, the gravitational torque stabilizes the upright orientation. This enables horizontal oscillations of the squirmer positions.
At higher squirmer numbers, the vorticity-induced mutual attraction leads to the formation of squirmer clusters, however only for sufficiently large gravitational torques. We quantify the cluster properties by the twodimensional density distribution, the mean cluster size, and the radial distribution function, from which we extract the mean cluster distance. Our simulations find that horizontal motion within the clusters decreases with increasing torque and the clusters become more compact. The cluster size has a maximum at a finite torque that is close but below the value where the vorticity-induced rotation of touching neighbors can be balanced. The maximum standard deviation of the cluster size at the same torque value shows the volatility of the clusters in the emerging cluster state. Increasing the torque further, squirmers in the clusters are less mobile. The clusters become more compact and ultimately the cluster size reaches a constant value for each density. The mean distance between the clusters also decreases with increasing torque and saturates. Interestingly, the areal density has little effect on the mean cluster distance.
Furthermore, we investigated how force and rotlet dipoles influence cluster formation. Puller squirmers require larger torques in order to form clusters, because their additional flow-field vorticity adds to the effective attraction and thus the stabilizing gravitational torque has to be larger. Their clusters are situated closer to the bottom wall where hydrodynamic wall interactions tilt squirmers towards each other. In contrast, pushers do not form noticeable clusters in our simulations. Squirmers with a rotlet dipole also form clusters for small dipolar strengths and due to their interaction with the noslip wall perform swirling motion. Finally, for increasing gravitational strength squirmers move closer to the bottom wall. The hydrodynamic wall interactions contribute to the vertical alignment of neutral squirmers, which also decreases the cluster size. In conclusion, the swimmer hydrodynamics has a profound effect on the cluster formation when the occurring hydrodynamic multipoles besides the stokeslet also contribute to the flow vorticity, which has to be balanced by the gravitational torque.
Gyrotactic cluster formation could be realized for microswimmers, such as Janus particles or Volvox algae, which are indeed bottom-heavy [45,54]. Janus particles also offer the possibility to tune the strength of the gravitational torque, while for Volvox algae it exhibits a natural variation. For an L-shaped particle the torque was controlled by changing the relative arm lengths and thus the effective lever arm [33]. However, the realization of the external torque is not decisive for our results. For example, phototaxis is an alternative way of inducing orientational alignment and the strength of the vertical reorientation can be controlled by light intensity [106,107]. But here, the size and orientation of the part of the swimmer, which is illuminated, also influences the torque, as was shown for a phototactic Janus particle [108]. Microswimmers casting shadows on each other, influence their collective dynamics. Furthermore, gradients in oxygen concentration also orient bacteria along the vertical, which was observed for B. subtilis swimming upwards towards an air-liquid interface at the top [59]. However, in general, any chemotaxis towards chemical fields, such as nutrients, can profoundly affect the dynamics of microswimmers as exemplified by the effect of trail avoidance in active emulsion droplets [109]. These droplets also form floating, spontaneously rotating clusters, when the concentration of the surfactant fuel is high enough [97]. Thus, while we have identified the generic hydrodynamic mechanism of gyrotactic cluster formation, the influence of chemical fields needs further study. Finally, we mention two possible extensions of our investigations. Biological microswimmers typically move in non-Newtonian fluids [110,111]. Therefore, it is worthwhile to investigate, how robust gyrotactic cluster formation is for squirmers in complex fluids. This requires to extent the method of MPCD to viscoelastic media [80]. Furthermore, current research also explores the dynamics of non-spherical swimmers [112][113][114], including squirmer rods [86,115]. The collective motion of rod-shaped particles is strongly influenced by their steric interactions. Under a gravitational torque, this suggests that the vertical orientation is even more stable in a dense suspension with smectic order, because collisions with neighboring rods counteract tilting against the vertical. Therefore, the interesting question arises how this affects the formation of gyrotactic clusters.

Author contribution statement
All the authors were involved in the preparation of the manuscript. All the authors have read and approved the final manuscript.

Declarations
Data Availability Statement This manuscript has data included as electronic supplementary material. The online version of this article contains supplementary material, which is available to authorized users.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.