P Degond

A Diez

A Walczak

Topological states and continuum model for swarmalators without force reciprocity

Keywords: individual-based model, macroscopic model, Fokker-Planck equation, BGK operator, self-organized hydrodynamics, synchronization, travelling-wave solution, index of a vector field, numerical simulations AMS Subject classification: 35F510, 35Q70, 35Q92, 37N25, 70F10, 82B40, 82C40

Swarmalators are systems of agents which are both self-propelled particles and oscillators. Each particle is endowed with a phase which modulates its interaction force with the other particles. In return, relative positions modulate phase synchronization between interacting particles. In the present model, there is no force reciprocity: when a particle attracts another one, the latter repels the former. This results in a pursuit behavior. In this paper, we derive a hydrodynamic model of this swarmalator system and show that it has explicit doubly-periodic travelling-wave solutions in two space dimensions. These special solutions enjoy non-trivial topology quantified by the index of the phase vector along a period in either dimension. Stability of these solutions is studied by investigating the conditions for hyperbolicity of the model. Numerical solutions of both the particle and hydrodynamic models are shown. They confirm the consistency of the hydrodynamic model with the particle one for small times or large phase-noise but also reveal the emergence of intriguing patterns in the case of small phase-noise.

Introduction

We investigate a new collective dynamics model based on the swarmalator concept described below. Collective dynamics refers to the coherent behavior of self-propelled agents subject to mutual interactions such as attraction, repulsion and alignment. Examples in nature are provided by birds [START_REF] Lukeman | Inferring individual rules from collective behavior[END_REF], fish [START_REF] Gautrais | Deciphering interactions in moving animal groups[END_REF][START_REF] Hemelrijk | Emergence of oblong school shape: models and empirical data of fish[END_REF], ants [START_REF] Boissard | Trail formation based on directed pheromone deposition[END_REF][START_REF] Couzin | Self-organized lane formation and optimized traffic flow in army ants[END_REF], bacteria [START_REF] Czirók | Formation of complex bacterial colonies via self-generated vortices[END_REF], sperm [START_REF] Creppy | Symmetry-breaking phase transitions in highly concentrated semen[END_REF], colloidal rollers [START_REF] Bricard | Emergent vortices in populations of colloidal rollers[END_REF]. How coherence emerges is still the subject of an intense literature in which mathematical models play a central role. In many models, self-propulsion is accounted for by imposing the agents a constant speed, identical for all agents, such as in the Vicsek model [START_REF] Vicsek | Novel type of phase transition in a system of self-driven particles[END_REF] and its many variants [START_REF] Aoki | A simulation study on the schooling mechanism in fish[END_REF][START_REF] Calovi | Swarming, schooling, milling: phase diagram of a data-driven fish school model[END_REF][START_REF] Cavagna | Flocking and turning: a new model for selforganized collective motion[END_REF][START_REF] Chaté | Collective motion of self-propelled particles interacting without cohesion[END_REF][START_REF] Costanzo | Spontaneous emergence of milling (vortex state) in a Vicsek-like model[END_REF][START_REF] Couzin | Collective memory and spatial sorting in animal groups[END_REF][START_REF] Hemelrijk | Schools of fish and flocks of birds: their shape and internal structure by self-organization[END_REF][START_REF] Hildenbrandt | Self-organized aerial displays of thousands of starlings: a model[END_REF] (see also the review [START_REF] Vicsek | Collective motion[END_REF]). By contrast, a constant speed is not assumed in the Cucker-Smale model [START_REF] Cucker | Emergent behavior in flocks[END_REF] and its variants [START_REF] D'orsogna | Self-propelled particles with soft-core interactions: patterns, stability, and collapse[END_REF][START_REF] Ha | From particle to kinetic and hydrodynamic descriptions of flocking[END_REF][START_REF] Ha | A simple proof of the Cucker-Smale flocking dynamics and mean-field limit[END_REF][START_REF] Motsch | A new model for self-organized dynamics and its flocking behavior[END_REF]. The synchronization of oscillator populations shares many similar features with collective dynamics. The paradigmatic model of oscillator synchronization is the Kuramoto model [START_REF] Kuramoto | Chemical Oscillations, Waves, and Turbulence[END_REF] which has stimulated an intense research activity (see e.g. the review [START_REF] Acebrón | The Kuramoto model: A simple paradigm for synchronization phenomena[END_REF]).

Among collective dynamics models, swarmalators have recently attracted increasing attention. Swarmalators are systems of agents that are simultaneously self-propelled particles and oscillators. In addition to position and velocity, they are endowed with a phase which may synchronize with the neighbors' phases. There is a two-way coupling between positions and phases: phase differences between neighboring particles modulate interaction forces and relative positions influence phase synchronization. The concept and terminology was introduced in [START_REF] O'keeffe | Oscillators that sync and swarm[END_REF] and further developed in [START_REF] Ha | Emergent behaviors of the swarmalator model for position-phase aggregation[END_REF][START_REF] Hong | Active phase wave in the system of swarmalators with attractive phase coupling[END_REF][START_REF] Hong | Coupling disorder in a population of swarmalators[END_REF][START_REF] Jiménez-Morales | Oscillatory behavior in a system of swarmalators with a shortrange repulsive interaction[END_REF][START_REF] Lee | Collective steady-state patterns of swarmalators with finite-cutoff interaction distance[END_REF][START_REF] Lizarraga | Synchronization and spatial patterns in forced swarmalators[END_REF][START_REF] O'keeffe | Collective behavior of swarmalators on a ring[END_REF][START_REF] O'keeffe | Ring states in swarmalator systems[END_REF] (see also the review [START_REF] O'keeffe | A review of swarmalators and their potential in bio-inspired computing[END_REF]). Swarmalator models have been applied to e.g. the collective swimming of nematode swarms [START_REF] Peshkov | Synchronized oscillations in swarms of nematode turbatrix aceti[END_REF] or the intercellular organization of multicellular organisms [START_REF] Japón | Intercellular communication and the organization of simple multicellular animals[END_REF]. Related earlier models associate phases with rotations of self-propulsion speeds [START_REF] Degond | Hydrodynamics of the Kuramoto-Vicsek model of rotating self-propelled particles[END_REF][START_REF] Levis | Activity induced synchronization: Mutual flocking and chiral self-sorting[END_REF][START_REF] Liebchen | Collective behavior of chiral active matter: Pattern formation and enhanced flocking[END_REF][START_REF] Paley | Oscillator models and collective motion[END_REF] but, in these models, there is only a one way coupling of phases on positions.

Topological states are solutions which have non-trivial topology quantified by some discrete topological index (such as a winding number). Topological states have appeared first in the quantum Hall effect and gradually in other applications such as topological insulators [START_REF] Hasan | Colloquium: topological insulators[END_REF][START_REF] Qi | Topological insulators and superconductors[END_REF]. Non-trivial topology endows the system with increased robustness against perturbations because breaking the topological state requires a jump of the topological index and thus, a finite amount of energy. This is the so-called 'topological protection' effect. In collective dynamics, topological states have been investigated in recent work [START_REF] Degond | Bulk topological states in a new collective dynamics model[END_REF][START_REF] Shankar | Topological sound and flocking on curved surfaces[END_REF][START_REF] Sone | Anomalous topological active matter[END_REF][START_REF] Souslov | Topological sound in active-liquid metamaterials[END_REF][START_REF] Zhang | Oscillatory chiral flows in confined active fluids with obstacles[END_REF] (see also the review [START_REF] Shankar | Topological active matter[END_REF]). Swarmalator models also support topological states, such as the "phase-wave states" of [START_REF] O'keeffe | Oscillators that sync and swarm[END_REF].

In the present paper, we propose a new swarmalator model and demonstrate that it possesses travelling-wave solutions having non-trivial topology. This model differs from the above cited swarmalator models because the phase-modulated force does not enjoy reciprocity: the force acted on a particle by another one is equal to the force acted by the latter on the former, and not the opposite as it should if reciprocity was enforced. The result is a pursuit behavior: when two particles interact, according to their relative phases one particle chases the other one or vice versa (see Fig. 1). Another difference with [START_REF] O'keeffe | Oscillators that sync and swarm[END_REF] is that our model is second order: the self-propulsion component of the velocity obeys a time-continuous version [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF] of the Vicsek model [START_REF] Vicsek | Novel type of phase transition in a system of self-driven particles[END_REF], namely, particles' self-propulsion velocities tend to align with their neighbors up to some noise. Lastly, the model includes noises in both self-propulsion velocity and phase, in contrast to most of the above cited literature.

Our methodology is based on studying a continuum version of the swarmalator model. In general, there are three stages of description of particle systems. The finest level of detail is provided by the particle system itself, which, for swarmalators, consists of a differential system for the positions, velocities and phases of all the particles involved. The next stage is given by the kinetic model where the system is described by the probability distribution of the particles in the space of positions, velocities and phases. The passage between particle and kinetic models requires letting N → ∞ where N is the number of particles. Its mathematical investigation has given rise to a large number of works, (see e.g. [START_REF] Cercignani | The Mathematical Theory of Dilute Gases[END_REF][START_REF] Degond | Macroscopic limits of the Boltzmann equation: a review[END_REF][START_REF] Spohn | Large Scale Dynamics of Interacting Particles[END_REF] in classical kinetic theory, [START_REF] Bolley | Mean-field limit for the stochastic Vicsek model[END_REF][START_REF] Briant | Cauchy theory and mean-field limit for general Vicsek models in collective dynamics[END_REF][START_REF] Diez | Propagation of chaos and moderate interaction for a piecewise deterministic system of geometrically enriched particles[END_REF] for the Vicsek model and [START_REF] Ha | A mean-field limit of the particle swarmalator model[END_REF] for the swarmalators). Kinetic models of collective dynamics have been proposed in e.g. [START_REF] Bertin | Hydrodynamic equations for self-propelled particles: microscopic derivation and stability analysis[END_REF][START_REF] Peruani | A mean-field theory for self-propelled particles interacting by velocity alignment mechanisms[END_REF]. For kinetic models derived from the Vicsek model, existence of solutions [START_REF] Figalli | Global well-posedness of the spatially homogeneous Kolmogorov-Vicsek model as a gradient flow[END_REF][START_REF] Gamba | Global weak solutions for Kolmogorov-Vicsek type equations with orientational interactions[END_REF], phase transitions [START_REF] Degond | Macroscopic limits and phase transition in a system of self-propelled particles[END_REF][START_REF] Degond | Phase transitions, hysteresis, and hyperbolicity for self-organized alignment dynamics[END_REF][START_REF] Frouvelle | Dynamics in a kinetic model of oriented particles with phase transition[END_REF], numerical methods [START_REF] Gamba | Spectral method for a kinetic swarming model[END_REF][START_REF] Griette | Kinetic equations and self-organized band formations[END_REF] and relations to models with no velocity normalization [START_REF] Bostan | Asymptotic fixed-speed reduced dynamics for kinetic equations in swarming[END_REF] have been studied.

The final stage which gives rise to the coarsest level of details but provides the highest tractability consists of continuum models. These models are obtained from kinetic ones through an asymptotic limit involving a small parameter representing the ratio of the microscopic scale (e.g. the range of particle interactions) and the macroscopic one (typically the size of the observed region). This approach originates from the kinetic theory of gases (see reviews in [START_REF] Cercignani | The Mathematical Theory of Dilute Gases[END_REF][START_REF] Degond | Macroscopic limits of the Boltzmann equation: a review[END_REF]). The continuum version of the Viscek model was derived via this asymptotic procedure for the first time in [START_REF] Degond | Macroscopic models of collective motion with repulsion[END_REF] (see also [START_REF] Frouvelle | A continuum model for alignment of self-propelled particles with anisotropy and density-dependent parameters[END_REF]) and gave rise to the new "Self-Organized Hydrodynamics (SOH)" model. This derivation provides closed formulas for the parameters of the hydrodynamic model as functions of those of the particle and kinetic model. Later, [START_REF] Jiang | Hydrodynamic limits of the kinetic selforganized models[END_REF] showed the mathematically rigorous validity of the asymptotic limit. Another, more intuitive approach [START_REF] Bertin | Boltzmann and hydrodynamic description for self-propelled particles[END_REF][START_REF] Bertin | Hydrodynamic equations for self-propelled particles: microscopic derivation and stability analysis[END_REF][START_REF] Toner | Flocks, herds, and schools: A quantitative theory of flocking[END_REF] was developed earlier and led to a different model, the "Toner-Tu" model. As interesting as it can be, the Toner-Tu model has no mathematically proven connection with the Vicsek particle model, and there are no rigorous formulas relating the parameters of the two models. There have been several extensions of [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF]. In relation to the present work, let us mention [START_REF] Degond | Macroscopic models of collective motion with repulsion[END_REF] which includes attraction-repulsion forces, and [START_REF] Degond | Body-attitude coordination in arbitrary dimension[END_REF][START_REF] Degond | A new flocking model through body attitude coordination[END_REF][START_REF] Degond | Quaternions in collective dynamics[END_REF][START_REF] Degond | Alignment of selfpropelled rigid bodies: from particle systems to macroscopic equations[END_REF] where alignment of body attitudes (in-stead of mere self-propulsion velocity) is considered and shown to support topological states [START_REF] Degond | Bulk topological states in a new collective dynamics model[END_REF]. Local existence of solutions for the SOH model was proved in [START_REF] Degond | Hydrodynamic models of selforganized dynamics: derivation and existence theory[END_REF][START_REF] Zhang | A local existence of viscous self-organized hydrodynamic model[END_REF] and numerical simulations can be found in [START_REF] Degond | Macroscopic models of collective motion with repulsion[END_REF][START_REF] Dimarco | Self-alignment driven by jump processes : Macroscopic limit and numerical investigation[END_REF][START_REF] Motsch | Numerical simulations of a nonconservative hyperbolic system with geometric constraints describing swarming behavior[END_REF].

In this paper, we first derive the continuum version of the swarmalator model under consideration. The derivation strongly relies on [START_REF] Degond | Macroscopic models of collective motion with repulsion[END_REF][START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF][START_REF] Frouvelle | A continuum model for alignment of self-propelled particles with anisotropy and density-dependent parameters[END_REF]] and details will be given in an appendix. The resulting model, called "Swarmalator Hydrodynamics (SH)" consists of three evolution equations for the particle density, mean self-propulsion velocity and mean phase. We will also consider the case where the phase noise is small which simplifies the system and leads to the "Noiseless Swarmalator Hydrodynamics (NSH)". We first analyze the conditions under which the NSH system is hyperbolic. We then specialize the SH system to two spatial dimensions and show the existence of doubly-periodic travellingwave solutions. These solutions present a non-trivial topology as evidenced by the index of the phase vector field being non-zero. We will then present numerical experiments which have three objectives: (i) the validation of the SH and NSH models as macroscopic descriptions of the swarmalators model; (ii) the numerical verification of the hyperbolicity conditions and (iii) the exploration and comparisons of the patterns obtained by the particle and SH or NSH models in relation with their topology. In forthcoming papers [START_REF] Degond | Topological ring waves in a swarmalator model without force reciprocity[END_REF][START_REF] Degond | Topological travelling waves of a swarmalator model within a strip[END_REF], we will explore other classes of travelling-wave solutions at both the particle and hydrodynamic level, investigate whether they enjoy topological protection and decipher the mechanisms of topological phase transitions when they occur.

The main innovations of this work are: (i) the introduction of a new swarmalator model involving force non-reciprocity; (ii) the derivation and hyperbolicity analysis of a system of continuum equations for this model, named the SH system; (iii) the derivation of a class of topologically non-trivial doubly-periodic travelling-wave solutions to the SH system (iv) the numerical validation of the SH model against the particle model and of its hyperbolicity, and the exploration of the patterns appearing with both models in relation with their topology.

The organization of this paper is as follows. In Section 2, we introduce the swarmalator system and derive the associated kinetic equations. Some technical points are deferred to Appendix 6. Section 3 presents the derivation of the hydrodynamic model, the analysis of its hyperbolicity and the determination of a special class of explicit travelling-wave solutions. Proofs of these results are collected in Appendices 7, 8 and 9. Numerical experiments are presented in Section 4 for the particle model and its hydrodynamic limit. Details on the numerical methods are deferred to Appendix 10. The videos of the outcome of the simulations can be found in the supplementary material. The list and description of the supplementary videos can be found in Appendix 11. Finally perspectives are drawn in Section 5.

Particle and kinetic models

In a first subsection, we present the particle swarmalator model on which this study is based. In a second subsection, we provide a kinetic formulation of this model in the limit of a large number of particles.

The particle model

