Collective response to local perturbations: how to evade threats without losing coherence

The most famous among SPP models of collective motion is the VM ([12]), which we now present in its continuous-time version. Consider a system made of N point-like active particles, each labeled by an index i, which move in a d-dimensional space following the equations of motion

$$\begin{align} \frac{\textrm{d}{\textbf{r}_i}}{\textrm{d}t} &= {\textbf{v}}_i \, , \phantom{000000000000} \end{align} \tag{ 1 }$$

$$\begin{align} \eta \frac{\textrm{d}{\textbf{v}}_i}{\textrm{d}t} &= - \frac{\partial{\mathcal{U}}}{\partial{\textbf{v}}_i}+ \lambda_i {\textbf{v}}_i + {\boldsymbol{\xi}}_i \, .\, \end{align} \tag{ 2 }$$

Each particle position $\textbf{{r}}_i\in \mathbb{R}^d$ evolves deterministically according to its corresponding velocity vector $\textbf{{v}}_i \in \mathbb{R}^{d_\textrm v}$ of fixed modulus $\lvert \textbf{{v}}_i \rvert = v_0$, while the latter undergoes a Langevin overdamped dynamics with a 'social' alignment force given by

$$\begin{align} \mathcal{U}(\{{\textbf{v}}_i\}) & = -\frac{J}{2v_0^2}\sum_{ij} n_{ij} {\textbf{v}}_i \boldsymbol{\cdot} {\textbf{v}}_j , \end{align} \tag{ 3 }$$

$$\begin{align} - \frac{\partial{\mathcal{U}}}{\partial{\textbf{v}}_i}& = \frac{J}{v_0^2}\sum_{j} n_{ij} {\textbf{v}}_j . \phantom{0000} \end{align} \tag{ 4 }$$

We recognize in equation (3) the Hamiltonian of the Heisenberg model, with the coefficient J > 0 setting the strength of the alignment interactions, and where the connectivity matrix n_ij specifies the interacting neighbors. However, in contrast to the standard Heisenberg model, here the connectivity matrix itself evolves in time through the particle positions, i.e. $n_{ij} = n_{ij}(\left\lbrace \textbf{{r}}_i(t) \right\rbrace_{i = 1}^N )$. Among the possible evolution rules, we will adopt in the following the paradigm of metric interactions: each particle interacts with those lying within a fixed interaction radius r_c , i.e. $n_{ij}(t) = \Theta(r_c- \lvert \textbf{{r}}_i(t)-\textbf{{r}}_j(t) \rvert )$.

Featuring in equation (2) are a viscosity coefficient η and a white Gaussian noise ${\boldsymbol{\xi}}_i$, which are linked by the Einstein-like relation

$$\begin{align} \langle{{\boldsymbol{\xi}}_i(t) \cdot {\boldsymbol{\xi}}_j (t^{^{\prime}}) }\rangle = 2\textrm d\eta T \delta_{ij} \delta(t-t^{^{\prime}}) \, , \end{align} \tag{ 5 }$$

where the temperature T sets the strength of the fluctuations. Finally, the parameter λ_i is a Lagrange multiplier which can be used to enforce the fixed speed condition $\textrm d\lvert \textbf{{v}}_i \rvert^2/\textrm dt = 0$ [18].

Due to the time dependence in $n_{ij}(t)$, it is known that the VM admits, even in $d\leqslant 2$, a T-driven transition from an ordered phase (flock), where the individual velocity vectors align and add up to a non-zero total velocity

$$\begin{align} {\textbf{V}} \equiv \frac{1}{N} \sum_{i=1}^{N} {\textbf{v}}_i \, , \end{align} \tag{ 6 }$$

to a disordered phase (swarm), where the mean group velocity $\textbf{{V}}$ is null. A convenient order parameter to describe the transition is the so-called polarization vector

$$\begin{align} {\boldsymbol{\Psi}} \equiv \frac{1}{N} \sum_{i=1}^{N} \frac{{\textbf{v}}_i}{|{\textbf{v}}_i|} = \frac{{\textbf{V}}}{v_0}\, , \end{align} \tag{ 7 }$$

and its corresponding (scalar) polarization $\Psi = |{\boldsymbol{\Psi}}| \in [0,1]$.

The VM successfully describes the large-scale behavior of several living and non-living active systems [13, 16]. However, the VM is not appropriate to explain collective turns in flocks of birds. Experiments indeed show that, in turning flocks, the directional information travels obeying a linear dispersion law, with a propagation speed only depending on the degree of order in the system [9]. The Vicsek dynamics does not reproduce this behavior. An intuition of why this occurs can be grasped by looking at equation (2): its structure is that of an overdamped Langevin equation for the velocity, and this kind of equations usually lead to a diffusive dispersion law [19]. The simplest way to obtain linearly dispersive solutions is to reinstate a second order derivative in equation (2), and this is precisely what the ISM does [10]. As in standard second order equations, it is convenient to write the model as a system of first order equations, by introducing appropriate conjugate variables. The ISM then reads

$$\begin{align} \frac{\textrm{d}{{\textbf{r}}_i}}{\textrm{d}t} &= {\textbf{v}}_i \, , \phantom{00000000000000000000000000} \end{align} \tag{ 8 }$$

$$\begin{align} \frac{\textrm{d}{{\textbf{v}}_i}}{\textrm{d}t} &= \frac{1}{\chi} {\textbf{s}}_i \times {\textbf{v}}_i \, , \phantom{000000000000000000000} \end{align} \tag{ 9 }$$

$$\begin{align} \frac{\textrm{d}{{\textbf{s}}_i}}{\textrm{d}t} &= \frac{{\textbf{v}}_i}{v_0} \times \left( - \frac{1}{v_0}\frac{\partial{\mathcal{U}}}{\partial{\textbf{v}}_i} -\frac{\eta}{v_0} \frac{\textrm{d}{\textbf{v}}_i}{\textrm{d}t} + {\boldsymbol{\xi}}_i \right) \, , \end{align} \tag{ 10 }$$

with the noise correlations given in equation (5). In contrast to the VM, in equation (10) the force term does not act directly on the velocity, but rather on its derivative, which is expressed in terms of the new variable $\textbf{{s}}_i$. The vectorial products enforce the constraint on the individual speeds. The particles can thus only change their directions: $\textbf{{s}}_i$ therefore plays the role of an internal angular momentum, regulating the rotations of the individual velocity vectors, and this is why it is called a 'spin'. The new parameter χ plays the role of a rotational inertia.

Let us now specialize the ISM to the planar two-dimensional case, where velocities lie on a plane (i.e. $d_\textrm v = 2$). This case is simpler to handle algebraically, and it gives a direct intuition of the physical and biological meaning of the spin and of inertial dynamics. Besides, it is also the relevant one for collective turns in flocks (as we will discuss later on), and for a variety of other biological groups. Generalization to the three-dimensional case is straightforward.

In the planar case, we can write

$$\begin{align} {\textbf{v}}_i(t) \equiv v_0 \left( \cos\varphi_i(t)\, , \, \sin\varphi_i(t) \right) \, , \end{align} \tag{ 11 }$$

where the phase ϕ_i specifies the angle of the velocity vector with respect to a reference direction. In terms of the phases, equations (9) and (10) assume a much simpler form, namely

$$\begin{align} \dot{\varphi}_i &= \frac{s_i}{\chi} \, , \phantom{00000000000000} \end{align} \tag{ 12 }$$

$$\begin{align} \dot{s}_i &= - \frac{\partial{\mathcal{U}}}{\partial{\varphi}_i} -\dfrac{\eta}{\chi} s_i + \xi_i \, .\,\, \end{align} \tag{ 13 }$$

In these equations the spin vector reduces to one single component, whose value s_i identifies the angular velocity of each particle. In the presence of a force, the spin acquires a non-zero value and, as a consequence, the particle turns. It can be shown that the instantaneous radius of curvature of the trajectory is proportional to the inverse value of the spin [10].

Equations (12) and (13) can be rewritten in terms of the phase only, giving

$$\begin{align} \chi \ddot{\varphi}_i+\eta \dot{\varphi}_i = - \frac{\partial{\mathcal{U}}}{\partial{\varphi}_i}+ \xi_i \, . \end{align} \tag{ 14 }$$

In the overdamped limit $\chi/\eta\to 0$ this second order equation reduces to a first-order one, which coincides with the planar Vicsek continuous model derived from equation (2). The ISM is therefore an inertial generalization of the VM, which reduces to the latter in the dissipative limit. The opposite limit in which η → 0 renders a deterministic equation with a Hamiltonian reversible structure (see equations (12) and (13)). In this case, since s_i represents the generator of the rotational symmetry of the velocity, its total value $S = (1/N)\sum_i s_i$ is a conserved quantity.

Let us now discuss how the ISM can appropriately describe the dispersion law observed in flocks of birds [9]. To do so, we recall first what is known from experiments [9, 20]:

• Bird flocks are highly ordered systems, with measured polarization values $\Psi \sim 0.9$;
• Turns are planar, i.e. each trajectory lies approximately on a 2d plane;
• Mutual distances do not change during turns, and individuals sweep equal radius paths;
• Turns have a localized spatial origin, and the signal propagates linearly in time through the group with speed $c_\textrm s \gg v_0$ and with negligible attenuation.

Given these premises, we can now focus on the planar version of the model, as in equation (14). This expression becomes even simpler if the system is in the deeply ordered phase, as flocks are. In this case, if we choose as a reference direction the one of the global velocity $\textbf{{V}}$, the individual phases are very small, i.e. $\varphi_i\ll 1$. One can then expand the potential $\mathcal{U}$ in equation (3) up to quadratic order—which is called the spin-wave approximation (SWA)—to find

