Collective Dynamics of Deformable Self-Propelled Particles with Repulsive Interaction

1.5 MB·English

by  Yu Itino

Typesetwithjpsj3.cls<ver.1.1> Letter
CollectiveDynamics ofDeformableSelf-Propelled Particles
withRepulsiveInteraction
YuItino,TakahiroOhkumaandTakaoOhta
∗
DepartmentofPhysics,GraduateSchoolofScience,KyotoUniversity,Kyoto606-8502,Japan
Weinvestigatedynamicsofdeformableself-propelledparticleswitharepulsiveinteractionwhosemag-
1 nitude depends on the relative direction of elongation of a pair of particles. A collective motion of the
1 particlesappearsintwodimensions.Howeverthisorderedstatebecomesunstablewhentheparticleden-
0 sity exceeds a certain critical threshold and the dynamics becomes disorder. We show by a mean field
2 analysisthatthisnoveltransitioncharacteristictodeformabilityoccursduetoasaddle-nodebifurcation.

n
a KEYWORDS: self-propelledparticle,collectivedynamics,saddle-nodebifurcation.
J

8
Dynamics of self-propelled objects have been studied for with s themagnitude.Thelasttermineq.(2)standsforthe
1 i
thesemorethanonedecade.Variouscollectivemotionshave interactionbetweenparticlesandisgivenby
] been found for interacting particles1–3) where the hydrody-
t N
f namicformulationhasbeendeveloped.4,5)Seeref.6)forare- f(i) = K F Q , (5)
o ij ij
s viewofexperimentalstudies.Effectsofnoiseshavealsobeen Xj=1
. investigated in various model systems to find a rich variety
t where K is a positive constant, N the total number of parti-
a ofcollectivedynamics.7–11) However,almostalloftheprevi-
m ousstudiesassumethatobjectsarerigiddespitethefactthat cles and Fij = −∂U(rij)/∂rij with rij = ri − rj. The form
of the potential U(r ) will be specified later. The other fac-
- coupling between shape deformation and migration in self- ij
d propulsioniscommoninbiologicalmicroorganisms12,13).In tor Qij expresses the dependence of the relative angle of a
n pairofelongatedparticlessuchthatwhentheinteractingtwo
thepresentpaper,weareconcernedwiththedeformationef-
o
particlesareparallel,therepulsiveinteractionisweaker.Itis
c fectonthecollectivedynamicsofself-propelledparticles.
definedby
[ Quiterecently,wehaveintroducedthemodelsystemwhere

the objectsare deformable.14) We have shown that deforma- Q
1 Q =1+ tr(S(i) S(j))2, (6)
ij
v tionproducesalotofinterestingdynamicsofindividualpar- 2 −
7 ticlessuchasrotatingmotion,zig-zagmotion,chaoticmotion whereQisapositiveparameter.
8 andhelicalmotion.14–17)Theseresultsclearlyindicatethatde- We have introduced the force (5) to take into considera-
3
formation is relevant to the dynamics. In the present paper, tion, approximately, of the excluded volume effect of elon-
3
. westudydynamicsofassemblyofdeformableself-propelled gatedparticles. Since the force dependson the variableS(i),
1 αβ
particleswherethereisarepulsiveshortrangeinteractionbe- itispossiblegenerallythatthereisacountertermineq.(3).
0
1 tweenapairofparticlessuchthatthemagnitudedependson However,forsimplicity,suchacontributionisnotconsidered
1 therelativedirectionofelongation.Ourmainconcernissta- inthepresentstudy.
: bilityofcollectivedynamicsduetodeformationofindividual Equation (3) implies that if the velocity is zero, v(i) = 0,
v α
i particles. theparticledoesnotdeformandtakesacircularshapeintwo
X We consider interacting particles in two dimensions. The dimensions. One of the characteristic features of the set of
r time-evolutionof the i-th particle obeysthe followingset of equations (1), (2) and (3) is that, when the interaction term
a equationsforthepositionofthecenterofmassrα(i),theveloc- isabsent,theseequationsdonotcontaintheparticlesizeex-
ityv(αi) andthedeformationtensorSα(iβ) plicitlyasaparameter.Thereforetheinteractionrangeinthe
potentialUistheonlyfundamentallengthscaleofthesystem.
d
r(i) = v(i), (1) Weintroduceanorientationalorderparameterdefinedby
dt α α
ddtv(αi) = γv(αi)− v(i)2v(αi)−aSα(iβ)v(βi)+ fα(i), (2) Φ=(cid:12)(cid:12)N1 N e2iθj(cid:12)(cid:12), (7)
(cid:12)(cid:12) (cid:12)(cid:12) (cid:12)(cid:12) Xj=1 (cid:12)(cid:12)
ddtSα(iβ) = −κSα(iβ) (cid:12)+b(cid:12) v(αi)v(βi)− 21 v(i)2δαβ!, (3) whereθj isdefinedthroug(cid:12)(cid:12)(cid:12)(cid:12)htherelatio(cid:12)(cid:12)(cid:12)(cid:12)n n(i) = (cosθi,sinθi).
(cid:12)(cid:12) (cid:12)(cid:12) In order to quantify the degree of coherentmotion, we may
whereγ 0,κ 0,aandbareconstants.(cid:12)Th(cid:12)issetofequa- defineanalternativeorderparameter1)
≥ ≥
tions is an extension of a single particle case.14) The tensor
S isrepresentedintermsoftheunitnormal nparalleltothe 1 N v(j) 1 N
longaxisofanelongatedparticleas Φv = N (cid:12)(cid:12)(cid:12)(cid:12)Xj=1 v(j) (cid:12)(cid:12)(cid:12)(cid:12)=(cid:12)(cid:12)(cid:12)(cid:12)N Xj=1eiφj(cid:12)(cid:12)(cid:12)(cid:12), (8)
Sα(iβ) = si n(αi)n(βi)− 12δαβ!, (4) wherewehaveputv(i)(cid:12)(cid:12)(cid:12)(cid:12)= v((cid:12)(cid:12)(cid:12)i) (c(cid:12)(cid:12)(cid:12)o(cid:12)(cid:12)(cid:12)(cid:12)sφ(cid:12)(cid:12)(cid:12)(cid:12)i,sinφi).In(cid:12)(cid:12)(cid:12)(cid:12) thefollowing,
weshalluseΦbutweha(cid:12)vev(cid:12)erifiednumericallythatthedif-
(cid:12) (cid:12)
(cid:12) (cid:12)
J.Phys.Soc.Jpn. Letter YuItinoTakahiroOhkumaandTakaoOhta
Fig. 3. (Coloronline)Hysteresisloopfor(a)N =512,(b)N =2048and
(c)N=8192
.
the dilute regime ρ = 0.0032 and the intermediate regime
Fig. 1. (Color online) Snapshots ofthe time-evolution ofparticles at(a)
ρ = 0.0852obtainedfor N = 512. Initially the velocityand
t = 256and(b)t = 512forρ = 0.075and N = 8192starting witha
thepositionoftheindividualparticlesarerandom.Whenthe
randominitialconfiguration.Thesefiguresdisplayonlythe1/4 1/4area
ofthesystem.Thecolorsdistinguishthedirectionofelongatio×nandthe density is low, the order parameter starts to grow at around
arrowsindicatethedirectionofmigration. t = 3000 and almost saturates at about t = 7000 as in Fig.
2(a). When the density is high, the transition to the ordered
stateismoreabrupt.Figure2(b)indicatesthatthereisanin-
cubationtimeandtheorderparameterstartstogrowataround
t=1500andcompleteswithinanarrowinterval.
Itseemsthatthecollectivemotionappearsasymptotically
in time even for a dilute limit since once an orderedmotion
appearsaftercollisions,thereisnomechanismtodestroyit.If
oneaddsnoisetermsinthetime-evolutionequations,asharp
transition from the disordered state in the dilute regime to
the ordered state is expected to occur. However, we do not
Fig. 2. (Color online) Growth of order for (a) ρ = 0.0032 and (b) ρ =
carryoutsuchnumericalsimulationsbecausetherehavebeen
0.0852,andfor N = 512.Thecurvegiven byeq.(10)isshownbythe
brokenlineinFig.(a). manystudiesforrigidparticleswithnoises7–11) and,further-
more,thedeformabilityisexpectedtobeirrelevantinthedi-
luteregime.
ferencebetweenΦ andΦv is negligiblysmall in theplotsin The time-dependence of the order parameter in Fig. 2(a)
Figs.2and3below.
canbefittedbythefollowingfunction
We shallsolveeqs.(1),(2)and(3)numerically.Hereafter
weshowtheresultsforthepotential Φ(t)=γ 1 e−α(t−t0)β , (10)
−
h i
r2 whereα=5.02 10 11,β=3,γ=0.955andt =2.25 103.
U(r )=exp ij . (9) × − 0 ×
ij −2σ2 Ift0isputtobezero,nosatisfactoryfitispossible.However,
Theparametersarechosenasa = 1andb = 0.5sothateach tthheisatbyrpueptocfhfaunngcetioshno.wTnhiisniFmigp.li2e(sbt)hcaatntnhoettbraenrseiptiroenseknitneedtibcys
particletendstoelongatetothedirectionparalleltotheveloc- is quitedifferentfromthe ordinarynucleationandgrowthin
ity.14) We fixtheparametersγ = κ = 1in whichtheparticle firstorderphasetransitionsinthermalequilibrium.Infact,we
undergoes straight motion.14) Unless stated explicitly, other have seen that once an ordered region appear in the present
parametersarechosenasσ = 1, K = 5andQ = 50.Theset self-propelled system, it spreads over the system almost in-
of equationsis solved by the Euler explicit methodwith the stantaneouslywhenthedensityismoderatelyhigh.
time incrementδt = 1/64.Particlesare confinedin a square Theorderedstatebecomesunstableforlargevaluesofden-
areawithsizeLandtheperiodicboundaryconditionsareim- sity.Wehaveinvestigatedthedisorderingprocessbyincreas-
posed.Inmostofthesimulations,wechooseN =512,butthe ing the density slowly. The compression and expansion are
cases N =1024, 2048 and 8192 are also examined to check madeataconstantrate
the numberdependence.The density ρ = N/L2 is varied by
dρ/dt= λ, (11)
changingthesystemsizeL.
±
Appearanceofthecollectivemotionofself-propelledpar- where
ticlesisshowninFig.1startingfromadisorderedinitialcon-
λ=(ρ ρ )2 n, (12)
dition.Theparticlesconstitutelocallyanapproximatehexag- 1 0 −
−
onal lattice such that one of the fundamental reciprocal lat- with typical values ρ = 0.07 and ρ = 0.14. The expo-
0 1
tice vectors is parallel to the propagatingdirection. Another nent n is varied. Figures 3 display the order parameter Φ as
orientationofthelatticesuchthatoneoftheprimitivelattice a function of ρ for n = 14. It is evident that the ordered
vectorsisparalleltothepropagatingvelocityhasnotbeenob- stateisspontaneouslybrokenwhenthedensityexceedsacer-
served.Wenotethatthepropagatinglatticeisnotcompletely tainthresholdvalueandthereisahysteresisinthetransition.
stationary.Eachparticleexhibitsasmalloscillation(orfluctu- Figure 3 shows the N-dependence of the hysteresis behav-
ation)initscenterofmassperpendicularlytothepropagating ior. The threshold ρ for the transition from the disordered
do
direction. state to the ordered state by decreasing the density is less
Figure 2 shows the growth of the order parameter Φ for sensitive to the system size. We note that the threshold for
J.Phys.Soc.Jpn. Letter YuItinoTakahiroOhkumaandTakaoOhta
Fig. 5. (Color online) Diffusive property of the disordered state for ρ =
0.14andN=512.(a)Self-intermediatescatteringfunctionforwavenum-
Fig. 4. (Coloronline)Theoder-disordertransitiondensityρodasafunction berskσ=1/16,kσ=1/8,kσ=1/4,kσ=1/2andkσ=1fromtopto
ofQforN=512.Thedottedlineisaguidetotheeye.Thevolumefraction bottom.Thedottedcurvesshowthefunctiongivenbyeq.(15).(b)Mean
ofthetotalparticlesdefinedbyφod=πσ2ρodisalsoindicatedasthescale squaredisplacementaveragedover16independentruns.
oftheverticalaxis.
creasesfromβ 1.09toβ 1.3byincreasingkσfrom1/16
k k
the ordered-to-disorderedstate by increasing the density oc- to1.Thisindica≈testhattheus≈ualdiffusionwiththewavenum-
curs at ρ 0.132 for N = 512 = 29, ρ 0.125 for
od od berdependenceβ ν 2occursforsmallvaluesofkwhereas
≈ ≈ k
N = 2048 = 211 andρ 0.115for N = 8192 = 213.This ≈
od a ballistic motion is dominant in the short length scale and
≈
indicatesthatthehysteresisregionρod ρdohasatendencyto clearly excludesthe stretched exponentialwith β 1.19) It
− k
decreaselinearlywithlnN.However,moreaccurateanalysis ≤
isalsoverified(notshown)thatthesimpleexponentialdecay
isnecessarytoobtainanydefiniteconclusion,whichisleftfor
independentofthevaluesofkislesssatisfactorytofitthenu-
afuturestudybecauseoftimeconsumingofnumericalsimu-
merical data. The mean square displacement is displayed in
lationsforlargersystem.Wehavealsocarriedoutnumerical
Fig. 5(b).Thishasbeenobtainedby the averageof16 inde-
simulationsforN =512forslowerchangeofthedensitywith
pendentruns.Itisapproximatelyproportionaltotimeasymp-
n=16andn=18.Thethresholddensityρ isapproximately
od toticallywiththemagnitudeofD 0.2.Itisevidentthatthere
givenbyρod = 0.132forn = 14,ρod = 0.131forn = 16and is no regime where C(t) has a pl≈ateau due to a cage effect.
ρod = 0.128forn = 18.Ontheotherhand,thethresholdρdo This also implies that the disordered state is different from
isgivenbyρdo = 0.0840forn= 14,ρdo = 0.0891forn= 16 theusualglassstate18) .
and ρ = 0.0902 for n = 18. Therefore, as expected, the
do We have shown that the coherent ordered state of the de-
hysteresisregiontendstoshrinkforslowerdensitychange. formableself-propelledparticlesbecomesunstable for suffi-
Thetransitionbetweentheorderedanddisorderedstatesin
cientlylargedensityandadisorderedstateappears.Wehave
the high density regime has not been observed in any rigid
verified that the same transition occurs not only for the soft
self-propelledsystemsandischaracteristictothedeformable
potentialgivenbyeq.(9)butalsofor
particles.Theorderedstateexistsirrespectiveofthevalueof
theanisotropicparameterQineq.(6).However,thestability U =exp( rij/σ). (16)
−
limit of the ordered state depends on it. The critical density
For example, the order-disorder transition occurs at around
decreases as the value of Q is increased as shown in Fig. 4. ρ = 0.02fordecreasingthedensityandρ = 0.06forincreas-
When the anisotropy is absent, Q = 0, the transition to the ingdensityfor N = 512keepingalltheparametersthesame
disorderedstatedoesnotoccur.
valuesasusedabove.
We have exploredthe propertiesin the disorderedstate at
Herewedescribeameanfieldanalysistoclarifytheorder-
the high density regime by evaluating numerically the self-
disordertransitioninthehighdensityregime.Theforceterm
intermediatescatteringfunctiondefineby
ineq.(2)canbewrittenas
N
1 N
Fs(k,t)= N *Xj=1e−ik·(rj(t)−rj(0))+. (13) v(i)· f(i) = Kv(i)Xj=1vˆ(i)·FijQij = Kv(i)NDvˆ(i)·FijQijE, (17)
Themeansquaredisplacementofthe individualparticlecan
wheretheunitvectorvˆ(i) = v(i)/ v(i) isregardedasastochas-
beobtainedfromF (k,t)as
s ticvariable.Theaverageisdefin(cid:12)eds(cid:12)uchthat
(cid:12) (cid:12)
∂2F (k,t) (cid:12) (cid:12)
C(t)= s = r(t) r(0)2 =4Dt, (14) 1 N
∂k2 (cid:12)(cid:12)k=0 D| − | E Xj = N Xj. (18)
whereDisthediffusion(cid:12)(cid:12) constantintwodimensions. D E Xj=1
Theintermediatescatteringfunctionisshownasafunction First, we make a decoupling approximationsuch that Q in
ij
of time in Fig. 5(a) for several values of k = k. It is clear
| | eq.(17)isreplacedbytheaverage Qij ,
that there is no multiple relaxation which is typical in glass
D E
dynamics.18) Actually, as shown in Fig. 5(a), the scattering v(i) f(i) = Kv(i)N vˆ(i) F Q . (19)
ij ij
functioncanbefittedby · D · ED E
By a dimensional analysis, we note that N vˆ(i) F might
Fs(k,t)=exp[−(t/τk)βk], (15) haveafactorσρandhaveverifiednumericallDytha·tthijiEsquan-
where τ 23.46 k ν with ν = 1.85.The exponentβ in-
k − k
≈ ×
J.Phys.Soc.Jpn. Letter YuItinoTakahiroOhkumaandTakaoOhta
transitionsduetoasaddle-nodebifurcationhavebeenstudied
theoreticallyinasingle-variableLangevinsystem.20)
Finally we discuss some apparently related phenomena.
Colloidalparticles(althoughnotself-propelled)exhibitglass
behavior for high density. In this glass dynamics, motion of
each particleis frozenor trappedin a cageanddoesnotun-
dergoausualBrownianmotion.Thisiscontrasttothepresent
result in Fig. 5 where particles have a finite diffusion con-
stantthoughtheoriginofthemotionisnotthermalbutcomes
from the self-propelledenergy.A jammingtransition occurs
Fig. 6. CorrelationχasafunctionofρforN=512.Themeanvaluesare ingranularmaterialsundershearflow.Thediffusionconstant
indicatedbytheblackdotswhereasthestandarddeviationsareshownby is finitein the jammedstate butitis definedperpendicularly
theerrorbars.Thisplotisobtainedbydecreasingthedensityslowly.The to the flow direction21) and inherently anisotropic since the
transitionoccursatρ=0.082.Theinsetshowsthatχisalsonegativefor
systemissubjectedtoshearflow.Thethirdrelatedsystemisa
theorderedstate.
trafficflowwhichalsoexhibitsajammingstatewhenparticles
are dense.However,the flowis generallydirectedin onedi-
tityisindeedproportionaltothedensityρsothatweput
mension22,23) or anisotropicallyone-dimensionalevenin the
N vˆ(i) F =ρσχ. (20) two-dimensionalspace24) andthereforethejammedstatehas
ij
·
D E nodirectanaloguewiththepresentisotropicdisorderedstate.
Theunknownquantityχisto be evaluatednumerically.An-
other decoupling approximation (S )2 = S 2 is also Acknowledgements
αβ αβ
employed.Secondly,weeliminateDSαβ byEputtDingdESαβ/dt = ThisworkwassupportedbytheGrant-in-Aidfortheprior-
0ineq.(3).Byusingtheseapproximations,eq.(2)becomes ity area“SoftMatterPhysics”andtheGlobalCOEProgram
1 d “TheNextGenerationofPhysics,SpunfromUniversalityand
w=γw gw2+χˆ√w[1+Qˆ(w v2 )2], (21) Emergence” both from the Ministry of Education, Culture,
2dt − −
D E Sports,ScienceandTechnology(MEXT)ofJapan.
where w = v2 0 with the superscript (i) omitted and
g=1+(ab)/(2κ),≥andQˆ = Qb2/(4κ2)andχˆ = Kρχσ.Putting
v2 = v 2 + (δv)2 with δv the fluctuating part, we iden-
1) T.Vicsek,A.Cziro´k,E.Ben-Jacob,I.Cohen,andO.Shochet:Phys.
h i
DtifyE v 2 with tDhe staEtionary stable solution of eq. (21). We Rev.Lett.75(1995)1226.
h i 2) J.TonerandY.Tu:Phys.Rev.Lett.75(1995)4326.
haveverifiednumericallythatχisnegativeasshowninFig.
3) N.Shimoyama, K.Sugawara, T.Mizuguchi, Y.Hayakawa, and M.
6. Therefore,we note that w = 0 is always a stable station-
Sano:Phys.Rev.Lett.76(1996)3870.
ary solution and that there is a pair of finite solutions, one 4) A.BaskaranandM.C.Marchetti:Phys.Rev.E77(2008)011920.
of which is unstable and the other larger solution is stable. 5) S.Ramaswamy:Annu.Rev.Condens.MatterPhys.1(2010)323.
For sufficiently large values of χ, this pair of solutionsdis- 6) S.J.EbbensandJ.R.Howse:SoftMatter6(2010)726.
| | 7) G.Gre´goireandH.Chate´:Phys.Rev.Lett.92(2004)025702.
appear. Thereforethere is a saddle-nodebifurcationat some
8) U.Erdmann,W.Ebeling,andA.Mikhailov:Phys.Rev.E71(2005)
valueofρ ρ beyondwhichonlythesolutionw=0exists
≡ od 051904.
andhencethisstatecorrespondstothedisorderedstate.Inthe 9) F.PeruaniandL.Morelli:Phys.Rev.Lett.99(2007)010602.
disorderedbranch,weshouldput v 2 =0ineq.(21)whereas 10) S.vanTeeffelen,U.Zimmermann,andH.Lo¨wen:SoftMatter5(2009)
h i
the stable finite stationary solution should be equated to v 4510.
h i 11) F.Ginelli,F.Peruani,M.Ba¨r,andH.Chate´:Phys.Rev.Lett.104(2010)
intheorderedbranch.Thefactthatχ < 0isessentialforthe
184502.
saddle-nodebifurcation.However,aquantitativecomparison
12) L.Li,S.F.Nørrelykke,E.C.Cox,andM.Hatakeyama:PLoSONE3
is impossible with the present crude approximation.For ex- (2008)e2093.
ample,we haveto requireχ = 0.57(χ = 0.73)to realize 13) Y.T.Maeda,J.Inose,M.Y.Matsuo,S.Iwaya,M.Sano,andM.Isalan:
− −
ρ = 0.12 (ρ = 0.082), which are substantially different PLoSONE3(2008)e3734.
od do
fromthevaluesinFig.6.Notealsothatthecoefficientofthe 14) T.OhtaandT.Ohkuma:Phys.Rev.Lett.102(2009)154101.
15) T.Ohta,T.Ohkuma,andK.Shitara:Phys.Rev.E80(2009)056203.
w5/2 termwhichisdominantforlargevaluesofwineq.(21)
16) T. Hiraiwa, K. Shitara, and T. Ohta: Soft Matter (2011) DOI:
has a factor ρQ. This implies that, if χ and (δv)2 are not 10.1039/c0sm00856g.
stronglydependenton Q,thethresholdρ deDcreaseEsas Qis 17) T.Hiraiwa,M.Y.Matsuo,T.Ohkuma,T.Ohta,andM.Sano:Europhys.
od
increased,whichisconsistentqualitativelywithFig.4. Lett.91(2010)20001.
18) K.BinderandW.Kob:Glassymaterialsanddisorderedsolids(World
Wementionthattheabovemeanfieldtheorydoesnotpre-
scientific,Singapore,2005).
dict the intrinsic hysteresis behavior of the transition in the
19) O. Dauchot, G. Marty, and G. Biroli: Phys. Rev. Lett. 95 (2005)
slow limit of the density change and the N limit. 265701.
→ ∞
We need to take account of the stochasticity which is self- 20) M.IwataandS.Sasa:Phys.Rev.E82(2010)011127.
generated from the nonlinear interaction among the self- 21) P.Olsson:Phys.Rev.E81(2010)040301.
22) O.J.O’Loan,M.R.Evans,andM.E.Cates:Europhys.Lett.42(1998)
propelledparticles.
137.
Weemphasizethatthiskindoforder-disordertransitionhas
23) A.Nakayama, M.Fukui, M.Kikuchi, K.Hasebe, K.Nishinari, Y.
notbeenobserved,toourknowledge,inanyrigidrod-shaped Sugiyama,S.Tadaki,andS.Yukawa:NewJ.Phys.11(2009)083025.
particlesandhenceitisacharacteristicfeatureofdeformable 24) O.Biham, A.Middleton, and D.Levine: Phys.Rev. A 46 (1992)
softparticles.Itismentionedherethatnonequilibriumphase R6124.