2.1. Single swimmer model

Micro-scale swimmers such as E. coli and Chlamydomonas use flagella, whip-like attachments to the cell body, to swim in a viscous fluid ambient. They differ in their swimming mechanism, however, and are classified as pushers and pullers, respectively. Pushers have a propulsive part (the posterior flagellar bundle in E. coli) which pushes the head of the organism through the fluid, leading to an extensile intrinsic dipole (see figure 2 a). Pullers have the propulsive portion (the pair of anterior flagella in Chlamydomonas) pulling the head through the fluid with the aid of a breast-stroke-like motion, leading to a contractile dipole (see figure 2 b). Typical pushers and puller velocity fields are shown in figure 2(e,f). We model these swimmers as (infinitely) slender rods which, given the large aspect ratios of typical pushers ( ${\sim}6$ for E. coli), is a fairly good approximation. The intrinsic force distribution is taken to be a line distribution of Stokeslets (Leal Reference Leal2007) with the Stokeslets being parallel (anti-parallel) to the swimming direction on the passive (propulsive) portion. The distribution is parameterized by ${\it\alpha}_{1}=L_{H}/L$ , where $L_{H}$ and $L$ are the head and total swimmer lengths, and is of the form $f(s)\boldsymbol{p}$ , with:

for pushers, where $F_{D}$ is the drag force on the head. Here, $\boldsymbol{p}$ denotes the swimmer orientation, and $s$ is the coordinate along the swimmer axis measured from the centre (see figure 2 c,d). The model above does not aim to elucidate the actual swimming mechanism, but merely assumes one to exist and that it leads to a non-trivial swimming velocity and force distribution. Thus, features of the velocity field specific to a given propulsion mechanism are neglected. For instance, in E. coli, the counter-rotation of the head and the propelling flagellar bundle leads to an additional torque dipole or a vortlet singularity (Ghose & Adhikari Reference Ghose and Adhikari2014), and the associated $O(1/r^{3})$ velocity field is neglected here (a reasonable approximation for slender swimmers since the singularity strength must scale with the much smaller head diameter). The puller model is the same as (2.1), but with the force distribution reversed. Again, our intention in modelling pullers is not to describe specific organisms such as Chlamydomonas, whose near-field velocity may be complicated by unsteady stroke kinematics and the presence of closed streamlines (Guasto, Johnson & Gollub Reference Guasto, Johnson and Gollub2010; Drescher et al. Reference Drescher, Dunkel, Cisneros, Ganguly and Goldstein2011), but to study the effect of the contractile nature of the dipole on collective motion. The main role for the puller simulations is to serve as a baseline for contrast.

Figure 2.

Schematic of the single swimmer model. (a,b) depict the Stokeslet distribution used to model a pusher and puller, respectively. The corresponding variations of the linear force density along the swimmer axis are shown in (c,d). (e,f) Show the dipole velocity field at large length scales; the dipole has an extensile character for the pusher and is contractile in nature for the puller.

The results presented here are for fore–aft symmetric swimmers $({\it\alpha}_{1}=1/2)$ with a Heaviside force density profile. The value of $F_{D}$ in (2.1) is obtained from equating the drag on a spheroidal E. coli head to that on a slender rod moving with the same speed. This leads to an aspect ratio of $6$ for the fore–aft symmetric swimmers simulated. Changes in ${\it\alpha}_{1}$ do not qualitatively alter the results. We have also considered a sinusoidally varying force density, and again find no qualitative differences.

2.2. Hydrodynamic interactions in a swimmer suspension

The swimmers, as they swim along, ‘stir’ the suspending fluid, and thereby, affect the motion of other swimmers. In the low Reynolds number regime, relevant to most micro-scale swimmers (for E. coli, $L\approx 12~{\rm\mu}\text{m}$ and $U\approx 10~{\rm\mu}\text{m}~\text{s}^{-1}$ , so that the Reynolds number, $Re={\it\rho}UL/{\it\eta}$ , is $O(10^{-4})$ for swimming in an aqueous medium. The $Re$ ’s characterizing the coherent structures, such as the vortices and jets, that dominate the collective dynamics, continues to remain small (Tuval et al. Reference Tuval, Cisneros, Dombrowski, Wolgemuth, Kessler and Goldstein2005)), the governing equations are the Stokes equations with an additional forcing term due to the swimmers and the continuity equation for an incompressible fluid. This set of equations, in dimensionless form, can be written as:

where $p$ and $\boldsymbol{u}$ are the pressure and velocity fields, and

is the forcing term due to the suspended swimmers. Here, $\boldsymbol{x}_{m}$ and $\boldsymbol{p}_{m}$ are the position and orientation of the $m$ th swimmer. We have used the following scales for non-dimensionalization: the length scale is $l_{c}=L/2$ ; the time scale is $t_{c}=L/(2U)$ ; and forces are scaled by $F_{c}={\it\mu}UL/2$ . Henceforth, we will be using dimensionless variables with the same notation as before. In dimensionless terms, the swimmer length is $2$ and its isolated swimming speed is unity.

The simulation domain is a cubic unit cell of side length $L_{box}$ having $N$ swimmers with centres at $\boldsymbol{x}_{1},\boldsymbol{x}_{2},\ldots ,\boldsymbol{x}_{N}$ , and with orientations $\boldsymbol{p}_{1},\boldsymbol{p}_{2},\ldots ,\boldsymbol{p}_{N}$ , supplemented by periodic boundary conditions. With this arrangement, the summation term in (2.4) can be written:

where the inner summation denotes a sum over the simulation box ( $\boldsymbol{l}=0$ ), and the outer summation involving the lattice vector $\boldsymbol{l}$ is over the periodic images. For numerical purposes, the swimmer line force distribution is specified using point forces at discrete points along the swimmer axis and the integrals involving the force distribution are computed using Gaussian quadrature. With each swimmer being discretized as $M$ quadrature points, (2.5) can be written as:

where $f(s_{j})\boldsymbol{p}_{i}$ is the $j$ th Stokeslet on the $i$ th swimmer in the simulation box, and the $w_{j}$ ’s are the quadrature weights. Absorbing the quadrature weights into the force density (Saintillan, Darve & Shaqfeh Reference Saintillan, Darve and Shaqfeh2005), the momentum equation, using (2.6), is given by:

