## Abstract

In this study, we consider suspensions of swimming microorganisms in situations where we might expect the promotion of two-dimensional flow, such as within thin fluid films. Given that two-dimensional, inertialess flows are notoriously long-ranged (although not afflicted by Stokes paradox in the case of self-motile bodies), this raises interesting questions around the care which must be taken with the semi-dilute assumption in such situations. Adopting the prototype squirmer as a model of a swimming microorganisms of the type previously considered, we find that although the flowfield decays algebraically with the characteristic separation distance between microorganisms, there remains a finite interaction between the squirmers even at asymptotically large distances. This finding is further borne out by asymptotic analysis, which confirms that the limiting form of the far-field interaction depends solely upon the relative orientation between the microorganisms. Those which swim in the same general direction are seen to experience very large lateral displacements (many times the size of the displacements experienced owing to interactions between less well-aligned swimmers). This clearly has potential implications for very dilute suspensions in which squirmers become broadly aligned in their swimming direction (e.g. during chemotaxis). We show that hydrodynamically enhanced cell spreading, previously reported for denser suspensions, can persist even at extreme dilutions. Moreover, we demonstrate that this induced spreading can continue in the presence of potentially decohering Brownian effects.

## 1. Introduction

The study of swimming microorganisms is an area that attracts a great deal of attention, no doubt in large part owing to the wonderful demonstrations of collective behaviour reported within the literature [1–5]. As might well be expected, many theoretical studies into these fascinating phenomena have ensued. One of the most well-studied and striking examples of self-organization in microorganism populations concerns bioconvection, where cells that have a tendency to swim upwards (owing to their low centres of mass) create unstable density stratifications, which eventually break down to produce patterns very reminiscent of those seen in thermal instabilities, e.g. the Rayleigh–Benard instability [6]. This same bottom-heaviness can also cause certain microorganisms in vertical pipeflow to aggregate at either the centre of the pipe, or its walls, depending upon the flow's direction [7]. In both cases, the patterns are expected to persist even in semi-dilute suspensions, where the cell concentrations are sufficiently high that the suspension can still be viewed as a continuum at the macroscopic level, but not so high that cell–cell hydrodynamic interactions are expected to be important. This is by virtue of the fact that in these examples the localization of microorganisms stem from each individual cell responding in a similar fashion to an external stimuli, such as gravity or flow shear, rather than via any fluid-mediated coordination between the cells. Under such circumstances, cell conservation equations can be used to great effect [7]. In these models, a diffusion term acknowledges that any population will contain a spectrum of cell physiologies and mobilities, whereas an advection term incorporates the effect of hydrodynamic forces upon an individual cell (for example, their possible reorientating influence). However, in general, no account is taken of cell–cell hydrodynamic interactions in these cell conservation approaches.

Because cells must be force- and torque-free in this inertialess environment, a force–dipole representation of the flow generated by a swimming cell has often been considered appropriate at reasonable distances from the cell body, and recent experimental measurements do seem to offer some support to this hypothesis (depending upon the type of microorganism) [8,9]. As well as generating propulsion, these flows exert hydrodynamic forces on their neighbours which, in turn, can generate additional cell transport. In this way, the motions of individual swimming cells in a colony can become coupled. Numerous studies have demonstrated that this can lead to superdiffusive transport in the cell population [10], including the seminal experimental work of Wu & Libchaber [4]. More recent experiments [11,12] and theoretical studies [13,14] have revealed a linear relationship between enhanced diffusion and the flux of swimming cells, with the strength of this relationship seen to be a function of distance to any nearby solid boundaries [14]. This is important given the propensity of many swimming microorganisms to aggregate close to surfaces. Interestingly, because tracer particles undergo closed loops when far from a swimmer on an infinite straightline trajectory, there is some discussion in the current literature as to the dominant mechanism by which net displacement can occur. Possible explanations include fluid entrainment by the swimmer [15,16], as well as finite length trajectories (e.g. owing to run-and-tumble swimming) [12]. In addition to the flow singularity models of swimming microorganisms predominantly used in these studies, there have also been computational studies that explicitly treat the fluid dynamics using a lattice Boltzmann scheme [17], as well as more rule-based simulations [18] that use models of flocking behaviour [19].

