Table Of Content

Elasticity-driven collective motion in active solids and active crystals Eliseo Ferrante Laboratory of Socioecology and Social Evolution, Katholieke Universiteit Leuven, Leuven, Belgium Ali Emre Turgut University of Turkish Aeronautical Association, Ankara, Turkey Marco Dorigo IRIDIA, Université Libre de Bruxelles, Brussels, Belgium Cristián Huepe 614 N. Paulina St., Chicago IL 60622-6062, USA 3 (Dated: January 15, 2013) 1 0 Weintroduceasimplemodelofself-propelledagentsconnectedbylinearsprings,withnoexplicit 2 alignmentrules. Belowacriticalnoiselevel,theagentsself-organizeintoacollectivelytranslatingor n rotatinggroup. Wederiveanalyticalstabilityconditionsforthetranslatingstateinanelasticsheet a approximation. We propose an elasticity-based mechanism that drives convergence to collective J motion by cascading self-propulsion energy towards lower-energy modes. Given its simplicity and 1 ubiquity, such mechanism could play a relevant role in various biological and robotic swarms. 1 ] Animal groups that move together, such as bacterial it was initially surprising that such systems could self- t f colonies, insect swarms, bird flocks, or fish schools [1–6], organize without them. While it can be argued that all o are all examples of biological systems displaying Collec- thesemodelsincludeatleastanimplicitalignmentinter- s . tive Motion (CM). In recent years, the dynamics of such action, it remains unclear if they are all driven to CM t a systems(referredtoheregenericallyasswarms)hasbeen by the same underlying mechanisms and to what extent m the subject of intense research [6–9]. A number of theo- agents must exchange orientation information, either ex- - reticalmodelshavebeenintroducedtostudyswarmsand plicitly or implicitly, to achieve CM. d to develop control rules that achieve similar coordinated In this Letter, we introduce a CM mechanism that is n o collectivedynamicsingroupsofautonomousrobots[8,9]. based on a very different paradigm: the emergence and c Despite this proliferation of CM algorithms, there is still growth of regions of coherent motion due to standard [ nocomprehensiveunderstandingoftheunderlyingmech- elasticity processes. We explore this mechanism by in- 1 anisms that can lead a group of self-propelled agents to troducing a simple two-dimensional Active Elastic Sheet v self-organize and move in a common direction. (AES)modelwithspring-likeinteractionsbetweenneigh- 0 boring agents and no explicit alignment, which describes ThecurrentCMparadigmhasbeenstronglyinfluenced 2 what we refer to as an active solid or an active crystal. 6 by the seminal work of Vicsek et al. [10], which intro- WedefinetheAESmodelasasystemofN agentsona 2 duced a minimal model for flocking, the Vicsek model, . that has become a referent in the field [7–9]. This model two-dimensional plane, where the position (cid:126)xi and orien- 1 tationθ ofeachagentifollowtheoverdampedequations 0 describes a group of point particles advancing at a fixed i of motion 3 common speed, only coupled through alignment interac- 1 tions that steer them towards the mean heading direc- (cid:126)x˙ =v nˆ +α (cid:104)(cid:16)F(cid:126) +D ξˆ(cid:17)·nˆ (cid:105) nˆ , (1) : i 0 i i r r i i v tion of all particles within a given radius [10–12]. In this (cid:104)(cid:16) (cid:17) (cid:105) i framework, a swarm can be viewed as a group of self- θ˙ =β F(cid:126) +D ξˆ ·nˆ⊥ +D ξ . (2) X i i r r i θ θ propelled spins with aligning interactions, described by r a a variation of the XY-model [13] where spins advance in Here, v0 is the forward biasing speed that induces their pointing direction rather than remaining affixed to self-propulsion (injecting energy at the individual par- a lattice. In the continuous limit, this system becomes ticle level), nˆ and nˆ⊥ are two unit vectors point- i i a fluid of self-propelled spins that follows the hydrody- ing parallel and perpendicular to the heading direction namic theory developed in [14–17]. More recently, other of agent i, and parameters α and β are the inverse modelsthatdonotrelyonexplicitalignmentinteractions translationalandrotationaldampingcoefficients,respec- havebeenintroduced. In[18], forexample, CMisdriven tively. The total force over agent i is given by F(cid:126) = i (cid:80) by escape-pursuit interactions only; in [19], by inelastic (−k/l )(|(cid:126)x −(cid:126)x |−l ), a sum of linear spring- j∈Si ij i j ij collisions between isotropic agents; and in [20] and [21], like forces with equilibrium distances l and spring con- ij byshort-rangeradialforcescoupledtoeachagent’sturn- stants k/l . Each set S contains all agents interacting ij i ing dynamics. Given that Vicsek-like algorithms rely on with agent i and remains fixed throughout the integra- explicit alignment rules to achieve CM [5, 10, 22, 23], tion. This system is thus akin to a spring-mass model 2 of elastic sheet [24] where masses are replaced by self- A B C propelledagentsthatturnaccordingtoF(cid:126) ·nˆ⊥ andmove 1 i i forward or backwards following F(cid:126) · nˆ and their self- i i propulsion. We include actuation noise (fluctuations of the individual motion) by adding D ξ to the heading θ θ angle, where D is the noise strength coefficient and ξ θ θ a random variable with standard, zero-centered normal probability distribution of variance 1. We include sens- 2 ing noise(errorsinthemeasuredforces)byaddingDrξˆr toF(cid:126) ,withD thenoisestrengthcoefficientandξˆ aran- i r r domly oriented unit vector. The degree of alignment in the system is monitored through the usual polarization order parameter (cid:13) (cid:13) 1 (cid:13)(cid:88)N (cid:13) 3 ψ = N (cid:13)(cid:13) nî(cid:13)(cid:13). (3) (cid:13) (cid:13) i=1 If all agents are perfectly aligned, we have ψ = 1. If agents are instead randomly oriented or rotating about the group’s barycenter, we have ψ =0. TheAESmodelwasdesignedtostudyCMundercon- ditions that are in many ways opposite to those in the FIG. 1. Active elastic sheet simulations of Eqs. (1) and (2). Vicsek model. While the only information that agents A:Hexagonalactivecrystalatt=0(A1),240(A2),and1700 exchange there is their relative heading angle, here they (A3). B: Rod-like active crystal at t=0 (B1), 400 (B2), and only sense their relative positions. While changing in- 1700(B3). C:ActivesolidatsametimesascolumnB;darker teracting neighbors over time has been shown to be nec- agents symbolize higher local alignment. essary there for achieving long-range order [11, 14, 25], here virtual springs connect the same agents throughout the dynamics. Furthermore, while both models de- lower-energy, translating solutions. However, we show scribe overdamped systems, these are implemented dif- one here to illustrate its dynamics, which cannot be at- ferently. In the Vicsek model, particles switch instanta- tained by the Vicsek model. neouslytothedesiredheadingangle;intheAEScase,we Column B displays an active elastic rod, comprised of integrate instead the overdamped angular equation (2), N = 118 agents arranged into three rows, for the same which turns out to be necessary here for achieving CM. noiseasincolumnA.Itisgenerated(B1)byplacingran- WeintegrateEqs.(1)and(2)numericallyusingastan- domly oriented agents with nearest-neighbor distances dardEulermethod. Allsimulationsbelowarecarriedout dB = 0.32 (within rows) and d∗B = 0.58 (between rows), with α=0.01, β =0.12, v =0.002, and dt=0.1. linkingallagentsseparatedbyd<1withspringsofnat- 0 Figure 1 presents three runs of the AES model. Col- urallengthdandspringconstantk/l=5/d. Hereagain, umnAdisplaysthedynamicsofanhexagonalactivecrys- largerandlargerregionsofcoherentdeformationemerge tal composed of N = 91 agents. At t = 0 (panel A1), (B2)untilCMisattainedandtherodstartsmoving(B3). agents are placed with random orientations on a per- Since the first bending mode has the largest final defor- fect hexagonal lattice, separated by d = 0.65. Nearest mation, a collective heading direction perpendicular to A neighbors are connected by springs with natural length the rod’s axis is favored. l = d and spring constant k/l = 5/0.65. Only sensing Column C shows N = 891 agents forming an active A noise is considered here (D = 0.5, D = 0), but results solid(giventheirregularagentpositions)withtwoholes. r θ remain qualitatively unchanged for other types of noise. Toconstructit,wedistributeagentsatrandom,homoge- Astimeadvances,growingregionsofcoherentmotionde- neously within the structure, connecting all agents sep- velop,eventuallydeformingthewholestructure(A2)un- arated by d < 1 with springs, following the same pro- tilthegroupstartstranslatingorrotatingcollectively. In cedure used in column B. Noise is set here to zero, but the case shown, the system converges to a rotating state equivalent dynamics are observed for small enough D r wheretheaxisofrotationandbarycenterdonotcoincide andD . Tohighlight orderedregions, eachagent’sdark- θ (A3), thusrotatingwhiletranslating. Notethatrotating ness is displayed proportional to the local order, defined states will always have higher elastic energy, since inner similar to ψ but summing only over the focal agent and and outer shells cannot move at the same v speed and otherslinkedtoit,insteadofthewholesystem. Initially, 0 must be sped up or slowed down by elastic forces. They most of the structure appears in light gray, since agents are less frequent and metastable, eventually relaxing to are randomly oriented (C1). As time advances, growing 3 1.0 11..00 AninterestingaspectoftheAESmodelisthatwecan 0.8 00..88 use a continuous elastic sheet approximation to perform analytical calculations. We follow this approach to carry 0.6 00..66 out a standard linear stability analysis [24] of the trans- 0.4 00..44 latingCMstateinthezeronoisecase. Webeginbywrit- 0.2 00..22 ing the elastic forces F(cid:126) =(F ,F ) that result from small x y 0 1.000 displacements (cid:126)u = (ux,uy) of points on the membrane with respect to their equilibrium positions 0.5 00..86 0 N=91, Dr=0.822 N=5=47=,= Djkrj=k0d.s8 06 000...645 N=91, Dθ=0.097 N=547, Dθ=0.095 Fx =(λ+2µ)∂∂2xu2x +µ∂∂2yu2x +(λ+µ)∂∂x2u∂yy, (4) -0.5 0 0.5 1 0 0.5 1 000...0243 0 0.5 1 0 0.5 1 Fy =(λ+2µ)∂∂2yu2y +µ∂∂2xu2y +(λ+µ)∂∂x2u∂xy, (5) where the elastic constants are given by the Lamé parameter λ and shear modulus µ [24]. We then linearize FIG. 2. (color online). Order parameter and Binder cumu- Eqs.(1)and(2)aroundanequilibriumsolutionwithun- lant vs. positional sensory noise Dr and angular actuation deformed membrane and all agents moving at speed v0 noise D for hexagonal active crystals with N =91 (same as θ in the xˆ direction, obtaining u˙ = αF , u˙ = v φ, and onFig.1,panelA1)andN =547agents. Toppanelsdisplay x x y 0 φ˙ = βF , where φ denotes perturbations to the θ = 0 the mean and local maxima of the distribution of ψ values y obtained in simulations. Insets show these distributions in equilibrium heading angle. Casting these expressions in the transition region. For large enough systems, both cases Fourier space with wavevector components (kx,ky), we display a first order transition with a bistable region. can write the perturbation dynamics in matrix form and compute its eigenvalues Λ to determine stability. These are found to satisfy the characteristic polynomial equa- regions of coherent motion emerge (C2), until the whole tion Λ3+C2Λ2+C1Λ+C0 =0, with structure starts moving when agents become sufficiently aligned (C3). C0 =αβµv0(λ+2µ)(cid:2)kx2+ky2(cid:3)2, (6) The AES model displays a discontinuous order- C1 =βv0(cid:2)µkx2+(λ+2µ)ky2(cid:3), (7) disorder transition similar to that in the Vicsek model. C =α(cid:2)(λ+2µ)k2+µk2(cid:3). (8) 2 x y Figure 2 examines this transition as a function of noise forthehexagonalactivecrystaloncolumnAofFig.1and Using Routh’s stability criterion (here given by C C > 1 2 foralarger(N =547)hexagonalconfigurationwithiden- C [27]), we find that the system is stable if 0 tical parameters. We performed 30 to 80 runs per noise αβv (λ+µ)2k2k2 >0,whichisalwaysverified. Wecon- 0 x y value, storing 2000 values of ψ per run (every 500 time- clude that translating CM solutions are always linearly stepsaftertheinitial106). Toppanelsshowthemeanand stable. This is not the case, however, for most variations localmaximaoftheψdistributionsandbottomones,the of the AES model. For example, if we consider a con- corresponding Binder cumulants G=1− 1(cid:10)ψ4(cid:11)/(cid:10)ψ2(cid:11)2 stantspeedalgorithmbysettingα=0,thecharacteristic [26]. As we increase either sensing noise D3 (left) or ac- polynomialbecomesΛ3+βv (cid:2)µk2+(λ+2µ)k2(cid:3)Λ=0, r 0 x y tuation noise D (right), ψ jumps from an ordered state which only has null or imaginary solutions. Linear per- θ where agents self-organize to a disordered state where turbations will therefore not dampen out, but produce they continue to point in random directions. The dis- instead permanent oscillations. Numerical simulations playedBindercumulantsandbimodaldistributionsshow confirmthat,evenforzeronoiseandstartingfromaper- that there is a region of bistability around the critical fectlyalignedinitialcondition, thegrouplosesorderand point in both cases. In the sensing noise case, G < 0 agents start rotating in place instead of aligning. close to the transition for N = 547, providing evidence We now characterize the nonlinear energy cascading forbistabilityinlargeenoughsystems. Fortheactuation mechanism that drives the AES model to self-organize. noise case, N =547 does not appear to be large enough Figure 3 presents the energy dynamics of an hexagonal toreachG<0,butthedipinthetransitionregiondrops activecrystalsimulationwithN =91andzeronoisethat further and further below G = 1/3 (the expected value converges to a translating solution. Top panels display for the disordered phase) as the system size is increased, the total kinetic and potential energies as a function of suggestingthatGwillreachnegativevaluesinlargersys- time. Panel C shows the spectral decomposition of the tems [12, 23, 26]. These numerical tests (and others we latter into its elastic normal modes, listed in order of performed with different noise types and agent configu- growingenergyandwithoutaccountingfordegeneracies. rations) show that the AES transition is first order and It is produced by first computing all 182 elastic normal has a bistable region. modes numerically (without considering agent orienta- 4 ³ 00..220 them. Self-propulsion thus feeds energy to lower and 0 1 E 0.15 lower modes, eventually reaching translational or rota- c00..110 A tional modes and achieving CM. Note that, despite this ite 0.05 mechanism, similar models may not converge to CM. n iK0.000 This will occur, for example, if agents inject too much 0.5 0 500 1000 1500 2000 2500 E 00..44 energy into high-energy modes while turning (as in the c 0.3 α=0 case described above) or if they overshoot the an- itsa 00..22 B gles that dampen high-energy modes by instantaneously lE 0.1 switching heading (as in [22, 23]) instead of integrating 0.00 180 0 500 1000 1500 2000 2500 Eq. (2). We have identified in this Letter an alternative, r elasticity-based mechanism that, in contrast to the Vic- eb m120 sek case, requires no exchange of heading information u edN C tinogancheiigehvbeoCrsMo.vIetratilmsoerteoquovireersconmoeswthitechMinergmoifni-nWteargancetr- o 60 theorem and achieve long-range order at non-zero noise M levels[11,14,25,28]. Instead,despiteincludingnofluid- likemixing,wefoundthattheAESmodeldisplaysnoloss 1 oflong-rangeorderforlargerstructures,althoughtheen- 0 500 1000 1500 2000 2500 Time ergy cascading mechanism takes longer to converge. We found only one other model, introduced in [20] to study FIG. 3. (color online). Kinetic energy (A), elastic energy the collective migration of tissue cells, that can display (B), and spectral decomposition of the elastic energy (C) as CM under similar conditions. In a version of this algo- a function of time for an hexagonal N = 91 active crystal rithm (designed to study active jamming) where agents simulation (with zero noise and same initial condition as in only have repulsive interactions and are confined to a panelA1ofFig.1)thatconvergestothealigned,translating state. BrighterpointsonCindicatehigherenergies. Afteran circular box, a similar elasticity-based mechanism was initial transient, A and B converge to their stationary values shown to be responsible for the dynamics of the jammed for collective translational motion. All modes display energy phase [21]. While that model describes a different situa- levelsthatoscillateastheydecay,withhighermodesdecaying tion,whereinteractionforcescandisplaceagentsperpen- faster. Elasticenergyflowstolowermodes,producingcoher- dicular to their heading, its detailed comparison to the ent motion that eventually reaches the lowest (translational AES model should yield a more complete understanding or rotational) modes. of the energy cascading mechanism. Given that some kind of attraction-repulsion interaction must be present for any group of moving agents to tions or self-propulsion) and then expanding the dynam- remain cohesive and avoid overlapping, we expect the ics into this basis. The initial condition is set as in panel elasticity-basedmechanismintroducedheretoplayarel- A1 of Fig. 1, with zero potential energy and kinetic en- evant role in the CM dynamics of a variety of systems. ergy E =Nv2/2=1.82×10−4 (setting the agent mass These could include microscopic biological or robotic k 0 to 1). As the membrane deforms during the initial tran- agents(orevencollectivelymigratingtissuecells[20,29]) sient, potential energy grows and becomes broadly dis- that can exert attraction-repulsion physical forces while tributed over all modes (as expected for disordered sys- being too simple to exhibit explicit aligning interac- tems), whilekineticenergydrops. Astimeadvances, the tions. We expect most animal groups to typically com- system rearranges itself into configurations with lower bineelasticity-basedandalignment-basedmechanismsin elasticenergyandhigherkineticenergy,eventuallyreach- theirCMdynamics,thuseffectivelybehavingasanactive ing again (now in the ordered, translating state) values viscoelastic material composed of aligning spins. An in- close to zero and E , respectively. The energy of each terestingopenquestionistheextenttowhicheachmech- k mode oscillates while decaying, as in an underdamped anism is responsible for the CM of specific real-world oscillator, with higher modes decaying faster than lower swarms,whichcandependonthetime-scaleanddynam- ones. This results from a combination of standard elas- ical state considered. Note that each mechanism results ticity, self-propulsion, and the coupling between elastic frominteractionsthatproducedifferentdynamicalsigna- forces and turning rate imposed by the AES model. In- tures, such as the properties of propagating waves or the deed, in standard damped elastic systems, higher energy response to perturbations. Now that there is a growing modes also decay faster than lower ones. Here, how- number of experiments that allow the precise tracking of ever,eachagentiscontinuouslyinjectingenergythrough different types of swarms [30, 31], these signatures could itsself-propulsionterm,somotioncannotbefullydamp- help determine their individual interaction rules based ened. Instead,modesdecaybysteeringagentsawayfrom only on observed properties of their collective dynamics. 5 TheworkofCHwassupportedbytheNationalScience Newman, The Theory of Critical Phenomena: An In- FoundationunderGrantNo. PHY-0848755. Theworkof troductiontoRenormalizationGroup (OxfordUniversity EF and AT was partially supported by the Vlaanderen Press, 1992). [14] J. Toner and Y. Tu, Phys. Rev. Lett. 75, 4326 (1995). Research Foundation Flanders (Flemish Community of [15] J. Toner and Y. Tu, Phys. Rev. E 58, 4828 (1998). Belgium) through the H2Swarm project. The work of [16] E. Bertin, M. Droz, and G. Grégoire, Phys. Rev. E 74, EF, AT and MD was partially supported by the Euro- 022101 (2006). pean Union’s ERC Advanced Grant (contract 246939). [17] A.Peshkov,S.Ngo,E.Bertin,H.Chaté, andF.Ginelli, We also acknowledge support by the Max Planck Insti- Phys. Rev. Lett. 109, 098101 (2012). tuteforthePhysicsofComplexSystemsinDresden,Ger- [18] P. Romanczuk, I. D. Couzin, and L. Schimansky-Geier, many, through the Advanced Study Group “Statistical Phys. Rev. Lett. 102, 010602 (2009). [19] D. Grossman, I. S. Aranson, and E. B. Jacob, New J. Physics of Collective Motion”, where part of this work Phys. 10, 023036 (2008). was conducted. [20] B. Szabó, G. J. Szo¨llo¨si, B. Go¨nci, Z. Jurańyi, D. Selmeczi, and T. Vicsek, Phys. Rev. E 74, 061908 (2006). [21] S. Henkes, Y. Fily, and M. C. Marchetti, Phys. Rev. E 84, 040301(R) (2011). [1] C. Breder, Ecology 35, 361 (1954). [22] G. Grégoire, H. Chaté, and Y. Tu, Physica D 181, 157 [2] A. Okubo, Adv. Biophys. 22, 1 (1986). (2003). [3] D.Gru¨nbaumandA.Okubo,FrontiersinTheoreticalBi- [23] G. Grégoire and H. Chaté, Phys. Rev. Lett. 92, 025702 ology,LecturenotesinBiomathematics, 100,296(1994). (2004). [4] J. A. Shapiro and M. Dworkin, eds., Bacteria as multi- [24] A.FetterandJ.Walecka,Theoretical Mechanics of Par- cellular organisms (Oxford University Press, 1997). ticles and Continua, Dover Books on Physics Series [5] I. D. Couzin, J. Krause, R. James, G. D. Ruxton, and (Dover, 2003). N. R. Franks, J. Theor. Biol. 218, 1 (2002). [25] M.AldanaandC.Huepe,J.Stat.Phys.112,135(2003). [6] D.J.T.Sumpter,CollectiveAnimalBehavior (Princeton [26] K. Binder, Rep. Prog. Phys. 60, 487 (1997). University Press, 2010). [27] K. Ogata, Modern Control Engineering (4th Edition) [7] A. Deutsch, G. Theraulaz, and T. Vicsek, Interface Fo- (Prentice Hall, 2001). cus 2, 689 (2012). [28] N.D.MerminandH.Wagner,Phys.Rev.Lett.17,1133 [8] T. Vicsek and A. Zafeiris, Phys. Rep. 517, 71 (2012). (1966). [9] M.Brambilla,E.Ferrante,M.Birattari, andM.Dorigo, [29] X. Trepat, M. R. Wasserman, T. E. Angelini, E. Millet, Swarm Intelligence (accepted for publication). D.A.Weitz,J.P.Butler, andJ.J.Fredberg,Nat.Phys. [10] T. Vicsek, A. Cziro´k, E. Ben-Jacob, I. Cohen, and 5, 426 (2009). O. Shochet, Phys. Rev. Lett. 75, 1226 (1995). [30] A. Cavagna, I. Giardina, A. Orlandi, G. Parisi, A. Pro- [11] A. Cziro´k, H. E. Stanley, and T. Vicsek, J. Phys. A: caccini,M.Viale, andV.Zdravkovic,AnimalBehaviour Math. Gen. 30, 1375 (1997). 76, 217 (2008). [12] H.Chaté,F.Ginelli,G.Grégoire, andF.Raynaud,Phys. [31] Y. Katz, K. Tunstrom, C. C. Ioannou, C. Huepe, and Rev. E 77, 046113 (2008). I.D.Couzin,P.Natl.Acad.Sci.USA108,18720(2011). [13] J. J. Binney, N. J. Dowrick, A. J. Fisher, and M. E. J.