$$\begin{align} \mathcal{U}=- \sum_{ij} n_{ij} \cos(\varphi_j -\varphi_i) \sim \sum_{ij} \Lambda_{ij} \varphi_j \varphi_i \ , \end{align} \tag{ 15 }$$

where we have introduced the discrete Laplacian Λ with

$$\begin{align} \Lambda_{ij} = -n_{ij} + \delta_{ij} \sum_k n_{ik} \, . \end{align} \tag{ 16 }$$

Since individuals do not change their mutual positions during turns, we can assume that the interaction network n_ij remains approximately the same. If the phases vary slowly from one individual to the next, we can further treat phases as a continuous field $\varphi_i\to \varphi(\textbf{{x}})$. The discrete Laplacian then becomes a continuous one, $\Lambda_{ij}\to -a^2 {\boldsymbol{\nabla}}^2$ (a being the mean nearest-neighbor distance), and equation (14) becomes

$$\begin{align} \chi \ddot{\varphi} + \eta \dot{\varphi} -a^2 J n_\textrm c \nabla^2 \varphi=\xi \, , \end{align} \tag{ 17 }$$

with $n_\textrm c = \sum_j n_{ij}$ (assumed to be constant on a regular network). While computations can be easily performed also in terms of the discrete variables, the continuous form is more convenient for visualizing the dispersion relation. Indeed, let us now set ξ = 0 into equation (17) (or equivalently focus on $\langle \varphi \rangle$), and Fourier transform in both space and time to get

$$\begin{align} \chi \omega^2 -i\eta \omega - a^2 J n_\textrm c k^2=0 \, . \end{align} \tag{ 18 }$$

From this relation we can immediately see that in the Vicsek case (χ = 0) the dispersion law becomes purely diffusive, i.e. $\omega \sim i k^2$, leading to strong attenuation and a non-linear relation between space and time. Conversely, when χ is different from zero a more complex dispersion law is obtained,

$$\begin{align} \omega_{\pm} = i\gamma \pm c_\textrm s \sqrt{k^2 - k_0^2} \, , \end{align} \tag{ 19 }$$

where we introduced the propagation speed $c_\textrm s$ and the effective dissipation γ as

$$\begin{align} c_\textrm s = a \sqrt{Jn_\textrm c/\chi} \, , \end{align} \tag{ 20 }$$

$$\begin{align} \;\;\,\,\gamma =\frac{\eta}{2\chi} \ , \phantom{000000} \end{align} \tag{ 21 }$$

and where $k_0\equiv \gamma/c_\textrm s$.

From equation (19) we can clearly evince that the two conditions observed in flocks of birds, i.e. linear dispersion law and no attenuation, are reproduced when $\gamma \ll 1$ and $k_0\ll k$. In this deeply underdamped regime, the ISM predicts a linear dispersion law with propagation speed given by equation (20). This is a highly non-trivial relation linking together the way directional information propagates through the system, and the degree of order (set by J). Remarkably, this prediction is very nicely verified in experimental data [9], thus supporting the idea that the ISM in the strongly underdamped regime provides a good description for collective turning in flocks.

More generally, according to the ratio of inertia and dissipation, and depending on the size of the system, the ISM interpolates between the dissipative Vicsek limit, and the efficient limit of flocks. For $k \ll k_0$ all the modes are overdamped, while for $k \gg k_0$ the system sustains linear propagation with speed $c_\textrm s$ and constant attenuation with damping time $\tau = \gamma^{-1}$. Since k₀ is the inverse of a length scale, it sets a limit on the size L of a flock through which a signal can linearly propagate: indeed, imposing $k_{\textrm {min}} = L^{-1} \gg k_0$ implies $L \ll c_\textrm s/\gamma$. Another equivalent way to understand this constraint comes by thinking in terms of timescales. Indeed, the time needed for a signal of speed $c_\textrm s$ to cross the entire flock of size L is

$$\begin{align} \tau_\textrm s = \frac{L}{c_\textrm s }=\frac{L}{a}\sqrt{\frac{\chi}{J n_\textrm c}} \, , \end{align} \tag{ 22 }$$

but the same signal gets attenuated over a time scale

$$\begin{align} \tau=\gamma^{-1} \ . \end{align} \tag{ 23 }$$

We can then identify two dynamical propagation regimes:

• Underdamped regime (inertial propagation): $\tau_\textrm s \lt \tau$. In this case directional changes travel linearly through the whole system before attenuation can deteriorate the signal. Flocks of birds observed in experiments [9] belong to the deep extreme of this regime, where attenuation is negligible for all the individuals. In terms of the parameters of the model, the condition defining this regime is $L^2 \eta^2/(4 a^2 J n_\textrm c \chi) \lt 1$.
• Overdamped regime (dissipative propagation): $\tau_\textrm s\gt\tau$. Propagation of information is inefficient because the signal gets dissipated before reaching the other end of the group. In this case, the signal travels unaffected up to certain scales $k^{-1}\lt L$, but it is damped on larger scales. This might result in a strong deformation of the group, or even in its fragmentation. It is likely that fish schools investigated in [4] belong to this regime.

The dispersion law and the way information is propagated through the system determine the dynamical behavior of the scalar polarization and of the correlation functions [18]. In this work, we are interested in understanding how collective turns might arise in finite groups, and we shall find that the dispersion law affects in a crucial way the occurrence of these events. To do so, we will study the response of finite systems to a local perturbation under different boundary and motility conditions, and we will explore the dynamical regimes numerically by tuning the various parameters of the model.

First of all, we consider a system in its polarized phase. Indeed, our aim is to investigate how a state of collective motion can change under external perturbations. We then perform the following approximations:

• Spin-wave approximation: the system is highly polarized, and the individual velocity vectors $\textbf{{v}}_i$ deviate weakly from the orientation of the global polarization $\textbf{{V}}$. If the latter is chosen to be initially aligned to the x-axis when the external field is applied, then $\varphi_i\ll 1$ in the initial phases of the collective turn.
• Fixed network: we assume that the adjacency matrix n_ij appearing in the equations of motion is no longer a function of time. This assumption is reasonable when the timescale over which the collective turn develops is much shorter than the typical reshuffling time of the interaction network n_ij . This condition of local equilibrium [23] is precisely what happens in flocks of starlings. A deviation from this condition is however expected to arise for sufficiently large values of v₀. This analysis will be the subject of section 5.

Under these approximations, the equations of motion for s_i and ϕ_i take the form

$$\begin{align} \dot{\varphi}_i &= \frac{s_i}{\chi} \, , \end{align} \tag{ 28 }$$

$$\begin{align} \dot{s}_i & \simeq - J \sum_{j} \Lambda_{ij} \varphi_j -\varphi_i h_i \cos \alpha_i + h_i\sin \alpha_i \!-\!\dfrac{\eta}{\chi} s_i + \xi_i \, , \end{align} \tag{ 29 }$$

where α_i denotes the direction of $\textbf{{h}}_i$, while $h_i = h_i(t)$ is its magnitude. We choose as the relevant observable the average polarization angle

$$\begin{align} \Phi(t) = \frac{1}{N} \sum_{i} \langle \varphi_{i} (t) \rangle \, . \end{align} \tag{ 30 }$$

This $\Phi(t)$ coincides, for small angles ϕ_i , with the angle formed by the polarization measured at time t (after the perturbation) with the initial polarization (i.e. the one maintained in the absence of perturbations). It therefore characterizes the way the system changes in time its collective flight direction.

The derivation of $\Phi(t)$ within the microscopic approach presented above is carried out in detail in appendix A.1 for a generic adjacency network n_ij ; an alternative derivation, for a regular lattice and in terms of coarse-grained fields $\varphi(\textbf{{x}},t)$, is presented in appendix A.2, leading to the same result. Here we simply present the final expression, which will be compared to numerical simulations in section 4.

As we observed in section 2.1, collective turns in real bird flocks generally present a well-localized origin, as if in response to some punctual external stimulus, and only then they propagate throughout the group. This indicates that the perturbing field should be chosen as local, i.e. $\textbf{{h}}_i\propto \delta_{ip}$, where p labels a specific particle chosen inside the flock. Applying a local field $\textbf{{h}}_p$ on the particle labeled by p induces on the mean polarization defined in equation (30) a response (see appendix A)

$$\begin{align} \Phi(t) = \frac{1}{N\eta} \int_{0}^{t} {\textrm{d}}{t^{^{\prime}}} h_p(t^{^{\prime}}) \sin \alpha_p(t^{^{\prime}}) \left[ 1 - e^{-2 \gamma (t-t^{^{\prime}})} \right] \, . \end{align} \tag{ 31 }$$

This equation has been derived under the assumption that the interaction network n_ij is symmetric (but generalizations are possible, see equation (A26)). Moreover, it should be considered valid within linear response for small $\textbf{{h}}$—however, its validity extends to higher orders by choosing $\alpha_p\equiv \pi/2$, as we show in appendix A. This is the choice we will adopt in our numerical simulations in the following sections.

In the following, we will choose a step-like perturbation in the form $h_p (t) = A_0 \Theta(t)$ applied at a fixed angle $\alpha_p(t)\equiv \alpha_p$, for which equation (31) reduces to

$$\begin{align} \Phi(t) = \frac{A_0 \sin \alpha_p}{\eta N} \left[ t - \dfrac{1 - e^{-2 \gamma t} }{2 \gamma } \right] \, . \end{align} \tag{ 32 }$$

Equation (32) tells us that, when the size $N$ of the system becomes very large (i.e. in the thermodynamic limit), the system will never change the direction of the global order parameter on observational time scales. This is what happens in a ferromagnet at equilibrium: the order is stable at low temperature, and a local field applied on one site cannot change the total magnetization. However, biological groups have much smaller sizes than a physical condensed matter system. In this case, times $t\sim \mathcal{O}(N)$ can in fact be small enough to be reachable in experiments, and the system can thus change its global direction on observational time scales.

