Open Access Article
This Open Access Article is licensed under a
Creative Commons Attribution 3.0 Unported Licence

Particle-resolved lattice Boltzmann simulations of 3-dimensional active turbulence

Dóra Bárdfalvya, Henrik Nordangera, Cesare Nardinib, Alexander Morozovc and Joakim Stenhammar*a
aDivision of Physical Chemistry, Lund University, Box 124, S-221 00 Lund, Sweden. E-mail:
bService de Physique de l'État Condensé, CNRS UMR 3680, CEA-Saclay, 91191 Gif-sur-Yvette, France
cSUPA, School of Physics and Astronomy, The University of Edinburgh, James Clerk Maxwell Building, Peter Guthrie Tait Road, Edinburgh, EH9 3FD, UK

Received 16th April 2019 , Accepted 31st July 2019

First published on 2nd August 2019

Collective behaviour in suspensions of microswimmers is often dominated by the impact of long-ranged hydrodynamic interactions. These phenomena include active turbulence, where suspensions of pusher bacteria at sufficient densities exhibit large-scale, chaotic flows. To study this collective phenomenon, we use large-scale (up to N = 3 × 106) particle-resolved lattice Boltzmann simulations of model microswimmers described by extended stresslets. Such system sizes enable us to obtain quantitative information about both the transition to active turbulence and characteristic features of the turbulent state itself. In the dilute limit, we test analytical predictions for a number of static and dynamic properties against our simulation results. For higher swimmer densities, where swimmer-swimmer interactions become significant, we numerically show that the length- and timescales of the turbulent flows increase steeply near the predicted finite-system transition density.

1 Introduction

An archetypical example of active matter is a suspension of synthetic or biological particles that possess the ability to convert the energy extracted from their surroundings into self-propulsion.1,2 This conversion of chemical energy into mechanical energy at the single-particle level and the resulting violation of detailed balance can result in rich displays of collective motion and dynamical self-assembly at larger lengthscales, such as strongly ordered bird flocks3 and dynamical clustering in suspensions of active colloids.4

One of the most well-studied active systems is a suspension of bacteria or algae that swim by beating or rotating a collection of flagella.5,6 At low densities, such systems show significantly enhanced diffusion of nonmotile particles compared to Brownian diffusion.7–11 At higher concentrations, though still rather dilute, suspensions of rear-actuated (pusher) bacteria exhibit a complex collective behaviour known as active turbulence, whereby the system starts exhibiting large-scale vortices and jets with higher fluid velocities than the velocity of the individual swimmer.12–15 For front-actuated (puller) swimmers such as Chlamydomonas, no such collective behaviour is observed in 3-dimensional suspensions,7,16 although instead a transition to a polar flocking state has been observed in simulations of puller stresslets confined to a 2-dimensional plane.17,18 We also note that, in the case of squirmers, which swim by an imposed slip flow along the spherical19–21 or elongated22,23 swimmer body, such a polar state is found for pullers also in 3 dimensions, while no sign of collective behaviour is found in the corresponding pusher suspensions.19,21 While squirmers is a more appropriate model for ciliated organisms such as Paramecium, these different collective behaviours highlight that the aspect ratio and seemingly subtle differences in the near-field flows can have large impacts on the non-equilibrium steady states.

The transition to active turbulence in pusher suspensions has been described as a hydrodynamic instability induced by the mutual reorientation of swimmers due to the long-ranged stresslet flow fields.1,24–30 Much of the theoretical understanding of active turbulence is based on continuum theories that describe the active suspension using effective equations of motion for the order parameter fields. These are typically derived from either a kinetic theory of a set of stresslet swimmers,1,16,24,28,29 or through the modification of the equations describing nematic liquid crystals through additional active stresses, resulting in so called active nematics models.31–33 For the former class of models, the microscopic parameters describing the transition to active turbulence can be determined through stability analysis of the linearised equations. For an unbounded suspension, the most unstable mode is the k = 0 one, and the ensuing analysis leads to the following prediction for the critical pusher number density nc required for collective motion:24,26,28,34,35

nc = 5λ/κ. (1)
Here, λ is the tumbling frequency by which individual swimmers randomise their swimming direction and κ is the stresslet magnitude, defined below. In order to probe the properties of the turbulent state itself, the linear theory is no longer accurate, and one instead needs to numerically integrate the equations of motion.31,34,36,37 This approach gives access to the full nonlinear behaviour of the turbulent state, although still resting on the approximations underlying the continuum equations. Particle-based simulation studies of active turbulence are more limited in number, partially due to the computational challenges of simulating large collections of hydrodynamically interacting particles over extended length- and timescales. These difficulties have so far prohibited a quantitative verification of the analytical predictions, such as that of the transition density. Hernandez-Ortiz et al.38,39 and Lushi and Peskin40 employed a model that describes each pusher swimmer as two connected spheres with a pair of embedded point forces, and observed a transition to a collectively flowing state resembling active turbulence. They characterised the coherent flows by calculating the properties of the fluid flows and the enhanced diffusion of tracer particles. Subsequently, Saintillan and Shelley41 and Krishnamurthy and Subramanian42 observed similar collective properties in suspensions containing up to N = 3 × 104 slender rod-like particles with pusher flow-fields in 3 dimensions.

