Soft Matter
PAPER
Hydrodynamic interactions of self-propelled swimmers
Cite this: Soft Matter, 2013, 9, 4923

John J. Molina,*a Yasuya Nakayamab and Ryoichi Yamamoto*a
The hydrodynamic interactions of a suspension of self-propelled particles are studied using a direct
numerical simulation method which simultaneously solves for the host ﬂuid and the swimming particles.
A modiﬁed version of the “Smoothed Proﬁle” (SP) method is developed to simulate microswimmers as
squirmers, which are spherical particles with a speciﬁed surface-tangential slip velocity between theparticles and the ﬂuid. This simpliﬁed swimming model allows one to represent diﬀerent types of
propulsion (pullers and pushers) and is thus ideal to study the hydrodynamic interactions among
swimmers. We use the SP method to study the diﬀusive behavior which arises due to the swimming
motion of the particles, and show that there are two basic mechanisms responsible for this phenomena:
the hydrodynamic interactions caused by the squirming motion of the particles, and the particle–particle
collisions. This dual nature gives rise to two distinct time- and length-scales, and thus to two diﬀusion
coeﬃcients, which we obtain by a suitable analysis of the swimming motion. We show that thecollisions between swimmers can be interpreted in terms of binary collisions, in which the eﬀective

Received 12th January 2013
Accepted 25th March 2013collision radius is reduced due to the collision dynamics of swimming particles in viscous ﬂuids. At short
time-scales, the dynamics of the swimmer is analogous to that of an inert tracer particle in a swimming

DOI: 10.1039/c3sm00140g

suspension, in which the diﬀusive motion is caused by ﬂuid-particle collisions. Our results, along withthe simulation method we have introduced, will allow us to gain a better understanding of the complex

www.rsc.org/softmatterhydrodynamic interactions of self-propelled swimmers.

1

Introduction

Swimming microorganisms, from bacteria, to algal cells, tospermatozoa, are ubiquitous in biological processes, and eventhough the specic propulsion mechanism can vary, the motionof these microswimmers is dened by two basic characteristics:(1) they are moving at very low Reynolds number, where viscous
forces are dominant (inertial forces can be neglected), and thus
(2) the net force on the body (including the cilia or agella used
to generate the motion) is zero.1,2 Although we have a clear
understanding of how these organisms can generate motion,
the physical properties of a suspension of such swimmers arestill not completely understood.3–7 In particular, the rheologicalproperties show a non-trivial dependence on the concentrationof swimmers,8–10 which deviates considerably from that of inert
force-free colloidal particles.11 In addition, hydrodynamic
interactions among swimmers have also been shown to give rise
to collective motion.12,13 However, theoretical studies on suchphenomena have either neglected the role of hydrodynamic
interactions, or used only far-eld approximations.7,14
In addition to the biological interest of understanding thehydrodynamic interactions of swimming microorganisms,
recent experimental work has shown how to design microa

Department of Chemical Engineering, Kyoto University, Kyoto 615-8510, Japan.
E-mail: john@cheme.kyoto-u.ac.jp; ryoichi@cheme.kyoto-u.ac.jpbDepartment of Chemical Engineering, Kyushu University, Fukuoka 819-0395, Japan

This journal is ª The Royal Society of Chemistry 2013

motors which can swim in the absence of external elds thanks
to interfacial phoretic eﬀects,15–19 which allow the particles totransform the local chemical or thermal energy into mechanicalenergy. Furthermore, the possibility of using these systems toperform work, whether it be to transport cargo20,21 or separate
particle mixtures,22 has led to many recent experimental andsimulation studies. Although a simple theoretical framework
for these swimming systems, which can be used to study theeﬀect of the shape and surface properties on the swimming
motion, has already been proposed,23,24 the collective dynamics,
and in particular the role of the hydrodynamic interactions
remains an open question. In this paper, we present a direct
numerical simulation method (DNS) for self-propelled particles
which attempts to overcome this diﬃculty, by solving theequations of motion for the host uid as well as the swimming
particles.
A detailed description of the swimming mechanism of real
microorganisms remains computationally prohibitive, particularly if one is interested in the collective dynamics of such
systems; therefore, we must look for a simple generic swimmer,
which can be easily modeled, but which manages to reproduce
the basic ow properties of actual swimmers, such that thehydrodynamic interactions are accurately reproduced. The“squirmer” model we use was introduced over 50 years ago25,26
to study the propulsion of spherical ciliate particles, which
generate motion through the synchronized beating of anenvelope of cilia at their surface. Fortunately, on a mesoscopic

Soft Matter, 2013, 9, 4923–4936 | 4923

Soft Matterscale, the eﬀect of this surface motion can be replaced by aspecied slip velocity between the particle and the uid. In
addition, the parameters of the model can be tuned to mimic
diﬀerent propulsion mechanism, allowing us to study thehydrodynamic interactions of diﬀerent (idealized) microswimmers within a unied framework. We present a modied“Smoothed Prole” method,27 previously developed to study
colloidal dispersions, which is capable of incorporating thissquirming motion by imposing the appropriate slip boundary
conditions at the surface of the particles. The advantage of thismethod is that it allows us to accurately and eﬃciently include
the many-body hydrodynamic interactions among the particles,
and it can be easily extended to complex host uids.
In this work we study the eﬀect of hydrodynamic interactions
on the dynamics of a suspension of self-propelled particles. The
swimming motion of the particles is known to give rise to adiﬀusive behavior,28 even in the absence of thermal uctuations,and a tracer particle placed in such a system will also undergo
diﬀusion.29,30 However, the nature of the diﬀusion is diﬀerent inboth cases, as is evidenced by the diﬀerent scaling behavior withthe concentration of swimmers.31 In the former case, the diﬀusion is caused primarily by particle–particle collisions, while
the latter is due to the hydrodynamic interactions caused by thesquirming motion of the particles. By a suitable analysis of theparticle displacements, and the decay in the velocity uctuations, we show that both phenomena are present in the motion of
the squirmers themselves, as should be expected. We are thus
able to extract the two underlying time- and length-scales of thesystem, corresponding to the two distinct diﬀusive motions, justfrom the particle trajectories. The importance of this cannot be
understated, as these two mechanisms (hydrodynamics andparticle collisions) are known to be fundamental in determining
the physical properties of swimming suspensions;32 additionally, measuring tracer diﬀusion in dense suspensions experimentally can be quite challenging, and the available simulation
methods for this are somewhat involved.31 Finally, using basic
results from kinetic theory, we analyze the collisions of theswimmers to derive eﬀective collision radii, which are shown to
be independent of concentration, and which depend only on thesquirming mode of the swimmers (i.e., the strength of thepusher/puller character). Somewhat surprising is the fact thatthese eﬀective particle sizes are considerably smaller than theactual size of the particles, which is a consequence of the collision dynamics of swimmers in viscous uids. This paper isorganized as follows: Section 2 introduces the mathematicalmodel of the swimmers and the simulation method, Section 3
validates the computational method against known analyticalresults, and Section 4 presents the results of our study.

2

The model and methods

2.1

Swimmer model: Blake's squirmers

We consider a simple model of self-propelled spherical swimmers, originally introduced by Lighthill25 and later extended byBlake,26 which move due to a self-generated surface-tangential
velocity us. This specic mechanism was proposed as a model
for an ideal ciliate particle, in which the synchronized beating

4924 | Soft Matter, 2013, 9, 4923–4936

Paperof the cilia at the surface gives rise to a net motion, in theabsence of any external elds. If one assumes that thedisplacements of this envelope of cilia are purely tangential,
then the eﬀective (time-averaged) slip velocity for thesesquirmers is given by26N X
us r^ ¼  
2rr^  e^ P n0 e^$^r
Bn e^$^n ðn þ 1Þ
n¼1(1)