We also note that, according to equation (32), the group changes its direction independently of whether the propagation of information is inertial or dissipative. This last feature is however specific to the fixed-network condition: in this case the perturbed individual does not move, and it is always interacting with the rest of group. Even though information arrives damped as it travels through the system, it is injected continuously at the perturbed site and sooner or later everyone will turn. As we will discuss later on, this is not what happens for finite groups of moving individuals, when—if information does not propagate unaltered and quickly enough—the perturbed individual will leave the group before everyone else can follow (see section 6).

The shape of the turn is described by equation (32): for times $t\ll \gamma^{-1}$ the polarization angle $\Phi (t)$ grows with a quadratic dependence, while at larger times a linear behavior is predicted. However, we remind that equation (32) is valid only for small phases ϕ_i , and it describes the polarization angle at most up to the time at which the turn is complete (i.e. the angles clearly do not increase indefinitely, but a saturation occurs, which is not captured by our equation). Therefore, it is convenient to introduce the timescale of the collective turn $\tau_\text{turn}$, defined as the time at which $\Phi(t)$ becomes of ${\mathcal{O}}(1)$, and which represents the observational time window of interest for our analysis. According to the value of $\tau_\text{turn}$ as compared to the typical timescales of the model, one might actually observe either a linear or a quadratic behavior before the saturation:

• Linear behavior. If $\tau_\text{turn}\gg \gamma^{-1}$, then we can disregard the exponential decay in equation (32) to find
  $$\begin{align} \Phi(t) &\simeq \frac{A_0 \sin \alpha_p}{\eta N} t \, , \end{align} \tag{ 33 }$$
  $$\begin{align} \tau_\text{turn} &\simeq \frac{\eta N}{A_0 \sin{\alpha_p}} \ . \end{align} \tag{ 34 }$$
  This condition is met for $\eta^2 N /(2\chi A_0 \sin{\alpha_p})\gg 1$, and it is typically what happens when dissipation is large with respect to inertia.
• Quadratic behavior. The opposite regime is found when, on the contrary, $\tau_\text{turn}\ll \gamma^{-1}$. In this case, we can expand the exponential in equation (32) for small arguments and we find
  $$\begin{align} \Phi(t) &\simeq \frac{A_0 \sin \alpha_p}{2 \chi N} t^2 \ , \phantom{000}\, \end{align} \tag{ 35 }$$
  $$\begin{align} \tau_\text{turn} &\simeq \left ( \frac{\chi N}{A_0 \sin\alpha_p}\right )^{1/2} \, . \end{align} \tag{ 36 }$$
  This conditions is met for $\eta^2 N /(2\chi A_0 \sin{\alpha_p})\ll 1$, and it is typically what happens for very underdamped systems, in which inertia dominates over dissipation. In this regime, however, when the angle α_p is different from $\pi/2$, the expression for $\Phi(t)$ is slightly more complex (see appendix A.1).

Finally, let us note that the conditions for the turn to appear linear or quadratic in the phase growth $\Phi(t)$ (which are stated above, and which involve the field amplitude A₀) are not exactly the same as the conditions for underdamped/overdamped propagation reported in section 2.4 (which involve instead the interaction strength J). One can thus envision realizations of the system in which propagation is underdamped, but the turn still appears linear. As we mentioned before, the linear/quadratic behaviors represent the two extreme limits of equation (32), which describes the polarization response in general; we distinguished these two behaviors mainly for practical reasons, since they are useful to check numerically the theoretical predictions, as we will do in the next sections.

The behavior of the polarization angle described by equation (32) is the immediate consequence of a very general feature of the potential $\mathcal{U}(\{{\textbf{v}}_i\})$, namely the rotational invariance of the velocity-velocity interactions. In the low-noise phase, the system polarizes in a well-defined direction, thus breaking the symmetry. As in standard O(n) models, the system remains however subject to soft Goldstone modes, i.e. low-energy excitations, in the subspace perpendicular to the polarization vector. In other terms, due to the presence of noise, the vectorial polarization freely fluctuates like a random walk within this zero-mode subspace (see appendix B). When an external field is applied, it provides a bias to this random walk, which results in the (almost) linear time dependence exhibited by equation (32).

However, the spontaneous fluctuations provide another reference time scale. In fact, even in the absence of a field, a finite system subject to fluctuations will eventually depart from its original direction if we wait for long enough. In appendix B we compute explicitly this wandering time $\tau_\textrm w$, which grows with the system size N and is regulated by the amplitude of the noise T. For times $t\gt\tau_\textrm w$, fluctuations have generally already changed (randomly) the direction of the polarization, and it becomes meaningless to speak of perturbation-response events. A condition to be satisfied in our analysis is therefore $\tau_\text{turn}\ll\tau_\textrm w$. Since the polarization response depends on the field amplitude, for this to happen we roughly need $A_0\gg T$ (see equations (34) and (36)), a requirement which is easily achieved for very ordered systems. More accurate estimates of this condition are provided in appendix B.

If the system lives in two dimensions (i.e. $d_{\textrm{v}} = 2\; \& \ d = 2$), when the network is kept fixed the model is analogous to the XY model. This is not the case for flocks (which perform planar turns but are really three-dimensional systems), but we will adopt this simplification in the numerical simulations discussed in the next sections. It is known that the XY model does not exhibit long-range order in the thermodynamic limit. However, for finite-sized systems the order parameter (the magnetization, equivalent to the polarization defined here) remains finite in the low-temperature phase, and it slowly decreases for $N\to\infty$. In particular, for sizes comparable to those considered here, the system is fully ordered at low temperature (see appendix C).

Here we describe our first simulations of collective turns. Hereafter, the parameters of the model are chosen so that the system is highly polarized (${\Psi} \gt 0.98$, as in real starling flocks [9]). This is useful for two reasons: first, it allows to obtain polarization curves that are not excessively noisy; second, by setting the flock in its deeply ordered state we are enforcing the conditions for the SWA, under which our predictions were derived (see section 3.2).

Let us introduce the following numerical protocol:

• (i)  
  The system of size N is initialized on a square lattice of side $L = \sqrt{N}$ in a disordered configuration. We wait for the system to thermalize and then we measure the direction of the order parameter, which is chosen as the reference value, i.e. $\Phi(t = 0)\equiv 0$;
• (ii)  
  We select a value for the field amplitude A₀;
• (iii)  
  We apply the external field on a single individual p inside the system, and we keep the field switched on for a time T_turn;
• (iv)  
  We switch off the field and let the system thermalize for a time T_term;
• (v)  
  To obtain more robust observations, we repeat steps (iii) and (iv) for n_turns times, and we average over the various realizations of the turning event.

We can first use our numerical simulation to check the linear and the quadratic behavior predicted for the polarization curves $\Phi(t)$, see equations (33) and (35). We keep the parameters J, T and η fixed, and we use the inertia χ as the tuning parameter to explore the different regimes of τ_turn vs. γ⁻¹. Let us start with the case $\tau_\textrm{turn}\gt \gamma^{-1}$, where our analytical derivation predicts a linear behavior of the polarization angle $\Phi(t)$—see equation (33). A typical mean polarization curve obtained in this regime is shown in figure 1(a) for two values of χ, together with a linear fit of the initial part of the curve (see inset). Recall that, due to the SWA, the agreement with the analytical prediction is expected to break down as the polarization angle saturates to its final value $\alpha_p = \pi/2$. To test the agreement of the prefactor in equation (33) with our numerical observations, we have repeated steps (i) to (iv) for several values of A₀ in the range $[5,30]$, and for each value of A₀ we have fitted the initial part of the mean curve $\Phi(t)$ with a linear function. Figure 1(b) shows our results for the same values of χ as in figure 1(a), which confirm the prediction of a linear dependence of the slope on A₀.

Zoom In Zoom Out Reset image size

Keeping the values of the parameters J, η and T fixed, while pushing the inertia χ towards higher values, brings the system more into the underdamped region where $\tau_\textrm{turn}\lt\gamma^{-1}$. Here the predicted time dependence for $\Phi(t)$ is quadratic, as described by equation (35). We note that an initial quadratic behavior is actually expected to occur in any regime, since there is always a time window (however short) in which the condition $\gamma t \ll 1$ is met—see section 3.2. However, as discussed in the previous section, the quadratic behavior holds for the whole duration of the turn only when $\tau_\text{turn}\ll\gamma^{-1}$. This occurs if $\eta^2 N /(2\chi A_0 \sin{\alpha_p})\ll 1$, which translates to $\chi \gg 5$ for the parameters used in figure 1(c), where we show two polarization curves $\Phi(t)$ for two distinct values of χ. We then repeat the same analysis as for the linear case: for each value of A₀ we fit the initial part of the polarization curve $\Phi(t)$ with a parabola (see inset in figure 1(c)), and check that the estimated prefactor grows linearly with A₀ (see equation (36)). This is precisely what happens, as displayed in figure 1(d).

Finally, a curious feature of the deeply underdamped regime is the occurrence of oscillations in the saturation region of the polarization angle (see the blue curve in figure 1(c)). This is intuitive once we recall that the parameter χ plays the role of a moment of inertia in equation (14). A quantitative estimate (see appendix E) predicts oscillations with frequency $\Omega = \sqrt{A_0/(\chi N)}$ and damping factor γ, in line with numerical observations.

