Stability of arrays of bottom-heavy spherical squirmers

TL;DR

This paper analytically studies the stability of dense arrays of spherical squirmers, revealing conditions for stability and discovering new dynamic structures, with results aligning well with numerical simulations.

Contribution

It introduces an analytical lubrication-based method to assess the stability of dense squirmer arrays, reducing reliance on computationally intensive simulations.

Findings

Vertical arrays are stable if gravitational torque is strong enough.

Puller squirmers exhibit more diverse dynamic structures than pushers.

The analytical approach aligns with numerical simulation results.

Abstract

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 uniformly-spaced, vertically-oriented squirmers is stable if the gravitational torque experienced due to bottom-heaviness is sufficiently strong, and an intercellular repulsive force is included. The predictions of instability and possible long time behavior is qualitatively similar 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.

Equations

$$u_{r}\big{|}_{r=a}=\sum_{n}A_{n}(t)P_{n}(\cos\theta),\qquad u_{\theta}\big{|}_{r=a}=\sin\theta\sum_{n}B_{n}(t)W_{n}(\cos\theta),$$

$$W_n(\cos \theta )=\frac{2}{n(n+1)}{P}_n^{\prime }(\cos \theta ).$$

$$\begin{array}{rl}u& =u_0+{ϵ}^{1/2}{u}_1+\mathcal{O}(ϵ),\\ v& =v_0+{ϵ}^{1/2}{v}_1+\mathcal{O}(ϵ),\\ w& ={ϵ}^{1/2}{w}_0+ϵw_1+\mathcal{O}({ϵ}^{3/2}),\\ p& ={ϵ}^{-3/2}{p}_0+{ϵ}^{-1}{p}_1+\mathcal{O}({ϵ}^{-1/2}).\end{array}$$

$$H=1+\frac{\lambda +1}{2\lambda }{\rho }^2+\mathcal{O}(ϵ),$$

$$q_{0}(\rho)=Q_{0}(\rho)\sum_{n}B_{n}W_{n}\big{(}-\bm{e}\cdot\bm{e}_{z}\big{)}\quad\text{and}\quad Q_{0}(\rho)=\frac{6\mu}{5a}\frac{\rho}{H^{2}}.$$

$$p_1(\rho ,\varphi )=f_{\text{p}}(\rho )+g(\rho )\cos 2\varphi .$$

$$F(\rho)=\frac{3}{4}\bigg{(}\frac{\lambda}{\lambda+1}\bigg{)}\frac{6H-1}{H^{2}},$$

$$F(\rho)\sum_{n}\bigg{[}B_{n}W_{n}\big{(}-\bm{e}\cdot\bm{e}_{z}\big{)}\bm{e}\cdot\bm{e}_{z}+\frac{1}{2}B_{n}W_{n}^{\prime}\big{(}-\bm{e}\cdot\bm{e}_{z}\big{)}(\bm{e}\cdot\bm{e}_{x})^{2}\bigg{]}=\frac{a}{\mu}f_{\text{p}}(\rho).$$

$$\displaystyle\begin{split}F_{x}^{(1)}&=-\frac{4}{5}\mu\pi a\ \bm{e}\cdot\bm{e}_{x}\frac{\lambda(\lambda+4)}{(\lambda+1)^{2}}\sum_{n}B_{n}W_{n}\big{(}-\bm{e}\cdot\bm{e}_{z}\big{)}\big{(}\log\epsilon+\mathcal{O}(1)\big{)},\\ F_{z}^{(1)}&=-9\mu\pi a\frac{\lambda^{2}}{(\lambda+1)^{2}}\sum_{n}\bigg{[}B_{n}W_{n}\big{(}-\bm{e}\cdot\bm{e}_{z}\big{)}\bm{e}\cdot\bm{e}_{z}+\frac{1}{2}B_{n}W_{n}^{\prime}\big{(}-\bm{e}\cdot\bm{e}_{z}\big{)}(\bm{e}\cdot\bm{e}_{x})^{2}\bigg{]}\big{(}\log\epsilon+\mathcal{O}(1)\big{)},\\ T_{y}^{(1)}&=\frac{16\lambda}{5(\lambda+1)}\mu\pi a^{2}\ \bm{e}\cdot\bm{e}_{x}\sum_{n}B_{n}W_{n}\big{(}-\bm{e}\cdot\bm{e}_{z}\big{)}\big{(}\log\epsilon+\mathcal{O}(1)\big{)},\\ T_{y}^{(2)}&=\frac{4\lambda^{2}}{5(\lambda+1)}\mu\pi a^{2}\ \bm{e}\cdot\bm{e}_{x}\sum_{n}B_{n}W_{n}\big{(}-\bm{e}\cdot\bm{e}_{z}\big{)}\big{(}\log\epsilon+\mathcal{O}(1)\big{)}.\end{split}$$

$${\bar{\mathbf{F}}}_1^{\prime }$$

$${\bar{\mathbf{T}}}_1^{\prime }$$

$${\bar{\mathbf{F}}}_2^{\prime }$$

$${\bar{\mathbf{T}}}_2^{\prime }$$

$$\textbf{A}=\left(\begin{array}[]{c c c}\log\epsilon&0&0\\ 0&\log\epsilon&0\\ 0&0&-\frac{3}{2\epsilon}+\frac{27}{20}\log\epsilon\\ \end{array}\right),\qquad\textbf{C}=\left(\begin{array}[]{c c c}0&-\log\epsilon&0\\ \log\epsilon&0&0\\ 0&0&0\\ \end{array}\right),\qquad\textbf{J}=\left(\begin{array}[]{c c c}\log\epsilon&0&0\\ 0&\log\epsilon&0\\ 0&0&0\\ \end{array}\right).$$

$${\mathbf{T}}_{\text{grav}}^i=-\frac{4}{3}\pi a^3{\rho }_{f}h{\mathbf{e}}_i\times \mathbf{g},$$

$$G_{bh}=\frac{2\pi {\rho }_{f}gah}{\mu B_1}.$$

$${\mathbf{F}}_{\text{rep}}={\kappa }_1{\kappa }_2\frac{\exp (-{\kappa }_2ϵ)}{1-\exp (-{\kappa }_2ϵ)}\frac{\mathbf{r}}{r}.$$

$${\mathbf{e}}_z^{\prime }$$

$${\mathbf{e}}_y^{\prime }$$

$${\mathbf{e}}_x^{\prime }$$

$${\mathbf{e}}_i\cdot {\mathbf{e}}_z^{\prime }={\mathbf{e}}_i\cdot \hat{\mathbf{r}},\text{and}{\mathbf{e}}_i\cdot {\mathbf{e}}_x^{\prime }=\sqrt{1-({\mathbf{e}}_i\cdot \hat{\mathbf{r}}{)}^2}.$$

$$\displaystyle\left(\begin{array}[]{c}\bar{\bm{F}}_{1}\\ \vdots\\ \bar{\bm{F}}_{n}\\ \hline\cr\bar{\bm{T}}_{1}\\ \vdots\\ \bar{\bm{T}}_{n}\\ \end{array}\right)^{\text{net}}$$

$$\text{M}\cdot \mathbf{x}=\mathbf{R},$$

$$\sum _{i=1}^n{V}_{i,x}=\sum _{i=1}^n{V}_{i,y}=\sum _{i=1}^n{V}_{i,z}=0.$$

$$\mathbf{r}=a\delta (\sin \varphi {\mathbf{e}}_x+\cos \varphi {\mathbf{e}}_z),\delta \ll 1.$$

$$\mathbf{e}=\sin \zeta {\mathbf{e}}_x+\cos \zeta {\mathbf{e}}_z,\zeta \ll 1.$$

$$\mathbf{x}={\mathbf{x}}_0+\zeta {\mathbf{x}}_1^{\text{r}}+\delta {\mathbf{x}}_1^{\text{t}}+\dots$$

$$={\text{M}}_0+\zeta {\text{M}}_1^{\text{r}}+\delta {\text{M}}_1^{\text{t}}+\dots$$

$$\mathbf{R}$$

$${\text{M}}_0\cdot {\mathbf{x}}_1^r={\mathbf{R}}_1^r\text{and}{\text{M}}_0\cdot {\mathbf{x}}_1^t={\mathbf{R}}_1^t.$$

$$x_1^{r,i}=a_i{G}_{bh}+b_i.$$

$$x_1^{t,i}=c_i\sin (\varphi -d_i)+e_i.$$

