Natural swarms in 3.99 dimensions

Abstract

The renormalization group is a key set of ideas and quantitative tools of statistical physics that allow for the calculation of universal quantities that encompass the behaviour of different kinds of collective systems. Extension of the predictive power of the renormalization group to collective biological systems would greatly strengthen the effort to put physical biology on a firm basis. Here we present a step in that direction by calculating the dynamical critical exponent z of natural swarms of insects using the renormalization group to order ϵ = 4 − d. We report the emergence of a novel fixed point, where both activity and inertia are relevant. In three dimensions, the critical exponent at the new fixed point is z = 1.35, in agreement with both experiments (1.37 ± 0.11) and numerical simulations (1.35 ± 0.04). Our results probe the power of the renormalization group for the quantitative description of collective behaviour, and suggest that universality may also play a decisive role in strongly correlated biological systems.

Main

Collective behaviour is found in a great variety of biological systems, from bacterial clusters and cell colonies, to insect swarms, bird flocks and vertebrate groups. A unifying ingredient, which also provides an insightful connection with statistical physics, is the presence of strong correlations. The correlation length, ξ, is often substantially larger than the microscopic scales^{1,2,3,4,5,6,7}, and in some instances ξ grows with the system size, giving rise to scale-free correlations^1,6,8. In the case of natural swarms of insects, a second key hallmark of statistical physics has been verified, namely dynamic scaling⁹. This is noteworthy, as dynamical scaling entangles spatial and temporal relaxation into one law, known as critical slowing down¹⁰: the collective relaxation time grows as a power of the correlation length, τ ~ ξ^z, thus defining the dynamical critical exponent, z. Strong correlations and scaling laws are the two essential prerequisites of the renormalization group (RG)^11,12. By coarse-graining short-wavelength fluctuations, the parameters of different systems flow towards few fixed points ruling their large-scale behaviour; RG fixed points therefore organize into universality classes the macroscopic behaviour of strongly correlated systems, thus providing parameter-free predictions of the critical exponents. The emergence of scale-free correlations and scaling laws also calls for an exploration of the RG path in collective biological systems.

In the broader field of active matter¹³, RG is already a key tool. The pioneering hydrodynamic theory of Toner and Tu^14,15 has been studied via the RG both in the polarized^16,17 and near-ordering phase^18,19, with applications in systems with nematic or polar order^20,21,22. RG has also been employed to study motility-induced phase separation^23,24, active membranes²⁵, bacterial chemotaxis²⁶ and cellular growth²⁷. Direct comparisons with experiments are few: although the exponent of giant number fluctuations in d = 2 (refs. ^15,17) has been confirmed in experiments on vibrated polar disks²⁸, in ref. ²⁹ the exponents of the Vicsek model in the ordered phase were found to be incompatible with those conjectured by Toner and Tu¹⁴. Other RG exponents have been checked in numerical simulations^30,31,32,33. Comparisons with biological experiments are scarcer. Experiments studying giant number fluctuations in swimming filamentous bacteria displaying long-range nematic order³⁴ found an exponent in disagreement with RG predictions of active nematic²¹ and polar¹⁷ systems. To the best of our knowledge, there is yet no successful test of an RG prediction against experiments on living active systems.

In this article we apply the RG to the dynamics of insect swarms. Swarms of midges in the field are near-critical, strongly correlated systems⁸ with short-range interactions⁴, obey dynamic scaling⁹ with an experimental exponent of z_exp = 1.37 ± 0.11, which is substantially smaller than the value z ≈ 2 of standard ferromagnets³⁵. This large gap indicates that fundamental new physics is required. Although the relation τ ~ ξ^z is not merely a dispersion law, smaller z nevertheless suggests that fluctuations are more swiftly transported across the system. Hence, the effort to match theory with experiments in natural swarms must look for more efficient mechanisms of information transfer. Activity is the first obvious candidate, as it allows fluctuations to propagate not only thanks to the inter-individual interaction, but also through the self-propelled motion of the particles¹⁴. Incompressible, near-critical active matter was first studied in ref. ¹⁸, where an RG analysis found that activity lowers z from 2 to 1.73. This was an important step towards bridging the gap between experiments and theory in natural swarms, although the chasm with the experimental exponent remains substantial.

The second ingredient known to foster information propagation is inertia. Behavioural inertia in the rotations of individual velocities was first introduced to explain the propagation of collective turns in bird flocks^36,37, but experiments also found clear evidence of underdamped inertial relaxation in natural swarms of midges⁹. At the general level, inertial dynamics stems from the existence of a reversible coupling between the primary field (playing the role of the generalized coordinate) and the generator of the symmetry (playing the role of the generalized momentum); in the case at hand, the symmetry is the rotation of the primary field, so we call ‘spin’ its generator. In the absence of explicit dissipation, this coupling leads to global conservation of the spin, a conservation law that is known—at equilibrium—to significantly decrease the dynamical exponent, from the z ≈ 2 of standard ferromagnets (model A in the classification of ref. ³⁵) to z = 1.5 in superfluid helium and quantum antiferromagnets (models E/F and G of ref. ³⁵). Note that, in model A, two-loop corrections to the one loop z = 2 value are very small³⁵, but in models E/F/G, z = d/2 is an exact result for all d ≤ 4 (ref. ³⁸). Overall, this scenario suggests that the combined effect of activity and inertia may account for the experimental exponent of natural swarms. In this article we perform an RG study of this theory, and find z = 1.34(8) in d = 3, a value in agreement with experiments on real swarms, z_exp = 1.37 ± 0.11, and also numerical simulations, z_sim = 1.35 ± 0.04. The RG result is a parameter-free prediction, with no input beyond the information that both activity and inertia must be part of the theory.

