Arbitrage Equilibrium, Invariance, and the Emergence of Spontaneous Order in the Dynamics of Bird-like Agents

The physics of active biological matter, such as bacterial colonies and bird flocks, exhibiting interesting self-organizing dynamical behavior has gained considerable importance in recent years. Current theoretical advances use techniques from hydrodynamics, kinetic theory, and non-equilibrium statistical physics. However, for biological agents, these approaches do not seem to recognize explicitly their critical feature: namely, the role of survival-driven purpose and the attendant pursuit of maximum utility. Here, we propose a game-theoretic framework, statistical teleodynamics, that demonstrates that the bird-like agents self-organize dynamically into flocks to approach a stable arbitrage equilibriumof equal effective utilities. This is essentially the invisible handmechanism of Adam Smith’s in an ecological context. What we demonstrate is for ideal systems, similar to the ideal gas or Ising model in thermodynamics. The next steps would involve examining and learning how real swarms behave compared to their ideal versions. Our theory is not limited to just birds flocking but can be adapted for the self-organizing dynamics of other active matter systems.

Flocking has been studied extensively from dynamical systems and statistical mechanics perspectives [10][11][12][13][14].Such analyses have contributed substantially to our evolving understanding of interesting emergent properties such as phase segregation, flock stability, etc.However, these approaches don't seem to model explicitly the critical feature of active biological agents, namely, the role of purpose and its attendant pursuit of maximum utility.Being biological agents, birds are innately purposeful, driven by the goal to survive and thrive in challenging environments as Darwin explained.We believe any comprehensive theory of active biological matter has to overtly account for this defining characteristic of the agents.
We address this need by using a novel game-theoretic framework, which we call statistical teleodynamics [15][16][17][18][19][20].The name comes from the Greek word telos, which means goal.Just as the dynamical behavior of gas molecules is driven by thermal agitation (hence, thermodynamics), the dynamics of purposeful agents * venkat@columbia.edu is driven by the pursuit of their goals and, hence, teleodynamics.Statistical teleodynamics may be considered as the natural generalization of statistical thermodynamics for purpose-driven agents in active matter.It is a synthesis of the central concepts and techniques of potential games theory with those of statistical mechanics towards a unified theory of emergent equilibrium phenomena in active and passive matter [20].
In this paper, we study a model of birds flocking that is inspired by the well-known Reynolds's Boids model, analytically and computationally.In the next two sections, we briefly review the dynamical systems and statistical mechanics perspectives, respectively.This is followed by our statistical teleodynamics formulation, which starts with a quick introduction to potential games.We conclude with a discussion of the main results and their implications.

II. DYNAMICAL MODELS OF FLOCKING
Flocking has been studied using dynamical models, which describe the time evolution of the position r i of the i th agent and its velocity v i using pre-specified rules for change.The state of the flock at any given time is specified by specifying r i and v i for all agents.When the agents move at a constant speed v 0 , the state of the system is then determined by the set of agents' positions and velocity directions or orientations {r i , s i } N i=1 , where N is the number of agents.Both the Reynolds's model [10] and the Vicsek model [11] describe the time evolution of an agent's velocity, but using different force models.
An agent i is said to be affected an agent j, if j is in the neighborhood of i, N i .The neighborhood N i of i is defined by a matrix whose elements are n ij , where The span of the neighborhood is specified in terms of the absolute distance between i and j, and a size parameter r 0 [11,21], such that As seen from Eq. 2, we consider agent i to be its own neighbor.One can also define a neighborhood in terms of a fixed topology of nearest neighbors [12][13][14], but we don't use this specification in this study.It follows that the number of neighbors of an agent i is given by n i = j n ij .
In the Reynolds's model, the agents (called boids) obey the following three rules as they fly around: In general, we can write the net effect of these forces on the i th boid by the equation (Eq.3), where a, b, and c are parameters corresponding to the rules of cohesion, separation, and alignment, respectively.Parameter η is the uncorrelated noise in the agent's velocity.The time-scale, ∆t, in Eq. 3 can be subsumed in a, b, c.The Vicsek model, similarly, updates velocity purely as a function of the alignment of the agent with its neighbors, though modifications have been proposed to include pairwise forces [13,21] (see also, Supplementary Information, Section 2).
It's important to note that in both the Reynolds model and the Viscek model, there is no concept of the "final" state, or an equilibrium state, of the system as time tends to infinity.In this regard, they are like molecular dynamics simulations of molecules where there is no concept of a final equilibrium state in the equations.They can only determine the immediate next move of the molecules, at any given time, not their final configurations.The final outcome is determined, a posteriori, after the simulations are run for a long time.We highlight this important point here as our theory differs conceptually in this regard.
While dynamical models of flocking have been studied extensively in literature, work has also been reported using statistical mechanics to analyze the flocking behavior.These methods typically use a maximum entropy formulation via an Ising-model inspired Hamiltonian of the boids' interaction [13,14].