Obtaining a continuum-level representation of the suspension dynamics, one that builds upon the earlier cell conservations models and can reproduce many of complexities observed in suspensions of microorganisms, is understandably a much-sought-after goal. One approach has been to incorporate additional active stresses into the flow equations, which account for the effect of the motile cells upon the fluid dynamics of the suspension. These flow equations are typically coupled to equations that describe the swimming orientation; some of these are phenomological in nature [20], others incorporate more specific orientation kinetics [21], whereas others invoke free-energy arguments [22,23]. Such models have been able to reproduce a number of observed phenomena, including the ability of motile cells to lower the effective viscosity of the suspension which, in turn, can degrade the stabilizing capabilities of the viscous solvent and allow perturbations to the suspension's configuration to grow in time. Such instabilities can lead to a *phase transition*, whereby a suspension of microorganisms with an isotropic distribution of swimming directions can switch to a nematic state with strong orientational order. It is thought that this is an important means by which the exotic types of collective motion, as seen in a number of experiments [3,5,24], are generated. Studies which exploit the analogies with (active) liquid crystals have demonstrated emergence of shear-banding and rolls, which transition to more disordered states as swimming becomes more vigorous [22]. This behaviour becomes more complex still as the number of degrees of freedom of the system increases [23], with periodic and chaotic structures observable. Importantly, the nature of these active stresses is seen to depend upon whether the microorganisms are *pulled* or *pushed* through the fluid by their swimming mechanisms, and different configurations are seen to be either stable or unstable, depending upon this distinction. These distinctions have also been reported to affect the bulk behaviour of a suspension, which can vary from solid-like behaviour, to that of a superfluid [25]. Some of these collective complexities have also been reproduced in studies that adopt the type of coarse-graining approach commonly used in non-equilibrium statistical mechanics [26,27,28]. These suggest that both steric and hydrodynamic effects are important in creating orientational order in dense suspensions of rod-like suspension. In addition, that long-range hydrodynamics are essential to the destabilization of orientational order.

However, the ability of microorganisms to alter the effective viscosity of a suspension is not restricted to rod-like bacterial microorganisms that align with the local flow. It has also been predicted in the so-called *squirmer* models of ciliated microorganisms [29,30], such as *Opalina*, which are largely spherical in shape. Squirmer models aim to emulate, in a mathematically concise manner, the means by which such microorganisms propel themselves through fluid through a synchronized beating of their surface coating of cilia. These use a mathematical description of the effective flow surface condition as formulated by Blake [31]. Such squirmer models take a time-average of the flow over multiple cilia beats, and have been shown to be a reasonable theoretical approximation to certain ciliated microorganisms [32]. In situations where the interaction time is comparable to this average time, however, it has been shown that unsteady effects owing to the beating could be important [33,34].

