Open Access Article

This Open Access Article is licensed under a

Creative Commons Attribution 3.0 Unported Licence

DOI: 10.1039/C9SM00774A
(Paper)
Soft Matter, 2019, Advance Article

Dóra Bárdfalvy^{a},
Henrik Nordanger^{a},
Cesare Nardini^{b},
Alexander Morozov^{c} and
Joakim Stenhammar*^{a}
^{a}Division of Physical Chemistry, Lund University, Box 124, S-221 00 Lund, Sweden. E-mail: joakim.stenhammar@fkem1.lu.se
^{b}Service de Physique de l'État Condensé, CNRS UMR 3680, CEA-Saclay, 91191 Gif-sur-Yvette, France
^{c}SUPA, 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 × 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.

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

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.

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_{s}p + U(r),
| (2) |

(3) |

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}

(4) |

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^{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. Each simulation was run for 2 × 10^{5} 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} ≡ ρ_{f}v_{s}l/μ, where ρ_{f} is the fluid density (set to unity in our LB simulations), and (ii) the nondimensional stresslet strength κ_{n} ≡ κ/(l^{2}v_{s}) = F/(μlv_{s}). The above parameter values yield Re_{s} = 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 v_{s} ≈ 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)a^{3}n, where n = N/L^{3} 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.

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

(5) |

(6) |

(7) |

(8) |

(9) |

(10) |

(11) |

(12) |

(13) |

(14) |

(15) |

(16) |

(17) |

(18) |

(19) |

(20) |

(21) |

(22) |

E_{k} = 4πk^{2}〈Û(k)·Û(−k)〉_{k=|k|},
| (23) |

(24) |

(25) |

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,^{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.

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^{−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,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^{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 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, 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. |

P(r) = 〈P_{1}(cosθ)〉_{|ri−rj|=r} = 〈cosθ〉_{r}
| (26) |

(27) |

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.

(28) |

(29) |

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

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.

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

## 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 λ^{−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 |