Dynamics of cooperative excavation in ant and robot collectives

The solution of complex problems by the collective action of simple agents in both biologically evolved and synthetically engineered systems involves cooperative action. Understanding the resulting emergent solutions requires integrating across the organismal behavior of many individuals. Here, we investigate an ecologically relevant collective task in black carpenter ants Camponotus pennsylvanicus: excavation of a soft, erodible confining corral. These ants show a transition from individual exploratory excavation at random locations to spatially localized collective exploitative excavation and escape from the corral. Agent-based simulations and a minimal continuum theory that coarse-grains over individual actions and considers their integrated influence on the environment leads to the emergence of an effective phase space of behaviors, characterized in terms of excavation strength and cooperation intensity. To test the theory over the range of both observed and predicted behaviors, we use custom-built robots (RAnts) that respond to stimuli to characterize the phase space of emergence (and failure) of cooperative excavation. Tuning the amount of cooperation between RAnts, allows us to vary the efficiency of excavation and synthetically generate the entire range of macroscopic phases predicted by our theory. Overall, our approach shows how the cooperative completion of tasks can arise from simple rules that involve the interaction of agents with a dynamically changing environment that serves as both an enabler and a modulator of behavior.


Introduction
Collective behavior is seen in organisms across many length scales, from the microscopic to the macroscopic (Nowak, 2006;Camazine et al., 2020;Gordon, 1999;Seeley, 2009;Couzin and Krause, 2003). These behaviours are often functional and serve as solutions to problems associated with tasks that cannot be solved efficiently at the individual level and range from brood care to foraging for food, protection from enemies and predation of prey, building complex architectures etc (Feinerman et al., 2018;Ocko and Mahadevan, 2014;Hölldobler and Wilson, 2009;Peleg et al., 2018;Rasse and Deneubourg, 2001). Since collective behavior involves multiple individuals, this necessarily involves some form of communication and/or cooperation that takes different forms across scales -from quorum sensing in unicellular bacterium and slime molds, to the waggle dance in bees, and various forms of physical signal propagation in animal societies and human organizations (Rasse and Deneubourg, 2001;Alcock, 2001;Pennisi, 2009;Nowak, 2006;Elster, 1998;Couzin and Krause, 2003).
The importance of environmental signals is particularly clearly seen in examples of collective task execution in social insects that have a long history of documented cooperative behavior (Hölldobler et al., 1990;Gordon, 1999;Perna and Theraulaz, 2017;Mikheyev and Tschinkel, 2004). Superorganisms made of individuals respond to local stimuli with stereotypical actions that leave their 'mark' on the environment, creating a spatio-temporal memory, commonly known as stigmergy (Hölldobler and Wilson, 2009). While stigmergy is usually associated with scalar pheromone fields, a broader definition might include the use of signaling via chemical, mechanical and hydrodynamic means (Buhl et al., 2005;Mikheyev and Tschinkel, 2004), as has been quantified in recent studies of bees (Ocko and Mahadevan, 2014;Peleg et al., 2018). To understand how collective task execution arises, we need to understand how individuals switch from local uncoordinated behavior to collective cooperation that translates to successful task execution in different social systems. From a biological perspective, this naturally involves understanding the neural circuits, physiology and ethology of an individual. A complementary perspective at the level of the collective is that of characterizing a 'crude view of the whole', which entails the quest for a small set of rules that are sufficient for task completion, along with the range of possible solutions that arise from these rules that might be tested experimentally. And finally, given the ability to engineer minimally responsive biomimetic agents such as robots (Rahwan et al., 2019), a question that suggests itself is that of the synthesis of effective behaviors using these agents. This allows us to explore regions of phase space that are hard to explore with social insects, and also to learn about the robustness of these behaviors using imperfect agents in uncertain and noisy physical environments, before looking for them in-vivo.
Here we use an ecologically relevant task in carpenter ants Camponotus Pennsylvanicus: excavation and tunneling, to quantify the dynamics of successful task execution by tracking individual ants, create a quantitative framework that takes the form of mathematical models for agent behavior, and finally synthesize the behavior using robots that can sense and act. Our work complements and builds on earlier studies on excavation (Buhl et al., 2005;Tschinkel, 2004;Deneubourg and Franks, 1995;Deneubourg et al., 2002) in social insects that looked at the effects of population size and the role of cooperation on the efficiency of digging, while developing 1-dimensional models to understand the excavation process. We go beyond these studies by (i) quantifying the collective behavior of ants by tracking them in space-time, following the dynamics of their interaction, and the process of excavation of the confining substrate, (ii) developing a theoretical framework that couples the change in ant density, substrate density and the rate of antennation in space and time to capture the collective execution of the task in terms of a few non-dimensional parameters that define the range of behaviors of the agents, (iii) synthesizing and recreating this collective task using custom-built robots that can respond to each other and the environment . An important outcome of our study is a phase diagram that shows the emergence of different collective behaviors associated with task completion as a function of just two dimensionless parameters that characterize the local rules underlying individual behavior and the nature of communication between agents such as ants and robots.

