Collective motion dynamics of active solids and active crystals

We introduce a simple model of self-propelled particles connected by linear springs that describes a semi-rigid formation of active agents without explicit alignment rules. The model displays a discontinuous transition at a critical noise level, below which the group self-organizes into a collectively translating or rotating state. We identify a novel elasticity-based mechanism that cascades self-propulsion energy towards lower-energy modes as responsible for such collective motion and illustrate it by computing the spectral decomposition of the elastic energy. We study the model's convergence dynamics as a function of system size and derive analytical stability conditions for the translating state in a continuous elastic sheet approximation. We explore the dynamics of a ring-shaped configuration and of local angular perturbations of an aligned state. We show that the elasticity-based mechanism achieves collective motion even in cases with heterogeneous self-propulsion speeds. Given its robustness, simplicity and ubiquity, this mechanism could play a relevant role in various biological and artificial swarms.


Introduction
Collective motion (CM) is observed in a broad range of biological systems, including bird flocks, fish schools, herds of quadrupeds, insect swarms and groups of bacteria [1][2][3][4][5][6]. In recent years, these systems (referred to here, generically, as swarms) have been the subject of intense research [7][8][9][10]. In theoretical studies, several models have been introduced to analyse their dynamical properties, their ability to process information collectively, and the components that are necessary to produce CM. In parallel, research in control theory and robotics has developed various decentralized control algorithms that mimic the CM observed in nature to achieve a similar level of coordinated collective behaviour in groups of autonomous robots [9][10][11].
From the perspective of physics, swarms represent a unique class of systems that combine the familiar setup of an ensemble of trajectories in space with new dynamical properties and concepts. Indeed, on one hand swarms can be modelled as a collection of moving particles with imposed kinetic energy and effective forces that determine their motion. On the other hand, these interactions can be non-central, non-additive, non-local and non-distance-dependent, and we can typically identify no conserved quantity. Swarms also share a property that is characteristic of many living systems: energy is injected at the smallest scale, at the level of each individual agent or particle. These properties, together with a growing number of experiments and potential applications, have made CM an exciting new field of research at the interface between physics and biology. However, despite intense recent research activity, there is still no comprehensive understanding of the underlying mechanisms that can lead groups of self-propelled agents to self-organize and move in a common direction.
The prevailing paradigm in the theory of CM has been strongly influenced by the seminal work of Vicsek et al [12]. This paper introduced a minimal model for flocking, the Vicsek model, that has become a referent in the field [8][9][10]. It describes a group of point particles advancing at a fixed common speed, only coupled through alignment interactions that steer each agent towards the mean heading direction of all particles within a given radius [12,13]. As the amount of noise is decreased, the system undergoes a dynamical phase transition at a critical noise level, below which particles self-organize and start moving in a common direction. In this framework, a swarm can be viewed as a fluid of self-propelled spins with aligning interactions, described by an extension of the XY-model [14] where spins advance in their pointing direction rather than remaining affixed to a lattice. In the same spirit of the Vicsek model, but using a continuous description rather than agent-based simulations, Toner and Tu [15,16] introduced a hydrodynamic theory of active fluids. In a different context, for swarm robotics studies, several control rules have been designed to achieve CM by implementing consensus algorithms that converge to a common heading direction [11].
While most of the work on CM has focused on systems with aligning interactions, some studies have considered agents that do not explicitly exchange information on their orientation. In [17], for example, CM is driven by escape-pursuit interactions only; in [18], by inelastic collisions between isotropic agents; and in [19,20], by short-range radial forces that are coupled to the agents' turning dynamics. In [21], CM is achieved due to local pairwise repulsive interactions between soft deformable self-propelled spherical particles. Given that Vicsek-like algorithms rely on explicit alignment rules to achieve CM [12,[22][23][24], it was initially surprising that other models could self-organize without including such rules. One can argue that these models must include an indirect effect that produces an effective, implicit aligning interaction between individual agents in order to achieve CM, thus reducing the dynamics again to reaching consensus in the heading direction. An example of implicit aligning dynamics is found in [18], where inelastic collisions tend to align post-collision trajectories due to the conservation of momentum. Despite these studies, it remains unclear if the same underlying mechanism is responsible for the emergence of CM in all these cases and whether or not the mechanism leading to CM must involve explicit or implicit aligning interactions.
In this paper, we introduce a mechanism for CM that is based on a novel paradigm: the emergence and growth of regions of coherent motion due to standard elasticity processes. We explore this mechanism by introducing a simple two-dimensional active elastic sheet (AES) model where the individual agent motion is determined by attraction-repulsion forces only and no orientation information is exchanged between agents. Rather than considering an active fluid where particles can flow with respect to each other, we consider here a two-dimensional active solid or active crystal by describing a configuration of self-propelled agents that act as an elastic membrane, where interacting neighbours remain coupled by linear elastic forces throughout the dynamics, regardless of the amount of strain in the system. The paper is organized as follows. In section 2, we introduce our AES model. In section 3, we characterize its typical dynamics and stationary solutions, focusing on the self-organization process that leads to CM and on the order-disorder transition observed at the critical noise. In section 4, we describe the energy cascading mechanism responsible for CM in the AES model and characterize its convergence dynamics. Section 5 presents an analytical linear stability calculation that provides a necessary condition for algorithms to sustain CM. Section 6 explores three examples of the AES dynamics: an elastic ring-shaped configuration with oscillating radius, the propagation of an angular perturbation on a group of aligned agents, and a variation of the AES model where each agent has a different self-propulsion speed. Finally, in section 7 we discuss potential applications and present our conclusions.

