1. Introduction
Flowing suspensions of non-Brownian particles often display a stochastic, diffusive-like motion known as hydrodynamic diffusion (Davis Reference Davis1996). Common examples include sheared suspensions, where particles driven by an externally imposed velocity gradient undergo chaotic and anisotropic displacements (Drazer et al. Reference Drazer, Koplik, Khusid and Acrivos2002; Pine et al. Reference Pine, Gollub, Brady and Leshansky2005), and sedimentation, where long-range velocity correlations result in tortuous and loopy falling paths (Guazzelli & Hinch Reference Guazzelli and Hinch2011). In both cases, the apparent random motions arise not from thermal fluctuations but rather from flow-induced multiparticle interactions, determined by the relative positions and orientations of many particles, leading to complex dynamics.
The dynamics may be more complex if the suspending particles are active, expending and dissipating energy at the microscale (Ramaswamy Reference Ramaswamy2017). In a seminal experiment, Wu & Libchaber (Reference Wu and Libchaber2000) (WL) studied particle diffusion in a freely suspended soap film containing E. coli bacteria. Adding trace amounts of micrometre-sized polystyrene spheres, they observed short-time superdiffusive and long-time diffusive dynamics of the tracers, with an effective diffusion coefficient 2–3 orders of magnitude higher than that of the background Brownian motion. This enhanced diffusion was attributed to transient formations of swirls and jets, or coherent structures, that increase in size and duration with the bacterial concentration. Many subsequent studies have investigated this mechanism, providing more detailed understanding of particle transport in active suspensions. For example, collective dynamics of microswimmers was reported in a number of experiments (Dombrowski et al. Reference Dombrowski, Cisneros, Chatkaew, Goldstein and Kessler2004; Sokolov et al. Reference Sokolov, Aranson, Kessler and Goldstein2007; Zhang et al. Reference Zhang, Be'er, Florin and Swinney2010; Wensink et al. Reference Wensink, Dunkel, Heidenreich, Drescher, Goldstein, Löwen and Yeomans2012; Lushi, Wioland & Goldstein Reference Lushi, Wioland and Goldstein2014; Li et al. Reference Li, Shi, Huang, Chen, Xiao, Liu, Chaté and Zhang2019) and simulations (Hernandez-Ortiz, Stoltz & Graham Reference Hernandez-Ortiz, Stoltz and Graham2005; Saintillan & Shelley Reference Saintillan and Shelley2007, Reference Saintillan and Shelley2008; Ishikawa & Pedley Reference Ishikawa and Pedley2008; Zöttl & Stark Reference Zöttl and Stark2014), suggesting that large-scale coherent motion could emerge purely from hydrodynamic interactions (Simha & Ramaswamy Reference Simha and Ramaswamy2002; Baskaran & Marchetti Reference Baskaran and Marchetti2009; Koch & Subramanian Reference Koch and Subramanian2011). On the other hand, even in dilute suspensions where swimmers’ dynamics is approximately uncorrelated, there is an enhanced diffusion proportional to the product of the swimmers’ density and their mean speed due to the advective flow they generate (Lin, Thiffeault & Childress Reference Lin, Thiffeault and Childress2011; Miño et al. Reference Miño, Mallouk, Darnige, Hoyos, Dauchet, Dunstan, Soto, Wang, Rousselet and Clement2011; Jepson et al. Reference Jepson, Martinez, Schwarz-Linek, Morozov and Poon2013; Pushkin & Yeomans Reference Pushkin and Yeomans2013; Morozov & Marenduzzo Reference Morozov and Marenduzzo2014). Besides bacteria which swim by pushing the fluid behind, the diffusion of passive particles was also measured in puller-type algal suspensions (Leptos et al. Reference Leptos, Guasto, Gollub, Pesci and Goldstein2009; Kurtuldu et al. Reference Kurtuldu, Guasto, Johnson and Gollub2011; Yang et al. Reference Yang, Peng, Liu, Tang, Xu and Cheng2016; von Rüling, Kolley & Eremin Reference von Rüling, Kolley and Eremin2021). An interesting observation in those systems is that the statistics of the particle displacements, though Gaussian at long times, are non-Gaussian transiently (Thiffeault Reference Thiffeault2015; Ortlieb et al. Reference Ortlieb, Rafaï, Peyla, Wagner and John2019).
Given the importance of hydrodynamic interactions, many theoretical or numerical models have been proposed to elucidate the rich hydrodynamic effects on particle diffusion in active suspensions. Graham and co-workers (Hernandez-Ortiz et al. Reference Hernandez-Ortiz, Stoltz and Graham2005; Underhill, Hernandez-Ortiz & Graham Reference Underhill, Hernandez-Ortiz and Graham2008) developed a minimal model, where self-propelled particles were treated as rigid dumbbells exerting force dipoles on the fluid. The dumbbells were made of point particles and interacted mainly via the fluid flow they generated. Their simulations showed collective dynamics similar to experimental observations and highlighted some differences between ‘pushers’ and ‘pullers’. Saintillan & Shelley (Reference Saintillan and Shelley2007, Reference Saintillan and Shelley2008) used a slender-body model and a kinetic theory to analyse the dynamics of self-propelling rods. Their results showed that, not only were aligned suspensions hydrodynamically unstable as predicted by an earlier theory (Simha & Ramaswamy Reference Simha and Ramaswamy2002), but also there was an instability in isotropic suspensions of pushers due to stress fluctuations. Ishikawa and co-workers (Ishikawa & Pedley Reference Ishikawa and Pedley2007; Ishikawa, Locsei & Pedley Reference Ishikawa, Locsei and Pedley2010; Kogure, Omori & Ishikawa Reference Kogure, Omori and Ishikawa2023) modified the Stokesian dynamics method (Brady & Bossis Reference Brady and Bossis1988) to simulate the diffusion of finite-size or Lagrangian tracer particles in suspensions of ‘squirmers’ (Pedley Reference Pedley2016). Their method provided an accurate description of the near-field hydrodynamic interactions, allowing them to consider more concentrated suspensions. Finally, many other methods have been proposed or adapted to study transport phenomena in active matter, such as the lattice-Boltzmann method (Llopis & Pagonabarraga Reference Llopis and Pagonabarraga2006; Stenhammar et al. Reference Stenhammar, Nardini, Nash, Marenduzzo and Morozov2017), multiparticle collision dynamics (Zöttl & Stark Reference Zöttl and Stark2014), dissipative particle dynamics (Chen et al. Reference Chen, Wei, Sheng and Tsao2016) and force-coupling method (Delmotte et al. Reference Delmotte, Keaveny, Climent and Plouraboué2018).
Despite the significant development and progress so far, a complete understanding of the underlying hydrodynamic interactions between active and passive particles has still not been established. For one, most previous studies were focused on suspensions of self-propelling particles, thus excluding active but individually immotile particles such as melanocytes in the skin (Simha & Ramaswamy Reference Simha and Ramaswamy2002; Baskaran & Marchetti Reference Baskaran and Marchetti2009), microtubule bundles powered by motor proteins (Sanchez et al. Reference Sanchez, Chen, DeCamp, Heymann and Dogic2012; Woodhouse & Goldstein Reference Woodhouse and Goldstein2012) or cells in crowded environments (Hallatschek et al. Reference Hallatschek, Datta, Drescher, Dunkel, Elgeti, Waclaw and Wingreen2023). Without separating these two aspects, one cannot delineate how much the observed collective dynamics is due to self-propulsion versus fluid-mediated hydrodynamic interactions. Furthermore, the passive particles employed earlier were often considered to be tracers that were merely advected by the flow produced by active swimmers. However, at high enough concentrations, the particles may develop certain spatial correlations (Saintillan & Shelley Reference Saintillan and Shelley2007; Li et al. Reference Li, Shi, Huang, Chen, Xiao, Liu, Chaté and Zhang2019), where the dynamics of active particles can be affected by their passive counterpart through both hydrodynamic and excluded-volume interactions. This is particularly relevant for characterizing realistic active systems, biological or synthetic, since they almost always contain some passive components in the form of dead cells or defects (Jepson et al. Reference Jepson, Martinez, Schwarz-Linek, Morozov and Poon2013).
Motivated by the aforementioned observations, we perform large-scale numerical simulations to study the hydrodynamic diffusion in apolar active suspensions. Specifically, we simulate suspensions of shakers, spherical squirmers that are active but not self-propelling, at volume fractions from 0.5 % to 55 %. The numerical method we use, an ‘active version’ of fast Stokesian dynamics (Fiore & Swan Reference Fiore and Swan2019), accounts for both far-field hydrodynamic interactions and near-field lubrication, as well as excluded-volume interactions amongst the particles (Elfring & Brady Reference Elfring and Brady2022). Our approach is similar to that of Ishikawa and Pedley, though the previous work did not consider shaker-type squirmers and the number of particles was often quite small (usually around 100). Here, extensive simulations of more than 1000 shakers show that the instantaneous and long-time translational dynamics are nearly identical between pullers and pushers, which vary non-monotonically with the volume fraction, in contrast to suspensions of self-propelling particles. The rotational dynamics tends to increase with the volume fraction as for self-propelling particles, though the detailed scalings are different. Remarkably, the translational motion of the shakers can be well described by the model of WL, despite the fact that it was originally proposed for passive beads in bacterial suspensions. We explain these results by providing detailed scaling and statistical analyses. Furthermore, we investigate spatial correlations in the suspension microstructure, and find weak local alignment due purely to hydrodynamic interactions.
The paper is organized as follows. In § 2, we describe the mathematical formulation of the squirmer model (§ 2.1) and the governing equations of the active fast Stokesian dynamics method (§ 2.2), followed by the solver verification (§ 2.3) and a brief discussion of pair interactions (§ 2.4). In § 3, we detail the numerical set-up (§ 3.1) and present the full simulation results, beginning with the instantaneous speeds and their distributions (§ 3.2), continuing to the hydrodynamic diffusion in both translational and rotational motion (§ 3.3) and ending with an analysis of spatial correlations in the suspension microstructure (§ 3.4). The discussion is focused on the effect of particle volume fraction, though we have also simulated binary suspensions at varying fractions of passive particles (see Appendix A). Finally, we conclude in § 4.
2. Models and methods
2.1. The squirmer model
The squirmer model, originally proposed by Lighthill (Reference Lighthill1952) and later extended by Blake (Reference Blake1971), has been widely employed to study the swimming mechanisms and collective dynamics of microorganisms at low Reynolds numbers (Pedley Reference Pedley2016). Mathematically, the tangential slip velocity on the surface of a spherical squirmer may be expressed as (Elfring & Brady Reference Elfring and Brady2022)
\begin{equation} {\boldsymbol u}_s (\kern0.7pt {\boldsymbol p} \boldsymbol{\cdot} \hat{{\boldsymbol r}}) = \sum_{n=1}^\infty \frac{2}{n(n+1)} P_n^\prime(\kern0.7pt {\boldsymbol p} \boldsymbol{\cdot} \hat{{\boldsymbol r}}) [B_n({\boldsymbol I}- \hat{{\boldsymbol r}}\hat{{\boldsymbol r}}) \boldsymbol{\cdot} {\boldsymbol p} + C_n {\boldsymbol p} \times \hat{{\boldsymbol r}}], \end{equation}
where 
${\boldsymbol p}$ is a unit vector defining the axis of a squirmer, 
$\hat {{\boldsymbol r}}$ is the unit normal vector on the surface, 
$P_n$ is the Legendre polynomial of degree 
$n$ and 
$P_n^\prime (x)=\textrm {d} P_n(x)/\textrm {d}x$. Note that the slip velocity is axisymmetric as 
${\boldsymbol u}_s$ is a function of 
${\boldsymbol p}\boldsymbol {\cdot } \hat {{\boldsymbol r}}$ only.
From (2.1) we can observe two modes of motion: the 
$B_n$ are the so-called squirming modes, which give rise to the longitudinal velocities about 
${\boldsymbol p}$ (such as the self-propulsion and self-straining); and the 
$C_n$ are associated with azimuthal slip velocities about 
${\boldsymbol p}$ (such as the self-spin). In the literature, it is customary to neglect all spinning modes while only keeping the two leading-order squirming modes (Ishikawa, Simmonds & Pedley Reference Ishikawa, Simmonds and Pedley2006; Elfring & Brady Reference Elfring and Brady2022). Figure 1 shows the local slip velocity due to 
$B_1$ and 
$B_2$, as well as the far-field stresslet flow due to 
$B_2$ (see the supplementary material of Mathijssen, Pushkin & Yeomans (Reference Mathijssen, Pushkin and Yeomans2015) for the far-field flow patterns of a few higher-order moments). Note the geometric asymmetry of the flow intensity in the axial and perpendicular directions due to mass conservation in three dimensions. If the particles are two-dimensional cylinders, however, the flow will be symmetric (i.e. zero velocity along 
$45^{\circ }$) and invariant for 
$\pm B_2$ up to a rotation.
(a) Slip velocity on a squirmer due to either 
$B_1$ (red) or 
$B_2$ (blue). (b) The far-field stresslet flow due to 
$B_2$ (arrows indicate direction and colour the magnitude). In both cases, the squirmer is oriented horizontally to the right.
Finally, as we describe in the next section, the prescribed slip velocity in our simulations is simply 
${\boldsymbol u}_s = {\boldsymbol U}_s + \boldsymbol{\mathsf{E}}_s \boldsymbol {\cdot } \hat {{\boldsymbol r}}$, where 
${\boldsymbol U}_s$ is the self-propulsion velocity and 
$\boldsymbol{\mathsf{E}}_s$ is the self-strain-rate tensor. Terms 
${\boldsymbol U}_s$ and 
$\boldsymbol{\mathsf{E}}_s$ are related to 
$B_1$ and 
$B_2$ as (Elfring & Brady Reference Elfring and Brady2022)
\begin{equation} {\boldsymbol U}_s ={-}\frac{2}{3}B_1{\boldsymbol p},\quad \boldsymbol{\mathsf{E}}_s ={-}\frac{3}{5} \frac{B_2}{a}\left(\kern0.7pt {\boldsymbol p}{\boldsymbol p}-\frac{\boldsymbol I}{3} \right), \end{equation}
where 
$a$ is the particle radius. Unlike the illustration in figure 1(a), this 
${\boldsymbol u}_s$ is not purely tangential. However, since the leading-order hydrodynamic interactions due to activity arise from 
$\boldsymbol{\mathsf{E}}_s$, our description can be considered as a minimal model that captures the hydrodynamic interactions between squirmers.
2.2. Active fast Stokesian dynamics
The dynamics of a suspension of squirmers can be formulated according to the active Stokesian dynamics framework (Elfring & Brady Reference Elfring and Brady2022), extended from the original Stokesian dynamics method (Brady & Bossis Reference Brady and Bossis1988). Assuming small, inertialess particles swimming in a viscous fluid at low Reynolds numbers, the external forces and torques on any particle must be balanced by their hydrodynamic counterparts, which are linearly related to the velocity moments on the particle through a hydrodynamic resistance tensor. In the absence of external flow and Brownian motion, this leads to
\begin{equation} {\boldsymbol U} = {\boldsymbol U}_s + \boldsymbol{\mathsf{R}}_{FU}^{{-}1} \boldsymbol{\cdot} ({\boldsymbol F}_{ext} - \boldsymbol{\mathsf{R}}_{FE}:\boldsymbol{\mathsf{E}}_s), \end{equation}
where 
${\boldsymbol U} \equiv (\boldsymbol U_1\ \boldsymbol U_2 \cdots \boldsymbol U_N)^{\rm T}$ denotes the array of velocities for 
$N$ particles (the same applies to 
$\boldsymbol U_s$, 
$\boldsymbol{\mathsf{E}}_s$ and 
${\boldsymbol F}_{ext}$) and 
$\boldsymbol{\mathsf{R}}_{FU}$ and 
$\boldsymbol{\mathsf{R}}_{FE}$ are the resistance tensors coupling the force and velocity moments of all particles. For example, if particle 
$i$ is active, both 
$\boldsymbol U_{s,i}$ and 
$\boldsymbol{\mathsf{E}}_{s,i}$ can be non-zero (
$\boldsymbol U_{s,i}=\boldsymbol{0}$ for a shaker); otherwise, 
$\boldsymbol U_{s,i}=\boldsymbol{\mathsf{E}}_{s,i} = \boldsymbol{0}$. In this compact notation, 
$\boldsymbol U_i$ includes both the linear and angular velocities of particle 
$i$; similarly, 
$\boldsymbol F_{{ext},i}$ includes both the external force and torque acting on particle 
$i$. Finally, because 
$\boldsymbol{\mathsf{R}}_{FU}$ and 
$\boldsymbol{\mathsf{R}}_{FE}$ only depend on the positions and orientations of the particles, the velocities can be computed quasi-statically to evolve the suspension.
Numerically, we solve a modified form of (2.3) using the fast Stokesian dynamics method (Fiore & Swan Reference Fiore and Swan2019); here, it is considered fast because the computational cost scales linearly with the number of particles. As in the conventional Stokesian dynamics method, the fast Stokesian dynamics method decomposes the hydrodynamic resistance into a far-field, many-body interaction and a near-field, pairwise contribution. Specifically, the far-field interaction is obtained by inverting a truncated multipole expansion of the Stokes flow induced by all particles, whereas the near-field interaction is the lubrication between two particles minus the duplicated parts in the far-field term. This leads to the following expressions of the resistance tensors:
\begin{equation} \boldsymbol{\mathsf{R}}_{FU}= \mathcal{B}^{\rm T} (\boldsymbol{\mathsf{M}}^{\,ff})^{{-}1} \mathcal{B} + \boldsymbol{\mathsf{R}}_{FU}^{nf},\quad \boldsymbol{\mathsf{R}}_{FE}= \mathcal{B}^{\rm T} (\boldsymbol{\mathsf{M}}^{\,ff})^{{-}1} \mathcal{C} + \boldsymbol{\mathsf{R}}_{FE}^{nf}, \end{equation}
where 
$\boldsymbol{\mathsf{R}}^{nf}_{FU}$ and 
$\boldsymbol{\mathsf{R}}^{nf}_{FE}$ are the near-field resistance tensors and 
$\boldsymbol{\mathsf{M}}^{ff}$ is the mobility tensor relating far-field force moments, 
$\mathcal {F}^{ff} \equiv ({\boldsymbol F}^{ff}\ {\boldsymbol S}^{ff})^{\rm T}$, to the velocity moments, 
$\mathcal {U} \equiv ({\boldsymbol U}\ {\boldsymbol E})^{\rm T}$, through mapping tensors 
$\mathcal {B}$ and 
$\mathcal {C}$:
\begin{equation} \mathcal{B} {\boldsymbol U} + \mathcal{C} {\boldsymbol E} = \mathcal{U}, \quad \mathcal{B}^{\rm T} \mathcal{F}^{ff} = {\boldsymbol F}^{ff}. \end{equation}After some algebra, (2.3) can then be reformulated as
\begin{equation} \begin{bmatrix} \boldsymbol{\mathsf{M}}^{\,ff} & \mathcal{B} \\ \mathcal{B}^{\rm T} & - \boldsymbol{\mathsf{R}}_{FU}^{nf} \end{bmatrix}\boldsymbol{\cdot} \begin{bmatrix} \mathcal{F}^{ff} \\ \boldsymbol{U}_* \end{bmatrix}= \begin{bmatrix} -\mathcal{C} \boldsymbol{\mathsf{E}}_s\\ \boldsymbol{F}_* \end{bmatrix}, \end{equation}with the shorthand notation
\begin{equation} \boldsymbol{U}_* = \boldsymbol{U} - \boldsymbol{U}_s,\quad \boldsymbol{F}_* ={-} {\boldsymbol F}_{ext} + \boldsymbol{\mathsf{R}}_{FE}^{nf}: \boldsymbol{\mathsf{E}}_s. \end{equation}This can be verified by substituting (2.4a,b), (2.5a,b), (2.7a,b) into (2.6) to recover (2.3).
Equations (2.6) and (2.7a,b) are the governing equations of our active fast Stokesian dynamics method. They can be solved efficiently using fast iterative methods (i.e. Krylov subspace methods) implemented on modern graphical processing units (Fiore & Swan Reference Fiore and Swan2019). For Krylov subspace methods to converge quickly, it is important to precondition the linear system (Benzi, Golub & Liesen Reference Benzi, Golub and Liesen2005); here, we can use exactly the same preconditioner as in the original fast Stokesian dynamics method since the matrix on the left-hand side of (2.6) does not depend on the particle activity. Once the solution 
${\boldsymbol U}$ is obtained, the particle positions and orientations are integrated forward in time using the standard second-order Runge–Kutta method; see the data availability statement for the source code of our implementation and Fiore & Swan (Reference Fiore and Swan2019) for more details of the numerical method.
2.3. Solver verification
The original fast Stokesian dynamics solver has already been verified in colloidal suspensions by Fiore & Swan (Reference Fiore and Swan2019) and in non-colloidal dense suspensions by Ge & Elfring (Reference Ge and Elfring2022). In the following, we present two additional verifications involving a pair of particles where exact or approximate analytical solutions are available.
First, we consider the relative motion of two identical passive particles in simple shear flow, for which an analytical solution was provided by Batchelor & Green (Reference Batchelor and Green1972). Figure 2 shows the simulated trajectories from different initial positions in comparison with the Batchelor & Green (Reference Batchelor and Green1972) solution. Here, the trajectories are fore–aft symmetric in the absence of non-hydrodynamic interactions or particle roughness. The interparticle gaps are minimal at 
$x=0$, when the two particles are aligned in the velocity gradient direction. The excellent agreement between the simulations and the analytical solutions confirms our numerical calculations of hydrodynamic interactions. Furthermore, the angular velocities of the particles, though not affecting the translational motion of passive spheres, are also verified. This is important for the simulation of active particles as their dynamics depend on the orientations.
(a) Relative trajectories of two passive spheres in the 
$xy$ plane and (b) their angular velocity 
$\varOmega _z$ in simple shear flow, where 
$a$ is the particle radius, 
$\dot \gamma$ the shear rate and 
$\varOmega _z^\infty =\dot \gamma /2$ the angular velocity of the flow. The different curves in (b) correspond to the three initial conditions 
$y_{min}/a=2.001$, 2.01 and 2.1. The analytical solutions (dashed lines) in (a,b) are due to Batchelor & Green (Reference Batchelor and Green1972).
Next, we consider the motion of a passive particle next to a shaker without any background flow. This is the most basic configuration from which the dynamics of a pair of active particles can be constructed, because of the linearity of the Stokes equation (Ishikawa et al. Reference Ishikawa, Simmonds and Pedley2006). In the purely hydrodynamic limit, i.e. no other interactions between the particles, the passive particle is simply advected by the flow generated by the shaker. Furthermore, if the distance between the two particles is relatively large, we may neglect the disturbance to the original flow field caused by the passive particle, because the disturbance flow decays faster with distance. Under this far-field approximation, the velocity of a force-/torque-free passive particle next to a shaker can be obtained by Faxén's law:
\begin{equation} {\boldsymbol U}_{p} = {\boldsymbol u} +\frac{a^2}{6} \nabla^2 {\boldsymbol u}, \quad {\boldsymbol \varOmega}_{p} = \frac{1}{2} \boldsymbol{\nabla} \times {\boldsymbol u}, \end{equation}
where 
${\boldsymbol U}_{p}$ and 
${\boldsymbol \varOmega }_{p}$ are, respectively, the translational and angular velocities of the passive particle and 
${\boldsymbol u}$ denotes the flow field generated by an isolated shaker (Ishikawa et al. Reference Ishikawa, Simmonds and Pedley2006):
\begin{equation} {\boldsymbol u}(\kern0.7pt {\boldsymbol p}, {\boldsymbol r}) = B_2 \left[\left(\frac{a^4}{r^4} - \frac{a^2}{r^2}\right)\left(\frac{3}{2} (\kern0.7pt {\boldsymbol p} \boldsymbol{\cdot} \hat{{\boldsymbol r}})^2 -\frac{1}{2}\right) \hat{{\boldsymbol r}} + \frac{a^4}{r^4} (\kern0.7pt {\boldsymbol p} \boldsymbol{\cdot} \hat{{\boldsymbol r}}) (\kern0.7pt {\boldsymbol p} \boldsymbol{\cdot} \hat{{\boldsymbol r}} \hat{{\boldsymbol r}} - {\boldsymbol p})\right]. \end{equation}
In the above, 
${\boldsymbol r}=r\hat {{\boldsymbol r}}$ is the distance vector from the centre of the shaker to a point in the fluid. We can see that 
$|{\boldsymbol u}| \sim 1/r^2$ to leading order and 
${\boldsymbol u}(\kern0.7pt {\boldsymbol p}, {\boldsymbol r}) = {\boldsymbol u}(-{\boldsymbol p}, {\boldsymbol r})$, both of which are expected for the dipolar flow generated by a shaker. Consequently, the hydrodynamically induced velocities follow the scalings 
$|{\boldsymbol U}_{p}| \sim 1/r^2$ and 
$|{\boldsymbol \varOmega }_{p}| \sim 1/r^3$ to leading order. One thus expects weaker rotations than translations in dilute suspensions of shakers.
Figure 3 shows the velocity of the passive particle as a function of distance (
$r/a$) at various angles (
$\theta$) away from a pulling shaker. Here, the translational velocity is decomposed into a radial (
$U_r$) and an azimuthal (
$U_\theta$) component, and the angular velocity is only non-zero in the perpendicular direction (
$\varOmega _z$). Note that by symmetry only 
$U_r \ne 0$ when 
$\theta =0$ or 
${\rm \pi} /2$, i.e. the passive particle is simply attracted towards or repelled from the active one; however, all three velocity components are non-zero when 
$\theta ={\rm \pi} /4$. Comparing our results with the far-field approximation shows that the two agree reasonably well even at small separations, particularly when 
$\theta ={\rm \pi} /4$, consistent with the simulations of Ishikawa et al. (Reference Ishikawa, Simmonds and Pedley2006). At large separations, we also observe the expected scaling laws estimated above.
Velocity of a passive particle as a function of distance (
$r/a$) next to a pulling shaker, in comparison with the far-field approximation (cf. (2.8a,b) and (2.9)). Here, the passive particle is at angle 
$\theta$ from the shaker, where 
$\cos \theta \equiv {\boldsymbol p}\boldsymbol {\cdot } \hat {{\boldsymbol r}}$; see the main text for details.
2.4. Pair interactions
We can integrate the particle velocities in time to obtain a physical picture of their pair interactions. Figure 4(a) shows the trajectories of a passive particle next to a shaker in a frame co-rotating with the reference shaker (for visualization, see supplementary movies available at https://doi.org/10.1017/jfm.2024.1071). For a puller-type shaker, the region 
$|x|/a \lesssim 1$ cannot be reached from 
$|x|/a \gtrsim 1$ due to excluded-volume interactions (see § 3.1 for details of the repulsive force we impose), implying hysteretic trajectories when reversing the boundary conditions on the shaker. This is confirmed by the trajectories of the passive particle next to a pusher-type shaker (see red dotted lines), which intersect the grey solid lines in the figure. Furthermore, these trajectories suggest that there are no stable configurations for a bound pair: depending on the initial position, the pair may temporarily approach each other but will always separate after the encounter, since all trajectories are open. We note that there are two saddle points located at the poles (
$x=\pm 2a$, 
$y=0$) of a pulling shaker, while there are infinitely many saddle points (
$x=0$) along the equatorial rings of a pushing shaker. The duration of the close encounters around those saddle points depends on the amount of angular noise in the dynamics. In concentrated suspensions, we expect the lifetime of such transient pairs to be short due to the effect of many-body interactions.
Relative trajectories between a reference shaker (pulling or pushing) and a passive particle (a) or another aligned shaker of the same kind (b). Grey solid lines are for pullers and red dotted lines for pushers. In both cases, the trajectory of the second particle is plotted in the frame co-rotating with the reference particle such that the reference particle is oriented horizontally.
Figure 4(b) shows the relative trajectories between two identical, initially aligned shakers (see supplementary movies for visualization). Similar to a shaker–passive pair, the trajectories are hysteretic (i.e. some regions are inaccessible) and all open (i.e. no bound pairs). Under purely hydrodynamic interactions, an initially aligned pair will remain aligned due to the symmetries of the flow induced by a force dipole, because the rotation of the first shaker due to the second one is exactly the same as the rotation of the second shaker due to the first one, which we observe numerically. However, we cannot expect all shakers to be aligned in a suspension because that implies nematic order, which is absolutely unstable at long wavelengths in apolar suspensions (Simha & Ramaswamy Reference Simha and Ramaswamy2002). Therefore, figure 4(b) provides some indication of how a pair of shakers interact, but does not represent a complete phase portrait.
3. Results
3.1. Simulation set-up
We simulate suspensions of shakers at various volume fractions (
$\phi$) ranging from dilute (
$\phi =0.5\,\%$) to dense (
$\phi =55\,\%$). At certain 
$\phi$, we have incorporated passive particles into the suspensions at varying ratios (
$\,f_p$): 
$f_p=0$ represents a fully active suspension, while 
$f_p>0$ represents binary active suspensions (
$\,f_p=1$ represents a fully passive suspension). Since we consider shaker-type squirmers, 
$B_1$ is always zero (thus the particles are apolar and invariant for 
${\pm }{\boldsymbol p}$), but 
$B_2$ can be either positive (pullers) or negative (pushers). Both 
$B_1$ and 
$B_2$ have the dimension of velocity, and at leading order the magnitude of 
$B_2$ determines the rate of the hydrodynamic interactions due to particle activity, as mentioned earlier. Throughout this work we fix 
$|B_2|$ and use it to define a reference time, 
$\tau _a \equiv a/|B_2|$, where 
$a$ is the particle radius (same for all particles in this work). This allows us to define the units in combination of 
$a$ and 
$\tau _a$ for all physical quantities involving length or time; i.e. 
$a/\tau _a$ for translational velocity, 
$1/\tau _a$ for rotational velocity, 
$a^2/\tau _a$ for translational diffusivity, etc. The results in this section are reported in these units unless otherwise stated. Finally, the values of the governing parameters are summarized in table 1.
Summary of the governing parameters.
Numerically, we simulate 
$N=1024$ particles in cubic domains with periodic boundary conditions; see figure 5 and supplementary movies for illustration. We have checked 
$N=2048$ and 4096, but observe no qualitative difference in the statistical results. The initial condition, in terms of the particle position and orientation, is random, and three realizations of each case (in 
$B_2$ and 
$\phi$) are run for 100
$\tau _a$ to reach the average steady state of the short-time dynamics (see § 3.2). For suspension simulations, it is customary to impose a short-range repulsive force to prevent the particles from overlapping. Here, we use the electrostatic repulsion typical for colloids, 
${F}_{ext} = F_0 \exp (-\kappa h)$, where 
${F}_0$ is the maximal repulsion, 
$\kappa$ the inverse Debye length and 
$h$ the surface gap between two particles (Mewis & Wagner Reference Mewis and Wagner2012). This introduces an additional time scale, 
$\tau _{e} \equiv 6{\rm \pi} \eta a^2/F_0$, which should be smaller than 
$\tau _a$ but greater than the numerical time step (
$\eta$ is the fluid dynamic viscosity). As in our previous work (Ge & Elfring Reference Ge and Elfring2022), we choose 
$\tau _{e}/\tau _a=0.02$ and 
$\kappa ^{-1}/a=0.01$ to model a strong and short-range repulsion. The numerical time step (
$\Delta t$) is 100 times smaller than 
$\tau _{e}$, i.e. 
$\Delta t/ \tau _a = 2\times 10^{-4}$.
(a–c) Suspensions of 1024 particles at different 
$\phi$ and 
$f_p$ (darker particles are passive).
3.2. Short-time dynamics: instantaneous speeds
We first examine the short-time suspension dynamics, quantified by the root-mean-square particle speed. Figure 6 shows a few examples of the translational (
$U_{rms}$) and rotational (
$\varOmega _{rms}$) root-mean-square speeds in time for dilute and dense suspensions with or without passive particles. When 
$f_p=0$, both 
$U_{rms}$ and 
$\varOmega _{rms}$ reach steady state in a relatively short time; however, the transients can be longer if the ratio of passive particles is higher or the volume fraction of the suspension is lower. The longer transients are better seen in the evolution of 
$\varOmega _{rms}$, which can also differ for shakers and passive particles. The difference stems from the non-reciprocity of hydrodynamic interactions. In our active suspensions of smooth spheres, the only mechanism to transfer angular momentum is via hydrodynamic interactions. From (2.8a,b), we can infer that the rotational motion of a passive particle due to a nearby shaker is larger than the converse, because the flow induced by a passive particle is of higher order. Therefore, unless the interactions cancel out, passive particles will have a larger 
$\varOmega _{rms}$. On the other hand, the translational momentum can be transferred through both hydrodynamic and excluded-volume interactions, where the latter tend to homogenize the momentum distribution. Consequently, the nearly equal 
$U_{rms}$ among shakers and passive particles suggests the importance of excluded-volume interactions for the short-time dynamics. This is observed in most of the cases that we simulated, except when 
$\phi (1-f_p)$ is small, in which the dynamics is slow because the concentration of active particles is low.
Instantaneous speeds of suspensions of pullers or pushers at 
$\phi =0.05$ or 0.5. In each case, a binary suspension with half passive particles (
$\,f_p=0.5$) is also shown.
The short-time dynamics can be analysed further by examining the distribution of particle speeds. Take suspensions of pulling shakers and passive particles for example (the results for pushers are similar). Figure 7 shows the probability density functions (p.d.f.s) of the translational (
$U$) and rotational (
$\varOmega$) speeds sampled over all particles at steady state. Here, we do not distinguish between active and passive particles because their distributions have nearly the same shape when sampled over a long time. Several interesting observations can be made from these p.d.f.s. First, as found in previous studies (Wu & Libchaber Reference Wu and Libchaber2000; Patteson et al. Reference Patteson, Gopinath, Purohit and Arratia2016), the translational speeds are always well fitted by the Maxwell–Boltzmann (MB) distribution, which describes the behaviour of fluids at equilibrium, even though our suspensions are highly dissipative and far from equilibrium. Second, although the angular speeds tend to be MB when 
$\phi (1-f_p)$ is high, they can also be bimodal in the opposite limit. The bimodal distribution suggests that there are two characteristic speeds in the suspension, yet both must result from hydrodynamic interactions, which only have one characteristic time scale (
$\tau _a$). Third, the relative magnitudes of the characteristic translational and rotational speeds seem to vary non-monotonically with 
$\phi$, though both of them tend to reduce with 
$f_p$ (see Appendix A).
Distribution of particle speeds at different 
$\phi$ and 
$f_p$ for suspensions of pulling shakers and passive particles. Lines are fits to the MB distribution.
We can explain these observations as follows. Consider any test particle in a suspension of shakers. Although no particle can move on its own, the flow induced by each active particle will generate a mean flow that advects all particles suspended in it. If the suspension is homogeneous, the average number of active neighbours per unit length at distance 
$r$ will be 
$n \sim 4{\rm \pi} r^2\phi _a$, where 
$\phi _a \equiv \phi (1-f_p)$. Furthermore, if the particle orientations are isotropic, the typical magnitude of the flow due to all particles at the same distance will scale as 
$\sqrt {n} /r^2 \sim \sqrt {\phi _a}/r$. From (2.8a,b) and (2.9), we thus expect the typical translational and rotational speeds due to particles at distance 
$r$ to scale as 
$\sqrt {\phi _a}/r$ and 
$\sqrt {\phi _a}/r^2$, respectively. As these are statistical estimates, the scalings are only valid at large distances. However, it is clear that the variances of both speeds are finite because the particle size is finite. Consequently, the distributions of the net velocities upon integrating the velocity contributions over all distances tend to be Gaussian per the central limit theorem. Furthermore, since there is no reason for any of the velocity components to differ from one another (otherwise the suspension would not be isotropic), their speeds necessarily display the MB distribution. These conditions are verified in Appendix B. On the other hand, if 
$\phi _a$ is small, the convergence towards Gaussian will be slow and both 
$U$ and 
$\varOmega$ can deviate from the MB distribution. We expect more deviation in 
$\varOmega$ because 
${\boldsymbol \varOmega }_p$ decays faster in 
$r$ than 
${\boldsymbol U}_p$. Hence, the bimodal distribution observed for 
$\varOmega$ is a statistical effect; the additional peak at low speeds in the p.d.f.s corresponds to particles without sufficient active neighbours within a large distance.
Figure 8 shows the root-mean-square speeds of suspensions of pullers or pushers at different volume fractions. The results are generally similar for pullers and pushers, though pullers tend to move slightly faster than pushers. When 
$\phi < 0.2$, both 
$U_{rms}$ and 
$\varOmega _{rms}$ increase with 
$\phi$, but the scaling exponents are less than the expected 
$1/2$ based on the statistical analysis valid at large distances. Therefore, the velocity contributions from nearby particles must have a weaker dependence on 
$\phi$, which we verify in § 3.4. Another observation is that these speeds do not increase monotonically with 
$\phi$. The non-monotonic behaviour of 
$U_{rms}$ is expected from the competition of hydrodynamic and excluded-volume interactions. In dilute suspensions, the mean kinetic energy of the translational motion increases with 
$\phi$ because hydrodynamic interactions increase with 
$\phi$. However, as 
$\phi$ approaches the jamming volume fraction, the dynamics necessarily slows down due to the reduction of free space between particles. If so, the persistence length (
$\ell _p$) of the ‘ballistic motion’ (i.e. a short run along a nearly straight path) may reduce with the volume fraction. The inset in figure 8(a) shows three estimations of 
$\ell _p$, which suggest a scaling of 
$\ell _p \sim \phi ^{-\mu }$, with 
$1/3 \leqslant \mu \leqslant 1/2$. At still higher 
$\phi$, the mean distance between neighbouring particles becomes comparable to two particle radii and 
$\ell _p$ reduces sharply. It is in this regime that 
$U_{rms}$ decreases with 
$\phi$. We defer the definitions of 
$\tau _r$, 
$t_{c}$ and 
$t_{u}$ to the next section (3.3) because they are related to the long-time dynamics; for now, it suffices to note that the different measures give similar results.
Typical instantaneous speeds of the translational (a) and rotational (b) motion in suspensions of pullers or pushers at different volume fractions. The inset in (a) shows the persistence length in suspensions of pullers, while the inset in (b) shows the perpendicular component of 
$\varOmega _{rms}$; see § 3.3 for details. In (a,b), 
$f_p=0$.
Finally, the slight non-monotonic behaviour of 
$\varOmega _{rms}$ is intriguing. From the reasoning above, we would expect 
$\varOmega _{rms}$ to increase monotonically with 
$\phi$, as excluded-volume interactions do not affect the rotation of spheres and the pairwise hydrodynamic interactions are stronger when particles are closer. Therefore, the slight decrease of 
$\varOmega _{rms}$ when 
$\phi > 0.4$ may suggest a partial cancellation of the hydrodynamic effect on particle rotation. Since hydrodynamic interactions are determined by the relative positions and orientations of all particles, this further implies that there may be certain correlations in the suspension microstructure, particularly in the dense regime. We analyse such correlations in detail in § 3.4.
3.3. Long-time dynamics: hydrodynamic diffusion
To characterize the long-time dynamics of our suspensions, we first calculate the mean square displacement (MSD) of the particles, defined as 
$\langle \Delta r^2(t) \rangle \equiv \langle [{\boldsymbol r}(\tau + t) - {\boldsymbol r}(\tau )]^2 \rangle _\tau$, where 
${\boldsymbol r}(t)$ is the position of a particle at time 
$t$ and 
$\langle \cdot \rangle _\tau$ is an average over all particles and reference times 
$\tau$. Figure 9(a) shows the temporal evolution of the MSDs for suspensions of pulling shakers at different volume fractions. The nearly linear growth of the MSDs signifies a diffusive process, which has indeed been observed in various active systems (Wu & Libchaber Reference Wu and Libchaber2000; Leptos et al. Reference Leptos, Guasto, Gollub, Pesci and Goldstein2009; Patteson et al. Reference Patteson, Gopinath, Purohit and Arratia2016; Peng et al. Reference Peng, Lai, Tai, Zhang, Xu and Cheng2016). For bacterial suspensions, WL developed a simple stochastic model to analyse the diffusion of immersed passive beads. Assuming a ‘collisional force’ that is exponentially correlated in time between the beads (i.e. the dynamics of the beads is not random but correlated), it can be shown that the MSD in three dimensions is governed by
\begin{equation} \langle \Delta r^2(t) \rangle = 6D [t - t_{c} (1- \exp({-t/t_{c}}))], \end{equation}
where 
$D$ is the effective diffusion coefficient and 
$t_{c}$ is the lifetime of the coherent structures in the flow. (In the original WL model, the MSD is expressed as 
$\langle \Delta r^2(t) \rangle = 4D t[1- \exp (-t/t_{c})]$, which is in two dimensions and neglects a short-time contribution. The full expression (3.1) fits our data better.) Fitting our data according to (3.1) (using 
$D$ and 
$t_c$ as fitting parameters) collapses all the results, as illustrated in figure 9(b), revealing that the dynamics is ballistic in short times and only becomes diffusive after 
$t>t_{c}$. This observation suggests that the interactions between the shakers may be effectively described by the WL model. However, it does not imply that the collective dynamics of shakers is equivalent to that of self-propelled particles or colloids undergoing so-called persistent Brownian motion (Howse et al. Reference Howse, Jones, Ryan, Gough, Vafabakhsh and Golestanian2007), as the scalings of 
$D$ and 
$\ell _p$ in suspensions of self-propelling versus non-self-propelling particles differ. This is an important distinction and is discussed later.
(a) The MSD in time at different volume fractions for suspensions of pulling shakers at 
$f_p=0$. (b) The same data rescaled according to (3.1), where the inset shows sample particle trajectories at 
$\phi =10\,\%$.
To further verify the WL model against our data, we fitted the velocity autocorrelation, 
$\langle {\boldsymbol U}(\tau ) \boldsymbol {\cdot } {\boldsymbol U}(\tau +t) \rangle _\tau$, according to
\begin{equation} \langle {\boldsymbol U}(\tau) \boldsymbol{\cdot} {\boldsymbol U}(\tau +t) \rangle_\tau = U_{rms}^2 \exp({-t/t_{u}}), \end{equation}
where 
${\boldsymbol U}(\tau )$ is the velocity of a particle at time 
$\tau$, 
$\langle \cdot \rangle _\tau$ is an average over all particles and reference times 
$\tau$ and 
$t_{u}$ is the characteristic time for the velocity decorrelation. The results of the fits are shown in figure 10(a). If the dynamics of our suspensions is correctly described by the WL model, the two time scales, 
$t_{c}$ and 
$t_{u}$, should be equal (Wu & Libchaber Reference Wu and Libchaber2000). Indeed, we observe 
$t_{c} \approx t_{u}$ at nearly all volume fractions except when 
$\phi$ is very small; cf. figure 8(a) inset. When 
$\phi \leqslant 0.01$, the dynamics is too slow for the velocity statistics to converge at the end of the simulations, thus 
$t_{u}$ is less reliable than 
$t_{c}$, as can also be seen in the quality of the fit in figure 10(a). Nevertheless, the scaling between the instantaneous and long-time dynamics is consistent, and we verify that 
$D \approx U_{rms}^2 t_{u}/3$ as expected from the Taylor–Green–Kubo relation (Graham Reference Graham2018) and the exponential decay of the velocity autocorrelation (figure 10b). We have checked that the scalings for pushers are the same. These comparisons thus cross-validate the WL model and our results. In the following, we use (3.1) to extract the diffusion coefficient (
$D$) and crossover time (
$t_c$) of the translational dynamics.
(a) Velocity autocorrelation in time at different volume fractions for suspensions of pulling shakers at 
$f_p=0$. (b) The scaling between 
$D$ and 
$U_{rms}^2 t_{u}$ for all simulated cases.
The long-time rotational dynamics can be quantified by the relaxation of the particle orientations. Specifically, we calculate the autocorrelation function, 
$\langle {\boldsymbol p}(\tau ) \boldsymbol {\cdot } {\boldsymbol p}(\tau +t) \rangle _\tau$, and fit it according to
\begin{equation} \langle \kern0.7pt {\boldsymbol p}(\tau) \boldsymbol{\cdot} {\boldsymbol p}(\tau +t) \rangle_\tau = \exp({-t/\tau_r}), \end{equation}
where 
${\boldsymbol p}(\tau )$ is the orientation of a particle at time 
$\tau$, 
$\langle \cdot \rangle _\tau$ is an average over all particles and reference times 
$\tau$ and 
$\tau _r$ is the characteristic time of the orientation decorrelation. Figure 11(a) shows that 
$\langle {\boldsymbol p}(\tau ) \boldsymbol {\cdot } {\boldsymbol p}(\tau +t) \rangle _\tau$ indeed decays exponentially in time at all volume fractions. For illustration, we also plot the ‘trajectories’ traced by the particle axes, which can be in any direction and change smoothly on a unit sphere. Physically, we expect 
$\tau _r$ to be similar to 
$t_c$, as both are related to the persistence of the ballistic motion. Figure 11(b) shows that, despite the scatter, 
$t_c \propto \tau _r$ and 
$t_c < \tau _r$. Therefore, the translational dynamics ceases to be ballistic before the particle orientations are fully randomized by the rotational diffusion, because the former is also affected by excluded-volume interactions.
(a) Orientation autocorrelation in time at different volume fractions for suspensions of pulling shakers at 
$f_p=0$; inset shows sample orientations traced by the particles on a unit sphere at 
$\phi =55\,\%$. (b) The scaling between 
$t_c$ and 
$\tau _r$ for all cases.
The fitting procedures presented above allow us to examine the long-time translational and rotational dynamics, characterized by 
$D$ and 
$\tau _r$, respectively, at different 
$\phi$ and 
$f_p$. Figure 12 shows that, as for the short-time dynamics, the results for pullers and pushers are similar (with pullers being slightly more diffusive) and 
$D$ is clearly non-monotonic in 
$\phi$. The dependence of 
$D$ on 
$\phi$ resembles that of 
$U_{rms}$ and can be explained by the same underlying mechanism (hydrodynamic versus excluded-volume interactions); however, we note that the peak 
$D$ occurs at a lower 
$\phi$ than the peak 
$U_{rms}$. This is due to the hydrodynamic coupling of the translational and rotational motion. Specifically, because 
$D \sim U_{rms}^2 t_u \sim U_{rms}^2 \tau _r$, its derivative 
${\partial D}/{\partial \phi } \sim D(({2}/{U_{rms}})({\partial U_{rms}}/{\partial \phi }) + ({1}/{\tau _r})({\partial \tau _r}/{\partial \phi }))$; and since 
$\tau _r$ decreases with 
$\phi$, the maximum 
$D$ is reached while 
$U_{rms}$ is still increasing. It is worth noting that somewhat similar non-monotonic diffusivities have also been observed in suspensions of ideally polarizable spheres, which, when subject to an electric field, generate an induced-charge electro-osmotic flow similar to the stresslet flow in figure 1 (Squires & Bazant Reference Squires and Bazant2004; Mirfendereski & Park Reference Mirfendereski and Park2019). The difference is that the induced-charge electro-osmotic flow aligns with the direction of the electric field, whereas the hydrodynamic interactions between shakers render no particular orientation.
Hydrodynamic diffusion. (a) Translational diffusion coefficient for suspensions of pullers or pushers at different volume fractions; inset shows the diffusivity estimated from short-time dynamics. (b) Rotational relaxation time for suspensions of pulling shakers; inset shows the rotational diffusivity 
$d_r$ versus 
$\phi$. In (a,b), 
$f_p=0$.
As for the rotational dynamics, the inverse of 
$\tau _r$ is sometimes used to define a rotational diffusion coefficient, 
$d_r \equiv \tau _r^{-1}$ (Saintillan & Shelley Reference Saintillan and Shelley2007). Figure 12(b) shows that 
$d_r$ increases sublinearly with 
$\phi$ when 
$\phi \lesssim 0.1$. This is roughly consistent with what we expect from the short-time dynamics, as can be verified with 
$\ell _p \sim U_{rms} t_c \sim \phi ^{-\mu }$ and 
$\tau _r \sim t_c$; cf. figures 8(a) and 11(b). We note that in the literature it is sometimes reported that 
$D \sim \phi ^{-1}$ and 
$d_r \sim \phi$ (Ishikawa & Pedley Reference Ishikawa and Pedley2007; Saintillan & Shelley Reference Saintillan and Shelley2007). A key difference between our system and those mentioned lies in the propulsion mechanism. The previous works considered self-propelling particles (
$U_{0} \sim$ const.), thus the persistence length is determined by the mean free path, 
$\ell _{0} \sim \phi ^{-1}$, leading to 
$D \sim U_{0} \ell _{0} \sim \phi ^{-1}$ and 
$d_r \sim \tau _r^{-1} \sim U_{0}/\ell _{0} \sim \phi$. In our case, there is no self-propulsion and the translational motion (
$U_{rms}$) is fully coupled to the particle rotation (
$\tau _r$); hence, the different relationship and scaling, even if the apparent motion of an individual particle may appear similar.
A comparison between figures 8(b) and 12(b) shows that 
$d_r$ does not reduce as in the case of 
$\varOmega _{rms}$ at higher 
$\phi$. One possible reason is that particles may rotate about their own axes, thus leaving their orientation unchanged. However, we have checked the perpendicular component of the angular speed and find it to exhibit the same 
$\phi$ dependence as the full speed; see the inset in figure 8(b). The translational and rotational velocities of the particles are also uncorrelated to their orientations in our suspensions; i.e. 
$\langle {\boldsymbol U} \boldsymbol {\cdot } {\boldsymbol p} \rangle \approx 0$ and 
$\langle {\boldsymbol \varOmega } \boldsymbol {\cdot } {\boldsymbol p} \rangle \approx 0$. We are uncertain about the source of this minor discrepancy.
Finally, we have also examined the long-time dynamics in binary suspensions; i.e. 
$f_p>0$. In general, both 
$D$ and 
$d_r$ tend to reduce with 
$f_p$ at a fixed 
$\phi$ as one would expect; however, the reduction of 
$D$ is nonlinear and, in some cases, it may even increase slightly with 
$f_p$ when 
$f_p \ll 1$. A brief discussion of the 
$f_p$ dependence is provided in Appendix A.
3.4. Microstructure: spatial correlations
From the short-time dynamics we inferred that certain spatial correlations, particularly over short distances, may be present in our suspensions. Therefore, in the reminder of the paper, we examine the near-field microstructure in both particle position and orientation.
We know that long-range nematic order is absolutely unstable in apolar active suspensions (Simha & Ramaswamy Reference Simha and Ramaswamy2002). For shakers, this was shown to result in an isotropic global orientation distribution (Evans et al. Reference Evans, Ishikawa, Yamaguchi and Lauga2011), which we confirm in our simulations; see Appendix B. To examine the local order, we calculate the pair correlation function:
\begin{equation} g({\boldsymbol r}) = \frac{\displaystyle\left\langle \sum_{i \ne j}^N \delta({\boldsymbol r}-{\boldsymbol r}_{ij})\right\rangle}{4{\rm \pi} r^2\rho}, \end{equation}
where 
$\rho$ is the average particle number density and 
$\langle \cdot \rangle$ is an average over all particles and time. Figure 13 shows a few examples of suspensions of pullers or pushers at different 
$\phi$. Here, the radial distribution, 
$g(r)$, displays the typical oscillatory behaviour (especially at higher 
$\phi$) with decaying peaks at 
$r/a \approx 2$, 4, 5.5 and so on, corresponding to the first few particle layers around any particle. However, the angular distributions in the first particle layer are anisotropic. Specifically, pulling shakers tend to find their nearest neighbours along the particle axis (i.e. 
$|{\boldsymbol p}\boldsymbol {\cdot }{\boldsymbol r}|/r \approx 1$), whereas pushing shakers tend to find them in the plane perpendicular to 
${\boldsymbol p}$ (i.e. 
${\boldsymbol p}\boldsymbol {\cdot }{\boldsymbol r}/r \approx 0$), both of which are consistent with what one would expect from their induced stresslet flows; see the insets in figures 13 and 1(b). Furthermore, since the same correlation applies also to the neighbours, adjacent shakers are approximately aligned: for pullers, the alignment is head-to-head; for pushers, it is side-by-side. However, we cannot visually identify any particle chains or sheets as the global microstructure remains disordered (see supplementary movies) and the angular distribution quickly becomes isotropic beyond the first particle layer.
Suspension microstructure. In each panel, the main plot shows the radial distribution function, 
$g(r)$, in suspensions of pullers or pushers at 
$\phi =0.05$ and 0.5, whereas the inset corresponds to the angular distribution at 
$r/a=2$, 4 and 5.5.
To quantify the local alignment, we compute the following correlation functions for the particle orientation (Saintillan & Shelley Reference Saintillan and Shelley2007):
\begin{align}
C_1(r) =
\frac{\displaystyle\left\langle \sum_{i \ne j}^N
(\kern0.7pt {\boldsymbol p}_i \boldsymbol{\cdot} {\boldsymbol p}_j)
\delta(r-r_{ij})\right\rangle}{\displaystyle\left\langle
\sum_{i \ne j}^N \delta(r-r_{ij})\right\rangle},\quad
C_2(r) = \frac{\displaystyle\left\langle \sum_{i \ne j}^N
\dfrac{1}{2}[3(\kern0.7pt {\boldsymbol p}_i \boldsymbol{\cdot}
{\boldsymbol p}_j)^2-1]
\delta(r-r_{ij})\right\rangle}{\displaystyle\left\langle
\sum_{i \ne j}^N \delta(r-r_{ij})\right\rangle},
\end{align}
where 
$\langle \cdot \rangle$ is an average over all particles and time. Here, 
$C_1(r) \in [-1,1]$ measures the polar order, while 
$C_2(r) \in [-\frac {1}{2},1]$ measures the nematic order; both tend to be zero if there is no order. Figure 14 shows that there is indeed a minor local nematic order, despite the absence of polar order at all distances. For pullers, particles are more aligned at lower volume fractions, where the alignment can persist to 
$r/a \approx 4$. For pushers, the alignment is generally weaker and appears to be insensitive to volume fraction. Both of these observations are consistent with what we inferred from the angular distributions (see the insets of figure 13). Comparing with suspensions of self-locomoting rods (Saintillan & Shelley Reference Saintillan and Shelley2007; Li et al. Reference Li, Shi, Huang, Chen, Xiao, Liu, Chaté and Zhang2019), the degree of alignment in our suspensions of spherical shakers is much weaker (in Saintillan & Shelley (Reference Saintillan and Shelley2007), 
$C_2(0) \approx 0.9$). This suggests that particle shape or self-propulsion can have a large effect on the local order.
Correlation functions for the particle orientations at 
$f_p=0$. (a) Polar order, 
$C_1(r)$. (b) Nematic order, 
$C_2(r)$. In both cases, the main plots correspond to suspensions of pullers, while the insets correspond to pushers (same legend).
Finally, we examine the scalings of the neighbour densities, proportional to 
$\rho g(r)$, at small and large distances. Figure 15 shows that, for both pullers and pushers, 
$\rho g(r)$ scales linearly with 
$\phi$ at 
$r/a=10$, consistent with what we expect in an ideal homogeneous suspension; however, the scalings are sublinear at 
$r/a=2$. Therefore, the slower growth of 
$U_{rms}$ and 
$\varOmega _{rms}$ with respect to 
$\phi$ (cf. figure 8) can be attributed to local correlations in the microstructure. Such local correlations may also explain the slight reduction of 
$\varOmega _{rms}$ or the weaker growth of 
$d_r$ when 
$\phi >0.4$, because particles tend to form layers at higher volume fractions and the effect of hydrodynamic interactions from aligned neighbouring particles on opposite sides cancels out. However, since we can only quantify the ordering and alignment statistically, any cancellation is also partial and statistical.
Neighbour densities at small and large distances for suspensions of pullers (a) or pushers (b). In both cases, 
$f_p=0$ and lines are fits of the data for 
$\phi \leqslant 0.2$.
4. Conclusion
In summary, we have presented a numerical study of the hydrodynamic diffusion of apolar active suspensions at volume fractions (
$\phi$) ranging from 0.5 % to 55 %. We model the active particles as squirmers, with zero self-propulsion but non-zero self-straining, and simulate their collective dynamics using a recently developed active fast Stokesian dynamics method. Our results show that both the instantaneous and long-time translational dynamics vary non-monotonically with 
$\phi$ due to the competition between hydrodynamic and excluded-volume interactions. Therefore, the suspension is most diffusive at a certain 
$\phi$, which is between 10 % and 20 % in our case. This is in contrast to suspensions of self-propelling particles, where the translational diffusivity typically reduces with 
$\phi$ due to increased rotational diffusion at higher 
$\phi$. Although the latter is also found in our system, the stronger hydrodynamic coupling between the translational and rotational motion herein changes the overall dynamics qualitatively.
The essential physical mechanism for the above observations lies in: (i) the coupling between the long-time and short-time dynamics, 
$D = U_{rms}^2 t_u/3$; (ii) the proportionality between the translational and rotational relaxation times, 
$t_u \approx t_c \sim \tau _r \equiv d_r^{-1}$; (iii) the scaling of the persistence length of the ballistic motion, 
$\ell _p \sim U_{rms}t_c \sim \phi ^{-\mu }$; and (iv) the scaling of instantaneous translational speed, 
$U_{rms} \sim \phi ^{\nu }$. From these relationships we only need to estimate the values of 
$\mu$ and 
$\nu$ to infer the scalings of 
$D$ and 
$d_r$ with respect to 
$\phi$, at least in the dilute regime. Our data suggest that 
$1/3 \leqslant \mu \leqslant 1/2$, well above the 
$\phi ^{-1}$ scaling for self-propelling particles, possibly due to the curved trajectories of shakers under hydrodynamic interactions (versus the more linear run-and-tumble motion). The estimation of 
$\nu$ is furnished by a statistical analysis valid at large distances, assuming the suspension is homogeneous and isotropic (which we have verified), while also taking into account the local correlations in the microstructure. In more concentrated suspensions, the last two scalings (iii and iv) change substantially because of crowding and more frequent collisions between particles.
Finally, we believe the present study can be extended in a number of ways to explore the dynamics of apolar active suspensions in more complex conditions. One open question concerns the effect of particle shape in dense suspensions. Although long-range nematic order cannot develop spontaneously in any apolar systems, non-spherical particles tend to have a stronger local alignment due to steric interactions at higher concentrations, and such alignment may affect the dynamics when the activities of nearby particles differ. Another interesting question is how the dynamics would be altered by a steady or time-dependent external shear. It is known that non-Brownian particles undergo diffusive motion akin to Taylor dispersion under simple shear, but can also self-organize into non-diffusive ‘absorbing states’ in oscillatory shear (Pine et al. Reference Pine, Gollub, Brady and Leshansky2005; Corté et al. Reference Corté, Chaikin, Gollub and Pine2008). Whether particles endowed with an internal activity can exhibit similar dynamics and, if so, what the ensuing rheology may be (see Ge et al. Reference Ge, Martone, Brandt and Minale2021; Ge & Elfring Reference Ge and Elfring2022) remain to be studied.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2024.1071.
Acknowledgements
The authors wish to thank S. Bagheri, L. Brandt, S.S. Ray, H. Diamant, J.-L. Thiffeault and S. Guo for useful discussions during the early phase of this work, as well as J.F. Brady for a discussion on the scaling while the manuscript was under review.
Funding
This research is supported by the Swedish Research Council via an International Postdoc Grant, no. 2021-06669VR.
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The source code that we used to generate the data of this study is openly available in GitHub at https://github.com/GeZhouyang/FSD.
Appendix A. Instantaneous and long-time dynamics in binary active suspensions
In §§ 3.2 and 3.3 of the main text, we presented the effect of volume fraction (
$\phi$) on the instantaneous and long-time dynamics in purely active suspensions (
$\,f_p =0$). In the following, we briefly describe the results when incorporating passive particles into the suspensions (
$\,f_p >0$).
Figure 16 shows the effect of 
$f_p$ on the instantaneous speeds at a few volume fractions. As discussed in § 3.2, if the suspension is homogeneous and isotropic, we would expect the typical flow generated by neighbouring active particles at a large distance 
$r$ to scale as 
$\sqrt {\phi _a}/r$, where 
$\phi _a \equiv \phi (1-f_p)$, leading to a flow proportional to 
$\sqrt {1-f_p}/r$ at a fixed 
$\phi$. Therefore, if the active and passive particles are randomly mixed at all distances, we would expect both 
$U_{rms}$ and 
$\varOmega _{rms}$ to scale as 
$(1-f_p)^{1/2}$. This is observed in the case of 
$U_{rms}$ but not 
$\varOmega _{rms}$, whose scaling exponent is less than 
$1/2$. Since 
${\boldsymbol \varOmega }_p$ decays faster in 
$r$ than 
${\boldsymbol U}_p$, the different scalings may be due to the local correlations in the microstructure, which have a larger effect on the final 
$\varOmega _{rms}$ than 
$U_{rms}$ (same for the 
$\phi$ dependence). Overall, it is clear that incorporating passive particles reduces the instantaneous particle speeds of the entire suspension.
Instantaneous speeds at different volume fractions (
$\phi$) and ratios of passive particles (
$\,f_p$). (a) Translational speed. (b) Rotational speed.
Figure 17 shows the long-time dynamics of binary suspensions at different volume fractions. Similar to the instantaneous speeds, both 
$D$ and 
$d_r$ tend to reduce with 
$f_p$; however, the reduction of 
$D$ is nonlinear near 
$f_p=0$ and, at 
$\phi =0.5$, it even increases slightly with 
$f_p$. To explain these behaviours, recall that the long-time and short-time dynamics are coupled, 
$D = U_{rms}^2 t_u/3$, and the rotational relaxation times are proportional, 
$t_u \approx t_c \sim \tau _r \equiv d_r^{-1}$. Since 
$U_{rms} \approx U_{rms,0}(1-f_p)^{1/2}$ and 
$t_u \approx t_{u,0}(1 + \alpha f_p)$ for small 
$f_p$, where 
$\alpha >0$ (verified but not shown), we have 
$D \approx U_{rms,0}^2 (1-f_p) t_{u,0}(1 + \alpha f_p)/3 = D_0 (1-f_p) (1+\alpha f_p)$. In these relations, 
$U_{rms,0}$, 
$t_{u,0}$ and 
$D_0$ correspond to 
$U_{rms}$, 
$t_{u}$ and 
$D$ at 
$f_p=0$, respectively. Therefore, if 
$\alpha >1$ (i.e. 
$t_u$ increases rapidly with 
$f_p$), we may observe an enhanced translational diffusion when incorporating passive particles in an active suspension.
Hydrodynamic diffusion at different volume fractions (
$\phi$) and ratios of passive particles (
$\,f_p$). (a) Translational diffusion coefficient. (b) Rotational diffusion coefficient.
Appendix B. Velocity distribution and nematic order
In § 3.2 of the main text, we mentioned that the velocity distribution tends to be Gaussian if there are sufficiently many interacting active particles. Furthermore, the distribution of each velocity components should be the same as we have no reason to expect any difference. Figure 18 verifies these statements. Specifically, we plot the three components of the instantaneous translational and rotational velocities, respectively, for 
$\phi =5\,\%$ or 50 % and 
$f_p=0$ or 0.5. In each case, the different velocity components collapse onto the same distribution, which is always nearly Gaussian except for 
$\varOmega$ when 
$\phi =5\,\%$. When the distribution of 
$\boldsymbol \varOmega$ deviates from Gaussian, its speed distribution also becomes bimodal as shown in figure 7. These trends are generally observed in all cases that we have simulated.
Distribution of particle velocities at different 
$\phi$ and 
$f_p$ for suspensions of pulling shakers and passive particles. Lines are Gaussian fits to the 
$x$ components.
The reason we do not expect the distribution of the velocity components to differ is the absence of nematic order. For a binary suspension of shakers and passive particles, we can calculate a nematic tensor, 
$\boldsymbol{\mathsf{Q}}$, given as
\begin{equation} {\mathsf{Q}}_{ij}=\tfrac{3}{2} \left( \langle p_i p_j \rangle_{s} - \tfrac{1}{3}\delta_{ij} \right), \end{equation}
where 
$\boldsymbol p$ is the orientation of the shaker, 
$\boldsymbol \delta$ the Kronecker delta and 
$\langle \cdot \rangle _s$ denotes the average over all shakers. By definition, 
$\boldsymbol{\mathsf{Q}}$ is symmetric and traceless, and its largest eigenvalue, 
$\lambda$, measures the nematic order: 
$\lambda =0$ if the suspension is isotropic; 
$\lambda =1$ if it is maximally nematic (Doi Reference Doi2013). Figure 19 shows the time evolution of 
$\lambda$ for suspensions of pullers and passive particles at 
$\phi =40\,\%$ and the steady-state values (
$\lambda _m$) at different 
$\phi$ and 
$f_p$. Clearly, there is no nematic order in all cases as 
$\lambda _m$ is always within one standard error, 
$\sigma _\lambda \sim 1/\sqrt {N(1-f_p)}$. Our data thus suggest that shakers remain isotropically oriented at all volume fractions and ratios of passive particles. As a consequence, the individual velocity components must be identically distributed.
Nematic order. (a) Time evolution of the nematic order parameter 
$\lambda$ for pullers and passive particles at 
$\phi =40\,\%$. (b) The mean nematic order 
$\lambda _m$ of all cases.
Open access