A hydrodynamic theory of active matter with reversible inertial couplings requires three fields—velocity, spin and density^39,40—making the calculation technically unfeasible. To make progress, following ref. ¹⁸, we eliminate the density field by imposing incompressibility, ∇ · v(x, t) = 0. Beyond its technical inevitability, this is a reasonable physical assumption. In compressible active systems the transition is first-order⁴¹, a framework that would make RG pointless and would rule out scaling. However, dynamic scaling is observed in natural swarms, suggesting a more complex scenario: in the absence of density fluctuations, the transition becomes second-order, and recent studies^42,43,44 suggest the existence of a crossover from a finite-size regime, where density fluctuations are weak and second-order physics is observed, to a infinite-size regime, where density fluctuations dominate, rendering the transition first-order. Weak density fluctuations and scaling laws in natural swarms⁹ suggest then that the dynamics of these finite-size systems can be studied through an incompressible second-order theory.

The polarization field, ψ, is defined by the relation, v(x, t) = v₀ψ(x, t), where v₀ is the microscopic speed. Having v₀ as an explicit parameter is useful so as to take the zero-activity limit, v₀ → 0, and compare it with equilibrium calculations^45,46,47. The generator of the rotations of ψ(x, t) is the spin, s(x, t). In d = 3 the spin is a vector; however, incompressibility requires ψ to have the same dimension as space, and because RG entails an expansion in powers of ϵ = 4 − d, we need a generic spin in d dimensions, which is a d × d anti-symmetric tensor. Reversible inertial dynamics arises from the Poisson brackets³⁵:

$${\{{s}_{\alpha \beta }\left({{{\bf{x}}}},\,t\right),{\psi }_{\gamma }\left({{{{\bf{x}}}}}^{{\prime} },\,t\right)\}=2g\,{{\mathbb{I}}}_{\alpha \beta \gamma \rho }{\psi }_{\rho }({{{\bf{x}}}},\,t){\delta }^{(d)}({{{\bf{x}}}}-{{{{\bf{x}}}}}^{{\prime} })},$$

(1)

stating that ${{s}_{\alpha \beta }\left({{{\bf{x}}}}\right)}$ is the generator of the rotational symmetry, thus leading to conservation of the total spin, ${{S}_{\alpha \beta }(t)=\int{\rm{d}}^{d}x\,{s}_{\alpha \beta }\left({{{\bf{x}}}},\,t\right)}$. The crucial constant g is the reversible coupling regulating this symplectic structure³⁵. Finally, ${{{\mathbb{I}}}_{\alpha \beta \gamma \nu }=({\delta }_{\alpha \gamma }{\delta }_{\beta \nu }-{\delta }_{\alpha \nu }{\delta }_{\beta \gamma })/2}$ is the identity in the space of s_αβ.

The dynamical field theory we study combines the irreversible off-equilibrium hydrodynamic approach of Toner and Tu^14,15,18 with the reversible conservative structure used to describe superfluid helium and quantum antiferromagnets (models E/F and G of ref. ³⁵). The equations of motion are given by

$${D}_{t}{\psi }_{\alpha }=-{{\varGamma }}\frac{\delta {{{\mathcal{H}}}}}{\delta {\psi }_{\alpha }}+g\,{\psi }_{\beta }\frac{\delta {{{\mathcal{H}}}}}{\delta {s}_{\alpha \beta }}-{\partial }_{\alpha }{{{\mathcal{P}}}}+{\theta }_{\alpha },$$

(2)

$${D}_{t}{s}_{\alpha \beta }=-{{{\varLambda }}}_{\alpha \beta \gamma \nu }\frac{\delta {{{\mathcal{H}}}}}{\delta {s}_{\gamma \nu }}+2g\,{{\mathbb{I}}}_{\alpha \beta \gamma \nu }{\psi }_{\gamma }\frac{\delta {{{\mathcal{H}}}}}{\delta {\psi }_{\nu }}+{\zeta }_{\alpha \beta },$$

(3)

where the material derivatives are defined as

$$\begin{array}{lll}{D}_{t}{\psi }_{\alpha }&=&{\partial }_{t}{\psi }_{\alpha }+{\gamma }_{v}{v}_{0}\,{\psi }_{\nu }{\partial }_{\nu }{\psi }_{\alpha },\\ {D}_{t}{s}_{\alpha \beta }&=&{\partial }_{t}{s}_{\alpha \beta }+{\gamma }_{s}{v}_{0}\,{\psi }_{\nu }{\partial }_{\nu }{s}_{\alpha \beta }.\end{array}$$

