Emergent collective dynamics of bottom-heavy squirmers under gravity

We present the results of hydrodynamic simulations using the method of multi-particle collision dynamics for a system of squirmer microswimmers moving under the influence of gravity at low Reynolds numbers. In addition, the squirmers are bottom-heavy so that they experience a torque which aligns them along the vertical. The squirmers interact hydrodynamically by the flow fields of a stokeslet and rotlet, which are initiated by the acting gravitational force and torque, respectively, and by their own flow fields. By varying the ratio of swimming to bulk sedimentation velocity and the torque, we determine state diagrams for the emergent collective dynamics of neutral squirmers as well as strong pushers and pullers. For low swimming velocity and torque we observe conventional sedimentation, while the sedimentation profile becomes inverted when their values are increased. For neutral squirmers we discover convective rolls of circulating squirmers between both sedimentation states, which sit at the bottom of the system and are fed by plumes made of collectively sinking squirmers. At larger torques porous clusters occur that spawn single squirmers. The two latter states can also occur transiently starting from a uniform squirmer distribution and then disappear in the long-time limit. For strong pushers and pullers only weak plume formation is observed.

Studying microswimmers under gravity is important because often they are not neutrally buoyant [24][25][26][27]. In such a setting non-equilibrium sedimentation has been observed [24,28,29] accompanied by polar order along Contribution to the Topical Issue "Motile Active Matter", edited by Gerhard Gompper, Clemens Bechinger, Roland G. Winkler, Holger Stark. Supplementary material in the form of 12 .mp4 files available from the Journal web page at https://doi.org/10.1140/epje/i2020-11949-8 a e-mail: ruehle@tu-berlin.de the vertical [30] and convection [31]. Numerical hydrodynamic studies also discovered two-dimensional Wigner fluids and swarming under strong gravity [32], as well as fluid pumps in a parabolic potential [33]. In experimental systems swimming under the influence of external fields generates intriguing and surprising phenomena such as the formation of thin layers of motile phytoplankton in coastal regions [34], algae in bound dancing states [35], and hovering rafts of active emulsion droplets [36]. Many microswimmers also perform gravitaxis, which is the ability to align (anti-)parallel to the direction of gravity. The gravitational torque to achieve this alignment can either result from hydrodynamic drag of a microswimmer with asymmetric morphology [37][38][39] or from bottom heaviness, i.e., when the center of mass is offset relative to the geometrical center [34,[40][41][42]. Gravitactic swimmers can induce an overturning instability when accumulating with higher density at the top boundary (reminiscent of the Rayleigh-Taylor instability), which then initiates various patterns of bioconvection [43][44][45][46]. However, also gyrotaxis clearly plays an important role in such settings [46][47][48][49][50][51]. There, the collective dynamics of microswimmers depends on the combined action of gravity and hydrodynamic flow [40,46]. The involvement of physiological aspects in biotic pattern formation has also been discussed [52][53][54].
In theory microswimmer systems under gravity have been investigated in the past using the versatile spherical squirmer model swimmer [27,31,32,[55][56][57]. Here, we simulate around 900 bottom-heavy squirmers under gravity with full hydrodynamics using the method of multiparticle collision dynamics (MPCD) [58,59]. Varying the ratio of swimming to bulk sedimentation velocity and the gravitational torque due to bottom heaviness, we determine state diagrams for neutral as well as strong pusher and puller squirmers. While for low swimming velocity and torque conventional sedimentation is recovered, the sedimentation profile becomes inverted when increasing their values. For neutral squirmers we discover a rich phenomenology between both states including a state where plumes consisting of collectively sinking squirmers feed convective rolls at the bottom of the system and dense clusters which spawn single squirmers. These two states can also occur transiently starting from a uniform squirmer distribution and then disappear in the long-time limit. For strong pushers and pullers only weak plume formation is observed. We thoroughly characterize all states by different quantities.
In the following, in sect. 2 we introduce the squirmer model swimmer and the simulation method of multi-particle collision dynamics. Then, in sect. 3 we present all our results. We start with the state diagram of neutral squirmers followed by a detailed characterization of the different states and also look at strong pushers and pullers. Finally, we finish with a summary and conclusions.

Squirmer model swimmer and simulation method 2.1 Spherical squirmer
Swimming on the micron scale is dominated by friction [17,60]. Hence, hydrodynamics is captured by the Stokes equations where u and p are the fluid velocity and pressure fields, respectively. Biological microswimmers often propel themselves by collective beating patterns of cilia, which create flow fields along their cell surfaces [3,61]. Also artificial microswimmers exist that either use phoretic self-propulsion mechanisms, such as diffusiophoresis and thermophoresis in the case of active colloids [3,62,63], or Marangoni stresses in the case of active emulsion droplets [11] in order to generate such surface flow fields. A simple and effective approximation to all these swimming mechanisms is offered by the spherical squirmer model [64,65], where an axisymmetric tangential flow field on the surface is prescribed u(r)| r=R = ∞ n=1 B n 2P n (e ·r) n(n + 1) [−e + (e ·r)r] .
Here, e is the swimmer orientation vector, R is the swimmer radius, P n is the nth Legendre polynomial, and P n means its first derivative. Typically, the expansion is truncated after the second term, leaving the two relevant modes B 1 and B 2 . Then, the flow field generated by the surface field of eq. (3) in the surrounding fluid is [65,66] where β = B 2 /B 1 is the squirmer-type parameter.

Free squirmer
The squirmer induces a hydrodynamic source dipole and for β = 0 also a force dipole, the far fields of which decay as 1/r 3 and 1/r 2 , respectively. Swimmers with β = 0 are called neutral squirmers, while β > 0 generates pullers and β < 0 pushers. Since free squirmers are force-free, a stokeslet term with a flow field decaying as 1/r is not allowed but apears in eq. (4) [64,65]. The reason is that for a moving squirmer also the swimming velocity v 0 e contributes to its surface velocity field, which is not included in eq. (3). Thus, following Pak and Lauga [66] the flow field of eq. (4) has to be interpreted as the pumping field generated by a squirmer held at a constant position by a force. This stalling force F a is given by the balance equation [66]: where here v is the swimming velocity of the freely moving squirmer. One can calculate the stalling force using Lamb's solution to the Stokes equations [66][67][68], F a = 4πηR∇(e · r)B 1 = 4πηRB 1 e, and arrive at the known relation v = 2 3 B 1 e with the swimming speed v 0 := 2 3 B 1 . In a freely translating squirmer, the Stokes flow field initiated by the pumping force is no longer present. Thus, in eq. (4) the stokeslet vanishes and the source-dipole term is modified leading to the flow field of a free squirmer [65,66],

Squirmer under gravity
Adding the gravitational force −mge z modifies the force balance of eq. (5) and yields for the total squirmer velocity, As in our previous publications [31,55] we introduce the velocity ratio α := v 0 /v sed to compare the self-propulsion to the bulk sedimentation velocity, v sed = mg/6πηR.
The gravitational force adds a Stokes flow field to the free-squirmer solution of eq. (6), which as usual contains a stokeslet and a source-dipole contribution: where we introduced the coordinate along the vertical z = r · e z . Due to their long-range nature it is important to always take stokeslet flow fields into account when they occur. This has been shown in experimental studies of the dancing motion of Volvox algae [35] or the Stokesian dynamics of swimmers in a harmonic trap [33].

Squirmer with bottom heaviness
In this article we assume the spherical squirmer to be bottom-heavy, i.e., its center of mass has an offset r 0 from the geometrical center [42] such that a torque mgr 0 (−e z × e) acts on the swimmer. Balancing external torque and rotational friction torque −8πηR 3 Ω, we find the angular velocity We will later use the dimensionless parameter to quantify the strength of the external torque. It compares -up to the factor 3/4-the characteristic time scale of self-propulsion, R/v 0 , to the characteristic time of reorientation by bottom heaviness, 8πηR 3 /(mgr 0 ). Sometimes it is called the gyrotactic orientation parameter [46,69]. The rotating squirmer generates the flow field of a rotlet [68], Note that the spatial decay of this flow field is the same as for pusher and puller squirmers but is more long-ranged compared to the neutral squirmer. However, u bh r vanishes when the squirmer is aligned with the vertical, meaning e = e z .

Multi-particle collision dynamics
In the following we present the algorithm for performing simulations with multi-particle collision dynamics that we have already used in the past [22,31,32]. Therefore, we only summarize it here. The algorithm is implemented in a massively parallelized code, which runs on a computer cluster. In addition, we also provide the parameters used in our simulations.

Algorithm
We numerically solve the Navier-Stokes equations at low Reynolds number using the mesoscale method of multiparticle collision dynamics (MPCD) [58,59,70,71]. Thermal noise is automatically included in this particle-based solver. Since the Reynolds numbers employed in our simulations are smaller than one, we effectively obtain solutions of the Stokes equations where any inertia is neglected.
In the MPCD method the fluid is composed of point particles of mass m 0 that are kept at temperature T 0 . A simulation step consists of the fluid particles performing consecutive streaming and collision steps. During the streaming step, which has a duration Δt, each fluid particle i moves ballistically with its velocity v i according to Δt is a simulation parameter that controls the fluid viscosity [70,72]. During the streaming step fluid momentum is advected and also transferred to swimmers or absorbed by walls. Furthermore, boundary conditions need to be applied to the squirmer surfaces and to bounding walls. For this we employ the bounce-back rule [73] to implement either the surface slip velocity field of a squirmer from eq. (3) or the no-slip boundary condition for walls. The dynamics of the squirmers themselves is also computed during the streaming step. We perform 20 molecular dynamics steps during each streaming step using Velocity-Verlet integration together with the gravitational force and steric interactions between squirmers [71].
The purpose of the collision step is to exchange momentum between fluid particles. To that end, the simulation box is divided into cubical cells of edge length a 0 . We also use this length to define the fluid particle density as the average particle number n fl per collision cell. Within each cell fluid-particle velocities are updated with the help of a collision operator, for which we use the MPC-AT+a rule [59,71]. Thus, a thermostat is set up and both linear and angular momentums are conserved [59]. We also apply a grid shift for each new collision step in order to enforce Galilean invariance [74]. During the collision step fluid and squirmers/bounding walls interact as well: if a collision cell overlaps with a boundary or a squirmer, the overlapping volume is filled with virtual fluid particles to ensure the fluid density n fl remains the same [71,75]. After the collision step the momentum gain of the virtual fluid particles is transferred to the involved squirmer where the virtual particles are located.
The flow fields calculated with the MPCD method are accurate on length scales larger than the mean free path of the fluid particles. Using large enough squirmer radii, the hydrodynamic behaviour of squirmer microswimmers is therefore well reproduced by the MPCD method [76,77]. Thus it has widely been used to simulate a variety of settings [22,31,55,[77][78][79].

Parameters
For the most part, we use the parameters presented below, in case of deviations we have stated them in the text. We use the duration Δt = 0.02a 0 m 0 /k B T 0 for the streaming step and a fluid particle density of n fl = 10. This implies a viscosity of η = 16.05 √ m 0 k B T 0 /a 2 0 [72,80]. For our squirmers we use a radius R = 4a 0 , therefore the translational and rotational thermal diffusivities in bulk respectively. We choose B 1 = 0.1 k B T 0 /m 0 and thus have the active Péclet number Pe = Rv 0 /D T = 330, which is comparable in order of magnitude to bacterial swimmers, such as E. Coli [81] or B. subtilis [82], with Péclet numbers around 100. Furthermore, the ballistic time scale of the squirmer is R/v 0 = 3000Δt. The Reynolds number Re = v 0 Rn fl /η = 0.17 is relatively high for creeping flow systems. Lowering it further would mean considerably additional computational cost. However, viscous effects are definitely dominant in our simulations.
Typically, the volume density of squirmers was held constant at 10%, simulating 914 squirmers in a box with a system size of 108a 0 × 108a 0 × 210a 0 , where the latter length is the box height. This means that in a crosssectional slab of width 2R perpendicular to the direction of gravity, the mean area fraction of squirmers is 15%. We use no-slip boundary conditions at the top and bottom walls and periodic boundary conditions in the horizontal plane. As initial condition we always choose uniformly distributed squirmers such that no squirmers overlap.
For our study of bottom-heavy squirmers we vary the velocity ratio α by changing the gravitational acceleration g, which in an experiment depends on the density mismatch between swimmer and fluid. We find that increasing α from zero up to around 7 captures the relevant features in our squirmer simulations. Values in this range can be obtained experimentally for example with Volvox algae in water [83] or active emulsion droplets in a mixture of H 2 O and D 2 O [36]. The rescaled torque r 0 /Rα is varied by changing both g and the center-of-mass offset r 0 . Note that the torque depends on real mass rather than buoyant mass. In order to account for this difference, in ref. [42] the parameter r 0 was redefined as r 0 m/Δm, with the buoyant mass Δm. This way we can use the same parameter α in both the torque and force equations. Alternatively, one could also assume that Δm ≈ m [42].
To analyze the data of our simulation runs, we save the position, orientation, as well as translational and angular velocities for each squirmer every 1000th time step of duration Δt.
Note for most densely packed squirmers a depletion of the MPCD fluid particles is observed between the squirmers due to finite compressibility [23]. In the states discussed in sects. 3.2.2 and 3.4 the squirmers are not most densely packed. Nevertheless, we ran one simulation for each of the states at higher fluid particle density n fl = 80 and decreased self-propulsion velocity B 1 = 0.01. This leaves the Péclet number almost constant [23] while lower-ing compressibility. In these simulations we could confirm that the phenomenology stays the same as for the lower fluid particle density, which we typically use.

Results
In our hydrodynamic simulations we explored the dynamics of bottom-heavy squirmers under gravity. Depending on the velocity ratio α := v 0 /v sed and the strength of the gravitational torque, we observed a variety of stable and transient states. In fig. 1(a) we show the state diagram in the parameter space α versus reduced gravitational torque r 0 /Rα, which we determined in simulations with neutral squirmers. We have marked four exemplary states by colored points, for which we show the density profiles in fig. 2 in the same colors and videos MOESM1-MOESM4 in the Electronic Supplementary Material (ESM). In the following, we introduce the main phenomenology of the observed states.
The horizontal dotted line at α = 1 in fig. 1(a) marks the upper limit, below which isolated squirmers sink down in a bulk fluid and settle at a finite distance from a lower bounding wall [55]. Hydrodynamically interacting squirmers show a non-equilibrium sedimentation state, which persists to α ≈ 3 for weak torques. For zero torque the non-equilibrium sedimentation was already observed in ref. [31]. A typical sedimentation profile at α = 1.5 and r 0 /Rα = 0.01 is depicted in fig. 2. In contrast, at large α and torques, the orientational bias of the squirmers leads to their enrichment at the top wall and inverted sedimentation occurs, which is even observable for small torques. The corresponding inverted sedimentation profile for α = 6.01 and r 0 /Rα = 0.02 is shown in fig. 2. For dilute systems of bottom-heavy active particles such states were already described in ref. [42]. Between sedimentation and inverted sedimentation interesting dynamic states occur, which we shortly introduce now.
In the region colored in gray in fig. 1(a) we observe collections of sinking squirmers, which we call plumes. Although oriented upwards on average, they can sink due to the reduced viscous friction of a squirmer cluster. The plumes supply a convective roll at the bottom of the system that is formed and kept running by the self-propelling squirmers. Solitary squirmer escape from the edges of the roll and swim upwards. The formation of rolls and plumes are reminiscent of bioconvection observed in experiments [46,82,84,85]. We show an example of this state at α = 2.31 and r 0 /Rα = 0.11 in video MOESM3 and also provide a snapshot in fig. 1(b), left. The density profile is non-monotonous with a broad maximum at the position of the bottom cluster, a minimum in the central region, where plumes pass through, and a sharp maximum at the top wall, where squirmers accumulate.
Notably, starting from an initially uniform distribution of squirmers we also observe plumes that form at the top wall, sink down, and then slowly evaporate. Likewise, convection rolls can be merely transient when the bottom cluster eventually disappears. The steady state for these  The profiles belong to conventional non-equilibrium sedimentation (red), inverted sedimentation (blue), stable plumes and convection rolls (orange), and the spawning cluster (green).
cases are inverted sedimentation profiles with strong layering at the top wall. The transient plumes and rolls occur at higher torques in the blue region of the state diagram of fig. 1 beyond the dashed line and to the right of the state of stable plumes and convection rolls. An interesting situation arises in the state diagram when α is situated in a narrow stripe above α = 1 for torques larger than a threshold value that we show by the lower solid line in fig. 1(a). The clearest representation of this spawning-cluster state arises for large torques, where the squirmer orientation is fixed to the upright direction. A big cluster of squirmers floats above the lower wall (see fig. 1(b), right). Hydrodynamic interactions between the squirmers increase their mobilities and thereby their sedimentation velocities, which can cancel the swimming ve-locity even for α ≥ 1. To be more specific, a squirmer is pulled downward by the stokeslet flow fields from surrounding squirmers (cf. eq. (8)). Taking 12 of them (hexagonal packing) at distances equal to the mean spacing in the spawning cluster gives a velocity roughly twice the swimming velocity v 0 . Hence, squirmers in a cluster sink downward until they reach and interact with the bottom wall. Figure 2 shows the sedimentation profile for α = 1.5 and r 0 /Rα = 0.5. In comparison to the convection roll the bottom cluster has a higher density visible by a more pronounced broad maximum and also the depletion in the middle of the cell is stronger by an order of magnitude. We call this state "spawning cluster" because individual squirmers occasionally escape from the pores within the cluster at high velocity. This can be seen in video MOESM4, as well as in the snapshot in fig. 1 In table 1 we summarize the characteristics of the different states of neutral squirmers, including the figures which describe them. In the following sections we discuss these states in more detail. We investigate conventional and inverted sedimentation of neutral squirmers in sect. 3.1 and the state with plumes and convection rolls in sect. 3.2. Transient plumes and rolls in the inverted sedimentation state are discussed in sect. 3.3 and in sect. 3.4 we address the spawning-cluster state. Finally, we show how the state diagram changes for pusher and puller squirmers in sect. 3.5.

Sedimentation
Collective sedimentation of squirmers under gravity has extensively been studied in refs. [27,31,32]. It also occurs, of course, for bottom-heavy squirmers. In the state diagram of fig. 1(a) the sedimentation regime at low torques in the middle extends beyond α = 1, although single squirmers with upright orientation can overcome gravity. For α = 1.5 and r 0 /Rα = 0.01 we already showed the exponential sedimentation profile in fig. 2. It occurs because the flow fields of nearby squirmers tilt the orientation of one squirmer away from the upright direction, which therefore sinks even for α > 1. The necessary flow vorticity ω = curl u/2 for the reorientation is provided only by the gravity-induced stokeslets (see eq. (8)) since the flow field of neutral squirmers has zero vorticity. The situation changes for large torques, where the squirmer orientation is always upright. At α < 1 squirmers are confined to clusters sitting on the bottom wall. When crossing α ≈ 1 they start to float, in the sense that the density at z = 0 develops a minimum. An exemplary density profile of these spawning clusters for α = 1.5 and r 0 /Rα = 0.5 was already discussed in connection with fig. 2.

Inverted sedimentation
In fig. 3 we show a set of density profiles for the inverted sedimentation state for α = 6.01. While at zero torque the profile is nearly uniform, increasing the torque from zero the sedimentation length of the inverted profile decreases and at high torques the inversion becomes strong enough that layers of squirmers form at the top wall. Note that for the three largest torque values transient plumes are observed. In concrete, a cluster separates from the top layers, sinks as a plume, and slowly evaporates (see video MOSEM8 in the ESM, and sect. 3.3). However, the sedimentation profiles are determined after reaching the steady state in the long-time limit.

Plumes and convective rolls
In the following we discuss squirmer plumes that constantly appear in the bulk and feed a convective roll at the bottom of the system.

Collective sinking and plumes
Plumes of squirmers sink with preferentially upright orientation. To characterize this motional state and contrast it with (inverted) sedimentation, we plot in fig. 4 the mean vertical squirmer orientation cos ϑ xy as a function of height z for all squirmers, which drift downwards: v z < 0. The average is taken over the horizontal xy plane. In addition, for each α the dashed line indicates the threshold value cos ϑ th = 1/α for the degree of upright orientation, which a single squirmer must not exceed in order to sink. In the conventional sedimentation state (red curve) the mean upright orientation is always below the threshold value, which explains the downward drifting of the squirmers. The same is true for the inverted sedimentation state (blue curve). However, while the threshold value cos ϑ th belongs to a small upright orientation, the almost constant mean orientation is negative here, thus squirmers with velocity v z < 0 swim downwards (rather than sink). In contrast, the mean orientation of squirmers in the plume state (orange curve) exceeds cos ϑ th in the region away from the walls, where plumes occur. Thus, their sinking cannot be explained by looking at single squirmers. Instead, it occurs since hydrodynamic friction in clusters of squirmers is reduced such that the mobility of each squirmer is increased. Thus, their sedimentation velocity, which acts against the upwards swimming, is larger compared to single squirmers and the whole cluster can sink. Indeed, for two squirmers the leading order flow field acting on the neighbor is given by the stokeslet contribution in eq. (8), which provides a flow in negative z-direction and thus reinforces the gravitational sinking. This hydrodynamically induced mobility increase has already been studied for passive colloids on the basis of Rotne-Prager mobilities, for example, in refs. [86,87] and for larger conglomerates in refs. [88,89].
Cluster velocities. To quantify the collective sinking in the plume state further, we show in fig. 5(a) the distributions of the vertical cluster velocity for different cumulated cluster sizes N cl . To determine them, every 1000th time step we monitor the clustering in the system by grouping squirmers that have a neighbor distance d < R/8 into the same cluster. This identifies clusters of different size N cl , where solitary squirmers have N cl = 1. We only consider squirmers in a region with height from 60a 0 < z < 160a 0 in order to avoid the densely accumulated squirmers close to the top and bottom walls. The distribution of cluster sizes (excluding solitary swimmers) is shown in the inset of fig. 5(a). Having grouped the squirmers into clusters of different sizes, we can determine the velocity distributions for all squirmers, solitary squirmers (N cl = 1), for squirmers in clusters with N cl ≥ 2, as well as those with N cl ≥ 10 (see fig. 5(a)). As references we have indicated three characteristic velocities by vertical lines: the maximum bulk velocity v = v 0 − mg/γ ∞ (right black line) and the zero mean velocity of all squirmers (left black line). The dashed orange line is the vertical velocity, which is calculated with the mean vertical orientation of all squirmers from the orange distribution with N cl ≥ 10. The total distribution of vertical squirmer velocities in the convective plume state has a broad shape symmetric about a mean very close to zero. This makes sense since in steady state there should not be a non-zero vertical flux of squirmers. Solitary squirmers (N cl = 1) have a broad distribution as well, but with a positive bias. This illustrates that solitary squirmers contribute more to upwards than to downwards motion. Interestingly, the velocities of solitary squirmers and small clusters can exceed the free-swimming limit v z = v 0 . We speculate that this comes from squirmers approaching the convective roll at the bottom, where they are pushed up by hydrodynamic flow fields originating from the convective rolls. Indeed, squirmers with v z > v 0 are hardly present in the velocity distribution of an inverted sedimentation state (not shown) where convective flows at the bottom do not exist.
Restricting ourselves to clusters of larger size, N cl ≥ k, the distributions are more and more shifted to negative velocities with increasing k and become narrower compared to both the single-particle and the overall distributions. They also lose the tail with positive drift velocity v z . The squirmer velocities in clusters are determined by the balance of self-propulsion and sedimentation velocities. The latter is increased in clusters of squirmers compared to solitary squirmers due to their hydrodynamic interacions, as discussed above. In contrast, when one calculates for N cl ≥ 10 the mean vertical orientation and balances the upwards swimming and sedimentation velocities using the single-squirmer mobility, one obtains a small positive velocity indicated by the dashed orange line. It hardly intersects the distribution at its far end. However, this velocity is close to the mean of the single-squirmer velocity distribution (red curve), because we find that the orientational distributions of squirmers swimming alone or in clusters are very similar. We attribute this to the missing vorticity in the flow field of a neutral squirmer. Finally, for the system displayed in fig. 5(a) at N cl ≥ 10 the mean value of the vertical velocity is close to −v 0 . It decreases further for even larger clusters and, of course, also depends on the velocity ratio α.
In fig. 5(b) the strong decrease of hydodynamic friction in clusters with increasing size is clearly visible. We estimated the effective friction coefficient γ eff = −mg/v eff by determining an effective sedimentation velocity v eff of a squirmer within a cluster. For this we averaged over the sedimentation velocities v z − v 0 cos ϑ of all squirmers within a cluster and subsequently took the mean over clusters with the same size that occur within a time window of 10 6 Δt. We normalize the resulting effective friction by the bulk value. In addition to the strong size dependence of γ eff , we realize how already for a single squirmer hydrodynamic interactions with its neighbors reduce the friction coefficient compared to γ ∞ of an isolated squirmer.
Flow vorticity. We have already mentioned how advection by hydrodynamic flow fields from neighboring squirmers (stokeslets in leading order) enhances the sinking velocity and thus the mobility of squirmers in a plume. However, squirmers in the neighborhood of a sinking plume are also reoriented by the vorticity of these flow fields. While sinking, a squirmer reorients and therefore swims towards the plume. It joins the plume and thereby contributes to its vertically extended shape. We note that such vorticities also play a role in the formation of a fluid pump by hydrodynamically interacting active particles moving in a harmonic trap potential [33]. In that study hydrodynamic torques only compete with rotational noise, whereas in the present case they have to balance the external gravitational torque acting on bottom-heavy squirmers.
At high external torques stable plume states do not exist in the state diagram of fig. 1(a), because bottom heaviness completely dominates the upright squirmer orientations. A reorientation by flow vorticity is not possible. However, when both hydrodynamic and gravitational torques are comparable, a neighboring squirmer tilts towards a plume but can only join it for a vertical alignment with cos ϑ 1/α, because otherwise it will swim upwards as explained before. This is why plumes persist at higher torques for decreasing α in the state diagram of fig. 1(a). The mechanism reported here resembles an instability described previously [46], where reorientations by flow vorticity also induce the formation of plumes composed of gyrotactic algae.
To illustrate the balance of stokeslet vorticity and rotations induced by bottom heaviness, we consider two squirmers at the same height with a distance r 12 . Taking the curl of eq. (8) and setting it equal to eq. (10) leads to The second squirmer is tilted towards the first one, then moves towards it so that they can form a plume. For a simple estimate in parameter space we choose the minimal distance r 12 = 2R and the orientation θ = π/2, where the bottom-heavy torque is maximal. The hyperbolic curve obtained from this estimate is shown in fig. 1(a) and locates the region of the plume and convective roll state fairly well.
Sizes of sinking clusters. We have already identified the plumes in fig. 4 via the mean vertical orientation of sinking squirmers. The clustering dynamics offers a further means to characterize and distinguish the stable plume state from the other states. In fig. 6 we show the mean number of squirmers N − in sinking clusters, i.e., where N − ≥ 2, over a period of 10 6 MPCD time steps Δt and within the same vertical region 60 < z/a 0 < 160 as in fig. 5. We can clearly recognize the plume state (orange line) by the high spikes that correspond to sudden events of collectively sinking squirmers passing through the region. Also shown is the transient plume state (purple line) that we discuss further in sect. 3.3. Here, the spikes disappear around 1.3·10 6 Δt when the plume has evaporated. In contrast, in the inverted sedimentation state (blue line), the average size of sinking clusters remains low. Some small spikes can be seen for the conventional sedimentation state (red line), where small clusters form and induce some convective dynamics, which has been observed before [31]. Likewise, some small spikes are visible for the spawning cluster state (green line) but much more rarely.

Convective roll
In the context of microswimmers, convection rolls are stationary rotational patterns, which are formed by gravitactic swimmers due to their self-propulsion and advection in the self-generated flow fields [36,50,82]. Here, we use this term for the recirculating motion inside a cluster sitting at the bottom wall. The cluster is visible in the left snapshot of fig. 1(b) and in video MOESM3. In video MOESM5, we show a 3D view of the system in the region 0 ≤ z/a 0 ≤ 100, where we have color-coded the vertical squirmer velocities. A snapshot from the video is depicted in fig. 7(a). A very dynamic situation is visible. Squirmers in a plume sink down on the right and join the cluster. Because of their negative velocities they are colored in blue (see color bar in fig. 7(c)). The plume is also visible on the left due to the periodic boundary condition. At the same time, individual squirmers colored in red swim upwards in a region with low density. In fig. 7(b) we see this depleted region more clearly from the top. All squirmers below z = 70a 0 are plotted. Furthermore, we recognize that the convective roll has an elongated shape. This will change for larger systems as we demonstrate below. The circulation pattern in the convective roll is visualized in fig. 7(c). Here, we plot all the trajectories of squirmers within a slice of thickness Δy = R spanning roughly the time interval 9 · 10 5 Δt. Each velocity vector of a squirmer is colored according to the value of its horizontal velocity component v x . We recognize two wavy patterns in the bottom region below z/a 0 = 100. They originate from squirmers moving down-and upwards while swimming to the left (blue lines) and from squirmers moving up-and downwards while swimming to the right (red lines). This generates the two circular patterns of the rolls. One is clearly visible in the right region, a second one on the left is not complete due to the periodic boundary condition. In the middle we see squirmers leaving the convective rolls.
In order to study the vertical motion in the system, we plot in fig. 8 the mean vertical current density of the squirmers j z (x, y), defined via j z (x, y) := ρ(r) v z (r) [31], where ρ is the squirmer density, . . . means average along the z-direction, and · · · average over time. Figure 8 shows the current density with the vertical average taken in three different regions: at the botton where the convective rolls are (0 ≤ z/a 0 ≤ 80), in the middle where the plumes predominately occur (80 ≤ z/a 0 ≤ 120), and in the bulk region above it (120 ≤ z/a 0 ≤ 160). At the bottom (left plot) we immediately recognize the low density region where single squirmers move up (red region), while the two counterrotating rolls meet where the squirmers drift downwards (blue region). In the middle Fig. 8. Heat map of the mean vertical current density jz(x, y) in the horizontal plane and for three regions z/a 0 ∈ [0, 80], [80,120], and [120, 160]. The time range for averaging is the same as in fig. 7. We have applied a low-pass filter (provided by python's scipy package) in order to smoothen the data. section (center plot) the sinking squirmers are spatially focussed due to the plume formation: squirmers swim towards each other and sink collectively. The plumes feed the two counterrotating rolls and, therefore, are found at the right edge. In contrast, rising squirmers move mainly individually, which is why their distribution is more spread out (center and right plot). In the middle height section (center plot) a weak plume is observable close to the center of the plane. In the videos MOESM3 and MOESM6 we see how such plumes move to the side while sinking. This exemplifies the strong hydrodynamic flows that uphold the convective roll.
While the two convective rolls convey the picture of a regular structure, the path of a single squirmer inside the dense cluster is irregular, as we show in fig. 9. We have reoriented the point of view in comparison to fig. 7(a) by 90 • , such that we now look at the long side of the rolls. Furthermore, taking into account the periodic boundary condition, we have slightly shifted the system box so that the squirmer does not leave the box. For special points in time we also display the squirmer's orientation vector in the figure. The trajectory starts at the blue dot. After joining the rolls at position I, the squirmer meanders inside the dense cluster and eventually escapes from an edge (IV) with an upright orientation. The plotted trajectory ends at the red dot. The squirmer attempts to leave the bottom cluster several times (II, III) but is unsuccesful because the orientation is tilted too strongly against the vertical. It even crosses the low density region close to position IV as the projected trajectory on the horizontal plane shows. The strong variation of the squirmer orientation against the normal is, of course, due to the vorticity of the generated hydrdynamic flow field. In our case, it results dominantly from the stokeslet due to the gravitational force acting on each squirmer: The source-dipole far-field of the neutral squirmer has zero vorticity, while wall-induced image fields decay with O(r −4 ) and therefore are weak [55,90]. Interestingly, the "minuet" dance of a pair of Volvox algae in the experiments of ref. [35] could also be reproduced by stokeslet interactions between the two microswimmers. The more chaotic meandering motion observed here inside the convective rolls then occurs due to the interactions with all the surrounding squirmers.
We already mentioned that the elongated convective rolls occupy the whole horizonzal plane of our simulation box. To avoid this finite-size effect, we double the edge length of the square cross section, keep the height constant, but reduce the squirmer volume fraction by a factor of two to 5%. As the top view of fig. 10 (left) demonstrates, several compact convective rolls with an island shape appear in contrast to a single elongated cluster. In the side view of fig. 10 (right) we recognize several sinking plumes. Videos MOESM6 and MOESM7 illustrates impressively the observed dynamics, in particular, how sinking plumes feed the convective rolls. Interestingly, rolls in the larger system have a characteristic distance. Increasing the density to 10% as in the small systems the islands become larger and partially touch each other, but the distance of their centers stays approximately the same. Such  regular structures are known in biological systems, where length scales can vary with system height or density, depending on the swimming mechanisms and vertical alignment [48,50]. However, a systematic variation of system size or density is beyond the scope of this article.
We clearly see the island shape of the convective rolls also in the mean vertical current density illustrated in fig. 11. The sinking plumes above the rolls and the sinking squirmers inside them are visible by the blue areas, while around these regions squirmers swim upwards. Figure 12 completes the picture of the convective roll. We show squirmer trajectories inside horizontal slices of thickness 4R either at the bottom (left) or top (right) of the convective rolls. The direction of the in-plane velocity is color-coded. At the bottom squirmers move radially outwards (away from the dense clusters) while at the top they move radially inwards (towards the centers of the dense clusters). Combined with the vertical squirmer current, this implies a toroidal flow pattern for the convective rolls, which is also visible in video MOESM7.

Transient plumes and rolls
Plumes and convective rolls as described in sect. 3.2 are not stable at high torques. Instead, the long-term steady state of the system is an inverted exponential sedimentation profile that develops after a long-lived transient. In the following, we distinguish between two different transient states, which we observed starting from an initially uniform distribution of squirmers: evaporating plumes forming at the top wall and unstable rolls.

Evaporating plumes
First, we consider systems with large α 5.6. In the state diagram of fig. 1(a) we have identified inverted sedimentation with transient plumes to the right of the dashed line. At such high velocity ratios α stable plumes cannot exist since the vorticity from the gravitational stokeslet is too weak to orient neighboring squirmers towards each other in order to form stable plumes while sinking. Instead, squirmers escape to the top wall already for external torques of medium strength. Hence, starting from a uniform distribution dense squirmer layers start to form at the top wall. Here, due to flow vorticity the orientations tilt and protrusions of squirmers form that eventually separate from the layers. The sinking squirmer cluster reaches a final height well above the bottom wall. It emits squirmers which join the layers at the top wall. Thus the plumes gradually evaporate and disappear. The whole process can be seen in video MOESM8. In steady state inverted sedimentation profiles with a few layers occur. We have already plotted some of them in fig. 3 for r 0 /Rα ≥ 0.04 and also show such a profile at the end of video MOESM8. The inability of a sinking cluster to attract more squirmers to support itself shows how gyrotactic structure formation fails for squirmers at high α. As a consequence the steady state here is always an inverted sedimentation for any torque. Note, gravitational detachment of protrusions from a top layer has been extensively studied in different theoretical and experimental settings in connection with bioconvection [43,44,46,82,85].

Transient convective rolls
We now consider lower values of α, where sinking plumes and convective rolls form, starting from the uniform initial distribution. Now, an increasing external torque dominates the orientational dynamics of squirmers meaning they experience a stronger vertical bias. This counteracts the hydrodynamic reorientation and the motion of squirmers towards a plume. Thus, at higher torques the plumes become thinner and eventually disappear in the long-time limit. As a result, the convective rolls at the bottom wall do not receive sufficient influx of squirmers and also evaporate. Again, the system reaches a steady state with an inverted sedimentation profile and with layering at the top wall. The process is visualized in video MOESM9. Both of the observed transient structures reveal themselves with clear signatures in the time evolution of the spatial squirmer distribution. In fig. 13 we plot the volume fraction η versus time, for both the lower and upper half of the system (z ∈ [0, H/2] and [H/2, H], where H is the system height). The curves for the steady states of inverted and conventional sedimentation (red and blue curves), as well as for stable plumes and convective rolls (orange curve) fluctuate around a constant value. However, transient plumes sinking from the top wall (purple curve) or transient convective rolls (green curve) have a steadily decreasing density in the lower half of the system. At the same time, the squirmer layers at the top wall grow at the expense of the shrinking plumes and rolls. Simulating the transients is very time-consuming. For example, the green curve belongs to a transient convective roll and has not equilibrated yet. Initially, the density is high in the lower half of the system and the roll takes a long time to dissolve.

Spawning clusters and transient hovering
In fig. 14 we show a spawning cluster viewed from the side (left) and from the top (right). The snapshots belong to a spawning cluster at α = 1.5 and large torque situated at the far right of the state diagram in fig. 1. We  Fig. 15. Heat map of the mean vertical current density jz(x, y) in the horizontal plane at α = 1.5 in the region 80 ≤ z/a 0 ≤ 120 and averaged over a time period of 9 · 10 5 Δt. The torque value r 0/Rα corresponds to sedimentation (left) and a spawning cluster (middle and right). We have applied a low-pass filter (provided by python's scipy package) in order to smoothen the data. also show this system in video M4. In the top view of the right snapshot we observe a porous structure of the cluster. In contrast to the convective rolls discussed above, holes strongly depleted by squirmers are visible. Thus the clusters are not compact objects. We colored the squirmers in the snapshots according to their vertical velocity v z using the same color code as in fig. 7. Hence, the green color of most squirmers shows that they move little in the vertical direction. Single squirmers perform a random walk or meander around within the cluster and when they reach the edge of a hole, the flow field of the neighboring squirmers strongly drifts them upwards with large velocities up to 3v 0 . They either rejoin the cluster or leave it. This is nicely visible in video MOESM4.
In the spawning-cluster state plumes and convective rolls no longer exist, as videos MOESM4 and MOESM7 demonstrate and when inspecting spawning clusters at different parameters. This becomes also clear from fig. 15, where we show the mean vertical current densities in the region z/a 0 ∈ [80, 120] above the spwaning clusters and for three different torques at α = 1.5.
The left plot for r 0 /Rα = 0.08 represents the sedimentation state. Similar to ref. [31] convective patterns in the region with an exponential sedimentation profile are visible, which consist of clearly separated areas with upward and downwards moving squirmers. Increasing the rescaled torque to 0.17 (see video MOESM10), these areas start to disintegrate and thereby mark the onset of the spawning-cluster state. Increasing the torque even further to 0.50, the mean vertical current density is zero everywhere except for some small patches. This is consistent with the very small density in the bulk region of our system as demonstrated by the corresponding density profile in fig. 2. Note, however, that compared to the sedimentation state the spawning-cluster state has a much higher density at the top wall, where squirmers leaving the cluster gather.
In the state diagram of fig. 1 we also mention "transient hovering". For very large external torques, well above the values used in the state diagram, we expect spawning clusters to dissolve with time at α > 1. However, close to α = 1 they dissolve very slowly.

Influence of squirmer type: pushers and pullers
For strong pullers and pushers the state diagrams, which we show in fig. 16, simplify considerably compared to neutral squirmers shown in fig. 1(a). First of all, we do not observe any stable or transient convective rolls nor spawning clusters. The main features in the state diagrams of fig. 16 are conventional and inverted sedimentation, where the separation line is shifted to higher torques and higher alpha compared to neutral squirmers. By visual inspection we observe plumes (see videos MOESM11 and MOESM12), where clusters of squirmers form in the upper region of the simulation cell, sink down, and dissolve. For pullers the plumes are more pronounced as we explain in sect. 3.5.2 and we indicate them by the grey shaded region as part of the sedimentation state.

Influence of hydrodynamics on sedimentation state
For pushers and pullers, which have a non-zero squirmer parameter β, the flow field of a force dipole is added to the total squirmer velocity field u(r) as documented by eq. (4). It decays like 1/r 2 and in contrast to the pure sourcedipole field of neutral squirmers possesses a non-zero vorticity, which acts on the orientation of nearby squirmers. Previous studies already investigated the consequences of this vorticity field for suspensions of microswimmers. They found that it weakens polar order in clusters, whereby pullers often retain a higher degree of polar order than pushers [33,[91][92][93]. Indeed, in the inset of fig. 17 we observe for the same parameters α = 6.01 and r 0 /Rα = 0.04 that the neutral squirmer has a stronger alignment along the vertical than pushers and pullers. As a result, pushers and pullers have a stronger tendency to sink under gravity. Therefore, while neutral squirmers show inverted sedimentation at α = 6.01 for all torque values, strong pullers (β = 5) exhibit a transition from conventional to inverted sedimentation with increasing r 0 /Rα. Examples for the respective sedimentation profiles are plotted in the main graph of fig. 17 together with a uniform profile right at the transition. Thus, the disturbance of the vertical alignment by the vorticity field of strong pushers and pullers shifts the transition line between conventional and inverted sedimentation to larger torques and swimming velocities. For the same reason the tendency to form layers at the upper wall is strongly reduced. This becomes obvious by comparing the sedimentation profiles for the same parameters α = 6.01 and r 0 /Rα = 0.04 for neutral squirmers in fig. 3 and strong pullers (β = 5) in fig. 17.

Plumes of pullers and pushers
Unlike for neutral squirmers we do not observe stable convective rolls for both strong pullers and pushers. This becomes obvious from the density profiles in fig. 18. The broad and high density peak near the bottom wall is missing for pullers and pushers. Note that we increased the rescaled torque r 0 /Rα for pullers and pushers to be clearly in the region where plumes occur. Convective rolls possess dense squirmer clusters with some polar order due to bottom heaviness. However, such clusters are hydrodynamically unstable for squirmers with a strong force-dipole contribution [91,92]. The same applies to spawning clusters, which we also do not observe. For neutral squirmers we observed plumes feeding convective rolls. For strong pullers and pushers we can also identify plumes by visual inspection, as demonstrated in videos MOESM11 and MOESM12, respectively. Pullers form visible plumes in the upper region that disband close to the bottom wall. Pusher plumes are very unstable and disintegrate already while they are sinking. Generally, the plume clusters are smaller compared to neutral squirmers. Thus, when plotting the mean number of sinking squirmers in a cluster, N − , versus time as in fig. 6, the plumes cannot clearly be identified by high spikes. Instead, for pullers we determined the mean squirmer number by also averaging over more than 5 · 10 5 time steps. The result is plotted in fig. 19 versus the rescaled torque for different swimming speeds α. For α = 1.5, 2.0, and 3.0 we roughly see a sigmoidal shape and locate a transition to a plume state at the inflection point. This is the meaning Fig. 19. Mean size of sinking puller clusters (β = 5) as a function of the rescaled torque for different α. The mean N− is calculated in the lower region 20 ≤ z/a 0 ≤ 120 and by averaging over more than 5 · 10 5 time steps. Furthermore, all sinking clusters with N − ≥ 2 are considered.
of the curved dashed line in the schematic state diagram of fig. 16, left. For pushers we do not see such a sigmoidal shape but rather a slow and steady increase. For this reason we did not identify a separate region for plumes but only indicate that we see them around the line separating conventional and inverted sedimentation from each other at higher torques.

Conclusions and outlook
In this article we have investigated the remarkable features of microswimmer suspensions under gravity using full hydrodynamic simulations of ca. 900 squirmer model swimmers, where we concentrated on the neutral squirmer but also looked at strong pushers and pullers. We have determined the respective state diagrams varying the ratio of swimming to bulk sedimentation velocity and the gravitational torque due to bottom heaviness. The general trend in all three cases reveals conventional sedimentation for low swimming velocity and torque, while the sedimentation profile becomes inverted when increasing both values.
In addition, for neutral squirmers we have discovered a rich phenomenology in between both sedimentation states. Squirmers sink collectively in plumes due to reduced hydrodynamic friction and feed fascinating convective roll patterns of elongated or toroidal shape that self-organize at the bottom of the system. The plume formation is supported by squirmer reorientation due to vorticity resulting from the stokeslet contribution to the flow field. The latter is induced by the gravitational force acting on each squirmer. The combination of upward swimming (gravitaxis) and reorientation by nearby flow fields (rheotaxis) is called gyrotaxis, a mechanism introduced and discussed in connection with bioconvection [46][47][48][49]51]. In our case plumes and convective rolls also form transiently at larger torques. When starting from an initially uniform squirmer distribution, a dense layering of squirmers forms at the top wall. It develops an instability due to the gyrotactic mechanism and then sinking plumes emerge. At increasing torque and moderate speed ratios we also observe dense but porous squirmer clusters that float above the bottom wall and spawn single squirmers. They also become transient for increasing speed ratio. For strong pushers and pullers the transition line between conventional and inverted sedimentation is shifted to higher torques and speed ratios. The reason is the non-zero vorticity of the additional force-dipole flow field, which reorients neighboring squirmer orientations from the vertical so that they can sink more easily. This is also the reason why strong pusher and pullers do not show such a rich phenomenology. Only weak plume formation is observed without any stable convective rolls occurring. In the case of pullers we could quantify it by the mean size of sinking puller clusters.
In systems with biological microswimmers plume formation is often traced back to the gyrotactic mechanism introduced above [46,48,49]. However, also the overturning instability of a dense layer of microswimmers at the top boundary, reminiscent of the Rayleigh-Taylor instability is discussed [44,94,95]. Indeed, in ref. [95] for plumes of the microorganism Tetrahymena emanating from a dense top layer, gyrotaxis is discarded. In contrast, in our case plumes for the convective rolls also form in bulk, where clearly gyrotaxis is the relevant mechanism. Even for the transient plumes of neutral squirmers emanating from a dense squirmer layer at the top boundary, we think that gyrotaxis is predominant. However, these comments also suggest that gyrotaxis of squirmers needs to be studied further in the future, in particular, how the characteristic length scales of plumes and convective rolls depend on system size and squirmer density.
In biological systems, microswimmers are often of pusher and puller type but nevertheless, in contrast to our simulations, can form stable plumes and convection cells [48,50,82]. Our investigations are performed for strong pushers and pullers, while sufficiently weak pusher and puller squirmers with β closer to zero will also show convective rolls. Indeed, estimates for the force-dipole moment of some biological microswimmers show a weak dipole strength [83,96,97]. Furthermore, in our simulations we look at a very generic system concentrating on gravitational force, bottom heaviness, and hydrodynamic interactions between the squirmers. Real microswimmers with a non-spherical shape also experience a drag torque [39,54] and their flagella might act on neighbors by steric forces.
Also, for the algae C.reinhardtii it has been shown that the flow field induced by the periodic beating pattern varies in time and during a short period is reminiscent of a pusher [98,99]. Thus, real microorganisms show a large variability and it was argued that different organisms could even be distinguished via their bioconvection patterns [50]. Generalizing our simulations in these directions provides opportunities for interesting future research.
Bottom heaviness is not the only source for orientational order in order to observe interesting pattern formation under gravity. It can also be induced by additional external fields, in which the microorganism performs taxis. For example, the alga C.reinhardtii relies on phototaxis [85,100], while the bacterium B. subtilis shows aerotaxis, where it aligns along a gradient of oxygen [50,82]. The response to chemical fields (chemotaxis) is fascinating on its own [101][102][103][104][105][106] and combining it with microswimmers [107] moving under gravity opens a new research direction.
Finally, microorganisms moving under gravity might also adapt their behavior to further external cues. For example, it has been argued that phytoplankton actively change their morphology to adjust their migration strategy to turbulent flow fields [39]. This connects to another new and fascinating research direction related to learning in active systems [108][109][110].
Open Access Funding provided by Projekt DEAL. We are grateful for stimulating discussions with and valuable input from J.-T. Kuhr, R. Kapral, K. Drescher, and A.J.T.M. Mathijsen, and thank the referees for their suggestions. This project was funded by Deutsche Forschungsgemeinschaft through the priority program SPP1726 (grant number STA352/11). The authors acknowledge the North-German Supercomputing Alliance (HLRN) for providing HPC resources that have contributed to the research results reported in this paper.
Publisher's Note The EPJ Publishers remain neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Appendix A. Table 2 provides an index of the videos referenced in this paper and made available in the ESM. It contains the system parameters and describes the states visualized by the videos.

Author contribution statement
All the authors were involved in the preparation of the manuscript. All the authors have read and approved the final manuscript.
Open Access This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.