Open Access Article
This Open Access Article is licensed under a Creative Commons Attribution-Non Commercial 3.0 Unported Licence

Collective sedimentation of squirmers under gravity

Jan-Timm Kuhr*, Johannes Blaschke, Felix Rühle and Holger Stark*
Institut für Theoretische Physik, Technische Universität Berlin, Hardenbergstr. 36, 10623 Berlin, Germany. E-mail:;

Received 14th June 2017 , Accepted 20th August 2017

First published on 2nd October 2017

Active particles, which interact hydrodynamically, display a remarkable variety of emergent collective phenomena. We use squirmers to model spherical microswimmers and explore the collective behavior of thousands of them under the influence of strong gravity using the method of multi-particle collision dynamics for simulating fluid flow. The sedimentation profile depends on the ratio of swimming to sedimentation velocity as well as on the squirmer type. It shows closely packed squirmer layers at the bottom and a highly dynamic region with exponential density dependence towards the top. The mean vertical orientation of the squirmers strongly depends on height. For swimming velocities larger than the sedimentation velocity, squirmers show strong convection in the exponential region. We quantify the strength of convection and the extent of convection cells by the vertical current density and its current dipole, which are large for neutral squirmers as well as for weak pushers and pullers.

1 Introduction

Microswimmers, whether biological or artificially produced, propel themselves forward without the help of any external force.1 Often however, external fields act on active particles.2 Examples are microswimmers in shear,3–6 Poiseuille,4,7–10 or swirling flow,11 in harmonic traps,12–15 and in light fields.16,17 The most natural example is gravity, which affects every swimmer that is not neutrally buoyant. In dilute active suspensions, where hydrodynamic interactions between swimmers can be neglected, both experimental18,19 and theoretical studies20,21 find exponential density profiles similar to that of passive colloids, but with a sedimentation length δ, which depends on activity and well surpasses that of passive colloids. Interestingly, in these dilute suspensions, analytical studies show the emergence of polar order,21 and if the microswimmers are bottom heavy, the sedimentation profile can even be inverted.22

For higher densities of microswimmers hydrodynamic interactions become important and collective behavior emerges.2,14,23–29 This includes motility-induced phase separation,30–34 swarming,16,35 and bioconvection,36–39 to name but a few phenomena. Furthermore, in real settings interactions with bounding surfaces are important,31,34,40–42 especially if the swimmers are not perfectly buoyant.43

In this article we consider systems with thousands of microswimmers under the influence of gravity. We simulate their full hydrodynamic flow fields using the method of multi-particle collision dynamics (MPCD)44,45 in order to include hydrodynamic interactions between swimmers as well as between swimmers and bounding walls. As a model microswimmer we use the squirmer,46–49 which is versatile enough to model the relevant swimmer types including pushers, pullers, and neutral swimmers.

In the following we concentrate on the case, where passive colloids would just strongly sediment to the bottom. We show how density or sedimentation profiles depend on the ratio of active to sedimentation velocity as well as on the squirmer type. During collective sedimentation squirmers develop densely packed layers in the bottom region of the simulation cell. In contrast, we observe an exponential density profile in the upper region, where squirmers form a more dilute active suspension. The mean vertical orientation of the squirmers depends strongly on their vertical position as well as on their swimmer type. For swimming velocities larger than the sedimentation velocity, we find that hydrodynamic interactions organize squirmers into convection cells. Importantly, both the strength of convection and the extension of the convection cells depend on the squirmer type.

The article is organized as follows. We first introduce the squirmer as our model microswimmer and then shortly address the simulation method of MPCD along with parameter settings and some details of our analysis in Section 2. In Section 3 we present our results of collectively sedimenting squirmers and analyze especially sedimentation and mean vertical orientation (in Section 3.1) and convection (in Section 3.2). Finally, in Section 4 we summarize our findings and conclude.

2 Methods and model

2.1 The squirmer as a model swimmer

In this work we use the squirmer46–49 as a model for a spherical microswimmer with radius R. It propels itself forward using a slip velocity field on its surface,
vs(rs) = B1(1 + βê·[r with combining circumflex]s)[(ê·[r with combining circumflex]s)[r with combining circumflex]sê], (1)
which generates a hydrodynamic flow field in the surrounding fluid. Here, rs is a vector, which points from the center of the squirmer to a point on its surface, [r with combining circumflex]s = rs/R is the corresponding unit vector, and the unit vector ê indicates the squirmer orientation. In the bulk of a quiescent fluid the squirmer orientation coincides with its swimming direction. Our model in eqn (1) only takes into account the first two terms introduced in ref. 46 and 47. They are sufficient to determine the swimming speed and swimmer type, by which microorganisms and artificial microswimmers like Janus particles50,51 or active droplets52–57 are typically characterized. Thus, the squirmer propels along its orientation vector ê with swimming speed v0 = 2/3B1 and creates fluid flow, the far field of which is controlled by the parameter β. The value of β therefore indicates the squirmer type. While β = 0 creates a neutral squirmer with the far field of a source dipole (∼r−3), β < 0 and >0 refer to pushers or pullers, respectively, the flow fields of which decay like a force dipole (∼r−2).58

