Smoluchowski Diffusion Equation for Active Brownian Swimmers

0.6 MB·English

by  Francisco J. Sevilla

Smoluchowski Diffusion Equation for Active Brownian Swimmers
Francisco J. Sevilla1,∗ and Mario Sandoval2,†
1Instituto de F´ısica, Universidad Nacional Auto´noma de M´exico,
Apdo. Postal 20-364, 01000, M´exico D.F., Mexico
2Department of Physics, Universidad Autonoma Metropolitana-Iztapalapa, Distrito Federal 09340, Mexico.
(ΩDated: February 3, 2015)
We study the free diffusion in two dimensions of active-Brownian swimmers subject to passive
5
1 fluctuations on the translational motion and to active fluctuations on the rotational one. The
0 Smoluchowski equation is derived from a Langevin-like model of active swimmers, and analytically
2 solved in the long-time regime for arbitrary values of the P´eclet number, this allows us to analyze
theout-of-equilibrium evolution of thepositions distribution of active particles at all time regimes.
n
Explicitexpressionsforthemean-squaredisplacement andforthekurtosisoftheprobabilitydistri-
a
butionfunction are presented,andtheeffectsof persistence discussed. Weshow throughBrownian
J
dynamics simulations that our prescription for the mean-square displacement gives the exact time
1
dependence at all times. The departure of the probability distribution from a Gaussian, measured
3
by the kurtosis, is also analyzed both analytically and computationally. We find that for P´eclet
numbers.0.1, thedistancefrom Gaussian increasesas∼t−2 atshorttimes,whileitdiminishesas
h] ∼t−1 in theasymptotic limit.
c
Keywords: ActiveParticles,DiffusionTheory,Telegrapher’sequation
e
m
- I. INTRODUCTION ral swimmers, and with the purpose of making devices
t
a able to deliver specialized drugs in precise region inside
st The study of self-propelled (active) particles moving our body, or to serve as micromachines able to detect,
t. at small scales has received much attention [1–9]. This anddiagnosediseases[15–17]. These microrobots,inthe
a is so since phenomena in nature like motion of plankton, same way as the smallest microorganisms, are subject
m
viruses andbacteria,have animportant rolein many bi- to thermal fluctuations or Brownian motion, which is
- ological processes as well as in industrial applications, an important issue to take into account, since thermal
d
n and more generally, because active particles are suitable fluctuations make the particles to lose their orientation,
o elementsofanalysisinthecontextonnonequilibriumsta- thus affecting the particles net displacement. Most clas-
c tisticalphysics. Moreover,confinedmattermadeby self- sical literature considering the effect of loss of orienta-
[ propelledparticleshasbeenobservedtobehaveinavery tion due to Brownian motion and activity on the par-

2 different way compared to matter conformed by passive ticles diffusion, has used a Langevin approach, since it
v ones, mainly due to the out-of-equilibrium nature of ac- can be considered as a direct way of finding the parti-
7 tive systems. For example, it has been observed that cles effective diffusion. Within this approach, isotropic
3 activeparticles,confinedinaboxwithsharpcorners,ac- self-propelledbodies subject to thermal forces [5, 18, 19]
2 cumulatemoreatthecornersofthecontainerratherthan and anisotropic swimmers [18] have been treated in the
7 spreading uniformly through the box. This effect origi- absenceandpresenceofexternalfields[20–23]. ASmolu-
0
nates that the generatedpressure on the walls do not be chowski approach to study the effective diffusion of ac-
.
1 uniformalongthecontainer[10]. Incontrast,mattercon- tive particles in slightly complex environments has also
0 formedbypassiveparticlesdevelopauniformpressurein been undertaken. For example, Pedley [24] introduced
5
the container. It has also been reported that interacting a continuum model to calculate the probability density
1
: activeparticleswithdifferentsizesandconfinedinabox, function for a diluted suspension of gyrotactic bacteria.
v formclustersleavingemptyregionsinsidethebox[11],a BearonandPedley[25]studied chemotacticbacteriaun-
Xi verydifferentscenariooccurswithpassiveparticlessince der shear flow and derived an advection-diffusion equa-
these particles would spread homogeneously in the con- tion for cell density. Enculescu and Stark [26] studied
r
a tainer. Thus wenotice thatthe interplayamongactivity the sedimentation of active particles due to gravity and
andconfinementmayprovidenovelpropertiestothelat- subject to translational and rotational diffusion. Simi-
ter systems, hence the necessity of studying them. larly, Pototsky and Stark [27] followed a Smoluchowski
From a technological point of view, the study of ac- approachandanalyzedactiveBrownianparticlesintwo-
tive particles is also very relevant. To mention an ex- dimensional traps. Saintillan and Shelley [28] used a ki-
ample, the bioengineering community is constructing netictheorytostudypatternformationofsuspensionsof
self-propelled micromachines [12–14] inspired by natu- self-propelled particles.
Self-propelledparticles in nature move ingeneralwith
a time-dependent swimming speed, like microorganisms
that tend to relax (rest) for a moment and then to con-
∗ fjsevilla@fisica.unam.mx tinue swimming [29]. In this work, we assume that par-
† [email protected] ticles move with a constant swimming speed, which is
2
a simplified model that has been supported by experi- mined from its own dynamics [32], and Langevin equa-
mental work [20, 30] and has also been adopted when tions foreachmay be written. Here we considerthe case
studying collective behavior [31]. Other more complex of a faster dynamics for the swimming speed such that
scenarios where for example active Brownian particles U (t)=U =const. In this way, the dynamics of this ac-
s 0
are subject to external flows, have been reported and tive Brownian particle, subject to passive-translational
solved following a Langevin approach [23]. and active-rotationalnoise, is determined by its position
In this paper we follow the approach of Smoluchowski x(t) and its direction of motion uˆ(t), computed from
and focus on the coarse-grained probability density of ϕ(t), that obey the following Langevin equations
finding an active Brownianparticle—that diffuses trans-
d
lationally and rotationally in a two-dimensional, un- x(t)=U uˆ(t)+ξ (t), (1a)
0 T
bounded space, and immersed in a steady fluid—at po- dt
sition x at time t without actually knowing its direction d
ϕ(t)=ξ (t). (1b)
R
of motion. We derive Smoluchowski’s equation for such dt
probability density and we solve it analytically.
In this way, the temporal evolution of the particle’s po-
WeareabletousetheFourieranalysistotranslatethe
sition [Eq.(1a)] is thus determined by two independent
Fokker-Planck equation for active particles, into an infi-
stochastic effects, one that corresponds to translational
nite set of coupled ordinary differential equations with
fluctuations due to the environmental noise, the other,
a hierarchy similar to the infinite BBGKY (Bogoliubov-
to the swimmer’s velocity whose orientation is subject
Born-Green-Kirkwood-Yvon) hierarchy that appears in
to active fluctuations [Eq. (1b)]. These same equations
the Boltzmanntheory,whichin principle canbe system-
have been considered in Ref.[33] to model the motion of
aticallysolved. Weclosethehierarchyatthesecondlevel
spherical Platinum-silica Janus particles in a solution of
and explicitly solve the first hierarchy equation, which
water and H O .
2 2
correspondsto the Smoluchowskiequation,andprovides From now on, we use D−1 and U D−1 as time and
Ω 0 Ω
the probability density function (p.d.f) we are interested
length scales respectively, such that D ≡ D D /U2
in. We transform back to the coordinate space, and we B B Ω 0
is the only free, dimensionless parameter of our analysis
areabletoexplicitlyfindthep.d.f. ofanactiveBrownian
which coincides with the inverse of thee so called P´eclet
particle in the short and long time regimes. The mean
number, which in this case corresponds to the ratio of
square displacement (msd) for a swimmer at short and
the advectivetransportcoefficient, U2/D , to the trans-
long times is also found, and classical results for the en- 0 Ω
lationaldiffusiontransportcoefficientD . FortheJanus
hanced diffusion of an active particle are recovered. We B
particles studied by Pallaci et al. in Ref. [34], we have
also find theoretical results for the kurtosis of the swim-
that D−1 ≈ 0.9 s−1, D ≈ 0.34 µm2/s and swimming
mer,henceadiscussionconcerningthenon-Gaussianbe- Ω B
speeds between 0.3 and 3.3 µm/s. These values give D
havior at intermediate times of an active swimmer p.d.f. B
between 0.035 and 4.2. Numerical results, on the other
is also offered. Finally and with comparison purposes,
Browniandynamics simulations were also performed ob- hand, are obtained by integrating Eqs. (1) using tehe
Euler-Crome explicit scheme with a time-step of 0.005.
taininganexcellentagreementamongourtheoreticaland
In Fig. 1 seven single-particle trajectories are presented
computational results at short and long times for both,
kurtosis and mean-square displacement of the swimmer. to show the effect of varying DB. The persistence effects
are conspicuous for small values of D (thick curves).
B
An equation for the one pearticle probability density
II. THE MODEL P(x,ϕ,t)≡hδ(x−x(t))δ(ϕ−ϕ(t))iceanbederivedbydif-
ferentiating this P(x,ϕ,t) with respect to time, namely
Consider a spherical particle of radius a, immersed
∂
in a fluid at fixed temperature T, that self-propels in P(x,ϕ,t)+U0uˆ·∇P(x,ϕ,t)=
∂t
a two-dimensional domain. The particle is subject to
∂
thermal fluctuations (modeled as white noise) in trans- − hξ (t)δ(x−x(t))δ(ϕ−ϕ(t))i
R
∂ϕ
lationand rotation,that is, ξ (t) andξ (t) respectively,
T R
where hξTi = hξRi = 0, hξi,T(t)ξj,T(s)i = 2DBδ(t−s) −∇·hξT(t)δ(x−x(t))δ(ϕ−ϕ(t))i, (2)
and hξ (t)ξ (s)i = 2D δ(t−s). Here D = k T/6πηa
R R Ω B B
where h·i denotes the average over translational and ro-
andD arerespectively,the translationalandrotational
Ω
tationalnoise realizationsand ∇=(∂/∂x,∂/∂y).By us-
diffusivities constants, with η the viscosity of the fluid.
The particle swimming velocity, U (t), is written ing Novikov’s theorem, we obtain the following Fokker-
s Planck equation for an active Brownian particle
explicitly as U (t)uˆ(t), where we denote by uˆ(t) =
s
[cosϕ(t),sinϕ(t)] (ϕ(t) being the angle between the di-
∂
rection of motion and the horizontal axis) the instan- P(x,ϕ,t)+U uˆ·∇P(x,ϕ,t)=
0
∂t
taneous unit vector in the direction of swimming, and
∂2
lUosc(itt)ytahloenignsuˆta(tn)t.anEeaocuhsomftahgensietuqdueanotfittihees smwaiymbmeindgetveer-- DB∇2P(x,ϕ,t)+DΩ∂ϕ2P(x,ϕ,t). (3)
3
D~~B = 0.001 ficient of the expansion pn(k,t), namely
D = 0.00316
~B
DDD~~BB === 000...0011316 ddtpn =−U20ike−DΩtee2nDΩte−iθpn−1+
DD~~BB == 01..3016 e (cid:2) e−2enDΩteiθpn+1 . (8)
B
0 Note that we have introduced the quantitiesekx±(cid:3)iky =
0 ke±iθ.
OurinterestliesonP (x,t)= 2πdϕP(x,ϕ,t),which
0 0
is given by the inverse Fourier transform of P (k,t) =
R 0
e−DBk2tp (k,t) as can be checked from Eq. (7). Let
0
us find p by solving Eq. (8) when the initialecondition
0
Pn(k,0)e= pn(k,0) = δn,0/2π is considered. This ini-
tial condeition corresponds to the case when the particle
deeparts fromethe origin in a random direction of motion
drawnfromauniformdistributionin[0,2π)andisequiv-
alent to the initial condition P(x,ϕ,0) = δ(2)(x)/2π,
FIG.1. (Coloronline)Singleparticletrajectoriesfordifferent
values of DeB = 0.001, 0.00316, 0.01, 0.0316, 0.1, 0.316, 1.0. where δn,m and δ(2)(x) denote, respectively, the Kro-
Data is generated by solving Eq. (1) during 104 time steps, necker delta and the delta function.
WefirstdealwiththesolutionofEq.(8)forlongtimes.
the initial position is chosen at the origin while the initial
To do so, we consider only the first three Fourier coeffi-
orientation is drawn from a uniform distribution in [0,2π).
cients p , i.e., those with n = 0, ±1, later on it will be
n
shownthattheothercoefficientsdonotcontributeinthis
regime. From Eq. (8) the three first Fourier coefficients
We now start to analytically solve Eq. (3). To do so, we e
satisfy
apply the Fourier transform to Eq. (3) and we obtain
d U
p =− 0ike−DΩt e−iθp +eiθp , (9a)
∂ 0 −1 1
P(k,ϕ,t)+iU uˆ·kP(k,ϕ,t)= dt 2
0
∂t e −D k2P(k,eϕ,t)+D ∂2 P(k,ϕ,t), (4) ddtp±1 =−eU20ik eDΩte∓iθp0+(cid:2) e−3DeΩte±iθp±e2(cid:3), (9b)
B Ω∂ϕ2 (cid:2) (cid:3)
whicheaftersomealgebraareceombinedintoaesingleequa-
e e tion for p , that is
where 0
d2 d U2
P(k,ϕ,t)=(2π)−1 d2xeik·xP(x,ϕ,t), (5) p e+D p =− 0k2p
dt2 0 Ωdt 0 2 0
Z
U2
and k =e(k ,k ) denotes the Fourier wave-vector. For e − e0k2e−4DΩtee2iθp2+e−2iθp−2 . (10)
x y 4
U = 0, Eq. (4) has as solution the set of eigenfunc-
0 (cid:0) (cid:1)
tions {e−(DBk2+DΩn2)teinϕ} with n an integer. Thus we At this stage we argue that in theelong timeeregime we
mayneglect the secondterm inthe rhs of Eq.(10). This
expand P(k,ϕ,t) on this set, expressly
approximationleads to the telegrapher’s equation
∞
P(k,ϕ,et)= 21π pn(k,t)e−(DBk2+DΩn2)teinϕ, (6) ddt22p0+DΩddtp0 =−U202k2p0, (11)
n=−∞
X
e e whichisrotationaellysymmetericinthek-sepaceandthere-
with k = |k|. This expansion allows to separate the ef-
fore, giving rise to rotationally symmetric solutions in
fects of translationaldiffusion contained in the prefactor
spatial coordinates if initial conditions with the same
e−DBk2t from the effects of rotational diffusion due to symmetry are chosen. One may check that Eq. (11) has
active fluctuations. the solution
The coefficients of the expansion are given by
D
p (k,t)=p (k,0)e−DΩt/2 Ω sinω t+cosω t ,
0 0 k k
2π 2ω
pn(k,t)=e(DBk2+DΩn2)t dϕP(k,ϕ,t)e−inϕ. (7) (cid:20) k (cid:21)(12)
Z0 witehω2 ≡U2ek2/2−D2/4. Therefore,byusingthisresult
k 0 Ω
AfteersubstitutingEq.(6)inEq.(4)eandusingthorthog- one can explicitly write that P0(k,t) = e−DBk2tp0(k,t),
onality of the eigenfunctions we get the following set of with e−DBk2t being the translational diffusion propa-
coupled ordinary differential equations for the n-th coef- gator. The solution in spatiael coordinates, P0(ex,t), is
4
obtained after taking the inverse Fourier transform of
P (k,t),whichisgivenbytheconvolutionofthetransla-
0
tional propagator with the probability distribution that
reetain the effects of active diffusion, i.e.,
P (x,t)= d2x′G(x−x′,t)p (x′,t), (13) FIG. 2. (Color online) Snapshots of the positions of 105
0 0
particles taken at different times, namely from left to right,
Z
DΩt=0.01, 0.1, 1.0, 10, and 100, for an inverse P´eclet num-
where p (x,t)is the inverseFourier transformofp (k,t) e
0 0 ber DB =0.001. Notice theformation of a rotationally sym-
and
metricring-likestructure(secondpanelfromleft)thatdevel-
e−x2/4DBt e ops from a Gaussian distribution (far left panel) due to the
G(x,t)= , (14) effectsofpersistenceinthemotionoftheparticles. Thestruc-
4πD t ture fades out as time passes (third and fourth panels from
B
left) toreach a Gaussian distribution at long times (far right
is the Gaussian propagator of translational diffusion in panel) .
two dimensions.
In the asymptotic limit (t → ∞), at which the coher-
ent wave-like behavior related to the second order time
formationof a rotationally symmetric ring-like structure
derivativeinEq.(11)canbeneglected(mainlyduetothe
at the intermediate time D t = 0.1 (second panel from
Ω
random dispersion of the particles direction of motion),
left), mainly due to the persistence effects on the parti-
p (x,t) tends to a Gaussian distribution with diffusion
0 cles’ motion. The structure fades out with time (third
constant U2/2D , namely
0 Ω and fourth panels from left) and starts turning into the
Gaussiandistribution[Eq.(16)],astimepasses(farright
e−x2/4(U02/2DΩ)t panel for D t=100). We must mention that in spite of
p (x,t)−−−→ . (15) Ω
0 t→∞ 4π(U02/2DΩ)t the good agreement in the short and asymptotic limits
givenby (18),thissolutiondoesnotprovideanadequate
SubstitutionofthislastexpressionintoEq.(13)andper- description of the particles distribution in the interme-
forming the integral, one gets in the long-time regime
diate time regime for all values of D , particularly for
B
the small ones, for which the effects of persistence are
P (x,t)= e−x2/4(DB+U02/2DΩ)t, (16) conspicuously prominent. The reasone of this failure lies,
0 4π(DB +U02/2DΩ)t not in the generalsolution (13) but in the approximated
equation(11),whichcorrespondstothe twodimensional
from which the classicaleffective diffusion constant D =
telegrapher’s. The solution to this equation [(12)] can
D +U2/2D is deduced [34].
B 0 Ω not be interpreted as a p.d.f. since becomes negative at
In the short time regime D t ≪ 1, we approximate
Ω shorttimes [35]. Thisundesiredfeature,resultsfromthe
p (k,t) ≈ p (k,0) = (2π)−1 and the p.d.f. in spatial
0 0 wake effect, characteristic of the solution of the two di-
coordinates corresponds to the Gaussian
mensional wave equation [9]. Notwithstanding this and
e e 1 e−x2/4DBt asisshownlateron,the mean-squaredisplacementcom-
P (x,t)≈ . (17) puted from our approximation coincides with the exact
0
2π 4πDBt one computed from numerical simulations, at all time
regimes.
Though it is a difficult task to obtain from Eq. (13) the
explicit dependence onx and t of P , a formula in terms
0
ofaseriesexpansionoftheoperator∇2appliedtoG(x,t)
can be derived, to say III. MEAN-SQUARE DISPLACEMENT
∞
1 1 D t 1 What are the effects of self-propulsionon the diffusive
P (x,t)= e−DΩt/2 1+ Ω ×
0 2π (2s)! 2 2s+1 behavior of the system is a question that has been ad-
Xs=0 (cid:18) (cid:19) dressed in several experimental and theoretical studies,
DΩt 2s 1+ 2U02∇2 s G(x,t), (18) and the msd, being a measure of the covered space as
2 D2 a function of time by the random particle, has been of
(cid:18) (cid:19) (cid:18) Ω (cid:19)
physical relevance in both contexts. Indeed, the msd is
Note that the last formula can be rewritten in terms of obtained in many experimental situations that consider
∂t instead of ∇2 by using that ∇2G=DB−1∂tG. active particles [18, 33, 36], hence a way of validating
In figure 2 snapshots of the positions of 105 particles existing theoretical approaches. On the other hand, it is
taken at different times from numerical simulations, and knownthatthediffusioncoefficientofactiveparticlesde-
for the inverse P´eclet number D = 0.001, are shown. pendsontheparticledensity,duemainlytoexcludedvol-
B
The first panel from left (D t = 0.01), corresponds to a ume effects of the particles. Here we assume an enough
Ω
p.d.f. close to the Gaussian giveen by Eq. (17). Note the diluted system to neglect those effects.
5
In what follows, an exact analytical expression for the >
2 ~t
msd will be obtained. From Eq. (11) one can show that ) 104
>
in the coordinate space P (x,t) satisfies
0 x
2
< 10
∂2 ∂ U2 -
P +D P = 0 +D D ∇2P x 0
∂t2 0 Ω∂t 0 2 B Ω 0 ( 10 ~
(cid:18) (cid:19) < D
+DB(cid:18)2∂∂t −DB∇2(cid:19)∇2P0. (19) 1100--42 B ~t 2 DDDD~~~~BBB ==== 0000....00000013131166
Iwfhlaosletesqpuacaet,iowneisombtualitnipfloierdthbeymx2sdanthienteeqguraattieodnoverthe 10-6 ~t DDD~~~BBB ===0 10..1.031623
B
d2 hx2(t)i+D dhx2(t)i=2U2+4D D , (20) 10-4 10-2 100 102 104
dt2 Ωdt 0 B Ω D t
Ω
whose solution can be easily found, namely
FIG. 3. (Color online) Mean-square displacement in units of
hx2(t)i=4DB 1−e−DΩt U02/DΩ2 asfunctionofthedimensionelesstimeDΩtfordifferent
DΩ valuesoftheinverseP´ecletnumberDB,namely0.001,0.0031,
2 2D D(cid:0) +U2 (cid:1) 0.01, 0.0316, 0.1, 0.3162, 1.0. Solid lines correspond to the
+ B Ω 0 D t− 1−e−DΩt , (21) analytical expression given byEq.(21),while thethin,dash-
(cid:0) DΩ2 (cid:1) Ω dottedlinescorrespondtotheresultsobtainedfromnumerical
(cid:2) (cid:0) (cid:1)(cid:3) simulations (see text).
where we have used that (d/dt)hx2(t)i = 4D in
t=0 B
this case.
(cid:12)
The linear dependence on time of the(cid:12)msd, expected We now use that
intheGaussianregime,ischeckedstraightforwardlyfrom
1 ∂ ∂
Eq. (21). In the long-time regime we get x2(t) =−∇2 P (k,t) =− k P (k,t) ,
k 0 0
k=0 k∂k ∂k k=0
(cid:12) (cid:12) (25)
hx2(t)i−−−−−→4 D + U02 t, (22) h(cid:10)ence,(cid:11)at short teimes, t(cid:12)(cid:12)he msd of the activee parti(cid:12)(cid:12)cle is
B
DΩt→∞ (cid:18) 2DΩ(cid:19) finally given by
which is a classical result indicating that the effective x2(t) ≈4D t+ 2D D +U2 t2, (26)
diffusion of a self-propelled particle is enhanced by its B B Ω 0
activity [3]. In the opposite limit. whichsho(cid:10)wsline(cid:11)arandquadr(cid:0)aticbehavioro(cid:1)fthemsd. In
the limit t →0, the above equation reduces to Eq. (23).
x2(t) −−−−−→4D t, (23)
B Thus wehavearrivedto the same resultsby two distinct
DΩt→0
methods, namely, one method that obtained the explicit
(cid:10) (cid:11)
which once again the latter equation shows, the linear expressionfortheprobabilitydensityandusedit,andthe
behavior with time of the msd. secondmethodthatextractedinformationfromEq.(19)
without solving it.
In order to validate our theoretical findings, we have
A. Comparison of the mean-square displacement also performed Brownian dynamics simulations. Fig. 3
showsacomparisonamongourtheoreticalpredictionfor
any time given by Eq. (21) (solid lines), and our Brown-
In this subsection we compare the results obtained for
ian dynamics simulations (dashed lines) implemented to
the msd at short and long times using Eq. (19) with re-
solve Eqs. (1) and averaging over 105 trajectories. Each
sults obtained following our new Fourier approach.
For long times and with the explicit expression for plottedlinecorrespondstodifferentvaluesofDB. Weob-
P (x,t) [Eq. (16)], we use the definition of the mean- serve and excellent agreement among theory and Brow-
0
squaredisplacement,hx2(t)i= d2xx2P (x,t)toeasily nian simulations. Additionally, Fig. 3 indiceates that a
0
transition of the msd from a linear to a quadratic be-
recover Eq. (22). For short times, we expand P (k,t) in
R 0
powers of k, and we keep only the quadratic terms in k havior is not always the case. For large values of DB
for the reasonthat will be clear immediately afeterwards, the linearbehaviorovertime dominates,whereasforDB
we have smaller than 0.3, a transition starts to occur. The aeb-
senceofatransitionfromalineartoaquadraticbehaveior
P (k,t)=1+ DΩ −D k2 t− DBDΩk2 + ωk2 t2. has also been reported by ten Hagen et al. [18]. They
0 2 B 2 2 suggest that at short times, the msd of a self-propelled
(cid:18) (cid:19) (cid:18) (cid:19)
(24) particleisaffectedbytheinitialorientationandthemag-
e
6
nitude ofitsforcepropulsionthatmakethatatransition where xT denotes the transpose of the vector x and Σ
occur or not. is the 2×2 matrix defined by the average of the dyadic
product (x − hxi)T · (x − hxi). In addition, it can be
shownthatfor circularlysymmetric distributions, as the
onesconsideredinthepresentstudy,Eq. (27)reducesto
hx4(t)i
IV. KURTOSIS κ=4 r, (28)
hx2(t)i2
r
Once we have obtained approximately P (x,t) from here h·i stands for the radial average over the ra-
0 r
Eq. (3), we wish to characterize its departure from a dial probability density distribution, namely h·i =
rad
∞
Gaussian p.d.f. as a function of time. This quantity drrP(r)(·). The analytical results obtained for the
0
has been measured experimentally is systems of Janus kurtosis will also be compared with those obtained from
R
particlesTothispurposewecalculateitskurtosisκ,given numerical simulations.
explicitly by [37] As we can see, Eq. (28) requires the calculation of the
fourth moment. After certain algebraic steps (see ap-
pendix B), we explicitly find that the fourth moment is
κ= (x−hxi)TΣ−1(x−hxi) 2 , (27) given by
D(cid:2) (cid:3) E
U4 D t D t 2 D t
hx4(t)i =25 0 e−DΩt/2 −3 Ω + Ω 1+4D (1+D ) cosh Ω +
r D4 2 2 B B 2
Ω (" (cid:18) (cid:19) (cid:16) (cid:17)# (cid:18) (cid:19)
D t e eD t 2 D t
Ω Ω Ω
3− (1+4D )+ (1+4D (1+D )) sinh . (29)
B B B
2 2 2
" (cid:18) (cid:19) # (cid:18) (cid:19))
e e e
The latter analytical expression is used in Fig. 4 to plot the values of κ corresponding to the snapshots shown in
the kurtosis, to be precise, the ratio of Eq. (29) to the Fig. 2). Inthis regimethe orientationalcorrelationssur-
square ofEq. (21). In this figure, kurtosis is plotted as a pass the effects of translational noise hence allowing the
function of the dimensionless time D t for different val- particles to move persistently outwardly from the origin
Ω
ues of D (dashed lines), while the corresponding exact and giving rise to circularly-symmetric, ring-like distri-
B
results from Brownian dynamics simulations are shown butions as the one shown in Fig. 2 (second from left)
in symbeols. Note that for these last ones, 4≤κ≤8. at D t = 0.1 for D = 0.001. Afterwards, this ring-
Ω B
We observe that the kurtosis shows a clear non- like distributions turns into a Gaussian distribution as
monotonic behavior for DB . 1, namely, it starts to the orientational coerrelations die out. Notice that the
diminish from κ = 8 as time passes, until it reaches a discrepancy between the results given by our analytical
minimum at a time arouned D−1 that depends on the in- approximationandtheexactresultsliesonthe wave-like
Ω
verse of the P´eclet number D . This behavior has been effects given by the second-order-time-derivativeterm of
B
observedinexperimentswith Janusparticlesinthreedi- Eq. (11) which, as shown in [9], gives 8/3 ≃ 2.6667 as
mensions [33] and theoreticaelly predicted in one quasi- the value for the kurtosis in the vanishing translational
one-dimensional channels [38]. In the limit of vanishing diffusion. For DB = 0.001, the minimum of the exact
translational diffusion, the kurtosis is a monotonic in- kurtosis is about 4.1352 at t ≈ 0.39D−1 while from our
Ω
creasing function of time, with κ = 4 as its minimum analytical approeximate κ≈3.0184 at t≈0.239D−1.
Ω
value [9]. At later times the kurtosis starts to increase In the next subsections we find asymptotic expression
reaching the value 8 in the asymptotic limit. This be- forthekurtosisforbothshortandlongtimesthatconfirm
havior becomes less and less evident as the translational theGaussianityofthedistributionfunctionforthesetime
diffusion coefficient surpassesthe effective diffusion coef- regimes.
ficientthatoriginatesintherotationaldiffusion, U2/D .
0 Ω
Due to this effect, we note a better agreement between
our analyticalapproximationandour simulationresults,
as D is increased. On the other hand, the persistence A. Kurtosis at short times
B
effects induced by orientationalcorrelations are revealed
in tehe diminishing regime of the kurtosis, as is shown in In ordet to find κ for t → 0, we expand P (k,t) in a
0
Fig. 4 for the cases in which D . 1 (see Table I for powerseriesink,usethefactthatω2 ≡U2k2/2−D2/4,
B k 0 Ω
e
e
7
B. Kurtosis at long times
TABLEI.KurtosisvaluesattimesDΩtthatareclosetothose
e
values chosen in Fig. 2 (DB =0.001).
In ordet to find κ for t → ∞, we proceed as before,
DΩt 0.01 0.1 1.01 10.245 100.08 we expand PP (k,t) in power series of k, use that ω2 ≡
κ 5.9598 4.3111 4.2729 6.7621 7.8713 0 k
U2k2/2−D2/4, but now keep all the terms in the series
0 Ω
expansion. Aefter applying Eq. (B1) we get that in the
long time limit
κ 9
25U4 4D D
hx4(t)i≈ 0 1+ B Ω
8 D4 U2
Ω (cid:20) 0
2 2
7 D~ +4 DBDΩ DΩt . (32)
B U2 2
(cid:18) 0 (cid:19) #(cid:18) (cid:19)
6
D = 0.001
B Finally we use Eq. (22) to get for long times that
D = 0.00316
B
5 D = 0.01
DB = 0.03162 hx4(t)i
4 DDBB == 00..131623 κ= 4 D + U02 t =8, (33)
DB = 1.0 B 2DΩ
B
3 (cid:16) (cid:17)
which once again, Eq. (33) shows that kurtosis has a
Gaussian behavior for long times. This result has also
2
-4 -2 0 2 4 been validated using Brownian dynamics and shown in
10 10 10 10 10
Fig.4,whereforlongvaluesofthe dimensionlesstime, κ
D t
Ω tends again to a value of 8.
FIG. 4. (Color online) Kurtosis of the particles distribution
V. THE ACTIVE FLUX: THE VELOCITY
as functeion of the dimensionless time DΩt for different val- FIELD V(x,t) OF ACTIVE PARTICLES
ues of DB, namely 0.001, 0.0031, 0.01, 0.0316, 0.1, 0.3162,
1.0. Dashed lines correspond to the exact analytical expres-
sion given by the quotient of Eq. (29) and Eq. (21), while The Fourier expansion given by Eq. (6) allows us to
the symbols mark the values of the kurtosis obtained from obtain the hydrodynamic fields. At the moment we
numerical simulations. Thin-dotted lines mark the values 8, have explicitly calculated P (k,t) that corresponds to
0
4, 8/3, that correspond respectively in two dimensions to: a
the Fourier transform of the density field ρ(x,t). The
Gaussian distribution, a ring-like distribution and wave-like
Fouriertransformof the comeponents of the velocity field
distribution [9].
V(x,t) are related to p (k,t) in the following way. By
±1
definition
e2π
and keep terms only of order k4 to get V (k,t)=U dϕ cosϕP(k,ϕ,t), (34a)
x 0
Z0
2π
D c4t5 c4t4 D c2D t4 Ve (k,t)=U dϕ sinϕPe(k,ϕ,t), (34b)
P (k,t)≈k4 Ω + + Ω B y 0
0 2 5! 4! 2 3! Z0
(cid:20)
e +c2DBt3 + DΩDB2t3 + DB2t2 , and afteer substitution of Eq. (6) ineEq. (34) we get
2 4 2 (cid:21) V (k,t)=U0 e−(DBk2+DΩ)t [p (k,t)+p∗(−k,t)],
=k4f(t). (30) x 2 1 1
(35a)
e e e
We now use Eq. (B1), that together with Eq. (26), one Vy(k,t)=U20 e−(DBk2+DΩ)t [p∗1(−k,t)−p1(k,t)],
finally gets for short times that
(35b)
e e e
where we have used that p (k,t) = p∗(−k,t). With
43f(t) −n n
κ≈ =8, (31) these expressions,it is straightforwardto verify that Eq.
(4DBt+2DBDΩt2+U02t2)2 (9a) corresponds to the coentinuity equaetion ∂tρ + ∇ ·
V(x,t)=0.Atthelongtimeregime,whereonlythefirst
threeFouriermodesareneeded,wehavethatp (k,t)can
which means that our derived p.d.f. has a Gaussian be- 1
be computed from Eq. (12) if we neglect the term that
havior at short times. This result has been validated
involves p , i.e., from
numerically and illustrated in Fig. 4, where kurtosis is 2 e
plotted versus time. It clearly can be seen that for short d U
times κ has a value of 8. p1 =− 0ike−iθeDΩtp0. (36)
dt 2
e e
8
Wenowusetheexplicitexpressionforp0,thatisEq.(12), κ 101
with p (k,0)=1/(2π)to integrateEq. (36). After some -
0 8 ~ t-1
steps we find that e 0 ~ t2
10
e ~
p (k,t)=−U0 (ik +k )eDΩt/2sinωkt. (37) DB
1 x y
4π ωk 10-1
If we iensert the latter into Eqs. (35a)-(35b) we obtain ~
D = 0.001
V (k,t)=−U02ik e−(DBk2+DΩ/2)tsinωkt, 10-2 DD~~~BBB == 00..0001316
x 4π x ωk D~B = 0.0316
-3 D = 0.1
Ve (k,t)=−U02ik e−(DBk2+DΩ/2)tsinωkt. 10 D~BB = 0.316
y y
4π ω
k
-4 -2 0 2 4
10 10 10 10 10
AfterinevertingtheFouriertransform,thispairofexpres-
D t
sions can be brought into the form
Ω
U2
V(x,t)=− 0∇F(x,t), (38) FIG.5. (Coloronline)ThedeparturefromGaussianbehavior
2 ofP0(x,t),measuredby8−κ,differsintheshort-timeregime
from the long-time regime. In the first case the departure
where the function F(x,t) depends on P0, explicitly captured by the power law (DΩt)−2 while in the long-time
regime such departureis well described simply by (DΩt)−1.
t
F(x,t)= dseDΩ(t−s) d2x′G(x−x′,t−s)P (x′,s).
0
Z0 Z
(39)
Last expression corresponds to a non-Fickian constitu-
tive relation that considers nonlocal effect in space and In summary, we have found a method and obtained
time. This relation, together with the continuity equa- withit,ananalyticalsolutiontotheSmoluchowskiequa-
tion, gives rise to a non-local (in both, space and time) tionofanactiveBrownianparticle self-propellingwith a
diffusion-like equation. constant velocity. We used a Fourier approach and ex-
ploited the circular symmetry of the distribution func-
In a similar way, one can systematically find higher
tion,toobtainaninfinitesystemofcoupledordinarydif-
Fourier modes for the probability density function,
ferentialequationsfortheFouriermodesoftheprobabil-
P(k,ϕ,t), in the wave vector space.
itydensity. Ourformalismshowedthatinordertohavea
wholedescriptionfortheparticlesdiffusion,weonlyneed
e
thefirstFouriercoefficient,sinceitwasfoundthathigher
VI. FINAL COMMENTS AND CONCLUSIONS
Fouriermodes donotcontribute to the mean-squaredis-
placement of the particle. With the explicit p.d.f. in
Though Eq. (21) has been derived from the long time hand, we calculated the particle mean-square displace-
approximation, it does give the correct time dependence ment for both long and short times, recovering classical
at all times as it is shown in Fig. 3, where our analyt- resultsofenhanceddiffusionduetoactivity. Wealsoval-
ical formula (solid lines) is compared against numerical idated these findings by performing Brownian dynamics
simulations (dashed lines). This fact can be understood simulations that showed and excellent agreementamong
directly from Eq. (10), since it can be shown that the theory and simulations. In addition, we also performed
inhomogeneous term does not contribute to the mean an asymptotic analysis for short and long times, to the
squared displacement (see appendix A), therefore agree- forth moment of our p.d.f. This calculation enabled us
ing with the expression obtained from the telegrapher’s to have analytical results for the kurtosis. The asymp-
Eq. (11). totic results for this measure were validated with Brow-
We want to point out that though the linear time de- niandynamicssimulationsthatshowedthatthe p.d.f. of
pendenceofthemsd,characteristicofnormaldiffusion,is an active particle always starts with a Gaussian behav-
reached at times DΩt∼1 (see Fig. 3), the Gaussian be- ior, and depending on the vale of the ratio DB, it may
havior of the distribution is reached at much later times evolve to a ring-like distribution, but it always returns
(see Fig. 4), suggesting that at finite times, DΩt > 1, to a Gaussian distribution for long times. Feinally this
normal diffusion does strictly corresponds to Gaussian work shed light on mathematical procedures to solve a
distributions. The kurtosis computed from our approxi- Smoluchowski equation for an active Brownian particle.
mation shows that it tends asymptotically to Gaussian,
κ = 8, as (8−D t)−1. In contrast, in the short time Future venues for this research, would be to solve an-
Ω
regime, the distribution departs from a Gaussian behav- alytically a Smoluchowski equation for confined active
ior with time as (8−D t)−2 (see Fig. 5). particles and/or particles interacting among them.
Ω
9
ACKNOWLEDGMENTS and therefore, using Eqs. (A1) we finally have that
d2 d
F.J.S. acknowledges support from DGAPA-UNAM dt2hx2(t)i+DΩdthx2(t)i0 =2U02+4DBDΩ, (A4)
through Grant No. PAPIIT-IN113114. M. S. thanks
CONACyT and Programa de Mejoramiento de Profeso- which coincides with Eq. (20), deduced from the tele-
rado (PROMEP) for partially funding this work. grapher Eq. (11) when considering only the first three
Fourier modes p ,p .
0 ±1
e e
Appendix B: Calculation of the Fourth moment of
Appendix A: Mean Square Displacement P0(x,t)
Ifweapplytheoperator−∇2 totherelationP (k,t)= The fourth moment hx4(t)i can be obtained through
k 0 r
e−DBk2tp (k,t) and evaluate at k = 0, one gets that the prescription
0
hx2(t)i=4D t+hx2(t)i ,whichleadsstraighteforwardly
to the paeir ofBrelations 0 hx4(t)i= 1 ∂ k ∂ 2 P (k,t) . (B1)
0
k∂k ∂k k=0
(cid:18) (cid:19) (cid:12)
d2 hx2(t)i= d2 hx2(t)i (A1a) Though it is possible to use this preescript(cid:12)(cid:12)ion directly on
dt2 dt2 0 P (k,t), we prefer to apply it to the series expansion in
0
dhx2(t)i=4D + d hx2(t)i , (A1b) powers of k of P0(k,t), this is due to the complicated
B 0
dt dt deependence of P on k. Such expansion is given by
0
e
pwcorhomebrpaeubhti·eliidt0yddimerneecoattsleyusrbteyh.eapTapvhleeyrianqgguea−onf∇ti(t2·y)twoddtiE22thqhx.p(201((t0x))i,0atn)cdaasnevtabhlee- P0(k,t)=e−eDΩt/2n,Xm∞=0Xjm=0(cid:18)mj(cid:19)(−DnB!t)n((i2tm)2m)!
k
uating at k=0, which results in e D t D2 m−j U2 j
× 1+ Ω − Ω 0 k2(n+j).
2(2m+1) 4 2
(cid:18) (cid:19)(cid:18) (cid:19) (cid:18) (cid:19)
(B2)
d2 d U2
dt2hx2(t)i0+DΩdthx2(t)i0 = 20 ∇2k(kx2+ky2)p0 k=0 where we have used the explicit dependence on k, of ω ,
k
U2 (cid:2) (cid:3) given in section II.
+ 40e−4DΩt ∇2k (kx−iky)2p−2+(kx+iky)2p2 ek=0 Since
(cid:8) (cid:2) (cid:3)(cid:9) (A2)
e e 1 ∂ ∂ 2
k k2(n+j) =[2(n+j)]2×
It is easy to show that the first term of rhs of the last k∂k ∂k
(cid:18) (cid:19)
equation is 2U02 while the next term, that involves the [2(n+j−1)]2k2(n+j−2), (B3)
modes p vanishes identically, thus we get
±2
the only terms that survives when evaluating at k = 0
e d2 d are those for which n+j = 2 with the condition that
hx2(t)i +D hx2(t)i =2U2, (A3) 0≤j ≤m. Thus we can write
dt2 0 Ωdt 0 0
U4 ∞ D t 2m+2 D t 2 D t 1 (m+2)(m+1)
hx4(t)i=44 0 e−DΩt/2 Ω Ω 1+ Ω +
D4 2 2 2 2m+5 2(2m+4)!
Ω m=0(cid:18) (cid:19) "(cid:18) (cid:19) (cid:18) (cid:19)
X
D t D t 1 m+1 D t 1 1
D Ω 1+ Ω +D2 1+ Ω (B4)
B 2 2 2m+3 (2m+2)! B 2 2m+1 2(2m)!
(cid:18) (cid:19)(cid:18) (cid:19) (cid:18) (cid:19) (cid:21)
e e
Thesumscanbewrittenintermsofhyperbolicfunction, and after some algebraic steps Eq. (29) is recovered.
[1] P.S.LovelyandF.W.Dahlquist,J.Theor.Biol.50,477 [2] T. J. Pedley and J. O. Kessler, Annu. Rev. Fluid Mech.
(1975).
10
24, 313 (1992). E. 84, 031105 (2011).
[3] H.C.Berg, Random walks in biology (Princeton Univer- [23] M. Sandoval, N. K. Marath, G. Subramanian, and
sity Press, Princeton, N.J., 1993). E. Lauga, J. Fluid Mech. 742, 50 (2014).
[4] T. Ishikawa and T. Pedley, J. Fluid Mech. 588, 437 [24] T. J. Pedley and J. O. Kessler, Journal of Fluid
(2007). Mechanics 212, 155 (1990), ISSN 1469-7645,
[5] J. R. Howse, R. A. L. Jones, A. J. Ryan, T. Gough, URL http://journals.cambridge.org/article_
R.Vafabakhsh,andR.Golestanian,Phys.Rev.Lett.99, S0022112090001914.
048102 (2007). [25] R.N.BearonandT.J.Pedley,Bull.Math.Biol.62,775
[6] P. Romanczuk, M. Bar, W. Ebeling, B. Lindner, and (2000).
L. Schimansky-Geier, Eur. Phys. J. Special Topics 202, [26] M. Enculescu and H. Stark, Phys. Rev. Lett. 107,
1 (2012). 058301 (2011), URL http://link.aps.org/doi/10.
[7] E. Lauga and T. Powers, Rep. Prog. Phys. 72, 096601 1103/PhysRevLett.107.058301.
(2009). [27] A. Pototsky and H.Stark, EPL 98, 50004 (2012).
[8] F.Ju¨licherandJ.Prost,TheEuropeanPhysicalJournal [28] D. Saintillan and M. J. Shelley, Phys. Rev. Lett.
E 29, 27 (2009). 100, 178103 (2008), URL http://link.aps.org/doi/
[9] F. J. Sevilla and L. A. Gomez Nava, Physical Review E 10.1103/PhysRevLett.100.178103.
90, 022130 (2014). [29] S. Babel, B. ten Hagen, and H. Lowen, J. Stat. Mech.
[10] Y.Fily,A. Baskaran, and M. F. Hagan, Soft Matter 10, P02011, 1 (2014).
5609 (2014). [30] S. Bazazia, J. Buhl, J. J. Hale, M. L. Anstey, G. A.
[11] X.Yang,M.L.Manning,andM.C.Marchetti,SoftMat- Sword, S. J. Simpson, and I. D. Couzin, Curr. Biol. 18,
ter10, 6477 (2014). 735 (2008).
[12] J.J.Abbott,K.E.Peyer,M.C.Lagomarsino, L.Zhang, [31] T. Vicsek, A. Czirok, E. Ben-Jacob, I. Cohen, and
L. Dong, I. K. Kaliakatsos, and B. J. Nelson, Int. J. O. Shochet, Phys.Rev.Lett. 75, 1226 (1995).
Robot. Res.28, 1434 (2009). [32] P. Romanczuk and L. Schimansky-Geier, Phys. Rev.
[13] G. Kosa, P. Jakab, G. Szekely, and N. Hata, Biomed. Lett. 106, 230601 (2011), URL http://link.aps.org/
Microdevices 14, 165 (2012). doi/10.1103/PhysRevLett.106.230601.
[14] W. Gao, R. Dong, S. Thamphiwatana, J. Li, W. Gao, [33] X. Zheng, B. ten Hagen, A. Kaiser, M. Wu, H. Cui,
L. Zhang, and J. Wang, ACS nano (2014). Z.Silber-Li,andH.Lowen,PhysicalReviewE88,032304
[15] T. E. Mallouk and A. Sen,Sci. Am. 300, 72 (2009). (2013).
[16] W. F. Paxton, S. Sundararajan, T. E. Mallouk, and [34] J.Palacci,C.Cottin-Bizonne,C.Ybert,andL.Bocquet,
A.Sen, Angew. Chem., Int.Ed. 45, 5420 (2006). Phys. Rev.Lett. 105, 088304 (2010).
[17] T. Mirkovic, N. S. Zacharia, G. D. Scholes, and G. A. [35] J.M.Porra,J.Masoliver,andG.H.Weiss,Phys.Rev.E
Ozin,ACS Nano4, 1782 (2010). 55, 7771 (1997), URL http://link.aps.org/doi/10.
[18] B. ten Hagen, S. van Teeffelen, and H. Lowen, J. Phys. 1103/PhysRevE.55.7771.
Condens. Matter 23, 194119 (2011). [36] D. Selmeczi, S. Mosler, P. H. Hagedorn, N. B. Larsen,
[19] V. Lobaskin, D. Lobaskin, and I. Kulic, Eur. Phys. J. and H.Flyvbjerg, Biophys. J. 89, 912 (2005).
Special Topics. 157, 149 (2008). [37] K.V.Mardia,Sankhya¯: TheIndianJournalofStatistics,
[20] S. Ebbens, R. A. L. Jones, A. J. Ryan, R. Golestanian, Series B pp.115 – 128 (1974).
and J. R. Howse, Phys.Rev.E 82, 015304(R) (2010). [38] B. ten Hagen, S. van Teeffelen, and H. Lowen, Cond.
[21] S.vanTeeffelen andH.Lowen,Phys.Rev.E78, 020101 Mat. Phys. 12, 725 (2009).
(2008).
[22] B.tenHagen,R.Wittkowski,andH.Lowen,Phys.Rev.