where ^e is the squirmer's xed swimming axis (i.e., we consider
that each squirmer carries with it a xed coordinate system which
determines its preferred swimming direction at each instant), ^r isa unit vector from the particle center to a point on the surface, Pn0 is
the derivative of the n-th order Legendre polynomial, and Bn is theamplitude of the corresponding mode. Neglecting all squirming
modes higher than three, Bn ¼ 0 (n $ 3), the following simple
expression for the surface tangential velocity, as a function of thepolar angle q ¼ cos1(^r $^e), is obtained


a
^
(2)
us ðqÞ ¼ B1 sin q þ sin 2q q
2where a ¼ B2/B1 determines whether the swimmer is a pusher(a < 0) or a puller (a > 0). An example of the former are spermatozoa and most bacteria, of the latter the unicellular alga
Chlamydomonas. A schematic representation of the ow prolegenerated by these two types of swimmers is given in Fig. 1.
Although the squirmer model we adopt forgoes describing thedetailed propulsion mechanism, it is capable of distinguishing
between pushers/pullers and provides an adequate approximation to the far-eld ow prole generated by these swimmers. For Newtonian uids, which is the only case considered
here, the swimming speed U of the squirmer is determineduniquely by the rst mode B1, irrespective of the size of theFig. 1 Schematic representation of the propulsion mechanism and ﬂow proﬁlesof a pusher and a puller, (a) and (b) respectively. These swimmers can be represented using Blake's squirming model, in which the detailed propulsion mechanism is replaced by a speciﬁed slip velocity at the surface of the particles, (c) and(d), for pushers and pullers, respectively.

This journal is ª The Royal Society of Chemistry 2013

PaperSoft Matterparticle, as U ¼ 2/3B1, while the second mode gives the strengthof the stresslet.33,34 The velocity eld generated by a single suchsquirmer, in the Stokes regime, was solved analytically by Ishikawa et al.33 and is given, in the laboratory frame (uid at rest
far away from the particle), by 
a2 a 2^uðrÞ ¼ B1 2e^ þ sin qq
r r 3
þ

a
2

 2


a
a2
^

1
3cos2 q  1 r^ þ 2 sin 2qq
2
r
r

(3)

where a is the radius of the particle. Notice that for neutralswimmers (a ¼ 0) the velocity eld decays as r3, while forpushers/pullers (a s 0) the velocity eld decays as r2. In
contrast, the velocity eld for a sedimenting particle (or a particle
experiencing a net body force) decays as r1.35 This will haveimportant consequences when considering the hydrodynamic
interactions of suspensions of swimmers. Finally, although thesquirmer model might seem overly simplistic, recent experimental work, focusing on spherical droplets which propel due toself-generated Marangoni stresses at their surface, has shownthat it is possible to create such swimmers articially.36

2.2

general mechanism used to obtain the squirming motion is thesame in all three cases: local conservation of momentum. For
the moment though, these DNS approaches have not been
extensively used to study these types of swimming systems; themost popular methods, which still account for the hydrodynamic interactions, have usually been based on Stokesian
Dynamics,8,28,31,33,47–49 and are thus limited to Newtonian uids
in the Stokes regime.
In what follows, we briey review the governing equations fora dispersion of inert colloids in a simple Newtonian-uid,
before considering how the equations must be modied for use
with the SP method for swimming particles with slip boundary
conditions. The formulation we present closely follows that ofref. 37, in which a more detailed description of the computational algorithm can be found. The motion of the host uid isdetermined by the Navier–Stokes equation with the incompressibility condition

This journal is ª The Royal Society of Chemistry 2013

(4)r(vt + uf$V)uf ¼ V$s(5)where r is the total mass density of the uid, uf is the host uid
velocity eld, and s is the stress tensor

Simulation method: basic equations

We propose a direct numerical simulation (DNS) procedure to
study a system of self-propelled squirmers based on the“Smoothed Prole” (SP) method,27,37 which allows one to eﬃciently solve both the Navier–Stokes equation (NS), for the uid
motion, and the Newton–Euler equations, for the colloids. Thismethod has been successfully used to study the diﬀusion,
sedimentation, electro-hydrodynamics, and rheology ofcolloidal dispersions in incompressible uids,38–41 and recent
work has shown how it can be extended to treat compressible
uids within a uctuating-hydrodynamics approach.42 A
detailed error analysis of the method can be found in ref. 43.
The basic idea is to replace the sharp boundary at the colloid–uid interface with a diﬀuse interface of nite thickness z. This
allows one to discretize the system using a xed Cartesian
coordinate grid, since the interface will always be supported bymultiple-grid points. Although a loss of accuracy at the surfaceof the particle is inevitable, we can easily impose periodic
boundary conditions (PBC), and use a Fourier spectral method
to solve for the uid equations of motion. The particles in the SP
method are not treated as boundary conditions for the hostuid, but rather as a body force in the NS equation. Thus, weavoid the mesh-reconstruction problems that plague mostcomputational uid dynamics methods for systems withmoving boundaries. We are aware of two alternative simulation
methods that aim to describe these squirmer suspensions at thesame level of description, the rst was developed by Ramachandran et al.44 using a Lattice Boltzmann (LB) model, andthe second was originally introduced by Downton and Stark45
within a multi-particle collision dynamics (MPC) framework,and later extended by Götze and Gompper46 to recover thecorrect rotational dynamics. Although the implementation
details are specic to each of the models (LB, MPC, SP), theV$uf ¼ 0

s ¼ pI + s0(6)

s0 ¼ h[Vuf + (Vuf)t](7)

with h the shear viscosity of the uid. Consider a mono-disperse
system of N-spherical particles, of radius a, mass Mp, andmoment of inertia Ip ¼ 2/5Mpa2I (with I the unit tensor). The
evolution of the colloids is given by the Newton–Euler
equations,50
R_ i ¼ V i
Mp V_ i ¼ F H þ F C þ F extiii

Q_ i ¼ skewðUi ÞQi_ i ¼ N H þ N ext
I p $U(8)iwhere Ri and Vi denote the center of mass position and velocityof particle i, respectively, Qi is the orientation matrix† and Ui theangular velocity, with skew (Ui) the skew-symmetric angularvelocity matrix
0
1
0
Uzi
Uyi
z
xA
@
Ui
0
Ui
(9)
skewðUi Þ ¼
Uyi
Uxi
0
The forces on the particles are comprised of hydrodynamic
contributions arising from the uid–particle interactions F H,
the colloid–colloid interactions due to the core potential of theparticles F C (which prevents particle overlap), and a possible
external eld contribution F ext (such as gravity). Likewise, thetorques on the particles can be divided into a hydrodynamic N H
and an external contribution N ext (for simplicity, the particle–
particle interactions are assumed to be given by a radialpotential). In what follows we consider buoyancy-neutral† For numerical stability we use quaternions, and not rotation matrices, torepresent the rigid body dynamics of the particles.

Soft Matter, 2013, 9, 4923–4936 | 4925

Soft Matter

Paper

particles, so that F ext ¼ N ext ¼ 0. Finally, conservation of
momentum between the uid and the particles implies thefollowing hydrodynamic force and torque on the i-th particleÐFH
i ¼ dSi$s

(10)

Ð
NH
i ¼ (x  Ri)  (dSi$s)(11)Ð
where dSi indicates an integral over the particle surface.
2.3

Simulation method: smoothed prole squirmers

We now present the computational algorithm used to simulate
the motion of spherical particles, with a given surface tangential
slip velocity us, using the SP method. We require that all eld
variables be dened over the entire computational domain(uid + particle). The concentration eld for the colloids is givenP
as f(x,t) ¼ N
i¼1fi(x,t), where fi ˛ [0, 1] is the smooth proleeld of particle i. This eld is dened such that it is unity within
the particle domain, zero in the uid domain, and smoothlyinterpolates between both within the interface regions. Details
on the specic denition and the properties of this prolefunction can be found in ref. 27. The particle velocity eld isdened in a similar fashion, as
fup ðx; tÞ ¼

N
X