2.2 Multi-particle collision dynamics

We investigate the behavior of many squirmers, which interact hydrodynamically with each other and with confining surfaces. We employ multi-particle collision dynamics (MPCD)44,45,59–61 to numerically solve the Navier–Stokes equations including thermal noise. They reduce to the Stokes equations at the low Reynolds numbers studied in this article.

In our MPCD simulations the fluid is modeled by ca. 2 × 107 point particles of mass m0. Their positions ri are updated in a streaming step using their velocities vi: ri(t + Δt) = ri(t) + viΔt. In the subsequent collision step fluid particles within cubic cells of linear extension a0 exchange momentum according to the MPC-AT + a rule.60 This conserves linear and angular momentum and also thermalizes velocities to temperature T. Further details of our implementation are described in ref. 31 and 34. We use here the parallelized version of ref. 34, which is suited to simulate large systems with many swimmers.

In the present work, we consider squirmers under gravity. So we have to add an acceleration term aΔt2/2 to the squirmers’ position in the streaming step, where the acceleration a is due to the gravitational force − mgez along the vertical with m the buoyant mass of a squirmer and g the gravitational acceleration. Since gravity does not induce a noticeable density change of the fluid on the micron length scale, we do not apply a gravitational acceleration to the fluid particles. If a fluid particle encounters a bounding wall or a squirmer, the particle's position and velocity are updated according to the “bounce-back rule”,62 which implements either the no-slip boundary condition or the surface flow field of eqn (1), respectively. During the streaming step momentum is transferred from the fluid particles to the squirmers, the velocities of which are updated by a molecular dynamics step. This includes steric interactions among squirmers and with bounding walls.

MPCD reliably reproduces the analytical results, including the flow field around passive colloids,59 the friction coefficient of a particle approaching a plane wall,63 the active velocity of squirmers,49 as well as the torque acting on them close to walls, where lubrication theory has to be applied.41 It also simulates correctly segregation and velocity oscillations in dense colloidal suspensions under Poiseuille flow.64,65 MPCD resolves flow fields on time and length scales large compared to the duration of the streaming step Δt and the mean free path of the fluid particles, respectively. Therefore, using a squirmer radius of R = 4a0, we expect to resolve hydrodynamic flow fields even when squirmers are close to each other.

2.3 Parameters

We simulate the behavior of N = 2560 squirmers of radius R = 4a0 under gravity in a cuboidal box of extensions Lx = Ly = 112a0 in the horizontal plane and Lz = 224a0 along the vertical. Thus, the mean volume fraction of squirmers amounts to ϕ ≈ 0.244. At z = 0 and z = Lz our system is bounded by walls, while periodic boundary conditions apply in the horizontal directions. For the duration of the streaming step we choose image file: c7sm01180f-t1.tif, which sets the shear viscosity to image file: c7sm01180f-t2.tif.66

In the following, an important parameter will be the ratio of active to bulk sedimentation velocity,

image file: c7sm01180f-t3.tif(2)
For spherical squirmers with buoyant mass m and gravitational acceleration g one has vg = mg/(6πηR). We will keep v0 fixed (it is set by v0 = 2/3B1 with image file: c7sm01180f-t4.tif), and choose three values of mg in order to study the cases α = 0.3, 1.0, and 1.5. In addition, we mention the sedimentation lengths of passive Brownian particles with the same buoyant masses m, δ0 = kBT/(mg), where kBT is thermal energy. For the values of α given above we calculate from this formula the respective values δ0 = 9.3 × 10−4R, 3.1 × 10−3R, and 4.7 × 10−3R, so that without activity squirmers would settle into a dense packing at the bottom of the simulation cell.

The Péclet number Pe = v0R/D, where D = kBT/(6πηR) is the translational diffusion coefficient, has the value Pe = 323 in all simulations, thus thermal translational motion is negligible. Furthermore, in all our simulations we have a Reynolds number of Re = v0Rnfl/η = 0.17, where nfl = 10 is the average number of fluid particles per collision cell. This implies Stokesian hydrodynamics where inertia can be neglected. Finally, with the thermal rotational diffusion coefficient Dr = kBT/(8πηR3), we introduce the persistence number Per = v0/(DrR). It measures the distance in units of particle radius R, where the squirmer moves persistently in one direction, before rotational diffusion changes its orientation. For all our simulations we have Per = 430. Thus, without gravity a single isolated squirmer would swim across the vertical extent Lz of the simulation cell in an almost straight line.