We consider N particles labeled k = 1, . . . , N , having position X k ∈ R n . We suppose these particles are self-propelled with constant and uniform self-propulsion speed c 0 and direction of self-propulsion v k ∈ S n-1 . Additionally we assume that each particle is endowed with a phase φ k ∈ R/(2πZ). Neighboring particles interact via alignment of their self-propulsion direction to a local average self-propulsion direction vk on the one hand and via an attractive-repulsive potential W depending on their positions and phases. The self-propulsion speed is also subject to Brownian noise. The phases of neighboring particle are subject to synchronization to a local average phase φk and some noise. Finally, the particles are subject to a confinement potential V . Specifically, the system for (X k , v k , φ k ) k=1,...,N is written

dX k dt = c 0 v k -γ∇ x W (X k , φ k ), (2.1
)

dv k = P v ⊥ k • ν vk -∇ x V (X k (t)) dt + √ 2D dB k t , (2.2
)

dφ k = -ν ′ sin(φ k -φk ) dt + √ 2D ′ dB ′ k t , (2.3)
the self propulsion speed c 0 , the alignment frequency ν, the synchronization frequency ν ′ , the velocity noise D and the phase noise D ′ being positive and given. The confinement potential V : R d → [0, ∞) with V → ∞ as |x| → ∞ is also given. The terms dB k t and dB ′ k t describe brownian noises on R n and R respectively. The symbol • in (2.2) means that the stochastic differential equation is interpreted in the Stratonovitch sense, a condition for the vector v k to remain on S n-1 [START_REF] Hsu | Stochastic Analysis on Manifolds[END_REF].

The attractive-repulsive potential W :

R d × R/(2πZ) × [0, ∞) → R takes the form W (x, φ, t) = 1 N N j=1 ω(|x -X j (t)|) sin(φ j (t) -φ), (2.4)
with a given sensing function ω: [0, ∞) → R such that x → ω(|x|) is smooth on R n and normalized, i.e. R n ω(|x|) dx = 1. The given constant γ can be either positive or negative and specifies the intensity of the attractive-repulsive force. The contribution of the j-th particle to the force -γ∇ x W (x, φ, t) depends on the phase difference φ j (t)φ. This contribution is in the direction of -∇ x [ω(x -X j (t))] if φ j (t) is slightly ahead of the phase φ and γ > 0 or if φ j (t) is slightly behind the phase φ and γ < 0. It is in the opposite direction in the converse cases. Thus, the attractive or repulsive character of this force depends on the relative phases and on the sign of γ.

In (2.2), for v ∈ S n-1 , P v ⊥ stands for the orthogonal projection onto {v} ⊥ and has expression P v ⊥ = Idv ⊗ v where Id is the n × n identity matrix and ⊗ stands for the tensor product. To define vk we first introduce the current

J k (x, t) = 1 N N j=1 ζ(|x -X j (t)|) v j (t), (2.5)
with again, a given sensing function

ζ: [0, ∞) → R. Then, we let vk = J k |J k | (X k (t), t), (2.6)
assuming that the denominator does not vanish. We proceed analogously to define φk . We first introduce

L k (x, t) = 1 N N j=1 η(|x -X j (t)|) e iφ j (t) , (2.7)
with a given sensing function η: [0, ∞) → R. Then, we define φk (t) ∈ R/(2πZ) by

e i φk (t) = L k |L k | (X k (t), t), (2.8)
again, assuming that the denominator does not vanish. Here again, dB

′ k t stands for independent brownian motions in R. System (2.1)-(2.
3) is an extension of the time-continuous version of the Vicsek model [START_REF] Vicsek | Novel type of phase transition in a system of self-driven particles[END_REF] proposed in [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF]. Indeed, let us temporarily assume that there is neither external potential (V = 0) nor phase noise (D ′ = 0) and that all particle phases are initially equal.

Then the attraction-repulsion potential vanishes (W = 0), the phases remain constant and the position and velocities follow the time-continuous version of the classical Vicsek model [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF]:

dX k dt = c 0 v k , (2.9
)

dv k = P v ⊥ k • ν vk dt + √ 2D dB k t . (2.10)
Eq. (2.9) describes motion in the direction of v k at constant speed c 0 as a consequence of self-propulsion. The first term at the right-hand side of (2.10) tends to align v k with the mean direction of the neighbors vk computed through (2.5), (2.6) at rate ν. The second term of (2.10) generates a brownian motion of v k on the sphere and models velocity noise with noise intensity D (see e.g. [START_REF] Degond | Hydrodynamics of the Kuramoto-Vicsek model of rotating self-propelled particles[END_REF][START_REF] Degond | Macroscopic models of collective motion with repulsion[END_REF][START_REF] Degond | A new flocking model through body attitude coordination[END_REF][START_REF] Degond | Quaternions in collective dynamics[END_REF][START_REF] Dimarco | Self-alignment driven by jump processes : Macroscopic limit and numerical investigation[END_REF][START_REF] Frouvelle | A continuum model for alignment of self-propelled particles with anisotropy and density-dependent parameters[END_REF] for references on this model and some of its variants). When added, the potential V biases the alignment direction in the direction of -∇ x V , as shown in (2.2). This takes into account for instance, external cues in the agents' navigation. With the assumption that V (x) → ∞ as x → ∞, the potential confines the particles in a bounded region of space. Now, if the phases are not constant and/or if phase noise is present (D ′ ̸ = 0), the attraction-repulsion potential W turns on. It may seem strange that W appears in (2.1) and not in (2.2) like V . We provide a justification of it in Appendix 6.1 through an overdamped limit. Now, W biases particle motion in the direction of -γ∇ x W . It endows the particles with a pursuit behavior illustrated in Fig. 1 in the case of a pair interaction. The dynamics of the phases follows (2.3) which is conceptually similar to (2.10) if phases φ k are associated with unit vectors e iφ k in the complex plane C. Indeed, (2.3) just states that this vector is subject to alignment with the mean phase vector of the neighbors e i φk at rate ν ′ and to noise with noise intensity D ′ .

This system describes agents that swarm (through velocity alignment) and which, at the same time, are oscillators subject to synchronization (through phase alignment). Moreover, swarming and synchronization are coupled through space as motion in space depends on velocity and phase while swarming and synchronisation depend on space through the target velocity and target phase. Thus, this system belongs to the class of "swarmalators", a term coined in [START_REF] O'keeffe | Oscillators that sync and swarm[END_REF]. Note that in the present system, the natural oscillator frequency is supposed to be zero. There is no conceptual objection to include a non-zero natural frequency but the patterns generated by the system and further explored could be different. This will be investigated in future work. Swarmalators have recently stimulated an intense research, see e.g. [START_REF] Ha | Emergent behaviors of the swarmalator model for position-phase aggregation[END_REF][START_REF] Hong | Active phase wave in the system of swarmalators with attractive phase coupling[END_REF][START_REF] Hong | Coupling disorder in a population of swarmalators[END_REF][START_REF] Jiménez-Morales | Oscillatory behavior in a system of swarmalators with a shortrange repulsive interaction[END_REF][START_REF] Lee | Collective steady-state patterns of swarmalators with finite-cutoff interaction distance[END_REF][START_REF] Lizarraga | Synchronization and spatial patterns in forced swarmalators[END_REF][START_REF] O'keeffe | Collective behavior of swarmalators on a ring[END_REF][START_REF] O'keeffe | Ring states in swarmalator systems[END_REF].

F 1 X 1 X 2 F 2 e iϕ2 e iϕ1
1 Figure 1: Pursuit dynamics between two particles located at X 1 and X 2 . The force

F 1 = -γ∇ x W (W 1 , φ 1 ) acting on X 1 is equal to the force F 2 = -γ∇ x W (W 2 , φ 2 ) acting on X 2 . Both are given by F 1 = F 2 = γ 2 ω ′ (|X 1 -X 2 |) X 2 -X 1 |X 2 -X 1 | sin(φ 2 -φ 1 )
and are depicted by black arrows. The phases of the two particles are symbolized by the unit vectors e iφ k , k = 1, 2 in the complex plane which are depicted in blue. The unit circle is drawn in red. The current situation which occurs for γω

′ (|X 1 -X 2 |) > 0 diplays the pursuit of X 2 by X 1 . The converse would occur if γω ′ (|X 1 -X 2 |) < 0.

The kinetic model

Our goal is now to derive a continuum version of the swarmalator model (2.1)-(2.3) in the form of hydrodynamic-type equations. To do so, it is convenient to introduce an intermediate description between the particle and hydrodynamic models, the so-called kinetic model.

We begin with introducing the empirical measure of the particles in the space R n × S n-1 × R/(2πZ) by

µ N (x, v, φ, t) = 1 N N j=1 δ (X j (t),v j (t),φ j (t)) (x, v, φ), (2.11)
where δ (X j (t),v j (t),φ j (t)) (x, v, φ) is the Dirac delta located at (X j (t), v j (t), φ j (t)). This is a random measure. In the limit N → ∞, under appropriate conditions which we will assume satisfied, µ N converges in the weak sense to a deterministic measure f (x, v, φ, t) dx dv dφ where f satisfies the following kinetic model posed on R n × S n-1 × R/(2πZ):

∂ t f + ∇ x • c 0 v -γ∇ x W f (x, φ) f + ∇ v • P v ⊥ ν vf -∇ x V (x) f +ν ′ ∂ φ sin( φf -φ)f = D∆ v f + D ′ ∂ 2 φ f, (2.12)
with W f given by

W f (x, φ, t) = R n ×S n-1 ×[0,2π] ω(|y -x|) sin(ψ -φ) f (y, w, ψ, t) dy dw dψ. (2.13)
To define vf and φf , we first define

J f (x, t) = R n ×S n-1 ×[0,2π]
ζ(|y -x|) w f (y, w, ψ, t) dy dw dψ, (2.14)

L f (x, t) = R n ×S n-1 ×[0,2π]
η(|y -x|) e iψ f (y, w, ψ, t) dy dw dψ.

(2.15)

Then, we let vf (x, t) = J f |J f | (x, t), e i φf (x, t) = L f |L f | (x, t). (2.16)
We note that the operators ∇ v , ∇ v •, ∆ v are respectively the spherical gradient, divergence and Laplace-Beltrami operators. Finally, by the fact that the empirical measure (2.11) is a probability measure, f satisfies the normalization condition

R n ×S n-1 ×[0,2π] f (x, v, φ, t) dx dv dφ = 1.
(2.17)

A hydrodynamic scaling detailed in Appendix 6.2, eventually leads to the following modified kinetic problem, which depends on the scaling parameter ε:

∂ t f ε + ∇ x • v -γ∇ x U f ε (x, φ) f ε -∇ v • P v ⊥ ∇ x V (x) f ε = 1 ε D ∇ v • -kP v ⊥ u f ε f ε + ∇ v f ε + D ′ ∂ φ -k ′ sin(α f ε -φ)f ε + ∂ φ f ε ,(2.18) U f (x, φ, t) = |ℓ f (x, t)| sin(α f -φ), (2.19)
u f = j f |j f | , j f (x, t) = S n-1 ×[0,2π] w f (x, w, ψ, t) dw dψ, (2.20)
e iα f = ℓ f |ℓ f | , ℓ f (x, t) = S n-1 ×[0,2π]
e iψ f (x, w, ψ, t) dw dψ,

with k = ν D and k ′ = ν ′ D ′ . The parameter ε ≪ 1 encodes the ratio of the microscopic scale, i.e. the typical distance or time over which particle response to interactions takes place, and the macroscopic scale, i.e. the typical size or duration of the experiment. The hydrodynamic model is obtained as the ε → 0 limit of this system. It describes how the macroscopic scale is influenced by the microscopic dynamics on average and is developed in the next section.

The hydrodynamic model

In a first subsection, we derive the hydrodynamic model by letting ε → 0 in the kinetic model (2.18)-(2.21). Details are given in Appendix 7. In a second subsection, we simplify the hydrodynamic model by assuming small phase noise, and we study the resulting model.

Derivation of the hydrodynamic model

We being with introducing the von-Mises distributions: for u ∈ S n-1 and α ∈ R/(2πZ), we let

M u (v) = 1 Z e kv•u , Z = S n-1 e kv•u dv, (3.1)
N α (φ) = 1 Z ′ e k ′ cos(φ-α) , Z ′ = [0,2π] e k ′ cos(φ-α) dφ. (3.2)
We note that Z does not depend on u but only on k and likewise, Z ′ does not depend on α but only on k ′ . Then, we have the formal theorem, whose proof can be found in Appendix 7.

Theorem 3.1 Suppose that there is a smooth solution f ε to the kinetic model (2.18) for all ε > 0 and that this solution converges smoothly as ε → 0 to a function f 0 . Then,

f 0 (x, v, φ, t) = ρ(x, t) M u(x,t) (v) N α(x,t) (φ), (3.3)
where ρ, u and α are functions from R n ×[0, ∞) to [0, ∞), S n-1 and R/(2πZ) respectively, which satisfy the following systems of partial differential equations:

∂ t ρ + ∇ x • ρ(c 1 u + bρ∇ x α) = 0, (3.4) ∂ t u + (c 2 u + bρ∇ x α) • ∇ x u + P u ⊥ ∇ x (Θ log ρ + κV ) = 0, (3.5)
ρ ∂ t α + (c 1 u + b ′ ρ∇ x α) • ∇ x α -Θ ′ ∇ x • ρ∇ x ρ = 0, (3.6)
where the coefficients c 1 , c 2 , b, b ′ , Θ, Θ ′ and κ are given in Appendix 7 and have the following properties:

c 1 : [0, ∞) → [0, 1), k → c 1 (k)
is an increasing function of k with c 1 (0) = 0 and lim k→∞ c 1 (k) = 1 [START_REF] Frouvelle | A continuum model for alignment of self-propelled particles with anisotropy and density-dependent parameters[END_REF].

c 2 : [0, ∞) → (-1, 1), k → c 2 (k) is a function of k such that c 2 (0) = 0 and lim k→∞ c 2 (k) = 1.
Furthermore, for small and large k, we have 0 < c 2 (k) < c 1 (k) [START_REF] Frouvelle | A continuum model for alignment of self-propelled particles with anisotropy and density-dependent parameters[END_REF]. Appendix 7). Θ ′ , b and b ′ are functions of k ′ and γ which have the opposite sign to γ (see Appendix 7).

