Tuning collective actuation of active solids by optimizing activity localization

Abstract

Active solids, more specifically elastic lattices embedded with polar active units, exhibit collective actuation when the elasto-active feedback, generically present in such systems, exceeds some critical value. The dynamics then condensates on a small fraction of the vibrational modes, the selection of which obeys non trivial rules rooted in the nonlinear part of the dynamics. So far, the complexity of the selection mechanism has limited the design of specific actuation. Here, we investigate numerically how localizing activity to a fraction of modes enables the selection of non-trivial collective actuation. We perform numerical simulations of an agent-based model on triangular and disordered lattices and vary the concentration and the localization of the active agents on the lattice nodes. Both contribute to the distribution of the elastic energy across the modes. We then introduce an algorithm, which, for a given fraction of active nodes, evolves the localization of the activity in such a way that the energy distribution on a few targeted modes is maximized – or minimized. We illustrate on a specific targeted actuation, how the algorithm performs as compared to manually chosen localization of the activity. While, in the case of the ordered lattice, a well-educated guess performs better than the algorithm, and the latter outperform the manual trials in the case of the disordered lattice. Finally, the analysis of the results in the case of the ordered lattice leads us to introduce a design principle based on a measure of the susceptibility of the modes to be activated along certain activation paths.

---

I. Introduction

A central goal of meta-material design is to realize multi-functionality, enabling a system to effectively perform a variety of tasks. Active solids, composed of elastically coupled active units that locally exert active forces, while being confined to the vicinity of well-defined reference positions, emerge as promising candidates for achieving this goal.

Correlated noise generated by an active bath is known to actuate nontrivial zero modes while suppressing harmonic modes to a degree dependent on temporal correlations.¹ This is the simplest evidence for the breakdown of equipartition in active solids. Active agents embedded in an elastic structure are further able to mobilize solid body motion² or a free-moving mechanism even in a topologically complex case.¹ Subsequently, the experimental and numerical evidence presented in ref. 3 revealed that the generic presence of a nonlinear elasto-active feedback of the elastic stress on the orientation of the active forces can induce selective and collective actuation of the solid: a collective oscillation of the lattice nodes around their equilibrium position emerges. Only a few elastic modes are actuated and crucially, they are not necessarily the lowest energy ones.

In the presence of several actuatable zero modes, whether trivially associated with solid body motion or more complex mechanisms, several dynamics coexist in phase space and a general formalism to describe the statistical evolution of collective motion has been derived.⁴ This coexistence can also hold in the case of mechanically stable solids, as illustrated experimentally by realizing a hysteretic tension-controlled switch between two actuation dynamics.⁵ Altogether active solids indeed offer a promising horizon for the design of multi-functional meta-materials. Furthermore, most of the active solids considered so far are ordered and have a spatially homogeneous distribution of active forces, leaving room for the large potential of alternative actuation strategies.

Exploring this opportunity is the main goal of the present work. To do so, we shall control the injection of energy in the lattice by taming the spatial distribution of activity in both ordered and disordered elastic lattices. More specifically we aim at answering the following question. To what extent, and with which guiding principles, can one design the spatial distribution of activity in the lattice, in order to activate some specific modes?

We perform numerical simulations of an agent based model on triangular and disordered lattices and vary the concentration and the localization of the active agents on the lattices’ nodes. We first show that, in sharp contrast with equilibrium solids, the distribution of energy in the elastic modes is very far from being equally distributed and can be controlled by the distribution of activity in the lattice. We then introduce an algorithm, which, for a given fraction of active nodes, evolves the localization of the activity in such a way that the energy distribution on a few targeted modes is maximized – or minimized. We illustrate on a specific targeted actuation, how the algorithm performs as compared to manually chosen localization of the activity. While in the case of the ordered lattice a well educated guess performs better than the algorithm, the latter outperform the manual trials in the case of the disordered lattice. Finally, the analysis of the results in the case of the ordered lattice leads us to introduce a simple design principle based on a measure of the susceptibility of the modes to be activated along certain activation paths.

The paper is organized as followed. In Section II, we introduce the agent models and the observables we will use throughout the paper. Section III is devoted to the characterization of the energy distribution among the modes for ordered and disordered networks of different sizes, varying the concentration of active nodes. Section IV explores the selection of actuated modes by tuning the distribution of active particles in a triangular lattice and leads us to propose an optimization algorithm in Section V, the performance of which is compared to the “manual” design, in the case of the ordered and disordered lattices. While the algorithm outperforms the manual trials in the case of the disordered lattice, a well educated guess performs better in the case of the ordered lattice, suggesting the existence of a simple design rule, which we propose in Section VI, before concluding.

II. Models and methods

A. Lattices

