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

Dora Bardfalvy; Henrik Nordanger; Cesare Nardini; Alexander Morozov; and Joakim Stenhammar

arXiv:1904.03069·cond-mat.soft·August 27, 2019

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

Dora Bardfalvy, Henrik Nordanger, Cesare Nardini, Alexander Morozov, and Joakim Stenhammar

PDF

TL;DR

This paper uses large-scale particle-resolved lattice Boltzmann simulations to study active turbulence in suspensions of microswimmers, revealing how collective chaotic flows emerge and evolve with swimmer density.

Contribution

It provides the first large-scale simulations of active turbulence with detailed hydrodynamic interactions, validating analytical predictions and exploring transition behavior at high densities.

Findings

01

Transition to turbulence occurs at predicted densities.

02

Flow length- and timescales increase steeply near transition.

03

Simulation results agree with analytical predictions in dilute limit.

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 \times 1 0^{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…

Equations63

n_{c} = 5 λ / κ .

n_{c} = 5 λ / κ .

\dot{r} = v_{s} p + U (r),

\dot{r} = v_{s} p + U (r),

\dot{p} = (I - pp) \cdot (\nabla U) \cdot p \approx (I - pp) \cdot \frac{U ( r ) - U ( r - p l )}{l} .

v_{s} = \frac{F}{6 π μ a} (1 - \frac{3 a}{2 l}) .

v_{s} = \frac{F}{6 π μ a} (1 - \frac{3 a}{2 l}) .

u_{i} (r) = \frac{( r ^{2} + 2 ϵ ^{2} ) δ _{ij} + r _{i} r _{j}}{( r ^{2} + ϵ ^{2} ) ^{3/2}} \frac{F _{j}}{8 π μ} .

u_{i} (r) = \frac{( r ^{2} + 2 ϵ ^{2} ) δ _{ij} + r _{i} r _{j}}{( r ^{2} + ϵ ^{2} ) ^{3/2}} \frac{F _{j}}{8 π μ} .

\overset{u}{^}_{i} (k) = (δ_{ij} - \hat{k}_{i} \hat{k}_{j}) F_{j} \frac{ϵ ^{2} K _{2} ( k ϵ )}{2 μ},

\overset{u}{^}_{i} (k) = (δ_{ij} - \hat{k}_{i} \hat{k}_{j}) F_{j} \frac{ϵ ^{2} K _{2} ( k ϵ )}{2 μ},

u_{i} (r; r_{0}, p) = \frac{1}{( 2 π ) ^{3}} \int d k e^{i k \cdot (r - r_{0})} [e^{- i k \cdot l p} - 1]

u_{i} (r; r_{0}, p) = \frac{1}{( 2 π ) ^{3}} \int d k e^{i k \cdot (r - r_{0})} [e^{- i k \cdot l p} - 1]

\times (δ_{ij} - \hat{k}_{i} \hat{k}_{j}) p_{j} \frac{F ϵ ^{2} K _{2} ( k ϵ )}{2 μ} .

U (r) = i = 1 \sum N u (r; r_{i}, p_{i}) .

U (r) = i = 1 \sum N u (r; r_{i}, p_{i}) .

c (t) = \frac{e ^{- λ t}}{V} \int d r \frac{1}{V ^{N}} \int d r_{1} \dots d r_{N} \frac{1}{( 4 π ) ^{N}} \int d p_{1} \dots d p_{N}

c (t) = \frac{e ^{- λ t}}{V} \int d r \frac{1}{V ^{N}} \int d r_{1} \dots d r_{N} \frac{1}{( 4 π ) ^{N}} \int d p_{1} \dots d p_{N}

\times (j = 1 \sum N u (r; r_{j}, p_{j})) \cdot (k = 1 \sum N u (r; r_{k} + v_{s} t p_{k}, p_{k})) .

c (t) = \frac{n}{2 π ^{2}} \int_{0}^{\infty} d k k^{2} I_{1} (k) (\frac{F ϵ ^{2} K _{2} ( k ϵ )}{2 μ})^{2},

c (t) = \frac{n}{2 π ^{2}} \int_{0}^{\infty} d k k^{2} I_{1} (k) (\frac{F ϵ ^{2} K _{2} ( k ϵ )}{2 μ})^{2},

I_{1} (k) = \int_{- 1}^{1} d x e^{- λ t - ik v_{s} t x} [1 - cos (k l x)] (1 - x^{2}) .

I_{1} (k) = \int_{- 1}^{1} d x e^{- λ t - ik v_{s} t x} [1 - cos (k l x)] (1 - x^{2}) .

c (t_{n}) = e^{- λ_{n} t_{n}} \frac{g _{1} ( t _{n} + 1 ) - 2 g _{1} ( t _{n} ) + g _{1} ( t _{n} - 1 )}{2 g _{1} ( 1 ) - \frac{5 π Δ ^{2}}{8}},

c (t_{n}) = e^{- λ_{n} t_{n}} \frac{g _{1} ( t _{n} + 1 ) - 2 g _{1} ( t _{n} ) + g _{1} ( t _{n} - 1 )}{2 g _{1} ( 1 ) - \frac{5 π Δ ^{2}}{8}},

g_{1} (x) = \frac{1 + x ^{2} Δ ^{2}}{x ^{2}} E (- \frac{1}{4} x^{2} Δ^{2}) - \frac{1 + \frac{1}{2} x ^{2} Δ ^{2}}{x ^{2}} K (- \frac{1}{4} x^{2} Δ^{2}),

g_{1} (x) = \frac{1 + x ^{2} Δ ^{2}}{x ^{2}} E (- \frac{1}{4} x^{2} Δ^{2}) - \frac{1 + \frac{1}{2} x ^{2} Δ ^{2}}{x ^{2}} K (- \frac{1}{4} x^{2} Δ^{2}),

c (R) = \frac{1}{V} \int d r \frac{1}{V ^{N}} \int d r_{1} \dots d r_{N} \frac{1}{( 4 π ) ^{N}} \int d p_{1} \dots d p_{N}

c (R) = \frac{1}{V} \int d r \frac{1}{V ^{N}} \int d r_{1} \dots d r_{N} \frac{1}{( 4 π ) ^{N}} \int d p_{1} \dots d p_{N}

\times U (r) \cdot U (r + R) .

c (R) = \frac{4 n}{( 2 π ) ^{3}}

c (R) = \frac{4 n}{( 2 π ) ^{3}}

\times \int d k e^{- i k \cdot R} [\frac{1}{3} + \frac{cos k l}{( k l ) ^{2}} - \frac{sin k l}{( k l ) ^{3}}] (\frac{F ϵ ^{2} K _{2} ( k ϵ )}{2 μ})^{2} .

c (R_{n}) = \frac{1}{c _{0} R _{n}} [g_{2} (R_{n} + 1) - \frac{10 Δ ^{2}}{4 + Δ ^{2} R _{n}^{2}} g_{2} (R_{n}) + g_{2} (R_{n} - 1)

c (R_{n}) = \frac{1}{c _{0} R _{n}} [g_{2} (R_{n} + 1) - \frac{10 Δ ^{2}}{4 + Δ ^{2} R _{n}^{2}} g_{2} (R_{n}) + g_{2} (R_{n} - 1)

g_{3} (R_{n} + 1) - \frac{10 Δ ^{2}}{R _{n} ( 4 + Δ ^{2} R _{n}^{2} )} g_{3} (R_{n}) - g_{3} (R_{n} - 1)],

g_{2} (x) = 10 x Δ^{2} [2 (5 + x^{2} Δ^{2}) E (- \frac{1}{4} x^{2} Δ^{2})

g_{2} (x) = 10 x Δ^{2} [2 (5 + x^{2} Δ^{2}) E (- \frac{1}{4} x^{2} Δ^{2})

- (4 + x^{2} Δ^{2}) K (- \frac{1}{4} x^{2} Δ^{2})],

g_{3} (x) = 2 (x^{4} Δ^{4} + 2 x^{2} Δ^{2} - 8) K (- \frac{1}{4} x^{2} Δ^{2})

- 4 (x^{4} Δ^{4} + 9 x^{2} Δ^{2} - 4) E (- \frac{1}{4} x^{2} Δ^{2}),

c_{0} = 15 Δ^{2} [8 (1 + Δ^{2}) E (- \frac{1}{4} Δ^{2})

c_{0} = 15 Δ^{2} [8 (1 + Δ^{2}) E (- \frac{1}{4} Δ^{2})

- 4 (2 + Δ^{2}) K (- \frac{1}{4} Δ^{2}) - \frac{5 π}{2} Δ^{2}] .

⟨ U^{2} ⟩ = \frac{κ ^{2} n}{480 π} \frac{c _{0}}{ϵ Δ ^{6}} .

⟨ U^{2} ⟩ = \frac{κ ^{2} n}{480 π} \frac{c _{0}}{ϵ Δ ^{6}} .

⟨ U^{2} ⟩ \approx \frac{21 n κ ^{2}}{8 ( 21 π l + 256 ϵ )} .

⟨ U^{2} ⟩ \approx \frac{21 n κ ^{2}}{8 ( 21 π l + 256 ϵ )} .

⟨ \hat{U} (k) \cdot \hat{U} (- k)⟩ = 4 n [\frac{1}{3} + \frac{cos k l}{( k l ) ^{2}} - \frac{sin k l}{( k l ) ^{3}}] (\frac{F ϵ ^{2} K _{2} ( k ϵ )}{2 μ})^{2} .

⟨ \hat{U} (k) \cdot \hat{U} (- k)⟩ = 4 n [\frac{1}{3} + \frac{cos k l}{( k l ) ^{2}} - \frac{sin k l}{( k l ) ^{3}}] (\frac{F ϵ ^{2} K _{2} ( k ϵ )}{2 μ})^{2} .

E_{k} = 4 π k^{2} ⟨ \hat{U} (k) \cdot \hat{U} (- k) ⟩_{k = ∣ k ∣},

E_{k} = 4 π k^{2} ⟨ \hat{U} (k) \cdot \hat{U} (- k) ⟩_{k = ∣ k ∣},

E_{k} = 4 π n κ^{2} [\frac{1}{3} + \frac{cos ( k l )}{( k l ) ^{2}} - \frac{sin ( k l )}{( k l ) ^{3}}] \frac{ϵ ^{4} k ^{2}}{l ^{2}} K_{2}^{2} (k ϵ) .

E_{k} = 4 π n κ^{2} [\frac{1}{3} + \frac{cos ( k l )}{( k l ) ^{2}} - \frac{sin ( k l )}{( k l ) ^{3}}] \frac{ϵ ^{4} k ^{2}}{l ^{2}} K_{2}^{2} (k ϵ) .

E_{k} \approx \frac{8 π}{15} n κ^{2},

E_{k} \approx \frac{8 π}{15} n κ^{2},

P (r) = ⟨ P_{1} (cos θ) ⟩_{∣ r_{i} - r_{j} ∣ = r} = ⟨ cos θ ⟩_{r}

P (r) = ⟨ P_{1} (cos θ) ⟩_{∣ r_{i} - r_{j} ∣ = r} = ⟨ cos θ ⟩_{r}

S(r)=\langle P_{2}(\cos\theta)\rangle_{|\mathbf{r}_{i}-\mathbf{r}_{j}|=r}=\Bigg{\langle}\frac{3\cos^{2}\theta-1}{2}\Bigg{\rangle}_{r}

S(r)=\langle P_{2}(\cos\theta)\rangle_{|\mathbf{r}_{i}-\mathbf{r}_{j}|=r}=\Bigg{\langle}\frac{3\cos^{2}\theta-1}{2}\Bigg{\rangle}_{r}

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

\makeFNbottom

[TABLE]

††footnotetext: a* Division of Physical Chemistry, Lund University, Box 124, S-221 00 Lund, Sweden. E-mail: [email protected]*††footnotetext: b* Service de Physique de l’État Condensé, CNRS UMR 3680, CEA-Saclay, 91191 Gif-sur-Yvette, France. ††footnotetext: c SUPA, School of Physics and Astronomy, The University of Edinburgh, James Clerk Maxwell Building, Peter Guthrie Tait Road, Edinburgh, EH9 3FD, United Kingdom. *††footnotetext: † 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/b000000x/

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 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. 21, 19 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,16, 29, 28, 24, 1 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:26, 28, 24, 34, 35

[TABLE]

Here, $\lambda$ is the tumbling frequency by which individual swimmers randomise their swimming direction and $\kappa$ 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. 34, 36, 37, 31 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 40 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 41 and Krishnamurthy and Subramanian 42 observed similar collective properties in suspensions containing up to $N=3\times 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 $\sim 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 $\pm F\mathbf{p}$ , where $\mathbf{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 $\kappa=\pm Fl/\mu$ where $\mu$ is the dynamic viscosity of the fluid, with $\kappa>0$ representing pushers and $\kappa<0$ pullers.

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

[TABLE]

Here, $\mathbf{U(r)}$ is the fluid velocity evaluated at the body position, $v_{s}$ is the constant swimming velocity of the individual swimmer, and $\mathbb{I}$ is the unit tensor. Equation (3) is the (discretised version of) Jeffery’s equation for infinite aspect ratio ( $\beta=1$ ); we checked our results also using a finite aspect ratio ( $\beta<1$ ), which yielded only small shifts for realistic values of $\beta$ . In addition to their reorientation as a result of the hydrodynamic interactions described by Eq. (3), swimmers also undergo Poisson-distributed random reorientations with an average frequency $\lambda$ . (Note that this mechanism cannot be cast as a continuous-time differential equation as that in Eq. (3). 28) This run-and-tumble motion results in a random walk with a persistence length $v_{s}/\lambda$ .

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

[TABLE]

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 $\delta$ 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\sim 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 $\Delta L=\Delta t=1$ ), we used the parameters $v_{s}=10^{-3}$ , $F=1.57\times 10^{-3}$ , $l=1$ , $\lambda=2\times 10^{-4}$ , and $\mu=1/6$ , with the latter value corresponding to the fluid relaxing to local equilibrium on each timestep. Each simulation was run for $2\times 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 $\sim 48$ hours on a single CPU, while the largest simulations ( $L=210$ ) took $\sim 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}$ .***An alternative way of non-dimensionalising the time units is to use the characteristic tumbling time $\lambda^{-1}$ , which would be the relevant timescale in a suspension of non-swimming stresslets (“shakers”).24 The relevant dimensionless numbers of the system are (i) the single-swimmer Reynolds number $\mathrm{Re}_{s}\equiv\rho_{\mathrm{f}}v_{s}l/\mu$ , where $\rho_{\mathrm{f}}$ is the fluid density (set to unity in our LB simulations), and (ii) the nondimensional stresslet strength $\kappa_{n}\equiv\kappa/(l^{2}v_{s})=F/(\mu lv_{s})$ . The above parameter values yield $\mathrm{Re}_{s}=6\times 10^{-3}$ , which is well below the Stokes-flow limit,47 and $\kappa_{n}\approx 9.4$ . The latter value can be compared with the corresponding value for E. coli for which $v_{s}\approx 22$ $\mathrm{\mu ms^{-1}}$ , $F=0.42$ pN, and $l=1.9$ $\mathrm{\mu m}$ , 48 yielding $\kappa_{n}\approx 11.2$ . Furthermore, an approximate volume fraction based on the spherical swimmer body can be calculated with the aid of Eq. (4) as $\phi=(4\pi/3)a^{3}n$ , where $n=N/L^{3}$ is the swimmer number density; with our parameters, Eq. (4) gives $a\approx 0.3$ . Using this estimate, the observed critical volume fraction needed for collective motion is $\phi_{c}\approx 0.02$ , which is clearly within the range of validity of our far-field hydrodynamic model. While the highest densities considered in this paper ( $\phi\approx 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 $\kappa$ is large enough, in accordance with Eq. (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 $\mathbf{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 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 24 but considerably more complicated; we therefore postpone the full analytical treatment of interacting swimmers to a forthcoming separate study.

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

[TABLE]

Here, $\epsilon$ is a factor describing the distance over which the regularisation acts and $r=|\mathbf{r}|$ ; in the limit $\epsilon\rightarrow 0$ this expression reduces to the usual Stokeslet. Its Fourier transform is given by

[TABLE]

where $K_{2}$ is the modified Bessel function of the second kind, and $\hat{k}_{i}=k_{i}/k$ , where $k=|{\bf 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 $\mathbf{r}_{0}$ and $\mathbf{r}_{0}+l\mathbf{p}$ , where $l$ is the dipolar length and $\bf p$ its orientation. The corresponding real-space velocity field of a single swimmer is thus

[TABLE]

The total velocity field at a position ${\bf r}$ created by a suspension of $N$ non-interacting swimmers with instantaneous positions ${\bf r}_{i}$ and orientations ${\bf p}_{i}$ , $i=1\dots N$ is, then,

[TABLE]

3.1 Temporal correlations

First we consider the velocity-velocity autocorrelation function $c(t)=\langle{\bf U}(0)\cdot{\bf U}(t)\rangle$ , which is given by the following ensemble average

[TABLE]

Here, we used the fact that the time-dependence of the velocity field only arises through changes of the swimmer positions; the factor $e^{-\lambda 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

[TABLE]

where

[TABLE]

Here, as before, $n=N/V$ , is the swimmer number density. Explicit integration gives

[TABLE]

where

[TABLE]

with $\mathbb{E}(x)$ and $\mathbb{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 $\lambda_{n}=\lambda l/v_{s}$ , and regularisation parameter $\Delta=l/\epsilon$ , 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)=\langle{\bf U}({\bf 0})\cdot{\bf U}({\bf R})\rangle$ is given by the following ensemble average

[TABLE]

Again, keeping only the ’self-correlation’ contributions, we obtain

[TABLE]

Direct evaluation of this integral yields

[TABLE]

where

[TABLE]

and

[TABLE]

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})\sim R_{n}^{-1}$ , in agreement with the results of Zaid et al. 50

3.3 Velocity variance

The velocity variance $\langle U^{2}\rangle$ can be obtained from Eq.(15) by setting $R=0$ . The result depends on $c_{0}$ and reads

[TABLE]

By combining the large- and small- $\Delta$ 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:

[TABLE]

3.4 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 ${\bf k}\rightarrow{-\bf k}$ due to our definition of the Fourier transform)

[TABLE]

The energy content associated with the velocity field at a lengthscale $k^{-1}$ is given by

[TABLE]

yielding

[TABLE]

In the double limit $\epsilon\rightarrow 0$ and $l\rightarrow 0$ , this expression reduces to

[TABLE]

which is independent of the wavevector $k$ .

4 Results and discussion

4.1 Fluid statistics

Figure 3 demonstrates how collective motion develops as a function of the concentration of microswimmers through the root-mean-square (RMS) fluid velocity $U_{\mathrm{RMS}}\equiv\langle\mathbf{U}^{2}\rangle^{1/2}$ for suspensions of pushers, pullers and noninteracting swimmers; for the latter simulations the terms containing $\mathbf{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 $\epsilon\approx 1.1$ (while fixing all other parameters to the values from the LB simulations), in good accordance with the interpolation length of the Peskin $\delta$ function used in the LB simulations. 46 At intermediate concentrations for interacting swimmers, there is a deviation from the square root dependence, with the RMS velocity increasing faster than $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.

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. (3.2)) for $R_{n}\approx 10$ , eventually falling below zero at $R_{n}\approx 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 $\xi$ as the distance where $c(R_{n})$ has decreased to $0.2$ (Fig. 4a), whereas the corresponding characteristic time $\tau$ 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 $\xi$ and $\tau$ 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\approx 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 $\tau$ curve, while the maximum in $\xi$ 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 caracteristic length or time-scales, while the non-monotonic behaviour of $\tau(n)$ was previously observed by Krishnamurthy and Subramanian. 42 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 $\xi$ and $\tau$ . 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\pi/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\pi/L$ , $4\pi/L$ and $8\pi/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 $\xi$ and $\tau$ 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\approx 100$ , corresponding to $N=3\times 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 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\rightarrow 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 $\sim 1.05$ inside the turbulent regime. This is in accordance with previous results 26, 28, 24, 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 $\theta$ 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 $\mathbf{p}$ of the central swimmer, i.e. within two cones centered around $\theta=0$ and $\theta=\pi$ . The polar and nematic order parameters $P(r)$ and $S(r)$ are furthermore defined as

[TABLE]

and

[TABLE]

where $P_{1}$ and $P_{2}$ 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\approx 5$ . In the case where the polar order parameter was calculated along $\mathbf{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\geq 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.

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.

[TABLE]

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 order in polar fluids. 51 Based on these order parameters, we define the characteristic lengthscales $\xi_{P}$ and $\xi_{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 $\xi_{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 $\xi$ from the fluid flows. The curve for $\xi_{S}$ shows a similar non-monotonic shape as the two other curves. It is, however, shifted towards slightly higher values of $\xi$ compared to the $\xi_{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 $\xi$ near the transition, which then plateaus to a finite value in the turbulent regime.

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\approx 3\times 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 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 $\xi$ and correlation time $\tau$ 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 $\xi_{P}$ derived from the local polar ordering between swimmers matches qualitatively that calculated from the flow field, while the corresponding lengthscale $\xi_{S}$ calculated from the nematic order parameter is significantly larger than $\xi_{P}$ , albeit with a similar shape, indicating a longer range of the local nematic order due to hydrodynamic interactions.

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

Conflicts of interest

There are no conflicts to declare.

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

Marchetti et al. 2013

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.

Bechinger et al. 2006

C. Bechinger, R. D. Leonardo, H. Löwen, C. Reichhardt, G. Volpe and G. Volpe, Rev. Mod. Phys., 2006, 88, 045006.

Bialek et al. 2012

W. Bialek, A. Cavagna, I. Giardina, T. Mora, E. Silvestri, M. Viale and A. M. Walczak, Proc. Natl. Acad. Sci. USA, 2012, 109, 4786.

Palacci et al. 2013

J. Palacci, S. Sacanna, A. P. Steinberg, D. J. Pine and P. M. Chaikin, Science, 2013, 339, 936–940.

Lauga and Powers 2009

E. Lauga and T. R. Powers, Rep. Prog. Phys., 2009, 72, 096601.

Cates 2012

M. E. Cates, Rep. Prog. Phys., 2012, 75, 042601.

Leptos et al. 2008

K. C. Leptos, J. S. Guasto, J. P. Gollub, A. I. Pesci and R. E. Goldstein, Phy. Rev. E, 2008, 77, 026709.

Jepson et al. 2013

A. Jepson, V. A. Martinez, J. Schwartz-Linek, A. Morozov and W. C. K. Poon, Phys. Rev. E, 2013, 88, 041002.

Semeraro et al. 2018

E. F. Semeraro, J. M. Devos and T. Narayanan, J. Chem. Phys., 2018, 148, 204905.

Kim and Breuer 2004

M. J. Kim and K. S. Breuer, Phys. Fluids, 2004, 16, 78.

Mino et al. 2013

G. L. Mino, J. Dunstan, A. Rousselet, E. Clément and R. Soto, J. Fluid Mech., 2013, 729, 423.

Creppy et al. 2015

A. Creppy, O. Praud, X. Druart, P. L. Kohnke and F. Plouraboué, Phys. Rev. E, 2015, 92, 032722.

Cisneros et al. 2007

L. H. Cisneros, R. Cortez, C. Dombrowski, R. E. Goldstein and J. O. Kessler, Exp. Fluids, 2007, 43, 737.

Wensink et al. 2012

H. H. Wensink, J. Dunkel, S. Heidenreich, K. Drescher, R. E. Goldstein, H. Löwen and J. M. Yeomans, Proc. Natl. Acad. Sci. USA, 2012, 109, 14308–14313.

Dunkel et al. 2013

J. Dunkel, S. Heidenreich, K. Dreschner, H. H. Wensink, M. Bär and R. E. Goldstein, Phys. Rev. Lett., 2013, 110, 228102.

Saintillan and Shelley 2007

D. Saintillan and M. J. Shelley, Phys. Rev. Lett., 2007, 99, 058102.

Pessot et al. 2018

G. Pessot, H. Löwen and A. Menzel, Mol. Phys., 2018, 116, 3401–3408.

Hoell et al. 2018

C. Hoell, H. Löwen and A. M. Menzel, J. Chem. Phys., 2018, 149, 144902.

Yoshinaga and Liverpool 2017

N. Yoshinaga and T. B. Liverpool, Phys. Rev. E, 2017, 96, 020603(R).

Pedley 2016

T. J. Pedley, IMA J. Appl. Math., 2016, 81, 488–521.

Alarcón and Pagonabarraga 2013

F. Alarcón and I. Pagonabarraga, J. Mol. Liq., 2013, 185, 56–61.

Theers et al. 2018

M. Theers, E. Westphal, K. Qi, R. G. Winkler and G. Gompper, Soft Matter, 2018, 14, 8590.

Theers et al. 2017

M. Theers, E. Westphal, G. Gompper and R. G. Winkler, Soft Matter, 2017, 12, 7372.

Stenhammar et al. 2017

J. Stenhammar, C. Nardini, R. W. Nash, D. Marenduzzo and A. Morozov, Phys. Rev. Lett., 2017, 119, 028005.

Yeomans et al. 2014

J. M. Yeomans, D. O. Pushkin and H. Shum, Eur. Phys. J. Special Topics, 2014, 223, 1771–1785.

Saintillan and Shelley 2013

D. Saintillan and M. J. Shelley, C. R. Physique, 2013, 14, 497.

Simha and Ramaswamy 2002

R. A. Simha and S. Ramaswamy, Phys. Rev. Lett., 2002, 89, 058101.

Subramanian and Koch 2009

G. Subramanian and D. L. Koch, J. Fluid. Mech., 2009, 632, 359.

Wolgemuth 2008

C. W. Wolgemuth, Biophys. J., 2008, 95, 1564.

Koch and Subramanian 2011

D. L. Koch and G. Subramanian, Annu. Rev. Fluid Mech., 2011, 43, 637–659.

Doostmohammadi et al. 2018

A. Doostmohammadi, J. Ignés-Mullol, J. Yeomans and F. Sagués, Nat. Commun., 2018, 9, 3246–3246.

Thampi et al. 2013

S. P. Thampi, R. Golestanian and J. M. Yeomans, Phys. Rev. Lett., 2013, 111, 118101.

Ramaswamy 2010

S. Ramaswamy, Annu. Rev. Fluid Mech., 2010, 1, 323–345.

Saintillan and Shelley 2008

D. Saintillan and M. J. Shelley, Phys. Rev. Lett., 2008, 100, 178103.

Hohenegger and Shelley 2010

C. Hohenegger and M. J. Shelley, Phys. Rev. E, 2010, 81, 046311.

Theillard et al. 2017

M. Theillard, R. Alonso-Matilla and D. Saintillan, Soft Matter, 2017, 13, 363–375.

Ezhilan et al. 2013

B. Ezhilan, M. J. Shelley and D. Saintillan, Phys. Fluids, 2013, 25, 070607.

Underhill et al. 2008

P. T. Underhill, J. P. Hernandez-Ortiz and M. D. Graham, Phys. Rev. Lett., 2008, 100, 248101.

Hernandez-Ortiz et al. 2009

J. P. Hernandez-Ortiz, P. T. Underhill and M. D. Graham, J. Phys.: Condens. Matter, 2009, 21, 204107.

Lushi and Peskin 2013

E. Lushi and C. P. Peskin, Comp. Struct., 2013, 122, 239.

Saintillan and Shelley 2012

D. Saintillan and M. J. Shelley, J. Royal Soc. Interface, 2012, 9, 571–585.

Krishnamurthy and Subramanian 2015

D. Krishnamurthy and G. Subramanian, J. Fluid Mech., 2015, 781, 422.

Nash et al. 2008

R. W. Nash, R. Adhikari and M. E. Cates, Phis. Rev. E, 2008, 77, 026709.

Saintillan and Shelley 2008

D. Saintillan and M. J. Shelley, Phys. Fluids, 2008, 20, 123304.

Nash et al. 2010

R. W. Nash, R. Adhikari, J. Tailleur and M. E. Cates, Phys. Rev. Lett., 2010, 104, 258101.

Peskin 2002

C. S. Peskin, Acta Numerica, 2002, 11, 479–517.

de Graaf and Stenhammar 2017

J. de Graaf and J. Stenhammar, Phys. Rev. E, 2017, 95, 023302.

Drescher et al. 2011

K. Drescher, J. Dunkel, L. H. Cisneros, S. Ganguly and R. E. Goldstein, Proc. Natl. Acad. Sci. USA, 2011, 108, 10940.

Cortez et al. 2005

R. Cortez, L. Fauci and A. Medovikov, Phys. Fluids, 2005, 17, 031504.

Zaid et al. 2011

I. M. Zaid, J. Dunkel and J. M. Yeomans, J. Royal Soc. Interface, 2011, 8, 1314–1331.

Böttcher 1973

C. Böttcher, Theory of electric polarization, Elsevier, Netherlands, 1973.

Bibliography51

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1Marchetti et al. 2013 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.
2Bechinger et al. 2006 C. Bechinger, R. D. Leonardo, H. Löwen, C. Reichhardt, G. Volpe and G. Volpe, Rev. Mod. Phys. , 2006, 88 , 045006.
3Bialek et al. 2012 W. Bialek, A. Cavagna, I. Giardina, T. Mora, E. Silvestri, M. Viale and A. M. Walczak, Proc. Natl. Acad. Sci. USA , 2012, 109 , 4786.
4Palacci et al. 2013 J. Palacci, S. Sacanna, A. P. Steinberg, D. J. Pine and P. M. Chaikin, Science , 2013, 339 , 936–940.
5Lauga and Powers 2009 E. Lauga and T. R. Powers, Rep. Prog. Phys. , 2009, 72 , 096601.
6Cates 2012 M. E. Cates, Rep. Prog. Phys. , 2012, 75 , 042601.
7Leptos et al. 2008 K. C. Leptos, J. S. Guasto, J. P. Gollub, A. I. Pesci and R. E. Goldstein, Phy. Rev. E , 2008, 77 , 026709.
8Jepson et al. 2013 A. Jepson, V. A. Martinez, J. Schwartz-Linek, A. Morozov and W. C. K. Poon, Phys. Rev. E , 2013, 88 , 041002.