Active elastic sheet (AES) model
We consider a system of N agents moving on a two-dimensional plane. The position x i and orientation θ i of each agent i are evolved according to the following set of overdamped dynamical equations: Here, v 0 is the forward biasing speed that induces self-propulsion (injecting energy at the individual particle level) and parameters α and β are inverse translational and rotational damping coefficients, respectively. Unit vectorn i points parallel to the heading direction of agent i and unit vectorn ⊥ i points perpendicular to it. The total force over agent i is given by the sum of linear spring-like forces with equilibrium distances l i j and spring constants k/l i j . We chose to define the stiffness of each spring as inversely proportional to its natural length l i j in order to mimic the elastic response of equivalent physical springs of different lengths. Each set S i contains all agents that interact with agent i. All S i sets are chosen at t = 0 from the local neighbourhood of agent i and remain constant throughout the integration. We note that the AES model is similar to a spring-mass model of an elastic sheet [25], with inert masses replaced by self-propelled agents that can only turn or move forward or backwards. For zero noise (D r = D θ = 0), each agent simply turns at a rate proportional to the projection of elastic forces perpendicular to its heading direction and moves forward or backwards driven by the projection of these forces parallel to its heading direction and by the self-propulsion term v 0 . We introduce sensing noise (errors in the measured forces) by adding D rξr to F i , where D r is the noise strength coefficient andξ r is a randomly oriented unit vector. We introduce actuation noise (fluctuations in the agent dynamics) by adding D θ ξ θ to the heading direction of each agent, where D θ is the noise strength coefficient and ξ θ a random variable with standard, zero-centred normal distribution of variance 1. We also tested cases where ξ θ followed a homogeneous distribution in the [−π, π] interval, finding equivalent results (data not shown).
We designed the AES model to explore CM under conditions that are very different from those in the Vicsek model. Firstly, in the Vicsek case agents only exchange information on their relative heading angle, while in the AES model they only know the relative positions of their neighbours. We define here as explicit alignment interactions those that depend directly on the relative angle between agents. We define as implicit alignment any other interaction that does not use information on relative heading angles but still tends to align the agents indirectly. Note that the elastic forces that produce alignment in the AES model are fundamentally different from ferromagnetic-like aligning interactions. Indeed, they are caused by elastic deformations that result from integrating misaligned trajectories over time. They will thus continue to act as long as the lattice remains deformed, regardless of the relative angle between agents. Therefore, no direct algebraic relationship can be established between the relative angle and the turning dynamics. Secondly, the Vicsek model displays internal flows that allow agents to interact with different neighbours over time. In fact, long range order cannot be achieved by the Vicsek dynamics if agents always interact with the same neighbours [13,15]. This result is established for equilibrium systems by the Mermin-Wagner theorem [26] and has been shown to extend to this non-equilibrium case: even relatively small systems do not display CM if interactions are only local and neighbours do not change over time [27]. In contrast, we show below that the AES model achieves CM despite having agents that interact with a fixed set of neighbours. Finally, while both models describe overdamped systems, in the AES model the angular equation of motion (2) gradually changes the heading angle instead of instantaneously switching it to the next desired heading direction [12,23,24]. We show below that this, apparently small, difference is essential for achieving CM in the AES model.

