# Supplemental Material for “Enhanced hydrodynamic diffusion in dense binary
active suspensions”

Zhouyang Ge1,2 Shervin Bagheri1 Gwynn J. Elfring2 1FLOW, Department of
Engineering Mechanics, KTH Royal Institute of Technology, 100 44 Stockholm,
Sweden  2Department of Mechanical Engineering and Institute of Applied
Mathematics, University of British Columbia, Vancouver V6T 1Z4, BC, Canada

###### Contents

  1. I Models and methods
    1. I.1 The squirmer model
    2. I.2 Fast Stokesian Dynamics
    3. I.3 Solver verifications
  2. II Results
    1. II.1 Simulation set-up
    2. II.2 Governing parameters
    3. II.3 Spatial and energy distributions
    4. II.4 Translational dynamics
    5. II.5 Global nematic order
    6. II.6 Spatial correlations
    7. II.7 Toy model

##  I Models and methods

###  I.1 The squirmer model

In the main article, we briefly introduced the squirmer model. Here, we
provide some additional details. The tangential slip velocity on the surface
of a spherical squirmer may be expressed as [1]

${\bm{u}}_{s}=\sum_{n=1}^{\infty}\frac{2}{n(n+1)}P_{n}^{\prime}({\bm{p}}\cdot{\bm{r}})\big{[}B_{n}({\bm{I}}-{\bm{r}}{\bm{r}})\cdot{\bm{p}}+C_{n}{\bm{p}}\times{\bm{n}}\big{]},$
(S.1)

where ${\bm{p}}$ is a unit vector along the particle axis, ${\bm{r}}$ the unit
normal vector on the surface, $P_{n}$ the Legendre polynomial of degree $n$,
and $P_{n}^{\prime}(x)=\partial P_{n}(x)/\partial x$. Note that, both
${\bm{u}}_{s}$ and $P_{n}$ are functions of ${\bm{p}}\cdot{\bm{r}}$.

(a)(b)

Figure 1: (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 while color the magnitude). In both cases, the squirmer is
oriented horizontally to the right.

From Eq. (S.1) we can observe two modes of motion: the $B_{n}$’s are the so-
called squirming modes, which gives rise to the longitudinal velocities about
${\bm{p}}$ (such as the self-propulsion and self-straining); and the $C_{n}$’s
are associated with azimuthal slip velocities about ${\bm{p}}$ (such as the
self-spin). As mentioned in the main article, it is customary to neglect all
$C_{n}$’s while only keeping the two leading-order squirming modes. Figure 1
shows the local slip velocity due to $B_{1}$ and $B_{2}$, respectively, as
well as the far-field stresslet flow due to $B_{2}$ (see the Supplementary
Information of Ref. [2] for the far-field flow patterns of a few higher-order
moments). Note the _asymmetry_ of the flow intensity in the axis and
perpendicular directions due to mass conservation in three dimensions (the
flow field is axisymmetric around ${\bm{p}}$); in two dimensions, the axis and
perpendicular directions are symmetric, i.e. the velocity is zero along 45
degrees. Because pullers and pushers only differ in the sign of $B_{2}$ (thus
the velocity direction in Figure 1b) in the squirmer model, the preferential
orientation of two nearby pulling shakers is different from that of two
pushing shakers: the former is more stable if they are aligned in the _axial_
direction (i.e. head-to-head), whereas the latter in the _perpendicular_
directions (i.e. side-by-side). This explains the anisotropic near-field
distributions of pullers and pushers observed in our simulations; c.f. Figure
12 in the supplemental material.

###  I.2 Fast Stokesian Dynamics

The results in the main article were obtained from large-scale hydrodynamic
simulations using the Fast Stokesian Dynamics (FSD) method [3], modified to
include particle activities according to the active Stokesian Dynamics (SD)
framework [1]. Below, we briefly describe the mathematical formulation of
active SD and its implementation in FSD. Full details of the respective
methods are available in Refs. [3, 1].

