Topological flocking models in spatially heterogeneous environments

Flocking models with metric and topological interactions are supposed to exhibit distinct features, as for instance the presence and absence of moving polar bands. On the other hand, quenched disorder (spatial heterogeneities) has been shown to dramatically affect large-scale properties of active systems with metric interactions, while the impact of quenched disorder on active systems with metric-free interactions has remained, until now, unexplored. Here, we show that topological flocking models recover several features of metric ones in homogeneous media, when placed in a heterogeneous environment. In particular, we find that order is long-ranged even in the presence of spatial heterogeneities, and that the heterogeneous environment induces an effective density-order coupling facilitating emergence of traveling bands, which are observed in wide regions of parameter space. We argue that such a coupling results from a fluctuation-induced rewiring of the topological interaction network, strongly enhanced by the presence of spatial heterogeneities. While the presence of spatial heterogeneities is known to affect emergent collective behavior in traditional active systems with metric interactions, its role in active systems with nonmetric interactions is still unclear. Here, the authors study the effect of quenched disorder in active matter with topological interactions, finding that topological interactions preserve long-range order in the presence of spatial heterogeneity, which also enhances the emergence of traveling bands.


F
locking is a fascinating self-organized phenomenon observed in a large number of artificial and biological systems 1,2 , including bacterial swarms [3][4][5][6] , fish schools 7 , and sheep herds 8 , among many other examples. The large-scale properties of these active systems crucially depend on the type of interaction neighborhood of the moving agents. Two fundamentally different types of interaction neighborhood have been explored, the so-called metric 9 and topological ones 10 .
In metric models, the neighborhood of a particle is defined via the Euclidean distance between the focal and the neighboring particles, and the number of neighbors of the focal particle scales with the local particle density. As a result of the competition between velocity alignment among neighbors and noise-induced decoherence, metric flocking models undergo spontaneous symmetry breaking 2 . In ideally homogeneous media, the order that emerges in these non-equilibrium systems in two dimensions is long-ranged (LRO) 11,12 , giant density fluctuations 12,13 are observed, and the phase transition is characterized by the presence of high-order, high-density bands that move across the system [14][15][16][17][18][19][20] . The presence of spatial heterogeneities or quenched disorder (e.g., obstacles, inhomogeneous substrates, etc.), ubiquitous in all experimental and real-world active systems 21 , dramatically affects the large-scale properties of metric flocking models. In scalar active matter, it was found that obstacles can lead to jamming, frozen states, and moving chains [22][23][24][25] . For vectorial active matter in the presence of quenched disorder, it was shown first numerically 26 and later analytically 27 that order becomes quasi-long-ranged (QLRO). It was also found, using a minimal model, there exists an optimal noise that maximizes collective motion 26 , a result later confirmed in more realistic simulations 28 . Furthermore, it was also predicted that spontaneous particle trapping leading to anomalous transport can occur 29 , a prediction in line with recent findings in bacteria 30 . In addition, it was also shown the existence of multiple attractors for flocks flying through the same realization of quenched disorder, meaning that the fate and history of the flock are strongly dependent on the initial condition 31 . Finally, it was found in models 32 and experiments 33 , that above a given density of spatial heterogeneities, polar bands vanish.
While metric flocking models have been successful in reproducing several real active systems, it has been suggested that animals interact with a specific number of neighbors, regardless of local density, and thus independently of the relative Euclidean distance between the individuals 10 . The large-scale properties of topological flocking models are believed to be fundamentally different from the ones of the metric counterparts. In particular, the phase diagram of these systems, so far only studied in homogeneous media, does not seem to possess a coexistence region characterized by the presence of polar, traveling bands [34][35][36] ; Fig. 1a, d. The absence of traveling bands has been attributed to an apparent lack of a density-order coupling. On the other hand, the impact of quenched disorder on the active matter with topological interactions has, so far, not been addressed.
Here, we address this open question in active matter theory by studying how the quenched disorder affects the emergent properties of topological flocking models using k-nearest neighbors (kNN) and Voronoi tessellation 37 . We find that topological models differ fundamentally from their metric counterparts by exhibiting long-range order even in the presence of heterogeneities. Furthermore, we observe that in topological models, spatial heterogeneities counter-intuitively facilitate the emergence of traveling, polar bands (Fig. 1b, e; Supplementary Movies 1 and 2), while such elongated structures are believed not to be present in homogeneous media 34,38 . Finally, we argue that band formation is related to the emergence of an effective coupling between local density and local order (Fig. 1c, f) due to local rewiring of the interaction network, which is strongly enhanced by the presence of spatial heterogeneities.
Our study provides a comprehensive characterization of the large-scale properties of topological flocking models in heterogeneous environments. The results reported here, together with those by Martin et al. 38 , strongly suggest that the established knowledge on topological flocking models needs to be fundamentally revised. Specifically, our analysis extends our understanding of topological interactions in active matter systems by showing that topological flocking models in complex environments behave as metric ones in homogeneous media.