Numerical dynamics and stationary solutions
In this section we present numerical simulations that characterize the typical dynamics and stationary solutions of the AES model.
We integrated equations (1) and (2) numerically using a standard Euler method, which yields the expressions where t is the numerical time-step. Note that D r and D θ are divided here by √ t in order to properly take account of the accumulation of noise over time [28]. The degree of alignment in the system is monitored by computing the usual polarization order parameter, defined as If all agents are perfectly aligned, we have ψ = 1; if they are instead randomly oriented or rotating about the group's barycentre, we have ψ = 0. All simulations in this paper were carried out using: α = 0.01, β = 0.12, v 0 = 0.002 and t = 0.1, unless otherwise noted. We also tested other parameter combinations without finding any significant qualitative difference in the resulting dynamics. Four examples of these explorations are presented in figure 2 below.  In row (C), each agent is coloured according to the degree of local alignment of its local neighbourhood, with ψ loc = 1 indicating full alignment and ψ loc = 0, no alignment. Simulation videos can be found in the supporting information, videos 1-3, available at stacks.iop.org/NJP/15/095011/mmedia. D θ and D r . As time advances, individual self-propulsion deforms the lattice, producing elastic forces that in turn affect the agent dynamics. Growing regions of coherent motion develop, eventually deforming the whole structure (A2 and A3) until the group starts translating or rotating collectively. The displayed case converges to a rotating state (A4) with an axis of rotation that does not coincide with its barycentre. The group therefore translates while rotating. In general, rotating solutions are often observed in our AES simulations. Note, however, that these states will always have higher elastic energy than translating ones, because AES structures must rotate like a solid body, where agent speeds grow linearly with the distance to the centre of rotation. If all agents have the same preferred speed v 0 , those in the inner and outer shells must be slowed down or sped up by elastic forces. Consequently, rotating states like the one on panel A4 have higher stored potential energies and are metastable. They are less frequently observed and, if integrated long enough and with high enough noise, they will eventually transition to a , 4000 (D3) and 6000 (D4). All simulations were started from the same initial condition (with random orientations and all virtual springs at their natural length), shown on panel A1. Each agent is coloured according to the degree of local alignment of its local neighbourhood, with ψ loc = 1 indicating full alignment and ψ loc = 0, no alignment. Simulation videos can be found in the supporting information, videos 4-7, available at stacks.iop.org/NJP/15/095011/mmedia. lower-energy, translating state. We chose to display a rotating state on panel A4, however, to illustrate its rich dynamics, which cannot be attained by the Vicsek model. We point out that rotating solutions of the AES system are very different from those in [29][30][31], where interacting agents can change and do not have to rotate collectively like a solid body. In those systems, no simple argument can show which state has higher potential energy and it cannot be determined, a priori, which is stable or metastable.
Row B presents the dynamics of an elastic rod comprised of N = 118 agents arranged into three rows. At t = 0 (B1), randomly oriented agents are positioned with distances to nearest neighbours d B = 0.32 (within each row) and d * B = 0.58 (between rows). All agents separated by a distance d < 1 are then linked by springs of natural length d and spring constant k/l = 1/d. Noise is the same as in row A. Here again, after the initial transient period, larger and larger regions of coherent deformation emerge (B2 and B3) until CM is reached and the rod starts moving (B4). Note that the first bending mode has here the largest final deformation, thus favouring a collective heading that is perpendicular to the rod's axis. This observation opens the interesting possibility of controlling the direction of self-organized CM by changing the shape of the agent formation.
Row C displays N = 891 agents forming an arbitrary square-like structure with two holes. We refer to this case as an active solid due to the lack of regularity in agent positions. Agents are initially distributed homogeneously within the predetermined shape of the structure, with random positions and orientations. All agents separated by d < 1 are then linked by springs of natural length d and spring constants k/l = 5/d. Noise is set here to zero, but we verified that the same qualitative dynamics is observed for small enough D r and D θ values. To highlight the ordered regions, we define the measure of local order where C(S i ) denotes the cardinality of set S i . If agent i and all the neighbours it interacts with are aligned, ψ loc i = 1; if they are isotropically oriented, ψ loc i = 0. We colour each agent in the elastic solid according to its ψ loc i value, following the scale displayed on the figure. At t = 0, agents are randomly oriented and ψ loc i values are typically small (C1). As time advances and the elastic sheet deforms, aligned regions of coherent motion (with ψ loc i close to 1) appear and grow (C2 and C3). Finally, the whole structure starts moving collectively when agents become sufficiently aligned (C3).
We carried out additional simulations to explore the effects of parameter changes on the dynamics described above. In order to reduce the number of parameters, we first write equations (1) and (2) in non-dimensional form by expressing them in terms of the characteristic units of length and time, l and l/v 0 , respectively (considering here cases where all springs have the same natural length l). We find with Here With this rescaling, the deterministic part of the dynamics depends only on the two non-dimensional parameters: A = αk/v 0 and B = βkl/v 0 . These control the extent to which elastic forces will result in translation and rotation, respectively. In figure 1  The AES model displays an order-disorder phase transition as a function of noise similar to that in the Vicsek model. Figure 3 examines this transition in the same hexagonal active crystal displayed on column A of figure 1 and in a larger (N = 547) hexagonal configuration with identical parameters. The leftmost column presents results as a function of sensing noise D r and the central column, as a function of actuation noise D θ . Top panels display the mean and local maxima of the distribution of ψ values computed for the last 10 6 time-steps (after discarding the initial 10 6 steps to ensure convergence to a statistical steady state) in 30 (or 80, in the transition region) equivalent runs per noise value. Bottom panels display the corresponding values of the Binder cumulant G = 1 − ψ 4 /3 ψ 2 2 [24]. As the level of either type of noise is increased, the system undergoes a discontinuous transition from an ordered state where agents self-organize, to a disordered state where they continue to point in random directions without achieving CM. This is evidenced by the discontinuous drop of the order parameter at the critical point and by the Binder cumulant, which is known to become negative in the transition region for first order transitions with bistable solutions. In the sensing noise case, we find that G reaches negative values for N = 547, indicating that the transition is discontinuous for large enough systems. In the actuation noise case, N = 547 does not appear to be large enough to reach G < 0, but the dip at the transition region drops further and further below G = 1/3 (the expected value in the disordered phase) as the system size is increased, which strongly suggests a discontinuous transition [24,32,33]. We confirm the presence of a bistable region (where ordered and disordered solutions coexist) for both cases by displaying on the rightmost column the distribution of ψ values in the transition region. It is apparent that these distributions become more bimodal as the system size is increased. Finally, we point out that we observed an equivalent discontinuous transition when using either Gaussian or homogeneous noise distributions and for other spatial configurations (data not shown).