In this paper, we significantly extend the above studies to system sizes reaching N > 106 microswimmers, which enables us to quantitatively study the properties of active turbulence with high numerical accuracy over extended length- and timescales. The simulations are based on an implementation of the lattice Boltzmann (LB) equation43 that allows us to describe each swimmer as a pair of point forces acting on the surrounding fluid, thus ignoring the effect of near-field hydrodynamics and excluded volume interactions, which can be justified by the relatively low number densities needed to reach the turbulent state. We show that the length- and timescales of the chaotic flows increase sharply at the transition to turbulence, before decreasing towards a plateau value in the turbulent state. The swimmer density where this increase is observed is in agreement with the theoretical prediction from kinetic theory, as long as we take into account the finite-size effects from using a finite box with periodic boundary conditions. We also show that studying these phenomena require very large system sizes: our results suggest that linear box dimensions of at least ∼100 times the swimmer length are required to eliminate finite-size effects. These large systems furthermore allow us to study the statistics of swimmer–swimmer correlations with unprecedented accuracy: the spatial two-body correlation functions show a transition to a state of strong orientational order, induced purely by far-field hydrodynamic interactions, but with no signs of significant density inhomogeneities. Our results provide a thorough characterisation of the hydrodynamically induced collective motion in pusher suspensions which should pave the way both for further experimental efforts and for testing analytical descriptions of active turbulence.

2 Model and methods

We consider a suspension of N swimmers represented by extended force dipoles (stresslets) moving in a three-dimensional fluid with periodic boundaries. The body and flagella exert two equal and opposite forces ±Fp, where p is the orientation of the swimmer, separated by a length l on the fluid as shown in Fig. 1. The swimmer is characterised by a dipole strength κ = ±Fl/μ where μ is the dynamic viscosity of the fluid, with κ > 0 representing pushers and κ < 0 pullers.
image file: c9sm00774a-f1.tif
Fig. 1 Schematic image of the pusher model, where F is the force, l is the swimmer length, vs the swimming speed, a the effective body radius and p is the orientation of the swimmer.

The position r and orientation p of each swimmer evolves according to the following equations of motion:24,44

= vsp + U(r), (2)
image file: c9sm00774a-t1.tif(3)
Here, U(r) is the fluid velocity evaluated at the body position, vs is the constant swimming velocity of the individual swimmer, and image file: c9sm00774a-t2.tif is the unit tensor. Eqn (3) is the (discretised version of) Jeffery's equation for infinite aspect ratio (β = 1); we checked our results also using a finite aspect ratio (β < 1), which yielded only small shifts for realistic values of β. In addition to their reorientation as a result of the hydrodynamic interactions described by eqn (3), swimmers also undergo Poisson-distributed random reorientations with an average frequency λ. (Note that this mechanism cannot be cast as a continuous-time differential equation as that in eqn (3).28) This run-and-tumble motion results in a random walk with a persistence length vs/λ.

Note that, in our model, the symmetry breaking in the single swimmer dynamics is only created by the self-propulsion term in the equation of motion (2): the single-swimmer flow-field is fully fore–aft symmetric, as there is no surface describing the swimmer body. However, an effective body radius a for the swimmer geometry in Fig. 1 can be calculated through the following relation between F, vs, and a:24

image file: c9sm00774a-t3.tif(4)
Here, the second term in brackets constitutes the leading-order correction to the Stokes–Einstein relation for the spherical swimmer body due to the disturbance flow created by the flagellar Stokeslet.

To simulate the hydrodynamic interaction between microswimmers, we used a D3Q15 BGK lattice Boltzmann (LB) method based on the point-force LB implementation of Nash et al.43,45 The key ingredient in this method is an algorithm for interpolating forces and velocities between the fluid and the off-lattice swimmers that employs the regularised version of the δ function due to Peskin.43,46 Since this function has a compact support of 2 lattice units, it effectively gives non-singular flow-fields that are distributed over distances comparable to the LB lattice spacing. Importantly, the method enables large-scale simulations of up to N ∼ 106 microswimmers at biologically relevant densities. In the simulations, cubic box sizes L3 ranging from (10)3 to (210)3 lattice sites were employed. In terms of LB units (where ΔL = Δt = 1), we used the parameters vs = 10−3, F = 1.57 × 10−3, l = 1, λ = 2 × 10−4, and μ = 1/6, with the latter value corresponding to the fluid relaxing to local equilibrium on each timestep. Each simulation was run for 2 × 105 timesteps, apart from in the transition region where longer simulations were necessary due to the slow dynamics. A typical simulation with L = 100 took ∼48 hours on a single CPU, while the largest simulations (L = 210) took ∼1 week in parallel on 5 CPUs.