Our numerical protocol additionally allows us to investigate and check the dispersion law predicted in section 2.4. The behavior of the polarization angle $\Phi (t)$ discussed above describes how the group as a whole rearranges its direction to align with the external field. It does not describe, however, how the perturbation—which is applied locally to a specific individual—is propagated from individual to individual through the system. To explore this issue, we need to consider the phases $\varphi_i(t)$ of the single individuals, and monitor how they change in time. In figures 2(a) and (b) we thus show these individual curves as a function of time, for two different sets of the model's parameters. The perturbed individual (blue curve) is the first to change its direction and to complete the turn. Other individuals display a very similar turning profile, but shifted in time with a certain delay $\Delta t_i$, which is larger the farther the individual is from the location of the perturbation (colored curves from left to right). This behavior indicates that the directional information spreads progressively through the system. Starting from these curves, we can quantify the dispersion law as follows. Let us identify the turning time of a given individual as the time when its direction reaches a threshold angle ϕ = 0.05 (dashed horizontal line in figures 2(a) and (b)). The set of individuals which start turning at a given time $\Delta t_i$ determine the turning front at that time, and their distance x_i from the perturbed individual identifies the location of the front, which travels obeying the dispersion law $x_i(\Delta t_i)$. We plot the latter in figure 2(c) for several simulations performed with different sets of parameters. The figure confirms the picture outlined in section 2.4, where two propagation regimes were identified according to whether the time for information to travel through the entire system, $\tau_\textrm s = L/c_\textrm s$, is smaller or larger than the time for the attenuation to damp the signal, $\tau = \gamma^{-1}$. When $\tau_\textrm s\ll \tau$ a linear dispersion law ($x\sim t$) is observed; conversely, when $\tau_\textrm s\gg \tau$ we find a dissipative behavior ($x\sim \sqrt{t}$). With the parameters used in figure 2, one has $\tau/\tau_\textrm s = 1$ for χ = 0.5.

Zoom In Zoom Out Reset image size

The numerical implementation of the off-lattice model is similar to that of the on-lattice case, but this time the individual positions $\textbf{{r}}_i$ evolve ballistically—according to the corresponding velocity vectors $\textbf{{v}}_i$—in a squared box with periodic boundary conditions. The interaction network $n_{ij}(t)$ must now be updated regularly: this is obtained by implementing the cell-list method (see [26–28] and appendix D). We remind that we chose a metric interaction rule, i.e. n_ij is different from zero if individuals i and j have a mutual distance lower than the interaction range $r_\textrm c$.

Studying the phase diagram of the model is made nontrivial by the presence of spatial aggregation effects. Even Vicsek-like models (i.e. first-order models) are known to exhibit strong spatial heterogeneities, for sufficiently low noise, in the case of additive interactions [25]. By this we mean an interaction term (in the evolution equation for the individual velocity vector) of the form

$$\begin{align} \sim \sum_{j \in \mathcal{N}_i} \sin ( \varphi_j - \varphi_i) \, , \end{align} \tag{ 38 }$$

which is not normalized by the total number $\mathcal{N}_i$ of particles interacting with the ith (as it happens instead in non-additive models). As a result, the interaction becomes stronger and stronger as new particles enter the interaction range of the ith, which might cause the formation of clusters. Since the equations of motion of our model feature an interaction term of this form (see section 3.2), we expect a similar mechanism to take place.

Clustering phenomena were indeed observed in our numerical simulations. We considered an initially uniform distribution of individuals; when a density fluctuation brings a few individuals close together, the interaction force increases and eventually dominates over the fluctuations, which are then unable to disrupt the cluster. A detailed analysis of the clustering formation goes beyond the scopes of this work, but one generally finds that these aggregation effects are more prominent at high activity, low temperatures (in the ordered phase of the system), and high values of the inertia (i.e. towards the underdamped regime). On general grounds, we expect the timescale of the clustering process to be mostly influenced by the inertia χ, which slows down the effective motility of the system.

The off-equilibrium phase diagram in figure 3(a) shows the scalar polarization Ψ of the model as a function of the temperature T, for various speeds v₀. In constructing the phase diagram, we selected a range of dynamical parameters $(\eta,\chi)$ for which the most severe aggregation effects are absent. As we noted, the density $\rho = N/L^2$ of the system and the interaction radius $r_\textrm c$ heavily affect clustering phenomena and, more in general, the dynamics of the model. Hereafter we work with ρ = 1 and $r_\textrm c = 1.5$, which ensure (within non-clustered, homogeneous systems) an average number of interacting neighbors $n_\textrm c \sim \pi r_\textrm c^2 \rho \sim 7$. This value reproduces the one measured in real starling flocks, where each bird is found to coordinate with its nearest 7–8 individuals [29].

Zoom In Zoom Out Reset image size

Note that, while numerical simulations of the 3d ISM have been previously reported in the literature [10, 30, 31], the ones presented in this work are the first numerical simulations of the 2d ISM.

After analyzing the unperturbed case, we now address the role of the activity in the presence of a local perturbation. We noted in section 4 that a turn always takes place in a finite-size, on-lattice system in response to a step-like local perturbation. However, once we bring in a nonzero activity, we are no longer guaranteed that a local perturbation will produce a similar behavior on the entire system. Indeed, two timescales are now expected to compete: that of information propagation (related to the speed $c_\textrm s$, see equation (20)), and the reshuffling time of the network $n_{ij}(t)$, which is the time taken by an individual to change its neighborhood. If the latter process is too fast, then a given individual may leave its interaction neighborhood before the signal is able to propagate. We thus generically expect a threshold activity $v_0^\text{lim}$ such that, for $v_0\lt v_0^\text{lim}$, the off-lattice system behaves similarly to the on-lattice case, and our analytical predictions for the time dependence of the polarization angle $\Phi(t)$ still apply.

We can estimate $v_0^{\text{lim}}$ by using a simple heuristic argument. Consider two distinct interacting individuals, i.e. at a relative distance smaller than the interaction range $r_\textrm c$. Let δϕ denote the mean angular deviation of their flight direction with respect to that of the flock $\textbf{{V}}$, as in the inset of figure 3(b). In the fully polarized limit $\Psi = 1$ one has $\delta \varphi = 0$, meaning that all the birds fly straight and never cross; conversely, in the presence of fluctuations δϕ, one can derive under the SWA the relation (see appendix B.2)

$$\begin{align} \delta \varphi \equiv \sqrt{\langle \varphi^2 \rangle} \sim \sqrt{2(1- \Psi)} \, . \end{align} \tag{ 39 }$$

In the worst-case scenario, the velocity vectors of the two individuals point outwards as in the inset of figure 3(b), so that they drift away from each other and stop interacting after a time t_r defined as

$$\begin{align} 2 v_0 \, \delta \varphi\, t_r \sim r_\textrm c \, . \end{align} \tag{ 40 }$$

The timescale $t_r \sim r_\textrm c/(2 v_0 \, \delta \varphi)$ thus estimates the reshuffling time of the connectivity matrix $n_{ij}(t)$. When an external field acts on a single individual, the information must reach its nearest neighbors before reshuffling happens, in order for a collective turn to take place in the way it occurs on a fixed network. Such information propagates across the inter-particle distance a within a time $t_a \sim a/c_\textrm s$, where $c_\textrm s$ is the information propagation speed introduced in equation (20). We conclude that information propagates similarly to the on-lattice case if $t_a \ll t_r$, so that the network can be seen as quasi-fixed; in turn, this implies

$$\begin{align} v_0 \ll \dfrac{r_\textrm c}{2a} \dfrac{c_\textrm s}{\delta \varphi} = \dfrac{r_\textrm c}{2} \sqrt{\dfrac{J n_\textrm c}{2 \chi (1-\Psi) } } \equiv v_0^{\text{lim}} \, . \end{align} \tag{ 41 }$$

Using the same protocol that we applied in the on-lattice case (see section 4), we now analyze the effects of a local perturbation. Again, the parameter χ can be tuned in order to explore the dynamical regimes of the model. The choice of the other parameters is instead dictated by two new requirements. First, we should make sure that the system is not in a clustered state. We then check beforehand that the time scale of clustering is much longer than the turning time $\tau_\text{turn}$ of the flock; this requires us to tune the temperature T in order to remain within a well-polarized state, where the system is spatially uniform. Secondly, numerical limitations (such as time discretization) would prevent us from exploring the region with $v_0 \gt v_0 ^{\text{lim}}$, if the latter were excessively high. All in all, this leads to the choice of parameters of figure 3(b), which corresponds to $\Psi \sim 0.98$ and $v_0 ^{\text{lim}} \sim 10$ according to equation (41).

We focus here as an example on the overdamped region $\tau_\text{turn}\gg \gamma^{-1}$, but similar results were observed also in the underdamped case. In this regime we know that the polarization angle $\Phi(t)$ should behave linearly (see equation (33)): in figure 3(b) we thus fit its initial part with a straight line, and we plot the so-obtained slope against A₀ for several values of v₀. Most of these lines overlap, which is expected since equation (33) is v₀-independent; conversely, the behavior starts to deviate from a straight line as we approach the region with $v_0 \sim v_0 ^{\text{lim}}$. This marks the breakdown of the on-lattice approximation.

We finally note that, in the entirety of this section, the activity v₀ of the system does not seem to affect the ability of the local field to trigger a collective reorientation, at least for speeds not much larger than $v_0 ^{\text{lim}}$ (we explored values up to $v_0 = 30$); indeed, the mean polarization angle always saturates to $\Phi (\infty) = \pi/2$, even in the overdamped regime where propagation of information is attenuated by damping. The reason why this happens is crucially related to the presence of periodic boundary conditions: even if the individual feeling the perturbation moves, it never leaves the system, and the perturbation therefore keeps acting on the latter indefinitely. This is not the case, on the contrary, when OBCs are considered. In the next section we will investigate this new setting, where—as we will discuss—the impact of the activity and the regime of information propagation are crucial in determining the occurrence of a turn.

