Pierre Degond
email: pierre.degond@math.univ-toulouse.fr

Angelika Manhart
email: a.manhart@ucl.ac.uk

Sara Merino-Aceituno

Diane Peurichard
email: diane.a.peurichard@inria.fr

Lorenzo Sala
email: lorenzo.sala@inria.fr

How environment affects active particle swarms: a case study

We investigate the collective motion of self-propelled agents in an environment filled with obstacles that are tethered to fixed positions via springs. The active particles are able to modify the environment by moving the obstacles through repulsion forces. This creates feedback interactions between the particles and the obstacles from which a breadth of patterns emerges (trails, band, clusters, honey-comb structures,...). We will focus on a discrete model first introduced in [17] and derived into a continuum PDE model. As a first major novelty, we perform an in-depth investigation of pattern formation of the discrete and continuum models in 2D: we provide phase-diagrams and determine the key mechanisms for bifurcations to happen using linear stability analysis. As a result, we discover that the agent-agent repulsion, the agent-obstacle repulsion and the obstacle's spring stiffness are the key forces in the appearance of patterns, while alignment forces between the particles play a secondary role. The second major novelty lies in the development of an innovative methodology to compare discrete and continuum models that we apply here to perform an in-depth analysis of the agreement between the discrete and continuum models.

Introduction

Understanding how patterns in collective motion arise from local interactions between individuals is an exciting and challenging endeavour that has drawn the attention of the scientific community [START_REF] Ben-Jacob | Cooperative self-organization of microorganisms[END_REF][START_REF] Bernoff | Biological aggregation driven by social and environmental factors: a nonlocal model and its degenerate Cahn-Hilliard approximation[END_REF][START_REF] Cavagna | Scale-free correlations in starling flocks[END_REF][START_REF] Boissard | Trail formation based on directed pheromone deposition[END_REF][START_REF] Degond | Macroscopic models of collective motion with repulsion[END_REF][START_REF] Helbing | Self-organized pedestrian crowd dynamics: Experiments, simulations, and design solutions[END_REF][START_REF] Noselli | Swimming Euglena respond to confinement with a behavioural change enabling effective crawling[END_REF][START_REF] Schoeller | From flagellar undulations to collective motion: predicting the dynamics of sperm suspensions[END_REF]. In many scenarios the environment plays a key role in the emergence of collective motion and of the resulting patterns [START_REF] Chepizhko | Optimal noise maximizes collective motion in heterogeneous media[END_REF][START_REF] Kamal | Enhanced locomotion, effective diffusion and trapping of undulatory micro-swimmers in heterogeneous environments[END_REF][START_REF] Lo | Cell movement is guided by the rigidity of the substrate[END_REF][START_REF] Majmudar | Experiments and theory of undulatory locomotion in a simple structured medium[END_REF][START_REF] Park | Enhanced Caenorhabditis elegans locomotion in a structured microfluidic environment[END_REF][START_REF] Tung | Fluid viscoelasticity promotes collective swimming of sperm[END_REF]. Examples are evacuation dynamics in the presence of obstacles [START_REF] Feliciani | Empirical analysis of the lane formation process in bidirectional pedestrian flow[END_REF][START_REF] Helbing | Self-organized pedestrian crowd dynamics: Experiments, simulations, and design solutions[END_REF][START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF], sperm dynamics in the seminal fluid [START_REF] Degond | Coupled self-organized hydrodynamics and Stokes models for suspensions of active particles[END_REF][START_REF] Tung | Fluid viscoelasticity promotes collective swimming of sperm[END_REF], swirl of fish under the presence of predators [START_REF] Carrillo | Particle, kinetic, and hydrodynamic models of swarming[END_REF], cells moving in a space filled with fibers [START_REF] Lo | Cell movement is guided by the rigidity of the substrate[END_REF] or over a substrate [START_REF] Noselli | Swimming Euglena respond to confinement with a behavioural change enabling effective crawling[END_REF],... In particular, we are interested by feedback interactions between self-propelled agents and their environment that they are able to modify. This happens, for example, (i) in the formation of paths in grass-land by active walkers [START_REF] Helbing | Modelling the evolution of human trail systems[END_REF][START_REF] Lam | Active walker models for complex systems[END_REF],(ii) in the modification of the extra-cellular matrix (fibers) by migratory cells [START_REF] Baricos | ECM degradation by cultured human mesangial cells is mediated by a PA/plasmin/MMP-2 cascade[END_REF], or (iii) in ant trail formation due to ant pheromone deposition [START_REF] Boissard | Trail formation based on directed pheromone deposition[END_REF]. In this paper, we will focus on the model introduced in [START_REF]Large-scale dynamics of self-propelled particles moving through obstacles: model derivation and pattern formation[END_REF] where collective motion happens in an environment filled with movable obstacles that are tethered to a fixed point via a spring. The authors in [START_REF]Large-scale dynamics of self-propelled particles moving through obstacles: model derivation and pattern formation[END_REF] showed that a variety of patterns are generated due to the feedback interactions between the obstacles and the self-propelled agents. Indeed, the capacity of the agents to modify their environment (i.e., to modify the position of the obstacles) is key for patterns to form.

Figure 1: Overview of the paper. It includes a summary of the scales, the models and the objects considered in this paper and introduced in [START_REF]Large-scale dynamics of self-propelled particles moving through obstacles: model derivation and pattern formation[END_REF] (first three grey lines). The blue boxes indicate the derivation of the different models and derivation assumptions. The main contributions in the paper appear in the last row corresponding to 'patterns' (at the discrete and continuum level and their correspondence) and the linear stability analysis (bottom right yellow box).

Figure 1 offers an overview of the ideas and messages of this paper. We will consider mostly two scales (marked in yellow). The reason for this is that understanding the emergent properties of collective dynamics requires to establish a link between the agent's interactions and the continuum dynamics that emerges at scales much larger than the size of the individual agents. As a consequence, it is natural to consider two different scales to investigate collective motion: a microscopic scale where the discrete dynamics of the agents can be described, and a macroscopic scale where the average/continuum behaviour of the large ensemble can be observed.

From a modelling perspective, it is natural to consider the microscopic scale, where individualbased models can describe individual-agent behaviour and their interactions. In the left column of Figure 1 we present key features of the individual-based model introduced in [START_REF]Large-scale dynamics of self-propelled particles moving through obstacles: model derivation and pattern formation[END_REF]. The model assumes that agents move trying to avoid obstacles via a repulsion force. Agents interact with each other following Vicsek-type dynamics [START_REF] Degond | An age-structured continuum model for myxobacteria[END_REF][START_REF] Ha | From particle to kinetic and hydrodynamic descriptions of flocking[END_REF][START_REF] Jiang | Hydrodynamic limits of the kinetic self-organized models[END_REF][START_REF] Vicsek | Novel type of phase transition in a system of self-driven particles[END_REF], i.e. they move at a constant speed trying to align their orientation of motion with the one of their neighbours, up to some noise, while repelling each other at short distances. The discrete system gives the time-evolution of the position of the obstacles (X i ) i=1,...,N tethered at fixed anchor points (Y i ) i=1,...,N via a spring and the position and orientation of the self-propelled agents (Z k , α k ) k=1,...,M , where X i , Y i , Z k ∈ R 2 and α k is a unit vector (see Eq. ( 1) for a full mathematical description of the system and Figure 1 for a list of the most relevant parameters). We will explore the variety of patterns that arise depending on the values of the model parameters.

However, the simulation of the discrete model becomes quickly computationally challenging for systems composed of millions of individuals. Therefore, for large-particle systems, continuum models are to be preferred since they provide information on the average behaviour and are computationally less costly (right column of Figure 1). Moreover, continuum models are the appropriate framework for studying large scale patterns and carry out mathematical analyses like linear stability analysis. The drawback is that, from a modelling perspective, they are harder to justify than individual based models. For this reason, one would like to derive the continuum dynamics from the discrete ones: this derivation validates the continuum models and provides understanding on the emergence of large-scale patterns. At the same time, during this derivation process, due to averaging and asymptotic analysis, some information on the discrete system can be lost. This rigorous derivation is precisely one of the purposes of kinetic theory. Kinetic theory has been successfully applied to the study of models like the Vicsek model [START_REF] Degond | An age-structured continuum model for myxobacteria[END_REF][START_REF] Ha | From particle to kinetic and hydrodynamic descriptions of flocking[END_REF][START_REF] Jiang | Hydrodynamic limits of the kinetic self-organized models[END_REF][START_REF] Vicsek | Novel type of phase transition in a system of self-driven particles[END_REF] and the Cucker-Smale model [START_REF] Aceves-Sanchez | Hydrodynamic limits for kinetic flocking models of Cucker-Smale type[END_REF][START_REF] Carrillo | Asymptotic flocking dynamics for the kinetic cucker-smale model[END_REF][START_REF] Cucker | On the mathematics of emergence[END_REF]. Tools from kinetic theory were applied in [START_REF]Large-scale dynamics of self-propelled particles moving through obstacles: model derivation and pattern formation[END_REF] to the discrete model described above, see second and third rows in Figure 1.

First, the authors derive the mean-field limit equation (large-particle limit N, M → ∞ for both agents and obstacles). This equation corresponds to a Kolmogorov-Fokker-Plank equation for the time-evolution of the distribution of the agents g = g(z, α, t) at position z ∈ R 2 and orientation α; and the time-evolution of the distribution of the obstacles f = f (x, y, t) at position x ∈ R 2 with anchor point at y ∈ R 2 .

Then, from the kinetic equations for these distributions, the authors in [START_REF]Large-scale dynamics of self-propelled particles moving through obstacles: model derivation and pattern formation[END_REF] obtained continuum equations for the system under some asymptotic assumptions on the parameters (right blue boxes in Figure 1). In particular, it is assumed a high stiffness of the obstacle springs, strong local agentagent repulsion and fast agent alignment. In this regime, it was shown in [START_REF]Large-scale dynamics of self-propelled particles moving through obstacles: model derivation and pattern formation[END_REF] that the obstacle density ρ f = ρ f (x, t) becomes a non-local function of the agent density ρ g = ρ g (x, t) and that the continuum model consists of a system of two non-linear non-local equations for ρ g and the local mean orientation of the agents Ω = Ω(x, t), see Eqs. [START_REF] Carrillo | Asymptotic flocking dynamics for the kinetic cucker-smale model[END_REF].

The main objective of this article is to investigate the influence of the tethered obstacles in pattern formation using the discrete and continuum models first introduced in [START_REF]Large-scale dynamics of self-propelled particles moving through obstacles: model derivation and pattern formation[END_REF]. The main contributions of this paper are listed below:

• We focus our study primarily on the continuum equations (which were analysed only in dimension one in [START_REF]Large-scale dynamics of self-propelled particles moving through obstacles: model derivation and pattern formation[END_REF]). Here we introduce 2D simulations of the continuum equations and an extensive phase diagram (Sec. 3.2) that shows the appearance of patterns depending on the value of the parameters (green box in Fig. 1). We carry out a linear stability analysis in 2D around uniform states and validate this analysis by comparing its predictions with the numerical simulations of the discrete and continuum models (right yellow box in Fig. 1).

• We document in which parameter regime the continuum equations capture the discrete patterns (bottom grey box in Fig. 1). To this aim, we propose a method to compare discrete and continuum simulations. This novel method provides an indicator of the distance between different patterns.

• Lastly, we also expand and greatly systematize the parameter exploration of the discrete model supported by a phase diagram. As a consequence, we detect two new patterns with respect to reference [START_REF]Large-scale dynamics of self-propelled particles moving through obstacles: model derivation and pattern formation[END_REF]: honeycombs structures and pinned agents states (left green box in Fig. 1).

Organisation of the paper. The paper is organized as follows: we first describe the models (discrete and continuum), including the derivation assumptions of the continuum model. Then we simulate both systems to construct two corresponding phase diagrams based on different values of the parameters. Next, to better understand pattern formation as function of the model parameters, we perform a linear stability analysis of the continuum equations around uniform states and identify bifurcation parameters controlling the formation of patterns. Finally, an innovative method is proposed to compare discrete and continuum simulations, which is used to determine in which parameter regime the continuum equations are in good accordance with the discrete dynamics. We conclude the paper with a discussion of the main results.

