D ec 2 01 9 Similarities between Insect Swarms and Isothermal Globular Clusters

Dan Gorbonos, Kasper van der Vaart, Michael Sinhuber, James G. Puckett, Nicholas T. Ouellette, Andrew M. Reynolds, and Nir S. Gov ∗ Department of Chemical Physics, The Weizmann Institute of Science, P.O. Box 26, Rehovot, Israel 76100 Department of Civil and Environmental Engineering, Stanford University, Stanford, California 94305, USA Department of Physics, Gettysburg College, Gettsyburg, Pennsylvania 17325, USA Rothamsted Research, Harpenden, AL5 2JQ, UK

Insect swarms are a canonical example of collective animal behavior [1,2], displaying group-level cohesion and stability even in the presence of environmental noise [3][4][5][6][7][8]. But while in many other forms of collective animal motion, such as flocking, schooling, and herding [2], the movement of individuals is coordinated, swarms are distinguished by their lack of globally aligned motion and the swarm state is not described by any order parameter.
It is thought that many swarming insect species such as Chironomus riparius, the midge species we consider here, interact predominantly via long-range acoustic sensing [9]. Indeed, a theoretical model that assumes that midges accelerate toward the sounds produced by other individuals has produced nontrivial results that are in good agreement with empirical observations [10,11]. The acoustic field produced by flying insects has a monopole component whose intensity falls off according to an inverse-square law [12], a scaling that is similar to the way the gravitational pull between objects falls off with distance. Dipole and other higher-order multipole components decay more rapidly, and are thus weaker than the monopole component. It is therefore tempting to speculate, as Okubo [3] and then Gorbonos et al. [10] did, that midge swarms are analogous to Nbody self-gravitating systems, where cohesion originates from the gravitational pull between the bodies compromising the system. A key contribution from [10] was to incorporate the adaptive gain of typical biological sensors, which leads to results that are different from Newtonian gravitation but are in good agreement with empirical measurements of swarms. The analogy between swarms and self-gravitating systems is appealing because it is well known that gravity can produce complex dynamical behavior from simple interactions-just as is thought to be the case for collective animal behavior [13]. Thus, by making a link between these two distinct systems, we can draw on the intuition built up from studying gravity to gain insight into collective animal behavior.
In this Letter, we examine this analogy more closely by studying the spatial variation of the number density and the velocity of individuals. We compare to the classic King model for the mass distribution in isothermal globular clusters [14] (see SI). Although we find some similarities with the observed swarms, there are also significant differences. These discrepancies reveal the limitations of a pure gravitational model for swarms. However, we find that the addition of adaptivity, which we can introduce only using numerical simulations, together with short-range repulsion can mitigate these limitations, producing mass distributions and dynamics that capture all the main features of the data. Thus, our results provide further support for adaptive gravity as a model for the swarm behavior.
We begin by considering the mass distributions predicted by the original King model. This model has been well studied previously. Chavanis et al. [15], for example, described the solution space for globular clusters by fixing the cluster mass and varying the other system parameters, such as temperature, cluster radius, energy, and so forth. To adapt this model for swarms, we use a fixed inverse temperature β, which describes isothermal globular clusters. In addition to being simple to implement (see SI), this choice is motivated by the observation that swarms of different sizes empirically display roughly the same amount of kinetic energy per midge [16,17]. For large globular clusters, a King model with fixed temperature gives a roughly constant kinetic energy per particle (independent of system size (Fig. S1)). With this assumption, we compute the density profiles predicted by the King model ( Fig. 1(a,b)), using Eqs. S15 and S16.
The King model predicts two distinct branches of solutions, which are clearly observed when we plot quantitative measures of the distributions such as the total mass M (Fig. 2(a), where we assume that all midges have the same mass), the density at the center ρ 0 (Fig. 2(b)), and the kurtosis (Fig. 2(c)) as functions of the overall swarm size R s . Here, R s is defined as the mean  1. (a, b) Two families of density profiles computed from the King model (Eqs. S15-S16) for different initial conditions, with shapes that are qualitatively similar to results from real swarms (c, d). (a) corresponds to the unstable branch, and (b) to the stable branch (see the discussion in the text). (c, d) Two families of density profiles measured in real swarms in two different laboratory setups (see text). (e) Density profiles of simulated swarms using the adaptive gravity model. Simulations were run using N = 12 and Rs = 16.2 (blue); N = 24 and Rs = 11.8 (red); N = 32 and Rs = 11.5 (green); and N = 48 and Rs = 12 (black). (f) Density profiles of simulated swarms including adaptivity and short-range repulsion. Simulations were run using N = 12 and Rs = 6 (blue); N = 24 and Rs = 6.7 (red); N = 32 and Rs = 8.2 (green); and N = 48 and Rs = 8.1 (black).
distance of a midge from the center of mass of the swarm: R s ≡ ∞ 0 r 3 ρ(r)dr/ ∞ 0 r 2 ρ(r)dr, where ρ(r) is the density profile. The two branches are termed "stable" and "unstable" (shown in blue and red, respectively, in Fig. 2(a-c)), and are distinguished by the sign of the heat capacity (positive or negative, respectively) of the cluster in the canonical ensemble [15]. It is found empirically that globular clusters reside on the unstable branch [15], presumably due to the strong destabilizing effects of "slingshots"-that is, anomalously high acceleration events that arise due to close encounters.
In Fig. 1(c,d) we plot the empirical density profiles measured for laboratory midge swarms. Details of the laboratory setup and measurement protocols are given in refs. [4] and [18]. The primary difference between the data in Fig. 1(c) and (d) is that the swarms in (c) were observed in a cubical laboratory enclosure measuring 91 cm on a side, while the enclosure for the swarms in (d) measured 122 cm on a side. Qualitatively, the density profiles from the swarms are similar in many aspects to those computed from the King model, though the agreement is not exact. To illuminate these similarities and differences further, we computed the same distribution measures for the swarms as for the King model, shown in Fig. 2(d-f). The density in the center of the midge swarm ρ 0 slightly decreases with increasing swarm size, which qualitatively fits the stable branch of the King model (where this decrease is much steeper). The total number of midges (as a proxy for the total mass) increases slowly with the swarm radius, a feature exhibited by the unstable branch of the King model (the stable branch has the opposite relation). The kurtosis of the midge swarms falls in a range of values between the King model branches, and is significantly larger (i.e., the swarms have heavier tails than they would if they were Gaussian) than the density profiles along the King model stable branch (which are very close to Gaussian).
To attempt to reconcile the model predictions with the empirical results for real swarms, we consider two modifications to normal Newtonian gravity. First, since particles in the simulations have a strong tendency to develop anomalously high accelerations due to slingshots and thus to evaporate from the cluster, it is common to include a softening parameter ǫ to the gravitational force [19], so that where C is a constant with units of (force × length 2 ), r i is the position vector for midge i, andr ij is the unit vector pointing from midge i to midge j. This "epsilon"-gravity modifies the force at short distances, up to | r i − r j | = O(ǫ).
In addition, we simulated our previously introduced adaptive-gravity model [10]. In this case, the effective force felt by midge i due to midge j is given by The logic underlying this model is that the strength of the signal received by each midge should be renormalized by the total buzzing noise due to all the other midges, to mimic the typical adaptive gain of biosensors [20]. R ad is the length scale over which this adaptivity occurs, and when r ij ≫ √ N R ad , such that the distance between a pair of midges is much larger than the adaptive range, the model reduces to epsilon-gravity (which itself reduces to Newtonian gravity when ǫ = 0).
Finally, we also simulated the adaptive-gravity model with an additional short-range repulsion between the midges. This effect was empirically observed and measured in real swarms [21], and was implemented in an earlier simulation scheme [10]. This type of effective interaction arises due to midges avoiding collisions while flying. Here we use a repulsive force with a cutoff to capture this effect, given by where we take L c ≈ 0.2 R ad . With these interaction laws, we considered several different swarm sizes R s , as defined above. We used R ad = 10 ≈ R s , placing us in the regime where adaptivity affects the vast majority of the midges in the swarm. To explore the spatial distribution of mass in the adaptive gravity model we turned to numerical simulations. We used a scheme for performing N -body dynamics originally developed by Aarseth (see Appendix 4.B of [22]), which allows for accurate numerical integration using fourthorder equations of motion. A complete description of the numerical method is given in ref. [19]. The acceleration of each mass in the simulation is computed by the direct summation of the forces due to the other N − 1 bodies following the form of the force law. The scheme was designed to work efficiently for up to N = 50, well within the range of typical midge swarms [4,18]. For further details about the simulation see the SI. Initially, particles were placed randomly in a simulation box of varying side lengths, in order to control the kinetic energy of the system. This was achieved (except for the epsilon-gravity system; see below) by varying the initial conditions of the simulations until the kinetic energy per particle was within 10% of the desired value in the long-time limit (Fig. S10). The initial velocities were zero, but under the influence of the forces the particles accelerated and eventually reached a quasi-stationary state (Figs. S2-S4) that exhibits the distinct non-Gaussian velocity statistics ( Fig. S5) observed in real swarms [4]. Note that we do not include explicit noise in the simulations, and the trajectories become ergodic due to the natural tendency of N-body gravitating systems to be chaotic.
We now compare the same quantitative features of the density profiles computed from the numerical simulations to those from the King model and the empirical swarm measurements. We start with the epsilon-gravity interaction (Eq. (1), using ǫ = 3.8, which was much smaller than the swarm size R s ; Fig. S12). Not surprisingly, some measures ( Fig. S12(b,c)) seem to agree well with the stable branch of the King model ( Fig. 2a-c), since close-encounter slingshots are greatly suppressed by the softened gravity. Other similarities include a decrease in the size R s and sharp decrease in ρ 0 for larger swarms. This force law therefore does not agree very well with the observed profiles of the midge swarms ( Fig. 2(d,e)), where the size increases with the number of particles and the density at the center shows a very weak decrease with size. Note that for epsilon-gravity, we were not able to keep the kinetic energy per particle constant when changing the overall number of particles in the swarm (Fig. S11).
Figures 1(e) and 2(g-i) show the results from adaptivegravity simulations (Eq. (2), using R ad = 10). Note that the swarms cover a much smaller range of sizes (R s ), compared to the real (Figs. 1(c,d) and 2(d-f)) or epsilongravity swarms (Fig. S12). This is due to the sharp increase in the density at the center when the number of particles increases combined with the fixed kinetic energy per particle, thereby maintaining a roughly constant R s . We find that ρ 0 decreases with increasing R s , similar to the observations in real swarms (as well as epsilon-gravity and the stable branch of the King model). The kurtosis is similar to the values seen in the real swarms. However, just as for epsilon-gravity and the stable branch of the King model, the biggest discrepancy with the real swarms is the relation between the number of particles (midges) and R s , which we find in the simulations to be a decreas-ing function (Fig. 2(g)) even though it is increasing in real swarms (Fig. 2(d)).
We therefore also tested the addition of a shortrange repulsion (Eq. 3) to the adaptive-gravity simulation (Figs. 1(f) and 2(j-l)). This additional ingredient makes ρ 0 roughly independent of R s while not changing the kurtosis significantly. Both measures are in reasonable agreement with the real swarm data, considering the small range of swarm sizes in the simulations. The primary change as a result of the short-range repulsion is the appearance of an overall increasing relation between the number of particles and swarm size (Fig. 2(j)). This feature, which is observed in the real swarms ( Fig. 2(d)), is absent from all the simulations that contain only longrange attractive interactions. This model therefore captures qualitatively all the main features of the density profiles of the midge swarms.
These results allow us to clearly delineate the properties of the midge swarm that are due to the long-range (adaptive) gravity interactions, namely the large kurtosis (heavy tails) and a slow decrease of the density at the swarm center for increasing swarm sizes. The shortrange repulsion is responsible for inflating the size (R s ) of swarms with increasing number of midges, which is counter to the behavior of purely long-range interactions (both epsilon and adaptive gravity).
So far we have considered the spatial density distribution of the swarm. We now turn to the spatial distribution of the velocities of the particles in the swarm, as another way to distinguish between the different models and compare to the observations from real swarms. We therefore computed the average speed of midges as a function of the distance r away from the center of mass of the swarm using the dataset from the larger midge enclosure [18]. As shown in Fig. 3(a), these speeds are essentially independent of position (as also noted in ref. [17]). Thus, the mean kinetic energy of a midge is also statistically uniform in space ( Fig. S8(a)).
The simulated speed profile for epsilon-gravity is shown in Fig. 3b. We see that the velocity decreases rapidly with increasing radius from the swarm center, in good agreement with the profiles calculated from the King model but contrary to the data for real swarms. By contrast, the velocity profile computed using adaptivegravity (Fig. 3(c)) is rather flat, with a small increase of speed for small r, presumably because the potential is too strongly softened in the swarm center. The velocity profile becomes even smoother when short-range repulsion is included (Fig. 3(d)), and is similar to the empirical results for real swarms. We find similar results for the standard deviation of the speed (Fig. S8).
Finally, we considered the relation between mean acceleration and speed for swarms of different sizes (Fig. S7) [17]. We find that epsilon-gravity fails to reproduce the observed collapse of the curves for all swarm sizes, while the adaptive gravity captures this relation extremely well. Taken together, our results provide further evidence that the interactions in insect swarms are well described by the adaptive-gravity framework, together with a short-range repulsion. This highlights that two key ingredients-long-range interactions and adaptivity-are essential. This model naturally captures many of the unusual properties of swarms, some of which have surprising similarities to globular clusters. Adaptivity, however, also induces additional effects. When the binding to the swarm is adaptive, the forces on one midge due to the other midges are not additive, and as a consequence the sum of the forces felt by all of the midge need not vanish, as it must for Newtonian gravity [10]. The center of mass of the swarm can therefore experience accelerations [23]. Such fluctuations have the potential to change fundamentally the characteristics of individual flight patterns; for example, Reynolds and Ouellette [23] showed that center of mass fluctuations allow for the emergence of Lévy flight patterns. Centre-of-mass movements may also help to stabilize insect swarms against environmental perturbations [24].
These results suggest that it may be fruitful in the future to consider modifications of the interactions within the King model as a general framework for describing the behavior of active systems with long-range interactions, whatever their form. The proper analytical treatment of a King model augmented with an adaptive-gravity in-teraction law remains as a challenge for future studies. Finally, future work on swarm modeling in particular should explore the effects of adding self-propulsion and stochasticity (noise), both of which are absent from the current models but certainly present for real insects.
We thank Tsvi Piran, Sverre Aarseth, and Tal Alexander for useful discussions.

The King Model
Most of the globular clusters are well-fitted by the King model [4], since it captures the basic structure and main observed features of them [2]. This model is one of a series of models that assume Maxwell-Boltzmann distribution of the velocities, although it is well-known that a self-gravitating system cannot be in thermodynamical equilibrium. This fact is a result of the long-range nature of the interactions, and the Boltzmann entropy which has no maximal value in an unbounded domain. Therefore self-gravitating systems, as part of their out-of-equilibrium state, have a strong tendency to evaporate over time. Since evaporation is a slow process, the system can be considered to be in a quasi-stationary state for times that are much shorter than the evaporation time, and the Maxwell-Boltzmann distribution can therefore serve as a good approximation.
In order to write a quasi-stationary expression for the distribution f ( r, v) (a function of the position and the velocity that does not depend on time) we start from the colisionless Boltzmann equation: where V = V (r) is a spherically symmetric gravitational potential. The density is defined by and it determines the gravitational potential through the Poisson equation where C is a constant with dimensions of mass · length 3 /time 2 which for gravitational force is the gravitational constant. A simple solution to this set of equations ((S1) and (S3)) is given by the Maxwell-Boltzmann distribution with a singular power-law behavior in the spatial part: where β is a positive constant which is identified with the inverse temperature. This solution is singular at the origin, unbounded and has an infinite mass since The singularity at the origin can be fixed by introducing regular boundary conditions: but it does not fix the asymptotic behavior. Numerical solutions of equations ((S1) and (S3)) with boundary conditions (S6) are solutions of isothermal spheres (since β is constant). Far from the center r ≫ (C ρ(0) β) − 1 2 the isothermal sphere solutions approach the singular solution (S4). Although the solution is not singular at the center, the solution is still unbounded and has an infinite mass (see for instance [5], p. 480-484).
The steady-state mass distribution that solves the Poisson-Boltzmann equations is singular, proportional to ρ(r) ∝ 1/r 2 . To find physical solutions, it is necessary to introduce a cut-off. The assumption in the King model is that a cluster cannot keep stars whose velocity exceeds a finite escape velocity. Therefore the velocity distribution should drop to zero at a finite limiting velocity. As a result, the cluster has a finite radius. Different cut-offs in the velocity distribution give different limiting radii. Observationally the finite boundary is set by the tidal forces of the galaxy with which the globular cluster is associated, that become dominant at a certain distance from the cluster center. The velocity distributions are simply the Maxwell-Boltzmann distributions minus a constant. Let us assume the following form at the origin: where α is a positive constant. This way the fastest particles v > v e (0) are subtracted from the distribution. The distribution of the velocities at r = 0 is defined by two parameters: β (determines the width of the distribution) and v e (0), the escape velocity.
According to Jeans' theorem any steady state solution of the collisionless Boltzmann equation is a function of the integrals of motion (see for instance [1], p. 220). Since energy is the integral of motion in a star cluster the distribution function can be written as a function of the energy E: We express the cut-off as an escape energy E e such that f = 0 for E ≥ E e . The escape energy corresponds to an escape velocity at a point r through: (S10) where α is a positive constant. Substituting the distribution function into (S2) and changing variables give us: (S13) where and R e is the radius of the cluster such that v e (R e ) = 0, and χ(R e ) = 0. In other words, any star that succeeds to reach R e , escapes the cluster and it is not a part of the cluster any more. When we substitute the density (S13) in the Poisson equation (S3) we get the fundamental equation of the King model with the boundary conditions where the rescaled dimensionless distance is ζ ≡ r(4 π C β ρ(0)) 1 2 . (S17) The set of King solutions for the density is a one parameter family of solutions which are conveniently parameterized by k,which is the rescaled potential at the center. Using the expression for the density (S13) and the fundamental equation (S15), we can write the mass in the following way: Note that here too, in the limit k → ∞, the singular solution (S4) is recovered (see [2]). In Fig. S1 we give the average kinetic energy (per unit mass) as a function of R s for various King solutions. The average kinetic energy (per unit mass) < E k > is given by: In the case of an isotropic distribution of velocities and spherical symmetry we get: Substituting (S10) and changing variables according to (S11) and (S12) we obtain: for r < R e . Integration by parts gives us: Repeating the steps in Eqs. (S21)-(S23) with additional integration over the spatial directions gives the average kinetic energy: The average velocity is given by and then we can write the dispersion normalized by the average velocity: FIG. S1. The average kinetic energy (per particle) as a function of Rs for a series of King solutions. The blue and the red lines stand for the stable and unstable branches in the canonical ensemble respectively. Note that in the limit of large clusters the average kinetic energy is approximately independent of the cluster size.

The Time Evolution of the Swarm in the Simulation
Here we show a typical oscillatory behavior as a function of time of the half-mass radius (Fig. S2) and the size of the swarm (mean distance from the center of the swarm) Fig. S3 for initial conditions that do not cause evaporation (over the time segment that is presented here). In Fig. S4 we show that the same swarm with different initial conditions evaporates over the same time segment (Fig. S4). Such cases were not considered as stable in this paper, and were ignored.

Rs Rs
FIG. S4. The size of the swarm (normalized by its mean value) as a function of simulation time for a swarm with adaptivity (R ad = 5) and 50 individuals. Here, in this example, different initial conditions create an evaporation at earlier time than the previous example in Fig. S3.
In the examples that we considered in the paper we did not include time segments similar to this one, that showed a significant evaporation.

The velocity and acceleration distributions
The distribution of a single component of the velocity in the simulation is given in Fig. S5 for various swarm sizes. In Fig. S6 we give the similar distribution of a single component of the acceleration. The structure of the tails is similar to the same distributions that were obtained from the laboratory data [3].

The relation between mean acceleration and speed
It was observed in [6] that if we take the mean acceleration of the right or left hemisphere of the swarm separately, it's absolute value is a monotonously increasing function of the speed. The experimental plot from [6] for various sizes of swarms is reproduced here in Fig. S7d where the speed is normalized with respect to its mean value. Comparing with the simulations results (for ǫ-gravity in Fig. S7a, adaptivity in Fig. S7b and adaptivity with short-range repulsion in Fig. S7c) for swarms of different sizes, we see that in the cases with adaptivity (with and without repulsion) all the graphs collapse at low velocities to have approximately the same acceleration, when in the ǫ-gravity case there are various values of acceleration at low velocities, and it does not agree with the experimental data in Fig. S7d. One can also see that with repulsion there is more noise at high velocities. The fitting that we see here between adaptivity and the experimental data suggests that adaptive-gravity naturally explains the dependence of the acceleration on the mean velocity, without the need to invoke an additional explicit velocity-dependent forces [6]. The kinetic energy of each midge in the laboratory is approximately constant and does not depend on the density or the size of the swarm. In the simulation, since we are interested only in the mutual forces, we do not require constant speed for the particles, and therefore the kinetic energy per individual is different for various numbers of particles and initial conditions. We chose initial conditions that give us approximately the same kinetic energy per individual in the case of simulations with adaptivity and adaptivity with repulsion (up to 10%). Without adaptivity we could not change significantly the kinetic energy by varying the initial conditions and it was determined mainly by the number of particles in the simulation.
The simulation with softened gravity ("epsilon-gravity") In order to avoid high accelerations due to close encounters ("slingshots") and thus evaporation of the cluster we included a softening parameter ǫ to the gravitational force (Eq. 1). The particles in the simulation still developed high accelerations in close encounters relative to simulations with adaptivity since adaptivity itself contributed to the decrease of accelerations by ǫ 2 /(R 2 ad + ǫ 2 ), and in order to compare the simulations in the two cases (with and without adaptivity) we have to reduce the accelerations in the "ǫ-gravity" case by this constant factor. The results in   were obtained this way. In Fig. S12(a) density profiles for various sizes of clusters are given (where ǫ 2 = 15) and in Fig. S12(b)-(d) we show three characteristic quantities of the clusters: the number of particles, the density at the center and the kurtosis. Here and in all the paper we truncated the tails at r ∼ 3 R s in the calculation of the kurtosis. The standard deviation of the speed normalized by the average speed in the swarm vs. the distance from the center of mass normalized by the swarm size Rs. We see that for r > 2 Rs the velocities become up to twice higher, the standard deviation of the velocities becomes twice lower and the density is less than 10 −5 of the density at the center (Fig. 1). Therefore this region probably describes particles that were repelled out of the swarm at high velocities (higher than the escape velocity). FIG. S10. The kinetic energy as a function of the size of the swarm Rs for cases with adaptivity (Fig. 2). Here the size of the swarm does not change much but the number of individuals change from 12 to 48 (Fig. 2).