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

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\times 10^6$) 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.


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 selfpropulsion. 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 3 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][8][9][10][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][13][14][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][20][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][25][26][27][28][29][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][32][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 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. 40,41 and Lushi and Peskin 42 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 38 and Krishnamurthy and Subramanian 39 observed similar collective properties in suspensions containing up to N = 3 × 10 4 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 6 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 43 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.

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 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 viscosity of the fluid, with κ > 0 representing pushers and κ < 0 pullers. The position r and orientation p of each swimmer evolves according to the following equations of motion: Here, U(r) is the fluid velocity evaluated at the body position, v s is the constant swimming velocity of the individual swimmer, and I is the unit tensor. Equation (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 Eq. (3), swimmers also undergo Poissondistributed random reorientations with an average frequency λ . 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: 24 Here, the second term in brackets constitutes the leadingorder 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,44 The key ingredient in this method is an algorithm for interpolating forces and velocities between the fluid and the offlattice swimmers that employs the regularised version of the δ function due to Peskin. 43,49 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 6 microswimmers at biologically relevant densities. In the simulations, cubic box sizes L 3 ranging from (10) 3 to (210) 3 lattice sites were employed. In terms of LB units (where ∆L = ∆t = 1), we used the parameters v s = 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.
All results will be presented in terms of the swimmer length

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 Eqs. (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 Fourierspace 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 24 but considerably more complicated; we therefore postpone the full analytical treatment of interacting swimmers to a forthcoming separate study.
Following Cortez et al., 47 we start from the expression for * 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 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 bŷ where K 2 is the modified Bessel function of the second kind, 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 0 and r 0 + lp, where l is the dipolar length and p its orientation. The corresponding real-space velocity field of a single swimmer is thus The 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,

Temporal correlations
First we consider the velocity-velocity autocorrelation function c(t) = U(0) · U(t) , which is given by the following ensemble average 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 v s is the swimming speed. Since the swimmers are statistically independent, only 'self-correlations' contribute to the average, yielding where Here, as before, n = N/V , is the swimmer number density. Explicit integration gives where with E(x) and 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.

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 average Again, keeping only the 'self-correlation' contributions, we obtain Direct evaluation of this integral yields c(R n ) = 1 c 0 R n g 2 (R n + 1) − 10∆ 2 4 + ∆ 2 R 2 n g 2 (R n ) + g 2 (R n − 1) where and Here, 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 −1 n , in agreement with the results of Zaid et al. 48

Velocity variance
The velocity variance U 2 can be obtained from Eq.(15) by setting R = 0. The result depends on c 0 and reads 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:

Energy spectrum
Eq.(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 −1 is given by yielding In the double limit ε → 0 and l → 0, this expression reduces to which is independent of the wavevector k.  Figure 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 2 1/2 for suspensions of pushers, pullers and noninteracting swimmers; for the latter simulations the terms containing U in Eqs. (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) Eq. (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. 49 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. 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.

Fluid statistics
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 (Eq. (12)), the corresponding spatial curves start deviating from the predictions (Eq. (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 Figure 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, 38 who did not observe any maximum at the transition density for neither the caracteristic length or time-scales, while the non-monotonic behaviour of τ(n) was previously observed by Krishnamurthy and Subramanian. 39 We attribute this difference to the significantly smaller system sizes used in earlier studies. We furthermore note that the predicted infinite-system critical density n c (Eq. (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 −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 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,50 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 5 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.
To further analyse the spatial structures of the flow, in Fig. 7 we calculate the Fourier space energy spectrum E k as defined in Eq. (23). For low densities, the spectrum (Fig. 7) is well described by the form predicted for uncorrelated swimmers in Eq. (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, E k decreases due to the short-range regularisation of the flow field; the slight discrepancy in be- tween 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, 38,39 , 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.

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 lines 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 1 (cos θ ) |r i −r j |=r = cos θ r (26) and S(r) = P 2 (cos θ ) |r i −r j |=r = 3 cos 2 θ − 1 2 r where P 1 and P 2 are the first and second Legendre polynomials. All the order parameters were calculated both for pushers (solid lines) and pullers (dashed lines). 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 lines) 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 sym- Fig. 9 Schematic image of the orientational correlation functions, showing the angles θ and ϕ used to sample the order parameters P(r) and S(r).
metry breaking due to self-propulsion is subdominant. 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 Eqs. (26)-(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 distancedependent Kirkwood G-factor employed to measure local order in polar fluids. 51 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 Figs. 10 and 11, and exhibits somewhat larger statistical fluctuations. Taken together, however, our three separate analyses (Figs. 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.

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 6 ) 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 noninteracting swimmers, we tested the method against a number  ξ n ξ P ξ S

Fig. 12
Characteristic lengthscales as measured from the orientational order between swimmers. (a) Cumulative polar order parameter G P , as defined in Eq. (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 from sub-averaging each simulation over four separate time intervals.
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 densitydependent 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 finitesize 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 shortranged 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.