For such spherical swimmers, it has been recently shown (through explicit full simulation) that hydrodynamic interactions alone are sufficient to generate orientational order [35]. Moreover, the effective viscosity of a suspension of bottom-heavy squirmers is seen to increase or decrease, depending upon whether the background shear flow is horizontal or vertical [36]. Certainly, a monolayer of a relatively modest (although still computationally challenging) number of such squirmers has been demonstrated to exhibit coherent structures [37]. Ishikawa *et al.* [37] do remark, however, that the two-dimensional nature of their monolayer configuration may be enhancing the strength of the correlations seen in their simulations, because the same degree of coherence was not observed in simulations of three-dimensional suspensions [38]. It is worth noting that in these simulations the monolayer exists in unbounded fluid, hence the flows generated by each squirmer are fully three-dimensional. However, in situations where a monolayer-like suspension comes about owing to close proximity to a nearby surface (as could be the case in the experiments of [3]), then it seems reasonable to expect that the flows generated by the squirmers might well be predominantly two-dimensional. Under such two-dimensional flow assumptions, there is evidence that the effective viscosity of a suspension of microorganisms can be lowered by their swimming motions [39] (albeit with a different type of microorganism model than that considered by Ishikawa *et al.* [29]). Of course, considering two-dimensional flow regimes brings with it certain theoretical advantages. Not least of which is the accessibility of complex-variable methods, and a number of studies have considered this two-dimensional regime. For example, such techniques have recently been used to determine the interaction of a squirming microorganism near a wall, where it is seen that squirming motions which would not produce propulsion in unbounded flow, can lead to motion owing to the presence of the squirmer's image system [40,41]. The motions predicted are seen to bear close resemblance to models and experiments of two-sphere-squirmers near to a solid surface [42,43], perhaps giving an indication that three-dimensional flow is suppressed to some extent by the wall's presence. Extensions to these studies, that exploit conformal mapping techniques, include consideration of squirmer swimming above a gap in a wall [44], beneath a free surface [45] and near to a corner [46]. There has also been some consideration of the effect of squirmer swimming on the effective viscosity of a suspension in a two-dimensional setting, although these have predominantly treated the collective interactions as being decomposable into a collection of pairwise interactions [13,47].

Of course, two-dimensional low Reynolds number flows are notoriously long-ranged (i.e. slow decay that scales with the inverse of distance from a moving force-free body). This, in turn, raises interesting questions around the modelling of some suspension types just discussed (e.g. suspensions within thin films bounded by a solid surface [3,11,14], or within suspended thin films [1,2,4]), where we might expect strongly two-dimensional flows. Hence, in this study, we examine the extent of the long-range two-dimensional hydrodynamic influence between squirmers, and discover that non-negligible cell–cell interactions have the potential to occur even at asymptotically large separations (§3*a*). However, the extent to which these hydrodynamic effects can influence bulk cell densities must be weighed against the decohering abilities of the intrinsic noise in the system owing to randomness in the swimming actions (which are greater in magnitude than pure thermal noise) [9]. Hence, we simulate squirmer suspensions, incorporating both hydrodynamic and stochastic effects, and find that in the presence of sufficiently many cells, hydrodynamic interactions do still leave a non-negligible imprint on the suspension (§3*b*).

## 2. Squirmer hydrodynamics