The dynamics of active suspensions can be formulated as follows. 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 can be 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{pmatrix}{\bm{U}}\\\
{\bm{\Omega}}\end{pmatrix}=\begin{pmatrix}{\bm{U}}_{s}\\\
{\bm{\Omega}_{s}}\end{pmatrix}+{\bm{R}}_{\text{FU}}^{-1}\cdot\begin{bmatrix}\begin{pmatrix}{\bm{F}}_{ext}\\\
{\bm{T}}_{ext}\end{pmatrix}-{\bm{R}}_{\text{FE}}:{\bm{E}}_{s}\end{bmatrix},$
(S.2)

where ${\bm{U}}$ and ${\bm{\Omega}}$ are, respectively, the linear and angular
velocities of the particle, ${\bm{U}}_{s}$, ${\bm{\Omega}}_{s}$ and
${\bm{E}}_{s}$ the respective self-propulsion, self-rotation and self-strain-
rate, ${\bm{F}}_{ext}$ and ${\bm{T}}_{ext}$ the external force and torque
(which may be zero), and ${\bm{R}}_{\text{FU}}$ and ${\bm{R}}_{\text{FE}}$ the
resistance tensors coupling the force and velocity moments. Note that, both
${\bm{R}}_{\text{FU}}$ and ${\bm{R}}_{\text{FE}}$ depend only on the positions
and orientations of the particles, thus the velocities can be computed quasi-
statically to evolve the suspension. Finally, the same relationship also holds
for passive particles; in that case, ${\bm{U}}_{s}$, ${\bm{\Omega}_{s}}$, and
${\bm{E}}_{s}$ are all zero.

Numerically, we solve a modified form of Eq. (S.2) using FSD [3]. As in the
conventional SD [4], FSD decomposes the hydrodynamic resistance into a far-
field and a near-field part,

${\bm{R}}_{\text{FU}}=\mathcal{B}^{T}(\mathcal{M}^{\text{ff}})^{-1}\mathcal{B}+{\bm{R}}_{\text{FU}}^{\text{nf}},\quad{\bm{R}}_{\text{FE}}=\mathcal{B}^{T}(\mathcal{M}^{\text{ff}})^{-1}\mathcal{C}+{\bm{R}}_{\text{FE}}^{\text{nf}},$
(S.3)

where $\mathcal{M}^{\text{ff}}$ is the mobility tensor relating far-field
force moments, $\mathcal{F}^{\text{ff}}\equiv({\bm{F}}^{\text{ff}}\
{\bm{T}}^{\text{ff}}\ {\bm{S}}^{\text{ff}})^{T}$, to the velocity moments,
$\mathcal{U}\equiv({\bm{U}}\ {\bm{\Omega}}\ {\bm{E}})^{T}$, through mapping
tensors $\mathcal{B}$ and $\mathcal{C}$,

$\mathcal{B}\begin{pmatrix}{\bm{U}}\\\
{\bm{\Omega}}\end{pmatrix}+\mathcal{C}{\bm{E}}=\mathcal{U},\quad\mathcal{B}^{T}\mathcal{F}^{\text{ff}}=\begin{pmatrix}{\bm{F}}^{\text{ff}}\\\
{\bm{T}}^{\text{ff}}\end{pmatrix}.$ (S.4)

After some algebra, Eq. (S.2) can then be reformulated as

$\begin{bmatrix}\mathcal{M}^{\text{ff}}&\mathcal{B}\\\
\mathcal{B}^{T}&-\bm{R}_{\text{FU}}^{\text{nf}}\end{bmatrix}\cdot\begin{bmatrix}\mathcal{F}^{\text{ff}}\\\
\bm{U}_{*}\end{bmatrix}=\begin{bmatrix}-\mathcal{C}\bm{E}_{s}\\\
\bm{F}_{*}\end{bmatrix},$ (S.5)

with the short-hand notation

$\bm{U}_{*}=\begin{pmatrix}\bm{U}-\bm{U}_{s}\\\
\bm{\Omega}-\bm{\Omega}_{s}\end{pmatrix},\quad\bm{F}_{*}=-\begin{pmatrix}{\bm{F}}_{ext}\\\
{\bm{T}}_{ext}\end{pmatrix}+\bm{R}_{\text{FE}}^{\text{nf}}:\bm{E}_{s}.$ (S.6)

