Stability of arrays of bottom-heavy spherical squirmers

The incessant swimming motion of microbes in dense suspensions can give rise to striking collective motions and coherent structures. However, theoretical investigations of these structures typically utilize either computationally demanding numerical simulations or simplified continuum models. Here we analytically investigate the collective dynamics of a dense array of steady, spherical squirmers. We first calculate the forces and torques acting on two closely separated squirmers, through solving the Stokes equations to second order in the ratio of mean spacing to squirmer radius. This lubrication analysis is then used to assess the stability of a dense, vertical, planar array of identical three-dimensional squirmers. The system of vertically oriented squirmers is unstable if there is no short-range repulsive force between them, even when there is a strong gravitational torque on them because they are bottom-heavy. When there is a repulsive force the positions of the squirmers are stable, but the orientations are unstable unless the bottom-heaviness parameter $G_{bh}$ is sufficiently large. The predictions of instability and possible long time behavior are qualitatively the same for monolayers confined between two parallel rigid planes as for unconfined monolayers. The predictions compare favorably with published numerical simulations, and reveal the existence of additional dynamic structures not previously observed; puller-type squirmers show a greater range of structures than pushers. The use of pairwise lubrication interactions provides an efficient means of assessing stability of dense suspensions usually tackled using full numerical simulations.

Figure 1

Geometry of the two squirmers. The origin of the coordinate system is located on the surface of sphere 2 closest to sphere 1. The vector ${\mathit{e}}_{\rho }$ points radially in the $x-y$ plane, and the vector ${\mathit{e}}_{\varphi }$ is the azimuthal direction.

Figure 2

First- and second-order components of the pressure distribution in the lubrication region. (a) $Q_0\left(\rho \right)$ and (b) $F\left(\rho \right)$ are shown as functions of $\rho$. Results are shown for $\lambda =0.2$, 0.5, 1, 2, 5, and $\infty$ (dotted).

Figure 3

(a) Boundary element simulations [38, 39] for $G_{bh}=100$ and $\beta =1$ show the development of stable monolayers of bottom-heavy squirmers, with an equilibrium spacing ${\epsilon }_0=0.002$. (b) Diagram showing the domain of spherical squirmers. The numbering scheme for a $d\times d$ diamond cell is shown. (c) The domain is subject to periodic boundary conditions, as shown. The direction in which the force of gravity acts is denoted by the vector $\mathit{g}$.

Figure 4

Plots showing the linearized solutions associated with all squirmers in the $10\times 10$ configuration. The results for the perturbed squirmer and its six nearest neighbors are shown in red and blue respectively. The green curves correspond to all other squirmers. Results have been computed with ${\epsilon }_0=0.002$ and $\beta =1$.

Figure 5

Direction of the tangential velocity at the surface of the squirmer, associated with the first two modes. Results are shown for (a) $B_1>0$ and (b) $B_2>0$. The boundary condition is a superposition of these modes, given by $u_{\theta }{|}_{r=a}=B_{1}sin\theta +B_{2}sin\theta cos\theta$. (c),(d) Plots showing the linearized solutions associated with the perturbed squirmer, for several different values of $\beta$. Results have been computed for a $10\times 10$ diamond configuration with an equilibrium spacing of ${\epsilon }_0=0.002$. Results are shown for $\beta =0, \pm 1$, and $\pm 5$. Positive and negative values of $\beta$ are shown in green and red respectively and the blue curves correspond to $\beta =0$.

Figure 6

Linearized solutions with repulsive force present and $\beta =1$. Plots showing the linearized solutions associated with all squirmers in the $10\times 10$ configuration. The results for the perturbed squirmer and its six nearest neighbors are shown in red and blue respectively. The green curves correspond to all other squirmers. Results have been computed with ${\epsilon }_0=0.002, \zeta =0.01, \delta ={\epsilon }_0/100, {\kappa }_1=1, {\kappa }_2={10}^3, \beta =1, G_{bh}=0$.

Figure 7

One single squirmer is perturbed in an otherwise uniform monolayer. (a) Figure showing the trajectories of all 64 squirmers over the course of the simulation. For the purposes of plotting, the positions have been scaled so that the radius of each squirmer is 0, even though the simulation was conducted with $a=1$. (b) Plots showing the orientation $\zeta$ and (c) corresponding standard deviation of all squirmers as a function of time. Parameters used are ${\epsilon }_0=2\times {10}^{-3}, \beta =1, {\kappa }_1=0, \varphi =3\pi /2, \delta ={\epsilon }_0/1000, \zeta =1/100$, and $G_{bh}=20$. Simulations were repeated with repulsive force present (${\kappa }_1=1$) for [(d),(e)] $G_{bh}=20$ and [(f),(g)] $G_{bh}=50$. Under these conditions, the lattice is translationally stable, but the squirmers still require sufficiently large value of $G_{bh}$ for orientational stability.

Figure 8

Orientation of every squirmer in an $8\times 8$ diamond, as a function of time. In each simulation, every squirmer is initially given a random perturbation to its position and orientation, with amplitude $\delta ={\epsilon }_0/100$ and $\zeta =1/100$ respectively. Results have been computed for $G_{bh}=35$, 40, 45, and 50 over the interval $t\in [0,100]$. Also shown is the standard deviation of $\zeta$ as a function of time. Additional parameters used are given by ${\epsilon }_0=2\times {10}^{-3}, {\kappa }_1=1, {\kappa }_2={10}^3$, and $\beta =1$.

Figure 9

Lubrication simulations of a squirmer monolayer situated between two plane parallel walls. (a) The mean angle from vertical, $M$, and (b) average variance for each squirmer, $S$, are shown across a broad range of $G_{bh}$ and $\beta$ values. (c)–(h) Representative results from different regions of parameter space highlight the qualitatively different long term dynamics. Parameters used include ${\epsilon }_0={\epsilon }_0^{\text{wall}}=2\times {10}^{-3}, {\kappa }_1={\kappa }_1^{\text{wall}}=1, {\kappa }_2={\kappa }_2^{\text{wall}}={10}^3$. The parameter combinations in (c)–(h) are depicted as circles in (a).

Figure 10

(a) Depiction of coherent structures (Ishikawa et al. [39]) which have formed within a monolayer of spherical squirmers. Results have been computed using $\beta =5$ and $G_{bh}=100$. (b) Schematic diagram showing this equilibrium configuration in which all squirmers are oriented at an angle of $\zeta =\pm {\zeta }_0$ from vertical. In this configuration, three different equilibrium spacings ${\epsilon }_{0}a, k_1{\epsilon }_{0}a$, and $k_2{\epsilon }_{0}a$ are permitted.

Figure 11

Equilibrium spacing in ensemble of uniformly tilted squirmers. Plots showing the values of $k_1$ (smooth) and $k_2$ (dashed) as functions of ${\zeta }_0$ for (a) $G_{bh}=20$ and (b) $G_{bh}=100$. Results have been computed with ${\epsilon }_0=2/1000, {\kappa }_1=1$, and ${\kappa }_2={10}^3$ for various $\beta$.