We begin by determining the flow generated by a single two-dimensional squirmer, of the type previously treated by Crowdy [40]. We will, however, allow for a more general form of squirming motion than has hitherto been considered. Working in the complex plane *z**=*x**+*iy** (asterisks denote dimensional quantities, and ** u***=(

*u**,

*v**), where

*u** and

*v** are the horizontal and vertical components of flow, respectively. We non-dimensionalize by scaling lengths on typical squirmer radius,

*U**

_{0}. This leads to a natural timescale,

*ψ*, to describe the two-dimensional Stokes flows generated by the squirmer (dropping asterisks on non-dimensional quantities):

*f*(

*z*) and

*g*(

*z*)

We initially choose to work in a swimmer frame of reference where the squirmer is positioned at the origin and orientated along the *x*-axis. It also proves convenient to define a local polar coordinate system (*r*,*θ*), with origin at the squirmer's centre. Hence, in this frame of reference, the swimming direction corresponds to *θ*=0. (Subsequently, determining flows generated by squirmers with arbitrary positions and orientations is then simply a matter of applying the relevant rotations and translations.) The squirming motion generated by each swimmer is specified by prescribing the radial, *U*_{r}, and azimuthal, *U*_{θ}, components of flow velocity at the surface of a squirmer in terms of the following Fourier decomposition:
*a*_{k}, *b*_{k}, *c*_{k} and *d*_{k} are real-valued constants. Expressing these in terms of a complex velocity (noting that *u*+*iv*=*e*^{iθ}(*u*_{r}+*iu*_{θ})) gives
*z*=*e*^{iθ}, and so *f*(*z*) and *g*(*z*), which produce a flow satisfying these surface conditions, and which tend to uniform flow in the far-field (in the swimmer frame of reference). These need to be chosen in a way which leaves the swimmer force- and torque-free (accomplished by the inclusion of algebraic terms only, i.e. no logarithmic terms). We therefore pose
*p*_{k} and *q*_{k} are complex-valued coefficients. Substituting into (2.1), we obtain
*U*+*iV* (which is precisely opposite to the uniform far-field flow in this frame of reference) is given by *U*+*iV* =*p*_{0}=−*γ*_{1}.

We have therefore determined the flow generated by the most general possible squirming motion, although, in what follows, we shall consider the simplest possible generated viscous flows that allows the microorganism to propel itself. These decay such as either 1/*r*, 1/*r*^{2} or 1/*r*^{3}. Specifically, we shall require only at most the first four coefficients in the tangential surface flow to be non-zero (i.e. only the first four Fourier modes in (2.3) are required). Under these circumstances, in the squirmer frame of reference, the complex velocity and vorticity take the following respective forms
*d*_{1}=*d*_{2}=2 (all other squirming coefficients zero), which combines elements of the self-motile potential swimmer considered by Thiffeault *et al.* [13], and the viscous treadmillers of Crowdy *et al.* [40,41,45] (which are not self-motile in an infinite fluid). This choice of squirming coefficients leads to a *puller* [48] (e.g. *Chlamydomonas*, as opposed to a *pusher*, which propels itself from behind).

In a suspension, each squirmer not only moves through self-propulsion, but is also transported by the flowfield generated by its fellow swimmers. Such coupling between the individual squirmers, via the fluid medium, is a well-known mechanism by which collective behaviour can occur [49]. The full treatment of the *N*-body hydrodynamic interaction poses a formidable computational challenge (Ishikawa & Pedley [36,37] have been able to simulate the full hydrodynamic interactions between 64 squirmers). Hence, in §3, we consider a simple dynamical system that aims to capture the consequences of the hydrodynamic coupling on a dilute squirmer suspension using some simplifying assumptions.

## 3. Squirmer interactions

We consider a suspension of *N* identical squirmers. (Of course, actual suspensions will consist of microorganisms with slightly varying sizes and motilities. These effects can be incorporated into the following formulation with relative ease, although we do not do so here.) In a laboratory frame of reference, we use *s*_{j}=*x*_{j}+i*y*_{j} and *θ*_{j} (*j*=1…*N*) to denote the cell centre and orientation, respectively, of the *j*th squirmer. Because we are focused in this study on suspensions in which individuals swim many body lengths apart, we approximate the flow generated by a collection of squirmers as the superposition of the flowfields generated by each individual cell. We also perform some supporting boundary element method (BEM) computations, to test the limits of this assumption.

#### (i) Dynamical system

In a laboratory frame of reference, each squirmer swims with velocity *θ*_{j} under its own motile force. It also generates flow velocity and vorticity fields, which are readily derived from (2.7) (i.e. velocity and vorticity in the squirmer frame). As such, squirmers are also transported by the flows generated by their neighbours. We will also allow for Brownian-like motion which stems from inherent randomness in the squirmer's swimming actions (and might also encompass thermal effects, although these are often considered less significant).

The squirmer motions can consequently be described by the following dynamical system:
*k*,*j*=1…*N*), where *w*=*W*−*γ*_{1} (see (2.7a)) and *σ*_{T} and *σ*_{R} are the (non-dimensional) translational and rotational diffusivities, respectively. The *x*-axis, into the laboratory frame of reference in which the squirmers are orientated with angles *θ*_{k} and *θ*_{j}, respectively.

The above system of stochastic differential equations (SDEs) is integrated using a fourth-order Runge–Kutta scheme (rather than an Euler–Mayuma scheme, which is only first-order accurate in time) to ensure that the dynamics are sufficiently resolved. The coefficients chosen for the scheme were those suggested by Burrage & Burrage [50], which give the standard Runge–Kutta scheme in the deterministic limit (and predict trajectories which agree with those obtained using Matlab's ODE45 solver). The form of this stochastic Runge–Kutta scheme has increments
*s*_{n}=(*s*_{1}(*n*Δ*t*),…,*s*_{N}(*n*Δ*t*)), *θ*_{n}=(*θ*_{1}(*n*Δ*t*),…,*θ*_{N}(*n*Δ*t*)) with Δ*t* the size of the timestep. The choice of coefficients for this particular Runge–Kutta scheme is given as [50]: *b*_{1}=−0.724, *b*_{2}=2.702, *b*_{3}=0.224, *b*_{4}=1.757, *b*_{5}=1, *b*_{6}=−2.918524118. Here, ** F**(

*s*_{n},

*θ*_{n})=(

*F*

_{1},…

*F*

_{N}) and

**(**

*G*

*s*_{n},

*θ*_{n})=(

*G*

_{1},…

*G*

_{N}), where

Finally, *α*_{0}≈7.8×10^{−9} [50])

Simulating full *N*-body interactions is nevertheless computationally expensive, even when simulating numbers of cells (e.g. *N*=100, 1000) which are orders of magnitude less than found in typical microorganism suspensions. In semi-dilute suspensions, however, it is usually accepted that the full *N*-body problem can be decomposed into a sequence of two-body interactions, by virtue of the large separations between squirmers, thereby offering substantial computational savings. In the next section, we therefore consider pairwise interactions between two squirmers (in the frame of reference of a squirmer). In doing so, however, we discover that when the flows can be expected to be largely two dimensional, reconstruction of the collective dynamics from pairwise interactions has only very limited applicability. Hence, we ultimately find ourselves needing to conduct full *N*-body simulations.

#### (ii) Boundary element computations

The superposition assumption supposes that the corrections to surface flow conditions from neighbouring squirmers has negligible influence at sufficient separation distances. In order to gauge the typical separations required for this to hold true, we also conduct some BEM simulations that take into account the full hydrodynamic interactions.

If the surfaces of the *N* squirmers at a given time *t* are denoted by *m*=1…*N*), then the flow on the surface of the *m*th squirmer can be expressed in boundary integral form as [52]
*l*,*m*=1…*N*, *i*,*j*,*k*=1,2), where ** f**=

**⋅**

*σ***is surface traction, and**

*n***is the known squirming velocity. Here,**

*U*

*V*_{m}and

*Ω*

_{m}are the translational and angular velocities of the

*m*th squirmer, respectively, that are induced by the hydrodynamic interactions with other cells. Here,

*S*

_{ij}and

*T*

_{ijk}are the two-dimensional steady Stokeslets and Stresslets, respectively [52]

*V*_{m}and

*Ω*

_{m}are determined through the fact that the squirmers remain force-free and torque-free, i.e.

*V*_{m}and angular velocities

*Ω*

_{m}of each squirmer at any instant in time. Their centre positions,

*s*_{m}=(

*x*

_{m},

*y*

_{m}), and orientations,

*θ*

_{m}, are then numerically evolved by time-stepping the following system

### (b) Deterministic pairwise squirmer interactions

We begin by considering the deterministic interaction of a squirmer pair, which will enable us to highlight some fundamental features of the hydrodynamic coupling. For pairwise interactions, the dynamical system given in (3.1) takes the form

The ability of this dynamical system to capture pairwise hydrodynamic interactions is explored in figure 2. Here, trajectories computed using the dynamical system (3.6, dashed lines) are compared with those made using BEM computations (3.5). It can be seen that for the squirmer motions primarily considered here, beyond distances of about 20 cell lengths the dynamical system performs well. For different squirmer motions, similar to those considered by Ishikawa *et al.* [29] (*d*_{1}=1, *d*_{2}=5), slightly greater separations are required.

We consider two ostensibly very different pairwise interactions. These are labelled interaction I and interaction II, and correspond to pair of squirmers that are initially *O*(10^{3}) and *O*(10^{4}) body lengths apart, respectively. Table 1 gives the specific initial values of squirmer locations and orientations. See also figure 3 for a sketch. We can see from figure 4, which plots the magnitude of the separation between the two squirmers |*s*_{2}−*s*_{1}| over time, that in both cases there remains an order of magnitude difference in squirmer separation throughout their swimming trajectory (interaction I:

Although the simulations in the laboratory frame of reference serve to illustrate the persistence of hydrodynamic interactions at large separations, the two interactions themselves appear quite different in form. However, there are some dynamical similarities between the pairwise interactions, which are most easily seen by returning to the frame of reference of one of the squirmers. Hence, we shall take squirmer 1 to be fixed at the origin, and orientated along the *x*-axis. We then solve for the relative displacement *s*=(*s*_{2}−*s*_{1})*e*^{−iθ1}, and relative orientation *θ*=*θ*_{2}−*θ*_{1} of the second squirmer (squirmer 2). See figure 6*b* for an illustration. Under these circumstances, we obtain from (3.6) the reduced dynamical system
*q*=*e*^{iθ}. The time evolution of both relative displacement and relative orientation is plotted in figure 7. The lines correspond to the predictions from the laboratory frame (3.6), whereas the markers denote the predictions obtained directly from the squirmer frame (3.7), with good agreement seen between the two. It is seen that while the relative orientation does not change significantly at such separations, relative displacement does.

In this squirmer frame of reference, it is now convenient to define the impact parameter, *b*(*t*), which allows us to uncover the above-mentioned dynamical similarities. The value of the impact parameter at any time is defined to be the distance of closest approach between the two squirmers, were the current instantaneous trajectory to be maintained (see figures 6*b* and 8 for illustrations), and is a commonly used quantity in kinetic gas theory [53]. A simple geometrical argument shows that it can be found by projecting the complex-valued trajectory, *s*(*t*), onto the unit vector perpendicular to the direction of motion, i.e.

If we now examine the time evolution of the impact parameter provided in figure 9, then we observe that the jump in its value, *b*(0)|≫1. The choice of sign corresponds to whether the initial impact parameter *b*(*t*=0) is positive or negative. This limit is overlaid as a dashed line in figure 9, where it is seen to agree with the numerically calculated jump, *θ*(0)=2*π*/3).