One has to, in effect, solve the Stokes equations for $N\times M$ interlacing simple cubic arrays of point forces. A solution for a single simple cubic array of point forces was first derived by Hasimoto (Reference Hasimoto1959). Hasimoto’s approach has been extended to the case of interlacing arrays in the context of simulating suspensions of slender fibres (Mackaplow & Shaqfeh Reference Mackaplow and Shaqfeh1998; Butler & Shaqfeh Reference Butler and Shaqfeh2002). Using a different approach, interlacing cubic arrays of finite-sized spheres corresponding to a random array in the original unit cell were analysed by Beenakker (Reference Beenakker1986), and this alternate approach has been adopted by Brady and co-workers in their Stokesian dynamics algorithm (Brady & Bossis Reference Brady and Bossis1988; Brady et al. Reference Brady, Phillips, Lester and Bossis1988). These two approaches are, of course, equivalent (Krishnamurthy Reference Krishnamurthy2014), and we follow Hasimoto (Reference Hasimoto1959).

As indicated above, the derivation of the disturbance velocity field in an infinite spatially periodic swimmer suspension, accounting for the re-normalization of the long-ranged hydrodynamic interactions, is now standard and follows the literature for passive suspensions (Brady et al. Reference Brady, Phillips, Lester and Bossis1988). Briefly, the velocity and pressure gradients are expanded in a Fourier series with the $\boldsymbol{k}=0$ terms being excluded due to the total force per unit cell being balanced by a mean pressure gradient. From a numerical point of view, one needs to exclude the $\boldsymbol{k}=0$ terms even for a suspension of force-free swimmers (these terms correspond to a specification of the particle stress in a homogeneous isotropic swimmer suspension). Once the $\boldsymbol{k}=0$ terms are interpreted correctly, an exact periodic solution to the Stokes equations, in terms of a slowly convergent series in wavenumber space, results. This series is only conditionally convergent at large $k$ , reflecting the fact that the Fourier transforms of the physical space singularities (force-multipoles) must be interpreted as generalized functions. The slow convergence in $k$ is resolved in the usual manner by using the Ewald summation technique (Ewald Reference Ewald1921; Hasimoto Reference Hasimoto1959; Beenakker Reference Beenakker1986), which involves splitting the original sum into real and Fourier space sums, with a cutoff parameter or Ewald coefficient, ${\it\alpha}$ , tuned such that the total cost of the computation is minimized. The final result for the disturbance velocity in terms of the real and Fourier sums is:

where the tensors $\unicode[STIX]{x1D63C}$ and $\unicode[STIX]{x1D63D}$ are given by:

$\unicode[STIX]{x1D644}$ being the identity tensor. The function ${\it\phi}_{{\it\nu}}$ in (2.9) and (2.10) is the incomplete Gamma function (Abramowitz & Stegun Reference Abramowitz and Stegun1972) and: The above expressions are the same as those given by Saintillan et al. (Reference Saintillan, Darve and Shaqfeh2005) in the context of sedimenting fibres or by Sierou & Brady (Reference Sierou and Brady2001) for suspensions of spheres. The main difference in our case is that the force distribution on each swimmer is the intrinsic one due to swimming, and is therefore known a priori, making the calculation of (2.8) an explicit one.

2.3. Swimmer kinematics

The coupling between the disturbance velocity and swimmer kinematics, for our slender swimmer model, is achieved through viscous slender-body theory (Batchelor Reference Batchelor1970; Leal Reference Leal2007). The swimmer force balance involves the zeroth moment of the disturbance velocity field along its length, and leads to the following swimmer translational velocity:

where $\boldsymbol{x}_{GC}$ denotes the swimmer centre, and $\boldsymbol{u}$ is the ambient field due to all other swimmers. The latter is obtained by subtracting the ‘self’ contribution from the total disturbance velocity field obtained above, since the intrinsic swimming velocity associated with this self-contribution is included separately. The torque balance involves the first moment of the ambient field, and leads to the following expression for the swimmer angular velocity:

The relations (2.13) and (2.14) are used to update the swimmer positions and orientations at each time step.

2.4. Tumbling and rotary diffusion models

A novel aspect of this work is the inclusion of intrinsic decorrelation mechanisms found in biological swimmers. Micro-scale swimmers in nature rarely swim in a straight line, even in isolation, orientation decorrelations being an integral feature of their motion. We study two such decorrelation mechanisms, tumbling and rotary diffusion, a combination of which is observed in bacterial motion.

2.5. Simulation algorithm

The initial condition is specified by initializing the swimmer positions ( $\boldsymbol{X}\equiv \{\boldsymbol{x}_{1},\boldsymbol{x}_{2}\ldots \boldsymbol{x}_{N}\}$ ) and orientations ( $\boldsymbol{P}\equiv \{\boldsymbol{p}_{1},\boldsymbol{p}_{2}\ldots \boldsymbol{p}_{N}\}$ ) in the simulation box. The positions are uniformly distributed in the simulation box and the orientations are uniformly distributed on the unit sphere. Note that the swimmers are assumed to be infinitely slender, and excluded volume considerations do not come into play in the choice of the initial configuration. In light of the known results for passive rod suspensions (Larson Reference Larson1988), we neglect their role in the dynamics also, since the volume fractions examined here are within the dilute regime. With the positions and orientations known, we can compute the velocity disturbance at the quadrature points along the length of each swimmer using (2.8). Once the velocity disturbance is known along the length of all swimmers, we can use (2.13) and (2.14) to calculate the translational and angular velocities. Formally, one may write the evolution equations as:

where $\mathscr{F}$ and $\mathscr{G}$ are known functionals of the velocity disturbance field. The positions and orientations of all the swimmers are updated using a second-order Runge–Kutta scheme (Wheatley & Gerald Reference Wheatley and Gerald1984), with a time step small enough to ensure converged results. For tumblers, the swimmers which tumble during the current time step are updated with a new orientation according to the model in § 2.4.1. For rotary diffusers, a small noise is added to the orientation vector of all the swimmers according to the model in § 2.4.2. The updated positions and orientations are used in the next time step, and so on.

