1. Introduction
Active-matter systems, such as self-propelled colloids and microorganisms, exhibit a wide range of collective behaviours, stemming from their inherently out-of-equilibrium nature (Ramaswamy Reference Ramaswamy2010; Elgeti et al. Reference Elgeti, Winkler and Gompper2015; Gompper Reference Gompper2025). In a liquid medium, self-propelled particles interact via hydrodynamic forces and steric effects. These interactions lead to rich collective dynamics, such as mobility-induced phase separation (MIPS) that produces aggregates (Fily & Marchetti Reference Fily and Marchetti2012; Cates & Tailleur Reference Cates and Tailleur2015), flocks with polar alignment (Ishikawa et al. Reference Ishikawa, Locsei and Pedley2008; Bricard et al. Reference Bricard, Caussin, Desreumaux, Dauchot and Bartolo2013) and chaotic vortex patterns known as bacterial turbulence (Dombrowski et al. Reference Dombrowski, Cisneros, Chatkaew, Goldstein and Kessler2004; Aranson Reference Aranson2022).
The impact of hydrodynamic interactions (HIs) on collective behaviour has been emphasised by comparing active Brownian particles that consider only the excluded-volume (steric) effect (Bechinger et al. Reference Bechinger, Di Leonardo, Löwen, Reichhardt, Volpe and Volpe2016; Zöttl & Stark Reference Zöttl and Stark2023) with squirmers that account for full hydrodynamics (Lighthill Reference Lighthill1952; Blake Reference Blake1971). A squirmer is a self-propelled particle with a surface slip velocity; it is called a ‘puller’ when the anterior thrust is stronger, and a ‘pusher’ when the posterior thrust is stronger (Ishikawa Reference Ishikawa2024). Mobility-induced phase separation has been widely reported in active Brownian particles, whereas studies of squirmers suggest that it is less likely to occur when HIs are fully taken into account (Theers et al. Reference Theers, Westphal, Qi, Winkler and Gompper2018; Zhou & Brady Reference Zhou and Brady2026). In contrast, MIPS-like clustering can emerge for strong puller squirmers (Ishikawa Reference Ishikawa2008; Kyoya et al. Reference Kyoya, Matsunaga, Imai, Omori and Ishikawa2015). This phenomenon suggests that the underlying mechanism differs from standard MIPS driven solely by motility and steric effects; rather, the clustering is mediated by HIs between strong puller squirmers. The effects of HIs require further clarification.
In Stokes flow, where inertia is negligible, the flow field induced by a microswimmer propagates over long distances, giving rise to HIs (Lauga & Powers Reference Lauga and Powers2009; Lauga Reference Lauga2020). Throughout this paper, we distinguish the flow field outside the lubrication region (referred to as the outer region) from that within the lubrication region (the inner region), with a crossover at a surface-to-surface distance of approximately
$0.5{-}1$
particle radii (Brady & Bossis Reference Brady and Bossis1988; Kim & Karrila Reference Kim and Karrila2005). This terminology is adopted to avoid confusion with the conventional use of the term far field in the active-matter literature, where it usually refers to the asymptotic force-dipole flow. The flow field in the outer region corresponds to the outer solution of lubrication theory and can be calculated with significantly greater accuracy than point-stresslet approximations. For spherical particles, this approach yields highly accurate predictions of particle motion down to surface-to-surface separations of the order of
$0.5{-}1$
particle radii. Due to the outer-region HIs, non-spherical swimmers align in the extensional direction, inducing collective swimming (Aditi Simha & Ramaswamy Reference Aditi Simha and Ramaswamy2002; Saintillan & Shelley Reference Saintillan and Shelley2008). When swimmers come into close proximity, HIs in the inner region, mediated by lubrication flows, become significant (Ishikawa et al. Reference Ishikawa, Simmonds and Pedley2006), together with steric effects (Kantsler et al. Reference Kantsler, Dunkel, Polin and Goldstein2013). Considerable debate has therefore persisted regarding whether HIs originating from the outer region or the inner region, or steric effects, are primarily responsible for the observed collective behaviour (Kim & Karrila Reference Kim and Karrila2005; Ishikawa et al. Reference Ishikawa, Locsei and Pedley2008; Saintillan & Shelley Reference Saintillan and Shelley2008; Drescher et al. Reference Drescher, Dunkel, Cisneros, Ganguly and Goldstein2011; Wensink Reference Wensink2012; Redner et al. Reference Redner, Hagan and Baskaran2013; Kyoya et al. Reference Kyoya, Matsunaga, Imai, Omori and Ishikawa2015; Bechinger et al. Reference Bechinger, Di Leonardo, Löwen, Reichhardt, Volpe and Volpe2016; Yoshinaga & Liverpool Reference Yoshinaga and Liverpool2018; Qi et al. Reference Qi, Westphal, Gompper and Winkler2022; Zöttl & Stark Reference Zöttl and Stark2023; Wu-Zhang et al. Reference Wu-Zhang, Fedosov and Gompper2024; Ishikawa Reference Ishikawa2025; Zhou & Brady Reference Zhou and Brady2026). In the Stokes flow regime, however, both contributions are intrinsically coupled and always present simultaneously, making it difficult to isolate their respective roles and draw definitive conclusions.
To date, numerous insights have been accumulated regarding swimming in non-Newtonian fluids (Wu et al. Reference Wu, Chen, Mukasa, Gao and Pak2020; Spagnolie & Underhill Reference Spagnolie and Underhill2023). Collective dynamics in non-plastic, non-Newtonian fluids has also been widely reported (Li & Ardekani Reference Li and Ardekani2016; Tung et al. Reference Tung, Lin, Harvey, Fiore, Ardon, Wu and Suarez2017; Aceves-Sanchez et al. Reference Aceves-Sanchez, Degond, Keaveny, Manhart, Merino-Aceituno and Peurichard2020; Liu et al. Reference Liu, Shankar, Marchetti and Wu2021; Liao & Aranson Reference Liao and Aranson2023; Kobayashi et al. Reference Kobayashi, Molina and Yamamoto2024; Zhang et al. Reference Zhang, Yu, Ouyang, Xia and Lin2024; Zhang et al. Reference Zhang, Yu, Xia, Zhu, Lin and Lu2025; Chamolly & Ishikawa Reference Chamolly and Ishikawa2026). In viscoelastic fluids, for instance, fluid elasticity modulates the local HIs between swimmers and promotes collective swimming (Li & Ardekani Reference Li and Ardekani2016; Tung et al. Reference Tung, Lin, Harvey, Fiore, Ardon, Wu and Suarez2017; Liao & Aranson Reference Liao and Aranson2023).
In contrast to non-plastic non-Newtonian fluids, viscoplastic fluids can restrict their spatial extent due to the presence of a yield stress; examples include mud, slurry, mucus and blood. Consequently, far-field HIs can be effectively eliminated in viscoplastic fluids, whereas they are unavoidable in Stokes flows of Newtonian fluids. This is because a microswimmer induces flow only within the yielded region, and the unyielded regions behave as solid rather than fluid (Balmforth & Hewitt Reference Balmforth and Hewitt2025). This property of viscoplastic fluids is ideal for eliminating only far-field HIs, offering a distinct setting that differs from previous studies. We note that, in a viscoplastic fluid, the stress field does not necessarily vanish outside the yielded region. However, the cutoff of far-field HIs can be understood in a dynamical sense, since HIs between swimmers are conventionally defined as flow-mediated couplings, whereby the velocity disturbance generated by one squirmer induces the translational and rotational motion of another. While the behaviour of a single swimmer (Hewitt & Balmforth Reference Hewitt and Balmforth2017; Eastham et al. Reference Eastham, Mohammadigoushki and Shoele2022; Hewitt & Balmforth Reference Hewitt and Balmforth2022; Nazari et al. Reference Nazari, Shoele and Mohammadigoushki2023; Xia et al. Reference Xia, Yu, Lin and Ishikawa2025c ) and two-swimmer interactions (Xia et al. Reference Xia, Yu, Hu, Zhu and Lin2025a ) in viscoplastic fluids have been examined, collective behaviour in viscoplastic fluids remains unexplored.
In this paper, using full hydrodynamic simulations, we investigate the collective dynamics of a monolayer of squirmers embedded in a viscoplastic medium. Strong pullers form MIPS-like clusters via HIs, but increasing the yield stress and inhibiting far-field HIs hinder their formation. Spherical neutral squirmers exhibit flocking with polar alignment, but increasing the yield stress suppresses its formation. The flocking is governed by the range of HIs and can be collapsed using yield stress and areal fraction of swimmers, underscoring the role of far-field HIs. In contrast, flocking in ellipsoidal neutral squirmers is only weakly affected by yield stress, indicating the dominance of near-field HIs and steric effects. This indicates a qualitatively different mechanism of collective behaviour between spherical and ellipsoidal swimmers. These findings are uniquely enabled by the use of a viscoplastic fluid, providing new insight into the long-standing debate on the relative importance of far- and near-field HIs on collective dynamics.
2. Model description
2.1. Basic equations
To model swimmers, we employ the squirmer model in spherical or ellipsoidal form (figure 1). The parameters
$b_x$
and
$b_z$
denote the semi-minor and semi-major axes, defining the aspect ratio
$a_r=b_z/b_x$
and effective radius
$a=\sqrt {b_xb_z}$
. The modified prolate spheroidal coordinate system that we use is given by (
$\zeta , \tau , \varphi$
) (the details are explained in Appendix A). The surface slip velocity of an ellipsoidal squirmer is given by (Ishikawa et al. Reference Ishikawa, Simmonds and Pedley2006; Theers et al. Reference Theers, Westphal, Qi, Winkler and Gompper2018)
\begin{align} \boldsymbol{u}_{{s}}=-B_1(1+\beta \zeta )\big (\boldsymbol{e}_\zeta \boldsymbol{\cdot }\boldsymbol{p}\big ) \boldsymbol{e}_\zeta , \end{align}
where
$\boldsymbol{p}$
is the unit orientation vector of the squirmer and
$\boldsymbol{e}_\zeta$
is the unit polar vector. The constant
$B_1$
determines the swimming velocity of a solitary squirmer
$v_0=B_1 \tau _0 (\tau _0- (\tau _0^2-1 ) \operatorname {coth}^{-1} \tau _0 )$
in the Stokes limit, with
$\tau _0=(1-a^{-2}_r)^{-1/2}$
(Keller & Wu Reference Keller and Wu1977; Theers et al. Reference Theers, Westphal, Gompper and Winkler2016). For
$b_z=b_x=a$
, the surface velocity for a spherical squirmer is recovered to
$\boldsymbol{u}_{s} = B_1\sin \theta ( {1 + \beta }\cos \theta )\boldsymbol{e}_\theta$
, where
$\theta$
represents the polar angle, while
${\boldsymbol{e}}_\theta$
is the unit tangent vector along the
$\theta$
direction (see figure 1
b). The active stress
$\beta$
helps in categorising the types of squirmer as pusher (
$\beta \lt 0$
), puller (
$\beta \gt 0$
) and neutral swimmer (
$\beta = 0$
). Here we focus specifically on the cases
$\beta = 0$
and
$\beta = \pm 3$
.
Model for (a) an ellipsoidal squirmer and (b) a spherical squirmer.
Figure 1. Long description
Panel A: Diagram of an ellipsoidal squirmer. The diagram shows an ellipsoid with a major axis labeled b_z and a minor axis labeled b_x. The center of the ellipsoid is marked as O. The diagram includes a vector labeled n pointing outward from the surface of the ellipsoid at an angle θ. The vector -e_ζ points downward from the center O along the major axis. The coordinate axes are labeled x’, y’, and z’, with the z’ axis aligned along the major axis of the ellipsoid. Panel B: Diagram of a spherical squirmer. The diagram shows a sphere with radius a. The center of the sphere is marked as O. The diagram includes a vector labeled n pointing outward from the surface of the sphere at an angle θ. The vector e_θ is perpendicular to n and lies in the plane of the sphere. The coordinate axes are labeled x’, y’, and z’, with the z’ axis passing through the center O of the sphere.
The full HIs between swimmers in a viscoplastic fluid are resolved by solving the Navier–Stokes equations using a fictitious domain method (Xia et al. Reference Xia, Yu, Lin and Ishikawa2025c
), while the steric interactions between swimmers are captured using a repulsive force model (Ishikawa et al. Reference Ishikawa, Simmonds and Pedley2006) to avoid the prohibitively small time step needed to overcome the overlapping particles. The basic concept of the fictitious domain method is to fill the interior of the particles with a ‘fictitious’ fluid and to locally apply a forcing term
$\boldsymbol f$
that enforces the boundary conditions on the particle surfaces and accounts for fluid–solid interactions. The flow field is simulated by solving the Navier–Stokes equations in the whole domain including the volume occupied by the particles:
\begin{align} \rho _f\left (\frac {\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}\right )= \boldsymbol{\nabla } \boldsymbol{\cdot }\left (-p \boldsymbol{I}+\boldsymbol{\tau }\right )+\boldsymbol f, \\[-28pt] \nonumber \end{align}
\begin{align} \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{u}=\boldsymbol{0}. \\[0pt] \nonumber \end{align}
Here,
$\boldsymbol u$
,
$p$
and
$\rho _f$
represent the fluid velocity, pressure and density, respectively. We employ a simple Bingham model, for which the shear-stress tensor reads
\begin{align} \begin{aligned} & \boldsymbol{\tau }=\left (2 \mu +\frac {\tau _y}{{\dot \gamma }}\right ) \boldsymbol{\dot \gamma } \quad \text{for}\quad \tau \gt \tau _y, \\ & \boldsymbol{\dot \gamma }=\boldsymbol{0} \quad \text{for}\quad \tau \leqslant \tau _y, \end{aligned} \end{align}
where
$\boldsymbol{\tau }$
denotes the stress tensor,
$\tau _y$
is the yield stress,
$\mu$
is the viscosity constant and
$\boldsymbol{ \dot \gamma }=\boldsymbol{\nabla }\boldsymbol{u}+\boldsymbol{\nabla }\boldsymbol{u}^T$
is the strain-rate tensor. Here
$\tau \equiv \sqrt {({1}/{2}) \boldsymbol \tau : \boldsymbol \tau }$
and
$\dot \gamma \equiv \sqrt {({1}/{2}) \boldsymbol{\dot \gamma }: \boldsymbol{\dot \gamma }}$
refer to the second invariants of the stress and strain-rate tensors, respectively. Swimmer velocity
$\boldsymbol U$
and angular velocity
$\boldsymbol \omega$
obey Newton’s equations of motion:
\begin{align} \rho _s V_p \frac {\mathrm{\,d} \boldsymbol{U}}{\mathrm{\,d} t}=\int _S \boldsymbol{n} \boldsymbol{\cdot }\left (-p \boldsymbol{I}+\boldsymbol{\tau }\right )\, \mathrm{d} s, \\[-28pt] \nonumber \end{align}
\begin{align} \frac {\mathrm{\,d}\left (\boldsymbol{J} \boldsymbol{\cdot }{\boldsymbol{\omega }}\right )}{\mathrm{d} t}=\int _S \boldsymbol{r} \times \left [\boldsymbol{n} \boldsymbol{\cdot }\left (-p \boldsymbol{I}+\boldsymbol{\tau }\right )\right ]\, \mathrm{d} s. \\[0pt] \nonumber \end{align}
Here,
$V_p$
and
$\boldsymbol J$
are the swimmer volume and moment of inertia tensor, respectively,
$\rho _s$
is the density of swimmers and
$\boldsymbol{r}$
is the position vector with respect to the particle centre of mass.
The position
$\boldsymbol{X}$
of each swimmer is determined from the following kinematic equation:
\begin{align} \frac {{{\textrm {d}}{\boldsymbol{X}}}}{{{\textrm {d}}t}} = {\boldsymbol{U}}. \end{align}
The orientation of swimmers is determined by
\begin{align} \frac {{\textrm {d}} \boldsymbol q}{{\textrm {d}}t} = \frac {1}{2}\left (\begin{array}{c@{\quad}c@{\quad}c@{\quad}c} q_4 & -q_3 & q_2 & q_1 \\[3pt] q_3 & q_4 & -q_1 & q_2 \\[3pt] -q_2 & q_1 & q_4 & q_3 \\[3pt] -q_1 & -q_2 & -q_3 & q_4 \end{array}\right )\left (\begin{array}{c} \omega '_x \\[3pt] \omega '_y \\[3pt] \omega '_z \\[3pt] 0 \end{array}\right ) \! . \end{align}
In (2.8),
$q_1$
,
$q_2$
,
$q_3$
and
$q_4$
are the components of the unit quaternion, which can be defined in terms of the Euler angles (Yu & Shao Reference Yu and Shao2007). The components of the angular velocity,
$\omega '_x$
,
$\omega '_y$
and
$\omega '_z$
, are measured in the body-fixed frame. With the quaternion, the rotation matrix from the space-fixed frame to the body-fixed frame reads
\begin{align} \boldsymbol{R}=2\left (\begin{array}{c@{\qquad}c@{\qquad}c} q_1^2+q_4^2-\frac {1}{2} & q_1 q_2+q_3 q_4 & q_1 q_3-q_2 q_4 \\[3pt] q_1 q_2-q_3 q_4 & q_2^2+q_4^2-\frac {1}{2} & q_2 q_3+q_1 q_4 \\[3pt] q_1 q_3+q_2 q_4 & q_2 q_3-q_1 q_4 & q_3^2+q_4^2-\frac {1}{2} \end{array}\right ) \! . \end{align}
Since the swimmer is rotationally symmetric, its orientation
$\boldsymbol{p}$
is fully determined by the direction of its symmetry axis, the
$z'$
axis, in the space-fixed frame. The orientation vector is given by
$\boldsymbol p = \boldsymbol{R} \boldsymbol{\cdot }(0,0,1)^T$
.
The fictitious domain method for modelling the swimming of a spherical or ellipsoidal squirmer in a Newtonian fluid has been validated in our previous studies (Lin & Gao Reference Lin and Gao2019; Xia et al. Reference Xia, Yu, Zhang, Lin and Ouyang2025d ). The accuracy of the method in simulating particle sedimentation (Yu & Wachs Reference Yu and Wachs2007) and swimming dynamics of a squirmer (Xia et al. Reference Xia, Yu, Lin and Ishikawa2025c ) in viscoplastic fluids has also been verified. Further details of the numerical implementation can be found in Appendix A.
2.2. Simulation set-up
The suspension consists of
$N$
swimmers in a computational box with size
$ [L,L,H ]$
. The centres and orientations of all swimmers are restricted to the mid-plane of the box, i.e. a monolayer. We use periodic boundary conditions in the
$x$
and
$y$
directions while the boundary conditions for velocity and pressure on the top and bottom boundaries are
$\partial \boldsymbol{u}/\partial z=\boldsymbol{0}$
and
$\partial p/\partial z=0$
, respectively. The influence of the boundaries in the
$z$
direction is negligible, since the fluid behaves as a solid outside the yielded region (figure 2). The areal fraction of swimmers is given by
$\phi _a = N\pi a^2/L^2$
. The dynamics of our system is governed by the Bingham number,
\begin{align} Bi = \frac {{{\tau _y}{a}}}{{\mu {B_1}}}, \end{align}
which measures the ratio between the yield stress and viscous stress.
In our previous research (Kyoya et al. Reference Kyoya, Matsunaga, Imai, Omori and Ishikawa2015), we investigated a monolayer of spherical or ellipsoidal squirmers in a Newtonian fluid, varying the aspect ratio from 1 to 3. The results show that the behaviour at an aspect ratio of 3 differed markedly from that at spherical shape, whereas the difference at an aspect ratio of 2 was less pronounced, indicating that the influence of aspect ratio becomes significant at higher values. Based on these findings, we focus on an aspect ratio of 3 for ellipsoidal swimmers in the present study.
Schematic representation of the model system. The red iso-surface represents the yield surface defined by
$\tau = \tau _y$
for the representative parameters
$Bi=0.1$
and
$\beta =0$
. (a) Three-dimensional perspective view, (b) cross-section in the
$x{-}y$
plane and (c) cross-section in the
$y{-}z$
plane.
The simulation box of size
$L=32a$
and
$H=16a$
is discretised in a three-dimensional uniform grid with mesh number
$N_x \times N_y \times N_z = 512 \times 512 \times 256$
. For the calculations presented in this work, 32 grid points are used across the particle’s diameter (
$2a/\Delta x=32$
) and the time step is set to be 0.001
$a/B_1$
; each is run for a duration exceeding
$1000a/B_1$
, and various parameters are monitored to ensure reaching statistically steady states. The numerical convergence of the method for swimming in viscoplastic fluids has been validated in our previous study (Xia et al. Reference Xia, Yu, Hu, Zhu and Lin2025a
,
Reference Xia, Yu, Lin and Ishikawac
) through systematic grid refinement. The swimmers were initially placed with random orientations. We choose a small but finite Reynolds number
$Re=\rho _f B_1a/\mu =0.1$
. The effect of the domain width
$L$
(as well as the swimmer number
$N$
) has been examined by using a larger width
$L = 64a$
in Appendix B, which indicates that the domain size does not influence the conclusions from our work.
2.3. Flow disturbance by a solitary squirmer
To establish a baseline of viscoplasticy-mediated HIs, we first examine the flow fields around a solitary puller, pusher or neutral swimmer in both Newtonian and viscoplastic environments. Since a detailed report has already been presented in Xia et al. (Reference Xia, Yu, Lin and Ishikawa2025c
), we provide only a brief overview here. Figures 3 and 4 compare the disturbance flow fields in the Newtonian fluid with those at
$Bi=0.1$
and
$1$
, respectively. As expected, the flow field at
$Bi=0.1$
closely resembles the Newtonian solution, whereas increasing the Bingham number progressively confines the disturbance to a finite yielded region. At
$Bi=1$
, viscoplasticity effectively screens the long-range disturbance generated by the squirmers. The magnitude of the flow inside the yielded region is also altered relative to the Newtonian case. These observations indicate that viscoplasticity does not simply ‘remove’ outer-region HIs. Rather, it accelerates the spatial decay of the disturbance flow.
Magnitude of the disturbance velocity induced by a single squirmer for (a)
$\beta =3$
, (b)
$\beta =-3$
and (c)
$\beta =0$
. In each panel, the upper half shows the Newtonian case (
$Bi=0$
), while the lower half shows the viscoplastic cases (
$Bi=0.1$
). The red contours delineate the yield surface in the viscoplastic fluid.
Magnitude of the disturbance velocity induced by a single squirmer for (a)
$\beta =3$
, (b)
$\beta =-3$
, and (c)
$\beta =0$
. In each panel, the upper half shows the Newtonian case (
$Bi=0$
), while the lower half shows the viscoplastic cases (
$Bi=1$
). The red contours delineate the yield surface in the viscoplastic fluid.
At a small Bingham number, the size of the yielded region can be estimated by balancing the stress induced by the squirmer with the yield stress. For a puller or pusher, we have
\begin{align} \tau (r) \sim \frac {S}{r^3} \sim \tau _y, \end{align}
where
$S \sim \mu v_0 a^2$
is the stresslet. This balance gives a characteristic yielded radius
$R_Y$
as
\begin{align} R_Y \sim \left (\frac {S}{\tau _y}\right )^{1/3}. \end{align}
Then we have
\begin{align} \frac {R_Y}{a} \sim Bi^{-1/3}. \end{align}
For a neutral swimmer, whose leading-order disturbance is a source dipole decaying as
$u \sim v_0(a/r)^3$
and the associated stress scales as
$\tau (r) \sim \mu v_0 a^3/r^4$
. Balancing this stress with the yield stress gives
\begin{align} \frac {R_Y}{a} \sim Bi^{-1/4}. \end{align}
In a homogeneous monolayer suspension, the average separation between swimmers scales as
$d/a = (\pi \phi _a)^{-1/2}$
. Then we expect that the yielded radius compared with the average particle separation scales as
$R_Y/d \sim Bi^{-1/3}\phi ^{1/2}_a$
(pushers/pullers) or
$Bi^{-1/4}\phi ^{1/2}_a$
(neutral swimmers).
The MIPS-like clustering in puller suspensions (
$\beta = 3$
). Snapshots of the system showing collective behaviour in (a–c) a Newtonian fluid and (d–f) a viscoplastic fluid with
$Bi=1$
: (a,d)
$\phi _a=0.25$
,
$a_r=1$
, (b,e)
$\phi _a=0.25$
,
$a_r=3$
and (c,f)
$\phi _a=0.45$
,
$a_r=1$
. All individual squirmers are coloured by the local bond-orientational order
$\psi _{6,k}$
, ranging between 0 (dislocations and disclinations) and 1 (hexagonal crystal). (g–l) Cluster-size distribution functions
$\mathcal{N}(n)$
in (g–i) a Newtonian fluid and (j–l) a viscoplastic fluid.
3. Results
3.1. Clustering of pullers
We begin by considering suspensions of puller squirmers (
$\beta = 3$
). In Newtonian suspensions, strong pullers tend to aggregate into clusters (Ishikawa Reference Ishikawa2008; Kyoya et al. Reference Kyoya, Matsunaga, Imai, Omori and Ishikawa2015). However, in a viscoplastic fluid, the yield stress can dissolve the aggregation (figure 5). To visualise the degree of packing inside clusters, we employ the bond parameter
$\psi _{6,k} \in [0,1]$
that measures local sixfold bond orientational order of swimmer
$k$
(Zöttl & Stark Reference Zöttl and Stark2014), where
$\psi _{6,k}\equiv ({1}/{6}) | \sum _{j \in N_6^k} \exp ({i 6 \vartheta _{k j}} ) |$
. Here the sum goes over the six nearest neighbours of swimmer
$k$
, and
$\vartheta _{kj}$
is the angle between the reference direction and the line connecting the particle
$k$
and its neighbouring particle
$j$
. We show the colour-coded bond parameter in figure 5(a–f). Specifically, the cluster of both spherical and ellipsoidal pullers with
$\phi _a=0.25$
is progressively destroyed by yield stress, where no long-lived clusters exist (see also supplementary movies available at https://doi.org/10.1017/jfm.2026.12036). In a denser suspension (
$\phi _a=0.45$
), the cluster exhibits pronounced hexagonal ordering (figure 5
c). As the Bingham number
$Bi$
increases, the resulting structure becomes progressively weaker and less coherent.
To quantify our findings, we identify the cluster structure using the cluster-size distribution function (Levis & Berthier Reference Levis and Berthier2014; Be’er et al. Reference Be’er, Ilkanaiv, Gross, Kearns, Heidenreich, Bär and Ariel2020; Qi et al. Reference Qi, Westphal, Gompper and Winkler2022). Clusters are identified using a distance criterion: a particle is assigned to a cluster if its nearest neighbour within that cluster satisfies
$\epsilon \leqslant a/4$
, where
$\epsilon$
is the surface-to-surface distance. The normalised cluster-size distribution function is defined as
\begin{align} \mathcal{N}(n)=\frac {\left \langle { n p(n)} \right \rangle }{N}. \end{align}
Here
$\left \langle \boldsymbol{\cdot }\right \rangle$
denotes temporal averaging in the steady regime. The normalised cluster-size distribution function gives the fraction of particles belonging to clusters of size
$n$
, where
$p(n)$
is the number of clusters containing
$n$
particles. In a Newtonian fluid, we observe a bimodal cluster-size distribution with a pronounced peak at large
$n$
, indicating phase separation. At
$Bi=1$
, the second peak of
$\mathcal{N}(n)$
diminishes and vanishes at
$\phi _a=0,25$
(figure 5
g,h), demonstrating the strong suppression of phase separation by yield stress on puller clustering. For denser cases (
$\phi _a = 0.45$
; figure 5
i,l), the second peak weakens with increasing yield stress, but remains visible even at
$Bi=1$
.
In order to quantify puller aggregation, the pair-correlation function is introduced as
\begin{align} g(r_g)=\frac {\left \langle \sum \limits _{i \neq j}^N \delta \left (r_g-r_{\textit{ij}}\right )\right \rangle }{2 \pi n r_g}, \end{align}
where
$r_{\textit{ij}}$
is the centre–centre distance of swimmer pairs and
$n=N/L^2$
is the average number density of swimmers in the plane.
The aggregation
$g(r_g)$
of spherical pullers with increasing Bingham number at
$\phi _a=0.3$
is shown in figure 6. In all cases, the first peak of
$g(r_g)$
locates at around
$r_g \simeq 2a$
, representing the strong tendency of attraction between squirmer pairs. The peak value of
$g(r_g)$
decreases with
$Bi$
, suggesting a weakening of effective attraction due to yield stress.
Pair-correlation function of spherical pullers (
$\beta = 3$
) at
$\phi _a=0.3$
with increasing
$Bi$
.
We introduce the expectation of
$\mathcal{N}(n)$
, i.e.
\begin{align} E_{\mathcal{N}}=\frac {\sum _{n=1}^N n \mathcal{N}(n)}{N} . \end{align}
Physically,
$E_{\mathcal{N}} \in [ {0,1} ]$
provides a measure of the degree of clustering: higher values indicate that particles are more concentrated in larger clusters, while lower values suggest that particles are more uniformly distributed. Here we focus on spherical pullers. Figure 7(a) shows
$E_{\mathcal{N}}$
plotted versus
$\phi _a$
for several values of
$Bi$
. The curves exhibit an increasing sigmoidal shape, where the steep region corresponds to the transition from the gas-like to the cluster state. In a Newtonian fluid, this transition is relatively sharp, whereas it becomes progressively smoother as
$Bi$
increases.
According to figures 7(a) and 5(g–i), we define the cluster state corresponding to
$E_{\mathcal{N}} \geqslant 0.5$
while
$\mathcal{N}(n)$
has a second peak. Figure 7(b) determines a diagram in the parameter space of
$\phi _a$
versus
$Bi$
for pullers. In general, the pullers exhibit a cluster state in the regime with higher values of
$\phi _a$
and lower values of
$Bi$
.
Cluster formation of pullers (
$\beta = 3$
). (a) The value
$E_{\mathcal{N}}$
as a function of the areal fraction at different
$Bi$
. (b) State diagram in the parameter space of
$\phi _a$
versus
$Bi$
. We identify the cluster state (squares) when a second peak appears in the cluster-size distribution function, whereas the gas-like state (circles) corresponds to a single-peak distribution. The pink and blue background regions indicate the clustered and gas-like phases, respectively.
Figure 7. Long description
Panel A: A line graph shows the value of E subscript N as a function of the areal fraction (phi subscript a) at different Bi values. The x-axis represents the areal fraction (phi subscript a) ranging from 0.05 to 0.45, and the y-axis represents E subscript N ranging from 0 to 1. Different colored lines with circle markers represent various Bi values (0, 0.1, 0.3, 0.6, and 1.0). The lines show an increasing trend in E subscript N as phi subscript a increases, with higher Bi values reaching higher E subscript N values more quickly. Panel B: A scatter plot with a background color gradient shows the state diagram in the parameter space of Bi versus phi subscript a. The x-axis represents phi subscript a ranging from 0.1 to 0.5, and the y-axis represents Bi ranging from 0 to 1. The plot uses circles and squares to indicate gas-like and cluster states, respectively. The background color gradient, ranging from blue to red, represents different E subscript N values. The blue region indicates the gas-like phase, while the red region indicates the clustered phase.
We demonstrate that MIPS-like clustering weakens with increasing yield stress as outer-region HIs are cut off, indicating that outer-region HIs play a key role in cluster formation. Correspondingly, the onset of clustering is shifted to higher area fraction (i.e. smaller average particle separation) and becomes more gradual as
$Bi$
increases. This behaviour can be interpreted as a competition between hydrodynamic attraction, which promotes aggregation, and decorrelating effects such as rotational diffusion. When viscoplasticity screens the disturbance flows, stronger crowding is required for hydrodynamic attraction to become dominant. The characteristic cluster size is also likely linked to the extent of the yielded region, since a smaller yielded zone reduces the effective range over which swimmers can interact and thereby suppresses the growth of large aggregates.
We note that the observed changes in collective behaviour should not be interpreted as arising solely from the direct removal of outer-region contributions. Rather, they likely result from a combined effect: the attenuation of long-range disturbances together with modifications of the flow field within the finite yielded region, which alter the local velocity and velocity-gradient fields experienced by nearby swimmers. The comparison with pusher suspensions presented in the following section, where no hydrodynamically mediated phase separation is observed, further highlights the crucial role of strong puller-induced HIs in driving the MIPS-like clustering reported here.
3.2. Whirl structures of pushers
We next examine the role of yield stress in the collective motion of pushers (
$\beta = -3$
). Previous studies in Newtonian fluids have established that elongated pushers tend to develop strong mesoscale spatiotemporal motions. We analyse the velocity correlation among pushers at
$\phi _a=0.3$
in figure 8. Here we use the surface-to-surface distance
$\epsilon$
to emphasise the positive correlation among nearby swimmers. The correlation is calculated only when a swimmer is present lateral to the reference squirmer, i.e.
$\pi /4 {-} 3\pi /4$
from its orientation:
\begin{align} I_V(\epsilon )= \frac { \left \langle \displaystyle \sum _{i\neq j} \boldsymbol U(\boldsymbol{X}_i)\boldsymbol{\cdot }\boldsymbol U(\boldsymbol{X}_j)\, \delta \!\left (\epsilon -\epsilon _{\textit{ij}}\right ) \right \rangle _{\tfrac {\pi }{4}\leqslant \theta _{\textit{ij}}\leqslant \tfrac {3\pi }{4}} }{ \left \langle \displaystyle \sum _{i\neq j} \delta \!\left (\epsilon -\epsilon _{\textit{ij}}\right ) \right \rangle _{\tfrac {\pi }{4}\leqslant \theta _{\textit{ij}}\leqslant \tfrac {3\pi }{4}} }, \end{align}
where
$\theta _{\textit{ij}}=\cos ^{-1} (\boldsymbol p_i \boldsymbol{\cdot }\boldsymbol (\boldsymbol X_j-\boldsymbol X_i)/|\boldsymbol X_j-\boldsymbol X_i|)$
represents the angle between the squirmer orientation and the relative position vector of the squirmer pair.
The profile of
$I_V$
is positive in the small-
$\epsilon$
region but negative in the large
$\epsilon$
region, which defines the whirl structure. Notably, these correlations are significantly more robust in ellipsoidal pushers than in their spherical counterparts. This highlights the importance of steric effects in whirl structures, considering that HIs are comparable for ellipsoidal and spherical squirmers with the same
$\beta$
. As
$Bi$
increases, the magnitudes of the positive correlations decrease for small
$\epsilon$
, indicating that HIs contribute to the formation of whirl structures. However, the persistence of whirl structures even at
$Bi = 1.0$
highlights the importance of steric effects for ellipsoidal swimmers (see also supplementary movies). In elongated pusher suspensions, the extensile stresslet generates flow fields that attract nearby swimmers in the lateral direction. Hydrodynamic attraction brings swimmers into close proximity, where steric interactions arising from their elongated shape reorient them, leading to enhanced alignment and positive velocity correlations. Increasing plasticity is expected to weaken this alignment, as the range of the attractive flow is reduced.
Effects of
$Bi$
on the whirl structures of pushers (
$\beta = -3$
) under
$\phi _a=0.3$
. Spatial correlation of the particle velocity for (a) spherical and (b) ellipsoidal (
$a_r=3$
) pushers.
Effect of
$Bi$
on the polar order of neutral swimmers (
$\beta = 0$
). Temporal evolutions of the global polar order parameter at different
$Bi$
for (a) spherical and (b) ellipsoidal (
$a_r=3$
) squirmers (
$\phi _a=0.3$
). (c) Mean polarisation as a function of
$Bi$
with various
$\phi _a$
; solid lines with circles correspond to spherical squirmers, while the dashed line with squares corresponds to ellipsoidal squirmers. (d) State diagram for spherical squirmers in the parameter space
$\phi _a$
versus
$Bi$
. We define the polar-ordered phase corresponding to
$\mathcal P \geqslant 0.7$
. The blue and pink background regions indicate the disordered and ordered phases, respectively. (e) The second invariant of strain-rate tensor versus the distance
$\epsilon$
relative to the surface of an isolated swimmer along the swimming direction. Inset: evaluated hydrodynamic range
$r_H$
as the distance at which the normalised strain rate becomes less than
$10^{-2}$
; the dashed line corresponds to
$r_H = 0.3Bi^{-1/2}$
.
3.3. Flocking of neutral swimmers
Finally, we consider suspensions of neutral squirmers (
$\beta =0$
), which are known to exhibit flocks with polar alignment in Newtonian suspensions (Ishikawa et al. Reference Ishikawa, Locsei and Pedley2008; Bricard et al. Reference Bricard, Caussin, Desreumaux, Dauchot and Bartolo2013). The strength of ordering can be described by the global polar order parameter
\begin{align} P(t) = \frac {\left | {\sum \nolimits _i^N {{\boldsymbol p_i}} } \right |}{N}. \end{align}
Quantitatively,
$P$
exhibits distinct variation curves depending on
$Bi$
(figure 9
a,b;
$\phi _a=0.3$
). For spherical swimmers, when
$Bi \leqslant 0.6$
, the system evolves towards a steady polar-ordered state (
$P \to 1$
). This ordering is purely hydrodynamic, as spheres possess no steric reorientation mechanism. Once
$Bi$
exceeds a critical threshold (
$Bi_c \approx 0.65$
), the system undergoes a transition from global order to global disorder (
$P \ll 1$
) (figure 9
a). Even for initially aligned configurations at
$Bi=0.7$
, the global polar order decays from perfect alignment, demonstrating that the ordered state is unstable at sufficiently large
$Bi$
. These results suggest that HIs play a crucial role in the emergence of polar order. These observations are also qualitatively consistent with our previous study of pairwise neutral squirmer interactions in viscoplastic fluids (Xia et al. Reference Xia, Yu, Hu, Zhu and Lin2025a
). In a Newtonian fluid, HIs between two neutral squirmers promote reorientation towards a more aligned configuration following close interaction, whereas increasing the Bingham number weakens this alignment-promoting reorientation. The reduced tendency for pairwise alignment provides a plausible microscopic mechanism for the weaker orientational ordering observed in the present collective simulations.
In contrast, ellipsoidal swimmers (
$a_r=3$
) exhibit pronounced fluctuations in
$P(t)$
(figure 9
b) consistent with previous observations (Kyoya et al. Reference Kyoya, Matsunaga, Imai, Omori and Ishikawa2015). This difference from spherical swimmers is mainly attributed to steric effects, suggesting that the observed fluctuations originate from steric interactions. Although they do not develop as strong an ordered state as spherical swimmers in a Newtonian medium, their ordering is markedly less sensitive to
$Bi$
. This reduced sensitivity stems from the steric effect to align with neighbouring ones, which dominates over HIs (Baskaran & Marchetti Reference Baskaran and Marchetti2008). We further examine the combined effects of
$Bi$
and
$\phi _a$
on the long-time mean polarisation
$\mathcal P= \langle P(t)\rangle$
. Figure 9(c) shows
$\mathcal P$
as a function of
$Bi$
for
$\phi _a=0.15$
to
$0.6$
. For spherical swimmers,
$\mathcal P$
drops sharply when
$Bi$
exceeds a critical value, signalling a clear transition from ordered to disordered states. In contrast, for ellipsoidal swimmers,
$\mathcal P$
decreases gradually with
$Bi$
, without exhibiting a well-defined transition, indicating a qualitative difference from spherical swimmers.
Our analysis produces a phase diagram (figure 9
d) in the
$(\phi _a, Bi)$
space that delineates the transition between ordered and disordered states for spherical swimmers. The red region indicates the presence of ordering among the swimmers, whereas in the blue region this ordering nearly vanishes. This polar order vanishes at high areal fraction. For
$\phi _a$
smaller than 0.65, a well-defined critical Bingham number
$Bi_c$
can be identified for each
$\phi _a$
. Interestingly,
$Bi_c$
exhibits a non-monotonic dependence on
$\phi _a$
. The value of
$Bi_c$
reaches its maximum at
$\phi _a \approx 0.3$
. In the semi-dilute regime (
$\phi _a \leqslant 0.3$
),
$Bi_c$
increases approximately linearly with
$\phi _a$
, whereas at higher areal fractions (
$\phi _a \geqslant 0.45$
),
$Bi_c$
drops sharply as
$\phi _a$
increases.
We now rationalise the observed non-monotonic dependence of
$Bi_c$
on
$\phi _a$
. Our results indicate that outer-region HIs are key in the development of polar order. In the semi-dilute regime, we assume that polar alignment is governed by HIs acting within an effective HI range
$r_H$
, over which HI-induced reorientation, governed by the local strain-rate tensor, remains significant. As shown in figure 9(e),
$\dot \gamma$
decays rapidly along the swimming axis, and its magnitude is strongly suppressed with increasing
$Bi$
, indicating a reduction of
$r_H$
. Defining
$r_H$
as the distance at which the normalised strain rate becomes negligible (
$(a/B_1)\dot \gamma \lt 10^{-2}$
), we extract
$r_H(Bi)$
and find an empirical scaling
$r_H\propto Bi^{-1/2}$
(inset of figure 9
e). Polar order can develop only when the typical nearest-neighbour distance
$r_{nn}\sim \phi _a^{-1/2}$
is smaller than
$r_H$
. Equating
$r_H=r_{nn}$
leads to the scaling
$Bi_c\sim \phi _a$
, consistent with the simulation results in figure 9(d), with
$Bi_c\simeq 2\phi _a$
for
$\phi _a\leqslant 0.3$
. In contrast, as
$\phi _a$
approaches the jamming areal fraction, the fluid throughout the middle plane is fully yielded, and the available free space between swimmers is strongly reduced, resulting in enhanced near-field HIs via lubrication forces. Because lubrication forces disrupt polar alignment (Ishikawa Reference Ishikawa2019), flocks become destabilised, leading to a decrease of
$Bi_c$
at high
$\phi _a$
.
4. Conclusion
We investigated the collective dynamics of microswimmers embedded in a viscoplastic medium. This system constitutes a distinct class in which HIs can be effectively cut off, shedding new light on the long-standing debate regarding the respective roles of far- and near-field HIs in collective behaviour. Specifically, our results highlight the importance of outer-region HIs for spherical squirmers in forming MIPS-like clusters, whirls and flocking. In contrast, steric interactions play a crucial role for ellipsoidal squirmers.
The mechanisms discussed in the present work should also be distinguished from the long-wavelength hydrodynamic instability of infinitely extended suspensions of straight swimmers described by continuum kinetic theories (Saintillan & Shelley Reference Saintillan and Shelley2007, Reference Saintillan and Shelley2008). Here, we focus on a confined microswimmer monolayer and investigate how viscoplasticity modifies the effective HI range governing the observed collective dynamics. The findings of this study provide general insight into the use of viscoplastic fluids as a means to tune collective behaviour in active suspensions.
In the absence of a yield stress, a strongly shear-thinning fluid (i.e. pseudo-plastic fluid) may also lead to a significant reduction of long-range HIs, and some of the observed features could potentially arise from this effect. A systematic comparison with shear-thinning fluids would be an interesting direction for future work.
Another important direction for future study is to extend the present monolayer configuration to fully three-dimensional suspensions. The HIs in a quasi-two-dimensional monolayer can differ fundamentally from those in an unconfined three-dimensional bulk system, where out-of-plane motion and three-dimensional structures may substantially modify collective motions. Investigating how viscoplasticity alters these genuinely three-dimensional hydrodynamic couplings would help clarify the universality of the mechanisms identified here and their relevance to active suspensions.
Supplementary movies
Supplementary movies are available at https://doi.org/10.1017/jfm.2026.12036.
Funding
The work was supported by the National Natural Science Foundation of China (grant nos. 12302340, 12132015 and 12332015), the JSPS KAKENHI (grant nos. 24KF0097, 21H04999 and 21H05308) and the JST CREST (grant no. JPMJCR25Q1).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Numerical implementation
A.1. The squirmer model
The swimmer considered has a prolate spheroidal shape. In a body-fixed Cartesian coordinate, the equation describing the swimmer’s surface reads
\begin{align} \frac {{{x^2} + {y^2}}}{{b_x^2}} + \frac {{{z^2}}}{{b_z^2}} = 1, \end{align}
where
$b_x$
and
$b_z$
are the semi-minor and semi-major axes, respectively. We denote half of the focal length by
$c = \sqrt {b_z^2 - b_x^2}$
, which yields the eccentricity
$e=c/b_z$
(
$0 \leqslant e \lt 1$
, with
$e = 0$
describing a sphere). The prolate spheroidal coordinates (
$\zeta , \tau , \varphi$
) are related to the Cartesian coordinates (
$x,y,z$
) via
\begin{align} \zeta = \frac {1}{{2c}}\left ( {\sqrt {{x^2} + {y^2} + {{(z + c)}^2}} - \sqrt {{x^2} + {y^2} + {{(z - c)}^2}} } \right ), \\[-28pt] \nonumber \end{align}
\begin{align} \tau = \frac {1}{{2c}}\left ( {\sqrt {{x^2} + {y^2} + {{(z + c)}^2}} + \sqrt {{x^2} + {y^2} + {{(z - c)}^2}} } \right ), \\[-28pt] \nonumber \end{align}
\begin{align} \varphi = \arctan \left ( {\frac {y}{x}} \right ). \\[-12pt] \nonumber \end{align}
In the equations,
$-1 \leqslant \zeta \leqslant 1$
,
$1 \leqslant \tau$
and
$0 \leqslant \varphi \lt 2\pi$
. The body-fixed Cartesian coordinates of a point on the swimmer surface
$\boldsymbol r_s = (x_s,y_s,z_s)^T$
are
$x_s=b_x \sqrt {1-\zeta ^2} \cos \varphi$
,
$y_s=b_x \sqrt {1-\zeta ^2} \sin \varphi$
and
$\,z_s=b_z \zeta$
. All points with
$\tau = \tau _0 \equiv e^{-1}$
lie on the ellipsoid’s surface. The unit vectors
$\boldsymbol{e}_\tau$
and
$\boldsymbol{e}_\zeta$
on the ellipsoid surface are related to those of the Cartesian coordinates via
\begin{align} \begin{aligned} & \boldsymbol{e}_\tau =\left (\frac {2 \tau \boldsymbol{\cdot }\sqrt {1-\zeta ^2}}{\sqrt {\tau ^2-1}} \boldsymbol{e}_x+\zeta \boldsymbol{e}_z\right ) \boldsymbol{\cdot }\frac {\sqrt {\tau ^2-1}}{\sqrt {\tau ^2-\zeta ^2}}, \\[5pt]& \boldsymbol{e}_\zeta =\left (\frac {2 \zeta \boldsymbol{\cdot }\sqrt {\tau ^2-1}}{\sqrt {1-\zeta ^2}} \boldsymbol{e}_x+\tau \boldsymbol{e}_z\right ) \boldsymbol{\cdot }\frac {\sqrt {1-\zeta ^2}}{\sqrt {\tau ^2-\zeta ^2}}. \end{aligned} \end{align}
In the squirmer model, the propulsion is generated by imposing surface squirming velocities. The surface velocity is assumed to be tangential, steady and axisymmetric. The general form of the surface velocity is expressed as (Blake Reference Blake1971)
\begin{align} \boldsymbol{u}_{s}=-\sum _{n=1}^{\infty } \frac {2 B_n}{n(n+1)} P_n^{\prime }(\boldsymbol{p} \boldsymbol{\cdot }\boldsymbol{n}) \boldsymbol{p} \boldsymbol{\cdot }(\boldsymbol{I}-\boldsymbol{n n}), \end{align}
where
$\boldsymbol{p}$
is the unit orientation vector,
$\boldsymbol{n}=\boldsymbol{e}_\tau$
is the unit normal vector and
$\boldsymbol{I}$
is the unit tensor. Here
$P'_n$
is the derivative of the
$n$
th-order Legendre polynomial. The coefficients
$B_n$
are often called polar modes. We consider the prescribed surface velocity (A6) with the first two polar modes (
$B_1$
and
$B_2$
); then it is simplified as (2.1).
A.2. Implementation of the fictitious domain method
We employ the augmented Lagrangian method proposed by Dean & Glowinski (Reference Dean and Glowinski2002) for the solution of the constitutive equations (2.4), which are rewritten as
\begin{align} \boldsymbol{\tau }=2 \mu \boldsymbol{\dot \gamma }+\tau _y \boldsymbol{B}. \end{align}
Tensor
$\boldsymbol{B}$
is a non-dimensional tensor which satisfies
$|{\boldsymbol{B}}| \leqslant 1$
,
${\boldsymbol{B}}:{\boldsymbol{\dot \gamma }} = {{\dot \gamma }}$
and the following relation (Dean & Glowinski Reference Dean and Glowinski2002):
\begin{align} \boldsymbol{B}=P_{\varLambda }\left (\boldsymbol{B}+\sigma \tau _y \boldsymbol{\dot \gamma }\right ) \quad \forall \,\sigma \gt 0, \end{align}
where
$P_{\varLambda }(\boldsymbol{q})$
is a symmetry-preserving orthogonal projection operator defined by
\begin{align} P_{\varLambda }(\boldsymbol{q})= \begin{cases}\boldsymbol{q} & \text{if}\quad |\boldsymbol{q}| \leqslant 1, \\[4pt] \boldsymbol{q} /|\boldsymbol{q}| & \text{if}\quad |\boldsymbol{q}|\gt 1 .\end{cases} \end{align}
Here
$\boldsymbol{q}$
represents a non-dimensional tensor. Equation (A8) implies that
$\boldsymbol{\dot \gamma } = \boldsymbol{0}$
when
$| \boldsymbol{B}+\sigma \tau _y \boldsymbol{\dot \gamma } | \leqslant 1$
, and
$\boldsymbol{B} = \boldsymbol{\dot \gamma }/{\dot \gamma }$
when
$| \boldsymbol{B}+\sigma \tau _y \boldsymbol{\dot \gamma } | \gt 1$
, and thus (A7) is equivalent to (2.4).
The body force
$\boldsymbol f$
on the right-hand side of (2.2) enforces a prescribed velocity boundary condition on swimmer surface
$S$
, given by
\begin{align} \boldsymbol{u}=(\boldsymbol{U}+\boldsymbol{\omega } \times \boldsymbol{r})+\boldsymbol{u}_{s}\;\; {\textrm {in}} \;\;{S}. \end{align}
To recover the tangential velocity (A10) on the surface of the particle, a solenoidal velocity field
$\boldsymbol{u}_a$
inside the squirmer is imposed (Xia et al. Reference Xia, Yu, Lin, Lin and Hu2025b
,Reference Xia, Yu, Zhang, Lin and Ouyang
d
). The imposed velocity leads to a zero translational and rotational velocity for the particle (i.e.
$\int _{P} \boldsymbol{u}_a \,{\textrm {d}} V=0$
and
$\int _{P} \boldsymbol{r} \times \boldsymbol{u}_a \,{\textrm {d}} V=0$
).
For convenience, we write the dimensionless quantities in the same form as their dimensional counterparts, unless otherwise specified.
The complete set of dimensionless governing equations comprises the following three parts.
(i) Equations of motion:
\begin{align} Re \left (\frac {\partial \boldsymbol{u}}{\partial t}+\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}\right ) =-\boldsymbol{\nabla }p+ \boldsymbol{\nabla} ^2 \boldsymbol{u}+Bi \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{B}+\boldsymbol{f} , \\[-28pt] \nonumber \end{align}
\begin{align} \boldsymbol{u}-(\boldsymbol{U}+\boldsymbol{\omega } \times \boldsymbol{r})-\boldsymbol{u}_a=\textbf{0}. \\[0pt] \nonumber \end{align}
(ii) Continuity equation:
\begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}=\textbf{0}. \end{align}
(iii) Constitutive equation:
\begin{align} \boldsymbol{B}=P_\varLambda \left (\boldsymbol{B}+\sigma Bi \boldsymbol{\dot \gamma }\right ) \quad \forall \,\sigma \gt 0 . \end{align}
The spatial derivatives in (A11) are evaluated using the second-order central-difference method. The convection and viscous terms are temporally integrated using the second-order Adams–Bashforth method and Crank–Nicholson method, respectively. Throughout this study, we set
$\sigma =Re/(2Bi^2)$
in (A14). This choice ensures numerical stability and has been widely adopted in previous simulations of Bingham flows (Dean & Glowinski Reference Dean and Glowinski2002). Problem (A11)–(A14) is solved in an iterative way with the fractional step time scheme (Yu & Wachs Reference Yu and Wachs2007). Particle quantities in (2.5)–(2.8) are updated using the first-order explicit Euler method. An efficient parallel algorithm based on message passing interface has been implemented to speed up computation.
Appendix B. Effect of the domain size
The effect of the domain width
$L$
(as well as the swimmer number
$N$
) has been examined by using a larger width
$L = 64a$
. Figure 10 shows the results for neutral swimmers and pullers. We can see that although the polar order evolution curves of neutral swimmers and the cluster size of pullers are affected by the domain size, the difference in the scaled results is negligible, indicating that the domain size does not influence the conclusions from our work. Similar trends were observed in viscoplastic fluids, and are therefore not shown here for brevity.
Effect of domain size on (a) temporal evolutions of the global polar order parameter
$P(t)$
for neutral swimmers (
$\beta =0$
,
$\phi _a=0.25$
) and (b,c) normalised cluster-size distribution function
$\mathcal{N}(n)$
for pullers (
$\beta =3$
,
$\phi _a=0.25$
) in Newtonian fluids. Insets in (b,c) show snapshots of the clusters.
Figure 10. Long description
Panel A: A line graph shows the temporal evolution of the global polar order parameter for neutral swimmers with domain sizes L = 32a and L = 64a. The x-axis represents B1t/a, and the y-axis represents P(t). The blue line corresponds to L = 32a, and the orange line corresponds to L = 64a. Both lines show an increase in P(t) over time, with the blue line reaching higher values more quickly than the orange line. Panel B: A histogram with an inset scatter plot shows the normalized cluster-size distribution function for pullers with domain size L = 32a. The x-axis represents the cluster size n, and the y-axis represents N(n). The inset scatter plot shows a snapshot of the clusters, with different colors representing different clusters. The histogram shows a peak around n = 80. Panel C: A histogram with an inset scatter plot shows the normalized cluster-size distribution function for pullers with domain size L = 64a. The x-axis represents the cluster size n, and the y-axis represents N(n). The inset scatter plot shows a snapshot of the clusters, with different colors representing different clusters. The histogram shows a peak around n = 200.
τ=τy
Bi=0.1
β=0
x−y
y−z
β=3
β=−3
β=0
Bi=0
Bi=0.1
β=3
β=−3
β=0
Bi=0
Bi=1
β=3
Bi=1
ϕa=0.25
ar=1
ϕa=0.25
ar=3
ϕa=0.45
ar=1
ψ6,k
N(n)
β=3
ϕa=0.3
Bi
β=3
EN
Bi
ϕa
Bi
Bi
β=−3
ϕa=0.3
ar=3
Bi
β=0
Bi
ar=3
ϕa=0.3
Bi
ϕa
ϕa
Bi
P⩾0.7
ϵ
rH
10−2
rH=0.3Bi−1/2
P(t)
β=0
ϕa=0.25
N(n)
β=3
ϕa=0.25
Open access