In order to examine the conditions required for (3.9) to be applicable, in figure 10, we plot the shifts in impact parameter, *θ*=0,2*π*. There are also steep gradients in the value of the asymptote near these parallel relative orientations. This is more clearly illustrated in the bottom images of figure 10, where profiles are taken through the impact parameter shift surface at *b*(0)=100 and *b*(0)=−100. Here, we also overlay the predictions of the asymptotic expression for impact parameter shift, as given by (3.9), which is seen to accurately predict this limiting shift over the full range of initial relative orientations. The implication is that squirmers which are initially orientated nearly parallel experience much greater hydrodynamic displacements than those that are less aligned. We shall see in the next section that this has implications for suspensions of squirmers that swim in broadly similar directions, for instance, as might be expected in the presence of nutrient gradients (i.e. during chemotaxis).

We also see in figure 11 that the long-range hydrodynamic effects are not restricted to flowfields which decay such as 1/|*z*| (as is the case when *d*_{1}=*d*_{2}=2). They also occur for squirming flowfields that decay like 1/|*z*|^{2} and 1/|*z*|^{3}.

### (c) *N*-body squirmer interactions

We saw in the previous section that the deterministic pairwise hydrodynamic interactions do not diminish with increased separation between the squirmers. We also found that two squirmers swimming in a similar direction can experience especially large displacements away from their non-interacting (i.e. straightline) trajectories. We might therefore wonder whether these effects can induce the kind of hydrodynamic spreading previously reported within suspensions [38,11], but at concentrations here perhaps previously considered too low for such effects to be relevant.