This can be verified by substituting Eqs. (S.3, S.4, S.6) into Eq. (S.5) to
recover Eq. (S.2).

Eqs. (S.5, S.6) are the governing equations of the active FSD. They can be
solved efficiently using the numerical procedure developed in the original FSD
on a graphical processing unit. Once ${\bm{U}}$ and ${\bm{\Omega}}$ are
obtained, the particle positions and orientations are integrated forward in
time using the standard or stochastic second-order Runge-Kutta method [5]. The
source code of our implementation (including bug fixes of the original FSD
solver) is available at https://github.com/GeZhouyang/FSD.

###  I.3 Solver verifications

The base solver of FSD was developed and verified to some extent by its
original developers, Fiore and Swan [3]. Here, we present two additional
verifications that are relevant to our physical problem. Specifically, we
simulated two cases involving a pair of particles where analytical or
numerical solutions were available previously. For validations of more complex
cases, such as the rheology of dense suspensions, see Ref. [6].

Figure 2: Relative trajectories of two passive spheres in the $xy$ plane
(left) and their angular velocity $\Omega_{z}$ (right) in simple shear flow,
where $a$ is the particle radius, $\dot{\gamma}$ the shear rate, and
$\Omega_{z}^{\infty}=\dot{\gamma}/2$ the angular velocity of the flow. The
different curves in the right figure correspond to the three initial
conditions, $y_{min}/a=2.001$, 2.01, and 2.1. The analytical solutions (dashed
lines) are due to Batchelor and Green [7].

Figure 2 shows the relative trajectories of two passive particles of radius
$a$ in simple shear flow, in comparison to the analytical solutions of
Batchelor and Green [7]. 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
simulation and the analytical solution verifies the numerical calculations of
hydrodynamic interactions (HIs), in both short and long ranges. Furthermore,
the angular velocities of the particles, though not affecting the
translational motion of passive spheres, are also validated. This is a crucial
check of the numerics, as the dynamics of active particles mainly depend on
their orientations.

Figure 3: Velocity of a passive particle as a function of distance, $r/a$,
next to an active particle, in comparison to the analytical and numerical
results of Ishikawa et al. [8]. Here, the passive particle is at 45 degree
from the active particle (left: puller with $\beta=1$; right: shaker with
$\beta=+\infty$); $U_{r}$ and $U_{\theta}$ are the radial and azimuthal
velocity components, while $\Omega_{z}$ is the only non-zero angular velocity;
cf. Figure 5 in Ref. [8].

Figure 3 shows the velocity of a passive particle as a function of distance
away from an active particle, in comparison to the analytical and numerical
results of Ishikawa et al. [8]. Here, the translational velocity is decomposed
into a radial and an azimuthal component, where both of them, as well as the
angular velocity, are normalized by the appropriate squirmer modes ($B_{1}$ or
$B_{2}$). As described in Ref. [8], the analysis of a passive particle next to
an active one is sufficient for obtaining the result of two active particles,
since the latter can be generalized from the former. We compare our simulation
results with the far-field approximations and boundary element simulations (if
available) of Ref. [8] in two cases: the active particle is (i) a puller with
$\beta\equiv B_{2}/B_{1}=1$ and (ii) a shaker with $\beta=+\infty$. In
general, we observe good agreement between our and previous results,
especially in terms of the angular velocity ($\Omega_{z}$) and in the shaker
case. The translational velocities we obtain are closer to the far-field
solutions, which are slightly lower than the results from boundary element
simulations. This deviation is expected since the entrainment effect due to
$B_{1}$ (which decays as $1/r^{3}$ in three dimensions) is neglected in the
current implementation to exploit the computational efficiency of FSD.
However, the qualitative behavior and the scaling of the velocity with respect
to distance are consistent with previous calculations.

##  II Results

###  II.1 Simulation set-up