Excavation in carpenter ants
We start with ants drawn from a mature colony of C. Pennsylvanicus that consist of a queen, the sole egg layer, and workers from three morphologically different castes -major, median, and minor (Hansen and Klotz, 2005). Although all ants perform different tasks like foraging, nest-keeping, and brood care to varied degrees, during excavation, major ants, equipped with their large mandibles, generally take the lead role, while median and minor ants transport the debris out of the nest. Ants communicate primarily through their antennae by using them to sense pheromones released by other ants and by touching other ants to identify their caste. It is this inter-organismal information exchange that enables the collective solution of complex tasks.
Our experiments consist of a dozen worker ants from the same colony that are anesthetized (using CO 2 ) and then brought into a confining ring-like corral made out of agarose (height 10mm, inner radius 35mm and outer radius 55mm) flanked above and below by two hard plastic sheets. To mimic their natural environment in a nest, we eliminated visible light and used infrared light to monitor the ants using video (see Figure 1(a)). We performed 4 experiments with a collective of 12 majors ants and 3 sets of experiments with a mixture of 4 major, 4 media and 4 minor ants. Once we introduce O 2 into the corral, the ants regain activity but stay still for a while before moving. They first exhibit wall-following until one or more of the ants initiates an exploratory excavation at a random location along the corral (ref Figure 2). After an initial exploratory phase the ants switch to an exploitative strategy in which they excavate a tunnel at a specific location and eventually break through the corral (see Video 1 and the sequence in Figure 1(b)). In contrast with the behavior of the 12 ant collective, when a single Major ant is introduced into the arena, the ant is unable to excavate through the agar barrier (see Video 1).
We can quantify this transition from rotationally isotropic exploration to localized excavation by considering both the behavior of individual ants or their effective density ϱa (r, ϕ, t) as a function of the polar coordinates (r, ϕ) and time t. We choose to use an effective coarse-grained density for two reasons: it is a more natural variable in the limit of large populations that vary in space and time, and is also amenable to building effective theories with fewer parameters that are easier to analyze and thus also compare to experiments. The ant density is obtained by averaging the position of the ants over a time window larger than the time taken for them to perform one task cycle , that starts with excavation at the boundary and ends with dropping debris in the interior of the corral (see Appendix 1 for further details). Over time, we see that the ant density becomes localized at a particular angle and location along the corral; here large-scale excavation eventually leads to excavation and escape from averaging along the radial direction. Ants start from an initially isotropic state and localize at an angle ϕ b along the boundary. T here is the excavation time. (e) Dynamics of the radial distribution of ant density P a r (r, t) as a function of radial distance, r obtained by averaging a sector of π/6 around the excavation site. We see that the ant density front propagates through the corral. The density is plotted for the same times as in (f) Evolution of the power spectrum |R(k, t)| 2 of first five Fourier modes capturing the number of tunnels formed during excavation R(ϕ, t) = ∑ kR (k, t)e ikϕ . Inset shows the real part of the Fourier coefficient, ℜ(R) at different time instants indicating that many modes are present in the boundary shape.
the corral (see Figure 2 and Appendix 1- figure 1 for the coarse-grained spatio-temporal evolution of the ant density, obtained by this averaging procedure). Simultaneously, collective excavation leads to an increase of the volume of excavated material, as shown in Figure 1(c) (see also Toffin et al., 2009). By averaging the ant density over radial positions, in Figure 1(d) we show the orientation distribution of the ant density P a ϕ (ϕ, t) =´ϱa(r, ϕ, t)dr is initially isotropic, and gradually starts to localize at a particular (arbitrary) value of the angle as time increases.
Averaging the density over the localized region, in Figure 1(e) we show the radial distribution of the ant density P a r (r, t) =´ϱa(r, ϕ, t)S(ϕ)dϕ (where S(ϕ) is a smoothing kernel localized around the excavation site) starts out by being initially uniform, and gradually propagates radially outwards as time increases. Consistent with localization and concomitant excavation (Figure 1(f) inset, Appendix 1-figure 2(c)), we see that the multiple azimuthal Fourier modes compete with each other initially before an elliptic mode (corresponding to a strongly localized state) is amplified as excavation progresses (shown in Figure 1(f), Appendix 1- figure 2(b)). All together, our quantitative observations show that an initially isotropic and homogeneous Figure 2. Evolution of the ant density field, ϱa(x, t) (in units of #/mm 2 ) as the tunneling progresses for experiments with 12 major ants. The density field is obtained by averaging the ant locations over 250 s during the tunneling process. In the second columns is the evolution of the boundary shape, R(ϕ) as a function of time where we see multiple excavation sites being explored before one of them succeeds. The darker spots in the image are the debris that the ants deposit as they excavate the boundary. distribution of ants in the corral induces exploration of multiple potential tunneling paths that transitions into the exploitative excavation of one specific location that eventually leads to an excavation route.

Model of cooperative excavation
In order to understand the dynamics of this cooperative excavation we first model the ants using discrete agents that mimic the microscopic behaviors of ants before turning to a coarse-grained field theoretic model for the evolution of the ant, pheromone and substrate density in space and time.
In the agent-based model each ant is represented as a circular disk of radius a with center position r j (t) and orientation p j (t) where j = 1 · · · n , n being the number of ants in the domain (see Figure 3(a)). We approximate the confining corral in the experiments using discrete boundary elements which the agents can pick and place in the interior of the domain (see Figure 3(b)). Initially, a random collection of agents engages in exploration within the corral in the absence of external gradients, consistent with observations (Trible et al., 2017) but their motion is rectified either by the presence of pheromone gradients or reinforcing antennating signals (Hölldobler et al., 1990;Reinhard and Srinivasan, 2009;Waters and Bassler, 2005;Gordon et al., 1993;Hillen and Painter, 2009;Toffin et al., 2009). Antennation involves information moving with the ants while pheromone gradients leads to information being laid down in the fixed environment. However, when ants move slowly relative to the time for the decay of the memory associated with antennation with other ants, the dynamics of both these processes is similar. Then the signals laid down (or transported) by ants increases locally at a rate proportional to their density (Gordon, 2021), and is subject to degradation and diffusion slowly. Accounting for these effects, we arrive at the following dynamical equations for the evolution of r j (t), θ j (t), c(x, t) as: Here, the orientation of the agent in Equation 1 is given by p j = (cos θ j , sin θ j ) with θ j being the heading angle, v o the characteristic speed of the agent, η j is a Gaussian white noise with correlation function The agents produce an antennating field at a rate k+ which decays at a rate k − centered around the agent, and captured by the function H(r j , a) = {1 if |x − r j | 2 − a 2 ≤ 0 , and vanishes otherwise}. We assume that the gradient in the antennating field along the local normal, on the right hand side of Equation 2, determines the rotation of the agents with G being the rotational gain. In order for the agents to initiate the excavation process, they can pick the elements from the boundary and drop them in the interior of the corral only when the local concentration of the antennating field is larger than a critical threshold c * , consistent with observations (Gordon, 2021;Gordon et al., 1993). Figure 3(b) shows snapshots (see Video 2 for a movie of the simulations) of the agent-based simulations following Equations 1-3 showing that the agents excavate successfully out of the corral when the gradient following behavior is strong (see Appendix 2 for details). Given this, we expect the time taken to  https://elifesciences.org/articles/79638/figures#video4 escape from the corral is a function of the number of agents. In Figure 3(c and d) we see that as we vary the number of agents from n = 1 − 100 , for very small or large number of agents in the corral, the agents are unable to escape over the time of simulations, Tstop (ref. Figure 3(c and d)), seen as saturation in the excavation time T/ts . In our agent-based simulations, we can encode the detailed behavior of individual ants and thus account for nuances and variations across the population. However, these simulations are computationally expensive as one needs to couple the dynamics of the antennating field (governed by a partial differential equation) with the motion of discrete agents while also evaluating the mutual interactions between all the agents in the corral. A complementary perspective that allows us to gain insights into the relevant parameters that govern the macroscopic dynamics of the collective is afforded by a theoretical framework that averages over the fast times and short length scale actions of the agents, considering spatial variations over scales much larger than a 'mean-free path' and 'collision time' associated with agent-agent interactions. Our effective theory attempts to couple three slowly-varying spatio-temporal fields: the ant density ϱa(x, t) , a communication field c(x, t) representing antennation and pheromone-based communication, and the corral density ϱs (x, t) , shown schematically in Figure 4(a). In the continuum picture, the agents' random motion is captured using diffusion of the density while the rectified motion due to pheromone gradients is captured through chemotaxis, in addition to being self-propelled with a velocity ua that is related to the local environment. Finally, motivated by observations of antennation (Gordon, 1999;Pagliara et al., 2018), we assume that when the ants are stimulated by the presence of the corral past a threshold of antennation, c * they start excavating. The rate of excavation is assumed to be proportional to the difference in the pheromone concentration relative to the threshold value (see further details). Accounting for these effects, we arrive at the following dynamical equations for the evolution of ϱa(x, t) and ϱs (x, t) that are coupled to Equation 3 for the evolution of the communication field: In Equation 4, the ant advection velocity is assumed to have the form ua = vo(1 − ϱs/ϱo)p where v o is the characteristic speed of the agents, and p is a unit vector pointing along the radial ( θ ) direction, and the term (1 − ϱs/ϱo) reflects the fact that excavating ants are slowed down by their labor; Da is the diffusivity of ants, χ is a chemotactic gain associated with the antennating-field-following behavior (related to the gain G in the agent-based model). Here is the average density of the ants defined by where is the domain size. This is a natural scale of the ant density as Equation 4 is in conservative form and the net density of the ants is preserved over the evolution. In Equation 5, k s is the rate of excavation of the corral and ϱ * a , c * are respectively the threshold concentration of ant density and antennating field required to initiate excavation. We assume that the behavioral switches have simple switch-like responses modeled here via the Heaviside function Θ(x) (or its regularization via hyperbolic or Hill functions). It is useful to note that in the absence of excavation dynamics, our framework reduces to the well known Keller-Segel model for chemotaxis (see Hillen and Painter, 2009 for a recent review) (also detailed in Appendix 2). The coupling of ant behavior to the dynamics of excavation introduces the all-important notion of functional collective behavior linking active agents, communication channels (the antennating and pheromone fields) and a dynamic, erodible corral that characterizes progress towards task completion.

Model parametrization and description
The evolution of the ant density in Equation 4 is a combination of three dynamical processes: ant migration, diffusion and biased motion due to antennating. There are three time-scales associated with these three processes: a diffusion time-scale τa ∼ l 2 /Da , a collective migration time-scale τv ∼ l/vo and a time-scale associated with taxis τx ∼ l 2 /χco , where l is a characteristic length-scale. This last scale can be either the width of the corral to be excavated L (which is assumed to be of same order as width of initial ant density profile l a ), the length-scale associated with the balance between antennating field diffusion and decay, l ∼ (Dc/k − ) 1/2 or the length-scale due to the advection of ant density and diffusion, l ∼ Da/vo . The antennating field in Equation 3 is governed by three processes, the generation of the antennating field, as well as its decay and diffusion. This leads to three more time-scales : an antennating field production time-scale τ+ ∼ co/(k+ϱo) , a diffusion time-scale τc ∼ l 2 /Dc , and a decay time-scale τ − ∼ 1/k − . Lastly, the dynamics of excavation from the corral which follows Equation 5 is governed by a characteristic time-scale τs ∼ 1/ks . The list of all seven time-scales and length-scales associated with the different processes in the model are in Appendix 2-table 2. In terms of the different time-scales (see Appendix 2 for a list along with their ranges), there are a total of six dimensionless parameters, of which two non-dimensional numbers are qualitatively important in capturing the etho-space of collective excavation: (i) the scaled cooperation parameter defined as C = τa/τx = χco/Da which determines the relative strength of antennation (gradient-following) to ant diffusion with c o being the maximum amplitude of the antennating field, (ii) the scaled excavation rate, E = τv/τs = ksl/vo . Here, l/vo is the characteristic time-scale of ant motion, with l ∼ min[(Dc/k − ) 1/2 , la] , where l a is the ant size (see Appendix 2 for details). The other four dimensionless parameters follow from the ratio of the time scale of ant motion and the diffusive time-scale as V = τx/τa = vol/Da . The ratio of the rate of production of pheromone and the rate of diffusion or decay, leading to the parameters k ± = τ − /τ+ = k+ϱo/(k − co) and Dc = τ − /τc = Dc/(l 2 k − ) so that the complete set of non-dimensional numbers that capture the dynamics of the ant collective is given by In terms of these parameters, the dynamics of the ant density, the antennating field and the corral density given by Equations 3-5 can be written in non-dimensional form as To complete the formulation of our model, we also need to specify some initial conditions and boundary conditions for the ant density, the pheromone density, and the location of the corral boundary which are detailed in the Appendix 2.

Linear analysis
Before we consider the different limits of the phase-space defined by the non-dimensional numbers, we show that the excavation process is an instability triggered by the scaled excavation parameter E in the system. Starting with the homogeneous state ϱ ss a = ϱ * a , c ss = c * = k+ϱo/k − , ϱ ss s = 1 which satisfies the Equations 6-8, and perturbing about this configuration using a plane wave ansatz (in 1D) we write: where we assume that ||ρa||, ||c(k)||, ||ρs(k)|| ≪ 1 . Then the linearized counterparts of the Equations 6-8 for the ant density, antennating field and the corral density read as: (Ωk 2 )ρa + ikVρsϱo = k 2 Cc ,: c =k ±ρa /(Ω + 1 + Dck 2 ) , Ωρs = −Eρs/2 . From this, we see that the growth rate Ω = −E/2 , is independent of all other parameters in the system, i.e. excavating begins when E > 0 , once the ants have created a sufficiently large spatially diffuse antennating field. To understand the dynamics of excavation of the corral and the different phases of collective behavior, we now explore the role of the other non-dimensional numbers.

Limits of phase-space
Next we discuss the different limits of the phase-space defined by the non-dimensional numbers {C, E, V,k ± , Dc} and the thresholds ϱ * a , c * .

Box 1. Ant behavior → Model → Robot behavior
Ants inside the corral move around, communicating with each other using their antennae before they cooperatively excavate the agarose corral. Though the detailed spatiotemporal dynamics of each ant's behavior is different at the microscopic level, we see that the cooperation between the ants results in a persistent density front (see Figure 1 (d, e) and Figure 2) that excavates the substrate. In the theoretical description of the collective's dynamics, the relevant behaviors are encoded through mutual interaction between the ants (via the antennating field) and the substrate. Such a description also inspires the robotic mimics that capture the ant collective's average behavior. We list below the comparison between relevant behaviors in ants and their analogous encoding in the theoretical model as well as in the robots. speed v b is stored internally and the rate of change Ω of the heading is a function of the cooperation parameter C, the photormone concentration c , and a stochastic process W (Brownian motion). A photormone threshold c * determines whether an object is grasped (with probability E) after it is detected by the distance sensor. (d) Orientation distribution of the RAnt density P r ϕ (ϕ, t) as a function of the azimuthal position ϕ is the orientation of the excavated tunnel. The density is plotted for different times. (e) Radial distribution of the RAnt density P r r (r, t) within a sector of π/2 centered around the position of the excavated tunnel as a function of distance from the center of the arena r . The density is plotted for the same times as in (d). (f) Confinement area A(t) as a function of time, normalized by initial circular confinement with radius Ro for different cooperation parameter C. (g) Normalized excavation time T as a function of cooperation parameter C, averaged over 5 experiments per cooperation parameter. Every experiment was run until the first RAnt excavated out or the experiment duration exceeded 15 min.
Small thresholds, when ϱ * ≪ ϱo and c * ≪ co When ϱ * a ≪ ϱo and c * ≪ co , we see the appearance of partial tunneling even with an initially inhomogeneous ant density ϱa , independent of the pheromone dynamics. However, depending on on the value of the ratio τs/τv , the ants can either excavate through the corral completely ( τv/τs ≪ 1 ) or partially ( τv/τs ≤ 1 ) (ref Appendix 2-table 3). If the ants are moving randomly, i.e. in the diffusiondominated regime, they can still tunnel through the corral if τc ∼ τs and partial tunnel through the corral if τc ≲ τs . In non-dimensional terms, this translates to the relations V ∼ O(1), C ≪ 1 or V, C ≪ 1 and E ∼ O(1) for the corral evolution. (Appendix 2- figure 1 shows the results of simulations of both the tunneling and the partial tunneling behavioral phases).
Cooperation dominated regimes when C ≫ 1 and E, V → 0 For efficient excavation, the ants need to work collectively by being localized and excavating fast. Spatial localization leads to cooperation via feedback from the antennating field (see Figure 4(b)) while for successful excavation, ants need to migrate towards the corral and tunnel through it, so that their effective speed v o needs to be non-zero. To quantify these behaviors, we first look at the dynamics of the ant density and the antennating field in the absence of migration i.e. V → 0 or corral evolution. This leads to three regimes: • Diffusion dominated regime: When the antennating field diffuses rapidly, i.e. Dc ∼k ± ≫ 1 , then the equations for the evolution of the antennating field and the ant density, Equations 3-4 simplify to • Chemotaxis dominated regime: When the chemocactic coefficient χ is large, i.e. in dimensionless terms C ≫ 1 , the ant collective gets jammed. To see this we linearize the Equation 11 about a uniform ant density ϱo and recognize that this leads to an effective negative diffusivity and thus a spatio-temporal focusing of the ant density; we leave a detailed analysis of the characteristics of this for future study.
To understand the balance between diffusion of the antennating field and its decay, we note the appearance of a natural length scale l ∼ (Dc/k − ) 1/2 which defines the zone of influence of the field and provides a measure of the non-dimensional tunneling rate indicated in Figure 8. All together, our analysis shows that the dynamics of the antennating field controls the aggregation or diffusion of ant density. But for efficient excavation, especially when the activation thresholds for excavation and localization ϱ * a , c * are large, we need both cooperation and finite velocity of migration. A catalog of the various regimes associated with partial tunneling, jamming, or diffusion as the dimensionless problem parameters are varied is listed in Appendix 2-table 3.
To understand how these different limits translate to the dynamics of excavation from the corral induced by the ants, we now consider the case when E, V ̸ = 0 , and solve the governing Equations 4; 5 in a one-dimensional setting (ref Appendix 2). We see that we can capture the two limits of the excavation behavior seen in experiments; for large excavation rates E > 1 and cooperation parameter, C > 1 , we see coordinated excavation (shown in Figure 4(b) and Figure 5 in a two-dimensional setting), while decreasing the cooperation parameter leads to disorganized excavation (shown in Figure 4(c)) (see Appendix 2- figure 1). While a direct comparison with the behavior of ants is not easy owing to the difficulty of inferring the dynamics of information transfer through antennation, Figure 8. Phases of cooperation Phase-diagram of cooperative task execution with different phases seen in ants and RAnts. In the robotic experiments we tune the Cooperation parameter C and the Excavation rate E while in the ant experiments we change the caste mixture. In the ant experiments we see the jammed and diffused phases transiently before the ants relax to cooperative excavation.
the minimal assumptions we have made about the antennating field dynamics suffice to capture the macroscopic behavior of the collective. All together, our agent-based model and the phase-field model shows the emergence of cooperativity without the need for a plan, optimization principle, or internal representations of the world; instead environmentally mediated communication between agents (Mataric, 1993) coupled to local behavioral rules suffice to realize robust excavation.

Robotic collective excavation
Although our quantitative observations of the collective behavior of the ants is qualitatively captured by both our agent-based and continuum models, a natural question we can ask is whether the coarse-grained averaging over of the communication field affects the emergence of the task in experiments, especially since we are unable to measure or directly control the microscopic behaviors of the ants. To go beyond our ability to explain the observations of ant behavior using our theoretical framework, we asked if we might be able to recreate the behavior in artificially engineered mimics, and probe a larger range of the parameters and phase-space spanned by the scaled excavation and cooperation parameters C,E , than our experiments allowed us to -see Table 1 for a list of the relevant variables across ants, models and robots Figure 5.
For this, we turn to a robotic platform to synthesize collective functional behaviors that arise from simple behavioral rules underlying individual programmable robots. Our custom designed robot ants (RAnts) are inspired by many earlier attempts to create artificial agents that are mobile and follow simple rules (Braitenberg, 1986;Brooks, 1991;Simon, 1996), can respond to virtual pheromone fields (Sugawara et al., 2004;Garnier et al., 2007) and are capable of robotic excavation (Aguilar et al., 2018). Our autonomous wheeled robots can exhibit emergent embodied behavior (Bricard et al., 2013), and are flexible enough to allow for a range of stigmergic interactions with the environment (Werfel et al., 2014;Petersen et al., 2019). This is made possible by having each RAnt equipped with an infrared distance sensor to detect obstacles and other RAnts, a retractable magnet that can pick up and drop wall elements with a ferromagnetic ring (shown in Figure 6(a)), and the ability to measure a virtual pheromone field generated by a light projected (from below) onto the surface of a transparent arena they operate in (see Figure 6(a, b), Theraulaz and Bonabeau, 1995;Sugawara et al., 2004;Garnier et al., 2007;Wang et al., 2021). The intensity of this 'photormone' field follows the antennating field Equation 2 and thus follows the dynamics of a field that is linked to the locations of the RAnts and diffuses and decays away from it. The photormone field is realized by a projected luminous field on the arena, which the robots can sense. This allows us to use a local form of Equations 4; 5 to define a robot's behavior in terms of an excavation rate E, a cooperation parameter C, and a threshold concentration for tunneling c * . This is encoded in the behavior-based rules (see Figure 6(c) and Appendix 3 for more details), that induces the following behavior: (i) follow gradient of projected photormone field; (ii) avoid obstacles and other RAnts at higher photormone locations; (iii) pick up obstacles from high photormone locations and drop them at low concentration levels. Since the robots have no symbolic representation of the different signals they sense (e.g. they cannot distinguish another RAnt from a wall element, since both merely produce a bump in the sensor signal), the observed behavior emerges from this simple sequence by depending on the current state of the environment and the robot.
Varying the parameter C ∈ [0, 1] allows us to tune the individual behavior from random motion ( C = 0 ) to tracking the photormone gradient ( C = 1 ) (see Video 3). Varying the non-dimensional excavation rate E changes the frequency at which the robots execute pick-and-drop behavior with detected objects, and serves to mimic what arises in ants as a function of their morphology and caste (see Appendix 1 for more details). For specific values of these parameters, we followed the collective behavior of RAnts by averaging their position over several pick-and-drop timescales to obtain the RAnt density field ϱr (r, ϕ, t) , just as for ants. When all the RAnts are programmed to have a cooperation parameter C = 1 , RAnts initially explore the region without picking the boundary element until the photormone concentration c ∼ c * , which happens once a particular location has enough visits by other RAnts. Just as for ants, we calculate the radially averaged RAnt density P r ϕ (ϕ, t) =´ϱr(r, ϕ, t)dr ; Figure 6(d) shows how RAnt density localizes at a (random) value of the azimuthal angle. As excavation progresses, the RAnt density propagates radially outwards as a density front just as in ants, shown in Figure 6(e) in terms of the quantity P r r (r, t) =´ϱr(r, ϕ, t)dϕ (also shown in Figure 7 for different trails when C = 1 ). Concommitantly, as excavation progresses, the corral area increases (Toffin et al., 2009); interestingly the scaled corral area A(t)/πR 2 0 is independent of the cooperation parameter C as shown in Figure 6(f) (all RAnts were programmed to have the same excavation rate).
However, cooperation does change the time for excavation; in Figure 6(g) we show the average excavation time (scaled by the characteristic time it takes for a rant to traverse the arena) and see that T/ts decreases with an increase in the cooperation parameter C . RAnts excavated out every time for C > 0.5 , but are unable to complete excavation for low values of the cooperation parameter (within a 15-min time window). Our results show that it is the localized collective excavation of RAnts mediated by photormone-induced cooperation that is responsible for efficient tunneling and excavation; for low values of C , tunneling is defocused and global, and thus not as effective (see Appendix 3- figure 2). When E → 0 (vanishing probability for a successful pick up) but C is large(see Figure 8 and Appendix 2 for theoretical predictions), the RAnts get jammed because they follow the photormone field they generate but are unable to tunnel through the boundary constriction. On the other hand, when E <1 and C <1 the agents do not cooperate and their diffusive behavior prevents successful tunneling. The range of strategies can be visualized in a twodimensional phase space spanned by the variables E and C shown in Figure 8. Low values of C and E lead to diffusive (and non-functional) behavior, while high values of these variables lead to coordinated excavation, with the other two quadrants corresponding to jammed states (large C , small E ) and partially tunneled states (large E , small C ). Interestingly, these states are also observed as transients in our ant experiments, for example in the initially diffused state that is characterized by random motion inside the corral, when transiently jammed states and partial tunneling occur (see Videos 1 and 4).

Discussion
Our analysis of collective behavior in a functional task, excavation, uses quantitative observations of ants to build theoretical and computational models to explain them, and recreate these behaviors using a swarm robotic system (see Video 4 for a summary). Our simple dynamical models involving individual agents as well as an effective continuum theory provide a phase diagram that shows how the transition from an individually exploratory strategy to an exploitative cooperative solution is mediated by the local chemical and mechanical environment. Our study suggested algorithms that we then deployed in an engineered system of robots that individually follow a minimal set of behavioral rules that mould the environment and are modulated by it.; the malleable environment serves both as a spatial memory as well as a computational platform (using the spatiotemporal photormone field and the corral). Our simulations of agent-based models and robotic experiments further suggest that a coarse-grained framework linking behavior, communication and a modulated environment is relatively robust to failure of and stochasticity in the behavior of individual agents (i.e. variations in initial conditions and number of agents), in the communication channels and in the corral geometry, in contrast to engineering approaches that aim to control all agents and optimize costs. Different strategies such as collective excavation, jamming, and diffusion then arise as a function of the relative strength of the cooperation (representing the ability to follow gradients and detect threshold values) and excavation parameters (representing the ability to move material), as manifested in a phase diagram, and the emergence of cooperation arises due to the relatively slow decay of an environmental signal (the pheromone/antennating/photormone field), coupled to a threshold excavation rate. Since the ability to solve complex eco-physiological problems such as collective excavation is directly correlated with a selective (functional) advantage in an evolutionary setting, perhaps collective behavior must always be studied in a functional context. L Mahadevan The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Supplementary files • MDAR checklist
Data availability All the data used to generate the figures in the article are available here: https://github.com/sgangaprasath/rantIFigData (copy archived at swh:1:rev:ba2c6291882cf2355c0fc5d27384a8ce0dc48cc5). The simulation code used in the article is also available in the same folder.

Appendix 1 Ant experiments Experimental setup -handling ants
We collected two queen-right mature colonies of Camponotus pennsylvanicus, established in logs of fallen trees, from the Middlesex Fells Reserve (42.45 °N,71.11 °W) in August 2019. Each mature colony consists of three morphologically distinct castes of worker ants: major, media and minor, with an average body length of 7 mm , 5 mm for media, and 4 mm respectively. We placed the collected wooden logs housing those colonies in two separate plastic "home" boxes. We coated the inner wall of each home box with ant-slip Fluon to prevent the ants from escaping the home box. Each home box was connected to a foraging box by a tube through which ants travelled to and fro. We kept the whole set up in the laboratory with a 12 hour light-dark cycle, 30°C temperature and 50-70% relative humidity. Before we moved to the next phase of the experiment, i.e. the data collection, we waited for the ants to resume foraging and excavation of woods (for expanding their galleries) inside their home wooden log; this took 3-5 days after the relocation.
Appendix 1-figure 1. Dynamics of ant density field, ϱa(x, t) (in units of #/mm 2 ) obtained by averaging the ant location and the boundary shape R(ϕ) when 4 ants each of major, media and minor types are confined inside the agar ring for different trials.
About 10 minutes prior to the experiments, we collected ants engaged in wood excavation from the surface of the nest log. We used insect aspirators for collecting the ants. Once we collected all ants needed for the experiment, we subjected the ants to Carbon dioxide anaesthesia for 1 minute. Next, we placed the anaesthetised ants in the agarose well in the experimental arena; we placed each ant at least 1 cm away from any other ant. Ants regained their activity in the next 5-10 minutes.

Experimental setup -confinement
For the next phase of the experiment, we needed to confine the ants in an excavatable enclosure. This is the corral that the ants need to bite through to free themselves. We used a ring-like confinement made of agarose gel, with a height of 10 mm, an inner radius of 35 mm and outer radius of 55 mm, making the ring 20 mm thick. To make a precise shape of the ring repeatedly, we custom-built a casting mold made of acrylic plastic. We started preparing for the Agar ring before we collected the ants. For making the ring, first, we thoroughly mixed 3 gm of Agar powder in 100 ml of tap water. We then warmed the solution in a microwave oven until the solution started bubbling and appeared clear. Next, we poured the solution in the plastic mold, and kept it in 30 C temperature for 25 minutes; the agarose gel solidified and become opaque during this time. Once the agarose turned solid, we placed the ring on top of plastic sheet in the arena. Next, we placed the ants inside the ring and put a perti dish lid on top of the agarose ring. Thus, we confined the ants -with a solid plastic floor and ceiling, and an excavatable agarose gel wall. A schematic of the set-up is shown in Figure 1(b).

Experimental setup -arena and video recording
The arena consists of a piece of white 3 mm thick plastic sheet as the substratum, illuminated with infrared back-light, and surrounded by a 1.5 cm high plastic wall coated with Fluon ant-slip. We placed a Point Grey (FLIR) Grasshopper3 GS3-U3-41C6NIR camera, fitted with a 65 mm macro lens, on top of the arena to capture the view of the whole ring. The camera recorded the videos with 30 fps recording speed and 1024×1024 pixels resolution.

Markerless tracking
Leveraging an open source, deep-learning based pose estimator package SLEAP (Pereira et al., 2020), we track 3 body parts in each ant -head, thorax, and abdomen (gaster). Sample results obtained from this tracking is shown in Appendix 1-figure 2(e) and in (f − h) we quantify the noise statistics of ant motion and its orientation using the tracking data. Ants initially move randomly in the confinement and one of the ants starts the excavation process after which several ants start excavating cooperatively at the same location. When the tunneling happens, all the ants are orientated along the tunnel. We see this through the progression of the orientation distribution of ants P(θ, t) in Appendix 1- figure 2(j). To characterize the localization in ant orientation as the excavation proceeds, we use a von-Mises distribution (the analog of a Gaussian distribution for a periodic variable, given by P(θ, t; µ(t), K(t)) = exp [K cos(θ − µ)]/2πI 0 (K) ) of the ants (where μ is the mean local orientation associated with location of tunnel along the boundary). In Appendix 1figure 2(k), we see that over time, K(t) increases, i.e. the variance decreases. During the excavation process, ants bite through the corral, carry the debris from the excavation site and drop it in the interior of the confinement. This happens over and over again until all the ants excavate out. We see this captured in the oscillations of the location of ants as shown in Appendix 1-figure 2(i).

Average dynamics
We have a total of 7 sets of experiments with four sets of experiments with a collective of 12 majors ants and 3 sets of experiments with a mixture of 4 major, 4 media and 4 minor ants. Using the recorded video of the ant excavation dynamics, we threshold the intensity to extract only the ant boundary and average the ant dynamics over 250 secs. This gives us a density field of ants representing the locations where the ants have been and the amount of time they spend. We found in our experiments that each ant bites the corral, picks the bitten piece and transports it into the interior of the confinement. This process takes approximately 60 secs (see Appendix 1-figure 2(i)) and we would like to average the ant dynamics over several 'turn over' time-scales. We chose 250 secs and the obtained density field is shown in Appendix 1-figure 2. We perform this averaging for the experiments with all major ants as well as the mixture of different castes. In all the experiments, an ant density front propagates through the corral as they excavate and gradually tunnel through.

Boundary tracking
From the recorded videos, we also track the locations in which the ants excavate for creating the tunnel. For that, we used a custom image processing Matlab script. First, we created a mask superimposing on the area encircled by the inner ring of the corral; we colored the mask with a shade different from the corral. When ants excavated the corral, the Matlab script could detect the difference in the shade/color of the excavated area. Using this contrast, we track the continuously changing boundary.

Appendix 2 Agent-based model of cooperative task execution
The results shown in Figure 3 are based on a numerical simulation where discrete agents operate in a continuum scalar communication field, and governed by Equation (1), Equations 2; 3 . Some additional behavioral rules have to be defined to model the interaction of agents with the substrate. We realize the substrate by discrete obstacles arranged in a circular ring. Agents will attach to an obstacle if they are within the detection range, l d , and if the measured communication field value is above the threshold, i.e. c ≥ c * hi . The agent will then reverse its direction of motion, by changing the sign of G in Equation (2). This results in a gradient descent behavior and the attached obstacle will be detached once the measured communication field value satisfies c < c * lo . After detachment, the sign of G is changed again. If agents encounter other agents or obstacles within the detection radius but c < c * , the agents will avoid the obstacle by turning randomly.
There are a few tuned behaviors we implemented to allow scaling the simulation to larger numbers of agents while maintaining the tunneling behavior. First, the gradient ∇ ⊥ c in Equation (2) is is multiplied by a sigmoid function to limit the turning rate of the agents. Second, the noise term in Equation (2) was set to zero for this simulation and the only source of randomness are the random turns during obstacle avoidance. Third, agents pause for t p1 when they encounter an obstacle and for t p2 when picking up an obstacle. This helps disrupting potential "pheromone traps" to be formed where agents are bound to a region of space due to a high field concentration.
The simulation parameters are described in the following table. All parameters are nondimensionalized by the corral size L as a natural length scale and the base speed of the agents, v 0 as a natural speed.

Continuum model of cooperative task execution
The dimensional equations for the ant-density ϱa(x, t) , antennating field c(x, t) and the corral ϱs (x, t) are given by, where velocity of the collective is ua = vo(1 − ϱs/ϱo)p , capturing the reduction in velocity as the ant collides with the corral. We approximate the Heaviside function, Θ(x) here using the hyperbolic function [1 + tanh(x)]/2 . In the coarse-grained picture describing the collective tunneling seen in experiments the relevant variables (shown schematically in Figure 4) are the density of ants, ϱa ; their velocity, ua ; amplitude of the antennating field, c ; the density of corral, ϱs . Here we discuss the limits of phase-space that are not described in the main text i.e. when E, ̸ = 0 and also the simulation details. Antennating field diffusion-decay τa ∼ l 2 /Da

Ant diffusion
Limits of phase-space when E ̸ = 0 Different phases of task execution/failure arise when the excavation parameter E and the cooperation parameter C are varied. In the cooperation dominated phase if the excavation rate of the agents is small, they get jammed and the analysis in the previous section holds true. When the cooperation among the agents is low, we have C ≪ 1 which results in diffusion dominated regime. Based on the strength of the excavation parameter E , the corral can be partially tunneled or just diffuse. Since we assume that the relevant length scale is of the same order as the width of the corral, L ∼ l , our analysis reduces to different phases based on whether E ≫ 1 (where we get partial-tunneling) or E ≪ 1 (we get diffusion). Based on this we get partial tunneling or diffused phase as listed in Appendix 2-table 3. In Appendix 2-figure 1 we show results from 1-D simulations highlighting the effect of different terms we have discussed from Equations 4 and 5 corresponding to different parts of the phase space of cooperative excavation. In the ant density diffusion dominated regime, i.e. C ≪ 1, E ≪ 1 , shown in Appendix 2- figure 1(b), there is little cooperation; rapid diffusion with slow excavation results in no tunneling. As we have seen in Figure 4(b, e), tunneling and partial tunneling are inferred through the ultimate state of the corral and the ant-density. In Appendix 2- figure 1(d, e) we show how the relative rate of the antennating field diffusion compared to decay, i.e. Dc ∼ O(1) leads to either tunneling or partial tunneling as we vary the cooperation parameter C . Decreasing C causes the maximum of ϱa(x, t) , c(x, t) to go down (ref Appendix 2-figure 1(c)),and the width of the initial ant density field increases. Increasing C leads to successful tunneling driven by the propagation of the location of maximum ant density x loc due to excavation of the corral. Furthermore, we see that the ants can be jammed either because the antennating field diffusion dominates, i.e. Dc ∼k ± ≫ 1 , or because of the same field decays rapidly, i.e. k ± ∼ O(1), Dc ≪ 1 . In both these cases however cooperation is what drives the aggregation. Lastly, we see that in order to achieve partial tunneling there are several routes depending upon the relative magnitudes of {C, E, V,k ± , Dc} listed in Appendix 2-

Partial-Tunneling
Appendix 2-table 3 Continued on next page
Appendix 2-table 3 Continued the wall element is slightly lifted from the ground (≈1 mm) for transportation due to the elevated position of the magnet relative to the ferromagnetic ring. RAnts have two typers of sensors; two light sensors (Adafruit ALS-PT19) located at the bottom left and right of the RAnt (relative to the direction of travel) and an infrared (IR) distance sensor (Everlight ITR20001 opto interrupter) capable of detecting objects within 3 cm in front of the RAnt. The chassis of the RAnt is 3D printed using acrylic styrene acrylonitrile (ASA) and supports all the internal components. The wheels of the RAnts were printed with the same material. Due to the design of the wheel arrangement, which was inspired by the zooid robots (Goc, 2016), we require two small steel caster balls of 3 mm in diameter that help stabilize the RAnt. The steel balls can be pressed into the bottom of the 3D printed chassis. A 3D printed case made of ASA encloses all internal components of the RAnts, except for a small switch to power the RAnt on or off. A small blue sticker of 6 mm in diameter was placed on the center top of the case and is used for tracking of the RAnt's position with the webcam mounted above the arena.

RAnt programming
The RAnt behavior is coordinated by the microcontroller which we programmed according to the pseudocode shown in Algorithm 1. The program is initialized with a variable d that encodes the direction of travel (1 for forward, -1 for backward), the cooperation parameter C ∈ [0, 1] , the RAnt's base speed v b , and the light intensity threshold c * . This threshold was set at 50% of the maximal light intensity that can be generated by the photormone field multiplied by the cooperation parameter, i.e. c * = 0.5 × cmax × C .
After initialization, the program enters a while loop which is running until the RAnt is switched off or the battery voltage drops below 3.5 V. The loop starts with setting the heading of the RAnt, which effectively sets the turning rate. The turning rate is a function of the cooperation parameter and a stochastic process W (Wiener process) which is integrated in the microprocessor. The turning rate follows the equation with c L and c R the photormone intensity measured in the left and right light sensors, respectively, cmax is the maximal photormone intensity measurable by the sensors, and b = 0.3 is a fixed amplitude. Using a sine function we map the stochastic process W to the range [−1, 1] to avoid getting stuck in constant rotation for large excursions of W . The first term in Equation 15 corresponds to phototaxis using the projected photormone and the second term to a random walk. We can tune the influence of either terms with the cooperation parameter C from pure phototaxis at C = 1 to a random walk at C = 0 . The turning rate is used to define the rotation speed of each wheel. One wheel is always turning at a base rate ω 1 = ω b = v b /R (with R the wheel's radius) and the other wheel at The assignment of ω 1 and ω 2 to the left and right wheel is flipped according to the sign of Ω . With this definition, at a value of Ω = ±1 a RAnt turns on the spot without any translation and at Ω = 0 the RAnt moves on a straight path without rotation. After the heading was defined and the turning rates sent to the motor driver, the distance sensor is checked for any obstacles that are present up to 3cm in front of the RAnt. At the same time, the light sensors are checked and compared to the threshold value c * . If an object is detected and the photormone value exceeds c * , the RAnt performs a fetching manoeuvre that consists of engaging the magnet with probability E , moving forward for a second with half the base speed v b then move backwards for the same amount of time. After the fetching manoeuvre, the direction parameter is inverted, i.e. d = −1 . If an object is picked up with the magnet after the fetching maanoeuvre, the distance sensor will report a detected object as long as it is attached to the magnet. Since d = −1 , the RAnt will perform the same type of gradient driven locomotion described in Equation 15 and Equation 16 but the sign of the signal sent to the motor driver will be inverted, resulting in a reverse motion of the RAnt. If an object is detected, but the photormone concentration in both sensors is lower than c * , an avoidance manoeuvre is performed which consists of a random rotation in place in any direction with the intent to turn away from the detected obstacle.
The next if-statement checks again if an obstacle is detected, but without the condition that the direction parameter is equal to one. If no obstacle is detected, the direction parameter d is set to one and the magnet is disengaged. This guarantees that if a fetching manoeuvre is performed but the wall element was not picked up or the object was another RAnt, the RAnt goes back to moving forward.
The last if-statement checks whether the RAnt is in the reverse mode d = −1 and if the photormone concentration dropped below the threshold c * . if both statements are true, the magnet is disengaged, dropping any potentially picked up wall elements, and the direction parameter is set back to d = 1 . In order to avoid the RAnt from picking up the just dropped element, it performs a random rotation in place in any direction before going back to the start of the main loop.

Experimental set-up
The photormone was projected with an Epson EX9200 projector onto an acrylic sheet with a translucent top, which served as the surface on which the RAnts are operating. The projector uses three-chip digital light processing (DLP) which is required for the light sensors in the RAnts to pick up the photormone field. Tests with single-chip DLP projectors generated large noise in the light sensors and phototaxis was not possible. The dynamics of the photormone field is a function of the RAnt's positions and is given by with c = c(x, t) the photormone concentration at position x = [x, y] and time t , D = 10 −5 m 2 s −1 the diffusion coefficient, k M = 1 s −1 the decay rate, k P = 6.5 s −1 the photormone production rate, n the number of RAnts detected in the arena, N (r i , Σ) a bivariate normal distribution with the position of the i th RAnt r i as the mean and covariance Σ with diagonal entries σ 2 = 10 −4 m 2 . The position of the RAnts are used as the centers of sources of photormone. If a RAnt is not moving, photormone is built up with rate k P at that location over time and diffuses out. When the RAnt moves to a new location, the built up photormone decays with rate k M . The reasoning for the parameter choices is as follows. The parameters were tuned to allow for a RAnt located at one position for one second to leave a detectable trace for 5 seconds. During that time, another RAnt moving at base speed v b ≈ 5cm/s can travel half the diameter of the arena. The diffusion length over the decay time scale is ≈ 3mm which may appear small, however, RAnts are not always moving at base speed but often located in a particular location for multiple seconds to even minutes. The parameter choice described here has shown to neither saturate the domain with photormone nor be too volatile, but allowing the photormone to act as a spatiotemporal memory for the RAnts over the course of an experiment.
The positions of the RAnts are tracked with a webcam mounted above the arena and evaluated in Matlab. Blue markers are attached on the centroid of the case's upper surface which allow to use a simple blob detection to identify the pixel position of the RAnts. The photormone concentration is then dynamically updated in the same Matlab script and displayed on the RAnt arena with the projector. The tracking and integration of the photormone field is executed in real time which restricted the update rate of the projected field to 8 Hz on average. The low refresh rate did not have any noticeable consequences for the conducted experiments but might have affected results for RAnts with a much larger base speed and more volatile photormone dynamics.
The set-up of the enclosure for the RAnts consisted of approximately 200 wall elements arranged in three concentric circles where the outermost circle had a diameter of 50 cm . The outermost circle was prevented from being pushed outward from their initial position by a thin plastic ring that was attached to the base of the arena. The plastic ring was thick enough to prevent wall elements from leaving the confinement, but thin enough for RAnts to roll over it to escape the arena. For every experiment we randomly placed the rants in the arena and waited for the first RAnt to excavate out or the time limit of 15 minutes to be reached. At that point, data was stored and the experiment ended. Most experiments required no intervention, but in case of an empty battery of a RAnt or any unexpected critical failure during the experiments, we had two RAnts standing by to replace the defective RAnt. Since all RAnts are identical and the main memory is communicated through the environment and the photormone concentration, a switch has no further statistical consequences on the outcome of the experiments. There was no leader and no dedicated roles, which makes every RAnt replaceable.
We conducted experiments for five cooperation parameters C = {0, 0.25, 0.5, 0.75, 1} at fixed excavation rate E = 1 and repeated experiments five times for each parameter. Every RAnt's software was updated before a new set of five experiments with the same cooperation parameter was conducted. For every experiment, we stored the webcam data and time stamps. The video frames were post-processed and locations of all RAnts and wall elements were stored as a function of time.
For the phase diagram experiments we used the previous data for cooperation parameters C = 0 and C = 1 for partial tunneling and tunneling, respectively. To induce jamming behavior and diffusion behavior the excavation rate had to be changed in the internal programming of the RAnts. By setting the excavation rate E = 0 the probability of the magnet engaging vanished which led to jamming for high cooperation parameters, and diffusion for low cooperation parameters. We only collected data for two trials of a few minutes each in the diffusion and jamming case as tunneling cannot be initiated with disengaged magnets which reduces the timescales over which the behavior occurs.

Cooperation parameter
We explored the effect of cooperation parameters on the excavation time and excavation performance as stated in the main text. Five cooperation parameters, i.e. C = {0, 0.25, 0.5, 0.75, 1} , were selected each of which was tested in five RAnt experiments with five RAnts. Appendix 3- figure 3 shows the final wall element distribution and the RAnt density averaged over the whole trial for all the conducted 25 experiments.
From the final wall element distribution one can deduce the degree of focus during the tunneling effort. For low cooperation parameters the initial three layers are excavated at multiple excavation sites. As the cooperation parameter increases, less excavation sites are visible and at C = 1 there is in general only one large excavation site. While the final wall distribution shows only a snapshot in time, the RAnt distribution is averaged over time and therefore displays where the RAnts were mostly located throughout the run. At low cooperation numbers, the RAnt density is generally distributed all across the arena. Localization of the density toward one region was observed for low cooperation parameters as excavated wall elements were forming a new boundary that confined the RAnt motion to that region (see e.g. C = 0 T4, C = 0.25 T4). As the cooperation number increases, more distinct localized density becomes apparent. Due to the photormone field Rants operating at higher cooperation parameter values are more likely to start excavating in locations where RAnts have previously been present. The location of that attracting field is not known a priori, but emerges spontaneously through the interaction with other RAnts. The location of the peak density field at higher cooperation numbers strongly correlates with the point of excavation in the wall. The difference in RAnt behavior as a function of the cooperation parameter is the degree of focus during excavation as represented by the von Mises parameter of the angular position of excavated boundary elements shown in Appendix 3- figure  2(a). A large value of the parameter indicates a high degree of concentration of the excavation effort, while low values indicate a scattered distribution of many digging sites. Another metric to assess the behavioral difference induced by the cooperation parameter is the traveled distance of the RAnts. Appendix 3-figure 2(b) displays the total travelled distance of a rant x normalized by the arena diameter D as a function of the normalized time t/ts , where ts = D/v b and v b the base speed, shows that RAnts travel a greater distance in the same amount of time at higher concentration parameters. The theoretical limit of the travelled distance is shown with the dashed line in the left-hand side figure, reflecting that RAnts do not constantly move at base speed, but are interrupted by other RAnts, obstacles, and fetching/dropping manoeuvres. As shown in Appendix 3-figure 2(c), RAnts travel at about a fifth of the base speed on average. An increase of the average speed is observed as a function of the cooperation parameter, which can be explained by the fact that obstacles are more scattered at lower cooperation parameters, effectively reducing the mean free path of a RAnt. A single RAnt can efficiently excavate a site if an initial photormone seed is present, but it is not robust. In fact, even though the RAnt in T2 managed to remove some elements in the last layer, it never excavated out but lost the photormone seed where it was digging and started diffusing again. Three rants were more successful in generating an initial photormone seed, but excavation occured at multiple sites even for C = 1 since the lower number of RAnts did not generate one dominating photormone field.

Phases of cooperation in RAnts
In the RAnt case we can infer the phase in which the RAnts operate by looking at the tunnel size, 1/K b (t) and the location along the boundary at which the RAnts are localized, ϱr(ϕ) as they execute their task. We find that in the jammed and diffused phase there exists no tunnel and the variance remains zero throughout the process. However the location along the boundary ϕ b at which the RAnts spend their time the most has a large peak around the jammed location due to high cooperation which in the case of diffusion remains widespread (see Appendix 3- figure 5).

Tunneling
Dif used Partial Tunneling f Appendix 3-figure 5. From the RAnt experiments in Figure 8, K b is the von Mises concentration parameter computed from the location of the boundary and ϱr is the angular distribution of the RAnts in the arena averaged over time. ϕ b is the time-averaged mean azimuthal location of the RAnts in the arena. RAnts present over longer periods in a particular sector of the arena will cause a peak in ϱr . One can infer the phase the RAnts are in by measuring these two quantities.
A successful tunnel, as we have already seen, has an initial increase in the variance that plateaus rapidly due to cooperation driven focus at a given location. As the RAnts are localized, focusing on their task, we again see peaks around the location of the tunnel. For a partial tunnel, due to low cooperation, the variance in the tunnel size is large and the location along the boundary the RAnts spend their effort is spread out. Thus the phase the RAnts operate in can be distinguished by using information about the environment, i.e. the tunnel size, in combination with agent dynamics, i.e. their location.