We consider two-dimensional elastic lattices at mechanical equilibrium, consisting of nodes with a well-defined reference configuration, connected by springs of stiffness κ, Fig. 1. The extremal nodes of the lattice are pinned to the lab frame. Both ordered and disordered lattices are considered. The ordered lattices are triangular lattices, composed of N = 1 + 3R(R − 1) nodes, located on R concentric hexagonal rings (see Fig. 1a for an example of this kind of lattice with R = 7). The disordered lattices are created by first generating a packing of soft discs at high density, following the protocol introduced in ref. 6 and 7 (see the Appendix). The center of each disc of the so obtained packing defines a node and two nodes i, j are connected by a spring whenever r_ij < R_j + R_i, where r_ij is the distance between nodes i and j and R_j is the radius of the disc j. To compare ordered and disordered lattices of the same sizes, large disordered lattices are generated, inside which we pin a hexagon of nodes at a distance d ≈ R + 1 from the center of the system in such a way to have approximately N(R) moving nodes within the pinned boundary (see Fig. 1b for an example of this kind of lattice with R = 7). All disordered networks used throughout this work have a high coordination number (z ≈ 6) and their density of states D(ω_k) – defined as the distribution of the frequencies ω_k – are statistically similar, for low frequencies, to the ordered case, as shown in Fig. 8 of Appendix A.

Fig. 1 Ordered and disordered lattices with randomly distributed driven nodes: (a) triangular lattice with N = 127 (R = 7) and 25 randomly distributed driven nodes, resulting in ϕ_dn = 0.2. (b) |A disordered lattice with N = 126 and ϕ_dn = 0.2. In both cases, the external layer of nodes is pinned, as indicated by dark grey points in the figure. White points represent empty nodes and red points thermally or activity driven nodes.

B. Dynamics

It was shown that the phenomenology of collective actuation is very well captured within the harmonic approximation, where the elastic forces are entirely encoded into a linear description F^el_i = −[Doublestruck M]_iju_j, with u_i the displacement of node i, and [Doublestruck M]_ij the dynamical matrix of the elastic lattice of interest. Here we follow the same scheme, but, in contrast with previous works, only a fraction ϕ_dn = n_dn/N of nodes are driven. We mostly consider active driving, but also introduce a thermally driven network for comparison.

In the active case, the driven nodes obey the self-aligning overdamped dynamics introduced in ref. 3. Each driven node is activated by a free to rotate self-aligning active particle that exerts a polar force on the node, while being reoriented by the total elastic force acting on that node. This nonlinear elasto-active feedback is the key ingredient responsible for the onset of collective actuation (see also ref. 8 for a review on self-aligning polar particles). The dynamical equations read:

u̇_i = Πn̂_i − [Doublestruck M]_iju_j, (1)

ṅ_i = n̂_i × u̇_i × n̂_i, (2)

where the unit vector n̂_i indicates the direction in which the particle sitting at node i exerts the active force. Π = l_e/l_a, the unique control parameter is the ratio of two lengths, the elastic length l_e = F_a/κ, which describes the elongation of a spring of stiffness κ under the action of an active force F_a and the alignment length l_a, which is the typical length a node i must be displaced to reorient the active particle sitting at this node. Eqn (2) describes the reorientation of particle i towards its displacement u̇_i according to the self-alignment mechanism. A noise term can be added to eqn (2), and we discuss its role in Appendix B. In this active setting, the non-driven nodes simply obey the same equation, with zero activity, that is Π = 0.

In the thermal case, it was chosen that the driven nodes obey the standard underdamped Langevin equation, since it estimates the thermodynamics properties better than the overdamped case,

where ξ_i(t) is a Gaussian white noise with zero mean 〈ξ_l(t)〉 = 0, 〈ξ_l(t)·ξ_k(t′)〉 = δ(t − t′)δ_lk. T_eff controls the noise amplitude. The non-driven modes simply respond elastically and transfer the elastic forces around the network obeying newtons law.

Numerics: In the case of active driving the dynamical equations were integrated using a Runge–Kutta 4th order method with dt = 0.01. For the thermal driving, the stochastic equations were integrated using the stochastic velocity Verlet algorithm⁹ (also dt = 0.01). The measures of interest are taken on the interval of t ∈ [250, 500] time steps (Δt = 250), disregarding the transient regime.

C. Observable

We are primarily interested in the distribution of the elastic energy among the vibrational modes of the lattices. These modes – also called normal modes (NM) – are the eigenvectors |φ_k〉 of the 2N × 2N symmetric dynamical matrix and form a complete orthonormal basis. The lattices considered here are all mechanically stable, they are associated with strictly positive eigenvalues ω_k². In the case of the triangular lattice, it is convenient to sort these modes by the four classes of rotational symmetries, which we will simply denote as the symmetry of class 1, 2, 3, and 4, that leave the lattice and its normal modes invariant. Examples of the normal modes classified by their classes of symmetries are shown in the Appendix E for a system with R = 7 (see also the ESI in ref. 3 for an explicit construction). The dynamics is then decomposed on these modes, by projecting the 2N-dimensional displacement field |u〉 = {u₁, u₂,…, u_N} on each mode: P_k(t) = 〈φ_k|u(t)〉. In practice, P_k(t) is averaged over time in the steady state and normalized as follows:where x̄ denotes a temporal average over 250 time steps (disregarding the first 250 steps).