where the pressure term p* in eqn (15) is determined by the
incompressibility condition V$u* ¼ 0. The remaining updatingprocedure applies to the slip condition at the particle boundary
and the rigidity constraint on the velocity eld.
We now consider the momentum change needed to maintain the slip velocity at the surface of each of the squirmers,
where the slip prole us is imposed with respect to the particle
velocities {Vi0 ;Ui0 }, using the previously updated positions andn+1
orientations {Rn+1
i ;Qi }. We note that at this point we do not yet
n+1
know the updated particle velocities {V n+1
i ;Ui }, which are thevalues that should be used when enforcing the surface slip
prole Vi0 ¼ V n+1(Ui0 ¼ Un+1ii ). Therefore, we adopt an iterative
solution, and as an initial guess, we use the particle velocities atthe previous time step, i.e., Vi0 ¼ Vni (Ui0 ¼ Uni). The updated total
velocity eld is now
2 t þh3nðu** ¼ u* þ 4
dsf sq 5(18)tn

3
2 t þh
nð
N
X
 0
0
4
dsf sq 5 ¼
4i V þ U  ri þ us  u*
tn

i¼1

þ
V i ðtÞ þ Ui ðtÞ  ri ðtÞ fi ðx; tÞ

(12)

i

i

i

N
X

h
fi ðdV i þ d Ui  ri Þ  Vpsq
r
i¼1

(19)

i¼1

with ri ¼ x  Ri, which allows one to dene the total uid
velocity eld asu(x,t) h (1  f)uf + fup(13)

where the incompressibility condition is satised on the entire
domain V$u ¼ 0. The evolution equation for u is then derived
assuming momentum-conservation between uid and particles27r(vt + u$V)u ¼ V$s + rffp + rfsq

(14)

where ffp represents the force density eld needed to maintainthe rigidity constraint on the particle velocity eld and fsq is theforce density eld generated by the particles' squirming motion.
We use a fractional step approach to update the total velocity
eld. Let un be the eld at time tn ¼ nh (h is the time interval).
We rst solve for the advection and hydrodynamic viscousstress terms, and propagate the particle positions (orientations)using the current particle velocities. This yields
*

n



tnð
þh

u ¼u þ

dsV$
tn

1 *
p I þ s0  uur

¼ Rni þ
Rnþ1
i

(15)

tnð
þh

dsV i

(16)

dsskewðUi ÞQi

(17)

tn

¼ Qni þ
Qnþ1
i

tnð
þh

tn

4926 | Soft Matter, 2013, 9, 4923–4936

The second term on the right hand side of eqn (19)
imposes the slip velocity prole us at the surface of each of
the squirmers; where 4i f (1  fi)|Vfi| is a smooth surface
prole function which is non-zero only within the interface
domain of the squirmer (normalized such that |Vfi| has amaximum value of one), and zero everywhere else. The third
term adds a counter-ow entirely within the particle domain,in such a way that local momentum conservation is obtained.
Assuming rigid-body motion, with velocities dVi and dUi, thisrequiresÐÐ
dxfi(dVi + dUi  ri) ¼  dx4i(Vi0 + Ui0  ri + usi  u*) (20)ÐÐdxri  fi (dVi + dUi  ri) ¼  dxri  4i(Vi0 + Ui0  ri
+ usi  u*)(21)

from which we can easily obtain the counter-ow terms dVi(dUi) from the particle velocities Vi0 (Ui0 ). A schematic representation of this procedure, used to enforce the specic slipboundary conditions for our model squirmers, is shown in
Fig. 2. Finally, the pressure term due to the squirming motion
psq is obtained from the incompressibility condition V$u** ¼
0. At this point, the momentum conservation is solved for the
total velocity eld.
The hydrodynamic force and torque exerted by the uid on
the colloids (which includes all contributions due to thesquirming motion) is again derived by assuming momentum
conservation. The time integrated hydrodynamic force andtorque over a period h are equal to the momentum exchangeover the particle domain

This journal is ª The Royal Society of Chemistry 2013

PaperSoft Matter
2.4

Fig. 2 Schematic representation of the updating scheme used to enforce
the slip boundary condition at the surface of the squirmers. Each particle isconsidered to exert a force on the ﬂuid at the surface, in order to maintainthe speciﬁed ﬂow proﬁle us (red arrows) for the squirming motion. To ensure
local momentum conservation, a counter-ﬂow is added within the particle
domain (blue arrows).

3
2 t þh
nð
ð




sq 5
H
4
¼ dx rfnþ1 u**  un
ds F þ F i
i

i

p

(22)

tn

3
2 t þh
nð
ð





5 ¼ dx rnþ1  rfnþ1 u**  un
4
ds N H þ N sq
i
i

i

i

p

(23)

tn

From this and any other forces on the colloids, the particle
velocities are updated as

Choosing the relevant reference frame