2.4 Determining sedimentation lengths

To determine the sedimentation lengths of the squirmers, we need to be sure that the system is in a steady state. In the beginning of the simulations we initialize the squirmers with random positions and orientations. We then observe that the collection of squirmers “collapses” towards the bottom wall by monitoring the mean squirmer height 〈z〉. In continuing the simulations, we ensure that ultimately 〈z〉 does not show any deterministic trend, but is only subject to fluctuations. We then simulate for a period of at least 104 MPCD time units (i.e. image file: c7sm01180f-t5.tif) and use this simulation data for our further analysis.

As explained in the results section, we determine the sedimentation or density profile ρ(z) of the squirmers, from which we identify some layering at the bottom wall of the simulation box, which is followed by a transitional and then an exponential region. In the latter we determine the sedimentation length δ using an exponential fit. The difficulty is to specify a range of heights [zb, zt], in which the exponential fit is performed. We have developed heuristic but robust criteria to identify this range. They ensure that neither layering at the bottom wall nor accumulation of squirmers at the top wall influences the fit values for δ. As a first constraint we demand that zt is at least a distance of 10R away from the top wall. For zb we require twice the height, which the squirmers would assume if they were all perfectly stacked in a hexagonal close packing. Within this first specification for the range [zb, zt] we then determine the final zt as the height, where the density ρ is minimal. For zb we take the smallest height z, where ρ(z) falls below 8% of the hexagonal-close-packed density. The value of 8% is a purely empirical value.

We use the data in the range [zb, zt] to obtain the sedimentation length δ from exponential fitting. In order to also estimate its error, we need to generate several estimates for δ from our data. Therefore, we split up the simulation time after reaching steady state into ten intervals. For each interval we perform exponential fits in four different ranges: (i) [zb, zt], (ii) [zb + 0.1Δz, zt], (iii) [zb, zt − 0.1Δz], and (iv) [zb + 0.1Δz, zt − 0.1Δz], where Δzztzb. We use the modified ranges as an additional measure to ensure that we are in the exponential regime. As an estimate for δ, we then take the mean of all 40 fits, while the corresponding standard deviation specifies the error.

3 Results

After a transient, the squirmers under gravity settle into a steady state, which we analyze in the following. Fig. 1 shows a snapshot for a simulation, where the ratio of swimming to sedimentation velocity was α = 1.0. In the lower part layers of squirmers have formed. In particular, the lowest layers display clusters of hexagonal packing (see inset of Fig. 1), which gradually dissolves when moving upward. After a transition region, where layering is not recognizable anymore, a dilute region of squirmers follows, where we will identify the exponential density profile. In the ESI, the two Videos V1 and V2 (for α = 1.5) illustrate very impressively how dynamic the whole sedimentation profile is, especially in the exponential regime. This is in stark contrast to passive particles. In the following we will investigate some features of sedimenting squirmers in more detail.
image file: c7sm01180f-f1.tif
Fig. 1 Snapshot of 2560 neutral squirmers (β = 0) moving under gravity in a steady state at α = 1. The volume fraction is ϕ ≈ 0.244. Three regions can be distinguished: layering at the bottom, followed by a transitional regime, and finally a region with an exponential density profile and convective motion of squirmers. Inset: Hexagonal clustering in the lowest layer of squirmers.

3.1 Sedimentation profile and vertical alignment

The sedimentation profile in Fig. 2 quantifies the observations from Fig. 1. Close to the bottom, layering is clearly observable in the volume fraction or density ρ(z) and indicated by peaks, the height of which gradually decreases within ca. 11 layers. After a transitional regime, ρ(z) decays exponentially and then is influenced by the upper bounding wall. Passive Brownian particles with buoyant mass m show an exponential sedimentation profile with sedimentation length δ0 = kBT/mg, as previously introduced. For very dilute suspensions of active particles, one can still derive2,18,20–22,67 and observe18,19 an exponential profile, however, with an increased sedimentation length δ > δ0. Even if passive particles all sink to the bottom due to their weight (δ0R), active particles with sufficiently large swimming speed can rise from the bottom with a sedimentation length δ > R = 4. In our simulations, squirmers strongly interact hydrodynamically by the flow fields they generate. Nevertheless, we observe an exponential decay of the density ρ(z) similar to lattice-Boltzmann simulations12,68 and experiments,18,19 which we find non-trivial. Thus from fits to the exponential part of the density profile, ρ(z) ∼ ez/δ, we extract the sedimentation length δ.
image file: c7sm01180f-f2.tif
Fig. 2 Semi-logarithmic plot of the sedimentation profile ρ(z) for the system in Fig. 1 (α = 1.0 and β = 0). Different regions are indicated. The green line is an exponential fit to extract the sedimentation length δ.