Modeling

Discrete dynamics

We consider as a starting point the model introduced in [START_REF]Large-scale dynamics of self-propelled particles moving through obstacles: model derivation and pattern formation[END_REF] for self-propelled particles undergoing collective motion in an environment filled with obstacles. Obstacles are tethered to a given fixed anchor point through a Hookean spring. They are characterised by their positions X i (t) ∈ R 2 over time t ≥ 0 and their anchor points Y i ∈ R 2 for i = 1, 2, . . . , N , where N is the total number of obstacles. The self-propelled particles are characterised by their positions Z k (t) ∈ R 2 and orientations α k (t) ∈ S 1 (unit circle) at time t ≥ 0, k = 1, 2, . . . , M , where M is the total number of agents. We assume that obstacles and agents interact through a given potential, as explained next.

The evolution for the obstacles (X i (t), Y i ) i=1,...,N and the agents (Z k (t), α k (t)) k=1,...,M over time is given by the following coupled system of stochastic differential equations:

dX i = - κ η (X i -Y i ) dt - 1 η 1 M M k=1 ∇ϕ (X i -Z k ) dt + 2d o dB i t , (1a)
dZ k =u 0 α k dt - 1 ζ 1 N N i=1 ∇ϕ (Z k -X i ) dt - 1 ζ 1 M M l̸ =k ∇ψ (Z k -Z l ) dt, ( 1b
)
dα k =P α ⊥ k • ν ᾱk dt + 2d s d Bk t , (1c)
where the mean direction ᾱk is defined via the mean flux J k as follows

ᾱk = J k |J k | , where J k = M j=1 |Z k -Zj |≤r A α j . (2)
Eq. (1a) gives the time-evolution for the obstacles' positions X i . The first term on the righthand side corresponds to the force generated by the Hookean spring anchored at position Y i with stiffness constant κ > 0. The tether positions Y i are given and do not change over time. The terms B i , i = 1, . . . , N are independent Brownian motions that introduce noise in the dynamics with intensity d 0 > 0. This term accounts for fluctuations in the dynamics. Finally, the second term on the right-hand side of Eq. (1a) is precisely the interaction force that couples the dynamics of the self-propelled agents with the ones of the obstacles. We assume that ϕ is an even and non-negative interaction potential. Typically we will assume ϕ to be a repulsive potential to model volume exclusion between obstacles and self-propelled particles. Now, Eq. (1b) gives the time-evolution for the position of the self-propelled agents Z k . The first term on the right-hand-side of (1b) expresses that agent k moves in the orientation α k at a fixed speed u 0 > 0. The second term is the force due to the interaction potential coupling the self-propelled agents and the obstacles, as we have seen before. Finally, the last term is a repulsive force between agents given by a potential ψ which is assumed to be non-negative and even. This force is added to the model to prevent agents clustering at a single point in space and represents volume exclusion interactions between the agents [START_REF] Degond | Macroscopic models of collective motion with repulsion[END_REF].

The last equation (1c) gives the time-evolution for the orientation of the agents and corresponds to the terms appearing in the Vicsek model [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF] which is a widely used model in collective motion. The right hand side of Eq. (1c) is the sum of two competing forces: a force that tries to align the orientation of the self-propelled agents with the mean orientation of their neighbours and a noise term that opposes this alignment. The noise is given by ( Bk ) k=1,...,M which are M independent Brownian motions (also assumed to be independent from B i , i = 1, . . . , N ) and the intensity of this noise is given by the parameter d s > 0. The operator P α ⊥ k represents the orthonormal projection onto α ⊥ k (where α ⊥ k is a vector orthogonal to α k ) and the symbol ′ • ′ indicates that the stochastic differential equation has to be understood in the Stratonovich sense [START_REF] Hsu | Stochastic analysis on manifolds[END_REF]. In particular, the projection ensures that, for all times where the dynamics are defined, α k (t) remains on the sphere, i.e., |α k | = 1. The alignment force is given by P α ⊥ k ν ᾱk where ν > 0 is a positive constant and ᾱk is the average orientation of the neighbouring agents that are at distance r A > 0 from agent k, as computed in Eq. ( 2). Indeed, this term corresponds to an alignment force since it can be rewritten as

P α ⊥ k ν ᾱk = ν∇ α k (α k • ᾱk )
, where ∇ α k denotes the gradient on the sphere. Therefore, this term is a gradient flow that relaxes α k towards the average orientation ᾱk at speed ν > 0.

Finally, notice that the discrete model (1) consists of first order equations: the model can be derived from second order equations in the overdamped (or inertialess) regime. This is the reason why the parameters η > 0 and ζ > 0 appear in the system: η corresponds to the obstacle friction and ζ to the agent friction. In an inertialess regime first-order equations give a good approximation of the dynamics and this regime appears in many biological applications, in particular involving micro-agents (like sperm cells) in highly viscous environments.

As we will see in later sections, the feedback interactions between agents and between agents and obstacles give rise to a variety of patterns depending on the value of the parameters.

Continuum dynamics

When the number of agents and obstacles becomes large, it is useful to derive equations that determine the average behaviour of the discrete system (1). These 'averaged' equations correspond to continuum equations, which were derived in [START_REF]Large-scale dynamics of self-propelled particles moving through obstacles: model derivation and pattern formation[END_REF] for the discrete system [START_REF] Aceves-Sanchez | Hydrodynamic limits for kinetic flocking models of Cucker-Smale type[END_REF]. In this section we summarise the results from this reference.

Main assumptions of the derivation

The derivation of the continuum equations in [START_REF]Large-scale dynamics of self-propelled particles moving through obstacles: model derivation and pattern formation[END_REF] is done under the following set of assumptions:

(a) Large-particle-system assumption. The number of obstacles and agents are assumed to tend to infinity, i.e., N → ∞, M → ∞.

Under this assumption, the authors derived formally equations for the evolution of obstacles and agent density (kinetic equations). Then, some of the parameters of the kinetic equation are scaled by a small factor ε ≪ 1 and the continuum equations are obtained in the limit ε → 0. We explain next the scaling assumptions considered.

(b) Scaling assumptions on the parameters. Three types of scaling assumptions are made: (i) the radius of alignment of the agents is supposed to be small and scaled as r A = O( √ ε); (ii) the agent-agent repulsion distance is supposed to be small and scales as r R = O(ε), but it is ensured that the potential stays of order 1 by setting

ψ(x)dx = µ < ∞; (3)
(iii) the agents alignment rate ν and orientational noise intensity d s in (1c) are supposed to be very large and scale as:

d s , ν = O( 1 ε ) with ds ν = O(1)
: this corresponds to fast agent-agent alignment and diffusion [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF].

(c) Uniform anchor density and stiff regime assumptions. It is assumed that the anchor density for the obstacles is constant (uniformly distributed) and that the obstacles' springs are very stiff (the parameter κ is very large). To this aim, we consider the ratio

γ = η κ ≪ 1 (4)
to be small. We suppose also a low obstacle noise regime, by considering the smallness of

δ = d o γ ≪ 1. (5)
The set of assumptions (a) is sufficient to derive continuum equations. The large-particle-limit or mean-field limit gives rise to kinetic equations for the obstacle density f = f (t, x, y) and the agent density g = g(t, z, α). The set of assumptions (b) and (c) are sufficient to obtain closed equations for the obstacle density ρ f = ρ f (t, x), the agent density ρ g = ρ g (x, t) and the meanagent orientation Ω = Ω(x, t). In particular, the scaling assumptions r A = O( √ ε) and r R = O(ε) imply that alignment and agent-agent repulsion forces become localized in space as ε → 0. The set of assumptions (c) is used to Taylor expand the function f with respect to γ and δ.

In summary, the continuum equations approximate a system with a very large number of agents and obstacles in the regime where the parameters of the system reach a given range of values, as described above, i.e., in the regime ε → 0 (by an asymptotic analysis) and γ ≈ 0, δ ≈ 0 (by a Taylor expansion approximation). These approximations will be taken into account when comparing discrete and continuum simulations, since they determine the range of validity of the continuum dynamics.

The continuum model

The authors in [START_REF]Large-scale dynamics of self-propelled particles moving through obstacles: model derivation and pattern formation[END_REF] obtain the following equations for the dynamics of the density of agents ρ g (x, t) ∈ R and their mean orientation Ω(x, t) ∈ S 1 at a point x ∈ R 2 at time t ≥ 0:

∂ t ρ g + ∇ • (U ρ g ) = 0, (6)
ρ g ∂ t Ω + ρ g (V • ∇) Ω + d 3 P Ω ⊥ ∇ρ g = γ s P Ω ⊥ ∆(ρ g Ω),
where

U = d 1 Ω - 1 ζ ∇ρ f - µ ζ ∇ρ g , V = d 2 Ω - 1 ζ ∇ρ f - µ ζ ∇ρ g ,
where ρ f (x, t) is the obstacle density given by:

ρ f /ρ A = 1 + 1 κ ∆ρ g + 1 κ 2 N (ρ g ) - η κ 2 ∂ t ∆ρ g + O η κ 3 , N (ρ g ) := detH(ρ g ), (7)
where ρ A is the distribution of the anchor points in space (assumed to be constant and here taken to be equal to 1 in the simulations and computations); H denotes the Hessian, 'det' denotes the determinant, and we have defined

ρ := ρ * ϕ, (8)
the convolution between ρ and ϕ, where ϕ is the repulsion kernel between agents and obstacles, Eq. [START_REF] Degond | Macroscopic models of collective motion with repulsion[END_REF]. In the numerical simulations we will drop the higher order terms in η/κ for ρ f . The model parameters are the friction constants ζ, η, the obstacle-spring constant κ, and the agentagent repulsion intensity µ given by Eq. ( 3).

The friction coefficient γ s reads

γ s = r 2 A 8 d s ν + c 2 . (9)
The constants d 1 , d 2 and d 3 are defined by

d i = u 0 c i , (10)
where u 0 is the agent speed, and c 1 , c 2 and c 3 are explicit constants that depend only on the fraction d s /ν:

c 1 = 2π 0 cos θ m(θ) dθ, (11a)
c 2 = π 0 sin 2 θ cos θ m(θ)h(θ) dθ π 0 sin 2 θ m(θ)h(θ) dθ , (11b)
c 3 = d s /ν, (11c)
where

m(θ) = 1 Z exp ν d s cos θ , Z := 2π 0 exp(ν cos θ/d s ) dθ.
and where the function h does not have a explicit form but it is the solution to a differential equation. Specifically, h(θ) = g(θ)/ sin(θ) where g is the unique solution (for the exact functional space in which this unique solution is defined, the reader is referred to [12, Lemma 2.3])

ν d s sin θ dg dθ + d 2 g dθ 2 = sin θ.
For an explanation on the meaning of these equations the reader is referred to [START_REF]Large-scale dynamics of self-propelled particles moving through obstacles: model derivation and pattern formation[END_REF]. We just point here that the system (6) for (ρ g , Ω g ) corresponds to the so-called Self-Organised Hydrodynamics with Repulsion (SOHR) [START_REF] Degond | Macroscopic models of collective motion with repulsion[END_REF] in the case where ∇ x ρf = 0 (i.e. when there is no influence from the obstacles). The SOHR is the continuum version of the Vicsek model with agent-agent repulsion [START_REF] Degond | Macroscopic models of collective motion with repulsion[END_REF].

Remark 1 (Approximation for ρ f and blow-up). The density ρ f may take negative values: in that case the continuum simulations will bestopped. Notice also that solutions may 'blow-up' in the sense that particle densities may concentrate at points in space. We here show some simulations of the discrete model [START_REF] Aceves-Sanchez | Hydrodynamic limits for kinetic flocking models of Cucker-Smale type[END_REF] to give an overview of the different types of patterns that emerge depending on the values of the parameters. Simulations are performed with N = M = 3000 agents and obstacles initially distributed uniformly in the periodic domain U = [0, 1]×[0, 1]. We also suppose that anchor points Y k for the obstacles are uniformly distributed in U , and fix the initial agent direction to π/4.