Results
Model. We consider active particles moving at constant speed v 0 in a two-dimensional, heterogeneous environment with periodic boundary conditions. The heterogeneous environment is modeled by a random distribution of "obstacles", which we also will refer to as quenched disorder or spatial heterogeneities. Each active particle interacts with its topological neighborhood (TN), which defines the particle's local environment. We use two definitions of TN: (i) the first k-nearest neighbor (kNN) objects, and (ii) all objects in the first shell by performing a Voronoi tessellation. Note that neighboring objects include other active particles, as well as obstacles. The behavior of particles is different for TN objects corresponding to active particles and obstacles: particles align their velocity to that of neighboring active particles and move away from obstacles. The equations of motion of i-th particle are given by: where dots on the left-hand side denote temporal derivatives, x i is the position of the particle, and θ i encodes the moving direction of the particle given by Vðθ i Þ ¼ ðcosðθ i Þ; sinðθ i ÞÞ. The first term in Eq. (2) describes the alignment of the particle with TN active particles (subset TN a ), while the second term describes repulsion from TN obstacles (subset TN o ). The symbol TN denotes the set of topological neighbors of particle i, including n b,i active particles and n o,i obstacles. The position of TN obstacles is given by y s , and α s,i denotes the angle, in polar coordinates, of the vector x i − y s . Note that "obstacles" are in fact areas that the active particles avoid by turning away from their center (y s ), which can be viewed as a soft-core repulsive interaction. Finally, γ is a constant and ξ i is a delta-correlated, dynamic noise such that 〈ξ i (t)〉 = 0 and hξ i ðtÞξ j ðt 0 Þi ¼ δ ij δðt À t 0 Þ; η is a constant that denotes the strength of the dynamic noise. We studied two options for g i (n o,i ) that lead qualitatively to the same results: (a) g i (n o,i ) = 1 for all values of n o,i , and (b) The latter option of g i (n o,i ) ensures that in the presence of obstacles, the active particle gives priority to obstacles, and move away from them, ignoring other active particles. Since results are easier to interpret with this rule, and are qualitatively the same as those obtained with g i (n o,i ) = 1, we illustrate the system behavior using the obstacle priority rule; results for g i (n o,i ) = 1 can be found in Supplementary Fig. S1. In the following, we fix γ = 1, v 0 = 1, dt = 0.1, and particle density ρ b = N b /L 2 = 1, with N b the number of active particles in the simulation box of linear size L (see "Methods" for further details). Note, that we have studied the dynamics of the above model recently also in the context of collective information processing 37 .
Dynamic noise vs. quenched disorder. The system considered here contains two sources of fluctuation that promote misalignment among the active particles: the dynamic noise and the quenched disorder (i.e., the obstacle field). For vanishing dynamic noise-i.e., in the limit of "cold" active matter-the initial condition and specific distribution of obstacles determine the temporal evolution of the flock, implying that the system is not ergodic 31 . By including a non-vanishing dynamic noise, the system remains strictly speaking non-ergodic; however, time average quantities over long time-intervals can become independent of the initial condition. Furthermore, we can expect that quenched disorder realizations sharing the same statistical properties-e.g., same density of randomly distributed obstacles-lead to similar time average quantities, as occurs for flocking models with metric interactions 26 .
To disentangle the level of fluctuation resulting from dynamic noise and quenched disorder, we compare the polar order parameter-defined as puted over different realizations of statistically identical disorders; Fig. 2 (see Supplementary Fig. S2 for the corresponding plots of kNN, k = 6). Note that the standard error of the mean, r, over disorder realizations (red vertical lines), is either smaller than or of the same order of the variance of the polar order r(t) over time for a single disorder realization (black curves). This strongly suggests that the large-scale properties of the system are highly similar among disorder realizations that share the same statistical properties. Finally, it is worth mentioning that for a given disorder realization in a finite system, though order can emerge in a large number of directions, not all of them exhibit the same probability.
Optimal noise and long-range order. As shown in Fig. 3a-c, the polar order parameter is a monotonically decreasing function of the noise strength η in homogeneous environments with vanishing obstacle density ρ o = 0, whereby ρ o = N o /L 2 and N o is the number of obstacles in the system. One of the most remarkable features of metric flocking models in complex environments, i.e., for ρ o > 0, is the non-monotonic functional form of the curve r vs. η that puts in evidence the presence of an optimal noise that maximizes collective motion 26 . This optimal noise is absent in topological flocking models with kNN interaction: the curve r vs. As observed for active particles with metric interactions in heterogeneous media 26,32 , we also find for Voronoi interaction, that small, yet finite values of dynamical noise facilitate exchange of directional information between clusters. At the optimal noise value, the different clusters merge into a band-like structure and the global orientational order becomes maximal. A further increase of dynamical noise leads then to a monotonous decrease in order. In short, the existence of an optimal noise in topological flocking models seems to be model-dependent. A fundamental difference between topological and metric flocking models in complex environments is observed at the level of the emergent order. Metric flocking models in homogeneous media display LRO, while in heterogeneous media, the order was shown, first numerically 26 and later by an RG argument 27 , to become QLRO: the polar order parameter r Fig. 1 The role of heterogeneity in flocking with topological interactions. Panels a-c corresponds to k-nearest neighbors interaction (kNN, k = 6 where k is the number of neighboring objects), and panels d-f are for Voronoi interaction. Snapshots indicate macroscopic configurations, a, d in obstacle-free environments, where obstacle density ρ o = 0, and b, e in complex environments, with obstacle density ρ o = 0.051. Black and red dots represent particles and obstacles, respectively. The insets in b, e show 1D band profiles, that is 1D particle density ρ L along the moving direction. c, f Local order (r ℓ ) vs local density (ρ ℓ ) measured in small cells of size l = 14 in a simulation box with linear size 140 corresponding to the systems in previous panels (see "Methods" for details). Black scatters are for homogeneous, and red scatters are for heterogeneous environments. Noise intensity is fixed at η = 0.65 here. decays algebraically with system size. On the other hand, topological models in homogeneous media and non-vanishing being the number of active particles, also exhibit LRO 34,35 . By keeping ρ b and ρ o constant, while increasing N b and N o , we provide solid numerical evidence indicating that the polar order parameter r converges towards a constant value in the thermodynamic limit for both kNN and Voronoi neighbors at low and high obstacle densities. Specifically, lim 1=N b !0 r ! r 1 , with r ∞ a non-vanishing constant; Fig. 4.
This result can also be obtained by studying h VðrÞ Vð0Þi, where VðrÞ refers to the local, average velocity of active particles in position r, that as expected for LRO converges to a non-zero value for |r| → ∞, see Supplementary Fig. S4. In addition, we have also confirmed the robustness of the observed LRO with respect to variation in the particle density ρ b by simulating systems with a larger and smaller density, ρ b = 1.5 and ρ b = 0.5 respectively (see Supplementary Fig. S5).
In short, topological flocking models in heterogeneous media exhibit LRO, in contrast to the QLRO reported for the metric counterpart. We note that as discussed further below, we observe the formation of large-scale bands for a wide range of parameters, in particular for the kNN model with k = 6. Thus the corresponding LRO results are obtained in the presence of such emergent spatial structures.
Traveling polar bands. In metric flocking models in homogeneous media, the emergence of polar bands has been explained as the result of a coupling between local polar order r ℓ and local density ρ ℓ 17 (see "Methods" for details regarding calculation of r ℓ and ρ ℓ ). On the other hand, topological flocking models have been introduced as active models that lead to large-scale collective motion independently of the local density of the active particles 10 . In short, it has been assumed that in topological flocking models the above-mentioned order-density coupling is not present. Thus, traveling polar bands are not expected to emerge, as illustrated in  In e, black arrows show the particles' instantaneous moving direction. In f the underlying (undirected) Voronoi interaction network is depicted. A link exists between two particles, if they are neighbors in the Voronoi tessellation, and if at least one of them does not have an obstacle in its neighborhood. The green circles point out clusters that are connected by long links, i.e., long-distance interactions. Error bars represent the standard deviation of the polar order parameter over different realizations of the obstacle field (see "Methods"). Note that error bars are often comparable or smaller than the symbol size.  (Fig. 1c, f). To quantify the emergence of traveling, polar bands we introduced a density modulation parameter β, defined via the amplitude of the largest  Fourier mode with finite wave number q > 0 of the Fouriertransformed coarse-grained density field (see "Methods" for details). Figure 5b-e indicates β at different noise values η and obstacle densities ρ o for k = 1 and k = 6. For k = 1, bands are observed only near transition point (orange and red regions). By increasing k-e.g., to k = 6-and ρ o , bands are observed for all η values such that η < η c,H , where η c,H is the critical value in homogeneous media (Fig. 5). We have confirmed that band structures emerge also for larger k (e.g., k = 12) over a wide range of parameters, in particular, different noise intensities also away from the order-disorder transition (see Supplementary Fig. S6).
A core finding is that for fixed values of η and k, bands become more pronounced as the obstacle density ρ o is increased; Fig. 5 (and Fig. 6 for Voronoi interaction). This means that counterintuitively the spatial heterogeneities promote band formation, while in metric models they hinder the formation of bands 26 . It is worth clarifying that this does not mean that spatial heterogeneities promote polar order, which decreases as ρ o increases. However, the presence of obstacles induces, as explained below, a coupling between local density and local (polar) order that leads to band formation. An important observation is that the speed of bands is independent of ρ o -i.e., of quenched disorder-and set by the amplitude η of dynamic noise (see Supplementary  Fig. S7a−d). As the density ρ o is increased, the number of active particles traveling in the bands diminishes, the disorder gas density increases, and as a result of this, the global polarization of the system decreases. One important lesson to draw is that it is not possible to reduce the impact of spatial heterogeneities to a re-normalized dynamic noise, since this would imply that the band speed depends on ρ o , which, we show, it does not.