3.1. Fluid velocity statistics

Figure 4(a–d) show typical fluid velocity fields and the associated number density fields in suspensions of pushers and pullers, with volume fractions of $0.05$ and $0.5$ , at a single time instant in the statistical steady state. At the lower volume fraction, the velocity fluctuations are only those due to individual uncorrelated swimmers which die off rapidly as one moves away from the particular swimmer (figure 4 a,b). In the suspension of pushers, as one moves to larger volume fractions, one begins to see coherent jets and vortices in the fluid on length scales larger than the swimmer size (see figure 4 c), similar to those observed in experiments (Mendelson et al. Reference Mendelson, Bourque, Wilkening, Anderson and Watkins1999; Dombrowski et al. Reference Dombrowski, Cisneros, Chatkaew, Goldstein and Kessler2004; Sokolov et al. Reference Sokolov, Aranson, Kessler and Goldstein2007), and in agreement with previous simulations (Saintillan & Shelley Reference Saintillan and Shelley2007, Reference Saintillan and Shelley2012). Pullers, however, continue to show an uncorrelated velocity field at larger volume fractions dominated by small-scale fluctuations corresponding to individual swimmers (figure 4 d). Although not shown here, the number density fluctuations for pushers and pullers at both volume fractions adhere closely to Poissonian statistics with $\langle N^{2}\rangle ^{1/2}\sim \langle N\rangle ^{1/2}$ ; we find no evidence of the giant number density fluctuations predicted/observed in other realizations of active matter (Simha & Ramaswamy Reference Simha and Ramaswamy2002; Narayan, Ramaswamy & Menon Reference Narayan, Ramaswamy and Menon2007) either in the stable or unstable regimes (for pushers).

The above qualitative observations are confirmed on examination of the probability distributions of the individual fluid velocity components (not shown) as well as the magnitude. Both measures indicate much larger fluctuations for pushers at higher volume fractions. In particular, the probability density for the velocity magnitude (figure 5 a,b) displays clear evidence of bulk fluid motion since, beyond a critical volume fraction, the distribution peaks at a non-zero speed of order unity. In contrast, the peak for pullers seems to occur at a much smaller value not resolved in the plots shown. The enhancement for pushers can also be seen in the plot of the mean fluid kinetic energy versus volume fraction in figure 5(c). The pusher kinetic energy curve diverges from the puller curve with increasing volume fraction, with larger boxes leading to larger kinetic energies and a steeper rate of increase with ${\it\nu}$ . This is consistent with linear stability theory in that the larger box sizes allow a greater range of unstable wavelengths, and thus a more intense fluid motion. Pullers do not display a noticeable trend with box size.

Figure 5.

Fluid velocity statistics. (a,b) Probability distribution function of fluid disturbance velocity magnitude ( $L_{box}=10L$ ). (c) Mean fluid kinetic energy plotted with respect to ${\it\nu}$ , for pushers and pullers, for three box sizes ( $L_{box}=5L,10L$ and $15L$ ). The volume fraction ${\it\nu}=1.0$ corresponds to 1000, 8000 and 27 000 swimmers for the box sizes of $5L$ , $10L$ and $15L$ , respectively.

Figure 6.

(a,b) The fluid velocity autocorrelation function, $F(t)$ , plotted for pushers and pullers at different volume fractions. (c) The correlation times for pushers and pullers plotted with respect to ${\it\nu}$ for $L_{box}=5L$ , $10L$ and $15L$ . The volume fraction ${\it\nu}=1.0$ corresponds to 1000, 8000 and 27 000 swimmers for the three box sizes $5L$ , $10L$ and $15L$ , respectively.

Further evidence of differences between pushers and pullers can be seen from plots of the Eulerian fluid velocity correlation, shown in figure 6(a,b), for pushers and pullers at different volume fractions. Figure 6(c) shows the correlation times, extracted from figure 6(a,b) by assuming an exponential decay, plotted as a function of the volume fraction (there is no rigorous basis for expecting an exponential decay since the correlation functions clearly change sign for intermediate times; the simple exponential fit nevertheless gives reasonable estimates). For pushers, the correlation time shows an initial increase with volume fraction reaching a maximum around ${\it\nu}=0.5$ . At these volume fractions, correlated fluid motions are persistent over time scales much larger than the single swimmer time scale of $L/2U$ ( $7$ times larger at ${\it\nu}=0.5$ for the $10$ -swimmer-length box). Further increase in the volume fraction leads to a decrease in the correlation times due to the dominant effect of a given point in space randomly sampling swimmer velocity fields at an increased frequency. The trend for pullers is clearly different in that puller suspensions show a monotonic decrease in the correlation times with volume fraction as seen in figure 6(b). This decrease may be attributed to the same effect which caused the decrease in pusher correlation times for higher volume fractions, mentioned above. The effect of box size can again be seen in figure 6(c) with larger simulation boxes ( $10$ and $15$ -swimmer-lengths) leading to much larger correlation times. Our focus, in trying to identify a stability threshold, will henceforth be restricted to ${\it\nu}\leqslant 0.6$ where the pusher correlation times increase with increasing volume fraction.

3.2. Swimmer statistics

The swimmer speed in a suspension deviates from its isolated value (of unity) owing to hydrodynamic interactions. To quantify this, in figure 7, the enhancement in the mean-squared swimming speed versus volume fraction, for pushers and pullers, is shown for three different box sizes. Beyond a critical volume fraction, there is a clear pusher–puller bifurcation. The pushers transition from a linear to a quadratic scaling with increasing volume fraction, while the pullers exhibit a linear increase over the entire range of volume fractions considered. The linear increase is consistent with the enhancement arising due to convection by the fluid velocity field which in turn may be estimated from an independent-swimmer argument. This latter observation is consistent with the findings of Hernandez-Ortiz et al. (Reference Hernandez-Ortiz, Underhill and Graham2009) for regularized dipole swimmers in a confined domain. Further, unlike pullers, the effect of box size on pusher speeds is noticeable in figure 7. Note that the increase in puller speed with volume fraction found here is unlike the weak decrease reported by Saintillan & Shelley (Reference Saintillan and Shelley2012), a discrepancy that may arise from the neglect of induced stresses.