We consider the following expressions for the agent-agent and agent-obstacle repulsion potentials:

ψ(x) = 6µ πr 2 R 1 - |x| r R 2 + , ϕ(x) = 3C ϕ 2πτ 1 - |x| τ 2 + , (12)
where

x 2 + = x if x 2 ∈ R + , 0 if x < 0.
Therefore, both potentials are compactly supported and act in a radius r R > 0 for agent-agent repulsion and a radius τ > 0 for agent-obstacle repulsion. Notice that the constants have been chosen such that µ = ψ(x)dx and C ϕ = |∇ϕ|(x)dx.

We fix a set of parameters as described in Table 1, and focus our study on the interplay between three parameters: the obstacle spring stiffness κ, the agent friction ζ and the agent-agent repulsion intensity µ.

Phase diagram

Fig. 2 shows the output of the simulations at time t = 10: at this time agents and obstacles patterns seem to have reached a steady state. In this figure, agents' positions and their orientations are represented with black arrows and obstacle's positions with blue dots. The output of the simulations are grouped into three panels: panel (A) corresponds to weak obstacle spring stiffness κ = 10, and panels (B) and (C) correspond to mild κ = 100 and strong κ = 1000 obstacle spring stiffness, respectively. Inside each panel, we arrange the simulations in a table: right-to-left columns correspond to increasing values of the agent-agent repulsion force µ, bottom-to-top rows correspond to increasing values of the friction coefficient ζ. Notice that the value for the agentagent repulsion force µ is not taken the same in all panels. Indeed, the values for µ selected are the ones that make different patterns appear in the simulations. We will justify further the particular choice of the parameters after the linear stability analysis of the continuum equations. Notice that the values for ζ are also different in panel (C). We refer the reader to the caption of Fig. 2 for the exact choices for the parameter values of µ and ζ. Finally, we point out that the figures marked with a red cross are the ones for which the videos can be found in the supplementary material (see Appendix A for more details).

From Fig. 2 we observe that a rich variety of agents' patterns emerges when varying the spring stiffness κ, the intensity of the agent-agent repulsion µ, and the friction coefficient ζ.

We classify these patterns into 4 main types and we outline the parameter regions corresponding to each with frames of different colors in Fig. 2: • Trails of agents (framed in red): agents organize into trails inside the obstacle pool. This behavior is mainly observed for weak and mild obstacle spring stiffness (κ = 10, panel (A) and κ = 100, panel (B) of Fig. 2, respectively)

• Honeycomb organization of the agents (framed in orange): For small obstacle spring stiffness κ = 10 (panel (A) of Fig. 2) and mild agent-agent repulsion µ > 0.1 (middle columns), we observe that the agents organize into fixed honeycomb structures, framing the obstacles which concentrate into aggregates of different sizes and shapes (not necessarily round). We point out that this pattern was not detected in the previous publication [START_REF]Large-scale dynamics of self-propelled particles moving through obstacles: model derivation and pattern formation[END_REF].

• Travelling bands of agents (framed in yellow): only observed for large values of the obstacle spring stiffness κ = 10 3 and large agent friction with the environment ζ = 5, here the agents organize into bands perpendicular to their direction of motion. The width of the bands increases with the agent-agent repulsion intensity µ (from left to right plots of the first row of panel (C)).

• Clusters of agents (framed in green): agents organize into clusters more or less round depending on the regime of parameters. Cluster formation appears in all regimes of obstacle spring stiffness κ = 10, 10 2 , 10 3 (all three panels), and the size of the clusters changes depending on the obstacle spring stiffness κ and on the agent-agent repulsion intensity µ but seems independent of the agent friction ζ. Particularly, we observe that the cluster sizes increase with µ, until a point is reached in which µ is so large that agent-agent repulsion counteracts all the other aggregation forces (right columns of Fig. 2). Moreover, the parameter µ acts as a phase transition parameter between different types of patterns. During the transition from clustered to near-homogeneous agent distributions with increasing repulsion intensity µ, we observe a passage to other pattern types such as trails (for weak κ = 10 or mild κ = 100 obstacle stiffness), or honeycomb organizations (for weak obstacle stiffness). Finally, we note that for large obstacle spring stiffness κ and small agent friction ζ (bottom row of panel (C), simulations marked with a green star) we observe the formation of 'pinned' clusters where the agents are grouped into very small clusters that do not move (see Supplementary material and Appendix A for access to the videos)

Each of these agent patterns is surrounded by obstacles that are kept at a given distance from the agents. This distance depends on the stiffness of the obstacles' springs κ and the agent-agent repulsion intensity µ. On one hand, if obstacles are loose enough (i.e., κ is small), the repulsion force between the agents and the obstacles may be large enough to keep them both at approximately the obstacle-agent repulsion distance τ (defined in the potential ϕ, Eq. ( 12)). On the other hand, agent-agent repulsion opposes this effect, by giving the agent population force to go against the pressure exerted by the obstacle pool. We indeed observe that increasing the agent-agent repulsion force µ (left-to-right columns of Fig. 2) decreases the typical distance between the agent structures and the obstacles.

Continuum dynamics

In this section we show numerical simulations of the continuum equations (6) using the numerical scheme detailed in Appendix C.

Simulation set up

We perform simulations of the continuum model on the periodic domain U = [0, 1] × [0, 1] discretized with space step ∆x ≈ 6.7 10 -3 (150 discretization points in each direction). The initial homogeneous agent direction Ω 0 is set to π 4 , and initial agent density ρ g is a small perturbation of a uniform distribution with Ω ρ g = 1. In order to compare the numerical results with the discrete model, we use the same parameters as for the discrete simulations presented in Section 3.1 (see Table 2).