In Fig. 3(a) we present sedimentation profiles for a larger swimming speed, α = 1.5, and different swimmer types β to explore the influence of pushers and pullers. At larger α, all profiles show an exponential regime with a larger sedimentation length compared to α = 1. Fig. 4 shows a parametric study of δ in units of squirmer radius R plotted versus β and for three values of α. Clearly, δ decreases with α and is only a fraction of R for α = 0.3, when the activity is too small for particles to swim upwards. However, already the ratio α = 1.5 is sufficient to have sedimentation lengths δ ≈ 10R. To address the robustness of our results, for the case β = 0 and α = 1.5, we reduced the volume density ϕ by decreasing the number of squirmers N to a value where layer formation does not occur, while keeping the height Lz of the simulation box fixed. This does not influence the sedimentation length δ significantly, as long as the squirmer density at the bottom of the box is similar to that of the top layers. For β = 0 we also reduced both N and Lz by a factor of two, which keeps ϕ constant. We find that the sedimentation length is reduced by about 30% and 40% for α = 1 and α = 1.5, respectively. This is not surprising, since hydrodynamic interactions with the top wall, which were not relevant before, push squirmers downwards and also turn them away from the wall,58,69 which makes them swim downwards.

image file: c7sm01180f-f3.tif
Fig. 3 Parameter study for varying β at α = 1.5. (a) Semi-logarithmic plot of the squirmer density ρ(z) as a function of height z. Vertical dashed lines in each plot indicate (from left to right) top of layering, followed by the start and end of the exponential regime, where the sedimentation length δ is extracted by an exponential fit (bold orange line). (b) Mean vertical orientation of squirmers as a function of height z. θ denotes the angle between the vertical and the squirmer orientation. (c) Vertical squirmer current density image file: c7sm01180f-t15.tif, averaged over the exponential regime and time, color-coded in the xy plane.

In Fig. 4 we realize that for weak pushers (β = −1) the sedimentation length is the largest and decreases for stronger pushers and also pullers. It has been reported in the literature that the interaction between parallel squirmers grows with |β|.48,70 We speculate that the collective interactions of many squirmers will therefore strongly randomize swimming directions and hinder squirmers with large |β| from reaching larger heights, as reflected in the sedimentation length. The larger sedimentation length for weak pushers as compared to neutral squirmers is, however, unexpected. A possible explanation comes from the shape of flow fields for β ≠ 0. Pushers, in their center-of-mass frame, have a stagnation point with vortical flow in front of them, while pullers have it at their back.23 Since squirmers are typically pointing up in the exponential regime (see below) and since their density ρ(z) decreases with height z, pullers reorient more nearby squirmers compared to pushers, which decreases δ. Interestingly, the trend of δ for varying β is inverted for small α with the minimum being at β ≈ 1. The reason for the inversion is not clear, but since δ < R, we assume that interactions with the densely packed squirmer layers are relevant.

image file: c7sm01180f-f4.tif
Fig. 4 Sedimentation length δ in units of R as a function of squirmer parameter β for different ratios of swimming to sedimentation velocities, α. Inset: Semi-logarithmic plot.

In Fig. 3(a) we also observe how the layer structure in the lower part of the system is influenced by β. For β = 0 and β = 1 the minima between successive layers are less pronounced, which implies less order. In contrast, especially for β = −1 and β = −2 layering is more pronounced. We ascribe this difference to the hydrodynamic interactions between neighboring squirmers, which depend on β.

Finally, in Fig. 5 we show the mean orientation of squirmers as a function of height z for the system illustrated in Fig. 1. The neutral squirmers in the bottom layer at z = 0 have an upright orientation due to hydrodynamic interactions with the bounding wall.31,41,42,48 The mean orientation then decreases to zero (see also Fig. 1) and drops to a negative value at the rim of the layering. This is simply because squirmers from above swim into the dense squirmer region and need some time to reorient and swim away. In the transitional region the orientation changes again rapidly to nearly upright and shows only small variations in the exponential regime.

image file: c7sm01180f-f5.tif
Fig. 5 Mean vertical orientation of squirmers, 〈cos[thin space (1/6-em)]θ〉, as a function of height z for the system in Fig. 1 (α = 1.0 and β = 0).