Figure 7.

Mean-square swimming speed minus unity plotted with respect to ${\it\nu}$ , for pushers and pullers, for $L_{box}=5L$ , $10L$ and $15L$ .

Next, we look at the polar pair orientation correlation defined as (Saintillan & Shelley Reference Saintillan and Shelley2007):

for both pushers and pullers. The polar correlation function above distinguishes between parallel ( $C_{1}=1$ ) and anti-parallel ( $C_{1}=-1$ ) swimmer alignments. We see from figure 8 that pushers and pullers show strong positive and negative correlations, respectively, at the volume fractions examined, and for separations of order the swimmer length. These aligning tendencies likely arise from pair interactions in the stable regime, and the decrease in the near-contact value of the correlation with increasing $nL^{3}$ suggests that higher volume fractions impede the effects of pair interactions for both pushers and pullers. Since our results compare well with Saintillan & Shelley (Reference Saintillan and Shelley2007) we conclude that the intrinsic stresses alone are sufficient to capture the main effects of the pair correlations, even in the stable regime. The general features of the orientation correlations, such as the opposing signs for pushers and pullers and the secondary role of the induced stresses, are consistent with what one might expect from viscous slender body theory at leading logarithmic order. A careful examination of figure 8 also reveals subtle signatures of the instability. Pushers show stronger correlations at a given separation, larger correlation lengths and importantly, the correlation length increases with increasing volume fraction, indicative of contributions arising from long-wavelength velocity fluctuations.

Another statistical measure of interest is the degree of alignment of a given swimmer along the local extensional axis of the ‘ambient’ rate-of-strain tensor resulting from the cumulative effect of the velocity disturbances due to all other swimmers. In the case of pushers, as discussed earlier, this alignment leads to an orientational anisotropy which reinforces the velocity disturbance, leading to an exponential instability (see Subramanian & Koch Reference Subramanian and Koch2009). In figure 9, we plot $\boldsymbol{E}\boldsymbol{ : }\boldsymbol{p}\boldsymbol{p}$ averaged over all the swimmers with respect to time. Both pushers and pullers show a weak parallel alignment at the lower volume fraction. At the higher volume fraction, the alignment is greater for pushers, consistent with the mutual reinforcement of velocity and orientation fluctuations that underlies the linear stability scenario (Subramanian & Koch Reference Subramanian and Koch2009; Subramanian & Nott Reference Subramanian and Nott2012).

Figure 8.

The polar pair-correlation function as a function of distance from the test swimmer for ${\it\nu}=0.1$ , ${\it\nu}=0.5$ and ${\it\nu}=1.0$ ; $L_{box}=10L$ ( ${\it\nu}=1.0$ corresponds to 8000 swimmers).

Figure 9.

$\langle \boldsymbol{E}\boldsymbol{ : }\boldsymbol{p}\boldsymbol{p}\rangle$ , a measure of the swimmer alignment with the local extensional axis, as a function of time; $L_{box}=10L$ ( ${\it\nu}=0.5$ corresponds to 4000 swimmers).

3.3. Tracer transport

In this section we study the transport of passive non-Brownian tracer particles convected by the swimmer disturbance fields, motivated by one of the earliest experiments on bacterial suspensions by Wu & Libchaber (Reference Wu and Libchaber2000) which involved the measurement of tracer diffusivities to characterize the presence of large-scale correlated motion in a freely standing soap film with swimming bacteria. Figure 10 shows the time evolution of tracer mean-square displacements. An initial ballistic regime is followed by long-time diffusive behaviour. The mean-square displacements in pusher and puller suspensions are comparable at low volume fractions ( ${\it\nu}=0.1$ ). However, there is a remarkable increase for pushers at higher volume fractions with the displacements exceeding those of pullers by around two orders of magnitude at ${\it\nu}=0.6$ . This is consistent with dramatic differences between the tracer trajectories in pusher and puller suspensions (figure 11). For pullers, the tracer particles execute a random walk with no significant persistence (figure 11 a,b). In striking contrast, the tracer trajectories in pusher suspensions show significant correlated motions, filling up the entire field of view (figure 11 c,d). The tracer diffusivities are plotted as a function of the swimmer volume fraction in figure 12, wherein the pusher curves show a clear bifurcation beyond a critical volume fraction with diffusivities that are much larger than those for pullers. Importantly, the simulation box size has a marked effect, with the bifurcation point shifting to lower volume fractions as the box size is increased. This is due to tracer diffusivities for pushers, at a given volume fraction (beyond the bifurcation point), increasing with increasing box size. The bifurcation point is in good agreement with that predicted by the linear theory (see (1.5)). The puller curves seem fairly box-size-independent, pointing to the absence of any instability and collective motion. While this is the case for the straight swimmers examined here, a non-trivial box-size-dependence will be seen to arise later in the presence of an intrinsic decorrelation mechanism (see §§ 4.2.2 and 4.3).

Among all statistical measures considered here, the tracer diffusivity is the most sensitive measure of the onset of collective motion. This can be seen by comparing figure 12 with earlier plots of the mean fluid kinetic energy (figure 5 c) and fluid velocity correlation times (figure 6 c), and is consistent with earlier simulations by Saintillan & Shelley (Reference Saintillan and Shelley2012); the diffusivity values and fluid velocity fluctuations in our simulations are consistently greater than those in Saintillan & Shelley (Reference Saintillan and Shelley2012) owing to our neglect of induced stresses. The tracer diffusivity scales as $D_{tracer}\sim \langle u^{2}\rangle {\it\tau}_{corr}$ . Both the fluid velocity magnitude ( $\langle u^{2}\rangle ^{1/2}$ ) and correlation time ( ${\it\tau}_{corr}$ ) show an increasing trend with volume fraction after the onset of collective motion, and the combined effect is to cause an even larger increase in the trace diffusivity for pushers. Underhill et al. (Reference Underhill, Hernandez-Ortiz and Graham2008) find similar results in their simulations with the tracer diffusivity in pusher suspensions exhibiting a pronounced box-size-dependence; in their case, however, the box-size-dependence is primarily due to a similar dependence of the correlation time. In what follows, we will use the tracer diffusivity to identify the box-size-independent threshold for instability in a suspension of pushers that undergo a run-and-tumble motion. As will be seen, the pusher–puller bifurcation in this case is much more subtle, and the choice of a sensitive statistical measure is indeed crucial.