All results will be presented in terms of the swimmer length l and the time scale l/vs. The relevant dimensionless numbers of the system are (i) the single-swimmer Reynolds number Resρfvsl/μ, where ρf is the fluid density (set to unity in our LB simulations), and (ii) the nondimensional stresslet strength κnκ/(l2vs) = F/(μlvs). The above parameter values yield Res = 6 × 10−3, which is well below the Stokes-flow limit,47 and κn ≈ 9.4. The latter value can be compared with the corresponding value for E. coli for which vs ≈ 22 μm s−1, F = 0.42 pN, and l = 1.9 μm,48 yielding κn ≈ 11.2. Furthermore, an approximate volume fraction based on the spherical swimmer body can be calculated with the aid of eqn (4) as ϕ = (4π/3)a3n, where n = N/L3 is the swimmer number density; with our parameters, eqn (4) gives a ≈ 0.3. Using this estimate, the observed critical volume fraction needed for collective motion is ϕc ≈ 0.02, which is clearly within the range of validity of our far-field hydrodynamic model. While the highest densities considered in this paper (ϕ ≈ 0.055) are large enough that near-field effects and other specific interactions would start to become significant, we would like to highlight that these strongly turbulent flows can also be achieved at small densities if the dipolar strength κ is large enough, in accordance with eqn (1). A relatively high-density suspension of weak dipoles can thus be viewed as a proxy for the turbulent properties of a dilute suspension of strong dipoles.

3 Analytical expressions for noninteracting swimmers

In this section, we will derive analytical results for the statistical and dynamical properties of a collection of noninteracting microswimmers, where we set the terms containing U in eqn (2) and (3) to zero. The swimmers still exert pairs of equal and opposite forces on the fluid, and generate long-ranged flow fields. Note that, in the non-interacting case, pushers and pullers are equivalent. In this limit, many average properties of the system can be obtained analytically with relative ease, since the swimmers are then statistically independent. In particular, we calculate the spatial and temporal correlation functions of the fluid, the fluid velocity variance, and the associated Fourier-space energy spectrum. These predictions will then be compared in the following sections to the corresponding LB simulations of non-interacting swimmers. The corresponding description taking into account swimmer–swimmer interactions is still achievable24 but considerably more complicated; we therefore postpone the full analytical treatment of interacting swimmers to a forthcoming separate study.

Following Cortez et al.,49 we start from the expression for the regularised flow field from a point force with magnitude F:

image file: c9sm00774a-t4.tif(5)
Here, ε is a factor describing the distance over which the regularisation acts and r = |r|; in the limit ε → 0 this expression reduces to the usual Stokeslet. Its Fourier transform is given by
image file: c9sm00774a-t5.tif(6)
where K2 is the modified Bessel function of the second kind, and [k with combining circumflex]i = ki/k, where k = |k|. The velocity field of an extended dipole is constructed from this expression by placing two point forces of equal magnitude F and opposite orientations at r0 and r0 + lp, where l is the dipolar length and p its orientation. The corresponding real-space velocity field of a single swimmer is thus
image file: c9sm00774a-t6.tif(7)
The total velocity field at a position r created by a suspension of N non-interacting swimmers with instantaneous positions ri and orientations pi, i = 1…N is, then,
image file: c9sm00774a-t7.tif(8)

3.1 Temporal correlations

First we consider the velocity–velocity autocorrelation function c(t) = 〈U(0)·U(t)〉, which is given by the following ensemble average
image file: c9sm00774a-t8.tif(9)
Here, we used the fact that the time-dependence of the velocity field only arises through changes of the swimmer positions; the factor eλt accounts for independent decorrelation events due to tumbling, and vs is the swimming speed. Since the swimmers are statistically independent, only ‘self-correlations’ contribute to the average, yielding
image file: c9sm00774a-t9.tif(10)
image file: c9sm00774a-t10.tif(11)
Here, as before, n = N/V, is the swimmer number density. Explicit integration gives
image file: c9sm00774a-t11.tif(12)
image file: c9sm00774a-t12.tif(13)
with [Doublestruck E](x) and [Doublestruck K](x) being the complete elliptic integrals of the first and second kind, respectively. We have furthermore introduced the dimensionless time tn = tvs/l, tumbling rate λn = λl/vs, and regularisation parameter Δ = l/ε, and we have normalised c(tn) so that c(0) = 1.

3.2 Spatial correlations

