Snowdrift game induces pattern formation in systems of self-propelled particles

Evolutionary games between species are known to lead to intriguing spatio-temporal patterns in systems of diffusing agent. However, the role of inter-species interactions is hardly studied when agents are (self-)propelled, as is the case in many biological systems. Here, we combine aspects from active-matter and evolutionary game theory and study a system of two species whose individuals are (self-)propelled and interact through a snowdrift game. We derive hydrodynamic equations for the density and velocity fields of both species from which we identify parameter regimes in which one or both species form macroscopic orientational order as well as regimes of propagating wave patterns. Interestingly, we find simultaneous wave patterns in both species that result from the interplay between alignment and snowdrift interactions - a feedback mechanism that we call game-induced pattern formation. We test these results in agent-based simulations and confirm the different regimes of order and spatio-temporal patterns as well as game-induced pattern formation.


I. INTRODUCTION
The composition and temporal evolution of microbial as well as other populations depends on a variety of factors including environmental conditions, the population structure, and the degree of mobility of each species [1][2][3]. Evolutionary game theory (EGT) [4,5] has been used as a mathematical framework to conceptually describe the evolutionary dynamics of such populations employing methods from nonlinear dynamics [6]. Two-player games, including the prisoner's dilemma and the snowdrift game, have played an important role in the development of the research field. In EGT, games are regarded as a population dynamics problem, where individuals who follow different strategies ("species") interact over many generations according to rules set by the respective game. The ensuing (nonlinear) dynamic process then determines which species survive, i.e. 'win the game'.
The mobility of individuals in particular has been shown to play a crucial role in the evolutionary success of a species. Previous theoretical studies have modeled mobility either as a diffusive process, in which individuals perform random walks (see e.g. Refs. [7][8][9]) or alternatively as a migration process, in which individuals move by winning a competition with other species that are spatially adjacent to them (see e.g. Refs. [10][11][12][13][14]). To the best of our knowledge, however, there are no studies that deal with the effect of self-propelled movement -a characteristic feature of many living beings -on the outcome of spatially extended evolutionary games. Since collectives of self-propelled particles (SPP) are known to generically exhibit intriguing self-organisation phenomena like flocking, clustering, and motility-induced phase separation(for reviews see e.g. Refs. [15][16][17][18][19]), this will affect the spatial proximity between individuals and thereby their competitive interaction (described as a game). Changes in the composition of the population could in turn affect the ordering processes of the SPPs.
Here, we investigate this interplay of collective effects due to self-propelled motion and nonlinear dynamics based on competitive (game theory) interactions. For specificity, we will address this question for the snowdrift game, one of the classical two-player games. Using both agent-based simulations and a hydrodynamic approach, we discover a phenomenon that we term game-induced pattern formation: When one of the two species starts to exhibit collectively moving patterns due to self-propelled motion and alignment of its agents, a pattern is induced in the second species. We attribute this to a non-linear feedback mechanism between the local snowdrift game interaction and the alignment between particles. Upon exploring the model parameters that determine the alignment and the interaction probability, we find that the system exhibits five distinct phases, each with qualitatively different collective behaviour. If the game interaction is sufficiently weak, we find that reducing the noise level leads to a phase transition from a disordered phase to an ordered phase in which one species induces a polar wave pattern in the second species. These patterns vanish when the noise strength is reduced, leading to a partially ordered phase in which one species shows uniform polar order while the second species is disordered. A further reduction of the noise strength then also leads to an ordering transition in the second species, which first induces a wave pattern in the previously uniformly ordered species and finally also shows uniform polar order. All of these phases are generalizations of the phases observed in single-species systems of SPPs, and we show that their existence depends on the interplay between alignment and inter-species interactions.
This paper is organized as follows: In the remainder of the introduction we give a concise overview over previous work on two-player games in a well-mixed environment and in spatially-extended systems. This is mainly intended to inform readers in the active matter field not familiar with this topic about those results relevant for the present study. In Section II we first summarize the Boltzmann approach for SPPs and the dynamics of the snowdrift game to then combine both models to find equations governing the dynamics of two SPP species with inter-species snowdrift game interaction. Section III takes these equations as a starting point to derive simpler hydrodynamic equations and analyses them using both analytical and numerical methods. In this section we obtain the phase diagram of the system in the space spanned by parameters of self-propelled motion and game interaction and discuss their properties. The results of the hydrodynamic approach and the overall phenomenology of the system are qualitatively validated in Section IV by agent-based simulations. We conclude with a summary and discussion of our results in Section V.
A. Two-player public good games Originally, the 'snowdrift' game (also known as 'hawkdove' game) was formulated as a strategic 'cooperatordefector' game [20]. It owes its name to its prominent explanatory framework: Two drivers are caught in a winter storm and a snowdrift is blocking the road. Both drivers now have the option of either remaining in the car or getting out and removing the snowdrift. Since both drivers want to continue their journey, but getting out is not really desirable, the best strategy for each driver is to do exactly the opposite of what the other driver does.
This strategic game can be transferred to the field of EGT and microbial systems such as yeast cells competing for resources and using different survival strategies [21]: In a metabolically costly process, part of the yeast cells (producer cells) convert sucrose to glucose, which it can use as an energy source. This process yields more glucose than can be metabolized, and the extra glucose diffuses into the surrounding medium. There it can be used by all yeast cells, regardless of whether they produce glucose or not. As long as there is an excess of glucose in the system, it is therefore advantageous to save the metabolic costs for its production and only consume it (non-producer cells). In contrast, when there is a lack of glucose, it becomes advantageous to produce glucose as this provides an energy source. This leads to the coexistence of producer (cooperator) and non-producer (defector) cells in the yeast population.
The snowdrift game is one out of four possible twoplayer ('cooperator-defector') games. These include the prisoner's dilemma, the coordination, and the harmony game. Their dynamics in a well-mixed environment (where diffusion is much faster than the reaction kinet-ics) can be formulated in terms of concepts from nonlinear dynamics [6]. The main difference between these four games is their evolutionary outcome, i.e. the composition of the population in the long run. While, for example, in the snowdrift game, cooperators and defectors coexist, defectors always dominate in the prisoner dilemma and cooperators go extinct. This is related to the rules of the prisoner's dilemma, which can be formulated as a public goods game where an individual can either cooperate by providing some kind of public good or defect by exploiting this public good. Under these conditions, defectors will prevail as they can benefit from the interaction with cooperators without paying any costs.