Discussion
How can spatial heterogeneities promote band formation? In the following, we argue that (local) rewiring of the underlying dynamical network leads to an effective density-order coupling. Our argument is based on the following observations: (i) local (orientational) order is strongly regulated by the level of (local) rewiring facilitating the fast exchange of orientational information between different sets of particles, (ii) obstacles induce local rewiring, and (iii) rewiring is strongly density-dependent, i.e., at high densities, it will occur highly localized in space. This, together with the two previous points results in an emergent coupling between local (particle) density and local order, a necessary condition for band formation.
The first assertion can be shown using a simple model. Assume a finite system of N s spins that when not connected to each other obey _ θ s ¼ ηξ s ðtÞ. At a rate, ν pick a pair i−j of spins and connect them for a finite time during which _ θ i ¼ γ sinðθ j À θ i Þ þ ηξ i ðtÞ and _ θ j ¼ γ sinðθ i À θ j Þ þ ηξ j ðtÞ; for details see "Methods". In this simple model, order-i.e., j∑ s expðiθ s Þ=N s j -increases with rewiring rate ν (Supplementary Fig. S8a). This non-spatial model serves to prove that local rewiring can promote order.
The next step of the argument is to understand that the spatial motion of agents implies rewiring. This is evident for diffusing spins with metric interactions, where the order is enhanced at larger densities or by using larger diffusion coefficients 39 . Here, both effects result in a faster exchange of interaction partners. In actual flocking models, however, the situation is more complex since particle velocity is coupled to θ i and it is not possible to control the rewiring rate ν-defined as the inverse of the average time an edge survives in the dynamical network-without affecting the dynamics of θ i . However, simulations performed with topological flocking models in small systems-Supplementary Fig. S8b and "Methods"-allow us to show that (local) polar order and the rewiring rate ν increase with the density of active particles ρ b . Furthermore, in the vicinity of obstacles, active particles are forced to modify their trajectories, which affects the distance to neighboring particles, and leads, in consequence, to rewiring. Figure 7a, d confirms that, as expected, ν increases with the density of obstacles ρ o . Here, we note that a finite obstacle density introduces naturally a characteristic maximal metric length scale of rewiring of the order of average obstacle distance 1= ffiffi ffi ρ p o . Finally, to quantify the coupling between local density ρ ℓ and local order parameter r ℓ we calculate the mutual information MI(ρ ℓ , r ℓ ) as a measure of the non-linear correlation between ρ ℓ and r ℓ , i.e., we quantify how much knowledge we gain on r ℓ by knowing ρ ℓ (see Fig. 7b, c, e, f). It is important to note that r ℓ is by Fig. 6 Bands in Voronoi interaction. a Typical macroscopic configurations observed in a system with Voronoi interaction, for different obstacle densities ρ o , and noise intensities η. Black dots are particles and red dots are obstacles. In obstacle-free environments (ρ o = 0), we observe rather homogeneous structures with no density segregation in the form of cluster or band. In heterogeneous environments, we observe cluster-like structures at low noise values, here η = 0.05, these clusters connect together to form bands at noise around η = 0.20 and lead to a maximum in polar order parameter (see Fig. 3c). Bands continue to form at higher noise values. b Density modulation parameter β in different regions of phase space specified with obstacle density ρ o and noise intensity η corresponding to the system in panel (a). It should be noted that since macroscopic structures for Voronoi interaction are different from those with k-nearest neighbors interaction (kNN, where k is the number of neighboring objects), with clusters at low noise values and bands at higher noise values, the same color (i.e., reddish regions) may refer to different types of structures when comparing to kNN interaction, and even for Voronoi when comparing high and low noise values. definition bounded to values smaller or equal to 1. For certain parameter choices, r ℓ saturates to almost 1 for all ρ ℓ values. This occurs for instance at low dynamic noise and low obstacles densities. In these situations, it is evident that r ℓ is independent of ρ ℓ , and thus MI(ρ ℓ , r ℓ ) adopts low values. Note that for r ℓ~1 particles are highly aligned and the relative distance between particles remains relatively constant implying a low rewiring rate. On the other hand, in the disordered state, we observe r ℓ ≈ 0 for all ρ ℓ . Overall, in order for r ℓ to be dependent on ρ ℓ , the system cannot be too disordered but also the level of polar order in the system should not be too high. This suggests, as actually observed, that sharper bands are observed at larger obstacle densities, where the level of global order is lower (see Supplementary Fig. S7e). In particular, for a fixed dynamic noise, MI(ρ ℓ , r ℓ ) is higher, indicating a stronger correlation between r ℓ and ρ ℓ , for larger obstacle densities ρ o , where the rewiring rate ν is also higher; see Fig. 7 and compare band snapshots in Figs. 5 and 6.
An interplay between r ℓ and ρ ℓ mediated by rewiring is arguably also present in spatially homogeneous systems. For fixed dynamic noise, it is expected that ν decreases with k and increases with particle density ρ b . Both trends are straightforward to understand under the assumption that the level of positional fluctuations of the particles is set by dynamic noise. For large values of k, positional fluctuations only lead to the replacement of the farthest neighbors, and this implies that most links (those corresponding to closer neighbors) are long-lived. In consequence, the average time an edge survives increases with k, and its inverse, ν, decreases. On the other hand, at high densities the average inter-individual distance between particles is small, and for the same level of positional fluctuations, a higher rewiring rate is expected. Simulations performed in an homogeneous medium, Fig. 8, are consistent with these arguments. In addition, all this suggests that the coupling between r ℓ and ρ ℓ in homogeneous media should be particularly strong for small k values in the vicinity of the onset of order. Figure 8c shows that in an homogeneous medium for k = 1 traveling bands robustly emerge, whereas they become quickly more diffuse with increasing k and at k = 6 are not observable in our simulations anymore (see bottom panels in Fig. 5a). This finding provides additional support for our arguments. At this point, we also would like to point the reader to a recent publication 38 , which we learnt about at the time of submission, providing an alternative explanation for band formation in homogeneous media of flocking models with topological interaction. We note that the different mechanisms are not mutually exclusive in facilitating band formation in active matter with topological interactions.