$$V_x=V_x^r\zeta +V_x^t\delta ,V_z=V_z^r\zeta +V_z^t\delta ,\Omega ={\Omega }^r\zeta +{\Omega }^t\delta .$$

Full text

Introduction

Motility is a pervasive feature among microorganisms, from the diurnal migration of marine phytoplankton Bollens et al. (2011) to the motion of bacteria in the gut Berg (2008). The acquisition of resources Blackburn et al. (1998), evasion from predators Kiørboe et al. (2014), and infection by pathogens Josenhans and Suerbaum (2002) all depend sensitively on organismal motility. Early microscopes dating back to the 18th century van Leeuwenhoek (1700) enabled glimpses into the dynamic nature of the microbial world. Since then, the role of cilia and flagella – ubiquitous, highly conserved propulsive appendages – has received considerable attention Sleigh (1962); Brennen and Winet (1977). Recent advances in imaging and microfluidic control offer new insights into the mechanics of cellular propulsion Son et al. (2015). The spatial distribution of cells in microbial consortia can influence nutrient cycling Smriga et al. (2016), horizontal gene transfer Moor et al. (2017) and fertilization processes Denissenko et al. (2012). Developing a quantitative framework for the collective dynamics of swimming microorganisms is therefore essential to understanding a vast array of biological processes. It has become clear that collective motions of many microorganisms can be very different from individual dynamics Elgeti et al. (2015). Striking examples of bacterial turbulence Dunkel et al. (2013), self-organization Wioland et al. (2013); Thutupalli et al. (2018) and coherent structures Saintillan and Shelley (2012) reveal the combined effects of confinement, hydrodynamic signatures, and steric interactions in determining emergent phenomena.

The squirmer model was first proposed by Lighthill in 1952 Lighthill (1952) and modified by his student Blake in 1971 Blake (1971), but its current wide applicability to a range of organisms was not initiated until relatively recently Pedley (2016); Gilpin et al. (2017); Shapiro et al. (2014); Magar and Pedley (2005); Michelin and Lauga (2011); Lin et al. (2011). Its elegance and simplicity enable modelling of cells in different environments Matas-Navarro et al. (2014), near air-liquid interfaces Wang and Ardekani (2013), or no-slip and repulsive walls Li and Ardekani (2014); Llopis and Pagonabarraga (2010); Lintuvuori et al. (2016). The conceptually simple model replaces an array of flagella with a single, no-slip, deformable surface, thereby linking discrete ciliary beating with an effective surface slip velocity. The model organism, Volvox carteri Goldstein (2015), renowned for its exquisite spherical symmetry, exemplifies the squirmer model, with strong agreement between predictions based on measured flagellar dynamics Pedley et al. (2016); Brumley et al. (2014, 2015) and the observed motion of freely-swimming colonies Drescher et al. (2010). Experimental Drescher et al. (2009) and theoretical evidence Delfau et al. (2016) hints at the importance of near-field interactions in determining collective properties of suspensions of squirmers.

Ishikawa et al. Ishikawa et al. (2006) investigated hydrodynamic interactions between two spherical squirmers, utilizing both lubrication theory and multipole expansions to model closely- and widely- separated squirmers respectively. Boundary element simulations of dense suspensions revealed stable collective states and intriguing oscillatory modes Ishikawa and Pedley (2008); Ishikawa et al. (2008), in which squirmers self organize into a densely packed lattice. Despite the conceptual simplicity of the squirmer model, it remains unclear what mechanisms are responsible for these states, and the precise conditions under which they are stable. In this paper, we analytically solve the Stokes equations between two bottom-heavy, spherical squirmers, in the limit of close separation. We use these results to predict the collective dynamics of a dense array of squirmers, and show that both orientational and translational stability are mediated through gravitational torques exerted on the cells, and a cell-cell repulsive force.

I Interactions between spherical squirmers

I.0.1 Interactions due to squirming motion

In order to calculate the forces and torques arising from the short-range interactions between two spherical swimming microorganisms, we will utilize the squirmer model. The single-squirmer model will be developed in the reference frame in which the center of the spherical squirmer is at rest, and the fluid at infinity has velocity given by $-U\bm{e}$. The value $U$ is the swimming speed of the sphere and $\bm{e}$ is its orientation vector – the unit vector along the axis of symmetry. The boundary conditions at the surface of the sphere are given by

[TABLE]

where $\theta$ is the angle measured from the anterior of the squirmer, $P_{n}$ is the $n^{\text{th}}$ Legendre polynomial, and $W_{n}$ is defined as

[TABLE]

The ultimate goal will be to consider the hydrodynamic interaction between two adjacent squirmers whose positions and orientations are arbitrary. Without loss of generality, consider the problem of two closely-separated spherical squirmers, as depicted in Fig. 1. The frame is chosen such that the orientation of squirmer 1 lies in the $x$-$z$ plane. ie $\bm{e}_{1}\cdot\bm{e}_{y}=0$. Any configuration in a laboratory frame can be mapped to the situation shown in Fig. 1 through a suitable linear transformation. By linearity of the Stokes equations, the problem involving two squirming spheres in a fluid that is at rest infinitely far away can be broken down into two distinct problems. The first has the squirming-sphere boundary condition on sphere 1 and zero velocity boundary condition on sphere 2. The second problem has zero velocity on sphere 1 and the squirming-sphere boundary condition on sphere 2. Only the former problem will be studied, since solving this will immediately yield the solution to the latter.

The radii of spheres 1 and 2 are given by $a$ and $\lambda a$ respectively, and the minimum separation between the spheres is taken to be $\epsilon a$ (with $\epsilon\ll 1$). The origin of the coordinate system is located at the surface of sphere 2, on the axis joining the centers of the two spheres. The $z$-axis passes through the spheres’ centers, so that spheres 1 and 2 lie in the regions $z>0$ and $z\leq 0$ respectively. The surfaces of spheres 1 and 2 are determined by $z=h_{1}$ and $z=h_{2}$, respectively. Let the two spheres, 1 and 2, have orientation vectors $\bm{e}_{1}$ and $\bm{e}_{2}$ and squirming sets $\bm{B}^{(1)}=(B_{1}^{(1)}(t),B_{2}^{(1)}(t),\ldots)$ and $\bm{B}^{(2)}=(B_{1}^{(2)}(t),B_{2}^{(2)}(t),\ldots)$, respectively. Squirmers with zero radial velocity on the sphere surface will be considered ($A_{n}(t)=0\;\forall\;n$). Although Fig. 1 depicts a configuration with two spheres, the following lubrication analysis can also be applied to the interaction between a sphere and a plane wall by considering the case where $\lambda\rightarrow\infty$.

The fluid velocity $\bm{u}=(u,v,w)$ and pressure $p$ in the gap between the squirmers, are expanded in powers of $\epsilon$:

[TABLE]

Similarly, the separation between the squirmers, $H$, non-dimensionalized by $\epsilon a$, can be written as a function of $\rho$ for $\epsilon\ll 1$,

[TABLE]

where $\rho$ is the distance from the $z$-axis (see Fig. 1). By expanding and solving the Stokes equations in powers of $\epsilon$, the leading order pressure distribution $p_{0}(\rho,\phi)=q_{0}(\rho)\bm{e}\cdot\bm{e}_{\rho}$ can be found (see SI Section S1 for detailed calculation), where

[TABLE]

Similarly, the second-order pressure is found to be of the form

[TABLE]

The component proportional to $\cos 2\phi$ disappears upon integration with respect to $\phi$ and so does not provide a net contribution to the force exerted between the spheres. It therefore suffices to consider

[TABLE]

where $F(\rho)$ is defined by

[TABLE]

The functions $Q_{0}(\rho)$ and $F(\rho)$ represent first- and second-order pressure increases due to the squirming motion of sphere 1. These are shown in Figs. 2(a-b) respectively.

Using Eqs. (5) and (7), the fluid velocity in the gap between the squirmers, Eq. (3), can be solved to first-order (see SI Section S1). These expressions enable the forces and torques acting on the two spheres to be calculated explicitly.

[TABLE]

The tangential and normal forces acting on sphere two, $F_{x}^{(2)}=-F_{x}^{(1)}$ and $F_{z}^{(2)}=-F_{z}^{(1)}$ respectively, are equal and opposite to the values on sphere one. By symmetry, the torque $T_{x}$ is precisely equal to zero. The torque in the $z$-direction can be evaluated; however it is found that $T_{z}=\mathcal{O}(\epsilon)$ so this need not be pursued. The torque exerted on sphere 2 in the $y$-direction has an extra factor of $\lambda$ compared to the results for sphere 1, arising from the discrepancy between their radii. It is also worth noting that for $\lambda=1$ (equally-sized spheres), the torque exerted on sphere 2 is one quarter times that exerted on sphere 1. Normal gradients in the fluid velocity are greater at the surface of the squirmer than they are at the boundary of the no-slip sphere.

I.0.2 Interactions due to motion of spheres

In addition to the effects of squirming, the two spheres will also experience forces and torques due to their linear and angular velocities, intercellular steric interactions, and gravity. The analysis so far has been performed in the frame shown in Fig. 1. The unit vectors appearing in Eq. (9) will now be given the primed coordinates to indicate that they are viewed in this frame. Suppose that in the reference frame $S^{\prime}$ depicted in Fig. 1, the spheres possess linear and angular velocity vectors $\bm{V}_{i}^{\prime}=(V_{x,i}^{\prime},V_{y,i}^{\prime},V_{z,i}^{\prime})$ and $\bm{\omega}_{i}^{\prime}=(\omega_{x,i}^{\prime},\omega_{y,i}^{\prime},\omega_{z,i}^{\prime})$ respectively. The subscript $i=1,2$ denotes either sphere 1 or sphere 2. Let $\bm{F}_{i}^{\prime}$ and $\bm{T}_{i}^{\prime}$ be the force and torques acting on sphere $i=1,2$ in this frame $S^{\prime}$. The forces and torques are scaled according to $\bar{\bm{F}}^{\prime}=\bm{F}^{\prime}/(\mu\pi a)$ and $\bar{\bm{T}}^{\prime}=\bm{T}^{\prime}/(\mu\pi a^{2})$. The following relationships can then be established Kim and Karrila (2005):

[TABLE]

where, correct to order $\mathcal{O}(\log\epsilon)$, the matrices are given by

[TABLE]

As expected, $\bar{\bm{F}}_{2}^{\prime}=-\bar{\bm{F}}_{1}^{\prime}$ and the forces and torques arising due to linear velocities are zero when $\bm{V}_{1}^{\prime}-\bm{V}_{2}^{\prime}=0$. Note that these results correspond to the case involving two equally-sized spheres ($\lambda=1$). Rotation of the spheres in the $z$-direction does not produce forces or torques that are singular as $\epsilon\rightarrow 0$, and the torque in the $z$-direction remains finite as the spheres become arbitrarily close together. It follows that the entries in the third row of C and J are all zero to order $\mathcal{O}(\log\epsilon)$.

I.0.3 Additional interactions

If the squirmers are bottom-heavy, there is an additional external torque acting on each sphere due to gravity. For species such as Volvox, this mechanism facilitates swimming in an upwards direction (negative gravitaxis). If the distance between the center of gravity and center of the squirmer is given by $h$, in the direction opposite to its swimming direction in an undisturbed fluid, the gravitational torque on the $i^{\text{th}}$ squirmer is given by

[TABLE]

where $\rho_{f}$ is the density and $\bm{g}$ is the acceleration due to gravity. The parameter introduced by Ishikawa et al. Ishikawa et al. (2006) is adopted here, comparing the gravitational and viscous torques:

[TABLE]

The non-dimensionalized gravitational torque can then be rewritten as $\bar{T}^{i}_{\text{grav}}=T^{i}_{\text{grav}}/(\mu\pi a^{2})=-\frac{2}{3\pi}B_{1}G_{bh}\sin\zeta_{i}$, where $\zeta_{i}$ is the angle of the squirmer from vertical. A repulsive force between spheres was included in the numerical simulations of ref Ishikawa et al. (2006):

[TABLE]

This was done to avoid the prohibitively small time step required to prevent squirmers from overlapping. The parameter $\epsilon$ is again the separation between squirmers, non-dimensionalized by the squirmer radius $a$. The parameter $\kappa_{1}$ represents the strength of the repulsion while $\kappa_{2}$ dictates the range at which this repulsion becomes significant. The values adopted will be the same as those used previously Ishikawa et al. (2006), namely $\kappa_{1}=1$ and $\kappa_{2}=10^{3}$. Importantly, this repulsive force can be “switched off” simply by choosing $\kappa_{1}=0$.

I.0.4 Generalization to multiple spheres

A larger system containing $n$ spheres will now be examined, and the total force and torque acting on each sphere will be calculated. Each sphere will interact hydrodynamically with neighbors that are sufficiently close, experience a gravitational torque due to bottom-heaviness, and be subject to the short-range repulsive force. The analysis so far has been performed in the frame shown in Fig. 1. The unit vectors will now be given the primed coordinates to indicate that they are viewed in this frame. In accordance with the preceding lubrication theory, the squirmer is required to be oriented in the $x^{\prime}$-$z^{\prime}$ plane. That is, $\bm{e}_{i}\cdot\bm{e}_{y}^{\prime}=0$. For any ordered pair of spheres $i$ and $j$ in the laboratory frame $S$, it is possible to transform to the reference frame $S_{ij}^{\prime}$ in which spheres $i$ and $j$ are positioned as squirmers 1 and 2 respectively in Fig. 1. Suppose that in the laboratory reference frame $S$, two squirmers $i$ and $j$ have position vectors $\bm{r}_{i}$ and $\bm{r}_{j}$ respectively and that their orientations are given by $\bm{e}_{i}$ and $\bm{e}_{j}$ respectively. Let $\bm{r}=\bm{r}_{i}-\bm{r}_{j}$. The coordinate system $S_{ij}^{\prime}$ is defined in the following way:

[TABLE]

By construction, this frame satisfies the condition that $\bm{e}_{i}\cdot\bm{e}_{y}^{\prime}=0$. Thus, the lubrication analysis presented in Section I.0.1 can be directly applied in this frame. In the calculation of the forces and torques due to squirming, it is also necessary to know the quantities

[TABLE]

The lubrication analysis is used in frame $S_{ij}^{\prime}$ for the case when sphere $i$ is a squirmer and sphere $j$ has the zero boundary condition. The complementary problem involving sphere $j$ as a squirmer and sphere $i$ as a sphere with zero boundary condition is considered separately in frame $S_{ji}^{\prime}$ since the lubrication analysis is only applicable when $\bm{e}_{j}\cdot\bm{e}_{y}^{\prime}=0$. In this fashion, each pair of squirmers will be considered twice when calculating the total force and torque on the system.

Equations (10)-(13) outline the forces and torques due to translational and rotational velocities of the two spheres, where everything is measured in the frame $S_{ij}^{\prime}$. However, it is desirable to find these quantities in the laboratory frame $S$. This is achieved by utilizing the appropriate transformation matrix, $\textbf{R}_{ij}$. A complete matrix-vector equation can be assembled as follows

[TABLE]

Since the fluid is considered to be at zero Reynolds number, the net force and torque on every squirmer must be zero. This condition is imposed simply by setting every entry on the left hand side of Eq. (54) to zero. The resulting matrix-vector equation can then be solved to find the linear and angular velocities corresponding to this condition. In particular, a system of the following form must be solved:

[TABLE]

where $\bm{x}$ contains the linear and angular velocities of the spheres. There are several important features of this equation that will now be discussed. Firstly, note that the matrix M depends only on the positions of the squirmers. It is completely independent of the squirming parameters, the strength of gravity, the orientations of the squirmers and the repulsive force. Matrix M can be assembled once the physical configuration of the suspension is known. The vector $\bm{R}$ in Eq. (55) depends on all of the parameters involved in the problem. Secondly, note that the $6n\times 6n$ matrix M has a rank which is precisely $6n-3$. In order to understand this, recall that the only hydrodynamic forces and torques are those arising from the lubrication regions, which depend on the relative motion of the squirmers compared to each other. In order to ensure that the matrix M is nonsingular, an arbitrary reference frame must be chosen. The frame in which the bulk velocity of the configuration of squirmers is zero is chosen. That is,

[TABLE]

II Uniform Monolayer of squirmers in an unbounded fluid