First, in figure 12, we examine the predictions for cell positions at time *t*=2000 within two different suspensions. One containing *N*=10 squirmers (left image), and another containing *N*=100 squirmers (right image). These have been calculated using the dynamical system (3.1), and using either successive pairwise interactions (crosses), or full *N*-body simulations (circles). As the number of squirmers increases, the accuracy of the successive pairwise interactions in this two-dimensional regime can be seen to reduce. Given the population sizes that we shall be considering shortly, we therefore choose to compute the full *N*-body interaction problem, rather than decomposing into a sequence of pairwise interactions.

We also examine the validity of using the dynamical system (3.1), by comparing its prediction with those of full BEMs for suspension of 6–12 squirmers. From figure 13, it can be seen that for squirmers of the type we have been considering (*d*_{1}=*d*_{2}=2), the dynamical system does a reasonable job of approximating the trajectories for suspensions where individuals are more than about 20 cell lengths apart (this number may vary slightly for different squirming motions, as was the case for pairwise interactions; cf. figure 2). Given that our main interest here is in suspensions of squirmers separated by distances much greater than some tens of body lengths, we feel confident in the reliability of the dynamical system description.

In figure 14, we show the positions of *N*=300 squirmers at time *t*=2000. These squirmers were initially uniformly distributed in a box 1000×1000 cell lengths in size. Two independent *N*-body simulations were computed for each set of initial conditions. One that included hydrodynamic interactions (dots), and another which did not include such interactions (crosses). The left-hand figure corresponds to simulations conducted with initial isotropic swimming directions (i.e. *θ*(0) uniformly distributed between 0 and 2*π*), which we can observe leads to little hydrodynamically enhanced spreading. However, given that in §3*a*, we observed much greater hydrodynamic displacements between pairs of swimmers when their swimming orientations were more closely aligned, we recompute the final positions from an initial state where squirmers have a swimming direction biased towards *θ*=0. In this situation, we now see from the right-hand figure that there is appreciable hydrodynamically induced spreading.

