Cohesive self-organization of mobile microrobotic swarms

Mobile microrobots are envisioned to be useful in a wide range of high-impact applications, many of which requiring cohesive group formation to maintain self-bounded swarms in the absence of confining boundaries. Cohesive group formation relies on a balance between attractive and repulsive interactions between agents. We found that a balance of magnetic dipolar attraction and multipolar repulsion between self-assembled particle chain microrobots enable their self-organization into cohesive clusters. Self-organized microrobotic clusters translate above a solid substrate via a hydrodynamic self-propulsion mechanism. Cluster velocity increases with cluster size, resulting from collective hydrodynamic effects. Clustering is promoted by the strength of cohesive interactions and hindered by heterogeneities of individual microrobots. Scalability of cohesive interactions allows formation of larger groups, whose internal spatiotemporal organization undergoes a transition from solid-like ordering to liquid-like behavior with increasing cluster size. Our work elucidates the dynamics of clustering under cohesive interactions, and presents an approach for addressing operation of microrobots as localized teams.

Untethered mobile microrobotic systems are envisioned to revolutionize our ability to manipulate the microscopic world with unprecedented flexibility. Mobile microrobots, actuated with external magnetic and acoustic fields, light, chemical reactions, and biological propellers, have been developed recently for precision micromanipulation, minimally invasive medical operations, and environmental applications [1][2][3][4][5][6][7][8] . However, translation of microrobots to most of these applications requires large numbers of microrobots (with sizes smaller than 50 µm) to work together to manipulate their macroscopic targets in much larger millimeter and centimeter scales.
Microrobot swarms have been introduced to address the need for collective functions and navigation of large numbers of microrobots in complex environments. Microrobot swarms of magnetic micro/nanoparticles have been utilized for enhancing functional output, and their reconfigurability improved multifunctionality and adaptability to dynamic environments [9][10][11][12][13] . For enhanced mobility and imaging, previous approaches relied on attractive interactions to keep the particles aggregated. However, in the absence of balancing attractive and repulsive interactions, these swarms are limited to disordered formations, lacking cohesive self-organization inherent in natural swarms [14][15][16] . Cohesive self-organization facilitates group formation in open spaces via attraction at large separations, and prevents jamming, overcrowding, and clumping at high densities via repulsion at small separations between agents 17 .
To emulate the bio-inspired self-organized cohesion in synthetic swarms, distancedependent interactions between microrobotic agents need to be controlled via physical forces. Here, we demonstrate self-organization of cohesive microrobotic teams emerging from their magnetic multipolar interactions in a liquid environment (Fig. 1). Formation, propulsion, and interactions of microrobots were controlled by a global precessing magnetic field. Microrobots were formed by the dynamic self-assembly of paramagnetic microparticles into anisotropic linear chains. The interactions between the microrobotic chains were controlled via their induced magnetic moments.
At specific opening angles of the precessing field, magnetic interactions between chains carry a slowly decaying dipolar attraction and a rapidly decaying multipolar repulsion. These counteracting effects give rise to a steady-state distance between pairs, where the sum of dipolar and multipolar forces equates to zero. Under cohesive interactions, chains self-organize into clusters by arranging themselves at steady-state distances from their neighbors. These clusters leverage collective synergies to move faster as cluster size increases, which is promoted by hydrodynamic interactions. Group formation and dissolution was mainly determined by the competition between cohesive interactions and intragroup heterogeneities. Internal organization of clusters ranged from solid-like ordering to liquid-like dynamic behavior, which depended on group size. Our approach addresses the operation of microrobots as localized teams, which will inspire researchers interested in active matter and microrobotics applications for designing advanced collective systems.