The present formulation facilitates calculation of the linear and angular velocities of all squirmers in any configuration where lubrication forces dominate. The consequences of perturbing a uniform monolayer of squirmers, subject to periodic boundary conditions, will now be explored. Only cells whose corresponding squirming sets are independent of time ($B_{n}(t)=B_{n}\;\forall\;n$. See Eq. (1)) will be studied. Consider the diamond-shaped configuration shown in Fig. 3(b), in which the equilibrium spacing between any two adjacent squirmers is given by $\epsilon_{0}a$. In the equilibrium state, all squirmers have an orientation vector $\bm{e}=\bm{e}_{z}$. For the time being, the motion of the squirmers is limited to the plane of the monolayer. Moreover, all translational and orientational perturbations are restricted to this plane, giving rise to what is essentially a three-dimensional system (2 translational $+$ 1 rotational).

II.1 Analytical approach

The consequences of perturbing the position and orientation of one squirmer will now be investigated. At this stage, time-evolution of the system will not be studied. The purpose of this section is to analytically address the behavior of the system in the small time limit. This corresponds to constructing and solving the matrix-vector equation in Eq. (55) once only. Without loss of generality, the perturbed cell is chosen to be squirmer 1, as depicted in Fig. 3(b). Let the origin of the coordinate system coincide with the center of this squirmer in its equilibrium position. A translational perturbation is initiated, of magnitude $a\delta$ in the direction $\phi$, such that the position vector of the squirmer is given by

[TABLE]

The orientation of the squirmer is also perturbed by $\zeta$, so that

[TABLE]

The matrix M and vector $\bm{R}$ depend on the small parameters $\delta$ and $\zeta$. A solution of the following form is sought

[TABLE]

where the superscripts ‘r’ and ‘t’ represent rotation and translation respectively. Other components of the matrix system can be linearized in the same way:

[TABLE]

These expressions are substituted into the original matrix-vector equation Eq. (55), and various orders of $\zeta$ and $\delta$ are equated. The vector $\bm{R}_{0}$ corresponds to the equilibrium configuration and is equal to $\bm{0}$. It follows that the leading-order solution is $\bm{x}_{0}=\bm{0}$. As one might expect, a suspension of evenly spaced squirmers, all pointing in the $z$-direction, do not experience a net force or torque. With this in mind, it is found that

[TABLE]

Recall that the matrix M depends only on the positions of the individuals cells. Thus, in the equilibrium configuration, $\textbf{M}_{0}$ depends only on the scaled equilibrium spacing, $\epsilon_{0}$. For a given suspension, the value of $\epsilon_{0}$ will be known. As such, $\textbf{M}_{0}$ can be constructed and inverted without knowing anything about the squirming parameters. The vectors $\bm{R}_{1}^{r}$ and $\bm{R}_{1}^{t}$ can be subsequently constructed. By definition, $\bm{R}_{1}^{r}$ must be independent of $\delta$ and thus $\phi$. From the form of $\bm{R}_{1}^{r}$, it follows that the elements of the solution $\bm{x}_{1}^{r}$ must be of the form

[TABLE]

The matrix $\textbf{M}_{0}$ cannot be inverted without assuming a particular value of $\epsilon_{0}$. As such, the coefficients in the above equation must be numerically fitted. The form of the solution arising through small perturbations in the position of the squirmer will now be examined. By definition, $\bm{R}_{1}^{t}$ is independent of $\zeta$ and thus $G_{bh}$. Terms in $\bm{R}_{1}^{t}$ involve either $\sin(\phi-\phi_{0})$ or $\cos(\phi-\phi_{0})$ for some $\phi_{0}$. Since $\textbf{M}_{0}$ is independent of $\phi$, a solution of the following form is sought:

[TABLE]

II.1.1 No repulsive force

Consider the situation in which the repulsive force between adjacent squirmers is absent. This is achieved by setting the value of $\kappa_{1}$ in Eq. (17) to be equal to zero. For a given $\epsilon_{0}$ and set of squirming parameters, the leading-order solution in the form of Eq. (59) is readily found:

[TABLE]

The term $\Omega^{r}$ represents the restoring effect that gravity has on the orientation of the squirmer, and is found to be directly proportional to $G_{bh}$ for the perturbed squirmer. When $\delta=0$ (ie. no translational perturbation), the orientation will be restored if $G_{bh}>0$ (since $\Omega^{r}<0$). For larger $\delta$, a correspondingly larger value of $G_{bh}$ is required to ensure that small perturbations to the orientation decay. In fact, the critical value of $G_{bh}$ is given by

[TABLE]

for some $k>0$. The ratio of the first two squirming modes, defined in Eq. (1), is given by

[TABLE]

The parameter $B_{2}$ is proportional to the stresslet of the squirmer, so $\beta<0$ and $\beta>0$ represent pushers and pullers respectively. Perturbing one squirmer in the configuration, will in general, affect all squirmers in the monolayer. Figure 4 summarizes the results for $\beta=1$. The central squirmer (red) is given either a rotational or a translational perturbation, as shown in Figure 4(a). The subsequent linear and angular velocities of all squirmers are then shown (Figs. 4(b-d) and Figs. 4(e-g) respectively). In each case, the blue curves represent the solutions for squirmers adjacent to the perturbed cell, and green represents the remaining cells in the monolayer. In the case of a rotational perturbation, the central squirmer will experience $V_{x}^{r}>0$ (Fig. 4(b)), indicating translational instability. The orientational perturbation will decay for that squirmer, but destabilizes the surrounding cells (Fig. 4(d)). Translational perturbations in the $x$ and $z$ directions will decay and grow respectively (see Figs. 4(e-f)), but at the same time, will destabilize the orientation of the central squirmer (Fig. 4(g))). Taken together, these results demonstrate that any perturbations to the uniformly spaced planar array will be linearly unstable, with rotational perturbations causing translational instability, and vice versa. Rotating the cell clockwise or anticlockwise will cause it to move right or left, respectively. Similarly, translating the squirmer right or left will cause it to move clockwise or anticlockwise, respectively.

For the case $\beta=-1$, the results (not shown) are extremely similar to those presented in Fig. 4. However, the sign of the red curves in panels (e) and (f) is reversed. In order to understand these results, it is helpful reconsider the mechanisms through which the squirming occurs. Figure 5(a-b) shows the direction of the tangential velocity for the first two modes of squirming. For $B_{2}>0$, the second squirming mode serves to draw fluid from the poles of the squirmer ($\theta=0,\pi$) to the equator ($\theta=\pi/2$). When the position of the squirmer is perturbed in the $x$-direction, this mode restores the position of the squirmer. Conversely, for perturbations in the $z$-direction, this active drawing of fluid away from the poles results in further destabilization from the equilibrium position. The results are reversed for $B_{2}<0$ (and therefore $\beta<0$).

The functions $V_{x}^{r}$, $V_{z}^{r}$, $\Omega^{r}$ and $\Omega^{t}$ associated with the perturbed squirmer do not vary with $\beta$. Any changes in the value of $\beta$ are manifested only in $V_{x}^{t}$ and $V_{z}^{t}$, the linear velocities associated with translational perturbations. Consider the plots in Fig. 5(c-d), which show these quantities for several different values of $\beta$. The linear velocity of the squirmer after a translational perturbation is directly proportional to $\beta$ (see red curves). It is emphasized again, that the angular velocity of the perturbed squirmer is independent of $\beta$. These observations however, are not in general true for the rest of the squirmers in the configuration.

Recall that there exists a critical value of $G_{bh}$, above which perturbations to the orientation of the squirmer will decay, regardless of the direction, $\phi$, of the translational perturbation. This critical value was shown to depend only on $\Omega^{r}$ and $\Omega^{t}$ (see Eq. (66)). It has just been found that these two functions associated with the perturbed squirmer are independent of the value of $\beta$ used. It thus follows that $G_{bh}^{\text{critical}}$ does not depend on the ratio $\beta=B_{2}/B_{1}$ of the squirming velocities. Importantly, the results obtained in this section are applicable only in the small time limit. The linear and angular velocities have been analyzed for particular squirmer configurations, but the time dependence of the problem has not yet been considered.

II.1.2 Repulsive force present

In the previous section it was found that perturbations in the position of the squirmer could be either unstable or stable, depending upon the direction of the perturbation $\phi$, and the squirming parameters. The previous analysis will now be repeated, but with the repulsive force outlined in Eq. (17) included. The velocities in the $x$ and $z$-directions can again be found, as well as the angular velocities for all of the squirmers. Figure 6 shows the linearized solutions for a $10\times 10$ configuration of squirmers. The reconstructed solutions $V_{x}$, $V_{z}$ and $\Omega$ for particular values of $\zeta$ and $\delta$ have been plotted, to demonstrate the significance of the repulsive force. It is evident that the position of the central squirmer will be stable subject to small perturbations in either the $x$ or $z$-directions. However, the functions $\Omega^{r}$ and $\Omega^{t}$ have not changed upon inclusion of the repulsive force. Consequently, the critical value of $G_{bh}$ required to eliminate small perturbations in $\zeta$ remains the same. That is, $G_{bh}^{\text{critical}}$ is independent of both the value of $\beta$ and the presence of the repulsive force.