B. Evolutionary games in spatially extended systems
The spatial extension of a population and the ensuing spatio-temporal arrangement of individuals can strongly influence the evolutionary dynamics [22][23][24][25]. It is well known that this plays an important role for microbial life, as many microbes on our planet are found in the form of dense microbial communities such as colonies and biofilms [23,26].
A key factor for the dynamics is the spatial clustering of individuals, which has already been argued by Hamilton to support cooperation in populations [27], because cooperating individuals are more likely to interact with each other and therefore are less likely to be exploited by defectors. To what degree this determines the composition of a population has been studied extensively on different kinds of lattices as well as for different kinds of interaction rules [8,[10][11][12][28][29][30][31][32][33][34]. Specifically, for the snowdrift game, it was shown that accounting for spatial extension can change population ratios compared to well-mixed systems [29,35].
While all the studies mentioned above consider diffusive dynamics, self-propulsion is an essential part of biological reality. In fact, there are a large variety of microbes that show directed mobility, often along chemotactic gradients [36][37][38][39] or due to chemical warefare [40]. Here, motility is an important factor of their strategy, e.g. for directed motion away from hostile environments, chemotaxis or self-organisation into flocks.

II. MODEL DESCRIPTION
A. Active matter model In order to implement propulsion of particles of both species, A and B, we employ a kinetic Boltzmann approach [41,42]. Here, particles with a diameter d 0 move ballistically on a two-dimensional substrate with constant speed v 0 and orientation θ. Ballistic motion is interrupted by stochastic self-diffusion events: Particles change their orientation by an angle η at a rate ω to where η is a stochastic variable (noise); see Fig. 1a. Hence, the particles move in a run-and-tumble-like fashion with a mean time ω −1 between tumbling events. Further, we assume that intra-species interactions are due to binary collisions between particles, that lead to alignment of their orientations, while inter-species interactions follow a snowdrift game (see the following section II B). We implement intra-species aligning collisions in terms of half-angle alignment rules; see Fig. 1b. The particle orientations after a collision event are given by where we allowed for additional noise terms η 1 and η 2 .
We have chosen these rules to keep our analysis conceptually simple. In actual systems one would need to account for the problem-specific interaction between the agents of the active system. This has, for example, been achieved for the actin-myosin motility assay [43][44][45]. For simplicity, diffusion and collision noise variables are drawn from a Gaussian distribution with the same standard deviation σ.
FIG. 1. Illustration of self-propelled dynamics: (a) Particles change their orientation θ by a stochastic angle η to θ = θ+η during self-diffusion events. (b) Binary collisions between particles of the same species lead to polar half-angle alignment of their orientations. We allow for additional stochastic noise terms η1 and η2 during the collision event.
For a dilute system of only one species of propelled particles, Bertin et al. [41,42] proposed a kinetic Boltzmann equation. It describes the dynamics of the oneparticle probability density function, f (r, θ, t), which indicates the probability that a particle is located at position r with the orientation θ at time t. The kinetic Boltzmann equation for self-propelled particles for one species reads [41,42] where I diff and I coll denote contributions resulting from diffusion and collision events, respectively, while e θ is the unit vector pointing in θ-direction. Their explicit expressions can be found in Appendix A 1. In the present context, the most important feature of I diff and I coll is that they depend on the parameter σ, which characterizes the noise strength.

B. Snowdrift game
In addition to the alignment interactions, we assume that particles during collisions also play a snowdrift game. The outcome of this interaction is specified by the payoff matrix (Fig. 2a), which is non-zero only if two different species meet. The fitness of a particular strategy, A or B, is defined as a constant background fitness (set to 1) plus the average payoff obtained from playing the game [6,46]: Here, a and b denote the densities of species A and B, respectively. In the following, we will refer to λ as the relative fitness parameter of particles of species B relative to A, and τ as the game strength.
The time evolution now follows 'survival of the fittest': The species with higher fitness value increases in population size, while the other species decreases in number. Specifically, the dynamics of the densities a and b are given by the replicator equations: These replicator equations govern the dynamics of the game. The fitness difference (f A − f B ) indicates the relative strength of the interaction partner, while the factor a b encodes the probability that particles of species A and B meet and interact with each other. Note that the overall density ρ = a+b is constant. Inserting the expressions for the fitness values f A and f B , Eqs. (4), into Eqs. (5) one can rewrite the replicator equations as The meaning of the terms relative fitness and game strength for λ and τ , respectively, now become evident: If b > λa, there is a fitness advantage for species A and its density grows relative to species B. In contrast, if b < λa, the growth of species B is favored relative to species A. Since the equations are anti-symmetric under the exchange of species A and B, we can, without loss of generality, impose that species A is always fitter than species B, i.e. constraining λ ∈ [0, 1]. By changing λ we can tune the relative fitness from the two species being equally competitive (λ = 1) to species A being strongly dominant over B (λ 1). The second game parameter, τ ∈ [0, ∞), can be understood as the game strength.
Since the total particle number is conserved (b =ρ − a), Eqs. (6) can be reduced to the temporal evolution of a single species (say A): Figure 2b shows the flow diagram of Eq. (7) with arrows pointing to the right for increasing a (∂ t a > 0) and to the left for ∂ t a < 0. There are three fixed points (stationary points with ∂ t a = 0): The fixed points a B and a A are unstable in the snowdrift game, while a 0 is stable. For a B = 0 and a A = 0 species A and B go extinct, respectively, and therefore the game ends (absorbing fixed points). The fixed point a 0 refers to a stable coexistence of both species with population ratio b/a = λ.
C. Active matter snowdrift game To describe our system of two interacting self-propelled particle species, we introduce the one-particle probability density functions α (r, θ, t) and β (r, θ, t), governing the temporal evolution of species A and B, respectively. We use the kinetic Boltzmann equation, Eq. (3), and extend it for each species separately by snowdrift interaction terms, I A game [α, β] and I B game [β, α], chosen in analogy to the replicator equations, Eqs. (6). This leads to the general form of the two-species Boltzmann equations We assume snowdrift interactions to be local and independent of particle orientation as in Eqs. (6). With the local densities a (r, t) and b (r, t) of species A and B, related to α and β by a (r, t) = π −π dθ α (r, θ, t) , the snowdrift interaction terms read The dependency on α and β in Eq. (11)a and Eq. (11)b, respectively, is chosen to match the form of the Boltzmann equation with its dependency on the one-particle densities α and β. The game strength τ tunes the relative weight of the snowdrift interaction compared to the alignment interaction (collision terms). If one sets τ = 0, the snowdrift interaction terms vanish and Eqs. (9) reduce to two independent equations for species A and B.
In the following, time, space, and densities are rescaled such that one measures time in units of 1/ω, position in units of v 0 /ω and spatial densities in units of ω/(v 0 d 0 ). This amounts to setting v 0 = d 0 = ω = 1.

III. HYDRODYNAMIC THEORY
In the previous section, we have formulated an extended version of two coupled kinetic Boltzmann equations for a dilute system of two species of self-propelled particles with a snowdrift interaction. It includes convection, self-diffusion, intra-species binary collision (with polar alignment), as well as a snowdrift interaction between particles. As solving the kinetic Boltzmann equations, Eqs. (9), analytically is generally impractical, one has to resort to numerical implementations [45,47,48]. However, as exemplified by studies of the kinetic Boltzmann equation for a single species [41,42,[48][49][50][51], one can also use Eqs. (9) to derive hydrodynamic equations for those collective variables that vary slowly in space and time.
In the next section, we first derive these hydrodynamic equations for the system of two species, Eqs. (9), and then study them analytically for spatially uniform systems. Then we investigate spatially non-uniform solutions using linear stability analysis of the uniform states and numerical solutions of the full hydrodynamic equations.

A. Derivation of hydrodynamic equations
As the probability densities α(r, θ, t) and β(r, θ, t) are periodic functions in θ, it is convenient to work with their Fourier modes, defined as α k (r, t) := π −π dθ e iθk α(r, θ, t) and β k (r, t) := π −π dθ e iθk β(r, θ, t). Note that the zeroth Fourier modes are the densities of the respective species, α 0 = a, β 0 = b, defined by Eqs. (10). Besides the densities of species A and B, the second set of slow, hydrodynamic variables are the polar order fields P A (r, t) and P B (r, t), which are related to the first Fourier modes Hence, the modes α 1 and β 1 are a measure for the polar order multiplied by the density at a given point in space and time of the system. If α 1 = β 1 = 0, the system is disordered. If α 1 or β 1 is nonzero, there is a preferred direction of motion, meaning the respective species moves collectively. We will refer to α 1 and β 1 in the following as velocity fields.
The expansion of the kinetic Boltzmann equations, Eqs. (9), in Fourier modes is given by (see Appendix A for details) where we introduced the abbreviations ∇ = ∂ x + i∂ y and ∇ * = ∂ x − i∂ y . The explicit expressions for the collision kernels I p,k (σ) are given in Appendix A. Equation (13) constitutes an infinite set of differential equations that couple lower-with higher-order Fourier modes. For k = 0, Eqs. (13) yield The first terms in Eqs. (14) can be understood as the temporal change of the densities due to spatial variations in the velocity fields, while the last terms can be interpreted as source or sink terms given by the replicator equations, Eqs. (6). These source terms couple the densities of species A and B and we will refer to them in the following as snowdrift terms. Note, that the sum of both snowdrift terms is zero at every point in space and time. Hence, the overall density ρ = a + b follows a continuity equation The spatial average of the density,ρ, is therefore a conserved quantity and a control parameter of the system. For a spatially uniform system, the spatial derivatives in the equations for the particle densities, Eqs. (14), vanish and Eqs. (14) reduce to the replicator equations. The corresponding steady states are given by the stable snowdrift fixed point with a (0) =ρ/(1 + λ) and b (0) =ρλ/(1 + λ). Hence, a state with uniform densities a = a (0) and b = b (0) and all higher Fourier modes vanishing is a solution to Eqs. (13). This solution describes a spatially uniform, disordered system with densities set to their snowdrift fixed points. Small spatially uniform perturbations δα k , δβ k to this solution evolve to linear order as The respectively first two terms in µ A k and µ B k are due to the self-propelled motion of both species and are familiar expressions from one-species SPP systems [41,42]. Here, we encounter an additional term due to the snowdrift game, which implies that the growth rate of every Fourier mode depends on the densities of both species. When we insert the stable solution for the densities, a (0) =ρ/(1+λ) and b (0) =ρλ/(1 + λ), into the expressions for µ A k and µ B k , the snowdrift terms vanish and we are left with the first two terms, which are then functions of σ,ρ and λ. Varying these parameters, one finds that all Fourier modes µ A k and µ B k with k ≥ 2 are negative and only the growth rates for the first Fourier mode, µ A 1 and µ B 1 , can become positive, meaning that small perturbations in the velocity fields α 1 or β 1 , respectively, will grow exponentially. For fixed overall density,ρ, and relative game strength, λ, the growth rates µ A 1 and µ B 1 are negative for high values of the noise strength σ and become positive when the noise strength is low enough. This defines threshold values for the noise strength, σ A t (λ,ρ) and σ B t (λ,ρ), at µ A 1 | σ A t ,ρ,λ = 0 and µ B 1 | σ B t ,ρ,λ = 0, respectively (see Fig. 3b). Below these thresholds (σ < σ A t or σ < σ B t ), the spatially uniform solution is unstable and small perturbations in the velocity fields, α 1 or β 1 , have positive growth rates.
Close to the threshold values σ A t and σ B t , a weakly nonlinear analysis yields further insights into the dynamics of the system and the resulting steady states. Here we follow Bertin et al. [41] and use the following truncation scheme: We assume that close to the onset of polar order the velocity fields are small, α 1 , β 1 1, and obey the scaling relations (a − a (0) ) ∼ α 1 , α k ∼ α |k| 1 (and analogously for β k ). Furthermore, we assume that spatial and temporal variations ∇, ∂ t ∼ α 1 are small. With these assumptions, one can truncate and close the Boltzmann equation in Fourier space, Eq. (13), at the third order term for the velocity field, α 3 1 . We refer the interested reader to Appendix A and Refs. [41,42,49] for a more detailed discussion on the truncation scheme. Following the described procedure, we obtain hydrodynamic equa-tions for the velocity fields of species A and B: For explicit expressions of the coefficients please refer to Appendix A. The essential difference between Eqs. (18) and the hydrodynamic equation for the velocity field of one-species SPP [41] lies in the coefficients µ A 1 and µ B 1 , which here couple the velocity fields to the density of the respective other species (see Eqs. (17)). Note that the snowdrift interaction term only appears in the term linear in the velocity fields in Eqs. (18). This is a consequence of the truncation scheme and is explained in more detail in Appendix A.
Taken together, Eqs. (14) and Eqs. (18) constitute a minimal set of hydrodynamic equations that describe the dynamics of the slow variables of our system, the densities and velocity fields of both species. Overall, we identify four control parameters: the noise strength σ, the overall densityρ, the game strength τ , and the relative fitness λ. For specificity, we use a fixed value for the average density,ρ = 1 (in units of ω/(v 0 d 0 )) and vary the remaining parameters. We observe the same phenomenology if we instead fix σ and useρ as control parameter [52] (this is due to the opposing roles of density and noise changes in the self-organization of active matter systems [16,17,19]).
In the next sections we will first determine solutions for the spatially uniform case of Eqs. (14) and Eqs. (18) and subsequently analyse their linear stability against spatially non-uniform perturbations.

B. Spatially uniform solutions
The stationary uniform solutions of the hydrodynamic equations, Eqs. (14) and (18), can be found by setting all temporal and spatial derivatives to zero. Then, the equations for the densities and velocity fields decouple. For the densities, Eqs. (14), one obtains These are precisely the replicator equations, Eqs. (6), with the stable solution (note we have setρ = 1) Hence, the spatially uniform densities of each species are determined by the relative fitness λ, and are independent of the game strength τ .
The hydrodynamic equations for the velocity fields, Eqs. (18), of a spatially uniform system are given by Hence, α 1 = β 1 = 0 (corresponding to a disordered state) is always a trivial solution. A nontrivial solution occurs when the terms in the brackets balance each other. Since ξ > 0 for all values of the control parameters (see Appendix A), such a solution for α 1 (β 1 ) exists for µ A 1 > 0 (µ B 1 > 0). We know from the last section that the growth rate µ A 1 (µ B 1 ) is negative at high values of the noise strength σ and becomes positive below the threshold value σ A t (σ B t ), which depends on the relative fitness λ, as shown in Fig. 3b. Note that the threshold values for the noise strength are different for species A and B because their densities are different in steady state for games with λ = 1. The 'weaker' species (here B) has a lower density and hence is more prone to the disordering effect of noise, i.e. σ B t ≤ σ A t . Moreover, all results are independent of the game strength τ , since the terms involving τ in µ A 1 and µ B 1 (see Eqs. (17)) cancel as the densities reach their spatially uniform steady state solutions a (0) and b (0) given by Eqs. (20). The game strength τ will, however, play an important role for the dynamics of spatially non-uniform systems.
For values of the noise strength below the respective threshold, the solution α 1 = β 1 = 0 is linearly unstable and additional, stable solutions appear, given by These solutions correspond to collective states which spontaneously break the rotational symmetry of the system by choosing random directions of macroscopic order θ a and θ b , respectively. Figure 3a shows the absolute values (amplitudes) |α  (20)). In this case the growth rates µ A 1 (a (0) , b (0) ) and µ B 1 (a (0) , b (0) ) are identical for all values of the noise strength σ (see Eqs. (17)). Thus, the amplitudes of the velocity fields, |α  When σ is decreased below σ A t = σ B t , the amplitudes of the velocity fields become non-zero, which means that both species exhibit spatially uniform macroscopic order. For λ < 1, the amplitudes |α  1 | differ. Since a relative game strength λ < 1 indicates that species A is 'stronger' than species B, its stationary density (fixed point) is larger than that of species B: a (0) > b (0) . As a consequence it has the larger growth rate, µ A 1 > µ B 1 , and thus polar order of A is less prone to noise resulting in a The threshold values for the noise strengths coincide if the relative game strength λ = 1.
larger threshold value, σ A t > σ B t , see Fig. 3b and Fig. 3a for λ = 0.6 and λ = 0.2. This underlines an important aspect of the game interactions. These determine the relative density of the species A and B, which directly affects the ability of these species to establish macroscopic polar order, i.e. the threshold values of the noise.
In summary, for our two-species system with game interactions, the spatially uniform solutions show a new aspect of macroscopic order that is absent in one-species SPP systems [41,42,49]. There is not only one transition between a polar-ordered and a disordered state at some threshold value for the noise strength but two successive transitions. The underlying reason is that the competition of the two species, as defined by the (snowdrift) game, implies that their stationary densities are different: the species with the lower fitness also has a lower density. As the competition between alignment, favoring polar order, and noise, favoring disorder, crucially depends on the particle density, the noise thresholds for the ordering transition of the two species are different. The more competitive species has a higher density and, therefore, orders at a higher noise strength than the less competitive species. This is summarized by the phase diagram shown in Fig. 3b. By varying the noise strength σ and the relative fitness λ there are three phases: a disordered phase above σ A t , an intermediate phase for where only the more competitive species, A, shows polar order, and a phase where both species exhibit polar order below σ B t .

C. Stability against spatially non-uniform perturbations
In this section we will test the linear stability of the spatially uniform solutions given by Eqs. (20) and Eqs. (22), against perturbations with a finite wavelength. To this end, we introduce spatially non-uniform perturbations around the uniform solutions Assuming that the perturbations are small, we insert the perturbed fields in the full hydrodynamic equations, Eq. (14) and Eq. (18), and keep only terms up to linear order in the perturbation. Note that one has to account for both perturbations δα 1 , δβ 1 and its complex conjugates δα * 1 , δβ * 1 . Introducing Fourier modes δa(r, t) = 1 2π ∞ −∞ dq e iq·r δa(q, t) and analogously for all other fields, one obtains a set of linear equations for with the Jacobian J q . It is a 6×6 matrix that depends on the wave vector q and the spatially uniform solution 1 ) and therefore on all control parameters (λ, τ, σ). An explicit expression for J q can be found in Appendix B. For simplification, we assume that the spatially uniform solutions of the velocity fields are parallel, α 1 . Numerical solutions of the hydrodynamic equations and agent-based simulations in the next sections will show that this assumption is justified.
The real parts of the eigenvalues of the Jacobian J q determine the linear growth rate of the perturbations: An eigenvalue with a positive real part indicates that the corresponding uniform solution is linearly unstable and the perturbation given by the respective eigenvector of the Jacobian will grow (exponentially) in time. Figure 4a shows examples of the real part of the eigenvalue with the largest real part, s q , as a function of q for different values of σ and fixed λ = 0.4 and τ = 0.2. For certain values of the noise strength σ there are long wavelength instabilities (i.e. for small values of the wave vector q) .
Our linear stability analysis shows that the spatially uniform ordered states are unstable with respect to spatial perturbations right at the onsets of both transitions to macroscopic order, σ A t and σ B t (Fig. 4b). In other words, the disordered phase (σ > σ A t ) is stable, while the phases with spatially uniform polar order in species A (σ A t > σ > σ B t ) and uniform polar order in both species (σ B t > σ) are both unstable towards spatial perturbations for σ close enough to the respective onset of or-  14) and (18). (a) Dispersion relation of the eigenvalue with the largest real part, sq, for different values of the noise strength σ, and fixed game strength τ = 0.2 and relative fitness λ = 0.4. sq is plotted for wave vector q pointing in the direction parallel to macroscopic order, q α1, β1 (denoted by q ). The spatially uniform state, Eqs. (22), is unstable when the maximum of the dispersion relation, smax, is larger zero. (b) Parameter regime (bands) of λ (relative fitness) and σ (noise strength) where the spatially uniform polar-ordered states are linearly unstable against spatial perturbations (i.e. smax > 0), for different values of the game strength τ (shades of blue) as indicated in the graph. The two bands are bounded from above by the threshold values σ A t and σ B t indicating the onset of spatially uniform polar order for species A and B, respectively (dashed lines). The width of the bands decreases with increasing game strength τ . White regions indicate parameter regimes where the spatially uniform ordered states are stable against spatial perturbations. (c), (d) Linearly unstable regions in (τ, σ)-space and (λ, σ)-space, respectively. The color code compares the two density field components (δa and δb) of the eigenvector corresponding to the eigenvalue with the largest positive real part: δa := δa/(δa + δb) which ranges from 0 (the largest linear growth rate corresponds to an eigenvector with δa-component equal to zero and a non-zero δb-component) to 1 (δb-component equal to zero and non-zero δa-component).
der. By computing the real part of the largest eigenvalue, s q , for different wave vectors q, we find that s q is largest for wave vectors parallel to the direction of macroscopic order (q α 1 ); i.e. instabilities are strongest when longitudinal. This suggests that the ensuing patterns will form along the direction parallel to macroscopic order, α (0) 1 and β (0) 1 , while the system remains homoge-neous along the direction perpendicular to macroscopic order. For τ = 0, i.e. in the absence of snowdrift interactions, we recover the same unstable regions in parameter space as previously obtained from one-species SPP systems [41,42] (see Fig. 4b). For a finite game strength, τ > 0, the upper boundaries of the unstable regions in parameter space are still σ A t and σ B t (as in the case of no snowdrift interaction). The width of the unstable regions, however, is smaller indicating that the snowdrift interaction has a stabilizing effect on the spatially uniform solutions.
The perturbation with the largest linear growth rate is given by the eigenvector δm q corresponding to the eigenvalue s q with the largest positive real part (and with wave vector q fixed such that s q is maximal); see Eq. (24). We compute this eigenvector as a function of the three control parameters σ, λ and τ in order to see how the interplay of snowdrift game and self-propelled motion affects which perturbations have the largest growth rate. In order to compare the growth rates of the fields of the two species, we first look at the two density field components of the eigenvector, δa and δb. We use the expression δ a := δa/(δa + δb) to measure the direction of the eigenvector projected onto the plane spanned by the two density fields. This value can range from 0 (δa = 0 and non-zero δb-component), meaning that only perturbations in the density of species B grow, to 1 (δb = 0 and non-zero δa-component) meaning that only the density of species A is unstable towards perturbations. Fig. 4c and Fig. 4d show the value of δ a as a function of the control parameters in (τ ,σ)-and (λ,σ)-space, respectively. Analysing the results in (τ ,σ)-space (Fig. 4c), we can read off how the game strength, τ , affects the relative growth rates of perturbations of the two densities, a and b: For large enough game strength, τ , perturbations in the densities of both species grow equally fast (δ a → 0.5). Decreasing τ changes the value of δ a dependent on the region in parameter space: In the unstable region with upper boundary σ = σ A t , perturbations in the density of species A are growing faster than perturbations in the density of species B such that δ a → 1 for small enough game strength τ (see Fig. 4c). In other words, the perturbation with the largest growth rate changes from an even orientation in the directions of both types to an orientation only in the a-direction. In the unstable region with upper boundary σ = σ B t , δ a → 0 for small enough τ , which means that only the density of species B is unstable. The inset of Fig. 4d shows the limiting case τ = 0 in the (λ, σ)-plane. From the color code we read off that, depending on the region, either δb = 0 (yellow region) or δa = 0 (blue region). In other words, for two non-interacting species (τ = 0), perturbations only grow in either species A or B. Note that in this case, instabilities in the density fields for each species occur at the respective threshold to macroscopic order: perturbations δa grow at σ < σ A t ; perturbations δb grow at σ < σ B t . There are four more components to the eigenvector, encoding the perturbations of the velocity fields. An analogous analysis as above for the velocity field components parallel to macroscopic order, δα 1 and δβ 1 , yields the same phenomenology as for the densities: For high values of τ , perturbations in the velocity fields of both species grow equally fast (δα 1 /(δα 1 + δβ 1 ) → 0.5). For τ → 0, on the other hand, perturbations δα 1 only grow in the unstable region in the vicinity of σ A t , while perturbations δβ 1 only grow in the unstable region at σ B t . These results are discussed in more depth in Appendix B 1. Components perpendicular to the macroscopic order are stable to perturbations.
In summary, the main insight gained in this section is that there are regions in parameter space where spatially uniform polar ordered states are linearly unstable against spatially non-uniform perturbations. These regions are located in the immediate vicinity of the threshold values for the noise strength, σ A t and σ B t , that mark the instability of a disordered state to states with different degrees of spatially uniform polar order. These threshold values of the noise strength are functions of the relative fitness, λ. The game strength, τ , affects the relative growth rates of perturbations of the two densities, a and b, and velocities parallel to macroscopic order, α 1 and β 1 : For high values of τ , perturbations in the fields of both species grow equally fast; for decreasing τ , growth rates differ between the species until for τ → 0 instabilities in the fields of A and B occur independently. In the next section we will study the spatial patterns that emerge in these unstable regions.