Conclusions
We performed a comprehensive study of flocking models with topological interactions in heterogeneous environments. We investigated two different types of topological interactions, kNN and Voronoi, which are the two most studied active topological models in the literature 10,[34][35][36]38 . Similarly to what occurs in an equilibrium system, where only few macroscopic details affect the emergent macroscopic behavior, here we found that the largescale properties of these systems do not depend on the details of the implementation of the model-e.g., on the choice of g i (n o,i )but on the nature of these interactions: i.e., topological (metricfree) interactions of polar symmetry. In that sense, our results are generic and we expect them to apply to other variations of topological interactions, as for example the spatially balanced kNN-model 40 . The two main results of our work on flocking models with topological interactions in heterogeneous environments are: (1) We found that in sharp contrast to metric models -where we observe quasi long-ranged order (QLRO) in heterogeneous environments-for topological models, according to our numerical study and up to the systems sizes investigated, the order is long-ranged (LRO) in the presence of spatial heterogeneities. (2) We showed that spatial heterogeneities promote the emergence of an effective density-order coupling that allows  active particles with topological interactions to form traveling polar bands, which share similar features to the bands observed in metric models. Importantly, bands are observed in parameter ranges in which metric models in homogeneous media do not develop bands. Furthermore, we argued that the counter-intuitive emergence of the density-order coupling for topological interactions is the result of the (local) rewiring of the underlying dynamical networks of active particles induced by the spatial heterogeneities. Finally, we expect that the numerical finding of LRO in heterogeneous media for nonmetric active modelsarguably related to the presence of long-ranged connections between distant clusters-will be supported by a RG calculation, as occurred in the past with the observation QLRO in metric models in the presence of quenched disorder 26,27 .
In summary, our results show that topological flocking models in the presence of spatial heterogeneities-which introduce a characteristic distance (1= ffiffiffiffi ffi ρ o p )-behave as metric ones in homogeneous media, an observation that invites to a reconsideration of "metric-free" interactions in active systems.