Convergence dynamics
We focus in this section on the convergence dynamics of the AES model. First, we show that the mechanism that leads to CM can be best understood by decomposing the energy of the system into its elastic modes. Then, we study how the convergence to the ordered state depends on the system size.

Energy cascading mechanism
We begin by computing the spectral decomposition of the energy into the elastic modes of the system. In order to do this numerically, we first determine the elasticity matrix K of the structure After t ≈ 1500, most elastic energy has been converted to kinetic energy, where it either dissipates or flows to lower modes, eventually reaching the zero (translational) mode. by perturbing the agent positions one by one with respect to an equilibrium configuration where all springs are at their natural length. We then find the eigenvalues and eigenvectors of the K matrix. As in standard elasticity, these eigenvectors define an orthogonal base over which we can decompose the displacements and velocities. The corresponding amplitudes yield the spectral decomposition of the potential and kinetic energy, respectively.
We display on figure 4 the dynamics of the total kinetic and potential energy (panels (A) and (B)), and of the spectral decomposition of the potential energy (panel (C)) for the same system simulated in row (A) of figure 1, but with zero noise (D r = D θ = 0) and for a run that converges to a translating solution. Note that the sum of potential plus kinetic energy is not conserved here due to the overdamped nature of the dynamics. For this system, we have 182 elastic modes, corresponding to 91 agents with two positional degrees of freedom per agent. The amplitude of these modes are displayed on figure 4(C) as a function of time, numbered in order of growing energy and smaller scales, without accounting for degeneracies. The initial condition is set with all agents randomly oriented and placed in an undeformed hexagonal lattice. At t = 0, the potential energy is therefore zero while the kinetic energy is equal to E * K = (1/2)N v 2 0 = 1.82 × 10 −4 , where we have set the agent mass to 1. At the beginning of the dynamics, kinetic energy drops and potential energy grows, due to elastic forces. Given the disorder in the system, this potential energy is initially broadly distributed throughout the different energy modes. As time advances, the system rearranges itself into states with decaying elastic energy and growing kinetic energy, until the former reaches values close to zero while the latter reaches again its E * K value. Figure 4(C) allows us to visualize the mechanism that leads to self-organization. After the initial transient, all modes are excited and their amplitudes oscillate while dampening, with higher modes decaying faster than lower ones. This results from a combination of a standard elasticity process and the coupling between elastic forces and heading angle imposed by the model. Indeed, it is well-known that higher energy modes dampen at a faster rate in elastic systems, since they are more rigid and have higher oscillation frequencies, which typically leads to faster dissipation. In this active system, however, each agent is continuously injecting energy at the individual level through its self-propulsion term, so motion cannot dampen out. Instead, elastic forces will steer agents away from higher modes more strongly than from lower modes. If while doing this agents do not re-excite higher elastic modes faster than these decay, the self-propulsion energy will be channeled to lower and lower modes until the first (rotational) mode or zero (translational) mode is reached and CM is achieved. If instead agents feed too much self-propulsion energy to higher modes while turning, these modes will always remain excited and no CM state will be reached. We conclude that not every active elastic system will achieve CM. For example, we consider in section 5 a constantspeed algorithm that does not display CM despite having the same angular dynamics as the AES model. Finally, we note that by writing equations for the flow of energy between the different modes, we can study analytically the conditions required for self-organization beyond the linear stability analysis described below. These calculations are beyond the scope of the current paper, however, and are left for future work.