To summarize, based on the fluid, swimmer and tracer statistics for straight-swimmer suspensions, our kinematic simulations appear to faithfully capture most of the salient features of the collective motion for pusher suspensions, as reported in earlier simulation efforts (Saintillan & Shelley Reference Saintillan and Shelley2007, Reference Saintillan and Shelley2012; Underhill et al. Reference Underhill, Hernandez-Ortiz and Graham2008), and this remains true for both the stable and unstable regimes.

Figure 10.

Mean-square displacements for the tracer particles plotted with respect to time for (a) pushers and (b) pullers; $L_{box}=10L$ .

Figure 11.

Simulation snapshots showing the (straight) swimmer positions (blue circles), orientations (red arrows) and the fluid disturbance velocity vectors (black arrows) on (a,c) at $t=100$ and tracer trajectories, integrated over the period $t=0{-}100$ , on (b,d). (a,b) Puller suspension and (c,d) pusher suspension; $L_{box}=10L$ with ${\it\nu}=0.5$ (4000 swimmers). The figures only display swimmers and tracers in a slab of size $20\times 20\times 3$ in the $x_{1}$ , $x_{2}$ and $x_{3}$ directions, respectively, centred around the midplane of the simulation box (along $x_{3}$ ). The fluid velocity vectors are projected onto this midplane.

Figure 12.

Tracer diffusivities as a function of ${\it\nu}$ , for pushers and pullers, for $L_{box}=5L$ , $10L$ and $15L$ . The vertical lines indicate the box-size-dependent theoretical predictions given by (1.5).

4.1. Theory for tracer diffusion in a suspension of tumblers

Prior to discussing the simulation results, we discuss existing theoretical work related to tracer transport in a dilute stable swimmer suspension. Lin, Thiffeault & Childress (Reference Lin, Thiffeault and Childress2011) considered tracer transport in a suspension of spherical squirmers that execute a run-and-tumble motion with constant run lengths ( $U{\it\tau}$ constant) using a kinetic theory approach first employed by Thiffeault & Childress (Reference Thiffeault and Childress2010). In the limit $U{\it\tau}\gg L$ , the dominant contribution to the tracer diffusivity was found to come from the $O(1/r^{2})$ stresslet velocity field, and scaled as $nUL^{4}$ , $n$ being the swimmer number density, and was therefore independent of ${\it\tau}$ . Interestingly, Pushkin & Yeomans (Reference Pushkin and Yeomans2013) showed that the tracer diffusivity in a suspension of swimming (point) force dipoles, that undergo a similar run-and-tumble motion, is, in fact, independent of ${\it\tau}$ for any $U{\it\tau}/L$ , on account of the self-similar nature of the $O(1/r^{2})$ velocity field associated with a point dipole, a conclusion also reached by Kasyap, Koch & Wu (Reference Kasyap, Koch and Wu2014). Miño et al. (Reference Miño, Dunstan, Rousselet, Clément and Soto2013) considered tracer diffusion arising in a suspension of straight swimmers, and surprisingly, the leading-order scaling for the diffusivity differed from the large- ${\it\tau}$ scaling of Lin et al. (Reference Lin, Thiffeault and Childress2011). The difference arises from the vanishing of the tracer displacement in the Eulerian approximation (fixed tracer location), due to the symmetry of the swimmer-induced displacements; indeed, even the classical calculation of the drift induced by a translating passive particle (Darwin Reference Darwin1953), with a fore–aft symmetric potential flow field, crucially relies on the Lagrangian asymmetry resulting from the convection of the tracer particle. As shown by Kasyap et al. (Reference Kasyap, Koch and Wu2014), and explained below, the discrepancy between Lin et al. (Reference Lin, Thiffeault and Childress2011) and Miño et al. (Reference Miño, Dunstan, Rousselet, Clément and Soto2013) points to the deterministic swimmer not being the same as a tumbler in the limit of a vanishingly small tumble frequency. In other words, the limit of vanishingly small stochasticity is a singular one as far as the tracer diffusivity scaling is concerned. More recently, Morozov & Marenduzzo (Reference Morozov and Marenduzzo2014), have again considered dipolar swimmers characterized by a constant path length, and find a tracer diffusivity $D=AnU{\it\sigma}^{4}$ similar to Lin et al. (Reference Lin, Thiffeault and Childress2011). Based on numerical simulations, the authors show that the prefactor $A$ goes to zero for large $U{\it\tau}/L$ . It was shown by Pushkin & Yeomans (Reference Pushkin and Yeomans2014) that this is due to the finite system size in the simulations, but a clear explanation for this non-trivial effect was not given, and neither was an expression for the tracer diffusivity in a periodic-box derived by rigorously accounting for hydrodynamic interactions with image swimmers.

The tracer diffusivity calculation most relevant to a suspension of slender-bodied swimmers that undergo a run-and-tumble dynamics, with the run-lengths governed by a Poisson distribution, is that of Kasyap et al. (Reference Kasyap, Koch and Wu2014). The diluteness assumption allowed the authors to neglect swimmer–swimmer and swimmer–tracer correlations, and each tracer particle is, at any given time, convected by the velocity field due to a single swimmer. As pointed out earlier in the context of the simulations by Graham and co-workers, a slender swimmer model, apart from being a more faithful representation of bacteria such as E. coli or B. subtilis, has the additional important advantage of avoiding a regularization length scale. Kasyap et al. (Reference Kasyap, Koch and Wu2014) found the tracer diffusivity to be of the form $D_{h}=nUL^{4}\widetilde{D}_{h}$ , where $\widetilde{D}_{h}$ is a monotonically increasing function of the non-dimensional mean-run-time $U{\it\tau}/L$ with $\widetilde{D}_{h}\rightarrow 1/(768{\rm\pi}M^{2})$ for $U{\it\tau}/L\rightarrow \infty$ ( $M$ here being the mobility of the swimmer head). For nearly straight swimmers ( $U{\it\tau}/L\gg 1$ ), the contribution to the tracer diffusivity comes from length scales much larger than the swimmer length, or equivalently, from asymptotically small wavenumbers in the tracer diffusion integral (see (4.1) below) where tumbling remains important. This peculiarity implies, as mentioned above, that the limit of a vanishing rate of relaxation (either tumbling or rotary diffusion) is a singular one. This also raises an interesting question regarding the tracer diffusivity in a periodic-box, relevant to our simulations, where the finite contribution due to distant tumbles will not occur once the swimmer trajectory exceeds the box dimensions.