Cohesive interactions of chain microrobots
Microrobots were formed by the dynamic self-assembly of superparamagnetic particles into linear chains under a global precessing magnetic field. Initially, a low concentration of monodisperse superparamagnetic particles (particle radius a was approximately 2.5 μm) in deionized water were dispersed in a microchannel and sedimented on the planar glass substrate. The particles were actuated with a precessing magnetic field (B0 = 10 mT, Ω = 18.8 rad/s, 3 Hz) (Supplementary Note 1). The precessing magnetic field is defined by two parameters: the precession angle Ψ, which is the angle between the instantaneous magnetic field vector B(t) and the axis of precession, and the tilt angle ϑ, which is the angle between the precession axis and the normal vector to the substrate plane (Fig. 1a). Under the magnetic field, particles interact with their induced magnetic dipoles and form chains by head-to-toe alignment of their induced dipoles. Formation of chains is a dynamic self-assembly process, which depends on the magnetic field strength, the field frequency, the fluid viscosity, and the size and magnetic susceptibility of particles 8 . Once assembled, chains synchronously rotate with B(t). By tilting the precession axis by ϑ, chains self-propel over the substrate, orthogonal to the direction of the tilt. Self-propulsion results from the hydrodynamic symmetry breaking mechanism due to the rotation of chains near the solid boundary, which has been described in detail elsewhere 8,18,19 .
We first investigated the pairwise dynamics of a homogeneous pair of chains for varying field parameters: 60 • ≤ Ψ ≤ 76 • at ϑ = 5 • (Figs. 1b, S1, Video S1). For Ψ < 65 • , the trajectories of the pairs diverged, during which the distance between the pairs steadily grew until there was not any discernible pairwise interaction (Figs. 1b, S1b). For 65 • ≤ Ψ ≤ 72 • , the pairwise distance converged to and oscillated around a steady-state value that persisted during the experiment duration. In this regime, pairs translated and rotated around their center, which led to the trochoidal trajectories ( Fig. 1b, S1b). Increasing Ψ further, the steady-state distance between the chains decreased until the pairs collapsed at Ψ = 74 • (Fig. 1b, S1b). Increasing the number of chains, chains self-organized into a cluster with a discernible order in their spatial organization, where each chain was distanced at an approximately equal steady-state distance from their nearest neighbors (Fig. 1c, Video S1). Clusters could be steered by changing the orientation of the tilt axis in the xy plane without changing Ψ and ϑ (Fig. 1c).

Pairwise magnetic and hydrodynamic interactions
Pairwise steady-state distances obtained from experiments are reported for different Ψ and categorized under "divergence", "cohesion", and "collision" states based on their pairwise behavior ( Fig. 2a). We attribute the emergence of a pairwise steady-state distance to dipolar and multipolar magnetic interactions between chains caused by their anisotropic shape magnetization. The magnetic interaction forces between chains were calculated numerically for all combinations of Ψ and ϑ, and for different number of microparticles per chain, n (Figs. 2a-b, S2, see Supplementary Note 2). Particles that have isotropic magnetizations, such as spherical paramagnetic particles (n = 1), interact solely with their dipoles. For ϑ = 0 • , time-averaged dipolar interactions are repulsive when Ψ is smaller than the magic angle 54.7º 20,21 , and attractive for Ψ > 54.7º. In contrast, chains (n > 1) additionally have higher-order magnetic moments due to spatially separated dipoles positioned at centers of particles that form the chains 20,22,23 . For doublets (n = 2), the dominant repulsive multipolar interaction is the hexapole-dipole for Ψ < 61.5º, which decays rapidly with r -6 compared to the r -4 decay rate of the attractive dipolar interaction forces 20 . Similar multipolar interactions are also in play for chains with n > 2 ( Fig. S2e-g). Due to the different decay rates of attractive dipolar and repulsive multipolar interactions, chains attract each other at large separations, and repel at small separations. The cross-over distance between the dipolar attraction and multipolar repulsion defines a dynamic steady-state distance r * between a pair of chains, and this distance depends on Ψ and n (i.e., the length of the chain L = 2an) (Figs. 2b, S2d-g). For ϑ = 0 • , pairwise magnetic interaction is axisymmetric about the precession axis, with a negligible anisotropy produced by a small tilt of ϑ < 5 • (Fig. S2d).
The pairwise motion of chains are mainly determined by the hydrodynamic interactions, especially when they are in the cohesive regime. Precessing chains generate fluid flows, which lead to their self-propulsion and to hydrodynamic interactions influencing the motion of their neighbors.
Simulations were used for calculating the velocity field around precessing chains (Supplementary Note 3). The chain rotation perpendicular to the substrate leads to a rotational flow in the xy plane, leading to rotation of pairs about their center (Figs. 2c, S3a). For ϑ > 0º, rotation has a component parallel to the substrate, leading to the self-propulsion of chains and the generation of a secondary flow field that increases the velocity of neighboring chains via advection (Fig. 2d,

S3b-c).
Simulations combining magnetic and hydrodynamic interactions of chains quantitatively captured the dynamics of pair motion (Fig. 2e-h, Video S2, Supplementary Note 4). In the experiments, translational and rotational velocities of pairs increased while pair distance decreased with Ψ, which was supported by simulations ( Fig. 2f-h). Furthermore, our model also captured pair dynamics for chains of different lengths. For increasing n, pair distance and translational velocity increased and rotational velocity decreased (Fig. S4).