Dependence on system size
We study here how the convergence dynamics depends on the system size. Figure 5 displays the global order parameter ψ as a function of time for hexagonal active crystals with N = 91, 547, 1027 and 5167 agents. Each simulation is started from a random initial condition, using the same parameters as in figure 4. Ten convergence curves are presented per system size. We observe that for N = 91 not all runs converge to the aligned (ψ ≈ 1) state. Instead, four runs reach the metastable rotating (ψ ≈ 0) state and remain trapped there until the end of our simulation time, which was set here at t = 10 4 (much longer than the displayed time frame). For the other sizes, however, no run ultimately converges to the rotating state. This is because the larger a rotating structure is, the faster outer agents must advance in order to maintain cohesion. For large enough systems, the drag introduced by these agents will be enough to destabilize the rotating solution. Figure 5 also shows that convergence times have a large variability. Despite this, the figure readily provides a rough estimate of how these times scale with the system size. Indeed, the time frame displayed on each panel is proportional to N 1/2 . Given that even with this rescaled temporal axis curves seem closer to the ordinate axis for higher N , it is apparent that the typical convergence time grows here slower than N 1/2 . The fact that the convergence dynamics is strongly non-monotonous suggests that its variability is a reflection of the complex dynamical landscape that the system navigates, where it can spend unpredictable amounts of time near local attractors and metastable states. Since the main metastable state is the rotating solution, Figure 5. Convergence dynamics of the order parameter ψ for four hexagonal active crystals of different sizes. All parameters are the same as in figure 4. Ten runs with different, randomly oriented initial conditions are displayed for each size. When ψ approaches 1, the system is converging to an aligned translating state and when it approaches 0, to a rotating (metastable) state. All system sizes display a broad variability of convergence times, with larger systems typically taking longer to converge. a systematic scaling analysis of convergence times will require not only a much larger set of simulation runs, but also structures that suppress the rotating state, such as strongly elongated shapes.