To implement a turn with OBCs, we first evolve the system off-lattice for a short time interval, until the scalar polarization reaches the stationary value $\Psi \sim 0.99$, and then we apply the field to an individual placed in the middle of the group. We then proceed, as in the previous sections, to measure the evolution of the mean polarization angle. Since (now) the group has finite boundaries, more stringent conditions must be met in order for a global turn to take place. Indeed, the perturbed individual will quickly follow the field direction and deviate from the group's mean velocity: if all the other individuals do not follow (and turn by the same angle) soon enough, a finite difference in orientation with respect to the turn initiator will persist when the latter reaches the boundary. This can result either in the flock remaining compact but turning by a smaller amount, or else in the group's fragmentation. For a full turn to occur we therefore need i) the information about the angular deviation to propagate intact to all individuals (i.e. underdamped inertial propagation), and ii) the turn to be complete before the initiator hits the boundary (i.e. not too large motility v₀).

To verify this picture quantitatively, let us start by considering the first condition. We choose a small value of v₀—for which we know the analytical predictions to hold even off-lattice—and we explore the various dynamical regimes of the model by appropriately varying the other parameters. Since in the off-lattice model large values of χ enhance clustering effects, we keep the inertia fixed and tune instead the dissipation η to span the underdamped/overdamped spectrum.

The resulting curves for the mean polarization angle $\Phi(t)$ are reported in figure 4. This figure shows that only for sufficiently low η (i.e. in the underdamped inertial regime) the flock is able to sustain a collective turn and follow the perturbation. At higher values of η (i.e. in the overdamped dissipative regime), the dissipation is stronger and the perturbed individual leaves the flock before the whole group can turn, so that the global flight direction is only partially affected by the perturbation. By definition, the mean polarization angle $\Phi(t)$ saturates to $\alpha_p = \pi /2$ for $t \to \infty$ only if the entire flock remains coherent and aligns with the field, whereas $\Phi(\infty) \lt \pi /2$ corresponds to a fragmentation of the group. Indeed, if the group splits into two components (one containing the perturbed individual and moving along the field direction, and another one going in another direction), then the resulting $\Phi(\infty)$ will be given by the weighted average over the two flight directions, with the weights represented by the sizes of the two clusters. The fewer individuals follow the perturbed one, the more $\Phi(\infty)$ will differ from $\pi/2$. In the inset of figure 4 we thus plot the limiting polarization angle $\Phi(\infty)$ as a function of $\tau/\tau_\textrm s$, i.e. the ratio between the dissipative and inertial timescales. As discussed in section 2.4, this ratio identifies whether the system is in the overdamped ($\tau/\tau_\textrm s\lt1$) or underdamped ($\tau/\tau_\textrm s\gt1$) regime of information propagation. This plot therefore shows that only when the system enters the underdamped inertial regime the turn becomes fully efficient and coherent (i.e. $\Phi(\infty) = \pi/2$).

Zoom In Zoom Out Reset image size

To further illustrate the mechanism described above, we also plot in figure 5 the single polarization angles $\varphi_i(t)$ of all the individuals in the group, for a single realization of the turning event. At small values of η in the underdamped regime, the flock remains coherent and the collective turn can take place. At intermediate values of η, only some individuals follow the perturbed one, and the group gets disrupted. Finally, at high η, the perturbed individual exits the flock while the flying direction of the rest of the group changes only slightly. Videos showing these turning events in real-time are included in the supplementary material, supporting the scenario we just described.

Zoom In Zoom Out Reset image size

Let us now discuss what is the role of the activity v₀ in determining the occurrence of a turn. To do so, we consider the system in its underdamped regime, where we know from the previous section that—at small values of v₀—the system is able to fully reorient, following the perturbed individual, in the direction of the applied field. We thus choose in the simulation the same parameters as in figure 4, with η = 0.3. Then we start increasing the activity of the individuals, to see whether and how much the turning performance is affected.

As discussed before, the ability of the flock to sustain a collective turn is determined not only by the efficiency of information propagation, but also by the time taken by the perturbed individual to exit the group. We can expect that the higher v₀, the harder will be for the flock to turn, since the perturbed individual will leave the system sooner. A simple argument can be used to predict the final direction of flight $\Phi(\infty)$ of the flock depending on its activity. At the beginning of the turn the flock has a large polarization Ψ, so we can imagine the system as if collectively moving straight; when an external step-like field is applied to a bird located in the middle of the flock, the perturbed individual changes its flight direction by $\pi /2$ within a negligible time span. The time taken by the perturbed bird to exit the flock is thus of the order of

$$\begin{align} t_\text{exit} \sim L / v_0 \, . \end{align} \tag{ 42 }$$

In general, we can assume that the system stops turning at $t\sim t_\text{exit}$, so that $\Phi(\infty) \sim \Phi(t_\text{exit})$ (at least in the absence of oscillations). For the parameters considered in the simulation, the system is in the underdamped propagation regime but has a linear growth (see figure 4). We can adopt the on-lattice prediction in equation (33) to estimate $\Phi(t_\text{exit})$, which gives (since $N = L^2$ and $\alpha_p = \pi/2$)

$$\begin{align} \Phi(t=t_\text{exit}) \sim \dfrac{A_0}{N\eta} \dfrac{L}{v_0} = \dfrac{A_0}{\eta L v_0}\, . \end{align} \tag{ 43 }$$

Indeed, our current choice of parameters gives $v_0 \ll c_\textrm s\simeq 2.12 \ll v_0^{\text{lim}}$, so that the on-lattice estimate in equation (33) is still reliable according to our discussion in section 5.2. Equation (43) suggests that, when v₀ increases, the final mean polarization angle should decrease. However, this estimate (and in particular the dependence on $1/v_0$) only holds when the exit time t_exit falls within the linear growth regime of the polarization angle. This is not the case for sufficiently low activity v₀, because we expect the system to perform a full turn in this limit, and $\Phi(t = t_\text{exit})$ must therefore approach $\pi /2$. This picture is confirmed by our numerical simulations. In figure 6(a) we display the mean polarization angle for different values of the activity v₀, and we show that indeed for sufficiently large values of v₀ the group is unable to complete the turn. In the inset, we plot the endpoint $\Phi(t_\text{exit})$ as a function of $1/v_0$: as expected, the linear dependence predicted by equation (43) gradually breaks down at small v₀, when $\Phi(t_\text{exit})$ approaches $\pi/2$.

Zoom In Zoom Out Reset image size

It is interesting to compare equation (43) with the response of the polarization $\Phi(t)$ to an impulse-like perturbation of finite duration $\mathcal{T}$. This is derived in appendix F for the on-lattice case with a calculation analogous to the one described in section 3.2, yielding

$$\begin{align} \Phi(t\to \infty) = \frac{A_0 \sin\alpha_p }{N \eta} \mathcal{T} \, . \end{align} \tag{ 44 }$$

The latter clearly agrees with the off-lattice prediction in equation (43) upon choosing $\mathcal{T} = t_\text{exit}$, i.e. the time during which the perturbed individual remains inside the group (see equation (42)).

Expression (43) also suggests that the final polarization angle, and therefore the efficiency with which the turn is performed, also depends on the size of the group itself. Ceteris paribus, large groups tend to be more in the overdamped regime, and they would take more time to perform a turn even on a fixed network. To better pinpoint these effects, we focus on the overdamped regime by using the same parameters as in figure 4 with η = 5, and then we vary the size $N \in [100,1600]$. The result is displayed in figure 6(b), where the behavior predicted by equation (43) is very well obeyed, including the linear dependence of the final flight direction on $1/ \sqrt{N}$ (see inset).

We start from the coupled Langevin equations (28) and (29), which we recast for convenience in the matrix form

$$\begin{align} \frac{\textrm{d}}{\textrm{d}t} {\boldsymbol{\Psi}} + \hat{L} {\boldsymbol{\Psi}} = {\textbf{F}}(t) \, , \end{align} \tag{ A1 }$$

upon defining the 2 N-dimensional vectors

$$\begin{align} {\boldsymbol{\Psi}} \equiv \begin{pmatrix} {\boldsymbol{\varphi}} \\ {\textbf{s}} \end{pmatrix} \, , \;\;\; {\textbf{F}} \equiv \begin{pmatrix} {\textbf{f}}_{(\varphi)} \\ {\textbf{f}}_{(s)} \end{pmatrix}= \begin{pmatrix} {\textbf{0}} \\ {\boldsymbol{\xi}}+{\textbf{h}}\sin{\boldsymbol{\alpha}} \end{pmatrix} \, . \end{align} \tag{ A2 }$$

In particular, the N-dimensional vector ϕ has components $\varphi_i(t)$. We also introduced the block matrix

$$\begin{align} \hat{L} = \begin{pmatrix} \mathbb{O} & -\mathbb{I}/\chi \\ J \tilde{\Lambda} & \eta \, \mathbb{I}/\chi \end{pmatrix} \, , \end{align} \tag{ A3 }$$

where $\tilde{\Lambda}$ is defined in terms of the discrete Laplacian $\Lambda_{ij}$ given in equation (16) as

$$\begin{align} \tilde{\Lambda} = \Lambda + \delta \Lambda \, , \;\;\;\;\; (\delta \Lambda)_{ij} = \delta_{ij} \frac{h_i \cos\alpha_i}{J} \, . \end{align} \tag{ A4 }$$

Equation (A1) admits a formal solution as

$$\begin{align} {\boldsymbol{\Psi}} (t) = e^{-t\hat{L}} {\boldsymbol{\Psi}}(0) + \int_0^{t} {\textrm{d}}{t^{^{\prime}}} \, e^{-(t-t^{^{\prime}})\hat{L} } {\textbf{F}}(t^{^{\prime}}) \, , \end{align} \tag{ A5 }$$