III. Energy distribution for a varying fraction of driven nodes

We first examine the distribution of energy among the modes when the driven nodes are randomly distributed in the lattice. We consider both disordered and ordered lattices and vary both the fraction of driven nodes ϕ_dn and the total number of nodes N. In the first order, the elastic energy of the system can be written in terms of the matrix [Doublestruck M]:

At equilibrium, equipartition dictates that each quadratic degree of freedom contributes equally to the system's energy, typically k_BT/2 per degree of freedom, where k_B is the Boltzmann constant. In the present context, this imposes that each term of the above decomposition contributes an equivalent amount of energy, so that P_k² ∝ ω_k⁻². This is precisely what is reported in Fig. 2a for the thermally driven case, when all nodes are driven.

Fig. 2 Distribution of energy in thermally and actively driven lattices: (a) Average squared amplitude of the modes vs. ω_k for a thermally driven, an actively driven ordered, and an actively driven disordered lattice as reported in the legend. Each point corresponds to one mode. The dashed line corresponds to the fit . While α = 2 in the thermal case, it is much larger in the active case. N = 127, ϕ_dn = 100%. (b) Dependence of α on the fraction of driven nodes ϕ_dn for N = 127. (c) Dependence of α on the system size N for ϕ_dn = 100%. The data are averaged over 50 initial conditions and 12 different disordered lattices were used (Π = 1.3).

In the case of active driving, equipartition does not hold. We start by simulating lattices of size N = 127, for which all nodes are actively driven (ϕ_dn = 100%). Quite remarkably, one still observes a power-law dependence P²_k ∝ ω^−α_k, for large enough ω_k, albeit with a much larger value of α, indicating a stronger condensation of the energy on the low energy modes, the effect being more pronounced in the case of ordered lattices. A strong condensation in the softest modes was also observed in the context of a jammed active system.¹⁰ These results are confirmed, whether varying the fraction of driven nodes (Fig. 2b), or the lattice size (Fig. 2c). The exponent α is very robust with respect to the fraction of driven nodes. Even for a fraction as low as ϕ_dn = 10%, α ≃ 4 remains twice larger than its equilibrium counterpart. We also note that in the case of an ordered lattice, increasing the system size amplifies the condensation, with α reaching values as large as close to 7, where it seemingly saturates. This is not the case for disordered lattices, where α ≃ 4.5 for all system sizes probed here. In all cases, the violation of equipartition reported above opens the path for manipulating the energy injection, in view of actuating modes preferentially.

IV. Mode selection by tuning the spatial distribution of active nodes

We first concentrate on the ordered lattice and investigate how the spatial distribution of the actively driven nodes condition the energy distribution amongst the modes, the modes being grouped according to their four classes of symmetries (see Appendix E). We consider four distinct spatial organization of the actively driven modes, corresponding to the four columns of Fig. 3: (a) all active, (b) the 4th ring only is active, (c) half of the 4th ring is active, and (d) a centerline is active. These choices are somehow arbitrary but illustrate well the role of the spatial localization of the active nodes. The results are shown for increasing values of the activity, as indicated by the dimensionless parameter Π.

Fig. 3 Selective actuation of modes by symmetry classes, using a different spatial distribution of active nodes: The top row indicates the localization of the active nodes. Black points represent pinned edges, red points active nodes and open circles are empty nodes. In (a) active particles are located in all the nodes of the lattice, (b) in the 4th ring, (c) in a semi-circle of the 4th ring and (d) in a horizontal line. The bottom row indicates the energy repartition among the modes grouped by symmetry classes, as indicated by the colors of the legend. The extreme right column (e) indicates this distribution when all nodes are thermally driven, as a reference. For each spatial distribution (column (a), (b), (c) and (d)), the distribution of the energy is represented for increasing values of Π. When Π is too small, collective actuation does not take place and the system remains frozen, as indicated by the blue overlay. The data shown here are averaged over 100 runs with random initial polarization; lattice size N = 127.