Case | $\phi$ | $N$ | Group 1 | Group 2 | $\beta_{1}$ | Per
---|---|---|---|---|---|---
1 | 40% | 4096 | S+ | P | $+\infty$ | $+\infty$
2 | 40% | 1024 | S+ | P | $+\infty$ | $+\infty$
3 | 40% | 1024 | S+ | P | $+\infty$ | 0.01
4 | 40% | 1024 | M | P | 0 | 0
5 | 40% | 1024 | S- | P | $-\infty$ | $+\infty$
6 | 40% | 4096 | S- | P | $-\infty$ | $+\infty$
7 | 40% | 1024 | S+ | S- | $+\infty$ | $+\infty$
8 | 50% | 1024 | S+ | S- | $+\infty$ | $+\infty$
9 | 50% | 1024 | S+ | P | $+\infty$ | $+\infty$
10 | 50% | 1024 | S- | P | $-\infty$ | $+\infty$
Table 1: Summary of the simulations. $\phi$ and $N$ are the total volume
fraction and number of particles, respectively. Each suspension contains two
groups of particles at various proportions: pulling shakers (S+), pushing
shakers (S-), neutral movers (M), and passive spheres (P). $\beta_{1}$ and Per
are, respectively, the squirmer parameter and Péclet number based on the
rotational noise of particles in group 1. $\beta_{2}$ can be defined
similarly.

We simulated binary suspensions in various compositions (in terms of the type
of particle) and proportions (in terms of their mixing fraction); see Table 1
for a summary. Although we mostly reported the results of active-passive
mixtures, (S+,P) and (S-,P), at $\phi=40\%$ and $N=4096$ (cases 1 and 6) in
the main article, the full simulations contain results of a smaller sample
size ($N=1024$), active-active mixtures (cases 7 and 8), and $\phi=50\%$
(cases 8–10). In all but two cases (3 and 4) the rotational dynamics of the
active particles were governed by HIs alone; in case 3 we imposed a strong
noise in the angular velocities (see the main article), and in case 4 we
suppressed the HIs due to particle activity entirely. Comparison between cases
1 and 2 as well as 5 and 6 showed no qualitative difference in the statistical
results, thus the number of particles were mostly set to $N=1024$.
Computationally, each simulation was performed in cubic domains with periodic
boundary conditions; different initial conditions (in terms of the particles’
position and orientation, see Figure 4) were also tested and had no effect on
the steady-state results.

Figure 4: Examples of three different initial conditions, where the color
represents the type of particle and the arrow denotes the particle
orientation.

###  II.2 Governing parameters