Starting from equation (21) in Kasyap et al. (Reference Kasyap, Koch and Wu2014) for the tracer diffusivity in an unbounded swimmer suspension, given by the following integral of the Fourier transformed swimmer velocity field:

the tracer diffusivity in a spatially periodic setting is given by:

where the integral over wavenumber space is simply replaced by a summation over discrete wavevectors normalized by the volume of the simulation box. Here, $\hat{\boldsymbol{U}}(\boldsymbol{k})=\int \text{d}\boldsymbol{x}\,\boldsymbol{u}^{\prime }(\boldsymbol{x})\text{e}^{-2{\rm\pi}\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}}$ , $\boldsymbol{u}^{\prime }(\boldsymbol{x})$ being the spatially periodic velocity field due to an infinite cubic array of swimmers. The tracer above is considered to be Brownian with a diffusivity $D$ . The expression for $\hat{\boldsymbol{U}}(\boldsymbol{k})$ is given by:

where $F(\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{p})$ is related to the swimmer (line) force distribution, and is given by Kasyap et al. (Reference Kasyap, Koch and Wu2014):

In (4.3) and (4.4), the wavevector is scaled by the factor $(2{\rm\pi}L_{box})^{-1}$ . Substituting (4.3) and (4.4) in (4.2), and rescaling the terms in the denominator, we obtain:

where we have rescaled ${\it\tau}$ by $L/(2U)$ , and $Pe=UL/D$ . For our simulations, ${\it\alpha}_{1}=1/2$ , and one obtains:

Substituting the above expression in (4.5), and considering $Pe\rightarrow \infty$ (a non-Brownian tracer, as is the case in the simulations), we find:

where we have denoted the ratio $L/L_{box}=L^{\ast }$ . Interchanging the integral and summation, and using a polar coordinate system aligned with $\hat{\boldsymbol{k}}$ , we have $\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{p}=k\cos {\it\theta}$ and $\text{d}\boldsymbol{p}=\sin {\it\theta}\,\text{d}{\it\theta}\,\text{d}{\it\phi}$ , and (4.7) can be written as:

with $m=\cos {\it\theta}$ , which is the final expression for the non-dimensional tracer diffusivity in a spatially periodic (dilute) swimmer suspension.

For arbitrary ${\it\tau}$ , the integral and the $\boldsymbol{k}$ -summation in (4.8) are performed numerically. The analytical forms for ${\it\tau}\rightarrow 0$ and ${\it\tau}\rightarrow \infty$ follow immediately from (4.8). Considering first the ${\it\tau}\rightarrow 0$ limit, and expanding the denominator $(\text{i}kmL^{\ast }{\it\tau}/2+1)$ for small ${\it\tau}$ , gives:

where we have used the fact that the integrand is an even function of $m$ . From (4.9), the dimensional tracer diffusivity for small ${\it\tau}$ scales as $O(nL^{3})U^{2}{\it\tau}$ . An analogous small- ${\it\tau}$ asymptote may be obtained for an unbounded domain from (4.1), and corresponds to $D_{h}=1/(576{\rm\pi}M^{2})(nL^{3})U^{2}{\it\tau}$ . In either case, the tracer is convected, with a probability of $O(nL^{3})$ by a bacterium velocity field of $O(U)$ that decorrelates after a (short) time of $O({\it\tau})$ . With decreasing $L^{\ast }$ , the periodic-box asymptote approaches the unbounded estimate, implying that periodicity is a regular effect for ${\it\tau}\rightarrow 0$ .

For ${\it\tau}\rightarrow \infty$ , we rewrite (4.8) as:

Expanding the denominator $[1+2/(\text{i}kmL^{\ast }{\it\tau})]$ for large ${\it\tau}$ , we obtain:

The integral in the first term is odd in $m$ and hence is zero. The second term, therefore, gives the leading-order contribution to the diffusivity as:

This leads to the counter-intuitive conclusion that the periodic-box tracer diffusivity vanishes in the limit ${\it\tau}\rightarrow \infty$ (nearly straight swimmers). This prediction (4.12) is in contrast to an unbounded suspension where a finite diffusivity emerges in this limit (Lin et al. Reference Lin, Thiffeault and Childress2011; Kasyap et al. Reference Kasyap, Koch and Wu2014). The mathematical reason for this can be deduced from (4.1) (see Kasyap et al. Reference Kasyap, Koch and Wu2014) where, in the limit of large ${\it\tau}$ , the integral requires a Cauchy principal-value interpretation, the dominant contribution coming from wavenumbers of $O(1/U{\it\tau})$ or smaller. Physically, this corresponds to a sufficiently large number of distant tumbles continuing to break the symmetry of the (deterministic) swimmer–tracer interaction at leading order. However, in a periodic-box, we see from the denominator of (4.8) that the smallest wavenumber is $O(L^{\ast })$ . For $1/{\it\tau}\ll L^{\ast }$ ( $(U{\it\tau})^{-1}\ll L_{box}^{-1}$ in dimensional terms), the dominant contribution in the unbounded case is not captured by the discrete sum over $\boldsymbol{k}$ . Thus, rather interestingly, the artefact of periodicity, inherent in any numerical simulation, renders the limit of vanishingly small stochasticity a regular one, in the sense of the large- ${\it\tau}$ limit being consistent with the zero tracer diffusivity obtained in the Eulerian fixed tracer approximation applied to a suspension of straight swimmers (Miño et al. Reference Miño, Dunstan, Rousselet, Clément and Soto2013; Pushkin, Shum & Yeomans Reference Pushkin, Shum and Yeomans2013). It may be shown that the leading-order tracer diffusivity vanishes even for ${\it\alpha}_{1}\neq 1/2$ , and the above comments therefore remain valid for swimmers with asymmetric force distributions (as must be the case, since all such swimmers appear as point force dipoles in the far-field). The dimensionless (non-Brownian) tracer diffusivity is plotted in figure 13 for a number of box sizes along with the result for an unbounded suspension. The latter increases monotonically with ${\it\tau}$ with an $O(nUL^{4})$ large- ${\it\tau}$ asymptote. The periodic-box diffusivity, however, is non-monotonic with a maximum of $O(nUL^{4})$ at a finite ${\it\tau}$ of $O(L/U)$ ; note that this maximum should spread out into a broad plateau in the range $L/U\ll {\it\tau}\ll L_{box}/U$ with increasing $L_{box}$ . The small and large ${\it\tau}$ asymptotes given by (4.9) and (4.12) are also shown.