The occurrence of polar order in the sedimentation profile of dilute suspensions has been predicted by theory and already occurs without any interactions (i.e. for dilute suspensions of active Brownian particles) just for kinetic reasons.2,21 Hydrodynamic interactions between squirmers obviously do not destroy the polar order. In our parametric study of Fig. 3(c), this is also confirmed for other squirmer types β. Differences occur in the layering and in the transitional region. For pullers the upright orientation in the bottom layer decreases for larger β, as expected by hydrodynamic interactions with the bottom wall in lubrication theory.41,42,48 In the adjacent layers hardly any polar order is visible in contrast to neutral squirmers. Weak pushers (β = −1 and −2) show a similar but weaker trend compared to neutral squirmers. For strong pushers (β = −5), again, there is hardly any polar order in the layering.

3.2 Convection

As Videos V1 and V2 in the ESI demonstrate, the squirmers in the exponential density region are very mobile. In fact, while their mean vertical velocity is zero, the steady-state distribution of vertical velocities for β = 0 and α = 1.5 can be well fitted by a Gaussian with a standard deviation comparable to v0 (see Fig. 6). For large |β| small deviations from the Gaussian form occur. Thus, also vertical squirmer speeds larger than ⅓v0 arise (the maximal vertical velocity a single bulk squirmer can have at α = 1.5). This means that squirmers are advected by flow fields set up by their neighbors. Indeed, for neutral squirmers and weak pullers/pushers, we see evidence for convection flow extending over the whole simulation cell.
image file: c7sm01180f-f6.tif
Fig. 6 Distribution of vertical squirmer velocity (blue) and Gaussian fit (black) for α = 1.5 and β = 0 in the exponential regime.

To quantify convection, we take the vertical squirmer current density and average it along the vertical in the exponential regime:

jz(x) = 〈ρ(x)〈vz(x) (3)
where x = (x,y) is a position in the horizontal plane, 〈⋯〉 means average along the vertical, and ρ is the squirmer density. In the following we always indicate a time average over some quantity q in the steady state by [q with combining macron]. We plot image file: c7sm01180f-t6.tif in Fig. 7 for neutral squirmers and α = 1.5 in the xy plane of the simulation box. While in the lower left region squirmers move upward, the vertical current goes downward in the upper right indicating a convection cell, which extends over the whole horizontal plane. In Fig. 3(c) we present image file: c7sm01180f-t7.tif for different squirmer types at the same α. For weak pullers (β = 1 and 2) the extent of the convection cell decreases and at β = 5 large-scale convection is no longer observable. For weak pushers (β = −1 and −2) the current density becomes weaker but still large-scale convection is visible, which then vanishes for β = −5. For the case β = 0, α = 1.5 we checked that large-scale convection is stable against a reduction in the squirmer density ϕ at constant Lz and for reduced Lz while keeping ϕ constant. In both settings a single convection cell extends across the simulation box.

image file: c7sm01180f-f7.tif
Fig. 7 Squirmer current density, image file: c7sm01180f-t16.tif, averaged over the exponential regime and time, color-coded in the xy plane for α = 1.5 and β = 0.

To further characterize the squirmer density current jz(x), we calculate its zeroth and first moment. The zeroth moment is the volume average of the vertical current density:

image file: c7sm01180f-t8.tif(4)
where A is the area of the xy plane of the simulation cell and 〈⋯〉 means average over the whole volume of the exponential region. The vertical current density 〈jz〉 strongly fluctuates in time (see Fig. 8) but in the steady state its temporal mean, image file: c7sm01180f-t9.tif, has to be zero because of particle conservation. Thus, we use 〈|jz|〉 to quantify how mobile the squirmers are in the exponential region (see Fig. 8). In Fig. 9 we show the time average image file: c7sm01180f-t10.tif versus β for different rescaled swimming velocities α. As expected, at α > 1 the squirmers are more mobile than for α ≤ 1 since they are able to move against gravity. Furthermore, for α = 1.5 the mean vertical current density image file: c7sm01180f-t11.tif decreases for large |β|, revealing again the importance of the advective flow set up by the different squirmer types.

image file: c7sm01180f-f8.tif
Fig. 8 Volume average of the squirmer current density, 〈jz〉, and its magnitude, 〈|jz|〉, plotted versus time for α = 1.5 and β = 0.

image file: c7sm01180f-f9.tif
Fig. 9 Time and volume average of the magnitude of the current density, image file: c7sm01180f-t17.tif, as a function of squirmer parameter β for different rescaled swimming speeds α.

We call the first moment of jz(x) current dipole,

image file: c7sm01180f-t12.tif(5)
where we choose the center of the xy plane as the origin of x. The current dipole serves to quantify the strength and horizontal extension of the convection cell by its magnitude jD ≔ |jD|. The cell’s orientation relative to the x-axis is given by the angle φD with cos[thin space (1/6-em)]φD = jD·ex/jD. Fig. 10 shows how jD strongly fluctuates in time, reflecting again the high mobility of the squirmers. The orientation angle φD also fluctuates and in the example of Fig. 10 assumes two mean orientations around 1.5π and π. Overall, we can record that the spatial arrangement of convection is subject to strong fluctuations and is strongly variable in time.