which can be made explicit by diagonalizing $\hat{L}$.

A.1.1. Spectrum of $\hat{L}$

Let us initially set $h_i = 0$ in the expression for $\hat{L}$, which is equivalent to replacing $\tilde{\Lambda}$ with Λ in equation (A3). We will first compute the unperturbed spectrum of $\hat{L}$ and only later add the effect of a small perturbing field h_i . Let us adopt for convenience the bra-ket formalism familiar from Quantum Mechanics [33], which allows us to express ${\boldsymbol{\varphi}} \rightarrow \vert \varphi\rangle$, $\textbf{{s}} \rightarrow \vert s\rangle$, and

$$\begin{align} {\boldsymbol{\Psi}} \rightarrow \vert \Psi\rangle = \vert \varphi\rangle \otimes \vert s\rangle \, , \end{align} \tag{ A6 }$$

where the symbol ⊗ denotes the tensor product operation. The components of $\vert \varphi\rangle$ and $\vert s\rangle$ along the position basis

$$\begin{align} \mathcal{B}_{\textrm {pos}}: \; \vert i\rangle \otimes \vert i\rangle \, , \end{align} \tag{ A7 }$$

where $\vert i\rangle$ is a N-dimensional vector, can be found by taking the scalar products

$$\begin{align} \varphi^i=\langle i|\varphi\rangle\,,\;\;\;s^i=\langle i|s\rangle\, . \end{align} \tag{ A8 }$$

Conversely, let us denote as

$$\begin{align} \mathcal{B}_{a}: \; \vert a\rangle \otimes \vert a\rangle \, , \;\;\;\;\; \Lambda \vert a\rangle = \lambda_a \vert a\rangle \end{align} \tag{ A9 }$$

the basis in which Λ is diagonal. The new components ϕ^a , labelled by the superscript a, are connected to the old ones (ϕⁱ ) via

$$\begin{align} \varphi^a&=\langle a|\varphi\rangle=\langle a\vert\sum_i\vert i\rangle\langle i\vert\varphi\rangle \nonumber\\ &=\sum_i \langle a|i\rangle\varphi^i\equiv \sum_i(U^{-1}){^a{}_i}\,\varphi^i \, , \end{align} \tag{ A10 }$$

where in the first line we inserted a decomposition of the identity $\mathbb{I}$, while in the second line we introduced the N × N matrix U, which enforces the change of basis

$$\begin{align} U \, : \; \mathcal{B}_{a} \rightarrow \mathcal{B}_{\textrm {pos}} \, . \end{align} \tag{ A11 }$$

When written in the basis $\mathcal{B}_{\textrm {pos}}$, the matrix U contains the eigenvectors of Λ as its columns; as a consequence, U diagonalizes Λ as

$$\begin{align} \text{diag}(\lambda_1,\,\lambda_2,\,\dots) = U^{-1} \Lambda U \, . \end{align} \tag{ A12 }$$

In the case of symmetric interactions $n_{ij} = n_{ji}$, the condition $\Lambda_{ij} = \Lambda_{ji}$ makes U unitary ($U^{-1} = U^{\dagger}$); note that the interaction network does not need to be regular (e.g. a square lattice) in order for this to occur.

Working in the basis $\mathcal{B}_{a}$, we now look for the eigenvectors of $\hat{L}$:

$$\begin{align} \hat{L}\vert l\rangle = \omega_l \vert l\rangle\, , \;\;\;\;\; \vert l\rangle = \vert l\rangle_\varphi \otimes \vert l\rangle_s \, . \end{align} \tag{ A13 }$$

Without loss of generality, we choose both $\vert l\rangle_\varphi$ and $\vert l\rangle_s$ to be proportional to $\vert a\rangle$, i.e.

$$\begin{align} \vert l\rangle = \vert a\rangle \otimes c(l)\vert a\rangle \;\;\; \xrightarrow{{\mathcal{B}_a}}\;\;\; \kappa \begin{pmatrix}1 \\ c(l) \end{pmatrix}\, , \end{align} \tag{ A14 }$$

where κ is some normalization constant. The functional form of c(l) has to be found by solving the eigenvalue equation (A13), which we can rephrase as

$$\begin{align} \begin{pmatrix}\omega_l & 1/\chi \\ - J \lambda_a & \omega_l - \eta/\chi \end{pmatrix} \begin{pmatrix}1\\ c(l) \end{pmatrix} = \begin{pmatrix}0\\ 0 \end{pmatrix} \, . \end{align} \tag{ A15 }$$

Solving for c(l) leads to

$$\begin{align} -J \lambda_a - \omega_l^2 \chi + \eta \omega_l = 0 \, , \end{align} \tag{ A16 }$$

so that, calling as usual $\gamma = \eta/(2\chi)$, we obtain

$$\begin{align} \omega_l = \gamma \pm i \omega_a \, , \;\;\;\;\; \omega_a \equiv \sqrt{J\lambda_a/\chi-\gamma^2} \, . \end{align} \tag{ A17 }$$

Using equation (A15) yields $c(l) = -\chi \omega_l$, which finally gives the eigenvectors

$$\begin{align} \vert l\rangle \xrightarrow{{\mathcal{B}_a}} \kappa \begin{pmatrix}1 \\ -\chi \omega_l \end{pmatrix} \, . \end{align} \tag{ A18 }$$

For each fixed value of a (there are N such values), equation (A17) shows as expected a twofold degeneracy for $\omega_l \equiv \omega_\pm (\lambda_a)$. We can thus represent

$$\begin{align} \hat{L} \xrightarrow[\text{fixed } a]{{\mathcal{B}_a}} L \equiv \begin{pmatrix}0 & -1/\chi \\ J \lambda_a & \eta/\chi \end{pmatrix} \, , \end{align} \tag{ A19 }$$

each of which is diagonalized by a matrix M such that

$$\begin{align} \text{diag}(\omega_+, \, \omega_-) = M^{-1}LM \, . \end{align} \tag{ A20 }$$

The latter is non-unitary (L is not symmetrical even when $\Lambda_{ij} = \Lambda_{ji}$), and reads

$$\begin{align} M&=\kappa \begin{pmatrix} 1&1\\ -\chi \omega_+&-\chi \omega_- \end{pmatrix} ,\;\nonumber\\ M^{-1}&=\frac{1}{2 i \chi \kappa\omega_a} \begin{pmatrix} -\chi \omega_-&-1 \\ \chi\omega_+&1\end{pmatrix}\, . \end{align} \tag{ A21 }$$

A.1.2. Temporal evolution

Equation (A5) takes a simple form in the basis $\mathcal{B}_a$. Indeed, by expressing

$$\begin{align} \vert\Psi\rangle =\sum_l \vert l \rangle\langle l \vert\Psi\rangle = \sum_l \Psi^l \vert l\rangle \, ,\;\;\; \vert F\rangle =\sum_l F^l \vert l\rangle \, ,\end{align} \tag{ A22 }$$

and by using the completeness relation

$$\begin{align} e^{-t\hat{L}} = \sum_l e^{-t\omega_l} \vert l \rangle \langle l \vert \, , \end{align} \tag{ A23 }$$

one obtains the set of 2 N decoupled equations

$$\begin{align} \Psi^l (t) = e^{-t\omega_l} \Psi^l(0) + \int_0^{t} {\textrm{d}}{t^{^{\prime}}} \, e^{-(t-t^{^{\prime}})\omega_l } F^l(t^{^{\prime}}) \, . \end{align} \tag{ A24 }$$

In order to translate this result back into the position basis, we need to:

• (i)  
  Use the matrices U and U⁻¹ on each of the two N-dimensional subspaces of the source term $\textbf{{F}}$ in equation (A2), so as to move to the basis in which the Laplacian is diagonal;
• (ii)  
  Apply the matrices M and M⁻¹ on each 2-dimensional subspace at fixed a, moving to the basis in which L is diagonal. This gives the components $F^l(t)$, which contain the source terms;
• (iii)  
  Use equation (A24) to compute the time evolution of each of the components $\Psi^l (t)$;
• (iv)  
  Apply the inverse of the two changes of basis used above (in reverse order), which finally gives $\varphi^i (t)$ and $s^i (t)$.

The above procedure is tedious but straightforward. Here we state the result for the variable $\varphi^i (t)$,

$$\begin{align} \varphi^i (t) & = \sum_{a,j=1}^N U^{\,i}_{a} (U^{-1}){^a{}_j} \int_0^t {\textrm{d}}{t^{^{\prime}}}\, e^{-\gamma (t-t^{^{\prime}})} \nonumber\\ & \quad \times \left\lbrace f^{\;j}_{(\varphi)}(t^{^{\prime}}) \left[ \frac{\gamma}{\omega_a}\sin\omega_a(t-t^{^{\prime}}) + \cos\omega_a(t-t^{^{\prime}}) \right] \right. \nonumber \\ & \qquad \quad + \left. \frac{1}{\chi\omega_a} f^{\;j}_{(s)}(t^{^{\prime}}) \sin\omega_a(t-t^{^{\prime}}) \right\rbrace \, , \end{align} \tag{ A25 }$$

where the second line vanishes in our case since $f^{\;j}_{(\varphi)}\equiv 0$. Equation (A25) is still written in terms of the matrix U, whose columns are the normalized eigenvectors of the discrete Laplacian $\Lambda_{ij}$: as such, it is valid for any time-independent interaction network n_ij , whose eigenvectors could in principle even be obtained numerically. In the following, we will see that we do not need to compute explicitly these eigenvectors in order to derive an expression for $\Phi(t)$ in the case of symmetric interactions.

If we now specialize to the case considered in the main text, where only a field coupled to the phases ϕ_i is considered, the term $f^{\;j}_{(\varphi)}$ is zero (see equations (26) and (27)) and we get

