††thanks: J.M. and M.O. contributed equally to this work.††thanks: J.M. and
M.O. contributed equally to this work.

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

Johanna Mayer1 Michael Obermüller1 Jonas Denk1,2,3 Erwin Frey1<EMAIL_ADDRESS>1Arnold Sommerfeld Center for Theoretical Physics (ASC) and Center for
NanoScience (CeNS), Department of Physics, Ludwig-Maximilians-Universität
München, Theresienstrasse 37, D-80333 München, Germany  2Department of
Physics, University of California, Berkeley, CA 94720, USA  3Department of
Integrative Biology, University of California, Berkeley, CA 94720, USA

###### Abstract

Evolutionary games between species are known to lead to intriguing spatio-
temporal patterns in systems of diffusing agents. 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.

###  I.1 Two-player public good games

Originally, the ‘snowdrift’ game (also known as ‘hawk-dove’ game) was
formulated as a strategic ‘cooperator-defector’ 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 two-player (‘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 kinetics) 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.

###  I.2 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 [10, 28, 11, 29, 30, 31, 12, 8, 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 [35, 29].

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

###  II.1 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 $\theta$. Ballistic motion is interrupted by
stochastic self-diffusion events: Particles change their orientation by an
angle $\eta$ at a rate $\omega$ to

$\theta^{\prime}=\theta+\eta,$ (1)

where $\eta$ is a stochastic variable (noise); see Fig. 1a. Hence, the
particles move in a run-and-tumble-like fashion with a mean time $\omega^{-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.2). 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

$\theta_{1}^{\prime}=\frac{\theta_{1}+\theta_{2}}{2}+\eta_{1},\qquad\theta_{2}^{\prime}=\frac{\theta_{1}+\theta_{2}}{2}+\eta_{2},$
(2)

where we allowed for additional noise terms $\eta_{1}$ and $\eta_{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 $\sigma$.

Figure 1: Illustration of self-propelled dynamics: (a) Particles change their
orientation $\theta$ by a stochastic angle $\eta$ to
$\theta^{\prime}=\theta+\eta$ 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 $\eta_{1}$ and $\eta_{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 one-particle probability density function,
$f\left(\mathbf{r},\theta,t\right)$, which indicates the probability that a
particle is located at position $\mathbf{r}$ with the orientation $\theta$ at
time $t$. The kinetic Boltzmann equation for self-propelled particles for one
species reads [41, 42]

$\partial_{t}f\left(\mathbf{r},\theta,t\right)+v_{0}\,\mathbf{e}_{\theta}\cdot\nabla
f=I_{\text{diff}}\left[f\right]+I_{\text{coll}}\left[f,f\right],$ (3)

where $I_{\text{diff}}$ and $I_{\text{coll}}$ denote contributions resulting
from diffusion and collision events, respectively, while $\mathbf{e}_{\theta}$
is the unit vector pointing in $\theta$-direction. Their explicit expressions
can be found in Appendix A.1. In the present context, the most important
feature of $I_{\text{diff}}$ and $I_{\text{coll}}$ is that they depend on the
parameter $\sigma$, which characterizes the noise strength.

###  II.2 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 [46, 6]:

$\displaystyle f_{A}$ $\displaystyle=1+\tau\,b,$ (4a) $\displaystyle f_{B}$
$\displaystyle=1+\tau\,\lambda\,a\,.$ (4b)

Here, $a$ and $b$ denote the densities of species $A$ and $B$, respectively.
In the following, we will refer to $\lambda$ as the relative fitness parameter
of particles of species $B$ relative to $A$, and $\tau$ as the game strength.

Figure 2:  Illustration of the snowdrift game. (a) Payoff matrix: There is a
non-zero payoff only if individuals belonging to different species meet. Then,
$A$ receives a payoff $\tau$, while $B$ receives a payoff $\lambda\tau$. The
relative fitness $\lambda$ and game strength $\tau$ serve as control
parameters. (b) Flow diagram of the replicator equation, Eq. (7), for species
$A$ with population density $a$: Unstable fixed points are marked with an open
circle, while the stable fixed point at $a_{0}$ is marked with a filled
circle. The arrows indicate the flow towards the stable fixed point.

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:

$\displaystyle\partial_{t}a$ $\displaystyle=\left(f_{A}-f_{B}\right)b\,a,$
(5a) $\displaystyle\partial_{t}b$
$\displaystyle=\left(f_{B}-f_{A}\right)a\,b.$ (5b)

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 $\rho=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

$\displaystyle\partial_{t}a$ $\displaystyle=\tau\left(b-\lambda a\right)a\,b,$
(6a) $\displaystyle\partial_{t}b$ $\displaystyle=\tau\left(\lambda
a-b\right)a\,b.$ (6b)

The meaning of the terms relative fitness and game strength for $\lambda$ and
$\tau$, respectively, now become evident: If $b\,{>}\,\lambda a$, there is a
fitness advantage for species $A$ and its density grows relative to species
$B$. In contrast, if $b\,{<}\,\lambda 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
$\lambda\,{\in}\,[0,1]$. By changing $\lambda$ we can tune the relative
fitness from the two species being equally competitive ($\lambda\,{=}\,1$) to
species $A$ being strongly dominant over $B$ ($\lambda\,{\ll}\,1$). The second
game parameter, $\tau\,{\in}\,[0,\infty)$, can be understood as the game
strength.

Since the total particle number is conserved ($b\,{=}\rho\,{-}\,a$), Eqs. (6)
can be reduced to the temporal evolution of a single species (say $A$):

$\partial_{t}a=\tau\left(b-\lambda
a\right)ab=\tau\left[\left(\rho-a\right)-\lambda
a\right]a\left(\rho-a\right).$ (7)

Figure 2b shows the flow diagram of Eq. (7) with arrows pointing to the right
for increasing $a$ ($\partial_{t}a\,{>}\,0$) and to the left for
$\partial_{t}a\,{<}\,0$. There are three fixed points (stationary points with
$\partial_{t}a=0$):

$a_{B}=0,~{}a_{A}=\rho,\textrm{ and }a_{0}=\frac{\rho}{1+\lambda}.$ (8)

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\,{=}\,\lambda$.

###  II.3 Active matter snowdrift game

To describe our system of two interacting self-propelled particle species, we
introduce the one-particle probability density functions
$\alpha\left(\mathbf{r},\theta,t\right)$ and
$\beta\left(\mathbf{r},\theta,t\right)$, 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_{\text{game}}^{A}\left[\alpha,\beta\right]$ and
$I_{\text{game}}^{B}\left[\beta,\alpha\right]$, chosen in analogy to the
replicator equations, Eqs. (6). This leads to the general form of the two-
species Boltzmann equations

$\displaystyle\partial_{t}\alpha{+}v_{0}\,\mathbf{e}_{\theta}\nabla\alpha$
$\displaystyle=I_{\text{diff}}\left[\alpha\right]{+}I_{\text{coll}}\left[\alpha,\alpha\right]{+}I_{\text{game}}^{A}\left[\alpha,\beta\right],$
(9a) $\displaystyle\partial_{t}\beta{+}v_{0}\,\mathbf{e}_{\theta}\nabla\beta$
$\displaystyle=I_{\text{diff}}\left[\beta\right]{+}I_{\text{coll}}\left[\beta,\beta\right]{+}I_{\text{game}}^{B}\left[\beta,\alpha\right].$
(9b)

We assume snowdrift interactions to be local and independent of particle
orientation as in Eqs. (6). With the local densities
$a\left(\mathbf{r},t\right)$ and $b\left(\mathbf{r},t\right)$ of species $A$
and $B$, related to $\alpha$ and $\beta$ by

$\displaystyle a\left(\mathbf{r},t\right)$
$\displaystyle=\int_{-\pi}^{\pi}d\theta\,\alpha\left(\mathbf{r},\theta,t\right),$
(10a) $\displaystyle b\left(\mathbf{r},t\right)$
$\displaystyle=\int_{-\pi}^{\pi}d\theta\,\beta\left(\mathbf{r},\theta,t\right),$
(10b)

the snowdrift interaction terms read

$\displaystyle I_{\text{game}}^{A}\left[\alpha,\beta\right]$
$\displaystyle=\tau\left(b-\lambda a\right)\,b\,\alpha,$ (11a) $\displaystyle
I_{\text{game}}^{B}\left[\beta,\alpha\right]$ $\displaystyle=\tau\left(\lambda
a-b\right)\,a\,\beta.$ (11b)

The dependency on $\alpha$ and $\beta$ 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 $\alpha$ and $\beta$. The game
strength $\tau$ tunes the relative weight of the snowdrift interaction
compared to the alignment interaction (collision terms). If one sets
$\tau\,{=}\,0$, the snowdrift interaction terms vanish and Eqs. (9) reduce to
two independent equations for species $A$ and $B$. On the other hand, in the
absence of convection ($v_{0}\,{=}\,0$), diffusion
($I_{\text{diff}}\,{=}\,0$), and collision ($I_{\text{coll}}\,{=}\,0$), one
recovers the replicator equations (6) by integrating Eqs. (9) over the angle
$\theta$.

In the following, time, space, and densities are rescaled such that one
measures time in units of $1/\omega$, position in units of $v_{0}/\omega$ and
spatial densities in units of $\omega/(v_{0}d_{0})$. This amounts to setting
$v_{0}\,{=}\,d_{0}\,{=}\,\omega\,{=}\,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 [47, 48, 45]. However, as exemplified by studies
of the kinetic Boltzmann equation for a single species [41, 49, 42, 48, 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.

###  III.1 Derivation of hydrodynamic equations

As the probability densities $\alpha(\mathbf{r},\theta,t)$ and
$\beta(\mathbf{r},\theta,t)$ are periodic functions in $\theta$, it is
convenient to work with their Fourier modes, defined as
$\alpha_{k}(\mathbf{r},t):=\int_{-\pi}^{\pi}d\theta\,e^{i\theta
k}\alpha(\mathbf{r},\theta,t)$ and
$\beta_{k}(\mathbf{r},t):=\int_{-\pi}^{\pi}d\theta\,e^{i\theta
k}\beta(\mathbf{r},\theta,t)$. Note that the zeroth Fourier modes are the
densities of the respective species, $\alpha_{0}\,{=}\,a,\beta_{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
$\mathbf{P}_{A}(\mathbf{r},t)$ and $\mathbf{P}_{B}(\mathbf{r},t)$, which are
related to the first Fourier modes

$\displaystyle
a\,\mathbf{P}_{A}=\left(\begin{array}[]{c}\text{Re}(\alpha_{1})\\\
\text{Im}(\alpha_{1})\end{array}\right),\quad
b\,\mathbf{P}_{B}=\left(\begin{array}[]{c}\text{Re}(\beta_{1})\\\
\text{Im}(\beta_{1})\end{array}\right).$ (16)

Hence, the modes $\alpha_{1}$ and $\beta_{1}$ are a measure for the polar
order multiplied by the density at a given point in space and time of the
system. If $\alpha_{1}\,{=}\,\beta_{1}\,{=}\,0$, the system is disordered. If
$\alpha_{1}$ or $\beta_{1}$ is nonzero, there is a preferred direction of
motion, meaning the respective species moves collectively. We will refer to
$\alpha_{1}$ and $\beta_{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)

$\displaystyle\partial_{t}\alpha_{k}+\frac{1}{2}$
$\displaystyle\bigl{(}\nabla\alpha_{k-1}+\nabla^{*}\alpha_{k+1}\bigr{)}=-\bigl{(}1-e^{-(k\sigma)^{2}/2}\bigr{)}\,\alpha_{k}$
$\displaystyle+\sum\limits_{p}\alpha_{k-p}\alpha_{p}\,\mathcal{I}_{p,k}+\tau(b-\lambda
a)\,b\,\alpha_{k},$ (17a) $\displaystyle\partial_{t}\beta_{k}+\frac{1}{2}$
$\displaystyle\bigl{(}\nabla\beta_{k-1}+\nabla^{*}\beta_{k+1}\bigr{)}=-\bigl{(}1-e^{-(k\sigma)^{2}/2}\bigr{)}\,\beta_{k}$
$\displaystyle+\sum\limits_{p}\beta_{k-p}\beta_{p}\,\mathcal{I}_{p,k}-\tau(b-\lambda
a)\,a\,\beta_{k},$ (17b)

where we introduced the abbreviations
$\nabla\,{=}\,\partial_{x}+i\partial_{y}$ and
$\nabla^{*}\,{=}\,\partial_{x}-i\partial_{y}$. The explicit expressions for
the collision kernels $\mathcal{I}_{p,k}(\sigma)$ are given in Appendix A.
Equation (17) constitutes an infinite set of differential equations that
couple lower- with higher-order Fourier modes. For $k\,{=}\,0$, Eqs. (17)
yield

$\displaystyle\partial_{t}a$
$\displaystyle=-\frac{1}{2}\left(\nabla\alpha_{1}^{*}+\nabla^{*}\alpha_{1}\right)+\tau\left(b-\lambda
a\right)\,b\,a,$ (18a) $\displaystyle\partial_{t}b$
$\displaystyle=-\frac{1}{2}\left(\nabla\beta_{1}^{*}+\nabla^{*}\beta_{1}\right)+\tau\left(\lambda
a-b\right)\,a\,b.$ (18b)

The first terms in Eqs. (18) 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 $\rho\,{=}\,a+b$ follows a continuity equation

$\partial_{t}\rho=-\frac{1}{2}\Bigl{[}\nabla(\alpha_{1}^{*}+\beta_{1}^{*})+\nabla^{*}(\alpha_{1}+\beta_{1})\Bigr{]}.$
(19)

The spatial average of the density, $\bar{\rho}$, 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. (18), vanish and Eqs. (18) reduce to the
replicator equations. The corresponding steady states are given by the stable
snowdrift fixed point with $a^{(0)}\,{=}\,\bar{\rho}/(1+\lambda)$ and
$b^{(0)}\,{=}\,\bar{\rho}\lambda/(1+\lambda)$. Hence, a state with uniform
densities $a\,{=}\,a^{(0)}$ and $b\,{=}\,b^{(0)}$ and all higher Fourier modes
vanishing is a solution to Eqs. (17). This solution describes a spatially
uniform, disordered system with densities set to their snowdrift fixed points.
Small spatially uniform perturbations $\delta\alpha_{k},\delta\beta_{k}$ to
this solution evolve to linear order as

$\displaystyle\partial_{t}\delta\alpha_{k}$
$\displaystyle=\mu_{k}^{A}(a^{(0)},b^{(0)})\,\delta\alpha_{k},$ (20a)
$\displaystyle\partial_{t}\delta\beta_{k}$
$\displaystyle=\mu^{B}_{k}(a^{(0)},b^{(0)})\,\delta\beta_{k},$ (20b)

with

$\displaystyle\mu^{A}_{k}(a,b)=$
$\displaystyle\bigl{(}e^{{-}\frac{\sigma^{2}}{2}}{-}1\bigr{)}+(\mathcal{I}_{0,k}{+}\mathcal{I}_{k,k})a+\tau(b{-}\lambda
a)b,$ (21a) $\displaystyle\mu^{B}_{k}(a,b)=$
$\displaystyle\bigl{(}e^{{-}\frac{\sigma^{2}}{2}}{-}1\bigr{)}+(\mathcal{I}_{0,k}{+}\mathcal{I}_{k,k})b+\tau(\lambda
a{-}b)a.$ (21b)

The respectively first two terms in $\mu^{A}_{k}$ and $\mu^{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)}\,{=}\,\bar{\rho}/(1+\lambda)$ and
$b^{(0)}\,{=}\,\bar{\rho}\lambda/(1+\lambda)$, into the expressions for
$\mu^{A}_{k}$ and $\mu^{B}_{k}$, the snowdrift terms vanish and we are left
with the first two terms, which are then functions of $\sigma$, $\bar{\rho}$
and $\lambda$. Varying these parameters, one finds that all Fourier modes
$\mu_{k}^{A}$ and $\mu_{k}^{B}$ with $k\geq 2$ are negative and only the
growth rates for the first Fourier mode, $\mu^{A}_{1}$ and $\mu^{B}_{1}$, can
become positive, meaning that small perturbations in the velocity fields
$\alpha_{1}$ or $\beta_{1}$, respectively, will grow exponentially. For fixed
overall density, $\bar{\rho}$, and relative game strength, $\lambda$, the
growth rates $\mu^{A}_{1}$ and $\mu^{B}_{1}$ are negative for high values of
the noise strength $\sigma$ and become positive when the noise strength is low
enough. This defines threshold values for the noise strength,
$\sigma^{A}_{t}(\lambda,\bar{\rho})$ and $\sigma^{B}_{t}(\lambda,\bar{\rho})$,
at $\mu_{1}^{A}|_{\sigma^{A}_{t},\bar{\rho},\lambda}\,{=}\,0$ and
$\mu_{1}^{B}|_{\sigma^{B}_{t},\bar{\rho},\lambda}\,{=}\,0$, respectively (see
Fig. 3b). Below these thresholds ($\sigma\,{<}\,\sigma^{A}_{t}$ or
$\sigma\,{<}\,\sigma^{B}_{t}$), the spatially uniform solution is unstable and
small perturbations in the velocity fields, $\alpha_{1}$ or $\beta_{1}$, have
positive growth rates.

Close to the threshold values $\sigma_{t}^{A}$ and $\sigma_{t}^{B}$, 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, $\alpha_{1},\beta_{1}\,{\ll}\,1$, and obey the
scaling relations $(a-a^{(0)})\,{\sim}\,\alpha_{1}$,
$\alpha_{k}\,{\sim}\,\alpha_{1}^{|k|}$ (and analogously for $\beta_{k}$).
Furthermore, we assume that spatial and temporal variations
$\nabla,\,\partial_{t}\,{\sim}\,\alpha_{1}$ are small. With these assumptions,
one can truncate and close the Boltzmann equation in Fourier space, Eq. (17),
at the third order term for the velocity field, $\alpha_{1}^{3}$. 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 equations for the velocity fields of species $A$ and $B$:

$\displaystyle\begin{split}\partial_{t}\alpha_{1}&=\left(\mu^{A}_{1}(a,b)-\xi(a)|\alpha_{1}|^{2}\right)\alpha_{1}+\nu(a)\,\nabla^{*}\nabla\alpha_{1}\\\
&\qquad-\gamma(a)\,\alpha_{1}\nabla^{*}\alpha_{1}-\kappa(a)\,\alpha_{1}^{*}\nabla\alpha_{1}-\frac{1}{2}\nabla
a,\end{split}$ (22a)
$\displaystyle\begin{split}\partial_{t}\beta_{1}&=\left(\mu^{B}_{1}(a,b)-\xi(b)|\beta_{1}|^{2}\right)\beta_{1}+\nu(b)\,\nabla^{*}\nabla\beta_{1}\\\
&\qquad-\gamma(b)\,\beta_{1}\nabla^{*}\beta_{1}-\kappa(b)\,\beta_{1}^{*}\nabla\beta_{1}-\frac{1}{2}\nabla
b.\end{split}$ (22b)

For explicit expressions of the coefficients please refer to Appendix A. The
essential difference between Eqs. (22) and the hydrodynamic equation for the
velocity field of one-species SPP [41] lies in the coefficients $\mu^{A}_{1}$
and $\mu^{B}_{1}$, which here couple the velocity fields to the density of the
respective other species (see Eqs. (21)). Note that the snowdrift interaction
term only appears in the term linear in the velocity fields in Eqs. (22). This
is a consequence of the truncation scheme and is explained in more detail in
Appendix A.

Taken together, Eqs. (18) and Eqs. (22) 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 $\sigma$, the overall
density $\bar{\rho}$, the game strength $\tau$, and the relative fitness
$\lambda$. For specificity, we use a fixed value for the average density,
$\bar{\rho}\,{=}\,1$ (in units of $\omega/(v_{0}d_{0})$) and vary the
remaining parameters. We observe the same phenomenology if we instead fix
$\sigma$ and use $\bar{\rho}$ as control parameter [52] (this is due to the
opposing roles of density and noise changes in the self-organization of active
matter systems [19, 16, 17]).

In the next sections we will first determine solutions for the spatially
uniform case of Eqs. (18) and Eqs. (22) and subsequently analyse their linear
stability against spatially non-uniform perturbations.

###  III.2 Spatially uniform solutions

The stationary uniform solutions of the hydrodynamic equations, Eqs. (18) and
(22), 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. (18), one obtains

$\displaystyle\partial_{t}a$ $\displaystyle=\tau(b-\lambda a)\,a\,b,$ (23a)
$\displaystyle\partial_{t}b$ $\displaystyle=\tau(\lambda a-b)\,a\,b.$ (23b)

These are precisely the replicator equations, Eqs. (6), with the stable
solution (note we have set $\bar{\rho}=1$)

$a^{(0)}=\frac{1}{1+\lambda},\,\quad b^{(0)}=\frac{\lambda}{1+\lambda}.$ (24)

Hence, the spatially uniform densities of each species are determined by the
relative fitness $\lambda$, and are independent of the game strength $\tau$.

The hydrodynamic equations for the velocity fields, Eqs. (22), of a spatially
uniform system are given by

$\displaystyle\partial_{t}\alpha_{1}$
$\displaystyle=\left(\mu^{A}_{1}-\xi|\alpha_{1}|^{2}\right)\alpha_{1},$ (25a)
$\displaystyle\partial_{t}\beta_{1}$
$\displaystyle=\left(\mu^{B}_{1}-\xi|\beta_{1}|^{2}\right)\beta_{1}.$ (25b)

Hence, $\alpha_{1}\,{=}\,\beta_{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 $\xi\,{>}\,0$ for all values
of the control parameters (see Appendix A), such a solution for $\alpha_{1}$
($\beta_{1}$) exists for $\mu^{A}_{1}\,{>}\,0$ ($\mu^{B}_{1}\,{>}\,0$). We
know from the last section that the growth rate $\mu^{A}_{1}$ ($\mu^{B}_{1}$)
is negative at high values of the noise strength $\sigma$ and becomes positive
below the threshold value $\sigma_{t}^{A}$ ($\sigma_{t}^{B}$), which depends
on the relative fitness $\lambda$, 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
$\lambda\,{\neq}\,1$. The ‘weaker’ species (here B) has a lower density and
hence is more prone to the disordering effect of noise, i.e.
$\sigma_{t}^{B}\,{\leq}\,\sigma_{t}^{A}$. Moreover, all results are
independent of the game strength $\tau$, since the terms involving $\tau$ in
$\mu_{1}^{A}$ and $\mu_{1}^{B}$ (see Eqs. (21)) cancel as the densities reach
their spatially uniform steady state solutions $a^{(0)}$ and $b^{(0)}$ given
by Eqs. (24). The game strength $\tau$ 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
$\alpha_{1}\,{=}\,\beta_{1}\,{=}\,0$ is linearly unstable and additional,
stable solutions appear, given by

$\displaystyle\alpha_{1}^{(0)}=\sqrt{\frac{\mu^{A}_{1}}{\xi}}\,e^{i\theta_{a}}\qquad\text{for}\,\,\sigma<\sigma_{t}^{A},$
(26a)
$\displaystyle\beta_{1}^{(0)}=\sqrt{\frac{\mu^{B}_{1}}{\xi}}\,e^{i\theta_{b}}\qquad\text{for}\,\,\sigma<\sigma_{t}^{B}.$
(26b)

These solutions correspond to collective states which spontaneously break the
rotational symmetry of the system by choosing random directions of macroscopic
order $\theta_{a}$ and $\theta_{b}$, respectively.

Figure 3: Solutions of the spatially uniform equations for the velocity fields
$\alpha_{1}$ and $\beta_{1}$ (see Eqs. (25)). (a) Amplitudes of the solutions
as functions of the noise strength $\sigma$ with the relative fitness set to
$\lambda\,{=}\,1$, $\lambda\,{=}\,0.6$, and $\lambda\,{=}\,0.2$, as indicated
in the graph. Solid lines are stable solutions, dashed lines unstable
solutions. If the noise strength is smaller than a threshold value,
$\sigma\,{\leq}\,\sigma_{t}^{A}\,(\sigma\,{\leq}\,\sigma_{t}^{B})$, the
disordered solution $|\alpha_{1}^{(0)}|\,{=}\,0$ $(|\beta_{1}^{(0)}|\,{=}\,0)$
becomes unstable and a new stable solution with $|\alpha_{1}^{(0)}|\,{>}\,0$
$(|\beta_{1}^{(0)}|\,{>}\,0)$ appears. (b) Bifurcation (phase) diagram for the
spatially uniform solution in control parameter space ($\lambda,\sigma$). At
$\sigma_{t}^{A}(\lambda)$, there is a transitions from a disordered state
(fixed point $|\alpha_{1}^{(0)}|\,{=}\,|\beta_{1}^{(0)}|\,{=}\,0$) to a broken
symmetry state with macroscopic polar order in $A$ only
($|\alpha_{1}^{(0)}|\,{>}\,0$, $|\beta_{1}^{(0)}|\,{=}\,0$). $A$ second
transition to a state with polar order in both species
($|\alpha_{1}^{(0)}|\,{>}\,0$ and $|\beta_{1}^{(0)}|\,{>}\,0$) occurs at
$\sigma_{t}^{B}(\lambda)\,{\leq}\,\sigma_{t}^{A}(\lambda)$. The threshold
values for the noise strengths coincide if the relative game strength
$\lambda\,{=}\,1$.

Figure 3a shows the absolute values (amplitudes) $|\alpha_{1}^{(0)}|$ and
$|\beta_{1}^{(0)}|$ as functions of the noise strength $\sigma$ for different
values of $\lambda$. For equally competitive species, $\lambda\,{=}\,1$, the
densities are the same for both species: $a^{(0)}\,{=}\,b^{(0)}\,{=}\,1/2$
(see Eqs. (24)). In this case the growth rates $\mu_{1}^{A}(a^{(0)},b^{(0)})$
and $\mu_{1}^{B}(a^{(0)},b^{(0)})$ are identical for all values of the noise
strength $\sigma$ (see Eqs. (21)). Thus, the amplitudes of the velocity
fields, $|\alpha_{1}^{(0)}|$ and $|\beta_{1}^{(0)}|$, are also identical as
shown in the upper panel of Fig. 3. When $\sigma$ is decreased below
$\sigma_{t}^{A}\,\,{=}\,\sigma_{t}^{B}$, the amplitudes of the velocity fields
become non-zero, which means that both species exhibit spatially uniform
macroscopic order. For $\lambda\,{<}\,1$, the amplitudes $|\alpha_{1}^{(0)}|$
and $|\beta_{1}^{(0)}|$ differ. Since a relative game strength
$\lambda\,{<}\,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,
$\mu_{1}^{A}\,{>}\,\mu_{1}^{B}$, and thus polar order of $A$ is less prone to
noise resulting in a larger threshold value,
$\sigma_{t}^{A}\,{>}\,\sigma_{t}^{B}$, see Fig. 3b and Fig. 3a for
$\lambda\,{=}\,0.6$ and $\lambda\,{=}\,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 $\sigma$ and the relative fitness $\lambda$ there are three
phases: a disordered phase above $\sigma_{t}^{A}$, an intermediate phase for
$\sigma_{t}^{A}\,{>}\,\sigma\,{>}\,\sigma_{t}^{B}$ where only the more
competitive species, $A$, shows polar order, and a phase where both species
exhibit polar order below $\sigma_{t}^{B}$.

###  III.3 Stability against spatially non-uniform perturbations

In this section we will test the linear stability of the spatially uniform
solutions given by Eqs. (24) and Eqs. (26), against perturbations with a
finite wavelength. To this end, we introduce spatially non-uniform
perturbations around the uniform solutions

$\begin{split}a(\mathbf{r},t)&=a^{(0)}\,{+}\,\delta
a(\mathbf{r},t),\;\alpha_{1}(\mathbf{r},t)=\alpha_{1}^{(0)}\,{+}\,\delta\alpha_{1}(\mathbf{r},t),\\\
b(\mathbf{r},t)&=b^{(0)}\,{+}\,\delta
b(\mathbf{r},t),\;\,\,\beta_{1}(\mathbf{r},t)=\beta_{1}^{(0)}\,{+}\,\delta\beta_{1}(\mathbf{r},t).\end{split}$
(27)

Assuming that the perturbations are small, we insert the perturbed fields in
the full hydrodynamic equations, Eq. (18) and Eq. (22), and keep only terms up
to linear order in the perturbation. Note that one has to account for both
perturbations $\delta\alpha_{1}$, $\delta\beta_{1}$ and its complex conjugates
$\delta\alpha_{1}^{*}$, $\delta\beta_{1}^{*}$. Introducing Fourier modes
$\delta
a(\mathbf{r},t)\,{=}\,\frac{1}{2\pi}\int_{-\infty}^{\infty}\mathrm{d}\mathbf{q}\,\mathrm{e}^{\mathrm{i}\mathbf{q}\cdot\mathbf{r}}\delta
a(\mathbf{q},t)$ and analogously for all other fields, one obtains a set of
linear equations for $\delta\mathbf{m}_{\mathbf{q}}\,{=}\,(\delta
a(\mathbf{q}),\delta
b(\mathbf{q}),\delta\alpha_{1}(\mathbf{q}),\delta\alpha_{1}^{*}(\mathbf{q}),\delta\beta_{1}(\mathbf{q}),\delta\beta_{1}^{*}(\mathbf{q}))^{T}$

$\partial_{t}\,\delta\mathbf{m}_{\mathbf{q}}=\mathbf{J}_{\mathbf{q}}\,\delta\mathbf{m}_{\mathbf{q}}$
(28)

with the Jacobian $\mathbf{J}_{\mathbf{q}}$. It is a $6{\times}6$ matrix that
depends on the wave vector $\mathbf{q}$ and the spatially uniform solution
$(a^{(0)},b^{(0)},\alpha_{1}^{(0)},\beta_{1}^{(0)})$ and therefore on all
control parameters $(\lambda,\tau,\sigma)$. An explicit expression for
$\mathbf{J}_{\mathbf{q}}$ can be found in Appendix B. For simplification, we
assume that the spatially uniform solutions of the velocity fields are
parallel, $\alpha_{1}^{(0)}\,{\parallel}\,\beta_{1}^{(0)}$. Numerical
solutions of the hydrodynamic equations and agent-based simulations in the
next sections will show that this assumption is justified.

Figure 4:  Stability analysis of the hydrodynamic equations Eqs. (18) and
(22). (a) Dispersion relation of the eigenvalue with the largest real part,
$s_{\mathbf{q}}$, for different values of the noise strength $\sigma$, and
fixed game strength $\tau\,{=}\,0.2$ and relative fitness $\lambda\,{=}\,0.4$.
$s_{\mathbf{q}}$ is plotted for wave vector $\mathbf{q}$ pointing in the
direction parallel to macroscopic order,
$\mathbf{q}\parallel\alpha_{1},\beta_{1}$ (denoted by
$\mathbf{q}_{\parallel}$). The spatially uniform state, Eqs. (26), is unstable
when the maximum of the dispersion relation, $s_{\text{max}}$, is larger zero.
(b) Parameter regime (bands) of $\lambda$ (relative fitness) and $\sigma$
(noise strength) where the spatially uniform polar-ordered states are linearly
unstable against spatial perturbations (i.e. $s_{\text{max}}>0$), for
different values of the game strength $\tau$ (shades of blue) as indicated in
the graph. The two bands are bounded from above by the threshold values
$\sigma_{t}^{A}$ and $\sigma_{t}^{B}$ 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 $\tau$. White
regions indicate parameter regimes where the spatially uniform ordered states
are stable against spatial perturbations. (c), (d) Linearly unstable regions
in ($\tau,\sigma$)-space and ($\lambda,\sigma$)-space, respectively. The color
code compares the two density field components ($\delta a$ and $\delta b$) of
the eigenvector corresponding to the eigenvalue with the largest positive real
part: $\delta_{a}\,{:=}\,\delta a/(\delta a+\delta b)$ which ranges from $0$
(the largest linear growth rate corresponds to an eigenvector with $\delta
a$-component equal to zero and a non-zero $\delta b$-component) to $1$
($\delta b$-component equal to zero and non-zero $\delta a$-component).

The real parts of the eigenvalues of the Jacobian $\mathbf{J}_{\mathbf{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_{\mathbf{q}}$,
as a function of $\mathbf{q}$ for different values of $\sigma$ and fixed
$\lambda\,{=}\,0.4$ and $\tau\,{=}\,0.2$. For certain values of the noise
strength $\sigma$ there are long wavelength instabilities (i.e. for small
values of the wave vector $\mathbf{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, $\sigma_{t}^{A}$ and $\sigma_{t}^{B}$ (Fig.
4b). In other words, the disordered phase ($\sigma>\sigma_{t}^{A}$) is stable,
while the phases with spatially uniform polar order in species $A$
($\sigma_{t}^{A}>\sigma>\sigma_{t}^{B}$) and uniform polar order in both
species ($\sigma_{t}^{B}>\sigma$) are both unstable towards spatial
perturbations for $\sigma$ close enough to the respective onset of order. By
computing the real part of the largest eigenvalue, $s_{\mathbf{q}}$, for
different wave vectors $\mathbf{q}$, we find that $s_{\mathbf{q}}$ is largest
for wave vectors parallel to the direction of macroscopic order
($\mathbf{q}\parallel\alpha_{1}^{(0)},\beta_{1}^{(0)}$); i.e. instabilities
are strongest when longitudinal. This suggests that the ensuing patterns will
form along the direction parallel to macroscopic order, $\alpha_{1}^{(0)}$ and
$\beta_{1}^{(0)}$, while the system remains homogeneous along the direction
perpendicular to macroscopic order. For $\tau\,{=}\,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, $\tau\,{>}\,0$, the upper boundaries of the
unstable regions in parameter space are still $\sigma_{t}^{A}$ and
$\sigma_{t}^{B}$ (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 $\delta\mathbf{m}_{\mathbf{q}}$ corresponding to the eigenvalue
$s_{\mathbf{q}}$ with the largest positive real part (and with wave vector
$\mathbf{q}$ fixed such that $s_{\mathbf{q}}$ is maximal); see Eq. (28). We
compute this eigenvector as a function of the three control parameters
$\sigma,\lambda$ and $\tau$ 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,
$\delta a$ and $\delta b$. We use the expression $\delta_{a}:=\delta a/(\delta
a+\delta b)$ to measure the direction of the eigenvector projected onto the
plane spanned by the two density fields. This value can range from $0$
($\delta a\,{=}\,0$ and non-zero $\delta b$-component), meaning that only
perturbations in the density of species $B$ grow, to $1$ ($\delta b\,{=}\,0$
and non-zero $\delta a$-component) meaning that only the density of species
$A$ is unstable towards perturbations. Fig. 4c and Fig. 4d show the value of
$\delta_{a}$ as a function of the control parameters in ($\tau$,$\sigma$)- and
($\lambda$,$\sigma$)-space, respectively. Analysing the results in
($\tau$,$\sigma$)-space (Fig. 4c), we can read off how the game strength,
$\tau$, affects the relative growth rates of perturbations of the two
densities, $a$ and $b$: For large enough game strength, $\tau$, perturbations
in the densities of both species grow equally fast ($\delta_{a}\rightarrow
0.5$). Decreasing $\tau$ changes the value of $\delta_{a}$ dependent on the
region in parameter space: In the unstable region with upper boundary
$\sigma\,{=}\,\sigma_{t}^{A}$, perturbations in the density of species $A$ are
growing faster than perturbations in the density of species $B$ such that
$\delta_{a}\rightarrow 1$ for small enough game strength $\tau$ (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
$\sigma\,{=}\,\sigma_{t}^{B}$, $\delta_{a}\,{\rightarrow}\,0$ for small enough
$\tau$, which means that only the density of species $B$ is unstable. The
inset of Fig. 4d shows the limiting case $\tau\,{=}\,0$ in the
$(\lambda,\sigma)$-plane. From the color code we read off that, depending on
the region, either $\delta b\,{=}\,0$ (yellow region) or $\delta a\,{=}\,0$
(blue region). In other words, for two non-interacting species
($\tau\,{=}\,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 $\delta a$
grow at $\sigma\,{<}\,\sigma_{t}^{A}$; perturbations $\delta b$ grow at
$\sigma\,{<}\,\sigma_{t}^{B}$.

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, $\delta\alpha_{1}^{\parallel}$ and
$\delta\beta_{1}^{\parallel}$, yields the same phenomenology as for the
densities: For high values of $\tau$, perturbations in the velocity fields of
both species grow equally fast
($\delta\alpha_{1}^{\parallel}/(\delta\alpha_{1}^{\parallel}+\delta\beta_{1}^{\parallel})\rightarrow
0.5$). For $\tau\rightarrow 0$, on the other hand, perturbations
$\delta\alpha_{1}^{\parallel}$ only grow in the unstable region in the
vicinity of $\sigma_{t}^{A}$, while perturbations
$\delta\beta_{1}^{\parallel}$ only grow in the unstable region at
$\sigma_{t}^{B}$. 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, $\sigma_{t}^{A}$ and $\sigma_{t}^{B}$, 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, $\lambda$. The game strength, $\tau$, affects the relative
growth rates of perturbations of the two densities, $a$ and $b$, and
velocities parallel to macroscopic order, $\alpha_{1}^{\parallel}$ and
$\beta_{1}^{\parallel}$: For high values of $\tau$, perturbations in the
fields of both species grow equally fast; for decreasing $\tau$, growth rates
differ between the species until for $\tau\rightarrow 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.

###  III.4 Numerical solutions of the hydrodynamic equations

The spatially uniform solutions, Eqs. (24) and (26), and their linear
stability analysis (Fig. 4) suggest five distinct phases. Especially in the
regimes where spatially uniform ordered states are unstable against spatially
non-uniform 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. (18) and (22), 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.

####  III.4.1 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.2 and III.3), 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, polar-
ordered 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.3). 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 [53, 42, 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.4.2. We have also investigated phase diagrams
for different values of $\tau$ and found the same topology.

Figure 5: Phase diagram in control parameter space ($\lambda,\sigma$) for game
strength $\tau\,{=}\,0.2$. Coloured regions indicate parameter regimes where
one finds phases that exhibit traveling polar waves, traveling polar wave
phase A (TWA) and traveling polar wave phase B (TWB). White regions indicate
regimes showing spatially uniform phases. Regions of linear instability are in
yellow-brown (depending on the value of the maximum of the largest real part
of the eigenvalue of the Jacobian, $s_{\text{max}}$). Bistable regions at the
transitions between phases are coloured in blue. Blue dots are calculated
numerically and indicate transitions between spatially uniform phases,
bistable regions and traveling wave phases. The blue star indicates the
vicinity of a small parameter regime, where the transition between the
partially polar-ordered uniform phase and the TWA phase becomes supercritical.

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 $\sigma_{t}^{A}$ as traveling polar wave phase A (TWA) and the one in
the unstable region at $\sigma_{t}^{B}$ 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, $|\alpha_{1}^{\parallel}|$ and
$|\beta_{1}^{\parallel}|$, are shown in Fig. 6; the superscript for the
velocity fields indicates that we consider the component parallel to the
direction of macroscopic order, $\alpha_{1}^{(0)}$and $\beta_{1}^{(0)}$, and
wave vector, $\mathbf{q}$ (as mentioned in section III.3, the stability
analysis showed that instabilities are largest for
$\mathbf{q}\parallel\alpha_{1}^{(0)},\beta_{1}^{(0)}$, 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. (24) and (26). 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 $\bar{a},\bar{b},|\bar{\alpha}_{1}|$ and
$|\bar{\beta}_{1}|$ are approximately also given by the corresponding uniform
fixed point values.

Figure 6:  Typical wave profiles for the density fields (solid lines), $a$ and
$b$, and the velocity fields (dashed lines), $|\alpha_{1}^{\parallel}|$ and
$|\beta_{1}^{\parallel}|$, in the traveling wave phases TWA (a) and TWB (b).
Control parameters are set to relative fitness $\lambda\,{=}\,0.4$, game
strength $\tau\,{=}\,0.2$, and noise strength $\sigma\,{=}\,0.67$ in (a), and
$\sigma\,{=}\,0.47$ in (b). The insets show how the density and velocity
fields are related to each other: At every point in space, the velocity fields
$\alpha_{1}$ and $\beta_{1}$ (dashed lines in the insets) are approximately
equal to the values
$|\alpha_{1}(\mathbf{r},t)|\,{\approx}\,|\alpha_{1}^{(0)}(\mathbf{r},t)|\,{=}\,\sqrt{\mu_{1}^{A}(a,b)/\xi(a)}$
and
$|\beta_{1}(\mathbf{r},t)|\,{\approx}\,|\beta_{1}^{(0)}(\mathbf{r},t)|\,{=}\,\sqrt{\mu_{1}^{B}(a,b)/\xi(b)}$
corresponding to spatially uniform solutions (see Eqs. (26)) for the
respective local densities $a({\bf r},t)$ and $b({\bf 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).

Moreover, we find that at every point in space, the velocity fields,
$|\alpha_{1}(\mathbf{r},t)|$ and $|\beta_{1}(\mathbf{r},t)|$, are well
approximated by the local values of the spatially uniform solutions, Eq. (26).
For local densities $a(\mathbf{r},t)$ and $b(\mathbf{r},t)$ with
$\mu_{1}^{A/B}\,{>}\,0$, these are given by
$|\alpha_{1}^{(0)}(\mathbf{r},t)|\,{=}\,\sqrt{\mu_{1}^{A}(a,b)/\xi(a)}$ and
$|\beta_{1}^{(0)}(\mathbf{r},t)|\,{=}\,\sqrt{\mu_{1}^{B}(a,b)/\xi(b)}$, else
we have $|\alpha_{1}|,|\beta_{1}|\,{=}\,0$. Hence, non-zero values for the
velocity fields $\alpha_{1}$ and $\beta_{1}$ requires that the densities are
large enough such that locally one is above the threshold for polar order,
i.e. both $\mu_{1}^{A}(a,b)$ and $\mu_{1}^{B}(a,b)$ are positive (see section
III.2). The only exception are patterns of the $B$ species in the TWA phase in
a regime where both $\lambda$ and $\tau$ are small enough such that the
density $b$ falls below the threshold value for polar order,
$\mu_{1}^{B}(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 ($|\alpha_{1}|>0$ or
$|\beta_{1}|>0$) and low density regions that are disordered
(($|\alpha_{1}|\,{=}\,0$ or $|\beta_{1}|\,{=}\,0$)).

####  III.4.2 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 $(\lambda,\tau)$ and $\sigma>\sigma_{t}^{A}$ where there is no
macroscopic order, neither in $A$ nor in $B$. Then, we reduce the noise
strength $\sigma$ 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
$\sigma$. Afterwards, we increase the noise strength again quasi-statically
back to its initial value. During this $\sigma$-sweep, we monitor the spatial
average of the velocity fields, $|\bar{\alpha}_{1}|$ and $|\bar{\beta}_{1}|$.
To identify the presence of spatial patterns, we use parameters that we call
variation parameters $\eta_{A}$ and $\eta_{B}$, which measure spatial
variations in the densities $a$ and $b$. They are defined as

$\eta_{A}:=\overline{\left(\frac{|\Delta a|}{\Delta
x}\right)},\quad\text{and}\quad\eta_{B}:=\overline{\left(\frac{|\Delta
b|}{\Delta x}\right)},$ (29)

where $\frac{|\Delta a|}{\Delta 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 $\eta_{A}\,{=}\,0$
($\eta_{B}\,{=}\,0$), species $A$ ($B$) is in a spatially uniform state, while
a non-zero value $\eta_{A}\,{>}\,0$ ($\eta_{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:  Bifurcations of the spatial average of the velocity fields,
$|\bar{\alpha}_{1}|$ and $|\bar{\beta}_{1}|$, and the spatial variation
parameters, $\eta_{A}$ and $\eta_{B}$, as a function of the noise strength
$\sigma$ for a fixed game strength $\tau\,{=}\,0.2$ and for different values
of the relative fitness $\lambda$: (a) $\lambda\,{=}\,1$, (b)
$\lambda\,{=}\,0.2$, and (c) $\lambda\,{=}\,0.5$. Blue and red lines consist
of sets of dots that each indicate a stable state calculated numerically. The
dashed lines indicate the analytical results for the spatially uniform
amplitudes of the velocity fields, Eq. (26); compare with Fig. 3. Regions
shaded in blue and yellow indicate parameter regimes, in which the system is
bistable or the spatially uniform states are linearly unstable with respect to
spatially heterogeneous perturbations, respectively.

Figure 7 shows the spatial variation parameters, $\eta_{A}$ and $\eta_{B}$,
and spatial averages of the velocity fields, $|\bar{\alpha}_{1}|$ and
$|\bar{\beta}_{1}|$, obtained from $\sigma$-sweeps for different values of the
relative fitness $\lambda$ and a fixed game strength $\tau\,{=}\,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 ($\lambda\,{=}\,1$), species $A$ and $B$ show identical
behavior (Fig. 7a): Quasi-statically decreasing (grey arrows) and then
increasing (black arrows) the noise strength $\sigma$ one observes a
hysteresis loop. Both, the variation parameters, $\eta_{A}$ and $\eta_{B}$,
and the average velocities, $|\bar{\alpha}_{1}|$ and $|\bar{\beta}_{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. (26). 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 $\lambda$ ($\lambda\,{=}\,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 ($|\alpha_{1}|,|\beta_{1}|>0$), there is now an intermediate
phase, in which only species $A$ shows polar order, but species $B$ does not
($|\alpha_{1}|>0,|\beta_{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 $\lambda$ ($\lambda\,{=}\,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 $\eta_{A}$,
$\eta_{B}$ and $|\alpha_{1}|$, $|\beta_{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.

####  III.4.3 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.3 already indicated a
possible role of the game interaction in pattern formation. There we found
that for finite game strength ($\tau\,{>}\,0$) spatial perturbations in the
density and velocity fields of species $A$ and $B$ grow simultaneously, while
in the absence of game interactions ($\tau\,{=}\,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 $\tau$. Figure 8 shows profiles of
the densities, $a$ and $b$, and the velocity fields,
$|\alpha_{1}^{\parallel}|$ and $|\beta_{1}^{\parallel}|$, for different values
of the game strength $\tau$. The other two control parameters are set to
$\sigma\,{=}\,0.68$ and $\lambda\,{=}\,0.4$, such that the system is in the
TWA phase. For $\tau\,{=}\,1$, traveling polar waves emerge in both species,
as expected (Fig. 8a). For smaller values of the game strength $\tau$, we
observe a decrease of wave amplitudes in $b$ and $|\beta_{1}^{\parallel}|$ but
an increase of wave amplitudes in $a$ and $|\alpha_{1}^{\parallel}|$ (Fig.
8b). In general, we find that the smaller the value of $\tau$, the lower the
amplitude of waves in species $B$ and the higher the amplitudes in species
$A$. In the absence of a snowdrift interaction ($\tau\,{=}\,0$), the wave
pattern in species $B$ completely disappears and species $B$ is in a
disordered, spatially uniform state with
$b\,{=}\,\frac{\lambda}{1+\lambda}\,{\approx}\,0.29$ and
$|\beta_{1}^{\parallel}|\,{=}\,0$ throughout the system (Fig. 8c). We conclude
that the pattern in $B$ observed for $\tau\,{>}\,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 $\tau$ in Section
III.3 (see Fig. 4c). There we found that the $b$\- and
$\beta_{1}^{\parallel}$-components of the eigenvector pointing in the
direction of the instability gradually decreased with the game strength
$\tau$, such that the $a$\- and $\alpha_{1}^{\parallel}$-components remained
as the only non-zero contributions for $\tau\rightarrow 0$.

Figure 8:  Wave profiles for noise strength $\sigma\,{=}\,0.68$, relative
fitness $\lambda\,{=}\,0.4$ and different game strengths $\tau$ as indicated
in the graphs. The system is in the TWA phase. (a) For $\tau\,{=}\,1$, species
$A$ and species $B$ exhibit polar waves in the density fields, $a$ and $b$,
with polar-ordered wave peaks (non-zero velocity fields
$|\alpha_{1}^{\parallel}|$ and $|\beta_{1}^{\parallel}|$). (b) For
$\tau\,{=}\,0.1$, we observe an increase of wave amplitudes in $a$ and
$|\alpha_{1}^{\parallel}|$ as well as a decrease in wave amplitudes in $b$ and
$|\beta_{1}^{\parallel}|$. This effect is stronger the smaller the value of
the game strength $\tau$ is. (c) For $\tau\,{=}\,0$, the polar wave pattern in
$B$ has completely disappeared and species $B$ is in a disordered, spatially
uniform state.

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, $\tau\,{>}\,0$, the game interaction induces a
shift of the density of the $B$ species towards the value
$b\,{=}\,\lambda\,a$, given by the fixed point of the snowdrift game (Eq.
(24)). The magnitude of this shift (flow in the phase space of the game)
depends the game strength $\tau$. The stronger the game interaction, the
stronger the mutual influence of the two species. For large $\tau$, one
expects an adiabatic limit where $a$ and $b$ are strictly correlated by
$b/a\,{=}\,\lambda$; 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 $\tau$, the pattern amplitudes in $a$ and
$|\alpha_{1}^{\parallel}|$ decrease while they increase in $b$ and
$|\beta_{1}^{\parallel}|$. For vanishing game strength ($\tau\rightarrow 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 $\tau$ 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 [55, 56, 42]. In addition, at the onset of
macroscopic order, particles were observed to form polar-ordered high-density
waves rather than spatially uniform polar order [56, 53].

###  IV.1 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\,{\times}\,L$ with periodic boundary conditions. Each
particle, located at position $\mathbf{r}_{i}$, is either of species A or B,
and moves with a velocity
$\mathbf{v}_{i}\,{=}\,v_{0}\,\mathbf{e}_{\theta}^{i}$; the speed $v_{0}$ is
assumed to be constant and the direction vector reads
$\mathbf{e}_{\theta}^{i}\,{=}\,(\cos\theta_{i},\sin\theta_{i})$, where
$\theta_{i}$ denotes the angle enclosed with the $x$-axis.

The positions of all particles are updated in discrete time steps $dt$

$\mathbf{r}_{j}^{t+dt}=\mathbf{r}_{j}^{t}+\mathbf{v}_{j}^{t}\,dt.$ (30)

Particles move ballistically unless they change their direction
$\mathbf{e}_{\theta}^{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 $\mathcal{A}_{j}^{t}$ at time $t$.
$\mathcal{A}_{j}^{t}$ is defined as the disk around $\mathbf{r}_{j}^{t}$ with
radius $d_{\text{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 $\eta$.
We assume that $\eta$ is a random variable that is uncorrelated in time and
space and given by a flat distribution over the interval $(-h\pi,h\pi]$ with
$h\,{\in}\,[0,1]$. The standard deviation $\sigma$ of a flat distribution is
then given by $\sigma\,{=}\,h\pi/\sqrt{3}$. We will refer to $\sigma$ as the
noise strength. The direction of each particle at time $t\,{+}\,dt$ is
therefore given by

$\theta_{j}^{t+dt}=\textrm{arg}\biggl{(}\sum_{k\in\mathcal{A}_{j}^{t}}e^{i\theta_{k}{t}}\biggr{)}+\eta^{t}_{j}\,.$
(31)

Previous studies have shown that if the noise strength $\sigma$ is
sufficiently low, species begin to form macroscopic polar order, i.e.
particles move collectively over distances much larger than their alignment
range $d_{\text{align}}$ [55, 56] In order to measure polar order of each
species we define the polar order parameters

$\alpha_{1}^{t}=\frac{1}{N_{A}^{t}}\Bigl{|}\sum_{k=1}^{N_{A}^{t}}e^{i\theta_{k}^{t}}\Bigr{|},~{}\;\text{and}\quad\beta_{1}^{t}=\frac{1}{N_{B}^{t}}\Bigl{|}\sum_{k=1}^{N_{B}^{t}}e^{i\theta_{k}^{t}}\Bigr{|},$
(32)

where $N_{A}^{t}$ and $N_{B}^{t}$ denote the total particle number of species
A and B, respectively, at time $t$; $\alpha_{1}^{t}$ and $\beta_{1}^{t}$ range
from $0$ (no order) to $1$ (perfect alignment of all particles). The temporal
averages of these polar order parameters are given by

$\langle\alpha_{1}\rangle=\frac{1}{T}\sum\limits_{t=1}^{T}\alpha_{1}^{t},~{}\;\text{and}\quad\langle\beta_{1}\rangle=\frac{1}{T}\sum\limits_{t=1}^{T}\beta_{1}^{t}.$
(33)

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
$\mathcal{G}_{j}^{t}$ of particle $j$ at time $t$ by a disk around particle
$j$ with radius $d_{\text{game}}$, which we will refer to in the following as
the game range.

All particles $N_{j}^{t}$ 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_{A,j}^{t}$ and $N_{B,j}^{t}$, respectively.
Consider the case that particle $j$ belongs to species A: This A-particle will
be at a disadvantage in the competition with the B-particles if the relative
abundance, $N_{B,j}^{t}/N_{A,j}^{t}$, lies below a threshold $\lambda$, 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\,{\to}\,B$) with a rate
$\tau(a_{j}^{t}\lambda\,{-}\,b_{j}^{t})$. Here,
$a_{j}^{t}\,{=}\,N^{t}_{A,j}/N_{j}^{t}$ and
$b_{j}^{t}\,{=}\,N^{t}_{B,j}/N_{j}^{t}$ 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 $\lambda$, 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

$p_{j}^{t}=\begin{cases}\tau(a_{j}^{t}\lambda{-}b_{j}^{t})dt&\textrm{if
$j{\in}\text{A}$ and $\tau(a_{j}^{t}\lambda{-}b_{j}^{t})>0$}\\\
\tau(b_{j}^{t}{-}a_{j}^{t}\lambda)dt&\textrm{if $j{\in}\text{B}$ and
$\tau(a_{j}^{t}\lambda{-}b_{j}^{t})<0$}\\\ \hfil
0&\textrm{otherwise.}\end{cases}$ (34)

As in Sec. II, $\tau\,{\geq}\,0$ is interpreted as the game strength: With
increasing $\tau$ 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
$\lambda\,{\in}\,[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, 56, 55] 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_{c}^{t}$ and $b_{c}^{t}$ and polar order
parameters $\alpha_{1,c}^{t}$ and $\beta_{1,c}^{t}$ 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
$\nu_{A}^{t}$, defined as twice the standard deviation of polar order
parameter $\alpha_{1,c}^{t}$ over all channels:

$\nu_{A}^{t}=2\sqrt{\frac{1}{C}\sum_{c=1}^{C}\left(\alpha_{1,c}^{t}-\langle\alpha_{1}^{t}\rangle\right)^{2}}.$
(35)

The factor of $2$ is introduced such that $\nu_{A}^{t}\,{\in}\,[0,1]$. The
velocity variation $\nu_{B}^{t}$ 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
$\Delta\nu_{A}\,{=}\,\langle\nu_{A,\perp}\rangle\,{-}\,\langle\nu_{A,\parallel}\rangle$,
with channels perpendicular and parallel to the direction of macroscopic
order, respectively. Further details can be found in Appendix D.

###  IV.2 Results of the agent-based simulation

Of the eight model parameters ($\lambda$, $\tau$, $\rho$, $\sigma$, $v_{0}$,
$d_{\text{align}}$, $d_{\text{game}}$, and $L$) we focus mainly on analyzing
the role of two: (i) the noise strength $\sigma$, and (ii) the relative
fitness $\lambda$. The noise strength $\sigma$ is used as a control parameter
to distinguish between the different phases, with the density kept fixed at
$\rho\,{=}\,N/L^{2}\,{=}\,1.5$. Varying the relative fitness $\lambda$ 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].

Figure 9: Regimes of different macroscopic order and spatio-temporal patterns
obtained from agent-based simulations for game strength $\tau\,{=}\,0.4$ and
two distinct values of the game theory interaction radius $d_{\text{game}}$ as
indicated in the graph. Yellow areas mark parameter regimes with traveling
polar wave patterns, where TWA and TWB mark the regions where species A and B,
respectively, promote wave formation. These two phases merge with increasing
values of the relative fitness $\lambda$. 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 phases — a 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.2, 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 $d_{\text{game}}$. The vertical orange lines in the right
panel indicate the three slices in parameter space displayed in Fig. 10.

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 $\tau$, relative fitness $\lambda$,
and noise strength $\sigma$, we observe up to five different phases. For
neutral games ($\lambda\,{=}\,1$), we observe, similar to previous studies
[42, 56, 55], three different phases: a disordered, a traveling polar wave,
and a uniformly polar-ordered phase. Figure 9 shows the diagram for a game
strength $\tau\,{=}\,0.4$, and two distinct game ranges
$d_{\text{game}}\,{\in}\,\\{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 the game range 111This is the case as
long as the game range is not chosen too small, i.e. that the average number
of particles participating in an interaction is insufficient to roughly
achieve the desired ratio given by $\lambda$.. In the following discussion,
the specific parameter values given refer to the case
$d_{\text{game}}\,{=}\,10$.

The agent-based simulations show that above some threshold value for the noise
strength ($\sigma\,{\gtrsim}\,0.7$) both species are in a spatially uniform,
disordered state, with the polar order parameters, $\alpha_{1}$ and
$\beta_{1}$, and the velocity variation $\Delta\nu$ 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 $\lambda$, as can be inferred from Fig. 10a
where $\alpha_{1}$ and $\Delta\nu_{A}$ are shown as functions of $\sigma$. 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 $\lambda$. 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=\lambda\cdot
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,
$\lambda\lesssim 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
$\lambda$ (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.4.3). A movie
displaying the dynamics in the TWA phase can be found in the supplementary
material (movie2).

Upon further decreasing the noise strength $\sigma$, there is a parameter
domain for $\lambda\lesssim 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 ($\alpha_{1}\,{>}\,0$), species B no longer
exhibits polar order ($\beta_{1}\,{\approx}\,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.

Figure 10:  Polar order (solid line) and velocity variation (dashed line)
against noise strength $\sigma$ for species A (top row) and species B (bottom
row) and a set of different values for the relative fitness $\lambda$ (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
$\lambda\,{=}\,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
($\alpha_{1}\,{>}\,0$) and finite velocity variations
($\Delta\nu_{A}\,{>}\,0$), while there is neither significant polar order in
species B ($\beta_{1}\,{\approx}\,0$) nor velocity variations $\Delta\nu_{B}$.
In contrast, the TWB phase is characterized by polar order as well as velocity
variation in both species A and B. The case $\lambda\,{=}\,0.28$ (b) is
phenomenologically similar to (a), but in the TWA phase there is now also
polar order ($\beta_{1}\,{>}\,0$) and finite velocity variations
($\Delta\nu_{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 $\lambda=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.

Upon decreasing the noise strength $\sigma$ 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 polar-ordered 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
agent-based 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 ($\lambda\,{=}\,0.26$).

The size of these demixed regions grows linearly with increasing game range
$d_{\text{game}}$, as is illustrated in Fig. 15. If demixing occurs, the
ensuing density ratio $b/a$ becomes larger than the relative fitness
$\lambda$. 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 $\tau$ 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.

Figure 11:  Density and velocity profiles in the TWB phase for noise strength
$\sigma\,{=}\,0.32$ and relative fitness $\lambda\,{=}\,0.26$. (a) Snapshot
from the agent-based simulation with particle positions of species A and B
marked in blue and red, respectively. Inset: blow-up showing species demixing
inside the wave. (b) 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 $\rho$ (upper plot) and polar
order $\alpha_{1}$, $\beta_{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.

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=\lambda$ 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 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 $\sigma$ than in the hydrodynamic approach,
and becomes even larger with increasing relative fitness levels $\lambda$; 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.

###### Acknowledgements.

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.

###  A.1 Extended kinetic Boltzmann equations

As discussed in Section II, we model the dynamics of two species $A$ and $B$
with self-propulsion and self-diffusion, 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
$\alpha\left(\mathbf{r},\theta,t\right)$ and
$\beta\left(\mathbf{r},\theta,t\right)$, for the species $A$ and $B$,
respectively:

$\displaystyle\partial_{t}\alpha\left(\mathbf{r},\theta,t\right)$
$\displaystyle+v_{0}\left(\mathbf{e}_{\theta}\cdot\nabla\right)\alpha\left(\mathbf{r},\theta,t\right)=I_{\text{diff}}\left[\alpha\right]+I_{\text{coll}}\left[\alpha,\alpha\right]+I_{\text{game}}^{A}\left[\alpha,\beta\right]$
(36a) $\displaystyle\partial_{t}\beta\left(\mathbf{r},\theta,t\right)$
$\displaystyle+v_{0}\left(\mathbf{e}_{\theta}\cdot\nabla\right)\beta\left(\mathbf{r},\theta,t\right)=I_{\text{diff}}\left[\beta\right]+I_{\text{coll}}\left[\beta,\beta\right]+I_{\text{game}}^{B}\left[\beta,\alpha\right]$
(36b)

with the terms describing self-diffusion, $I_{\text{diff}}\left[f\right]$, and
alignment interaction, $I_{\text{coll}}\left[f,f\right]$

$\displaystyle I_{\text{diff}}\left[f\right]$ $\displaystyle=-\omega
f\left(\theta\right)+\omega\int_{-\pi}^{\pi}d\theta^{\prime}\int_{-\infty}^{\infty}d\eta_{0}\,f\left(\theta^{\prime}\right)\
\mathcal{P}_{0}\left(\eta_{0}\right)\,\delta_{2\pi}\left(\theta^{\prime}+\eta_{0}-\theta\right),$
(37a) $\displaystyle I_{\text{coll}}\left[f,f\right]$
$\displaystyle=-f\left(\theta\right)\int_{-\pi}^{\pi}d\theta^{\prime}\,\mathcal{R}_{\theta,\theta^{\prime}}f\left(\theta^{\prime}\right)+\int_{-\pi}^{\pi}d\theta_{1}\int_{-\pi}^{\pi}d\theta_{2}\,\mathcal{R}_{\theta_{1},\theta_{2}}f\left(\theta_{1}\right)f\left(\theta_{2}\right)\int_{-\infty}^{\infty}d\eta\,P\left(\eta\right)\,\delta_{2\pi}\left(\tfrac{\theta_{1}+\theta_{2}}{2}{+}\eta{-}\theta\right),$
(37b)

and the terms for the game interaction

$\displaystyle I_{\text{game}}^{A}\left[\alpha,\beta\right]$
$\displaystyle=\tau\left(b-\lambda a\right)b\,\alpha\left(\theta\right),$
(38a) $\displaystyle I_{\text{game}}^{B}\left[\alpha,\beta\right]$
$\displaystyle=\tau\left(\lambda a-b\right)a\,\beta\left(\theta\right)\,.$
(38b)

The collision kernel
$\mathcal{R}_{\theta_{1},\theta_{2}}\,{=}\,4\,d_{0}v_{0}\left|\sin\left(\tfrac{1}{2}(\theta_{1}-\theta_{2})\right)\right|$
characterizes the collisions between the agents and

$\delta_{2\pi}\left(x\right):=\sum_{m=-\infty}^{+\infty}\delta(x+2\pi m).$
(39)

denotes the $2\pi$ periodic delta distribution. We assume that both random
variables, $\eta_{0}$ and $\eta$, describing diffusion and collision noise,
respectively, are Gaussian distributed:

$\displaystyle\mathcal{P}_{0}(\eta_{0})$
$\displaystyle=\frac{1}{\sqrt{2\pi\sigma_{0}^{2}}}\,\exp\biggl{[}-\frac{\eta_{0}^{2}}{2\sigma_{0}^{2}}\biggr{]},$
(40) $\displaystyle\mathcal{P}(\eta)$
$\displaystyle=\frac{1}{\sqrt{2\pi\sigma^{2}}}\,\exp\biggl{[}-\frac{\eta^{2}}{2\sigma^{2}}\biggr{]}.$
(41)

For a more detailed derivation of $I_{\text{diff}}$ and $I_{\text{coll}}$, we
refer to Ref. [42]. By using the rescaling

$\begin{split}\tilde{t}:=&\,t\cdot\omega,\qquad\tilde{\mathbf{r}}:=\,\mathbf{r}\cdot\omega\,v_{0}^{-1},\\\
\tilde{\tau}:=&\,\tau\cdot\omega\,d_{0}^{-2}v_{0}^{-2},\qquad\tilde{f}:=\,f\cdot
d_{0}v_{0}\omega^{-1}.\end{split}$ (42)

we can set $\omega$, $v_{0}$, and $d_{0}$ equal to $1$ without loss of
generality. In the following, we will omit the tilde and just write $t$,
$\mathbf{r}$, $\tau$, and $f$ to simplify notation.

###  A.2 Fourier transformed Boltzmann terms

As a first step, we perform the Fourier transformation of the extended kinetic
Boltzmann equations (36) in $\theta$ with the Fourier representation

$\displaystyle f\left(\mathbf{r},\theta,t\right)$
$\displaystyle=\frac{1}{2\pi}\sum\limits^{\infty}_{k=-\infty}e^{-ik\theta}\,\hat{f}_{k}(\mathbf{r},t),$
(43) $\displaystyle\hat{f}_{k}(\mathbf{r},t)$
$\displaystyle=\int\limits^{\pi}_{-\pi}d\theta~{}e^{ik\theta}\,f\left(\mathbf{r},\theta,t\right),$
(44)

where $f$ stands for $\alpha$ or $\beta$. With this convention for the Fourier
transformation the following identities hold:
$\int_{-\pi}^{\pi}e^{ik\theta}=2\pi\delta_{k,0}$ and
$\sum_{k=-\infty}^{\infty}e^{-ik\theta}=2\pi\delta(\theta)$. In the following,
we omit spatial and temporal dependencies, and, if not specified differently,
all integrals over $\theta$ have boundaries $-\pi$ to $\pi$, all integrals
over the noise have boundaries $\pm\infty$, and sums over wave numbers $k$ run
from $-\infty$ to $\infty$. In order to calculate the Fourier transform of the
convection term, we use
$\mathbf{e}_{\theta}=\left(\cos\theta,~{}\sin\theta\right)^{T}$. In Fourier
space, the convection term reads

$\displaystyle\int\\!d\theta~{}e^{ik\theta}\,\mathbf{e}_{\theta}\cdot\partial_{\mathbf{r}}\alpha(\theta)=\int
d\theta~{}e^{ik\theta}\Bigl{[}\cos\theta\,\partial_{x}+\sin\theta\,\partial_{y}\Bigr{]}\alpha(\theta)$
$\displaystyle\qquad=\int\\!\mathrm{d}\theta\,\mathrm{e}^{\mathrm{i}k\theta}\,\tfrac{1}{2}\left[\left(\mathrm{e}^{\mathrm{i}\theta}{+}\mathrm{e}^{-\mathrm{i}\theta}\right)\partial_{x}-\mathrm{i}\left(\mathrm{e}^{\mathrm{i}\theta}{-}\mathrm{e}^{-\mathrm{i}\theta}\right)\partial_{y}\right]\alpha(\theta)$
$\displaystyle\qquad=\tfrac{1}{2}\partial_{x}\left(\hat{\alpha}_{k+1}+\hat{\alpha}_{k-1}\right)-\tfrac{i}{2}\partial_{y}\left(\hat{\alpha}_{k+1}-\hat{\alpha}_{k-1}\right)$
$\displaystyle\qquad=\tfrac{1}{2}\bigl{(}\nabla\hat{\alpha}_{k-1}+\nabla^{*}\hat{\alpha}_{k+1}\bigr{)},$
(45)

where $\nabla\,{:=}\,\partial_{x}{+}i\partial_{y}$ and
$\nabla^{*}\,{:=}\,\partial_{x}{-}i\partial_{y}$. The Fourier transform of the
diffusion term can be calculated as

$\displaystyle\int\\!d\theta~{}e^{ik\theta}I_{\text{diff}}\left[\alpha\right]=-\int\\!d\theta~{}e^{ik\theta}\alpha(\theta)$
$\displaystyle\quad+\int\\!d\theta~{}e^{ik\theta}\int\\!d\theta^{\prime}\\!\int\\!d\eta~{}\mathcal{P}_{0}(\eta)\,\alpha(\theta^{\prime})\,\delta_{2\pi}(\theta{-}\theta^{\prime}{-}\eta)$
$\displaystyle\quad=-\hat{\alpha}_{k}+\int
d\theta^{\prime}\int\\!d\eta~{}e^{ik(\theta^{\prime}+\eta)}\mathcal{P}_{0}(\eta)\,\alpha(\theta^{\prime})$
$\displaystyle\quad=-\hat{\alpha}_{k}+\hat{\alpha}_{k}\hat{\mathcal{P}}_{0,k}=(\hat{\mathcal{P}}_{0,k}-1)\hat{\alpha}_{k}$
(46)

with $\hat{\mathcal{P}}_{0,k}=$ exp$(-k^{2}\sigma_{0}^{2}/2)$ the $k$-th
Fourier mode of the Gaussian distribution $\mathcal{P}_{0}(\eta)$. The
diffusion term consists of a gain and a loss term, given by the first and the
second term in Eq. (37)b, respectively. We denote the Fourier transforms of
the gain and the loss terms as $(I)$ and $(II)$, respectively, $\int
d\theta~{}e^{ik\theta}I_{\text{coll}}\left[\alpha,\beta\right]=(I)+(II)$:

$\displaystyle(I)=$ $\displaystyle-\int
d\theta\int\limits_{\theta-\pi}^{\theta+\pi}d\theta^{\prime}~{}e^{ik\theta}~{}\mathcal{R}_{\theta,\theta^{\prime}}\,\alpha(\theta)\,\beta(\theta^{\prime})$
$\displaystyle=-\int
d\theta\int\limits_{\theta-\pi}^{\theta+\pi}d\theta^{\prime}~{}e^{ik\theta}\mathcal{R}_{\theta,\theta^{\prime}}\sum\limits_{q,p}\frac{\hat{\alpha}_{q}\hat{\beta}_{p}}{(2\pi)^{2}}e^{-iq\theta}e^{-ip\theta^{\prime}}$
$\displaystyle=-\int d\theta\int
d\vartheta\sum\limits_{q,p}\frac{\hat{\alpha}_{q}\hat{\beta}_{p}}{(2\pi)^{2}}\mathcal{R}_{|\vartheta|}\,e^{i\theta(k-q-p)}e^{-ip\vartheta}$
$\displaystyle=-\frac{1}{2\pi}\int
d\vartheta~{}\sum\limits_{q,p}\hat{\alpha}_{q}\hat{\beta}_{p}\mathcal{R}_{|\vartheta|}\,\delta_{q,k-p}\,e^{-ip\vartheta}$
$\displaystyle=-\frac{1}{2\pi}\sum\limits_{p}\hat{\alpha}_{k-p}\hat{\beta}_{p}\int
d\vartheta~{}\cos(p\vartheta)\mathcal{R}_{|\vartheta|}.$ (47)

Analogously, $(II)$ can be calculated as

$\begin{split}(II)=&\int d\theta~{}e^{ik\theta}\int
d\theta_{1}\int\limits_{\theta_{1}-\pi}^{\theta_{1}+\pi}d\theta_{2}\int\\!d\eta~{}\alpha(\theta_{1})\beta(\theta_{2})\mathcal{R}_{\theta_{1},\theta_{2}}\mathcal{P}(\eta)\,\delta_{2\pi}\left(\theta{-}\tfrac{\theta_{1}{+}\theta_{2}}{2}{-}\eta\right)\\\
=&\int d\theta\int
d\theta_{1}\int\limits_{\theta_{1}-\pi}^{\theta_{1}+\pi}d\theta_{2}\int
d\eta~{}\sum\limits_{q,p}e^{ik\theta}e^{-iq\theta_{1}}e^{-ip\theta_{2}}\frac{\hat{\alpha}_{q}\hat{\beta}_{p}}{(2\pi)^{2}}\mathcal{R}_{\theta_{1},\theta_{2}}\mathcal{P}(\eta)\,\delta_{2\pi}\left(\theta{-}\tfrac{\theta_{1}{+}\theta_{2}}{2}{-}\eta\right)\\\
=&\int d\theta_{1}\int
d\vartheta~{}\sum\limits_{q,p}e^{i\theta_{1}(k-q-p)}e^{i\vartheta\left(k/2-p\right)}\frac{\hat{\alpha}_{q}\hat{\beta}_{p}}{(2\pi)^{2}}\mathcal{R}_{|\vartheta|}\hat{\mathcal{P}}_{k}=\frac{1}{2\pi}\int
d\vartheta~{}\sum\limits_{q,p}\delta_{q,k-p}e^{i\vartheta\left(k/2-p\right)}\hat{\alpha}_{q}\hat{\beta}_{p}\mathcal{R}_{|\vartheta|}\hat{\mathcal{P}}_{k}\\\
=&\frac{1}{2\pi}\sum\limits_{p}\hat{\alpha}_{k-p}\hat{\beta}_{p}\hat{\mathcal{P}}_{k}\int
d\vartheta~{}\cos\Bigl{(}\bigl{(}\tfrac{k}{2}-p\bigr{)}\vartheta\Bigr{)}\mathcal{R}_{|\vartheta|},\end{split}$
(48)

In total, we get for the collision term:

$\displaystyle\int
d\theta~{}e^{ik\theta}I_{\text{coll}}\left[\alpha,\beta\right]=\sum\limits_{p}\hat{\alpha}_{k-p}\hat{\beta}_{p}\mathcal{I}_{p,k},$
(49)

with

$\displaystyle\mathcal{I}_{p,k}:=\int\frac{d\vartheta}{2\pi}\mathcal{R}_{|\vartheta|}\Bigl{[}\hat{\mathcal{P}}_{k}\,\cos\bigl{(}(\tfrac{k}{2}-p)\vartheta\bigr{)}-\cos(p\vartheta)\Bigr{]}.$
(50)

The last remaining part to complete the Fourier representation of Eq. (36) is
$I_{\text{game}}$:

$\int d\theta~{}e^{ik\theta}I^{A}_{\text{game}}=\tau(b-\lambda
a)b\,\hat{\alpha}_{k}$ (51)

and $I^{B}_{\text{game}}$ analogously. The Fourier transformation of the
equation for species B, Eq. (36b) works analogous to the transformation for
species A shown here. In summary, the Fourier representations of Eqs. (36)
read

$\displaystyle\partial_{t}\hat{\alpha}_{k}+\frac{1}{2}$
$\displaystyle\left(\nabla\hat{\alpha}_{k-1}+\nabla^{*}\hat{\alpha}_{k+1}\right)$
(52a)
$\displaystyle=(\hat{\mathcal{P}}_{0,k}-1)\hat{\alpha}_{k}+\sum\limits_{p}\hat{\alpha}_{k-p}\hat{\alpha}_{p}\mathcal{I}_{p,k}+\tau(b-\lambda
a)b\,\hat{\alpha}_{k}$ $\displaystyle\partial_{t}\hat{\beta}_{k}+\frac{1}{2}$
$\displaystyle\left(\nabla\hat{\beta}_{k-1}+\nabla^{*}\hat{\beta}_{k+1}\right)$
(52b)
$\displaystyle=(\hat{\mathcal{P}}_{0,k}-1)\hat{\beta}_{k}+\sum\limits_{p}\hat{\beta}_{k-p}\hat{\beta}_{p}\mathcal{I}_{p,k}-\tau(b-\lambda
a)a\hat{\beta}_{k}.$

where we used the fact that the 0-th Fourier mode is the local density of the
respective species, since

$\displaystyle\alpha_{0}(\mathbf{r},t)$
$\displaystyle=\int_{-\pi}^{\pi}\mathrm{d}\theta\,e^{ik\cdot
0}\alpha(\mathbf{r},\theta,t)$ (53)
$\displaystyle=\int_{-\pi}^{\pi}\mathrm{d}\theta\,\alpha(\mathbf{r},\theta,t)=a(\mathbf{r},t)$
(54)

and $\beta_{0}(\mathbf{r},t)=b(\mathbf{r},t)$ analogously.

###  A.3 Uniform disordered solution

One solution of (52) can be identified by using the fact that
$(\hat{\mathcal{P}}_{0,k}-1)=0$ and $\mathcal{I}_{n,k}=0$ for $k=0$ for all
$n$. Assuming spatial homogeneity, we set
$\nabla\alpha_{k}\,{=}\,\nabla\beta_{k}\,{=}\,0$. The kinetic Boltzmann
equations for $k\,{=}\,0$ then read

$\partial_{t}a=\tau\left(b-\lambda
a\right)a\,b,\quad\partial_{t}b=-\tau\left(b-\lambda a\right)a\,b.$ (55)

Substituting the control parameter $\bar{\rho}\,{=}\,a+b$, we get the three
solutions

$\displaystyle(a^{(0)},b^{(0)})=$ $\displaystyle\hskip
18.77882pt(0,\bar{\rho})$ (56a) $\displaystyle(a^{(0)},b^{(0)})=$
$\displaystyle\hskip 18.77882pt(\bar{\rho},0)$ (56b)
$\displaystyle(a^{(0)},b^{(0)})=$
$\displaystyle(\bar{\rho}\,\tfrac{1}{1+\lambda},\bar{\rho}\,\tfrac{\lambda}{1+\lambda})$
(56c)

and $\alpha_{k}\,{=}\,\beta_{k}\,{=}\,0$ for all $k\,{>}\,0$. Eqs. (55) 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 (52) around the solution Eq. (56c) and calculate the linear
growth rate of small perturbations $\delta\alpha_{k}$ and $\delta\beta_{k}$.
The linearized dynamics of $\delta\alpha_{k}$ and $\delta\beta_{k}$ is given
by (for $k\,{>}\,0$)

$\displaystyle\partial\delta\alpha_{k}$
$\displaystyle=\left[\hat{\mathcal{P}}_{0,k}-1+(\mathcal{I}_{k,k}+\mathcal{I}_{0,k})\,\bar{\rho}\,\frac{1}{1+\lambda}\right]\delta\alpha_{k}$
$\displaystyle=\mu_{k}^{A}(a^{(0)},b^{(0)})\delta\alpha_{k},$ (57)
$\displaystyle\partial\delta\beta_{k}$
$\displaystyle=\left[\hat{\mathcal{P}}_{0,k}-1+(\mathcal{I}_{k,k}+\mathcal{I}_{0,k})\,\bar{\rho}\,\frac{\lambda}{1+\lambda}\right]\delta\beta_{k}$
$\displaystyle=\mu_{k}^{B}(a^{(0)},b^{(0)})\delta\beta_{k},$ (58)

with $\mu_{k}^{A}$ and $\mu_{k}^{B}$ defined as in Eqs. (21). The solution Eq.
(56c) is therefore unstable when the linear growth rates $\mu_{k}^{A}$ or
$\mu_{k}^{B}$ have positive real parts. Because $\hat{\mathcal{P}}_{0,k}$ is
the Fourier transformed distribution of the noise, we have
$\hat{\mathcal{P}}_{0,k}\,{\leq}\,1$, thus
$\hat{\mathcal{P}}_{0,k}-1\,{\leq}\,0$. The linear growth rates $\mu_{k}^{A}$
and $\mu_{k}^{B}$ can therefore only be positive when
$\mathcal{I}_{k,k}+\mathcal{I}_{0,k}\,{>}\,0$. Calculating
$\left(\mathcal{I}_{k,k}+\mathcal{I}_{0,k}\right)$ for different values of
$k\,{>}\,0$ shows that it is positive only for $k\,{=}\,1$ (for collision
noise values $\sigma<1$) and negative for all $k\,{>}\,1$. Hence, small
perturbations of the modes $\alpha_{k},\beta_{k}$ for $k\,{>}\,1$ will decay
and only perturbations in $\alpha_{1}$ and $\beta_{1}$ may grow (if
$\mu_{k}^{A}\,{>}\,0$ or $\mu_{k}^{B}\,{>}\,0$, respectively).

The definition of the Fourier series yields for the first mode, $k\,{=}\,1$:

$\displaystyle\alpha_{1}(\mathbf{r},t)$
$\displaystyle=\int_{-\pi}^{\pi}d\theta\,e^{i\theta}\alpha(\mathbf{r},\theta,t),$
(59) $\displaystyle\beta_{1}(\mathbf{r},t)$
$\displaystyle=\int_{-\pi}^{\pi}d\theta\,e^{i\theta}\beta(\mathbf{r},\theta,t).$
(60)

The first modes can be related to the polar order fields $\mathbf{P}_{A}$ and
$\mathbf{P}_{B}$, which measure the polar order of the system, through

$\displaystyle
a\,\mathbf{P}_{A}=\left(\begin{array}[]{c}\text{Re}(\alpha_{1})\\\
\text{Im}(\alpha_{1})\end{array}\right)\,,\quad
b\,\mathbf{P}_{B}=\left(\begin{array}[]{c}\text{Re}(\beta_{1})\\\
\text{Im}(\beta_{1})\end{array}\right).$ (65)

Hence, the real and imaginary part of $\alpha_{1}$ and $\beta_{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 $\alpha_{1}$ and $\beta_{1}$ as
velocity fields of the species $A$ and $B$, respectively. This means that when
$\mu_{k}^{A}$ or $\mu_{k}^{B}$ 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.2, we find critical values for $\sigma$ (as functions
of $\lambda$) where either the real part of $\mu_{k}^{A}$ or $\mu_{k}^{B}$ is
zero, indicating a transition from zero to nonzero polar order.

###  A.4 Scaling ansatz

Following [41, 42, 49] we employ a scaling ansatz to derive hydrodynamic
equations for the density, $a$ and $b$, and velocity fields $\alpha_{1}$ and
$\beta_{1}$, of the species $A$ and $B$, respectively: Specifically, we
introduce the small parameter $\varepsilon$ and make the ansatz

$\begin{split}a=\mathcal{O}(\varepsilon)~{},\quad\alpha_{k}=&\mathcal{O}(\varepsilon^{|k|})\quad\forall
k\neq 0\\\
b=\mathcal{O}(\varepsilon)~{},\quad\beta_{k}=&\mathcal{O}(\varepsilon^{|k|})\quad\forall
k\neq 0\end{split}$ (66)

Furthermore, we assume the scaling

$\partial_{t}~{}\sim~{}\partial_{x}~{}\sim~{}\partial_{y}~{}\sim~{}\varepsilon.$
(67)

We can now expand the Fourier representation of the kinetic Boltzmann
equations (52) up to the lowest order which gives non-trivial solutions, which
is $\varepsilon^{3}$ in our case and truncate higher orders. The remaining
terms are:

$\displaystyle k=0:$
$\displaystyle\partial_{t}a+\tfrac{1}{2}(\nabla\alpha_{1}^{*}+\nabla^{*}\alpha_{1})=\tau(b-\lambda
a)a\,b$ (68a) $\displaystyle k=1:$
$\displaystyle\partial_{t}\alpha_{1}+\tfrac{1}{2}(\nabla\alpha_{0}+\nabla^{*}\alpha_{2})=(\mathcal{P}_{0,1}-1)\alpha_{1}+(\mathcal{I}_{0,1}+\mathcal{I}_{1,1})a\,\alpha_{1}+(\mathcal{I}_{2,1}+\mathcal{I}_{-1,1})\alpha_{1}^{*}\alpha_{2}+\tau(b-\lambda
a)\alpha_{1}b$ (68b) $\displaystyle k=2:$
$\displaystyle\partial_{t}\alpha_{2}+\tfrac{1}{2}\nabla\alpha_{1}=(\mathcal{P}_{0,2}-1)\alpha_{2}+(\mathcal{I}_{0,2}+\mathcal{I}_{2,2})a\,\alpha_{2}+\mathcal{I}_{1,2}\,\beta_{1}^{2}$
(68c)

We only show the equations for $a$, $\alpha_{1}$ and $\alpha_{2}$ here. The
equations for $b$ and $\beta_{1}$ are analogue. Because we only take into
account $\mathcal{O}(\epsilon^{3})$ in our final equations for $a$, $b$,
$\alpha_{1}$, and $\beta_{1}$, we can drop the terms $\partial_{t}\alpha_{2}$
in (68c) and solve (68c) for $\alpha_{2}$. By defining

$\nu(a)=-\frac{1}{4[P_{0,2}-1+(\mathcal{I}_{0,2}+\mathcal{I}_{2,2})a]}$ (69)

we can write
$\alpha_{2}=4\,\nu(a)\mathcal{I}_{1,2}\,\beta_{1}^{2}-2\,\nu(a)\nabla\alpha_{1}$.
Substituting this into (68b), we get

$\begin{split}\partial_{t}\alpha_{1}=&-\tfrac{1}{2}\nabla
a-2\,\nabla^{*}\nu(a)\left(\mathcal{I}_{1,2}\,\alpha^{2}_{1}-\tfrac{1}{2}\nabla\alpha_{1}\right)-2\,\nu(a)\left(2\mathcal{I}_{1,2}\,\alpha_{1}\nabla^{*}\alpha_{1}-\tfrac{1}{2}\nabla^{*}\nabla\alpha_{1}\right)+(\mathcal{P}_{0,1}-1)\alpha_{1}\\\
&+(\mathcal{I}_{0,1}+\mathcal{I}_{1,1})a\,\alpha_{1}+4\,\nu(a)\mathcal{I}_{1,2}(\mathcal{I}_{2,1}+\mathcal{I}_{-1,1})\alpha_{1}^{*}\alpha_{2}^{2}-2\,\nu(a)(\mathcal{I}_{2,1}+\mathcal{I}_{1,1})\alpha_{1}^{*}\nabla\alpha_{1}+\tau(b-\lambda
a)b\,\alpha_{1}\end{split}$ (70)

Neglecting order $\mathcal{O}(\varepsilon^{4})$ and rearranging the equation,
we arrive at

$\begin{split}\partial_{t}\alpha_{1}=&\big{[}(\mathcal{P}_{0,1}-1)+(\mathcal{I}_{0,1}+\mathcal{I}_{1,1})a+\tau(b-\lambda
a)b+4(\mathcal{I}_{2,1}+\mathcal{I}_{-1,1})|\alpha_{1}|^{2}\nu(a)\mathcal{I}_{1,2}\big{]}\alpha_{1}\\\
&-\tfrac{1}{2}\nabla
a+\nu(a)\nabla^{*}\nabla\alpha_{1}-4\nu(a)\mathcal{I}_{1,2}\alpha_{1}\nabla^{*}\alpha_{1}-2\nu(a)(\mathcal{I}_{2,1}+\mathcal{I}_{-1,1})\alpha_{1}^{*}\nabla\alpha_{1}.\end{split}$
(71)

Following the same steps for $\beta_{1}$ and with the definitions

$\displaystyle\mu_{1}^{A}$
$\displaystyle=\,\left(\mathcal{P}_{0,1}{-}1\right){+}\left(\mathcal{I}_{0,1}{+}\mathcal{I}_{1,1}\right)a{+}\tau\left(b{-}\lambda
a\right)b,\quad$ (72a) $\displaystyle\mu_{1}^{B}$
$\displaystyle=\,\left(\mathcal{P}_{0,1}{-}1\right){+}\left(\mathcal{I}_{0,1}{+}\mathcal{I}_{1,1}\right)b{+}\tau(\lambda
a{-}b)a$ (72b) $\displaystyle\nu(x)$
$\displaystyle=-\tfrac{1}{4}\left[\mathcal{P}_{0,2}-1+(\mathcal{I}_{0,2}+\mathcal{I}_{2,2})x\right]^{-1},$
(72c) $\displaystyle\gamma(x)$ $\displaystyle=4\,\nu(x)\mathcal{I}_{1,2}$
(72d) $\displaystyle\kappa(x)$
$\displaystyle=2\,\nu(x)(\mathcal{I}_{2,1}+\mathcal{I}_{-1,1}),$ (72e)
$\displaystyle\xi(x)$
$\displaystyle=-4(\mathcal{I}_{2,1}+\mathcal{I}_{-1,1})\nu(x)\mathcal{I}_{1,2},$
(72f)

where $x$ stands for either $a$ or $b$ and

$\displaystyle\nabla$
$\displaystyle=\partial_{x}+i\partial_{y},\quad\nabla^{*}=\partial_{x}-i\partial_{y},$
(73a) $\displaystyle\mathcal{P}_{0,k}$
$\displaystyle=\int\\!d\eta\,e^{ik\eta_{0}}\,\mathcal{P}_{0}(\eta_{0})=\exp\bigl{[}-\tfrac{1}{2}\sigma_{0}^{2}k^{2}\bigr{]},$
(73b) $\displaystyle\mathcal{P}_{k}$
$\displaystyle=\int\\!d\eta\,e^{ik\eta}\,\mathcal{P}(\eta)=\exp\bigl{[}-\tfrac{1}{2}\sigma^{2}k^{2}\bigr{]},$
(73c) $\displaystyle\mathcal{I}_{n,k}$
$\displaystyle=\int\\!d\psi\,\mathcal{R}_{|\psi|}\left[\mathcal{P}_{k}\cos\bigl{(}(n{-}\tfrac{k}{2})\psi\bigr{)}{-}\cos(n\psi)\right].$
(73d)

we get the final set of hydrodynamic equations for the densities $a$ and $b$
and the polar order fields $\alpha_{1}$ and $\beta_{1}$

$\displaystyle\partial_{t}a$
$\displaystyle=-\frac{1}{2}\left(\nabla\alpha_{1}^{*}+\nabla^{*}\alpha_{1}\right)+\tau(b-\lambda
a)b\,a$ (74a) $\displaystyle\partial_{t}b$
$\displaystyle=-\frac{1}{2}\left(\nabla\beta_{1}^{*}+\nabla^{*}\beta_{1}\right)+\tau(\lambda
a-b)a\,b$ (74b) $\displaystyle\partial_{t}\alpha_{1}$
$\displaystyle=\left(\mu_{1}^{A}(a,b)-\xi(a)|\alpha_{1}|^{2}\right)\alpha_{1}+\nu(a)\,\nabla^{*}\nabla\alpha_{1}$
(74c)
$\displaystyle-\gamma(a)\,\alpha_{1}\nabla^{*}\alpha_{1}-\kappa(a)\,\alpha_{1}^{*}\nabla\alpha_{1}-\tfrac{1}{2}\nabla
a$ (74d) $\displaystyle\partial_{t}\beta_{1}$
$\displaystyle=\left(\mu_{1}^{B}(a,b)-\xi(b)|u|^{2}\right)\beta_{1}+\nu(b)\,\nabla^{*}\nabla\beta_{1}$
(74e)
$\displaystyle-\gamma(b)\,\beta_{1}\nabla^{*}\beta_{1}-\kappa(b)\,\beta_{1}^{*}\nabla\beta_{1}-\tfrac{1}{2}\nabla
b.$ (74f)

The set of coupled partial differential equations (74) describes the evolution
of the two densities $a$ and $b$ and the corresponding polar order fields
$\alpha_{1}$ and $\beta_{1}$ for the species $A$ and $B$, respectively. When
we set the game parameter $\tau=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. (74) is a result of the truncation scheme: The
snowdrift terms in Eqs. (52) have the form $\tau(b-\lambda
a)b\,\hat{\alpha}_{k}$ in the dynamic equation for $\alpha_{k}$ and
$-\tau(b-\lambda a)a\,\hat{\beta}_{k}$ in the equation for $\beta_{k}$. Since
we assume a scaling
$a=\mathcal{O}(\varepsilon),b=\mathcal{O}(\varepsilon),\alpha_{k}=\mathcal{O}(\varepsilon^{|k|}),\beta_{k}=\mathcal{O}(\varepsilon^{|k|})$
for all $k\neq 0$, the snowdrift term is of order $\epsilon^{3}$ for $k=0$ and
$k=1$. For all $k\,{>}\,1$, the snowdrift term has a order of at least
$\epsilon^{4}$ and is therefore cut off in the truncation scheme. Hence, the
only snowdrift terms possible are the ones linear in the density and velocity
fields, hence $\pm\tau(b-\lambda a)b\,a$, $\tau(b-\lambda a)b\,\alpha_{1}$ and
$\tau(b-\lambda a)a\,\beta_{1}$, which all appear in the hydrodynamic Eqs.
(74).

##  Appendix B Jacobian of the hydrodynamic equations

To analyze the stability of our solutions against spatial perturbations, we
employ a plane wave ansatz:

$\displaystyle f(\mathbf{r},t)$
$\displaystyle=\frac{1}{2\pi}\int\limits_{-\infty}^{\infty}d\mathbf{q}~{}e^{i\mathbf{qr}}f(\mathbf{q}),$
(75a) $\displaystyle f(\mathbf{q})$
$\displaystyle=\frac{1}{2\pi}\int\limits_{-\infty}^{\infty}d\mathbf{r}~{}e^{-i\mathbf{qr}}f(\mathbf{r}),$
(75b) $\displaystyle f^{*}(\mathbf{r},t)$
$\displaystyle=\frac{1}{2\pi}\int\limits_{-\infty}^{\infty}d\mathbf{q}~{}e^{-i\mathbf{qr}}f^{*}(\mathbf{q}),$
(75c) $\displaystyle f^{*}(\mathbf{q})$
$\displaystyle=\frac{1}{2\pi}\int\limits_{-\infty}^{\infty}d\mathbf{r}~{}e^{i\mathbf{qr}}f^{*}(\mathbf{r}),$
(75d)

where $f$ stands for $\delta a$, $\delta b$, $\delta\alpha_{1}$ and
$\delta\beta_{1}$. Linearizing the hydrodynamic Eqs. (18) and Eqs. (22) around
the solutions given by Eqs. (24) and Eqs. (26) and substituting (75a), we
arrive at the following linear equations for $\delta a(\mathbf{q})$, $\delta
b(\mathbf{q})$, $\delta\alpha_{1}(\mathbf{q})$,
$\delta\alpha_{1}^{*}(-\mathbf{q})$, $\delta\beta_{1}(\mathbf{q})$ and
$\delta\beta_{1}^{*}(-\mathbf{q})$:

$\displaystyle\partial_{t}\delta a(\mathbf{q})=$
$\displaystyle\left[\tau(b-2\lambda a)b\right]\delta
a(\mathbf{q})+\left[\tau(2b-\lambda a)a\right]\delta
b(\mathbf{q})+\left[-\tfrac{1}{2}(iq_{x}+q_{y})\right]\delta\alpha_{1}(\mathbf{q})+\left[-\tfrac{1}{2}(iq_{x}-q_{y})\right]\delta\alpha_{1}^{*}(-\mathbf{q})$
(76a) $\displaystyle\partial_{t}\delta b(\mathbf{q})=$
$\displaystyle\left[-\tau(b-2\lambda a)b\right]\delta
a(\mathbf{q})+\left[-\tau(2b-\lambda a)a\right]\delta
b(\mathbf{q})+\left[-\tfrac{1}{2}(iq_{x}+q_{y})\right]\delta\beta_{1}(\mathbf{q})+\left[-\tfrac{1}{2}(iq_{x}-q_{y})\right]\delta\beta_{1}^{*}(-\mathbf{q})$
(76b)
$\displaystyle\begin{split}\partial_{t}\delta\alpha_{1}(\mathbf{q})=&\left[\left(\frac{\partial\mu_{1}^{A}}{\partial
a}-\frac{\partial\xi}{\partial
a}|\alpha_{1}|^{2}\right)\alpha_{1}-\frac{1}{2}(iq_{x}-q_{y})\right]\delta
a(\mathbf{q})+\left[\frac{\partial\mu_{1}^{A}}{\partial
b}\alpha_{1}\right]\delta b(\mathbf{q})\\\
&+\left[\mu_{1}^{A}-2\xi|\alpha_{1}|^{2}-\nu(q_{x}^{2}+q_{y}^{2})-\gamma\alpha_{1}(iq_{x}+q_{y})-\kappa\alpha_{1}^{*}(iq_{x}-q_{y})\right]\delta\alpha_{1}(\mathbf{q})+\left[-\xi\alpha_{1}^{2}\right]\delta\alpha_{1}^{*}(-\mathbf{q})\end{split}$
(76c)
$\displaystyle\begin{split}\partial_{t}\delta\alpha_{1}^{*}(-\mathbf{q})=&\left[\left(\frac{\partial\mu_{1}^{A}}{\partial
a}-\frac{\partial\xi}{\partial
a}|\alpha_{1}|^{2}\right)\alpha_{1}^{*}-\frac{1}{2}(iq_{x}+q_{y})\right]\delta
a(\mathbf{q})+\left[\frac{\partial\mu_{1}^{A}}{\partial
b}\alpha_{1}^{*}\right]\delta b(\mathbf{q})\\\
&+\left[-\xi(\alpha_{1}^{*})^{2}\right]\delta\alpha_{1}(\mathbf{q})+\left[\mu_{1}^{A}-2\xi|\alpha_{1}|^{2}-\nu(q_{x}^{2}+q_{y}^{2})-\gamma\alpha_{1}^{*}(iq_{x}-q_{y})-\kappa\alpha_{1}(iq_{x}+q_{y})\right]\delta\alpha_{1}^{*}(-\mathbf{q})\end{split}$
(76d)
$\displaystyle\begin{split}\partial_{t}\delta\beta_{1}(\mathbf{q})=&\left[\frac{\partial\mu_{1}^{B}}{\partial
a}\beta_{1}\right]\delta
a(\mathbf{q})+\left[\left(\frac{\partial\mu_{1}^{B}}{\partial
b}-\frac{\partial\xi}{\partial
b}|\beta_{1}|^{2}\right)\beta_{1}-\frac{1}{2}(iq_{x}-q_{y})\right]\delta
b(\mathbf{q})\\\
&+\left[\mu_{1}^{B}-2\xi|\beta_{1}|^{2}-\nu(q_{x}^{2}+q_{y}^{2})-\gamma\alpha_{1}(iq_{x}+q_{y})-\kappa\beta_{1}^{*}(iq_{x}-q_{y})\right]\delta\beta_{1}(\mathbf{q})+\left[-\xi\beta_{1}^{2}\right]\delta\beta_{1}^{*}(-\mathbf{q})\end{split}$
(76e)
$\displaystyle\begin{split}\partial_{t}\delta\beta_{1}^{*}(-\mathbf{q})=&\left[\frac{\partial\mu_{1}^{B}}{\partial
a}\beta_{1}^{*}\right]\delta
a(\mathbf{q})+\left[\left(\frac{\partial\mu_{1}^{B}}{\partial
b}-\frac{\partial\xi}{\partial
b}|\beta_{1}|^{2}\right)\beta_{1}^{*}-\frac{1}{2}(iq_{x}+q_{y})\right]\delta
b(\mathbf{q})\\\
&+\left[-\xi(\beta_{1}^{*})^{2}\right]\delta\beta_{1}(\mathbf{q})+\left[\mu_{1}^{B}-2\xi|\beta_{1}|^{2}-\nu(q_{x}^{2}+q_{y}^{2})-\gamma\beta_{1}^{*}(iq_{x}-q_{y})-\kappa\beta_{1}(iq_{x}+q_{y})\right]\delta\beta_{1}^{*}(-\mathbf{q})\end{split}$
(76f)

This system of linear equations for the Fourier modes
$\delta\mathbf{m}_{\mathbf{q}}\,{=}\,(\delta a(\mathbf{q}),\delta
b(\mathbf{q}),\delta\alpha_{1}(\mathbf{q}),\delta\alpha_{1}^{*}(\mathbf{q}),\delta\beta_{1}(\mathbf{q}),\delta\beta_{1}^{*}(\mathbf{q}))^{T}$
can be written as

$\partial_{t}\,\delta\mathbf{m}_{\mathbf{q}}=\mathbf{J}_{\mathbf{q}}\,\delta\mathbf{m}_{\mathbf{q}}$
(77)

where we introduced the Jacobian matrix $\mathbf{J}_{\mathbf{q}}$, which can
be read off from Eq. (76). As an example, for the homogeneous, isotropic fixed
point with $\alpha_{1}=\alpha_{1}^{*}=\beta_{1}=\beta_{1}^{*}=0$,
$a=\bar{\rho}/(1+\lambda)$ and $b=\bar{\rho}\lambda/(1+\lambda)$, $\mathbf{J}$
is given by:

$\mathbf{J}_{\mathbf{q}}=\begin{pmatrix}-\rho^{2}\tau\lambda^{2}/(1+\lambda)^{2}&\rho^{2}\tau\lambda/(1+\lambda)^{2}&-(iq_{x}+q_{y})/2&(q_{y}-iq_{x})/2&0&0\\\
\rho^{2}\tau\lambda^{2}/(1+\lambda)^{2}&-\rho^{2}\tau\lambda/(1+\lambda)^{2}&0&0&-(iq_{x}+q_{y})/2&(q_{y}-iq_{x})/2\\\
-(iq_{x}-q_{y})/2&0&\mu_{1}^{A}-\nu q^{2}&0&0&0\\\
-(iq_{x}+q_{y})/2&0&0&\mu_{1}^{A}-\nu q^{2}&0&0\\\
0&-(iq_{x}-q_{y})/2&0&0&\mu_{1}^{B}-\nu q^{2}&0\\\
0&-(iq_{x}+q_{y})/2&0&0&0&\mu_{1}^{B}-\nu q^{2}\end{pmatrix}$ (78)

with $q^{2}=q_{x}^{2}+q_{y}^{2}$. The eigenvalues of this matrix encode the
stability of the fixed point: If for some vector $\mathbf{q}$, the Jacobian
$\mathbf{J}_{\mathbf{q}}$ has an eigenvalue with positive real part, wave-like
perturbations with the wave-vector $\mathbf{q}$ will grow exponentially,
indicating the formation of patterns. On the other hand, if all eigenvalues
have negative real part, for all $\mathbf{q}$, spatial perturbations decay and
the spatially uniform state is (linearly) stable.

###  B.1 Velocity components of the eigenvector

Figure 12: Stability analysis of the hydrodynamic equations in
($\tau,\sigma$)-space. The color code measures how fast perturbations in the
velocity field of species $A$ ($\delta\alpha_{1}$) grow compared to
perturbations in the density of species $B$ ($\delta\beta_{1}$):
$\delta\alpha_{1}/(\delta\alpha_{1}+\delta\beta_{1})$ which ranges from $0$
(only perturbations $\delta\beta_{1}$ grow) to $1$ (only perturbations
$\delta\alpha_{1}$ grow).

In Sec. III.3 we compute $\delta a/(\delta a+\delta b)$, as a function of
$\tau$ (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, $\delta\alpha_{1}$, $\delta\alpha_{1}^{*}$, $\delta\beta_{1}$ and
$\delta\beta_{1}^{*}$ , encoding the perturbations of the velocity fields. An