image file: c7sm01180f-f10.tif
Fig. 10 Magnitude jD and orientation angle φD of the current dipole plotted versus time for α = 1.5 and β = 0.

Finally, in Fig. 11 we plot the time average image file: c7sm01180f-t13.tif versus squirmer type β for different α. Convection is largest for large α and neutral squirmers. The current dipole we define in eqn (5) can be interpreted in analogy with a charge dipole in electrostatics. Its magnitude changes when either the current density (the separated “charges”) changes in magnitude or when the distance between regions of positive vs. negative vertical speed is altered (the distance of “charge separation”). In accordance with Fig. 3(c) we find that image file: c7sm01180f-t14.tif for α = 1.5 decreases when the magnitude of the current density jz(x) decreases (weak pusher, β < 0) or the extension of the convection cell becomes smaller (weak puller, β > 0). Strong pushers/pullers with |β| = 5 only show weak convection. For small rescaled swimming speed α ≤ 1 convection is generally small.

image file: c7sm01180f-f11.tif
Fig. 11 The time average of the magnitude of the current dipole, image file: c7sm01180f-t18.tif, as a function of squirmer parameter β for different rescaled swimming speeds α.

4 Conclusions

In this work we addressed the collective sedimentation of squirmers in a gravitational field concentrating on the relevant case, where the corresponding passive particles would completely sediment. We showed that the sedimentation profile can be divided into three distinct regions; closely packed layers near the bottom, a transition region, and a rather dilute active suspension at the top. Here, the density depends exponentially on height, as for passive colloids, and the profile is very dynamic. The exponential dependence is non-trivial given the strong hydrodynamic interactions between the squirmers due to their flow fields. We also identified a strongly height-dependent mean orientation of the swimmers. While the mean orientation is strongly varying across the closely packed layers, in particular for small |β|, it varies far less in the exponential region. From the latter we extracted sedimentation lengths and showed that these not only grow with the ratio of active to sedimentation velocity, but also depend on the squirmer type. We argue that neutral squirmers or weak pushers and pullers are more persistent when swimming upwards and thereby show larger sedimentation lengths.

Furthermore, sedimenting squirmers create strong convective currents due to their hydrodynamic interactions. The spatial extension and the strength of convection are again determined by the rescaled swimming speed and by the squirmer type. Neutral squirmers as well as weak pushers and pullers show the strongest convectional flow. In particular, for swimming speeds larger than the sedimentation velocity pronounced convection cells occur, which extend over the whole simulation box. Finally, as another signature of the highly dynamic sedimentation profile in the exponential region, we identified strong temporal fluctuations of the convective currents.

What we could not resolve in our current simulations due to limited computational resources is the question of what determines the lateral extent of the convection cells. Is there an intrinsic length scale, which sets it? For this we would need to increase the simulation cell in the lateral directions. In future work, we plan to include bottom heaviness of the squirmers in our simulations and study in detail the inversion of the sedimentation profile, which was discussed in ref. 22 for very dilute swimmer suspensions. Preliminary results in denser systems show that instabilities occur due to hydrodynamic interactions.71 In a harmonic trapping potential a similar instability leads to the formation of fluid pumps by breaking the rotational symmetry of the trap.12,14 In the present case we expect convectional patterns to occur. Thereby, we will connect to the phenomenon of bioconvection.36–39

Conflicts of interest

There are no conflicts to declare.