Figure 13.

Theoretical predictions for the tracer diffusivity in an unbounded suspension and in a periodic array of swimmers ( $L_{box}=10L,15L$ and $25L$ ). In the former case, the diffusivity monotonically increases to a ${\it\tau}$ -independent limit (solid black curve). The periodic-box results show a non-monotonic trend, with the small ${\it\tau}$ and large ${\it\tau}$ asymptotes being $O({\it\tau})$ and $O(1/{\it\tau})$ , respectively.

The non-trivial role of periodicity above is evidently relevant to our simulations. The $O(1/{\it\tau})$ scaling of the periodic-box tracer diffusivity, leading to the aforementioned non-monotonicity, results from two approximations in the above analysis – the approximation of a fixed tracer location and the neglect of hydrodynamic interactions between swimmers. As already pointed out, accounting for the tracer displacement during a swimmer–tracer interaction will lead to a non-zero tracer diffusivity as ${\it\tau}\rightarrow \infty$ , even in a spatially periodic setting, and such displacements occur naturally in the simulations discussed below; this effect is, however, small in the limit of slender swimmers owing to the weak disturbance field that convects the tracer. More importantly, in the presence of inter-swimmer hydrodynamic interactions, the decrease in the periodic-box tracer diffusivity, for large ${\it\tau}$ , will be cutoff once the rate of swimmer orientation decorrelation due to tumbles becomes comparable to that arising from hydrodynamic interactions. In other words, for large enough ${\it\tau}$ , the tracer diffusivity must approach a plateau controlled by the hydrodynamically induced orientation decorrelation arising from (weak) swimmer interactions. In the dilute regime, and for large aspect ratios ( ${\it\gamma}\gg 1$ ), this decorrelation arises from the cumulative effects of small $O(\ln {\it\gamma})^{-1}$ orientation changes associated with successive pair interactions, and may be characterized by a rotary diffusivity (Rahnama, Koch & Shaqfeh Reference Rahnama, Koch and Shaqfeh1995). This hydrodynamic rotary diffusivity $(D_{rh})$ for slender-bodied swimmers was calculated by Subramanian & Koch (Reference Subramanian and Koch2009), who found it to be $O[nUL^{2}/(\ln {\it\gamma})^{2}]$ for an unbounded domain. Now, from (4.12), the large- ${\it\tau}$ asymptote of the periodic-box tracer diffusivity is $O(nUL^{4}{\it\tau}^{-1})$ , or in dimensional terms $O(nL^{5}/{\it\tau})$ , and this decrease with ${\it\tau}$ will continue until ${\it\tau}\sim O(D_{rh}^{-1})$ when the tracer diffusivity curve should turn around, eventually saturating at a value of $O(nL^{5}D_{rh})\sim ((nL^{3})/(\ln {\it\gamma})^{2})O(nUL^{4})$ for ${\it\tau}D_{rh}\gg 1$ (on use of the above infinite-domain estimate for $D_{rh}$ valid for $L_{box}>O(UD_{rh}^{-1})$ ). Thus, the non-monotonic box-size-dependence predicted above arises essentially due to the hydrodynamic orientation decorrelation being much weaker than the intrinsic decorrelation (due to tumbling or rotary diffusion) for the slender swimmers considered here, corresponding to $nL^{3}\ll O(\ln {\it\gamma})^{2}$ for large enough box sizes. In the other extreme, where the hydrodynamic interactions are strong and induce a decorrelation on a time scale much shorter than the $O(L_{box}/U)$ scale corresponding to periodicity, the rate of orientation decorrelation is given by $[1/{\it\tau}+2D_{rh}]^{-1}$ , and an estimate of the tracer diffusivity is obtained as $(nL^{3})U^{2}[1/{\it\tau}+2D_{rh}]^{-1}$ with ${\it\tau},D_{rh}^{-1}\ll L_{box}/U$ . This expression implies a monotonic increase of the tracer diffusivity with ${\it\tau}$ with a large- ${\it\tau}$ plateau of $O(nL^{3}U^{2}/D_{rh})$ . The decorrelation due to strong hydrodynamic interactions is not a diffusive process in orientation space, and $D_{rh}^{-1}$ here must merely be regarded as the inverse of a hydrodynamic decorrelation time ( ${\it\tau}_{hyd}$ ). Thus, for strong hydrodynamic interactions, the periodic-box tracer diffusivity must rise monotonically to its large- ${\it\tau}$ asymptote. It is clear then that the behaviour of the tracer diffusivity, as a function of ${\it\tau}$ , in a spatially periodic swimmer suspension, depends on the strength of the inter-swimmer hydrodynamic interactions; the interactions may be termed weak or strong depending on the magnitude of the hydrodynamic orientation decorrelation time in relation to the time scale associated with periodicity. The strength of the inter-swimmer hydrodynamic interactions is a sensitive function of the swimmer model with earlier simulations having considered both weakly (Saintillan & Shelley Reference Saintillan and Shelley2007, Reference Saintillan and Shelley2012) and strongly (Underhill et al. Reference Underhill, Hernandez-Ortiz and Graham2008) interacting swimmers. Towards the end of this section, we show, in our simulations, that there is indeed a transition of the tracer diffusivity curves from a non-monotonic to a monotonic variation with ${\it\tau}$ with decreasing hydrodynamic decorrelation time.

Figure 14.