D. Numerical solutions of the hydrodynamic equations
The spatially uniform solutions, Eqs. (20) and (22), and their linear stability analysis (Fig. 4) suggest five distinct phases. Especially in the regimes where spatially uniform ordered states are unstable against spatially nonuniform perturbations (coloured regimes in Fig. 4b,c,d), we expect the formation of spatio-temporal patterns. In order to study the nonlinear dynamics of the hydrodynamic equations, Eqs. (14) and (18), beyond linear stability, we solve these equations numerically. Specifically, we use a finite difference scheme for a two-dimensional grid with periodic boundary conditions. We employ an explicit iterative process, the 4th order Runge-Kutta method, to determine the time evolution; for details on the implementation please see Appendix C.

Phase diagram and phase transitions
These numerical solutions of the hydrodynamic equations confirm the existence of five distinct phases (Fig. 5). We find three spatially uniform phases, a disordered phase, an intermediate phase where species A shows polar order while species B does not, and a phase where both species exhibit polar order. The parameter regimes determined numerically are in accordance with our analytical calculations (described in sections III B and III C), where we found uniform solutions to be linearly stable against non-uniform perturbations (the three white regions in Fig. 4b,c,d).
In the parameter regime where the spatially uniform solutions are linearly unstable against spatial perturbations, we observe the formation of spatial patterns. These patterns are characterized by density-segregated, polarordered bands in both species that concurrently move through the system; see Movie movie1 in the supplementary material for an illustration (an explanation and the parameter set for the movie can be found in Appendix E). We will refer to these traveling bands as polar waves. The wave vector of the traveling patterns always points in the same direction as the macroscopic polar order, as suggested by linear stability analysis (see Section III C). When viewed from a co-moving frame, the patterns are uniform in the direction perpendicular to the wave vector. Similar density-segregated traveling wave states are known from one-species SPP systems [42,53,54]. Interestingly, the domains where traveling wave patterns exist, extends beyond the parameter regime where the spatially uniform polar-ordered states are unstable to spatial perturbations of finite wavelength. In the blue shaded regime in Fig. 5 the system exhibits bistabilty and subcritcal behavior as discussed in detail below; see section III D 2. We have also investigated phase diagrams for different values of τ and found the same topology.
We will see in the next section that the mechanism underlying pattern formation varies according to the region in the parameter space. A distinction is therefore made between two different phases with spatial patterns. We refer to the traveling polar wave state emerging in the unstable region in parameter space at σ A t as traveling polar wave phase A (TWA) and the one in the unstable region at σ B t as traveling polar wave phase B (TWB). For both of these phases, examples of the wave profiles for the density fields, a and b, and the velocity fields, |α 1 | and |β 1 |, are shown in Fig. 6; the superscript for the velocity fields indicates that we consider the component parallel to the direction of macroscopic order, α 1 and β (0) 1 , and wave vector, q (as mentioned in section III C, the stability analysis showed that instabilities are largest for q α 1 , meaning that patterns form along the direction of macroscopic order). Our numerical solutions show that the polar waves for species A and B are strongly correlated spatially, meaning that they follow each other through the system. These findings are consistent with our results from linear stability analysis, where we found that instabilities always occur simultaneously for interacting species A and B.
In the three spatially uniform phases, all fields show values that are close to their uniform fixed point values given by Eqs. (20) and (22). Although in the TWA and TWB phase, the fields exhibit spatial patterns and deviate significantly from their uniform fixed point values at different points in space, the spatial averagesā,b, |ᾱ 1 | and |β 1 | are approximately also given by the corresponding uniform fixed point values.
Moreover, we find that at every point in space, the velocity fields, |α 1 (r, t)| and |β 1 (r, t)|, are well approximated by the local values of the spatially uniform solutions, Eq. (22). For local densities a(r, t) and b(r, t) with µ A/B 1 > 0, these are given by |α , else we have |α 1 |, |β 1 | = 0. Hence, non-zero values for the velocity fields α 1 and β 1 requires that the densities are large enough such that locally one is above the threshold for polar order, i.e. both µ A 1 (a, b) and µ B 1 (a, b) are positive (see section III B). The only exception are patterns of the B species in the TWA phase in a regime where both λ and τ are small enough such that the density b falls below the threshold value for polar order, µ B 1 (a, b). Nevertheless we see order emerge in the high-density wave peak (see lower inset of Fig. 6a). In summary, the pattern can segregate the system in regions of high density with macroscopic order (|α 1 | > 0 or |β 1 | > 0) and low density regions that are disordered ((|α 1 | = 0 or |β 1 | = 0)).
corresponding to spatially uniform solutions (see Eqs. (22)) for the respective local densities a(r, t) and b(r, t) (green solid lines in the insets). The waves in species B in the TWA phase seem to be an exception to this observation (inset in the bottom left figure).

Bistability and subcriticality
Previous studies on SPP systems indicate that patterns are stable even in regimes where uniform solutions are linearly stable [47,49]. To identify the full regime in which pattern solutions exist, we use the following approach ('hysteresis loop'): We first compute the stationary solution for fixed parameters (λ, τ ) and σ > σ A t where there is no macroscopic order, neither in A nor in B. Then, we reduce the noise strength σ quasi-statically in small steps, thereby sweeping through the different phases in parameter space. By 'quasi-static' we mean that the system is given enough time to equilibrate between successive adjustments in σ. Afterwards, we increase the noise strength again quasi-statically back to its initial value. During this σ-sweep, we monitor the spatial average of the velocity fields, |ᾱ 1 | and |β 1 |. To identify the presence of spatial patterns, we use parameters that we call variation parameters η A and η B , which measure spatial variations in the densities a and b. They are defined as where |∆a| ∆x is the absolute value of the numerical spatial derivative of the density field a and the bar indicates the average over all grid points (see App. C). For η A = 0 (η B = 0), species A (B) is in a spatially uniform state, while a non-zero value η A > 0 (η B > 0) indicates that the system exhibits some kind of spatial patterns in A (B). By combining the information about average velocities and spatial variations we can infer the phase of the system, as we discuss next. Figure 7 shows the spatial variation parameters, η A and η B , and spatial averages of the velocity fields, |ᾱ 1 | and |β 1 |, obtained from σ-sweeps for different values of the relative fitness λ and a fixed game strength τ = 0.2. A general observation from our numerical analysis is that phases with traveling polar waves exist in parameter regimes wider than the domains predicted by the linear stability analysis of spatially uniform states; compare the yellow-and blue-shaded regimes in Fig. 7. While in the yellow-shaded parameter regime only the traveling waves are stable solutions, in the blue-shaded parameter regime both these waves and spatially uniform states are stable solutions, i.e. the system shows a bistable dynamics and subcritical behavior.
For a neutral game (λ = 1), species A and B show identical behavior (Fig. 7a): Quasi-statically decreasing (grey arrows) and then increasing (black arrows) the noise strength σ one observes a hysteresis loop. Both, the variation parameters, η A and η B , and the average velocities, |ᾱ 1 | and |β 1 |, show discontinuous changes (up and down arrows) at the boundaries of the bistable regime at high noise values (right blue-shaded region). Hence, in the bistable regime, the disordered phase and the traveling wave phase (TWA) are possible metastable states. In contrast, in the bistable regime at low noise strength (left blue-shaded region) only the spatial variation parameter exhibits a hysteresis loop while the average velocities change continuously. Moreover, the spatial averages of the numerically determined values for the velocities closely follow the results calculated analytically for a spatially uniform system (dashed lines in Fig. 7a)); compare Eq. (22). This has two implications. First, a spatially uniform polar-ordered state and a traveling wave state are metastable solutions (in the left bistable regime). Second, the spatially uniform and spatially non-uniform state have the same average velocities. We therefore call this transition quasi-continuous.
For small enough values of the relative fitness λ (λ = 0.2 in Fig. 7b), one observes all of the five different phases (compare also Fig. 5) with bistable regions as well as hysteresis loops at each transition. In addition to the disordered phase at high noise strength and the spatially uniform polar-ordered phase at low noise strength, in which both species show polar order (|α 1 |, |β 1 | > 0), there is now an intermediate phase, in which only species A shows polar order, but species B does not (|α 1 | > 0, |β 1 | = 0). Similar as for a neutral game, the spatial mean of both velocities is fairly well approximated by the analytical results we obtained for a spatially uniform system (dashed lines in Fig. 7b), even in the regime where the system actually shows phase separation into polar bands of high density and disordered spatial regions of low density.
For large enough values of the relative fitness λ (λ = 0.5 in Fig. 7c), the two traveling polar wave phases, TWA and TWB, merge; compare also Fig. 5. The phase behavior is similar as in the neutral case with the difference, that the travelling polar waves in species A and B differ in shape and degree of polarization (see the different values of η A , η B and |α 1 |, |β 1 | in Fig. 7c).
Last, in a small parameter regime (close vicinity around the blue star in Fig. 5), the transition between the uniform partially polar-ordered phase and the TWA phase becomes supercritical, i.e. we do not find jumps in the variation parameters and average velocities, but a continuous change between zero and non-zero values. Moreover, in this regime the dynamics become extremely slow, so that the system takes a long time to reach its spatially non-uniform, stationary state, indicating critical slowing down.

Game-induced pattern formation
In the last section we found that the two traveling wave phases, TWA and TWB, exhibit coupled traveling polar wave patterns. In a next step, we now investigate the relationship between the interaction mediated by the snowdrift game and the formation and coupling of polar wave patterns. This will shed light on the differences between the TWA and TWB phases and on the underlying mechanism of pattern formation.
The linear stability analysis performed in Section III C already indicated a possible role of the game interaction in pattern formation. There we found that for finite game strength (τ > 0) spatial perturbations in the density and velocity fields of species A and B grow simultaneously, while in the absence of game interactions (τ = 0) instabilities in the fields of A and B occur independently. This suggests that the games interaction is responsible for the coupling between the polar wave patterns of species A and B.
To pinpoint this idea, we systematically study the formation of polar wave patterns as a function of the game strength τ . Figure 8 shows profiles of the densities, a and b, and the velocity fields, |α 1 | and |β 1 |, for different values of the game strength τ . The other two control parameters are set to σ = 0.68 and λ = 0.4, such that the system is in the TWA phase. For τ = 1, traveling polar waves emerge in both species, as expected (Fig. 8a). For smaller values of the game strength τ , we observe a decrease of wave amplitudes in b and |β 1 | but an increase of wave amplitudes in a and |α 1 | (Fig. 8b). In general, we find that the smaller the value of τ , the lower the amplitude of waves in species B and the higher the amplitudes in species A. In the absence of a snowdrift interaction (τ = 0), the wave pattern in species B completely disappears and species B is in a disordered, spatially uniform state with b = λ 1+λ ≈ 0.29 and |β 1 | = 0 throughout the system (Fig. 8c). We conclude that the pattern in B observed for τ > 0 in the TWA phase is due to the snowdrift interaction. This behaviour is consistent with our analysis of the direction of the linear instability as a function of τ in Section III C (see Fig. 4c). There we found that the b-and β 1 -components of the eigenvector pointing in the direction of the instability gradually decreased with the game strength τ , such that the a-and α 1 -components remained as the only non-zero contributions for τ → 0.
Our interpretation of these observations is the following. The alignment interactions between the particles of species A cause polar wave patterns in the TWA phase. This is consistent with previous studies on SPP systems, that find traveling polar waves at the onset of ordered motion.
For a finite game strength, τ > 0, the game interaction induces a shift of the density of the B species towards the value b = λ a, given by the fixed point of the snowdrift game (Eq. (20)). The magnitude of this shift (flow in the phase space of the game) depends the game strength τ . The stronger the game interaction, the stronger the mutual influence of the two species. For large τ , one expects an adiabatic limit where a and b are strictly correlated by b/a = λ; compare Figs. 8a-c. We call this mechanism game-induced pattern formation.
There is an analogue mechanism in the TWB phase. Our simulations show that with decreasing game strength τ , the pattern amplitudes in a and |α 1 | decrease while they increase in b and |β 1 |. For vanishing game strength (τ → 0), patterns in species A fully disappear and one is left with spatially uniform polar order in species A and polar wave patterns in species B. Hence, in the TWB phase, species B induces the traveling wave pattern in species A.
In parameter regions where the TWA and TWB phases overlap, both species form polar waves independently, due to alignment interactions. Game interactions cause synchronization of these wave patterns. Interestingly, when we tune the game strength τ to very small values, in a parameter regime where the densities of both species are high enough to exhibit pattern formation, we observe "oscillating" solutions where a faster moving traveling wave periodically passes a slower wave. Our explanation for this phenomenon is that the active properties responsible for the wave formation dictate different wave speeds for both species and the game strength is too weak to synchronize both patterns.

IV. AGENT-BASED SIMULATION
As a complementary approach and to test the robustness of the results obtained with the kinetic Boltzmann approach, we implemented agent-based simulations of a system of two species that exhibit intra-species alignment and inter-species interactions mediated by the snowdrift game. For a single species, previous agent-based simulations showed that aligning collisions can lead to the formation of macroscopic order when the density of particles is high enough or orientational noise is weak [42,55,56]. In addition, at the onset of macroscopic order, particles were observed to form polar-ordered high-density waves rather than spatially uniform polar order [53,56].
A. Agent-based model Following Bertin et al. [42], we use a generalized Vicsek model [55] for self-propelled particles and extend it to two species, which interact through a snowdrift game. We consider a system of N particles on a two-dimensional square plane of area L × L with periodic boundary conditions. Each particle, located at position r i , is either of species A or B, and moves with a velocity v i = v 0 e i θ ; the speed v 0 is assumed to be constant and the direction vector reads e i θ = (cos θ i , sin θ i ), where θ i denotes the angle enclosed with the x-axis.
The positions of all particles are updated in discrete time steps dt Particles move ballistically unless they change their direction e i θ in random 'diffusion' or 'collision' events. A 'collision' event is defined in the spirit of a Vicsek model [55]: Each particle aligns its direction to the average direction of all particles of the same species in its alignment vicinity A t j at time t. A t j is defined as the disk around r t j with radius d align , the alignment range, corresponding to d 0 in the kinetic model. Independent of alignment, each particle is subject to 'diffusion' events, in which its direction changes by a random amount η. We assume that η is a random variable that is uncorrelated in time and space and given by a flat distribution over the interval (−hπ, hπ] with h ∈ [0, 1]. The standard deviation σ of a flat distribution is then given by σ = hπ/ √ 3. We will refer to σ as the noise strength. The direction of each particle at time t + dt is therefore given by Previous studies have shown that if the noise strength σ is sufficiently low, species begin to form macroscopic polar order, i.e. particles move collectively over distances much larger than their alignment range d align [55,56] In order to measure polar order of each species we define the polar order parameters e iθ t k , (28) where N t A and N t B denote the total particle number of species A and B, respectively, at time t; α t 1 and β t 1 range from 0 (no order) to 1 (perfect alignment of all particles). The temporal averages of these polar order parameters are given by In addition to intra-species alignment interactions, particles also play an inter-species snowdrift game each time they 'encounter' each other. Analogously to the alignment vicinity, we define a snowdrift game vicinity G t j of particle j at time t by a disk around particle j with radius d game , which we will refer to in the following as the game range.
All particles N t j within this range are assumed to participate in the game, similar as in an urn model [6]. The corresponding particle numbers for species A and B are denoted by N t A,j and N t B,j , respectively. Consider the case that particle j belongs to species A: This A-particle will be at a disadvantage in the competition with the Bparticles if the relative abundance, N t B,j /N t A,j , lies below a threshold λ, which defines their relative fitness. We assume, in analogy to the replicator equations Eqs. (6), that a disadvantaged particle will switch species, i.e. change strategy (here A → B) with a rate τ (a t j λ − b t j ). Here, a t j = N t A,j /N t j and b t j = N t B,j /N t j denote the particle densities of species A and B, respectively (within the game range of particle j). On the other hand, if a particle is at an advantage in the game it plays, i.e. the relative abundance is greater than the threshold λ, we assume that it keeps its strategy. An analogous argument holds if j happens to a B-particle.
Taken together, we define the probability for switching species as As in Sec. II, τ ≥ 0 is interpreted as the game strength: With increasing τ the strength of the snowdrift interaction compared to the alignment interaction increases. Because the species are antisymmetric under relabeling A and B, we can restrict the relative fitness to λ ∈ [0, 1], i.e. species A is the 'fitter' species. We assume that particles which change species lose their directional information and move in a random direction afterwards until they align with partners of their new species. This assumption, which we refer to as 'switch randomization', prevents transmission of directional information between species, such that all inter-species interactions are purely due to the snowdrift game.
To be able to identify spatial structure as observed in previous agent-based SPP models [42,55,56] we divide the system into C channels perpendicular to the direction of macroscopic order. Each channel has area L 2 /C. We then calculate the densities a t c and b t c and polar order parameters α t 1,c and β t 1,c for each channel c. The same partitioning is made parallel to the direction of macroscopic order. If waves are present in the system, we expect a difference in the magnitude of the order parameter for different channels, whereas for an uniformly ordered phase the order parameters for each of the channels should be almost identical. We measure this difference with the velocity variation parameter ν t A , defined as twice the standard deviation of polar order parameter α t 1,c over all channels: The factor of 2 is introduced such that ν t A ∈ [0, 1]. The velocity variation ν t B for species B is defined analogously. To see if waves are present over an extended period of time, we use the differences between the time averages ∆ν A = ν A,⊥ − ν A, , with channels perpendicular and parallel to the direction of macroscopic order, respectively. Further details can be found in Appendix D.

B. Results of the agent-based simulation
Of the eight model parameters (λ, τ , ρ, σ, v 0 , d align , d game , and L) we focus mainly on analyzing the role of two: (i) the noise strength σ, and (ii) the relative fitness λ. The noise strength σ is used as a control parameter to distinguish between the different phases, with the density kept fixed at ρ = N/L 2 = 1.5. Varying the relative fitness λ tunes the ratio of species B to A. The remaining parameter values are listed in Appendix D. We have tested that the same phenomenology found for these specific parameter values is also observed in a broader parameter range, similar to the kinetic Boltzmann approach [57].
In our agent-based simulations for the 'active matter game' we observe a similar phenomenology as in the kinetic Boltzmann approach: Depending on the parameters, especially the game strength τ , relative fitness λ, and noise strength σ, we observe up to five different phases. For neutral games (λ = 1), we observe, similar to previous studies [42,55,56], three different phases: a disordered, a traveling polar wave, and a uniformly polar-ordered phase. Figure 9 shows the diagram for a game strength τ = 0.4, and two distinct game ranges d game ∈ {5, 10}. The topology of the diagram is similar to the one obtained from the hydrodynamic equations (Fig. 5), and it is independent of the specific value for Grey shaded areas indicate parameter regimes where the system is in either a disordered or uniformly polar-ordered state, as indicated in the figure legend. For all transitions between phases there are transition regions in which both phasesa uniform phase (either ordered or disordered, depending on the transition) and a polar wave pattern -can dynamically change (light blue regions). Note that, as explained at the end of Section IV B, due to demixing of species inside the waves, the relative size of the parameter regime covered by the TWA and merged wave phase widens with increasing dgame. The vertical orange lines in the right panel indicate the three slices in parameter space displayed in Fig. 10. the game range [58]. In the following discussion, the specific parameter values given refer to the case d game = 10.
The agent-based simulations show that above some threshold value for the noise strength (σ 0.7) both species are in a spatially uniform, disordered state, with the polar order parameters, α 1 and β 1 , and the velocity variation ∆ν for both species approximately zero; see upper-most light grey area in Fig. 9 and Fig. 10, as well as Fig. 14a. Below this threshold noise value, species A shows polar waves independent of the value for λ, as can be inferred from Fig. 10a where α 1 and ∆ν A are shown as functions of σ. In analogy to the corresponding phase observed in the kinetic Boltzmann approach, we call this phase traveling wave phase A (TWA). How these polar waves in species A affect the degree of order in species B depends on the strength of the relative fitness λ. Within the polar waves, the density of A is enhanced, which in turn implies an increased density for species B due to the game interaction; recall that at the mean-field level we have b = λ · a. Whether this induces polar order depends on whether the density of species B within the band exceeds the threshold for polar order (in the absence of the game). Indeed, we find that for low values of the relative fitness, λ 0.26, there is only a density wave in species B but no (or only weak) polar order (see Fig. 10a and Fig. 14b). For larger values of λ (Fig. 10b,c), however, the density of species B becomes large enough within the polar bands of species A such that the threshold for polar order is crossed and species B exhibits polar bands as well. Since this polar order is mediated by the polar waves of species A through the game interaction and not through the alignment interaction of species B per se, it is a game-induced phenomenon similar as in the hydrodynamic approach (Sec. III D 3). A movie displaying the dynamics in the TWA phase can be found in the supplementary material (movie2).
Upon further decreasing the noise strength σ, there is a parameter domain for λ 0.29 where we find an intermediate phase (see the grey-shaded region in Fig. 9 enclosed between the TWA and TWB phases). In this region of the diagram, there are no longer any wave patterns in species A or B; see central white regions in Figs. 10a,b. Furthermore, while the polar order of species A is maintained (α 1 > 0), species B no longer exhibits polar order (β 1 ≈ 0). In other words, species A is now in a uniformly polar-ordered state, while species B is in a disordered state; see also Fig. 14c.
Upon decreasing the noise strength σ even further, we observe wave patterns again. Now, both species A and B display traveling waves with polar order; see Figs. 10a,b and Fig. 14d. Species B forms polar waves by itself due to alignment interactions. The polar waves in species A, however, have a different origin: Without the game interaction species A would already be in a uniform polarordered phase. The role of the interaction mediated by the snowdrift game is to spatially modulate the density of species A such that it mirrors the wave pattern of species B; for spatial profiles of the density and the polar order see Fig. 11b. We call this phase traveling wave phase B (TWB); the dynamics are displayed in movie3. At even lower noise strength, the waves disappear again, and both species exhibit uniform polar order (dark grey area in Fig. 9 and movie4). Transitions between different phases do not occur at exact parameter values, instead there is a region with stochastic switching, where the system alternates between the two adjacent phases.
While the diagrams obtained from the agent-based simulations (Fig. 9) and using the hydrodynamic theory (Fig. 5) are qualitatively fairly similar, the agentbased simulation reveal interesting additional phenomena which are absent in the Boltzmann approach. For phases in which both species show long-range polar order, i.e. polar waves or uniform polar order, we observe a demixing of species A and B, meaning clusters or alternating bands of only one species or the other species. Figure 11a illustrates this demixing in the TWB phase for polar waves parallel to the direction of wave propagation. Note that the bands are wider in species A , as it has a fitness advantage compared to species B (λ = 0.26).
The size of these demixed regions grows linearly with increasing game range d game , as is illustrated in Fig. 15.  λ (a-c), indicated in the graph. There are parameter regimes with spatially uniform states (white regions, corresponding to grey shaded areas in Fig. 9), wave patterns (yellow shaded region), and bistable behavior where both wave patterns and uniform (dis-)order (blue shaded region) stochastically alternate. We observe three qualitatively different regimes. For λ = 0.23 (a) there are distinct TWA and TWB phases, separated by an intermediate phase. The TWA phase shows non-zero polar order for species A (α1 > 0) and finite velocity variations (∆νA > 0), while there is neither significant polar order in species B (β1 ≈ 0) nor velocity variations ∆νB. In contrast, the TWB phase is characterized by polar order as well as velocity variation in both species A and B. The case λ = 0.28 (b) is phenomenologically similar to (a), but in the TWA phase there is now also polar order (β1 > 0) and finite velocity variations (∆νB > 0) for species B. The intermediate phase in which species A is uniformly ordered and species B disordered (center white strip) is much smaller than in (a). For λ = 0.31 (c), the TWA and TWB phases merge to one broad, single travelling wave phase which exhibits polar order and finite velocity variation in both species.
If demixing occurs, the ensuing density ratio b/a becomes larger than the relative fitness λ. Both demixing and the density shift occur only if the relative motion between particles is slow (e.g. in the polar-ordered phases) compared to the rate τ of snowdrift game interactions. In the agent-based simulations used in this study, this requirement is met only for the traveling wave phase induced by species B (TWB) and the uniformly polar ordered phase, the only phases in which both species A and B display sufficient polar order. While the demixing manifests itself as bands in the TWB phase, the species form clusters in the uniformly polar ordered phase; see Fig. 14e. The observed demixing of species is an interesting feature of the snowdrift game and occurs even in the absence of alignment interactions (see left panel of Fig. 15), as long as the diffusive mixing of the system is sufficiently slow.
We attribute the formation of such clusters and bands to the finite interaction range of the snowdrift game. In a well-mixed population the density ratio is given by b/a = λ and, hence, all particles have equal fitness. Upon demixing within the interaction range of the game, however, the overall fitness can be increased as alike species cluster together. This tendency for demixing dominates Corresponding wave profile as a function of position x averaged over the y-direction. For the analysis we partitioned the system into 256 channels along the x-axis, counted the particles per channel and evaluated the channel densities a, b, and ρ (upper plot) and polar order α1, β1 (lower plot). There is a clear correlation between the density wave and polar order; outside the wave the system shows a low-density disordered phase.
over 'entropic' mixing effects due to random particle motion, i.e. for slow movement of particles relative to each other (in comparison to the game strength).
In the agent-based approach we observe that the TWB phase is stable over a much wider range of noise strength σ than in the hydrodynamic approach, and becomes even larger with increasing relative fitness levels λ; cf. Fig. 9 and Fig. 5. We hypothesize that this stabilization of travelling waves is due to the interplay between demixing and the loss of directional information during switch randomization: Particles that change species at the rear end of a wave and thereby lose their directional information counteract the formation of uniform polar order.

V. DISCUSSION AND CONCLUSIONS
We investigated an active matter system consisting of two different types of self-propelled particles interacting both mechanically and through birth-death processes described by a game. For the mechanical interaction we used standard alignment of the particles' velocities implemented either as pairwise collisions with half-angle alignment (following Bertin et al. [41]) or alignment in a certain neighborhood (in the spirit of a Vicsek model [55]). In addition, we considered competition between the different particle types ('species') mediated through a snowdrift game, again played either upon collision or in some neighborhood. Since these interactions can change the (local) density of particles -a critical parameter for symmetry-breaking and pattern-formation in active matter -they are an essential and interesting extension of classical active matter models. We would even suggest that they may be seen as a simple classic example of a more general class of interactions that affect the density of the self-propelled particles rather than their orientation.
We have derived hydrodynamic equations for the densities and velocities of the particles, starting from a kinetic Boltzmann equation for polar particles that we extended by the replicator dynamics of the snowdrift game. Using linear stability analysis we have studied the system's dynamics as a function of two key control parameters, the noise amplitude of the polar alignment interactions and the relative fitness of the two species in the snowdrift game. These studies revealed different regimes of macroscopic polar order and pattern formation for one or both species. Below a critical noise strength, we observe that particles of the species with the higher fitness form polar order and evolve into polar traveling wave patterns. Interestingly, these patterns-mediated by snow drift interactions between the two species-induce polar travelling waves in the other species as well. Furthermore, we find that in a regime of even lower noise values, the species with the lower fitness forms polar wave patterns, and these patterns again induce polar traveling waves in the other species. We referred to this phenomenon as game-induced pattern formation. Although we have discovered this phenomenon for the specific case of polar active particles interacting through a snowdrift game, we would argue that it is a fairly general phenomenon that applies to a broad class of active systems where the local density of particles is affected by (nonlinear) reaction kinetics. We have complemented the kinetic approach with agent-based simulations in which interactions between particles take place in some extended neighbourhood rather than in collisions. These simulations qualitatively confirm the phase diagram and the various regimes of macroscopic order and pattern formation, including the game-induced patterns. This shows that our findings are robust and do not depend on the particular implementation of the active matter system. Moreover, our numerical simulations of the hydrodynamic theory and the agent-based model consistently show that the transitions between the different types of collective behavior are characterised by discontinuous changes in the velocity and spatial density variation. In particular, they also exhibit hysteresis and bistability. This shows that the ensuing phase transitions are generically first order and the formation of patterns is subcritical. Only by fine-tuning parameters could we find a regime where the transition becomes continuous. Finally, our agent-based simulations show the formation of species separated clusters (disks and bands) in the ordered phases due to the finite range of game interactions. The simulations suggest that these spatial structures can lead to a stabilization of game-induced wave pattern leading to a broader parameter regime of these waves.
From a broader perspective, this study suggests a possibly general mechanism in active matter systems where patterns in one species induce patterns in another species.
Here, since we focused on two species with qualitatively identical pattern phenomenology (traveling waves), the patterns induced by game interactions are qualitatively the same (traveling waves). More generally, it would be interesting to ask what could happen if the individual species showed qualitatively different types of collective behavior and different pattern phenomenology. For instance, different symmetries of intra-species alignment, e.g. due to different shapes of the species' agents, could lead to density regimes of patterns with different symmetries for the two species (e.g. polar clusters and nematic patterns) [59]. Based on our study, patterns in one species induce density variations in the second species (game-induced patterns). These density variations in the second species could then locally enter regimes of pattern formation and induce patterns in the second species. As a result there would be patterns of different symmetries, which undergo mutual feedbacks and possibly lead to a rich pattern phenomenology with novel collective states in active matter.
Our study could also, perhaps only on a conceptual level, relate directly to 'real life' systems, where motile organisms undergo chemical interactions with each other. Examples for biological systems whose spatiotemporal dynamics are shaped by the interplay of motility and chemical interactions include motile bacteria that compete for shared resources [60] and bacterial systems that interact through chemotaxis and other types of extracellular signals. For a direct comparison with these systems, it will certainly be necessary to consider a more complex interaction scenario than the snowdrift game considered here and, for example, explicitly include the spatiotemporal dynamics of the extracellular chemicals and more complex chemical reaction kinetics [1]. Other systems, where we expect similar phenomenology, may include systems with different chemical species engineered such that these species are self-propelled and show chemical reactions.

ACKNOWLEDGMENTS
We would like to thank Lorenz Huber for stimulating discussions, and Fridtjof Brauns and Laeschkir Würthner for their help in implementing the numerical simulations. This research was supported by the Excellence Cluster ORIGINS, which is funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy EXC-2094-390783311, and the funding initiative "What is life"of the Volkswagenstiftung.
Appendix A: Derivation of the hydrodynamic equations The kinetic Boltzmann Eqs. (9) set the mathematical framework for the investigation of our active particle system with snowdrift game interaction. Following [42,49], we use these equations as a starting point to derive hydrodynamic equations that describe the dynamics of the density and polar order at least in a parameter regime close to the onset of macroscopic polar order.

Extended kinetic Boltzmann equations
As discussed in Section II, we model the dynamics of two species A and B with self-propulsion and selfdiffusion, that undergo intraspecies polar alignment, and interspecies game interactions. All agents are assumed disc-shaped with diameter d 0 and velocity v 0 . We therefore assume the following dynamics for the one-particle probability density functions α (r, θ, t) and β (r, θ, t), for the species A and B, respectively: with the terms describing self-diffusion, I diff [f ], and alignment interaction, and the terms for the game interaction The collision kernel R θ1,θ2 = 4 d 0 v 0 sin 1 2 (θ 1 − θ 2 ) characterizes the collisions between the agents and denotes the 2π periodic delta distribution. We assume that both random variables, η 0 and η, describing diffusion and collision noise, respectively, are Gaussian distributed: For a more detailed derivation of I diff and I coll , we refer to Ref. [42]. By using the rescaling we can set ω, v 0 , and d 0 equal to 1 without loss of generality. In the following, we will omit the tilde and just write t, r, τ , and f to simplify notation.

Fourier transformed Boltzmann terms
As a first step, we perform the Fourier transformation of the extended kinetic Boltzmann equations (A1) in θ with the Fourier representation where f stands for α or β. With this convention for the Fourier transformation the following identities hold: π −π e ikθ = 2πδ k,0 and ∞ k=−∞ e −ikθ = 2πδ(θ). In the following, we omit spatial and temporal dependencies, and, if not specified differently, all integrals over θ have boundaries −π to π, all integrals over the noise have boundaries ±∞, and sums over wave numbers k run from −∞ to ∞. In order to calculate the Fourier transform of the convection term, we use e θ = (cos θ, sin θ) T . In Fourier space, the convection term reads dθ e ikθ e θ · ∂ r α(θ) = dθ e ikθ cos θ ∂ x + sin θ ∂ y α(θ) where ∇ := ∂ x +i∂ y and ∇ * := ∂ x −i∂ y . The Fourier transform of the diffusion term can be calculated as withP 0,k = exp(−k 2 σ 2 0 /2) the k-th Fourier mode of the Gaussian distribution P 0 (η). The diffusion term consists of a gain and a loss term, given by the first and the second term in Eq. (A2)b, respectively. We denote the Fourier transforms of the gain and the loss terms as (I) and (II), respectively, dθ e ikθ I coll [α, β] = (I) + (II): Analogously, (II) can be calculated as In total, we get for the collision term: with The last remaining part to complete the Fourier representation of Eq. (A1) is I game : and I B game analogously. The Fourier transformation of the equation for species B, Eq. (A1b) works analogous to the transformation for species A shown here. In summary, the Fourier representations of Eqs. (A1) read where we used the fact that the 0-th Fourier mode is the local density of the respective species, since and β 0 (r, t) = b(r, t) analogously.

Uniform disordered solution
One solution of (A17) can be identified by using the fact that (P 0,k − 1) = 0 and I n,k = 0 for k = 0 for all n. Assuming spatial homogeneity, we set ∇α k = ∇β k = 0. The kinetic Boltzmann equations for k = 0 then read (A20) Substituting the control parameterρ = a + b, we get the three solutions and α k = β k = 0 for all k > 0. Eqs. (A20) are the replicator equations described in Section II. As argued there, considering the dynamics of the equations , only the third solution is linearly stable against perturbations. In order to investigate the linear stability of this solution, we linearize the equations (A17) around the solution Eq. (A21c) and calculate the linear growth rate of small perturbations δα k and δβ k . The linearized dynamics of δα k and δβ k is given by (for k > 0) with µ A k and µ B k defined as in Eqs. (17). The solution Eq. (A21c) is therefore unstable when the linear growth rates µ A k or µ B k have positive real parts. BecauseP 0,k is the Fourier transformed distribution of the noise, we haveP 0,k ≤ 1, thusP 0,k − 1 ≤ 0. The linear growth rates µ A k and µ B k can therefore only be positive when I k,k + I 0,k > 0. Calculating (I k,k + I 0,k ) for different values of k > 0 shows that it is positive only for k = 1 (for collision noise values σ < 1) and negative for all k > 1. Hence, small perturbations of the modes α k , β k for k > 1 will decay and only perturbations in α 1 and β 1 may grow (if µ A k > 0 or µ B k > 0, respectively).
The definition of the Fourier series yields for the first mode, k = 1: dθ e iθ β(r, θ, t). (A25) The first modes can be related to the polar order fields P A and P B , which measure the polar order of the system, through Hence, the real and imaginary part of α 1 and β 1 are the x-and y-components of the overall polar order of the system at a given point in space and time. Since we have set the particle velocity v 0 to one, the polar order fields are identical to the velocity fields of both species.
In the following, we will refer to α 1 and β 1 as velocity fields of the species A and B, respectively. This means that when µ A k or µ B k have positive real parts, the absolute value of the polar order of species A or B, respectively, will grow. Indeed, as discussed in section III B, we find critical values for σ (as functions of λ) where either the real part of µ A k or µ B k is zero, indicating a transition from zero to nonzero polar order.

Scaling ansatz
Following [41,42,49] we employ a scaling ansatz to derive hydrodynamic equations for the density, a and b, and velocity fields α 1 and β 1 , of the species A and B, respectively: Specifically, we introduce the small parameter ε and make the ansatz Furthermore, we assume the scaling We can now expand the Fourier representation of the kinetic Boltzmann equations (A17) up to the lowest order which gives non-trivial solutions, which is ε 3 in our case and truncate higher orders. The remaining terms are: k = 1 : ∂ t α 1 + 1 2 (∇α 0 + ∇ * α 2 ) = (P 0,1 − 1)α 1 + (I 0,1 + I 1,1 )a α 1 + (I 2,1 + I −1,1 )α * 1 α 2 + τ (b − λa)α 1 b (A29b) k = 2 : ∂ t α 2 + 1 2 ∇α 1 = (P 0,2 − 1)α 2 + (I 0,2 + I 2,2 )a α 2 + I 1,2 β 2 1 (A29c) We only show the equations for a, α 1 and α 2 here. The equations for b and β 1 are analogue. Because we only take into account O( 3 ) in our final equations for a, b, α 1 , and β 1 , we can drop the terms ∂ t α 2 in (A29c) and solve (A29c) for α 2 . By defining ν(a) = − 1 4[P 0,2 − 1 + (I 0,2 + I 2,2 )a] (A30) we can write α 2 = 4 ν(a)I 1,2 β 2 1 − 2 ν(a)∇α 1 . Substituting this into (A29b), we get Neglecting order O(ε 4 ) and rearranging the equation, we arrive at ∂ t α 1 = (P 0,1 − 1) + (I 0,1 + I 1,1 )a + τ (b − λa)b + 4(I 2,1 + I −1,1 )|α 1 | 2 ν(a)I 1,2 α 1 − 1 2 ∇a + ν(a)∇ * ∇α 1 − 4ν(a)I 1,2 α 1 ∇ * α 1 − 2ν(a)(I 2,1 + I −1,1 )α * 1 ∇α 1 . (A32) Following the same steps for β 1 and with the definitions where x stands for either a or b and we get the final set of hydrodynamic equations for the densities a and b and the polar order fields α 1 and β 1 The set of coupled partial differential equations (A35) describes the evolution of the two densities a and b and the corresponding polar order fields α 1 and β 1 for the species A and B, respectively. When we set the game parameter τ = 0, we recover (for each species) the hydrodynamic equations derived in [41,42]. The fact that the snowdrift interaction term only appears in the term linear in the velocity fields in Eqs. (A35) is a result of the truncation scheme: The snowdrift terms in Eqs. (A17) have the form τ (b − λa)bα k in the dynamic equation for α k and −τ (b − λa)aβ k in the equation for β k . Since we assume a scaling a = O(ε), b = O(ε), α k = O(ε |k| ), β k = O(ε |k| ) for all k = 0, the snowdrift term is of order 3 for k = 0 and k = 1. For all k > 1, the snowdrift term has a order of at least 4 and is therefore cut off in the truncation scheme. Hence, the only snowdrift terms possible are the ones linear in the den-sity and velocity fields, hence ±τ (b−λa)b a, τ (b−λa)b α 1 and τ (b − λa)a β 1 , which all appear in the hydrodynamic Eqs. (A35). by: with q 2 = q 2 x + q 2 y . The eigenvalues of this matrix encode the stability of the fixed point: If for some vector q, the Jacobian J q has an eigenvalue with positive real part, wave-like perturbations with the wave-vector q will grow exponentially, indicating the formation of patterns. On the other hand, if all eigenvalues have negative real part, for all q, spatial perturbations decay and the spatially uniform state is (linearly) stable. In Sec. III C we compute δa/(δa + δb), as a function of τ (see Fig. 4c) in order to see how the snowdrift game affects which density perturbations grow in time. There are four more components to the eigenvector, δα 1 , δα * 1 , δβ 1 and δβ * 1 , encoding the perturbations of the velocity fields. An analogue analysis as in Sec.III C for the velocity field components parallel to macroscopic order, δα 1 and δβ 1 , yields the same phenomenology as for the densities and is shown in Fig. 12: For high values of τ , perturbations in the velocity fields of both species grow equally fast (δα 1 /(δα 1 + δβ 1 ) → 0.5). For τ → 0 on the other hand, perturbations δα 1 only grow in the unstable region in the vicinity of σ A t , while perturbations δβ 1 only grow in the unstable region at σ B t . The components per-pendicular to the macroscopic order, δα ⊥ 1 and δβ ⊥ 1 , are always zero. Hence, perturbations of the velocity field components perpendicular to macroscopic order decay in time.

Appendix C: Numerical methods -hydrodynamic simulation
To investigate the full non-linear dynamics of the hydrodynamic Eqs. (14) and Eqs. (18), we numerically implemented a finite-difference approximation (4th order Runge-Kutta method). Specifically, we model our system by a 2D-grid with periodic boundary conditions and replace all spatial derivatives by their finite difference approximations.
with h a small number. These approximations are the same for partial derivatives. In our case f represents the density fields a and b or the velocity fields α 1 and β 1 (then f is a complex function) and z is a space variable. To solve our equations by finite-differences, we start by defining a square grid in space with grid points The fields are then represented by time dependent N ×N matrices with complex entries in the case of α 1 and β 1 . Using the approximations for spatial derivatives we get the follow-ing replacement rules for the field equations and the corresponding rules for b, α 1 and β 1 . Overall we get 4N 2 ODEs (one ODE for every entry of the four N × N matrices) of the general form ∂ t a i,j = G(a k,l , b k,l , α 1;k,l , β 1;k,l ) with k, l ∈ 1, 2, 3, ..., N . We assume now that we know the state of the system at a time t = 0. To iteratively solve for later times, the 4th order Runge-Kutta (RK4) method worked very well for our system. The general concept of RK4 is the following: Consider the initial value problemḟ where the initial value f 0 and the function F are given. As a first step we discretise the time t n = ∆t n, n ∈ 1, 2, 3, ..., M with ∆t a small number. Hence, time dependent functions are disctretised as f n := f (t n ). The RK4 method now states that f n+1 = f (t n+1 ) can be approximated by In our case, f is an array of the four relevant fields and F = F (f ) is given by the hydrodynamic equations (14,18). We chose initial conditions (a 0,i,j , b 0,i,j , α 1;0,i,j , β 1;0,i,j ) for the fields at every grid point (i, j) and find the solution for t = ∆t, 2∆t, 3∆t, ... via the RK4 method outlined above. We chose periodic boundary conditions for our system. In order to discretise the fields we had to select the size of the 2D grid (the matrices) and the spacing between grid points, h. For our simulations we chose gridsize = 160 × 160, h = 0.4, ∆t = 0.25 * h 2 .
A finer discretization in time and space did not change the outcome of the calculated stationary states of the system. For the simulation results discussed in section III D, the simulation started at homogeneous fields (a 0,i,j , b 0,i,j , α 1;0,i,j , β 1;0,i,j ) = (ρ/2, ρ/2, 0 ,0) for fixed control parameters ρ = 1, σ, τ , λ with small spatial random fluctuations.
We tested the validity of our Runge-Kutta implementation by using COMSOL Multiphysics finite element analysis software and found very well agreement with our implementation.
Appendix D: Agent-based simulations

Implementation and parameters
The agent-based simulation is designed, as described in Sec. IV, as a system of particles in a 2D-box of size L × L with periodic boundary conditions. Each of the N particles in the system has a position r, velocity v 0 = v 0 e(θ), species type S = {A, B}, and fitness f . The simulation is set up in discrete time steps and interactions are evaluated using a modified Verlet algorithm [61]. The update for each time step dt is performed as follows: First, the game theoretical interspecies interaction is evaluated and the particle species are updated according to Eq. (30). In the next step, the intraspecies collisions are processed, the particle directions adjusted, and noise η is added according to Eq. (27). Then, the particle positions are updated according to Eq. (26). As a last step, a snapshot of the system is saved into a hdf5 database for later analysis. These steps are repeated sim steps = T /dt times, where T is the total time of the simulation. Where possible, these steps were implemented utilizing parallel programming.
Parameters used for all shown data are given by dt = 1, L = 256, N = 98304, d align = 1, v 0 = 0.5 (except for Fig. 15: v 0 = 0.015), T = 10 5 and Table I. The hdf5 database generated by the simulation was analyzed using Python scripts; the movies were generated with ffmpeg.

Detailed explanation of the velocity variation parameter
In order to identify polar wave patterns in the system, we introduce the velocity variation parameter ν (Eq. (31)). We therefore choose to divide the system into C channels along the average direction of particle propagation and perpendicular to it (see Fig. 13).  Fig. 11 0.26 0.32 5 Fig. 14a 0.23 0.71 10 Fig. 14b, movie2 ( Each particle in the system is assigned to the corresponding channel based on its location. The channel densities a t c and b t c are calculated at time t for each channel c as well as the corresponding channel polar order parameter α t 1,c and β t 1,c . The velocity variation nu t A is then defined as twice the standard deviation of polar order parameter α t 1,c over all channels: The factor of 2 is introduced such that ν t A ∈ [0, 1]. The velocity variation ν t B for species B is defined accordingly. The temporal averages are denoted by ν A and ν B .
If waves are present in the system, the difference between ν A,⊥ and ν A, for channels perpendicular and parallel to the average direction of motion is nonzero, while for systems without spatial structure they cancel each other. We denote this by ∆ν A = ν A,⊥ − ν A, .
(D2) The channel subdivisions and corresponding densities a t c and b t c and polar order parameters α t 1,c and β t 1,c are also used to generate the wave profile plots such as Fig. 11b.

Five phases
For small enough values of relative fitness λ, we observe the five distinct phases (disorder in Fig. 14a, TWA (see Fig. 14b), partial order (only one species polar ordered, Fig. 14c), TWB (Fig. 14d), and full polar order (Fig. 14e) in the agent-based simulations. The travelling wave phases induced by species A and B (TWA and TWB), as well as the phase of full polar order are shown in the attached movies (see Section E).

Size and consequences of the demixing patterns
A distinct feature of the agent-based simulation is the formation of demixing patterns in environments with slow relative particle movement, compared to the game strength τ . These patterns are a consequence of the range-based interactions in the implementation of the snowdrift game and can also occur in the absence of an alignment interaction, as long as the movement of particles is slow (v 0 /τ 1, s.t. the number of games played by the same set of agents in a row is much larger than one.) This can be seen from the demixing patterns in Fig. 15a, which shows a snapshot of a system where instead of self-propelled particles we simulated particles that show slow diffusive motion. The size of these patterns is proportional to the game range d game (see Fig. 15b). On the other hand, the spatially averaged densities of A and B do not depend on the game range d game . For simulations with alignment interaction, the system can also display demixing patterns in the form of bands (see Fig. 14d).
FIG. 14. Snapshots (top row), density (middle row) and polar order (bottom row) of a set of different noise strength σ as indicated in the graph; the relative fitness is fixed to λ = 0.23 and the interaction range to dgame = 10. (Remaining parameters can be found in Table I.) The profiles for the densities and the polar order are displayed along the direction of polar order as described in Appendix D 2. In the density plot (middle row), the density of species B is adjusted by a factor of inverse relative fitness λ −1 to emphasize the relation between the snowdrift game and local densities of both species; note that at the mean-field level one has a = b/λ. Appendix E: Movie Descriptions movie1.mp4: Simulations of the hydrodynamic equations show the formation of spatial patterns. These spatial patterns are characterised by density-segregated polar ordered bands in both species that move concurrently through the system. When viewed from a comoving frame, the patterns are uniform in the direction perpendicular to the wave vector. The movie shows the profiles of the two density fields a and b over time. Control parameters are set to σ = 0.68, λ = 0.4 and τ = 1.0, such that the system is in its traveling wave phase A (TWA).
Movies 2-4 show simulations of the agent-based approach. Parameter values for movie2-movie4 are given in Table I.
movie2.mp4: The formation of a traveling wave phase induced by species A (TWA). The wave is clearly visible as a broad density pattern and travels in positive y-direction. The dominance of species A can be inferred from the relative abundance of species A compared to B. Only species A shows polar order inside the wave.
movie3.mp4: The formation of a traveling wave phase induced by species B (TWB). The sharply contrasted wave travels in positive x-direction with polar order in both species. Species B is in a parameter regime where it would display wave pattern on its own. Species A on the other hand would be uniformly ordered in this parameter regime and the observed wave patterns are induced by the snowdrift interaction between the species. Furthermore there is clearly visible demixing, as the species form alternating bands inside the wave due to the finite range of the snowdrift interaction.
movie4.mp4: During the formation of order the two species initially cluster together and start to display uniform polar order in both species A and B. The species then spread out over the whole system while maintaining high polar order. In the end small clusters with internal demixing can be observed moving thorugh the system parallel to each other with high polar order.