Programming active cohesive granular matter with mechanically induced phase changes

Interactions of simple robots can be mapped to a discrete algorithm that accurately predicts the ensemble’s emergent behaviors.


Introduction
Self-organizing collective behaviors are found throughout nature, including shoals of fish aggregating to intimidate predators [1], fire ants forming rafts to survive floods [2], and bacteria forming biofilms to share nutrients when they are metabolically stressed [3].Inspired by such systems, researchers in swarm robotics and programmable active mat-ter have used many approaches towards enabling ensembles of simple, independent units to cooperatively accomplish complex tasks [4][5][6].Both control theoretic and distributed computing approaches have achieved some success, but often rely critically on robots computing and communicating complex state information, requiring relatively sophisticated hardware that can be prohibitive at small scales [7,8].Alternatively, statistical physics approaches model swarms as systems being driven away from thermal equilibrium by robot interactions and movements (see, e.g., [9,10]).Tools from statistical physics such as the Langevin and Fokker-Planck equations can then be used to analyze the mesoscopic and macroscopic system behaviors [11].Current approaches present inherent tradeoffs, especially as individual robots become smaller and have limited functional capabilities [12,13] or approach the thermodynamic limits of computing and power [14].
To apply to a general class of micro-or nano-scale devices with limited capabilities, we focus on systems of autonomous, self-actuated entities that utilize strictly local interactions to induce macroscale behaviors.Two behaviors of interest are dynamic free aggregation, where agents gather together without preference for a specific aggregation site (see Section 3.2.1 of [5]), and dispersion, its inverse.These problems are widely studied, but most work either considers robots or models with relatively powerful capabilities -e.g., persistent memory for complex state information [15,16] or long-range communication and sensing [17][18][19] -or lack rigorous mathematical foundations explaining the generality and limitations of their results as sizes scale [20][21][22].Recent studies on active interacting particles [23] and inertial, self-organizing robots [24] employ physical models to treat aggregation and clustering behaviors, but neither prove behavior guarantees that scale with system size and volume.Supersmarticle ensembles [25] are significantly more complex, exhibiting many transient behavioral patterns stemming from their many degrees of freedom and chaotic interactions, making them less amenable to rigorous algorithmic analysis.
Here we take a two-pronged approach to understanding the fundamental principles of programming task-oriented matter that can be implemented across scales without requiring sophisticated hardware or traditional computation that leverages the physics of local interactions.We use a theoretical abstraction of self-organizing particle systems (SOPS), where we can design and rigorously analyze simple distributed algorithms to accomplish specific goals that are flexible and robust to errors.We then build a new system of deliberately rudimentary active "cohesive granular robots" (which, to honor granular physics pioneer Robert Behringer, we call "BOBbots" for Behaving, Organizing, Buzzing robots) to test whether the theoretical predictions can be realized in a real-world damped driven system.Remarkably, the lattice based equilibrium model quantitatively captures the aggregation dynamics of the robots.With a provable algorithmic model and even simpler BOBbots capturing the algorithm's essential rules, we next explore how contact stress sensing -a capability that is readily available in the robotic platform but interestingly not easily computable by a strictly local, distributed algorithm -can enhance aggregation performance, as suggested by insights from the theoretical model.This complementary approach demonstrates a new integration of the fields of distributed algorithms, active matter, and granular physics that navigates a translation from theoretical abstraction to practice, utilizing methodologies inherent to each field.

Aggregation algorithm
While many systems use interparticle attraction and sterical exclusion to achieve system-wide aggregation and interparticle repulsion to achieve dispersion, these methods typically use some long-range sensing and tend to be nonrigorous, lacking formal proofs guaranteeing desirable system behavior.To better understand these collective behaviors, the abstract model of self-organizing particle systems (SOPS) allows us to define a formal distributed algorithm and rigorously quantify long-term behavior.Particles in a SOPS exist on the nodes (or vertices) of a lattice, with at most one particle per node, and move between nodes along lattice edges.Each particle is anonymous (unlabeled), interacts only with particles occupying adjacent lattice nodes, and does not have access to any global information such as a coordinate system or the total number of particles.
In earlier work, Cannon et al. [26] analyzed a distributed SOPS algorithm for aggregation and dispersion under the assumption that the particle system remained simply connected (i.e., the system forms a single connected cluster with no holes).This SOPS algorithm defines a finite Markov chain with local moves that connect the state space of all simply connected configurations of particles.Moves are defined so that each particle, when activated by its own Poisson clock (i.e., after a delay chosen at random from a Poisson distribution with constant mean), chooses a random neighboring node and moves there with a probability that is a function of the number of neighbors in the current and new positions provided the node is unoccupied and the move satisfies local conditions that guarantee the configuration stays simply connected.In particular, for configurations σ and τ differing by the move of a single particle p along a lattice edge, the transition probability is defined as P (σ, τ ) ∝ min(1, λ n −n ), where λ > 0 is a bias parameter that is an input to the algorithm, n is the number of neighbors of p in σ and n is the number of neighbors of p in τ .These probabilities arise from the celebrated Metropolis-Hastings algorithm [27,28] and are defined so that the Markov chain converges to a unique Boltzmann distribution π such that π(σ) is proportional to λ E(σ) , where E(σ) is the number of nearest neighbor pairs in σ (i.e., those pairs that are adjacent on the lattice).
It was shown in [26] that the connected SOPS ensemble provably aggregates into a compact conformation when λ > 3.42 and expands to a conformation with nearly maximal (linear) perimeter when λ < 2.17 with high probability, i.e., with a probability of failure that is exponentially small in N , the number of particles.However, despite rigorously achieving both aggregation and dispersion, this distributed algorithm has two notable drawbacks that make it infeasible for direct implementation in a physical system of simple robots: the connectivity requirement that tethers the particles together and the "look ahead" requirement used to calculate transition probabilities ensuring convergence to the desired Boltzmann distribution.
To address these issues, we define a modified aggregation and dispersion algorithm M AGG where particles can disconnect and moves rely only on the current state.Here, particles occupy nodes of a finite region of the triangular lattice, again moving stochastically and favoring configurations with more pairs of neighboring particles.Each particle has its own Poisson clock and, when activated, chooses a random adjacent lattice node.If that node is unoccupied, the particle moves there with probability λ −n , where n is the number of current neighbors of the particle, for bias parameter λ > 0. Thus, rather than biasing particles towards nodes with more neighbors, we instead discourage moves away from nodes with more neighbors, with larger λ corresponding to a stronger ferromagnetic attraction between particles (Figure 1A).This new chain M AGG converges to the same Boltzmann distribution π(σ) ∝ λ E(σ) over particle system configurations σ as the original SOPS algorithm.Details of the proofs can be found in the Materials and Methods.
where k 0 is a scaling constant, P M C is the number of particles on the periphery of the largest component, and N is the total number of particles.This phase change is qualitatively invariant to the system's size.
Let Ω be the set of configurations with N particles within our bounded lattice region.We will use the following definition to quantify aggregation for particles that can be disconnected, capturing both the size and compactness of aggregates.Definition 1.For β > 0 and δ ∈ (0, 1/2), a configuration σ ∈ Ω is (β, δ)-aggregated if there is a subset R of lattice nodes such that: 1.At most β √ N edges have exactly one endpoint in R; 2. The density of particles in R is at least 1 − δ; and 3. The density of particles not in R is at most δ.
Here, β is a measure of how small the boundary between R and its complement R must be, measuring the compactness of the aggregated particles, and δ is a tolerance for having unoccupied nodes within the cluster R or occupied nodes outside of R. We say that a configuration is dispersed if no such (β, δ) exist.
By carefully analyzing the stationary distribution of M AGG , which is just the desired Boltzmann distribution, we establish conditions that provably yield aggregation when the particles are confined to a compact region of the triangular lattice (Figure 1B).The proof uses arguments from [29]; see the Materials and Methods for details.Theorem 2. Let configuration σ be drawn from the stationary distribution of M AGG on a bounded, compact region of the triangular lattice, when the number of particles N is sufficiently large.If λ > 5.66, then with high probability there exist β > 0 and 0 < δ < 1/2 such that σ will be (β, δ)-aggregated.However, when 0.98 < λ < 1.02, the configuration σ will be dispersed with high probability.
Varying values of λ in simulation gives strong indication that dispersion persists for larger values of λ and the aggregation algorithm undergoes a phase transition whereby the macroscopic behavior of the system suddenly changes from dispersion to aggregation (Figure 1C-D, Movie S1), mimicking the fixed magnetization ferromagnetic Ising model which motivated our Markov chain algorithm.Nonetheless, our proofs demonstrate that our system has two distinct phases of behavior for different ranges of λ for any system with a sufficiently large number of interacting particles, which is enough for our purposes.