Although we have not included thermal uctuations in oursystem, the swimming motion of the particles is known to giverise to diﬀusive behavior at suﬃciently high particle concentrations.28 However, due to the self-propelled nature of theparticles, the largest contribution to their displacement will
naturally come from the their inherent swimming. As such, onemust wait a very long time before the particles exhibit any typeof diﬀusive behavior, and even then, it is diﬃcult to establish
what role the hydrodynamic interactions among neighboring
particles are playing. All particles will be swimming in the oweld generated by their neighbors, and the interactions among
them can give rise to motion perpendicular to the particle's
preferred swimming direction, as well as provide a momentary
impulse that can increase/decrease the velocity parallel to theswimming axis, or even change its orientation in space. In order
to better understand this phenomena, we analyze the particle
motion with respect to the frame of reference of the moving
squirmers, as was proposed by Han et al.51 to study the Brownian motion of ellipsoidal particles. We begin by decomposing
the particle trajectories in terms of displacements between
successive time intervals, which we take to be equally spaced. If
R(tn) denotes the position of a tagged particle at time step n (tn ¼ndt, with dt a suitably small time interval), we can express thetime evolution of its position {R(ti)} as
)
(
i
n
X
X
Dj ; .; Rðt0 Þ þDn(28)Rðt0 Þ; Rðt0 Þ þ D1 ; .; Rðt0 Þ þ
j¼1

j¼1

3
3
2 t þh2 t þhnðnð where Di ¼ R(ti)  R(ti1). Using this notation, the translational


sq 5
nþ1
n
H
C
ext 5
1 4
1 4
þ Mp
ds F i þ F i
ds F i þ F i
V i ¼ V i þ Mp
diﬀusivity is given as52
tn
tn
*"
#2 +
n
E
X
1 D
1
(24)
2
DT ðtn Þ ¼
Dn
¼
(29)
Rðtn Þ  Rðt0 Þ
6ndt
6ndt
i¼1
2 t þh2 t þh
3
3
nð
nð


sq 5
H
1 4
1 4
n
5
þ
I
Unþ1
¼
U
þ
I
$
ds
N
þ
N
$
ds N ext
To study the coupling between translation and rotation of
p
p
i
i
i
i
i
the particles, we will also consider the rotational diﬀusivity,tntn(25) dened asE
1 D 
2
We recall that we have imposed the slip prole us withDR ðtn Þ ¼(30)
w tn  wðt0 Þ
6ndt
respect to the primed velocities {Vi0 ;Ui0 }, which need not be
equal to the nal velocities of the particle at step n + 1. Toðtnmaintain consistency, we iterate over eqn (18)–(25) untilwðtn Þ ¼ dsU ðsÞ(31)convergence in the velocities is achieved. Finally, the resulting
t0
particle velocity eld fn+1un+1
is enforced on the total velocitypwhere w is the unbounded rotational displacement of theeld asparticle.
2 t þh3
In order to separate the “random” contribution given by thenð(26) surrounding conguration from the particle's own swimming
dsf f p 5
unþ1 ¼ u** þ 4
motion, we consider the displacements within the body frametnof reference,
2 t þh3nðh^ i ¼ Qti1Dic(32)
4
dsf f p 5 ¼ fnþ1 unþ1
 u**  Vpp
(27)
p
r
tn
where tildes are used to denote quantities with respect to thewith the pressure due to the rigidity constraint obtained from coordinate-system attached to each of the squirmers (^e1,^e2,^e3),
the incompressibility condition V$un+1 ¼ 0. The total pressure with ^e3 the preferential swimming axis (see Fig. 3). The advantage of this approach is illustrated in Fig. 4, where the trajectoryeld is then given by p ¼ p* + pp + psq.
This journal is ª The Royal Society of Chemistry 2013

Soft Matter, 2013, 9, 4923–4936 | 4927

Soft Matter

Paperof a single particle, from a suspension of pushers a ¼ +2 at f ¼
0.1, is given in both representations: with respect to the lab- andbody-space displacements. One can immediately see where thediﬃculty in analyzing the particle motion comes from, as thelength scales for the directed (swimming) and uctuating
motion diﬀer by an order of magnitude. We dene an addi~ T, in terms of thesetional eﬀective hydrodynamic diﬀusivity D
perpendicular and parallel displacements, as
#2 +
*"
n
X
1
t
t
~ ðtn Þ ¼
~
(33)
D
D
T
i
4 ndt
i¼1
~k ðtn Þ ¼
D
T

1
2 ndt

*"

n
X
~ k  ndtU
D
i

#2 +
(34)

i¼1

 k

~ ðtn Þ þ D~t ðtn Þ~T ðtn Þ ¼ 1 2D
D
T
T
3

(35)

where the second term inside square brackets in eqn (34) isrequired to remove the (average) contribution to the particle
displacements from the inherent swimming motion (U ¼ hV$^e3i isthe average velocity along the particle's swimming axis). These
diﬀusivities provide a direct measure of the strength of thehydrodynamic interactions among the squirmers. The corresponding diﬀusion coeﬃcients can be obtained from the long-time
limit (assuming a plateau has been reached and the limit exists)
DT ¼ lim DT ðtÞ(36)

DR ¼ lim DR ðtÞ(37)~T ðtÞ
~T ¼ lim D
D

(38)

t/N

t/N

t/N

Since the parallel and perpendicular eﬀective diﬀusion coef~ kT and D~tcients, D
T , exhibit no deviation from isotropic behavior(although this could be expected to change if persistent longrange order appears), we only consider the average eﬀective~ T. A similar analysis can be performed for
diﬀusion coeﬃcient D

Fig. 4 Trajectory of a single particle in a suspension of pullers a ¼ +2 at f ¼ 0.1.
(a) The real lab-space trajectory of the particle, and (b) the body-space trajectory.
The latter is constructed by transforming the individual particle displacements tothe frame of reference of the particle, and expressing them in terms of^3. Indisplacements parallel and perpendicular to the particle's swimming axis e
addition, the (parallel) motion due to the average swimming speed, Dx(t) ¼^3it, has been removed. Thus, (b) gives the random motion induced on thehV$eparticle by the local ﬂuctuations in the surrounding ﬂow ﬁeld. The length scale isgiven by the particle radius a/D ¼ 5 and the projections of the trajectories onto
the bottom plane have been added as a visual guide.

the rotational diﬀusivities, in order to study the rotations about
the three particle axes independently, but it provides no additional information and will not be presented here.
Finally, we also consider the velocity auto-correlation
functions52
CV(t) ¼ hV(t)$V(t0)i(39)

CU(t) ¼ hU(t)$U(t0)i(40)

and in particular, the correlation functions for the velocity
components parallel and perpendicular to the swimming axis
Fig. 3 Analysis of the particle displacements with respect to the particle's
moving coordinate frame. The particle displacement D during a given time
^ t to^ k and perpendicular cinterval is decomposed into its components parallel cthe swimming axis. For a single isolated squirmer c^ t ¼ 0 for any given timeinterval; for a suspension of swimmers the ﬂow induced by the neighboring
particles gives rise to non-zero perpendicular displacements.

4928 | Soft Matter, 2013, 9, 4923–4936

This journal is ª The Royal Society of Chemistry 2013k

CV(t) ¼ hVk(t)$Vk(t0)i(41)

k
C~V(t) ¼ hV~k(t)$V~k(t0)i

(42)

~t
~t
C~t
V (t) ¼ hV (t)$V (t0)i(43)

PaperSoft Matter

Although eqn (41) and (42) both measure the correlations of
the parallel velocity components, the former does so within thexed lab frame, while the latter uses the moving body frame.
Therefore, the rst will be sensitive to changes in the magnitude
and direction of the velocity vector, while the second will only
register changes in the magnitude. At short times, we can expect
the decay in correlations of the parallel velocity components to
be determined by the local conguration of the suspension(where the ow generated by the nearby particles can act toenhance or suppress the swimming motion), and at long times
by the particle collisions, which will reorient the particle's
swimming direction. This means the initial decay should be the~ kV, since the particle has not had time tosame for both CkV and C
reorient itself; at long times however, correlations measured
within the lab frame of reference should decay to zero, whereas
those measured within the body reference should reach a nite
value (determined by the average swimming speed of thesquirmers). For the decay in correlations of the perpendicular~tvelocity components C
V we expect a behavior analogous to theshort-time decay of the parallel components, since it is due
entirely to the ow eld generated by the neighboring squirmers. We note that a similar analysis has been successfully used
to study the interplay between hydrodynamic and Brownianuctuations on the motion of (inert) sedimenting colloids.53

3

Validation

The rst obvious test of our simulation method is to make surethat an isolated swimmer will move at the expected velocity U ¼
2/3B1, regardless of the particle size or the value of the secondsquirming mode B2. We performed simulations for a single
squirmer, inside a periodically replicated cubic simulation box
of dimension L ¼ 128D (D is the grid spacing), for various
particles sizes and squirming modes. Fig. 5 shows the results
obtained for a neutral squirmer (a ¼ 0), for particle sizes a/D ¼
4, 5, 6, 8, and for various pullers (a > 0), with a particle size ofa/D ¼ 6. The particle Reynolds number Re ¼ rUa/h and thewidth of the diﬀuse interface z/D were the same for all thesimulations, 0.01 and 2, respectively. In all cases, the particle's(average) swimming velocity shows excellent agreement withthe theoretical predictions. For the neutral squirmers, theswimming velocity is within x2% of the exact value, regardlessof the particle diameter, although the speed shows a small
decrease with increasing particle size. The small oscillations
exhibited by the velocity are due to discretization errors, as thenumber of grid points on the two hemispheres of the particle
surface will vary depending on the relative position of theparticle center within the computational bins. Similar agreement is obtained for the pullers (pushers), although the discretization error increases with increasing rar (at xed particle
size a/D). As such, for the system sizes we have considered, weare limited to moderate values of rar ( 5.
To verify that the velocity eld generated by our SP squirmers
is correct, we compare with the exact (at Re ¼ 0) analyticalexpression given in eqn (3). The steady-state velocity elds
generated by a single squirmer (L/D ¼ 128, A/D ¼ 6, z/D ¼ 2, a ¼+2, Re ¼ 0.01), along with the corresponding stream line plot (inThis journal is ª The Royal Society of Chemistry 2013

Fig. 5 Swimming speed of an isolated squirmer in a periodically replicated cubic
simulation box of length L/D ¼ 128 at Re ¼ 0.01 as a function of time (in simulation units). (a) Neutral swimmer at various particle sizes a/D. (b) Puller of size
a/D ¼ 6 for various swimming modes 0 # a # 5. Velocities are scaled by thetheoretical value for the swimming speed U ¼ 2/3B1.

the lab and particle reference frames), are shown in Fig. 6.
Excellent agreement with the analytical results is obtained,
although diﬀerences in the stream lines arise at large distances
r/a x 4 for q x p/4 (with respect to the swimming axis +^z). Even
though the uid velocity within these regions is vanishingly
small (compared to the swimming speed of the squirmer), a clear
systematic deviation is observed in the direction of the streamlines. This is due to our use of periodic boundary conditions,which causes the stream lines to close in on themselves, in orderto match at the boundaries of the simulation cell. A similar
deviation is observed for the velocities along the swimming
direction, but these occur at much larger distances.
Following Götze and Gompper,46 we consider the interactions
between two xed squirmers at a distance r from each other, withparallel orientations. Fig. 7 shows the results we have obtained
for the force parallel to the displacement vector between two
pullers (perpendicular to their swimming axes), normalized by
the Stokes force FS ¼ 6phaU for an inert particle moving with thesame velocity (U ¼ 2/3B1). The simulations were carried out at
Re ¼ 0.01 for a box size of L/D ¼ 128, with a particle radius ofa/D ¼ 6, and a swimming mode a ¼ +2. The functional form forthis perpendicular force has been given by Ishikawa et al.33
Fnear f log(3)(44)

Soft Matter, 2013, 9, 4923–4936 | 4929

Soft Matter

Paper

Fig. 7 Perpendicular force felt by two (ﬁxed) parallel squirmers (a ¼ +2), of diameter s/D ¼ 10 and interface thickness z/D ¼ 2, as a function of distance r, for Re ¼ 0.01
and a system size of L/D ¼ 128. The dashed (yellow) line gives the ﬁt to the expectedfunctional form F/Fs ¼ Aln(r/s  1) of the near-ﬁeld force, with A ¼ 0.327, while thedotted (violet) line gives the (exact) far-ﬁeld force. All forces are normalized by the
Stokes force for an inert sphere with the same speed (U ¼ 2/3B1).

to r ( 1.5s, while the far-eld force is approached asymptotically for r T 3s. The force is observed to be proportional to theswimming mode, such that F f a and F(a) ¼ F(a).

Fig. 6 Azimuthally averaged steady-state ﬂuid velocities for a single puller a ¼+2 of size a/D ¼ 6 swimming along the ^z-axis, within a periodically replicatedcubic simulation box of size L/D ¼ 128. (a) Fluid velocity vectors within thelaboratory frame, the red and blue arrows show the simulation and analyticalresults, respectively. (b) Fluid velocity stream lines within the laboratory frame, thecolored lines (color-coded with respect to the magnitude of the velocity) represent the analytical solution, while the gray lines show the simulation results. (c)Same as (b) but within a reference frame moving with the particle. Due to discretization errors, streamlines can begin/end at the ﬂuid–particle interface.
Length scales have been scaled by the particle radius.

where 3 ¼ r  s is the minimum separation distance between thesurface of the particles (with s ¼ 2a the particle diameter). Additionally, they also obtained exact far-eld expressions for the force
between the two squirmers from a generalization of Faxen's laws,
Ffar ¼

3 ps3 a B1
16r2

(45)

At short distances the repulsive force experienced by the twosquirmers is seen to follow the expected scaling relationship up

4930 | Soft Matter, 2013, 9, 4923–4936

4

Results and discussions

4.1

Diﬀusion coeﬃcients

We have studied the diﬀusive behavior of a semi-dilute
suspension of identical non-buoyant squirmers, for various
concentrations f # 0.124 (with f ¼ 4pa3N/3V the packing
fraction) and squirming modes a ¼ 0, 1, 2. We work inreduced units, in which the density and viscosity of the hostuid are unity r ¼ m ¼ 1. All simulations were performed forsquirmers of radius a/D ¼ 5 (interface width z/D ¼ 2), in a cubic
simulation box of length L/D ¼ 64, for a particle Reynolds
number of Re ¼ rUa/m ¼ 0.05 (U is the swimming speed of an
isolated squirmer). Although our systems are not very large L/a 12, nite size eﬀects have been shown to be small,28 and since
we will be mainly focused on studying the short-range hydrodynamic interactions among particles, they can be safely
ignored. In what follows, all quantities are presented in nondimensionalized form, using as characteristic length, speed,
and time units the particle radius a, swimming speed U ¼ 2/3B1,
and the time required for an isolated squirmer to move a
distance equal to its radius T ¼ a/U.
~ T(t), dened in eqn (29) and (35),The diﬀusivities DT(t) and D
for a system of pullers (a ¼ +2), at various concentrations, areshown in Fig. 8. The standard diﬀusivities DT(t) show a ballistic
regime which extends to very long times t x 100, aer which theslope starts to decrease, indicating a transition towards adiﬀusive regime. However, purely diﬀusive motion (slope ofzero) is only obtained for the highest concentrations and alphavalues. This behavior has been analyzed in detail in ref. 28. In~ T(t) reaches the diﬀusive regime at much shorter
contrast, Dtimes, for all concentrations and all non-zero values of alpha.
We note that both quantities are measuring relatedThis journal is ª The Royal Society of Chemistry 2013

Paperphenomena, but the latter does so directly, since the motion
due to the inherent particle swimming has been removed from~ T shows a clear
the analysis. As a function of concentration, D
increase, which is consistent with the interpretation of thediﬀusive motion being caused by interactions with neighboring
particles: as the concentration increases, so does the number ofneighbors, and thus the strength of the interactions. In
contrast, we see that DT decreases with concentration. The
reason for this is simple, since DT mainly measures the eﬀect ofthe particles' own swimming, an increase in the diﬀusive
behavior can only hinder this motion. The rotational diﬀusivities DR, also shown in Fig. 8, present the same basic features as~ T, the onset of the diﬀusive regime is obtained within the sameDtime interval and they exhibit a similar concentration dependence. Apart from a diﬀerence in scale, there is no clear
diﬀerentiating factor between these two quantities. Similar
results are obtained for all other non-zero values of a considered; where, as was pointed out by Ishikawa and Pedley,28 theeﬀect of increasing (decreasing) a is the same as that ofincreasing (decreasing) the concentration, at least with respect

Fig. 8 Translational and rotational diﬀusivities for a system of pullers (a ¼ +2) atvarious concentrations. (a) Translational diﬀusivities obtained from the meansquared displacements in the ﬁxed-lab reference frame DT(t), and the body frame~ T(t) (with the displacement due to the intrinsic swimming motion suitablyDremoved), dashed and solid lines, respectively; and (b) the rotational diﬀusivities
DR(t). Lines of slope one, corresponding to purely ballistic motion, have beendrawn for comparison. The arrows show the increase/decrease of the various
diﬀusivities as a function of concentration f.

This journal is ª The Royal Society of Chemistry 2013

Soft Matterto the diﬀusive motion of the squirmers, since they both lead toan increase (decrease) in the strength of the hydrodynamic
interactions among particles.

4.2

Velocity uctuations and correlations

To have a better understanding of the how the diﬀusive motion
arises, we now consider the single particle velocity correlation
functions dened in eqn (39)–(43). We present in detail the
results obtained for a system of pullers (a ¼ +2) at f ¼ 0.124,
shown in Fig. 9, but the same analysis applies for the othervalues of f(a) we have studied. The standard denition of thevelocity auto-correlation CV(t) results in a function that shows a
two-time decay process, as evidenced by the shoulder thatappears around 1 < t < 10. The initial decay (t x 1) is due to thechanging conguration of neighboring particles (and their selfgenerated ow eld), while the second, long time decay (t x 10)
arises due to the reorientation of the particles' swimming
direction, which is a consequence of the particle–particle
collisions. Conrmation for this interpretation is provided by
the three other correlation functions we have dened. Considerrst the decay in correlations of the parallel velocity compo~ kV (particle frame). Both functions
nents, CkV(t) (lab frame) and C
show the same initial decay, caused by local uctuations in thevelocity eld, as the particle has not experienced any signicantchange in its swimming direction. However, at longer times,
these two functions show a completely diﬀerent behavior. The
parallel correlations measured within the xed lab system CkV(t)exhibit a slow (long-time) decay before eventually collapsing
onto the full correlation function CV(t); the deviations at longtimes t > 50 are due to statistical uncertainties in our data. In
contrast, the correlations measured within the particle's frameof reference reach a plateau aer the initial decay. This meansthat while the particle is constantly reacting to the ow eldgenerated by its neighbors, by changing it's swimming speed,these uctuations decay very fast, typically within the time ittakes for the particle to travel its diameter (the characteristic

Fig. 9 Single-particle velocity time-correlation functions for a system of pullers(a ¼ +2) at f ¼ 0.124. The two-time relaxation approximation to the velocity
correlation function C(2)V (t)/CV(0) (eqn (46)), with the amplitudes and relaxation
times obtained from the simulation, is also shown (dash-dot red line).

Soft Matter, 2013, 9, 4923–4936 | 4931

Soft Matter

Paper

length over which the ow eld varies signicantly). As such,~ kV(t) is due primarily tothe long-time decay in both CV(t) and C
changes in the orientation of the swimming direction, and not
to changes in the magnitude of the swimming speed. As~texpected, the correlations in the perpendicular velocities CV (t)show a very fast decay, over a characteristic time equal to that of
the initial short-time decay of the parallel velocity components.
The decay in correlations of the angular velocity CU(t) (not~tshown) presents the same characteristic behavior as C
V (t).
Finally, the characteristic times ss and sl for these twoprocesses, governed by the hydrodynamic interactions and theparticle collisions, can be obtained by assuming an exponential~t
decay exp(t/s) for the appropriate correlation function: CV (t)kand CV(t) for ss and sl, respectively (see Fig. 9). The two-timescale relaxation of the velocity uctuations can then be
approximated as

particle–particle collisions of the swimmers. From the kinetic
theory of gases,54 we know that the mean free path l (or average
distance between collisions) should be proportional to theinverse of the concentration of particles. In reduced units,
we havepﬃﬃﬃ
2 2 2 1
sc f
l¼
(47)
3
with f the packing fraction, and sc the collision or cross-section
diameter, which is not necessarily equal to the actual physicaldiameter of the particle. This diﬀerence is due to the selfgenerated ow prole around each particle, and as such should
only depend on the swimming mode a, not on the concentration of particles. For the moment though, let us assume that weare dealing with hard-spheres sc ¼ 2, the average time betweenHScollisions is simply tHS
c ¼ l /Uc, which gives

ð2ÞCV ðt; f; aÞ ¼ cs ðf; aÞet=ss ðf;aÞ þ cl ðf; aÞet=sl ðf;aÞ
cs ¼ hdV$dVi
cl ¼ Uk 2(46)

with Uk ¼ hV$^e3i and where we have explicitly shown thedependence of the amplitudes c and decay times s on theconcentration and swimming mode of the particles. Fig. 9shows excellent agreement between this two-time relaxation
process (with all four parameters obtained from the simulationdata) and the full correlation function CV(t) dened in eqn (39).
Furthermore, in what follows we show how the long-time decay
can be simplied and interpreted in terms of a binary collision
process between the swimmers.

4.3

Scaling

To conrm the interpretation we have given for the diﬀusive
motion of the squirmers, and the emergence of the two timescales, we perform a scaling analysis for the relevant parameters(diﬀusion coeﬃcients and collision time scales), as a functionof concentration f. Although a similar analysis has been presented in ref. 28 by Ishikawa and Pedley, they have not~ T, which means that only theconsidered the eﬀective diﬀusion D
long-time behavior of the system, given by the particle collisions, was analyzed (the hydrodynamic interactions which giverise to the rapidly decaying velocity uctuations are completely
masked by the swimming motion). In a subsequent study,31 theauthors considered the diﬀusion of uid particles, as well as
inert colloidal particles, in a suspension of swimmers, and in
this case they were able to observe the emergence of a second(shorter) time-scale. Our results, which are in agreement with
their scaling arguments, provide a complementary view of thediﬀusive motion of these squirmer systems. The benet of theanalysis we propose lies in the fact that all relevant time andlength scales can be obtained just from the motion of thesquirmers, i.e., there is no need to consider the motion of theuid or to introduce inert (non-swimming) particles. This
results in simpler simulations and will also be relevant whentrying to compare with experimental data.
4.3.1 Correlation times. First, lest us consider the longtime dynamics of the system, which is determined by the4932 | Soft Matter, 2013, 9, 4923–4936

tHS
c ¼

pﬃﬃﬃ
2
ðUc fÞ1
6(48)

where Uc is the average velocity of the particles. In thermal
systems, this velocity would be determined by the reservoir andbe independent of concentration, so that one obtains the wellknown f1 dependence for the collision time tHS
f f1.
c
Although the velocity of our squirmers is not concentration
independent, to rst order, we can safely assume that theparticles are all swimming at an average velocity which is near
the swimming velocity of an isolated squirmer, Uc ¼ 1. Therefore, the collision time for our systems should show the same
concentration dependence predicted by kinetic theory, tc f f1.
To obtain the collision times tc from the decay times scomputed from the simulations, we use the Enskog
approximation52,54
tc ¼ 2s/3(49)

Our results, shown in Fig. 10, conrm the scaling predictions
for the long-time decay tlc ¼ 2sl/3. Of special interest is thediﬀerence between tlc for the diﬀerent squirming modes, withthe collision time decreasing with increasing a. This can be
explained by an increase in the eﬀective size of the particle, as
the mean free path l is inversely proportional to the collision
diameter sc. We also show the mean collision time tHS
c for anequivalent system of hard-spheres. The fact that the collisiontimes for the squirmers are signicantly larger than the corresponding hard-spheres values is surprising, since it implies thatthe squirmers have a cross-section diameter which is smaller
than the hard-sphere diameter of the particle (tc f sc2). The
reason for this lies in the collision dynamics of the squirmers.
The scaling of the fast time-scale ss (or tsc) is harder toelucidate, since it shows only a very weak concentration
dependence. Given that ss is related to the velocity uctuations
caused by the short-range hydrodynamic interactions between
particles, it should be related to the time it takes for the ow
region around a particle to change: this is essentially the timenecessary for a particle to swim its diameter t x 2. Our results,
also shown in Fig. 10, agree with this rough estimate and also
indicate a power-law behavior with an exponent of x1/6.

This journal is ª The Royal Society of Chemistry 2013

PaperSoft Matter

Fig. 10 (a) Scaling of the collision times, tlc and tsc, with concentration f. Also
shown is the collision time tHS
c for an equivalent system of hard spheres, with a(scaled) velocity of Uc ¼ 1, given by the kinetic theory of gases. (b) Scaling of the~ T and rotational DR diﬀusion coeﬃcients with concentration. (c)translational D
Scaling of the velocity ﬂuctuations with concentration.

Although the decrease with concentration seems clear, it is very
diﬃcult to accurately measure such small variations and wehave no suitable explanation for the value of this exponent(given by a t to the data). These results are analogous to thoseobtained by Ishikawa et al.31 for the time-scale of diﬀusing uid
particles in a suspension of squirmers.
4.3.2 Diﬀusion coeﬃcients. For the scaling of the translational diﬀusion coeﬃcients, a simple dimensional analysis yields
DT f Uc2tc
This journal is ª The Royal Society of Chemistry 2013(50)

where Uc and tc are the characteristic velocity and time scales of
the collision process which gives rise to the diﬀusive motion.31
For the standard diﬀusion coeﬃcient DT (eqn (36)), the reorientation of the particles is caused by the particle–particle
collisions, thus tc ¼ tlc and DT f f1 (assuming a constantswimming velocity). This behavior has been analyzed in detailin ref. 28. More relevant to our study is the scaling of the~ T (eqn (38)), in which the motion
eﬀective diﬀusion coeﬃcient D
due to the particle swimming has been removed. Our results,
shown in Fig. 10, indicate a linear dependence with concen~ T f f. This is precisely the scaling behavior reported bytration D
Ishikawa et al.31 for the diﬀusion of inert/uid particles, and thesame concentration dependence that has been observed
experimentally for the diﬀusion of tracer particles in swimming
suspensions (for both pushers and pullers).29,30 This is not
surprising, since we can consider that this diﬀusive motion
arises due to collisions between the particle (swimmer or not)and the uid. By denition, it is clear that the velocity and timescales cannot be the same swimming speed or collision timeused to dene DT. In this case, the characteristic velocity scale
will be set by the velocity uctuations dv around the averageswimming velocity, which can be directly measured bycomputing the velocity components perpendicular to theswimming axis of the particles Vt. For these uctuations, oursimulations indicate a square-root dependence with concenpﬃﬃﬃ
tration dvf f (see Fig. 10). The characteristic time-scale forthis process will be the time necessary for the uid owsurrounding a given particle to show considerable variations,and this can be estimated by the time it takes the particle totravel its own diameter tc f 1. We thus recover the linear
dependence given by the simulations. The same scaling
behavior is obtained for the rotational diﬀusion coeﬃcient
DR f f.
4.3.3 Collision diameters. So far, we have not considered
the eﬀect of the squirming mode a on the diﬀusive properties ofthe swimmers. While it is clear that increasing the magnitude ofa will give rise to stronger hydrodynamic interactions, and thusincrease the diﬀusion of the particles (as can be seen in Fig. 10),
it is not clear how this dependence can be quantied. We
propose that a, which denes the strength/range of the selfgenerated uid-ow around a squirmer (eqn (3)), can be directly
related to the eﬀective collision diameter of the swimmers.
From the expression of the mean free path, eqn (47), we obtainthe following for the collision radius rc ¼ sc/2
 pﬃﬃﬃ
1=2
rc ¼ 3 2Uc tc f

(51)

Assuming tc ¼ f(a)f1, with f(a) a function only of thesquirming mode a, and considering the swimming velocity to
be constant (Uc ¼ 1), we arrive at a concentration independent
collision radius
rc f [ f(a)]1/2(52)

As f (a) is a decreasing function of a, the eﬀective collision
radius of the particle should increase with the magnitude of thesquirming mode. Our results for the collision radius, given by

Soft Matter, 2013, 9, 4923–4936 | 4933

Soft Mattereqn (51), with Uc ¼ Uk ¼ hV$^e3i and tc ¼ 2/3sl obtained directlyfrom the simulations, are shown in Fig. 11. Allowing for thelarge uncertainties in measuring the decay times at low
concentrations, we obtain very good agreement with theprevious scaling analysis: the radii show only a small variation
with f, and a clear distinction is observed as rar is increased.
We note however that the pushers seem to present a larger
collision radius than the pullers, something which is diﬃcult to
explain with the simple collision model we have presented.
Also shown in Fig. 11 is the mean collision radius rc as afunction of a. Although more data is needed to accurately
specify the functional form of f (a), a simple linear dependence
predicts a value of rc (a ¼ 0) which is an order of magnitude
smaller than the hard-sphere radius of the particles rc x 0.1.
This is consistent with our simulations, from which we wereunable to obtain reliable estimates for sl, because the velocity
correlation functions exhibited no substantial decay over thetime scales we studied t x 102 to 103. However, we were able to~ T ofobtain estimates for the eﬀective diﬀusion coeﬃcients Dthese neutral squirmers, which are an order of magnitude
smaller than the corresponding values for a ¼ +1. This meansthat the diﬀerence in collision times, with respect to a corresponding system of hard-spheres (see Fig. 10), will be evenPaperlarger for the neutral squirmers than it is for the swimmers with
non-zero a values. This is consistent with an interpretation interms of hard-spheres when analyzing the diﬀerences amongsquirmers (larger a equivalent to larger collision radius rc), but
it can appear contradictory when comparing to actual hardsphere systems. Aer all, we could expect the systems for a ¼ 0,
in which the ow prole due to the squirming motion is most
localized (U(r) f r3), to be closest to a system of hard-spheres.
Yet our results indicate that the collision radius of the neutral
squirmers is an order of magnitude smaller than the actual
radius of the particles. As we have already mentioned, this isdue to the diﬀerence in the collision dynamics of the particles,
which diﬀer dramatically from that of hard-spheres.33,55
Although the exact collision process will depend not only on therelative velocity of the particles, but also on their relativeorientations, the squirming motion results in deection angles
which are (on average) smaller than the corresponding hardsphere values. In essence, this means that the ow elds
generated by the particles allows them to swim past each otherwith a relatively small change in their swimming direction(again, as compared with an actual hard-sphere collision).
Finally, using eqn (51), the velocity correlation function for asuspension of squirmers (eqn (46)) can be reduced to thefollowing functional formpﬃﬃ
CV ðt; f; aÞ ¼ cs et=ss þ Uk 2 e 2sc Uk ft(53)
where the long-time decay process is expressed only in terms of
the average swimming speed Uk(f, a) and the hard-spherecollision diameter sc(a) of the squirmers. The former has a weak
linear dependence on concentration Uk ¼ 1  A(a)f, while thelatter is concentration independent, as was shown in Fig. 11. The
scaling behavior for the short time process is still not understood, but a t to our simulation results suggest a power law
behavior for the concentration dependence of the amplitude anddecay time of the velocity uctuations, with an exponent of 1/2
and 1/6, respectively; i.e., cs f f and ss f f1/6. A more detailed
analysis is required to rmly establish these scaling relationships, as well as to determine the dependence on the swimming
parameter a. Work along these lines is in progress.

5

Fig. 11 Eﬀective collision radius rc of the squirmers as a function of (a)concentration f and (b) swimming mode a. Results were obtained by using themean free path expression for a system of hard-spheres given by the kinetic
theory of gases eqn (51), with the velocity and collision time given by the simulation data.

4934 | Soft Matter, 2013, 9, 4923–4936

Conclusions

We have investigated the hydrodynamic interactions ofsuspensions of squirmers using a modied version of thesmoothed prole method (SP) for particle dispersions. The SP
method allows one to fully resolve the hydrodynamic interactions in many particle dispersions in an accurate and eﬃcientmanner, and we have shown how it can be extended to systems
with slip boundary conditions, such that it is possible todescribe squirmers (active swimmers which move due to selfgenerated surface tangential velocities). The validity of themethod was conrmed by comparing the simulation data withexact results for the case of a single swimmer, for which werecover the correct swimming speed and are able to accurately
reproduce the uid ow generated by the squirming motion,and for two aligned swimmers at a xed distance, for which werecover the expected hydrodynamic force. The advantage of theThis journal is ª The Royal Society of Chemistry 2013

PaperSP method for swimming particles, as opposed to Stokesian
Dynamics (which has been successfully, and extensively, used to
study these systems) is its applicability to particle dispersions incomplex uids. This is relevant in the case of swimming microorganisms, since the role of the nutrient and the possibility of
having a non-Newtonian host uid must be considered whencomparing with experiments.
In this paper we have analyzed the eﬀect of the hydrodynamic interactions on the motion of semi-dilute squirmer
suspensions, up to volume fractions of f ¼ 0.124, for various
swimming modes rar # 2. Although we have no yet included
thermal uctuations in our description, the swimming motionof the particles gives rise, over suﬃciently long time scales, to adiﬀusive regime. In order to distinguish between the contributions due to the hydrodynamic interactions, caused by thesquirming motion, and those due to the particle–particle collisions, which are the two basic mechanisms responsible for thediﬀusive motion, we have analyzed the particle dynamics interms of movement due to the inherent swimming of theparticles, and that due to the (hydrodynamic) interactions
among them. This is easily done by looking at the motion fromthe particle's own frame of reference, i.e., decomposing the
motion parallel and perpendicular to the swimming axis. This
analysis has allowed us to demonstrate the appearance of two
distinct time scales within our system, one related to the time
between particle–particle collisions, the other to the uid–
particle interactions. This two-time scale nature of the particle
interactions in swimming suspensions can be clearly seen in
the two-time relaxation of the velocity correlation function. We
are thus able to dene an eﬀective hydrodynamic diﬀusion
coeﬃcient (corresponding to the short-time uid–particle
interactions), in which the self-motion of the particle has been
removed, which shows a linear scaling with concentration. In
contrast, the standard diﬀusion coeﬃcient is inversely proportional to the concentration of swimmers. This is in agreementwith simulation and experimental results on tracer diﬀusion inswimming suspensions. Additionally, since the long-time
dynamics of the system is related to the particle–particle collisions, we have used the well-known results from kinetic theory
to deduce an eﬀective, concentration independent, collision
radius for our swimmers. Due to the complex collision
dynamics of these particles, this collision radius is actuallysmaller than the hard-sphere radius of the particle, andincreases with increasing a.

Acknowledgements
The authors would like to express their gratitude to the Japan
Society for the Promotion of Science for nancial support(Grants-in-Aid for Scientic Research KAKENHI no 23244087).

References
1 E. M. Purcell, Am. J. Phys., 1977, 45, 3–11.
2 E. Lauga and T. R. Powers, Rep. Prog. Phys., 2009, 72, 096601.
3 M. Manghi, X. Schlagberger, Y.-W. Kim and R. R. Netz, So
Matter, 2006, 2, 653–668.

This journal is ª The Royal Society of Chemistry 2013

Soft Matter
4 R. Golestanian, J. M. Yeomans and N. Uchida, So Matter,
2011, 7, 3074–3082.
5 J. S. Guasto, R. Rusconi and R. Stocker, Annu. Rev. Fluid
Mech., 2012, 44, 373–400.
6 S. Ramaswamy, Annu. Rev. Condens. Matter Phys., 2010, 1,
323–345.
7 D. L. Koch and G. Subramanian, Annu. Rev. Fluid Mech., 2011,
43, 637–659.
8 T. Ishikawa and T. J. Pedley, J. Fluid Mech., 2007, 588, 399–
435.
9 A. Sokolov and I. S. Aranson, Phys. Rev. Lett., 2009, 103,
148101.
10 S. Rafaı̈, L. Jibuti and P. Peyla, Phys. Rev. Lett., 2010, 104,
098102.
11 G. K. Batchelor, J. Fluid Mech., 1970, 41, 545–570.
12 A. Sokolov, I. S. Aranson, J. O. Kessler and R. E. Goldstein,Phys. Rev. Lett., 2007, 98, 158102.
13 K. Drescher, K. C. Leptos, I. Tuval, T. Ishikawa, T. J. Pedley
and R. E. Goldstein, Phys. Rev. Lett., 2009, 102, 168101.
14 J. Toner, Y. Tu and S. Ramaswamy, Ann. Phys., 2005, 318,
170–244.
15 J. L. Anderson, Annu. Rev. Fluid Mech., 1989, 21, 61–99.
16 W. F. Paxton, S. Sundararajan, T. E. Mallouk and A. Sen,
Angew. Chem., Int. Ed., 2006, 45, 5420–5429.
17 W. F. Paxton, A. Sen and T. E. Mallouk, Chem.–Eur. J., 2005,
11, 6462–6470.
18 J. R. Howse, R. A. L. Jones, A. J. Ryan, T. Gough,R. Vafabakhsh and R. Golestanian, Phys. Rev. Lett., 2007,
99, 048102.
19 H.-R. Jiang, N. Yoshinaga and M. Sano, Phys. Rev. Lett., 2010,
105, 268302.
20 L. Baraban, M. Tasinkevych, M. N. Popescu, S. Sanchez,
S. Dietrich and O. G. Schmidt, So Matter, 2012, 8, 48–52.
21 G. V. Kolmakov, V. V. Yashin, S. P. Levitan and A. C. Balazs,
So Matter, 2011, 7, 3168–3176.
22 W. Yang, V. R. Misko, K. Nelissen, M. Kong andF. M. Peeters, So Matter, 2012, 8, 5175–5179.
23 R. Golestanian, T. B. Liverpool and A. Ajdari, New J. Phys.,2007, 9, 126.
24 J. F. Brady, J. Fluid Mech., 2011, 667, 216–259.
25 M. J. Lighthill, Annu. Rev. Fluid Mech., 1969, 1, 413–446.
26 J. R. Blake, J. Fluid Mech., 1971, 46, 199–208.
27 Y. Nakayama and R. Yamamoto, Phys. Rev. E: Stat.,
Nonlinear, So Matter Phys., 2005, 71, 036707.
28 T. Ishikawa and T. J. Pedley, J. Fluid Mech., 2007, 588, 437–
462.
29 M. J. Kim and K. S. Breuer, Phys. Fluids, 2004, 16, L78.
30 K. Leptos, J. Guasto, J. Gollub, A. Pesci and R. Goldstein,Phys. Rev. Lett., 2009, 103, 198103.
31 T. Ishikawa, J. T. Locsei and T. Pedley, Phys. Rev. E: Stat.,
Nonlinear, So Matter Phys., 2010, 82, 021408.
32 A. Sokolov and I. S. Aranson, Phys. Rev. Lett., 2012, 109,248109.
33 T. Ishikawa, M. P. Simmonds and T. J. Pedley, J. Fluid Mech.,
2006, 568, 119–160.
34 L. Zhu, E. Lauga and L. Brandt, Phys. Fluids, 2012, 24,051902.

Soft Matter, 2013, 9, 4923–4936 | 4935

Soft Matter
35 W. B. Russel, D. A. Saville and W. R. Schowalter, Colloidal
Dispersions, Cambridge University Press, Cambridge, 1stedn, 1992.
36 S. Thutupalli, R. Seemann and S. Herminghaus, New J. Phys.,
2011, 13, 073021.
37 Y. Nakayama, K. Kim and R. Yamamoto, Eur. Phys. J. E: So
Matter Biol. Phys., 2008, 26, 361–368.
38 T. Iwashita, Y. Nakayama and R. Yamamoto, J. Phys. Soc.
Jpn., 2008, 77, 074007.
39 A. Hamid and R. Yamamoto, J. Phys. Soc. Jpn., 2013, 82,
024004.
40 K. Kim, Y. Nakayama and R. Yamamoto, Phys. Rev. Lett.,
2006, 96, 208302.
41 H. Kobayashi and R. Yamamoto, J. Chem. Phys., 2011, 134,
064110.
42 R. Tatsumi and R. Yamamoto, Phys. Rev. E: Stat., Nonlinear,So Matter Phys., 2012, 85, 066704.
43 X. Luo, M. R. Maxey and G. E. Karniadakis, J. Comput. Phys.,2009, 228, 1750–1769.
44 S. Ramachandran, P. B. Sunil Kumar and I. Pagonabarraga,
Eur. Phys. J. E: So Matter Biol. Phys., 2006, 20, 151–
158.

4936 | Soft Matter, 2013, 9, 4923–4936

Paper
45 M. T. Downton and H. Stark, J. Phys.: Condens. Matter, 2009,21, 204101.
46 I. Götze and G. Gompper, Phys. Rev. E: Stat., Nonlinear, So
Matter Phys., 2010, 82, 041921.
47 J. F. Brady and G. Bossis, Annu. Rev. Fluid Mech., 1988, 20,
111–157.
48 J. W. Swan, J. F. Brady, R. S. Moore and ChE 174, Phys. Fluids,
2011, 23, 071901.
49 T. Ishikawa and T. J. Pedley, Phys. Rev. Lett., 2008, 100,
088103.
50 J. V. José and E. J. Saletan, Classical Dynamics: A Contemporary
Approach, Cambridge University Press, New York, 1st edn,
1998.
51 Y. Han, A. M. Alsayed, M. Nobili, J. Zhang, T. C. Lubensky
and A. G. Yodh, Science, 2006, 314, 626–630.
52 J.-P. Hansen and I. R. McDonald, Theory of Simple Liquids,
Academic Press, New York, 3rd edn, 2008.
53 J. Padding and A. Louis, Phys. Rev. Lett., 2004, 93, 220601.
54 D. A. McQuarrie, Statistical Mechanics, University Science
Books, Sausalito, 1st edn, 2000.
55 I. Llopis and I. Pagonabarraga, J. Non-Newtonian Fluid Mech.,
2010, 165, 946–952.

This journal is ª The Royal Society of Chemistry 2013