Notice that the agent-agent alignment distance at the continuum level is chosen to be r A = 0.15 whereas for the discrete simulations it was 0.1. This choice corresponds to having rescaled r A approximately by a scaling factor ε = 0.5, i.e., r ′ A = √ εr A , where r ′ A = 0.1 is the parameter used in the discrete simulation (see the scaling assumption (b) in Sec. d s /ν is relevant for the continuous model, independently of their individual values. We therefore just ensure that this ratio is kept the same as for the discrete simulations and use d s /ν = 0.01.

Phase diagram

We present the output of the continuum simulations. To facilitate the comparison with the discrete system, we adopt the same representation as the one presented in Fig. 2. In particular, the continuum densities are discretized as follows: at a simulation time t we distribute randomly N = 3000 agent points in the domain according to the distribution ρ g (•, t), and similarly for the obstacle points using ρ f (•, t). In Fig. 3 we show the simulation results at the final time of the simulation, corresponding either to the time before blow-up or appearance of negative density for the obstacles (see Rem. 1) or to t = 10, as for the discrete simulations. As in Fig. 2, the simulations are separated in three panels: panel (A) is obtained for weak obstacle stiffness κ = 10, and panels (B) and (C) are for κ = 100 and κ = 1000 respectively. In each panel, we organize the simulations in tables for which bottom-to-top rows correspond to increasing values of the friction coefficient ζ, while left-to-right columns correspond to increasing values of the agent-agent repulsion force intensity µ. See the legend of Fig. 3 for more details on the parameter values considered for ζ and µ.

In Fig. 3 we observe different patterns for the agents, each framed using the same color code as for the discrete simulations: cluster formation (framed in green, present in all three panels), travelling bands (framed in yellow, panel (C), trails (framed in red, panel (B)), near-honeycomb structures (framed in orange, panel (A)), uniform distributions (unframed) and in-between states. Here again, increasing the obstacle spring stiffness κ (from top to bottom panels) decreases the distance between agents and obstacles (i.e., the white area around the agents is reduced with increasing κ). We also observe that increasing the agent-agent repulsion intensity µ increases the size of the agent clusters and this parameter again serves as a transition parameter between clusters and uniform distribution of the agents, passing through honeycomb structures (first row of panel (A)), trails (third row of panel (B)) or travelling bands (first three rows of panel (C)). The effect of the friction parameter ζ becomes more relevant for large values of κ. For example, in panel (C) the parameter ζ serves as a transition parameter between clusters, trails and uniform states.

Comparing phase diagrams. We compare the two phase diagrams from the discrete simulations in Fig. 2 and the continuum simulations in Fig. 3. Note, though, that there is not an exact correspondence of the values for the parameter ζ used in panel (C) for the two cases.

It is noteworthy that the patterns observed with the continuum simulations are similar to the patterns of the discrete simulations (Fig. 2) for strong and mild obstacle spring stiffness (compare panels (B) and (C) of Figs. 3 and2), while the two models lead to different types of behavior in the weak obstacle stiffness regime (panel (A)). These are expected results since the continuum model has been obtained in a strong obstacle spring stiffness regime (1/κ ≈ 0). As a result, the continuum model seems to be unable to produce the rich variety of patterns offered by the discrete model when considering loose obstacles. Also, we do not observe the pinned state with the continuum model, which appeared with the discrete dynamics when considering large obstacle spring stiffness κ and small agent friction ζ. Even though pinned-states are observed for large values of κ, they correspond to states where agents collapse into a very small cluster and then the numerical simulations of the continuum equations blow-up due to a high concentration of the agent density ρ g (see Rem. 1).

Linear stability of uniform states 4.1 Analysis of the continuum model

Continuum equations are amenable to linear stability analysis around constant solutions or uniform states. This is useful because the presence of instabilities signals the formation of patterns. In this section we obtain an explicit condition for the stability of uniform states.

Before stating the main result, we introduce the following notation: denote by φ the Fourier transform of ϕ defined as, for k ∈ R 2 :

φk := φ|k| = φ(k) = R 2 e -ik•x ϕ(x)dx ∈ R.
Notice that ϕ is assumed rotationally invariant, therefore φ is real and rotationally invariant (so we abused notation and wrote φ|k| instead of φk ).

Theorem 1 (Linear instability). Consider fixed constant values ρ 0 > 0 and Ω 0 ∈ S 1 . Then, the linearized system of (6) around (ρ 0 , Ω 0 ) is unstable if and only if there exists z > 0 such that z 2 ( φz ) 2 > µκ.

(

) 13
The proof of the theorem is given later. First, we derive sufficient conditions for the system to be stable: Corollary 1 (Conditions for stability). Suppose that ϕ is absolutely continuous, rotationally invariant, and ϕ, ϕ ′ ∈ L 1 . Then, it holds that

c 0 := max z∈R+ z 2 ( φz ) 2 < ∞ (14)
and if µκ > c 0 , then the continuum equations (6) are linearly stable. Moreover, if ϕ is given by [START_REF] Degond | Macroscopic models of collective motion with repulsion[END_REF], define c ′ 0 = c 0 /C 2 ϕ . It holds that the constant c ′ 0 is independent on the obstacle-agent repulsion radius τ and the intensity C ϕ and the system is stable whenever

µκ C 2 ϕ > c ′ 0 .
Proof. Since by assumption ϕ is absolutely continuous and ϕ, ϕ

′ ∈ L 1 , we have that | φ′ (k)| = |k|| φ(k)|. Moreover, since ϕ ′ ∈ L 1 , then φ′ is bounded. Therefore, |k| 2 | φ(k)| 2 is
bounded and c 0 is finite. In this case, for µκ > c 0 the instability condition [START_REF] Degond | An age-structured continuum model for myxobacteria[END_REF] does not hold, so the system is stable.

In the particular case where ϕ takes the shape given in [START_REF] Degond | Macroscopic models of collective motion with repulsion[END_REF], one can check that the following self-similarity condition holds

|k| φk = τ|k| φ(1) (τ|k|),
where ϕ (1) corresponds to ϕ when taking τ = 1. Therefore, it holds that

C 2 ϕ c 0 = max k |k| 2 ( φk ) 2 = max k (τ|k|) 2 ( φ(1) (τ|k|)) 2 = max y |y| 2 ( φ(1) (|y|)) 2 ,
and so c 0 is independent of τ. The rest of the corollary follows: the value of c ′ 0 is also clearly independent of C ϕ as it is just a multiplicative factor of ϕ.

Remark 2 (Limiting case of pillar obstacles). In the case where the obstacles are fixed pillars, i.e., the case where κ → ∞, then the uniform distribution of agents and pillars is always a stable solution. The effect of this limiting case is that the equations for the agents on (ρ g , Ω) become decoupled from the obstacles' density ρ f = ρ A , which is just constant (take the formal limit κ → ∞ on the continuum equations (6)). Therefore, there is an abrupt behavioural change between static obstacles and obstacles that can move a bit (anchored at a fixed point via a very stiff spring). This shows that, in this particular set up, the fact that the agents are able to modify their environment is crucial for interesting patterns to emerge.

The role of the parameters. From the instability condition (13), we observe that the main drivers of the formation of instabilities are: the shape of the agent-obstacle repulsion potential ϕ, the obstacle-spring stiffness κ, and the agent-agent repulsion intensity µ. High agent-agent repulsion -high values of µ -has a stabilising effect while high agent-obstacle repulsion -high values of C ϕ -has a destabilising effect, and vice-versa. Also, high values of the spring constant κ have a stabilising effect and small values have the opposite effect.

From Cor. 1, in the case when ϕ is given in [START_REF] Degond | Macroscopic models of collective motion with repulsion[END_REF] the ratio given by

b p := µκ C 2 ϕ c ′ 0 ( 15
)
is the single value that acts as a bifurcation parameter. However, the obstacle-agent repulsion radius τ plays a role in determining the size of the patterns (see Fig. 3). Also, from the instability condition (13) and corollary 1, for typical shapes of the potential ϕ, we expect to have stability for small and large values of the wave vector k but instabilities can appear at intermediate values whenever c 0 > µκ (where c 0 is given in ( 14)).

The rest of this section is devoted to the proof of Th. 1.

Proof of Th. 1. We start by linearising the continuum equations ( 6) around (ρ 0 , Ω 0 ) by expanding the solution using a small perturbation parameter β

ρ g = ρ 0 + βρ 1 + O(β 2 ), Ω = Ω 0 + βΩ 1 + O(β 2 ), |Ω| = 1. ( 16
)
Dropping the higher order terms, we obtain the linearised system (where the over-script bar notation is defined in Eq. ( 8)):

∂ t ρ 1 + d 1 Ω 0 • ∇ρ 1 + d 1 ρ 0 ∇ • Ω 1 = μρ 0 ∆ρ 1 + ρ 0 λ ∆ 2 ρ1 -γ∆ 2 ∂ t ρ1 , (17a)
ρ 0 ∂ t Ω 1 + ρ 0 d 2 (Ω 0 • ∇) Ω 1 + d 3 P Ω ⊥ 0 ∇ρ 1 = γ s ρ 0 P Ω ⊥ 0 ∆Ω 1 , (17b)
Ω 0 • Ω 1 = 0, ( 17c
)
where ∆ 2 is the bi-Laplacian, i.e., ∆ 2 ρ = ∆(∆ρ) and P Ω ⊥ 0 is the orthogonal projection on Ω ⊥ 0 . Note also that μ = µ/ζ, γ = η/κ and λ = ρ A /(κζ) (we assume ρ A = 1).

We now define the functions F, G : R + → R by:

F (z) := z 2 ρ 0 ζ 1 κ z 2 ( φz ) 2 -µ , (18)
G(z) := 1 + ρ 0 η κ 2 ζ z 4 ( φz ) 2 > 0,
and given k ∈ R 2 , we denote by k 0 , k 1 the quantities

k 0 = (k • Ω 0 ), k 1 = (k • Ω ⊥ 0 ), (19)
where Ω ⊥ 0 is the image of Ω 0 by the rotation of angle π/2. Th. 1 is then a direct consequence of the following proposition. Proposition 1. System (17) allows for non-trivial plane wave solutions, i.e. solutions of the form

ρ 1 (x, t) = ρe ik•x+αt , Ω 1 (x, t) = Ωe ik•x+αt , (20)
where k ∈ R 2 is the wave vector, α ∈ C, ρ ∈ C and Ω ∈ C 2 , and (ρ, Ω) ̸ = (0, 0) if and only if α and k fulfil the following dispersion relations:

Case A: k ∥ Ω 0 Option 1: ρ ̸ = 0, Ω = 0, α = α 1 (k) := -i d 1 k 0 G(|k 0 |) + F (|k 0 |) G(|k 0 |) . ( 21
)
Option 2: ρ = 0, Ω ̸ = 0.

α = α 2 (k) := -id 2 k 0 -|k| 2 γ s . ( 22
)
Case B: k ∦ Ω 0 . Then, α is a root of the following polynomial of degree 2:

α 2 G + α G|k| 2 γ s -F + ik 0 (Gd 2 + d 1 ) (23)
+ d 1 (ρ 0 d 3 k 2 1 -d 2 k 2 0 ) -|k| 2 γ s F + i d 1 k 0 |k| 2 γ s -d 2 k 0 F = 0.
The real parts of α are negative if and only if the following holds:

G(k)|k| 2 γ s -F (k) > 0, (24)
and

H(k) := G(k)|k| 2 γ s -F (k) 2 d 1 d 3 k 2 1 ( 25
) -γ s F (k)|k| 2 (d 1 -d 2 G(k)) 2 k 2 0 + G(k)|k| 2 γ s -F (k) 2 > 0.
Proof of Prop. 1. Substituting the plane wave ansatz into the equation yields

ρα + iρd 1 (Ω 0 • k) + iρ 0 d 1 Ω • k = -|k| 2 μρ 0 ρ + |k| 4 λρ 0 ρ( φk ) 2 (1 -γα) , (26a)
ρ 0 α Ω + iρ 0 d 2 Ω (Ω 0 • k) + iρd 3 P Ω ⊥ 0 k = -|k| 2 ρ 0 γ s Ω, (26b)
Ω 0 • Ω = 0. (26c) or (if Ω = ωΩ ⊥ 0 ) (G(|k|)α -F (|k|) + id 1 k 0 )ρ + iρ 0 d 1 k 1 ω = 0, id 3 k 1 ρ + ρ 0 (α + id 2 k 0 + |k| 2 γ s )ω = 0.
This is a homogeneous linear system in (ρ, ω) which has a non-trivial solution if and only if the determinant of the system is 0, i.e.:

(G(|k|)α -F (|k|) + id 1 k 0 ) α + id 2 k 0 + |k| 2 γ s + d 1 d 3 k 2 1 = 0. ( 27
)
If k 1 = 0, there are two roots corresponding to either bracket being zero. This leads to [START_REF] Helbing | Self-organized pedestrian crowd dynamics: Experiments, simulations, and design solutions[END_REF] or [START_REF] Helbing | Modelling the evolution of human trail systems[END_REF]. If k 1 ̸ = 0, we can recast [START_REF] Lo | Cell movement is guided by the rigidity of the substrate[END_REF] in [START_REF] Hsu | Stochastic analysis on manifolds[END_REF].

To determine the sign of the real part of α, we use the Routh-Hurwitz criterion for polynomials with complex coefficients [START_REF] Gantmacher | Applications of the Theory of Matrices[END_REF][START_REF] Morris | The routh and routh-hurwitz stability criteria: Their derivation by a novel method using comparatively elementary algebra[END_REF]. In our case the Routh-Hurwitz criterion states that the Re(α) < 0 for all solutions α if and only if expressions ( 24) and ( 25) hold.

With Prop. 1 we conclude the proof of Th. 1 as follows. Suppose (13) holds and let z 0 > 0 be such that

z 2 0 φ2 z0 > µκ. Let k = z 0 Ω 0 . Then k 0 = z 0 and k 1 = 0. So F (|k|) = F (z 0 ) > 0 and α = α 1 (k) is such that Re (α) > 0.
Hence, the linearized system is unstable.

Suppose now [START_REF] Degond | An age-structured continuum model for myxobacteria[END_REF] does not hold, i.e., z 2 φ2 z < µκ, for all z ∈ R + . Then, F (|k|) < 0, for all k ∈ R 2 . It results that Re(α 1 (|k|)) < 0, Re(α 2 (|k|)) < 0. Furthermore ( 24) and ( 25) are obviously satisfied for all k ∈ R 2 . Hence the system is stable.

Numerical validation of the linear stability analysis

In this section we compare the pattern predictions given by the linear stability analysis with the results obtained from numerical simulations. This way we check that the linear stability analysis truly captures pattern formation, i.e., that nonlinear effects are of second order and most of the patterns characteristics are captured by linear effects.

Predictions from the theoretical analysis and qualitative agreement with the macroscopic simulations

We start by giving insights on the size and shape of the expected patterns based on the theoretical predictions offered by the stability analysis. To this aim, we consider perturbations introduced in the stability analysis (see Prop. 1), around the homogeneous density ρ 0 = 1 and in constant direction Ω 0 ∈ S 1 . As we are particularly interested in characterizing the patterns corresponding to clusters or bands, we will focus on the theoretical values for wave vectors parallel to Ω 0 and parallel to Ω ⊥ 0 :

k th ∥ = argmax k∥Ω0 Re(α(k)), k th ⊥ = argmax k∥Ω ⊥ 0 Re(α(k)),
where α(k) corresponds to case A (Eq. ( 21)), α(k) corresponds to case B (larger root of Eq. ( 23), computed numerically) and the symbol 'Re' indicates the real part. With these wave vectors maximizing the real part of α, we define the quantities:

S th 1 = 2π |k th ∥ | , S th 2 = 2π |k th ⊥ |
.

These quantities give the size of the expected patterns in each direction. We will also compute the maximal growth rates of the perturbations in these two directions:

α ∥ max = max k∥Ω0 Re( α(k)), α ⊥ max = argmax k∥Ω ⊥ 0 Re(α(k)).
Equipped with these quantifiers, we now study the influence of the model parameters on the expected pattern shapes and sizes. As predicted by the stability analysis, patterns can be expected if the bifurcation parameter b p (Eq. ( 15)) is below 1. We fix C ϕ = 5 and use ϕ as in Eq. ( 12) giving c 0 ≈ 5.6 independent on τ as shown in the proof of corollary 1 (definition of c 0 in Eq. 14). We vary b p by changing the values of the agent-agent repulsion intensity µ and aim to study the influence of the friction constant ζ, the obstacle spring stiffness κ and the agent-obstacle repulsion distance τ. For each subsection, we compare qualitatively these predictions based on the linear stability analysis with simulations of the macroscopic model presented in Fig. 3 Influence 4 (left panel) that we indeed recover the critical value 1 of the bifurcation parameter, below which perturbations grow (Re(α) > 0) and after which they are damped, independently on the value of ζ. This shows that b p is indeed a relevant bifurcation parameter. Moreover, one can observe that perturbations grow faster for smaller values of the friction constant ζ (compare the blue and red curves in the left panel). From the right panel of Fig. 4, we note first that the size of the clusters increases when increasing the bifurcation parameter (here, by increasing the agent-agent repulsion µ). These are expected results as stronger agent repulsion leads to higher pressure in the agent population, leading to larger clusters. Secondly, we observe that the size of the patterns is independent on the friction constant ζ, but the parameter zone in which patterns are of travelling band type (i.e S th 1 > 0 and S th 2 = 0) is larger for larger values of ζ -compare the yellow and blue dashed curves in the right panel-. Thus, high friction substrates seem to favor the formation of travelling bands compared to low friction environments, provided the bifurcation parameter is large enough (large obstacle spring stiffness and/or large agent-agent repulsion compared to agent-obstacle repulsion).

Qualitative comparison with the macroscopic simulations. The influence of the agent friction ζ for τ = 0.15 and κ = 1000 can be observed in the macroscopic simulations presented in Fig. 3 panel (C), comparing the rows together (from bottom to top for increasing values of ζ). We first note that in the simulations of the three panels of Fig. 3, the values considered for the product µκ were always the same, i.e., µκ ∈ {0.02, 0.2, 0.5, 1, 4, 6}, and the value of C ϕ = 5 was kept constant, corresponding to the following values for the bifurcation parameter:

b p ∈ {0.0036, 0.0357, 0.0893, 0.1785, 0.7142, 1.0713}.

We then observe that in each panel of Fig. 3, patterns are indeed observed in the first 5 columns of the tables while the last column displays a homogeneous distribution of agents. This validates the fact that patterns are observed only when the bifurcation parameter b p is below 1. Moreover, focusing on the last panel (for which κ = 1000), we recover most of the observations predicted by the stability analysis: (a) the pattern size increases when increasing the bifurcation parameter (increasing µ: compare simulations from left to right in Fig. 3 Influence of the obstacle spring stiffness κ. Here we adopt the same representation as in the previous paragraph, but fixing the agent friction constant ζ = 0.5 and playing on the obstacle spring stiffness κ (we keep the agent-obstacle distance τ = 0.15). Fig. 5 shows the values of α ∥ max and α ⊥ max (left panel) and of S th 1 and S th 2 (right panel), as functions of the bifurcation parameter b p and for κ = 10 (blue curves), κ = 100 (red curves) and κ = 1000 (yellow curves). From Fig. 5 (right), we can observe a similar evolution of the pattern size playing on the obstacle spring stiffness as when changing the friction constant ζ: increasing the obstacle spring stiffness κ slightly increases the zone of parameters favoring the formation of bands of agents (compare yellow and red dashed curves in the right panel). One can particularly note (blue curve of Fig. 5 (right)) that environments composed of loose obstacles (κ = 10) will only promote agent clusters the size of which is independent of the value of the bifurcation parameter. Finally we note from Fig. 5 (left) that the growth rate of perturbations does not evolve monotonically with the spring stiffness κ: faster perturbations are observed for κ = 100 compared to κ = 10 or κ = 1000 (compare red with blue and yellow curves in the left panel).

Qualitative comparison with the macroscopic simulations. The influence of the obstacle spring stiffness κ for τ = 0.15 and ζ = 0.5 can be observed in the macroscopic simulations presented in Fig. 3, comparing the second rows (starting from the bottom) in each panel (panel (A) for κ = 10, panel (B) for κ = 100, panel (C) for κ = 1000).

Again, we obtain a very good agreement with the theoretical predictions: (a) the pattern sizes increase when increasing the bifurcation parameter (by increasing µ: compare simulations from left to right in each panel), (b) the increase in pattern size as function of µ seems less important for κ = 10 (panel (A)) than for larger obstacle spring stiffness (panels (B) and (C)), and (c) travelling bands are only observed for κ = 1000 (panel (C)).

Influence of the agent-obstacle repulsion distance τ

Finally we aim to document the role of the agent-obstacle repulsion distance τ. We adopt the same methodology as in the two previous paragraphs: we fix ζ = 0.5 and κ = 1000 and show in Fig. 6 the values of α ∥ max and α ⊥ max (left panel) and of S th 1 and S th 2 (right panel), as functions of the bifurcation parameter b p and for τ = 0.15 (blue curves), τ = 0.2 (red curves) and τ = 0.3 (yellow curves). We first observe that increasing the value of τ slows down the growth of the perturbation modes (compare blue red and yellow curves of Fig. 6 (left)). Moreover, as predicted by the stability analysis, the critical value of µ for which patterns appear does not depend on τ: patterns are once again only observed as long as the bifurcation parameter b p does not exceed the value 1. Secondly, Fig. 6 (right) shows that the agent-obstacle distance τ has a strong impact on the size of the clusters: larger τ leads to larger agent clusters (compare for instance yellow and blue curves in Fig. 6 (right)), and agent-obstacle repulsion distance does not impact the shape of the patterns (clusters or bands types).

As the simulations of Fig. 3 have been generated only for τ = 0.15, we are not able at this point to compare qualitatively the predictions of the stability analysis with the simulations of the macroscopic model as functions of this parameter. We will however assess the influence of τ via a quantitative comparison between the model and the theory in the next section.

Altogether, these results show that agent-agent repulsion favors the spreading of the agents while agent-obstacle repulsion tends to aggregate the agents (and consequently clusters obstacles together). Travelling bands of agents seem to be favored in low friction environments composed of stiff obstacles, and the size of agent clusters seem to be controlled primarily by the agent-obstacle distance and the bifurcation parameter (ratio between the agent-agent repulsion intensity and the agent-obstacle repulsion intensity).

Quantitative agreement between the macroscopic simulations and the stability analysis

Here we provide a quantitative assessment of the pattern sizes computed numerically on the simulations of the macroscopic model and the ones predicted by the stability analysis. To this aim, we first compute numerically the pattern sizes using the 2D Discrete Fourier transform of the agent density at equilibrium F [ρ g ] = F [ρ g ](k), and extract the frequency of the two maximal modes k ∥ , k ⊥ ∈ R 2 aligned in the direction of Ω 0 and Ω ⊥ 0 , respectively:

k ∥ = argmax k∥Ω0 | F [ρ g ](k)|, k ⊥ = argmax k∥Ω ⊥ 0 | F [ρ g ](k)|,
where | • | is the modulus of a complex number. Then, the theoretical quantifiers S th 1 and S th 2 will be compared with

S 1 = 2π |k ∥ | , S 2 = 2π |k ⊥ | .
In Fig. 7, we show the values of S 1 and S 2 (dotted curves) and S th 1 , S th 2 (plain curves), for three different values of the obstacle spring stiffness κ = 10 (panel (A)), κ = 100 (panel (B)) and κ = 1000 (panel (C)), and three different value of τ: τ = 0.15 (blue curves), τ = 0.2 (orange curves) and τ = 0.3 (yellow curves). Note that here ζ = 0.5 so that theoretical predictions correspond to Fig. 6. Simulations are completed for Ω 0 = cos π 4 , sin π 4 and ρ 0 = 1. As one can observe, we obtain a fairly good agreement between the values computed on the numerical solution and the ones predicted by the linear stability analysis as presented in Fig. 7. As predicted, the size of the repeating patterns increases as τ increases, (compare blue, red and yellow curves), and as the agent-agent repulsion intensity µ increases while staying below the critical threshold µ * (above which the homogeneous steady-state profile is stable), corresponding to b p = 1. For κ = 1000 (panel (C) of Fig. 7), we also recover the regime of travelling bands predicted for µ = 4 • 10 -3 and here S 1 > 0 and S 2 = 0, i.e., patterns (the travelling bands) are only in the direction Ω 0 .

Quantitative assessment of the continuum model

In this section, we aim to compare quantitatively the continuum and discrete models. As our goal is to compare continuous density profiles (continuum model) with clouds of points representing individual positions (discrete model), a method to quantify the 'proximity' between these two different types of solutions has to be devised. A first natural choice would be to use the quantifiers defined in the previous section, i.e. to compute the maximal eigenmode of the Fourier transform of the agent distributions from the discrete simulations. This would enable to construct a spaceindependent quantifier which could give an insight into the main structures of the discrete model. However, as one can observe in Fig. 2, the agent and obstacle structures that emerge from the discrete dynamics are not necessarily regularly spaced in the domain, which makes the use of the Fourier transform imprecise for the discrete simulations. In the following section, we propose a new method to compare discrete point clouds and continuum densities which does not require some spatial regularity of the patterns.

Methodology to compare discrete and continuum simulations

In Table 3, we summarize the steps of the method we propose to compare discrete and continuum simulations. After generating two simulations (one with the continuum model and one with the discrete dynamics, step 1), we first aim to find the optimal Cartesian mesh on which (i) we interpolate the continuum solution and (ii) we compute the density of the point clouds using a Particle In Cell (PIC) method (step 2, Section 5.1.1). At the end of this step, both solutions (continuum and discrete) are projected on the same Cartesian mesh. In step 3 (section 5.1.2), we then compute a Wasserstein-type distance based on the histograms of the two density distributions.

STEP 1

Simulation of the continuum equations [START_REF] Carrillo | Asymptotic flocking dynamics for the kinetic cucker-smale model[END_REF] Simulation of the particle dynamics (1a)-(1b) ⇓ ⇓

STEP 2

Discretization of the output density on a grid Π ∆x U using the PIC method: ρ ∆x mac Approximate the particle density on the grip Π ∆x U using the PIC method: ρ ∆x mic → Compute the optimal grid size ∆x using the ℓ 2 distance (Sec. 5.1.1)

↘ ↙ STEP 3
Compare ρ ∆x mic and ρ ∆x mac with the distance W (ρ ∆x mic , ρ ∆x mac ) → Computed with the EMD method (Sec. 5.1.2) Table 3: Diagram of the methodology used to compare the simulations for the continuum equations ( 6) and the simulation of the discrete dynamics (1a)-(1b). PIC: Particle-In-Cell; EMD: Earth Movers Distance; W : Wasserstein distance.

Discretization of the particle density

A natural choice for comparing point clouds and continuum densities is to choose a Cartesian grid for both models, and compute the density of the individual agents using for instance a Particle-In-Cell (PIC) method [START_REF] Degond | Asymptotic-preserving particle-in-cell methods for the vlasov-maxwell system near quasi-neutrality[END_REF]. However, the choice of grid points spacing is critical, as it depends on the profile of the distribution as well as on the number of particles present in the computational domain: highly clustered agent distributions require fine meshes to enable to capture the characteristics of the small and concentrated agent clusters, while more homogeneous agent distributions require coarser grids to allow the capture of larger patterns (see Fig. 8). To be efficient, the grid spacing must, therefore, account for the characteristic size of the continuum structures that can be captured with a finite number of individual points. As we want to compare a continuum model with a discrete one, we will use the continuum simulations as a reference. Our goal here is to find the optimal Cartesian mesh on which a continuum density ρ ∆x (x, t) would be best represented by a cloud of N points, for N given. Note that the continuum density ρ ∆x (x, t) is itself already discretized on a Cartesian mesh with spacing ∆x, because it corresponds to a solution of the discretized continuum model.

Given a continuum density profile ρ ∆x (x, t) -discretized on a Cartesian mesh Ω ∆x ⊂ Ω with grid spacing ∆x = 1

Nx in each direction -we first throw N individual points (y 1 , . . . , y N ) ∈ Ω according to the distribution ρ ∆x (x, t). We now denote by ρ h P IC (y 1 , . . . , y N ) the density of the individual points (y 1 , . . . , y N ) computed on a Cartesian mesh of spacing h > 0 using a PIC method, and Π ∆x ρ h P IC its linear interpolation on the initial mesh Ω ∆x . We aim at finding the optimal grid spacing h minimizing the L 2 distance between the initial continuum density ρ ∆x (x, t) and its approximation by N individual points:

h = argmin h ||ρ ∆x -Π ∆x ρ h P IC (y 1 , . . . , y N ) || ℓ 2 (Ω∆x) ,
where ||.|| ℓ 2 (Ω h ) denotes the discrete l 2 norm on a Cartesian mesh Ω h :

||ρ h || ℓ 2 (Ω h ) = h 2 Nx i=1 Nx j=1 |ρ h (x i , y j )| 2
The optimal h is computed numerically. It therefore corresponds to the best grid spacing one can hope for approximating a density ρ ∆x with a set of N points. We therefore will use this quantity to compare a simulation of the continuum model with one of the discrete model performed with N agents. The discretized macroscopic density will be denoted by ρ h mac = ρ h P IC (y 1 , . . . , y N ), and the approximation from the discrete particle simulation will be denoted by ρ h mic (computed via the PIC method on a grid with spacing h). In the following section, we describe how to compare ρ h mac with ρ h mic .

Comparing discretized and discrete dynamics

The comparison between the discretized and the discrete dynamics will be done in several steps:

Step 1) Choosing the right distance to compare the micro-and macro-simulations: we want to construct a quantifier enabling to compute the distance between the two distributions ρ h mac and ρ h mic described in the previous section (solutions of the continuum and discrete models projected on a Cartesian mesh with spacing h). The first natural choice would be to use the discrete L 2 norm as both quantities are defined on the same meshes. However, we need a quantifier independent on space translations, as there is no reason for the patterns of the discrete model to match exactly the locations of those of the continuum model at a given time. For example if the discrete and continuum simulations produce band patterns with same width and speed but not at the same positions, we still want to consider that the two solutions are very close to each other. Therefore, we propose here to use a Wasserstein-like distance. Inspired from [START_REF] Rubner | The earth mover's distance as a metric for image retrieval[END_REF], we choose to work with the Earth Movers Distance (EMD). The EMD is based on the minimal cost that must be paid to transform one distribution into the other and relies on the solution to a transportation problem issued from linear optimization. As solving the transport problem in 2 dimensions is very costly, we 'compress'/approximate the density distributions using their signatures (histograms).

Step 2) Construction of the signatures of the distributions: given a density profile on a grid containing N h = 1 h points in each direction (ρ ij ), i = 1 . . . N h , j = 1 . . . N h , the signature of ρ, P [ρ] = {(p 1 , ω 1 ), . . . , (p m , ω m ))} is defined as:

p k = kM n b , ω k = N h i=1 N h j=1 1 [p k-1 ,p k ] ρ ij , k = 1 . . . n b , (28)
where M = ||ρ|| ∞ = max i,j ρ ij and the number of bins n b has been chosen using the Freeman Diaconis rule, for which the bin width corresponds to 2 IQR n 3/2 , where IQR is the interquartile range of the data and n is the number of observations (in our case the number of grid points, n = 1 h 2 ). We give in Fig 9 a visual representation of computing the signature of a toy distribution with 4 bins and in Fig. 10 an example of the histograms of two simulations of the continuum model. Note that when computed on density distributions, the points p k in each cluster correspond to local density values and the corresponding weights ω k are the number of grid (spatial) points in which the density is comprised between the values p k-1 and p k .

Step 3) Definition of the EMD between two signatures: following the lines of [START_REF] Rubner | The earth mover's distance as a metric for image retrieval[END_REF], we apply the following linear programming problem: Let P = {(p 1 , ω 1 ), . . . , (p m , ω m ))} and Q = {(q 1 , v 1 ), . . . , (q n , v n ))} be two signatures with m and n clusters represented by their representatives p k , q ℓ and their respective weights ω k , v ℓ for k = 1 . . . m, ℓ = 1 . . . n. We want to find a flow F = (f kℓ ) minimizing the overall cost:

W (P, Q, F ) = m k=1 n ℓ=1 d kℓ f kℓ ,
where d kℓ is the ground distance matrix between clusters p k and q ℓ :

d kℓ = |p k -q ℓ |.
The minimization is made under the following set of constraints:

f ij ≥ 0, 1 ≤ i ≤ m, 1 ≤ j ≤ n, (29)
n j=1 f ij ≤ ω pi , 1 ≤ i ≤ m (30)
m i=1 f ij ≤ ω qj , 1 ≤ j ≤ n (31)
m i=1 n j=1 f ij = min( m i=1 ω pi , m i=1 ω qj ). (32)
If we look at the signatures P and Q as a set of goods at given locations (represented by p and q) each with a given amount (represented by the weights ω and v), the EMD can be seen as a transportation problem consisting in finding the least expensive flow of goods from the suppliers to the consumers, where the cost of transporting a single unit of goods is given. Then, constraint [START_REF] Maury | A macroscopic crowd motion model of gradient flow type[END_REF] expresses that 'supplies' can be transported from P to Q only, while constraints ( 30), ( 31) limits the amounts of supplies that can be given by P to Q and that can be received from Q to P , respectively. The final constraint [START_REF] Noselli | Swimming Euglena respond to confinement with a behavioural change enabling effective crawling[END_REF] expresses the fact that the total amount of mass transported must be optimal. Once this transportation problem is solved, the EMD between signatures P and Q, EM D(P, Q) is then defined as:

EM D(P, Q) = m i=1 n j=1 d ij f ij m i=1 n j=1 f ij .
Rubner et al. proved in [START_REF] Rubner | The earth mover's distance as a metric for image retrieval[END_REF] that when the ground distance is a metric and the total weights of the two signatures are equal, the EMD is a true metric. Therefore, by considering the Euclidean distance as ground distance we can use the EMD as a valid dissimilarity measure between signatures. However, as two different density distributions may have the same signature, the EMD with (28) as signatures is a pseudo-metric. However, as shown in Fig. 10, the histograms between band like patterns and clustered state are very different distributions, making this pseudo metric suitable for measuring the dissimilarity between pattern types. Moreover, we check carefully in the next paragraph the validity of the EMD when it can be compared to the classical L 2 distance.

Step 4) Validation of the pseudometric EMD. In order to check the validity of the pseudometric constructed in this section, we aim to compare the efficiency of the EMD in cases where it can be compared to the classical L 2 distance. More specifically, we use it to measure the dissimilarity between the density profile of the continuum simulation (high density clustered simulation of Fig. 8 left column) and its approximation by a cloud of N points reconstructed on a grid, using the procedure described in the previous section. We show these dissimilarity measures in Fig. 11, as a function of the number of grid points for the PIC method (N P IC , horizontal axis of Figs. 11) and different numbers N of discrete particle (different curves), using the EMD distance based on histograms (left panel) or the L 2 distance based on point values (right panel). As one can observe in Fig. 11, the EMD and the L 2 norm are in good accordance. As previously observed in Sec. 5.1.1, both metrics show that for each number of particles used to approximate the continuum density distribution, there exists an optimal number of grid points for the PIC method which minimizes the distance between the initial density and its approximation by particles. As expected, this optimal value increases as the number of particles increases, suggesting that using a larger number of agents allows the use of finer grids which enables us to better capture the fine structures of the continuum density distribution. Moreover, this figure shows that the Wasserstein distance based on the EMD between density signatures seems to be a valid tool to compare density distributions.

In the next section we present the numerical comparison between the discrete and continuum models.

Results

We aim to compare quantitatively the steady-states of the discrete and continuum models in different regimes of the parameters, and study the influence of the number of agents for the discrete model N as well as the scaling parameter ϵ. We recall that the assumptions for the derivation of the continuum equations are given in Sec. 2.2.1. In particular, some of the parameters are scaled by a factor ε ≪ 1 in the following way (denoting by a tilde the values used for discrete simulations):

rR = ϵr R , rA = √ ϵr A , ds = d s ϵ , ν = ν ϵ . (33)
For all simulations, we consider the same number of agents and obstacles and set M = N , and we fix the values of C ϕ = 5 (leading to c 0 = 5.6) and ζ = 0.5. For each set of parameters, we use the method previously described in Sec. 5.1 to compare discrete and continuum simulations.

Mild obstacle spring stiffness

In Fig. 12, we show the simulations obtained for mild obstacle spring stiffness κ = 100. The left panel is obtained for µ = 2.10 -3 (corresponding to a bifurcation parameter b p ≈ 0.036), and the right panel is for µ = 4.10 -2 (corresponding to b p ≈ 0.7, close to the stability threshold 1). Top figures show the EMD between the continuum and discrete solutions as a function of ϵ, for different number of agents used for the discrete simulations N : N = 500 (blue curve), N = 1000 (red curve), N = 3000 (yellow curve) and N = 5000 (purple curves). The corresponding simulations are shown below in tables: for each, the left column shows the simulations of the continuum model, and the next columns are simulations of the discrete model for different values of ϵ: ϵ = 0.1 (second column), ϵ = 0.5 (third column), ϵ = 0.8 (fourth column), ϵ = 1 (last column). The different rows of the tables correspond to different number of agents for the discrete simulations as well as for the discretization of the continuum density (from top to bottom: N = 500, N = 1000, N = 3000, N = 5000). Fig. 12 suggests that the discrete and continuum models are in quite good agreement in the case of week agent-agent repulsion (b p ≪ 1, left panel), where both models are able to reproduce agent clusters, while their correspondence is more tenuous for stronger agent-agent repulsion (b p close to the instability threshold, right panel), where the discrete dynamics seems to produce more trail-like patterns than the continuum model. For both regimes however, we can observe a significant improvement of the discrete-continuum correspondence as ϵ decreases, suggesting that the continuum model becomes a good approximation of the discrete dynamics as ϵ goes to zero. Indeed, for weak agent-agent repulsion (left panel), we observe that decreasing ϵ is accompanied by an increase in the cluster sizes and a decrease of the distance between the boundary of the clusters and the obstacles, getting closer to the cluster types observed with the macroscopic dynamics. For stronger agent-agent repulsion (right panel), the clusters thicken as ϵ decreases and get closer to the continuum structures.

These observations are confirmed by the measurements of the EMD between the discrete and continuum agent distributions (top plots of Fig. 12). Indeed, one notes in the left panel that the distance between the two distributions decreases as the scaling parameter ϵ decreases, independently on the number of agents. Moreover, the top plot on the right panel shows that the discrete-continuum distance is larger for stronger agent-agent repulsion (b p close to 1) compared to the case where b p ≪ 1 (left panel). From the right figure, we also observe a strong dependency of the discrete-continuum distance as a function of the number of agents used in the discrete model. When the agent-agent repulsion is strong (or equivalently when b p is close to 1), it becomes crucial to use a large number of individuals for the discrete simulations, while the number of agents does not seem to significantly impact the discrete-continuum agreement in regimes favoring the apparition of small and dense clusters (small agent-agent repulsion or equivalently small b p ).

These first observations tend to suggest that the choice of the number of agents in the discrete setting seem to depend both on the choice of ϵ and on the regime of parameters. In order to give more insights on the influence of N and b p on the discrete-continuum match, we plot in Fig. 13 the EMD between the discrete and continuum models as a function of b p (by changing the value of µ for fixed κ = 100), having fixed ϵ = 0.1 and for different N :N = 500 (blue curve), N = 1000 (red curve) and N = 3000 (yellow curve) and N = 5000 (purple curve).

As one can see in Fig. 13, the discrete-continuum distance increases with b p independently on the number of agents N , suggesting indeed that the discrete and continuum models are closer far from the instability threshold. As the agent-agent repulsion increases (increasing values of b p ), the number of agents used in the discrete simulations has increasing influence on the match between the discrete and continuum simulations. These results suggest that large agent clusters with low density are better captured by a large number of agents.

Strong obstacle spring stiffness

Here we aim to study the discrete-continuum agreement for strong obstacle spring stiffness κ = 1000. In Fig. 14, the left panel is obtained for µ = 2.10 -4 (corresponding to a bifurcation parameter b p ≈ 0.036), and the right panel is for µ = 4.10 -3 (corresponding to b p ≈ 0.7, close to the stability threshold 1). Top figures show the EMD between the continuum and discrete solutions as a function of ϵ, for different number of agents used for the discrete simulations N : N = 500 (blue curve), N = 1000 (red curve), N = 3000 (yellow curve) and N = 5000 (purple curves). As in the previous section, the corresponding simulations are shown below in tables: for each, the left column shows the simulations of the continuum model, and the next columns are simulations of the discrete model for different values of ϵ: ϵ = 0.05 (second column), ϵ = 0.1 (third column), ϵ = 0.5 (fourth column) ϵ = 0.8 (fifth column) and ϵ = 1 (last column). As before, the different rows of the tables correspond to different number of agents for the discrete simulations as well as for the discretization of the continuum density (from top to bottom: N = 500, N = 1000, N = 3000, N = 5000).

For strong obstacle spring stiffness κ = 1000, we again observe that the discrete and continuous models are in good agreement far from the instability threshold (left panel), where both models reproduce clusters, while the agreement between the two models worsen for stronger agent-agent repulsion (right panel), where the discrete system fails to reproduce the travelling bands patterns observed with the continuum model. Again, the discrete-continuum agreement improves as ϵ decreases: for low agent-agent repulsion the discrete pattern sizes converge to the ones of the continuum model as ϵ decreases (from right to left in the left panel), and for strong agent-agent repulsion (last rows of the right panel), decreasing ϵ induces a phase transition between clustered states and trail-like agent patterns, closer to the formation of bands.