Linear stability analysis
One of the interesting aspects of our AES model is that we can use a continuous elastic sheet approximation to carry out analytical calculations. We follow this approach and perform a standard linear stability analysis [25] of the zero noise case to investigate which specific dynamical rules can sustain translating CM solutions. We begin by writing the elastic forces F = (F x , F y ) that result from small displacements u = (u x , u y ) of points on the membrane with respect to their equilibrium position. These are given by the standard elasticity equations where the elastic constants are the Lamé parameter λ and shear modulus µ [25]. By linearizing the equations of motion (1)-(2) around an equilibrium solution with undeformed membrane and all agents moving at speed v 0 in thex direction, we find the following expressions for u x , u y , and the perturbation field φ of the heading angle: Replacing equations (13) and (14) into (15) and casting the resulting expression in Fourier space with wavevector components k x and k y , we can express the perturbation dynamics in matrix form and compute its eigenvalues to determine the linear stability of the system. We find that satisfies the characteristic equation 3 Using Routh's stability criterion, here given by C 1 C 2 > C 0 [34], we find that the system will be linearly stable if αβv 0 (λ + µ) 2 k 2 x k 2 y > 0, which is always verified. We conclude that the translating CM solution is linearly stable for all parameter values.
Analytical calculations like those presented above allow us to determine, a priori, which elasticity-based equations of motion will be able to sustain CM. We found that most variations of the AES model cannot support stable aligned solutions. Consider, for example, an algorithm where the heading angle is determined by equation (2) but the speed is fixed to v 0 , by setting α = 0. Given that the angular dynamics remains unchanged, one could naïvely think that this system will align like the AES model. We will now show, however, that this is not the case. The characteristic equation for α = 0 becomes 3 which has solutions = 0 and 2 = −v 0 β[µk 2 x + k 2 y (λ + 2µ)]. Since ∈ i , linear perturbations will not dampen out, but produce instead permanent oscillations. We confirmed through numerical simulations that, even for zero noise and starting from an aligned initial condition, the group will lose order as agents rotate in place. After testing several elasticitybased active systems, we found only one other example that can sustain CM: a variation of the model introduced in [19] that we will further discuss in section 7.