BOBbots: a model active cohesive granular matter system
Next, to test whether the lattice-based equilibrium system can be used to control a real-world swarm in which there are no guarantees of detailed balance or Boltzmann distributions, we introduce a collective of active cohesive granular robots which we name BOBbots (Figure 2A-C, S1) -Behaving, Organizing, Buzzing robots -whose design physically embodies the aggregation algorithm.Driven granular media provide a useful soft matter system to integrate features of the physical world into the toolkit for programming collectives.This builds upon three decades of work understanding how forced collections of simple particles interacting locally can lead to remarkably complex and diverse phenomena, not only mimicking solids, fluids, and gasses [30,31] -e.g., in pattern formation [32,33], supercooled and glassy phenomena [34,35], and shock waves [36] -but also displaying phenomena characteristic of soft matter systems such as stress chains [37] and jamming transitions [38,39].While cohesive granular materials are typically generated in situations where particles are small (powders, with interactions dominated by electrostatic or even van der Waals interactions) or wet (with interactions dominated by formation of liquid bridges between particles) [40,41], we generate our cohesive granular robots using loose magnets which can rotate to always achieve attraction.
The movement and interactions between BOBbots were designed to capture the salient features of the abstract stochastic algorithm while replacing all sensing, communication, and probabilistic computation with physical morphology and interactions.Each BOBbot has a cylindrical chassis with a base of elastic "brushes" that are physically coupled to an off-center eccentric rotating mass vibration motor (ERM).The vibrations caused by the rotation of the ERM are converted into locomotion by the brushes (Figure 2C).Due to asymmetry in our construction of this propulsion mechanism, the BOBbots traverse predominantly circular trajectories [42] that are randomized through their initial conditions but -unlike the SOPS particles -are inherently deterministic with some noise and occur at a constant speed per robot distributed as v 0 = 4.8 ± 2.0 cm/s.See the Materials and Methods for further details.
Analogous to the modified transition probabilities in the aggregation algorithm that discourage particles from moving away from positions where they have many neighbors, each BOBbot has loose magnets housed in shells around its periphery that always reorient to be attractive to nearby BOBbots (Figure 2C).The probability that a BOBbot detaches from its neighbors is negatively correlated with the attractive force from the number of engaged magnets, approximating the movement probabilities given by the algorithm which scale inversely and geometrically with the number of neighbors.We subsequently verify this assertion experimentally (see Section S5 of the Supplementary Materials for details).The strength of the magnets F M 0 determines whether the system aggregates or disperses in the long run, analogous to λ in the algorithm.
To allow for study of larger BOBbot ensembles and more comprehensive sweeps of parameter space, we also performed Discrete-Element Method (DEM) simulations of the BOBbots (see Figure 2D-F and the Materials and Methods for more details).The motion of an individual BOBbot is modeled as a set of overdamped Langevin-type equations governing both its translation and rotation subject to its diffusion, drift [43], magnetic attraction, and sterical exclusion with other BOBbots.The translational drift corresponds to the speed from the equilibrium of the drive and drag forces while the rotational drift corresponds to the circular rotation.Similar methods have been used to understand macroscale phenomena emerging from collectives of microscopic elements [11] and to model particle motion in active matter [44].
Mitigating the effects of the arena's fixed boundaries in both experiments and simulations presented a significant design challenge.BOBbots can persist along the boundary or in corners, affecting system dynamics by, for example, enabling aggregates to form where they would not have otherwise or hindering multiple aggregates from integrating.
To address these issues, uniform airflow was employed to gently repel BOBbots away from the boundary and similar effects were implemented in simulation.More details about the experimental apparatus and protocol can be found in the Materials and Methods.

Clustering dynamics explained by algorithm analysis
Since the critical elements of the SOPS algorithm can be physically embodied by robots as simple as our BOBbots, to test if the SOPS model could quantitatively capture collective dynamics, we next investigated the degree to which collectives of BOBbots aggregate as a function of their peripheral magnet strength F M 0 in both robotic experiments and DEM simulations.(For convenience, F M 0 is normalized by the gravity of Earth g = 9.81 m/s 2 when using the unit of gram.)The experimental protocol begins with placing magnets of a particular strength F M 0 into the BOBbots' peripheral slots.The BOBbots are positioned and oriented randomly in a rectangular arena and are then actuated uniformly for a fixed time during which the BOBbots' positions and the size of the largest connected component are tracked (Figure 3A-C).These trials are conducted for several F M 0 values with repetition.We followed the same protocol in simulations.In experiment and DEM simulation, we observe an abrupt, rapid rise and then saturation in the size N M C of the largest connected component as the magnetic attraction F M 0 increases (Figure 3D).These curves resemble those in Figure 1D, with the magnetization F M 0 playing a role analogous to the bias parameter λ.Given this correspondence, we explored whether the equilibrium SOPS model could be used to make quantifiable predictions in the robot experiments.First, we designed a test to examine how force and λ scale.Recall that in the SOPS algorithm, the force acting on each particle is proportional to λ n , where n is the particle's current number of neighbors.In the experiments, BOBbots cannot count their neighbors, but the magnets are expected to provide a similar force that also increases geometrically when more magnets are engaged.
To estimate the relationship between force and λ, we investigate the rate at which a BOBbot loses or gains neighbors over a fixed amount of time.Viewing a BOBbot's completion of half its circular motion as analogous to a particle moving to a new lattice node in the SOPS algorithm and using this time interval to evaluate the transition, simulation data shows that a BOBbot's transition probability from having a higher number of neighbors n to a lower number n closely follows the algorithm's P (σ, τ ) ∝ min(1, λ n −n ) transition probabilities (Figure 4A, S10).Further, we evaluated the BOBbots' effective bias parameter λ eff as a function of F M 0 and found an exponential relation λ eff = exp(βF M 0 ), where β is a constant representing inverse temperature (Figure 4B).The BOBbots' transition probabilities can then be approximated as , where β is the inverse temperature of the system and n = n • F M 0 can be interpreted as the energy contributed by a BOBbot's n neighbors.With the relation between F M 0 and λ eff established, we next compare the aggregation behaviors exhibited by the SOPS algorithm and the BOBbot ensembles.Figure 4C shows the fraction of particles/BOBbots in the largest component N M C /N observed in both the SOPS algorithm and BOBbot simulations after converting with respect to λ eff ; the algorithm does indeed capture the maximum cluster fraction observed in the simulations.Notably, the aggregated and dispersed regimes in λ-space established by Theorem 2 provide a rigorous understanding of these BOBbot collective behaviors.For instance, the proven dispersed regime 0.98 < λ < 1.02 gives a clear explanation for why agents will not aggregate even in the presence of mutual attraction.Further, it also helps establish the magnitude of attraction needed to saturate the aggregation.We additionally test the SOPS prediction that the maximum cluster should not only be large but also compact, occupying a densely packed region.The results from [29] that we apply here for aggregation suggest the following relationship between the size of the largest component N M C and its perimeter P M C .In dispersed configurations, P M C should scale linearly with N M C , meaning that most BOBbots lie on the periphery of their components.In aggregated configurations, however, P M C should scale as N