It is noteworthy that for small agent-agent-repulsion (b p ≪ 1, left figure), the agreement between the discrete and continuum dynamics seem to be better when using less agents in the discrete model, independently on the value of ϵ (compare purple and blue curves on the left panel), while close to the instability threshold (b p close to 1, right figure) the choice of N seems to be related to ϵ: the discrete-continuum error decreases when using larger N for small ϵ, smaller N for larger ϵ.

Summary of observations

We conclude that the continuum equations are a good approximation of the discrete dynamics in the limit of small rescaling parameter ε, as long as the agent-agent repulsion µ is small enough (i.e in a parameter regime far from the instability threshold, b p ≪ 1). On the contrary, the trend is less apparent when µ gets closer to the instability threshold µ * (corresponding to b p = 1). In particular, when µ ≈ µ * , the rescaling factor ε can act as a phase transition parameter between different types of patterns (right panel, Fig. 14). This phase transition is due to the fact that the instability condition is given by [START_REF] Degond | An age-structured continuum model for myxobacteria[END_REF]. Indeed, the presence of µ in this formula hints to the fact that at the discrete level the agent-agent repulsion potential ψ plays a key role in determining the patterns that emerge. Therefore, it is no wonder that by rescaling the value of agent-agent repulsion radius rR = εr R (and therefore changing the value of ψ) the shape of the patterns also changes. However, the smaller the µ the less relevant the role of ψ, thus the predictions of the continuum simulations become more robust.

There is also another important factor to take into account: the continuum dynamics just gives averaged behaviour of the discrete dynamics. If there is a wide variability in the discrete dynamics, due to its intrinsic stochasticity, then the average behaviour will not be able to represent well particular realizations of the discrete dynamics. It seems that closer to the boundary of the instability region (µ ≈ µ * ) this variability is larger.

Discussion

In this article we have investigated a model for collective dynamics in an environment filled with obstacles that are tethered to a fixed point via a spring. The model was first introduced in [START_REF]Large-scale dynamics of self-propelled particles moving through obstacles: model derivation and pattern formation[END_REF]. In particular, the paper has presented the following novelties: (i) phase diagram of the continuum equations in dimension 2; (ii) a linear stability analysis of constant solutions; (iii) method to discriminate between different types of patterns that has been used to compare quantitatively the relation between discrete and continuum simulations; (iv) and, finally, a more extensive phase diagram of the discrete dynamics that has allowed to identify two new types of patterns with respect to [START_REF]Large-scale dynamics of self-propelled particles moving through obstacles: model derivation and pattern formation[END_REF] (honey comb structures and pinned cluster states).

The continuum description captures well the behavior of the system when it is comprised of a large number of agents and obstacles, and involves huge computational savings compared with the simulation of the discrete system. Comparing discrete and continuum simulations is in general not straightforward. We have proposed a method to compare the two types of solutions to investigate in which parameter regime they are in good correspondence. This parameter regime includes the assumptions made for the derivation of the continuum equations in Section 2.2.1: the spring stiffness must be large κ ≫ 1, the number of agents and obstacles must be large N, M ≫ 1, the scaling parameter ε ≪ 1 (see [START_REF] Oldroyd | On the formulation of rheological equations of state[END_REF] for the rescaled parameters) must be small. However, we require one more condition to have a good correspondence between discrete and continuum dynamics: the agent-agent repulsion intensity µ must be much smaller than the critical value µ * , which is at the threshold of the instability condition [START_REF] Degond | An age-structured continuum model for myxobacteria[END_REF]. For values closer to µ * the intrinsic variability of the system is too large to be described just with the averaged behaviour that captures the continuum equations.

This work has also showcased the impact of the environment in pattern formation in collective dynamics. The phase diagrams of both discrete and continuum dynamics show that the feedback interactions between agents and obstacles give rise to a rich variety of patterns. In particular, we have observed that trails, travelling bands, moving clusters, uniform configurations and other in-between patterns emerge. The fact that agents can modify their environment by moving the obstacles is fundamental to this pattern emergence. This can be clearly seen in the linear stability analysis where the instability condition (13) depends crucially on the agent-obstacle repulsion force ϕ which is the only interaction force between agents and obstacles, and on the spring stiffness κ which indicates the degree of mobility of the obstacles around their tethered positions.

As a prospective work, we would like to use the models investigated here to study the impact of the environment in collective dynamics under different set up. One of these set up is collective motion in a complex fluid. To investigate this, the idea is to couple the current model with a fluid model. Then the environment in which collective motion takes place will be the combination of the fluid with the obstacles. The idea of representing a complex fluid in this manner is similar to other existing models in the literature, such as the Oldroyd-B model that describes the visco-elasticity of fluids filled with spring dumbbells [START_REF] Oldroyd | On the formulation of rheological equations of state[END_REF]. The coupling of the current discrete model with a fluid model will require a new derivation of the continuum equations and a new linear stability analysis to understand how the presence of the fluid impacts the dynamics and pattern formation. Another extension of this work will investigate the impact in collective dynamics of an environment filled with a different type of obstacles (i.e., obstacles of a different nature than the ones considered in this work). For example, one can consider solid obstacles that are movable but that are not tethered or that have a particular shape (like elongated fibers). , and the second step (third and fourth rows) computes the approximated density using a PIC method with different spacing from the point distributions (third row: using a coarse grid with spacing h = 0.1, fourth row: using a finer grid h = 0.02). As one can see in Fig. 8, while the number of points to throw to approximate the continuum density does not play a major role for high density clustered distributions, it becomes critical for approximating more homogeneous distributions (compare left and right columns). Moreover, high density clusters require the use of a fine enough grid to correctly recover the initial distribution (compare third and fourth rows in the left column), while smoother distributions are better approximated using a large number of agents and coarse grids (third row of the left columns). These first results highlight the necessity for adapting the numerical grid used to compute the density of agents from the discrete model if one hopes to have a consistent quantifier to compare with the continuum model.

Increasing swimmer-swimmer repulsion Increasing friction coefficient

Continuous model code

The MATLAB code corresponding to the simulations of Fig. 3 (Continuous model) can be found at the following link: https://doi.org/10.6084/m9. figshare.19939409.v1 This supplementary material contains the following files: C Numerical method for the continuum model

One of the difficulties in solving the nonlinear model ( 6) is the geometric constraint |Ω| = 1, and the resulting non-conservativity of the model arising from the presence of the projection operator P Ω ⊥ . We rely on a method proposed in [START_REF] Motsch | Numerical simulations of a nonconservative hyperbolic system with geometric constraints describing swarming behavior[END_REF][START_REF] Degond | Macroscopic models of collective motion with repulsion[END_REF] where the SOH model is approximated by a relaxation problem consisting of an unconstrained conservative hyperbolic system supplemented with a relaxation operator onto vector fields satisfying the constraint |Ω| = 1. The relaxation model writes:

∂ t ρ ϵ g + ∇ • U ϵ ρ ϵ g = 0 ∂ t (ρ ϵ g Ω ϵ ) + ∇ • ρ ϵ g V ϵ ⊗ Ω ϵ + d 3 ∇ρ ϵ g -γ s ∆ ρ ϵ g Ω ϵ = ρ ϵ g ϵ (1 -|Ω ϵ | 2 )Ω ϵ , (34)
where

U ϵ = d 1 Ω ϵ - 1 ξ ∇ ρϵ f - µ ξ ∇ρ ϵ g , V ϵ = d 2 Ω ϵ - 1 ξ ∇ρ ϵ f - µ ξ ∇ρ ϵ g .
In the limit ϵ → 0, one can shows that (34) converges towards [START_REF] Carrillo | Asymptotic flocking dynamics for the kinetic cucker-smale model[END_REF]. The main idea is based on the fact that the right-hand side of ( 34) is parallel to Ω ϵ , and the proof is similar to [START_REF] Motsch | Numerical simulations of a nonconservative hyperbolic system with geometric constraints describing swarming behavior[END_REF].

The numerical method is adapted from [START_REF] Motsch | Numerical simulations of a nonconservative hyperbolic system with geometric constraints describing swarming behavior[END_REF], using the so-called splitting scheme in two steps (dropping the ϵ for clarity):

• Step 1: Solve the conservative part:

∂ t ρ g + ∇ • U ρ g = 0 ( 35
)
∂ t (ρ g Ω) + ∇ • ρ g V ⊗ Ω + d 3 ∇ρ g -γ s ∆ ρ g Ω = 0, (36)
• Solve the relaxation part:

∂ t ρ g = 0 ∂ t (ρ g Ω) = ρ g ϵ (1 -|Ω| 2 )Ω Writting Q =   ρ g ρ g Ω x ρ g Ω y 
 , and Q = ρ g * Φ * Φ, system (35) can be written as:

∂ t Q + ∂ x F (Q, Q x , ∆Q x ) + ∂ y G(Q, Q y , ∆Q y ) = 0, ( 37
)
where we have denoted by ∆Q x (resp. ∆Q y ) the derivative in x (resp. in y) of the Laplacian of Q, and the fluxes write

F (Q, Q x , ∆ Qx ) =    d 1 Q(2) -1 ξ γρ A η Q(1)∆ Qx * Φ * Φ -µ ξ Q(1)Q x (1) d 2 Q(2) 2 Q(1) -1 ξ γρ A η Q(2)∆ Qx * Φ * Φ -µ ξ Q(2)Q x (1) + d 3 Q(1) -γQ x (2) d 2 Q(2)Q(3) Q(1) -1 ξ γρ A η Q(3)∆ Qx * Φ * Φ -µ ξ Q(3)Q x (1) -γQ x (3),    and G(Q, Q y , ∆ Qy ) =    d 1 Q(3) -1 ξ γρ A η Q(1)∆ Qy * Φ * Φ -µ ξ Q(1)Q y (1) d 2 Q(2)Q(3) Q(1) -1 ξ γρ A η Q(2)∆ Qy * Φ * Φ -µ ξ Q(2)Q y (1) -γQ y (2) d 2 Q(3) 2 Q(1) -1 ξ γρ A η Q(3)∆ Qy * Φ * Φ -µ ξ Q(3)Q y (1) + d 3 Q(1) -γQ y (3).

  

The explicit time-discretization for Eq. ( 6) writes:

Q * i,j = Q n i,j - ∆t ∆x F n i+ 1 2 ,j -F n i-1 2 ,j - ∆t ∆y G n i,j+ 1 2 -G n i,j-1 2 ,
where

F i+ 1 2 ,j = F (Q i,j ) + F (Q i+1,j ) 2 - 1 2 P 2 ( ∂F ∂Q ( Qi,j , Qxi,j , ∆ ˜x Q i,j ) Q i+1,j -Q i,j , with Qi,j = Q i,j + Q i+1,j 2 , Q xi,j = Q i+1,j -Q i,j 2 , Qxi,j = Q xi,j + Q xi+1,j 2 ,

  of the friction constant ζ. First we fix κ = 1000 and τ = 0.15, and show in Fig. 4 the values of α ∥ max and α ⊥ max (left panel) and of S th 1 and S th 2 (right panel), as functions of the bifurcation parameter b p and for different values of the agent friction constant ζ: ζ = 0.1 (blue curves), ζ = 0.5 (red curves), ζ = 1 (yellow curves). One can first observe in Fig.

  panel (C)), (b) the zone of parameters showing travelling bands increases when increasing the agent friction ζ (compare bottom to top rows of panel (C)). Therefore, we obtain a very good qualitative agreement between the simulations of the macro model and the tendencies predicted by the stability analysis as function of ζ.