As a first observation, we note that the active driving leads to very different distributions of the elastic energy across the modes of different symmetry classes. When all nodes are driven (a), the modes of class 1 and 2 vastly dominate the dynamics, as compared to the passive case, where the equipartition of energy favors the modes of class 1 and 3. Second, one sees that the distribution of energy is only slightly changed when reducing the activation to a unique line on the fourth ring (b). Conversely the distribution of energy among the four classes of symmetry is strongly altered when only half of the ring is driven (c), or when only the central line is driven (d) pointing at an important role of the symmetry of the spatial distribution of the driven nodes. We also note that, depending on the geometrical organization of the driven nodes, the dependence on Π can be minimal, as in case (a), (b), and (d), or pretty strong as in case (c).

These observations drive us to search for spatial distributions of the active nodes that enhance the projection of the dynamics in a desired class of symmetry or, even, in a few specific numbers of modes.

V. Optimization algorithm

Can one determine a spatial distribution of active particles that effectively amplifies the dynamics within a desired class of symmetry or a specific normal mode? To address this question, we propose an algorithm that combines molecular dynamics simulations of the dynamics and the Metropolis Monte Carlo method to evolve the active configurations of the lattice.

A spatial distribution of actively driven nodes is denoted by a vector |σ〉 = {0, 1, 1, 0,…} of size N, where the one's indicate the nodes that are actively driven, and the zero's the non-driven nodes.

This configuration is evaluated by running the dynamics during 2Δt = 200 time steps and computing the cost function , where the temporal averaged runs over the Δt = 100 time steps composing the second half of the simulation time window (the choice of this time window is discussed in Appendix D). C_k^|σ〉 evaluates the ability of the configuration |σ〉 to concentrate the energy of the system in a given mode k. Every 2Δt time steps, a Monte-Carlo move is proposed from the configuration |σ〉 to another configuration |σ′〉, by changing the location of one active node, hence keeping the overall fraction of active nodes constant (see Fig. 4). The goal being to maximize, and respectively minimize, the cost function, depending on whether one wants to increase, and respectively decrease, the projection on a given mode, and the new configuration |σ′〉 is accepted with probability:

where T_e is an effective temperature which allows for the exploration of the configurations space. If the new configuration |σ′〉 is accepted, the next step starts from it; if not, the original configuration |σ〉 is restored and another move is proposed. A Monte Carlo Step (MCS) is defined as N trials configurations. We use T_e = 10⁻³ and test the dependence of our results on T_e in Appendix C, Fig. 11.

Fig. 4 Monte Carlo move: schematic example of two configurations, |σ〉 and |σ′〉, which differ by one Monte Carlo step, where an active node is made inactive in favor of another node, keeping the fraction ϕ_dn constant.

The algorithm is evaluated on its ability to find a configuration, for which the spatial distribution of 24 active nodes among 127 nodes maximizes the condensation of the dynamics on the mode k = 3, as indicated by the value of . The mode k = 3 belongs to the class of symmetry 2 and is therefore not predominantly actuated in a typical configuration. We recall that mode 3 is shown in Fig. 7(a) and (b). The performance of the configuration optimized by the algorithm after 1200 Monte Carlo steps (last column on the right) is compared to configurations, where the same number of active nodes are localized in preset geometries (first six columns) as indicated by the drawings in red, on the top row of Fig. 5: randomly chosen around all the lattice (uniform random); localized on two layers r = {1,6}, r = {2,5}, and r = {3,4}; on two separated areas, |x| > 3 which do not respect the rotational symmetries of the lattice; with the 4th layer fully occupied, r = 4. The top, respectively the bottom, row displays the distributions obtained for the ordered, respectively the disordered, lattices for 30 different initial conditions. In the case of the disordered network, we generated 30 different lattices – except for the optimized case where we used 10 – and, for each one, 30 initial conditions were simulated.

Fig. 5 Performance of the optimization algorithm: Violin plot showing the distributions of the final projections onto mode 3, denoted , obtained after 1200 Monte Carlo steps (last column on the right). These distributions are compared across 6 preset configurations with spatial distributions of the active nodes, indicated in red on the top row. The white points represent the mean of the distribution, while the black bars indicate the dispersion around the mean. The density of active nodes, ϕ_dn = 24/127 is identical across all cases. The results are shown for both ordered hexagonal lattices (top) and disordered lattices (bottom). In the ordered networks, 30 distinct initial conditions (ICs) were used for each configuration, except for the optimization case, for which were used 350 distinct ICs. In the disordered networks, 30 ICs were simulated for 30 different networks in all configurations, except for the optimized cases, which involved only 10 networks. R = 7, N = 127 and Π = 1.3.

In general, the algorithm performs much better than the preset configurations. In the disordered lattices this is always the case. In the case of the ordered hexagonal lattice, the activation of the 4th layer, which precisely matches the location of maximal polarization inside the mode k = 3, performs better. Note that the dynamics being deterministic, a fraction of the initial condition leads to complete failure of the activation of the third mode, when the spatial distribution of the active nodes is preset, while the optimization algorithm allows for an adaptation of the spatial distribution of the active nodes to the randomly chosen initial condition.