(4)

Activity breaks Galilean invariance¹⁵, so that the couplings γ_v and γ_s can be different from 1 and from each other. ${{{\mathcal{H}}}}$ is the classic Landau–Ginzburg coarse-grained Hamiltonian¹⁴:

$${{{\mathcal{H}}}}={\int{\rm{d}}^{d}x\left[\frac{1}{2}{\partial }_{\beta }{\psi }_{\alpha }{\partial }_{\beta }{\psi }_{\alpha }+\frac{r}{2}{\psi }^{2}+\frac{u}{4}{\psi }^{4}+\frac{1}{2}{s}^{2}\right]}.$$

(5)

The ψ-dependent part of ${{{\mathcal{H}}}}$ is the standard alignment interaction, while s²/2 is the ‘kinetic’ term³⁵. Because natural swarms have scale-free correlations⁸, we will perform the calculation on the critical manifold, r = 0.

The terms proportional to g in equations (2) and (3) are the reversible forces, characterizing inertial dynamics. Instead of acting directly on the polarization, the alignment force ${\delta {{{\mathcal{H}}}}/\delta \psi}$ acts on the spin, which in turns rotates ψ (ref. ³⁷). On the other hand, the terms proportional to the kinetic coefficients Γ and Λ represent the irreversible forces, responsible for relaxation. The pressure ${{{\mathcal{P}}}}$ enforces incompressibility, which in k-space translates into k_αψ_α = 0, or ${P}_{\alpha \beta }^{\perp }{\psi }_{\beta }={\psi }_{\alpha }$, with the projector ${P}_{\alpha \beta }^{\perp }={\delta }_{\alpha \beta }-{k}_{\alpha }{k}_{\beta }/{k}^{2}$. To respect this constraint, one must project the dynamic equation for ψ, equation (2), and its noise correlator¹⁸:

$${\langle {\theta }_{\alpha }\left({{{\bf{k}}}},\,t\right){\theta }_{\beta }\left({{{{\bf{k}}}}}^{{\prime} },\,{t}^{{\prime} }\right)\rangle ={(2\uppi )}^{d}2\tilde{{{\varGamma }}}\,{P}_{\alpha \beta }^{\perp }\,{\delta }^{(d)}\left({{{\bf{k}}}}+{{{{\bf{k}}}}}^{{\prime} }\right)\delta \left(t-{t}^{{\prime} }\right)}$$

where ${\tilde{{{\varGamma }}}\ne {{\varGamma }}}$ out of equilibrium. Notice that, if we write a general anisotropic form of this kinetic coefficient, ${{{\varGamma }}}_{\alpha \beta }={{{\varGamma }}}^{\perp }{P}_{\alpha \beta }^{\perp }+{{{\varGamma }}}^{\parallel }({\mathbb{I}}-{P}_{\alpha \beta }^{\perp })$, only Γ^⊥ survives the projection of equation (2), so that effectively Γ = Γ^⊥ in our notation (and similarly for ${\tilde{{{\varGamma }}}}$). Further non-diagonal forms of the kinetic coefficient may be conceivable, but we do not study them here, as they would make the calculation too intricate and go beyond the scope of this article. On the other hand, because equation (3) is not projected, anisotropic corrections to Λ in k-space are generated by the RG⁴⁷, so it is convenient to assume from the outset a general anisotropic form:

$${{{\varLambda }}}_{\alpha \beta \gamma \nu }={\left({\lambda }^{\perp }{{\mathbb{P}}}_{\alpha \beta \gamma \nu }^{\perp }+{\lambda }^{\parallel }{({\mathbb{I}}-{{\mathbb{P}}}^{\perp })}_{\alpha \beta \gamma \nu }\right){k}^{2}},$$

(6)

where ${{\mathbb{P}}}_{\alpha \beta \gamma \nu }^{\perp }$ is the projector in the anti-symmetric space⁴⁷. The fact that this kinetic coefficient is proportional to k² means that the irreversible terms also conserve the total spin. Notably, thanks to the Poisson structure, this k² term is generated by the RG if one tries neglecting it⁴⁶. Finally, the spin noise has correlator

$${\langle {\zeta }_{\alpha \beta }\left({{{\bf{k}}}},\,t\right){\zeta }_{\gamma \nu }\left({{{{\bf{k}}}}}^{{\prime} },\,{t}^{{\prime} }\right)\rangle ={(2\uppi )}^{d}4{\tilde{{{\varLambda }}}}_{\alpha \beta \gamma \nu }{\delta }^{(d)}\left({{{\bf{k}}}}+{{{{\bf{k}}}}}^{{\prime} }\right)\delta \left(t-{t}^{{\prime} }\right)}$$

where ${\tilde{{{\varLambda }}}}$ has the same structure as Λ, although out of equilibrium we may have, ${\lambda }^{\perp ,\,\parallel }\ne {\tilde{\lambda }}^{\perp ,\,\parallel }$. In addition to the terms discussed here, new interactions compatible with symmetries are generated by the RG transformation. To be self-consistent and closed, the RG calculation must take into account these new terms (Supplementary Section IC4).