The control parameter in each case is the fraction of particles in the second
group, termed $f_{p}$ in the main article since we were focusing on the effect
of passive particles. $f_{p}$ was varied between 0 and 1 with smaller
increments near 0. In principle, one could also vary the strength of the
second squirming mode, $B_{2}$, to examine the potential rate dependence of
the global dynamics; c.f. Ref. [9]. However, we fixed $|B_{2}|$ to reduce the
number of parameters. Numerically, it is customary to impose a short-range
repulsive force to prevent the particles from overlapping. Here, we used the
electrostatic repulsion typical for colloids, viz.
${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 [10]. These allow us to define two time scales, $\tau_{a}\equiv
a/|B_{2}|$ and $\tau_{e}\equiv 6\pi\eta a^{2}/F_{0}$. Throughout this work we
fixed $\tau_{e}/\tau_{a}=0.02$ and $\kappa^{-1}/a=0.01$ to model a strong,
short-range repulsion. This choice is consistent with our previous work on
passive suspensions [6].

###  II.3 Spatial and energy distributions

In the main article, we showed the spatial and energy distributions of the
(S+,P) suspension at $f_{p}=0.5$ (case 1), claiming that the former is random
and the latter equally shared among the two particle groups. Here, we provide
more details in support of those claims.

Figure 5 shows the probability $\rho_{k}^{(ij)}$ of a particle in group $i$ to
find its first $k$ nearest neighbors to be all in group $j$, as introduced in
the main article. Specifically, the top row is the (S+,P) suspension with
$f_{p}=0.5$ (case 1), the bottom row the (M,P) suspension with $f_{p}=0.5$
(case 4); the left column shows the temporal evolutions of these
probabilities, the right column their steady-state values. The behaviors for
(S-,P) suspensions as well as those at $\phi=50\%$ are similar to case 1. We
note that there is a slight bias for passive particles (group 2) to be closest
to active ones (more so for pushers than pullers), as $\rho_{k=1}^{(21)}>0.5$.
In the shaker case, this bias may be attributed to the attractive stresslet
flow induced by an active particle in its axial or perpendicular direction
(c.f. Figure 1b), which is absent between passive spheres. In the mover case,
the bias is due to the directed motion – it is more likely for a passive
particle to be hit by (thus close to) an active mover than another passive
particle. However, the main difference between the top and bottom rows of
Figure 5 lies in $\rho_{k}^{(ij)}$ for $k\geq 3$. Clearly, a passive particle
is less likely to be surrounded by a self-propelling mover than what would be
expected in a random configuration, whereas the probability remains random
with increasing $k$ in the shaker case.

Figure 5: Probability $\rho_{k}^{(ij)}$ of a particle in group $i$ to find its
first $k$ nearest neighbors to be all in group $j$; top row: the (S+,P)
suspension with $f_{p}=0.5$ (case 1); bottom row: the (M,P) suspension with
$f_{p}=0.5$ (case 4). Dotted lines correspond to random mixtures.

Figure 6: Root mean square translational ($U_{rms}$) and angular
($\Omega_{rms}$) velocities in time, where $\tau_{a}=a/B_{2}$. Left column:
(S+,P) suspensions (case 1); right column: (S-,P) suspensions (case 6). Top
row: $f_{p}=0.1$; middle row: $f_{p}=0.5$; bottom row: $f_{p}=0.9$. Dotted
lines are the mean.

Figure 6 shows the root mean square translational ($U_{rms}$) and angular
($\Omega_{rms}$) velocities in (S+,P) and (S-,P) suspensions under three
different $f_{p}$’s. As mentioned in the main article, both $U_{rms}$ and
$\Omega_{rms}$ are approximately the same in the two particle groups,
regardless of the composition and mixing fraction. Note that, suspensions of
pulling shakers are more energetic than pushing shakers, consistent with their
long-time diffusivities; c.f. Figure 2(b) in the main article. This is
probably because pulling shakers tend to be aligned along their own axes,
where the induced stresslet flow is stronger than in the perpendicular
directions; then, replacing a shaker with a passive particle will lead to
larger net flow (c.f. the toy model in the main article). Note also that the
results are noisier if they are averaged over less particles (red lines in the
top row and black lines in the bottom row). Nevertheless, their mean values
are roughly the same between the two groups, and we consistently observe
$U_{rms}>\Omega_{rms}a$, indicating the violation of the equipartition
theorem.

###  II.4 Translational dynamics

In the main article, we showed the rescaled mean square displacements (MSD),
$\langle\Delta r^{2}(t)\rangle/6Dt_{c}$, based on the model of Wu and
Libchaber [11]; c.f. Figure 2(a) in the main article. Figure 7(a) shows the
raw MSDs, which crossover from ballistic to diffusive dynamics around time
$t_{c}$. Although the ballistic dynamics persist for more than a few
$\tau_{a}$, the distributions of the displacements all have finite variances,
thus the dynamics will eventually become diffusive according to the central
limit theorem; see Figure 7(a) inset. Figure 7(b) shows the cumulative
probability density function (CDF) of the particle displacements in the $z$
direction over a short time $t=\tau_{a}$ in the ballistic regime (the CDFs in
$x$ and $y$ are similar). These CDFs are generally well fitted by the error
functions, except for the case at $f_{p}=0.99$, where the suspension is mostly
passive. However, even in that case, the CDF, $\mathcal{C}$, quickly converges
to 1 for large displacement, $\Delta s$, as $1-\mathcal{C}(\Delta s)\sim\Delta
s^{-2}$ at least (thus the probability density function, PDF, decays as
$\Delta s^{-3}$ at least); see the inset of Figure 7(b). Indeed, direct
inspection of the PDFs shows that the data are well fitted by either the
Gaussian distribution or an alternative form proposed by Leptos et al. [12],

$\mathcal{P}(\Delta x)=\frac{1-f}{\sqrt{2\pi}\delta_{g}}e^{-\Delta
x^{2}/2\delta_{g}^{2}}+\frac{f}{2\delta_{e}}e^{-|\Delta x|/\delta_{e}},$ (S.7)

where $\delta_{g}$ is the standard deviation of the Gaussian core,
$\delta_{e}$ the exponential decay length, and $f$ the fraction of the
enhanced displacements. The variance is finite in either case. However, we
note that the exponential decay at $f_{p}=0.99$ is probably a transient
behavior due to the poor statistics (not enough active particles in the
sample) [13].

(a)(b)(c)(d)

Figure 7: Statistics of the particle displacements in (S+,P) suspensions (a–c)
and (S-,P) suspensions (d) at $\phi=40\%$, $N=4096$. (a) MSD in time, showing
crossover from ballistic (slope 2) to diffusive (slope 1) dynamics, with
crossover time $t_{c}$ (c.f. Figure 2(a) in the main article). Inset shows the
same data on linear scales. (b) Cumulative probability density function (CDF)
of the particle displacements in the $z$ direction (similarly in $x$ and $y$)
over $t=\tau_{a}$. The full lines are from the simulations and the dashed
lines are their Gaussian fits. Inset shows the decay of $1-\mathcal{C}(\Delta
s)$, where $\mathcal{C}$ denotes the CDF and $\Delta s$ is the scalar net
displacement. (c,d) Probability density function (PDF) of the particle
displacements in $x$ direction over $t=\tau_{a}$, normalized by the standard
deviation of the Gaussian core, $\delta_{g}$; (S+,P) suspensions in (c) and
(S-,P) suspensions in (d). Data at $f_{p}=0.99$ are better fitted by Eq. (S.7)
as in Leptos et al. [12], though the exponential tail is probably a transient
phenomenon [13].

(a)(b)

Figure 8: Relative short-time dynamics (a) and long-time diffusivity (b) of
the binary active suspensions. $\sigma^{2}$ is the MSD per unit time (in
$\tau_{a}$) and $D$ is the diffusion coefficient extracted from the MSD as
described in the main article; their values at $f_{p}=0$ are denoted by the
subscript zero. In both plots, lines are least squares fits of the form
$y=ax^{\nu_{1}}(1-x)^{\nu_{2}}+1-x$. The inset in (b) shows four
representative cases on linear scales.

Figure 9: Correlations between the crossover time $t_{c}$ and the diffusivity
$D$ as well as $U_{rms}$.

Figure 8 shows the translational dynamics in short and long times of the
binary suspensions for all cases. Here, we define the relative short-time
dynamics, $\sigma^{2}/\sigma_{0}^{2}$, as the ratio of the MSDs per unit time
at a finite $f_{p}$ and $f_{p}=0$, and the relative long-time diffusivity,
$D/D_{0}$, same as in the main article. Clearly, the dynamics in both short
and long times can be significantly enhanced in the presence of a small
fraction of passive particles. This applies to both (S+,P) and (S-,P)
suspensions, as well as the (S+,S-) suspensions, as Pe${}_{r}\to\infty$. We
note that the dynamics enhancement is stronger in the short time than in the
long time (nevertheless, they are positively correlated; c.f. Figure 2(c) in
the main article). This implies that there is only one timescale in our
dynamics, which originates from the HIs. Furthermore, the correlation between
$t_{c}$ and $D$ (and also $U_{rms}$) must be negative, since on dimensional
ground one expects $D\sim U_{rms}^{2}t_{c}$. Indeed, we observe $t_{c}\sim
D^{-\nu_{1}}$ and $t_{c}\sim U_{rms}^{-\nu_{2}}$, with $1/3\leq\nu_{1}\leq
2/3$ and $1/2\leq\nu_{1}\leq 1$, see Figure 9. Despite the large scatter of
the data (likely due to uncertainties in estimating $D$ and $t_{c}$), the
lower bounds of these exponents are consistent with the observed $D\sim
U_{rms}^{1.5}$ scaling, c.f. Figure 2(c) in the main article. Overall, the
dynamics enhancement in binary suspensions of shakers is robustly observed,
i.e., it is not sensitive to the volume fraction, sample size, or type of
shakers.

###  II.5 Global nematic order

In the main article, we mentioned that the orientations of the shakers at
steady state are isotropic in our binary suspensions. Below, we provide
explicit evidence supporting this claim.

Since shakers have ${\bf n}\to-{\bf n}$ symmetry, we calculate the nematic
tensor, $\bm{Q}$, defined as

$Q_{ij}=\frac{3}{2}\bigg{(}\langle
n_{i}n_{j}\rangle_{1}-\frac{1}{3}\delta_{ij}\bigg{)},$ (S.8)

where $\langle\cdot\rangle_{1}$ denotes the average over all particles in the
first group. By definition, $\bm{Q}$ is symmetric and traceless, and its
largest eigenvalue, $\lambda$, measures the nematic order: $\lambda=0$ if the
system is isotropic; $\lambda=1$ if it is maximally nematic [14]. Figure 10
shows the time evolution of $\lambda$ and its mean value $\lambda_{m}$ at
steady state. Clearly, there is no nematic order in all cases as $\lambda_{m}$
is always negligible regardless of $f_{p}$, the combination of particles, and
their initial conditions. Note that, the increase of $\lambda_{m}$ as
$f_{p}\to 1$ is due to the finite sample size. By central limit theorem,
$|\lambda_{m}-\bar{\lambda}|\sim\sigma_{\lambda}/\sqrt{N(1-f_{p})}$, where
$\bar{\lambda}$ and $\sigma_{\lambda}$ are the true mean and standard
deviation of $\lambda_{m}$, respectively. If the suspension was isotropic
(i.e. $\bar{\lambda}=0$), we would expect $\lambda_{m}\sim
1/\sqrt{N(1-f_{p})}$, which is indeed the case as shown in Figure 10 (see the
inset of Figure 10a). This is also true for the second particle group, if it
is composed of active particles. Therefore, our data show that shakers remain
isotropically oriented in the presence of passive particles.

S+(a)(b)

Figure 10: Global nematic order. (a) Time evolutions of the nematic order
$\lambda$ of (S+,P) suspensions at $\phi=40\%$ ($N=4096$); inset shows the
orientation distribution (Aitoff projection) of $f_{p}=0$ at steady state. (b)
Mean nematic order $\lambda_{m}$ of all cases. The shade corresponds to the
standard error in the estimation of $\lambda_{m}$ given a limited sample size.

###  II.6 Spatial correlations

Figure 11: Orientation correlations in (S+,P) suspensions at different $f_{p}$
($\phi=40\%$, case 1). Left column: radial distributions of
$\langle|{\bm{n}}_{0}\cdot{\bm{n}}_{r}|\rangle$. Right column: probability
density functions of ${\bm{n}}_{0}\cdot{\bm{n}}_{r}$. Top row: $f_{p}=0$.
Middle row: $f_{p}=0.1$. Bottom row: $f_{p}=0.5$.

Figure 12: Angular distributions at three representative distances for
suspensions of pulling (left) and pushing (right) shakers; $f_{p}=0.5$,
$\phi=40\%$ (cases 1 and 6). Top row: $r/a=2.0$; middle row: $r/a=4.0$; bottom
row: $r/a=6.0$; c.f. Figure 11.

In the main article, we showed an activity-induced local ordering of shakers;
c.f. Figure 4 therein. As can be expected, the ordering is strongest at
$f_{p}=0$ (i.e., when there are only shakers in the suspensions) and will
decrease with $f_{p}$, since the effect of passive particles on the
orientation of other particles are weaker. Figure 11 shows the orientation
correlation, $\langle|{\bm{n}}_{0}\cdot{\bm{n}}_{r}|\rangle$, in (S+,P)
suspensions at $f_{p}=0$, 0.1 and 0.5. Here, the averaging is only over active
particles as the orientation of passive spheres has no effect on the dynamics.
Clearly, the distance over which the orientation of the shakers are correlated
reduces with $f_{p}$. However, nearly touching shakers remain correlated even
at $f_{p}=0.5$ (further increasing $f_{p}$ leads to large fluctuations in the
results as the sample size becomes too small). This can also be seen in Figure
12, where we plot the angular distributions of pulling and pushing shakers in
their own body frame at $f_{p}=0.5$. Note that in these maps the reference
particle is always active, whereas the neighbors can be either active or
passive (since in the toy model we consider the effect of passive particles on
an active lattice). The same patterns are observed at $f_{p}=0.5$ and
$f_{p}=0$ (c.f. Figure 4(a) in the main article), thus the same local order
and the same mechanism for enhanced hydrodynamic diffusion by passive
particles are at play.

###  II.7 Toy model

Figure 13: Mean square velocity in arbitrary units as a function of $f_{p}$
obtained from 100 periodic lattices of 100 particles; the active particles are
aligned in the main plots (black) but oriented randomly in the insets (red).
Left: hexagonal lattices. Right: square lattices with noise in the particle
position (similarly if noisy in the orientation).

In the main article, we illustrated the effect of passive particles on the
dynamics of binary suspensions of shakers using a two dimensional square
lattice, but mentioned that the same effect would be observed in any lattices
with reflection and twofold symmetries, such as hexagonal or cubic lattices,
because they preserve the symmetry of the stresslet flow. As an example,
Figure 13 (left) shows the results on two dimensional hexagonal lattices,
where alignment leads to nonmonotonic dynamics, with a peak at $f_{p}=0.5$,
and random orientations lead to monotonically decreasing dynamics with
$f_{p}$. This is the same as in square lattices (see Figure 5(c) in the main
article). The reason that the dynamics reduce with the amount of passive
particles in randomly oriented lattices is because the induced flow therein is
_additive_ (i.e. it is more likely for the superimposed length of two randomly
oriented vectors to be greater than any of the original one in two and higher
dimensions). In Figure 13 (right) we also show the results on a square
lattice, but with noise added in the particle position (i.e. the particles do
not perfectly sit on the lattice points). Clearly, the dynamics are
nonmonotonic if the particles are aligned, and they resemble the actual
observations better, since $\langle U^{2}\rangle\neq 0$ at $f_{p}=0$ and the
peak occurs at a fraction $f_{*}<0.5$; c.f. Figure 2(b) in the main article.
Finally, we note that the results are essentially the same if a small amount
of noise is added to the particle orientation.

## References

  * Elfring and Brady [2022] G. J. Elfring and J. F. Brady, J. Fluid Mech. 952, A19 (2022).
  * Mathijssen  _et al._ [2015] A. J. Mathijssen, D. O. Pushkin, and J. M. Yeomans, J. Fluid Mech. 773, 498 (2015).
  * Fiore and Swan [2019] A. M. Fiore and J. W. Swan, J. Fluid Mech. 878, 544 (2019).
  * Brady and Bossis [1988] J. F. Brady and G. Bossis, Ann. Rev. Fluid Mech. 20, 111 (1988).
  * Iniesta and Garcia de la Torre [1990] A. Iniesta and J. Garcia de la Torre, J. Chem. Phys. 92, 2015 (1990).
  * Ge and Elfring [2022] Z. Ge and G. J. Elfring, Phys. Rev. E 106, 054616 (2022).
  * Batchelor and Green [1972] G. Batchelor and J.-T. Green, J. Fluid Mech. 56, 375 (1972).
  * Ishikawa  _et al._ [2006] T. Ishikawa, M. Simmonds, and T. J. Pedley, J. Fluid Mech. 568, 119 (2006).
  * Ge  _et al._ [2021] Z. Ge, R. Martone, L. Brandt, and M. Minale, Phys. Rev. Fluids 6, L101301 (2021).
  * Mewis and Wagner [2012] J. Mewis and N. J. Wagner,  _Colloidal suspension rheology_ (Cambridge University Press, 2012).
  * Wu and Libchaber [2000] X.-L. Wu and A. Libchaber, Phys. Rev. Lett. 84, 3017 (2000).
  * Leptos  _et al._ [2009] K. C. Leptos, J. S. Guasto, J. P. Gollub, A. I. Pesci, and R. E. Goldstein, Phys. Rev. Lett. 103, 198103 (2009).
  * Thiffeault [2015] J.-L. Thiffeault, Phys. Rev. E 92, 023023 (2015).
  * Doi [2013] M. Doi,  _Soft matter physics_ (Oxford University Press, USA, 2013).