In Fig. 6, we present the convergence of the algorithm for the case where the goal is to maximize the projection onto mode 3. We also analyze the distribution of active particles in the lattices as a function of the algorithm's convergence. Fig. 6 displays the evolution of the performance of the algorithm as a function of the number of Monte Carlo steps for three different lattices: the ordered hexagonal one, a typical disordered lattice, and the disordered lattice, for which the algorithm obtains the best performance. For each of them the figure shows the mean performance, averaged over 100 initial conditions, while the worst and best cases delimitate the colored areas. From these evolutions one sees that disordered networks can achieve a better performance than the ordered one, in the sense that the optimization algorithm identifies a spatial distribution of the active nodes that favors a stronger condensation of the dynamics on the selected mode of interest. From the dynamical evolution of the localization of the active nodes during the optimization process it appears that the algorithm has a hard time in identifying the best configuration in the case of the ordered lattice. More specifically, it only very slowly condensates the active nodes toward the 4th layer of the lattice, in contrast with the disordered case, where this condensation takes place in the early steps of the optimization. We interpret this behavior as a sign that the symmetries of the ordered lattice favors the nearly degenerescence of many configurations with respect to their evaluation by the cost function.

Fig. 6 Convergence of the algorithm performance: (a) evolution of the cost function for the projection on the mode k = 3, , with the number of Monte Carlo steps in three different networks: the ordered hexagonal lattice (green and white dotted line), a randomly chosen disordered lattice (yellow and black dashed line) and the disordered lattice, for which the algorithm obtains the best performance (orange and black dashed line). The lines indicate the mean performance over 100 random initial conditions, while the colored areas are delimited by the worst and best cases. (b) Spatial distribution of the active nodes as a function of their distance to the central one at different moment of the optimization process (as indicated in the legend) for the three networks (as labeled on the right of each panel); (R = 7, N = 127; Π = 1.3).

VI. Optimization rule for the selective actuation of hexagonal lattices

The analysis of the results from the optimization algorithm suggests a possible simple strategy for optimizing the spatial distribution of the active nodes in order to achieve the condensation of the dynamics on some specific modes. We observed that the optimal configuration to amplify the projection of the dynamics on the third mode of the hexagonal lattice, for a network size with R = 7, corresponds to a localization of the active nodes on the fourth layer of the lattice. As noted already, this coincides with the regions of largest magnitude of the displacement field in mode k = 3 (see Appendix E).

We propose to generalize this observation and organize the spatial distribution of the active nodes in light of the displacements field geometry of the mode of interest. The underlying hypothesis is that localizing the active nodes in the regions of high displacement of a mode should favor the coupling of the activity to that mode. To do so we first define an activation path L on the lattice, composed of |L| adjacent nodes and define S_k(L) the activation susceptibility of that path in mode k as the local projection of the polarization field on the tangent to this path:

where l indexes the nodes along the path L, ê_l is the unitary vector tangent to the path at node l, and (φ_k)_l is the displacement of mode k at node l. This susceptibility is defined such that it is maximal when the path L runs parallel to the local maxima of the displacement field of the mode of interest. The configurations of active lines transport the alignment information along the path backward and forward, favoring oscillations to occur parallel to the path.

Fig. 7 illustrates the application of the above ideas in the case of an ordered hexagonal lattice with R = 14 layers, for which two types of active paths are tested: (i) concentric active rings of size r(L(r)), with |L(r)| = 6r and ê_l = ê_θ(l), see Fig. 7c and (ii) linear radial paths that do not cross the center (L(θ), with |L(θ)| = R − 2 and ê_l = ê_r(l), see Fig. 7d). Fig. 7(c) and (d) show the susceptibility S_k(r), for the concentric ring path of radius r and S_k(θ) for the radial path with orientation θ for the five modes k = 3, 5, 9, 13, 14 displayed on the top rows. The core observation is that the susceptibility strongly depends on the combination of the path and the mode. For instance, S₉(r) is negligible for all concentric paths, while S₉(θ) is systematically large for all radial path. The dependence on the path can be relatively simple, as it is in the case for S₃(r), where one recovers the observation made earlier that the mode k = 3 has a maximal susceptibility when the distribution of active nodes concentrate on a ring of radius r = R/2. But it can also be less obvious for higher modes with less symmetric displacements, such as mode k = 5 or k = 14. One can also note that, once a path is chosen, there are several modes with susceptibility greater than zero. This implies that these modes would be activated if the nodes were excited along that path. Therefore, optimizing for a single mode is generally likely to be impossible.