Within the RG analysis it is possible to define a set of effective parameters and couplings that are independent of the field dimensions and upon which physical predictions uniquely depend (Supplementary Section IC5). Because the most important factors are activity and inertia, we focus here on their effective coupling constants:

$${{c}_{v}={v}_{0}\ {\gamma }_{v}\ \frac{{\tilde{{{\varGamma }}}}^{1/2}}{{{{\varGamma }}}^{3/2}}\ \,\,,\,\,f={g}^{2}\ \frac{{\tilde{\lambda }}^{\parallel }}{{{\lambda }^{\parallel }}^{2}{{\varGamma }}}}.$$

(7)

c_v is the effective coupling regulating activity, which vanishes for v₀ → 0, while f quantifies the effective reversible coupling giving rise to inertial dynamics, and we will therefore refer to it as the inertial coupling constant. The scaling dimension of all effective couplings is proportional to 4 − d, indicating that the upper critical dimension is d_c = 4 and that an expansion in powers of ϵ = 4 − d is appropriate. A momentum shell RG calculation at one loop^11,12 produces 65 Feynman diagrams (full details are provided in Supplementary Sections I and II). A rich fixed-point structure emerges, as shown in Fig. 1a.

Fig. 1: RG flow.

a, Flow in the conservative case. The novel fixed point (red circle), with non-zero off-equilibrium activity and non-zero inertial coupling, is the only stable one, with a dynamical critical exponent z = 1.35 in d = 3. The equilibrium non-inertial fixed point (black square), z = 2.0, corresponds to standard ferromagnets (model A of ref. ³⁵). The equilibrium inertial fixed point (blue diamond), z = 1.5, corresponds to superfluids and quantum antiferromagnets (models E and G of ref. ³⁵). Finally, the active non-inertial fixed point (green triangle), z = 1.73, corresponds to active matter without reversible coupling between velocity and spin¹⁸. This last fixed point is not connected to the active inertial one onto this plane. b, Flow with spin dissipation. Spin dissipation, η, is a relevant parameter that brings the flow out of the conservative plane. Because η grows up to infinity with the RG iterations, it is convenient to use the reduced dissipation ${\hat{\eta }}={\eta /(1+\eta )}$ to represent the flow. If we perturb the active inertial fixed point, z = 1.35, with some dissipation, the RG flow leaves the ${\hat{\eta }=0}$ plane, until it eventually reaches the active overdamped fixed point for ${\hat{\eta }=1}$ (green pyramid), where z = 1.73. When ${\hat{\eta }\ne 0}$ it is better to represent the flow through the reduced inertial coupling, ${\hat{f}=(1-\hat{\eta })f}$, instead of f, so that in the overdamped limit, ${\hat{\eta }=1}$, we have one less parameter, as the inertial coupling drops out of the calculation. The overdamped fixed point, z = 1.73, is best seen as belonging to the overdamped ${\hat{\eta }=1}$ line, rather than to the conservative but non-inertial line, η = 0, f = 0. Even though the value of ${\hat{f}}$ is the same on the two lines, only the first corresponds to the correct overdamped limit. All flow lines are actual numerical solutions of the RG equations.

The simplest fixed point corresponds to zero activity and zero inertial coupling, ${c}_{v}^{* }={f}^{* }={0}$ (black square, Fig. 1a). This equilibrium non-inertial fixed point describes non-active systems, such as classical ferromagnets, where the polarization is not coupled to the spin; here z = 2 at one loop (model A of ref. ³⁵). Incompressibility is merely a solenoidal constraint on ψ, leading to the universality class of dipolar ferromagnets⁴⁸. If we perturb this fixed point by adding an inertial coupling we reach the equilibrium inertial fixed point (blue diamond), which still has zero activity, but non-zero inertial coupling, ${{c}_{v}^{* }={0},\,{f}^{* }\ne 0}$. This fixed point describes equilibrium superfluids and antiferromagnets (models E/F and G of ref. ³⁵) and it has z = d/2, hence z = 1.5 in d = 3. Here, too, incompressibility is a solenoidal constraint on ψ, which changes the static universality class, but not the dynamical one⁴⁷. This fixed point is unstable against activity, which leads the RG flow towards a novel active inertial fixed point (red circle), where both ${{c}_{v}^{* }\ne 0}$ and f^* ≠ 0. The combined effect of activity and inertia significantly lowers the dynamical critical exponent; in d = 3 we find z = 1.34(8). This fixed point is stable against perturbations of all the parameters considered so far. As we shall discuss more thoroughly later on, we believe this to be the fixed point describing natural swarms.

Finally, there is a fourth fixed point (green triangle), which has non-zero activity, ${{c}_{v}^{* }\ne 0}$, but zero inertial coupling, f* = 0, corresponding to z = 1.73 in d = 3. Here, the inertial reversible terms are absent from the dynamics, so the polarization is decoupled from the spin¹⁸. This active non-inertial fixed point is stable against activity fluctuations, but as soon as we perturb it with an inertial coupling, f ≠ 0, the RG flow diverges (shaded area). There is a sound reason for this: the correct way to attain non-inertial dynamics is not to kill the reversible coupling between coordinate and momentum, but to introduce dissipation and let it take over in the overdamped limit. This is the consistency check our calculation must pass next.