Nonetheless, before drawing any conclusions with regard to hydrodynamic spreading in very dilute suspensions, we should also perhaps investigate whether such hydrodynamic interactions are sufficiently strong to overcome the decohering influence of Brownian effects. In order to determine whether this is the case, we consider suspensions of squirmers uniformly distributed within a domain that is 1000×1000 cell lengths in size. In line with the findings above, we also concentrate on the case where the squirmers are swimming in a preferential direction, perhaps owing to a chemotactic gradient. Hence, swimming orientations are initially normally distributed about *θ*=0, with a standard deviation of 0.1. Consequently, we solve the stochastic differential equations (3.1) for population sizes of *N*=100, *N*=500 and *N*=1000. We assume that the (non-dimensional) rotational diffusivity of these squirmers is similar to that measured in the experiments of Vladimirov *et al.* [54,38] on *Chlamydomonas nivalis*, giving us *t*=2000 and *t*=4000, in order to gauge how any hydrodynamic influences evolve over time. From figure 15, we note that for *N*=100, Brownian effects do seem to dominate over the hydrodynamic spreading noted in the earlier deterministic simulations (figure 14). However, the long-range hydrodynamic interactions are additive, and so as the number of squirmers increases, hydrodynamic spreading appears to reassert itself. For a suspension of *N*=500 squirmers, a clear difference in the distributions of squirmers subject to hydrodynamic interactions emerges. This is a trend which strengthens still further when the population is increased to *N*=1000 squirmers. Hydrodynamic interactions are seen to retard the progress of the squirmers, leading to a greater tail in the spatial profile. (Given that flows generated by singularity solutions may be expected to decay more rapidly depending upon the particular swimming regime [55], in section 3 of the electronic supplementary material, we also consider the impact of different squirmer motions, with flowfields that decay such as 1/|*z*|^{2} and 1/|*z*|^{3}. In addition, in all of these simulations, we have made the usual assumption that rotational diffusivity is more significant than translational diffusivity, which we have taken to be negligible. However, to demonstrate that the above observations are insensitive to the introduction of significant translational diffusivity, in section 3 of the electronic supplementary material we also consider the situation where rotational and translational diffusivities are comparable.)