Exploring AES dynamics
We carry out in this section an initial exploration of three different dynamical setups of the AES model that allow us to better understand its typical behaviour and modelling possibilities. We first consider a rotating ring solution, then the propagation of perturbations on an aligned rodlike configuration, and finally a variation of the AES model where each agent has a different self-propulsion speed v 0 .  (clockwise). The system is then integrated forward in time with zero noise. Top panels display snapshots at t = 5135, 7655, 10 5 and 4 × 10 5 ; the bottom panel shows the mean radius of the ring as a function of time. Initially, the ring expands, increasing elastic forces until a maximum radius is reached (panel 1) and the ring starts contracting. The contraction speed then increases until it reaches a maximum (panel 2) and agents start turning outwards until they move again tangentially, reaching a minimum radius. This breathing mode continues to oscillate with decaying amplitude until it fully dampens out. The circular configuration then loses stability, exciting higher modes that deform it (panel 3). This state survives with different levels of deformation until the end of our integration time (panel 4). Figure 7 shows the propagation dynamics of a local perturbation of the heading angle on an aligned CM state. We set up a three-row rod-like structure similar to that in figure 1(B), with Agents are initially placed on a three-row rod-like structure similar to that in figure 1(B) but longer (containing N = 499 agents), with all agents heading in the same direction (θ = 0) and aligned to the rods main axis. Simulations are carried out with the same parameters used in figure 1(B), but for zero noise. At t = 0, the frontmost agent is rotated by π/18 (left column) or π/2 (right column). The plots display the heading angle of all 167 agents on the central row, numbered from back to front. The large perturbation propagates faster and produces a wake of longer wavelength than the small one. the same parameters but in a longer configuration with N = 499 agents. All agents are initially placed aligned and pointing on the same direction as the rod axis, which we define as θ = 0. At t = 0, we perturb the orientation of the frontmost agent, rotating it by π/18 (small perturbation) or π/2 (large perturbation). The system is then integrated forward in time with zero noise. We plot the heading angles of all agents on the central row of the rod-like structure at four different moments in time, indexed in order of their position from back to front. Both small and large perturbations display here a wake of persistent angular oscillations behind them. Note, however, that other preliminary simulations that we have carried out on less elongated structures show this wake rapidly decaying after the passage of the initial perturbation. For the long rod-like case presented here, we observe that small angular perturbations propagate faster and leave a wake of shorter wavelength and smaller amplitude than large ones. These results illustrate the rich dynamics exhibited by the propagation of perturbations in the AES model. Their study will require further systematic analyses that are left for future work.

Heterogeneous self-propulsion speeds
We now consider a variation of the AES model where all agents have different preferred speeds. Figure 8 displays the dynamics of an active elastic hexagon with the same parameters Snapshots of a hexagonal AES simulation with heterogeneous preferred speeds at t = 0 (1), 7300 (2), 12 000 (3) and 22 000 (4). Each agent's self-propulsion speed v 0 is chosen at random between 0 and 0.04. After a very long integration time, the system reaches a quasi-ordered state of CM where the orientation of faster agents oscillates broadly so they can remain cohesive with their slower neighbours. Agents are coloured based on their degree of local alignment ψ loc , as in figure 1(C). The simulation video can be found in the supporting information, video 9, available at stacks.iop.org/NJP/15/095011/mmedia. as in figure 1(A), but where instead of fixing all self-propulsion speeds to v 0 = 0.02 we select a different v 0 for each particle, at random, from the interval 0 v 0 0.04. Agents are coloured based on their local alignment, as in figure 1(C). Interestingly, even for this highly heterogeneous system agents manage to self-organize and achieve CM, albeit after a very long relaxation time. As time advances, growing regions of coherent motion emerge (panel 2), but a localized part of the hexagonal structure remains persistently disordered (the lower-left quadrant of the hexagon on panel 3). This is the area that displays the highest differences of preferred speeds, which makes CM harder to reach. Eventually, the system finds a way to fully self-organize and the groups starts moving collectively. In the case shown, it advances in a curved trajectory due to the random accumulation of faster agents at one side of the structure. A salient feature of the final quasi-ordered state is that it can never reach the stationary solution where all agents are fully aligned. This is because the orientation of faster agents must oscillate strongly so they can remain cohesive with their slower neighbours. The ability displayed here by the AES model to self-organize even in highly heterogeneous situations opens the possibility of constructing a rich variety of active solids that produce diverse collective dynamics by assembling groups of agents with different individual characteristics.