Dynamics and organization of homogeneous and heterogeneous clusters
Increasing the number of chains N, we observed their self-organization into clusters (Fig. 3a, b).
For homogeneous clusters, each chain was distanced evenly from their nearest neighbors (Fig. 3c Video S3). The pairwise distance decreased with Ψ (Fig. 3c). The pairwise distance also decreased from N = 2 to N = 3, but remained relatively unchanged with further increases in N (Fig. 3c). We ascribe this to the increased mutual dipolar attraction when multiple chains are present, which compacts the cluster. Cluster velocity increased with Ψ and N (Fig. 3d, Video S3). The latter effect can be attributed to an enhancement of collective hydrodynamic advection resulting from the superposition of flows generated by N chains.
To understand the effects of intragroup variations on the collective dynamics, we investigated clusters formed by heterogeneous members (Fig. 4). Dynamic self-assembly allows tuning the size distribution in small clusters (N = 7), enabling systematic investigation of heterogeneity. Group stability was evaluated based on two order metrics: rotational order parameter, R, and connectivity (Supplementary Note 6). R measures the degree of coherence of rotational motion of chains around the cluster center. R → 1 for perfectly coherent rotational motion, and R → 0 for no rotation 24 . Connectivity is based on total magnetic interaction potential; thus, quantifies the strength of cohesive interactions that hold the chains together.
Three clusters of varying degrees of heterogeneity were formed (Fig. 4a). Cluster heterogeneity was measured as the standard deviation of the number of particles in chains (σn), while the average ̅ was fixed at 3. Homogeneous clusters actuated at Ψ = 68 • and ϑ = 5 • remained as a single cluster, in which chains preserved their ordered spatial organization over time with R ~ 1 (Fig. 4b, c). However, when heterogeneities were introduced, clusters were more likely to break down, with decreasing R and connectivity over time as the chain interactions weakened (Fig. 4b, c, Video S4). Decreasing ϑ to 3 • reduced the cluster velocity and re-established the R and connectivity of the heterogeneous clusters (Fig. 4b, c, Video S4). Overall, heterogeneities cause a variance in chain velocities (Fig. S5) and weakens cohesive forces, resulting in cluster dissolution.
By decreasing the velocity, the balance of these competing effects is shifted towards cohesion.
To elucidate the group formation and dissolution dynamics, we investigated the chain trajectories after subtracting the mean cluster translation and rotation, revealing the internal motion of individual chains (Fig. 4d, Video S4, Supplementary Note 6). For stable clusters, chain displacements were constrained around their mean internal positions, where they perform small oscillations, as quantified by their mean-squared displacement (MSD) curves (Fig. 4e). We observed that the amplitude of these positional fluctuations increased with increasing heterogeneity and ϑ (Fig. 4f). Positional fluctuations were bounded to a finite amplitude when clusters exhibited structural ordering, and grew indefinitely when the clusters were breaking apart (i.e., Het2, Ψ68 • -ϑ5 • ) (Fig. 4f).

