Pairwise scattering and bound states of spherical microorganisms

Suspensions of active matter feature prominently in biological and engineered systems, and have generated great interest among fluid dynamicists and soft matter physicists [1]. Notable examples include suspensions of swimming microorganisms [2], collective motion of active filaments [3] and a range of chemically- [4], magnetically- [5, 6] and acoustically-driven [7] particles. Underpinning striking experimental observations such as bacterial turbulence [8], stable bound states [9], and clustering of particles [10], is an interplay between hydrodynamic and steric interactions, and responses to external forces and torques [11, 12]. Theoretical works of active suspensions have typically centred around continuum models – either in the dilute limit [13] or through minimal models [14] – and numerical simulations which resolve interactions across multiple length scales [15, 16]. Recent work has also shown that summation of pairwise lubrication interactions between closely-separated spherical squirmers is sufficient to capture the bulk rheological properties of concentrated suspensions, across a range of packing fractions and external field strengths [17, 18].

In the pursuit of an understanding of active matter collective dynamics, experimental works have systematically examined the fluid flows generated by propulsive appendages [19, 20, 21], swimming microorganisms [22, 23, 24, 25], the interactions between pairs of cells [9] and active droplets [26], and the role of nearby boundaries [27, 28, 29]. Colonies of Volvox exhibit striking bound states [30] mediated by gravity and a bounding wall, while cells of Chlamydomonas reflect in non-trivial ways from no-slip surfaces [31]. Detailed numerical studies have carefully calculated the scattering angles between pairs of squirmers in quiescent Stokes flow [15, 32, 33], and examined the role of inertia [34], density stratification [35], and unsteady swimming [36]. For axisymmetric interactions, it is possible to find closed form analytical solutions for the interaction between squirmers [37]. However, for general interactions, works involving closely-separated squirmers typically involve computationally-intensive simulations to calculate the fluid flows and resolve the organisms’ trajectories [15]. The success of lubrication theory in predicting selected hydrodynamic bound states [9] as well as bulk suspension rheology [17] offers a tantalising prospect that it may be used more generally to calculate pairwise dynamics of hydrodynamically coupled microorganisms. This paper utilises analytical expressions for forces and torques between closely-separated spherical squirmers whose orientations are coplanar, calculating the hydrodynamic interactions between them and resolving their full two-dimensional trajectories. This enables rapid calculation of non-trivial scattering dynamics between pairs of arbitrarily positioned squirmers, across a range of squirming parameters (pushers vs. pullers) and for varying degrees of bottom-heaviness.

The principal goal of this paper is to calculate the scattering dynamics of pairs of colliding spherical microorganisms of equal radius $a$. The celebrated squirmer model [38] is used to capture the motility of the cells, which move due to radial and tangential velocity boundary conditions,

respectively, 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

The swimming speed of a solitary squirmer is related to the primary squirming mode, $V_{s}=2B_{1}/3$, and $\beta=B_{2}/B_{1}$ captures the stresslet strength. We consider squirmers with zero radial velocity ($A_{n}(t)=0\;\forall\;n$), and although the analysis can be completed to arbitrary order in $B_{n}$, we will truncate the tangential squirming modes to $n=1,2$. The overall velocity boundary condition is therefore given by

The cartesian coordinates of the spheres’ centres $\bm{x}_{i}$ and orientations $\bm{p}_{i}$ (for $i=1,2$) are considered to lie in the $x$-$y$ plane. The spheres may be bottom heavy, so that when the swimming direction, $\bm{p}_{i}$, is not vertical, the squirmer experiences a gravitational torque,

where $v$ and $\rho$ are the cell volume and average density, respectively; $h$ is the displacement of the centre of mass from the geometric centre; and $\bm{g}$ is the gravitational vector with $|\bm{g}|=g$. At each instant in time, the squirmers each experience forces and torques due to their intrinsic swimming motion (propulsion); hydrodynamic drag due to the background fluid; mutual hydrodynamic interactions; and the action of gravity. We also include a repulsive interparticle force to prevent the spheres from overlapping,

where $\hat{\bm{r}}$ is the unit vector along the line of centres and $\epsilon$ is the minimum dimensionless clearance between the two cells, $\epsilon a=|\bm{x}_{1}-\bm{x}_{2}|-2a$. We apply $\kappa_{1}=10$ and $\kappa_{2}=100$.

To calculate the hydrodynamic interactions between two squirmers, we utilise linearity of the Stokes equations and solve for the interaction between a squirmer and a no-slip sphere. Figure 1(a) shows the arrangement of the squirmer-sphere pair used to derive the forces and torques due to squirming in the lubrication limit. The solution in this frame can be mapped to any relative position/orientation of squirmers through a linear transformation. Without loss of generality, we apply the squirming-sphere boundary condition on sphere 1 and zero velocity boundary condition on sphere 2. The forces and torques acting on the two spheres were previously calculated explicitly [18], and are given by

The above expressions are valid for all tangential squirming modes, but henceforth we will consider only $n=1,2$. We also consider spheres of equal radii, $\lambda=1$.

Consistent with previous works [18, 39], the total dimensional force $\bm{F}_{i}$ and torque $\bm{T}_{i}$ on each squirmer are scaled according to $\bar{\bm{F}}_{i}=\bm{F}_{i}/(\eta_{0}\pi a)$ and $\bar{\bm{T}}_{i}=\bm{T}_{i}/(\eta_{0}\pi a^{2})$, so that they each have units of velocity. We introduce a dimensionless number $G_{bh}$, defined as

which measures the relative importance of gravitational torques compared to viscous torques for the squirmers. The gravitational torque on each squirmer can then be written as

Each squirmer is subject to a propulsive force of constant magnitude, parallel to its orientation vector, according to

which ensures propulsion of an isolated squirmer, at the prescribed swimming speed $V_{s}$. The complete resistance formulation for the two spheres in Stokes flow is then given by

where the resistance matrix is $\bm{R}=\bm{R}^{sq-sq}+\bm{R}^{drag}$. The matrix $\bm{R}^{sq-sq}$ incorporates the hydrodynamic forces and torques acting on both spheres, arising from their linear and angular velocities. These expressions correspond to squeezing, shearing, and rotation of no-slip spheres in the lubrication limit, and can be found in full in [18] and [40]. The diagonal matrix $\bm{R}^{drag}$ represents the hydrodynamic drag due to the translation and rotation of an isolated sphere in Stokes flow, and is given by

where $\bm{I}_{n}$ is the $n\times n$ identity matrix. On the right hand side of Eq. (10), the forces, $\bar{\bm{F}}^{sq}$, and torques, $\bar{\bm{T}}^{sq}$, arise due to the squirming motion of each sphere. We refer the reader to [18] for detailed expressions noting that a change of reference frame must be made for each pair, in order for the expressions to apply. The resistance formulation in Eq. (10) can be written in the form

where $\mathbf{x}$ is a $6\times 1$ vector containing the position and orientation of both spheres (confined to lie in the $x$-$y$ plane), and both $\bm{R}$ and $\bm{S}$ depend on $\mathbf{x}$. This enables us to calculate the linear and angular velocities of both squirmers based on their respective positions and orientations. However, we note that interparticle interactions – mediated through $\bm{R}^{sq-sq}$, $\bar{\bm{F}}^{sq}$, $\bar{\bm{T}}^{sq}$ and $\bar{\bm{F}}^{rep}$ – occur only in the lubrication limit. The squirming forces and torques exhibit a logarithmic dependence [18] on the interparticle spacing, $\epsilon=(|\bm{x}_{1}-\bm{x}_{2}|-2a)/a$, so $\epsilon=1$ is a natural choice for the limit for lubrication interactions.

For widely-separated squirmers ($\epsilon>1$), the resistance matrix is diagonal, and therefore Eq. (10) can be solved explicitly to yield

As expected, solitary squirmers will swim at a speed $V_{s}$ in a direction $\bm{p}_{i}$ which can change due to the presence of gravitational torques.

In the absence of gravitational torques, widely-separated squirmers swim with a constant velocity, $V_{s}\bm{p}_{i}$. Deviations from these trajectories only occur if the spacing is sufficiently close ($\epsilon\leq 1$) so as to facilitate hydrodynamic interactions. In order to characterise the pairwise scattering dynamics of non-bottom heavy squirmers, it is therefore sufficient to consider pairs whose initial non-dimensionalized separation is exactly $\epsilon=1$. The cases where squirmers begin with a larger separation will either have the squirmers fail to come into close contact with one another and continue in straight lines, or reach a point where $\epsilon=1$ so that the remaining path is described by our results. We consider the initial positions of the spheres to lie on the $x$-axis, as shown in Fig. 1(a), with orientations $\phi_{i}$ measured anticlockwise from the positive $x$-axis. Equation 12 is integrated directly over $t\in[0,200]$, to find the position and orientation of both squirmers as functions of time. The duration of these simulations is long enough to ensure that squirmers had either departed from each other’s vicinity ($\epsilon>1$) and thereafter drifted away from one another in straight lines, or entered a state that they would stay in indefinitely, such as a steady state, closed orbiting pattern or periodic motion. At the end of each simulation, the final angles $\Phi_{i}$ were recorded for each squirmer $i$. We also recorded the duration of the interaction, $\tau$, defined as the length of time for which the two squirmers interact hydrodynamically ($\epsilon\leq 1$). Unless otherwise stated, all positions and distances will be considered as being non-dimensionalized with respect to the sphere radius $a$.

Although it might appear natural to express the outgoing squirmer angles $\Phi_{i}$ in terms of incoming angles $\phi_{i}$, we note that half of the possible initial conditions result in squirmers immediately swimming apart from one another (e.g., for $\text{mod}(\phi_{2},2\pi)<\text{mod}(\phi_{1},2\pi)$). Furthermore, these configurations do not correspond to initial collisions between straight-swimming squirmers, since reversing the direction of time would result in hydrodynamic interactions. When characterising the scattering properties of the squirmers, it is therefore most appropriate to consider the outgoing angles $\Phi_{i}$ in terms of the difference, $\Delta=(\phi_{2}-\phi_{1})/\pi$, and sum, $\Sigma=(\phi_{2}+\phi_{1})/\pi$, of the incoming angles. Figure 1(c) classifies the different dynamics observed in the numerical simulations.

Many microorganisms, including Volvox [9], Chlamydomonas [41] and Heterosigma akashiwo [42], exhibit a displacement between their geometric centre and centre of mass, resulting in an effective bottom-heaviness. We now consider the effect of an external gravitational torque, as outlined in Eq. (4), on the pairwise dynamics of spherical squirmers. To exhibit representative dynamics, we focus on the case where $\alpha=0$, so that the direction of gravity is perpendicular to the line joining the spheres’ centres (see Fig. 1), though in principle this could be varied. An isolated squirmer whose initial orientation $\bm{p}$ is not parallel to the gravitational field $\bm{g}$, will experience an external torque which balances the viscous torques of rotation, to gradually reorient the squirmer until $\bm{p}$ is antiparallel to $\bm{g}$. However, for pairs of squirmers, the hydrodynamic interactions provide additional forces and torques which can facilitate richer dynamics.

For sufficiently strong $G_{bh}$, the gravitational torque is expected to dominate the dynamics, resulting in upward swimming motion of the squirmer pair. Figures 8(a-b) show representative results for $G_{bh}=10$, with $\beta=1$ and $\beta=0$ respectively. The hydrodynamic interactions temporarily compete with gravitational torques, but ultimately the squirmers’ swimming direction is righted against gravity. Figure 8(c) reveals an interesting bound state in which both stronger gravity ($G_{bh}=50$) and greater stresslet ($\beta=-5$) result in fluctuations in the spacing between the squirmers. With the squirmers splayed away from one another, they interact solely through their posterior regions. The pusher boundary condition on each squirmer provides a torque which opposes that due to bottom-heaviness. Based on their orientations, the cells have a tendency to swim away from one another. However, this is balanced by the fact that pushers ($\beta<0$) draw fluid away from their equator, and therefore attract one another in this region.

For smaller values of $G_{bh}$, an interplay between the hydrodynamic interactions and gravitational torques can give rise to non-trivial bound states. For strong pushers, pairs of squirmers may exhibit orbiting motion, with the pair moving diagonally against gravity. Figure 8(d) shows an example of one such trajectory, which we refer to as a “drifting orbit”. The orbiting direction depends on the initial angles of the two squirmers. These results can be reconciled with Fig. 4(b,d), in which strong non-bottom heavy pushers may orbit one another. If gravity is sufficiently weak, it may bias the mean direction of the orbiting pair without disrupting the perpetual bound state.

To investigate the parameter combinations giving rise to drifting orbits, we calculated the squirmer trajectories over a range of $G_{bh}$, $\beta$, $\Delta$ and $\Sigma$ values. For given $G_{bh}$ and $\beta$, we calculated the proportion of initial orientations that resulted in drifting orbits, the results of which are summarised in Fig. 8(g). In Section III.2, we saw that for $\beta=-5$, $G_{bh}=0$, approximately 23% of incoming angles result in permanent trapping (red shaded region of Fig. 4(d)), effectively all of which were orbiting dynamics. By increasing $G_{bh}$ from zero, we see that the overall probability of orbiting decreases, up to a value of $G_{bh}\approx 2.5$, where orbiting no longer occurs, and the squirmers separate.

Highly specific symmetries can give rise to hydrodynamic bound states of downward swimming squirmers. The squirmers in Fig. 8(e), oriented symmetrically about the $y$-axis with $\Delta=1.7$ and $\Sigma=1$, are subject to gravitational torques which would tend to rotate the squirmers towards one another. However, this is balanced by the viscous torques generated by pullers interacting through anterior contact regions. Importantly, the results are qualitatively different to those in Fig. 8(c), since pullers draw fluid from their poles to the equator, and therefore ‘repel’ one another when aligned close to parallel. This perfectly balances the propulsive force, giving rise to an equilibrium state. As is the case for an isolated squirmer, this downward swimming state is unstable, with a small asymmetry in the initial swimming angles ultimately resulting in a U-turn (see Fig. 8(f) for $\Delta=1.75$, $\Sigma=1.05$). However, we note that the bottom-heavy squirmer pair is capable of swimming downwards for a considerable amount of time before gravitational torques eventually reorient the pair. Specifically, the bound pair in Fig. 8(f) swim 4.5 times further downward than they would as isolated squirmers, revealing that hydrodynamic trapping can provide ‘inertia’ from the influence of external torques.

The main focus of this paper has been to calculate the dynamics of colliding spherical squirmers, over a range of incoming angles, squirming parameters and external gravitational fields. Recent work has demonstrated that the bulk rheology of a concentrated suspension of spherical squirmers can be quantitatively captured by considering only lubrication interactions [17] – that is, hydrodynamic interactions only between cells and their nearest neighbours. Here we have explicitly calculated the swimming trajectories of pairs of squirmers subject to these short-range hydrodynamic interactions, revealing good agreement with simulations in which hydrodynamic interactions are fully resolved using boundary element methods [15]. Our results reveal a rich range of dynamics, including transient scattering, stalling motion, pairwise linear and circular orbits.

We have examined the pairwise dynamics for squirmers with initial angles $\phi_{1}$ and $\phi_{2}$, measured anticlockwise from the positive $x$-axis. However, it is more appropriate to consider the difference, $\Delta=(\phi_{2}-\phi_{1})/\pi$, and sum, $\Sigma=(\phi_{2}+\phi_{1})/\pi$, of the incoming angles over the range $[0,2]$, since this excludes initial configurations in which squirmers swim directly away from one another. The simulations were run for sufficiently long so as to determine whether cells would interact transiently, and swim apart with scattered angles, $\Phi_{1}$ and $\Phi_{2}$, or remain bound indefinitely. In all cases, we calculated the squirmers’ change in orientation $\delta_{i}=(\Phi_{i}-\phi_{i})/\pi$ due to the interaction. Neutral squirmers exhibit deflections of magnitude up to $\approx 2\pi/3$ ($\delta=\pm 0.64$, see Fig. 2). The results exhibit strong discontinuities in scattering angles across the line which determines whether squirmers slide past one another or rotate to swim as a pair. Importantly, all interactions are transient, except for highly specific symmetric cases ($\Delta=0$, $2$ and $\Delta=\Sigma=1$).

Introducing a stresslet through a non-zero value of $\beta$ results in qualitatively different results. For strong pullers ($\beta=5$), squirmers may exhibit pairwise bound states, in which cells swim together with a constant velocity. For both $\beta=1$ and $\beta=5$, Fig. 3 reveals strong discontinuities in the outgoing angle emerging with respect to incoming angles. It is also worth noting that as $\beta$ increases, the sign of the deflection switches for many of the initial configurations. Pairs of pushers exhibited the richest behaviour, with a large fraction of the initial conditions resulting in hydrodynamic bound states (see Fig. 4). For strong pushers ($\beta=-5$), we observed pairwise linear motion, or circular orbiting motion in a direction that depends on the squirmers’ initial orientations (clockwise for $\Sigma>1$, anticlockwise for $\Sigma<1$).

Including the effect of bottom-heaviness results in external torques which must balance those from the squirming motion and the relative motion of the two spheres. Many initial configurations ultimately converge towards a pair of upward swimming cells, but other observations include periodic orbiting with drift, downward swimming, or oscillating bound states.

The good agreement between lubrication theory results (LT) and boundary element method (BEM) highlights the utility of the present method in efficiently calculating the scattering dynamics. Indeed, Ishikawa et al. showed that BEM and LT predict similar forces and torques for squirmer separations as large as $\epsilon\sim 1$ [15]. The ability for lubrication theory to capture hydrodynamic bound states and orbiting dynamics of Volvox carteri – an experimental realisation of a spherical squirmer – near no-slip boundaries has been well-documented [9, 30], but the scattering dynamics have not been systematically explored. Near-field hydrodynamics are also known to dominate the collective motions of ellipsoidal squirmers [43]. Other works involving individual squirmers swimming near a wall [44, 45] have demonstrated hydrodynamic oscillations. For that case, in which the ratio of minimum clearance to squirmer radius was $\sim 0.2$, lubrication theory interactions were not expected to be very accurate [44]. However, in the present model of colliding squirmers, the minimum separation is significantly smaller (e.g., $\epsilon\sim 0.028$ for $\beta=0$). For pairs of squirming spheres subject to axisymmetric interactions, Papavassiliou & Alexander were able to use the reciprocal theorem to find exact solutions for the dynamics [37]. However, this was only a partial solution, since the non-axisymmetric components have not been determined. Moreover, the authors did not consider the effects of external torques, for example due to gravity. The present lubrication theory scheme provides greater generality across squirmer configurations and external torques, albeit with some compromise of accuracy, and the boundary element method represents an accurate benchmark for assessing solutions.

There has been considerable work examining the interactions of other kinds of microswimmers, e.g., phoretic particles [46], active droplets [26], active particles [10] and mixed suspensions [47]. However, in addition to hydrodynamic interactions, chemical trails can be far-reaching, and therefore affect the far-field dynamics. By examining only the near field hydrodynamic interactions, we have overlooked the structure of the flow field generated by individual squirmers. This could be important for strong pullers, which generate closed streamlines in their wake. A nearby squirmer could experience flow in the opposite direction to the tangential boundary condition. This likely contributes to the discrepancy between LT and BEM methods for $\beta=5$ in the specific case of orthogonal approach, Fig. 7. However, since organisms generally approach one another and collide with anterior hemispheres, this is likely to be important only in limited cases.

Although the focus of this paper has been on the interaction between identical squirmers whose orientations are confined to the same plane, this work could be generalised to non-identical squirmers colliding with arbitrary orientation vectors. Another possible extension to this work involves considering axial rotation of squirmers while swimming, with particular application to Volvox [39]. Finally, the classification of pairwise collisions presented in this paper could be used to predict the collective motion of non-dilute suspensions of squirmers, in which interactions are dominated by pairwise collisions. This could be accomplished either using look-up-tables (e.g. results of Figs. 2-4), or through developing a collision integral formulation [48] akin to the study of rarefied gases.