A
1/2 M C , approximating the minimal perimeter for the same number of BOBbots by at most a constant factor.We test these scaling relationships in simulations with 400 BOBbots (Figure 5A) and find that the theory's predictions hold in the dispersed regime; however, the 0.66±0.07sublinear scaling power for the aggregated case is slightly higher than the theory's prediction of 0.5.This discrepancy may in part be due to boundary and finite-size effects -in fact, DEM simulations with periodic boundaries show a scaling power of 0.59 ± 0.18 that is closer to the SOPS theory (Figure S17) -but is also affected by non-reversibility inherent in the BOBbots' circular trajectories.To make a quantitative comparison that captures when components are both large and compact, we track where k 0 is a scaling constant defined such that AGG M C = 1 when the system is optimally aggregated, achieving the minimum possible perimeter.Physically, AGG M C is reminiscent of surface tension for which energy minimization leads to a smaller interface (in our setting, smaller perimeter P M C ), yielding an AGG M C closer to 1.We obtain agreement between the SOPS and DEM simulations with respect to this metric as well (Figure 4D), further validating the theory's prediction, though the DEM simulations yield slightly smaller AGG M C than the SOPS algorithm for large λ.
We noticed the size of the largest component N M C grows roughly proportional to t 1/2 over time (Figure 3C).Since the perimeter of the largest cluster P M C scales proportional to N 0.66 5A), this implies the length scale grows like t 1/3 .This is reminiscent of coarsening in a broad class of systems described by Cahn-Hilliard equation ∂u/∂t = ∇ 2 (Φ (u) − γ∇ 2 u) where order parameter u takes continuous values in (−1, 1) where −1 and 1 are analogous to empty and occupied nodes in the SOPS lattice, respectively.To bridge the SOPS algorithm with the Cahn-Hilliard equation, we first observe that the SOPS algorithm with bias parameter λ can be exactly mapped to an Ising model with fixed magnetization [45,46] with coupling strength J = 1 2β log λ, where β is inverse temperature (see Section S7 of the Supplementary Materials for details).As shown by Penrose [47], the fixed magnetization Ising model with coupling strength J can be mapped to the surface tension γ of the Cahn-Hilliard equation as γ = βJ.Thus, the SOPS and BOBbot ensemble behaviors map to the Cahn-Hilliard equation with γ = 1 2 log λ ∝ F M 0 .This suggests that in the limit, the SOPS and BOBbot aggregation behavior should display a second-order phase transition at a critical λ c corresponding to the critical surface tension γ in the Cahn-Hilliard equation.The corresponding critical value λ c = e 2/7 ≈ 1.33 on the hexagonal lattice lies within the λ c ∈ (1.02, 5.66) range proven by the SOPS theory (see Section S8 of the Supplementary Materials for details).Thus, we obtain agreement between the SOPS theory for a finite lattice system and the Cahn-Hilliard equation for an active matter system at the continuum limit.This mapping gives further confirmation of the universality of our results and provides another perspective for "programming" active collectives.

Enhancing clustering via local stress sensing
We have demonstrated that the BOBbot ensembles mimic a lattice model that can provably aggregate for large enough λ, corresponding physically to highly attractive interaction that favors large components with small perimeter.We now ask whether we can achieve rudimentary collective intelligence determining, for example, how robots could tune their responses to enhance or dampen aggregation, thereby achieving a more tightly clustered or dispersed state.In particular, we explore whether such tuning can help counteract some ways the system deviates from the theory, such as variations in the BOBbots' speeds and magnetic attraction, improving the fidelity to the original algorithm.While the BOBbots remain unable to count neighbors or estimate the Gibbs probabilities directly as prescribed by the algorithm, we take advantage of physical effects of the BOBbot ensembles to "program" desirable behavior without using any traditional computation.
The first effect relies on observations that for a fixed magnet strength, the size of the largest component N M C decreases with increasing BOBbot speed v 0 (Figure S9); a full investigation of the behavior of BOBbot collectives at varying uniform speeds will be the subject of a separate study.We further observe that N M C scales linearly with z, the average number of neighbors per BOBbot at equilibrium (Figure 6A, inset).Thus, BOBbot speed v 0 is inversely correlated with the average number of neighbors per BOBbot z.This arises from v 0 being a proxy for β −1 in the effective attraction λ eff .Consequently, we can mimic enhanced aggregation via increased magnet strength by reducing a BOBbot's speed as a function of its number of neighbors.
Without adapting a BOBbot's speed based on its number of neighbors, a BOBbot collective actuated uniformly at a speed v converges to an average of z std (v) neighbors per BOBbot at equilibrium (Figure 6A, red); any point in speedneighbor space deviating from z std (v) is transient.To enhance aggregation, we engineer reduced speeds v eng (z) that a BOBbot with z neighbors should adapt to (Figure 6A, blue).These slowed speeds allow the collective to reconverge to a new steady-state with a larger number of average neighbors per BOBbot (Figure 6A, arrows).This feedback between the engineered speeds v eng and the steady-state average number of neighbors z std iterates until reaching the fixed point in speed-neighbor space where the steady-state and engineered behaviors meet as z = z std (v eng (z)).The distribution of a BOBbot's number of contacts over six 10-minute experiments using F M 0 = 3 g.Each sample is an average of number of contacts over 1 second.Inset: Simulation results using the same conditions as the experiment.(B) A simulation demonstrating enhanced aggregation in an ensemble of 400 BOBbots using a weak magnet strength of F M 0 = 7 g.Each BOBbot's speed decreases from 6 cm/s to 1.2 cm/s as its stress s 0 = j∈neighbors s j /F M 0 ≈ z increases from 0 to 6, where z is its current number of neighbors.BOBbots in an aggregate's interior experience the most stress (dark gray) and thus have the slowest speeds, enabling larger aggregates to form.Without adapting speed in response to stress, the cluster sizes remain the same magnitude as in the 0 minute snapshot (left).
While adapting speeds based on numbers of neighbors would be relatively straightforward to implement in more complex robots capable of counting neighbors (e.g., optically as in [15,16,48,49]), implementing such a scheme in the deliberately simple BOBbots is challenging given their lack of such sensing.Here we utilize a second physical effect: inspired by the correlation of particle density and stress on individual particles in granular systems [50], we propose that monitoring local contact stress can function as a proxy for counting numbers of neighbors.An immediate benefit of such a scheme is that it can be implemented on the existing robots via custom, low-cost, analog surface stress sensors (see Figure 6B and the Materials and Methods for details).The implemented stress sensors function such that for sufficiently large stress (e.g., when in a cluster), motor speed is decreased by 70% (Figure 6C).
We implemented this "physical algorithm" on BOBbot ensembles with weakly attractive magnets (Movie S6).In experiments with ensembles of 10 BOBbots in a circular arena, adapting BOBbot speeds in response to stress sensing significantly increases the average number of neighbors per BOBbot (Figure 7A).Further, there is a quantitative match in the final average number of neighbors per BOBbot between the experiments and the fixed points predicted in Figure 6A, validating our control strategy for enhancing aggregation.Simulations using the same arena and stress-mediated response reproduce the experimental results (Figure 7A, inset).In simulations of 400 BOBbots with F M 0 = 7 g, we observe that BOBbots with more neighbors experience higher stress and thus have the slower speeds (Figure 7B).This stress-mediated decrease in speed enables large aggregates to form that would not have existed otherwise in the weakly attractive regime.The use of stress sensing opens an interesting avenue for collectives of rudimentary robots to incorporate higher-order information without complex vision systems; further, contact stress provides insights (e.g., closeness to a jamming transition) that could be valuable in densely packed clusters [51].