Fig. 7 Optimization of the selective actuation of hexagonal lattices: (a) and (b) five specific modes of the triangular lattice, of size N = 547, R = 14, (k = 3, 5, 9, 13, 14), with (a) the amplitude and (b) the vector field of the local polarization of the mode (red, respectively blue, regions indicate large, respectively small, displacement in the mode polarization). (c) The susceptibility S_k(r) as a function of r and (d) S_k(θ) as a function of θ for the five selected modes of interest as indicated in the legend. (e)–(h) Respective condensation of the energy on modes k, colored by their class of symmetry, as indicated in the legend of panel (f), for four configurations of active nodes, as indicated in each panel by the sketch of the activation path in red; (e) a ring on r = 3; (f) a ring on r = 7; (g) two segments on θ = 0° and θ = 180°.

Fig. 7(e)–(g) show the actual distribution of energy among the modes for four different paths. For the concentric ring path, one verifies very clearly the strong selection of mode k = 13 by the ring of radius r = 3 (Fig. 7e) and the even stronger one for the mode k = 3, when active nodes are distributed along the ring of radius r = 7 (Fig. 7f). The case of the radial paths confirms that the mode k = 9 is excited for both path configurations. More interestingly one sees that the radial path along θ = 0°, 180° is unable to select differentially the modes k = 5 and the mode k = 14, while the path along θ = 0°, 180° does select the mode k = 14, without activating the mode k = 5.

Altogether, the use of the activation susceptibility is thus a good design principle, although it clearly also reveals the limitation of the selectivity that can be reached. Nevertheless, we observe that the larger the system is, the more selective the activation design can be, especially for modes of high energy, because of the larger specificity of their polarization geometry and the possibility of combining several activation paths.

VII. Conclusions

In this paper, we explored different strategies for injecting energy into elastic lattices by exciting their nodes with active agents. We investigated both ordered triangular lattices and disordered lattices with a coordination number of approximately z ≈ 6, similar to that of the ordered lattice. By distributing the active nodes with various spatial organization, we demonstrate the possibility of tuning, at least to some extent, the energy partition amon the models. This of course contrasts with the classical scenario of equipartition imposed by thermal equilibrium.

When all lattice nodes are active (ϕ_dn = 1), a pronounced concentration of energy is observed on lower-frequency modes, as previously observed both numerically and experimentally.^1,3 More specifically, the amplitude distribution of the energy among the modes exhibits a power-law decay at large enough frequency, , with an exponent α significantly larger than 2, the value indicative of equipartition. When reducing the fraction of active nodes, the distribution of energy among the modes is heavily influenced by the spatial distribution of the active nodes inside the lattice. While complete control over energy concentration in specific modes remains elusive, we nevertheless demonstrate that there is room for optimization.

In the case of an ordered lattice, a well educated guess is possible and we propose a simple design principle, which identifies optimal paths for the distribution of the active nodes, on the basis of the geometry of the spatial structure of the modes, characterizing the path by their activation susceptibility. In the case of a disordered lattice, specific spatial distribution of the active nodes, combined with specific realization of the disorder, can better condensate the energy on a class of modes of interest, as compared to an ordered lattice with the same size. This suggests that the optimization algorithm has greater potential for evolution toward an “optimal distribution” in the presence of disorder. An interesting perspective would be to jointly optimize for the disorder of the network and the spatial distribution of the active nodes.

When considering disordered lattices, our focus has primarily been on cases where the coordination number is high, approximately z ≈ 6 and the lattice is largely hyperstatic and mechanically stable. It would be valuable to extend our investigation to lattices with lower coordination numbers, approaching z → 4. In this limit, it is known that the lattices are on the verge of losing their mechanical stability,¹¹ leading to an abundance of low-frequency normal modes compared to crystalline structures. The properties of these modes have been extensively studied in disordered solids,^12,13 with correlations drawn to solid rearrangements.^14,15 They also have been proposed as elemental defects controlling the flow of disordered solids.¹⁶ It would be intriguing to investigate whether active particles distributed in specific lattice regions could selectively excite the desired low-frequency mode intervals or unveil unpredictable features absent from the current work.

Conflicts of interest

There are no conflicts of interest to declare.

Appendices

Appendix A: Setting the disordered lattices

Disordered networks are created by first generating packings of soft discs at high density. To generate these packings, we employ a protocol used previously which is here summarized. We initiate by setting the density of particles ρ = N_p/V = 0.5 (where V represents the volume of the system and N_p stands for the number of particles), followed by conducting a high-temperature equilibration (T = 1.0) of the system under the influence of the potential energy , where r_ij is the distance between the centers of particles i and j, R_j is the radius of the particle j and H(x) is the Heaviside function. Subsequently, we employ the FIRE algorithm¹⁷ to minimize the potential energy U. Throughout this minimization process, we maintain a constant pressure of p = 1.0 using a Berendsen barostat¹⁸ with a time constant τ_Ber = 10.0. The minimization procedure is stopped when the interparticle force falls below 10⁻¹.