Θ = Θ(k) = k -1 , κ = κ(k) = n-1 k + c 2 (k) (see
Remark 3.1 Numerical simulations in dimension n = 2 [71] suggest that c 2 : [0, ∞) → [0, 1), k → c 2 (k) is an increasing function of k and that 0 < c 2 (k) < c 1 (k), ∀k ∈ (0, ∞).
But a rigorous proof of these properties is still lacking. Remark 3.2 Because of (2.17), ρ satisfies the normalization condition

R n ρ(x, t) dx = 1. (3.7)
In the remainder of this paper, we will assume that k is small or large enough so that the property 0 < c 2 (k) < c 1 (k) is guaranteed. Formula (3.3) relies on the following lemma, whose proof can be found in [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF]. Defining

Q(f ) = D ∇ v • -kP v ⊥ u f f + ∇ v f + D ′ ∂ φ -k ′ sin(α f -φ)f + ∂ φ f , (3.8)
we have :

Lemma 3.2 (i) We can write Q(f ) = D ∇ v • M u f ∇ v f M u f + D ′ ∂ φ N α f ∂ φ f N α f . (3.9)
(ii) Defining

D(f ) = S n-1 ×[0,2π] Q(f ) f M u f N α f dv dφ, (3.10)
we have

D(f ) = - S n-1 ×[0,2π] M u f N α f D ∇ v f M u f N α f 2 + D ′ ∂ φ f M u f N α f 2 dv dφ ≤ 0. (3.11) (iii) Let f = f (v, φ).
Then, the following assertions are equivalent:

(a)

Q(f ) = 0, (b) D(f ) = 0 (c) ∃(ρ, u, α) ∈ [0, ∞) × S n-1 × R/(2πZ) such that f (v, φ) = ρ M u (v) N α (φ). (3.12)
We note that Eq. (2.18) can be written

T (f ε ) = 1 ε Q(f ε ), (3.13)
with

T (f ) = ∂ t f + ∇ x • v -γ∇ x U f (x, φ) f -∇ v • P v ⊥ ∇ x V (x) f , (3.14)
being the transport operator. Lemma 3.2 and Formula (3.13) suggest that we could modify the collision operator (3.9) and introduce the relaxation (BGK-type) operator

Q R (f ) = D ρ f M u f N α f -f , (3.15)
with u f and α f given by (2.20), (2.21) and ρ f by

ρ f (x, t) = S n-1 ×[0,2π] f (x, w, ψ, t) dw dψ, (3.16)
and consider the analogous perturbation problem to (3.13), namely

T (f ε ) = 1 ε Q R (f ε ). (3.17)
It can be easily shown that Lemma 3.2 (ii) and (iii) still holds with (3.11) replaced by

D(f ) = -D S n-1 ×[0,2π] M u f N α f ρ f - f M u f N α f 2 dv dφ ≤ 0. (3.18)
Then, the following theorem holds Theorem 3.3 Theorem 3.1 still holds for Eq. (3.17), but with different expressions of the constants involved. Their expressions can be found in Appendix 7.

Kinetic Eq. (3.17) can be given an interpretation in terms of an interacting particle system following a jump process, aka a Piecewise Deterministic Markov Process (PDMP). This interpretation is given in Appendix 6. [START_REF] Bertin | Boltzmann and hydrodynamic description for self-propelled particles[END_REF].

We now make comments on the hydrodynamic model (3.4)- (3.6). Again, let us temporarily assume that there is no external potential (V = 0), no phase noise (Θ ′ = 0, see Appendix 8.1) and that the phases are initially constant (α| t=0 independent of x). Then, α remains constant in time (and independent of x) and the equations for ρ and u reduce to

∂ t ρ + ∇ x • (c 1 ρu) = 0, (3.19)
∂ t u + c 2 (u • ∇ x )u + Θ ρ P u ⊥ ∇ x ρ = 0, (3.20) |u| = 1, (3.21)
where we have highlighted in (3.21) the fact that u (the average self-propulsion direction) is a normalized vector. This model is the continuum version of the Vicsek model derived in [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF], referred to as the Self-Organized Hydrodynamics (SOH). Eq. (3.19) is the mass conservation (or continuity) equation for the fluid density ρ. It shows that the fluid velocity c 1 u establishes along the average self-propulsion direction u and has norm c 1 . Hence, the fluid speed is less than the particle speed (whose value is 1 after scaling) and this is because the direction of particle velocities is spread around u according to the von Mises distribution (3.1). Eq. (3.20) is akin to the momentum conservation equation in the isothermal Euler equation of gas dynamics. Indeed, the first two terms correspond to the material derivative of u and are balanced by a pressure force -Θ∇ x ρ. However, we can spot several differences. The first one is the normalization condition (3.21) which has no counterpart in the Euler equation and which gives rise to the projection operator P u ⊥ in factor of the pressure term. Indeed, this projection is needed to ensure consistency with the constraint (3.21). Another difference is that the material derivative does not involve the fluid velocity c 1 u but a different velocity c 2 u since c 2 ̸ = c 1 . This feature makes the model non Galilean-invariant, but this is no surprise because the particle model itself is not Galilean-invariant: there is a preferred frame where the particle speed is 1. The quantity c 2 is the speed at which information propagates among the agents for them to update their velocity in response to density gradients. We noted that, fairly generally, we have c 2 < c 1 . The fact that c 2 ̸ = c 1 will be key to the existence of topological travelling-wave solutions discussed in the next section.

Compared with the SOH model, the full model involves the additional equation (3.6) for the average phase α(x, t). Naturally, only spatial gradients of α influence the dynamics of ρ and u (as only differences of phases influence the particle dynamics) and their influence is the same in Eqs. (3.4) and (3.5): it adds the same term bρ∇ x α to both the fluid velocity c 1 u and the information velocity c 2 u. This reflects the influence of the attractionrepulsion potential which, at the particle level, adds to the self-propulsion velocity a vector whose orientation depends on the differences (i.e. gradients) of the phases and on the sign of γ (which controls that of b as stated in the theorem). The phase equation (3.6) is a balance equation similar to the velocity equation (3.5). The material derivative of α combines c 1 u and the contribution of ρ∇ x α as in (3.4) but intriguingly, this contribution is weighted by a different coefficient b ′ compared to Eqs. (3.4), (3.5). Finally the material derivative is balanced by a diffusive term of the density ρ. However, note that Θ ′ may have either signs as well as γ. Finally, the external potential appears in the velocity equation (3.5) and biases the mean alignment direction u (or its material derivative) consistently with what it does in the particle model. The full model (3.4)-(3.6) will be referred to as "Swarmalator Hydrodynamics" (SH).

Small noise limit in the phase equation

In this section, we neglect the noise in the phase equation, i.e., we make k ′ → ∞. It has been numerically observed that the most interesting patterns are obtained in this regime. Indeed, phase gradients increase coherence between the particles through the pursuit mechanism (see Fig. 1). On the opposite, phase noise, which contributes to equalizing the phases, destroys this coherence.

We have the following lemma, the proof of which can be found in Appendix 8.1.

Lemma 3.4 In the limit k ′ → ∞, we have b → -γ, b ′ → -γ, Θ ′ → 0. (3.22)
In this limit, the macroscopic system becomes (assuming ρ ̸ = 0):

∂ t ρ + ∇ x • ρ(c 1 u + bρ∇ x α) = 0, (3.23) ∂ t u + (c 2 u + bρ∇ x α) • ∇ x u + P u ⊥ ∇ x (Θ log ρ + κV ) = 0, (3.24) ∂ t α + (c 1 u + bρ∇ x α) • ∇ x α = 0, (3.25)
We note that the contributions of ρ∇ x α are now weighted by the same coefficient b in all three material derivatives (there was a different coefficient b ′ for the phase equation in the general SH (3.4) -(3.6)). Consequently, in (3.25), the phase is transported by the fluid velocity as defined from (3.23), i.e. c 1 u + bρ∇ x α.

In the remainder of this paper, we will focus on this system, further referred to as the "Noiseless Swarmalator Hydrodynamics" (NSH).

Taking the gradient of (3.4) and introducing z = ∇ x α, we find that (ρ, u, z) satisfies the following system:

∂ t ρ + ∇ x • ρ(c 1 u + bρz) = 0, (3.26) ∂ t u + (c 2 u + bρz) • ∇ x u + P u ⊥ ∇ x (Θ log ρ + κV ) = 0, (3.27) ∂ t z + ∇ x (c 1 u + bρz) • z = 0, (3.28) ∇ x ∧ z = 0, (3.29)
where ∇ x ∧ z is the exterior derivative of z if z is identified to a differential form, i.e. ∇ x ∧ z is the antisymmetric matrix with entries (

∇ x ∧ z) ij = ∂ x i z j -∂ x j z i . Eq. (3.29)
is a structural constraint that says that z is a gradient, namely that there exists α such that z = ∇ x α. Eq. (ii) Conversely, let (ρ, u, z) be a solution of System (3.26)- (3.29). Suppose z 0 = z| t=0 satisfies ∇ x ∧ z 0 = 0. Let α 0 be a solution of ∇ x α 0 = z 0 and let α be a solution of

∂ t α + (c 1 u + bρz) • z = 0, α| t=0 = α 0 . Then, (ρ, u, α) is a solution of the NSH System (3.23)-(3.25).
We investigate the hyperbolicity of System (3.26)-(3.29) in the case V = 0. Consider a uniform steady state of density ρ 0 , velocity u 0 and phase gradient z 0 (note that for this study, we ignore the normalization condition (3.7)). In view of Lemma 3.5, such a steady-state corresponds to a solution (ρ 0 , u 0 , α 0 ) of the NSH System (3.23)-(3.25) with

α 0 = α 0 (x, t) = z 0 • x -(c 1 u 0 + bρ 0 z 0 ) • z 0 t. ( 3
δ ∈ [-π, π), θ ∈ [0, π], ϕ ∈ [-π, π) such that u 0 = e 1 , z 0 = |z 0 | (cos δ e 1 + sin δ e 2 ), τ = sin θ cos ϕ e 1 + sin θ sin ϕ e 2 + cos θ e 3 . (3.31)
In other words, δ is the polar angle of z 0 in the frame (e 1 , e 2 ) and (θ, ϕ) are the spherical angles of τ in (e 1 , e 2 , e 3 ). This geometric setting is illustrated in Figure 2. Now, we have Lemma 3.6 In this Lemma, we make

V = 0. (i) If z 0 = 0 or z 0 ∥ u 0 and if c 1 -c 2 ̸ = 2|b|ρ 0 |z 0 |, System (3.26)-(3.29) is hyperbolic about (ρ 0 , u 0 , z 0 ).
λ 1 = (c 2 + bρ 0 |z 0 |) sin θ cos ϕ,
is of multiplicity at least n -2. The other three eigenvalues are simple and are the real roots of a cubic equation. In some special cases, one of these roots can coincide with λ 1 , in which case it has multiplicity n -1. In the general case, λ 1 has multiplicity equal to n -2.

If |z 0 | = 0, w.l.o.g. we can choose θ = π/2.
Then, the eigenvalues are (generically): * λ 1 = c 2 cos ϕ of multiplicity n -2 corresponding to propagation of the components of the velocity normal to both u 0 and ξ, * λ 2 = c 1 cos ϕ of multiplicity 1 corresponding to propagation of phase perturbations, * a pair of simple eigenvalues

λ ± = 1 2 (c 1 + c 2 ) cos ϕ ± (c 1 -c 2 ) 2 cos 2 ϕ + 4c 1 Θ sin 2 ϕ ,
corresponding to the intertwining of the density perturbation and the component of the velocity perturbation along

P u ⊥ 0 ξ. (ii) If z 0 ̸ = 0 and z 0 ̸ ∥ u 0 , there are two constants C 1 and C 2 with 0 < C 1 < C 2 , depending on c 1 , c 2 and Θ, such that for all ρ 0 |b||z 0 | ∈ (0, C 1 ) ∪ (C 2 , ∞), there exists values of the angle δ between u 0 and z 0 such that System (3.26)-(3.29) is not hyperbolic about (ρ 0 , u 0 , z 0 ).
The proof of this Lemma can be found in Appendix 8.2. (ii) The condition that c 1 -c 2 ̸ = 2|b|ρ 0 |z 0 | is purely technical. In this case, there is a double real eigenvalue which is not explicit. Thus, the dimension of the corresponding eigenspace is not known and we cannot guarantee hyperbolicity. However, this is an isolated value of |b|ρ 0 |z 0 | and so, we may state that generically, the problem is hyperbolic when z 0 ∥ u 0 .

(iii) As soon as z 0 ̸ ∥ u 0 and ρ 0 |b||z 0 | ie either small or large, there are configurations (ρ 0 , u 0 , z 0 ) that lead to hyperbolicity loss. The proof in Appendix 8.2 shows that this loss happens when the angle between u 0 and z 0 is close to π/2. So, there is no contradiction with the hyperbolicity result when z 0 ∥ u 0 . Although we do not know if hyperbolicity is lost at intermediate values of ρ 0 |b||z 0 |, we may expect to observe instabilities as soon as the angle between u 0 and z 0 is sizeable.

In the next section, we will consider the NSH model and derive several classes of explicit solutions in dimension n = 2.

Doubly periodic travelling-wave solutions in dimension n = 2

In this section, we restrict ourselves to dimension n = 2. We let (x 1 , x 2 ) be the cartesian coordinates of a point x ∈ R 2 and (e 1 , e 2 ) be the cartesian coordinate basis. We denote by (u 1 , u 2 ) the two coordinates of the self-propulsion velocity u in this basis. We recall that u is a normalized vector, i.e.

u 2 1 + u 2 2 = 1. (3.32)
We assume a spatial domain Ω = (0, 1) 2 with periodic boundary conditions. Then, we have the following

ρ = 1, (3.33) u = U, (3.34) α = 2π (px 1 + mx 2 ) -λt + α 0 , (3.35) λ = 2πc 1 (pU 1 + mU 2 ) + 4π 2 b(p 2 + m 2 ), (3.36)
where U = U 1 e 1 + U 2 e 2 and α 0 ∈ R is an arbitrary constant.

(ii) Stationary solutions with non-constant phase: let (p, m) ∈ Z 2 \ {(0, 0)} be such that

p 2 + m 2 ≤ c 1 2π|b| . (3.37)
Then there exist two vectors (if the inequality in (3.37) is strict) or a unique vector (if there is equality in (3.37)) U ∈ S 1 such that The proof of this theorem is given in Appendix 9. These periodic solutions can easily be generalized to arbitrary dimensions. They are also solutions to the SH model (3.4)

pU 1 + mU 2 = - 2πb c 1 (p 2 + m 2 ). ( 3
-(3.6) if in Prop. 3.7, we change b into b ′ .
We note that when x 1 or x 2 increases by one (i.e. one period), the phase α is increased by an integer multiple of 2π. So, as an element of R/(2πZ), α is Z-periodic, as ρ and u are. An alternate view consists of introducing the unit vector field e iα(x,t) where the plane is identified with C. In this representation, the integers p and m are the indices of this unit vector field when, at a given time t, one moves along the x 1 -and x 2 -axes by one period respectively (see Fig. 3). Thus, the pair (p, m) is a topological index of this solution and obviously, remains unchanged with time. The vectors e iα(x,0) at the points (x 1 , x 2 ) = ( k 8 , ℓ 8 ), for k, ℓ in {0, 1, . . . , 8} are shown. The color code corresponds to the angle α with red, pink, yellow, green, light blue, blue, magenta, purple corresponding to α = 0, π 4 , π 2 , 3π 4 , π, 5π 4 , 3π 2 , 7π 4 respectively. The indices of this vector field about the origin are 1 and 2 when one moves along the x 1and x 2 -axes respectively. These travelling-wave solutions are nothing but the restrictions to a periodic spatial domain of the uniform solutions of System (3.26)-(3.29), with the additional constraint that the phase itself must be periodic (hence the restriction to a discrete set of values of z 0 , namely z 0 ∈ 2πZ 2 \ {(0, 0)}). The stability of these solutions has been studied in Section 3.2. Thus, we may apply Lemma 3.6 and conclude that, as long as U = (p, m)/ p 2 + m 2 , these solutions are stable. Also, we may state that otherwise, if |b| p 2 + m 2 is either small or large, these solutions are unstable. In particular, the stationary solutions which are guaranteed to be stable are those for which we have equality in (3.37). In the unstable case, one may wonder what happens to the topological index (p, m). Since, it can only change by integer values, only two possibilities may occur: either the solution transitions to another solution with the same topological index, or, a phase discontinuity occurs somewhere. To distinguish between these two scenarios, numerical simulations are needed and will be developed in future work.

Numerical experiments

In this section, we present numerical simulations of both the particle model (2.1)-(2.3) and the hydrodynamic SH model (3.4)- (3.6). The goals and main conclusions of this section are summarized below.

In Section 4.1, we provide a validation of the SH model as a macroscopic description of the particle model by showing that the particle model follows the behavior predicted by the SH model. For this, we use the explicit doubly periodic solutions of Section 3.3 as a baseline and we provide a quantitative analysis of the convergence of the particle model to the SH model as the number of particles goes to infinity.

In Section 4.2 we study the long-time behavior of the particle system and the stability of the doubly-periodic travelling-wave solutions. We investigate the role of the noise in the phase equation (i.e. of the parameter k ′ ) and of the direction of the initial velocity u 0 . In particular, we demonstrate that the doubly-periodic solution is all the more stable than the level of phase-noise is high (i.e. than the value of k ′ is low). For low levels of noise in the phase equation (i.e for large values of k ′ ), the behavior of the particle system after a certain simulation time is characterized by the emergence of a strong segregation phenomenon between regions of constant phase separated by thin boundaries of very-low density. By contrast, for large levels of noise in the phase equation (i.e. for k ′ small), the doubly-periodic travelling wave solution is very stable starting from any configuration. Moreover, in this case, topologically non-trivial states can emerge even starting from an initially disordered state. Finally, these experiments also give a numerical verification of the hyperbolicity condition stated in Lemma 3.6. In particular we observe that when u 0 and ∇ x α 0 are normal and for low phase-noise levels, the destabilization of the corresponding travelling-wave solution is faster than when u 0 and ∇ x α 0 are aligned. In addition, we also observe that the solution is more stable when u 0 and ∇ x α 0 are negatively aligned (i.e. in opposite directions) than when they are positively aligned (i.e. in the same direction).

In Section 4.3, we present simulations of the SH model. The segregation behavior observed in the particle simulations is not clearly observed in the simulation of the SH system although we still observe the formation of thin regions of low-density for low levels of noise in the phase equation. Finally, the particle and SH simulations both show that the doubly-periodic travelling wave solution is more stable for large levels of noise in the phase equation and when the initial velocity is negatively aligned with the phase gradient.

Validation of the hydrodynamic limit

In order to numerically validate the derivation of the Swarmalator Hydrodynamics model (3.4)-(3.6), we simulate the particle system (2.1)-(2.3) with a set of parameters chosen accordingly to the scaling conditions presented in Appendix 6.2 and under which the mean-field and hydrodynamic limits are taken. Then we confront statistical quantities measured at the particle level with their prediction given by the hydrodynamic model.

Scaling

Following the notations of Appendix 6.2, we first fix the macroscopic space and time scales x 0 = 1 and t 0 = 1. Then we choose a large number of particles N and a small radius of interaction R. In a spatially homogeneous setting, each particle interact in average with N neigh = πR 2 N neighboring particles. As a rule of thumb, the mean-field regime is attained when a particle typically interacts with at least a few tens of other particles. We will therefore choose the parameters N and R such that N neigh ∼ 10 2 .

Then, as explained in Appendix 6.2, in order to take the hydrodynamic limit, the radius of interaction R is the only scaling parameter from which all the other parameters of the particle simulations can be defined. Namely, we choose a linear potential for the phase attraction-repulsion force

ω(r) = x 4 0 t 0 R 2 ω r R ,
where ω(r) = 3 π (1r)1 r≤1 and r = r/R. Note that with x 0 = 1 and t 0 = 1 then it holds that R 2 ω(|x|) dx = 1 and for y ∈ R 2 ,

∇ x ω(|y -x|) = C y -x |y -x| 1 |y-x|≤R , (4.1)
with C =

x 4 0 t 0 R 3 . Moreover, we choose ζ(r) = η(r) = x 2 0 R 2 1 r≤R .
The other parameters are also chosen depending on R as

D = x 0 t 0 R D, ν = x 0 t 0 R ν, D ′ = x 0 t 0 R D′ , ν ′ = x 0 t 0 R ν′ ,
so that D, ν, D′ , ν′ are dimensionless parameters. Note that the parameter γ is already a dimensionless parameter.

Setting of the experiment

In order to check the behavior of the particle system as N → +∞ and R → 0, we use as a test case the doubly periodic travelling wave solution derived in Section 3.3. Given arbitrary alignment parameters D, ν, D′ , ν′ , we initialize the particle system by drawing N particles uniformly on the torus (0, 1) 2 and the initial velocities with respect to the von Mises distribution M u 0 with u 0 = (0, -1) T . The phase of a particle at position (x 1 , x 2 ) is sampled from the von Mises distribution N 2πx 2 (corresponding to p = 0 and m = 1 in Proposition 3.7). Then, we choose the intensity of the phase attraction-repulsion force such that b

′ = 1 2π (1 + c 1 ), i.e. γ = - c 1 + 1 2πc ′ 1 ( 1 k ′ + c ′ 2 )
, so that the theoretical travelling wave speed is equal to

λ = -2πc 1 + 4π 2 b ′ = 2π,
where we recall that b ′ is given by (7.35). We recall that with the potential (4.1), a particle chases the particles which are slightly ahead of phase when γ < 0 and flees them when γ > 0. In the present case, with γ < 0, p = 0 and m = 1, it means that the phase attraction-repulsion force is opposite to the mean self-propulsion velocity of the particles (see e.g. the term inside the divergence operator in (3.4)) and λ = 2π means that the travelling-wave is moving upward and is 1-periodic. At the particle level, in order to observe a travelling wave and measure its speed, we need to define a suitable statistical indicator. First, we define the "mass" of a particle with phase φ as m(φ) = 1 + cos φ.

Then we consider the x 2 -coordinate of the "center of mass" (in the torus) at time t of the particle system, given by the formula:

x N 2 (t) = 1 2π arg 1 N N k=1 m(φ k )e 2iπX k,2 (t) , (4.2)
where X k,2 is the x 2 -coordinate of the position of particle k. The formal mean-field limit when N → +∞ leads to

1 N N k=1 m(φ k )e 2iπX k,2 (t) -→ N →+∞ 1 0 2π 0 (1 + cos φ) e 2iπx 2 N α(t,x 2 ) (φ) dx 2 dφ,
where we recall that N α denotes the von Mises distribution (3.2) with the parameter α given by the macroscopic model (3.35). Namely in this experiment α(t, x 2 ) = 2πx 2 -λt+α 0 where we chose α 0 = π. A direct computation shows that

x N 2 (t) → 1 2 + t mod 1. (4.3)

Results

We run the particle simulation with the parameters described above and we measure the quantity (4.2). First, Figure 4 illustrates the behavior of the particle simulation when N goes to infinity with a fixed small interaction radius. We measure the quantity x N 2 (t) for various values of N and we compare it to the theoretical prediction (4.3) when N → +∞. For small values of N , in a regime where the particles have only a few neighbors to interact with (N neigh ∼ 1), only a noisy behavior is observed although the center of mass is clearly biased to move in the expected direction. As N grows and since R is fixed, each particle interacts in average with a larger number of other particles. When the number of neighbours crosses ∼ 100, the system can be considered as being in a mean-field regime and we observe a very good qualitative agreement with the macroscopic prediction (4.3). Having shown a qualitative agreement between the measured and predicted behaviors, we then perform a quantitative analysis of the convergence of the particle scheme when N → +∞. Since the quantity x N 2 (t) qualitatively and theoretically converges as N → +∞ towards a straight line, we compute the slope of this line (using a standard linear regression method) and compare it to the predicted travelling wave speed, theoretically equal to 1. The results are shown in Figure 5.

In Figure 5a, we observe an excellent quantitative agreement between the measured and predicted travelling wave speeds when N is large. Namely, when N > 3 • 10 6 , the absolute value of the difference between the two quantities is, in average, of order 10 -3 with a standard deviation of the order 10 -2 . Moreover, Figure 5a gives an indication on the behavior of the particle scheme for small and moderate values of N . In particular, it should be noted that, when N is small, the measured speed is actually larger than the predicted one. To understand this phenomenon, let us first point out that with the chosen parameters, there are two competing effects: the phase gradient produces a positive force in the (0, 1) T direction (i.e. pointing upward) while the self-propulsion velocity is in the direction (0, -1) T (i.e. pointing downward). However there is a difference in the nature of these two forces. The self-propulsion velocity is subject to noise and different particles have independent noises. On the contrary, the force exerted by the phase gradient is just computed by taking an average over all (neighboring) particles and is thus less sensitive to the individual noises affecting the particles self-propulsion velocities and phases. Even when N is small, we can infer that the phase attractionrepulsion force has a smaller variance than the self-propulsion velocity. Since the norm of the self-propulsion velocity is constant equal to 1, a larger variance means that in average, the self-propulsion velocity opposite to the phase attraction-repulsion force has a norm smaller than 1. The phase attraction-repulsion force is thus winning over the two, which can explain the positive bias observed. Note also that when the self-propulsion velocity and the phase attraction-repulsion forces are positively aligned, we have observed (not shown here) that the measured speed is slower than the one predicted, which is expected since due to the noise, the two forces are not perfectly aligned and the norm of their sum is smaller than the theoretical value obtained when they are perfectly aligned.

A second important observation that we have made, but which is not directly shown in Figure 5, is that for larger values of R (R ∼ 0.01) and for very large values of N (N > 10 6 ), there is a perceptible negative bias in the measured travelling-wave speed, meaning that it is slower than expected (although it is still larger for N small). For R = 0.01 and for all values of N > 10 6 , the measured speed is about 0.97 instead of 1.

Since the theoretical speed is proportional to c 1 (k ′ ) this can be explained by the fact that when taking the hydrodynamic limit (see Section 7), this latter quantity appears as the limit when R → 0 of (y,t) dy under the assumption that the particle distribution is equal to the equilibrium distribution f (y, w, ψ, t) = ρ(y, t)M u(y,t) (w)N α(y,t) (ψ) with ρ(y, t) ≡ 1 and u, α are given by the macroscopic model. The kernel η R ≡ η R (|x|) := R -n η(|x|/R) is non negative, has integral one and tends to a Dirac delta at 0 when R → 0. Consequently, since η R has integral one, by the triangle inequality, the last term on the right-hand side is smaller than c 1 (k ′ ) for any R > 0. However, this effect is not perceptible for the value R = 0.005 as shown on Figure 5.

|L R | = R n ×[0,2π] η R (|y -x|)e iψ N α(y,t) (ψ) dy dψ = c 1 (k ′ ) R n η R (|y -x|)e iα
Finally, Figure 5b shows that the standard deviation around the average measured travelling-wave speed decreases as an inverse power law as N is increasing. The exponent of this inverse power law is computed using a standard linear regression and is equal to approximately 0.48. This experimental value has a simple theoretical explanation. In the mean-field limit, the particles can be shown to behave as independent and identically distributed random variables (thanks to the so-called propagation of chaos property). Consequently, by the central limit theorem, the standard deviation around the average (4.2) should behave like N -1/2 , which is consistent with the observations.

Segregation phenomena

In a macroscopic regime, the doubly periodic travelling wave observed in the particle simulations is stable during a certain time but may eventually transition towards another state. This destabilization phenomenon may be understood as a natural consequence of the numerical noise induced by the scheme (which combines both the finite size effect and the inherent stochasticity of the particles). It can also be seen as a practical way to study the stability of the different solutions of the macroscopic system (3.4)- (3.6). In this section we investigate how the particle system departs from the doubly periodic travelling wave configuration depending on the initial direction of the velocity and on the value of the parameter k ′ . This will give a numerical confirmation of the results of Lemma 3.6. In addition, we show that for low levels of noise in the phase equation (i.e. for large values of k ′ ), the particle system transitions towards configurations characterized by a strong segregation between populations of particles with constant phase separated by thin boundaries of low density.

Parameters and setting of the experiments

All the simulations in this section take place in an intermediate regime where the interaction radius is sizeable compared to the dimension of the domain. Consequently, this setting is slightly farther from the hydrodynamic limit (which requires R → 0) than that of Section 4.1 and we note a small but perceivable departure of the observed travellingwave speed from the value predicted by the hydrodynamic model. We choose this setting because simulations are quicker to run without qualitatively altering the results. All the simulations presented below use the following parameters:

N = 10 6 , R = 0.01, ν = 5, D = 1, γ = 0.2.
The particles' positions are sampled independently uniformly in the torus. The initial phase of a particle at position (x 1 , x 2 ) is sampled from the von Mises distribution N 2πx 2 . With this choice, the phase attraction-repulsion force drags the particles downward (i.e. in the direction (0, -1) T ). Given an initial velocity v 0 ∈ S 2 , the particles' velocities are sampled according to the von Mises distribution M v 0 .

In order to test the influence of the phase attraction-repulsion interaction, we vary the parameters ν′ and D′ and the choice of the initial velocity v 0 . First, we consider four choices for the parameters ν′ and D′ which correspond to four different levels of noise in the phase equation.

1. Very-low noise: k ′ ≡ k ′ hi+ = 10000 (ν ′ = 10, D′ = 0.001).

Low noise: k

′ ≡ k ′ hi = 200 (ν ′ = 10, D′ = 0.05). 3. Medium noise: k ′ ≡ k ′ med = 10 (ν ′ = 10, D′ = 1). 4. Large noise: k ′ ≡ k ′ lo = 3 (ν ′ = 3, D′ = 1).
For smaller values of k ′ , only a noisy behavior is observed as it can be expected. Then, for each value of k ′ , we will consider three different choices for the initial velocity v 0 .

1. When v 0 = (0, -1) T we say that the velocity and the phase attraction repulsion force are Positively Aligned (PA).

2. When v 0 = (0, 1) T we say that the velocity and the phase attraction repulsion force are Negatively Aligned (NA).

3. When v 0 = (1, 0) T we say that the velocity and the phase attraction repulsion force are Orthogonal (OT).

As a control system, for each value of k ′ , we also consider the system where the positions, velocities and phases of the particles are initially sampled independently uniformly respectively in (0, 1) 2 , S 1 and [0, 2π]. This configuration is referred as UF in the following.

Results

The results for all the simulations are shown in Videos 1 to 18 in Appendix 11.1. The main observations are summarized below.

1. In the very-low noise case (k ′ hi+ = 10000) and for any choice of the direction of the initial velocity v 0 , the particles immediately segregate into small regions of equal phase separated by very-low density thin boundaries.

Starting from a NA configuration, after a transition period during which many small constant-phase regions are forming, merging or expanding, the system stabilizes into a configuration where a succession of constant-phase regions with a band-like shape move in the direction v 0 = (0, 1) T (see Figs. 6a and6b and Video 1). The global direction of motion of the system v 0 remains constant throughout the simulation. Note that since there is no phase gradient within each region, the particles in each region are not subject to the phase attractionrepulsion force.

Starting from a PA configuration, the situation is initially analogous to the NA case except that the particles are moving in the direction (0, -1) T . However, as times grows, unlike the NA case, the system does not reach a stable configuration within the 40 units of time of the simulation. Although we still observe the formation of band-like regions of constant phase moving at a constant speed, these regions do not have a constant shape and are perpetually subject to destruction and recombination (see Video 2).

Starting from a OT configuration, unlike the NA and PA cases, the global velocity of the system does not remain constant. During the first units of time and simultaneously to the formation and recombination of constant-phase regions, the global velocity of the particles (and direction of motion) transitions from (1, 0) T to (0, 1) T . After 40 units of time, the systems reaches a configuration similar to the one starting from a NA configuration but no bands stable over a long time can be clearly identified (see Video 3).

In the low noise case (k ′

hi = 200), the dynamics is similar to the very-low noise case with a few exceptions. First, the final configuration starting from the NA and OT configurations is still composed of band-like constant-phase regions moving in the direction (0, 1) T but the size of the regions is increased (see Figs. 6c and6d and Videos 5 and 7) and the final state is more stable starting from a OT configuration. We also note that the phase in each band is not constant over time but slowly evolve. Secondly, starting from the PA configuration, the analogous configuration is much more unstable and the system finally ends up in a flocking phase with all the phases equal and an arbitrary direction of motion (see Video 6).

3. In the medium noise case (k ′ med = 10). For all choices of v 0 , the initial structure is preserved during a longer time (approximately 5 units of time). Then we observe the formation of thin elongated low-density regions. Unlike the previous cases, they do not clearly delimitate segregation regions and we de not observe the formation of constant-phase clusters. As time grows, the system finally ends up in a situation similar to the previous cases where band-like structures separated by thin lowdensity boundaries are moving along the x 2 -axis. In the OT case, the velocity transitions from (1, 0) T to (0, 1) T . We note that the final outcome of the system is still characterized by a segregation phenomenon into band-like structures but unlike the cases where k ′ is larger, these structures are larger and although they are clearly separated by thin low-density boundaries, there is an inner gradient of phase in each band and the phase is not preserved over time in each band as it was in the previous cases (see Figs. 6e and 6f and Videos 9 to 11). 4. In the large noise case (k ′ lo = 3) and for all choices of v 0 the initial doubly periodic travelling wave is stable throughout the simulation (see Figs. 6g and6h and Videos 13 to 15). Finally, starting from a uniformly disordered UF state with random velocities and phases, for very-low to medium levels of noise (k ′ = k ′ hi+ to k ′ = k ′ med ), the system always ends up in a flocking phase with all the phases and velocities equal. Similarly to what is observed starting from a NA, PA or OT configuration, we initially observe clusterization and phase-segregation phenomena which are all the more important that k ′ is large (see Videos 4, 8 and 12). The situation is more complex in the large noise scenario (k ′ = k ′ lo ) as independent experiments with the same parameters can lead to different outcomes. In addition to the flocking phase (see Video 16), we have also observed cases where the system ends up in various topologically non-trivial states. These states are still characterized by a continuous gradient of phase which is topologically constrained by the periodicity of the domain. However, unlike the doubly-periodic travelling wave solution, a wide range of complex structures can emerge characterized by non constant densities and velocities (see Videos 17 and 18).

Discussion

Regardless of the level of noise k ′ , all the experiments confirm that the NA configuration is the most stable one. Even though the initial doubly periodic travelling wave is not preserved, the only stable final configuration correspond to a case where the velocity is pointing in the direction (0, 1) T and the phase is constant along the x 1 -axis and piecewise constant and increasing along the x 2 axis. It can be understood as a more general version of the NA case where in the latter the phase is linear and increasing along the x 2 -axis. As predicted by Lemma 3.6, the OT configuration is the most unstable one and except for low values of k ′ , it never persists and quickly transitions towards a more stable NA case. Note that Lemma 3.6 does not make a distinction between the stability of the NA and PA cases. Numerically the former is the most stable.

An important observation is the ability of the particle model to produce segregation. We were not able to predict this behavior using the SH model (3.4)-(3.6), although we cannot exclude that it also corresponds to (possibly singular) solutions that remain to be identified. The simulation of the macroscopic model presented in the next section may also support this idea. It is also worth mentioning that the emergence of band-like structures is a well-known phenomenon in the Vicsek model [START_REF] Chaté | Collective motion of self-propelled particles interacting without cohesion[END_REF]. However, the phenomenon observed is quite different on many aspects. First, the emergence of bands in the Vicsek model can be observed only in a very specific range of parameters and in very-low density regimes. The conditions of the presented simulations are much more general. Moreover, the bands observed here have a very different profile from the ones observed in the Vicsek model. In the Vicsek model, bands have an asymmetric profile characterized by a sharp front edge and an exponentially decaying tail whereas the bands observed here are larger and have a symmetric profile (see Fig. 7). Regarding the doubly periodic travelling wave solution that we theoretically identified in Section 3.3, the particle simulations tend to indicate that the noise in the phase equation has a stabilizing effect. In the large noise case and for any choice of the initial velocity, we indeed do not observe any segregation but rather a stable doubly periodic travelling wave as expected, with a moving speed close to the theoretical speed predicted by the macroscopic model.

Finally, these experiments show the robustness of topological states. Even if for lower levels of noise in the phase equation, the theoretical doubly periodic travelling wave does not persist, the system still transitions towards a state characterized by a nontrivial topology. Such scenario is never observed starting from the topologically trivial UF state. The topological robustness increases with the level of noise in the phase equation and in such regime, more general topological states can even emerge form the UF configuration, which shows that the flocking state is not the only global attractor of the system. The doubly periodic travelling wave may be seen as the simplest topological state and further investigations are needed to determine whether the topological states observed at the particle level also correspond to solutions of the hydrodynamic model.

Simulation of the hydrodynamic NSH and SH models

Following the methodology introduced in [START_REF] Motsch | Numerical simulations of a nonconservative hyperbolic system with geometric constraints describing swarming behavior[END_REF] for the SOH model, a finite volume discretization of the SH model is also presented. The details of the numerical method are described in Appendix 10.2. The results for all the simulations are shown in Videos 19 to 21 in Appendix 11.2. The main observations are summarized below.

NSH case

Simulating the NSH system is computationally easier, in particular because it does not require to discretize the term ∇ x • (ρ∇ x ρ) in the phase equation. Starting from a doubly periodic initial condition and regardless of the direction of the initial velocity, we observe a very stable travelling wave moving in the direction of the phase gradient at the speed predicted by the model. However, starting from a slightly perturbed initial condition, we immediately observe the formation of shocks with strong local variations of the density and of the phase. Due to the CFL condition which becomes too stringent in this situation, we were not able to continue the simulation further. The same difficulty happens in the very-low noise case k ′ hi+ . Note however that this observation is consistent with the particle simulations which also show this behavior.

SH case

By adding the terms corresponding to a nonzero k ′ ̸ = 0 (i.e. by adding noise in the phase equation), the situation becomes more stable even starting from perturbed initial conditions. The behavior is not immediately comparable to the behavior of the particle simulations, but shares some of the main features. For the experiments discussed below, in each scenario, the initial state is perturbed by adding a small random uniform noise for the density, phase and velocity independently for each cell. For each cell, we add to the theoretical density (equal to 1) a uniform random variable in [-0.25, 0.25]. For the two angles which define the phase and velocity we add to the theoretical value a uniform random variable in the interval [-0.75, 0.75] (in radians).

In the low noise case k ′ hi = 200, starting from a perturbed OT state we observe the formation of thin low-density regions very reminiscent of the ones observed in the particle simulations and a transition towards a more stable state close to a NA configuration (see Video 21). Unlike the particle simulations, such behavior is not observed starting from either perturbed NA or PA states. In the PA case, even though the density does not stay uniform and we can observe the formation of thin regions of lower density, they quickly fade away and do not degenerate as in the OT case or in the particle simulations (see Video 20). In the NA case, the initially perturbed density, phase and velocity are quickly restored to their theoretical unperturbed state (see Video 19), which is another numerical confirmation of the increased stability of this state. In all cases, the simulation reaches a doubly periodic travelling wave solution with a travelling wave speed correctly predicted by the model (3.36).

In the medium noise and large noise cases k ′ med = 10 and k ′ lo = 3, regardless of the initial condition (OT, NA or PA), unlike the particle simulations, we do not observe the formation of inhomogeneities but rather a stable doubly periodic travelling wave solution (see Videos 22 to 27). This behavior is similar to the one observed in the particle simulations in the large noise scenario and confirm the stabilizing effect of the noise in the phase equation.

The simulations of the SH and particle models agree well for short times or when the phase-noise level is large enough. At long time scales when the phase-noise is small, the doubly-periodic solutions are more unstable with the particle simulations than with the SH model. In the particle simulations, we observe the formation of regions of constant phase separated by thin low density regions. By contrast, simulations of the hydrodynamic model always maintain continuous gradients in phase.

Conclusion and perspectives

In this paper, we have presented a new swarmalator model without force reciprocity and derived its hydrodynamic limit. We have studied the hydrodynamic model in the limit of small phase noise and determined its hyperbolicity regime. Then, we have derived a class of explicit doubly-periodic travelling-wave solutions in two spatial dimensions. These solutions have non-trivial topology quantified by the index of the phase vector over a period in either dimension. Solutions with index values larger than one are possible. Then numerical simulations of these doubly-periodic travelling-wave solutions with both the particle and hydrodynamic models have been presented. They confirm that the hydrodynamic model is an accurate approximation of the particle one for short time or large phase noise. They also provide a validation of the hyperbolicity result. However, for long times and small values of the phase noise, the two models differ but both give rise to topological solutions. In forthcoming papers [START_REF] Degond | Topological ring waves in a swarmalator model without force reciprocity[END_REF][START_REF] Degond | Topological travelling waves of a swarmalator model within a strip[END_REF], we will pursue the investigation of topological states in this swarmalator model by deriving and studying classes of travelling-wave solutions in other geometries such as strips or annulae.

Further studies can be envisioned. For instance, numerical simulations suggest that in some parameter ranges, periodic-in-time solutions are generated. Their mathematical investigation is still open. Another direction is to explore other phase spaces that would generate solutions with more complex topologies such as higher order homotopy groups of spheres. Appendices 6 Particle and kinetic models: additional material 6.1 Overdamped limit for potential W System (2.1)-(2.3) is the limit as ε → 0 of the following system

dX k dt = c 0 v k + w k , (6.1
)

dv k = 1 ε ε P v ⊥ k • ν vk -∇ x V (X k (t)) dt + √ 2D dB k t , (6.2
)

dw k dt = - 1 ε w k + γ∇ x W (X k , φ k ) , (6.3)
complemented with (2.3) (here, the parameter ε has a different meaning compared to the next section). This can be interpreted as follows. Forgetting ε for the time being, from (6.1), we see that the particle velocity is decomposed in two terms: one c 0 v k stemming from self-propulsion which has constant norm, and a complementary one w k . The total force acting on the particle is thus the sum of a component acting on the self-propulsion velocity given by (6.2) and of one acting on its complement given by ( 6.3). The first component is unchanged from the original system (see (2.2)) while the second one is a relaxation force. It describes the competition between the external force ∇ x W and a friction force which hypothetically results from the surrounding medium (note that a similar friction could also be included in (6.2) but would vanish anyway because P v ⊥ k v k = 0). Now, 1/ε in factor of (6.2), (6.3) is the friction coefficient and is very large. We see that, in the limit ε → 0, we recover (2.1)-(2.3) under the condition that the forces involved in (6.2) are of order ε (hence the multiplication by ε of all the terms involved). So, System (2.1)-(2.3) is obtained as the overdamped limit of the unrelaxed system under the assumption that the alignment force, the noise and the exterior potential are very small, of the same order as the inverse of the friction coefficient.

Scaling of the kinetic model

We first non-dimensionalize the kinetic model (2.12). We let x 0 , t 0 be space and time units which we relate to each other by x 0 = c 0 t 0 . we note that v, φ, v, φ are already dimensionless. We introduce the change of variables x = x/x 0 , t = t/t 0 and functions

f (x, v, φ, t) = x n 0 f (x 0 x, v, φ, t 0 t), Ṽ (x) = t 0 x -1 0 V (x 0 x), W (x, φ, t) = t 0 x -2 0 W (x 0
x, φ, t 0 t). We also assume that there exists R > 0 and functions ω, ζ and η such that, for all r ∈ [0, ∞):

ω(r) = x n+2 0 t 0 R n ω r R , ζ(r) = x n 0 R n ζ r R , η(r) = x n 0 R n η r R .
We define dimensionless constants

D = Dt 0 , D′ = D ′ t 0 , k = ν D , k ′ = ν ′ D ′ , R = R x 0 .
In these new variables, the kinetic model reads (after dropping the tildes and bars for simplicity):

∂ t f + ∇ x • v -γ∇ x W f (x, φ) f -∇ v • P v ⊥ ∇ x V (x) f = D ∇ v • -kP v ⊥ vf f + ∇ v f + D ′ ∂ φ -k ′ sin( φf -φ)f + ∂ φ f , (6.4) W f (x, φ, t) = R n ×S n-1 ×[0,2π] 1 R n ω |y -x| R
sin(ψφ) f (y, w, ψ, t) dy dw dψ, (6.5)

J f (x, t) = R n ×S n-1 ×[0,2π] 1 R n ζ |y -x| R w f (y, w, ψ, t) dy dw dψ. (6.6) L f (x, t) = R n ×S n-1 ×[0,2π] 1 R n η |y -x| R e iψ f (y, w, ψ, t) dy dw dψ, (6.7) vf (x, t) = J f |J f | (x, t), e i φf (x, t) = L f |L f | (x, t). (6.8)
We now make the following scaling assumptions:

R = ε → 0, D = O 1 ε , D ′ = O 1 ε , k = O(1), k ′ = O(1), γ = O(1).
Thus, introducting D and D′ such that D = D/ε and D ′ = D′ /ε, we may assume that D and D′ are constants. After this scaling, the problem is written (again dropping the tildes for simplicity):

∂ t f ε + ∇ x • v -γ∇ x W ε f ε (x, φ) f ε -∇ v • P v ⊥ ∇ x V (x) f ε = 1 ε D ∇ v • -kP v ⊥ vε f ε f ε + ∇ v f ε + D ′ ∂ φ -k ′ sin( φε f ε -φ)f ε + ∂ φ f ε , (6.9) W ε f (x, φ, t) = R n ×S n-1 ×[0,2π] 1 ε n ω |y -x| ε sin(ψ -φ) f (y, w, ψ, t) dy dw dψ, (6.10)
J ε f (x, t) = R n ×S n-1 ×[0,2π] 1 ε n ζ |y -x| ε w f (y, w, ψ, t) dy dw dψ. (6.11)
L ε f (x, t) = R n ×S n-1 ×[0,2π] 1 ε n η |y -x| ε e iψ f (y, w, ψ, t) dy dw dψ, (6.12)
vε f (x, t) = J ε f |J ε f | (x, t), e i φε f (x, t) = L ε f |L ε f | (x, t). (6.13)
Now, expanding expressions (6.10)-(6.13) in powers of ε, we get:

W ε f = U f + O(ε 2 ), U f (x, φ, t) = S n-1 ×[0,2π]
sin(ψφ) f (x, w, ψ, t) dw dψ, (6.14)

J ε f (x, t) = ζ 0 j f + O(ε 2 ), j f (x, t) = S n-1 ×[0,2π] w f (x, w, ψ, t) dw dψ, (6.15)
L ε f (x, t) = η 0 ℓ f + O(ε 2 ), ℓ f (x, t) = S n-1 ×[0,2π]
e iψ f (x, w, ψ, t) dw dψ, (

)

vε f = u f + O(ε 2 ), u f (x, t) = j f |j f | (x, t), (6.17)
φε f = α f + O(ε 2 ), e iα f (x, t) = ℓ f |ℓ f | (x, t), (6.18)
with

ζ 0 = R n ζ(|x|)
dx, and a similar definition for η 0 . Furthermore, expanding sin(ψφ) in (6.14), we note that U f is given by (2.19).

Finally, introducing expansions (6.14)-(6.18) into (6.9) and neglecting the resulting O(ε) terms which will have no influence on the final result, we are led to (2.18).

Particle system associated with the BGK operator

We use the same notations as Section 2.1. Each particle k ∈ {1, . . . , N } is associated with an increasing sequence of random numbers T 1 k , T k 2 , . . . , T k n , . . . which are subject to the condition that the interval between two consecutive numbers are independent random variables following a Poisson process with intensity D. At time T k n the (velocity, phase) pair of the k-th particle jumps between from

(v k , φ k )(T k n -0) to (v k , φ k )(T k n + 0), while X k is continuous (i.e. X k (T k n + 0) = X k (T k n -0)). For t ∈ (T k n , T k n+1
), the triple (X k , v k , φ k ) evolves according to the following differential system:

dX k dt = c 0 v k -γ∇ x W (X k , φ k ), dv k dt = -P v ⊥ k ∇ x V (X k (t)), dφ k dt = 0. with initial condition (X k , v k , φ k )(T k n +0) and (X k , v k , φ k )(T k n+1 -0)
are the values obtained by the solution of this system at time T k n+1 . Finally at jump time

T k n , the pair (v k , φ k )(T k n + 0) is drawn according to the von Mises distribution M vk (T k n -0) N φk (T k n -0)
where vk (T k n -0) and φk (T k n -0) are computed by (2.6) and (2.8) in which t is taken equal to T k n -0. This type of jump process is known as a Piecewise Deterministic Markov Process (PDMP).

In [START_REF] Diez | Propagation of chaos and moderate interaction for a piecewise deterministic system of geometrically enriched particles[END_REF], it is proved that, in the limit N → ∞, the empirical measure of this process (see Section 2.2) converges to the following kinetic equation:

∂ t f + ∇ x • c 0 v -γ∇ x W f (x, φ) f -∇ v • ν P v ⊥ vf f = D f -M vf N φf ,
with vf and φf given by (2.16). The scaling developed in Section 6.2 can be developed analogously here. In particular, it results in the localisation of vf and φf which are then replaced by u f and α f given by (2.20) and (2.21) respectively. After scaling and neglect of higher order terms in ε, the kinetic model reduces to (3.17) with the BGK type collision operator (3.15). BGK-type models of Vicsek-type dynamics have been investigated in [START_REF] Degond | Body-attitude coordination in arbitrary dimension[END_REF][START_REF] Degond | Phase transitions and macroscopic limits in a BGK model of body-attitude coordination[END_REF][START_REF] Degond | Bulk topological states in a new collective dynamics model[END_REF][START_REF] Degond | Alignment of selfpropelled rigid bodies: from particle systems to macroscopic equations[END_REF][START_REF] Dimarco | Self-alignment driven by jump processes : Macroscopic limit and numerical investigation[END_REF].

7 Limit ε → 0: proofs In this section, we prove Theorems 3.1 and 3.3. The proofs are identical for the two theorems. We develop it for Theorem 3.1 and only point out what is different for Theorem 3.3 when necessary. The proof follows a certain number of steps.

Step 1: f is given by (3.3) with the functions ρ, u and α to be determined. Indeed, if f ε → f 0 as ε → 0 smoothly, then f 0 satisfies Q(f 0 ) = 0, which in view of Lemma 3.2 (iii) means that, at any given point (x, t), f is given by (3.12). At a different point (x ′ , t ′ ), the equilibrium (3.12) may be different. This means that ρ, u and α are functions of (x, t), still to be determined, and that f is given by (3.3).

Step 2: derivation of the mass conservation equation (3.4). This is simply done by integrating (3.13) with respect to (v, φ) and using that for any smooth functions f , Q(f ) dv dφ = 0, (in this discussion, we omit the integration domain S n-1 × [0, 2π) any time the context is clear). This cancels the 1/ε singularity and leads to

T (f ε ) dv dφ = 0.
Letting ε → 0, we finally get

T (f 0 ) dv dφ = 0.
We note that the ∇ v term in the expression (3.14) of T cancels in the integration with respect to v. The time and space derivatives commute with the integrals in v and φ and we get

∂ t ρ + ∇ x • ρ v -γ∇ x U f 0 (x, φ) M u (v) N α (φ) dv dφ = 0. ( 7.1)
Now, we have from (2.19):

U f 0 (x, φ) = c ′ 1 ρ sin(α -φ), with c ′ 1 = 2π 0 cos φ e k ′ cos φ dφ 2π 0 e k ′ cos φ dφ . (7.2)
So, we get Step 3: Computation of the generalized collision invariants (GCI). To find equations for u and α is a not as straightforward, as there are no collision invariants, i.e. function χ(v, φ) such that for all smooth functions f Q(f ) χ dv dφ = 0, other than constant functions. In [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF], it was shown that this difficulty can be solved by the concept of generalized collision invariant (GCI). We summarize the approach here and refer to [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF][START_REF] Frouvelle | A continuum model for alignment of self-propelled particles with anisotropy and density-dependent parameters[END_REF] for details. To define the GCI concept, we first introduce the following operators if the collision operator is Q given by (3.8), then

∇ x U f 0 (x, φ) = c ′ 1 ∇ x ρ sin(α -φ) + ρ cos(α -φ)∇ x α . (7.3)
Q(f ; u, α) is defined for any (u, α) ∈ S n-1 × R/(2πZ) by Q(f ; u, α) = D ∇ v • M u ∇ v f M u + D ′ ∂ φ N α ∂ φ f N α . (7.6)
if the collision operator is Q R given by (3.15), then

Q R (f ; u, α) = D ρ f M u N α -f . (7.7)
We note that

Q(f ) = Q(f ; u f , α f ), (7.8)
and we have a similar relation between Q R and Q R . Then, we have the:

Definition 7.1 Given (u, α) ∈ S n-1 × R/(2πZ), the function χ u,α : S n-1 × R/(2πZ) → R, (v, φ) → χ u,α (v, φ) is a GCI for Q associated to (u, α)
if and only if the following holds:

Q(f ; u, α) χ u,α dv dφ = 0, ∀f such that P u ⊥ f v dv dφ = 0 and f sin(φ -α) dv dφ = 0. (7.9)
We have a similar definition of a GCI for Q R by replacing Q by Q R .

We have the Proposition 7.2 The set G u,α of GCI χ u,α is a vector space given by

G u,α = A h cos(φ-α) sin(φ-α)+h(u•v) P u ⊥ v •B +C | A, C ∈ R, B ∈ {u} ⊥ , (7.10)
where the functions h and h: (-1, 1) → R are given as follows:

Case of operator Q: h and h are such that the functions

g(φ) = h cos φ sin φ, g(θ) = h(cos θ) sin θ, (7.11)
are the solutions of the following equations:

-∂ φ e k ′ cos φ ∂ φ g = sin φ e k ′ cos φ , (7.12)
for g and

-∂ θ sin n-2 θ e k cos θ ∂ θ g + n -2 sin 2 θ sin n-2 θ e k cos θ g = sin θ sin n-2 θ e k cos θ (7.13)
for g; these solutions are unique in the spaces H 1 (0, π) for g and

g | sin n 2 -2 θ g ∈ L 2 (0, π), sin n 2 -1 θ g ∈ H 1 (0, π) ,
for g; the functions g and g are nonnegative on [0, π] and g can be extended into a smooth odd function on [-π, π],

Case of operator Q R : h and h are given by

h = 1, h = 1. (7.14)
Proof of Proposition 7.2. Case 1: collision operator Q. In [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF][START_REF] Frouvelle | A continuum model for alignment of self-propelled particles with anisotropy and density-dependent parameters[END_REF], it is shown that (7.9) is equivalent to saying that

∃A ∈ R, ∃B ∈ {u} ⊥ such that Q * (χ u,α ; u, α) = A sin(φ -α) + B • v, (7.15)
where

Q * (•; u, α) is the formal L 2 -adjoint of Q(•; u, α).
Computing this adjoint, we find that χ = χ u,α is a GCI if and only if it satisfies the following problem

D 1 M u ∇ v • M u ∇ v χ + D ′ 1 N α ∂ φ N α ∂ φ χ = A sin(φ -α) + B • v. (7.16)
For given A and B, by Lax-Milgram theorem, it can be shown that this problem has a unique solution in the subspace Ḣ1 of H 1 (S n-1 × R/(2πZ)) consisting of functions g satisfying g dv dφ = 0. Furthermore, any solution to (7.16) in H 1 is equal to this special solution up to an additive constant. Denoting by χ A,B the unique solution of (7.16) in Ḣ1 corresponding to the pair (A, B), we have χ A,B = χ A,0 + χ 0,B with χ A,0 only depending on φ and solution to

D ′ 1 N α ∂ φ N α ∂ φ χ A,0 = A sin(φ -α), (7.17)
and χ 0,B only depending on v and solution to

D 1 M u ∇ v • M u ∇ v χ = B • v. (7.18)
These two problems have been solved in [START_REF] Degond | Continuum limit of self-driven particles with orientation interaction[END_REF][START_REF] Frouvelle | A continuum model for alignment of self-propelled particles with anisotropy and density-dependent parameters[END_REF] (note that (7.17) is a special case of (7.18) corresponding to the dimension n = 2). Their solutions are given as follows:

χ A,0 = -A h cos(φ -α) sin(φ -α), χ 0,B = -h(u • v) B • v, (7.19)
with h and h given by (7.11). This leads to (7.10).

Case 2: collision operator Q R . Eq. (7.15) stands but now (7.16) is changed into

D χ(w, ψ) M u (w) N α (ψ) dw dψ -ψ(v, φ) = A sin(φ -α) + B • v. (7.20)
It is clear that any two solutions of this equation differ by a constant. So, we can single out a solution by requiring that χ(w, ψ) M u (w) N α (ψ) dw dψ = 0. It follows that G u,α is given by (7.10) with h = 1 and h = 1.

We will now write

χ α (φ) = h cos(φ -α) sin(φ -α), ⃗ χ u (v) = h(u • v) P u ⊥ v, so that G u,α = A χ α (φ) + ⃗ χ u (v) • B + C | A, C ∈ R, B ∈ {u} ⊥ . (7.21)
We remark that ⃗ χ u is a vector (perpendicular to u) further referred to as the vector GCI. The main use of the GCI is as follows. We first note that

P u ⊥ f f v dv dφ = P u ⊥ f (|j f |u f ) = 0, f sin(φ -α f ) dv dφ = |ℓ f | sin(α f -α f ) = 0.
Therefore, f satisfies the conditions of (7.9) for (u, α) = (u f , α f ). Thanks to (7.8), we deduce that

Q(f ) χ u f ,α f dv dφ = Q(f ; u f , α f ) χ u f ,α f dv dφ = 0.
From this, we obtain:

T (f ε )χ u f ε ,α f ε dv dφ = 1 ε Q(f ε )χ u f ε ,α f ε dv dφ = 0.
In particular, we get

T (f ε )⃗ χ u f ε dv dφ = 0, T (f ε )χ α f ε dv dφ = 0.
But it is clear that χ α is continuous with respect to α and ⃗ χ u with respect to u. We can then let ε → 0 and obtain

T (f 0 )⃗ χ u f 0 dv dφ = 0, T (f 0 )χ α f 0 dv dφ = 0. (7.22)
These are the two missing equations of the fluid model as shown in the next step.

Step 4: Explicit forms of the equations for u and α. In this step, we are making the two equations (7.22) explicit. We have

T (f 0 ) = (∂ t + v • ∇ x )(ρM u N α ) -γ∇ x • ∇ x U ρMuNα ρM u N α -∇ v • P v ⊥ ∇ x V (x) ρM u N α . (7.23)
Using the notations D = ∂ t + v • ∇ x and v ⊥ = P u ⊥ v as well as (7.3), computations give

D(ρM u N α ) = M u N α Dρ + kρv ⊥ • Du + k ′ ρ sin(φ -α)Dα , (7.24)
∇ x • ∇ x U ρMuNα ρM u N α = -c ′ 1 ρ∆ x ρ -ρ 2 |∇ x α| 2 + |∇ x ρ| 2 sin(φ -α) -3ρ∇ x ρ • ∇ x α + ρ 2 ∆ x α cos(φ -α) +kρv ⊥ • ∇ x ρ sin(φ -α) -ρ∇ x α cos(φ -α) • ∇ x u +k ′ ρ sin(φ -α) ∇ x ρ sin(φ -α) -ρ∇ x α cos(φ -α) • ∇ x α M u N α , (7.25)
-∇ v • P v ⊥ ∇ x V (ρM u N α ) = ρM u N α (n -1) ∇ x V • v -k P v ⊥ ∇ x V • P v ⊥ u . (7.26) Using the decomposition v = (v • u)u + v ⊥ , we can write T (ρ M u N α ) = M u N α T ee + T eo + T oe + T oo ,
where T ee is even with respect to both v ⊥ and sin(φα), T eo is even with respect to v ⊥ and odd with respect to sin(φα), T oe is odd with respect to v ⊥ and even with respect to sin(φα) and finally T oo is odd with respect to both arguments. We have

T eo = k ′ ρ sin(φ -α) ∂ t + (v • u) u • ∇ x α + γc ′ 1 ρ∆ x ρ -ρ 2 |∇ x α| 2 +|∇ x ρ| 2 sin(φ -α) -k ′ ρ 2 |∇ x α| 2 sin(φ -α) cos(φ -α) , (7.27)
T oe = v ⊥ • ∇ x ρ + kρv ⊥ • ∂ t u + (v • u) (u • ∇ x )u + (n -1)ρv ⊥ • ∇ x V +ρk(u • v)v ⊥ • ∇ x V -γc ′ 1 kρ 2 cos(φ -α) v ⊥ • (∇ x α • ∇ x )u , (7.28)
while the other terms will not be needed in the forthcoming computations. We now consider the first Eq. (7.22). Since χ u f 0 is odd with respect to v and even with respect to α, the only term in T (ρ M u N α ) which will not vanisy in the integration by imparity will be that corresponding to T oe . By Lemma 4.1 of [START_REF] Degond | Nematic alignment of self-propelled particles: From particle to macroscopic dynamics[END_REF], we have, for any function k(v • u),

S n-1 k(v • u) v ⊥ ⊗ v ⊥ dv = 1 n -1 S n-1 k(v • u) (1 -(v • u) 2 ) dv P u ⊥ .
Using this, the first Eq. (7.22) leads to (7.29) with

C 1 P u ⊥ ∇ x ρ + kC 1 ρ∂ t u + kC 2 ρ(u • ∇ x )u + (n -1)C 1 + kC 2 ρP u ⊥ ∇ x V -γkc ′ 1 2 C 1 ρ 2 (∇ x α • ∇ x )u = 0,
C i = 1 n -1 π 0 cos i-1 θ h(cos
Θ = 1 k , κ = n -1 k + c 2 . (7.32)
and b is given by (7.5). We proceed similarly for the second Eq. (7.22). In this case, χ α f 0 is even with respect to v and odd with respect to φ. So, the only term of T (ρ M u N α ) which remains is that corresponding to T eo . This leads to

k ′ C ′ 1 ρ∂ t α + c 1 C ′ 1 k ′ ρ(u • ∇ x )α + γc ′ 1 C ′ 1 (ρ∆ x ρ -ρ 2 |∇ x α| 2 +|∇ x ρ| 2 ) -γc ′ 1 k ′ C ′ 2 ρ 2 |∇ x α| 2 =
′ = -γc ′ 1 1 k ′ + c ′ 2 , (7.35)
c ′ 2 = C ′ 2 C ′ 1 = 2π 0 cos φ h(cos φ) e k ′ cos φ sin 2 φ dφ 2π 0 h(cos φ) e k ′ cos φ sin 2 φ dφ , (7.36)
Θ ′ = - γc ′ 1 k ′ . (7.37)
Finally, the properties of the coefficients listed at the end of the theorem statement are direct consequences of [START_REF] Frouvelle | A continuum model for alignment of self-propelled particles with anisotropy and density-dependent parameters[END_REF]. This ends the proof of Theorem 3.1.

8 Small noise limit in the phase equation 8.1 Expression of the system: proof of Lemma 3.4 We recall that b, b ′ and Θ ′ are given by (7.5), (7.35), (7.37) respectively. Now, it is proved in [START_REF] Frouvelle | A continuum model for alignment of self-propelled particles with anisotropy and density-dependent parameters[END_REF]

that c ′ 1 → 1, c ′ 2 → 1 as k ′ → ∞.
Thus, (3.22) follows immediately.

8.2 Hyperbolicity: proof of Lemma 3.6

Let (ρ 0 , u 0 , z 0 ) ∈ (0, ∞) × S n-1 × R n be given. Then, (ρ, u, z) = (ρ 0 , u 0 , z 0 ) is a spatially uniform stationary solution of System (3.26)- (3.29). The linearization of System (3.26)-(3.29) about this equilibrium solution is then given by

∂ t ρ + (c 1 u 0 + 2bρ 0 z 0 ) • ∇ x ρ + c 1 ρ 0 ∇ x • u + bρ 2 0 ∇ x • z = 0, (8.1) ∂ t u + Θρ -1 0 P u ⊥ 0 ∇ x ρ + (c 2 u 0 + bρ 0 z 0 ) • ∇ x u = 0, (8.2)
∂ t z + b|z 0 | 2 ∇ x ρ + ∇ x u z 0 + ∇ x z (c 1 u 0 + 2bρ 0 z 0 ) = 0, (8.3) ∇ x ∧ z = 0. (8.4) u • u 0 = 0. (8.5)
In (8.3), ∇ x u denotes the gradient matrix of u i.e. (∇ x u) ij = ∂ x i u j , for all i, j ∈ {1, . . . , n}.

The expression ∇ x u z 0 refers to the multiplication of the matrix ∇ x u and the vector z 0 . Similar definitions apply to ∇ x z. The constraint (8.5) expresses that the first order variation of a normalized vector is orthogonal to that vector. We take the partial Fourier transform of this system with respect to x and denote the resulting unknown by (ρ, û, ẑ)(ξ, t), with ξ the Fourier dual variable to x. We assume ξ ̸ = 0 (the case ξ = 0 corresponds to constants, which we already know are solutions of the linearized system). We recall that τ = ξ/|ξ|. From (8.4) we deduce that ẑ = zτ , where z ∈ R. Also, due to (8.5), we can project (8.2) on {u 0 } ⊥ , and the component u 1 of u on e 1 is always zero.

Using these remarks, the Fourier transform of System (8.1)-(8.4) leads to

1 i|ξ| ∂ t ρ + τ • (c 1 u 0 + 2bρ 0 z 0 )ρ + c 1 ρ 0 τ • û + bρ 2 0 z = 0, 1 i|ξ| ∂ t û + Θρ -1 0 P u ⊥ 0 τ ρ + (c 2 u 0 + bρ 0 z 0 ) • τ û = 0, 1 i|ξ| ∂ t z + b|z 0 | 2 ρ + û • z 0 + (c 1 u 0 + 2bρ 0 z 0 ) • τ ẑ = 0.
Choosing a reference frame such that u 0 , z 0 and τ are expressed by (3.31), and denoting by V = (ρ, û2 , û3 , û4 , . . . , ûn , z) T , we can write

1 i|ξ| ∂ t V + AV = 0,
where A is given in dimension n ≥ 3 by:

A =            X 1 c 1 ρ 0 sin θ sin ϕ c 1 ρ 0 cos θ 0 . . . 0 bρ 2 0 Θρ -1 0 sin θ sin ϕ X 2 0 0 . . . 0 0 Θρ -1 0 cos θ 0 X 2 0 . . . 0 0 0 0 0 X 2 .
0 | 2 |z 0 | sin δ 0 . . . . . . 0 X 1            , (8.6)
with

X 1 = sin θ c 1 cos ϕ + 2bρ 0 |z 0 | cos(ϕ -δ) , X 2 = sin θ c 2 cos ϕ + bρ 0 |z 0 | cos(ϕ -δ) .
In dimension n = 2, we can choose θ = π/2 and A is given by:

A =   X 1 c 1 ρ 0 sin ϕ bρ 2 0 Θρ -1 0 sin ϕ X 2 0 b|z 0 | 2 |z 0 | sin δ X 1   .
We recall that the system is hyperbolic about (ρ 0 , u 0 , z 0 ) if and only if A is diagonalizable with real eigenvalues for all values of θ and ϕ.

In the case n ≥ 3, developing the determinant with respect to the last column, we find:

det(A -λId) = (X 1 -λ)(X 2 -λ) n-3 X 1 -λ c 1 ρ 0 sin θ sin ϕ c 1 ρ 0 cos θ Θρ -1 0 sin θ sin ϕ X 2 -λ 0 Θρ -1 0 cos θ 0 X 2 -λ +(-1) n+2 bρ 2 0 n k=3 (-1) 2k-1 (X 2 -λ) Θρ -1 0 sin θ sin ϕ X 2 -λ b|z 0 | 2 |z 0 | sin δ = (X 2 -λ) n-2 (X 1 -λ) (X 1 -λ)(X 2 -λ) -L + M R -M 2 (X 2 -λ) =: (X 2 -λ) n-2 P (λ), (8.7)
with

L = c 1 Θ(cos 2 θ + sin 2 θ sin 2 ϕ), M = bρ 0 |z 0 |, R = Θ sin θ sin ϕ sin δ. (8.8)
After rearranging, we get

P (λ) = Aλ 3 + Bλ 2 + Cλ + D, (8.9)
with

A = -1, B = 2X 1 + X 2 , C = -2X 1 X 2 -X 2 1 + L + M 2 , (8.10) D = X 2 1 X 2 -LX 1 + M R -M 2 X 2 . (8.11)
In the case n = 2, a direct computation shows that (8.7) with P given by (8.9), (8.10) (8.11) is still true provided we make θ = π/2 in (8.8).

In the case n ≥ 3, we see that X 2 is an eigenvalue of A with multiplicity at least n -2. We show that the associated eigenspace has dimension at least n -2. First, by inspection of (8.6), it is clear that the space V = {(0, 0, 0, v 4 , . . . , v n , 0) | (v 4 , . . . , v n ) ∈ R n-3 } is a subspace of ker(A -X 2 Id). Now, we may remove the lines and columns of A -X 2 Id of indices comprised between 4 and n and check the following determinant

Y c 1 ρ 0 sin θ sin ϕ c 1 ρ 0 cos θ bρ 2 0 Θρ -1 0 sin θ sin ϕ 0 0 0 Θρ -1 0 cos θ 0 0 0 b|z 0 | 2 |z 0 | sin δ 0 Y , with Y = X 1 -X 2 = sin θ (c 1 -c 2 ) cos ϕ + M cos(ϕ -δ) .
Developing with respect to the last column, we easily realize that this determinant is equal to 0, showing that there is a one-dimensional complement space to V in ker(A -λId) and consequently, that this eigenspace is at least of dimension n -2.

According to classical results about the cubic equation [START_REF] Jacobson | Basic Algebra I[END_REF], P has three distinct real roots if and only if its discriminant ∆ is positive, i.e. ∆ =:

1 27 4(B 2 + 3C) 3 -(2B 3 + 9BC + 27D) 2 > 0.
Some algebra leads to

B 2 + 3C = Y 2 + 3(L + M 2 ), 2B 3 + 9BC + 27D = -2Y 3 -9Y (L -2M 2 ) + 27M R,
and finally we get ∆ as a polynomial of Y :

∆ = 4Y 4 M 2 + 4Y 3 M R + Y 2 (L 2 + 20LM 2 -8M 4 ) +18Y M R(L -2M 2 ) + 4(L + M 2 ) 3 -27M 2 R 2 .(8.12)
Proof of (i). We first consider the case where δ = 0 or δ = π. In this case (8.8) shows that R = 0 and ∆ reduces to

∆ = 4Y 4 M 2 + Y 2 (L 2 + 20LM 2 -8M 4 ) + 4(L + M 2 ) 3 = 4 M 2 (Y 2 -M 2 ) 2 + L 3 + 3L 2 M 2 + 3LM 4 + Y 2 (L 2 + 20LM 2 ).
Since L ≥ 0, all terms in the last expression are nonnegative so we get ∆ ≥ 0. If L ̸ = 0, we see that ∆ > 0 and so, the three roots of P are real and distinct. If L = 0, we have

∆ = 4M 2 (Y -M ) 2 (Y + M ) 2 , so that ∆ = 0 ⇐⇒ Y = ϵ 1 M,
with ϵ 1 = ±1. On the other hand, from (8.8), we have

L = 0 ⇐⇒ θ = π 2
(1 -ϵ 1 ϵ 2 ϵ 3 )M = ϵ 1 ϵ 2 (c 1 -c 2 ).
Since c 1c 2 > 0 (see Theorem 3.1), we must have (1ϵ 1 ϵ 2 ϵ 3 ) = 2. Hence ∆ can only be zero if 2|M | = c 1c 2 , which is ruled out in the assumptions of the theorem. Thus, it follows that P has three real and distinct roots. Now, we study whether one of these roots coincides with the root X 2 found earlier. Indeed, in this case, we have P (X 2 ) = 0. From (8.9), (8.10), (8.11) we readily get that P (X 2 ) = LY . So, P (X 2 ) = 0 implies L = 0 or Y = 0.

If L = 0, injecting (8.13) into (8.6) shows that the eigenspace associated to the eigenvalue X 2 contains the subspace

W = {(0, v 2 , . . . , v n , 0) | (v 2 , . . . , v n ) ∈ R n-1 }. Since Y 2 -M 2 =
0 is not allowed by the assumptions of the theorem, we see that A -X 2 Id has rank 2 which shows that ker(A -X 2 Id) has dimension exactly equal to n -1. Since the roots of P are distinct, X 2 is a simple root of P so that, as an eigenvalue of A, X 2 has multiplicity n -1. So, in this case, A is diagonalizable with real eigenvalues.

If Y = 0, injecting it into (8.6) readily shows that A-X 2 Id has rank 2. We conclude similarly as in the previous case. Now, in the general case where no root of P coincides with X 2 , the eigenspaces associated with the roots of P are one-dimensional. Then, X 2 is an eigenvalue of multiplicity exactly equal to n -2 and since the associated eigenspace has dimension at least n -2 as previously shown, its dimension is exactly n -2. It follows that A is diagonalizable with real eigenvalues.

This shows that the model is hyperbolic about (ρ 0 , u 0 , z 0 ) when z 0 ∥ u 0 . The case where z 0 = 0 is obvious and left to the reader.

Proof of (ii), case |M | large. We notice that the condition that ρ 0 |b||z 0 | large just means that |M | is large. If we introduce T and Z such that

Y = T + M Z, T = (c 1 -c 2 ) sin θ cos ϕ, Z = sin θ cos(ϕ -δ),
and insert it in (8.12), we get that ∆ is a polynomial in M of degree 6 which is written:

∆ = 4(Z 2 -1) 2 M 6 + 16T Z(Z 2 -1)M 5 + 8T 2 (3Z 2 -1) + 12L + 20LZ 2 + 4Z 3 R -36ZR M 4 + l.o.t,
where "l.o.t" stand for "lower order terms". Suppose Z = ±1 =: ϵ. Then, sin θ = 1 and cos(ϕδ) = ϵ. This means that θ = π/2 and ϕ = δ

(if ϵ = 1) or ϕ = δ + π (if ϵ = -1). This implies L = c 1 Θ sin 2 δ, R = ϵΘ sin 2 δ, T = ϵ(c 1 -c 2 ) cos δ.
In this case, ∆ reduces to

∆ = 16 T 2 + 2(L -ϵR) M 4 + l.o.t = 16 (c 1 -c 2 ) 2 cos 2 δ -2(1 -c 1 )Θ sin 2 δ M 4 + l.o.t
Thus, whenever δ is such that

0 ≤ cot 2 δ < E, E = 2Θ(1 -c 1 ) (c 1 -c 2 ) 2 ,
(we notice that E > 0 thanks to Theorem 3.1), the leading order of ∆ as a polynomial in M is a negative coefficient times M 4 . Thus, ∆ is negative for large enough M . The coefficients of the lower order terms are bounded by constants that only depend on c 1 , c 2 and Θ. Thus, there exists a constant C 2 > 0 which only depends on c 1 , c 2 and Θ such that for |M | > C 2 and δ such that cot 2 δ < E/2, then ∆ < 0, showing that the System is not hyperbolic about a state of corresponding (ρ 0 , u 0 , z 0 ).

Proof of (ii), case |M | small. Suppose now θ = π/2 and ϕ such that L = c 1 Θ sin 2 ϕ = e ≪ 1. Let us also choose δ = ϕ + π/2 so that, for instance:

sin ϕ = ϵ e c 1 Θ , cos ϕ = sin δ = 1 - e c 1 Θ = 1 + O(e),
with ϵ = ±1. Then, we get

R = ϵ eΘ c 1 1 + O(e) , Y = c 1 -c 2 + O(e).
Finally, let us choose M = αe with α a constant to be chosen later. Inserting these assumptions into (8.12), we check that the leading order term when e → 0 is coming from the second term of (8.12), so that

∆ = 4ϵ(c 1 -c 2 ) 3 Θ c 1 αe 3/2 1 + o(1) ,
as e → 0. Taking α = -ϵ, we see that this leading order term is negative. By the same arguments as in the previous case, this shows that there exists a constant C 1 > 0 which only depends on c 1 , c 2 and Θ such that for all |M | < C 1 , we can find δ (depending on |M |) such that ∆ < 0. This proves that the System is not hyperbolic about a state of corresponding (ρ 0 , u 0 , z 0 ) either.

-λ + (c 1 U 1 + 2πbp)2πp + (c 1 U 2 + 2πbm)2πm = 0,
which, after rearrangement, is nothing but (3.36). This gives the travelling-wave solutions. Now, such solutions are stationary if λ = 0. In such a case, supposing (p, m) ̸ = (0, 0), U must satisfy (3.38). Letting X = pe 1 + me 2 , this equation is written

U • X |X| 2 = - 2πb c 1 . (9.1)
Using Cauchy-Schwarz inequality and the fact that |U | = 1, we find that

|X| ≤ c 1 2π|b| , (9.2)
is a necessary condition for the existence of U . But this condition is exactly (3.37). It is also a sufficient condition. If the inequality (9.2) is satisfied, we can write

U = - 2πb|X| c 1 X |X| + σ 1 - 2πb|X| c 1 2 1/2 X |X| ⊥ , = - 2πb c 1 p m + σ 1 p 2 + m 2 - 2πb c 1 2 1/2 -m p (9.3)
where ( X |X| ) ⊥ is the vector obtained by rotating X |X| by an angle of π/2 and σ = ±1. This gives two solutions U except if the factor of σ is zero, which is the case where (9.2) is an equality. In this case, the solution is given by

U = -Sign(b)X/|X| = -Sign(b) 1 p 2 + m 2 p m .
Finally, the case (p, m) = (0, 0) is obvious, which ends the proof.

Numerical methods

In this section, we give additional details on the numerical methods used to produce the simulations shown in Section 4. The code is freely available on the GitHub page of the second author at https://github.com/antoinediez/Swarmalators

The particle scheme is written in Python and the finite volume scheme in Julia.

Particle scheme

Simulating mean-field particle systems is relatively easy though computationally expensive when the number of particles becomes large. In order to simulate the particle system with up to 3.5 millions particles, we rely on the highly-efficient GPU framework introduced in the SiSyPHE library [START_REF] Diez | SiSyPHE: A Python package for the Simulation of Systems of interacting mean-field Particles with High Efficiency[END_REF] which is based on the KeOps library [START_REF] Charlier | Kernel operations on the GPU, with autodiff, without memory overflows[END_REF]. The SiSyPHE library is a versatile Python library designed for the simulation of collective dynamics models which already includes classical models such as the Vicsek model. Thanks to the objectoriented implementation of the library and since the present model is an elaboration of the Vicsek model, only a simple extension of the base class Vicsek of the SiSyPHE library is needed in order to incorporate the new phase variable and its contribution to the dynamics. From a methodological point of view, the (stochastic) particle system (2.1)-(2.3) is discretized using a first-order Euler-Maruyama scheme. For a given interaction radius R and with the notations of Section 4.1.1, the time-step is taken equal to ∆t = 10 -2 / max(Cγ, ν, D, ν ′ , D ′ ). All the particle simulations have been run using an Nvidia GTX 2080 Ti GPU chip on the GPU cluster of the Department of Mathematics at Imperial College London.

Finite volume scheme

Following the methodology introduced in [START_REF] Motsch | Numerical simulations of a nonconservative hyperbolic system with geometric constraints describing swarming behavior[END_REF] for the SOH model, the finite volume scheme is based on the following formulation of the system (3.4)-(3.6) as the relaxation limit ε → 0 of a system written in conservative form.

∂ t ρ ε + ∇ x • [ρ(c 1 u ε + bρ ε ∇ x α ε )] = 0, (10.4)
∂ t (ρ ε u ε ) + ∇ x • [ρ ε u ε ⊗ (c 2 u ε + bρ ε ∇ x α ε )] + Θ∇ x ρ ε = ρ ε ε (1 -|u ε | 2 )u ε , (10.5)
∂ t (ρ ε cos α ε ) + ∇ x • [ρ ε cos α ε (c 1 u ε + bρ ε ∇ x α ε )] = -sin α ε (b -b ′ )|ρ ε ∇ x α ε | 2 + Θ ′ ∇ x • (ρ ε ∇ x ρ ε ) , (10.6) ∂ t (ρ ε sin α ε ) + ∇ x • [ρ ε sin α ε (c 1 u ε + bρ ε ∇ x α ε )] = cos α ε (b -b ′ )|ρ ε ∇ x α ε | 2 + Θ ′ ∇ x • (ρ ε ∇ x ρ ε ) , (10.7)
where for numerical stability reasons, we use the variable (cos α, sin α) instead of just α.

We first solve the conservative part using a custom HLLE scheme [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF][START_REF] Einfeldt | On Godunov-type methods near low densities[END_REF] and we use a splitting method for the source terms as outlined below.

1. Using a dimensional splitting, solving the conservative part

∂ t ρ + ∇ x • [ρ(c 1 u + bρ∇ x α)] = 0, ∂ t (ρu) + ∇ x • [ρu ⊗ (c 2 u + bρ∇ x α)] + Θ∇ x ρ = 0, ∂ t (ρ cos α) + ∇ x • [ρ cos α(c 1 u + bρ∇ x α)] = 0 ∂ t (ρ sin α) + ∇ x • [ρ sin α(c 1 u + bρ∇ x α)] = 0,
reduces to solving two 1D equations. In order to compute the numerical flux between two datas (ρ ℓ , u 1ℓ , u 2ℓ , cos α ℓ , sin α ℓ ) and (ρ r , u 1r , u 2r , cos α r , sin α r ) we first need to approximate the phase gradient (in the x 1 -direction). A finite difference approximation can be computed from (cos α ℓ , sin α ℓ ) and (cos α r , sin α r ) only by taking the argument of the following complex number:

∂ x 1 α ≃ z ℓr := arg cos α r cos α ℓ + sin α r sin α ℓ + i(sin α r cos α ℓsin α ℓ cos α r ) /∆x, (10.8) where ∆x is the space discretization step. Then, using the change of variable (q 0 , q 1 , q 2 , q 3 , q 4 ) = (ρ, ρu 1 , ρu 2 , ρ cos α, ρ sin α) we are led to the computation of the Jacobian matrix of the flux function f (q 0 , q 1 , q 2 , q 3 , q 4 ) :=       c 1 q 1 + bz ℓr q 2 0 c 2 q 2 1 /q 0 + Θq 0 + bz ℓr q 0 q 1 c 2 q 1 q 2 /q 0 + bz ℓr q 0 q 2 c 1 q 3 q 1 /q 0 + bz ℓr q 0 q 3 c 1 q 4 q 1 /q 0 + bz ℓr q 0 q 4      

.

A direct computation shows that the Jacobian matrix of f has four eigenvalues, one with multiplicity 2: ν 1 = c 1 q 1 /q 0 + bz ℓr q 0 , and three with multiplicity 1: ν 2 = c 2 q 1 /q 0 + bz ℓr q 0 , ν + = 1 2 2c 2 q 1 /q 0 + 3bz ℓr q 0 + √ ∆ , ν -= 1 2 2c 2 q 1 /q 0 + 3bz ℓr q 0 -√ ∆ ,

where ∆ = 4c 2 (c 2c 1 )(q 1 /q 0 ) 2 + 4c 1 Θ + 4bz ℓr (c 1c 2 )q 1 + (bz ℓr q 0 ) 2 . Using these values, the computation of the numerical flux using a HLLE scheme is explained in [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF][START_REF] Einfeldt | On Godunov-type methods near low densities[END_REF]. Note that since the eigenvalues depend on z ℓr and thus on ρ and ∇ x α (which are not ensured to be bounded), we have to use an adaptive time step ∆t in order to guarantee the CFL condition ∆t ∆x max(|ν 1 |, |ν 2 |, |ν + |, |ν -|) ≤ 1 at each iteration. For better stability, in the experiments, ∆t is chosen so that the CFL number does not exceed 0.1.

The relaxation part reads

∂ t ρ ε = 0, ∂ t (ρ ε u ε ) = ρ ε ε (1 -|u ε | 2 )u ε , ∂ t (ρ ε cos α ε ) = 0, ∂ t (ρ ε sin α ε ) = 0.
It can be solved explicitly but as shown in [START_REF] Motsch | Numerical simulations of a nonconservative hyperbolic system with geometric constraints describing swarming behavior[END_REF], when ε → 0 it reduces to a mere normalization of the velocity. We solve this part using an explicit Euler scheme. We use a finite difference approximation of the spatial derivatives on the right-hand side, with a classical five-point discretisation of the Laplacian term and using the same method as before (10.8) for the gradient in α. The system reaches a kind of doubly periodic travelling wave solution but with a much more complex shape.

Simulations of the SH system

The following videos supplement the results presented in Section 4.3.2.

Video 19 (Low noise NA) Simulation of the SH system starting from a perturbed NA configuration with k ′ = k ′ hi = 200 and the other parameters given in Section 4.2.1. The initial noise quickly resorbs and the system finally reaches a stable doubly periodic travelling wave with no perceptible inhomogeneities.

Video 20 (Low noise PA) Simulation of the SH system starting from a perturbed PA configuration with k ′ = k ′ hi = 200 and the other parameters given in Section 4.2.1. The initial noise quickly resorbs and the system finally reaches a stable doubly periodic travelling wave although small inhomogeneities in the density are perceptible.

Video 21 (Low noise OT) Simulation of the SH system starting from a perturbed OT configuration with k ′ = k ′ hi = 200 and the other parameters given in Section 4.2.1. The initial noise quickly resorbs but small inhomogeneities in the density are perceptible and eventually degenerate into thin very-low density regions. After a transition period during which the global velocity transitions towards a limit value close to (0, 1) T , these structures disappear and the system reaches the corresponding stable doubly periodic travelling wave solution (although although small inhomogeneities in the density remain perceptible).

Video 22 (Medium noise NA) Simulation of the SH system starting from a perturbed NA configuration with k ′ = k ′ lo = 10 and the other parameters given in Section 4.2.1. The initial noise quickly resorbs and the system finally reaches a stable doubly periodic travelling wave with no perceptible inhomogeneities.

  (3.28) preserves this structural constraint as the exterior derivative of a gradient is always zero. Now, the relation between the NSH System (3.23)-(3.25) and System (3.26)-(3.29) is expressed in the following lemma, the proof of which is obvious. Lemma 3.5 (i) Let (ρ, u, α) be a solution of the NSH System (3.23)-(3.25). Then, (ρ, u, z) with z = ∇ x α is a solution of System (3.26)-(3.29).

. 30 )

 30 In other words, a uniform steady-state solution for System (3.26)-(3.29) is a travellingwave solution of the NSH System (3.23)-(3.25). We take the spatial Fourier transform of the linearized system to (3.26)-(3.29) and denote by ξ the Fourier variable and τ = ξ/|ξ| its direction. Using frame indifference, we introduce a reference frame (e 1 , . . . , e n ) such that u 0 ∈ Span{e 1 }, z 0 ∈ Span{e 1 , e 2 } and ξ ∈ Span{e 1 , e 2 , e 3 }. Then, we introduce the angles

1 Figure 2 :

 12 Figure 2: Geometric setting of the hyperbolicity analysis

Remark 3 . 3

 33 (i) In the case z 0 = 0, the eigenvalues λ 1 and λ ± are the eigenvalues of the SOH model (3.19), (3.20)[START_REF] Degond | Hydrodynamic models of selforganized dynamics: derivation and existence theory[END_REF]. The phase equation only adds the extra eigenvalue λ 2 . In this case, no eigenmode depends on the coupling constant b.

Proposition 3 . 7

 37 (i) Travelling-wave solutions: let (m, p) ∈ Z 2 and let U ∈ S 1 be arbitrary. Then, the following is a periodic travelling-wave solution of the NSH system (3.23)-(3.25) (with V = 0) in Ω satisfying the normalization condition (3.7):

. 38 )

 38 For these two choices of U (respectively unique choice of U ), then λ = 0 and the solution (3.33)-(3.35) is a stationary solution. (iii) Stationary solutions with constant phase: If (p, m) = (0, 0), then λ = 0 and α = α 0 is a constant. Then the solution (3.33)-(3.35) is a stationary solution for any choice of U ∈ S 1 .

1 Figure 3 :

 13 Figure3: Plot of the phase field α(x, 0) (when α 0 = 0) given by(3.35) in the case (p, m) = (1, 2). The vectors e iα(x,0) at the points (x 1 , x 2 ) = ( k 8 , ℓ 8 ), for k, ℓ in {0, 1, . . . , 8} are shown. The color code corresponds to the angle α with red, pink, yellow, green, light blue, blue, magenta, purple corresponding to α = 0, π 4 , π 2 , 3π 4 , π, 5π 4 , 3π 2 , 7π 4 respectively. The indices of this vector field about the origin are 1 and 2 when one moves along the x 1and x 2 -axes respectively.

Figure 4 :

 4 Figure 4: Twenty independent measurements of the x 2 -coordinate of the center of mass of the particle system over 0.1 units of time for three different values of N . The thick dashed red line is the theoretical prediction derived from the macroscopic model. For each value of N , twenty independent experiments are displayed in colored plain thin lines. For each value of N , the number N neigh = πR 2 N is the average number of neighbours in a ball of radius R. Parameters: R = 0.005, ν = 3, D = 1, ν′ = 5, D′ = 1.

  (a) Speed for various N (b) Standard deviation around the mean value

Figure 5 :

 5 Figure 5: Convergence of the particle scheme as N → +∞. (a) For different values of N , twenty independent experiments are run and for each of them the measured travellingwave speed is computed as the slope coefficient of x N 2 (t) using a standard linear regression. The blue points correspond to the mean value over the twenty experiments for each value of N and the error bar in orange shows the standard deviation. The theoretical speed is equal to 1 as indicated by the dotted black horizontal line. The inset shows the absolute value of the difference between the mean value of the measured speed averaged over the 20 simulations and the theoretical speed in log-scale. (b) The blue points show the standard deviation around the mean value of the measured travelling-wave speed for each value of N . The orange line is the regression line (excluding the first point) and has a slope coefficient equal to 0.48 which is close to the value 0.5 predicted by the central limit theorem. Parameters: R = 0.005, ν = 3, D = 1, ν′ = 5, D′ = 1.

Figure 6 :

 6 Figure 6: Final state after 40 units of time starting from a NA configuration for four values of k ′ . For each value of k ′ , the domain is discretized into a uniform grid with 10 4 cells of size 0.01. (a),(c),(e),(g) The density of particles is obtained by counting the proportion of particles in each cell. (b),(d),(f),(h) The phase in each cell is the average phase of the particles in this cell. It is arbitrarily set to 0 when the cell is empty. The corresponding videos can be found in the supplementary material: (a)-(b) Video 1 (c)-(d) Video 5 (e)-(f) Video 9 (g)-(h) Video 13. Parameters: N = 10 6 , R = 0.01, ν = 5, D = 1, γ = 0.2.

Figure 7 :

 7 Figure 7: Number of particles in the strips [0, 1] × [ k 100 , k+1 100 ] for k ∈ {0, . . . , 99} after 40 units of time starting from the NA configuration with k ′ = k ′ med . Same parameters as Fig. 6e.

Inserting ( 7 . 3 ) 1 = π 0 cos θ e k cos θ sin n-2 θ dθ π 0 e k cos θ sin n- 2

 731002 into (7.1) leads to (3.4) with c

  θ) e k cos θ sin n θ dθ π 0 e k cos θ sin n-2 θ dθ , i = 1, 2. (7.30) Dividing (7.29) by kC 1 , we get (3(cos θ) e k cos θ sin n θ dθ π 0 h(cos θ) e k cos θ sin n θ dθ , (7.31)

  1 φ h(cos φ) e k ′ cos φ sin 2 φ dφ 2π 0 e k ′ cos φ dφ i = 1, 2. (7.34) Dividing by C ′ 1 k ′ , we get (3.6) with b

  and (ϕ = 0 or ϕ = π).(8.13) In this case, we have Y = ϵ 2 (c 1c 2 ) + ϵ 3 M with ϵ 2 = cos ϕ = ±1 and ϵ 3 = cos δ = ±1. Therefore, we have ∆ = 0 if and only if

9

  Doubly periodic travelling-wave solutions: proofs We look for a solution of the form (3.33), (3.34), (3.35). Since ρ, u and ∇ x α = 2π(pe 1 + me 2 ) are constant in space and time, Eqs (3.23) and (3.24) (with V = 0) are trivially satisfied. The only equation left to verify is (3.25); It leads to

3 .

 3 The other source terms reads∂ t ρ = 0, ∂ t u = 0 ∂ t cos α =sin α (bb ′ )ρ|∇ x α| 2 + Θ ′ (∆ρ + |∇ x √ ρ| 2 ) , ∂ t sin α = cos α (bb ′ )ρ|∇ x α| 2 + Θ ′ (∆ρ + |∇ x √ ρ| 2 ) .

Video 5 (

 5 Low noise NA) Particle simulation starting form a NA configuration with k ′ = k ′ hi = 200 and the other parameters given in Section 4.2.1. The behavior is the same as the one described in the caption of Video 1 but the bands are larger and their phase slowly varies in time. See Figs.6c and 6d. Video 6 (Low noise PA) Particle simulation starting form a PA configuration with k ′ = k ′ hi = 200 and the other parameters given in Section 4.2.1. The behavior is initially the same as the one described in the caption of Video 2 but the bands are not stable and the systems finally reaches a flocking phase. Video 7 (Low noise OT) Particle simulation starting form a OT configuration with k ′ = k ′ hi = 200 and the other parameters given in Section 4.2.1. The behavior is initially the same as the one described in the caption of Video 3 but the systems finally reaches a stable state which is the same state as the one starting from the NA configuration. Video 8 (Low noise UF) Particle simulation starting form a UF configuration with k ′ = k ′ hi = 200 and the other parameters given in Section 4.2.1. The behavior is the same as the one described in the caption of Video 4 but the segregation regions are larger. Video 9 (Medium noise NA) Particle simulation starting form a NA configuration with k ′ = k ′ med = 10 and the other parameters given in Section 4.2.1. The initial doubly periodic travelling wave persists during about 5 units of time. Then thin low-density regions emerge and delimitate band-like structures. Unlike the cases with less noise in the phase equation, the bands are larger and an inner gradient in phase can be identified inside each band. Video 10 (Medium noise PA) Particle simulation starting form a PA configuration with k ′ = k ′ med = 10 and the other parameters given in Section 4.2.1. The behavior is similar to the one presented in the caption of Video 9 except that the bands are moving in the opposite direction. Moreover, the bands shape is less stable. Video 11 (Medium noise OT) Particle simulation starting form a OT configuration with k ′ = k ′ med = 10 and the other parameters given in Section 4.2.1. The behavior is initially similar to the one presented in the caption of Video 9 except that the global velocity of the particles transition from (1, 0) T to (0, 1) T and the systems finally reaches the same stable state as the one starting from a NA configuration. Video 12 (Medium noise UF) Particle simulation starting form a UF configuration with k ′ = k ′ med = 10 and the other parameters given in Section 4.2.1. The behavior is the same as the one described in the caption of Video 8. Video 13 (Large noise NA) Particle simulation starting form a NA configuration with k ′ = k ′ lo = 3 and the other parameters given in Section 4.2.1. The initial doubly periodic travelling wave is stable and persists during the 40 units of time of the simulation. The theoretical travelling wave speed is λ ≃ 0.03. The measured speed is approximately equal to 0.05. Video 14 (Large noise PA) Particle simulation starting form a PA configuration with k ′ = k ′ lo = 3 and the other parameters given in Section 4.2.1. The initial doubly periodic travelling wave is stable and persists during the 40 units of time of the simulation. The theoretical travelling wave speed is λ ≃ -1.75. The measured speed is approximately equal to -1.72. Video 15 (Large noise OT) Particle simulation starting form a OT configuration with k ′ = k ′ lo = 3 and the other parameters given in Section 4.2.1. The initial doubly periodic travelling wave is stable and persists during the 40 units of time of the simulation. The theoretical travelling wave speed is λ ≃ -0.86. The measured speed is approximately equal to -0.84. Video 16 (Large noise UF 1) Particle simulation starting form a UF configuration with k ′ = k ′ lo = 3 and the other parameters given in Section 4.2.1. The system reaches a flocking state (although two phases can be identified) Video 17 (Large noise UF 2) Particle simulation starting form a UF configuration with k ′ = k ′ lo = 3 and the other parameters given in Section 4.2.1. The system reaches a state close to a doubly periodic travelling wave but there are also waves of velocities. Video 18 (Large noise UF 3) Particle simulation starting form a UF configuration with k ′ = k ′ lo = 3 and the other parameters given in Section 4.2.1.

Video 23 (

 23 Medium noise PA) Simulation of the SH system starting from a perturbed PA configuration with k ′ = k ′ lo = 10 and the other parameters given in Section 4.2.1. Same observations as Video 22. Video 24 (Medium noise OT) Simulation of the SH system starting from a perturbed OT configuration with k ′ = k ′ lo = 10 and the other parameters given in Section 4.2.1. Same observations as Video 22. Video 25 (Large noise NA) Simulation of the SH system starting from a perturbed NA configuration with k ′ = k ′ lo = 3 and the other parameters given in Section 4.2.1. Same observations as Video 22. Video 26 (Large noise PA) Simulation of the SH system starting from a perturbed PA configuration with k ′ = k ′ lo = 3 and the other parameters given in Section 4.2.1. Same observations as Video 22. Video 27 (Large noise OT) Simulation of the SH system starting from a perturbed OT configuration with k ′ = k ′ lo = 3 and the other parameters given in Section 4.2.1. Same observations as Video 22.

PD holds a visiting professor association with the Department of Mathematics, Imperial College London where part of this research was conducted. Part of this research was conducted when AD was supported by an EPSRC-Roth scholarship cofunded by the Engineering and Physical Sciences Research Council and the Department of Mathematics at Imperial College London. No new data were collected in the course of this research.

List of supplementary videos

The videos can be found using the following link: https://figshare.com/projects/Topological_states_and_continuum_model_for_swarmalators_ without_force_reciprocity/139912

The videos show the outcome of the simulations discussed in Sections 4.2 and 4.3. For each video, the left panel shows the spatial density and the right panel the average phase. These quantities are computed using a spatial discretization on a uniform grid with constant step ∆x = 0.01. For the particle simulations, the value of the density in a cell is computed as the proportion of particles in this cell and the phase is their average phase (it is arbitrarily set to 0 if the cell is empty). Similarly the average velocity is computed on a spatial grid with step size ∆x = 0.05 and depicted by black arrows. For the particle simulations, three particles are represented on the right panel by three disks (of arbitrary radius) colored according to the phases of these particles.

Particle simulations

The following videos supplement the results presented in Section 4.2.2 and discussed in Section 4.2.3. Note that the phases of the particles range from 0 to 2π but they have been rescaled so that the range of the colorbar of the right-panel is between 0 and 1.

Video 1 (Very low noise NA) Particle simulation starting form a NA configuration with k ′ = k ′ hi+ = 10000 and the other parameters given in Section 4.2.1. The particles segregate into small regions of equal phase separated by thin low-density regions. The systems eventually reaches a stable state where well-separated band-like structures characterized by a constant phase travel at a constant speed in the direction opposite to the phase gradient. See Figs. 6a and6b.

Video 2 (Very low noise PA) Particle simulation starting form a PA configuration with k ′ = k ′ hi+ = 10000 and the other parameters given in Section 4.2.1. The particles segregate into small regions of equal phase separated by thin low-density regions. After 40 units of time, well-separated band-like structures characterized by a constant phase can be identified. They travel at a constant speed in the same direction as the phase gradient but their shapes are not stable.

Video 3 (Very low noise OT) Particle simulation starting form a OT configuration with k ′ = k ′ hi+ = 10000 and the other parameters given in Section 4.2.1. The particles segregate into small regions of equal phase separated by thin low-density regions. The global velocity of the particles transitions from (1, 0) T to (0, 1) T . Although band-like structures can be identified, they are not as stable as in the NA case.

Video 4 (Very low noise UF) Particle simulation starting form a UF configuration with k ′ = k ′ hi+ = 10000 and the other parameters given in Section 4.2.1. The particles segregate into small regions of equal phase separated by thin low-density regions. These regions are very dynamic with a lot of merging and mixing. After about 10 units of time, the system reaches a flocking state in which all the phases and velocity equal.