Formation of large clusters
Scalability of cohesive interactions enabled formation of large clusters (up to N = 53), which were formed by introducing more particles during the cluster formation (Fig. 5, Video S5). Global motion and internal fluctuations of such clusters were measured experimentally (Fig. 5a, b, Video   S6). We observed that the mean neighbor distance did not vary considerably with N (Fig. S6). On the other hand, cluster velocity increased (Fig. 5c, Video S6), continuing the trend observed for smaller clusters (Fig. 3d). We ascribe this observation to scaling effects associated with the collective hydrodynamic interactions with increasing cluster size. To assess this argument, we The internal motion of clusters reveals a tendency towards large positional fluctuations of chains with increasing N (Fig. 5b, Video S6). In small clusters (N ≤ 11), chains perform small fluctuations around their mean internal positions, which remain relatively fixed over long times, indicative of a solid-like order (Fig. 5b, d). As cluster size grew (N > 19), chains started moving inside the cluster, while remaining confined within the cluster radius RC (Fig. 5b, d). As such, amplitude of positional fluctuations grew proportional to RC when N > 19 (Fig. 5d). MSD data showed that the chain motion exhibits ballistic ~ t 2 behavior at short times (Fig. 5e). For t > 10 s, chain motion varied from small bounded oscillations evidenced in solid-like structuring in clusters for N ≤ 11, to liquid-like diffusive ~ t behavior confined to the cluster radius for N ≥ 15 (Fig. 5e).
We ascribe the tendency towards larger fluctuations with increasing N to the different distance-dependent decay rates of magnetic and hydrodynamic interactions. Magnetic dipolar interactions decay with r -4 , whereas the far-field of rotlet flow parallel to the wall decays with r -2 25,26 . Being short-ranged, magnetic forces holding chains together in a solid-like order rapidly reach saturation as the number of neighboring chains increases. On the other hand, hydrodynamic forces displacing chains keep increasing due to their longer range. It can be expected that hydrodynamic effects compete magnetic effects at a certain cluster size, and the cluster transitions to a liquid-like internal state.

Discussion
Collective motion manifests itself at all length scales relevant to biology 27 . Examples range from cytoplasmic transport in plant cells 28 , maze-solving slime moulds 29 , cooperative foraging in social ant groups 30 , migratory flocking of white storks 31 , and human motion in crowded environments 32 .
Inspired from natural systems, robotic swarm systems are being developed to address complex tasks such as collective construction and search operations 33,34 . A similar trend is prevalent in the microscale robotic swarms, with the aim of enhancing functional throughput, multitasking capabilities, and to impart microrobots with reconfigurability to enhance their adaptability to environmental constraints. However, a direct transfer of algorithmic approaches from the macroscale to microrobotic swarms present significant challenges due to challenges in miniaturization and powering of analogous components. Instead, microrobotic systems currently need to rely on micron-scale physical interactions for local coupling that generate global collective behaviors. Here, we demonstrated a microrobotic system where pairwise magnetic interactions can guide self-organization of cohesive clusters, which leads to a synergistic enhancement of the collective mobility over individual microrobots via hydrodynamic interactions.
Cohesive interactions play an important role in biological swarms (imagine a herd of sheep), where the group formation needs to be maintained without confining boundaries [15][16][17] . As opposed to purely attractive interactions that lead to the collapse of agents into tightly packed clusters, inclusion of short-range repulsive interactions can promote inter-agent spacing. In microrobotic swarms, several works have achieved cohesive organization in non-propulsive systems by combining repulsive capillary/hydrodynamic and attractive magnetic forces, however only at the liquid-air interfaces 35,36 . There is a great interest to form self-organizing swarms in fully liquid environments, due to their relevance to biomedical applications. In this work, we have achieved cohesive self-organization of microrobotic swarms in a fully immersed liquid environments by taking advantage of multipolar interactions resulting from anisotropic shape-magnetization of chains actuated with time-varying fields. Programming multipolar interactions holds promise for the design of advanced collective motion and manipulation of microrobotic swarm systems 6,8,22,23,37,38 , as we have shown here.
Swarm heterogeneity can enhance resilience of social collectives against random noise 39 , or can allow collaborative task division in robotic collectives 40 . On the other hand, introduction of large variances in population characteristics (e.g., velocity, interaction strength) may have deleterious effects on the order and cohesiveness of flocking swarms that interact locally, which can be regulated or alleviated by self-sorting and mixing mechanisms 41-43 . Despite its importance, heterogeneity has remained relatively unexplored in the field of synthetic active matter swarms. In the present system, heterogeneities mainly contributed to disordering of cohesive clusters due to the variance of individual mobilities. However, collective order can be restored by slowing down the swarm, which re-adjusts the competing effects of cohesiveness and individuality.
Self-organized cohesion can be further scaled to form large swarms. In our system, we observed that increasing group size enhances swarm mobility, akin to hydrodynamic cooperation observed in collectives of sperm cells 44 and of active colloidal rollers 26,45,46 . On the other hand, different scaling of magnetic (decaying with r -4 , short range) and hydrodynamic (decaying with r -2 , long range) interactions affect the spatiotemporal organization of chains as swarms grow larger.
We observed a solid-like spatial organization in small clusters, where magnetic cohesion is dominant over hydrodynamics. Two effects became comparable in larger groups, resulting in liquid-like dynamics, where chains displaced with respect to their neighbors. These states are highly similar to the flying crystal and the moving droplet formations in flocks interacting with cohesive alignment rules 16 .
In conclusion, we have presented mobile microrobotic swarms cooperating to generate a cohesive organization much larger than an individual microrobotic unit, which introduce synergistic advantages and display rich spatiotemporal organizations arising from collective dynamics. Our approach addresses the operation of microrobots as localized teams, which could inspire researchers in active matter and microrobotics fields for designing advanced collective systems for future applications in biomedicine, precision manipulation and manufacturing, and environmental sensing and remediation.