$$\begin{align} \varphi^i (t) & = \sum_{a,j=1}^N U^{\,i}_{a} (U^{-1}){^a{}_j} \int_0^t {\textrm{d}}{t^{^{\prime}}} \, e^{-\gamma (t-t^{^{\prime}})} \nonumber\\ & \quad \times \frac{1}{\chi\omega_a} [h_j(t^{^{\prime}})\sin\alpha_j+\xi_j(t^{^{\prime}})] \sin\omega_a(t-t^{^{\prime}}) \, . \end{align} \tag{ A26 }$$

A.1.3. Corrections to the spectrum of $\hat{L}$

The first-order correction to the unperturbed spectrum when a small external field $\textbf{{h}}$ is switched on can be computed by standard perturbation theory [33],

$$\begin{align} \tilde{\lambda}_a = \lambda_a + \left\langle a |\delta \Lambda |a\right\rangle \, , \;\;\;\;\; \vert \tilde{a}\rangle = \vert a\rangle + \sum_{b\neq a} \frac{\langle{b}\vert{\delta \Lambda}\vert{a}\rangle}{\lambda_a - \lambda_b} \vert b\rangle \, , \end{align} \tag{ A27 }$$

where λ_a are the eigenvalues of the discrete Laplacian Λ and $\vert a\rangle$ are its eigenvectors. This leads to

$$\begin{align} \tilde{\lambda}_a &= \lambda_a + \sum_{i} (U^{-1}){^a{}_i} \frac{h_i\cos\alpha_i}{J} U^{\,i}_{a} \, , \phantom{0}\,\, \end{align} \tag{ A28 }$$

$$\begin{align} \tilde{U}{^i{}_a} &= U^{\,i}_{a} + \sum_{b\neq a} \sum_{j} \frac{(U^\dagger){^b{}_j} (h_j\cos\alpha_j) U{^j{}_a}}{J(\lambda_a - \lambda_b)} U{^{\,i}{}_b} \, ,\nonumber\\ \end{align} \tag{ A29 }$$

which can be used into equation (A25) to replace $U^{\,i}_{a}$ and λ_a , the latter being contained in ω_a via equation (A17). Note that the corrections in equations (A28) and (A29) vanish by choosing $\alpha_i\equiv \pi/2 \; \forall i$.

A.1.4. Average polarization angle

A remarkable property of the discrete Laplacian $\Lambda_{ij}$, which can be inferred from its definition in equation (16), is

$$\begin{align} \sum_j \Lambda_{ij} = \sum_j \left( -n_{ij} + \delta_{ij} \sum_k n_{ik} \right) = 0 \, , \end{align} \tag{ A30 }$$

which means that $\Lambda_{ij}$ always admits a zero mode (i.e. a $\lambda_0 = 0$ eigenvalue) with constant right eigenvector. This implies that the matrix U, which contains the normalized eigenvectors as its columns, must be such that

$$\begin{align} U{^{\,i}{}_0} = 1/\sqrt{N} \;\;\; \forall i \, , \end{align} \tag{ A31 }$$

or equivalently $(U^{\dagger}){^0{}_i} = 1/\sqrt{N}$, since the transposed of the matrix U contains the same (right) eigenvectors as its rows. Note that in general these do not coincide with the left eigenvectors of Λ: this is only the case when $\Lambda_{ij}$ is symmetric, because then the orthogonality condition reads

$$\begin{align} \mathbb{I} = U^{-1}U = U^\dagger U \, . \end{align} \tag{ A32 }$$

In the symmetric case we then have in particular

$$\begin{align} \sum_i U^{\,i}_{a} = \delta_{a0}\sqrt{N} \, , \end{align} \tag{ A33 }$$

and summing over i in equation (A25) in order to obtain $\Phi(t)$, only the term with a = 0 survives. Using equation (A2) (or equation (A26)) with $\langle {\boldsymbol{\xi}} \rangle = 0$, we thus find

$$\begin{align} \Phi (t) = \frac{1}{N}\sum_{j=1}^N \int_0^t {\textrm{d}}{t^{^{\prime}}} \, e^{-\gamma (t-t^{^{\prime}})} \frac{h_j(t^{^{\prime}})\sin\alpha_j}{\chi \omega_0} \sin \omega_0(t-t^{^{\prime}}) \, , \end{align} \tag{ A34 }$$

where we called $\omega_0 = \omega_{a = 0} = i \gamma +{\mathcal{O}}(h)$ (see equation (A17)). Note that equation (A34) is already of ${\mathcal{O}}(h)$, so that in general including the corrections to the spectrum of $\hat{L}$ due to the external field does not modify the leading-order contribution in perturbation theory. The validity of this approximation extends to higher orders by choosing $\alpha_j \equiv \pi/2$, because then we have seen that the first order corrections in equations (A28) and (A29) vanish. This is the choice we adopt in our numerical simulations.

By choosing a local field $h_i\propto \delta_{ip} = A_0\Theta(t)$, we finally recover the expression for $\Phi(t)$ reported in the main text, equation (31). When $\alpha_p\neq \pi/2$, however, one has to be careful in some regions of the parameter space. Indeed, the full explicit expression of ω₀ is (including corrections) $\omega_0 = \sqrt{A_0 \cos\alpha/(N\chi)-\gamma^2}$. Usually the second term dominates, but for strong underdamping (i.e. $\eta^2 N/(2\chi A_0\cos\alpha)\ll 1$) it does not, and $\Phi(t)$ acquires a more complex structure.

The expression of the response $\varphi^i(t)$ we derived in equation (A25) is general and it applies to any interaction network n_ij , even when it is not symmetric. However, we can get to equation (31) of the main text (which refers to the on-lattice case) with much less effort if we work in terms of the coarse-grained fields $\varphi(\textbf{{x}},t)$ and $s(\textbf{{x}},t)$.

The coarse-grained counterpart of our dynamical equations can be easily derived from equations (12) and (13) and read [22]

$$\begin{align} \dot{\varphi} = \frac{\delta{\mathcal{U}}}{\delta{s}} \, , \end{align} \tag{ A35 }$$

$$\begin{align} \quad\;\;\; \dot{s} = - \frac{\delta{\mathcal{U}}}{\delta{\varphi}}- \eta s + \xi \, , \end{align} \tag{ A36 }$$

where the continuous version of the potential is

$$\begin{align} \mathcal{U} = \int \frac{{\textrm{d}}^{d}{x}}{a^d} \bigg[ & \frac{Jn_\textrm c a^2 }{2} \lvert {\boldsymbol{\nabla}}{\varphi} \rvert^2 +\frac{1}{2}\cos\alpha\, h({\textbf{x}},t)\varphi^2 \nonumber \\ &- \sin\alpha\, h({\textbf{x}},t)\varphi \bigg]. \end{align} \tag{ A37 }$$

Note that the component of the external field $\textbf{{h}}$ parallel to the initial polarization $\textbf{{V}}(t = 0)$ plays the role of a mass term for the field ϕ. We stress that $\mathcal{H}$ above has been derived under the SWA, hence it has to be considered as a low-temperature expansion of the ISM. In fact, a field theory can still be constructed if one relaxes this assumption, but it will in general be more complicated and possibly contain new, non-Gaussian terms [31]. From equations (A36) and (A35) we get

$$\begin{align} \dot{\varphi} &= \frac{s}{\chi} \, , \nonumber\\ \dot{s} &= Jn_\textrm ca^2 {\nabla^2}{\varphi} - h \varphi \cos\alpha + h \sin\alpha - \frac{\eta}{\chi}s + \xi \, , \end{align} \tag{ A38 }$$

and deriving again the first equation we obtain

$$\begin{align} \chi \ddot{\varphi} + \eta \dot{\varphi} -Jn_\textrm ca^2 {\nabla^2}{\varphi} + h \varphi \cos\alpha = h \sin\alpha + \xi \, . \end{align} \tag{ A39 }$$

In analogy with the microscopic derivation, we now discard the term proportional to h on the left-hand side, as we expect it to produce higher-order corrections for small h (note that such term actually vanishes for $\alpha = \pi/2$). This leads us to the linear problem

$$\begin{align} \mathcal{D}_{{\textbf{x}},t} \varphi({\textbf{x}},t) = f({\textbf{x}},t) \end{align} \tag{ A40 }$$

upon defining

$$\begin{align} \mathcal{D}_{{\textbf{x}},t} &\equiv \chi \partial_{tt}^2 + \eta \partial_t - Jn_\textrm c a^2 {\nabla^2}_{{\textbf{x}}} \, , \end{align} \tag{ A41 }$$

$$\begin{align} f({\textbf{x}},t) &\equiv \sin\alpha\, h({\textbf{x}},t) + \xi({\textbf{x}},t) \, .\phantom{00}\,\, \end{align} \tag{ A42 }$$

The propagator of the differential operator in equation (A41) can be computed by standard methods [18] to give, in the time-momentum domain,

$$\begin{align} G({\textbf{k}},t) = e^{-\gamma t}\frac{\sin(\omega_{{\textbf{k}}} t)}{\chi\omega_{{\textbf{k}}}} \, , \end{align} \tag{ A43 }$$

where $\omega_{{\textbf{k}}}\equiv c_\textrm s \sqrt{k^2-k_0^2}$ is precisely the one appearing on the right-hand side of equation (19). This yields the average response of the system as

$$\begin{align} \langle \varphi({\textbf{k}},t) \rangle = \int_0^t {\textrm{d}}{t^{^{\prime}}} \, G({\textbf{k}},t-t^{^{\prime}}) \langle {\,f}({\textbf{k}},t^{^{\prime}}) \rangle \, . \end{align} \tag{ A44 }$$