So far, our theory has conserved the total spin, thanks to the Poisson structure generating the reversible terms in the dynamics and to the fact that the irreversible kinetic coefficient Λ is zero at k = 0. Although the Poisson structure has no reasons to change, Λ could—within natural swarms we cannot exclude that some spin dissipation exists, not as a result of a violation of the rotational symmetry, but because individual midges might exchange spin with the environment in a way that is unaccounted for in the equations of motion. Spin dissipation is produced by a k-independent term η in the kinetic coefficient, Λ ~ (λ^∥ + λ^⊥)k² + η, and similarly in ${\tilde{{{\varLambda }}}}$. The RG calculation (reported in Supplementary Section IE) shows that the scaling dimension of η is always positive, so if we perturb the active inertial fixed point with η ≠ 0, the RG flow moves out of the conservative plane and eventually reaches a fixed point at η = ∞, where polarization decouples from the spin and z = 1.73 (Fig. 1b). This is the correct way to obtain the overdamped limit in which inertia becomes irrelevant, so we call this the active overdamped fixed point (a hopefully clarifying map of the theory is depicted in Fig. 2). Yet, if the overdamped fixed point is the only asymptotically stable one, why should we be interested in the inertial fixed point?

Fig. 2: Map of incompressible active theory.

The coarse-grained dynamical equations may either have or not have reversible terms giving rise to inertial coupling between the polarization ψ (that is, the generalized coordinate) and the spin s (that is, the generalized momentum). In the first case (g ≠ 0) we have an inertial theory, with a Poisson structure expressing the fact that s is the generator of the rotational symmetry, thus leading to conservation of the total spin. In the second case (g = 0) we recover the non-inertial theory of ref. ¹⁸, where polarization is decoupled from the spin and the symmetry does not entail any Poisson structure (the equation for s becomes irrelevant). In this case z = 1.73. On the other hand, in the inertial theory, the irreversible kinetic coefficient of the spin may be either conservative or non-conservative. In the conservative case there is no spin dissipation (η = 0), which produces the inertial-conservative fixed point with z = 1.35. In the non-conservative case, the kinetic coefficient contains a dissipative term (η ≠ 0), although the impact of dissipation depends on how strong that is compared to the system size. In the underdamped regime, ηL^z ≪ 1, collective fluctuations are still ruled by the inertial-conservative fixed point, so z = 1.35. This is the regime of natural swarms. Conversely, in the overdamped regime, ηL^z ≫ 1, the Poisson structure is washed out, the spin drops out of the calculation, and collective fluctuations are ruled by the fully non-conservative fixed point, hence giving z = 1.73.

The answer is that finite-size systems can be ruled by a partially stable RG fixed point if the physical parameters are close enough to it. Consider a ferromagnet slightly above its critical temperature T_c. The stable fixed point is T* = ∞, and yet, if the temperature is close enough to T_c, the critical fixed point governs the physics as long as the system’s size is smaller than the correlation length, L ≪ ξ. This is a general mechanism: when the physical couplings are close to a mixed-stability fixed point, the RG flow remains for many iterations in the vicinity of it, and because iterating RG corresponds to observing at larger and larger scales, the flow of a finite-size system may never get out of that basin of attraction. This balance is always regulated by a crossover length scale, ${{{\mathcal{R}}}}$, which is in general a more complicated quantity than ξ, but the upshot is the same: as long as ${L\ll {{{{\mathcal{R}}}}}^{\kappa }}$ (where κ is the crossover exponent), the metastable fixed point rules the system^42,46. How is this relevant for natural swarms? The underdamped shape of the dynamic correlation functions in natural swarms (Fig. 3a) is solid experimental evidence that spin dissipation is weak. On the other hand, rewiring of the interaction network in swarms occurs over the same timescale as velocity relaxation (Fig. 3b); that is, activity is strong. Hence, the RG flow starts close to the conservative plane, η = 0, but far from the equilibrium plane, c_v = 0. As a result, RG rapidly leads the system in the vicinity of the active inertial fixed point, z = 1.35, lingering there for many iterations, before flowing to the overdamped fixed point (Fig. 3c,d). We find ${{{\mathcal{R}}}}={\sqrt{{\lambda }^{\parallel }/\eta }}$ and κ = 2/z (Supplementary Section IE2), so that for ${L\ll {({\lambda }^{\parallel }/\eta )}^{1/z}}$ a finite-size system is ruled by the active inertial fixed point. Given that λ^∥ is finite along the flow, we conclude that as long as ηL^z ≪ 1, the underdamped inertial scenario must hold. Because experimental relaxation is underdamped (Fig. 3a), we conclude that the dynamical critical exponent in natural swarms is z = 1.35.

Fig. 3: RG crossover.

