Order–disorder transitions in a minimal model of active elasticity

We introduce a new minimal model for self-propelled agents that attract, repel, and align to their neighbors through elastic interactions. This model has a simple mechanical realization and provides an approximate description of real-world systems ranging from active cell membranes to robotic or animal groups with predictive capabilities. The agents are connected to their neighbors by linear springs attached at a distance R in front of their centers of rotation. For small R, the elastic interactions mainly produce attraction-repulsion forces between agents; for large R, they mainly produce alignment. We show that the agents self-organize into collective motion through an order–disorder noise-induced transition that is discontinuous for small R and continuous for large R in finite-size systems. In large-scale systems, only the discontinuous transition will survive, as long-range order decays for intermediate noise values. This is consistent with previous results where collective motion is driven either by attraction–repulsion or by alignment forces. For large R values and different parameter settings, the system displays a novel transition to a state of quenched disorder. In this regime, lines of opposing forces are formed that separate domains with different orientations and are stabilized by noise, producing locally ordered yet globally disordered quenched states. Video Abstract: Order–disorder transitions in a minimal model of active elasticity

Several minimal models have been introduced to help understand the self-organizing principles and dynamics of collective motion [18]. These models are typically defined by a set of self-propelled agents that interact with their neighbors. Many of them include explicit alignment interactions through which each agent tends to align to the heading direction of its neighboring agents. We will refer to these as velocity-based models, since they require exchanging information on the relative velocities between agents (or, at least, their relative heading directions). One of the first and most influential minimal models of collective motion, the Vicsek model [19], is in fact purely velocity-based. We expect, however, that real-world systems will generally have interactions that involve both the relative velocities and the relative positions of neighbors, as it has been shown, for example, in fish schooling experiments [6,20,21]. Attraction-repulsion forces based on relative positions have thus been added to multiple models, although contrast to the AE case, however, the AEA model can be fully mapped to a mechanical system, as explained below.

Model description
Let us imagine each agent as a self-propelled cart that can only turn about its center of rotation and move forward or backward along its heading direction (no sideways displacements are allowed), as sketched in figure 1. This restricted motion can be achieved by imposing an additional constraint that stops the cart from moving sideways (i.e. a force that restricts motion along n ⊥ i ). This would correspond, for example, to the dynamics of a cart mounted on two no-slip side-by-side wheels (represented by the dashed lines in the figure). In such a system, the no-slip constraint allows no sideways motion, the linear speed of the center point between the wheels is (s l + s r )/2, and the angular speed about it is 2(s r − s l )/D (where s l and s r are the tangential speeds of the left and right wheels, respectively, and D is the distance between them). We introduce self-propulsion forces such that, if no external forces are applied, the cart advances in a straight line at a fixed preferred speed v 0 , alongn i . We introduce interactions by connecting each cart to its neighbors through linear springs attached to an arm, a distance R in front of its center of rotation. As the springs stretch or contract, they apply elastic forces to the arm, affecting the speed and heading direction of the cart according to a standard decomposition into radial and tangential forces. The total force F over agent i is thus decomposed into a displacement force F ·n i and torque Rn i × F. Figure 1 presents a diagram of the mechanical system described above for a pair of interacting agents, i and j, represented by the gray discs with positions x i and x j (also indicating their centers of rotation), and heading directionsn i andn j , respectively. If we consider an overdamped system of N of agents moving on a two-dimensional arena, each agent i will obey the following dynamical equationṡ Here, x i represents the position of the center of rotation of each agent i, the unit vectorn i points in its heading direction, with θ i being its angle, and the unit vectorn ⊥ i points perpendicular to it. The parameter v 0 controls the preferred self-propulsion speed. The term D θ ξ θ introduces fluctuations in the heading direction, where D θ controls the noise strength and ξ θ is a delta-correlated random variable with standard, zero-centered normal probability distribution. Note that, in contrast to the original AE model, here we only add noise to the angles and not to the positions. It was shown in [29,30], however, that both types of noise produce very similar results.
The sum of elastic forces over the arm applied by all springs connected to neighboring agents is represented by F i . The F i ·n i term in equation (1) thus gives the projection of these forces along the direction of motion, which will affect translation, whereas the F i ·n ⊥ i term in equation (2) is the perpendicular projection, which will affect rotation. Parameters α and β control the coupling between these forces and the resulting translation and rotation, respectively.
In the physical representation of this overdamped system, α is determined by the friction associated to the cart translation along the heading directionn i . The value of β depends in turn on the friction associated to the cart rotation and will grow linearly with the length R of the arm, as it affects the torque. In what follows, however, we will set the value of β independently of R, in order to be able to consider the effect of R = 0 while keeping β > 0, and thus recover the original AE model.
As in the AE model, F i is given by the sum of all elastic forces on agent i produced by the linear spring connections to its neighbors, that is, Here, l 0 is the natural length and k/l 0 the elasticity constant of all springs. The set S i contains the indexes of the neighbors that interact with agent i. In all simulations below, each agent i interacts with its six nearest neighbors, which remain constant throughout the simulation. The system thus forms an elastic hexagonal lattice. Up to this point, the AE and AEA models are defined by identical equations, since we have not yet specified how r ij is computed. In the AE case, r ij = d ij = x j − x i because the purpose was to define a model where interactions depended exclusively on the relative positions of neighbors, and not on their relative headings. This implied that the AE model could not be mapped to a mechanical system, however, since the decomposition into translational and rotational forces would still require the springs to be connected to an arm. Instead, in the AEA model, r ij is defined as the vector that connects the two points where the spring linking agents i and j would be attached in a real-world mechanical system (see figure 1), and is thus given by For R = 0, we recover the AE model, where the interactions depends only on relative positions and not on relative headings. For R > 0, the larger the value of R, the more the interactions depends on relative headings. Finally, in the simulations below we will measure the degree of alignment order using the standard polarization order parameter With this definition, p = 1 if all agents are perfectly aligned and p = 0 if they are heading in different directions. Note that, as we will see in the following sections, the latter can correspond to cases with random headings or with internally aligned domains that point in different directions.