III. STATISTICAL TELEODYNAMICS OF FLOCKING -A GAME-THEORETIC FORMULATION
The Reynolds and Vicsek models specify bottom-up agent-level dynamical behavior, but they don't provide an analytical framework to predict the dynamics of the entire flock.This is determined only computationally via agent-based simulations.That is, there is no analytical framework to derive the behavior of the whole from the behaviors of the parts.
On the other hand, the statistical mechanics formulation [13] is a top-down approach that starts with the specification of the Hamiltonian of the flock and then imposes the maximum entropy distribution on it.It is not clear why maximum entropy, which is obviously relevant for passive matter systems, would be applicable for survival-driven birds.The typical statistical mechanics approach uses the superficial similarity between spins in magnetic systems (e.g., the Ising model) and the orientation of the birds to apply maximum entropy methods.The deeper question of why this is conceptually relevant for the birds is not addressed.Most importantly, all these approaches don't seem to recognize explicitly that active agents such as birds act instinctively to improve their survival prospects.
We address these challenges using our statistical teleodynamics framework [15][16][17][18][19][20].In this theory, the fundamental quantity is an agent's effective utility which is a measure of the net benefits pursued by the agent.Every agent behaves strategically to increase its effective utility by switching states and exploiting arbitrage opportuni-ties.
In our theory of flocking, we propose that birds are arbitrageurs that always maneuver to increase their effective utilities, which determine their survival prospects, dynamically in flight.The effective utility of a bird depends on its position, speed, and alignment with the rest of the members in its neighborhood.
Thus, we interpret the three rules of engagement for Reynolds's boids not as externally imposed forces on the agents but as innate, self-actualizing, properties of the agents acquired over millions of years of Darwinian evolution.These instinctive characteristics enable the agents to incessantly search for better effective utilities in order to improve their survival chances.
Hence, we believe that the proper formulation of flocking ought to start with a model of effective utility that a bird uses to make such decisions dynamically in flight.Seen from this perspective, we suggest that birds do not fly randomly (as statistical mechanics-based formulations implicitly assume), but maneuver strategically to improve their utilities.We exploit this critical insight to model the dynamical behavior of birds in flight by using the concepts and techniques from potential games.
In potential games, there exists a single scalar-valued global function, called a potential (φ(x)) that has the necessary information about the payoffs or the utilities of the agents.The gradient of the potential is the utility, h i , of the ith agent [15,[22][23][24][25]. Therefore, we have where x i = N i /N and x is the population vector.A potential game reaches equilibrium, called Nash equilibrium (NE), when the potential φ(x) is maximized.Furthermore, this Nash equilibrium is unique if φ(x) is strictly concave [23].At Nash equilibrium, all agents enjoy the same effective utility, i.e., h i = h * .In fact, the equality of effective utilities in all states is the fundamental criterion of game-theoretic equilibrium for active matter.It is an arbitrage equilibrium [19] where the agents don't have any incentive to switch states anymore as all states provide the same effective utility h * .Thus, the maximization of φ and h i = h * are exactly equivalent criteria, and both specify the same outcome, namely, an arbitrage equilibrium.
There is a deep and beautiful connection between potential game theory and statistical mechanics as discussed by Venkatasubramanian [17,20].Since an elaborate discussion about this would take us afar from the objectives of this paper, we refer interested readers to [17,20] .

A. Garud's Utility: Position Dependence
Our goal here is to develop a simple model of the effective utility (h i ) of our boid-like agent, called garud (after the legendary king of birds, Garuda, in Indian mythology).
We want the model to be an appropriate coarse-grained description of the system that can make useful predictions not restricted by system-specific nuances.We have tried, deliberately, to keep the model as simple as possible without losing key insights and relevance to empirical phenomena.One can add more complexity as and when desired later on.What we are aiming for is the equivalent of the ideal gas model or the Ising model for birds flocking.
We develop our teleodynamical model using Reynolds's model as the start, but our approach is not restricted to this example alone; it is applicable to other models as well.
We consider a garud's position in the frame of reference of the center of mass of its neighborhood.We then apply the rule of cohesion and rule of separation to formulate the model for utility.
The rule of cohesion requires the garuds to come together and hence an i th garud's utility increases as it has more neighbors, n i .However, the increased utility comes at the cost of congestion, the disutility of congestion (corresponding to the rule of separation).The trade-off between the two terms, the benefit-cost trade-off, results in an inverted-U profile, which, following Venkatasubramanian [17], can be parameterized as, where r is the position component of the utility for the i th α, β > 0. Note that the positional dependence is accounted for in the computation of n i .Given a configuration of {r i }, the neighborhood of the i th garud is defined by the parameter r 0 , where if j th garud is within this radius then it is considered a neighbor.
This in turn identifies a direction of increased utility, given by, = α ∂n i ∂r i − 2βn i ∂n i ∂r i (6) is dependent on the garuds in the perimeter of the neighborhood of the reference garud i.

B. Garud's Utility: Velocity Dependence
The utility of a garud is also dependent on the velocity of its neighbors in that the garud attempts to match the orientation with its neighboring garuds.This utility component (h v ) can be written as, The utility of the i th garud, then, depends on the orientation of the other garuds in its neighborhood, i.e., s i • s j where j is a neighbor of garud i.This gives the alignment utility for the i th garud as, where n ij shows if the j th garud is a neighbor of garud i, i.e., If each garud is perfectly aligned with its neighbors, this utility component is maximal, whereas if they are oriented in the opposite direction this is minimal.Therefore, the garuds prefer to be aligned.This gives the i th garud an arbitrage opportunity to adjust its velocity vector towards this direction to increase its utility.This opportunity for increasing its utility generates a self-propelled force on the i th garud.If the i th garud is not aligned with its neighbors, this direction of increased utility is given by,

C. Garud's Effective Utility
There is one other utility component remaining to be considered.This is not stated explicitly in the three rules of the boids.However, it is implied because it is assumed that the boids have to be moving constantly.
So, as a garud incessantly moves and jockeys for better positions and orientations, its ability to do so is hampered by the competition from other garuds in the neighborhood that are also trying to do the same.As Venkatasubramanian explains [17], this disutility of competition can be modeled as −δ ln n i .
This term, when integrated to obtain the potential φ(x), leads to entropy in statistical mechanics.Thus, maximizing the potential φ(x) is equivalent to maximizing entropy under certain conditions.For more details, the reader is referred to Venkatasubramanian [17,20].Now, by combining all these components, we arrive at the effective utility for the i th garud given by where l i = 1 ni j n ij s i • s j is the average alignment of agent i.Without any loss of generality, δ can be assumed to be 1, and will be done so for the rest of this paper.
When α, β, γ = 0, the garuds don't have any preferences and hence fly randomly.This is what is captured by the remaining − ln n i term, which we call entropic restlessness.
Statistical teleodynamics, via potential game theory, proves that the self-organizing dynamics of the garuds would eventually result in an arbitrage equilibrium where the effective utilities of all the garuds are the same, i.e., h i = h * .
In the next section, we discuss our simulation results that confirm this prediction.

IV. RESULTS
Effective utility and its derivative as a function of the number of neighbors ni, for different values of alignment li (α, β, γ, δ = 0.5, 0.005, 0.25, 1).There are two locations where the derivative of the effective utility is zero for different alignments.
For the simulation details, the reader is referred to the methodology section VII below.If the garuds are flying randomly, without any rules of behavior, then this base case corresponds to α, β, γ, = 0; δ = 1.This result is discussed in Section S4.
For the other cases, the effective utility function in Eq. 12 is plotted in Fig. 1 in terms of the number of neighbors of garud i (n i ), for different alignments (l i ) for a given set of α, β, γ.We see that there are two values of n i (n − and n + ) where the gradient of utility, for a given value of alignment, is zero.These values are determined by (α (see, Supplementary Information, Section 1).In Fig. 1, at the lower value (n − ), any deviation in n i increases the utility of the garud, and hence leads to an unstable point.However, for the higher n i point (n + ), we see that any deviation reduces the garud's utility.Therefore, this leads to a stable point as any deviation would bring a garud back to the higher utility state.Therefore, this is the point a garud will try to reach to maximize its utility.For example, for the red curve, this would correspond to the point where n i = 73.6.
However, despite this point's stability, a garud will not be able to stay there indefinitely as the other garuds in its neighborhood are constantly changing their positions and orientations in their flights.Therefore, the i th garud would be fluctuating around this point.
In Fig. 2 and 3, we show the simulation results of both the Reynolds's Boids model and our utility-driven Garud model dynamics.We show the results of the Reynolds's model for illustrative purposes only as our objective is not to mimic the Reynolds's model exactly.We just want to show that the utility-driven model's collective behavior is very similar to that of the Reynolds model.
The results are shown for different parameter values of (a, b, c) and (α, β, γ).From the simulations, we obtain a set of position and velocity values {r i , v i } of each agent i at every time step.Once this is obtained, we extract the features n i and l i for all the agents.This in turn, is used to compute the average number of neighbors n and average alignment l for the entire population for all time points (Fig. 2).
Figure 2 shows the snapshots at different time points of the evolution in 3D-space and the n i -l i phase-space.While the exact dynamics, the exact configuration of the population, and the time-scale of evolution of these two models cannot be the same (and that is not the aim, either), we observe that the qualitative patterns of collective behavior are very similar in both cases.In particular, we see that all agents gravitate towards a certain region in the n i − l i phase space for both models.They both start at the lower far right point (where the average alignment is near zero as the agents are all randomly oriented initially, and closely packed) and evolve towards the upper center-right region in black.We notice a qualitative match of both the trajectories.
Furthermore, we can also see a quantitative similarity between the two models for specific parameters (Fig. 3).Fig. 3a shows the phase-space for the Reynolds's boids model, and Fig. 3b shows the same for the utility-driven model.
In both models, we notice similar features of evolution towards the arbitrage equilibrium states, starting from the lower right point at time t = 0 to ending in the colored regions, where the average number of neighbors and the average alignment fall in similar corresponding regions.The plots in the right show the average alignment and average number of neighbors of the agents in the last 100 time steps.
While both models exhibit similar collective behaviors, it is not apparent, however, from the three rules of the Reynolds's model, that its dynamics would result in an equilibrium state in the n i − l i phase space.This is different for the utility-driven model, however.Since its potential game formulation predicts an arbitrage equilibrium outcome, it is clear right from the beginning where in the phase space the system is going to end up in.We can make a quantitative prediction about the average n * i and h * values at equilibrium.This ability to predict the final outcome of the collective behavior of the population, given the individual agent-level properties captured in the utility function h i , is an important defining feature and strength of the statistical teleodynamics framework.An additional characteristic is the ability to prove the stability of the final outcome.
We also ran the simulations for different time-step sizes of ∆t = 0.01, 0.1, 0.5 in Eq. 3 to understand the dynamics of the evolution better.Note in Figure 4, at the start (t = 0), the utilities of all the garuds are spread out, with many having negative utility values, and the average utility ( h) low.
But as the dynamics evolves, every garud tries to increase its utility by maneuvering to a better neighborhood and better orientation, the distribution becomes narrower, the average utility keeps increasing, and reaches a near-maximum value ( h = 22.17 ± 2.90 in Fig. 4a) and fluctuates around it.Note that this is around the maximum theoretical value of about 23.8 (given by Eq. 12), where the histogram peaks.This suggests that nearly all the garuds have similar effective utilities asymptotically, approaching the maximum.This, of course, is the arbitrage equilibrium outcome predicted by the theory (see also, Supplementary information, Fig. S3).
The garuds do not converge exactly on h * but fluctuate around it because of the stochastic dynamics.This is also seen in Table I where nearly the top 10 % of the garuds at a particular time-step are very close to the maximum utility value.In fact, the top 50% of the garuds have an average utility of greater than 23.This arbitrage equilibrium state is unique only if the potential φ(x) is strictly concave [23].For garuds, this is not the case in general as the concavity would depend on α, β, and γ having some particular values.So, for the typical case where φ(x) is not concave, there could be multiple equilibrium configurations of the garuds.Thus, instead of an equilibrium point in the n i − l i phase space, one has an equilibrium region, in general.In other words, invoking a terminology from chaos and nonlinear dynamics, there would be a basin of attraction in the phase space where the garuds finally settle in and fly around.This is what we see in Fig. 2 and 3 in the colored regions.

A. Stability of the Arbitrage Equilibrium
We can ascertain the stability of this equilibrium by performing a Lyapunov stability analysis [17].A Lyapunov function V is a continuously differentiable function that takes positive values everywhere except at the equilibrium point (i.e., V is positive definite), and decreases (or is nonincreasing) along every trajectory traversed by the dynamical system ( V is negative definite or negative semidefinite).A dynamical system is locally stable at equilibrium if V is negative semidefinite and is asymptotically stable if V is negative definite.
Following Venkatasubramanian [17], we identify our Lyapunov function V (n i ) where φ * is the potential at the Nash equilibrium (recall that φ * is at its maximum at NE) and φ(n i ) is the potential at any other state.Note that V (n i ) has the desirable properties we seek: (i) V (n * ) = 0 at NE and V (n i ) > 0 elsewhere, i.e., V (n i ) is positive definite; (ii) since φ(n i ) increases as it approaches the maximum, V (n i ) decreases with time, and hence it is easy to see that V is negative definite.Therefore, the arbitrage equilibrium is not only stable but also asymptotically stable.
Our simulation results confirm this theoretical prediction (see Figure 5).After the garuds population reached equilibrium, we disturbed the equilibrium state by randomly changing the positions and/or velocities of the garuds.The simulation is then continued from the new disturbed far-from-equilibrium state.We conducted experiments with three kinds of disturbances: • Disturbance 1:Velocity disturbance, where each garud's velocity is changed to a random orientation and magnitude.
• Disturbance 3: Position and velocity disturbance, where both position and velocity vectors are changed.
As seen in Figure 5, after the 100 th time step when the population had reached equilibrium, we introduced these disturbances.The 101 st time step shows in red color the new far-from equilibrium states, where the average utility has dropped considerably.
In all cases, the population recovers quickly, typically in another 100 time steps or so to reach the original equilibrium region (shown in green).The figure shows the disturbance (red) and recovery (green) in both the 3-D space and the phase space.
In Fig. 5a, as the velocities are randomized at the 101 st time step, the alignment goes down to 0, but recovers to the original equilibrium quickly.In Fig. 5b, as the new configuration corresponds to a similar value of average number of neighbors as before, the disturbance is not that much.Note that the drop in average utility is small.In Fig. 5c, we see that this disturbance is huge, pushing the configuration close to the original random state.But still, the population is able to recover to the arbitrage equilibrium quickly.
This shows that the arbitrage equilibrium region is not only stable, but asymptotically stable.That is, the garuds flocking configuration is resilient and self-healing.Given the speed of the recovery, it could possibly be exponentially stable but we have not proved this analytically here.
The asymptotic stability of this arbitrage equilibrium is similar to that of the income-game dynamics as discussed by Venkatasubramanian [15,17] using a similar Lyapunov stability analysis.

V. CONCLUSION
For three centuries we have known that there are constants of motion, such as energy and momentum, for passive matter.Nevertheless, it comes as a surprise to discover that the dynamics of active matter populations could also have an invariant, namely, the effective utility.However, the role of invariance here is different from its role in dynamics.The constants of motion such as energy and momentum are conserved, but effective utility is not.
The role of this invariance is more like that of set point tracking and disturbance rejection in feedback control systems.These are called the regulation and servo prob-lems, respectively, in control theory [26,27].The system, i.e., the garuds population, adjusts itself dynamically and continually, in a feedback control-like manner, to maintain its overall effective utility.
It is important to emphasize, however, that this control action is decentralized as opposed to the typical centralized control system in many engineering applications.The agents individually self-organize, adapt, and dynamically course-correct to offset the negative impact on their effective utilities by other agents or other external sources of disturbance.The population as a whole stochastically evolves towards the stable basin of attraction in the phase space in a self-organized and distributed-control fashion.This is essentially Adam Smith's invisible hand mechanism of economics.As Smith observed [28], "It is not from the benevolence of the butcher, the brewer, or the baker, that we expect our dinner, but from their regard to their own interest.We address ourselves, not to their humanity but to their self-love, and never talk to them of our own necessities but of their advantages."Thus, every garud is pursuing its own self-interest, to increase its own h i , and a stable collective order spontaneously emerges via such self-organization.
Invariants are quite rare in physics, rarer still in biology and economics.That's why it is exciting to see them as their presence usually signals something deep, something fundamental.In physics, their existence has revealed deeply fundamental symmetries of the cosmos, as Emmy Noether showed.Therefore, it is important to understand the implications of our discovery in sufficient depth and breadth.
We do realize, of course, that our bird-like agents are not real birds.Our model and simulations are not real biological systems.Nevertheless, our results suggest intriguing possibilities for real biological entities that need to be explored carefully.
Therefore, the interesting and surprising result, seen both analytically and in simulations, that the emergent arbitrage equilibrium is asymptotically stable is an important one with potentially far-reaching consequences, particularly in biological, ecological, and economic contexts.For example, this could be an important mechanism of pattern formation and pattern stability in biological systems.Populations of cells could self-organize, under different spatial and temporal conditions and constraints, driven by their incessant and instinctive hunt for better utilities, to settle into various stable basins of attraction -i.e., into different types of stable emergent order -to form stable organized structures.Their asymptotic stability property bestows upon them the resilient self-healing feature found so commonly in many biological systems.This process could be a core mechanism behind the design, control, and optimization of stable biological systems via self-organization.
This mechanism is applicable to different length and time scales, from molecular to macroscopic to planetary scales.These results raise several interesting questions about populations of biological active matter competing with one another.For instance, consider all the different kinds of microbial populations in the human body, or for that matter, in any living organism.Not only are all the individual microbes competing with one other strategically for resources to increase their effective utilities, to improve their survival and growth fitness, different microbial populations are also competing with one another at a higher-scale.
So, is there a hierarchy of arbitrage equilibria?That is, the microbes in a given population at some arbitrage equilibrium among themselves (say, level-1 equilibrium), and such populations themselves are in equilibrium with one another at a higher level (say, level-2), and so on.Is there a planet-scale equilibrium?That is, are all the living species on our planet, along with the environment, at some arbitrage equilibrium?Or are we evolving towards one, just like the garuds population did in this study.What happens when this equilibrium is upset by either internal disturbances (such as climate change) or external shocks (such as asteroid impact)?
As one can see, our theory is not limited to just birds flocking.It is also applicable to the self-organizing dynamics and evolution of a wide variety of systems in physics, biology, sociology, and economics.As Venkatasubramanian et al. [20] showed the emergence of the exponential energy (i.e., Boltzmann) distribution for gas molecules can be modeled by the effective utility Similarly, they showed, [20] as examples of biological systems, bacterial chemotaxis can be modeled by and the emergence of ant craters by The same study showed how the Schelling game-like segregation dynamics in sociology can be modeled by and the income game in economics by What we have is for ideal systems, similar to the ideal gas or Ising model in thermodynamics.Just as real gases and liquids don't behave exactly like their ideal versions in statistical thermodynamics, we don't expect real biological systems (or economic or ecological systems) to behave like their ideal counterparts in statistical teleodynamics.Nevertheless, the ideal versions serve as useful starting and reference points as we develop more comprehensive models of active matter systems.The next steps would involve examining and learning how real-world biological systems behave compared to their ideal versions.This would, of course, necessitate several modifications to the ideal models.
We note, from equations 14-18, a certain pattern in the structure of the effective utility functions in different domains.Thus, we see that the same mathematical and conceptual framework is able to predict and explain the emergence of spontaneous order via self-organization to reach arbitrage equilibrium in dynamical systems in physics, biology, sociology, and economics.
This kind of universality is particularly striking, prompting us to conclude with a quote a from the inimitable Richard Feynman that seems apropos here: "Nature uses only the longest threads to weave her patterns, so that each small piece of her fabric reveals the organization of the entire tapestry."It appears that the emergence of spontaneous order via self-organizing stable arbitrage equilibria is such a thread.

VI. ACKNOWLEDGEMENTS
This work was supported in part by the Center for the Management of Systemic Risk (CMSR), Columbia University.This manuscript was written when the corresponding author (VV) was the Otto Monsted Distinguished Visiting Professor at the Danish Technical University (DTU) as well as a resident of Nyhavn 18, Copenhagen, as a guest of the Danish Central Bank in the summer of 2022.It is with great pleasure that VV acknowledges the generous support of these institutions and the warm hospitality of his colleagues in the Chemical Engineering Department at DTU.

VIII. METHODOLOGY
We created a simulation of 1000 garuds in a periodicbox of dimensions 20×20×20, where each garud's neighborhood is a sphere with radius r 0 = 3.Each garud starts at a random location and orientation inside a 10×10×10 block.Speed of each garud is limited between 0.5 and 1.The update algorithm works similar to the Reynolds's garuds update, except the force is driven by the numerical estimates of direction of increased utility (additively based on position and velocity).An additional noise is also added to the velocity update strategy similar to the Reynolds's model to capture the erroneous strategies of velocity update for each garud.This is given by a noise parameter (0.01, unless specified) times the magnitude of the velocity.The noise indicates that a garud does not make perfect choices in updating its velocity.

Venkat Venkatasubramanian *
Complex Resilient Intelligent Systems Laboratory, Department of Chemical Engineering, Columbia University, New York, NY 10027, U.S.A.

S1. UTILITY MODEL FOR BIRDS FLOCKING
We can see that the number of neighbors for a garud n i = j n ij .The garuds try to maximize the effective utility, given by . We define the average alignment of each garud as This results in the alternative formulation of the utility in terms of the alignment, The first three terms are the utility of cohesion, disutility of congestion, and the utility of alignment.The last term is the disutility due to entropic restlessness.

S2. OPTIMUM BASED ON THE UTILITY FORMULATION
At equilibrium, all garuds have the same utility, i.e., Now, h i is maximum at two n i values where the gradient is zero, given by, Note that n − is an unstable point as any deviation in the number of neighbors would result in increasing utility, thereby causing the garud to move away from there.On the other hand, n + is a stable point because any deviation would decrease the utility, thereby causing the garud to return to its original state.

S3. DYNAMICAL MODELS OF FLOCKING
Regarding the discussion in Section II, the net effect of the three forces in the Reynolds model on the velocity of the ith boid is modeled by the equation, where a, b, and c are parameters corresponding to the rule of cohesion, rule of separation, and rule of alignment, respectively, v c,i is the average velocity of the neighbors of i, and r c,i is the center of the neighborhood as perceived by the agent i [1].These are given by Parameter η is the uncorrelated noise in the agent's velocity.Substituting the average velocity of the neighbors and the center of the neighborhood as perceived by the agent i, Eq.S1 can be simplified to give, In general, we can write the above equation as, * venkat@columbia.edufor a time-step ∆t.The time-scale ∆t in Eq.S3 can be subsumed by the parameters to give Eq.S2 The Vicsek model is a similar model where the velocity update is purely a function of the alignment of an agent with it's neighbors.The constant velocity dynamics is sometimes modified to include other pair-wise attractionrepulsion forces f ij and is written as [2], where Θ scales the dynamics to a unit vector to ensure the constant speed of the individual agents.
S4. BASE CASE: α, β, γ = 0 Fig. S1 shows the case where α, β, γ are set to zero, so the garuds are entropically driven (only − ln n i component is driving the motion of the garuds).In this case, the garuds start with a high density and zero alignment (Fig. S1a) initial configuration and finally settle in a configuration where they have random positions and velocities (Fig. S1b).If each garud is randomly located in 3D space, the number of neighbors for each garud on average is given by, ni = N L 3 where L = 20 is the length of the domain and r 0 = 3 is the size of the neighborhood.We subtract 1 in Eq.S4, as the i th garud itself is not considered in the average number of neighbors.With a total of 1000 boids, this gives an estimate of the average number of neighbors of the i th boid as 13.1.The random alignment of the garud also gives the average alignment with neighbors as zero, as there is no incentive to increasing the alignment.The simulation results confirm the theoretical expectations.

S5. NOISE TREND
The noise parameter is also varied in the simulations.To reiterate, the noise-parameter dictates the magnitude of randomness in the velocity vector which is added to the ideal "direction" of increasing utility.As the noise keeps increasing, we see that the system tends more and more towards randomness characterized by random positions and velocities of agents (dashed lines in Fig. S2).

S6. DYNAMICS VARIATION FOR DIFFERENT TIME STEP SIZES
We also ran the simulation for different time step sizes of ∆t = 0.01, 0.1, 0.5 in Eq.S3.Fig. S3 shows that

1 . 2 . 3 .
Rule of cohesion: A boid steers to move towards the average position of local flockmates Rule of separation: A boid steers to avoid collision and crowding of local flockmates Rule of alignment: A boid steers towards the average heading of local flockmates i − βn 2 i + γn i l i − δlnn i ) li = -1.0li = 0.0 li =
FIG. S2.Neighbors (a) and alignment (b) dependence on noise parameter for utility parameters α, β, γ = 0.5, 0.005, 0.25.We see that increasing the noise in decisionmaking results in the system tending more towards random behavior (dashed lines in a and b)

TABLE I .
Utility of different percentiles of the garuds at the 1000 th time step corresponding to Figure4Population Time step size, ∆t Average utility