a, Weak dissipation. Given the velocity correlation function, C(t), and its relaxation time, τ, we can define the shape function as ${h(t/\tau )}={-\log C(t/\tau )/(t/\tau )}.$ In the limit t/τ → 0, h(t/τ) → 1 for overdamped exponential dynamics, and h(t/τ) → 0 for inertial underdamped dynamics⁹. Experiments on natural swarms (orange line) and numerical simulations of near-critical ISM (purple line) both show underdamped inertial relaxation. The Vicsek model, on the other hand, belongs to the overdamped class (green line; data from ref. ⁹). b, Strong activity. Shown are velocity (full line) and network (dashed line) dynamical correlation functions for natural swarms (orange) and near-critical ISM simulations (purple). The network correlation function measures the fraction of particles remaining within the n_c nearest neighbours after time t (Supplementary Section III), so it quantifies how quickly the interaction network reshuffles with time. In both natural swarms and ISM, the network decorrelates on the same timescale as the velocity, so they are strongly active systems. Here, n_c = 18, which is the mean number of interacting neighbours in simulations. In Supplementary Section III and Supplementary Fig. 15 we show that in natural swarms the two timescales are the same over all spatial scales. c, Crossover of the flow. A close up of the RG flow around the active inertial fixed point shows that when the flow starts at weak dissipation and strong activity, it first rapidly approaches the active inertial fixed point, staying in its neighbourhood for many RG iterations, and then it crosses over to the overdamped regime (red circles). d, Crossover of the critical exponent. Shown is the RG evolution of the dynamical critical exponent z and of the reduced dissipation, ${\hat{\eta }=\eta /(1+\eta )}$, along the crossover flow line depicted in c. The RG crossover from underdamped to overdamped fixed points corresponds to an actual crossover in real space, such that z = 1.35 for ${L\ll {{{{\mathcal{R}}}}}^{2/z}}$ and z = 1.73 for ${L\gg {{{{\mathcal{R}}}}}^{2/z}}$.

Critical slowing down in natural swarms was first experimentally observed in ref. ⁹, but the spatio-temporal span of the events in that study was somewhat too limited to obtain an accurate determination of z, as the largest swarm had N = 278 individuals. Here we added eight new swarming events to the experimental dataset, notably including a swarm of 780 insects. The relaxation time τ versus correlation length ξ is reported in Fig. 4a. In ref. ⁹ the exponent was determined through least-squares (LS) linear regression of ${\log \tau}$ versus ${\log \xi}$. However, LS works under the hypothesis that the independent variable is perfectly determined and that all experimental uncertainty is in the dependent variable. When this hypothesis is violated, LS systematically underestimate the slope⁴⁹. In our experiments, errors certainly impact both τ and ξ, so LS is not appropriate, and this is why z was unfortunately underestimated in ref. ⁹. Reduced major axis (RMA) regression⁴⁹, on the other hand, treats fluctuations over x and y on the same statistical footing (Methods and Supplementary Section IV). When applied to our dataset, RMA gives z_exp = 1.37 ± 0.11 (Fig. 4). The substantial error bar should make us cautious about the agreement between experiments and theory, also considering the rather uncontrolled approximations our calculation made, most notably incompressibility and the first-order perturbative expansion in powers of ϵ, with ϵ = 1. For this reason, we make a final sanity check of our RG calculation through numerical simulations.

Fig. 4: Experimental and numerical results.

a, Experiments. The logarithm of relaxation time τ is plotted versus the logarithm of the correlation length ξ in natural swarms (logarithms are in base 10). The critical exponent z_exp is the slope of the linear fit. Because experimental uncertainty affects both τ and ξ, standard LS regression (which assumes no uncertainty on the abscissa) systematically underestimates the exponent. RMA regression treats uncertainty on the two variables in a symmetric way by minimizing the sum of the areas of the triangles formed by each point and the fitted line. The inset shows, for a sample point, the construction of the triangle whose area is given by |ΔxΔy|/2 (see Methods). RMA regression gives z_exp = 1.37. The physical ranges for τ and ξ are 80 ms < τ < 610 ms and 50 mm < ξ < 250 mm (Table I of Supplementary Section VI). b, Numerical simulations. Plot of logτ versus ${\log \xi}$ in the ISM. Numerical errors are so small that LS and RMA give the same result, z_sim = 1.35, and the LS error is ±0.04. c, Experimental resampling. To estimate the error bar on the experimental exponent we use a resampling method. We randomly draw 10⁷ subsets with half the number of points and in each subset we determine z using RMA. We report here three such random subsets (orange, selected point; grey, unselected point). Rare experimental fluctuations under resampling can produce an unphysical value of z smaller than 1; this, however, happens in only 0.002% of the data subsets. d, Final comparison. Probability distribution of the experimental critical exponent (orange) from the resampling method of c. The standard deviation of this distribution gives the error on the experimental exponent, ±0.11. The vertical band (purple) indicates the position and error of the numerical critical exponent, and coloured symbols indicate the various RG fixed-point values of z.

The field theory we have studied is the coarse-grained expression of the Inertial Spin Model (ISM³⁷), in which the particles’ velocities are rotated by the spins, and the spins are acted upon by the social alignment forces:

$$\begin{array}{lll}\frac{{\rm{d}}{{{{\bf{v}}}}}_{i}}{{\rm{d}}t}&=&\frac{1}{\chi }{{{{\bf{s}}}}}_{i}\times {{{{\bf{v}}}}}_{i},\\ \frac{{\rm{d}}{{{{\bf{s}}}}}_{i}}{{\rm{d}}t}&=&{{{{\bf{v}}}}}_{i}\times \frac{J}{{n}_{i}}\mathop{\sum}\limits_{j}{n}_{ij}(t)\ {{{{\bf{v}}}}}_{j}-\frac{\eta }{\chi }{{{\bf{{s}}}_{i}}}+{{{{\bf{v}}}}}_{i}\times {{{{\mathbf{\upzeta}}}}}_{i},\\ \frac{{\rm{d}}{{{{\bf{r}}}}}_{i}}{{\rm{d}}t}&=&{{{{\bf{v}}}}}_{i},\end{array}$$

(8)

with noise correlator ${\langle {{{{\mathbf{\upzeta }}}}}_{i}(t)\cdot {{{{\mathbf{\upzeta }}}}}_{j}({t}^{{\prime} })\rangle =2dT\ \eta \ {\delta }_{ij}\delta (t-{t}^{{\prime} })}$. χ is the generalized turning inertia, J is the alignment strength, η is the microscopic spin dissipation (with a small abuse of notation we called it as its mesoscopic counterpart), and T is the noise amplitude (or temperature). The adjacency matrix n_ij(t) is defined by a metric interaction radius, r_c. We want to compare the numerical results with the incompressible RG calculation. Hence, even though we do not impose incompressibility in the simulation, we employ a normalized alignment strength, J/n_i = J/∑_kn_ik, a prescription known to make alignment-based models less prone to phase separation⁵⁰; moreover, we monitor each simulation to be sure that phase separation does not occur. We run three-dimensional simulations in the near-ordering scale-free regime, where ξ ~ L (Supplementary Section VB). In the overdamped limit, η → ∞, and the ISM converges to the non-inertial Vicsek model³⁷, exactly as our dynamical field theory converges to the non-inertial theory of ref. ¹⁸. However, our aim is to check that in the underdamped regime the dynamics of a finite-size ISM simulation is ruled by the inertial fixed point; hence, the dissipation η has been chosen small enough to yield inertial relaxation, as in natural swarms (Fig. 3a). On the other hand, the speed v₀ = |v_i|, has been chosen large enough to be in the active regime, namely to have a network relaxation time of the same order as the velocity relaxation time (Fig. 3b). Full details of the simulation are reported in Methods and in Supplementary Section V. A plot of relaxation time versus correlation length is shown in Fig. 4b. The numerical errors are quite small, so LS and RMA give the same value of z, and we can therefore calculate the error simply through LS. The result is z_sim = 1.35 ± 0.04, in remarkable agreement with the RG theoretical prediction. This consistency also validates the idea that the incompressible theory can indeed be used to describe finite-size compressible systems, as long as density fluctuations are not strong.

A final comment is in order. Of the two keystones of the RG, rescaling and coarse-graining, only the latter produces anomalous critical exponents, giving rise to non-trivial collective behaviours. The technical fingerprint of coarse-graining is the presence of Feynman diagrams: this is the case of the present calculation, which therefore probes the core element of RG. The consistency between theory, simulations and experiments attained here strongly supports the idea that the RG—and its most fruitful consequence, universality—may have an incisive impact also in biology.

Methods

Experiments

Data were collected in the field by acquiring video sequences using a multi-camera system of three synchronized cameras (IDT-M5) shooting at 170 frames per second. Two cameras (the stereometric pair) were at a distance between 3 m and 6 m depending on the swarm and on the environmental constraints. A third camera, placed at a distance of 25 cm from the first camera, was used to solve tracking ambiguities. We used Schneider Xenoplan 50-mm f = 2.0 lenses. Typical exposure parameters were aperture f = 5.6 and an exposure time of 3 ms. Recorded events had a time duration between 0.88 and 15.8 s (Table I of Supplementary Section VI). More details are availabile in ref. ⁴. To reconstruct the 3D positions and velocities of individual midges, we used the tracking method described in ref. ⁵¹. Our tracking method is accurate even on large moving groups and produces very low time fragmentation and very few identity switches, therefore allowing for accurate measurements of time-dependent correlations.

Fit of the dynamic critical exponent

