Collective Movement and Collective Information Acquisition With Signaling

We consider a population of mobile agents able to make noisy observations of the environment and communicate their observation by production and comprehension of signals. Individuals try to align their movement direction with their neighbors. Besides, they try to collectively find and travel towards an environmental direction. We show that, when the fraction of informed individuals is small, by increasing the noise in communication, similarly to the Vicsek model, the model shows a discontinuous order-disorder transition with strong finite-size effects. In contrast, for a large fraction of informed individuals, it is possible to go from the ordered phase to the disordered phase without passing any phase transition. The ordered phase is composed of two phases separated by a discontinuous transition. Informed collective motion, in which the population collectively infers the correct environmental direction, occurs for a high fraction of informed individuals. When the fraction of informed individuals is low, the misinformed collective motion, where the population fails to find the environmental direction, becomes stable as well. Besides, we show that an amount of noise in the production of signals is more detrimental for the inference capability of the population and increases temporal fluctuations, the density fluctuations, and the probability of group fragmentation, compared to the same amount of noise in the comprehension.

INTRODUCTION
Many species, from bacteria [1][2][3] and cells [4,5] to insects [6,7], large animals [8][9][10][11] and humans [12] show collective motion, an intriguing phenomena in which individuals in a population move in ordered groups, presumably due to local interactions [13]. Such a collective motion is suggested to endow many advantages, such as avoiding predators [14], or enhancing the information acquisition capability of the population in noisy environments [15,16]. Besides the biological examples, similar phenomena and similar challenges to address these phenomena can arise in the collective motion of artificial agents [17][18][19][20].
Although the collective motion has been subject to intense study [3,13,[21][22][23][24][25][26][27][28][29][30][31], an important point not considered in the relevant literature, is that in many cases, the information individuals reach from others in the population is provided through communication by exchanging signals [5,[32][33][34][35][36]. Besides, in many cases, collectively moving populations try to find and travel to a preferred goal, such as a nutrient source or a migration root [8,12,16,24]. They may do so, while only a fraction of individuals have information about the preferred root and try to lead the group [8,12,16,24]. These considerations raise the important question that how collective motion in a population of individuals who exchange their social information by production and comprehension of signals is formed, and how the noise inherent in the communication system affects the collective information acquisition capability of the moving population?
To answer this question, we consider a population of mobile agents trying to collectively find and travel towards an environmental direction. Individuals communicate information about their movement by production and comprehension of signals. Individuals can be uninformed or informed. Uninformed individuals make decisions over their direction of motion based on the information they receive from their neighboring individuals. Informed individuals, on the other hand, can make noisy observations of an environmental direction and make decisions over their direction of motion by combining this information with the information they receive from the direction of motion of their neighboring individuals. We show that an ordered phase, where individuals move towards the same direction, emerges for a low level of noise in communication, and a disordered phase, where individuals travel to different directions, emerges for high levels of noise. The ordered phase is composed of two phases. Informed collective motion, in which the population collectively moves towards the environmental direction, emerges if the noise in communication is low enough and the fraction of informed individuals is high enough. As the fraction of informed individuals decreases (for small communication noise), the model shows a discontinuous phase transition to a misinformed collective motion phase in which the population moves collectively, but towards a direction other than the environmental direction. The fraction of informed individuals needed for the population for collectively travel to the environmental direction is larger for higher noise levels.
By increasing the noise in communication, for a low fraction of informed individuals, the model shows a discontinuous phase transition to the disordered phase in which individuals move towards random directions. As is the case in Vicsek-like models [21][22][23], the order-disorder transition suffers from strong finitesize effects that make its discontinuous nature apparent only in large system sizes. On the other hand, when the fraction of informed individuals is high, such that the system is well into the informed collective motion phase, by increasing the communication noise, contrary to the Viscek-like models, the population moves gradually from the ordered phase to the disordered phase without any phase transition. This shows how the amount of information about the environment contained in the population can change the nature of the order-disorder transition observed in Vicsek-like models.
Finally, we show an amount of noise in the production of signals is more detrimental for the collective information acquisition capability of the population compared to the same amount of noise in the comprehension of signals. This result is in keeping with a recently found comprehension-production asymmetry in a model of collective decision making, where individuals residing on a network try to form a consensus about an environment that can be found in a finite number of possible states [37]. While Ref. [42] provides theoretical evidence for a comprehension-production asymmetry in a case where a population of immobile individuals has access to a Potts variable as their belief variable, our finding here provides evidence for a similar asymmetry in the case of collective motion, where agents can move on space and need to make decisions over their direction of motion. Besides, we show that the production noise increases the density fluctuations and the probability of group fragmentation, and temporal fluctuations compared to comprehension noise. Thus, in this regard, by extending results found in Ref. [37], our findings suggest asymmetry between signal comprehension and production is a fundamental characteristic of biological communication systems.

THE MODEL
A schematic representation of the model is provided in Figure 1.
We consider a population of N mobile agents, moving with constant speed v 0 on a L × L two-dimensional surface with periodic boundaries. There is a favorable direction of motion ε, called the environmental direction. As a direction in two dimensions can be repressed with an angel, ε is an angle in a two-dimensional surface. In this manuscript, we measure angles with respect to the x-axis. Individuals try to find and travel towards ε. Individuals decide about their direction of motion based on their personal observation (if available) and social information. To model observation, each individual is equipped with a noisy information channel R, and to communicate, individuals are equipped with a signal production G, and signal comprehension channel C, which they use for exchanging signals. Below we describe these channels in turn.
Observation. Personal information is acquired by noisy observation of the environmental direction. We assume observation is made through a noisy channel R(r|ϵ), such that the result of an observation in environment ε is in the interval [r − dr/2, r + dr/2] with probability R(r|ϵ)dr. We take R(r|ϵ) to be uniformly distributed in the interval [ϵ − η R , ϵ + η R ]. η R can be thought of as the noise level in observation. To implement this, we simply set the result of an observation made by an individual according to r ϵ + ξ R , where ξ R is a uniformly distributed noise term drawn from the interval [−η R , η R ]. We assume only a fraction h of the individuals, called informed individuals, are able to observe the environment.
Communication. Social information is acquired by production and comprehension of signals. To communicate its direction of motion b, each individual produces a signal σ, using its signal production channel G(σ|r b). That is a signal in the interval [σ − dσ/2, σ + dσ/2] is produced when an individual intends to signal direction b with probability G(σ|r b)dσ. Here, dσ is a differential element, and b and σ, referring to a direction in two dimensions, satisfy b, σ ∈ [0, 2π). We note that both b and σ are angles in two dimensions. We take G(σ|b) to be uniformly distributed in the interval [b − η G , b + η G ]. η G is a measure of noise in signal production. This can be implemented by setting σ b + ξ G , where ξ G is a uniformly distributed noise term drawn from the interval [−η G , η G ].
Signals are transmitted up to a distance l. That is all the individuals in a circle of radius l centered around the transmitter receive the signal. The receivers, comprehend a signal to refer to a direction of travel r ′ ∈ [0, 2π), according to their comprehension channel C(r ′ |σ). That is, signal σ is comprehended as referring to a direction in the interval [r ′ − dr ′ /2, r ′ + dr ′ /2] with probability C(r ′ |σ)dr ′ . We take C(r ′ |σ) to be uniformly distributed in the interval [σ − η C , σ + η C ]. η C is a measure of noise in signal comprehension. This can be implemented by setting r ′ σ + ξ C , where ξ C is a uniformly distributed noise term drawn from the interval [−η C , η C ]. In this way, we have taken the fact that signal production and comprehension are in general subject to noise [37][38][39][40][41][42] into account. This will allow us to study the effect of noise in communication on collective motion and collective information acquisition of the population.
Decision making. As a result of signals an individual receives, it reaches a set of representations r ′ . We note that bold letters show a set, and thus, r ′ is a set. This set is composed of all the directions an individual receives by comprehending signals in its l-neighborhood. Besides, informed individuals make a personal observation, r. Each informed individual, α, makes a decision about its direction of motion b α , based on its observation r α and social information r ′ α , by a weighted averaging rule. That is: Here, r ′ α is the number of representations (directions) individual α has received, and ω is the self-confidence of individuals. Thus, is the average direction of the motion of neighboring individuals of the individual α. ω lies between 0 and one and determines how individuals weigh their personal observation compared to their social information in decision making. We note that, here and in the following, to calculate summation over angles it is necessary to take periodicity into account. For this purpose, to average a set of angles (θ 1 , .., θ n ), we first decompose the angles into their x and y components, (x θ1 , .., x θn ) and (y θ1 , .., y θn ), and define the average angle based on its X and Y-components, X n i 1 x θi and Y n i 1 y θi , as The same care is taken in calculating weighted average between two angles as in Eq. 1.
An uninformed individual, β, makes a decision about its direction of motion by simply averaging the social information it receives: As in the Vicsek model [21], the dynamic is synchronous. At each time step, informed individuals make an observation using R(r|ϵ), and all the individuals transmit a signal σ based on their direction of motion b using G(σ|b). The signals are received by all the individuals up to a distance l of the transmitter. Receivers, comprehend the signals as referring to direction r ′ using C(r ′ |σ). Finally, individuals make a decision about their direction of travel using their decision-making rule and update their direction accordingly. The simulations start with a random distribution of the position and direction of individuals. The base parameter values used in the simulations (unless otherwise stated) are l 1, In the following we use two variables to distinguish different phases of the system. The absolute value of the normalized average velocity of the population is defined as m This variable is commonly used in the Viscek-like models to distinguish ordered from the disordered motion. As a second variable we define the angular deviation of the average population direction from the environmental direction as δθ N α 1 θα−ϵ N . We will use this variable to distinguish informed collective motion Each agent is equipped with two information channels to communicate. Production channel G(σ|b) is used for signal production, such that signal σ is produced for a direction b with probability G(σ|b). In the same way, the Comprehension channel C(r′|σ), is used for comprehension of the signals. In addition to these, a fraction q of the individuals are informed individuals, who are equipped with an observation channel R(r|ϵ). (B) and (C): The dynamics of the model. Agents move on a two-dimensional space with constant speed v 0 and towards different directions b i . In each time step, each agent i produces a signal σ i based on its direction b i to communicate its direction of motion (B). The signal travels up to a distance l, such that all the agents up to a distance l from a transmitter receive its signal. Receivers comprehend signals using their comprehension channel (C). Here, agents one and three receive a signal only from agent 2, but agent 2, having two individuals in its l-neighborhood, receive signals from both one and3. Besides, agent 3, who is an informed individual, makes an observation using its observation channel. Each agent i, makes a decision based on the information they receive from communication r′ i , and the result of observation r i (in the case of informed individuals), using the weighted averaging rule Eq. 1. Here, r′ from misinformed collective motion. Informed collective motion and misinformed collective motion can also be distinguished using average angular deviation of the direction of motion of individuals from the environmental direction, defined as δθ ′ N α 1 |θα−ϵ| N . As we will see in the following δθ and δθ ′ lead to similar pictures. For a population with observation noise, production noise, and comprehension noise, equal respectively to, η R , η G and η C , and with a fraction h of informed individuals, we show these variables by,

Collective Motion Phases
We begin by studying the behavior of the system as a function of the noise level and the fraction of informed individuals. To do this, in Figure 2A, we plot the absolute value of the normalized average velocity of the population in the case that an amount of noise η is in the comprehension and production is noiseless, m(0.75, 0, η, h). The case that the noise is in production is qualitatively similar. For low noise level, m is close to 1. This shows that collective motion emerges, and individuals travel in the same direction. By increasing the noise level, the model shows a transition to a disordered phase in which individuals fail to align their motion and travel in random directions. Thus, m takes a small value.
Interestingly, the ordered phase is composed of two distinct phases. To see this, in Figure 2B, we plot the angular deviation of the average direction of motion from the environmental direction, δθ(0.75, 0, η, h). In the ordered phase, for large h, the population possesses enough informed individuals to be able to collectively find the environmental direction. Consequently, δθ(0.75, 0, η, h) takes a value close to zero. We call this phase informed collective motion phase. On the other hand, for small h, misinformed collective motion, in which the population moves collectively in a wrong direction, becomes possible as well, and the system becomes bistable (see below). Consequently, δθ(0.75, 0, η C , h) takes a large value. We call this phase misinformed collective motion phase.
In 2 (c), we plot average angular deviation of the direction of motion of individuals from the environmental direction δθ ′ (0.75, 0, η, h). This measure shows similar behavior to that of δΘ for small h. However, the two measures differ for large h and large noise levels. The reason is that, for large noise levels, ordered motion does not emerge. However, as long as h is large, individuals travel on average towards the environmental direction. Although, due to noise in observation, a large fluctuation around the environmental direction exists. For δΘ, which gives the average direction of motion of the population, these fluctuations, being distributed symmetrically around ε cancel. On the other hand, in δΘ ′ , which gives the average deviation of the individuals' direction of motion from the environmental direction, these fluctuations show up. Consequently, the two measures show different behavior for large h and large η.

Phase Transitions
The nature of the order-disorder transition in the Vicsek and related models has been the subject of intense research [13,[21][22][23]. It is well known that the order-disorder transition is discontinuous in many cases (depending on the parameter values and the way noise is implemented in the model) [13]. However, due to strong finite-size effects, the discontinuous nature of the transition shows up only in very large sizes [22]. Our model shows a similar phenomenology only for small h. However, the situation is different for large h. This can be seen in Figure 3A and Figure 3B for respectively, the production and comprehension noise, where the probability distribution of the order parameter for different noise levels and constant (large) h is plotted. Here, the distribution is derived from a single time series of the system of length T 50000, after discarding the first 1,000 time steps. As can be seen, by increasing the noise level, the order parameter gradually decreases. This excludes a discontinuous transition. Besides, the distribution remains peaked at a single value and no broadening of the distribution resulting from large fluctuations characteristic of a continuous transition occurs [43]. This suggests it is possible to go from the informed collective motion phase to the disordered phase without passing any phase transition. In Figure 3C and Figure 3D, by plotting the distribution of the order parameter for fixed h and different noise levels, we study the discontinuous nature of the order-disorder transition for small h.
In Figure 3C the case of production noise is considered. Here, L 100, v 0 0.6, ω 0.25, η R 0.75 and the density ρ 0.1. The distributions are derived from a single time series of the system of size T 1.2 × 10 6 , after discarding the first 10 5 time steps. As can be seen, the distribution shows distinct peaks corresponding to the ordered and disordered phases. As the noise level increases, the peaks corresponding to the ordered phase decrease, while that corresponding to the disordered phase increases. This phenomenology is characteristic of a discontinuous transition [22,44]. In Figure 3D, the case of comprehension noise is considered, where the same bimodality which suggests a discontinuous transition is observed. Here, L 200, v 0 0.1, ω 0.25, η R 0.75 and ρ 0.1. We note that the finite-size effects are much stronger for comprehension noise compared to production noise, such that the discontinuous nature of the transition becomes apparent for larger population size in the former. Besides, for production noise, the finite-size effects are stronger for smaller velocities, while for comprehension noise, they are stronger for larger velocities. We have not been able to conclusively infer the discontinuous nature of the transition for large velocities when noise is in the comprehension. Figure 3C, we see that for the case of production noise the distribution of the order parameter has two major peaks in the ordered phase. This happens only for small enough h, i.e., when the ordered phase is bistable, and the system can be found in both informed and misinformed collective motion phases. The two major peaks of m correspond to these two phases. This can be seen in Figure 3E, where m together with δΘ, as a function of time, are plotted. As can be seen, m shows intermittency between two values corresponding to the two peaks of the ordered phase in Figure 3C. When m is very large, corresponding to the rightmost peak in Figure 3C, δΘ is very small, indicating that the average direction of motion coincides with the environmental direction. On the other hand, when m takes the smaller value, δΘ becomes large, indicating the average population direction differs from the environmental direction. The reason why the value of m in the misinformed collective motion phase is smaller than that in the informed collective motion phase is that when noise is in production (as is the case here), the probability of group fragmentation is high, such that in many times, the population is composed of different groups each collectively moving towards a different direction independently of others. Occasionally a group is decomposed into smaller groups, or different groups can merge to form a larger group (see the Supplementary Video for a visual manifestation). In the informed collective motion phase, all the groups head towards the  Figure 3C suggests. For comprehension noise, m does not show a similar bimodality that indicates intermittency between informed and misinformed collective motion. The reason is that, contrary to the production noise, with comprehension noise, the population rarely is decomposed into different dense groups with different directions of travel, and the strong fission-fusion dynamics observed for production noise is absent in the case of comprehension noise (see the Supplementary Video). We will shortly return to this difference between comprehension and production noise.

Returning to
The bi-stability associated with a discontinuous transition results in hysteresis, which provides an alternative way to test the nature of the informed-misinformed phase transition [44]. This is shown in Figure 3F, where the hysteresis loop for the case of comprehension noise is shown. Production noise shows similar hysteresis effects. Here, we run a simulation beginning with h 6 400 , which lies in the informed consensus phase. We gradually decrease h down to h − 6 400 (a negative h results from reversing the environmental direction) and then increase it back to the initial value. The resulting hysteresis loop results from the memory effects and indicates a discontinuous transition.

Noise in Signal Production Increases the Density Fluctuations, Temporal Fluctuations, and Degrades Collective Information Acquisition Capability
We have already seen that production noise increases the probability of group fragmentation and leads to a strong fission-fusion dynamic absent for the comprehension noise. This can be shown more quantitatively. For this purpose, we define the (relative) asymmetry in density fluctuations as δρ(η R ,η,0,h)−δρ(η R ,0,η,h) δρ(η R ,η,0,h)+δρ(η R ,0,η,h , where, the density fluctuation δρ, is defined as the standard deviation of spatial density of the population. In a situation where the individuals are distributed uniformly in space, which is often the case for comprehension noise, this takes a small value. On the other hand, when the group is decomposed into independently traveling dense groups, as is often the case for the production noise, this takes a large value. Consequently, the asymmetry in density fluctuations is always positive, as can be seen in Figure 4A. This shows, compared to comprehension noise, the production noise increases the density fluctuations and the probability of group fragmentation. There is another asymmetry between the comprehension and production of signals. This asymmetry arises in the collective information acquisition capability of the population. To see this, we define the asymmetry in the inference capability of the population as the difference between the angular deviation of the average population direction from the environmental direction when an amount of noise is in signal production compared to the case when the noise is in signal comprehension, δΘ(η R , η, 0, h) − δΘ(η R , 0, η, h). To justify this choice, we note that in the case that production noise is more detrimental for the collective information acquisition of the population, it leads to a higher angular deviation compared to the case that the same amount of noise is in the comprehension. Consequently, this quantity becomes positive. The asymmetry in the inference capability of the population, δΘ(η R , η, 0, h) − δΘ(η R , 0, η, h), for η R 0.75, in the η − h plane is plotted in Figure 4B.b. As can be seen, this quantity is always positive. This indicates that, compared to the same amount of error in the comprehension, an amount of error in production leads to a poor collective inference, and thus, a larger angular deviation from the environmental direction.
As a second measure of asymmetry in the inference capability of the population, we consider the difference between the average deviation of the individuals' direction of motion from the environmental direction, when an amount of noise is of production type compared to the case where the noise is of comprehension type δΘ(η R , η, 0, h) − δΘ(η R , 0, η, h). This is plotted in Figure 4C in the η − h plane. As can be, this We note that the positivity of the measures of asymmetry in collective information acquisition results from two shifts in the phase transitions. First, production noise shifts the order-disorder transition to smaller noise levels. This shows production noise is more detrimental to the ordering of the population. Second, production noise shifts the monostable informed collective motion phase to larger values of h. This means with production noise, a larger fraction of informed individuals is necessary for the population to successfully infer the correct environmental direction. As shown in the Supplementary Material, both asymmetries in the density fluctuation and collective information acquisition are robust features of the model, valid for all the parameter values.
Finally, we note that noise in signal production also increases temporal fluctuations. To see this, in Figures 5A,B, we plot the normalized average velocity, m, and average angular deviation of the individuals from the environmental direction, δΘ ′ , for, respectively, the cases where an amount of noise, η, is in signals production and signal comprehension. As can be seen, the same amount of noise is in signal production leads to higher temporal fluctuation in both m and δΘ ′ . We note that the same picture holds for δΘ. To see production noise leads to higher temporal fluctuation in the entire phase diagram, we calculate the standard deviation of the time series of m and δΘ ′ , denoted respectively as Σ m and Σ δΘ′ . Subtracting the fluctuations in m when an amount of noise is in signal production from that when the same amount of noise is in signal comprehension, Σ m (0.75, η, 0, h) − Σ m (0.75, 0, η, h), we arrive at a measure of asymmetry in temporal fluctuations when an amount of noise is of production type compared to a case where the same amount of noise is of comprehension type. Similarly, we define asymmetry in the temporal fluctuations of the average angular deviation of the direction of motion of individuals from the environmental direction, by first calculating the standard deviation of the time series of δΘ ′ (0.75, η, 0, h) and δΘ ′ (0.75, 0, η, h), denoted respectively as Σ δΘ (0.75, η, 0, h) and Σ δΘ (0.75, 0, η, h).
Subtracting the first term from the second one, we arrive at a measure of asymmetry in the temporal fluctuation of the collective information acquisition capability when an amount of noise is of production type from that when the same amount of noise is of comprehension type. These two measures are plotted in Figures  5C,D in η − h plane. As can be seen, both measures are nonnegative in the entire η − h plane. This shows that an amount of noise in signal production increases temporal fluctuations compared to the same amount of noise in signal comprehension.

DISCUSSION
We have introduced a model of collective movement in which individuals in a population try to collectively find and travel in a preferred direction. To do so, a fraction of individuals, informed individuals, can use private information provided by the noisy observation of the environment, and social information, provided by communication between individuals by exchanging signals. Others, uninformed individuals, can only use social information provided by communication. By analysis of the model, we showed that signal production noise not only decreases the inference capability of the population but also increases density fluctuations and the probability of group fragmentation. Besides, by identifying two different phases of collective motion (informed and misinformed) separated by a discontinuous transition, we have identified the mechanism by which a fraction of informed individuals able to only make noisy observations of the environment can lead the group. Finally, we have cast the muchstudied order-disorder transition in Vicsek-like models into a broader context by showing that the nature of this transition depends on the amount of information about the environment available in the population, in other words, the fraction of informed individuals and the accuracy of their observation.
A similar comprehension-production asymmetry in the inference capability and similar phase transitions had been recently observed in a model of collective decision making in a population of immobile agents [37]. In this model, a population Frontiers in Physics | www.frontiersin.org of agents who reside on a fixed communication network and live in an environment that can take one out of a finite number of possible states, try to collectively infer the environmental state [37,45]. Similarly to the environmental state, in this model, individuals can take one out of a finite number of beliefs and are equipped with a production and a comprehension probability transition matrices to produce and comprehend a set of a finite number of signals. It is shown that such a model of collective decision making with signaling shows a similar comprehension-production asymmetry, according to which an amount of noise in signal production is more detrimental for the collective capability of the population to infer the environmental state, compared to the same amount of noise in signal production [37]. Our work here extends this study in different ways. In previous models of biological communication by production and comprehension of signals [37,42,45], or more generally, models of a proto-language [46][47][48][49], individuals try to communicate a finite number of states, which lack a distance measure. In other words, in that case, the state or belief variable is a categorical variable. In contrast, here, we have explored a case when individuals can form beliefs over a continuous space and produce and comprehend a continuous set of signals, which possesses a distance measure i.e., an ordinal variable. In this regard, in the limit of zero movement speed, our model reduces to a model of biological communication in a population of immobile agents, who can form beliefs and produce and comprehend signals over a continuous space of possible beliefs which possesses a distance measure. Second, by introducing movement into the model, we have been able to study collective movement in a communicating population. Our finding here shows that a similar asymmetry in the collective inference capability of biological populations is at work in this case as well. Furthermore, by showing that the production noise also increases temporal fluctuations, the density fluctuations, and the probability of group fragmentation, our analysis reveals new ways in which noise in signal production can be more detrimental than noise in the comprehension. This theoretical study thus strengthens the case for the existence of a comprehension production asymmetry and calls for empirical investigations to shed light on the question that whether such an asymmetry exists in biological populations.
Our study also provides insight on the question that how a fraction of informed individuals can lead the group on the move. Previous work has shown that enough fraction of informed individuals can indeed lead the collective movement of the group [24,50]. Our study extends this line of research in two ways. First, our model introduces the more realistic possibility that informed individuals have only noisy and imperfect information of a preferred goal, such as a nutrient source. Second, our model studies how different types of noise in communication can affect the ability of the informed individuals to lead the group and the group's ability to find the preferred direction of motion. In this regard, our finding reveals that depending on the fraction of informed individuals and noise in observation and communication, a collectively moving population can be found in different ordered or disordered phases, separated by phase transitions. Besides, our analysis reveals that the fraction of informed individuals needed to successfully lead the group increases by increasing noise in communication. This can be interpreted to result from the fact that a higher noise level in communication disrupts the flow of information in the population, and thus, a higher fraction of informed individuals and a higher information flow from the environment to the population through informed individuals is necessary for finding the environmental direction.
Finally, we note that the similarity between the phenomenology of the two apparently different systems, that is, the model of collective movement introduced here and the model of collective decision making mentioned before [37,45,49], can be understood by noting that movement decision involves choosing a direction between a continuous set of possible directions. This observation highlights the similarity of collective movement with collective decision making, which can provide insights into the physics of collective movements such as the nature of the phase transitions observed in such models or the mechanisms by which such populations can optimize their information acquisition capabilities [45], which can be subject to future researches.

The Simulations
All the simulations are started with a random distribution of positions and directions of motions. The parameter values used in the simulations differ for each figure and are given in the figure captions. To derive the hysteresis loop in Figure 3.f, we consider a population of size N 400, with a fraction of informed individuals equal to h 6 400 (other parameter values are presented in the figure caption). Starting from random distribution of positions and directions of motion, the simulation is run for T 10000 time steps, after which h is decreased by a value 1 400 , and the procedure is repeated until h − 1 400 is reached. As mentioned before, a negative h results from reversing the environmental direction (here ϵ π/2 is changed to ϵ −π/2). Then h is increased again in steps of δh 1 400 , each time after performing the simulation for T 10000 time steps, until h reaches its initial value of h 6 400 . The angular deviation and its error bars for each value of h are calculated based on the last T 2000 steps (after the system equilibrates) of the simulation for each h value.

Density Fluctuations
To calculate the density fluctuation, we first define the density field, ρ(x, y) as the number of individuals per unit area. To define this quantity, we divide the space into square bins of linear size l 0 1. The density fluctuation is then defined as the standard deviation of this quantity δρ [〈[ρ(x, y) − 〈ρ(x, y)〉] 2 〉] 1/2 . Where averages, 〈.〉, indicate an average over space.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.


Characterizing The Limits of Linear Modeling of Non-Linear Swarm Behaviors

We study the limits of linear modeling of swarm behavior by characterizing the inflection point beyond which linear models of swarm collective behavior break down. The problem we consider is a central place object gathering task. We design a linear model which strives to capture the underlying dynamics of object gathering in robot swarms from first principles, rather than extensively relying on post-hoc model fitting. We evaluate our model with swarms of up to 8,000 robots in simulation, demonstrating that it accurately captures underlying swarm behavioral dynamics when the swarm can be approximated using the mean-field model, and when it cannot, and finite-size effects are present. We further apply our model to swarms exhibiting non-linear behaviors, and show that it still provides accurate predictions in some scenarios, thereby establishing better practical limits on linear modeling of swarm behaviors.

Introduction
Swarm Robotics (SR) is the study of the coordination of large numbers of simple robots (Ş ahin, 2005). SR systems can be homogeneous, i.e., single robot type and identical control software, or heterogeneous, i.e., multiple robot types and control software (Prorok et al, 2017;Ramachandran et al, 2020). The main differentiating factors between SR systems and multi-agent robotics systems stem from the mechanisms on which SR systems are based. Historically, these were principles of biological mimicry or problem solving techniques inspired by natural systems of agents such as bees, ants, and termites (Labella and Dorigo, 2006), though modern SR systems typically incorporate more mathematically rigorous elements of conventional multi-robot system design (Castello et al, 2016).
In this work, we study the central place foraging problem in which robots gather objects (blocks) across a finite operating arena and bring them to a central location (nest) under various conditions and constraints. Foraging is one of the most studied applications of SR, due to its straightforward mapping to real-world applications (Hecker and Moses, 2015); for an extensive discussion of the state of the art, see (Lu et al, 2020). The complexity of the foraging task frequently gives rise to non-linear behavioral dynamics caused by inter-robot interactions.
It has been established that above a certain swarm density ρ S (ratio of swarm size to arena size) (Sugawara and Sano, 1997;Hamann, 2013), interactions between robots can produce non-linear behaviors which are not predictable from the individual components, e.g., from the swarm control algorithm (Cotsaftis, 2009;Georgé and Gleizes, 2005;Hunt, 2020;De Wolf and Holvoet, 2005). Hence, the collective performance of a swarm S of N robots each running an identical control algorithm κ can be a non-linear function of the behavior of a system of N independent robots (Harwell et al, 2020). Below this density threshold, swarm behavior can be well approximated using linear models.
Many real problems across scales can be tackled using swarms in low density environments independently of swarm size, in which each robot is responsible for a large area on the order of 100m 2 ; these include indoor warehouses (≈ 64m 2 ), outdoor search and rescue, precision agriculture, and field monitoring (600,000m 2 ≈ 150 acres). Recent results summarizing the challenges of moving swarms into the real world argue that directing research towards low density swarm applications is critical (Tarapore et al, 2020). Some researchers hold that swarms with low ρ S are not properly characterized as swarms, and are instead systems of independent robots because they lack the high level of inter-agent interaction which characterizes natural swarms. We argue that a high level of inter-agent interaction, and therefore potentially emergent self-organization, is only one of the defining properties of swarms. Low density swarms exhibiting the properties of scalability, flexibility, and robustness can properly still be considered swarms (Harwell et al, 2020). The exact value of ρ S at which a given linear model of swarm behavior breaks down is influenced by many factors, including the control algorithm κ, the number of robots N, and characteristics of the problem being solved, so in general it cannot be determined a priori.
In this work we seek to characterize, as a function of ρ S , the practical limits of using linear models to model non-linear behaviors of a swarm engaged in real-world foraging. Specifically, we will show that non-linear behavioral dynamics in the foraging problem can be captured using linear modeling which is derived directly from problem features, such as swarm density, swarm size, and robot control algorithm.

Motivation and related work
Many SR systems have been designed around imitating natural systems, and use heuristic decision making (Castello et al, 2016) rather than combining natural principles with a strong mathematical grounding (Talamali et al, 2020). Nevertheless, heuristic approaches to swarm control have been effective for robots that operate with incomplete information and limited computing power. SR researchers taking this approach average large numbers of simulation runs to develop more accurate models of swarm behavior and obtain empirical insights into real-world problems (Harwell and Gini, 2019). The emphasis on empirical rather than mathematical models has been a chief impediment to a wider use of SR systems (Lerman et al, 2004). Systematically varying individual agent parameters to study their effect on collective swarm behavior is impractical, even in simulation. Mathematical characterization of collective swarm behavior is more difficult, but provides the means to precisely predict it a priori-without the need of repetitive experiments.
Given the complexity of SR systems, and the frequently non-linear ways in which behaviors can arise, it is difficult to obtain precise bounds on the collective behavior of a swarm S of N robots each running a control algorithm κ. Robots might need to respond to environmental signals that arrive at unpredictable times; such systems are well-modeled as asynchronous, and therefore difficult to predict precisely. However, if we conceptualize S as a differentiable, continuous quantity, its dynamics can be modeled with Ordinary Differential Equations (ODEs) whose variables are the population counts associated with the different roles. We can apply a macroscopiccontinuous ODE modeling approach for the average behavior of S in the steady state (Berman et al, 2007), with the caveat that when using such a model to determine the behavior of S, actual system behaviors could be far from the average (Lerman et al, 2004). Usually, the larger the system, the smaller the fluctuations; in small systems the fluctuations resulting from finite size effects can be of order N, resulting in models which are accurate at asymptotically large scales but not at small scales. The master equation (Kampen, 2007), which is typically used to model the expected average behavior of systems, can be used to calculate the deviation from the average, but such calculations are often intractable or algebraically difficult.
A promising mathematically rigorous methodology utilizing macro-and microscopic models for group dynamics and individual behavior over time has been developed (Lerman and Galstyan, 2002;Berman et al, 2007;Galstyan et al, 2005;Sugawara and Sano, 1997), which sidesteps the difficulty in modeling the average behavior of the swarm, by instead modeling the change in the average behavior of the swarm, which is much easier. It uses (1) differential equations to model the behavior of the average number of robots in S in a given state, (2) discrete difference equations to model the stochastic transitions between robot states, and (3) stochastic simulation of discrete difference equations to compute state transition rates for all robots. It draws on implementations of the stochastic master equation in chemistry and statistical physics (Kampen, 2007). Through the usage of rate constants and population fractions in each state, it is possible to assess the general behavior of a system under a variety of stochastic circumstances. Most importantly, this approach has predictive control and performance guarantees-crucial components to translating laboratory models into viable real-world solutions without needing behavioral characterization of robot swarms via simulation experiments. Models in this differential equation paradigm operate on both the forward problem, i.e., predicting collective behavior from features of the control algorithm each robot runs (Lerman and Galstyan, 2002), and the inverse problem, i.e., incorporating design constraints into algorithm design in order to produce a desired collective behavior (Berman et al, 2007;Hsieh et al, 2008). However, the models are not overly robust, for two reasons. First, they make large simplifying assumptions such as homogeneous agent distributions, homogeneous environments (e.g., no obstacles and/or a completely visible arena), and Markov/semi-Markov scenarios (Berman et al, 2007). Second, they require many free parameters and extensive post-hoc model fitting which are specific to the implementation of the algorithm under study. Nevertheless, many notable applications of various forms of this methodology have appeared in the literature, demonstrating its practical utility. A few noteworthy examples include the stick pulling experiment (Ijspeert et al, 2001), foraging of green/red pucks using agent memory (Lerman and Galstyan, 2003), the "house hunting model" (Hsieh et al, 2008;Berman et al, 2007), and ant-inspired models that collaborate with or without communication (Sugawara and Sano, 1997).

Background in ODE modeling and diffusion
To model swarm collective behavior using ODEs, we consider the individual robot Finite State Machine (FSM) (Fig. 1). Our FSM is identical to previous work (Lerman et al, 2001;Lerman and Galstyan, 2002). Each state maps directly to a single robot behavior or a set/sequence of robot behaviors which together make up the robot controller for executing a foraging task. It is a coarse-grained model of robot behavior, which omits controller details such as sensing and actuation and contains the minimum number of states needed to describe the system dynamics for the problem. Such an approach is generally more mathematically tractable. If the results of the analysis do not sufficiently agree with the observed collective behavior, other states can be added (Lerman and Galstyan, 2002).  Another benefit of the ODE modeling approach is that each of the states in Fig. 1 maps directly to an ODE describing this change, where the transition rates for states become the ODE terms, and all ODEs can therefore be directly written down from robot control algorithm characteristics. We refer the reader to (Lerman et al, 2001;Kampen, 2007) for the theoretical underpinnings of the validity of this translation, which can be applied to any SR system, not just those doing a foraging task. The quantities modeled in Fig. 1 are listed in Table 1, and the ODE model of these quantities from previous work (Lerman et al, 2001;Lerman and Galstyan, 2002) is summarized in Eqs. (1) to (4). In our model, shown later in Section 4.3, we simplify the original equations and replace some parameters with mathematical derivations. Eq.
(1) describes the change in the number of robots in the searching state, which decreases as searching robots pickup blocks or encounter other robots and switch to collision avoidance.
The rate at which robots leave the searching state and switch to collision avoidance, α r N s (t)[N s (t) + N], can be understood as follows. When a searching robot encounters another searching robot, both switch to the collision avoidance state, decreasing the N s (t) by two; when a searching robot encounters either a homing robot or a robot already in the collision avoidance state, N s (t) decreases by one. The total decrease is then proportional to 2N s (t)+N h (t)+N s av (t)+N h av (t) = N s (t) + N, where N is the total number of robots in the swarm. Finally, searching robots encounter other robots at rate α r , regardless of the state of the other robot. Eq. (1) increases as homing robots deposit blocks in the nest or as searching robots exit the collision avoidance state. Eq.
(2) describes the change in the number of robots in the homing state, which increases as robots acquire and pickup blocks or leave the collision avoidance state, and decreases as robots enter the collision avoidance state or deposit their block in the nest.
can be understood analogously to the one in Eq.
(1). The model assumes that the rate at which homing robots encounter other robots is different than for searching robots, reasoning that there will be more congestion near the nest, and therefore uses a separate parameter α r to account for this. Eq.
(3) describes the change in the number of robots avoiding collision with other robots.
Eq. (4) describes the change in the number of free blocks available for robots to find; it decreases whenever a searching robot acquires a block (i.e., once deposited in the nest, blocks are not re-distributed in the arena).
In addition to using ODEs to model collective behavior, we also draw on diffusion theory. In diffusing systems composed of homogeneous particles undergoing normal diffusion, the average particle displacement is proportional to the diffusion time, i.e., a linear relationship; this assumes the particle is moving in an infinite, structureless medium close to equilibrium. However, there are important circumstances in which the relationship between average particle displacement and time is non-linear, referred to as anomalous diffusion (Oliveira et al, 2019;Metzler et al, 2014;Vlahos et al, 2008). For instance, in biological systems, interactions with other particles or membranes could influence the observed diffusion for macro-proteins as they traverse through biological media (Santamaria-Holek and Vainstein, 2009;Weiss et al, 2003;Nicolau Jr et al, 2007). Anomalous subdiffusion arises due to crowding in a concentrated system, which can make it heterogeneous and disordered (Surya K Ghosh and Metzler, 2016). Crowded environments can obstruct the path or make certain spaces inaccessible to the diffusing particle and prompt the rise of this behavior.

Changes to previous work
The ODE model described in Eqs.
(1) to (4) from (Lerman and Galstyan, 2002) has some limitations which we address in our model.

Uniformity assumptions
In previous work S is assumed to be uniformly distributed in 2D space, and blocks to be scattered randomly; therefore, the model is only accurate on scenarios meeting these criteria. In Fig. 2 we show some of those cases. Random (RN) block distributions (Sugawara and Sano, 1997; Hecker and Moses, 2015) ( Fig. 2d), are appropriate in scenarios such as order fulfillment in a warehouse, but many scenarios cannot be modeled by such distributions (Campo and Dorigo, 2007;Sugawara and Sano, 1997;Castello et al, 2016;Hecker and Moses, 2015;Lerman et al, 2001). For example, transferring material from one side of a building to another location requires a single source (SS) ( . For other scenarios, such as evacuation of civilians from a disaster zone, the block distribution cannot be inferred a priori, and a power law (PL) distribution ( Fig. 2c), in which blocks are clustered in groups of various sizes, is appropriate (Hecker and Moses, 2015;Harwell et al, 2020). By considering the distribution of blocks in a general way using spatial clusters, we guarantee that the model itself describes underlying characteristics of the environment and is robust enough to be applied to a wide range of foraging scenarios.

Invalid steady-state assumptions
The model assumes that S has reached steady-state, but at t = 0 assumes a finite number of objects for the swarm to find (Eq. (4)). This does not model accurately long-running SR systems or those with a small number of blocks relative to N. We formulate our ODE model to handle both steady-state and non steady-state behaviors by considering two types of operating environments. First, those in which any block deposited in the nest is re-distributed immediately in the operating area, so the overall number of blocks available for robots to find in the arena remains the same over time. Second, those in which any block deposited in the nest is not re-distributed, resulting in fewer blocks available for robots to find over time, which is the paradigm used in (Lerman and Galstyan, 2002). Our model treats non steady-state foraging environments as a special case of steady-state behaviors by setting to zero the number of areas where blocks can be distributed after nest deposition (see A j later in Eq. (8)). In this work, we are interested in steady-state behaviors.

Many free parameters
Free parameters do not invalidate the theoretical basis of the model, but limit its reuse by requiring iterative parameter refinement. The need to do posthoc model fitting requires experimental simulations, which makes it difficult to use the model as a predictive tool. We instead derive analytical expressions for all parameters except τ av from scenario characteristics, greatly reducing the need for post-hoc model fitting and reducing the number of free parameters from five to one. τ av is directly computable from the robot controller characteristics, or otherwise obtainable via experiments with a single robot. We introduce C ca (m), to characterize the inter-robot interference in scenario m. Although it is a free parameter, it has a strong semantic context, making fitted values more meaningful and suggesting the possibility of future derivation.

Preliminaries
To make accurate predictions of collective swarm behavior in the steady-state for each of the block distributions shown in Fig. 2, our model makes the following assumptions: 1. The swarm density ρ S is relatively low, so the behavior of N robots can be approximated by a linear function of the behavior of a single robot. An important consequence of this is that we can ignore the effect of a robot avoiding collision and encountering another robot during its avoidance maneuvers. This is not necessarily a limiting assumption, as many real-world applications require sparsely distributed swarms (Tarapore et al, 2020). We relax it as part of our model evaluation experiments in Section 6. 2. N is sufficiently large such that S can be approximated using the mean field model. We relax also this assumption as part of our model evaluation experiments Section 6. 3. S is homogeneous, and the control algorithm used by the robots is a Correlated Random Walk (CRW) (Renshaw and Henderson, 1981;Harwell et al, 2020), which is a random walk where the direction of the next step is biased based on the direction of the previous step. The bias angle θ is drawn from a probability distribution f (θ ) (Codling et al, 2010). Robots do a CRW at an average velocity v s until they acquire a block, which they then transport to the nest using phototaxis, i.e., motion in response to light, at a constant speed v h . Robots have no memory. 4. To reduce congestion in the nest, robots do not return to the exact nest center to drop their carried object, i.e., they employ a reactive avoidance strategy. Instead, they choose a random point along their trajectory from where they enter the nest to treat as the nest center. This shortens the homing distance, and can be computed as shown later in Section 4.4.
5. S has reached steady-state at some t 0. This maps naturally to swarms with long-running autonomy, such as those deployed on large-scale agriculture and industrial automation applications.
None of the assumptions we made reduces the utility of our model in its application to real-world problems; i.e., they are consistent with the constraints imposed by the problems themselves and commonly used robot hardware. We will relax some of them during the evaluation of our model, showing that while they are helpful during derivation, they are not essential for practical applications of our model.

ODE model
Let M = {SS,DS,PL,RN} be the set of scenarios based on the block distributions shown in Fig. 2. In each m ∈ M, let the area where blocks can be distributed be a subset of the overall area A of the arena. Let j = 1, . . . , J be the sub-areas within A in which blocks can be distributed, each described by a tuple (A j ,c j ,d j ,ρ j ). A j is the area occupied by the sub-area j, c j is the center of the sub-area, d j is its dimensions, and ρ j is the mean steady-state block density within the sub-area. The area within A where blocks can be found is the union of these disjoint subsets: The value of ρ j varies across sub-areas in our model, so that the block encounter rate α b can be captured accurately even in extreme non-homogeneous block distributions. As an example, consider a degenerate case in which there are two areas, both represented by a circular sector around the nest. Robots would likely forage only from the inner circle and rarely reach the outer circle, and so a single ρ j is not sufficient to capture collective dynamics in such cases. We require both ρ j and the robot spatial occupancy distribution in order to accurately compute α b in a general way (see Section 4.5).
(1) to (3) in previous work described earlier in Section 3.
Eq. (8) uses the described block modeling method to capture the underlying behavioral dynamics of the swarm. We made the following additional assumptions about the block distribution for its derivation. First, whenever a block is redistributed all j are selected with probability proportional to the fraction of distributable area they contain. Second, blocks are distributed uniformly within a given j. Third, every j can hold any number of blocks, allowing two blocks to occupy the same location (i.e., stacking). We note that under our-steady state assumption, Eq. (8) can be solved analytically given B(0) and does not need to be solved numerically. We also observe that from Eq. (8) the number of blocks in the arena available for robots to find can increase and decrease, addressing one of the weaknesses in the model from previous work.
Next we derive analytical models for τ h ,α b , and α r from the arena geometry, number of blocks, block distribution, etc. We do not derive τ av , because it depends intrinsically on the interference avoidance strategy employed by S and therefore cannot be derived independently from κ without additional assumptions.

Derivation of homing time τ h
We build the intuition behind the block acquisition probability for a robot at location x in the arena as follows, modeling the nest as a single point x n . Since searching begins from the nest, the density of N s (t) must be greater near the nest because robots perform biased random walks which originate from a common point. Consequently, this non-uniform swarm spatial distribution means that the mean distance from the nest at which a searching robot encounters a block is not the same as the mean distance of a block from the nest. From (Codling et al, 2010), we have that the spatial occupancy distribution from the central point, as a result of a biased random walk with bias distribution Uni f orm(−θ ,θ ), falls off linearly, and we would intuitively expect the following: 1. It is moderated by ρ j = B j (t)/A j , because the rate of decay of the mean block acquisition distance as a function of distance from the nest within a given A j is slower for low ρ j . Consequently, a j sparsely populated with blocks will have minimal effect on the overall swarm block acquisition probability distribution, while a j densely packed with blocks will create an area of higher acquisition probability. 2. ρ j would play an exponentially moderating role only when the block acquisition location x is close to x n , such as for RN or PL block distributions. For SS and DS distributions, where the mean distance from a block to the nest is large, the effect of ρ j on block acquisition locations should be minimal.
Using the results of (Codling et al, 2010) and these intuitions, we formulate Eq. (9) as a close approximation of the occupancy distribution of a single random walker performing a correlated random walk starting from a central point x n . C is a normalization constant to ensure p acq j (x) integrates to 1 over all j.
Having defined the probability density function, we now derive the expected acquisition location by finding the expected values of the marginal density functions in x by integrating Eq. (9) and summing over all j: and similarly for y. We now write an expression for τ 1 h : v h is the robot phototaxis velocity, specified in the input configuration, and E[x acq ] − x n is the expected distance an acquired block will be from the center of the nest. d cr is the distance the homing path is shortened due to the employed congestion reduction strategy, and is straightforward to calculate from the arena geometry. We provide an example of computing d cr for a square arena, which in this work corresponds to the RN and PL scenarios; similar calculations can be done for the other arena shapes. Using a square nest of length L, without loss of generality we find the average distance to the origin of a randomly selected point in a region R := {(x,y) : 0 ≤ y ≤ x ≤ L/2} that is a triangle 1/8 of the nest. The uniform density on this region is then: Then, the average distance d cr is computed via: Finally, to derive τ h , we note that under our assumption of low to moderate ρ S , the homing time increases linearly with N according to the expected value of time lost due to inter-robot interference (Lerman and Galstyan, 2002), averaged across all robots: 4.5 Derivation of block acquisition rate α b Using our mean-field assumption, S can be approximated as a fluid composed of robot particles, and be considered to obey many of the same laws; in the long-term limit, the governing equation for the biased random walk used in this work is the advectiondiffusion equation (Codling et al, 2010). Given sufficiently simple robots, this approximation gives good results (Codling et al, 2010;Pang et al, 2019). Using simple robots does not detract from the utility of the system, which mimics the self-organizing behavior of real swarm populations such as social insects. Using this intuition, we obtain α b by computing the mean time it takes a robot "particle" starting at x n to "diffuse" within the operating area A to the expected acquisition location E[x acq ]. Viewing E[x acq ] − x n as the Root Mean Square (RMS) displacement distance and assuming a linear relationship between displacement and diffusion time, we obtain: where F(N) is the diffusion constant for a swarm of N robots, and α b is the expected time to diffuse from x n to E[x acq ]. We note that while the calculation for α b does not directly consider the distribution of blocks in the arena and their densities within a given block cluster, it depends on E[x acq ], which does.
An exact calculation of F(N) is out of the scope of this paper, so we approximate it as shown in Eq. (16) using the results of the RMS diffusion for CRW (Codling et al, 2010) and our intuition that smaller θ will result in quicker swarm diffusion (more straight line motion): where C d f (m) is a per-scenario parameter characterizing the linearity of the diffusion rate, and D xy is defined as follows from (Codling et al, 2010): where v s is the robot searching velocity and a scenario parameter. Under "normal" circumstances with inert particles, Eqs. (15) and (16) give good results because system diffusion varies linearly with time, which these equations require. However, many foraging environments are heterogeneous and have non-uniform distributions of blocks. As a result, the swarm may experience non-linear anomalous diffusion (Hasnain et al, 2018;Woringer et al, 2020), i.e., crowding in some areas of the environment, giving rise to the need for C d f (m), which implicitly captures these artifacts, as shown in Eq. (16).
To calculate C d f (m) from first principles, we take the following approach. First, we calculate the environment heterogeneity H(m) in reference to an environment with perfectly uniform block distribution (homogeneous case). For the square RN and PL scenarios, we compute the variance of the average inter-cluster distance. For the rectangular SS and DS scenarios with J ≤ 2, the average inter-cluster distance is not meaningful, so we calculate H(m) as the mean cluster distance from the "ideal", which is defined as an environment in which all c j are incident with x n . Modulating H(m) in Eqs. (18) and (19) by the expected level of crowdedness within each cluster we obtain an expression for C d f (m) which accounts for the number of block clusters (J) and the area cumulatively covered by all clusters (A d ): In Fig. 3, we plot the calculated values for C d f (m) for all scenarios to verify we have restored linearity of diffusion, and see strong correlations between C d f (m) and F(N) for all block distributions. As expected we see the highest heterogeneity associated with the lowest diffusion rate (PL), and the lowest heterogeneity associated with the highest diffusion rate (RN). With more observed heterogeneity within the environment, the scaling factor necessary to restore the linearity of the diffusion rate increases to accommodate for this occurrence. We also note that the differing x intercepts depend on the block distribution, and that as the size of the arena increases the crowding observed is not as pronounced (less anomalous diffusion). This might be due to the decreased likelihood of robots being "trapped" (Chepizhko and Peruani, 2013) by obstacles in the environment, leaving the robot more able to continue moving through the arena.

Derivation of robot encounter rate α r
To derive α r , we use our assumption of low to moderate ρ S , to model the robot encounter rate as a function of α 1 r and F(N). We view Fig. 1 as a queueing network, where robots are either performing collision avoidance maneuvers or not. The input/output transition rates for a state are summed to form the arrival and service rates for the collision avoidance queue Q ca (λ ca = α 1 r and µ ca = 1/τ av , respectively). Modeling Q ca as a M/M/1 queue, i.e., at most one robot exits collision avoidance per ∆t, which is reasonable if ∆t is small, we can write the following using Little's Law (Šeda et al, 2017) as: The second term in Eq. (20) is a corrective factor accounting for robots that experience interference due to encountering arena walls, not other robots, which is simply the scaled rate at which a single robot experiences interference. We cannot use N av (t) directly, as α r needs to be computed a priori, but we can estimate it as a function of α 1 r and N 1 av (t), using our intuition regarding swarm diffusion: where α 1 r and N 1 av (t) can be computed from κ using the results of (Codling et al, 2010). We increase the influence of θ in Eq. (21) by introducing another D θ in the denominator; smaller θ will result in more interrobot interference due to straight line motion. C ca (m) characterizes the sub-or super-linearity of Eq. (21) for a scenario m in relation to the random (RN) scenario.

Experimental setup
We use the ARGoS simulator (Pinciroli et al, 2012) 1 with a dynamical physics model of the marXbot robot in a 3D space for maximum fidelity (robots are still restricted to motion in the XY plane).
Robots run the FSM shown in Fig. 1, use short range proximity sensors for collision avoidance, and ground sensors to check if they are currently on top of a block. Blocks are abstracted as black squares on the arena floor. We ran four sets of experiments characterized by different swarm density, and evaluate our model across the scenarios shown in Fig. 2, relaxing some core assumptions made when deriving our model as we proceed.
We ran our experiments with N = 1 . . . 6,006 robots in environments with SS and DS source block distributions, and with N = 1 . . . 11,837 robots for RN and PL block distributions. Our scenarios are obstaclefree in the sense that they do not contain external obstacles; however, robots act as obstacles to each other. We use θ = π 36 for CRW and compute τ av from the details of our robot control algorithm κ (details omitted). C ca (m) is computed post-hoc for each experiment. For all experiments we performed 32 runs of T = 200,000 seconds, with a timestep of t = 0.2 seconds, i.e., the robot control algorithm runs at 5 Hz.

Results
For all the scenarios, we omit predictions for N s (t) since it can be computed from the conservation of robots via N s (t) = N − N av (t) − N h (t). We present our results with a non-scaled X-axis, with one data point per experiment, and a logarithmic Y-axis for larger N, to improve readability.
6.1 Constant ρ S , large N These are the "ideal" conditions under which to evaluate our model, because the linearity and mean-field assumptions we made when deriving our model hold.
In Fig. 4, we see strong agreement between the predictions of our model and experimental results for all scenarios, providing compelling evidence that our model is capturing the underlying dynamics of interrobot interference and searching accurately and that our mean-field and linearity assumptions are accurate. PL scenarios are the least favorable of all foraging environments, since they are asymmetrical and do not contain easily exploitable block clusters. Our model struggles to predict N h (t), N av (t) within the 95% confidence interval for PL scenarios, but does track the overall trend reasonably well, showing that our underlying diffusion model and assumptions about linearity of α r are generally accurate even in this difficult case. The divergence that does exist between predictions and experiments for PL scenarios suggests that Eq. (9) is moderately inaccurate; this is further supported by slight differences between experimental data and predictions for N h (t) in the RN and SS scenarios (a) Comparison of predictions of N h (t) (homing) behavior with experimental data for SS and DS scenarios.
(b) Comparison of predictions of N av (t) = N h av (t) + N s av (t) (avoidance) behavior with experimental data for SS and DS scenarios.
(c) Comparison of predictions of N h (t) (homing) behavior with experimental data for RN and PL scenarios.
(d) Comparison of predictions of N av (t) = N h av (t) + N s av (t) (avoidance) behavior with experimental data for RN and PL scenarios. environments. In Fig. 4b, when our predictions of N av (t) are inaccurate for some DS scenarios (in terms of confidence intervals), the inaccuracy may not be of practical concern as the predictions only differ from experimental results by < 1 robot/6,000.

Constant ρ S , small N
We relax the mean-field assumption, and evaluate our model's ability to capture finite-size effects when S cannot be reasonably approximated using the meanfield model (with N ≤ 50).
In Fig. 5, we see strong agreement between experimental results and the predictions of our model for all tested scenarios with small N. Even when this crucial assumption is relaxed, our model accurately captures finite-size effects, suggesting that the mean-field assumption may be valid at very small N, depending on the control algorithm used. The greatest discrepancy between predictions and results is shown in Fig. 5c for RN scenarios, which are inherently the most non-linear because of the placement of the nest relative to the objects to be gathered. In comparison (a) Comparison of predictions of N h (t) (homing) behavior with experimental data for SS and DS scenarios.
(b) Comparison of predictions of N av (t) = N h av (t) + N s av (t) (avoidance) behavior with experimental data for SS and DS scenarios.
(c) Comparison of predictions of N h (t) (homing) behavior with experimental data for RN and PL scenarios.
(d) Comparison of predictions of N av (t) = N h av (t) + N s av (t) (avoidance) with experimental data for RN and PL scenarios. with results from the previous section, stronger agreement between predictions and results with constant ρ S is seen when N is small, likely due to the inherent stochasticity of swarms which comes into play more strongly at larger N.

Variable ρ S , small N
We relax the linearity and mean-field assumptions, and evaluate our model's ability to capture finite-size effects when S exhibits non-linear behaviors with ρ S = 0.01 − 0.1 and N ≤ 50.
In Fig. 6 we see that our model generally performs well even when we violate its low ρ S assumption. We see strong agreement between the model and experimental data for all scenarios for ρ S = 0.01 − 0.05 for both N h (t) and N av (t). In Fig. 6d we see a striking divergence between our model and results at ρ S = 0.09, which may be a numerical anomaly or the point at which our ODE solver begins to struggle with the non-linear behavioral inputs.
(a) Comparison of predictions of N h (t) (homing) behavior with experimental data for SS and DS scenarios.
(b) Comparison of predictions of N av (t) = N h av (t) + N s av (t) (avoidance) behavior with experimental data for SS and DS scenarios.
(c) Comparison of predictions of N h (t) (homing) behavior with experimental data for RN and PL scenarios.
(d) N av (t) = N h av (t) + N s av (t) (avoidance) with experimental data for RN and PL scenarios.

Variable ρ S , large N
We restore the mean-field assumption to evaluate our model's ability to capture collective dynamics across scales on the order of 1,000s of robots when S exhibits non-linear behaviors. We believe that ρ S = 0.1 is an appropriate upper limit for real-world swarms of practical utility, as it corresponds to 1 robot/10m 2 which is an extremely high density in the context of feasible real-world applications.
In Fig. 7, we see that our model generally does not model swarm behavior at large scales well. In Fig. 7 we see divergence between model predictions and results for N av (t), N h (t) for all scenarios beyond ρ S = 0.02, and beyond ρ S = 0.01 for RN, PL scenarios, and our model fails to provide meaningful predictions, due to the highly non-linear swarm behaviors present, which is as we expected. For RN and PL scenarios, we see a strong correlation between the ρ S at which our model's predictions for both N h (t) and N av (t) break down, which happens at ≈ ρ S = 0.01.
(a) Comparison of predictions of N h (t) (homing) behavior with experimental data for SS and DS scenarios.
(b) Comparison of predictions of N av (t) = N h av (t) + N s av (t) (avoidance) behavior with experimental for SS and DS scenarios.
(c) Comparison of predictions of N h (t) (homing) behavior with experimental data for RN and PL scenarios.
(d) Comparison of predictions of N av (t) = N h av (t) + N s av (t) (avoidance) behavior with experimental data for RN and PL scenarios.

Discussion
Our results advance the state of the art in ODE modeling for foraging swarms (Lerman et al, 2001;Lerman and Galstyan, 2003) by demonstrating accurate predictions of swarm behaviors across scales and block distributions when the assumptions made in Section 4.2 are valid. To the best of our knowledge, our model is the first presented in the literature capable of capturing the collective dynamics of swarms when the mean-field approximation is valid, i.e., large N, but also when the mean-field approximation is not valid and finite-size effects are present. i.e., small N. Building on our results, accurate predictions of swarm performance and SR system properties, such as scalability and emergent self-organization using (Harwell et al, 2020), are possible.
As an example, we can predict swarm performance with reasonable accuracy using our model, as shown in Fig. 8. We define a performance measure P as the (a) SS and DS: N = 1 . . . 50.
(b) RN and PL: N = 1 . . . 50. Fig. 8 Predictions of swarm performance from first principles for different scenarios rate of block collection by the swarm S. We then predict it via P = α b /µ h , where µ h is the rate of robots leaving the homing queue Q h (analogous to Q ca ). Under our assumption of Q h as a M/M/1 queue, we set µ h = 1. In Fig. 8 we see strong agreement between our predictions and experimental results for SS and PL scenarios across all swarm sizes, and reasonable agreement for DS scenarios. Our predictions are less accurate for RN scenarios, likely due to our approximation of Q ca as an M/M/1 queue not being correct; i.e., our model does not account for multiple robots dropping an object in the nest on a single t. When we consider the swarm behavior with constant ρ S at large N (Fig. 4), we see our model struggle to predict behavior for all portions of the curve for all scenarios. We see that while our model predicts behavior at N ≥ 100 well, it struggles to predict behavior at N < 100 for SS, DS, and RN scenarios, underscoring the difficulty in creating a single linear model which can capture both mean-field and finite-size behaviors simultaneously. PL scenarios are the most challenging to model, and we see that while our model follows the general trend of behavior well, its predictions are usually outside of the confidence interval range. However, our model is accurate to within ≈ 16 robots/11,837, which may be sufficient in many cases.
Our linearity assumption, strictly speaking, is clearly not valid in all contexts across scales or across scenarios at similar N. Nevertheless, it is valid in many contexts (e.g., Figs. 5 and 6, parts of Fig. 4), and thus our model is of practical utility in designing SR systems to tackle problems across scales. Furthermore, our findings suggest that N always has a non-linear effect on behavior, even at extremely low densities, which is a function of the scenario itself. For more favorable scenarios, this non-linearity may be negligible even at large N, but for less favorable scenarios, such as RN and PL, it is not, and non-linear modeling must be applied for accurate prediction across scales.
When we relax our linearity assumption and increase ρ S to introduce non-linear behaviors in the swarm via inter-robot interference, our model still produces comprehensively accurate predictions of behavior at small N for all scenarios (Fig. 6). While our model will undoubtedly break down at some ρ S > 0.1, our results demonstrate that the threshold is above any realistic scenario, because densities > 0.1 are infeasible in general for real-world problems.
Capturing non-linear behaviors using linear models is much more successful if N is small; this aligns well with our intuition. If we imagine that the amount of time each robot spends avoiding collision with other robots at each timestep (t ca ) is roughly of the form t ca = N ρ S then it is clear that t ca will increase more quickly with N, even if ρ S is low.
Based on this discussion, we can confidently characterize the limits of linear modeling, at least for ODE based approaches, and say that it is possible to capture non-linear behaviors of swarms using linear models across scales. Furthermore, for N ≤ 50, doing so is also practical, as our model still provided predictions accurate within 95% confidence in all scenarios. However, scenarios exhibiting comparatively more nonlinear behaviors via N av (t), such as RN, have large 95% confidence intervals, reflecting the greater variability in behavior as a result of non-linearity. As a result, our model's predictions of the swarm average behavior may be less useful in such scenarios, since the model provides predictions of the average case. We note that while the results of this work suggest that linear modeling of non-linear swarm behaviors is possible using other modeling approaches, they do not necessarily generalize to other modeling paradigms.
Our model still relies on one parameter: C ca (m). While this is a free parameter in the strictest sense, it has a strong semantic context rather than just being a fitting parameter to ODE terms to improve prediction accuracy. This is a substantial improvement over previous work which used five context-free fitting parameters (Lerman and Galstyan, 2002). Thus, while our model is not (yet) strictly fundamental, the semantic context of C ca (m) strongly suggests that it can be determined analytically, and first-principle modeling of general foraging behavior is possible.
Finally, it is important to emphasize the presented results were obtained in simulation, commensurate with other mathematical approaches to modeling SR systems (Guerrero-Bonilla et al, 2020;Tarapore et al, 2017). Even if our model does not produce equivalently accurate results to those shown in Figs. 5 and 6 with systems of real robots, we believe that it will still provide a strong starting point for modeling systems of real robots, thereby accelerating the deployment of swarm solutions to real-world problems, for two reasons. First, the robot FSM used in this work (Fig. 1) is simple: it has minimal sensor usage, and does not use memory. This simplicity will make real-world systems of such robots especially amenable to mathematical modeling. Second, our model has been shown to be accurate with swarms of non-trivial size (e.g., up to 46 robots in Fig. 5, up to 51 robots in Fig. 6, parts of Fig. 4 with up to 11,837 robots). In such swarms, nondeterministic transient behaviors will reliably arise even in simulation, e.g., floating point representation errors; similar errors will arise even in small swarms of real robots.

Conclusions and future work
We have investigated how to characterize the limits of linear modeling of swarm behaviors. We have developed a robust ODE model for foraging swarms, and shown it is accurate in a wide variety of scenarios in swarms which exhibit linear behaviors, and are large enough to be approximated using mean-fields, demonstrating its utility as an effective modeling tool. Furthermore, in swarms for which either assumption is not true, our model in many cases is still accurate. Future work will derive the block acquisition density function using random walks theory and refine our model with experiments using real robots.