We would like to thank C. Cottin-Bizonne and A. Zöttl for stimulating discussions. This project was funded by Deutsche Forschungsgemeinschaft through the research training group GRK 1558 and priority program SPP 1726 (grant number STA352/11). Simulations were conducted at the “Norddeutscher Verbund für Hoch- und Höchstleistungsrechnen” (HLRN), project number bep00050.


  1. E. Lauga and T. R. Powers, Rep. Prog. Phys., 2009, 72, 096601 Search PubMed.
  2. H. Stark, Eur. Phys. J.: Spec. Top., 2016, 225, 2369–2387 Search PubMed.
  3. A. Sokolov and I. S. Aranson, Phys. Rev. Lett., 2009, 103, 148101 Search PubMed.
  4. S. Rafa, L. Jibuti and P. Peyla, Phys. Rev. Lett., 2010, 104, 600 Search PubMed.
  5. H. M. López, J. Gachelin, C. Douarche, H. Auradou and E. Clément, Phys. Rev. Lett., 2015, 115, 121 Search PubMed.
  6. E. Clément, A. Lindner, C. Douarche and H. Auradou, Eur. Phys. J.: Spec. Top., 2016, 225, 2389–2406 Search PubMed.
  7. A. Zöttl and H. Stark, Phys. Rev. Lett., 2012, 108, 218104 Search PubMed.
  8. S. Uppaluri, N. Heddergott, E. Stellamanns, S. Herminghaus, A. Zöttl, H. Stark, M. Engstler and T. Pfohl, Biophys. J., 2012, 103, 1162–1169 Search PubMed.
  9. A. Zöttl and H. Stark, Eur. Phys. J. E: Soft Matter Biol. Phys., 2013, 36, 3 Search PubMed.
  10. R. Rusconi, J. S. Guasto and R. Stocker, Nat. Phys., 2014, 10, 212–217 Search PubMed.
  11. M. Tarama, A. M. Menzel and H. Löwen, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 90, 032907 Search PubMed.
  12. R. W. Nash, R. Adhikari, J. Tailleur and M. E. Cates, Phys. Rev. Lett., 2010, 104, 258101 Search PubMed.
  13. A. Pototsky and H. Stark, Europhys. Lett., 2012, 98, 50004 Search PubMed.
  14. M. Hennes, K. Wolff and H. Stark, Phys. Rev. Lett., 2014, 112, 238104 Search PubMed.
  15. A. M. Menzel, A. Saha, C. Hoell and H. Löwen, J. Chem. Phys., 2016, 144, 024115 Search PubMed.
  16. J. A. Cohen and R. Golestanian, Phys. Rev. Lett., 2014, 112, 365 Search PubMed.
  17. C. Lozano, B. ten Hagen, H. Löwen and C. Bechinger, Nat. Commun., 2016, 7, 12828 Search PubMed.
  18. J. Palacci, C. Cottin-Bizonne, C. Ybert and L. Bocquet, Phys. Rev. Lett., 2010, 105, 088304 Search PubMed.
  19. F. Ginot, I. Theurkauff, D. Levis, C. Ybert, L. Bocquet, L. Berthier and C. Cottin-Bizonne, Phys. Rev. X, 2015, 5, 011004 Search PubMed.
  20. J. Tailleur and M. E. Cates, Europhys. Lett., 2009, 86, 60002 Search PubMed.
  21. M. Enculescu and H. Stark, Phys. Rev. Lett., 2011, 107, 058301 Search PubMed.
  22. K. Wolff, A. M. Hahn and H. Stark, Eur. Phys. J. E: Soft Matter Biol. Phys., 2013, 36, 9858 Search PubMed.
  23. A. A. Evans, T. Ishikawa, T. Yamaguchi and E. Lauga, Phys. Fluids, 2011, 23, 111702 Search PubMed.
  24. M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao and R. A. Simha, Rev. Mod. Phys., 2013, 85, 1143–1189 Search PubMed.
  25. D. Saintillan and M. J. Shelley, C. R. Phys., 2013, 14, 497–517 Search PubMed.
  26. F. Alarcón and I. Pagonabarraga, J. Mol. Liq., 2013, 185, 56–61 Search PubMed.
  27. W. Yan and J. F. Brady, Soft Matter, 2015, 11, 6235–6244 Search PubMed.
  28. A. Zöttl and H. Stark, J. Phys.: Condens. Matter, 2016, 28, 253001 Search PubMed.
  29. C. Krüger, C. Bahr, S. Herminghaus and C. C. Maass, Eur. Phys. J. E: Soft Matter Biol. Phys., 2016, 39, 103 Search PubMed.
  30. M. E. Cates and J. Tailleur, Annu. Rev. Condens. Matter Phys., 2015, 6, 219–244 Search PubMed.
  31. A. Zöttl and H. Stark, Phys. Rev. Lett., 2014, 112, 118101 Search PubMed.
  32. R. Matas-Navarro, R. Golestanian, T. B. Liverpool and S. M. Fielding, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2014, 90, 032304 Search PubMed.
  33. R. M. Navarro and S. M. Fielding, Soft Matter, 2015, 11, 7525–7546 Search PubMed.
  34. J. Blaschke, M. Maurer, K. Menon, A. Zöttl and H. Stark, Soft Matter, 2016, 12, 9821–9831 Search PubMed.
  35. N. Oyama, J. J. Molina and R. Yamamoto, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2016, 93, 043114 Search PubMed.
  36. H. Wager, Philos. Trans. R. Soc., B, 1911, 201, 333 Search PubMed.
  37. T. J. Pedley and J. O. Kessler, Annu. Rev. Fluid Mech., 1992, 24, 313–358 Search PubMed.
  38. C. R. Williams and M. A. Bees, J. Exp. Biol., 2011, 214, 2398–2408 Search PubMed.
  39. C. R. Williams and M. A. Bees, J. Fluid Mech., 2011, 678, 41–86 Search PubMed.
  40. I. Llopis and I. Pagonabarraga, J. Non-Newtonian Fluid Mech., 2010, 165, 946–952 Search PubMed.
  41. K. Schaar, A. Zöttl and H. Stark, Phys. Rev. Lett., 2015, 115, 038101 Search PubMed.
  42. J. S. Lintuvuori, A. T. Brown, K. Stratford and D. Marenduzzo, Soft Matter, 2016, 12, 7959–7968 Search PubMed.
  43. F. Rühle, J. Blaschke, J.-T. Kuhr and H. Stark, New J. Phys., submitted Search PubMed.
  44. A. Malevanets and R. Kapral, J. Chem. Phys., 1999, 110, 8605 Search PubMed.
  45. G. Gompper, T. Ihle, D. M. Kroll and R. G. Winkler, Advances in polymer science 221, Springer, Berlin, 2009, pp. 1–87 Search PubMed.
  46. M. J. Lighthill, Comm. Pure Appl. Math., 1952, 5, 109–118 Search PubMed.
  47. J. R. Blake, J. Fluid Mech., 1971, 46, 199 Search PubMed.
  48. T. Ishikawa, M. P. Simmonds and T. J. Pedley, J. Fluid Mech., 2006, 568, 119–160 Search PubMed.
  49. M. T. Downton and H. Stark, J. Phys.: Condens. Matter, 2009, 21, 204101 Search PubMed.
  50. R. Golestanian, T. B. Liverpool and A. Ajdari, Phys. Rev. Lett., 2005, 94, 53 Search PubMed.
  51. S. Ebbens, D. A. Gregory, G. Dunderdale, J. R. Howse, Y. Ibrahim, T. B. Liverpool and R. Golestanian, Europhys. Lett., 2014, 106, 58003 Search PubMed.
  52. S. Thutupalli, R. Seemann and S. Herminghaus, New J. Phys., 2011, 13, 073021 Search PubMed.
  53. M. Schmitt and H. Stark, Europhys. Lett., 2013, 101, 44008 Search PubMed.
  54. M. Schmitt and H. Stark, Eur. Phys. J. E: Soft Matter Biol. Phys., 2016, 39, 062901 Search PubMed.
  55. M. Schmitt and H. Stark, Phys. Fluids, 2016, 28, 012106 Search PubMed.
  56. C. C. Maass, C. Krüger, S. Herminghaus and C. Bahr, Annu. Rev. Condens. Matter Phys., 2016, 7, 171–193 Search PubMed.
  57. C. Jin, C. Krüger and C. C. Maass, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, 5089–5094 Search PubMed.
  58. S. E. Spagnolie and E. Lauga, J. Fluid Mech., 2012, 700, 105–147 Search PubMed.
  59. J. T. Padding and A. A. Louis, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 74, 031402 Search PubMed.
  60. H. Noguchi, N. Kikuchi and G. Gompper, Europhys. Lett., 2007, 78, 10005 Search PubMed.
  61. R. Kapral, Advances in Chemical Physics 140, John Wiley & Sons, Hoboken, 2008, pp. 89–146 Search PubMed.
  62. J. T. Padding, A. Wysocki, H. Löwen and A. A. Louis, J. Phys.: Condens. Matter, 2005, 17, S3393 Search PubMed.
  63. J. T. Padding and W. J. Briels, J. Chem. Phys., 2010, 132, 054511 Search PubMed.
  64. P. Kanehl and H. Stark, J. Chem. Phys., 2015, 142, 214901 Search PubMed.
  65. P. Kanehl and H. Stark, Phys. Rev. Lett., 2017, 119, 018002 Search PubMed.
  66. H. Noguchi and G. Gompper, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 78, 016706 Search PubMed.
  67. J. Tailleur and M. E. Cates, Phys. Rev. Lett., 2008, 100, 218103 Search PubMed.
  68. M. E. Cates, K. Stratford, R. Adhikari, P. Stansell, J.-C. Desplat, I. Pagonabarraga and A. J. Wagner, J. Phys.: Condens. Matter, 2004, 16, S3903 Search PubMed.
  69. A. P. Berke, L. Turner, H. C. Berg and E. Lauga, Phys. Rev. Lett., 2008, 101, 23 Search PubMed.
  70. I. O. Götze and G. Gompper, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 82, 041921 Search PubMed.
  71. M. Hennes, K. Wolff and H. Stark, unpublished results.


Electronic supplementary information (ESI) available: 2 video files. See DOI: 10.1039/c7sm01180f

This journal is © The Royal Society of Chemistry 2017