Figure 2 :

 2 Figure 2: Simulations of the discrete model for the parameters indicated in table 1. Agents are represented as black arrows giving their direction of motion, obstacles are represented as blue circles. Panel (A): for weak obstacle stiffness κ = 10, panel (B): for mild obstacle stiffness κ = 100, panel (C): for large obstacle stiffness κ = 1000. In each panel, the vertical axis represents different values of the friction coefficient ζ (from bottom to top: ζ = 0.2, 0.5, 1, 2 for panels (A) and (B) and ζ = 0.2, 1, 2, 5 for panel (C); and the horizontal axis represents different values of the agent-agent repulsion µ: panel (A): µ ∈ {0.002, 0.02, 0.05, 0.1, 0.4, 0.6} , panel (B): µ ∈ {0.0002, 0.002, 0.005, 0.01, 0.04, 0.06}, panel (C): µ ∈ {2.10 -5 , 2.10 -4 , 6.10 -4 , 2.10 -3 , 4.10 -3 , 6.10 -3 }.

Figure 4 :

 4 Figure 4: Prediction of the linear stability analysis. Left: values of the maximal growth rate of the plane wave perturbations in the direction of Ω 0 (continuous lines) and in the orthogonal direction Ω ⊥ 0 (dashed lines) as functions of the bifurcation parameter b p , for different values of the agent friction ζ: ζ = 0.1 (blue curves), ζ = 0.5 (red curves), ζ = 1 (yellow curves). Right: same representation for the size of the perturbations in the two directions S th 1 and S th 2 .

Figure 5 :

 5 Figure 5: Left: values of the maximal growth rate of the plane wave perturbations in the direction of Ω 0 (continuous lines) and in the orthogonal direction Ω ⊥ 0 (dashed lines) as functions of the bifurcation parameter b p , for different values of the obstacle spring stiffness κ: κ = 10 (blue curves), κ = 100 (red curves), κ = 1000 (yellow curves). Right: same representation for the size of the perturbations in the two directions S th 1 and S th 2 .

Figure 6 :

 6 Figure 6: Left: values of the maximal growth rate of the plane wave perturbations in the direction of Ω 0 (continuous lines) and in the orthogonal direction Ω ⊥ 0 (dashed lines) as functions of the bifurcation parameter b p , for different values of the agent-obstacle distance τ: τ = 0.15 (blue curves), τ = 0.2 (red curves), τ = 0.3 (yellow curves). Right: same representation for the size of the perturbations in the two directions S th 1 and S th 2 .

Figure 7 :Figure 8 :

 78 Figure 7: Values of the maximal eigenmodes of the Fourier transform of the numerical solution (dotted curves) and predicted by the stability analysis (plain curves) in direction Ω 0 = (cos π/4, sin π/4) (left figures) and Ω ⊥ 0 (right coloured frames). Three different values of the obstacle spring stiffness are considered: κ = 10 (panel (A)), κ = 100 (panel (B)) and κ = 1000 (panel (C)), and three different agent-obstacle repulsion force distances τ = 0.15 (blue curves), τ = 0.2 (orange curves) and τ = 0.3 (yellow curves). Right column: examples of simulations for parameters reported on the graphs: simulations with diamond symbol match panel (B) (yellow frame: clusters; red and blue frame: trails) and the simulation with the circle symbol matches panel (C) (travelling bands).

Figure 9 :

 9 Figure 9: Example of the signature of a distribution with 4 bins, M = 85 and N x = 13. The different colors represent the different compartments of the signature.

Figure 10 :

 10 Figure 10: Histograms for simulations of Fig. 8 as defined in Eq. (28).

Figure 11 :Figure 12 :

 1112 Figure 11: Error between the density profile of the continuum simulation (high density clustered simulation of Fig. 8 left column) and its approximations using the procedure described in Sec. 5.1.1, as a function of the number of grid points for the PIC method N P IC (horizontal axis) and different number N of discrete particles (see insert for correspondance between curve color and N ), using the EMD distance based on histograms (left panel) or the L2 distance based on point values (right panel)

Figure 13 :

 13 Figure 13: EMD between the approximated continuum density and the discrete one as a function of b p for κ = 100, ζ = 0.5 and ϵ = 0.1, and for different values of the number of agents N : N = 500 (blue curve), N = 1000 (red curve) and N = 3000 (yellow curve) and N = 5000 (purple curve).

Figure 14 :

 14 Figure 14: Comparison between the discrete and continuum simulations for strong obstacle spring stiffness κ = 1000 and agent friction ζ = 0.5. Left figures: for weak agent agent repulsion µ = 2 10 -4 , right figures, for µ = 4 10 -3 . Top figures: EMD between the approximated continuum density and the discrete one as function of ϵ for different values of the number of agents N : N = 500 (blue curve), N = 1000 (red curve) and N = 3000 (yellow curve) and N = 5000 (purple curve). Bottom tables: simulations of the continuum model (left column), and of the discrete one for different values of ϵ: ϵ = 0.05 (second column), ϵ = 0.1 (third column), ϵ = 0.5 (fourth column) ϵ = 0.8 (fifth column) and ϵ = 1 (last column). The different rows correspond to different number of agents for the discrete simulations as well as for the discretization of the continuum density (from top to bottom: N = 500, N = 1000, N = 3000, N = 5000)

•

  run Swimmer demo.m : file which allows to fix the parameter values (set to reproduce Fig 3 of the paper), create folders for registering the datafiles generated by the model, and run the simulations • SwimmerSOH.m : MATLAB function that allows, given a set of parameters provided as entry, to run an entire simulation and register regularly the files in the specified data folder • The other files (get constants.m, get eigenvalues F.m, get eigenvalues G.m, get F.m, get G.m, get fluxes.m, get Jacobian F.m, get Jacobian G.m, mmat.m, convol.m, make plot.m) are intrinsic functions of the model, the description of which is contained in the corresponding files.

Table 1 :

 1 2.2.1). Note that only the ratio Parameters used for the discrete simulations of Fig.2. The various values considered for µ, ζ, κ are specified in the caption of Fig.2.

	Parameters	Value	Description
	N		3000	number of obstacles
	M		3000	number of agents
	u 0		1	agent speed
	r R		0.075	agent-agent repulsion distance
	r A		0.1	agent-agent alignment distance
	ν		2	agent-agent alignment intensity
	τ		0.15	agent-obstacle repulsion distance
	C ϕ		5	agent-obstacle repulsion intensity
	d s		0.02	noise in the agents' orientation
	η		1	obstacle friction
	d 0		0	obstacle positional noise
	µ	various	agent-agent repulsion intensity
	ζ	various	friction constant of the agents
	κ	various	spring constant coefficient
	Parameters	Value	Description
	h	≈ 6.7 • 10 -3	step-size spatial discretization
	v 0	1		agent speed
	r A	0.15	agent-agent alignment distance
	ds ν	0.01	parameter coming from alignment forces
	τ	0.15	agent-obstacle repulsion distance
	C ϕ	5		agent-obstacle repulsion intensity
	η	1		obstacle friction
	γ s	28 • 10 -4	viscosity coefficient
	µ	various	agent-agent repulsion intensity
	ζ	various	agent friction constant
	γ	various	γ = η/κ

Table 2 :

 2

Parameters used for the simulations of the continuum equations

[START_REF] Carrillo | Asymptotic flocking dynamics for the kinetic cucker-smale model[END_REF]

shown in Fig.

3

. The constants d 1 , d 2 , d 3 only depend on ν/d s and are obtained by computing the expressions

[START_REF] Cucker | On the mathematics of emergence[END_REF]

and

[START_REF] Degond | Asymptotic-preserving particle-in-cell methods for the vlasov-maxwell system near quasi-neutrality[END_REF]

.

Acknowledgements

The authors wish to thank Prof. E. Keaveny, Department of Mathematics, Imperial College London, for stimulating discussions. PD holds a visiting professor association with the Department of Mathematics, Imperial College London, UK. Part of this research was done when PD was affliated with this Department and supported by the Engineering and Physical Sciences Research Council (EPSRC) under grant no. EP/P013651/1. SMA was supported by the Vienna Science and Technology Fund (WWTF) with a Vienna Research Groups for Young Investigators, grant VRG17-014. The research by SMA was partially supported by the Austrian Science Fund (FWF) through the project F65.

A Supplementary material: Videos of IBM simulations

In the following paragraphs, we give some details on the videos of the IBM simulations available online as supplementary material. Each video is composed of two simulations representative of the different types of patterns shown in Fig. 2 and highlighted by a red cross. In each movie, agents are represented by black arrows and obstacles by colored points. The colors indicate the distance of the obstacles to their tethered points (from blue (close to their attachment site) to red (stretched springs)). Otherwise stated, the parameters used for the simulations are the ones indicated in table 1.

S1 -Trails

Link: https://doi.org/10.6084/m9.figshare.19615599.v2

In this video, two simulations are shown : The left movie is obtained for weak obstacle spring stiffness κ = 10, agent friction ζ = 2 and agent-agent repulsion µ = 0.02 and the right movie is obtained for mild obstacle spring stiffness κ = 100, agent friction ζ = 1 and agent-agent repulsion µ = 0.04. Both simulations show the spontaneous formation of trails of agents in a more or less deformable field of obstacle. For weak obstacle spring stiffness (left movie), agents easily repulse the obstacles as they move, creating large trails empty of obstacles that can merge or evolve over time. For larger obstacle spring stiffness (right movie), the agents also organize in trails that push the obstacles, creating tunnels with strong walls, more robust over time. It is noteworthy that trail structures are only observed for weak or mild obstacle spring stiffness.

S2 -Travelling bands -strong obstacles

Link: https://doi.org/10.6084/m9.figshare.19615884.v3 Here, both simulations feature strong obstacle spring stiffness κ = 1000 and agent friction ζ = 2. The left movie is obtained for weak agent-agent repulsion µ = 2.10 -5 while the right movie is obtained for µ = 10 -3 . In both situations, agents end up organizing in travelling bands on the long run, but we can observe a first phase when agents try to organize in trails. This suggest that the trail-like formation is only stable when agents have enough strength to push the obstacles as they move. Comparing the left and right movie, we also observe that larger agent-agent repulsion leads to larger clusters of agents.

S3 -Honneycomb structures Link: https://doi.org/10.6084/m9.figshare.19615116.v2

We consider here the case of weak obstacle spring stiffness κ = 10. The left movie is obtained for agent friction ζ = 0.2 and mild agent-agent repulsion µ = 0.05 while the right movie is obtained for ζ = 1 and larger agent-agent repulsion µ = 0.1. As one can observe, when the agent-agent repulsion is large enough in an easily deformable obstacle field, the agent phase wins over the obstacle phase, creating regularly spaced islands of obstacles in the form of honneycomb structures.

S4 -Moving Clusters

Link: https://doi.org/10.6084/m9.figshare.19615878.v2 Here, the left movie features weak obstacles κ = 10, agent friction ζ = 0.2 and agent-agent repulsion µ = 0.002, while the right movie is for κ = 100, ζ = 1 and µ = 0.004. When agentagent repulsion is low enough, agents spontaneously organize into more or less round clusters surrounded by obstacles. This cluster formation happens very fast, and agent clusters then move more or less fast depending on their environment. Large and slow clusters of agents are observed in the left movie, where obstacles are very loose and agent-agent repulsion is a bit larger, while more numerous, smaller and more stable agent clusters are observed in the right movie (where agent-agent repulsion is a bit larger but obstacles are stronger).

S5 -Pinned Clusters Link: https://doi.org/10.6084/m9.figshare.19615866.v1

In these movies, we consider a small agent friction ζ = 0.2 and agent-agent repulsion µ = 2.10 -4 , for mild obstacle spring stiffness κ = 100 (left pannel) and for strong obstacles κ = 1000 (right pannel). We classify these clusters as 'pinned' as the agents organize very fast into small and highly concentrated clusters that do not move (right pannel) or move very slow in the case of mild obstacles. We can also observe merging of clusters in the case of mild obstacles (left video). ∆ Qi,j ≈ Qi-1,j + Qi+1,j -2 Qi,j ∆x 2 + Qi,j-1 + Qi,j+1 -2 Qi,j ∆y 2

and where P 2 ( ∂F ∂Q ) is a second polynomial of a matrix at the intermediate state of ( Qi,j , Qxi,j , ∆ ˜x Q i,j ) and ( Qi+1,j , Qxi+1,j , ∆ Q xi+1,j ) computed following [START_REF] Degond | Polynomial upwind schemes for hyperbolic systems[END_REF]. Terms in G are computed the same way.

Convoluted terms are computed using a double fast-Fourier transform as:

where we used the fast Fourier transform.