The ordered and disordered networks can be compared in terms of the frequency distribution of their normal modes: Fig. 8 shows the density of states D(ω), defined as the distribution of the frequencies ω, that are the square root of the dynamical matrix ([Doublestruck M]) eigenvalues, for both ordered and disordered networks. For small frequencies (ω < 1), both types of lattices have similar distributions.

Fig. 8 Density of states of the triangular lattice compared to the disordered case for several lattice sizes, showing that they are all statistically the same.

Additionally, the definition of the dynamical matrix does not account for any tension in the lattice. Fig. 9 displays the average projection of the dynamics onto the normal modes of the lattice, categorized into four rotational symmetries. This is shown for the case where active particles occupy all the lattice nodes, as depicted in the figure's inset. The key difference between the two cases is the presence of tension in the lattice, which has been demonstrated to be a relevant control parameter in promoting different types of dynamics in active solids.⁵ In Fig. 9a, the case with high tension is presented, showing that regardless of the value of Π, the dynamics project almost entirely onto symmetry type 1, recovering the results shown in ref. 3. When the tension is reduced, additional symmetries emerge in the dynamics, as shown in Fig. 9b.

Fig. 9 Average projection of the dynamics on the normal modes of the lattice, divided into 4 rotational symmetries, as indicated in the legend, for a ranging value of π. In (a) the system is tensioned and in (b) the tension is zero. The blue line indicates the percentage of moving systems for each π. In both cases, the system size is R = 4 and all the nodes of the lattice are full with active particles. The black points represent pinned edges, the red points active nodes and the open circles are empty nodes. The data shown here are an average over 100 runs of dynamics starting with different polarity directions.

Appendix B: Noise influence

To understand the implications of non-determinism on our system, we have also tested the convergence of the time series projection with the addition of noise. The noise is taken into account in eqn (2) (while eqn (1) remains unchanged), where it becomes:for D = αγ²/τ²κ², where α is the noise amplitude, γ is the drag coefficient, τ is the characteristic alignment time of the agent, κ is the spring stiffness and ξ is a unitary normal random variable.

In Fig. 10 are shown the average values of P₃² which were considered in a system with size R = 7 and a density of active nodes of 24/127 – the same as in Fig. 5 of the manuscript – and the distribution of the 24 driven active nodes in the ring r = 4 for 30 distinct initial conditions. In Fig. 10 we observe that the average projection of the dynamics on the third mode is robust to the addition of noise, even for noise with an order of one.

Fig. 10 Influence of noise: average projections over the 3rd mode for 30 initial conditions, the system has size R = 7, and the active distribution (with ϕ_dn = 0.19) is concentrated on r = 4, for Π = 1.3.

Appendix C Robustness of effective temperature, T_e

To adjust the temperature like parameter T_e used in the Monte Carlo rules for the optimization algorithm (eqn (6)), we explored three logarithmic scales of T_e. Fig. 11 displays the convergence of the algorithm for T_e = 10⁻³, T_e = 10⁻², and T_e = 10⁻¹. While too large T_e = 10⁻¹ ruins the optimization dynamics, we don't see the strong effect of T_e in the range [10⁻³, 10⁻²].

Fig. 11 Effective temperature, T_e, used on the optimization algorithm. On the graph with points and solid lines are represented the average of 100 optimizations and the shade areas are limited by the best (upper) and worst (lower) samples.

Appendix D Optimization algorithm time window

To test the robustness of the projection with respect to the time interval of the measurement, we performed the following test. We considered a case with R = 7 and a density of active nodes of 24/127, where the 24 active particles are randomly distributed in the system. We attempted to change the location of an active particle every 2Δt = 200 (molecular dynamics time steps) to optimize the projection of the dynamics onto mode 3 and summarized the results in Fig. 12.

Fig. 12 Optimization time series: projections over the 3rd mode for 5 initial conditions, the system has size R = 7, and random active distribution (with ϕ_dn = 24/127 = 0.19), for Π = 1.3. The blue intervals show the time window for measurements Δt = 100 and the white refers to intervals where no measurements are performed.

Each red plot in this figure corresponds to the projection of the dynamics onto mode 3 for one initial condition (IC), while the black line represents the average over 5 ICs. The white regions indicate periods during which no measurements were taken, and the blue regions represent the measurement intervals of Δt = 100 time steps. Note that the dynamics quickly reach a steady state shortly after the simulation begins. Since only small changes are made every 2Δt = 200 (with only one particle being moved while velocity and polarization parameters are maintained from one trial to the next), the dynamics converge rapidly, justifying the use of such measurement intervals.

Appendix E Modes of the ordered lattices with their class of symmetries