Similar to the temporal correlation function, the spatial velocity–velocity correlation function c(R) = 〈U(0U(R)〉 is given by the following ensemble average
image file: c9sm00774a-t13.tif(14)
Again, keeping only the ‘self-correlation’ contributions, we obtain
image file: c9sm00774a-t14.tif(15)
Direct evaluation of this integral yields
image file: c9sm00774a-t15.tif(16)
image file: c9sm00774a-t16.tif(17)
image file: c9sm00774a-t17.tif(18)
image file: c9sm00774a-t18.tif(19)
Here, Rn = R/l, and we have again normalised c(Rn) such that c(0) = 1. For large Rn, c(Rn) ∼ Rn−1, in agreement with the results of Zaid et al.50

3.3 Velocity variance

The velocity variance 〈U2〉 can be obtained from eqn (15) by setting R = 0. The result depends on c0 and reads
image file: c9sm00774a-t19.tif(20)
By combining the large- and small-Δ asymptotics of this expression, we obtain the following uniform approximation, which interpolates well between the two regimes and is significantly easier to use than the full expression:
image file: c9sm00774a-t20.tif(21)

3.4 Energy spectrum

Eqn (15) can be interpreted as a Fourier transform of the velocity–velocity spatial correlation function, and we therefore can identify (after the substitution k → −k due to our definition of the Fourier transform)
image file: c9sm00774a-t21.tif(22)
The energy content associated with the velocity field at a lengthscale k−1 is given by
Ek = 4πk2Û(kÛ(−k)〉k=|k|, (23)
image file: c9sm00774a-t22.tif(24)
In the double limit ε → 0 and l → 0, this expression reduces to
image file: c9sm00774a-t23.tif(25)
which is independent of the wavevector k.

4 Results and discussion

4.1 Fluid statistics

Fig. 3 demonstrates how collective motion develops as a function of the concentration of microswimmers through the root-mean-square (RMS) fluid velocity URMS ≡ 〈U21/2 for suspensions of pushers, pullers and noninteracting swimmers; for the latter simulations the terms containing U in eqn (2) and (3) are set to zero, so that the swimmers do not interact with each other through the fluid. (Note again that, for the noninteracting case, pushers and pullers are equivalent.) For noninteracting swimmers, the RMS fluid velocity is accurately described by (the square root of) eqn (21), yielding an n1/2 dependence over the whole concentration range. Using that expression, we fit the value ε ≈ 1.1 (while fixing all other parameters to the values from the LB simulations), in good accordance with the interpolation length of the Peskin δ function used in the LB simulations.46 At intermediate concentrations for interacting swimmers, there is a deviation from the square root dependence, with the RMS velocity increasing faster than n1/2 for pushers and slower for pullers: this is a signature of the build-up of long-ranged orientational correlations underlying the collective behaviour.24 At concentrations of n > 0.2 the turbulent state is fully developed for pushers, as is clearly visible in the snapshots in Fig. 2 and videos available as ESI. In the following sections, we will focus mainly on pusher suspensions and their transition from disordered swimming to active turbulence.
image file: c9sm00774a-f2.tif
Fig. 2 Snapshots of the fluid velocity field in a 2-dimensional lattice plane of the 3-dimensional simulation box (L = 150) at densities before the transition to turbulence (n = 0.05), close to the transition (n = 0.15), and in the turbulent regime (n = 0.3). Vectors show the velocity field in the xy-plane and the colours indicate its z-component. The length of the velocity vectors have been rescaled for clarity. See also the corresponding videos available as ESI.

image file: c9sm00774a-f3.tif
Fig. 3 Root-mean-square fluid velocity URMS, as a function of the swimmer number density n. The results were obtained using a linear box dimension L = 150 for pushers and L = 100 for noninteracting swimmers and pullers, due to the significant finite-size effects in the turbulent regime of the pusher suspensions. The dashed curve indicates a fit using eqn (21), yielding ε = 1.1.

To characterise the length- and timescales of the chaotic flows, we now turn to their spatial and temporal correlation functions. In Fig. 4, we show both the equal-time velocity correlation functions and the velocity autocorrelation functions, computed for four different densities. While the temporal correlation function c(tn) in suspensions of noninteracting swimmers show reasonable agreement with the theoretical predictions (eqn (12)), the corresponding spatial curves start deviating from the predictions (eqn (16)) for Rn ≈ 10, eventually falling below zero at Rn ≈ 80. We attribute this poor agreement at intermediate and large Rn to the significant effect of periodic boundary conditions on the long-ranged part of the flow fields as well as the presence of higher multipoles in the LB swimmer flow fields. All the correlation functions become increasingly long-ranged when going from n = 0.05 to n = 0.10. Their range, especially that of the autocorrelation function, then increases significantly near the transition around n = 0.15, followed by a slight decrease inside the turbulent regime (n = 0.3). To quantitatively characterise the length- and timescales of the chaotic flow, in accordance with ref. 32, we define the characteristic length ξ as the distance where c(Rn) has decreased to 0.2 (Fig. 4a), whereas the corresponding characteristic time τ is defined as the point when c(tn) has decayed to 0.4 (see Fig. 4b). The difference in the two threshold values is due to the slow decay of c(tn) around the transition density: it reaches 0.2 only after prohibitively long times, resulting in poor statistics. The results, showing ξ and τ as a function of n for a wide range of densities are shown in Fig. 5. In accordance with the results in Fig. 4, there is a very sharp increase of both quantities at n ≈ 0.15, followed by a gradual decrease towards a plateau value. The difference between the transition region and the turbulent region is most pronounced in the τ curve, while the maximum in ξ is somewhat broader and has its peak at n = 0.2 rather than n = 0.15. The clearly non-monotonic curves reported here are different from the results reported previously by Saintillan and Shelley,41 who did not observe any maximum at the transition density for neither the characteristic length or time-scales, while the non-monotonic behaviour of τ(n) was previously observed by Krishnamurthy and Subramanian.42 We attribute this difference to the significantly smaller system sizes used in earlier studies.

image file: c9sm00774a-f4.tif
Fig. 4 Build-up of long-ranged velocity correlations due to collective motion. (a) The spatial velocity correlation function c(Rn), and (b) velocity autocorrelation functions c(tn) for four different densities as indicated. The dashed curve shows the value we used to determine the characteristic length- and timescales ξ and τ. The results were obtained using L = 210. Insets show comparisons between LB simulations of noninteracting swimmers (symbols) and the analytical predictions of eqn (16) and (12), using ε = 1.1 (lines).

image file: c9sm00774a-f5.tif
Fig. 5 The characteristic length ξ (black) and time τ (red) as a function of density n, obtained from systems with L = 150. Error bars indicate the estimated standard deviations obtained from dividing the simulation into four equally sized time intervals. The dashed curves indicate, from left to right, the predicted transition densities nc for an unbounded system (kc = 0, eqn (1)), and for the finite wavenumbers kc = 2π/L, 4π/L, and 8π/L, as discussed in the text.

We furthermore note that the predicted infinite-system critical density nc (eqn (1)) falls somewhat below the observed increases in ξ and τ. In order to qualitatively explain this, we first consider the effect of a finite box size, which shifts the critical wavevector from kc = 0 to kc = 2π/L. In addition, the use of periodic boundary conditions will effectively screen the stresslet flow-fields, yielding them more short-ranged than the r−2 spatial decay in an infinite fluid. While this screening is gradual, it becomes significant already at length scales between L/4 and L/2. Since the long-wavelength instability leading to active turbulence is an effect of the r−2 decay of the flow field, this screening will shift the instability to even shorter lengthscales than the finite-box reasoning alone. Thus, in Fig. 5, we also plot the values of nc corresponding to the critical wavenumbers kc = 0, 2π/L, 4π/L and 8π/L, obtained through a linear stability analysis as detailed elsewhere.35,44 However, a more in-depth analysis of the effect of PBCs together with a more rigorous numerical treatment to numerically localise the transition would be necessary to confirm these qualitative arguments.

In order to further investigate the system size dependence discussed above, in Fig. 6 we show ξ and τ plotted as a function of the linear box size L for the same four densities studied in Fig. 4. In both panels, we observe significant finite-size effects until L ≈ 100, corresponding to N = 3 × 105 for n = 0.3, although these appear to persist to even larger systems in the transition region. These results again highlight the importance of using large-scale simulations when studying collective motion in microscopic models of microswimmers.

image file: c9sm00774a-f6.tif
Fig. 6 Finite-size effects of the characteristic length- and timescales. (a) ξ and (b) τ plotted as functions of the linear system size L. Error bars indicate the estimated standard deviations obtained either from four independent runs with different initial conditions (n = 0.15) or from dividing the full simulation into four equally sized time intervals.

To further analyse the spatial structures of the flow, in Fig. 7 we calculate the Fourier space energy spectrum Ek as defined in eqn (23). For low densities, the spectrum (Fig. 7) is well described by the form predicted for uncorrelated swimmers in eqn (24) with a flat shape at intermediate k, corresponding to a superposition of the r−2 flow fields of uncorrelated swimmers. At the lowest accessible values of k, the spectra decrease slightly compared to the infinite-system prediction, again likely due to the effect of PBCs. At high k, corresponding to the length-scale of individual swimmers, Ek decreases due to the short-range regularisation of the flow field; the slight discrepancy between the data and the analytical prediction in this regime is due to the different forms of regularisation used in the two treatments. In the turbulent regime, most of the kinetic energy is localised at scales much bigger than l, in accordance with previous studies,41,42 even though the collective motion is driven by energy injected at small lengthscales. Another feature of the energy spectrum, which was absent in previous studies due to finite-size effects, is the peak in the spectrum which evolves for high swimmer densities, again indicating a characteristic, finite length-scale of the flow field. This should be contrasted with the curve corresponding to the transition density n = 0.15 (light green curve in Fig. 7), which monotonically increases as k → 0.

image file: c9sm00774a-f7.tif
Fig. 7 Energy spectra Ek, as defined in eqn (23) for different pusher densities n: note the evolving peak for finite k in the turbulent regime. Results were obtained using a system with L = 150. The dashed line shows the low-density prediction of eqn (24), using the dialled value of κ and fitted ε = 1.0.

4.2 Swimmer statistics

In order to further characterise the properties of the suspension we now turn to the local ordering of the swimmers. First of all, there are no significant density inhomogeneities in any of the systems: the pair correlation function g(r) (Fig. 8) has a maximum peak height of ∼1.05 inside the turbulent regime. This is in accordance with previous results24,26,28,34,35 showing that the transition to collective motion is orientational in nature. We thus turn to analyse the orientational order between swimmers, as quantified by the angle θ between the orientations of two swimmers separated by a distance r (Fig. 9). The polar and nematic order parameters were calculated as a function of swimmer-swimmer separation in two different ways: first, for all swimmers a distance r from the central swimmer (solid curves in Fig. 10a and 11), and secondly, for only those swimmers that lie along the axis of the swimming direction p of the central swimmer, i.e. within two cones centered around θ = 0 and θ = π. The polar and nematic order parameters P(r) and S(r) are furthermore defined as
P(r) = 〈P1(cos[thin space (1/6-em)]θ)〉|rirj|=r = 〈cos[thin space (1/6-em)]θr (26)
image file: c9sm00774a-t24.tif(27)
where P1 and P2 are the first and second Legendre polynomials. All the order parameters were calculated both for pushers (solid curves) and pullers (dashed curves). Looking at the full polar order parameter P(r) (Fig. 10a), we observe a weak local polar alignment for pushers and antialignment for pullers, which however converge rather quickly to zero around r ≈ 5. In the case where the polar order parameter was calculated along p (Fig. 10b), the curves for pushers instead fall below zero, showing that the swimmers are weakly aligned in opposite directions along their swimming direction for r ≥ 3. Looking at the corresponding S(r) curves in Fig. 11, for pushers (solid curves) we observe a significantly more pronounced increase in both the magnitude and range of the order parameter when going into the turbulent state. Furthermore, unlike the case of P(r), no significant anisotropy in the nematic order parameter is observed (data for S|| not shown). These observations indicate that far-field hydrodynamics alone is sufficient to induce significant nematic ordering between swimmers, while the effect of the polar symmetry breaking due to self-propulsion is subdominant.

image file: c9sm00774a-f8.tif
Fig. 8 Radial distribution function g(r) for pusher suspensions of four different densities n; for pullers, the corresponding curves (not shown) are completely flat at all densities within the statistical uncertainty.

image file: c9sm00774a-f9.tif
Fig. 9 Schematic image of the orientational correlation functions, where p is the orientation of a swimmer, showing the angles θ and φ used to sample the order parameters P(r) and S(r).

image file: c9sm00774a-f10.tif
Fig. 10 Distance-dependent polar order parameter P(r) as defined in eqn (26) (a) sampled isotropically, and (b) sampled in a cone of angle φ = 20 degrees in front of and behind the swimmer as shown in Fig. 9. Solid curves denote pushers and dashed curves pullers. The apparent kinks at small r are due to poor sampling at small separations.

image file: c9sm00774a-f11.tif
Fig. 11 Distance-dependent nematic order parameter S(r) as defined in eqn (27) for pushers (solid curves) and pullers (dashed curves).

Another approach to calculate the characteristic lengthscales based on the polar and nematic order parameters of the swimmers is to consider the integrated form of the order parameters P(r) and S(r) in eqn (26) and (27), i.e.

image file: c9sm00774a-t25.tif(28)
image file: c9sm00774a-t26.tif(29)
These forms of the order parameter measure the range of polar and nematic order around a single swimmer, in a manner equivalent (modulo a constant shift of unity) to the distance-dependent Kirkwood G-factor employed to measure local orientational order in polar fluids.51 Based on these order parameters, we define the characteristic lengthscales ξP and ξS as the values of R where GP(R) and GS(R) take on their maximum values (see Fig. 12a). The resulting lengthscale curves are shown in Fig. 12b as a function of density. A direct comparison between ξP and the characteristic lengthscale calculated from the fluid (Fig. 5) shows a striking similarity between the curves, modulo a shift in the y direction attributable to the somewhat arbitrary cutoff value used to calculate ξ from the fluid flows. The curve for ξS shows a similar non-monotonic shape as the two other curves. It is, however, shifted towards slightly higher values of ξ compared to the ξP curve, in accordance with the observations made in Fig. 10 and 11, and exhibits somewhat larger statistical fluctuations. Taken together, however, our three separate analyses (Fig. 6 and 12b) of the emerging lengthscales based on either fluid flows or the swimmer orientation provides a consistent picture showing a sharply increasing ξ near the transition, which then plateaus to a finite value in the turbulent regime.

image file: c9sm00774a-f12.tif
Fig. 12 Characteristic lengthscales as measured from the orientational order between swimmers. (a) Cumulative polar order parameter GP, as defined in eqn (28), together with the definition of the corresponding lengthscale ξP; ξS is defined analogously. The y-axis has been shifted by unity for visualisation purposes. (b) ξP and ξS as a function of the microswimmer density n: note the similarity between ξS and ξ as measured from the fluid velocity correlations (Fig. 5). Error bars indicate the estimated standard deviations obtained dividing the simulation into four equally sized time intervals.

5 Conclusions

In this study, we have provided an analysis of the structure and dynamics in suspensions of swimming microorganisms using lattice Boltzmann simulations of model microswimmers. The model, which accurately includes the effect of far-field stresslet flow fields, enables us to simulate large enough systems (N ≈ 3 × 106) to quantitatively study the spatio-temporal properties of the ensuing long-range fluid flows in the turbulent regime, something which has not been possible using more complex models due to the high computational costs. For non-interacting swimmers, we tested the method against a number of analytical predictions for the structural and dynamical observables. In the semidilute regime, where swimmer–swimmer interactions become significant, we showed that the correlation length ξ and correlation time τ of the flows diverge steeply near a swimmer density close to the one predicted from kinetic theory, after a qualitative inclusion of the effect of periodic boundary conditions. Beyond the transition, both these quantities relax to significantly smaller values, indicating the emergence of flows with a finite, well-defined length- and timescale. This was further confirmed in the Fourier-space energy spectra of the flow field, which shows the development of a peak at low k in the turbulent regime. The statistics of the swimmer-swimmer correlations are consistent with the above analysis: the density-dependent lengthscale ξP derived from the local polar ordering between swimmers matches qualitatively that calculated from the flow field, while the corresponding lengthscale ξS calculated from the nematic order parameter is significantly larger than ξP, albeit with a similar shape, indicating a longer range of the local nematic order due to hydrodynamic interactions.

Apart from the results described above, our study highlights the need for employing very large systems when studying collective behaviours in microswimmer suspensions: our finite-size studies indicate that linear system sizes of at least 100 times the swimmer length is necessary to capture the properties of the chaotic flows. This computational efficiency comes at the cost of ignoring effects of near-field hydrodynamic and short-ranged steric interactions between microswimmers, which will become significant at high microswimmer densities. On the other hand, the simplicity of the model provides us with full control of how the different system parameters affect the collective behaviour, and enables a direct comparison with kinetic theories that employ microswimmer models with long-range stresslet flow fields.24 The model can also be extended to include the effect of external gradients, system boundaries, and short-ranged interactions, thus providing further insight into the interplay between microscopic interactions and dynamics, external perturbations, and collective behaviour.

Conflicts of interest

There are no conflicts to declare.


Joost de Graaf, Davide Marenduzzo and Rupert Nash are kindly acknowledged for helpful discussions during the early phases of this project. The work was funded (DB, JS) by the Swedish Research Council (grant ID 2015-05449) and the Crafoord Foundation (grant ID 20170678). CN acknowledges the support of an Aide Investissements d'Avenir du LabEx PALM (ANR-10-LABX-0039-PALM). All simulations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at LUNARC.

Notes and references

  1. M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao and R. A. Simha, Rev. Mod. Phys., 2013, 85, 1143 CrossRef CAS.
  2. C. Bechinger, R. D. Leonardo, H. Löwen, C. Reichhardt, G. Volpe and G. Volpe, Rev. Mod. Phys., 2006, 88, 045006 CrossRef.
  3. W. Bialek, A. Cavagna, I. Giardina, T. Mora, E. Silvestri, M. Viale and A. M. Walczak, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 4786 CrossRef CAS PubMed.
  4. J. Palacci, S. Sacanna, A. P. Steinberg, D. J. Pine and P. M. Chaikin, Science, 2013, 339, 936–940 CrossRef CAS PubMed.
  5. E. Lauga and T. R. Powers, Rep. Prog. Phys., 2009, 72, 096601 CrossRef.
  6. M. E. Cates, Rep. Prog. Phys., 2012, 75, 042601 CrossRef CAS PubMed.
  7. K. C. Leptos, J. S. Guasto, J. P. Gollub, A. I. Pesci and R. E. Goldstein, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 77, 026709 CrossRef PubMed.
  8. A. Jepson, V. A. Martinez, J. Schwartz-Linek, A. Morozov and W. C. K. Poon, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 88, 041002 CrossRef PubMed.
  9. E. F. Semeraro, J. M. Devos and T. Narayanan, J. Chem. Phys., 2018, 148, 204905 CrossRef PubMed.
  10. M. J. Kim and K. S. Breuer, Phys. Fluids, 2004, 16, 78 CrossRef.
  11. G. L. Mino, J. Dunstan, A. Rousselet, E. Clément and R. Soto, J. Fluid Mech., 2013, 729, 423 CrossRef.
  12. A. Creppy, O. Praud, X. Druart, P. L. Kohnke and F. Plouraboué, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2015, 92, 032722 CrossRef PubMed.
  13. L. H. Cisneros, R. Cortez, C. Dombrowski, R. E. Goldstein and J. O. Kessler, Exp. Fluids, 2007, 43, 737 CrossRef.
  14. H. H. Wensink, J. Dunkel, S. Heidenreich, K. Drescher, R. E. Goldstein, H. Löwen and J. M. Yeomans, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 14308–14313 CrossRef CAS PubMed.
  15. J. Dunkel, S. Heidenreich, K. Dreschner, H. H. Wensink, M. Bär and R. E. Goldstein, Phys. Rev. Lett., 2013, 110, 228102 CrossRef PubMed.
  16. D. Saintillan and M. J. Shelley, Phys. Rev. Lett., 2007, 99, 058102 CrossRef PubMed.
  17. G. Pessot, H. Löwen and A. Menzel, Mol. Phys., 2018, 116, 3401–3408 CrossRef CAS.
  18. C. Hoell, H. Löwen and A. M. Menzel, J. Chem. Phys., 2018, 149, 144902 CrossRef PubMed.
  19. N. Yoshinaga and T. B. Liverpool, Phys. Rev. E, 2017, 96, 020603(R) CrossRef PubMed.
  20. T. J. Pedley, IMA J. Appl. Math., 2016, 81, 488–521 CrossRef.
  21. F. Alarcón and I. Pagonabarraga, J. Mol. Liq., 2013, 185, 56–61 CrossRef.
  22. M. Theers, E. Westphal, K. Qi, R. G. Winkler and G. Gompper, Soft Matter, 2018, 14, 8590 RSC.
  23. M. Theers, E. Westphal, G. Gompper and R. G. Winkler, Soft Matter, 2017, 12, 7372 RSC.
  24. J. Stenhammar, C. Nardini, R. W. Nash, D. Marenduzzo and A. Morozov, Phys. Rev. Lett., 2017, 119, 028005 CrossRef PubMed.
  25. J. M. Yeomans, D. O. Pushkin and H. Shum, Eur. Phys. J.: Spec. Top., 2014, 223, 1771–1785 Search PubMed.
  26. D. Saintillan and M. J. Shelley, C. R. Phys., 2013, 14, 497 CrossRef CAS.
  27. R. A. Simha and S. Ramaswamy, Phys. Rev. Lett., 2002, 89, 058101 CrossRef PubMed.
  28. G. Subramanian and D. L. Koch, J. Fluid Mech., 2009, 632, 359 CrossRef.
  29. C. W. Wolgemuth, Biophys. J., 2008, 95, 1564 CrossRef CAS PubMed.
  30. D. L. Koch and G. Subramanian, Annu. Rev. Fluid Mech., 2011, 43, 637–659 CrossRef.
  31. A. Doostmohammadi, J. Ignés-Mullol, J. Yeomans and F. Sagués, Nat. Commun., 2018, 9, 3246 CrossRef PubMed.
  32. S. P. Thampi, R. Golestanian and J. M. Yeomans, Phys. Rev. Lett., 2013, 111, 118101 CrossRef PubMed.
  33. S. Ramaswamy, Annu. Rev. Fluid Mech., 2010, 1, 323–345 Search PubMed.
  34. D. Saintillan and M. J. Shelley, Phys. Rev. Lett., 2008, 100, 178103 CrossRef PubMed.
  35. C. Hohenegger and M. J. Shelley, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2010, 81, 046311 CrossRef PubMed.
  36. M. Theillard, R. Alonso-Matilla and D. Saintillan, Soft Matter, 2017, 13, 363–375 RSC.
  37. B. Ezhilan, M. J. Shelley and D. Saintillan, Phys. Fluids, 2013, 25, 070607 CrossRef.
  38. P. T. Underhill, J. P. Hernandez-Ortiz and M. D. Graham, Phys. Rev. Lett., 2008, 100, 248101 CrossRef PubMed.
  39. J. P. Hernandez-Ortiz, P. T. Underhill and M. D. Graham, J. Phys.: Condens. Matter, 2009, 21, 204107 CrossRef PubMed.
  40. E. Lushi and C. P. Peskin, Comput. Struct., 2013, 122, 239 CrossRef.
  41. D. Saintillan and M. J. Shelley, J. R. Soc., Interface, 2012, 9, 571–585 CrossRef PubMed.
  42. D. Krishnamurthy and G. Subramanian, J. Fluid Mech., 2015, 781, 422 CrossRef CAS.
  43. R. W. Nash, R. Adhikari and M. E. Cates, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 77, 026709 CrossRef CAS PubMed.
  44. D. Saintillan and M. J. Shelley, Phys. Fluids, 2008, 20, 123304 CrossRef.
  45. R. W. Nash, R. Adhikari, J. Tailleur and M. E. Cates, Phys. Rev. Lett., 2010, 104, 258101 CrossRef CAS PubMed.
  46. C. S. Peskin, Acta Numerica, 2002, 11, 479–517 CrossRef.
  47. J. de Graaf and J. Stenhammar, Phys. Rev. E, 2017, 95, 023302 CrossRef PubMed.
  48. K. Drescher, J. Dunkel, L. H. Cisneros, S. Ganguly and R. E. Goldstein, Proc. Natl. Acad. Sci. U. S. A., 2011, 108, 10940 CrossRef CAS PubMed.
  49. R. Cortez, L. Fauci and A. Medovikov, Phys. Fluids, 2005, 17, 031504 CrossRef.
  50. I. M. Zaid, J. Dunkel and J. M. Yeomans, J. R. Soc., Interface, 2011, 8, 1314–1331 CrossRef PubMed.
  51. C. Böttcher, Theory of electric polarization, Elsevier, Netherlands, 1973 Search PubMed.


Electronic supplementary information (ESI) available: Videos showing the swimmers and the fluid in the disordered regime, transition regime and in the turbulent regime. See DOI: 10.1039/c9sm00774a
An alternative way of non-dimensionalising the time units is to use the characteristic tumbling time λ−1, which would be the relevant timescale in a suspension of non-swimming stresslets (“shakers”).24

This journal is © The Royal Society of Chemistry 2019