Object transport in the aggregated phase
Encouraged by the close connections between the physical system and the underlying theoretical model along with the successful control scheme for enhanced aggregation using stress sensing, we sought to test whether aggregated BOBbots could collectively accomplish a task.In particular, could an aggregated BOBbot collective "recognize" the presence of a non-robot impurity in its environment and cooperatively expel it from the system?Typically, such collective transport tasks -e.g., the cooperative transport of food by ants [52,53] -either manifest from an orderdisorder transition or rely heavily on conformism between agents for concerted effort and alignment of forces.With our BOBbot collectives, we instead aim to accomplish transport via simple mechanics and physical interactions emergently controlling global behavior without any complex control, communication, or computation.
By maintaining a high magnetic attraction F M 0 , we remain in the aggregated regime where most BOBbots connect physically and can cumulatively push against untethered impurities (e.g., a box or disk) introduced in the system (Figure 8A, Movie S7).The BOBbot collective's constant stochastic reconfiguration grants it the ability to envelop, grasp, and dislodge impurities as their individual forces additively overcome the impurities' friction, leading to large displacement in the aggregated regime (Figure 8B, right) with a median displacement of 7.9 cm over 12 minutes.On the contrary, we find that systems with weak magnetic attraction (i.e., those in the dispersed regime) can typically only achieve small impurity displacement (Figure 8B, left) with a median displacement of 0.9 cm over 12 minutes (see Figure S11 for distributions).We observe infrequent anomalies in which dispersed collectives achieve larger displacement than aggregated ones, but these outliers arise from idiosyncrasies of our rudimentary robots (e.g., an aggregated cluster of BOBbots may continuously rotate in place without coming in contact with an impurity due to the BOBbots' individual orientations in the aggregate; see Movie S7).
Characterizing the impurity's transport dynamics as mean-squared displacement over time r 2 (τ ) = vτ α reveals further disparities between the aggregated and dispsered BOBbot collectives (Figure 9A).On a log-log plot, the intercept indicates log(v), where v is the characteristic speed of the impurity's transport; we observe that in all but one fringe case the strongly attractive collectives achieve transport that is orders of magnitude faster than those of the weakly attractive ones (Figure 9B).The slope of each trajectory indicates the exponent α that characterizes transport as subdiffusive (α < 1), diffusive (α = 1), or superdiffusive (α > 1).While all the strongly attractive collectives immediately achieve nearly ballistic transport (with α = 1.85 ± 0.11 for τ < 20 s) indicating rapid onset of cluster formation and pushing, the weakly attractive collectives initially exhibit mostly subdiffusive transport (with α = 0.89 ± 0.56 for τ < 20 s) caused by intermittent collisions from the dispersed BOBbots (Figure 9C).When the slight heterogeneous distribution of the dispersed BOBbots remains unchanged for a sufficiently long time, the accumulation of displacement in a persistent direction can cause a small drift, leading to ballistic transport at a longer time scale.These results align with the predictions of a simple model combining subdiffusive motion with small drift (Figure S12).Nonetheless, the transport speeds achieved by the dispersed collectives are two orders of magnitude smaller than those of the strongly attractive ones.Simulations of impurity transport (see S6 for details) reproduce the experimental results (Figure 9B, inset, Movie S7), including the rare anomalies.Seven of the 100 simulations of weakly attractive collectives succeeded in transporting the impurity to the arena boundary at slow speeds while 76 of the 100 simulations of strongly attractive collectives did so ballistically.The remaining 24 simulations of attractive collectives that did not achieve ballistic transport consistently formed an aggregate that never came into contact with the impurity.We found that disaggregating established aggregates by introducing time periods with no attraction enabled them to dissolve and reform for another attempt at transport.Using different disaggregating sequences, the attractive collectives achieved ballistic transport in 15-20% more simulations than without disaggregating (Figure S13).Physically and interestingly, in the Cahn-Hilliard picture, impurity transport can be interpreted as the expulsion of an obstacle in a continuum mixture with sufficiently high surface tension to yield phase separation.If the obstacle occupies a position that is later occupied by the solid phase, the obstacle is expelled due to sterical exclusion; when its position is unvisited by the solid phase during the process of coarsening, however, it remains stagnant, similar to the anomalies for attractive collectives.In this interpretation, disaggregating effectively repeats the coarsening process to that the probability any given position is unvisited by the solid phase is significantly diminished.

Discussion and Conclusion
In this paper, we use mathematical ideas from distributed computing and statistical physics to create task-oriented cohesive granular media composed of simple interacting robots called BOBbots.As predicted by the theory, the BOBbots aggregate compactly with stronger magnets (corresponding to large bias parameter λ) and disperse with weaker magnets (or small λ).Simulations capturing the physics governing the BOBbots' motions and interactions further confirm the predicted phase change with larger numbers of BOBbots.The collective transport task then demonstrates the utility of the aggregation algorithm.
There are several noteworthy aspects of these findings.First, the theoretical framework of the underlying SOPS model can be generalized to allow many types of relaxations to its assumptions, provided its dynamics remain reversible and model a system at thermal equilibrium.For example, noting that the probability that a robot with n neighbors detaches may not scale precisely as λ −n as suggested by the Boltzmann weights, we can generalize the SOPS model to be more sensitive to small variations in these weights: the proofs establishing the two distinct phases can be shown to extend to this setting, provided the probabilities p n of detaching from n neighbors satisfy The robustness of the local, stochastic algorithms makes the macro-scale behavior of the collective resistant to many types of idiosyncrasies inherent in the BOBbots, including bias in the directions of their movements, the continuous nature of their trajectories, and nonuniformity in their speeds and magnet strengths.Moreover, our algorithms are inherently self-stabilizing due to their memoryless, stateless nature, always converging to a desired system configuration -overcoming faults and other perturbations in the system -without the need for external intervention.In our context, the algorithm will naturally continue to aggregate, even as some robots may fail or the environment is perturbed.
We find agreement not only between the BOBbot ensembles and the discrete SOPS model, but also with continuum models of active matter.The SOPS algorithm for aggregation and dispersion was initially defined as a distributed, stochastic implementation of a fixed magnetization Ising model.In addition to showing that our experimental system follows guarantees established by the analysis of a discrete model, we also observe that the growth of its largest component matches the power-law derived for the Cahn-Hilliard equation, a continuous analog of the Ising model [47].This mapping provides an intuitive understanding of how the SOPS bias parameter λ, the physical inter-BOBbot attraction F M 0 , and the surface tension γ in the Cahn-Hilliard equation correspond; thus, as γ controls the phase change in the Cahn-Hilliard equation, so do λ and F M 0 in their respective settings.This observation buttresses our confidence that the SOPS model provides a useful algorithmic framework capable of producing valid statistical guarantees for ensembles of interacting robots in continuous space.
Moreover, we find that the nonequilibrium dynamics of the BOBbots are largely captured by the theoretical models that we analyze at thermal equilibrium, which is in agreement with the findings of Stenhammar et al. [54].For example, in addition to visually observing the phase change as the magnetic strengths increase, we are able to test precise predictions about the size and perimeter of the largest connected components based on the formal definitions of aggregation and dispersion from the SOPS model.We additionally use simulations to study the transition probability of a BOBbot from having n neighbors to having n neighbors to see if the magnetic interactions conform to the theory, and indeed we see a geometric relation decrease in the probability of moving as we increase the number of neighbors, as predicted.The resultant correspondence between the magnetic attraction and effective bias in the algorithm confirms a quantitative connection between the physical world and the abstract algorithm.
In summary, the framework presented here using provable distributed, stochastic algorithms to inspire the design of robust, simple systems of robots with limited computational capabilities seems quite general.It also allows one to leverage the extensive amount of work on distributed and stochastic algorithms, and equilibrium models and proofs in guiding the tasks of inherently out of equilibrium robot swarms.Preliminary results show that we likely can achieve other basic tasks such as alignment, separation (or speciation), and flocking through a similar principled approach.We note that exploiting physical embodiment with minimal computation seems a critical step in scaling collective behavior to encompass many cutting edge settings, including micro-sized devices that can be used in medical applications and cheap, scalable devices for space and terrestrial exploration.Additionally, we plan to further study the important interplay between equilibrium and nonequilibrium dynamics to better solidify these connections and to understand which relaxations remain in the same universality classes.

Details of the SOPS algorithm and proofs
The SOPS algorithm M AGG for aggregation and dispersion is given in Algorithm 1.The algorithm is presented as a Markov chain, but could easily be modified to function as a distributed algorithm executed by each particle independently and concurrently as shown in [26,29].
Algorithm 1 Markov chain M AGG for aggregation and dispersion in SOPS Beginning at any configuration of N particles in a bounded region, fix λ > 1 and repeat: 1: Choose a particle P uniformly at random; let be the lattice node it occupies.2: Choose an adjacent lattice node and q ∈ (0, 1) each uniformly at random.3: if is empty and q < λ −n , where n is the number of neighbors P has at then Recall that Theorem 2 analyzes the stationary distribution π of the Markov chain M AGG for aggregation and dispersion.In particular, Theorem 2 was shown in [29] to hold for π(σ) ∝ λ −b(σ) ∝ λ E(σ) , where b(σ) is the number of "boundary edges" of the lattice that have exactly one endpoint occupied by a particle.So it remains to show that M AGG converges to this stationary distribution π.

Lemma 3. The unique stationary distribution of
Proof.Let σ and τ be any two SOPS configurations with σ = τ such that Pr(σ, τ ) > 0, implying that τ can be reached from σ by a single move of some particle P .Suppose P has n neighbors in σ and has n in τ .We must show the detailed balance condition holds with respect to the transition probabilities: The algorithms in [26,29] were designed using the Metropolis-Hastings algorithm [28] which specifies transition probabilities Pr(σ, τ ) = min{π(τ )/π(σ), 1} to capture the ratio between stationary weights of the current and proposed configurations.So we have that π(τ )/π(σ) = λ n −n .It is then easy to see that this ratio is unchanged by the modified transition probabilities where Pr(σ, τ ) = λ −n and Pr(τ, σ) = λ −n , and thus detailed balance is satisfied: Therefore, since π satisfies detailed balance and M AGG is an ergodic finite Markov chain, we conclude that π is the unique stationary distribution of M AGG .
We conclude by outlining the proof of Theorem 2 that shows M AGG achieves aggregation when λ is large enough and dispersion when λ is close to one.Our proof is a series of information-theoretic arguments about the stationary distribution π.We use ideas similar to Peierls arguments, which are often used in statistical physics to study phase changes in behavior space for infinite systems [55].In [29] it was shown that, for finite systems, particles of two different colors could either separate into monochromatic clusters or integrate, indifferent to color.This separation algorithm can be applied to the setting where a bounded region of the lattice is completely filled with particles that move by "swapping" places with their neighbors.By viewing particles of one color as "empty space" and particles of the other color as our particles of interest, the swap moves in the separation algorithm correspond to particle moves within a bounded area.These are precisely the moves used in our aggregation algorithm, where separation corresponds to aggregation and integration corresponds to dispersion.Thus, it is straightforward to leverage the arguments for separation and integration in [29] to show aggregation and dispersion in a bounded region.
For large enough bias λ, we prove aggregation occurs with high probability as follows.Using techniques introduced in [56], we define a map from any configuration without an aggregate to a configuration with an aggregate by (i) choosing some scattered particles in a systematic way and (ii) rearranging them as an aggregate in a carefully chosen location.We then show that no aggregate configuration has too many preimages under this map because of the careful way we remove scattered particles.On the other hand, we show that applying this map to a dispersed configuration leads to a large increase in its stationary probability.Provided λ is large enough that the probability gain outweighs the number of preimages, these two facts imply that aggregated configurations are much more likely to occur in the stationary distribution than dispersed ones.More formally, the above argument shows that the stationary probability of being in a dispersed configuration is at most (c 1 /λ) c2 √ N , where c 1 , c 2 > 0 are constants that depend on the map described above.Thus, provided λ is large enough, this probability of being in a dispersed configuration is very small, proving that aggregation is achieved with high probability.
When the bias λ is close to one, we can prove that dispersion occurs with high probability.We show that there exist polynomially many events such that if aggregation occurs, then at least one of these events must also occur.These events correspond to certain regularly-shaped subregions of the lattice being almost entirely occupied by particles.We then use a Chernoff-type bound to show that each of these events is exponentially unlikely when λ is close to one.This implies that the stationary probability for aggregated configurations is at most the sum of polynomially many terms that are each exponentially small, so dispersion must occur with high probability for this range of λ.

BOBbot design
The BOBbot mechanical design was developed in SolidWorks, and its skeleton was 3D printed in ABS plastic by a Stratsys UPrint SE Plus printer at a layer resolution of 0.010 inches and sparse density (Figure S1).Each BOBbot contains a lithium ion polymer battery (Adafruit Industries) that is equipped with Qi wireless charging for recharging between experiments (Adafruit Industries).The brushbot design is implemented using an ERM (BestTong) for vibrations and two Pienoy dog toothbrush heads as feet, yielding noisy circular trajectories (Movie S2).The BOBbot's motor circuitry was assembled on a Printed Circuit Board (PCB) designed in EagleCAD (Figure S2).The PCBs were printed at the Georgia Tech Interdisciplinary Design Commons makerspace and outsourced from JLCPCB.This circuitry is switched and modulated by a phototransistor (Adafruit Industries), which acts as a proportional controller for motor speed.Grade N42 neodymium magnets (K&J Magnetics) are housed in the BOBbot chassis for inter-robot attraction, and can be swapped for magnets of different strengths to modulate the BOBbots' cohesion.A complete list of BOBbot components can be found in Table S1.
To achieve stress sensing, each BOBbot is equipped with four triggers that mechanically deform and close the circuit upon collisions to sense the locally exerted stress (Figure 6B, S2).These triggers are positioned radially in front of the permanent magnets in the chassis.The stress sensors function such that a robot decreases its motor speed for sufficiently large stress (Figure 6C).The analog circuit is designed to reduce the motor's current in a manner proportional to the total number of contacts, starting with roughly 70% reduction for a single triggered sensor (Figure 6C, top).When multiple sensors are triggered, a BOBbot's speed is practically negligible.

Simulations
To simulate the SOPS, we execute the algorithm on a hexagonal lattice.The size of the lattice in Figure 1 is chosen to be sufficiently large so that boundary effects are mitigated.The size of the lattice for Figure 4 is chosen to match the area density and the number of agents in the physical evolution and algorithm.To determine the constant k 0 in the aggregation metric 2 • 6 and perimeter P M C = 6 , setting k 0 so that AGG M C = 1.This yields k 0 = 1/ 8 √ 3.
Beyond the information described in the main text, the DEM simulations faithfully represent the spherical loose magnets with exponentially decaying force housed in each BOBbot's chassis slots, resulting in patchy magnetic interaction as the magnets move freely in their slots.Attraction between two simulated BOBbots is calculated based on these magnetic spheres' strength and the minimum physical separation between any interacting pair, which depends on the relative position and orientation of the two BOBbots.
To calibrate our DEM simulations, we measure the BOBbots' physical parameters and use these values for the simulated BOBbots (Table 1).Most parameters such as the mass and dimensions of each BOBbot are directly measured.
For others, we use a series of experiments designed to isolate individual parameters.For instance, to avoid possible system errors such as in-plane friction when measuring the magnetic force, we measured the minimum force needed to overwhelm the magnetic force in vertical direction (Figure S4).Other indirect measurements involve the translational and rotational drag (Figure S5, S6).The key ingredient in these experiments is to use a known force (Earth's gravity) to calibrate these intricate forces.Details can be found in Section S2 of the Supplementary Materials.
The DEM simulations use the Euler-Maruyama method with a time step of 1 ms to integrate the following Newton equations: As the agents are in the overdamped regime where |m ¨ r| η| ˙ r|, the Newton equations are equivalent to the Langevin equations for active Brownian particles by taking the limit m, I → 0. ˙ r = v 0 û + F env ( r, ϕ)/η + ξ(t)/η φ = ω 0 + τ env ( r, ϕ)/η ϕ + ξ ϕ (t)/η ϕ As we see from the reduced equations, in the steady state, a BOBbot will perform a circular motion with a saturated speed v 0 = F D /η and a frequency of ω 0 = τ D /η ϕ .This suggests that we can control a BOBbot's speed v 0 by changing its motor vibration strength, varying F D .
The initial placement of the BOBbots is achieved by greedy rejection sampling, sequentially placing BOBbots in random positions that do not overlap with the previously placed BOBbots.A cell list search method is used to speed up the simulation's computation by subdividing the simulated arena into square cells so that, when integrating forces for a given BOBbot, only consider interactions with BOBbots from the same or adjacent cells.The size of the cells is chosen such that the relative error caused by this approximation is within 10 −3 .and theoretical perspectives.Data and Materials Availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials.Additional data related to this paper may be requested from the authors.translational drag coefficients; the one that produces the trajectory most closely matching those in the experiments is chosen as the simulation η (Figure S5).

Supplementary Materials
The measurement of the rotational drag η ϕ exploits its balance with the driving torque.To measure the rotational torque exerted on a BOBbot, a very light rigid straw is attached across the diameter of a BOBbot (Figure S6).We then let the BOBbot use the straw to push objects at various arm lengths.For a given obstacle to push, the rotational torque is obtained by finding the largest torque of friction on an obstacle to balance.We decrease the arm length from a large value to a point the BOBbots can just push the obstacle.Given the measured saturated angular velocity ω 0 , the rotational drag can be inferred as η ϕ = τ D /ω 0 .

S3. Sensitivity to initial conditions
The BOBbots move in deterministic but noisy circular trajectories.Although the noise only causes small deviations in the individual trajectories, the combination of all the forces in the ensemble makes the system behavior very sensitive to initial conditions.Simulations with the same random seed started from almost exactly the same initial conditions except for a discrepancy of 0.5 mm of one robot's initial position deviate from each other significantly after 1 minute (Figure S8).We therefore regard this system as ergodic, yielding reasonable statistical sampling on longer time scales.

A B
0.5 mm

S4. Dependence of maximum cluster size on BOBbot speed and curvature
To investigate the effect of the BOBbots' individual speeds on the size of the maximum cluster, we run simulations for F M 0 = 3, 7, 10 g with 20 repetitions for 8 speeds equally spaced in range v 0 = 1-8 cm/s.Figure S9 shows how the maximum cluster size N M C decreases as the BOBbots' individual speed is increased.We also find that N M C increases with larger radii of curvature corresponding to decreased torque.

S5. Detachment probability
In Figure 4A we showed the probability of a BOBbot detaching from three neighbors.Figure S10 shows the complete set of detachment probabilities for various magnetic attraction strengths F M 0 .We find that the λ eff increases with the magnetic attraction.For a particular attraction strength, λ eff only varies within a small range such that the fitted exponential lines are nearly parallel to each other.

S6. Object transport by BOBbot collectives
Figure S11 shows the displacement of the box impurity over time for weakly attractive and strongly attractive BOBbot collectives.While nearly all experimental runs with strongly attractive collectives exhibit rapid and large displacement, some of the weakly attractive collectives exhibit two-stage transport dynamics that start with very little displacement and eventually achieve larger displacement.We posit that the weakly attractive collectives' two-stage dynamics are composed of subdiffusion arising from the BOBbots' collisions and a small drift caused by the persistent heterogeneity of the BOBbots around the impurity.To validate this hypothesis, we developed a toy model in MATLAB where the subdiffusion r H (t) with mean-squared displacement |r ) is generated by the fractional Brownian motion generator (from the MATLAB Wavelet Toolbox) and is added to a drift motion r D (t) = v D t.The relative magnitude difference between the subdiffusion and drift is chosen to match the experiments.When the drift is small (i.e., |v D | = 0.005), we observe two-stage transport dynamics consistent with the experiments (Figure S12, magenta).On the other hand, when the drift is dominant over the subdiffusion as in the strongly attractive collectives (i.e., |v D | = 0.5), the toy model reproduces the nearly ballistic trajectories observed in experiment (Figure S12, blue).In fact, the mean-squared displacement of this composed motion MSD(r H + r D ) is where (r H (t + τ ) − r H (t)) t vanishes due to the isotropy of subdiffusion.The final two equations demonstrate how the subdiffusive power is dominated by the ballistic power 2 when the drift speed |v D | is large.
DEM simulations of the impurity transport task use the same obstacle parameters (e.g., rectangular shape and a mass of 60 g) as in the experiments and approximates the friction coefficient by 0.25.The friction of the box is integrated all over the contact area.The friction on each point is anti-parallel to the instantaneous velocity.Sample animations of the simulations can be found in Movie S7.
We additionally performed DEM simulations with intermittent time periods with no attraction in order to dissolve and disaggregate formed aggregates.Movie S7 demonstrates how aggregates initially missing the obstacle to transport can find it successfully after disaggregating.All three disaggregating sequences we investigated result in more successful transports to the boundary when compared to the base attractive case without disaggregating.

S7. From SOPS to the fixed-magnetization Ising model
In this section, we prove that the SOPS algorithm can be mapped to an fixed-magnetization Ising model with coupling strength J = 1 2β log λ.For a given SOPS configuration of particles in a bounded region of the lattice, construct a corresponding Ising lattice gas where the spin σ of an occupied (resp., unoccupied) node in the SOPS is +1 (resp., −1) in the gas.The SOPS algorithm has a transition probabilities P SOPS = λ −∆HSOPS , where H SOPS = −n move is its Hamiltonian and n move is the number of neighbors that the moving particle is leaving.The Ising model has transition probabilities P Ising = exp(−β∆H Ising ), where H Ising = − 1 2 J (i,j) σ i σ j is its Hamiltonian and σ i ∈ {−1, +1} is the spin of site i.Lemma 4. Let H SOPS be the Hamiltonian of the SOPS algorithm, H Ising be the Hamiltonian of a fixed-magnetization Ising model, and J be the coupling strength of the Ising model.Then for any particle move in the SOPS algorithm and the corresponding spin updates in the Ising model, we have ∆H Ising = 2J∆H SOPS .
Proof.Consider a particle moving from node i to node j in the SOPS algorithm and let n i (resp., n j ) be the number of neighbors the particle has at node i (resp., node j).It is easy to see that ∆H SOPS = n j − n i for this move.To ∆H Ising , observe that the corresponding spin changes in the Ising model are σ i : +1 → −1 and σ j : −1 → +1.Let z be the coordination number (i.e., degree) of the lattice.Consider all sites k adjacent to i and j; we have four cases: 1. k is occupied and adjacent to i, so σ i σ k : +1 → −1.There are n i such sites.
2. k is occupied and adjacent to j, so σ j σ k : −1 → +1.There are n j such sites.
3. k = j is unoccupied and adjacent to i, so σ i σ k : −1 → +1.There are z − n i − 1 such sites.
To calculate ∆H Ising , we simply sum the spin changes of these cases; all other spins remain the same and thus cancel in the difference.We have: Figure S14 shows examples of particle moves and their corresponding changes to the Ising Hamiltonian to illustrate the relationship established by Lemma 4. Lemma 4 shows that it is exactly when J = 1 2β log λ that we achieve P SOPS = P Ising , which completes the mapping from the SOPS algorithm to the fixed-magnetization Ising model.The SOPS theory predicts that strongly attractive ensembles should produce aggregates that are both large and compact, namely, that the perimeter of the largest connected component P M C should scale with its size N M C with a 1/2 power.However, data from BOBbot experiments follow a slightly larger 0.66 ± 0.07 exponent (Figure S17, blue).This discrepancy is due in part to boundary and finite-size effects, so in simulation we investigated periodic boundary conditions (Figure S17, red).These simulations have a 0.59 ± 0.18 exponent that better aligns with the predictions of the SOPS theory.M C with periodic boundary conditions.Scaling between the largest component's size N M C and perimeter P M C in number of BOBbots for simulated systems of 100-400 BOBbots with F M 0 = 19 g using fixed boundary conditions (blue) and periodic boundary conditions (red).The fixed boundary conditions achieve a scaling power of 0.66 ± 0.07 while periodic boundary conditions achieve a scaling power of 0.59 ± 0.18.

Figure 1 :
Figure 1:The theoretical selforganizing particle system (SOPS).(A) A particle moves away from a node where it has n neighbors with probability λ −n , where λ > 0. Thus, moves from locations with more neighbors are made with smaller probability than those with fewer (e.g., in the insets,p 1 = λ −3 < p 2 = λ −2 < p 3 =1).(B) Time evolution of a simulated SOPS with 1377 particles for λ = 7.5 showing progressive aggregation (Movie S1).The bulk of the largest connected component is shown in blue and its periphery is shown in light blue.(C) Time evolution of N M C , the size of the largest connected component, showing dispersion for λ = 1.5 and aggregation for λ = 12.The simulations use 400 particles.(D) Phase change in λ-space for the aggregation metric AGG M C = N M C /(k 0 P M C √ N ), where k 0 is a scaling constant, P M Cis the number of particles on the periphery of the largest component, and N is the total number of particles.This phase change is qualitatively invariant to the system's size.

Figure 2 :
Figure 2: BOBbots and their collective motion.(A) Schematic of experimental setup.BOBbots are placed in a level arena with airflow gently repelling them from the boundaries.(B) A closeup of the experimental platform.(C) Mechanics of the BOBbots.Loose magnetic beads housed in the BOBbots' peripheries can reorient so BOBbots always attract each other.The vibration of the ERM motor and the asymmetry of bristles lead to the directed motion.The light sensor activates the motion.(D) Discrete element method simulation setup.(E) The BOBbotboundary interactions: airflow repulsion f A , BOBbot-boundary friction f BW , and normal force F BW,n .(F) The inter-BOBbot interactions: attraction between magnetic beads F M , inter-BOBbot friction f BB , and sterical exclusion F BB,n .

15 BFigure 3 :
Figure 3: Evolution of BOBbot clusters.(A) Time evolution snapshots of both experiment (Movie S3) and (B) simulation (Movie S4) for a system of 30 BOBbots with different magnet strengths: F M 0 = 5 g (left) where we observe dispersion, and F M 0 = 19 g (right) where we observe aggregation.Experimental images have been processed with a low-pass filter for better visual clarity.(C) Time evolutions of the size of the largest component N M C in experiment and simulation for a system of 30 BOBbots with F M 0 = 5 g (magenta) and F M 0 = 19 g (blue).(D) Scaling of cluster size vs. magnetic strength for a system of 30 BOBbots showing an increase in N M C as the magnet strength F M 0 increases.The yellow plot line shows the mean and standard deviation of N M C in the 150 simulation runs for each magnetic strength F M 0 between 1-35 g, with a step size of 1 g.Experimental data is shown in red with error bars showing the standard deviation of largest cluster size N M C and the uncertainty of F M 0 due to empirical measurement.

Figure 4 :
Figure 4: Algorithmic interpretation of BOBbot clustering.(A) A diagram showing how the effective bias parameter λ eff is evaluated from the DEM simulation.(B) The dependence of λ eff on the magnetic attraction force F M 0 .The (C) maximum cluster fraction N M C /N and (D) aggregation metric AGG M C for different values of λ in both the SOPS algorithm (blue) and physical simulations (red).The green and blue shaded regions show the dispersed and aggregated regimes proved from theory, respectively.

Figure 5 :
Figure 5: Perimeter scaling of BOBbot clusters.(A) Log-log plot showing the scaling relationship between the largest component's size N M C and perimeter P M C in number of BOBbots for simulated systems of 400 BOBbots with F M 0 = 5 g (magenta) and 19 g (cyan) for fixed boundary conditions.Each data point is the average of 20 simulations.While the SOPS predicts a scaling power of 0.5 for the aggregated case (cyan), the data shows a slightly larger -but still sublinearpower of 0.66 ± 0.07.(B) Final snapshot of the collective motion of 400 BOBbots with F M 0 = 5 g (left) and 19 g (right).BOBbots shown in black belong to the largest connected component; those outlined in red are on its periphery.

Figure 6 :
Figure 6: Design and implementation of stress sensing for enhanced aggregation.(A) The effect of the engineered, adaptive speeds (blue) on the steady-state average number of neighbors per BOBbot (red) for F M 0 = 3 g.Without adapting speeds, BOBbots actuated at a given speed v i would obtain an average of z i neighbors per BOBbot at equilibrium (initial point i).With the adaptive speeds, an average of z i neighbors per BOBbot causes the average speed to slow (i → 1) which in turn enables convergence to the steady-state response with more neighbors per BOBbot (1 → 2).This feedback iterates until the steady-state and engineered responses coincide at final point f = (v f , z f ), where v f < v i and z f > z i .Inset: The mapping between maximum cluster size N M C and the average number of neighbors per BOBbot z indicates the stress-sensing control strategy will increase component sizes.(B) A BOBbot equipped with a stress sensor and the schematic top-view sketch of the triggered and not-triggered states.(C) The BOBbot's response to stress.Top: the speed of a BOBbot when its sensor is and is not triggered.Bottom: The rate of sensor triggering as function of the stress applied.

Figure 7 :
Figure 7: Adapting speed via stress sensing enhances aggregation.(A)The distribution of a BOBbot's number of contacts over six 10-minute experiments using F M 0 = 3 g.Each sample is an average of number of contacts over 1 second.Inset: Simulation results using the same conditions as the experiment.(B) A simulation demonstrating enhanced aggregation in an ensemble of 400 BOBbots using a weak magnet strength of F M 0 = 7 g.Each BOBbot's speed decreases from 6 cm/s to 1.2 cm/s as its stress s 0 = j∈neighbors s j /F M 0 ≈ z increases from 0 to 6, where z is its current number of neighbors.BOBbots in an aggregate's interior experience the most stress (dark gray) and thus have the slowest speeds, enabling larger aggregates to form.Without adapting speed in response to stress, the cluster sizes remain the same magnitude as in the 0 minute snapshot (left).

Figure 8 :
Figure 8: Object transport using aggregation.(A) Schematic of the experimental setup.(B) Time evolution snapshots of box transport by a system of 30 BOBbots with magnet strength F M 0 = 5 g and 19 g (Movie S7).The box has a mass of 60 g.The final panel shows the object's complete trajectory, where D denotes the Euclidean distance of the final displacement.

Figure 9 :
Figure 9: Object transport using aggregation.(A) Mean-squared displacement of the box over time in log-log scale for collectives with F M 0 = 5 g (magenta) and 19 g (blue).(B) Distribution of the average speed, calculated as the final displacement D (as shown in Fig. 8B) divided by total time.Inset: Simulation results for the overall transport speed.The two peaks for F M 0 = 19 g correspond to pushing to the edges and corners.(C) Distributions of the mean-squared displacement exponent α at short time scale τ < 20 s.

4 :P
moves to .5: else P remains at .

Figure S1 .
Figure S1.Cross-sectional views of the BOBbot mechanical design.FigureS2.BOBbot circuitry.FigureS3.Experimental platform design and details.FigureS4.Calibration experiment for calculating magnet force F M 0 .FigureS5.Calibration experiment for calculating translational drag coefficient η.FigureS6.Calibration experiment for calculating rotational drag coefficient η ϕ .FigureS7.Boundary airflow effects in experiment and simulation.FigureS9.Dependence of maximum cluster size N M C on BOBbot speed v 0 and curvature R C .FigureS10.Probability of detachment for various magnetic attraction.FigureS11.Object transport trajectories.FigureS12.Toy model for object transport.Figure S13.Transport enhanced by disaggregating.Figure S14.Examples showing ∆H Ising = 2J∆H SOPS .Figure S15.Critical surface tension γ and bias parameter λ.Figure S16.Pattern formation below and above critical λ.Figure S17.Approaching P M C ∝ N

Figure S6 .
Figure S1.Cross-sectional views of the BOBbot mechanical design.FigureS2.BOBbot circuitry.FigureS3.Experimental platform design and details.FigureS4.Calibration experiment for calculating magnet force F M 0 .FigureS5.Calibration experiment for calculating translational drag coefficient η.FigureS6.Calibration experiment for calculating rotational drag coefficient η ϕ .FigureS7.Boundary airflow effects in experiment and simulation.FigureS9.Dependence of maximum cluster size N M C on BOBbot speed v 0 and curvature R C .FigureS10.Probability of detachment for various magnetic attraction.FigureS11.Object transport trajectories.FigureS12.Toy model for object transport.Figure S13.Transport enhanced by disaggregating.Figure S14.Examples showing ∆H Ising = 2J∆H SOPS .Figure S15.Critical surface tension γ and bias parameter λ.Figure S16.Pattern formation below and above critical λ.Figure S17.Approaching P M C ∝ N

Figure S13 .
Figure S1.Cross-sectional views of the BOBbot mechanical design.FigureS2.BOBbot circuitry.FigureS3.Experimental platform design and details.FigureS4.Calibration experiment for calculating magnet force F M 0 .FigureS5.Calibration experiment for calculating translational drag coefficient η.FigureS6.Calibration experiment for calculating rotational drag coefficient η ϕ .FigureS7.Boundary airflow effects in experiment and simulation.FigureS9.Dependence of maximum cluster size N M C on BOBbot speed v 0 and curvature R C .FigureS10.Probability of detachment for various magnetic attraction.FigureS11.Object transport trajectories.FigureS12.Toy model for object transport.Figure S13.Transport enhanced by disaggregating.Figure S14.Examples showing ∆H Ising = 2J∆H SOPS .Figure S15.Critical surface tension γ and bias parameter λ.Figure S16.Pattern formation below and above critical λ.Figure S17.Approaching P M C ∝ N

Figure S16 .
Figure S1.Cross-sectional views of the BOBbot mechanical design.FigureS2.BOBbot circuitry.FigureS3.Experimental platform design and details.FigureS4.Calibration experiment for calculating magnet force F M 0 .FigureS5.Calibration experiment for calculating translational drag coefficient η.FigureS6.Calibration experiment for calculating rotational drag coefficient η ϕ .FigureS7.Boundary airflow effects in experiment and simulation.FigureS9.Dependence of maximum cluster size N M C on BOBbot speed v 0 and curvature R C .FigureS10.Probability of detachment for various magnetic attraction.FigureS11.Object transport trajectories.FigureS12.Toy model for object transport.Figure S13.Transport enhanced by disaggregating.Figure S14.Examples showing ∆H Ising = 2J∆H SOPS .Figure S15.Critical surface tension γ and bias parameter λ.Figure S16.Pattern formation below and above critical λ.Figure S17.Approaching P M C ∝ N

Figure S1 :
Figure S1.Cross-sectional views of the BOBbot mechanical design.Figure S2.BOBbot circuitry.Figure S3.Experimental platform design and details. Figure S4.Calibration experiment for calculating magnet force F M 0 .Figure S5.Calibration experiment for calculating translational drag coefficient η.Figure S6.Calibration experiment for calculating rotational drag coefficient η ϕ .Figure S7.Boundary airflow effects in experiment and simulation.Figure S9.Dependence of maximum cluster size N M C on BOBbot speed v 0 and curvature R C . Figure S10.Probability of detachment for various magnetic attraction.Figure S11.Object transport trajectories.Figure S12.Toy model for object transport.Figure S13.Transport enhanced by disaggregating.Figure S14.Examples showing ∆H Ising = 2J∆H SOPS .Figure S15.Critical surface tension γ and bias parameter λ.Figure S16.Pattern formation below and above critical λ.Figure S17.Approaching P M C ∝ N 1/2 M C with periodic boundary conditions.Table S1.List of BOBbot components.Movie S1.Aggregation dynamics in a self-organizing particle system (SOPS).Movie S2.Individual BOBbot dynamics.Movie S3.Aggregation and dispersion in BOBbot collectives.Movie S4.Aggregation and dispersion in simulated BOBbot collectives.Movie S5.BOBbot collective dynamics in large simulated systems.Movie S6.Enhanced BOBbot aggregation using mechanical stress sensing.Movie S7.Object transport by BOBbot collectives.

Figure S3 :𝑟Figure S4 :
Figure S3: Experimental platform design and details.(A) The experimental platform is composed of a T-slot base supporting the foam board, aluminum, and particle board layers.(B) Levelling screws in the T-slot framing allow for incline adjustment.(C) A leaf blower with a multi-pronged tygon tubing attachment provide airflow to the PVC pipe boundary to mitigate boundary effects.

Figure S6 :Figure S7 :
Figure S6: Calibration experiment for calculating rotational drag coefficient η ϕ .(A) The experimental setup and (B) the corresponding force diagram, where f max denotes the largest frictional torque that the driving torque τ D can balance.

Figure S8 :
FigureS8: Ensemble sensitivity to initial conditions.(A) To measure the ensemble's sensitivity to initial conditions, we simulated three runs with the same random seeds and nearly identical initial conditions, with the exception of a spatial difference of 0.5 mm for the initial position of one BOBbot, shown in blue, yellow, and red.(B) The x-position over time of the perturbed BOBbot for different starting positions, where the curve color corresponds to the starting position in (A).The trials only coincide for roughly the first 30 s, and after one minute they diverge significantly.

Figure S9 :
Figure S9: Dependence of maximum cluster size N M C on BOBbot speed v 0 and curvature R C .The statistics shows the ensemble average from 20 simulations for each data point.

Figure S10 :
Figure S10: Probability of detachment for various magnetic attraction.Colored lines show the simulation data and black lines show the exponential fits.

Figure S13 :
Figure S13: Transport enhanced by disaggregating.(A) Different patterns for magnetic strength over time used in simulations.The last three interleave periods of strong magnetic attraction with periods without any magnetic attraction (disaggregating).(B) The percentage of simulation runs achieved by weakly attractive collectives (red), strongly attractive collectives (green), and strongly attractive collectives with disaggregating (blue).The statistics use 100 simulations for each scenario.An animation of simulations using disaggregting sequence C can be found in Movie S7. (C) Comparison of transport time to the boundary among collectives that are weakly attractive, strongly attractive, and strongly attractive with disaggregating.

1 −Figure S14 :Figure S16 : 2 M
Figure S14: Examples showing ∆H Ising = 2J∆H SOPS .A particle moves to the right with a decrease of (A) 1, (B) 2, and (C) 0 neighbors, illustrated by the graphs showing the local configuration before and after the move.The numbers in black between sites i and j are the value of σ i σ j .The number of neighbors of the moving particle is shown in white.

Table 1 :
List of parameters used in physical simulations.