Motivation and applications
In addition to its mechanical interpretation, an important motivation for the AEA model is that it allows us to study interactions that are mostly dependent either on the relative positions of neighboring agents or on their relative headings. Indeed, for R = 0 we recover the original AE model and the interactions depend only on relative positions. For larger R values, changes in the relative headings will directly change r ij , thus introducing alignment forces. The nondimensional ratio R/l 0 can therefore be viewed as a parameter that interpolates between purely position-based and mostly alignment-based interactions, when taking values between 0 and 1, respectively. For R/l 0 > 1, we note that the collective dynamics can display additional complications, since the point where the forces are applied on the arm of each agent would be beyond the position of its nearest neighbors, which can result in a various frustrated or 'tangled' states that have no direct physical interpretation. In the following sections, we will therefore focus on analyzing systems with 0 R l 0 .
Further motivation for the AEA model is given by its potential application as a minimal approximate description of other active systems. It could capture, for example, essential dynamical features of any group of self-propelled components with short range, contact-mediated interactions with a fix set of nearest neighbors, ranging from migrating cells to self-propelled rods [36,37]. (Note that the fixed neighbors condition can be satisfied even without attraction forces if the components remain densely packed throughout the dynamics.) Although in these systems the contacts between nearest neighbors typically produce a combination of central forces and torques that are more complex than those in the AEA model, they share similar properties and could produce equivalent collective states. Finally, the AEA model could also capture the dynamics of agents that react to attraction-repulsion social forces over their estimated future position. Indeed, since the forces are applied at a distance R in front of each agent, they can be interpreted as interactions over a linear extrapolation of its future position after a time τ , when advancing at constant velocity (R/τ )n. In regimes where they are moving collectively at the same approximate speed, the AEA model could thus provide a simple description of the dynamics of agents with predictive capabilities, such as animal groups or robot swarms. Snapshots of an AEA model simulation showing N = 100 self-propelled agents in their self-organizing process. Each black dot represents the center of rotation of an agent, the blue line is its forward arm (with respect to its heading direction), and the dotted gray lines are the linear springs that mediate elastic interactions. At t = 0 (a) the system is initiated with random headings. At t = 3000 (b) several regions display alignment. At t = 6000 (c) most agents are heading in the same direction, as the system achieves collective motion. (Simulation videos in the supplementary material (https://stacks.iop.org/NJP/23/023019/ mmedia)).

Simulation setup
The initial conditions for all simulations were set up by placing N agents, each connected to its six nearest neighbors, forming a perfect hexagonal lattice in a periodic two-dimensional arena. Their initial orientations were either random or perfectly aligned, as specified below. All simulation arenas had an aspect ratio of √ 3/2, in order to contain a homogeneous hexagonal distribution of √ N × √ N agents. In the following sections, we will investigate the dynamics and order-disorder transition of the AEA model as a function of the noise level D θ for R > 0 and different parameter values. In order to reduce the parameter space and to compare with previous simulations of the original AE model, we will fix the self-propulsion speed to v 0 = 0.002, the spring constants to k = 5, and their natural lengths to l 0 = 1 in all the runs presented below. In sections 3-5 we will explore the detailed AEA dynamics by simulating a relatively small system of N = 10 × 10 = 100 agents, and in section 6 we will extend these analyses to larger systems.

Self-organizing dynamics
We begin our analysis by describing the self-organizing dynamics of the AEA model for α = 0.01, β = 0.12, an intermediate arm length R = 0.3, and noise level D θ = 0.2. Figure 2 displays snapshots of an AEA simulation with N = 100 agents linked by linear springs in a 10 × 5 √ 3 periodic box where all can be at their equilibrium distances when forming a perfect hexagonal lattice. Each snapshot displays the positions of the centers of rotation of the agents (black points), their forward arms (blue lines), and the springs that mediate their elastic interactions (gray dotted lines). Panel (a) presents the initial condition at time t = 0, in which we place all agents with random orientations, filling the plane homogeneously (their centers of rotation forming a prefect hexagonal lattice). Panel (b) displays the system at t = 3000, after the self-propulsion and elastic forces start aligning the system. Finally, panel (c) shows the self-organized statistical steady state at t = 6000, after achieving collective motion, in which most agents head in the same approximate direction and the polarization p value is close to one.
We found in our numerical explorations that the AEA simulations appear to always rapidly self-organize (from initial random headings to the aligned state of collective motion) for any R 1 and D θ below a critical noise value, as long as the structure of the hexagonal lattice is not compromised (i.e. if each agent remains within the hexagonal cell defined by the springs connecting its immediate neighbors). For high values of v 0 , the dynamics becomes more complicated, as agents can overcome the elastic forces and the hexagonal lattice structure can become 'crumpled' (with some agents crossing the boundaries of their local hexagonal cells). These crumpled regions tend to persist for a very long time, often until the end of our simulations, and in most cases the system remains disordered. In a few test cases where we performed extremely long simulations, however, the system always eventually reached the aligned state of collective motion, below a critical noise, even as the crumpled regions remained. This suggests that the self-organizing dynamics is robust even for high v 0 , but given its slow relaxation and additional complexity when the hexagonal lattice structure is lost, this regime requires further systematic analyses that will be left for future work.
In the small v 0 regime with 0 < R 1 that is the focus of our attention, the aligned state of collective motion of the AEA model does not display the persistent oscillations observed in the original AE model for all v 0 , even without noise. This type of persistent oscillations were identified analytically in [38] in a minimal two-agent system with zero noise that is reduced to only two degrees of freedom by imposing symmetric initial conditions.
In order to understand the dynamics of small perturbations in the AEA case, we carried out the same two-agent linear stability analysis about the aligned state performed in [38]. The detailed calculations are presented in the appendix and yield the following eigenvalues of the stability matrix For R = 0, we recover the result previously obtained for the original AE model, with purely imaginary λ ± eigenvalues and thus persistent oscillations. For R > 0 and v 0 < βkR 2 /2, both eigenvalues will be real and negative, and any oscillation will be overdamped. This is consistent with the lack of oscillations observed in the simulations displayed on figure 2, where v 0 = 0.002 and βkR 2 /2 = 0.027. For v 0 > βkR 2 /2, the eigenvalues have a negative real part and an imaginary part, so we expect to see some damped oscillations as the system reaches the aligned state. This transition from overdamped to damped oscillations is apparent when comparing two additional videos that we include in the supplementary material, where we set the preferred self-propulsion speed to v 0 = 0.08 or v 0 = 0.008 while keeping all other parameters unchanged.
We observe there that the v 0 = 0.08 case displays oscillations whereas the v 0 = 0.008 does not, even in the presence of noise.

Transitions to dynamic disorder
We will now study the order-disorder transition to a state of dynamic disorder as a function of noise, examining its continuous or discontinuous nature for different values of R. We will focus in this section on the transition to dynamic disorder, where agents are persistently rotating in random directions without showing any significant local alignment, which is the usual disordered state in collective motion transitions.
To avoid the newly identified quenched disorder state that will be discussed in section 5, we found through numerical explorations that we can use smaller α and larger β values than those in section 3. We will therefore set α = 0.001 and β = 0.3 for all simulations in this section, which moves the quenched disorder transition to higher R and D θ values. Figure 3 displays the polarization as a function of noise D θ for different values of R, ranging from 0 to 0.4. Each point is the result of 20 simulations of 10 6 timesteps in a system of N = 100 agents, half initialized with random orientations and the other half with all agents aligned. We selected the last 10 5 timesteps of every simulation and computed their instantaneous polarization p values per frame, using equation (5). Each point in the figure is a local maximum of the distribution of all p values collected for a given R and D θ . In cases where the system is bistable, the plot will thus show two polarization values at a given noise level. It is apparent in this figure that the nature of the transition changes as a function of R. It is discontinuous for R = 0, as previously reported for the original AE model, and remains discontinuous for a low R = 0.01 value. For R 0.02, however, the transition becomes continuous.  The results presented in this section are consistent with the hypothesis that the AEA model switches from a discontinuous to a continuous order-disorder transition because the system switches from a position-based to a velocity-based self-organizing mechanism and has a constant number of neighbors, as discussed in the introduction. Indeed, for small R values, the interactions will mostly depend on relative positions, as in [29,30], whereas for large R, they will be strongly affected by relative headings (aligned here with the relative velocities), as in [32]. We can therefore expect the transition to be discontinuous for small R and continuous for large R, which we confirmed above. Note that we do not expect this transition to survive in very large systems, however, since the interacting neighbors do not change over time. Indeed, as we also mentioned in the introduction, velocity-based models with a fixed interaction network of nearest neighbors cannot achieve long-range order for nonzero noise. To explore this, we will study in section 6 how the transition changes in systems of up to N = 40 000 agents.

Transition to quenched disorder
In addition to the usual transition to dynamic disorder described above, we identified a novel type of order-disorder transition in the AEA model. For larger values of R and α, the polarization drops sharply as D θ reaches a critical noise where the system falls into a state of quenched disorder, in which the positions and headings of the agents remain approximately fixed over time, in average, and are thus 'quenched' or 'frozen'. This is in contrast to the state of dynamic disorder presented in the previous section, which is often referred to as 'annealed' disorder. In this section, we will describe the properties and underlying dynamics of the transition to quenched disorder in the AEA model.
The solid curves on figure 5 present the polarization as a function of noise for combinations of α = 0.001 or 0.002, and R = 0.4 or 1.0. All other simulation parameters and procedures were kept unchanged from the previous section. We display here the mean polarization values because the distribution maxima plotted in figure 3 would not properly represent the converged quenched disordered states, since these can reach a slightly different low-polarization value in each run, as we will explain below. The curve with α = 0.001 and R = 0.4 corresponds to one of the standard cases of continuous transition to dynamics disorder included in figure 3. For all the other parameter combinations presented here, we find a discontinuous transition with an abrupt reduction of the polarization level occurring at a critical noise, corresponding to the transition to quenched disorder. Note that, in contrast to the discontinuous transition to dynamic disorder identified in section 4 for small R, the quenched disorder transition appears here for large R values.
In order to further analyze this phenomenon, we will focus on the interactions that have opposing self-propulsion forces, in which the self-propulsion of both agents is acting more strongly in opposite directions than in parallel, that is, where the difference between their heading angles is greater than 90 degrees. These forces will thus mainly contribute to compressing, extending, or rotating the spring between the two agents, depending on their relative positions, rather than to translating it. Indeed, for antiparallel heading directions, for example, the self-propulsion forces will tend to compress the spring if the agents are head-to-head, extend it if they are tail-to-tail, or rotate it if they are side-by-side. Instead, for parallel headings, the self-propulsion forces will tend to translate the spring. More specifically, we define as opposing self-propulsion forces all interactions between agents i and j that have c ij < 0, where The dashed lines with dots on figure 5 display the fraction of opposing self-propulsion forces for each state, computed as the total number of these forces over the total number of springs. In the α = 0.001 and R = 0.4 case, where the system transitions continuously to a state of dynamic disorder, this fraction increases smoothly as the polarization decreases. For all the other combinations of α and R in the figure, the discontinuous transition to quenched disorder coincides with a sharp increase in the fraction of opposing self-propulsion forces. This fraction then displays significant fluctuations above the critical noise, because each simulation stabilizes at a different, approximately constant number of opposing forces when the state of quenched disorder is reached.
The differences between the dynamic disorder and quenched disorder transition curves described above can be understood by examining figures 6 and 7, and the corresponding simulation videos in the supplementary material. These figures display the final states of 12 of the simulations used to generate figure 5 near the transition point. In each snapshot, we labeled with a red cross the midpoint of any spring subject to opposing self-propulsion forces, that is, with c ij < 0.
When the simulations are performed below the critical noise, they all reach an aligned state with almost no opposing self-propulsion forces, as shown in panels (a) and (d) of both figures.
As the noise is increased in the α = 0.001 and R = 0.4 case, the system transitions to a state of dynamic disorder where agents lose alignment. We observe increasing fluctuations in the orientations, which result in the presence of more points of opposing self-propulsion forces that are persistently appearing and disappearing at random locations homogeneously distributed throughout the simulation arena, as shown in figures 6(b) and (c), and in the corresponding supplementary material simulation videos.  For all other combinations of α and R in these figures, the system transitions instead to a quenched disorder state. At the critical noise level, we observe the formation of lines of opposing self-propulsion forces that percolate through the periodic box, as shown in figures 6(e) and 7(b), and (e). These lines then stabilize into persistent structures that become the boundaries between relatively fixed domains of locally aligned agents, each trying to advance in different directions but unable to do so because they are blocked by the others. As the noise level is increased, additional lines of opposing forces appear and their structures become more complex, as shown in figures 6(f) and 7(c), and (f). This explains the variations in the low-polarization values obtained for different runs above the critical noise, which we referred to in the context of figure 5, since the sizes and orientations of the domains that form in each specific run will determine the exact p values. Finally, in the supplementary material simulation videos, we note that the agent orientations, domains, and lines of opposing forces still display significant temporal fluctuations around a mean, due to the strong finite size effects that develop in these relatively small systems. Indeed, we observe much smaller fluctuations in the larger simulations presented in section 6.
We will now further describe the state of quenched disorder by considering the minimal two-agent symmetric dynamics detailed in the appendix. We begin by noting that the opposing self-propulsion forces in our simulations would correspond to the head-to-head or tail-to-tail unstable equilibrium solutions in the two-agent system. Although each point of opposing forces is unstable, when a line is formed at nonzero noise we can argue heuristically that the orientations of each pair of opposing agents is stabilized by its neighboring pairs. Indeed, for a line of opposing forces to disappear, agents at each side of the line would have to turn in groups, so that they can align with those at the other side without being stopped by the spring forces of their other neighboring agents. The presence of noise can frustrate this coordinated turning, however, thus stabilizing the lines of opposing forces. This heuristic explanation is consistent with the observation that the quenched disorder state appears when the first line of opposing forces is stabilized as it percolates across the periodic box.
The linear perturbation analysis of the two-agent head-to-head and tail-to-tail equilibrium states of opposing forces presented in the appendix produces the same stability matrix and characteristic eigenvalues, given by Since λ + is always positive, this shows that both states will always be unstable under linear perturbations, as expected. Although this minimal analysis cannot capture the collective noise-driven mechanism that stabilizes the lines of opposing forces, discussed above, it can show how this stability may depend on model parameters. Indeed, for larger α and smaller β or v 0 , the positive eigenvalue becomes smaller, so the level of instability decreases. This implies that the critical noise will also decrease, since the effect of the noise-driven stabilizing mechanism does not need to be so strong. We confirmed this result in figure 5, where the critical noise value decreases for larger α, and in other numerical explorations where we found that the critical noise increases with β and v 0 (data not shown). In addition, we note that the λ + and λ − expressions do not depend on R, which is also consistent with our numerical results since the critical noise is not strongly affected by the different values of R in figure 5(b).

Large system analysis
In previous sections, we analyzed the finite-size order-disorder transitions in the AEA model by focusing on relatively small systems. This allowed us to demonstrate that they can be clearly identified even for small N values, to examine their detailed agent and interaction dynamics, and to perform extensive parameter scans. In this section, we will study the properties of the transitions described above in larger systems. Figure 8 presents the three types of order-disorder transitions detailed in sections 4 and 5 for systems with N = 100, 1600, 10 000, and 40 000 agents, with equal mean density (in periodic boxes of size 10 × 5 √ 3, 40 × 20 √ 3, 100 × 50 √ 3, and 200 × 100 √ 3, respectively). The model parameters were set here to l 0 = 1, α = 0.001, β = 0.3, and v 0 = 0.002. Each point shows the mean polarization of the last 10 5 timesteps of at least four independent simulations lasting 2 × 10 6 timesteps.
The resulting curves clearly illustrate how the three types of transitions previously obtained in smaller systems change with system size. For a small R = 0.01 value, the transition is discontinuous. As in the original position-based AE model, we observe that it displays a clear asymptotic behavior akin to a first order phase transition in the thermodynamic limit [29]. Above the critical noise, we find the same state of dynamic disorder described in section 4.
For an intermediate R = 0.2, the transition appears as continuous for all system sizes, converging to a critical noise D θ ≈ 0.6. As shown in section 4 for smaller systems, this is equivalent to the situation in a velocity-based Vicsek model with topological or metric-free interactions [32]. In the N = 40 000 simulations, however, the ordered state loses polarization for D θ > 0.3 and we can extrapolate that it will disappear, with the continuous transition, for large enough systems. This result is consistent with the discussion in previous sections that links the continuous or discontinuous nature of the finite-size transition to the presence of velocity-or position-based self-organizing mechanisms, respectively. Indeed, as discussed in the introduction, the alignment forces in velocity-based models cannot achieve long-range order at nonzero noise levels if the interaction network is fixed in time (i.e. if the same neighbors interact throughout the simulation), as in our AEA model. Therefore, as R grows and the alignment forces become dominant, we observe a continuous transition that will disappear for large systems. A discontinuous transition will appear instead at D θ ≈ 0.3, since the elasticity-based self-organizing mechanism can still lead to an ordered state at lower noise values in large-scale systems, as shown for the AE model in [29].
Finally, for a large R = 1, the transition in large-scale systems is again discontinuous and leads to the state of quenched disorder, as in section 5. Although the critical noise converges to a well-defined value (here D θ ≈ 1.0), the ordered phase loses polarization in larger systems (here for D θ > 0.6), as in the R = 0.2 case. We thus expect this order-disorder transition to also disappear in large enough systems, although the state of quenched disorder will still be reached for high D θ values. Figure 9 displays the final state of some of the simulations used to generate figure 8 near the transition points. Videos of the same simulations are included in the supplementary material. For each curve, we selected two cases; one just above and one just below the critical noise value. In order to simplify these larger system snapshots, we do not show here the linear springs (all neighboring agents are connected, as in previous simulations) and each agent is simply represented by an unit vector pointing in its heading direction that starts at the agent position and has a fixed length (regardless of the value of R). The large overlaid blue arrows serve to guide the eye, showing the approximate mean heading direction of groups of agents. As in figures 6 and 7, the crosses indicate the midpoint between interacting neighbors with opposing self-propulsion forces.
These larger simulations show that the three types of transition have the same characteristics described in sections 4 and 5 for smaller systems. The ordered states in the left column present a similar aligned configuration for all values of R, with only a few points of opposing forces appearing and disappearing at random. The states of dynamic disorder and quenched disorder also display the same features identified in previous sections, which become more apparent in this larger system. For R = 0.01 and R = 0.2, the dynamic disorder state in panels (b) and (d) is characterized by the presence of multiple points of opposing forces that are continuously appearing and disappearing at random, homogeneously throughout the system. Instead, for R = 1 the quenched disorder state displayed in panel (f) does not evolve in time, with the mean orientation of each agent remaining constant even as the headings continue to fluctuate about that mean at every timestep due to noise. We can also deduce from this panel that the quenched disorder state was reached when the horizontal line of opposing self-propulsion forces first percolated through the system. This line divides locally ordered regions with different heading directions (indicated by the large overlaid arrows) that remain fixed in time, stabilizing this state of low global polarization. Finally, the corresponding simulation videos in the supplementary material show that here the lines of opposing forces display almost no change over time, in contrast to the smaller cases where finite size effects resulted in significant fluctuations.
The analysis of the equal-time correlation functions of the heading directions as a function of the distance between agents can provide additional information for understanding the different types of ordered and disordered states. We include the corresponding curves in the supplementary material for runs equivalent to those in figure 9. These allow us to directly distinguish between the different states: all ordered states display long-range correlations, in the dynamic disorder states the correlations decay quickly to zero, and in the quenched disorder states they oscillate from positive to negative at the scale of the orientation domains (showing that headings tend to be antiparallel at the lines of opposing forces). In future work, these correlation functions can be used to draw the complete phase diagrams as a function of all parameters in finite-size and large-scale systems.

Discussion and conclusions
In this paper, we introduced a new mechanical AEA model system of self-propelled components with elastic interactions that has a direct physical realization and can also be used to approximately describe a variety of active dynamics. This system self-organizes into collective motion for low levels of the noise D θ and displays a finite-size order-disorder transition, as D θ is increased, that can be of three different types, depending on the parameter values. For very large system sizes, only one of these transitions will remain, since in the other two the ordered state relies on velocity-based alignment forces that cannot achieve long-range order because the interaction network is fixed.
For R/l 0 ≈ 0, the interactions depend mainly on the relative positions between neighboring agents. We thus recover the discontinuous order-disorder transition to a state of dynamic disorder previously observed in the original AE model (equivalent to the R = 0 case), since the self-organization process appears to result from the same elasticity-based mechanism described in [29,30]. This phase transition appears to remain for very large systems, in the thermodynamic limit.
For larger values of R, the interactions depend strongly on the relative orientations. We thus obtain the continuous finite-size order-disorder transition into a state of dynamic disorder that had been previously observed in velocity-based models with alignment interactions and fixed topological neighbors [35], because the self-organization process appears to result from the distributed consensus mechanism on the velocities or heading angles described in [24,25]. We have observed that the critical noise will grow with R values, as the alignment elastic forces become stronger. This appears to be true for, both, the finite-size continuous transition and the large-scale discontinuous transition. For very large systems, the velocity-based ordered phase that leads to the continuous transition will disappear and only a discontinuous transition from the position-based ordered phase will survive.
For nonzero values of R and in certain parameter combinations, we also find a novel state of quenched disorder, in which the mean heading direction of each agent remains fixed in time. This state appears after a finite-size discontinuous order-disorder transition in small and medium-sized systems, and we expect it to emerge in a transition from dynamic disorder to quenched disorder in very large systems. At the transition point, lines of opposing self-propulsion forces are stabilized when they become long enough to percolate through the system. These lines divide regions of locally aligned agents that point in different directions from those in other regions, which results in low global polarization. Although the exact mechanism for reaching the state of quenched disorder is still unclear, the lines of opposing forces seem to nucleate spontaneously due to random fluctuations and to be stabilized by noise.
It is interesting to try to connect this AEA quenched disorder state to other frozen disorder states found in physical systems. When comparing to the frustrated quenched disorder states in spin glasses, for example, we note that no direct analogy can be established because the AEA model presents no geometrical frustration. Indeed, in the AEA case local alignment always corresponds to a simple energy minimum that can be reached by reducing the noise. Some similarities can be found when comparing instead to the jammed disordered states in granular media, since large-scale displacements that could lead to lower energy are also stopped by localized stable structures: force arches in granular media and lines of opposing forces in the AEA model. We note, however, that a fundamental difference with both of these cases is that the AEA quenched disorder state is only reached at high noise levels.
Much work remains to be done to further investigate this novel state of quenched disorder. For example, it would be relevant to analyze the nucleation and dynamics of the lines of opposing forces for a broad range of parameter combinations. Interestingly, some of these studies could be carried out analytically through a combination of stochastic calculations and elasticity theory. It would also be important to explore if this type of quenched disorder transition can be found in other models or in experimental systems. Indeed, since the dynamics and interactions of the AEA model are very general, all three types of transition described in this work could be observable, in principle, in a variety of biological or technological systems, such as active cell membranes or decentralized robot swarms.
Finally, given that the AEA model has a direct mechanical interpretation, it would be interesting to build it as a simple table-top experiment to explore its real-world features and capabilities. This could be done, for example, by using bristlebot-like agents connected by springs. It could even be constructed using 3D-printed passive components that advance through a similar ratchet effect, which would be linked into an elastic structure and placed on a shaken surface.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files). Figure A1. Polarization as a function of noise D θ in simulations of an AEA model with R = 0 (equivalent to the original AE model studied in [30]) for a system of N = 100 agents homogeneously distributed in a 10 × 10 hexagonal lattice with periodic boundary conditions. The blue line displays the mean, and the black crosses the local maxima, of the distribution of polarization p values in 20 independent simulations. We observe a discontinuous transition that is almost identical to that obtained in [30] for an hexagonal-shaped formation with free boundaries moving in an infinite arena.

Appendix A. Simulation timestep and validation
The AEA model was simulated by integrating equations (1) and (2) in time using a standard first-order forward Euler timestep, given by The numerical timestep was set for all simulations in this paper to Δt = 0.1. In preliminary analyses, we tested that our results did not significantly change when using smaller Δt values or a second-order Runge-Kutta method [39].
In contrast to previous studies of the original AE model [15,29,30] where agents were placed in a hexagonal-shaped formation with free boundaries in an infinite arena, here we filled the periodic arena with the agents. In order to verify that the original AE model behaves similarly in both cases, we simulated the AEA model with R = 0 for N = 100 agents in a 10 × 5 √ 3 periodic arena and compared them to the N = 91 hexagonal formation results in [30], using the same parameters: α = 0.01, β = 0.12, v 0 = 0.002, l 0 = 0.65, and k = 5. Figure A1 displays the resulting polarization values as a function of noise D θ . To generate this figure, we performed 20 simulations of 10 6 timesteps for each noise value, collecting the polarization p defined in equation (5) of the last 10 5 timesteps (to guarantee convergence to a final statistical steady state). We then computed the mean (blue line) and local maxima (black crosses) of the resulting distributions. Half of the simulations were initialized with random orientations and half perfectly aligned. All agents were initially placed at their equilibrium distances. By comparing this figure to figure 3 in [30], we can verify that the transition is almost identical, with the current case displaying only a slightly higher critical noise, which can be attributed to the stabilizing effect of the periodic boundary conditions.

Appendix B. Linear stability analysis
We present here the linear stability analysis of a minimal two-agent system that follows the AEA dynamics with zero noise. We will reduce the degrees of freedom to only two variables by imposing mirror symmetry to the stationary configurations and their perturbations. This will allows us to compute simple analytical expressions that can capture the dependency of the expected dynamics on the model parameters.
We note that one can develop a different linear stability analysis in the continuous limit, where agent positions are described by a displacement field, orientations by a spin field, and forces by the linear elasticity equations. However, this approach does not directly capture the transient oscillatory nature of some of the Figure B1. Schematic representation of a minimal two-agent system with symmetric configuration. As in figure 1, each agent is represented by a gray disc (with positions x 1 and x 2 ) that has a blue 'arm' projecting forward (with headingsn 1 andn 2 ). For the stability analysis, the equilibrium distance is defined as r 0 , its perturbations as Δr, the equilibrium angles as φ 0 , and their perturbations as Δφ/2. Note that we depict here arbitrary agent headings to define the variables generically; the actual equilibrium headings are always parallel or antiparallel. solutions described below (as shown in [29,30] for the AE model) and cannot easily describe the unstable equilibrium solutions that we will also analyze here. Since the purpose of these calculations is only to provide an intuitive approximate description of the stationary states, we will leave further exploration of the continuous approach for future work. Figure B1 presents the two-agent system with mirror symmetry and its corresponding variables. As in section 2, each agent i (gray circles) is located at x i and advances with preferred self-propulsion speed v 0 in its heading directionn i . The forward-projecting arms are represented by the blue bars but we do not display the linear spring connecting their tips, in order to simplify the diagram. The equilibrium positions and orientations computed below are defined by the distance r 0 and angles φ 0 , respectively, which are displayed for an arbitrary configuration in the figure. Their corresponding symmetric perturbations are given by Δr and Δφ/2.
In the symmetric stationary states that we consider, the agents can be either static or translating vertically (with respect to the figure) with constant speed. The dynamical equations (1) and (2)  As in section 2, here R is the arm length, k is the spring constant, l 0 its natural length, andn ⊥ 1 a unit vector perpendicular ton 1 .
When we substitute these expressions into equations (B.1) and (B.2), we find four possible equilibrium solutions. Two of these correspond to having both agents advancing in parallel upwards or downwards on figure B1, with K y = v 0 and φ 0 = +π/2, or K y = −v 0 and φ 0 = −π/2, respectively, and at a distance r 0 = l 0 from each other. The other two correspond to both agents oriented horizontally and immobile (K y = 0), either pushing or pulling in opposite directions (head-to-head or tail-to-tail), with φ 0 = 0 and r 0 = l 0 + 2R − v 0 αk , or φ 0 = π and r 0 = l 0 − 2R + v 0 αk , respectively. We now investigate the dynamics of perturbations over the distance Δr and angle Δφ about the four equilibrium solutions computed above. Since we are imposing symmetric motion, we can use equations (1) and (2) to express the dynamics of Δr and Δφ only in terms of F 1 andn 1 , obtaininġ Here, each θ i denotes the angle of its correspondingn i whereasx is a unit vector pointing along the horizontal axis in figure B1. The force and unit vectors associated to agent 1 for the perturbed variables are given by By replacing these expressions into equations (B.6) and (B.7) for each equilibrium solution, we can thus find its corresponding perturbation dynamics. In what follows, we will use this procedure to carry out the linear stability analysis for each equilibrium state. We first consider the two equivalent cases where the agents are advancing in parallel, either upwards or downwards, by substituting φ 0 = π/2 and r 0 = l 0 into the equations above to obtaiṅ From here, we can compute the eigenvalues of the resulting stability matrix, which are given by In section 3 we use this result to evaluate the stability of the polarized equilibrium solution as a function of the model parameters.
We now turn our attention to the two cases where the agents are oriented horizontally and are immobile, either pushing head-to-head or pulling tail-to-tail on each other. For the head-to-head case, we substitute φ 0 = 0 and r 0 = l 0 + 2R − v 0 αk , into equations (B.8)-(B.10) to obtaiṅ On the other hand, for the tail-to-tail case, we replace φ 0 = π and r 0 = l 0 − 2R + v 0 αk into the same equations, which giveṡ When we linearize both sets of expressions for small Δx and Δφ perturbations, we find that they produce the same system of equations, given by d dt From here, we can again compute the eigenvalues of the resulting stability matrix, which take the following two values Since λ + > 0, this equilibrium solution is always unstable and the value of λ + provides a measure of its degree of instability. In sections 4 and 5, we use the linear stability eigenvalues computed above for a minimal two-agent mirror system to deduce the dependence on the parameters of the states and dynamics for larger systems.