Stability of milling patterns in self-propelled swarms on surfaces
Jason Hindes, Victoria Edwards, Sayomi Kamimoto, George Stantchev, Ira B. Schwartz
SStability of milling patterns in self-propelled swarms on surfaces
Jason Hindes , Victoria Edwards , Sayomi Kamimoto , George Stantchev , and Ira B. Schwartz U.S. Naval Research Laboratory, Washington, DC 20375, USA and Department of Mathematics, George Mason University, Fairfax Virginia, 22030, USA
In some physical and biological swarms, agents effectively move and interact along curved surfaces.The associated constraints and symmetries can affect collective-motion patterns, but little is knownabout pattern stability in the presence of surface curvature. To make progress, we construct ageneral model for self-propelled swarms moving on surfaces using Lagrangian mechanics. We findthat the combination of self-propulsion, friction, mutual attraction, and surface curvature producemilling patterns where each agent in a swarm oscillates on a limit cycle, with different agents splayedalong the cycle such that the swarm’s center-of-mass remains stationary. In general, such patternsloose stability when mutual attraction is insufficient to overcome the constraint of curvature, and weuncover two broad classes of stationary milling-state bifurcations. In the first, a spatially periodicmode undergoes a Hopf bifurcation as curvature is increased which results in unstable spatiotemporaloscillations. This generic bifurcation is analyzed for the sphere and demonstrated numerically forseveral surfaces. In the second, a saddle-node-of-periodic-orbits occurs in which stable and unstablemilling states collide and annihilate. The latter is analyzed for milling states on cylindrical surfaces.Our results contribute to the general understanding of swarm pattern-formation and stability in thepresence of surface curvature, and may aid in designing robotic swarms that can be controlled tomove over complex surfaces.
I. INTRODUCTION
Swarming occurs when emergent spatiotemporal pat-terns are produced from the interaction of large collec-tives of coupled mobile agents with limited dynamicsand simple rules. Examples have been discovered overa wide range of space and time scales in nature includ-ing: colonies of bacteria and bees[1, 2], swarms of antsand locusts[3, 4], schools of fish[5–7], flocks of starlingsand jackdaws[8–10], and crowds of people[11]. As a con-sequence, much work in biophysics has focused on study-ing simple models to understand general swarm pattern-formation[12–15], and the nonequilibirum statistics ofactive-matter systems[16–19].In general, natural swarms are robust to the loss ofindividual agents and communication links, and are dy-namically responsive to complex, changing environments.Because of these useful properties, there is great inter-est in engineering multi-robot systems that can imitatenature’s ability to produce robust collective behaviorfrom simple components. Research in robotic swarm-ing has focused on developing control laws that gov-ern individual robot motion within a group so that sta-ble group-dynamics emerges[20–25]. In addition to pat-tern formation, robotic swarms have been proposed forsolving a variety of cooperative-dynamics tasks includ-ing mapping[26], resource allocation [27–29], and leaderfollowing[30].Swarms in three-spatial dimensions can be subject tothe geometry and curvature of surfaces on which theymove. Examples range from embryonic development[31]and corneal growth[32] to the collective dynamics of ac-tive particles confined to droplets[33, 34] and roboticconsensus and control on manifolds[35] (e.g., track-ing and sensing on a surface with known geometry).Yet, most of the known theoretical results concerning swarms on surfaces pertain to novel steady-state pat-terns due to curvature, equilibrium statistics, and re-laxation properties[36–41]. Despite extensive work ongeneral swarm stability[42–44], to our knowledge littleis known about the detailed bifurcation structure of col-lective swarming patterns on various curved surfaces.Here, we analyze stationary, milling patterns of self-propelled swarms on surfaces with simple attractive inter-actions. We find two types of bifurcations that appear assurface curvature is increased: a generalized Hopf bifur-cation where spatially-periodic modes introduce unstableoscillations at a definite frequency, and a saddle-node-of-periodic-orbits (SNpo) bifurcation where stable and un-stable milling states coalesce. Our results shed light onhow surface curvature destabilizes swarming patterns inboth generic and geometry-specific ways, and provides in-sight into how patterns emerge on different surfaces. Inaddition to general understanding, the techniques used toreveal bifurcation structure may aid in designing roboticswarms to move cooperatively over complex surfaces andterrains.To begin, let us consider a swarm of self-propelledagents interacting through conservative forces – the sim-plest model for position-dependent interactions. In ad-dition to such forces, each agent is acted upon byactivation-dissipation forces that depend on its velocity, A l = µ ( v l ) v l /v l , where v l is the velocity of the l th agent, v l is its speed, and µ ( v l ) is a nonlinear function to bespecified[45–47]. The latter is a consistent feature of ac-tive matter systems, since the “friction coefficient”, µ ( v l ),can be both positive and negative, due to input energyfrom an agent’s environment[45–48] (or onboard motorin the case of mobile robots[49, 50]).We are interested in studying a swarm’s dynamicsin generalized curvilinear coordinates, q kl , in which La- a r X i v : . [ n li n . AO ] J u l grange’s equations-of-motion take the covariant form ddt ∂L∂ ˙ q kl − ∂L∂q kl + ∂F∂ ˙ q kl = 0 , (1)where L is the swarm’s Lagrangian, F = (cid:80) l (cid:82) v l µ l ( v l ) dv l is the Rayleigh dissipation (activation) function[47], anddots denote time derivatives. Equation(1) is standardfor Lagrangian mechanics with dissipation – the impor-tant distinction being the possibility of effectively neg-ative dissipation for active agents. The Lagrangian ap-proach is expected to significantly simplify the descrip-tion and analysis of collective-motion states in swarms,especially in the presence of symmetries and generalizedconstraints, as in classical mechanics problems. More-over, it does not assume an over-damped limit for theagent dynamics[36, 40].For simplicity, in this paper we take F = − (cid:88) n v n (cid:16) α − β v n (cid:17) , (2) L = m (cid:88) n v n − a N (cid:88) n,j | r n − r j | , (3)as in many works on swarms of self-propelled agents inCartesian coordinates[45, 46, 48, 51, 52] (which we re-cover as a special case). In Eqs. (2-3) m is the mass ofeach agent, α is a self-propulsion constant, β is a fric-tion constant, a is a coupling constant, N is the numberof agents, and r j = ( x j , y j , z j ) is the position-vector inCartesian coordinates for the j th agent in three spatialdimensions. Note that Eq.(3) implies that the pairwiseinteraction force between two agents, n and j is linearin r n − r j and global, though these assumptions can berelaxed[50, 53]. In addition, note that the spring-likepotential in Eq.(3) is defined in terms of the pairwisedistance between agents in Euclidean space, rather thanalong, e.g., geodesics [54]. This is done for the sake ofsimplicity and to allow for a proof of principle in theweak mean-curvature limit. Beyond theory, we note thatCartesian versions of Eqs.(2-3) have been implementedin experiments with several robotics platforms includingautonomous cars, boats, and quad-rotors[49, 50].In what follows we study the stability of swarms onsurfaces using the above formalism. In Sec.II, we definethe milling patterns of interest as swarm limit-cycles, andcompute their properties on several surfaces. In Sec.III,we analyze the stability of limit cycles by tracking whichFloquet multipliers cross the unit circle under parametervariation. Two destabilization patterns are observed assurface curvature is increased. In the first, two milling-state multipliers cross in a Hopf bifurcation of spatially-periodic modes. In the second, many multipliers cross si-multaneously, resulting in a SNpo bifurcation. The firstis discussed in Sec.III A, while the second is discussed inSec.III B. Throughout, predictions are compared to sim-ulations. Sec.IV provides a summary and further discus-sion. II. LIMIT-CYCLE MILLING
We are interested in the dynamics of Eqs.(1-3) whenagents are constrained to surfaces, e.g., when a partic-ular curvilinear coordinate is constant. Classic exam-ple surfaces considered in detail are: the sphere of ra-dius r ( x l , y l , z l ) s = r (sin θ l cos φ l , sin θ l sin φ l , cos θ l ); thecylinder of radius ρ , ( x l , y l , z l ) c = ( ρ cos φ l , ρ sin φ l , z l ),and the torus of radii b and c , ( x l , y l , z l ) t =(( c + b cos θ l )cos φ l , ( c + b cos θ l )sin φ l , b sin θ l ). However,our approach can be generally applied to any coordinatesurface.Substituting these curvilinear coordinates into Eqs.(1-3) gives the example equations-of-motion for swarmingon a sphere,¨ φ l sin θ l − ˙ φ l sin θ l (cid:104) α − βr (cid:0) ˙ θ l + ˙ φ l sin θ l (cid:1)(cid:105) − φ l ˙ θ l cos θ l − aN (cid:88) j sin θ j sin( φ j − φ l ) = 0 , (4)¨ θ l − ˙ θ l (cid:104) α − βr (cid:0) ˙ θ l + ˙ φ l sin ˙ θ l (cid:1)(cid:105) − ˙ φ l sin θ l cos θ l − aN (cid:88) j (cid:104) sin θ j cos θ l cos( φ j − φ l ) − cos θ j sin θ l (cid:105) = 0 , (5)and on a cylinder,¨ z l = ˙ z l (cid:104) α − β (cid:0) ρ ˙ φ l + ˙ z l (cid:1)(cid:105) + aN (cid:88) j ( z j − z l ) , (6)¨ φ l = ˙ φ l (cid:104) α − β (cid:0) ρ ˙ φ l + ˙ z l (cid:1)(cid:105) + aN (cid:88) j sin (cid:0) φ j − φ l (cid:1) , (7)where we set m = 1 (since it entails an overall time-constant only). Similar equations for swarming on a toruscan be found in the App.V A. We note that the sameequations-of-motion can be derived from the Cartesianversion of Newton’s law for the model considered, by sub-stitution, though this approach is more cumbersome[44].There are two tendencies in the swarm’s dynamicsworth noting based on the physics of Eqs.(1-3). First,when a = 0 the activation-dissipation forces drive agentsto follow independent geodesics on the surface withequal speed, but no particular direction of motion (seeApp.V B). Second, agents tend to cluster in space, suchthat the distance between agents is minimized. Thelatter is implied by the spring-like interaction forces inEq.(3). Together the combined forces produce a varietyof collective-motion states depending on parameter val-ues and initial conditions.However, in this work we concern ourselves with sta-tionary milling solutions of Eqs.(1-3) where each agentmaintains periodic motion, while the center-of-mass ofthe swarm, R ≡ (cid:80) j r j /N , remains constant. A primaryreason for studying milling states is that they emergefrom broad initial conditions, e.g., spatially uniform dis-tributions of agents over a surface with random initial ve-locities. In contrast, flocking states require initial align-ment of velocities. The periodic motion of milling statesoccurs on limit cycles (LCs), with agents splayed approx-imately uniformly at different points along the same LC.The top panels in Fig.1 show snap-shots in time of ex-ample milling states for large swarms.In general, we can only calculate exact LCs for specialcases such as on the sphere. However, we can computeLCs on other surfaces, without resorting to large multi-agent simulations, by using a self-consistency approach inthe limit N (cid:29)
1. The approach entails trading the sums-over-agents that appear in Eqs.(4-7), with time-averagesover a single-particle LC. The former requires simulat-ing a 4N-dimensional system, the latter a 4-dimensionalsystem with a properly chosen R , R = (cid:90) T ( R )0 r ( LC ) ( t ; R ) dtT ( R ) . (8)The trajectory r ( LC ) ( t ; R ) is a LC of the single-particlesystem, and T ( R ) is its period. A single-particle systemcan be found by replacing R , which appears in the in-teraction sums in Eqs.(4-7), with Eq.(8). For example,substituting Eq.(8) into Eqs.(6-7) we find¨ z = ˙ z (cid:104) α − β (cid:0) ρ ˙ φ + ˙ z (cid:1)(cid:105) − az, (9)¨ φ = ˙ φ (cid:104) α − β (cid:0) ρ ˙ φ + ˙ z (cid:1)(cid:105) − a sin φ (cid:104) cos φ (cid:105) , (10)where (cid:104) cos φ (cid:105) = (cid:90) T cos φ ( t ) dtT . (11)is the time average of cos φ over a LC period. Note thatin Eqs.(9-11) we have chosen (cid:104) z (cid:105) = (cid:104) sin φ (cid:105) = 0, since theseconstants simply shift and rotate LCs along the cylinder,and we have dropped the agent subscript, l . See App.V Afor the analogous equations on a torus.Limit cycles can be computed numerically by: inte-grating the single-particle system[55], r ( t ), with an initialchoice of R , determining the period of the LC after tran-sients, and updating the choice of R based on a quasi-Newton evaluation of Eq.(8). Such an algorithm can beimplemented using numerical-integration software com-bined with a fixed-point-solver. Example comparisonsare shown in Fig.1 with excellent agreement. In the toppanels we plot the LCs in red on each surface comparedto large multi-agent simulations shown with magenta cir-cles. The initial conditions for the large multi-agent simu-lations were uniformly-random spatial distributions overeach surface with random velocities. In the bottom pan-els, we plot solutions of Eq.(8) compared to averages-over-particles in the swarm, shown with blue circles, fora wide range of parameters.We note that in general when milling states are sta-ble in swarms with attractive interactions, for a given R there are two possible stable LCs that are simply re-flections of one another, e.g., φ → − φ given (cid:104) z (cid:105) = (cid:104) y (cid:105) = 0.Both LCs can be seen clearly in the top panels of Fig.1(b-c). As a consequence, a general milling state on a surface is built from some fraction of the agents moving on oneLC, while the remaining move on the reflected LC– thefraction depending on initial conditions. The existence oftwo LCs also occurs in Cartesian swarms in 2d[48, 51, 52],where reflection simply reveres the angular velocity.For the special case of milling on the sphere, LCs aresimply circular orbits. Assuming a swarm milling with auniform distribution of constantly rotating φ j ’s all in thesame direction, and oriented such that the polar angle isconstant, we find θ j ( t ) = sin − (cid:32)(cid:114) αβar (cid:33) , φ j ( t ) = 2 π ( j − N + √ at. (12)Equation (12) is easy to check by direct substitution, andcompares well to simulations, Fig.1(a). Note that for con-sistency, a LC does not exist on the sphere if α/ [ βar ] > φ j = −√ a . III. STABILITY
As an estimate, we expect milling states to change sta-bility when the arc length of the LC is roughly equalto the inverse of the mean surface curvature. If thesetwo quantities are very different, then a periodic-solutionshould be hard to realize on the surface. The naturalperiod of oscillations for milling scales as ∼ / √ a , e.g.Eq.(12), while the average speed of an agent is (cid:112) α/β ;the latter is the speed at which self-propulsion and fric-tion forces cancel. See App.V B for more details. Alto-gether we expect instability to arise approximately when H (cid:112) α/ ( βa ) ∼
1, where H is the mean surface curvature.Beyond this crude estimate, we would like to under-stand and analyze the stability of milling states quanti-tatively, and determine if there are any differences in thebifurcations between surfaces. In our stability analysisbelow, we focus on milling states that correspond to asingle LC, where all agents rotate in the same direction,for three reasons: this case persists when repulsive forcesare added, the stability of any given milling state has onlya weak dependence on the number of agents on each LC,and it is analytically tractable. We begin by analyzingthe linear stability of milling states on the sphere, wherea complete bifurcation picture can be derived, and com-pare to numerical Floquet analysis for other surfaces.To determine the local stability on the sphere weneed to understand how small perturbations to Eq.(12)change in time. Our first step is to substitute a gen-eral perturbation, θ j ( t ) = sin − (cid:104)(cid:112) α/ ( βar ) (cid:105) + B j ( t ) and φ j ( t ) = 2 π ( j − /N + √ at + A j ( t ), into Eqs.(4-5) andcollect terms to first order in A j ( t ) and B j ( t ) (assuming | A j | , | B j | (cid:28) ∀ j ). The result is the following linear sys-tem of differential equations with constant coefficients –the latter property is a consequence of our transforma- (a) (b) (c) xyzxyz xyz -2 2-2 2-11-1 11-1-11 -0.8 0.60.5-0.5-0.52 0.5 1 1.5 20.920.960.940.98 0.6 0.8 1.2 1.4 1.610.960.970.980.99 FIG. 1. (Color online) Milling states on coordinate surfaces. Top panels show time-snapshots from simulations with N = 600agents on the (a) sphere ( r = 1 . , α = 0 . , β = 5 . , a = 0 . ρ = 0 . , α = 0 . , β = 1 , a = 2), and (c) torus( b = 0 . , c = 1 . , α = 0 . , β = 1 . , a = 1 . tion into the proper coordinate system and is what allows for an analytical treatment for milling:¨ B l + αβr B l − (cid:32) αβ √ ar (cid:113) βar α − (cid:33) ˙ A l = αβr N (cid:88) j (cid:34)(cid:32) (cid:16) βar α − (cid:17) cos (cid:16) π ( j − l ) N (cid:17)(cid:33) B j − sin (cid:16) π ( j − l ) N (cid:17)(cid:113) βar α − · A j (cid:35) , (13)¨ A l + 2 α ˙ A l + 2 √ a (cid:113) βar α − · (cid:2) ˙ B l + αB l (cid:3) = aN (cid:88) j (cid:34) sin (cid:16) π ( j − l ) N (cid:17)(cid:113) βar α − · B j + cos (cid:16) π ( j − l ) N (cid:17) A j (cid:35) . (14) A. Spatially periodic modes
Given the periodicity implied by the splayed phasevariables in Eq.(12), it is natural to look for eigen-solutions of Eqs.(13-14) in terms of the discrete Fouriertransforms of A j ( t ) and B j ( t ). In fact, by inspection wecan see that only the first harmonic survives the summa-tions on the right-hand sides of Eqs.(13-14), because ofthe sine and cosine terms, and hence we look for the par-ticular solutions: A j ( t ) = A exp( λt + 2 πi ( j − /N ) and B j ( t ) = B exp( λt + 2 πi ( j − /N ). Substitution gives thefollowing equation for the stability exponent, λ , of the spatially-periodic modes on the sphere: (cid:34) λ + 2 αλ − a (cid:35) · (cid:34) λ + a (cid:16) αβar − (cid:17)(cid:35) − a (cid:16) − αβar (cid:17) · (cid:34) λ − i √ a (cid:35) · (cid:34) i √ a − α − λ (cid:35) = 0 . (15)Note that the complex-conjugate of λ is also a solutionof Eqs.(13-14).In general, the milling state on the sphere is linearlystable if there are no solutions to Eq.(15) with Re [ λ ] >
0. By dividing λ by α in Eq.(15), we can see that thespectrum for the spatially periodic mode is effectively afunction of two positive parameters only, µ ≡ a/α and ν ≡ α/ ( βar ), i.e. λ/α = L ( µ, ν ) and has four solutions.Moreover, in the relevant parameter region ν < ν = 1 in Sec.III B.Generically, a Hopf bifurcation occurs when λ = ± iω ,with ω (cid:54) = 0. It is easy to verify that this condition[59] issatisfied in Eq.(15) when ν H = 38 (16)Importantly, the Hopf bifurcation uncovered is acurvature-induced instability; it does not occur in flatspace. For instance, if one takes r → ∞ , Eq.(16) can onlybe satisfied in the trivial infinite-speed or zero-couplinglimits. Also, it is worth mentioning that since Eqs.(13-14)describe the linearized dynamics in co-moving referenceframes, the Hopf bifurcations entailed by Eq.(16) are gen-eralized Hopf bifurcations of LCs (leading to motion ona high-dimensional torus) . However, for brevity and toavoid double usage of “torus”, we use the short-hand,Hopf, throughout.Example Hopf curves are shown in Fig.2(a) withdashed lines for several values of α . For each value of α , the milling state is predicted to be stable for all largervalues of a (i.e., moving to the right of the dashed-lines atfixed r ). The circles denote simulation-determined tran-sition points: the smallest a ( r ) for which a swarm of 600agents, initially prepared in the milling state with a smallrandom perturbation (i.e., independent and uniformlydistributed A j and B j over [ − − , − ]), returns to amilling state after an integration time of t = 10000. Pre-dictions from Eq.(16) show excellent agreement with sim-ulation results. Similarly determined transition pointsfor a milling state in which half the agents rotate in onedirection, and half rotate in the opposite direction, areshown with squares. We can see that the Hopf-transitionline still gives a good approximation for this more generalcase, especially for larger values of r .So far we have demonstrated that the spatially peri-odic modes that determine milling stability on the sphereare simple plane waves. On other surfaces the modes areperiodic but not plane waves. Still, Hopf bifurcations oc-cur, when such modes have complex-valued Floquet mul-tipliers (FMs) that cross the unit circle. Floquet multi-pliers are the LC-analogs of the eigenspectrum for smallvariation around fixed-points in dynamical systems[56–58, 60]. Unlike the sphere, however, the FMs and theirHopf-points have to be determined numerically. In fact,for the Hopf bifurcations on the cylinder and torus wemust compute the FMs for the whole multi-agent system;See App.V C for further details. Examples for the cylin-der are shown in Fig.3. In panel (a), we plot the stabilitydiagram for several values of α . The dashed-lines are nu-merically determined Hopf bifurcations. For each valueof α , the milling state is predicted to be stable for alllarger values of a (i.e., moving to the right of the dashed-lines at fixed ρ ). The circles are simulation-determined I II (a) (b)
FIG. 2. (Color online) Milling stability on the sphere. (a)Hopf bifurcation curves are drawn with dashed lines for α = 0 .
05 (green, bottom), α = 0 .
20 (blue, middle), and α = 0 .
40 (red, top). In each case β = 5 .
0. Points denotesimulation-determined stability changes for N = 600: millingwith all agents rotating in the same direction (circles), andmilling with half the agents rotating in the opposite direction(squares). (b) Stability diagram. Region (I, red) has no un-stable modes, and region (II, blue) has two unstable modes.(I) and (II) are separated by the Hopf bifurcation, Eq.(16),drawn with a solid black line. The milling state exists to theleft of the dashed-line. transition points– computed in the same manner as forthe sphere. As with the sphere, the LC linear-stabilityanalysis agrees well with simulations. An example peri-odic mode that changes stability in a Hopf bifurcationis shown in panel (b). The mode is an eigenvector cor-responding to one of the two FMs that crosses the unitcircle. Plotted are the real and imaginary parts of thevector’s z l and φ l -components. Similar results are shownfor Hopf bifurcations on the torus in Sec.V A. (a) (b) (c) -0.6 -0.2 0.2 0.6-0.40-0.4 NN/2 (i)(ii)(iii)(iv) FIG. 3. (Color online) Milling stability on the cylinder.(a) Hopf bifurcation curves are drawn with dashed lines for α = 4 . α = 2 . α = 0 . α = 2 . β = 1 .
0. Points denote simulation-determined stability changes for N = 600, and follow the samecolor scheme. (b) Spatially periodic mode that changes stabil-ity in the Hopf bifurcation. Plotted are the mode componentsversus the agent number, l when a = 0 . ρ = 1 . α = 2, and β = 1: Re [ z l ] (solid blue, iii), Im [ z l ] (dotted blue, ii), Re [ φ l ](dashed red, iv), and Im [ φ l ] (dotted-dashed red, i). (c) Sta-ble (solid red, outer) and unstable (dashed blue, inner) limitcycles that collide in the saddle-node-of-periodic-orbits bifur-cation. Parameters are a = 2, ρ = 0 . α = 0 .
5, and β = 1. B. Repeated-eigenvalue modes
We have seen that the Hopf bifurcation of spatiallyperiodic modes is qualitatively the same for millingstates on different surfaces. A natural question concernswhether other types of stability changes occur, not dueto a single complex pair of FMs, but the rest of the spec-trum. To make progress, it is again useful to consider thesphere and the remaining solutions of Eqs.(13-14), first.After which, we can take a closer look at the cylinder,where a different bifurcation is possible.As mentioned previously, higher Fourier-modes havesums in Eqs.(13-14) that vanish. Consequently, Eqs.(13-14) have N − s ,satisfying sα (cid:34)(cid:16) sα (cid:17) + aα αβar (cid:35)(cid:34)(cid:16) sα (cid:17) + 2 (cid:35) + sα aα (cid:34) − αβar (cid:35)(cid:34)(cid:16) sα (cid:17) + 1 (cid:35) = 0 . (17)Just as with Eq.(15), there are two effective parametersand four solutions to Eq.(17)– giving a total of 4( N − Re [ s ] ≤ ν = α/ [ βar ] < ν = α/ ( βar ) = 1, a degenerate Hopf bifurcationoccurs – degenerate, since 2( N −
1) FMs cross the unitcircle simultaneously . Putting the whole spectrum of themilling state on the sphere together, we can now drawthe complete stability diagram, shown in Fig.2(b). Thered region (I) has no unstable modes, and the blue region(II) has two unstable modes. The Hopf bifurcation forthe periodic mode is drawn with a solid black line, whilethe degenerate Hopf bifurcation of the repeated part ofthe spectrum is drawn with a dashed black line. It isclear from Fig.2(b), that the milling state on the spherechanges stability only through the Hopf bifurcation ofspatially-periodic modes. Accounting for the rest of thespectrum does not change the stability picture.However, this is not true in general for milling onother surfaces. In particular, on the cylinder we findthat N − real FMs can approach unity in a degener-ate SNpo bifurcation[56–58]. An example FM-spectrumnear bifurcation is shown in Sec.V C for reference.To track the SNpo numerically, it is sufficient to com-pute the FMs for a LC of the single-particle system,Eqs(9-11) – making the computation significantly fasterthan for the Hopf bifurcation of spatially-periodic modes.Examples of SNpo curves are shown in Fig.3(a) with solidlines. Simulation-determined transitions points are plot-ted with colored circles. As with the Hopf bifurcations,the milling state is predicted to be stable for all largervalues of a (i.e., moving to the right of the solid lines atfixed ρ ).
1. Stable and unstable cycle pairs
An important hallmark of SNpo bifurcations is the col-lision and annihilation of two periodic orbits[56–58]. Asa consequence, we expect a stable milling-state LC tocoalesce with an unstable LC, for swarms on the cylin-der, both of which have N − z ( t ) ( LC ) ≈ Z sin( ωt + γ ) and φ ( t ) ( LC ) ≈ ( Z/ρ ) sin ωt , withthe assumption that higher harmonics make only smallcontributions. Substituting our ansatz into Eqs(9-11),multiplying the equations by sin ωt (and cos ωt respec-tively), and integrating over a period, we obtain the fol-lowing three fixed-point equations: ω − a = 14 βZ ω sin 2 γ, (18) Z ω = 4 αβ (5 + cos 2 γ ) , (19) Zω ρa (cid:104) ω − aω (cid:105) = (cid:90) π cos (cid:16) Zρ sin σ (cid:48) (cid:17) dσ (cid:48) π · (cid:90) π sin (cid:16) Zρ sin σ (cid:17) sin σdσ π , (20)which can be solved for the approximate parameters ofboth stable and unstable milling states on the cylinder.Equations (18-20) have two non-trivial solutions thatmeet at a critical curvature, corresponding to the SNpobifurcation, as expected. Example solutions are shown inFig.4 for large (red), medium (blue), and small (green)coupling. The stable amplitude is shown with a solid lineand the unstable with a dashed line. The agreement be-tween the first-harmonic approximation and large-agentsimulations (points) is fair in general, and increases withthe coupling.It is clear from Fig.4 that as ρ is increased above theSNpo bifurcation, the amplitudes converge quickly totheir asymptotic values. The limiting solutions corre-spond to γ → , π/ Z = 2 α βa , αβa . (21)as ρ → ∞ . The asymptotic results imply an importantproperty of milling on the cylinder: unstable states ex-ist even when mean curvature is weak . The persistenceof such dynamically unstable orbits makes swarming on FIG. 4. (Color online) Saddle-node-of-periodic-orbits dia-gram for milling states on the cylinder. Predictions for thelimit-cycle amplitudes, Eqs.(18-20), are shown for the stable(solid) and unstable (dashed) cycles with a = 20 (red, lower), a = 5 (blue, middle), and a = 2 (green, top). Points denotesimulation amplitudes for a swarm of N = 600 agents. Param-eters are α = 0 . β = 1 . cylindrical surfaces qualitatively different from, e.g., thesphere, where no such orbits exist. IV. CONCLUSION
In this work we studied the stability of stationary,milling patterns in self-propelled swarms constrained tomove on surfaces. Using Lagrangian mechanics, wefound that with simple attractive forces between agentssuch patterns correspond to surface-dependent limit cy-cles where agents in a swarm are splayed along differentpoints on the same cycle. We showed that the constraintof curvature can destabilize milling swarms, and does soin both generic and geometry-specific ways. In the for-mer, a spatially periodic mode of the milling state be-comes unstable in a generalized Hopf bifurcation. Thisbifurcation was demonstrated and analyzed for swarm-ing on the sphere, cylinder, and torus, as examples. Onthe other hand, we found that on cylindrical surfaces un-stable milling patterns exist, even in the limit of smallcurvature, and can merge with stable milling states in asaddle-node-of-periodic-orbits bifurcation. Our analyti-cal and numerical results were verified in detail with largemulti-agent simulations.To our knowledge, these are the first formal bifurca-tion results for swarming patterns in general geometries,and they hint at a more general theory for swarming sta-bility on Riemannian manifolds, especially in the limitof weak mean-curvature. Extending the Lagrangian for-mulation of the swarm dynamics in terms of geodesicpotentials on arbitrary surfaces, as well as a geodesicallydefined center-of-mass that is constrained to the surface,will allow for a more accurate analysis of milling pat-terns on surfaces with higher sectional curvature vari-ability. We expect our methods to be useful for ana-lyzing patterns in other swarming problems where boththe patterns themselves and the constraints imposed are expressible in terms of symmetries and generalized coor-dinates. For instance, an important future generalizationof the swarms studied here is the inclusion of repulsiveforces, since such forces are known to create a variety ofmilling patterns on surfaces beyond the simple splayedlimit cycles that we analyzed.Considering the patterns of swarm dynamics on spe-cific classes of surfaces used in practical applications isanother possible direction of future research. For in-stance, piece-wise polynomial surfaces, such as splines,or more generally rational splines can be of interest forground robotics applications. Parametrically control-ling swarms of mobile robots to cooperatively move overcomplex surfaces and terrains in the future will requirestudying the effects of communication topology, both dy-namic and incomplete, as well as uncertainty and noise.Of course, many open questions remain. Yet, our ap-proach shows that general stability analysis for large self-propelled swarms on surfaces is possible and we demon-strate how it might be done.JH and IBS were supported by the U.S. Naval Re-search Laboratory funding (N0001419WX00055) andthe Office of Naval Research (N0001419WX01166),(N0001419WX01322), and (N0001420WX01237). TEwas supported through the U.S Naval Research Labo-ratory Karles Fellowship. SK was supported throughthe GMU Provost PhD award as part of the IndustrialImmersion Program. GS was supported through theOffice of Naval Research funding (N0001420WX01237).
V. APPENDIXA. Swarming on the torus
As defined in Sec.II, the torus has two radial pa-rameters: b and c . The former is the tube ra-dius, while the latter is the radius from the cen-ter of the hole to the center of the tube. Us-ing the coordinate transformation, ( x l , y l , z l ) t =(( c + b cos θ l )cos φ l , ( c + b cos θ l )sin φ l , b sin θ l ), where θ isan azimuthal angle and φ is the angle inside of the tube,we find the swarming Eqs.(1-3):¨ φ l = ˙ φ l (cid:34) α − β (cid:16) b ˙ θ l + ( c + b cos θ l ) ˙ φ l (cid:17)(cid:35) + aN ( c + b cos θ l ) (cid:88) j ( c + b cos θ j )sin( φ j − φ l ) +2 b sin θ l ( c + b cos θ l ) ˙ φ l ˙ θ l , (22)¨ θ l = ˙ θ l (cid:34) α − β (cid:16) b ˙ θ l + ( c + b cos θ l ) ˙ φ l (cid:17)(cid:35) + aN (cid:88) j sin( θ j − θ l ) − ( c + b cos θ l ) sin θ l b ˙ φ l + a sin θ l N b (cid:88) j (1 − cos( φ j − φ l ))( c + b cos θ j ) . (23)Furthermore the single-particle system, used to com-pute LCs associated with milling, can be found by sub-stituting Eq.(8) into Eqs.(22-23). By choosing a solutionwith (cid:104) z (cid:105) = (cid:104) y (cid:105) = 0, the single-particle system only dependson the time-average of x , (cid:104) x (cid:105) = (cid:104) ( c + b cos θ )cos φ (cid:105) , or¨ φ = ˙ φ (cid:34) α − β (cid:16) b ˙ θ + ( c + b cos θ ) ˙ φ (cid:17)(cid:35) +2 b ˙ φ ˙ θ sin θ − a (cid:104) x (cid:105) sin φc + b cos θ , (24)¨ θ = ˙ θ (cid:34) α − β (cid:16) b ˙ θ + ( c + b cos θ ) ˙ φ (cid:17)(cid:35) + ab (cid:0) c −(cid:104) x (cid:105) cos φ (cid:1) sin θ − ( c + b cos θ ) sin θb ˙ φ . (25)Equations (24-25) were used for the theory comparisonsin Fig.1(c).Note that if we take c → θ j → θ j − π/ ∀ j , inEqs.(22-23), we get the spherical system, Eq.(4-5), backwith r = b . Consequently, for relatively small c we expectto qualitatively reproduce the stability picture for millingstates on the sphere. Namely, such states will changestability through Hopf bifurcations of spatially-periodicmodes, as discussed in Sec.III A for the sphere and cylin-der. Three examples of Hopf-curves that demonstratethis prediction are shown in Fig.5, and agree well withlarge multi-agent simulations. B. Uncoupled dynamics
In the uncoupled limit ( a = 0) Eq.(1) becomes: ddt (cid:34) v ∂v∂ ˙ q k (cid:35) = v ∂v∂ ˙ q k (cid:16) α − βv (cid:17) . (26)Note that we have dropped the agent subscript, l . Fromgeneric initial conditions, we find that the forces on the FIG. 5. (Color online) Milling stability on the torus. Hopfbifurcation curves are drawn with solid lines for α = 0 . α = 0 . α = 0 . b = 1 . β = 5 .
0. Points denote simulation-determined stabilitychanges for N = 600, and follow the same color scheme. right-hand-side of Eq.(26) vanish as t → ∞ , resulting ina constant speed v = (cid:112) α/β . By Hamilton’s principle, weknow that ddt (cid:2) v ∂v∂ ˙ q k (cid:3) = 0, describes trajectories that ex-tremize the integral S [ q , ˙q ] = (cid:82) v dt , which is the actionfor a free-particle with speed v . Since v = (cid:112) α/β = dsdt ,where ds is the differential arc length along a surface, S = 12 (cid:114) αβ (cid:90) ds. (27)Equation (27) implies that, in the absence of interactionsand in the long-time limit, the integral which is extrem-ized along an agent’s trajectory is proportional to the arclength. Hence, such trajectories are surface geodesics. C. Numerical Floquet multipliers
Computing FMs numerically amounts to solving aneigenvalue problem for a linear, discrete-dynamical sys-tem. We start by considering the repeated part of theFloquet spectrum for milling states, since we can restrictourselves to the single-particle systems, e.g., Eqs.(6-7)and Eqs.(24-25). Recall that this part of the spectrum isrelevant for the SNpo bifurcations discussed in Sec.III B.Let us define a 4-dimensional phase-space vector (co-ordinates plus velocities) for the single-particle system, P = (cid:0) q , ˙ q , q , ˙ q (cid:1) . Note that the superscript is a labelfor the generalized coordinates (not an exponent), and wedrop the agent subscript l . Our goal is to compute thecharacteristic dynamics for small variations around LCs, (cid:15) ( t ) = P ( LC ) ( t ; R ) − P ( t ), where ( LC ) denotes the phase-space coordinates evaluated on a LC. Since the dynamicsare assumed to occur near LCs, for sufficiently small (cid:15) ( t ),we can employ Floquet’s theorem at lowest order in (cid:15) ( t ),which states that (cid:15) ( t + T ) = M (cid:15) ( t ) [56–58, 60], where M is a constant, non-singular monodromy matrix of dimen-sion 4x4 and T is the period of the LC.The first step is to compute the LC as described inSec.II for given parameters. The second step is to com- N 2N 3N 4N-101 (a) mode m u l t i p li e r N 2N 3N 4N01 (b) mode m u l t i p li e r FIG. 6. (Color online) Example Floquet multipliers formilling states on the cylinder, the full multi-agent system,Eqs.(6-7). (a) Multipliers near a Hopf bifurcation, a = 0 . ρ = 1 . α = 2, and β = 1. The real and imaginary parts of themultipliers are plotted with red and blue circles, respectively.The two multipliers that cross the unit circle at bifurcationare plotted with asterisks. (b) Multipliers near a SNpo bifur-cation, a = 3, ρ = 0 . α = 0 .
5, and β = 1. The repeated(real) multipliers that approach 1 at bifurcation are plottedwith a dashed line. pute M . The repeated FMs of the milling state aresimply the eigenvalues of M . A straightforward way todetermine M is to compute four perturbations, (cid:15) ( T ), (cid:15) ( T ), (cid:15) ( T ), and (cid:15) ( T ). Each can be found by integrat-ing the single-particle system over one period from theinitial conditions (cid:15) (0) = (cid:15) (1 , , , (cid:15) (0) = (cid:15) (0 , , , (cid:15) (0) = (cid:15) (0 , , , (cid:15) (0) = (cid:15) (0 , , , (cid:15) >
0. The integrated perturbations fillout the columns of M ; in particular M jn = [ (cid:15) n ( T )] j /(cid:15) ,where j, n ∈ { , , , } . For reference, (cid:15) = 10 − in all ofthe computations in this work. Lastly, to compute SNPobifurcation-points on the cylinder for fixed ρ , we simplyreduce a in small steps from some starting value until thesecond largest eigenvalue of M is equal to 1 (within somesmall error tolerance).On the other hand, for the Hopf bifurcations wemust compute the FMs for the full multi-agent LC,since the spatially-periodic modes that change stabil- ity entail collective oscillations of all agents. Keep inmind that the repeated part of the full spectrum willequal the single-particle FMs discussed in the previ-ous paragraph for large N . For reference, N = 300in all FM-computations in this work, which is largeenough to produce negligibly small finite-size effects.Our approach is essentially the same as for the single-particle system, except that monodromy matrix in thiscase is 4 N x4 N dimensional. Following the same pro-cedure, let us define the full phase-space vector X = (cid:0) q , ˙ q , q , ˙ q , q , ˙ q , q , ˙ q , ..., q N , ˙ q N , q N , ˙ q N (cid:1) . The firststep is to construct the complete milling state from thecomputed single-particle LC. The assumption of Eq.(8)implies that agents are splayed uniformly in time alonga LC, and therefore we take X ( t ; R ) ( LC ) = (cid:0) P ( LC ) ( t + [1 − T /N ; R ) , P ( LC ) ( t + [2 − T /N ; R ) , ..., P ( LC ) ( t + [ N − T /N ; R ) (cid:1) . (28)Next, we define the perturbation δ ( t ) = X ( t ; R ) ( LC ) − X ( t ), which has dynamics δ ( t + T ) = D δ ( t ), for small δ by Floquet’s theorem. The matrix D can be determinedby computing 4 N integrated perturbations, δ l ( T ), where δ l ( t = 0) = (cid:15) (0 , , ..., l , , ... ) with l ∈ { , , ..., N } . Incontrast to the preceding paragraph, these perturbationsare integrated in the full system, Eqs.(1-3). The inte-grated perturbations fill out the columns of D , such that D jl = [ δ l ( T )] j /(cid:15) , where j ∈ { , , ..., N } . As before, theFMs for the full system are the eigenvalues of D .Finally, to compute Hopf-bifurcation points on thecylinder or torus for fixed curvature, we simply reduce a in small steps from some starting value until a pair ofeigenvalues of D , not in the repeated part of the spec-trum, have unit absolute value (within some small errortolerance). Examples of the FMs for both milling-statebifurcations on the cylinder are shown in Fig.6 for thefull multi-agent system. As mentioned above, the N realFMs that approach 1 in panel (b), are identical to thesingle-particle FMs. [1] A. Polezhaev, R. Pashkov, A. I. Lobanov, and I. B.Petrov, Int. J. Dev. Bio. , 309 (2006).[2] J. Li and A. H. Sayed, EURASIP Journal on Advancesin Signal Processing , 18 (2012).[3] G. Theraulaz, E. Bonabeau, S. C. Nicolis, R. V.Sol´e, V. Fourcassi´e, S. Blanco, R. Fournier, J.-L. Joly, P. Fern´andez, A. Grimal, P. Dalle,and J.-L. Deneubourg, Proceedings of the Na-tional Academy of Sciences , 1 (2012).[5] K. Tunstrm, Y. Katz, C. C. Ioannou, C. Huepe, M. J.Lutz, and I. D. Couzin, PLOS Computational Biology , 1 (2013). [6] D. S. Calovi, U. Lopez, S. Ngo, C. Sire, H. Chat´e, andG. Theraulaz, New Journal of Physics , 015026 (2014).[7] A. Cavagna, L. Del Castello, I. Giardina, T. Grigera,A. Jelic, S. Melillo, T. Mora, L. Parisi, E. Silvestri,M. Viale, and A. M. Walczak, Journal of StatisticalPhysics , 601 (2015).[8] G. F. Young, L. Scardovi, A. Cavagna, I. Giardina, andN. E. Leonard, PLOS Computational Biology , 1 (2013).[9] M. Ballerini, N. Cabibbo, R. Candelier, A. Cav-agna, E. Cisbani, I. Giardina, V. Lecomte,A. Orlandi, G. Parisi, A. Procaccini, M. Viale,and V. Zdravkovic, Proceedings of the Na-tional Academy of Sciences T. A., and N. T. Ouellette, R. Soc. B. (2019),10.1098/rspb.2019.0865.[11] K. Rio and W. H. Warren, Transportation Research Pro-cedia , 132 (2014), the Conference on Pedestrian andEvacuation Dynamics 2014 (PED 2014), 22-24 October2014, Delft, The Netherlands.[12] T. Vicsek and A. Zafeiris, Phys. Rep. , 71 (2012).[13] M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B.Liverpool, J. Prost, M. Rao, and R. A. Simha, Rev.Mod. Phys. , 1143 (2013).[14] M. Aldana, V. Dossetti, C. Huepe, V. M. Kenkre, andH. Larralde, Phys. Rev. Letts. , 095702 (2007).[15] S. Chandra, M. Girvan, and E. Ott, Phys. Rev. X ,011002 (2019).[16] A. Solon, Y. Fily, A. Baskaran, M. E. Cates, Y. Kafri,M. Kardar, and J. Tailleur, Nature Phys. , 673 (2015).[17] E. Fodor, C. Nardini, M. E. Cates, J. Tailleur, P. Visco,and F. van Wijland, Phys. Rev. Lett. , 038103 (2016).[18] F. G. Woodhouse, H. Ronellenfitsch, and J. Dunkel,Phys. Rev. Lett. , 178001 (2018).[19] E. Woillez, Y. Zhao, Y. Kafri, V. Lecomte, andJ. Tailleur, Phys. Rev. Lett. , 258001 (2019).[20] J. P. Desai, J. P. Ostrowski, and V. Kumar, in IEEETransactions on Robotics and Automation , Vol. 17(6)(2001) pp. 905–908.[21] A. Jadbabaie, Jie Lin, and A. S. Morse, IEEE Transac-tions on Automatic Control , 988 (2003).[22] H. G. Tanner, A. Jadbabaie, and G. J. Pappas, in , Vol. 2 (2003) pp. 2016–2021Vol.2.[23] H. G. Tanner, A. Jadbabaie, and G. J. Pappas, in , Vol. 2 (2003) pp. 2010–2015Vol.2.[24] V. Gazi, IEEE Transactions on Robotics , 1208 (2005).[25] H. G. Tanner, A. Jadbabaie, and G. J. Pappas, IEEETransactions on Automatic Control , 863 (2007).[26] R. K. Ramachandran, K. Elamvazhuthi, and S. Berman,“An optimal control approach to mapping gps-denied en-vironments using a stochastic robotic swarm,” in RoboticsResearch: Volume 1 , edited by A. Bicchi and W. Bur-gard (Springer International Publishing, Cham, 2018) pp.477–493.[27] H. Li, C. Feng, H. Ehrhard, Y. Shen, B. Cobos, F. Zhang,K. Elamvazhuthi, S. Berman, M. Haberland, and A. L.Bertozzi, in (2017) pp. 4341–4347.[28] S. Berman, A. Halasz, V. Kumar, and S. Pratt,in
Proceedings 2007 IEEE International Conference onRobotics and Automation (2007) pp. 2318–2323.[29] M. A. Hsieh, ´A. Hal´asz, S. Berman, and V. Kumar,Swarm Intelligence , 121 (2008).[30] J. Wiech, V. A. Eremeyev, and I. Giorgio, ContinuumMechanics and Thermodynamics , 1091 (2018).[31] P. J. Keller, A. D. Schmidt, J. Wittbrodt, and E. H. K.Stelzer, Science , 1065 (2008).[32] J. Collinson, L. Morris, A. Reid, T. Ramaesh,M. Keighren, J. Flockhart, R. Hill, S. Tan, K. Ramaesh,B. Dhillon, and J. West, Dev. Dynam. , 432 (2002). [33] F. Keber, E. Loiseau, T. Sanchez, S. DeCamp, L. Giomi,M. Bowick, M. Marchetti, Z. Dogic, and A. Bausch,Science , 1135 (2014).[34] R. Zhang, Y. Zhou, M. Rahimi, and J. J. de Pablo, Nat.Commun. , 13483 (2016).[35] J. Markdahl, J. Thunberg, and J. Gonalves, IEEE Trans-actions on Automatic Control , 1664 (2018).[36] R. Sknepnek and S. Henkes, Phys. Rev. E , 022306(2015).[37] W. Li, Sci. Rep. , 13603 (2015).[38] S. Praetorius, A. Voigt, R. Wittkowski, and H. L¨owen,Phys. Rev. E , 052615 (2018).[39] L. M. C. Janssen, A. Kaise, and H. L¨owen, Sci. Rep. ,13603 (2015).[40] L. Apazaa and M. Sandoval, Soft Matter , 9928 (2018).[41] P. Castro-Villarreal and F. J. Sevilla, Phys. Rev. E ,052605 (2018).[42] V. Gazi and K. M. Passino, Swarm Stability and Opti-mization (Springer, New York, 2003).[43] G. Albi, D. Balagu, J. A. Carrillo, and J. von Brecht,SIAM J. Appl. Math. , 794 (2014).[44] J. Hindes, V. Edwards, S. Kamimoto, I. Triandaf, andI. B. Schwartz, Phys. Rev. E , 042202 (2020).[45] H. Levine, W. J. Rappel, and I. Cohen, Phys. Rev. E , 017101 (2000).[46] U. Erdmann, W. Ebeling, and A. S. Mikhailov, Phys.Rev. E , 051904 (2005).[47] E. Minguzzi, European Journal of Physics , 035014(2015).[48] M. R. D’Orsogna, Y. L. Chuang, A. L. Bertozzi, andL. S. Chayes, Phys. Rev. Lett. , 104302 (2006).[49] V. Edwards, P. deZonia, M. A. Hsieh, J. Hindes,I. Triandof, and I. B. Schwartz.[50] K. Szwaykowska, I. B. Schwartz, L. Mier-yTeran Romero, C. R. Heckman, D. Mox, and M. A.Hsieh, Phys. Rev. E , 032307 (2016).[51] E. Forgoston and I. B. Schwartz, Phys. Rev. E ,035203(R) (2008).[52] L. M. y Teran-Romero, E. Forgoston, and I. B. Schwartz,IEEE Transactions on Robotics , 1034 (2012).[53] J. Hindes, K. Szwaykowska, and I. B. Schwartz, Phys.Rev. E , 032306 (2016).[54] J. Osborne and G. Hicks, Notices of the American Math-ematical Society , 544 (2013).[55] In practice, we find that randomly generated initial con-ditions are sufficient for computing limit cycles in thesingle-particle system.[56] Y. A. Kuznetsov, Elements of Applied Bifurcation Theory (Springer, Berlin, 2004).[57] S. Strogatz,
Nonlinear Dynamics and Chaos: With Appli-cations to Physics, Biology, Chemistry, and Engineering (Westview Press, 2015).[58] S. Wiggins,
Introduction to Applied Nonlinear DynamicalSystems and Chaos (Springer, New York, 2003).[59] Note the agreement between Eq.(16) and the estimategiven in the first paragraph of Sec.III. The estimate be-comes exact by replacing 1 with 3 /
8. Similarly, if we re-place r with the inverse, mean curvature of the cylinder,2 ρ , in Eq.(16), we get an accurate approximation for theHopf bifurcation on the cylinder.[60] G. Moore, SIAM J. Numer. Anal.42