Tracer diffusivities for both straight swimmers and tumblers with ${\it\tau}=1$ and ${\it\tau}=3$ , plotted as a function of volume fraction; $L_{box}=10L$ .

For slender swimmers, the above discussion shows that there exists an upper limit of ${\it\tau}\sim O(L_{box}/U)$ above which box size effects lead to a non-monotonic dependence of $D_{h}$ on $U{\it\tau}/L$ . However, the threshold criterion discussed earlier implies that ${\it\tau}$ must be greater than $O(nUL^{2})^{-1}$ for instability in a suspension of run-and-tumble swimmers. Thus, for the ${\it\tau}$ -varying simulations, considered in §§ 4.2.2 and 4.2.3, one only has a finite ${\it\tau}$ -window, given by $(nUL^{2})^{-1}<{\it\tau}<L_{box}/U$ , available for examining the transition to instability representative of an unbounded suspension. In particular, if the volume fraction chosen is too small (with the intent of rendering swimmer–swimmer correlations vanishingly small in the stable dilute regime), the lower bound would begin to approach $O(L_{box}/U)$ , and effects related to periodicity would become dominant.

4.2. Simulations of suspensions of tumblers

Since the parameter $nUL^{2}{\it\tau}$ is predicted to govern the transition to collective motion in a suspension of run-and-tumble swimmers (Subramanian & Koch Reference Subramanian and Koch2009), one may cross-over to the unstable regime by increasing $nL^{3}$ at a fixed ${\it\tau}$ or vice versa. We discuss the volume-fraction-varying and ${\it\tau}$ -varying simulations in sequence.

4.3. Simulations of suspensions of rotary diffusers

Our approach to studying the effects of rotary diffusion parallels that of § 4.2.2 for tumbling. We use the tracer diffusivity as the measure to estimate the stability threshold in simulations that vary the rotary diffusivity at a fixed volume fraction ( ${\it\nu}=0.05$ ). Again, we contrast the tracer diffusivities in suspensions of pushers and pullers via sets of ‘interactions-off’ and ‘interactions-on’ simulations, respectively. The results for both sets are presented in figure 24 in a format similar to that of figure 18. The results correspond to box sizes of $10$ -, $15$ - and $30$ -swimmer-lengths, with each row corresponding to a given box size, and the left and right columns corresponding to interactions being switched off and on, respectively. Focusing first on the ‘interactions-off’ results, for the smaller $D_{r}$ values, we see a trend similar to that seen in the case of tumblers in that the pusher and puller diffusivities are nearly coincident and decrease with decreasing $D_{r}$ . Although an analytical calculation is not available in this case, one expects the dependence of the tracer diffusivity on $D_{r}^{-1}$ , in an unbounded swimmer suspension, to parallel that found by Kasyap et al. (Reference Kasyap, Koch and Wu2014) for tumblers. That is to say, the tracer diffusivity is expected to be a monotonically increasing function of $D_{r}^{-1}$ , asymptoting to a value of $O(nUL^{4})$ in the deterministic limit. Although, in light of the singular nature of this limit, the numerical pre-factor will likely differ from that found for tumbling (the tracer diffusivity, in the limit of vanishing stochasticity, is expected to be of the general form $f({\it\tau}D_{r})nUL^{4}$ where $f$ takes different limiting values, both of order unity, for ${\it\tau}D_{r}\rightarrow 0$ and ${\it\tau}D_{r}\rightarrow \infty$ ). Thus, the decrease in the tracer diffusivity for large $D_{r}^{-1}$ , and the resulting non-monotonicity observed in the ‘interactions-off’ simulations, are again due to the effects of spatial periodicity. In contrast, for large $D_{r}$ ’s, the tracer diffusivity will be limited by diffusive decorrelation of the swimmer trajectories. This is evident from the scaling of the small- ${\it\tau}$ estimate which involves the swimmer translational diffusivity of $(U^{2}{\it\tau})/3$ arising from random tumbles. Since the long-time translational diffusivity arising from rotary diffusion is $U^{2}/(6D_{r})$ , the tracer diffusivity for large $D_{r}$ is obtained directly from (4.9) by replacing ${\it\tau}$ in the diffusivity scale with $1/(2D_{r})$ (Sandoval et al. Reference Sandoval, Navaneeth, Subramanian and Lauga2014). However, from figure 24(a,c,e), it is clear that the ‘interactions-off’ simulations do not conform to the $O(D_{r}^{-1})$ scaling expected. The reasons for this discrepancy are not entirely clear. A possible reason is again the excluded volume around a swimmer leading to short-range steric interactions that appear to modify the tracer transport much more strongly for rotary diffusers than for tumblers. More importantly, the tracer diffusivities for both pushers and pullers are almost coincident at the intermediate $D_{r}$ ’s, for the two larger box sizes (see figure 24 c,e), which is the range where the pusher–puller bifurcation occurs on turning on the swimmer interactions. Indeed, when interactions are switched on, a clear bifurcation can be seen between pushers and pullers (figure 24 b,d,f) despite the artefacts of excluded volume for large $D_{r}$ . Beyond a critical value of $1/D_{r}$ , pushers suspensions show consistently higher tracer diffusivities when compared to pullers. Significantly, the pusher tracer diffusivity, at a fixed volume fraction, shows a clear increase with increasing box size, once again pointing to hydrodynamic instability, rather than local pair interactions, as the underlying cause. In light of the large- $D_{r}$ numerical artefacts, the critical $1/D_{r}$ is identified as the value next to the cross-over point of the pusher and puller curves. For the $10$ and $30$ -swimmer-length boxes, the resulting value, $1/D_{r}\approx 10$ , is also the point at which the pusher diffusivities first exceed the puller ones by an amount greater than ${\it\nu}$ , consistent with (4.14). For the $15$ -swimmer-length box, application of (4.14) gives $1/D_{r}\approx 5$ , but this is marginal since the difference between the diffusivities only just exceeds the threshold at this point. Given the finite resolution along the $D_{r}$ axis, a reasonable estimate of the bifurcation corresponds to $1/D_{r}|_{crit}\approx 10$ . The linear stability for rotary diffusers is a critical $nUL^{2}/D_{r}$ , and based on the above critical $D_{r}$ , with ${\it\nu}=0.05$ , one obtains: