Understanding the thermodynamic properties of insect swarms

Sinhuber et al. (Sci Rep 11:3773, 2021) formulated an equation of state for laboratory swarms of the non-biting midge Chironomus riparius that holds true when the swarms are driven through thermodynamic cycles by the application external perturbations. The findings are significant because they demonstrate the surprising efficacy of classical equilibrium thermodynamics for quantitatively characterizing and predicting collective behaviour in biology. Nonetheless, the equation of state obtained by Sinhuber et al. (2021) is anomalous, lacking a physical analogue, making its’ interpretation problematic. Moreover, the dynamical processes underlying the thermodynamic cycling were not identified. Here I show that insect swarms are equally well represented as van der Waals gases and I attribute the possibility of thermodynamic cycling to insect swarms consisting of several overlapping sublayers. This brings about a profound change in the understanding of laboratory swarms which until now have been regarded as consisting of non-interacting individuals and lacking any internal structure. I show how the effective interactions can be attributed to the swarms’ internal structure, the external perturbations and to the presence of intrinsic noise. I thereby show that intrinsic noise which is known to be crucial for the emergence of the macroscopic mechanical properties of insect swarms is also crucial for the emergence of their thermodynamic properties as encapsulated by their equation of state.


I: No evidence for long-range interactions in laboratory insect swarms
Tetrahedra are the minimum configuration required to describe the 3-dimensional movements of midges (since 3 midges will always lie in a plane). Here I show that the statistical properties of these tetrahedra mirror expectations for Gaussian, independent individual positions; thereby bolstering the analysis of Puckett et al. [2014] who reported the acceleration measurements show a clear short-range repulsion but do not provide conclusive evidence of long-range interactions between individuals and their nearest neighours.
Following Biferale et al. [2005] the dynamics of the tetrahedra can be analysed by introducing the relative velocities: 1 = ( 2 − 1 )/√2, 2 = (2 3 − 2 − )/√6, 3 = (3 4 − 3 − − )/√12. The evolution of a tetrahedron is given by The velocities of the tetrahedra tend to have one extensional component and  one compressional component: exactly mirroring expectations for Gaussian, independent individual positions (Fig. S3). Note also that the alignment of the eigenvectors associated with g1 and k1 mirrors expectations for Gaussian, independent individual positions; indicating that strong velocity differences are not preferentially associated with intense elongations approximately in the same direction.
Following Biferale et al. [2005] the dynamics of the tetrahedra were further examined by determining the average times, 〈 ( )〉 for the eigenvalues, g1, g2, g3 to double (Figs.S4a and S5a). Midge swarms exhibit a scaling regime, 〈 ( )〉~1 . The stochastic model of Reynolds et al. [2017] for the trajectories of non-interacting individuals, within which velocities and positions evolve jointly as Markovian process, predicts such a scaling regime, albeit with a scaling exponent of 2 rather than 1. The presence of the scaling regime and the self-similarity are more evident after applying simple multiplicative factors on the g-axis ( Fig. S4b and Fig.   S5b). The presence of a range where the doubling times for different eigenvalues are the same is equivalent to stating that the typical shape of the tetrahedra is preserved while its size increases. Stochastic modelling suggests that the difference between the observed and predicted scaling exponents cannot be attributed to short-range repulsion as reported by Puckett et al [2014]. Nor can it be attributed to interspecific variability of model parameters.
These do not change the prediction for the scaling exponent. It may instead to be attributed to having presupposed that the acceleration autocorrelation approaches a δ function at the origin, corresponding to an uncorrelated (white noise) component in the acceleration, and hence a Markov process. When the white-noise in the model of Reynolds [2017] is replaced by noise with exponential autocorrelation, i.e., when the white-noise is replaced by an OUprocess, the predicted scaling exponent can in brought into line with observations ( Fig. S6).
In such models, accelerations, velocities and positions are jointly Gaussian and evolve collectively as a Markov process . Predicted scaling exponents depend on model parameters, suggesting that the scaling is species dependent.
'Perceived velocity gradients' were examined following the approach reviewed in Pumir and 2 . This is because the eigenvalues of M are given by the roots of 3 + + = 0 . The quantity Δ = 4 3 + 27 2 separates purely real solutions ( < 0 ), physically interpreted as the flow being locally non-rotational, from one real and a complex conjugate pair of roots ( > 0), corresponding to a flow that is locally rotating The distribution of the two invariants, R and Q, mirror predictions from stochastic models for noninteracting individuals (Fig. S7). In accordance with observations, the model predicts that the flow is locally rotating 88% of the time. The ≠ 0 indicating that the flows are compressible.
The above approaches can in principle be extended to include accelerations, e.g. via and via perceived gradients in accelerations, −1 . It is, however, potentially problematic to compare observations with expectations for independent individuals derived from 1 st -order stochastic trajectory models. This is simply because accelerations per se, as opposed to average accelerations, are not well-defined model quantities. Nonetheless, plots of 1 ′ against 2 ′ etc. (ala Fig. S3) where 1 ′ ≥ 2 ′ ≥ 3 ' are the eigenvalues of closely match model predictions.
. Figure S1. The tetrahedra within the midge swarms tend to have almost twodimensional, strongly elongated geometries mirroring expectations for Gaussian, independent individual positions. Distributions of the dimensionless shape parameters = / 2 where = √ 1 + 2 + 3 is the radius of gyration. Distributions are shown for swarms containing on average 94 (Ob1) and 19 (Ob18) individuals. Data are taken from Sinhuber et al. [2019]. Shown for comparison are simulation data for Gaussian, independent individual positions. Covariances 〈 〉 also closely mirror expectations for Gaussian, independent individual positions.

II: Indirect speculative evidence of confinement to horizontal slabs
Slab confinement can be assessed quantitively following the methodology of Pumir et al.
[2000] and Biferal et al. [2005] since 4 points along an individual's trajectory will form a tetrahedron. At short times, individual trajectories are, as expected, colinear (as 1 ≫ 2 , 3 ) but at longer times they are coplanar (as 1 ≈ 2 ≫ 3 ); and as expected mirror expectations for Gaussian, independent positions (Figs. S8). The tendency to form almost two-dimensional structures therefore has mostly an "entropic" origin: indeed, there are many pancake-like tetrahedra (very small I3) when positions are decorrelated. At long times the proportions of midge trajectories elongated in the x-, y-and z-directions are approximately 5:5:1, indicating that confinement is mostly on horizontal planes.
In this regard, it is interesting to note that Cavagna et al. [2017] showed that wild swarms appear to belong to a novel dynamical universality class: as correlations exhibit scaling behaviour that is atypical 3-dimensional systems. The observed scaling is distinctly different from the scaling that characterizes the 3-dimensional Vicsek et al. [1995] model of swarming self-propelled particles and much more like the scaling that characterizes the 2-dimensional form of this model [Baglietto and Albano 2008]. That is, the apparent novelty of the scaling may be attributed to slab confinement.

Figure S8. Individual trajectories are colinear at short times and coplanar at long times.
Distributions of the dimensionless shape parameters = / 2 where = √ 1 + 2 + 3 is the radius of gyration. Data are taken from Sinhuber et al. [2019]: swarm Ob1 which on average contains 94 individuals.

for arbitrarily defined density profiles -associated mean accelerations are non-physical
The density profiles of laboratory swarms are skewed perhaps because of the energetic demands of flying upwards and/or because they are bounded at ground level [Kelley and Ouellette 2013]. Such density profiles could for argument sake be crudely represented as Here I show that the number of Gaussians in the best representations of an such arbitrarily defined density profile does not increase monotonically with swarm size, contrary to the analysis of natural swarms (Fig. 3). I also show that the mean accelerations associated with such mono-modal distributions are non-physical and incompatible with observations. These results are consistent with the claims in the main text that the multi-Gaussian fits to the swarm density profiles are indicative of actual structure and with swarm growth proceeding slab by slab.
As for natural swarms (see main text) the best representation of the putative gamma distributions in terms of multi-Gaussians was found using maximum likelihood method and the best multi-Gaussian distribution was identified objectively using the Akaike information criterion. Typically, the number of Gaussians in the best representation does not increase monotonically with swarm size (   Velocities are taken to be Gaussian distributed with mean zero and unit variance.

Model 1
In this model, as the population size of the swarm increases, existing slabs move apart to accommodate more individuals and in doing so become less stable (become noisier). This noise creates new centres of attraction, spawning 'embryonic' slabs and is thereby constructive. The swarm population then redistributes, re-establishing equilibrium, and the intrinsic noise diminishes. The cycle then repeats. In this regard it is interesting to note that swarms are analogous to self-gravitating systems [Okubo 1986, Reynolds 2019a and references therein] and that the model of slab formation is analogous to Hoyle's Steady State Model. Near constancy of average density is a hallmark of laboratory swarms [Kelley and Ouellette 2014]. Motivation for the approach comes from: the association of tensile strength with intrinsic noise [Reynolds 2019b]; the subsequent identification of intrinsic noise in quiescent swarms ; the presence of an embryonic-like slab located about midway between two more prominent slabs in laboratory data ( Fig. 1b main text).
A swarm with equi-spaced, equi-size slabs can be represented in simple stochastic models by a periodic potential ( ) = 0 ( ) (S1) since the associated vertical density profile is unjustified (and although putative phases were not identified), the above analysis again raises the spectre of near criticality because the number of slabs is predicted to be near to a doubling transition; more generally near a slab-creation transition.

Model 2
In this model, increases in acoustic interference that accompany increases in population size are offset by increases acoustic clustering [Aldersley et al. 2017 )) (S6) (Fig. S12). Intrinsic noise, which here acts on the slabs rather than on the confining potential, causes the swarm size to grow exponentially with the number of slabs (Figs. S13); thereby providing a mechanism for the regularization of the average density. By inducing acoustic clustering, acoustic interference can thus cause explosive growth of the swarm.

Adaptation
Finally, it is worth remarking that 'adaptation' may also be implicated in slab formation. As noted by Gorbonos et al. [2016] midges interact acoustically and the perception of sound may not be fixed but instead adapts to the local sound intensity so that acoustic intensity drops when there is a strong background noise; preventing damage and saturation of the sensory organs. Gorbonos et al. [2016] utilized the similarity in form between the decay of acoustic and gravitational sources to construct an 'adaptive gravity' model of swarming that agrees well with experimental observations of laboratory swarms. In the outskirts of a swarm, where the background noise is weak, the effects of adaptivity on the cohesive forces that bind the swarm together are weak. But in the centres of large swarms where the background noise is strong, adaptivity can sufficiently weaken the cohesive force that holds the swarm together. This led Gorbonos et al. [2016] to speculate on the appearance of a maximal swarm size, beyond which they become unstable and split. Here, this possibility is examined within the context of stochastic models of the trajectories of swarming insects.
Midges are, to good approximation, bound to the centre of the swarm by an effective mean acceleration (mean force) that grows linearly with distance from the swarm centre [Okubo 1986, Kelley and. This is consistent with the results of stochastic model which predicts that the effective mean acceleration is given by = − In this case the equilibrium density profiles is given by where N is a normalization factor. Adaptation thereby results in a near flat density profile within the core of the swarm (Fig. S14a). This increases the likelihood that the swarm will split because fluctuations in number density are now more likely to result in the spontaneous appearance of clusters of individuals away from the centre of the swarm and so result in the spontaneous appearance of new centres of attraction (Fig. S14b).
Note also that according to Eqn. S8, the physical size of the swarm increases as the strength of the adaption increases, i.e., increases as the population size and so the background noise increases. This is consistent with  who reported that adaption

Figure S14a
). Adaptation is predicted to flatten density profiles within the cores of swarms. b). Adaptation increases the potential for splitting. The likelihood that there are more individuals located outside the core, | | > 1 than within it as a function of Λ 2 . Results are shown for a swarm of size = 1 containing 100 (red line) and 1000 (blue line) individuals.
Positions were drawn at random from the distribution, Eqn. S8.

V: Outline of stochastic model formulation
The formation of stochastic models for the trajectories of swarming insects was originally presented in  and is, here, outlined for completeness. The simplest such models, like Eqn. 4 in the main text, assume that the positions, z, of individual insects can be described by the stochastic differential equation

VI: Interpreting the pressure-density relation for the 'vapour phase'
Laboratory insect swarms consist of a core 'condensed' phase surrounded by a dilute 'vapour' phase [Sinhuber and Ouellette 2017]. These two phases coexist in equilibrium but have distinct macroscopic properties. In the condensed phase, pressure, P, is approximately proportional to density, n, (i.e., its 'isothermal') whereas in the vapour phase ∝ with 1/ 2   [as predicted by the model of Reynolds et al. 2018a]. The vapour phase relationship does not have a thermodynamic analogue, but it does have a gravitational analogue, corresponding to descriptions of interstellar dust clouds [Shu et al. 1972, Viala 1972, as prefigured by . The core is analogous to the globular cluster [Reynolds 2018b], as anticipated by Gorbonos et al. [2016]. The association of the vapour phase with interstellar dust clouds adds to a long-standing ] and productive analogy between insect swarms and self-gravitating systems [Reynolds 2018b, Gorbonos et al. 2020. In this regard it is interesting to note the smallest swarms (with less than 10 individuals) are predicted to lack condense cores and be entirely analogous to interstellar dust clouds ].

VII: Critical Damping
Here I show that laboratory swarms of Chironomus riparius midges are critically damped and I show that critical damping arises freely in simple 3-dimensional stochastic models.
Laboratory swarms of Chironomus riparius midges [Sinhuber et al. 2019] are found to be either weakly overdamped (as the negative lobe is barely evident) or weakly underdamped (as there is no evidence of a secondary positive lobe on the velocity autocorrelation functions, i.e., the midges just fail to overshoot) (Fig. S15). Moreover, estimates for the frequencies of oscillatory modes are typically close to zero. This suggests that the laboratory swarms are close to being critically damped, poised between being over-and underdamped; contrary to the observations of laboratory swarms of Anarete midges ] and observations of wild swarms of Anopheles gambiae [Butail et al. 2013]. Consequently, the Chironomus riparius midges tend to return to their equilibrium positionthe centre of the swarmin the minimum time; typically, just failing to overshoot and not making single oscillations. A precursor of this critical damping can be found in  who showed that mean-squared displacements saturates when 〈 2 〉 ≈ 2 ; i.e., when the average midge has reached the edge of the swarm.
Critical damping seemingly requires fine tuning. It occurs in 1-dimensional models for swarms with Gaussian position and velocity statistics when 2 2 = 2 ]. 3-dimensional models do not appear to be not analytically tractable. Nonetheless, the results of numerical simulations indicate that overshooting (resulting in a secondary positive lobe in the velocity autocorrelation) is generally absent or weak in the 3-dimensional versions of  model [Reynolds et al. 2017] when 'spin' (a preferred sense of rotation of individual trajectories) is the absent or is sufficiently weak (Fig. S16). [This is also true of second-order models wherein accelerations, velocities and positions rather than just velocities and positions evolve jointly as Markovian processes. And it is true of simulated spinless swarms as they relax back to their equilibrium positions after being displaced]. Stochastic models with spin produce looping trajectories [Borgas et al. 1997]. Spin has been detected in laboratory swarms of Chironomus riparius midges [Reynolds 2019, see also  but it may be stronger in laboratory swarms of Anarete midges ] and in wild swarms of Anopheles gambiae [Butail et al. 2013]. Spin suppresses dispersion, thereby promoting cohesiveness but impedes a swarms' return to its' equilibrium form after being perturbed. This may account for the observations of Butail et al. [2013] who reported that the frequencies of oscillation increased with mean wind speed: environmental disturbances may excite the spin degrees of freedom because spin pushes the swarm into a more robust state. Robustness may be further enhanced by the formation of transient, local order (synchronized subgroups) [Butail et al. 2013, Shishika et al. 2014].
Critical damping (or nearly so) also arises freely in 2-dimensional stochastic models and is evident at the population level in two species of zooplankton, Daphnia and Temora, that swarm around light shafts [Banas et al. 2003]. Individual velocities decorrelate almost entirely in less than one period of the harmonic attractive force.

Figure S16. Near critical dampening arises freely in 3-dimensional stochastic models
when spin is sufficiently weak , Reynolds et al. 2017]. It occurs for many values of the model parameters. Ouellette [2016], they expand rather than contract in the perpendicular direction. This suggests that swarms do not become inherently weaker when stretched because they thicken in response to a tensile force, rather than becoming thinner like positive-Poisson-ratio material.

VIII: Anticipated emergence of auxetic material-like properties
Auxetic material-like properties are predicted to arise because stretching a swarm in one direction will, by virtue of diluting the population, tend to increase the fluctuations in the resultant centrally attractive force. Such fluctuations (intrinsic noise) will tend to further increase the overall size of swarm (see main text).

IX: An illustrative example of a further application: modelling of honeybee drone congregation areas and flyways
Males of many insect species form dense, lek-like aerial swarms, above visual cues known as swarm markers. These often maintain a relatively stable size and shape even as individuals leave and others arrive, leading Sullivan [1981] to hypothesise that individual males move between adjacent swarms. Woodgate et al. [2021] provide the best evidence for a mating strategy in which individuals, honeybee drones, travel between multiple aerial leks whose locations are fixed. Woodgate et al. [2021] used harmonic radar technology [Riley et al. 1996] to record the flight paths of individual drones. They found clear evidence that drones favour certain areas in which to congregate, and that these areas (drone congregation areas, DCAs) were stable over two sequential years. Surprisingly, drones often visit multiple potential lekking sites within a single flight and take shared flight paths, 'flyways' between them. Flights between such sites, and between such sites and the hive are relatively straight. Once inside a drone congregation area, drones display convoluted, looping flight patterns those velocity and acceleration statistics closely resemble that of midges   [Okubo 1986, Kelley and Ouellette 2013] have to good approximation Gaussian density profiles. Consequently, the overall density profile for 2 DCAs can be represented by a bi-Gaussian distribution, 2 2 )) (S13) Woodgate et al. [2021] also reported that velocities within the DCA's are to good approximation Gaussian distribution. A 1-dimension stochastic trajectory simulation model that is consistent with a bi-Gaussian density profile and Gaussian velocity statistics is given in  ] that takes account of intrinsic noise may, for example, explain why the DCAs appear to be in tension [Woodgate et al. 2021] and why DCAs, like some midge swarms [Poda et al. 2019], are not located directly over swarm markers [Woodgate et al. 2020].
One big question raised by the study of Woodgate et al. [2021] is why do drones switch between DCAs when lekking vertebrates are faithful to a single location? Although the switching between DCAs is here predicted to be accidental it may nonetheless be advantageous: either because when queens are rare in the environment, a regular patrol between multiple locations gives a better chance of finding a queen than staying in a single place; or because gatherings of drones have properties that are attractive to queens and that drones hunt around for a particularly good DCA to join. This suggests that there could be selection for maintaining switching rather than selection for switching. Nonetheless, switching is only predicted to occur when the DCAs are sufficiently close together. This is consistent with Feugère et al. [2021] who reported that remote mosquito swarms are acoustically isolated.