## 4. Discussion

We have examined the hydrodynamic interactions in a dilute suspension of squirmers for a regime where the flow can be considered to be predominantly two-dimensional. In doing so, we have uncovered the unexpected result that the pairwise interaction between two squirmers does not become negligible even at asymptotically large squirmer separations. This is in spite of the fact that the flowfield actually decays with distance. As the suspension is made more dilute, the relative distance a squirmer is displaced by a hydrodynamic interaction with any single neighbouring squirmer becomes less significant. However, the hydrodynamically induced displacements accumulate, because, in this regime, no squirmer is too far away to exert a hydrodynamic influence. This therefore seems to suggest that the importance of hydrodynamic interactions in this very dilute regime depends not only upon average separation between microorganisms (e.g. the suspension concentration), but also the population size. For similar reasons, we have seen that even in very dilute suspensions, the *N*-body squirmer interactions cannot necessarily be decomposed into a succession of pairwise interactions.

Population size was also seen to be a factor in whether hydrodynamically enhance spreading in very dilute suspensions can persist in the presence of Brownian effects (both biological and thermal). For relatively small populations, rotational diffusivity dominates over any hydrodynamic spreading. However, as the number of squirmers increases, so too do the cumulative hydrodynamic interactions. In sufficiently large populations, they are able to leave their imprint upon the bulk density of a suspension. For the highest population suspensions considered here (*N*=1000, which is still quite modest compared with cell populations in actual suspensions), the hydrodynamics induce an enhanced cell spreading even in the presence of Brownian noise. Although the rotational diffusivities considered here lie towards the lower end of those measured experimentally [54], we expect that higher populations (greater than we are currently able to feasibly compute) will likely lead to an hydrodynamic imprint even for microorganisms with greater intrinsic rotational diffusivities. We have also assumed that the rotational diffusivity is isotropic, taking into account a potential swimming bias only through the initial conditions. However, it may be more rational to incorporate this swimming bias into the random behaviour, so that it persists over time. In that case, we might expect the hydrodynamics to be more significant still, given that the pairwise interactions seem to suggest that greater alignment leads to greater long-range hydrodynamic interactions.

It is worth noting that although Drescher *et al.* [9] have recently argued that in dense suspensions the radius of influence of hydrodynamics about a microorganism is actually very small (less than 10 *μm*), and hence that self-organization may be collision-dominated and steric effects potentially important [27], their scaling argument assumes much shorter interaction times than is the case in this very dilute, two-dimensional regime. In addition, the scaling argument only considers the instantaneous flowfield. The squirmer displacements, however, are integrated quantities. Our full stochastic simulation seems to indicate that our dilutions (*ϕ*≈10^{−3}) and interaction times are such that hydrodynamics remain important in this regime. It is important to note, however, that we are not suggesting here that these long-range hydrodynamic interactions are capable of producing the organizational instabilities and phase transitions seen in much denser microorganism suspensions. The effects on bulk density predicted here are on altogether longer spatial and temporal scales.

In terms of future directions, there is an interesting current debate as to the whether advection or entrainment are more important to the enhanced diffusion seen within suspensions, with a suggestion being made that the answer may depend upon whether the flow is predominantly two-dimensional, or three-dimensional [12,15,16]. It might be interesting to pursue this question further using the type of model squirmers used in this study. In addition, given the increasing computational difficulty as the population size increases (each *N*=1000 stochastic realization simulated here took over a week's computational time across six CPUs), following Baskaran *et al.* [28] and others, a fruitful future approach may be to incorporate these hydrodynamic expressions into a coarse-grained continuum description for the suspension.

## Acknowledgements

We are grateful to Prof. John Sader (University of Melbourne) and the referees for some very useful suggestions. Simulations were run on the NZ eScience Infrastructure (NeSI) facility.

- Received August 5, 2013.
- Accepted April 8, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.