II.2 Numerical approach

In Section II.1, the analytical form of the linear and angular velocities associated with squirmers in a large uniform suspension were studied. In particular, the case where one squirmer was subject to small translational and angular perturbations was considered. However, these results were not able to address the long-term behavior of the suspension following a perturbation from the equilibrium. A numerical study into the dynamics of the monolayer will now be undertaken, using the formulation presented in Section I.0.4. The configuration is the same as the one presented in Fig. 3(b) and is again subject to periodic boundary conditions, as depicted in Fig. 3(c).

To begin with, the consequences of perturbing one squirmer in the uniform monolayer will be explored. For small perturbations, this system was studied analytically in Section II.1. It was found that the position of the squirmer was stable, provided the repulsive force between adjacent squirmers was included. It was also found that there exists a critical value of $G_{bh}$, above which perturbations to the orientation of the squirmer will decay, regardless of the direction, $\phi$, of the translational perturbation. Recall that this value was independent of the presence of the repulsive force. The functional form of this critical value is given in Eq. (66), and was derived under the assumption that the neighboring squirmers were all left unperturbed. The results in Fig. 7 show the consequences of perturbing one squirmer in an $8\times 8$ system. The equilibrium spacing is again considered to be $\epsilon_{0}=2\times 10^{-3}$ and the perturbation is given by $\zeta=1/100$ and $\delta=\epsilon_{0}/1000$ with $\phi=3\pi/2$. From the analysis in Section II.1, it is known that for these parameters, to ensure $\Omega<0$ for the perturbed squirmer requires $G_{bh}>1.84$. The value $G_{bh}=20$ is used, which is known to be well beyond this critical value.

From the previous analytical work, it is known that the value of $G_{bh}$ used here guarantees $\Omega<0$ for the first time-step. However, the only way the long-time behavior can be assessed is through these numerical simulations. For $t\ll 1$, the system without intercellular repulsive force ($\kappa_{1}=0$) appears to be stable, with the orientation of the perturbed squirmer beginning to be restored (Fig. 7(b)). However, as time progresses, the other squirmers in the system begin to move (Fig. 7(a)). Indeed, the system becomes unstable as time progresses, with the orientation of all squirmers growing in magnitude. The perturbations that develop in the surrounding squirmers will act to destabilize the central squirmer.

It was found earlier that in the absence of the repulsive force, for $\beta>0$, each squirmer is stable and unstable to translational perturbations in the $x$ and $z$-directions respectively, with the converse true for $\beta<0$. Figure 7 demonstrates this phenomenon clearly, with the squirmers drifting towards one another in the $z$-direction. The repulsive force outlined in Section II.1 is now reinstated. This prevents the squirmers from coming too close together, since translational perturbations are quickly eliminated. Consider the plots in Fig. 7(d-g), which show the orientation of 64 steady squirmers over the interval $t\in[0,50]$, for two different values of $G_{bh}$. Even with the orientation of only one single squirmer perturbed, the whole system eventually becomes unstable for $G_{bh}=20$. However, the system is stable for large $t$ when $G_{bh}=50$.

Since the repulsive force quickly restores the position of the squirmers to their equilibrium value, where $\delta\rightarrow 0$, the value of $G_{bh}$ used is expected to be well above the critical value derived earlier. Nevertheless, instability is observed among the orientation of the squirmers. In the early stages of the simulation, the orientations of the neighbors become perturbed, causing the original perturbed squirmer to become further destabilized. This results in an increase in the corresponding value of $G_{bh}$ required to eliminate all angular perturbations. Since the restoring force quickly returns the squirmers to their equilibrium positions, the critical value of $G_{bh}$ required for angular stability does not depend strongly on the translational perturbations initially given to the squirmers. The stability of the orientation of the squirmers is dictated by the angular velocities associated with small perturbations to the orientations rather than positions. The interactions between $\Omega^{r}$ for various squirmers are key in determining the critical value of $G_{bh}$ required for stability. In addition to the fact that perturbations to the orientations grow when $G_{bh}<G_{bh}^{\text{critical}}$, another interesting feature of Fig. 7(d) is the splitting of these orientations in a dichotomous fashion. As time progresses, the squirmers rotate away from vertical in a coordinated manner. It will be shown later that this phenomenon also occurs in other configurations in which instabilities develop.

We investigated the effects of perturbing the position and orientation of all squirmers in the uniform monolayer. The squirmers were given a random perturbation to both their orientation and position, with amplitudes $\zeta=1/100$ and $\delta=\epsilon_{0}/100$ respectively. The stabilizing effect that gravity has on the suspension is evident (see Fig. 8). For $G_{bh}=35$ and 40 the system is unstable, while for $G_{bh}=50$ the system is stable for large $t$. Although not shown here, the results for $G_{bh}=0$ yield $\text{std}(\zeta)\rightarrow\pi/\sqrt{3}$ for large $t$, corresponding to a uniform distribution in which there is no preferred orientation.

III Monolayer of squirmers between vertical rigid walls

In the preceding sections, the positions of the squirmers were restricted to lie in the $x$-$z$ plane. This condition can be relaxed in order to permit out of plane motion. Under these conditions, the monolayer is unstable, with perturbations in the $y$-direction growing (see Supplementary Information Section S2 for detailed analysis). However, a uniform monolayer of squirmers in an otherwise unbounded fluid is, in any case, an unrealistic situation. One method of maintaining a monolayer of spherical squirmers is to use a Hele-Shaw cell which is sufficiently thin. To this point, the forces and torques on the squirmers arising due to either sphere-sphere interactions or the effects of gravity have been considered. It is straightforward to extend to the case where the uniform monolayer of steady, spherical squirmers is situated between two plane parallel walls. The two planes are defined by $y=\pm a(1+\epsilon_{0}^{\text{wall}})$ so that the minimum clearance between the squirmers and the wall in the equilibrium configuration is $\epsilon_{0}^{\text{wall}}a$. Gravity is still considered to act in the negative $z$-direction.

As the lubrication analysis presented earlier applies to two spheres, each of arbitrary radius, the forces and torques acting on a squirmer interacting with the planes can easily be found by taking the limit $\lambda\rightarrow\infty$. The forces acting on the squirmers due to their translational and rotational motion must also be considered. For this, the results presented in ref Kim and Karrila (2005) are applied. A short-range repulsive force between the spheres and the walls is also incorporated into the model, as in Eq. (17), with the parameters $\kappa_{1}^{\text{wall}}$ and $\kappa_{2}^{\text{wall}}$. Equation (54) is modified by the inclusion of extra terms to account for the walls (see SI Eq. (S72)).

To explore the influence that squirming strength and bottom heaviness have on the monolayer stability, we performed 702 simulations across a range of values for $G_{bh}$ and $\beta=B_{2}/B_{1}$ (see Fig. 5(a-b) for schematic). This enables us to explicitly investigate the differences between pushers ($\beta<0$) and pullers ($\beta>0$). In each simulation, the squirmers’ orientations were subject to random initial conditions, and the long time dynamics were observed. The effect of the repulsive force is to stabilize the positions of the squirmers, retaining the lattice like structure. It is therefore sufficient to consider the orientation from vertical, $\theta_{i}$, of squirmers in each monolayer. Several qualitatively different dynamics emerge, depending on the parameter combination $(\beta,G_{bh})$.

In order to quantify the dynamics for various parameter choices $(\beta,G_{bh})$, we analyze the time-dependent angle $\{\theta_{i}(t)\}$ for all squirmers $i=1,\ldots,N$ in a given simulation. Firstly, we define $M=\langle|\theta_{i}(t)|\rangle_{t,i}$, averaged over time and all squirmers in the monolayer. The parameter $M=0$ if and only if all squirmers converge to a vertical orientation at large $t$. However, in the case of $M\neq 0$, this parameter is unable to distinguish between steady states (e.g. Fig. 9(e)) and chaotic results (Fig. 9(h)). We therefore define a second parameter, $S$, calculated by taking the variance of each time-dependent signal $\theta_{i}(t)$, and subsequently averaging over the squirmer population.

[TABLE]

