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

Abstract

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 × 10⁶) 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 flocks³ and dynamical clustering in suspensions of active colloids.⁴

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 spherical^19–21 or elongated^22,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 n_c required for collective motion:^{24,26,28,34,35}

n_c = 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 Peskin⁴⁰ 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 Shelley⁴¹ and Krishnamurthy and Subramanian⁴² observed similar collective properties in suspensions containing up to N = 3 × 10⁴ 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 > 10⁶ 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) equation⁴³ 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.

Fig. 1 Schematic image of the pusher model, where F is the force, l is the swimmer length, v_s 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

ṙ = v_sp + U(r), (2)

Here, U(r) is the fluid velocity evaluated at the body position, v_s is the constant swimming velocity of the individual swimmer, and 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).²⁸) This run-and-tumble motion results in a random walk with a persistence length v_s/λ.

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, v_s, and a:²⁴

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 ∼ 10⁶ microswimmers at biologically relevant densities. In the simulations, cubic box sizes L³ ranging from (10)³ to (210)³ lattice sites were employed. In terms of LB units (where ΔL = Δt = 1), we used the parameters v_s = 10⁻³, F = 1.57 × 10⁻³, l = 1, λ = 2 × 10⁻⁴, and μ = 1/6, with the latter value corresponding to the fluid relaxing to local equilibrium on each timestep. Each simulation was run for 2 × 10⁵ 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/v_s.‡ The relevant dimensionless numbers of the system are (i) the single-swimmer Reynolds number Re_s ≡ ρ_fv_sl/μ, where ρ_f is the fluid density (set to unity in our LB simulations), and (ii) the nondimensional stresslet strength κ_n ≡ κ/(l²v_s) = F/(μlv_s). The above parameter values yield Re_s = 6 × 10⁻³, which is well below the Stokes-flow limit,⁴⁷ and κ_n ≈ 9.4. The latter value can be compared with the corresponding value for E. coli for which v_s ≈ 22 μm s⁻¹, F = 0.42 pN, and l = 1.9 μm,⁴⁸ 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)a³n, where n = N/L³ 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 achievable²⁴ but considerably more complicated; we therefore postpone the full analytical treatment of interacting swimmers to a forthcoming separate study.

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

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 bywhere K₂ is the modified Bessel function of the second kind, and k̂_i = k_i/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 r₀ and r₀ + lp, where l is the dipolar length and p its orientation. The corresponding real-space velocity field of a single swimmer is thusThe total velocity field at a position r created by a suspension of N non-interacting swimmers with instantaneous positions r_i and orientations p_i, i = 1…N is, then,

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 averageHere, 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 v_s is the swimming speed. Since the swimmers are statistically independent, only ‘self-correlations’ contribute to the average, yieldingwhereHere, as before, n = N/V, is the swimmer number density. Explicit integration giveswherewith [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 t_n = tv_s/l, tumbling rate λ_n = λl/v_s, and regularisation parameter Δ = l/ε, and we have normalised c(t_n) so that c(0) = 1.

3.2 Spatial correlations

Similar to the temporal correlation function, the spatial velocity–velocity correlation function c(R) = 〈U(0)·U(R)〉 is given by the following ensemble averageAgain, keeping only the ‘self-correlation’ contributions, we obtainDirect evaluation of this integral yieldswhereandHere, R_n = R/l, and we have again normalised c(R_n) such that c(0) = 1. For large R_n, c(R_n) ∼ R_n⁻¹, in agreement with the results of Zaid et al.⁵⁰

3.3 Velocity variance

The velocity variance 〈U²〉 can be obtained from eqn (15) by setting R = 0. The result depends on c₀ and readsBy 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:

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)The energy content associated with the velocity field at a lengthscale k⁻¹ is given by

E_k = 4πk²〈Û(k)·Û(−k)〉_k=|k|, (23)

yieldingIn the double limit ε → 0 and l → 0, this expression reduces towhich 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 U_RMS ≡ 〈U²〉^1/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 n^1/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.⁴⁶ At intermediate concentrations for interacting swimmers, there is a deviation from the square root dependence, with the RMS velocity increasing faster than n^1/2 for pushers and slower for pullers: this is a signature of the build-up of long-ranged orientational correlations underlying the collective behaviour.²⁴ 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.

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.†

Fig. 3 Root-mean-square fluid velocity U_RMS, 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(t_n) 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 R_n ≈ 10, eventually falling below zero at R_n ≈ 80. We attribute this poor agreement at intermediate and large R_n 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(R_n) has decreased to 0.2 (Fig. 4a), whereas the corresponding characteristic time τ is defined as the point when c(t_n) has decayed to 0.4 (see Fig. 4b). The difference in the two threshold values is due to the slow decay of c(t_n) 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,⁴¹ 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.⁴² We attribute this difference to the significantly smaller system sizes used in earlier studies.

Fig. 4 Build-up of long-ranged velocity correlations due to collective motion. (a) The spatial velocity correlation function c(R_n), and (b) velocity autocorrelation functions c(t_n) 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).

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 n_c for an unbounded system (k_c = 0, eqn (1)), and for the finite wavenumbers k_c = 2π/L, 4π/L, and 8π/L, as discussed in the text.

We furthermore note that the predicted infinite-system critical density n_c (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 k_c = 0 to k_c = 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⁻² 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⁻² 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 n_c corresponding to the critical wavenumbers k_c = 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 × 10⁵ 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.

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 E_k 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⁻² 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, E_k 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.

Fig. 7 Energy spectra E_k, 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 results^{24,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) = 〈P₁(cos θ)〉_|ri−rj|=r = 〈cos θ〉_r (26)

andwhere P₁ and P₂ 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.

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.

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).

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.

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.

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.⁵¹ Based on these order parameters, we define the characteristic lengthscales ξ_P and ξ_S as the values of R where G_P(R) and G_S(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.

Fig. 12 Characteristic lengthscales as measured from the orientational order between swimmers. (a) Cumulative polar order parameter G_P, 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 × 10⁶) 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.²⁴ 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.

Acknowledgements

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

---

Footnotes

† 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 λ⁻¹, which would be the relevant timescale in a suspension of non-swimming stresslets (“shakers”).²⁴

---