Experimental Setup
External magnetic fields required for self-assembly and actuation of chain microrobots were generated using a custom five-coil magnetic setup integrated on an inverted optical microscope (Axio Observer A1, Carl Zeiss) (Fig. S8). The magnetic coil system was arranged to generate up to 20 mT in in-plane directions and 10 mT in out of plane z-direction 7  precessing magnetic field with pre-determined tilt and precession angles is applied to assemble and actuate chain microrobots, which is detailed in Supplementary Note 1.

Dynamic Model and Simulations
Magnetic interactions between two chains were calculated by averaging dipolar interaction forces between particles over a cycle of magnetic field precession. We defined a characteristic magnetic interaction force, 0 = 12 0 ( ) 2 , for two chains separated by a distance L. Detailed methods for calculating magnetic interactions are described in Supplementary Note 2. where interaction forces between induced magnetic dipoles of particles ( ), particle-particle ( )

Supplementary Note 1: Applied magnetic fields
Precessing magnetic field is defined by two angles. Precession angle Ψ is the angle between the axis of precession (w) and the magnetic field vector (B), and the tilt angle ϑ is the angle between the precession axis and the normal vector to the planar substrate surface (i.e., z-axis) (Fig. 1a).
Applied magnetic field is varied in time by revolving the magnetic field vector about the precession axis, with angular velocity Ω via the following mathematical operation: where 0 = ‖ ‖ is the magnetic field magnitude, t is time, I is the identity matrix, [] × is the crossproduct operator, w = (0, sin(ϑ), cos(ϑ)), and n0 = (0, sin(ϑ + Ψ), cos(ϑ + Ψ)) is the magnetic field vector at t = 0.

Supplementary Note 2: Magnetic interactions between chains
For calculation of magnetic interactions forces between two chains (Figs. 2a-b and S2), we considered that each chain consists of n paramagnetic particles and precess by following the magnetic field given by Eq. S1. Distance vector pointing from the j th particle to i th particle is denoted by Rij = Rj -Ri. The induced magnetic dipole moment mi of each particle is given by = / 0 , where vp is the particle volume, χ is the volumetric magnetic susceptibility and 0 is the vacuum permeability. The interaction force between two magnetic dipoles is calculated with the following equation 1  Following, the time-averaged magnetic force acting between two chains can be calculated by summing the interaction forces between particles and averaging over a precession cycle, where the summation is performed over each particle pair i, j belonging to chains. Observing Eqs.
S2 and S3, the strength of the magnetic force between two chains varies with the following proportionality: where r is the distance between two chains. Therefore, for a pair of chains separated by a distance of one chain length (r = 2na), a characteristic magnetic interaction force can be defined as:

Supplementary Note 3: Numerical model for simulating flow fields around a single chain
Flow field generated by the precession of a chain near a wall was calculated with simulations using where interactions between particles i and j arise from magnetic dipole-dipole forces (fm), particleparticle (fb) and particle-wall (fw) excluded volume forces, and gravitational (fg) forces. Magnetic dipolar interactions between particles are calculated via Eq. S2. The grand mobility tensor couples the velocities of particles (̇) to the forces acting on each particle through contributions of self and pair hydrodynamic mobility tensors that account for no-slip boundary conditions at the substrate surface 3 . Simulations implement the mathematical expressions for the grand mobility tensor that were provided by Swan and Brady 3 . Following the approach presented by Sing et al. 4 , particle-particle and particle-wall excluded volume forces were modeled with modified Lennard-Jones force terms, ) (Eq. S10) with σ = 0.1a, h is the distance of a particle from the wall, and ϵ is sufficiently small to neglect attractive terms. fg = -Δρvpg where Δρ is the buoyant density of particles, and g is the gravitational acceleration. Eq. S8 was integrated with an explicit Euler scheme to obtain the trajectories of each particle in chains. Each cycle of chain rotation was divided into 5×10 5 time steps, and simulations were performed for 270 cycles (~30 seconds of real time experiments).

Supplementary Note 5: Reduced-order discrete chain model
For simulating the dynamics of clusters consisting of many chains, calculating the motion of each particle separately is computationally intensive. For this reason, we developed a reduced order simulation that models the dynamic of the collection of chains that constitute the cluster, in which each chain is treated as a discrete point ( Fig. S6 and Video S7). This model accounts for the time-averaged magnetic interactions and near-wall hydrodynamic self-propulsion and interactions between chains. The equation of motion for each chain is given by where the velocity of i th chain (̇), is the sum of its self-propulsion velocity, v0,i, the velocity of the hydrodynamic flow generated by all its j th neighbours at the position of i th chain, vh,ij, and the displacement velocity due to the magnetic forces imposed by its neighbors, vm,ij.
Self-propulsion of a chain arises from its self-advection under the hydrodynamic flow generated by its precessing motion. Similarly, a chain is also advected by the flow generated by its neighbors. The essential features of hydrodynamic flows were captured by modeling each chain as a rotlet singularity above a solid wall. Specifically, a rotlet solution provides the hydrodynamic velocity field generated by a point torque at Stokes regime. We consider that each chain has an effective hydrodynamic radius ah, and its rotation is given with the angular velocity vector Ω. The axis of rotation is tilted from the z-axis (i.e., normal to the plane of the substrate) by angle ϑ, therefore, the angular velocity vector can be decomposed into two components which are parallel and perpendicular to the substrate, such that Ω = (0, Ω|| , Ω⊥) where Ω|| = Ω.sin(ϑ) and Ω⊥ = Ω.cos(ϑ).
Chains are located at a distance h from the wall at r = (x, y, h). The solid wall imposes hydrodynamic no-slip boundary conditions, which can be accounted through (Eq. S17) In the simulations, we use the functional form given by Eq. S17, which lets specifying r * as an input parameter. Eq. S15 tells us that, if A and r* are known, then the only remaining unknown is k, which tunes the stiffness of the repulsive multipolar interaction. We found that setting k to different values in the range of 6 to 8 produce qualitatively similar results in our simulations.
Lastly, vm,ij is calculated by multiplying the magnetic interaction force with the effective mobility of a chain: (Eq. S18) The reduced-order model successfully captures the essential trends observed in experiments: Chains form cohesive clusters with a characteristic steady-state distance between neighbors, clusters performed rotation and translation, and the cluster velocity increased with number of chains. An example of simulated trajectories is displayed in Fig. S6c and Video S7. In order to achieve a quantitative fit between the model and the experiments for translation and angular velocity of clusters, we need to tune two parameters ah and h (Figs. S6a-b). An experimental determination of ah and h is difficult due to the anisotropic shape of chains. Also