The parameter $S$ will be zero if every squirmer converges to a constant orientation, regardless of its value. This parameter therefore provides great utility in distinguishing between equilibrium structures and other results. For each of the 702 simulations, the parameters $M$ and $S$ were calculated, the results of which are summarized in Fig. 9(a) and Fig. 9(b) respectively.

Across the range of parameters studied, 5 different characteristic behaviors were observed for $t\gg 1$. The simplest case is that in which all squirmers eventually orient vertically, $\theta_{i}\rightarrow 0$ (Case II in Fig. 9, see for example Fig. 9(d)), and is identified as when both $M=0$ and $S=0$. This is precisely the equilibrium structure initially observed by Ishikawa et al. Ishikawa and Pedley (2008); Ishikawa et al. (2008) for $\beta=1$, $G_{bh}=100$, and investigated analytically in Section II.1. From the random initial conditions studied here, the system can converge to this vertical state for either pushers ($\beta<0$) or pullers ($\beta>0$), provided $G_{bh}$ is sufficiently large. For pushers ($\beta<0$), only one other type of behavior is possible, in which all squirmers converge to a non-zero equilibrium orientation (Case I in Fig. 9, see for example Fig. 9(c)). These dynamics occur when the second-order squirming mode is large enough to destabilize the vertically oriented monolayer.

For pullers ($\beta>0$), a richer set of dynamics is possible. For a given value of $G_{bh}$, increasing $\beta$ beyond a critical value results in an abrupt transition from $M=0$ to $M>0$. The system adopts a bistable state, in which squirmers possess a finite and constant tilt angle (Case III in Fig. 9, see Fig. 9(e)), qualitatively similar to Case I. Increasing $\beta$ further results in oscillations about these values (Case IV, see also Fig. 9(f,g)). For sufficiently large $\beta$, the entire system becomes unstable (Case V). Figure 9(h) illustrates these unstable dynamics, with the orientation of one squirmer shown in red.

IV Investigation of tilted structures and oscillatory states

The numerical simulations of Section III revealed the existence of stable states in which all squirmers adopt a non-zero mean orientation from vertical, either constant in value or oscillating in time. By symmetry, the configuration in which all squirmers in the lattice are vertically oriented is a steady state, and the linear and angular velocities of all squirmers in the periodic lattice will be zero. However, the conditions under which the “tilted equilibrium” can occur are not immediately clear. Ishikawa et al. also discovered stable coherent structures in which the squirmers do not orient themselves in a vertical direction, even in the presence of strong bottom-heaviness (see Fig. 10). The only difference between Figs. 3(a) and 10(a) is that the value of $\beta$ has been increased from 1 to 5. This corresponds to a shift in parameters equivalent to moving from Case II to Case IV in Fig. 9. The structure in which all squirmers possess some orientation of magnitude $\zeta_{0}$ from vertical, as shown in Fig. 10(b), will now be studied. It will be the goal of this section to understand the nature of this equilibrium state.

Until now, the equilibrium spacing between adjacent cells has been considered to be uniform throughout. Three different values are now permitted, namely $\epsilon_{0}a$, $k_{1}\epsilon_{0}a$ and $k_{2}\epsilon_{0}a$. In the case where $k_{1}=k_{2}=1$, it is found that the net force on the squirmers is zero provided $\zeta_{0}=0$ or $B_{1}=0$. The former case has already been studied in detail, whilst the latter case involves squirmers that would not be able to swim in an unbounded fluid (swimming speed = $2B_{1}/3$). If the columns in Fig. 10(b) are evenly spaced, with $k_{1}=k_{2}\neq 1$, then the same conditions are required. The equilibrium configuration depicted in Fig. 10(a) cannot be achieved with $k_{1}=k_{2}$ unless $B_{1}=0$. Such a configuration would, by symmetry, be independent of the inter-particle repulsive force.

Consider now the case where $k_{1}\neq k_{2}$. In this scenario, it is immediately obvious that the net force experienced by each squirmer as a consequence of the repulsive force presented in Eq. (17) will be non-zero, and so, the existence of an equilibrium configuration will depend on the presence of this force. Nevertheless, the analysis is continued in an attempt to account for the results in Fig. 9 and boundary element simulations of Ishikawa et al. (Fig. 10(a)). By specifying the values of $\epsilon_{0}$, $\kappa_{1}$ and $\kappa_{2}$, it is possible to calculate the values of $k_{1}$ and $k_{2}$ for any $\beta$, $G_{bh}$ and equilibrium orientation $\zeta_{0}$. For a given experiment, $\beta$ and $G_{bh}$ will be known a priori and so the spacings $k_{1}\epsilon_{0}a$ and $k_{2}\epsilon_{0}a$ will be functions of the equilibrium orientation $\zeta_{0}$. Figures 11(a) and 11(b) show the values of $k_{i}$ associated with $G_{bh}=20$ and $G_{bh}=100$ respectively. The curve in Fig. 11(b) corresponding to $\beta=5$ incorporates exactly the same parameters as in Fig. 10(a). In order to determine the equilibrium orientation $\zeta_{0}$, an additional piece of information is required. It would be possible, for instance, to demand that the mean equilibrium spacing between adjacent cells is equal to $\epsilon_{0}$. That is, $(k_{1}+k_{2})/2=1$. The corresponding values are given by $k_{1}=1.07134$, $k_{2}=0.92866$ and $\zeta_{0}=1.12638$. This value of $\zeta_{0}$ is very similar to that observed in Fig. 10(a). The advantage of choosing values of $k_{i}$ as close to 1 as possible is that it minimizes the effects associated with the repulsive force. For given values of $\beta$ and $G_{bh}$ it is possible to find the equilibrium spacings $k_{1}$ and $k_{2}$, and orientation $\zeta_{0}$. Although this reveals the existence of a tilted equilibrium configuration, it does not assess the stability of the monolayer in that case. Depending on the parameter configuration $(\beta,G_{bh})$, the monolayer may be bistable (Case III), oscillatory (Case IV), or completely unstable (Case V).

V Discussion and Conclusions

While fully resolved boundary element simulations of interacting spherical squirmers revealed stable hexagonal lattice configurations Ishikawa and Pedley (2008), the model complexity prevented a simple understanding of the mechanisms behind this stability. Moreover, the computational intensity precluded a broad exploration of parameter space relevant to a range of biological and synthetic micro-swimmers. Here we have developed a semi-analytical framework to predict the dynamics of dense suspensions of spherical squirmers. We began by solving the Stokes equations to second order between closely-separated squirmers. This followed similar steps to ref Ishikawa et al. (2006), but to higher order, as required to calculate the normal force. These analytical expressions were then utilized in a ‘lubrication simulation’, to assess the global dynamics of a dense monolayer of squirmers. This revealed that pairwise lubrication interactions, in conjunction with a short-range repulsive force, were sufficient to account for the stable states observed in previous studies. This framework therefore provides a computationally inexpensive means of investigating the dynamics of dense suspensions of swimming microorganisms.

Initial studies of the monolayer restricted the motion of the spheres to lie in a plane, even though the fluid was unbounded and three-dimensional. Further analysis of this monolayer confirms the intuitive result that it is unstable subject to small out-of-plane perturbations (see SI Fig. S1 for further information). The inclusion of nearby plane parallel walls, as in the case of a rigid Hele-Shaw cell, maintains the structure of the monolayer (see SI Fig. S2), with orientational perturbations again eliminated for sufficiently large values of $G_{bh}$. For every value of $\beta$ and $G_{bh}$ studied, suspensions of pushers ($\beta<0$) were stable for large $t$, with all squirmers converging either to vertical, or a finite tilt angle. Conversely, pullers ($\beta>0$) exhibited a range of qualitatively different states (see Fig. 9), with orientations being completely unstable for sufficiently large $\beta$.

In the present framework, we have neglected any density difference between the squirmers and the fluid, which would lead to sedimentation Drescher et al. (2009). The inclusion of a Stokeslet term would modify the flow through the monolayer, and therefore potentially influence the stability calculations. This is the subject of future work.