We now turn to the polarization angle in equation (30), whose continuous-space version is

$$\begin{align} \Phi(t) = \frac{1}{\mathcal{V}} \int_\mathcal{V} {\textrm{d}}^{d}{x}\, \langle \varphi({\textbf{x}},t) \rangle = \frac{1}{\mathcal{V}} \langle \varphi({\textbf{k}}={\textbf{0}},t) \rangle \, . \end{align} \tag{ A45 }$$

Here we called $\mathcal{V}$ the volume occupied by the entire flock (i.e. $\mathcal{V}\sim N$ on a square lattice), and we adopted the continuum prescription

$$\begin{align} \sum_{{\textbf{k}}} f({\textbf{k}}) \simeq \frac{\mathcal{V}}{(2\pi)^d} \int {\textrm{d}}^{d}{k} \,f({\textbf{k}}) \, . \end{align} \tag{ A46 }$$

Using equations (A43) and (A44) then simply gives

$$\begin{align} \Phi(t) &= \frac{1}{\mathcal{V}} \int_{-\infty}^t {\textrm{d}}{t^{^{\prime}}} \, G({\textbf{k}}={\textbf{0}},t-t^{^{\prime}}) \sin\alpha H({\textbf{k}}={\textbf{0}},t^{^{\prime}}) \nonumber \\ &= \frac{\sin\alpha}{\mathcal{V} \eta} \int_{-\infty}^t {\textrm{d}}{t^{^{\prime}}}\, H({\textbf{k}}={\textbf{0}},t^{^{\prime}}) \left[ 1- e^{-2\gamma (t-t^{^{\prime}})}\right]\, . \end{align} \tag{ A47 }$$

To make contact with equation (31), it is sufficient to choose a local field and give it a stepfunction-like time dependence,

$$\begin{align} H({\textbf{x}},t) = A_0 \delta({\textbf{x}}-{\textbf{x}}_0) \Theta(t) \;\; \rightarrow \;\; H({\textbf{k}},t) = A_0 e^{i {\textbf{k}} \cdot {\textbf{x}}_0} \Theta(t) \, . \end{align} \tag{ A48 }$$

Let us note, however, the difference with respect to the microscopic approach of appendix A.1: in that case, we have studied the effect of a perturbation applied on one of the microscopic variables $\textbf{{v}}_i$. Here, conversely, the coarse-graining procedure has washed out the identity of the single microscopic degrees of freedom $\textbf{{v}}_i$, and the perturbation in equation (A48) is thus applied locally at position $\textbf{{x}}_0$. Moving off-lattice by introducing a large activity v₀ (see section 5), the former approach is expected to fail because of the time dependence in the interaction matrix $n_{ij}(t)$, while the latter approach fails both because the propagator $G(\textbf{{k}},t)$ in equation (A43) was computed on-lattice, and because it erroneously assumes that the perturbation remains fixed in space at $\textbf{{x}} = \textbf{{x}}_0$.

When the system is in its ordered phase, it is useful to decompose each velocity vector into its components parallel and perpendicular to the global polarization $\textbf{{V}}$,

$$\begin{align} {\textbf{v}}_i = v_i^{\,L} \boldsymbol{\hat{n}} + \boldsymbol{\pi}_i \, , \end{align} \tag{ B1 }$$

where $\boldsymbol{\pi}_i$ is a $(d_\textrm v-1)$-dimensional vector, and by construction one has

$$\begin{align} \sum_i {\boldsymbol{\pi}}_i = 0 \, . \end{align} \tag{ B2 }$$

In the planar model, π_i is simply a scalar,

$$\begin{align} \pi_i = v_0 \sin\varphi_i \, . \end{align} \tag{ B3 }$$

Now let us call $\boldsymbol{\hat{n}}$ the direction assumed by the polarization $\textbf{{V}}$ at time t = 0, but allow $\textbf{{V}}(t)$ to change in time. We can decompose

$$\begin{align} {\textbf{V}} = \frac{1}{N} \sum_i {\textbf{v}}_i = ({\textbf{V}}\cdot \boldsymbol{\hat{n}}) \boldsymbol{\hat{n}} + \boldsymbol{\delta V}^{\perp} \, , \end{align} \tag{ B4 }$$

where in the planar case

$$\begin{align} &({\textbf{V}}\cdot \boldsymbol{\hat{n}})(t) = \frac{1}{N} \sum_i v_i^L(t) \, ,\phantom{000000000000} \end{align} \tag{ B5 }$$

$$\begin{align} \boldsymbol{\delta V}^{\perp}(t) & = \frac{1}{N} \sum_i {\boldsymbol{\pi}}_i = \frac{v_0}{N} \sum_i \sin\varphi_i(t) \boldsymbol{\hat{\pi}} \nonumber\\ & \simeq \frac{v_0 \boldsymbol{\hat{\pi}}}{N} \sum_i \varphi_i(t), \end{align} \tag{ B6 }$$

the latter being valid within the SWA.

The SWA decomposition can be used to derive a relation between the scalar polarization, i.e. $\Psi = \lvert \textbf{{V}} \rvert/v_0$, and the average fluctuation δϕ of a single velocity vector $\textbf{{v}}_i$ around the global flight direction $\textbf{{V}}$. Indeed, using equations (B1) and (B2) we have

$$\begin{align} \Psi & = \dfrac{1}{v_0 N} \Big| \sum_i \mathbf{v_i} \Big| = \dfrac{1}{N} \Big| \sum_i \cos \varphi_i \Big|\nonumber\\ & \simeq 1 - \dfrac{1}{2N} \sum_i \varphi_i ^2 \, , \end{align} \tag{ B7 }$$

where in the last step we used the SWA. Defining the mean fluctuation as

$$\begin{align} \delta \varphi \equiv \sqrt{\langle \varphi^2 \rangle} = \dfrac{1}{N} \sum_i \varphi_i ^2 \, , \end{align} \tag{ B8 }$$

the relation stated in equation (39) follows immediately from equation (B7).

We define the wandering time $\tau_\textrm w$ by the condition

$$\begin{align} \langle \lvert \boldsymbol{\delta \boldsymbol{\Psi}}^{\perp}(t\sim \tau_\textrm w) \rvert^2 \rangle \sim {\mathcal{O}}(1) \, , \end{align} \tag{ B9 }$$

where the brackets denote the average over the noise, and $\textbf{{V}} = v_0{\boldsymbol{\Psi}}$. Using equation (B5), we find under the SWA

$$\begin{align} \langle \lvert \boldsymbol{\delta {\boldsymbol{\Psi}}}^{\perp}(t) \rvert^2 \rangle &= \frac{1}{v_0^2 N^2} \sum_{ik} \langle {\boldsymbol{\pi}}_i \cdot {\boldsymbol{\pi}}_k \rangle \nonumber\\ &\simeq \frac{1}{N^2} \sum_{ik} \langle \varphi^i(t)\varphi^k(t) \rangle \, . \end{align} \tag{ B10 }$$

The evolution of $\varphi^i(t)$ is given by equation (A25) upon setting the external source field $\textbf{{h}} = 0$. This gives, for a symmetric interaction network,

$$\begin{align} \varphi^i (t) & = \sum_{a,j=1}^N U^{\,i}_{a} (U^{\dagger}){^a{}_j} \int_0^t {\textrm{d}}{t^{^{\prime}}}\, \frac{e^{-\gamma (t-t^{^{\prime}})}}{\chi\omega_a} \xi^j(t^{^{\prime}})\nonumber\\ & \quad \times \sin\omega_a(t-t^{^{\prime}}), \end{align} \tag{ B11 }$$

where it is useful to isolate the contribution of the a = 0 mode

$$\begin{align} \varphi^i (t) & = \frac{1}{N} \sum_{j=1}^N \int_0^t {\textrm{d}}{t^{^{\prime}}}\, \frac{\xi^j(t^{^{\prime}})}{\eta} \left[ 1 - e^{-2\gamma (t-t^{^{\prime}})} \right] \; \nonumber\\ & \quad + \; (a\neq 0) \, . \end{align} \tag{ B12 }$$

Indeed, thanks again to the property in equation (A33), the modes a ≠ 0 give no contribution; using the noise variance in equation (5) we thus obtain

$$\begin{align} \langle | \delta {\boldsymbol{\Psi}}^{\perp}(t)|^2 \rangle = \dfrac{2T}{N \eta } \left[t - \dfrac{1-e^{-2 \gamma t} }{\gamma} + \dfrac{1-e^{-4 \gamma t} }{4 \gamma} \right]. \end{align} \tag{ B13 }$$

Specializing this expression for the two regimes $t \gg \gamma^{-1}$ and $t \ll \gamma^{-1}$ and using the condition in equation (B9), we finally get

$$\begin{align} \tau_\textrm w = \left\{\begin{array}{ll} \frac{\eta N}{2 T} \, , & \text{for} \; t \gg \gamma^{-1} \, , \\ \left(\frac{3}{16}\frac{N \eta}{T \gamma^2}\right)^{1/3} \, , & \text{for} \; t \ll \gamma^{-1} \, . \end{array}\right. \end{align} \tag{ B14 }$$

We stress again that equation (B14) can be generically interpreted as the persistence time of the average polarization in the XY model subject to an underdamped Langevin dynamics. Its origin, as discussed above, is related to the presence of the zero modes of the Laplacian, i.e. to the original rotational invariance of the potential function $\mathcal{U}$ of the velocities (see equation (3)). This is a common feature of O(n) models, and it is independent of the kind of dynamics adopted to let the system evolve. Indeed, the wandering of the order parameter and its persistence time have been computed in a variety of works, see e.g. [34, 35] where a microcanonical dynamics for the XY model was considered, or the more recent [21, 36] where the ISM and the VM were analyzed.