Collective Sedimentation of Squirmers under Gravity

TL;DR

This study investigates how active spherical microswimmers called squirmers sediment under gravity, revealing layered structures, orientation patterns, and convection behaviors influenced by swimming speed and squirmer type.

Contribution

It introduces a detailed simulation of large-scale squirmer sedimentation under gravity, highlighting the effects of swimming velocity and squirmer type on collective behavior.

Findings

Layered sedimentation profiles with close-packed bottom layers.

Exponential density distribution in the upper regions.

Strong convection observed for swimming velocities exceeding sedimentation speed.

Abstract

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 close 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.

Equations

$${\mathbf{v}}_s({\mathbf{r}}_s)=B_1\left(1+\beta \hat{\mathbf{e}}\cdot {\hat{\mathbf{r}}}_s\right)\left[\left(\hat{\mathbf{e}}\cdot {\hat{\mathbf{r}}}_s\right){\hat{\mathbf{r}}}_s-\hat{\mathbf{e}}\right],$$

$$\alpha :=\frac{v_0}{v_g}.$$

$$j_z(\mathbf{x})=⟨\rho {⟩}_{\Vert }(\mathbf{x})⟨v_z{⟩}_{\Vert }(\mathbf{x})$$

$$⟨j_z⟩=\frac{1}{A}{\int }_A{j}_z(\mathbf{x})d^{2}x=⟨\rho ⟩⟨v_z⟩,$$

$${\mathbf{j}}_D=\frac{1}{A}{\int }_A\mathbf{x}{j}_z(\mathbf{x})d^{2}x,$$

Full text

I Introduction

Microswimmers, whether biological or artificially produced, propel themselves forward without the help of any external force Lauga and Powers (2009). Often however, external fields act on active particles Stark (2016). Examples are microswimmers in shear Sokolov and Aranson (2009); Rafaï et al. (2010); López et al. (2015); Clément et al. (2016), Poiseuille Rafaï et al. (2010); Zöttl and Stark (2012); Uppaluri et al. (2012); Zöttl and Stark (2013); Rusconi et al. (2014), or swirling flow Tarama et al. (2014), in harmonic traps Nash et al. (2010); Pototsky and Stark (2012); Hennes et al. (2014); Menzel et al. (2016), and in light fields Cohen and Golestanian (2014); Lozano et al. (2016).

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 experimental Palacci et al. (2010); Ginot et al. (2015) and theoretical studies Tailleur and Cates (2009); Enculescu and Stark (2011) find exponential density profiles similar to that of passive colloids, but with a sedimentation length $\delta$, which depends on activity and well surpasses that of passive colloids. Interestingly, in these dilute suspensions analytical studies show the emergence of polar order Enculescu and Stark (2011), and if the microswimmers are bottom heavy, the sedimentation profile can even be inverted Wolff et al. (2013).

For higher densities of microswimmers hydrodynamic interactions become important and collective behavior emerges Evans et al. (2011); Marchetti et al. (2013); Saintillan and Shelley (2013); Alarcón and Pagonabarraga (2013); Hennes et al. (2014); Yan and Brady (2015); Zöttl and Stark (2016); Stark (2016); Krüger et al. (2016). This includes motility-induced phase separation Cates and Tailleur (2015); Zöttl and Stark (2014); Matas-Navarro et al. (2014); Navarro and Fielding (2015); Blaschke et al. (2016), swarming Cohen and Golestanian (2014); Oyama et al. (2016), and bioconvection Wager (1911); Pedley and Kessler (1992); Williams and Bees (2011a, b), to name but a few phenomena. Furthermore, in real settings interactions with bounding surfaces are important Llopis and Pagonabarraga (2010); Schaar et al. (2015); Zöttl and Stark (2014); Lintuvuori et al. (2016); Blaschke et al. (2016), especially if the swimmers are not perfectly buoyant Rühle et al. .

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) Malevanets and Kapral (1999); Gompper et al. (2009) in order to include hydrodynamic interactions between swimmers as well as between swimmers and bounding walls. As a model microswimmer we use the squirmer Lighthill (1952); Blake (1971); Ishikawa et al. (2006); Downton and Stark (2009), 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 II. In Section III we present our results of collectively sedimenting squirmers and analyze especially sedimentation and mean vertical orientation (in Section III.1) and convection (in Section III.2). Finally, in Section IV we summarize our findings and conclude.

II Methods and Model

II.1 The squirmer as a model swimmer

In this work we use the squirmer Lighthill (1952); Blake (1971); Ishikawa et al. (2006); Downton and Stark (2009) as a model for a spherical microswimmer with radius $R$. It propels itself forward using a slip velocity field on its surface,

[TABLE]

which generates a hydrodynamic flow field in the surrounding fluid. Here, $\bm{r}_{s}$ is a vector, which points from the center of the squirmer to a point on its surface, $\hat{\bm{r}}_{s}=\bm{r}_{s}/R$ is the corresponding unit vector, and the unit vector $\hat{\bm{e}}$ indicates the squirmer orientation. In the bulk of a quiescent fluid the squirmer orientation coincides with its swimming direction. Our model in eq. (1) only takes into account the first two terms introduced in Lighthill (1952); Blake (1971). They are sufficient to determine the swimming speed and swimmer type, by which microorganisms and artificial microswimmers like Janus particles Golestanian et al. (2005); Ebbens et al. (2014) or active droplets Thutupalli et al. (2011); Schmitt and Stark (2013, 2016a, 2016b); Maass et al. (2016); Jin et al. (2017) are typically characterized. Thus, the squirmer propels along its orientation vector $\hat{\bm{e}}$ with swimming speed $v_{0}=2/3B_{1}$ and creates fluid flow, the far field of which is controlled by the parameter $\beta$. The value of $\beta$ therefore indicates the squirmer type. While $\beta=0$ creates a neutral squirmer with the far field of a source dipole ($\sim r^{-3}$), $\beta<0$ and $>0$ refer to pushers or pullers, respectively, the flow fields of which decay like a force dipole ($\sim r^{-2}$) Spagnolie and Lauga (2012).

II.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) Malevanets and Kapral (1999); Padding and Louis (2006); Noguchi et al. (2007); Kapral (2008); Gompper et al. (2009) to numerically solve the Navier-Stokes equations including thermal noise. They reduce to the Stokes equations at low Reynolds numbers studied in this article.

In our MPCD simulations the fluid is modeled by ca. $2\cdot 10^{7}$ point particles of mass $m_{0}$. Their positions $\bm{r}_{i}$ are updated in a streaming step using their velocities $\bm{v}_{i}$: $\bm{r}_{i}(t+\Delta t)=\bm{r}_{i}(t)+\bm{v}_{i}\Delta t$. In the subsequent collision step fluid particles within cubic cells of linear extension $a_{0}$ exchange momentum according to the MPC-AT+a rule Noguchi et al. (2007). It conserves linear and angular momentum and also thermalizes velocities to temperature $T$. Further details of our implementation are described in Refs. Zöttl and Stark (2014); Blaschke et al. (2016). We use here the parallelized version of Ref. Blaschke et al. (2016), 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 $\bm{a}\Delta t^{2}/2$ to the squirmers position in the streaming step, where the acceleration $\bm{a}$ is due to the gravitational force $-mg\bm{e}_{z}$ 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” Padding et al. (2005), which implements either the no-slip boundary condition or the surface flow field of eq. (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. It includes steric interactions among squirmers and with bounding walls.

MPCD reliably reproduces analytical results, including the flow field around passive colloids Padding and Louis (2006), the friction coefficient of a particle approaching a plane wall Padding and Briels (2010), the active velocity of squirmers Downton and Stark (2009), as well as the torque acting on them close to walls, where lubrication theory has to be applied Schaar et al. (2015). It also simulates correctly segregation and velocity oscillations in dense colloidal suspensions under Poiseuille flow Kanehl and Stark (2015, 2017).

MPCD resolves flow fields on time and length scales large compared to the duration of the streaming step $\Delta t$ and the mean free path of the fluid particles, respectively. Therefore, using a squirmer radius of $R=4a_{0}$, we expect to resolve hydrodynamic flow fields even when squirmers are close to each other.

II.3 Parameters

We simulate the behavior of $N=2560$ squirmers of radius $R=4a_{0}$ under gravity in a cuboidal box of extensions $L_{x}=L_{y}=112a_{0}$ in the horizontal plane and $L_{z}=224a_{0}$ along the vertical. Thus, the mean volume fraction of squirmers amounts to $\phi\approx 0.244$. At $z=0$ and $z=L_{z}$ our system is bounded by walls, while periodic boundary conditions apply in the horizontal directions. For the duration of the streaming step we choose $\Delta t=0.02a_{0}\sqrt{m_{0}/k_{B}T}$, which sets the shear viscosity to $\eta=16.05\sqrt{m_{0}k_{B}T}/a_{0}^{2}$ Noguchi and Gompper (2008).

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

[TABLE]

For spherical squirmers with buoyant mass $m$ and gravitational acceleration $g$ one has $v_{g}=mg/(6\pi\eta R)$. We will keep $v_{0}$ fixed (it is set by $v_{0}=2/3B_{1}$ with $B_{1}=0.1\sqrt{k_{B}T/m_{0}}$), and choose three values of $mg$ in order to study the cases $\alpha=0.3$, $1.0$, and $1.5$. In addition, we mention the sedimentation lengths of passive Brownian particles with the same buoyant masses $m$, $\delta_{0}=k_{B}T/(mg)$, where $k_{B}T$ is thermal energy. For the values of $\alpha$ given above we calculate from this formula the respective values $\delta_{0}=9.3\cdot 10^{-4}R$, $3.1\cdot 10^{-3}R$, and $4.7\cdot 10^{-3}R$, so that without activity squirmers would settle into a dense packing at the bottom of the simulation cell.

The Péclet number $\mathrm{Pe}=v_{0}R/D$, where $D=k_{B}T/(6\pi\eta R)$ is the translational diffusion coefficient, has the value $\mathrm{Pe}=323$ in all simulations, thus thermal translational motion is negligible. Furthermore, in all our simulations we have a Reynolds number of $\mathrm{Re}=v_{0}Rn_{\text{fl}}/\eta=0.17$, where $n_{\text{fl}}=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 $D_{r}=k_{B}T/(8\pi\eta R^{3})$, we introduce the persistence number $\mathrm{Pe}_{r}=v_{0}/(D_{r}R)$. It measures the distance in units of particle radius $R$, which the squirmer moves persistently in one direction, before rotational diffusion changes its orientation. For all our simulations we have $\mathrm{Pe}_{r}=430$. Thus, without gravity a single isolated squirmer would swim across the vertical extent $L_{z}$ of the simulation cell on an almost straight line.

II.4 Determining sedimentation lengths

To determine the sedimentation lengths of the squirmers, we need to be sure that the system is in 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 $\langle z\rangle$. In continuing the simulations, we ensure that ultimately $\langle z\rangle$ does not show any deterministic trend, but is only subject to fluctuations. We then simulate for a period of at least $10^{4}$ MPCD time units (i.e. $10^{4}a_{0}\sqrt{m_{0}/k_{B}T}$) and use this simulation data for our further analysis.

As explained in the results section, we determine the sedimentation or density profile $\rho(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 $\delta$ using an exponential fit. The difficulty is to specify a range of heights $[z_{b},z_{t}]$, 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 $\delta$. As a first constraint we demand that $z_{t}$ is at least a distance of $10R$ away from the top wall. For $z_{b}$ 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 $[z_{b},z_{t}]$ we then determine the final $z_{t}$ as the height, where the density $\rho$ is minimal. For $z_{b}$ we take the smallest height $z$, where $\rho(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 $[z_{b},z_{t}]$ to obtain the sedimentation length $\delta$ from exponential fitting. In order to also estimate its error, we need to generate several estimates for $\delta$ 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) $[z_{b},z_{t}]$, (ii) $[z_{b}+0.1\Delta z,z_{t}]$, (iii) $[z_{b},z_{t}-0.1\Delta z]$, and (iv) $[z_{b}+0.1\Delta z,z_{t}-0.1\Delta z]$, where $\Delta z:=z_{t}-z_{b}$. We use the modified ranges as an additional measure to ensure that we are in the exponential regime. As an estimate for $\delta$, we then take the mean of all 40 fits, while the corresponding standard deviation specifies the error.

III Results

After a transient, the squirmers under gravity settle into a steady state, which we analyze in the following. Figure 1 shows a snapshot for a simulation, where the ratio of swimming to sedimentation velocity was $\alpha=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 supplemental material†, the two videos V1 and V2 (for $\alpha=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.

III.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 $\rho(z)$ and indicated by peaks, the height of which gradually decreases within ca. 11 layers. After a transitional regime, $\rho(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 $\delta_{0}=k_{B}T/mg$, as previously introduced. For very dilute suspensions of active particles, one can still derive Tailleur and Cates (2008, 2009); Palacci et al. (2010); Enculescu and Stark (2011); Wolff et al. (2013); Stark (2016) and observe Palacci et al. (2010); Ginot et al. (2015) an exponential profile, however, with an increased sedimentation length $\delta>\delta_{0}$. Even if passive particles all sink to the bottom due to their weight ($\delta_{0}\ll R$), active particles with sufficiently large swimming speed can rise from the bottom with a sedimentation length $\delta>R=4$. In our simulations, squirmers strongly interact hydrodynamically by the flow fields they generate. Nevertheless, we observe an exponential decay of the density $\rho(z)$ similar to lattice-Boltzmann simulations Cates et al. (2004); Nash et al. (2010) and experiments Palacci et al. (2010); Ginot et al. (2015), which we find non-trivial. Thus from fits to the exponential part of the density profile, $\rho(z)\sim e^{-z/\delta}$, we extract the sedimentation length $\delta$.

In Fig. 3a) we present sedimentation profiles for a larger swimming speed, $\alpha=1.5$, and different swimmer types $\beta$ to explore the influence of pushers and pullers. At larger $\alpha$, all profiles show an exponential regime with larger sedimentation length compared to $\alpha=1$. Fig. 4 shows a parametric study of $\delta$ in units of squirmer radius $R$ plotted versus $\beta$ and for three values of $\alpha$. Clearly, $\delta$ decreases with $\alpha$ and is only a fraction of $R$ for $\alpha=0.3$, when activity is too small for particles to swim upwards. However, already the ratio $\alpha=1.5$ is sufficient to have sedimentation lengths $\delta\approx 10R$. To address the robustness of our results, for the case $\beta=0$ and $\alpha=1.5$, we reduced the volume density $\phi$ by decreasing the number of squirmers $N$ to a value where layer formation does not occur, while keeping the height $L_{z}$ of the simulation box fixed. This does not influence the sedimentation length $\delta$ significantly, as long as the squirmer density at the bottom of the box is similar to that of the top layers. For $\beta=0$ we also reduced both $N$ and $L_{z}$ by a factor of two, which keeps $\phi$ constant. We find that the sedimentation length is reduced by about 30% and 40% for $\alpha=1$ and $\alpha=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 Spagnolie and Lauga (2012); Berke et al. (2008), which makes them swim downwards.

In Fig. 4 we realize that for weak pushers ($\beta=-1$) the sedimentation length is largest and decreases for stronger pushers and also pullers. It has been reported in literature that the interaction between parallel squirmers grows with $|\beta|$ Ishikawa et al. (2006); Götze and Gompper (2010). We speculate that the collective interactions of many squirmers will therefore stronger randomize swimming directions and hinder squirmers with large $|\beta|$ 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 $\beta\neq 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 Evans et al. (2011). Since squirmers are typically pointing up in the exponential regime (see below) and since their density $\rho(z)$ decreases with height $z$, pullers reorient more nearby squirmers compared to pushers, which decreases $\delta$. Interestingly, the trend of $\delta$ for varying $\beta$ is inverted for small $\alpha$ with the minimum being at $\beta\approx 1$. The reason for the inversion is not clear, but since $\delta<R$, we assume that interactions with the densely packed squirmer layers are relevant.

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

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 Ishikawa et al. (2006); Zöttl and Stark (2014); Schaar et al. (2015); Lintuvuori et al. (2016). 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.

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 Enculescu and Stark (2011); Stark (2016). Hydrodynamic interactions between squirmers obviously do not destroy the polar order. In our parametric study of Fig. 3c), this is also confirmed for other squirmer types $\beta$. Differences occur in the layering and in the transitional region. For pullers the upright orientation in the bottom layer decreases for larger $\beta$, as expected by hydrodynamic interactions with the bottom wall in lubrication theory Ishikawa et al. (2006); Schaar et al. (2015); Lintuvuori et al. (2016). In the adjacent layers hardly any polar order is visible in contrast to neutral squirmers. Weak pushers ($\beta=-1$ and $-2$) show a similar but weaker trend compared to neutral squirmers. For strong pushers ($\beta=-5$), again, there is hardly any polar order in the layering.

III.2 Convection

As videos V1 and V2 in the supplemental material† 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 $\beta=0$ and $\alpha=1.5$ can be well fitted by a Gaussian with standard deviation comparable to $v_{0}$ (see Fig. 6). For large $|\beta|$ small deviations from the Gaussian form occur. Thus, also vertical squirmer speeds larger than $\frac{1}{3}v_{0}$ arise (the maximal vertical velocity a single bulk squirmer can have at $\alpha=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.

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

[TABLE]

where $\bm{x}=(x,y)$ is a position in the horizontal plane, $\langle\ldots\rangle_{\|}$ means average along the vertical, and $\rho$ is the squirmer density. In the following we always indicate a time average over some quantity $q$ in the steady state by $\overline{q}$. We plot $\overline{j_{z}(\bm{x})}$ in Fig. 7 for neutral squirmers and $\alpha=1.5$ in the $xy$ plane of the simulation box. While in the lower left 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. 3c) we present $\overline{j_{z}(\bm{x})}$ for different squirmer types at the same $\alpha$. For weak pullers ($\beta=1$ and $2$) the extent of the convection cell decreases and at $\beta=5$ large-scale convection is no longer observable. For weak pushers ($\beta=-1$ and $-2$) the current density becomes weaker but still large-scale convection is visible, which then vanishes for $\beta=-5$. For the case $\beta=0$, $\alpha=1.5$ we checked that large-scale convection is stable against a reduction in the squirmer density $\phi$ at constant $L_{z}$ and for reduced $L_{z}$ while keeping $\phi$ constant. In both settings a single convection cell extends across the simulation box.

To further characterize the squirmer density current $j_{z}(\bm{x})$, we calculate its zeroth and first moment. The zeroth moment is the volume average of the vertical current density:

[TABLE]

where $A$ is the area of the $xy$ plane of the simulation cell and $\langle\ldots\rangle$ means average over the whole volume of the exponential region. The vertical current density $\langle j_{z}\rangle$ strongly fluctuates in time (see Fig. 8) but in steady state its temporal mean, $\overline{\langle j_{z}\rangle}$, has to be zero because of particle conservation. Thus, we use $\langle|j_{z}|\rangle$ to quantify how mobile the squirmers are in the exponential region (see Fig. 8). In Fig. 9 we show the time average $\overline{\langle|j_{z}|\rangle}$ versus $\beta$ for different rescaled swimming velocities $\alpha$. As expected, at $\alpha>1$ the squirmers are more mobile than for $\alpha\leq 1$ since they are able to move against gravity. Furthermore, for $\alpha=1.5$ the mean vertical current density $\overline{\langle|j_{z}|\rangle}$ decreases for large $|\beta|$, revealing again the importance of the advective flow set up by the different squirmer types.

We call the first moment of $j_{z}(\bm{x})$ current dipole,

[TABLE]

where we choose the center of the $xy$ plane as the origin of $\bm{x}$.

The current dipole serves to quantify the strength and horizontal extension of the convection cell by its magnitude $j_{D}:=|\bm{j}_{D}|$. The cell’s orientation relative to the $x$-axis is given by the angle $\varphi_{D}$ with $\cos\varphi_{D}=\bm{j}_{D}\cdot\bm{e}_{x}/j_{D}$. Figure 10 shows how $j_{D}$ strongly fluctuates in time, reflecting again the high mobility of the squirmers. The orientation angle $\varphi_{D}$ also fluctuates and in the example of Fig. 10 assumes two mean orientations around $1.5\pi$ and $\pi$. Overall, we can record that the spatial arrangement of convection is subject to strong fluctuations and strongly variable in time.

Finally, in Fig. 11 we plot the time average $\overline{j_{D}}$ versus squirmer type $\beta$ for different $\alpha$. Convection is largest for large $\alpha$ and neutral squirmers. The current dipole we define in Eq. (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. 3c) we find that $\overline{j_{D}}$ for $\alpha=1.5$ decreases when the magnitude of the current density $j_{z}(\bm{x})$ decreases (weak pusher, $\beta<0$) or the extension of the the convection cell becomes smaller (weak puller, $\beta>0$). Strong pushers/pullers with $|\beta|=5$ only show weak convection. For small rescaled swimming speed $\alpha\leq 1$ convection is generally small.

IV 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; close 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 close packed layers, in particular for small $|\beta|$, 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 argued 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. Wolff et al. (2013) for very dilute swimmer suspensions. Preliminary results in denser systems show that instabilities occur due to hydrodynamic interactions Hennes et al. . In a harmonic trapping potential a similar instability leads to the formation of fluid pumps by breaking the rotational symmetry of the trap Nash et al. (2010); Hennes et al. (2014). In the present case we expect convectional patterns to occur. Thereby, we will connect to the phenomenon of bioconvection Wager (1911); Pedley and Kessler (1992); Williams and Bees (2011a, b).

Acknowledgements.

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.