The equilibrium spacing between adjacent squirmers, $\epsilon=0.002$, was chosen to match the stable value emerging from full boundary element simulations Ishikawa and Pedley (2008). Under these conditions, the logarithmic singularities (see Eq. (9)) dominate the expressions for the hydrodynamic forces and torques. This paper focusses on the collective dynamics of monolayers of spherical squirmers, but the framework could be readily extended to model fully 3D concentrated suspensions. Although the present analysis could in principle also be applied to polydisperse suspensions, it is likely that substantial variations in the separation $\epsilon$ would limit applicability of the lubrication approximations. The colonial alga Volvox carteri is a very good realization of Lighthill’s spherical squirmer Lighthill (1952) (with $\beta<0$), but there are significant experimental challenges in preparing a monodisperse suspension of Volvox. Experimental investigation of the present system is therefore most likely to be achieved for large suspensions of identical synthetic microswimmers Thutupalli et al. (2018) situated in a vertical Hele-Shaw cell.

Acknowledgements

The authors thank T. Ishikawa, R.E. Goldstein, and M. Polin for useful discussions, and The University of Melbourne’s High Performance Computer Spartan. This work was supported by a Gates Cambridge Scholarship, a Human Frontier Science Program Cross-Disciplinary Fellowship, and a Discovery Early Career Researcher Award DE180100911 (D.R.B.).

S1 Hydrodynamic interactions between spherical squirmers

This analysis follows similar steps to that presented in ref Ishikawa et al. (2006), but with some corrections, and completed to higher order. The spatial coordinates are scaled according to $\epsilon^{1/2}aX=x$, $\epsilon^{1/2}aY=y$, $\epsilon aZ=z$ and thus $\epsilon^{1/2}a\rho=\rho^{*}$. The boundaries of spheres 1 and 2 within the lubrication region can then be written as follows:

[TABLE]

where $\rho=\sqrt{X^{2}+Y^{2}}$. By linearity of the Stokes equations, the problem involving two squirming spheres in a fluid that is at rest infinitely far away can be broken down into two distinct problems. The first has the squirming-sphere boundary condition on sphere 1 and zero velocity boundary condition on sphere 2. The second problem has zero velocity on sphere 1 and the squirming-sphere boundary condition on sphere 2. Only the former problem will be studied, since solving this will immediately yield the solution to the latter. For a solitary squirmer immersed in a fluid which is at rest at infinity, one can express the fluid velocity field as

[TABLE]

where $\bm{e}$ is the swimming direction, $\bm{r}$ is the position vector and $r=|\bm{r}|$. By writing the position vector of sphere 1 as $\bm{r}_{1}=\bm{r}-(1+\epsilon)a\bm{e}_{Z}$ where $\bm{e}_{Z}=(0,0,1)$ and performing a Taylor series expansion of W_{n}\big{(}\tfrac{\bm{e}\cdot\bm{r}_{1}}{r_{1}}\big{)} about the point $-\bm{e}_{1}\cdot\bm{e}_{z}$, the fluid boundary condition on the surface of sphere 1 can be written in the form

[TABLE]

where the functions $\bm{u}_{A}$, $\bm{v}_{A}$ and $\bm{w}_{A}$ are expressed as infinite series over the squirming modes. From the boundary conditions presented in Eq. (S4), it seems logical to attempt to express the velocity and pressure as power series in $\epsilon^{1/2}$. Consider the following expansions for the fluid velocity $\bm{u}=(u,v,w)$ and pressure $p$:

[TABLE]

With these in mind, the $x$, $y$ and $z$-components of the Stokes equations to various orders in $\epsilon$ can be extracted.

S1.0.1 The first-order solution

From the leading-order $z$-component of the Stokes equations, it is evident that $p_{0}=p_{0}(X,Y)$. Furthermore, on spheres 1 and 2, the fluid velocity is given by $\bm{u}_{A}$ and $\bm{0}$ respectively. Subject to these boundary conditions, the leading-order $x$ and $y$-components of the Stokes equations are integrated to find the following expressions:

[TABLE]

where $H=H_{1}-H_{2}=1+\frac{\lambda+1}{2\lambda}\rho^{2}+\mathcal{O}(\epsilon)$. Integrating the leading-order component of the continuity equation and utilizing Eqs. (S8) and (S9) yields the Reynolds equation:

[TABLE]

where we have made use of the fact that \bm{\rho}\cdot\bm{u}_{A}=-\sum_{n}B_{n}W_{n}\big{(}-\bm{e}\cdot\bm{e}_{z}\big{)}\bm{e}\cdot\bm{\rho}. Recall that $p_{0}=p_{0}(X,Y)$, or equivalently $p_{0}(\rho,\phi)$. A solution to Eq. (S10) of the form $p_{0}(\rho,\phi)=q_{0}(\rho)\bm{e}\cdot\bm{e}_{\rho}$ is found, where

[TABLE]

Figure 2(a) shows how $Q_{0}(\rho)$ varies as a function of $\rho$ for various values of $\lambda$.

S1.0.2 The second-order solution

In a fashion completely analogous to the first-order case, the second-order fluid velocities can be derived. These are equivalent to the expressions in Eqs. (S8) and (S9) but with $\bm{v}_{A}$ replacing $\bm{u}_{A}$.

[TABLE]

In addition, the corresponding Reynolds equation for the next order pressure contribution, $p_{1}$, is as follows:

[TABLE]

The complete solution to Eq. (S14) is of the form

[TABLE]

where the particular integral is

[TABLE]

and $g(\rho)$ is the solution to a second-order differential equation (not shown). It is not necessary to solve for $g(\rho)$ since it will not contribute to the overall force on the spheres anyway. A new function, $F(\rho)$, is defined so that

[TABLE]

where

[TABLE]

Figure 2(b) shows the dependence of $F(\rho)$ on $\rho$ for various values of $\lambda$. Importantly, these results correspond to the contribution to the pressure which is independent of $\phi$. The contribution which is proportional to $\cos 2\phi$ will disappear upon integration with respect to $\phi$.

Before commencing the evaluation of the forces acting on sphere 1, it will be necessary to calculate the fluid velocity in the gap. Equations (S8) and (S9) contain the cartesian components of the fluid velocity in the $x$ and $y$-directions, correct to first-order, with analogous solutions for the second-order case (Eqs. (S12)-(S13)). The fluid velocity in the $\rho$ and $\phi$-directions at each order in $\epsilon$ is required, since the subsequent analysis will involve finding the rate-of-strain tensor in cylindrical coordinates.

S1.0.3 First-order velocities

From Eqs. (S8), (S9) and (S11), the leading-order fluid velocity in the $\rho$ and $\phi$-directions can be written respectively as

[TABLE]

The components $u_{0,\rho}$, $u_{0,\phi}$ and $u_{0,z}$ are defined so that $u_{0,\rho}(\rho,z)\bm{e}\cdot\bm{e}_{\rho}=U_{0}$, $u_{0,\phi}(\rho,z)\bm{e}\cdot\bm{e}_{\phi}=V_{0}$ and $u_{0,z}(\rho,z)\bm{e}\cdot\bm{e}_{\rho}=\epsilon^{1/2}W_{0}$. The new function $v_{M}$ is defined so that

[TABLE]

The fluid velocities can thus be expressed in the following way:

[TABLE]

S1.0.4 Second-order velocities

The second-order pressure distribution is given by $p_{1}(\rho,\phi)=f_{\text{p}}(\rho)+g(\rho)\cos 2\phi$, where $f_{\text{p}}(\rho)$ and $g(\rho)$ are known functions. The second-order fluid velocity in the $\rho$ and $\phi$-directions can be written as

[TABLE]

S1.0.5 Tangential force

The component of the force acting tangential to the surface of sphere 1 in the gap region is given by dF_{x}=\bm{e}_{x}\cdot\big{[}\bm{\sigma}\cdot\bm{n}\big{]}dA, where $\bm{\sigma}=-p\textbf{I}+2\mu\bm{\varepsilon}$ is the stress tensor and the particle surface is defined as $A$. It can be seen that

[TABLE]

and so

[TABLE]

In the region between the spheres, $\pi-\theta\ll 1$. Since $\sin\theta=\rho^{*}/a=\epsilon^{1/2}\rho$, it is appropriate to change variables according to $\rho=\epsilon^{-1/2}\sin\theta$, and so $d\rho=\epsilon^{-1/2}\sqrt{1-\epsilon\rho^{2}}d\theta$. In addition, $\cos\theta=-1+\rho^{2}\epsilon/2+\mathcal{O}(\epsilon^{2})$. It follows that

[TABLE]

The required components of the rate-of-strain tensor, $\bm{\varepsilon}$, can be evaluated using the expressions for the fluid velocities. The force element in the $x$-direction is found to be

[TABLE]

and so

[TABLE]

Since the form of $v_{0,\rho}$ and $v_{0,\phi}$ are known, the leading-order force can be calculated:

[TABLE]

where $\rho_{0}$ denotes the extent of the lubrication region. The three integrals can be evaluated analytically. The inner solution must now be matched with the outer solution in order to find a suitable value for the boundary of the lubrication region, $\rho_{0}$. Using $\rho_{0}=D/\epsilon^{1/2}$ where $D\sim\mathcal{O}(1)$, it follows that $\log\rho_{0}^{2}\sim-\log\epsilon$. With this in mind, the tangential force can be determined:

[TABLE]

This expression in Eq. (S31) represents a small correction to the results presented in ref Ishikawa et al. (2006). However, the results for identically sized squirmers ($\lambda=1$) – the situation of most practical interest – are unchanged. The process of calculating the tangential force on sphere 2 is almost identical to the former case, with only a few small modifications to the analysis being necessary. Since $p_{0}=p_{0}(\rho,\phi)$ is independent of $Z$, the pressure at the surface of sphere 2 is the same as at the surface of sphere 1. Analysis reveals that the tangential force exerted on sphere 2 is $F_{x}^{(2)}=-F_{x}^{(1)}$.

S1.0.6 Normal force

Consider the normal force component acting on sphere 1, given by

[TABLE]

Although the above equation has a dependence on $\theta$, the formulation will continue in cylindrical coordinates. The force element can be written as

[TABLE]

where

[TABLE]

The change of variables defined by $\rho=\epsilon^{-1/2}\sin\theta$ is again utilized. Since $\cos\theta=-1+\rho^{2}\epsilon/2+\mathcal{O}(\epsilon^{2})$, the dependence of Eq. (S33) on $\theta$ can be removed, and it is found that the contributions to the force at various orders in $\epsilon$ are given by

[TABLE]

Since $p_{0}$ is proportional to $\cos\phi$, the integral in Eq. (S36) is identically zero. Consider now the next-order contribution to the normal force, given in Eq. (S37). The function $p_{1}$ has already been found and is given in Eq. (S15). The term in Eq. (S15) which is proportional to $\cos 2\phi$ will be zero upon integration with respect to $\phi$. Thus, Eq. (S37) can be expressed solely in terms of $f_{\text{p}}(\rho)$:

[TABLE]

where $f_{\text{p}}(\rho)$ is defined in Eq. (S16). The corresponding integral can be evaluated analytically and subsequently expanded asymptotically for $\rho_{0}\gg 1$. Performing the same matching procedure as in the preceding analysis for the tangential forces, where $\rho_{0}=D/\epsilon^{1/2}$, it is found that the total force exerted on sphere 1 in the $z$-direction is

[TABLE]

The force exerted on sphere 2 in the $z$-direction can also be found, and is given by $F_{z}^{(2)}=-F_{z}^{(1)}$. The normal force contribution was overlooked in ref Ishikawa et al. (2006), but will form an important component of our subsequent analysis.

S1.0.7 Torque

The torque element exerted on sphere 1 in the $y$-direction is given by

[TABLE]

Noting that $\varepsilon_{\rho^{*}z}=\varepsilon_{z\rho^{*}}$ and remembering that $\rho=\epsilon^{-1/2}\sin\theta$ on sphere 1, the leading-order contribution to $T_{y}$ can be found to be

[TABLE]

which can be evaluated and simplified to yield

[TABLE]

By symmetry, the torque $T_{x}$ is precisely equal to zero. The torque in the $z$-direction can be evaluated, however it is found that $T_{z}=\mathcal{O}(\epsilon)$ so this need not be pursued. The torque exerted on sphere 2 in the $y$-direction is readily computable using the results of the preceding results. Most of the working is the same as before, and it is found that

[TABLE]

There exists an extra factor of $\lambda$ compared to the results for sphere 1, arising from the discrepancy between their radii. It is also worth noting that for $\lambda=1$, the torque exerted on sphere 2 is one quarter times that exerted on sphere 1. Normal gradients in the fluid velocity are greater at the surface of the squirmer than they are at the boundary of the no-slip sphere, giving rise to this somewhat counterintuitive result.

S2 Monolayer of squirmers in an unbounded fluid: unrestricted motion

To this point, the position and orientation vectors of the squirmers have been restricted to lie in the $x$-$z$ plane. Importantly, it was found that for certain combinations of parameters, there exists a critical value of $G_{bh}$ above which small perturbations to the orientations and positions of the squirmers will decay, and the equilibrium monolayer configuration is restored. In particular, this stability depended on the presence of the repulsive force. Without it, the monolayer was unstable, even for large $G_{bh}$. Perturbations to the position and orientation vectors will now be permitted to be out of the plane of the equilibrium configuration depicted in Fig. 3. For this the full matrix-vector equation presented in Eq. (55) is constructed and solved for the linear and angular velocities. The coordinates and orientation of each squirmer can be represented by $\bm{r}_{i}=(x_{i},y_{i},z_{i})$ and $\bm{e}_{i}=(\sin\theta\cos\phi,\sin\theta\sin\phi,\cos\theta)$ respectively. The time-stepping is done in the same fashion as before, using the Runge-Kutta-Fehlberg method (RKF45).

Firstly, consider the effects of perturbing the positions of the squirmers in any of the three spatial dimensions. For these simulations, the initial orientations of the squirmers were left unperturbed, pointing in the $z$-direction. The positions were given a perturbation in a random direction, of magnitude $\epsilon_{0}a/100$. The results corresponding to $G_{bh}=50$ and $\beta=1$ are presented, since this parameter combination previously gave rise to a stable monolayer.

It is clear from Fig. S1(a) that the small perturbations to the squirmers’ positions quickly grow, with the monolayer beginning to drift apart in the $y$-direction. At larger times, the squirmers move with constant velocity in the $y$-direction. Although the initial perturbation is only in the squirmers’ positions, Fig. S1(b) illustrates that their orientations quickly deviate from their equilibrium value. Results have also been computed for the case where only the orientations of the squirmers in the equilibrium monolayer are perturbed. This gives rise to a set of results which are qualitatively the same as those presented in Fig. S1. When motion of the squirmers was limited to lie in the plane of the monolayer, the configuration was stable (see green curves in Fig. 8(b)). It is clear, however, that the extra degrees of freedom in the squirmers’ position and orientation serve to destabilize this monolayer. This result is in accordance with previous findings Ishikawa and Pedley (2008) where elaborate coherent structures formed in 2D were not observed in full 3D simulations.

S3 Monolayer of squirmers between rigid walls

It has already been shown that the uniform monolayer of squirmers in an unbounded fluid is unstable when subjected to small translational or rotational perturbations in the direction perpendicular to their plane (SI Section S2). This out-of-plane motion is now suppressed by including the two walls. To achieve this, Equation (54) is modified by the inclusion of additional forces and torques given by:

[TABLE]

The matrices $\textbf{M}_{i}^{\text{wall}}$ are block-diagonal since the forces and torques acting on any sphere are independent of the linear and angular velocities of any other sphere. By again demanding that the squirmers are all force- and torque-free, it is possible to construct a matrix-vector equation of the same form as Eq. (55). In the absence of the two walls, a reference frame had to be chosen with some arbitrary velocity (see Eq. (56)). Upon inclusion of the walls, the choice is no longer arbitrary, since the forces and torques arising through interaction with the planes depend on their velocities relative to the squirmers.

Consider the consequences of initiating [small] perturbations to the positions and orientations of the squirmers in the monolayer between the walls. The results corresponding to $\beta=1$, $\epsilon_{0}=2\times 10^{-3}$ and $\epsilon_{0}^{\text{wall}}=5\times 10^{-3}$ are presented, and the differences between the cases involving $G_{bh}=0$ and $G_{bh}=50$ examined. From Figs. S2(a) and S2(c) it is clear that the walls have the effect of stabilizing the positions of the squirmers. In this respect, they help to maintain the structure of the monolayer. This phenomenon is insensitive to the value of $G_{bh}$ used. Figures S2(b) and S2(d) demonstrate the effect that changing $G_{bh}$ has on the system. For $G_{bh}=0$, the orientation of the squirmers grows with time while for $G_{bh}=50$ the squirmers in the monolayer are ultimately restored to their equilibrium orientation. This is the same behavior that was exhibited in Fig. 8 where a critical value of $G_{bh}$ was required to maintain the equilibrium configuration of the monolayer. The difference, now, is that the nearby plane walls hold the monolayer together in the $y$-direction instead of having to ignore motion in that direction.