Methods
Simulation details. The model was implemented in C++ programming language. Stochastic differential equations were solved using the Euler−Maruyama method with an integration time step of dt = 0.1. Topological interactions including kNN and Voronoi implemented using the CGAL computational geometry algorithm library 41 . For kNN, k-d tree algorithm is used, where in order to account for periodic boundary conditions, the main simulation box has been repeated in different directions. In order to find particles in the first shell of Voronoi neighbors, the dual graph of Voronoi diagram i.e., (periodic implementation of) Delaunay triangulation is used.
All the other calculations and data processing have been done using Python and dependent libraries, in particular Numpy 42 and Scipy 43 . Furthermore, we used NetworkX library 44 for network related calculations.
Local density ρ ℓ and local order r ℓ . Density-order coupling plots of Figs. 1 and 7 have been obtained by superimposing the r ℓ and ρ ℓ of 60 snapshots taken from three time windows of simulations. In order to find these local quantities, the simulation box is divided into small cells of linear size ℓ. Accordingly, local density is defined by n ' ' 2 , where n ℓ is the number of particles in the cell. And, local order is defined by r ' ¼ j 1 n ' ∑ i expðiθ i Þj, where the summation is over the n ℓ particles of the cell. For the simulation box of size 140, we have used a cell size ℓ = 14.
Quantification of bands 1D band profile. In order to obtain 1D band profile, the density field of particles is smoothed using the kernel density estimation algorithm, then integrated over the direction perpendicular to the moving direction. Profiles represented in Fig. 1 are the result of averaging over 200 snapshots taken every 10 time steps.
Band speed. Speed of band is obtained by measuring the displacement of the peak of 1D band profile during a fixed time period.
Bandwidth. Considering 1D band profile, bandwidth is obtained from the difference between two points on the horizontal axis where the height of the profile is equal to a quarter of the maximum value.
Density modulation parameter β. In order to quantify bands occurring in different regions of the parameter space, i.e., different k, ρ o , and η, we cannot rely on bandwidth obtained from 1D band profiles. Since, in addition to single bands, we also observe band-like density modulations or multiple bands, some of which are merging and splitting during the course of simulations. Therefore, obtaining a bandwidth, which is representative of configurations of all the time steps, is in general not possible. In order to address this problem, we use the Fourier transformation of the coarse-grained density field and identify the maximal amplitude of the resulting Fourier spectrum for a finite spatial frequency q > 0 (wave number). The density modulation parameter (maximum amplitude), β, is obtained after averaging over the power spectrum of 200 snapshots taken every 10 time steps. Please note that non-zero values of β may also indicate other density modulation besides traveling bands, as for example formation of dense clusters in the Voronoi model for small dynamical noise (see Fig. 6). However, β ≳ 0.5, typically indicates band formation.
Order parameter fluctuations and error bars. There are two kinds of fluctuations that affect the value of the polar order parameter in our system. One stems from different realizations of obstacles in the environment, the so-called quenched disorder, the other is due to fluctuations in particles orientation, that is dynamic noise. The variation of the polar order r(t) due to dynamic noise can be measured through its standard deviation, which will be correlated with the intensity of the dynamical noise. In the context of heterogeneous environments, the error of the (time-averaged) polar order parameter due to different realizations of the quenched disorder is the important quantity. The error bars in Figs. 3 and 4, are calculated from four and five different realizations of random obstacle fields per parameter point, respectively.
From the non-spatial rewiring model to rewiring in small systems. In order to show that rewiring can enhance (orientational) order, we consider a simple nonspatial model. A system of small number, N s = 10, of spins with initial random orientations is considered. The system is such that at each time step there is only one link connecting two spins, i and j. These two spins align with these rules, _ θ i ¼ γ sinðθ j À θ i Þ and _ θ j ¼ γ sinðθ i À θ j Þ, while there is a random contribution ηξ s (t) to the orientation of all the spins (s) in the system. The link between i, j remains for m time steps, then another two spins are selected randomly to interact. With this simple model, we show that smaller m, in other words, larger rate ν = 1/t m (t m = dt ⋅ m) of rewiring a single link results in a higher polar order r for the system (see Supplementary Fig. S8a). However, in flocking models, rewiring is associated to the relative motion of the particles, which in turn is related to the level of order. To verify that rewiring is correlated to the local level of order in the flocking model, we performed a series of small size simulations that clearly show such correlation between the rewiring rate ν-which increases with the local density of particles as well as the density of obstacles -and the level of order r; see Supplementary Fig. S8b.

Data availability
The data-sets generated during the current study are available from the corresponding author on request.