Discussion and conclusions
We have identified in this paper an alternative, elasticity-based mechanism for achieving CM and introduced the AES model to illustrate it. Up to now, the only existing theoretical framework for explaining how a variety of systems self-organize to achieve CM was based on the Vicsek model [12] and on the active hydrodynamic theory first introduced by Toner and Tu [15]. Our work develops a very different theoretical framework, providing a simple alternative mechanism for CM that is based on elasticity instead of alignment consensus or momentum transfer and requires no exchange of orientation information. We found only one other system that can display CM under similar conditions: the model introduced in [19] to study the collective migration of tissue cells. In a version of this model designed to study active jamming, where agents are confined in a circular box and only have repulsive interactions, a similar elasticitybased mechanism was shown to be responsible for the dynamics of the jammed phase [20].
The relevance of aligning interactions for achieving CM has been a long-standing issue in the field. Only a few studies have considered systems where agents align and move collectively without exchanging orientation information, and their underlying mechanism remained unexplained. In fact, seminal work by Grégoire and co-workers [23,24] suggested that CM was not possible in minimal models without aligning interactions. Our studies show that their model did not to converge to CM when aligning interactions were turned off because of two reasons. Firstly, it considered agents with constant speed. As shown above, this is enough to prevent the elasticity-based mechanism from achieving CM in the AES case. Secondly, their agents switch to the next desired heading direction in one time-step, instead of integrating equation (2). This instantaneous relaxation can stop energy from smoothly flowing to lower energy modes and from developing growing regions of coherent motion.
To the best of our knowledge, this work is the first to combine the theories of elasticity and CM. It could find applications in any system of self-propelled agents or active particles that hold an approximately rigid formation. The simple idea that the slower dissipation of lower, more coherent elastic modes can lead to self-organization in active matter provides a unifying framework for studying disparate systems such as cell migration in tissue development (where elasticity has been shown to play a relevant role) [19,21,35], flocking models or swarm robotic control algorithms without explicit alignment [36], or even the recently developed active gels, which have no self-propulsion but could also display elasticity-driven self-organization [37].
An important application of this work could be found in experimental studies. There is a growing number of experiments that study CM in natural and artificial systems. One of their current objectives is to verify if the properties of their propagating waves coincide with the predictions of the hydrodynamic theory of CM [15,16]. Our work suggests that these efforts may be incomplete, since they only search for alignment-based waves, whereas waves with different properties would result from active elastic interactions. This suggests the interesting possibility of determining the type of interactions present in a given system based only on collective properties such as the propagation of perturbations. We expect most natural groups to move collectively as a result of a combination of both types of mechanisms, alignment-based and elasticity-based, effectively behaving as an active viscoelastic material of aligning components. The relevance and role of each mechanism could change at different time-scales and in different regimes. However, given that some type of cohesion is required to form a cohesive group, the elasticity-driven mechanism should play a role in many cases. Note that the AES model does not intend to describe the details of specific experimental systems, but to illustrate instead the possible relevance of the elasticity-based mechanism in their convergence to CM. In order to further quantify the dynamics of a specific system, a more detailed model should be developed (potentially combining different aligning mechanisms), with its corresponding parameters and characteristic scales obtained from direct measurements.
Finally, an appealing aspect of the AES approach is that it is well-suited for analytical studies. In addition to the stability calculations described in this paper, we can envision carrying out a stochastic differential equation analysis to determine the critical noise level required for the transition. We could also develop energy-cascading arguments that describe the self-organization dynamics as requiring a net energy transfer from higher to lower elastic modes. These analytical approaches could help improve our fundamental understanding of CM and of the more general class of non-equilibrium, self-organizing systems where energy is injected at the smallest scales, which are highly relevant in the study of biological systems and new active materials.