Dynamic scaling states that the relaxation time τ_k at wavelength k and correlation length ξ are linked by the relation τ_k = ξ^zΩ(kξ), where Ω is a scaling function. To infer the value of z from experimental data, we measured the relaxation time τ_k of the mode at wavelength k = ξ⁻¹ in different swarming events. Experimental evaluation of τ_k and ξ is discussed in Supplementary Section IV, and follows ref. ⁹. Dynamic scaling in this case reduces to τ ~ ξ^z (where ${\tau \equiv {\tau }_{k = {\xi }^{-1}}}$); hence, logτ = zlogξ + c. In ref. ⁹, z was fitted through a standard LS regression, which gave z_exp = 1.12 ± 0.16 on the dataset of ref. ⁹, and z_exp = 1.16 ± 0.12 on the current larger dataset. The problem with LS, though, is that it assumes that experimental uncertainty is only present in the dependent variable y, which is not true for our experimental data, as both τ and ξ are subject to experimental uncertainty. When using it on a dataset where the error also affects x, LS systematically underestimate the slope⁴⁹. Therefore, LS is not a good method in our case. RMA regression, on the other hand, is a method that works under the hypothesis that both x and y are affected by uncertainties^52,53. RMA fits a set of Gaussian variables x_i and y_i with homogeneous variance ${\sigma }_{x}^{2}$ and ${\sigma }_{y}^{2}$ to a regression line, y = f(x) = αx + β, and it determines α and β through minimization of the sum of the areas of the triangles formed between each point and the regression line with sides parallel to the axis (Fig. 4a, inset). For each point, the area of this triangle is given by, ${\left\vert {{\Delta }}{x}_{i}{{\Delta }}{y}_{i}\right\vert /2}$, where

$${{\Delta }}{x}_{i}={x}_{i}-{f}^{-1}({y}_{i})={\frac{\alpha {x}_{i}+\beta -{y}_{i}}{\alpha }}$$

(9)

$${{{\Delta }}{y}_{i}={y}_{i}-f({x}_{i})={y}_{i}-\alpha {x}_{i}-\beta }.$$

(10)

The function to be minimized is therefore

$${{{\Sigma }}\left(\alpha ,\,\beta \right)=\frac{1}{N}\mathop{\sum }\limits_{i=1}^{N}\frac{{\left({y}_{i}-\alpha {x}_{i}-\beta \right)}^{2}}{\left\vert \alpha \right\vert }}.$$

(11)

The minimization equations, ∂_αΣ = ∂_βΣ = 0, give

$${\alpha }^{2}=\frac{{\mathbb{E}}\left[{y}^{2}\right]-{\mathbb{E}}{\left[y\right]}^{2}}{{\mathbb{E}}\left[{x}^{2}\right]-{\mathbb{E}}{\left[x\right]}^{2}}\,\,,\,\,\beta ={\mathbb{E}}\left[y\right]-\alpha {\mathbb{E}}\left[x\right],$$

(12)

where ${{\mathbb{E}}\left[g(x,\,y)\right]=\frac{1}{N}\mathop{\sum }\nolimits_{i = 1}^{N}g({x}_{i},\,{y}_{i})}$. The sign of α is the same as the sign of the correlation between x and y. A further benefit of RMA compared to other methods, such as LS or effective variance (EV) (both discussed in Supplementary Section IVB), is that the fit is invariant under an interchange of variables, y(x) versus x(y). Moreover, RMA is also invariant under any scale change of the variables, so is not sensitive to the values of σ_x and σ_y, at variance with other methods. RMA is the only method, among those in which the fitted coefficient can be expressed in terms of elementary regression coefficients, that obeys both properties above⁵³. In Supplementary Section IVB we also describe the EV regression method, which requires the experimental errors δτ and δξ as an input; we use EV with two different estimates of the (most problematic) experimental error δτ, and obtain results compatible with those of RMA (z_exp = 1.32 ± 0.18 and z_exp = 1.34 ± 0.18). Given the substantial difficulties in assigning a univocal experimental error on τ to each swarm (Supplementary Section IVA3), we prefer to quote the RMA result z_exp = 1.37 ± 0.11—which is error-neutral—as our most confident determination of the exponent.

Numerical simulations

Equations (8) of the microscopic ISM are numerically implemented using the RATTLE algorithm to enforce the constraint |v_i(t)| = v₀. Simulations are performed in d = 3 in cubic boxes with periodic boundary conditions; the average density is fixed at ρ₀ = 1 and the sizes explored are N = 256, 384, 512, 729, 1,024 and 2,048, with N the number of particles. The effective inertia is χ = 1 and the alignment strength is J = 18. The metric interaction range is r_c = 1.6, corresponding (at this density) to an average number n_c ~ 18 of interacting neighbours. The microscopic spin dissipation is η = 1 and the speed is v₀ = 2; this choice of η and v₀ ensures that the dynamics is clearly underdamped and active (see main text). The temperature (or noise strength) T serves as control parameter of the order–disorder transition; we explored the interval T ∈ [1: 8]. The time step of integration was chosen as dt = 10⁻⁴. For each size N we identify a finite-size ‘critical’ point, T_c(N). At every value of T we run five independent samples, initialized with polarized configurations, of length 6 × 10⁵ to 8 × 10⁵ steps, increasing with system size, and we compute the susceptibility χ and take as T_c(N) the point where this quantity reaches a maximum. We measure the correlation length with the inverse of the wavenumber k = 1/ξ, where the static correlation function at this temperature peaks, and we calculate the relaxation time τ of the velocity’s C(k, t), in complete analogy with the analysis of experimental data. This procedure ensures that the systems are always in a near-critical scaling regime, with the correlation length scaling linearly with the linear system’s size L (Supplementary Section VB).