The triangular lattice with hexagonal boundaries respects symmetries over rotations of angle π/3, done by the operator Θ, and reflections Σ (e.g., across the axis y = 0). The eigenvectors of the dynamical matrix, [Doublestruck M], can be either direct eigenvectors of those symmetry operators or bases for inner spaces that respect those symmetries, with reflection eigenvalues given by σ = ±1 and rotation eigenvalues given by θ = exp(ikπ/3) for k ∈ {−2,…,3}.

The eigenvectors corresponding to the complex eigenvalues of the rotational operator are also complex and occur in degenerated pairs, |φ_±〉 relative to θ_± = exp(±inπ/3) for some n. These paired modes can be combined into two real modes |φ_l〉 and |φ_m〉, which also have the same energy as |φ_±〉. Although |φ_l〉 and |φ_m〉 (eigenvectors of [Doublestruck M]) are not eigenvectors of Θ, the 2-dimensional space they span remains invariant under rotations of this kind. The effect of Θ on these modes is characterized by 〈φ_l|Θ|φ_l〉 = 〈φ_m|Θ|φ_m〉, which corresponds to the real part of the eigenvalue of |φ_±〉. Therefore, the complete symmetry of a normal mode |φ_k〉 is determined by two real numbers:

〈φ_k|Θ|φ_k〉 ∈ {1, 1/2, −1/2, −1} and 〈φ_k|Σ|φ_k〉 ∈ {1, −1}.

In this work the normal modes of the dynamical matrix are characterized by the real part of the rotational symmetry eigenvalues they relate to. More precisely, the symmetries 1, 2, 3 and 4, used throughout the whole text, are related respectively to Real(θ) = 0.5, 1, −0.5 and −1. More details on the symmetry class derivation can be found in the ESI of ref. 3.

In Fig. 13 are presented the first 24 normal modes for the ordered triangular lattice of size R = 7, the colors identify the class of rotational symmetry they follow.

Fig. 13 Normal modes displacement fields and rotational symmetries for the modes from 1 to 24.

Data availability

The data that support the findings of this study are available from the corresponding author, DL, upon reasonable request.

Acknowledgements

We thank Paul Baconnier for insightful discussions and the help provided on simulations. DL and CB thanks CNPq and CAPES for partially financing this study. We thank the supercomputing laboratory at New York University (NYU-HPC), where part of the simulations were run, for the computer time.

References

1. F. G. Woodhouse, H. Ronellenfitsch and J. Dunkel, Phys. Rev. Lett., 2018, 121, 178001.
2. E. Ferrante, A. E. Turgut, M. Dorigo and C. Huepe, Phys. Rev. Lett., 2013, 111, 268302.
3. P. Baconnier, D. Shohat, C. H. López, C. Coulais, V. Démery, G. Düring and O. Dauchot, Nat. Phys., 2022, 18, 1234.
4. C. Hernáandez-López, P. Baconnier, C. Coulais, O. Dauchot and G. Düring, Phys. Rev. Lett., 2024, 132, 238303.
5. P. Baconnier, D. Shohat and O. Dauchot, Phys. Rev. Lett., 2023, 130, 028201.
6. C. Brito, E. Lerner and M. Wyart, Phys. Rev. X, 2018, 8, 031050.
7. G. Kapteijns, W. Ji, C. Brito, M. Wyart and E. Lerner, Phys. Rev. E, 2019, 99, 012106.
8. P. Baconnier, O. Dauchot, V. Démery, G. Düring, S. Henkes, C. Huepe and A. Shee, arXiv, 2024, preprint,  DOI:10.48550/arXiv.2403.10151.
9. N. Grønbech-Jensen, Mol. Phys., 2020, 118, e1662506.
10. S. Henkes, Y. Fily and M. C. Marchetti, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 84, 040301.
11. S. Alexander, Phys. Rep., 1998, 296, 65.
12. M. Wyart, L. Silbert, S. Nagel and T. Witten, Phys. Rev. E: Stat., Nonlinear, Soft Matter. Phys., 2005, 72, 051306.
13. E. Lerner, G. Düring and E. Bouchbinder, Phys. Rev. Lett., 2016, 117, 035501.
14. C. Brito and M. Wyart, J. Chem. Phys., 2009, 131, 149.
15. A. Widmer-Cooper, H. Perry, P. Harrowell and D. Reichman, Nat. Phys., 2008, 4, 711.
16. M. Manning and A. Liu, Phys. Rev. Lett., 2011, 107, 108302.
17. E. Bitzek, P. Koskinen, F. Gähler, M. Moseler and P. Gumbsch, Phys. Rev. Lett., 2006, 97, 170201.
18. H. J. C. Berendsen, J. P. M. Postma, W. F. van Gunsteren, A. DiNola and J. R. Haak, J. Chem. Phys., 1984, 81, 3684.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4sm00868e

---