Numerical treatment of the Boltzmann equation for self-propelled particle systems

Kinetic theories constitute one of the most promising tools to decipher the characteristic spatio-temporal dynamics in systems of actively propelled particles. In this context, the Boltzmann equation plays a pivotal role, since it provides a natural translation between a particle-level description of the system's dynamics and the corresponding hydrodynamic fields. Yet, the intricate mathematical structure of the Boltzmann equation substantially limits the progress toward a full understanding of this equation by solely analytical means. Here, we propose a general framework to numerically solve the Boltzmann equation for self-propelled particle systems in two spatial dimensions and with arbitrary boundary conditions. We discuss potential applications of this numerical framework to active matter systems, and use the algorithm to give a detailed analysis to a model system of self-propelled particles with polar interactions. In accordance with previous studies, we find that spatially homogeneous isotropic and broken symmetry states populate two distinct regions in parameter space, which are separated by a narrow region of spatially inhomogeneous, density-segregated moving patterns. We find clear evidence that these three regions in parameter space are connected by first order phase transitions, and that the transition between the spatially homogeneous isotropic and polar ordered phases bears striking similarities to liquid-gas phase transitions in equilibrium systems. Within the density segregated parameter regime, we find a novel stable limit-cycle solution of the Boltzmann equation, which consists of parallel lanes of polar clusters moving in opposite directions, so as to render the overall symmetry of the system's ordered state nematic, despite purely polar interactions on the level of single particles.


I. INTRODUCTION
Developing a deeper understanding of active matter [1][2][3][4] has been the major focus of a considerable amount of theoretical work over the last decades [5][6][7][8][9][10][11][12][13][14][15][16][17][18][19][20]. In recent years, kinetic theory has gained considerable popularity to assess the ordering behavior in systems of actively propelled particles [21][22][23][24][25][26][27][28][29][30][31][32][33][34][35][36]. In this context, the Boltzmann equation provides a particularly compelling approach to active matter systems. Apart from its inherent limitations, which are largely due the assumptions of binary particle interactions and molecular chaos (cf. Refs. [28,34,37]), the following advantages of this framework are manifest: (i) The structure of the Boltzmann equation, relating convection and collision processes on the level of the one-particle distribution function, is ideally suited to explicitly implement a microscopic picture of particle dynamics. (ii) Due to its mesoscopic character, the Boltzmann equation provides an immediate connection to the system's hydrodynamic variables, which naturally arise in the form of the various moments of the oneparticle distribution function. Therefore, the Boltzmann equation sets up a direct link between the microscopic dynamics of the system's constituent particles and the corresponding physics emerging on hydrodynamic length and time scales. Moreover, since the Boltzmann equation keeps track of all parameters used to formulate a particlelevel description of the system under study, the resulting hydrodynamics can be studied explicitly in terms of those microscopic parameters.
The Boltzmann equation and related kinetic frame-works have been employed in a number of previous studies on two-dimensional systems of self-propelled particles [19, 22, 26-28, 30, 31, 33, 38], mainly to derive a hydrodynamic description. To this end, an expansion technique is used to derive an infinite hierarchy of coupled equations of motion for the Fourier modes of the angular degrees of freedom of the one-particle distribution function. Since the zeroth, first, second, and higher order Fourier modes are directly related to the local average particle number, polar, nematic, and higher symmetry order parameters, the Fourier space description lends itself as a starting point to extract a hydrodynamic description. In general, this is achieved by truncating the hierarchy of Fourier space equations by means of a perturbation ansatz for the Fourier modes, which itself typically depends on system's symmetries. Following a similar strategy, a numerical approach based on an Enskog-like kinetic theory has recently been proposed in Ref. [38].
While such Fourier space approaches have been successfully applied to investigate the onset of collective motion, the validity of the underlying truncation schemes remains largely elusive for parameters deeper inside the ordered region. In this work, we propose a Boltzmann equation based approach which we will refer to as SNAKE (solving numerically active kinetic equations). SNAKE focuses on a direct numerical solution of the Boltzmann equation in real space, without any explicit recourse to the aforementioned coupled hierarchy of equations in Fourier space. Such a direct numerical approach offers a number of advantages in the study of active systems, which can be summarized as follows. First of all, the real space Boltzmann equation constitutes a closed description for the dynamics of the one-particle distribution function. Consequently, the range of validity of this direct numerical strategy is restricted only by inherent limitations of the Boltzmann equation itself but is not principally confined to a parameter region around the onset of collective motion. Secondly, SNAKE provides full access to the flexibility of the Boltzmann equation in the implementation of specific model systems which are based on a concrete microscopic picture of particle dynamics. This, thirdly, equips this approach with particularly powerful capabilities to incorporate physical boundary conditions, which can be formulated on the basis of actual particle-wall interactions rather than on the level of macroscopic, hydrodynamic field variables. Finally, since all pertinent hydrodynamic fields are easily computed from the one-particle distribution function by taking appropriate averages, this numerical approach allows to fully exploit the Boltzmann equation's potential to mediate between the physics on microscopic and macroscopic scales.
To demonstrate the benefits of SNAKE, we will present a detailed study of an archetypical model for active particles subject to binary, polar particle interactions, which has been proposed originally by Bertin et al. [22,26] and which we will henceforth refer to as "BDG model". As has been established previously [26], the BDG model transitions from a spatially homogeneous, fully isotropic state at low densities / high noise levels to a collectively moving state at high densities / low noise levels. In the vicinity of this phase transition, the polar order parameter field was shown to be susceptible to a longitudinal instability [26], which in turn has been linked to the emergence of stable solitary wave patterns propagating on an otherwise isotropic low density background [34,38]. Well beyond the transition line to collective motion, this longitudinal instability disappears from the Boltzmann equation [26] and a spatially homogeneous state of collective motion has been shown to emerge [34].
Our work is structured as follows. In section II, we discuss a discretization scheme for the Boltzmann equation of self-propelled particles in two spatial dimensions. The central result of this section is an update rule, Eq. (21), which lends itself to direct numerical implementation. This section also contains a discussion of the general strategy to implement virtually arbitrary boundary conditions and geometries. Readers not interested in the details of the actual discretization procedure underlying the SNAKE algorithm may skip sections II B -II D and proceed to section III.
There, we will exploit SNAKE 's capabilities to address the stationary states of the BDG model over a broad range in parameter space (section III A). In particular, we will show that our implementation of SNAKE consistently reproduces previously established analytical results on this model [26]. In section III B we focus on the onset of collective motion, which is accompanied by the formation of density segregated patterns and hystere-sis effects, as previously reported in numerical studies of the Vicsek model [8,14,38]. Additionally, in section III C, we extend previous results on the flocking transition in polar systems [26,28,38] well beyond the onset of order. Specifically, we demonstrate the presence of an additional first order phase transition from a spatially inhomogeneous, density segregated phase to a spatially homogeneous polar phase, taking place well inside the ordered, broken-symmetry parameter regime (section III D). Overall, we find striking similarities between the flocking transition to spatially homogeneous polar order, on the one hand, and the features of liquid-gas phase transitions in equilibrium systems, on the other hand. Analogous observations have been made previously in the context of a lattice model of active Ising spins, where a discrete up-down symmetry is spontaneously broken in two spatial dimensions [20]. In section III E, we shift our focus to a deeper discussion of density segregated patterns and report on the emergence of previously unseen "cluster-lane patterns" which seem to occupy the same parameter region as the familiar solitary wave patterns and to constitute a stable limit-cycle solution of the underlying Boltzmann equation. These patterns consist of parallel lanes of polar clusters moving in opposite directions; cf. Figs. 6(b,c) and Figs. 11(b,c) for snapshots. Intriguingly, the stability of such patterns provides a mechanism to establish states of nematic macroscopic order, arising from perfectly polar particle interactions on microscopic scales. Finally, in section III F, we investigate the coarsening dynamics for both, wave and cluster lane patterns within the density segregated phase.
Section IV, then summarizes our numerical approach, discusses our main findings, and gives an outlook of potential applications of SNAKE. In particular, this section gives a comprehensive and self-contained discussion of the most pertinent results established in the more technical sections II and III, and can be read independently of these technical parts. Section IV has been written to provide the reader with a direct and non-technical accessibility to the core results of this work.

II. DISCRETIZATION OF THE BOLTZMANN EQUATION
Here and in the following, we focus on a kinetic description for two-dimensional systems of actively propelled particles in an over-damped environment. As a simplifying assumption, we consider particles moving at constant speed v but variable orientation θ, such that each particle's velocity vector v can be written v = vê θ (ê θ : unit vector). For particles with nearly constant propelling forces, this is a reasonable approximation in the "Aristotlean regime" where the dynamics is governed by a balance between propelling and dissipative forces and inertial effects can be neglected. Moreover, for sufficiently dilute systems, binary particle interactions are believed to dominate over higher order interactions and the spatio-temporal dynamics of the one-particle distribution function f (x, θ, t) can be captured by means of the Boltzmann equation [39].

A. Boltzmann equation for active particles
The most general form of the Boltzmann equation reads The density of particles in the phase space element ω = [x, x + dx] × [θ, θ + dθ] is being convected due the particles' active propulsion (left hand side), and is, simultaneously, subject to sudden changes due to particle collisions, chemical reactions, etc. which is captured on the right hand side of Eq. (1a) in the form of source terms I source . To make contact to previous works on the Boltzmann equation for active systems [22,26,30], we consider two separate contributions to the source term, which take into account rotational diffusion due to a fluctuating background and / or noisy propelling mechanism, I sd , and binary particle collisions, I c , respectively. To be specific, rotational diffusion is implemented using a constant "tumbling" rate λ, and drawing the corresponding angular shifts from a given probability distribution p 0 . The self-diffusion integral I sd then takes the following form (omitting any explicit reference to space and time coordinates for brevity): where we introduced the short hand notation to take into account the 2π-periodicity of f (θ). In Eq. (1c), the two terms on the right hand side quantify the diffusion-related loss and gain of particles in the phase space element ω, respectively.
To asses the effect of binary particle collisions, we assume that the one-particle distribution function f (x, θ, t) varies on time and length scales which are large compared to the typical time between particle collisions and typical inter-particle separations, respectively. Therefore, on the scale of the Boltzmann equation (1a), the impact of collisions on the distribution of angles θ can be formulated locally in space. To capture the effect of binary collisions, we resort to a collision rule, i.e. a model-specific mapping θ w between pre-and post-collisional angles: Further, we assume that the collision process is subject to noise, which causes a rotation of the post-collisional angles θ (1/2) w by some random amount drawn from a probability distribution p. We then arrive at the following form of the collision integrals I c (again, omitting any explicit reference to space and time coordinates): where Γ(|φ 1 − φ 2 |) is a phase space factor ("Boltzmann Stoßzylinder"), quantifying the rate of collisions between particles with orientations φ 1 and φ 2 . The first term on the right hand side of Eq. (1d) gives the rate at which the phase space element ω loses particles due to binary collisions. Similarly, the second term quantifies the rate at which particles are scattered into ω. For the sake of greater generality, we further introduced the parameter ψ ∈ [0, π] in Eq. (1d), to account for systems with limited angular interaction ranges [34]. Eqs. (1) have been investigated in a number of previous studies [22,26,30,31,33,34] by mapping the mesoscopic Boltzmann equation (1a) to a set of hydrodynamic equations of motion for the first few moments of the oneparticle distribution function f (x, θ, t). As discussed in the introduction, this mapping is based on a truncation of terms in Fourier space representation of the Boltzmann equation which itself rests on a perturbation ansatz for the Fourier modes around an isotropic base state. To which extent such an approach yields reliable results in the context of ordered systems remains an open question; cf. Ref. [38]. In the following sections, we take a different route and focus on a kinetic description for active systems by a direct numerical solution of the Boltzmann equation (1).
Before we derive a discretized version of Eqs. (1), one last remark is in order: Since Eqs. (1b)-(1d) imply π −π dθ I source [f ; θ] = 0, Eq. (1a) gives rise to a conservation law for the local particle density ρ(x, t) = π −π dθ f (x, θ, t): where denotes the local momentum density. In what follows, we will restrict ourselves to the discussion of closed systems in which the particle flux across its boundaries vanishes [40]. Eq. (2) then entails conservation of the overall particle density where Ω denotes the system's two-dimensional spatial domain, and ||Ω|| its volume.

B. Dimensionless Boltzmann Equation
We start by rewriting Eqs. (1) in dimensionless form, in order to reduce the number of parameters. To this end, we measure time in units of the self-diffusion time scale λ −1 , distances in units of the "ballistic flight length" v λ −1 , and the one-particle distribution function f in units of the (constant) overall particle density ρ 0 [cf. Eq. (4)]. Finally, the phase space factor Γ quantifies the rate particle collisions and is, therefore, proportional to the particle speed v and the particle extension d. Dimensional analysis then dictatesΓ = d v. We thus arrive at the following convention of units which will be used throughout this work. The Boltzmann equation (1a) then takes the dimensionless form Here the expressions for I sd and I c are given in Eqs. (1c) (with λ = 1) and (1d), respectively, with f and Γ in their dimensionless forms. In the present system of units, the dimensionless density parameter ρ 0 measures the rate of collisions [ρ 0 d v] relative to the rate of self-diffusion [λ]. As an aside we note that in a regime of high densities and / or large ballistic flight lengths with ρ 0 1, the characteristic scales introduced above are no longer meaningful. In such cases, the characteristic time scale would be set by the typical time between subsequent collisions, (ρ 0 d v) −1 , the characteristic length scale by the typical inter-particle distance (ρ 0 d) −1 . In this limit, self diffusion can be neglected and the problem becomes independent of the overall particle density.

C. Discretization scheme
Having established a dimensionless form for the Boltzmann equation for active, two-dimensional systems, Eq. (6), this section is devoted to the derivation of a discretized version of Eq. (6). In order to keep the notation compact, we use "co-and contravariant indices" such that summation over repeated co-and contravariant index pairs is implied. Throughout this work, we reserve greek letters for spatial indices, and latin letters for angular indices. Moreover, we use the symbol || . . . || to denote the measure of spatial domains and angular intervals.
To discretize the spatial domain Ω, we consider a general tessellation by means of an arbitrary set of connected domains {s α } α=1,...,M such that they are a partition of Ω, i.e. [δ α β : Kronecker delta] s α ⊂ Ω, (7a) for all α, β ∈ {1 . . . M }. Similarly, to partition the angular range A = (−π, π], we consider tessellations by means of a set of simply connected intervals {a n } n=1,...,K , such that a n ⊂ A, (8a) for all n, m ∈ {1, . . . , K}; cf Fig. 1. The particular choices for the spatial tessellations are of no importance for the general development of our discretization scheme. However, some care must be taken in the construction of the angular tessellations. To avoid creating an artificial angular bias, the angles {θ n } n=1,...,K , used to represent the partition {a n } n=1,...,K , should sum up to zero. In addition, the angular partition should be constructed such that any special lattice found in the spatial tessellations (e.g. when space-filling, regular partitions are being used) should be reflected in any particular choice of angular tessellations. In what follows, we will use a n = θ n − π K , θ n + π K , θ n = 2πn K , and choose K ∈ 4N such that (i) the angular partition {a n } n ⊃ {0, π, ±π/2} and (ii) for each angular domain a i ∈ {a n } n there exists a domain a j ∈ {a n } n such that θ i = θ j + π (modulo 2π). This then implies We proceed by defining indicator functions for the spatial partition of Ω and the angular partition of (−π, π], F n (θ) ≡ 1, if θ ∈ a n , 0, else, The system's spatial domain Ω is discretized using an arbitrary tessellation {s α }α obeying Eqs. (7). A set of angular channels according to the tessellation (9) is attached to each spatial cell s γ . Spatial transport: Convection in the bulk of the system along the direction θn (indicated by the dashed arrow) amounts to translating each spatial cell s γ by the convection length τ . The convective transport operator T αk nγ , Eq. (16b), then redistributes the current value of f γ n among all neighboring cells according to their intersection with the convected cell.
which can be used to write Here, the time dependent density matrix f α k (t) ≥ 0 constitutes a discrete representation of the continuous oneparticle distribution function f (x, θ, t). To assess the temporal evolution of the density matrix f α k , it suffices to establish a transformation law of this matrix under small time increments. This transformation law then provides us with an "update rule" to compute the temporal evolution f α k (t) iteratively, starting from some initial distribution f α k (0), and hence to numerically solve the Boltzmann equation (6). In the following, we restrict ourselves to considering the time evolution in the bulk of the system, postponing a detailed discussion of boundary conditions to section II D.

Convective transport operator
Consider a small translation in the (rescaled) time coordinate, t → t+τ , where τ 1. To assess the convective contribution to the time transformation law for the density matrix f α k , we replace the convective derivative on the left hand side of the Boltzmann equation (6) by finite differences. Using a forward difference to substitute for the partial time derivative, we get Further, noting thatê θ · ∇f is the directional derivative of f alongê θ , we writê where now a backward difference is used to replace the spatial derivative in order to achieve cancellation of the terms ∝ f (x, θ, t). We chose equal step lengths in Eqs. (14a) and (14b) to account for the fact that the particles move at unit speed v = 1 according to our units convention. To derive the time transformation law for the density matrix f α k due to particle convection, we add Eqs. (14a) and (14b), substitute the discretization ansatz (13) for the one-particle distribution function f (x, θ, t), and project the resulting expression onto the discretized phase space volume ω α n = s α ⊗ a n . Introducing the projection operator we get where we introduced the convective transport operator which quantifies the time propagation of the density matrix f α n , due to the convective transport of particles. From Eq. (16b) we can extract a direct geometrical interpretation for the implementation of convective transport processes via T αk nγ ; cf. Fig. 1: The convective transport operator is computed by "convecting" the indicator function F γ for the spatial domain s γ along θ n , and integrating the resulting function over the domain s α . Geometrically, this amounts to computing the area of intersection between s α and a convected copy of s γ divided by the volume of the phase space element ω α n . The operator element T αk nγ thus measures the number particles being transported from ω γ n to ω α n during (t − dt, t). More precisely, T αk nγ gives the increase in particle density in the phase space element ω α n relative to the particle density in the phase space element ω γ n at time t − dt due to convective coupling between these elements in the time interval (t − dt, t).
A few remarks on the properties of the convective transport operator are in order. First of all, since, in the absence of self-diffusion and collisions, the particles maintain their current orientation, the convective transport operator is diagonal with respect to the angular indices for convective processes in the bulk of the system. This is, however, no longer true for convection in the vicinity of a (non-periodic) boundary, where particles generally reorient due to collisions with boundary walls. A corresponding redefinition of the convective transport operator will be given in section II D, where we discuss the implementation of non-periodic boundary conditions. Secondly, since T αk nγ reduces to the identity operator in the limit τ → 0, the convective transport operator becomes "quasi diagonal" in the spatial indices for sufficiently small time steps τ . More precisely, while the diagonal elements T αn nα are O(1), off-diagonal elements are O(τ ), i.e. convective coupling between neighboring elements becomes increasingly small as τ → 0. This, in turn, has important consequences for the computation of the angular self-diffusion and collision contributions to the time evolution of the density matrix f α n , which themselves are O(τ ) as will be shown below. Thus, up to linear order in τ , convective processes decouple from angular diffusion and collision, such that the latter processes can be computed independently inside each spatial domain s α .

Self diffusion operator
To discretize the self-diffusion integral I sd , we proceed along the same lines as in the preceding section and substitute the discrete representation of f (x, θ, t), Eq. (13), and project the resulting expression onto the phase space volume ω α n . Using the analytical form of the self-diffusion matrix, Eq. (1c), we find where the matrix elements of the self-diffusion operator quantify the net gain of particles (per unit time) in the angular element a n , by diffusive scattering from the angular element a k .

Collision operator
Finally, the discretized version of the collision integral I c defines the collision operator C kl n , which can be calculated using the standard procedure outlined in the previous two subsections. Starting loss term in Eq. (1d), we obtain: where the matrix elements measure the loss of particles (per unit time) in the angular element a n due to collisions with particles in the element a k . Similarly, considering the gain term of Eq. (1d), we find: where the matrix elements count the number of particles (per unit time) being scattered into the angular element a n as a result of binary collisions between particles with angles φ 1 ∈ a k and φ 2 ∈ a l . Combining Eqs. (18) and (19), we find with the collision operator C kl n given by:

Discrete time transformation law and summation rules
Having discretized the various terms in the Boltzmann equation, Eq. (6), we can now combine our findings from the previous three subsections, to arrive at the following transformation law, capturing the evolution of the density matrix f α n (t) under a small transformation in time, with the various operators, acting on the density matrix f α n (t), given in Eqs. (16b), (17b), and (20b). Eq. (21) defines a recursive relation, which allows us to numerically compute the solutions of the Boltzmann equation, Eq. (6), given initial conditions f α n (0). It, therefore, constitutes a discrete version of the Boltzmann equation for active systems, as introduced in section II A.
We conclude the current section by demonstrating that the particle conservation property, inherent in the original Boltzmann equation , Eq. (6), is being preserved by the transformation law, Eq. (21). To this end, we use the following operator summation rules, which themselves follow directly from the respective definitions of the corresponding operators: In particular, Eq. (22a) expresses the fact that convective streaming in the bulk conserves the volume of the phase space elements ω k γ = s γ ⊗ a k . Eqs. (22b) and (22c) capture the particle conserving properties of the self diffusion and collision operators, by asserting that the total gain of particles in each angular element is balanced by the total loss of particles in all remaining angular elements, and vice versa. Now, using Eqs. (22), we find i.e. conservation of the total number of particles. Finally, we note that, without further simplification, the computation time for an algorithm based on Eq. (21) is of order O(K 3 ) [41], and is thus effectively set by the number of angular slices K. In the absence of collision noise, on the other hand, the transition probabilities in the definition of the collision operator [cf. Eq. (19b)] reduce to delta peaks about θ (1/2) w , rendering the implementation of Eq. (21) of order O(K 2 ). Unless stated otherwise, we will use K = 32 in what follows. To test whether this is an appropriate choice for the angular resolution, we computed the spatially homogeneous solution of the density matrix f α n for a range of values K ∈ {32, . . . , 128} and a number of different density levels, using a 1 × 1 spatial grid. In fact, we observe no significant quantitative improvement upon increasing the number of angular slices in the presence of collision noise; cf. Fig. 2(a). The stationary, spatially homogeneous solution for the one-particle distribution function f (θ) as depicted in Fig. 2(a) has previously been observed numerically [21] and predicted analytically [23].

D. Boundary conditions
The analytical form of Eq. (21) applies equally well in the bulk and at the boundaries of the system. For non-periodic boundaries, however, the convective transport operator T αk nγ differs from its bulk form, Eq. (16b), and has to be modified. In this section, we start by discussing these modifications in abstract terms, and will later apply the generalized form of the convective transport operator to the specific case of reflective boundaries.
In the following, we shall confine ourselves to the discussion of fixed boundaries with a model specific reflection rule. Specifically, we assume particles to undergo a spontaneous reorientation upon collisions with the boundary, where θ and θ BC (θ) denote the particle's orientation before and after the boundary contact, and where θ BC can be either a deterministic or a random variable. For the sake of clarity, we will discuss deterministic particleboundary interactions. The treatment of stochastic boundary conditions follows along similar lines. To the considered order, self-diffusion processes and particle-particle collisions are evaluated at a fixed point in space and time in Eq. (21). To account for the presence of boundaries, we are, therefore, led to reconsider the derivation of the convective transport operator T αk nγ ; cf. section II C 1. For small but finite values of the time shift τ , we have to redefine the convective derivative, Eqs. (14), for all Here B τ denotes the set of all points in space which are within "convective reach" of the boundary δΩ of the system's domain Ω within the next time interval [t, t + τ ]. An appropriate redefinition of the convective transport operator can then be given using "inverse propagation functions": c −1 s (x, θ; τ ) ≡ particle position at time t, given that the particle is located at x and oriented alongê θ at time t + τ , (26a) c −1 a (x, θ; τ ) ≡ particle orientation at time t, given that the particle is located at x and oriented alongê θ at time t + τ .
Using these inverse propagation functions, we redefine the finite difference approximation of the convective derivative (14) as follows: which reduces to the bulk form, Eq. (14), upon substituting the inverse bulk propagation functions c −1 s (x, θ; τ ) = x − τê θ , and c −1 a (x, θ; τ ) = θ. Using this generalized finite difference representation to replace the convective derivative, we can proceed along the same lines of section II C 1 to arrive at the following, generalized form of the convective transport operator The generalized convective transport operator, Eq. (28) lends itself to implement virtually arbitrary boundary conditions, which enter the Eq. (28) by an appropriate definition of the inverse propagation functions (26). To illustrate the general concept, we conclude this section with a brief discussion of the implementation of reflective boundary conditions.
After the collision of a particle with a reflective boundary, the particle's orientation changes according to where φ b (x) denotes the local orientation of the system's boundary, and where θ and θ denote the particle's orientation before and after the collision with the reflective boundary, respectively. Moreover, if a particle's state is described by (x, θ) at time t and is known to have undergone a collision with the system's boundary in the time interval (t − τ, t], its spatial position x at time t − τ is given by where R φ b denotes a reflection with respect to the boundary's local tangent of orientation φ b . In summary, the inverse propagation functions in the case of reflective boundary conditions read: and else, (30b) where a, b denotes the line segment joining the points a and b.
As a final remark we mention that Eqs. (30) can be implemented straightforwardly using arbitrary boundary geometries. Although a detailed study on the impact of confining geometries on the behavior of active systems lies beyond the scope of our present work, SNAKE lends itself as a convenient starting point to address this important question. For sample simulations of active polar systems in circular and ring-like confinements, we refer the interested reader to the Supplemental Material [42].

III. POLAR SYSTEMS: PHASES AND TRANSITIONS
Having established a general framework to discretize the Boltzmann equation for active particles, Eq. (6), we now discuss actual applications of SNAKE in the context of active polar systems. We will focus our considerations on a simple model system with a polar collision rule and full angular interaction range. For the self-diffusion and collision noise probability functions p and p 0 , we choose zero-mean Gaussian probability distributions with variances σ 2 0 and σ 2 , respectively. This model has previously been proposed and studied by Bertin et al. [22,26] using mainly analytical methods. Our purpose here is twofold. First, we can resort to a number of well established analytical results [26] for this model, as well as a number of previous studies of closely related systems [8,14,28,38], to test for the accuracy and reliability of our computational scheme. Secondly, we can use SNAKE to go beyond the hydrodynamic picture of the onset of collective motion, as discussed in Refs. [22,26], and study the nature of the various phases, as well as the ensuing transitions between them.
In the following sections, we will use the system of units introduced section II B; cf. Eqs. (5). Moreover, we will use spatial tessellations in the form of regular twodimensional arrays of L×L rectangular grid sites with linear extensions |s α | = const ≡ x , and periodic boundary conditions. Some care must be taken for the appropriate choice of the size of one "Boltzmann element", x : The spatial patterns typically formed by these systems at the onset of polar order, are known to develop rather steep wave fronts, with density profiles varying on length scales comparable to the collision length (i.e. the typical length a particle travels between successive collisions) [14,38]. To be able to resolve such patterns with sufficient accuracy, we need x ∼ ρ −1 0 . In all subsequent sections, we will use densities ρ 0 = O(0.1 − 1), such that x = 5 is an appropriate choice. Unless stated otherwise, x = 5 and L = 100 will be used throughout this work.

A. Phase diagram of stationary states
The simple model system referred to above, Eq. (31), has been shown to exhibit a bifurcation, which separates a parameter regime at low density / high noise, from one at high density / low noise levels. In the high noise (low density) regime, an isotropic, spatially homogeneous phase (IHP) is observed. This homogeneous and isotropic state can be seen as a dynamical attractor for the system's dynamics at low densities, which we will refer to as "IHP fixed point" in the following. For low enough noise levels (and high enough densities), this IHP fixed point loses its stability and the system undergoes a transition toward polar order. The location of the corresponding bifurcation can be calculated directly from the Boltzmann equation, Eq. (6), by performing an expansion in terms of the moments f k = π −π dθ e ikθ f (θ) of the one-particle distribution function f (θ), and studying the resulting dynamical equation for the polar order parameter f 1 (cf. Ref. [22,26] for details): Here the various coefficients are functions of ρ 0 , σ, and σ 0 , and the dots indicate higher order and gradient terms. Although a complete description of the dynamics of f 1 , Eq. (32), contains couplings to an infinite number of higher order moments f n , the linear coefficient µ 1 (ρ 0 , σ, σ 0 ) can be calculated exactly within the Boltzmann equation approach and does not depend on any approximation scheme used to truncate these couplings [26,34]. The condition therefore provides us with an exact benchmark for the location of the bifurcation surface in parameter space. It is obtained by solving the above equation for ρ t [26]: For simplicity, we choose σ = σ 0 throughout this work. To test our implementation of SNAKE, we prepared different systems with varying values for the density (0.0125 ≤ ρ 0 ≤ 0.85) and noise levels (0.1 ≤ σ = σ 0 ≤ 0.95), starting from random initial conditions [43]. We performed the computations on a regular rectangular spatial grid with 100 × 100 grid sites and K = 32 angular slices, and evaluated the systems' (stationary) states after a total computation time of T = 15 000. The resulting phase diagram is shown in Fig. 2(c) along with the analytical result for the bifurcation line, Eq. (34). In accordance with previous analytical [26] and numerical results [14], we observe a phase transition between an isotropic, homogeneous phase (IHP) and a brokensymmetry density segregated phase (DSP) upon crossing the bifurcation line ρ t (σ), Eq. (34). A typical system snapshot inside the DSP is shown in Fig. 2(b). Deeper inside the symmetry-broken parameter domain, systems asymptotically reach a spatially homogeneous, polarized phase (SHP) after a transient episode of polar clustering and subsequent coarsening; cf. section III F. In the subsequent sections, we will give a more detailed analysis of the various observed phases and the corresponding transitions in between.

B. First-order transition from disorder to polar order
The transition between the isotropic homogeneous phase toward the formation of polar order (IHP → DSP) is probably one of the most studied phase transitions in active polar systems. Despite prior beliefs that this transition is of second order [5], whereby polar order was assumed to build up continuously, a growing number of more recent agent based [8,14,44] and kinetic studies [28,38] clearly indicate that the formation of density seg- The results indicate that the transition is first-order. (b) Hysteresis study. Average momentum against overall density ρ0 for σ = σ0 = 0.35 and two system sizes L ∈ {50, 100}. Solid symbols (path indicated by solid arrows) indicate the system's stationary state when starting from random initial conditions. Open symbols (path indicated by open arrows) indicate the system's stationary states upon lowering the density "quasi-statically" starting from a stationary wave-like pattern at ρ0 > ρt. The system exhibits hysteresis, corroborating the assertion that the underlying ordering transition is first order.
regated patterns at large enough system sizes render this phase transition first order.
To check for consistency of our numerical Boltzmann approach with these previous results, we studied the transition of the IHP toward polar for systems of varying linear sizes L ∈ {25, 50, 100}, subject to two different noise levels σ = σ 0 ∈ {0.35, 0.6}. For each choice of the system size and noise level, we prepared a set of systems at different overall densities ρ 0 ∈ [0.05 < ρ t (σ), 0.64 > ρ t (σ)] using random initial conditions. To quantify the phase transition, we computed the local polar order parameter ϕ α from the numerical solution of the density matrix f α n , From this, we then computed the spatially averaged polar order parameter, after the systems have settled into a stationary state (which, here and in the following, shall include moving patterns of "stationary" shape, as observed inside the DSP).
The results of this study are shown in Fig. 3(a). Our results indicate a jump discontinuity of the polar order parameter ϕ x at the transition to polar order, which is accompanied by the formation of spatial heterogeneities, typically in the form of traveling wave patterns (cf. Fig. 2(b); see also section III E). To within our chosen density resolution, ∆ρ = 0.01, the size of the observed jump-discontinuity is virtually independent of the system size for L 50, but depends on the noise levels in the system. For the smallest system sizes considered, L = 25 and smaller (not shown), the size of the jumpdiscontinuity shrinks down to zero, eventually rendering the observed phase transition continuous for system sizes below a critical size L * . These observations are in agreement with previous studies on this phase transition in the context of the Vicsek model [8,14] and Enskog-like kinetic theories [28,38]. For a more detailed discussion of the critical system size L * , we refer the reader to Refs. [28,38].
To scrutinize the subcritical nature of the bifurcation leading to the formation of polar order in large enough systems, we further checked for the existence of hysteresis effects in this model system. To this end, we numerically computed the "stationary" solution of the density matrix f α n inside the DSP, and then quasi-statically reduce the overall density ρ 0 in small steps. Here and in the following, the term "quasi-static" change of the overall density refers to the fact that the system is given a sufficient amount of time to equilibrate in between successive adjustments in ρ 0 . Two typical outcomes of such hysteresis experiments are shown in Fig. 3(b) for two different system sizes. There, closed symbols and arrows indicate the path a system follows in the ρ 0ϕ x plane, when starting from random initial conditions, exhibiting the familiar jump discontinuity at ρ 0 = ρ t . Open symbols and arrows indicate the inverse situation, where systems evolve from an initial state inside the ordered DSP parameter regime. In this latter case, our results indicate that the corresponding dynamical fixed point is stable down to values of ρ 0 = ρ i < ρ t well below the actual transition density ρ t , thus corroborating the subcritical nature of the corresponding dynamical bifurcation. Similar hysteresis effects have been reported in Refs. [14,31,38].

C. Density segregated polar phase
Since we have chosen to numerically solve the Boltzmann equation in real space, rather than in Fourier space, we expect SNAKE to be applicable even for parameters well inside the ordered phase, where Fourier space approaches would have to include an exceedingly large number of modes. Therefore, we continue our study of polar systems by presenting a more detailed analysis of the DSP region in parameter space. In anticipation of our discussion in section III E, we note that we frequently observe two distinct types of patterns within the DSP: traveling wave patterns and "cluster lane patterns". While the former pattern shape is by now well known to occur in polar active systems with metric interactions [14,26,34,38], to our knowledge the emergence of cluster lane patterns has not yet been reported in the literature so far. We devoted the separate section III E to these novel patterns, and restrict our present discussions to the well-known case of traveling waves. After initial transients, during which "polar droplets" form and coalesce out of random initial conditions, the resulting "stationary" traveling wave patterns (when viewed in a co-moving frame) become translationally invariant in the lateral direction x ⊥ (perpendicular to the wave's velocity vector). The structure of these wave patterns can, therefore, be described by means of the longitudinal coordinate x (along the direction of the wave's velocity vector). Figure 4(a) shows typical density profiles of such polar waves as a function of the distance from the transition line at ρ t (σ). The observed wave profiles exhibit a characteristic front-rear asymmetry, similar to those reported in previous numerical studies [14], whereby the asymmetry becomes more pronounced as the parameters are chosen deeper inside the ordered region.
The waves themselves can be viewed as high-density structures, traveling on an isotropic, low-density background [cf. appendix A and Fig. 10]. To quantify the density segregation between the wave fronts and the isotropic background as function of the distance from the system's phase boundary, we introduce the following characteristic density scales: as measures for the density levels within the highdensity, moving phase ("HD phase"), and the lowdensity, isotropic background ("LD phase"), respectively. Figure 4(b) shows the characteristic density scales, and the density separation parameter, as a function of the system's overall density ρ 0 at fixed noise levels. The data shown in Fig. 4 has been com-puted using the following simulation protocol: For a system well inside the DSP parameter regime, we followed its time evolution starting from a random initial condition until it reached a stationary state with a traveling wave pattern. The density density matrix f α n corresponding to this stationary wave pattern was then used as initial condition in subsequent numerical computations, during which the density level was "quasi-statically" decreased or increased until the density segregation parameter dropped back to zero, i.e. until a spatially homogeneous state is reached.
This protocol allows us to explore the parameter range over which the "DSP fixed point" of the system's dynamics remains locally attracting. Before we embark on a more detailed discussion on the stability of the various dynamical fixed points and observed phase coexistence in parameter space in section III D, we take a little detour and briefly comment on the observed wave profiles as the overall density is varied. First, we observe that the density level of the LD phase, ρ min , is virtually independent of the overall density ρ 0 across the entire region of existence of the DSP. Any "surplus density" arising from an increase in ρ 0 must therefore be accommodated within the HD phase, and, consequently, implies a change in the wave profiles. At the lowest densities (ρ 0 < ρ t ) where wave patterns are still stable, we observe wave profiles of relatively weakly pronounced front-rear asymmetry [cf. Fig. 4(a), ρ 0 = 0.2]. As the overall density is increased slightly above ρ t , surplus density is being accommodated primarily by increasing the characteristic density ρ max , leaving the front-rear asymmetry virtually unchanged [cf. Fig. 4(1), ρ 0 = 0.225]. For still larger values of the overall ρ 0 , the characteristic density ρ max quickly saturates and surplus density is being accommodated by increasing the width of the solitary waves. This, in turn, leads to a drastic increase in front-rear asymmetry.
The latter regime of wave broadening can be understood intuitively by employing a simple mean field picture which assumes constant density levels within the HD and LD phases, respectively. To this end, consider a system of arbitrary volume ||Ω|| in the DSP wave broadening parameter regime. Since both characteristic densities, ρ min and ρ max , are constant by assumption, conservation of particle number N V hd ρ max + V ld ρ min (V hd/ld : volumes occupied by HD and LD phases, respectively) dictates: Figure 4(c), where the surface fraction V hd /||Ω|| is recorded as a function of ρ 0 confirms this linear dependence, and thus corroborates the above mean field picture. Moreover, since V hd ≤ ||Ω||, Eq. (38) implies a crossover scale ρ h = O(ρ max ), above which the LD phase gets depleted and the system becomes spatially homogeneous. We checked that the mean field picture conveyed by Eq. (38) is actually independent of the system size L, and becomes quantitatively robust against changes in the system size for L 50; cf. Fig. 4(c). Interestingly, the composition of the system in terms of HD and LD phases as a function of of the overall density ρ 0 exhibits qualitatively analogous behavior in the two-dimensional lattice model of active Ising spins of Ref. [20], where a Z 2 (up-down) symmetry is spontaneously broken. This observation strongly suggests that the qualitative features of the flocking transition toward spatially homogeneous polar order in active systems might actually be universal across different symmetry classes.

D. Transition to spatially homogeneous polar state and pattern selection
In the previous section, we investigated the evolution of systems inside the DSP regime, as the overall density ρ 0 is varied. We observed that the DSP parameter regime, which can be defined by non-zero values of the density separation parameter η, is hallmarked by a linear variation of the surface fraction of the HD phase, V hd /||Ω||, with ρ 0 . Importantly, the DSP regime is delimited by sharp jumps in η and V hd /||Ω|| at both of its boundaries, i.e. at the transition toward spatially homogeneous isotropic state (low densities) and at the transition toward a spatially homogeneous polar ordered state (high densities). As has been discussed before (e.g. Refs. [14,38] and section III B), such non-regular behavior at the low density limit of the DSP regime can be attributed to the first-order nature of the ordering transition between the IHP and DSP parameter regimes. However, for the transition between the DSP regime and the SHP regime comparable results are not available, yet. The present section is devoted to fill this gap, and to illuminate the nature of the phase transition between DSP and SHP regimes. Similar to the case of the IHP ↔ DSP phase transition, we will observe hysteresis effects which are highly suggestive of a first-order transition between the DSP and SHP regimes.
To assess the nature of the DSP ↔ SHP transition, we chose density and noise parameters inside the SHP regime (ρ 0 = 0.415, σ = σ 0 = 0.5), and let the system evolve until it has reached a state of spatially homogeneous polar order. The corresponding stationary solution for the density matrix f α n was then used as initial condition in subsequent computations. In contrast to the computation protocol in section III B (where parameters were adjusted quasi-statically), the system's overall density ρ 0 was then quenched to successively lower values ρ 0 ∈ {0.41, 0.405, . . . , 0.25}, and the resulting pattern at asymptotically long times was recorded; for details of the numerical protocol cf. Appendix C.
In Fig. 5(a) the corresponding results of this study are shown in terms of the surface fraction V hd /||Ω|| occupied by the HD phase. There, open symbols and arrows indicate data points corresponding to initial conditions inside the SHP regime. Closed data points refer to initial conditions inside the DSP regime (where the overall density was then quasi-statically increased). The data show that the transition between the DSP and SHP regimes is accompanied by hysteresis effects, supporting our initial claim that the corresponding transition is first-order. Phase space coexistence between both regimes can be observed for densities ρ dh (σ = σ 0 = 0.5) ≈ 0.345 ≤ ρ 0 ≤ ρ h (σ = σ 0 = 0.5) ≈ 0.39 [45].
and find excellent agreement. Summarizing the results from section III, we can now extend the picture conveyed by the phase space diagram in Fig. 2(c), which has been recorded using random initial conditions at different points in parameter space. Our discussions in this and previous sections disclosed regions in parameter space, where phase selection by the system is not uniquely related to the specific choice of parameters (i.e. density and noise levels). Instead, phase selection is sensitive to initial conditions in these parameter regions, and IHP / DSP coexistence at the low-density end, and DSP / SHP coexistence at the high-density end of the DSP regime is observed. We condensed these findings by sketching the corresponding bifurcation diagram in the ρ 0 -η plane in Fig. 5(b). There, solid lines indicate the position of stable fixed points of the system's nonlinear dynamics, as measured by actual computations. Dashed lines are drawn "by hand" and indicate an estimate of the position of the corresponding unstable fixed points. Figure 5(b) corresponds to a horizontal cut through the phase diagram shown in Fig. 2(c) at σ = σ 0 = 0.5. Comparing both figures, we see that the phase diagram in Fig. 2(c) displays the DSP only in the relatively narrow window ρ t ρ 0 ρ dh since initial conditions very close to the η = 0 fixed point are chosen. Due to the subcritical character of both transitions, IHP ↔ DSP and DSP ↔ SHP, density segregated patterns can actually be observed over the much broader parameter range ρ i < ρ 0 < ρ h . Interestingly, an analogous extension of the existence region of density segregated patterns has recently been reported in the context of nematic ordered active systems [31,35].

E. Phase separated patterns
In previous works on polar systems with identical [26], or similar model properties [14,38], the emergence of density segregated states has been discussed in terms of traveling wave patterns with "infinite" extent in the lateral direction [46]; cf. snapshot in Fig. 2(b). However, in our numerical studies we find an additional, genuinely distinct type of density segregated patterns, which is fundamentally different from the familiar solitary wave patterns. These patterns consist of a number of regularly shaped, moving high density regions of finite lateral extent, where each such region carries a large net polarity. To emphasize their finite spatial extent, we will refer to these high-density polar regions as (polar) "clusters". Intriguingly, the polarities of these clusters are oriented along a common broken-symmetry axis with orientations distributed such that the overall symmetry of the system is nematic. More precisely, clusters are arranged in lanes, where each lane hosts a train of equally spaced clusters of parallel polarities. In our numerical studies, we observe two distinct types of cluster lane patterns: First, we find antisymmetric cluster lane solutions, where all cluster lanes have identical geometrical properties (same lane width and number and shape of polar clusters), except for a reflection of the cluster polarity from one lane to the next, thus rendering the overall symmetry of the system purely nematic. Typical snapshots of this first type of antisymmetric cluster lane patterns are shown in Figs. 6(b,c) and in Fig. 11(b). Second, we observe asymmetric cluster lane patterns where neither the number of clusters per lane, nor the morphological properties of clusters in different lanes are symmetric; cf. Fig. 11(c). We note that the lane widths of this latter type of asymmetrical clusters were subject to a very slow but continuous drift during the computation times considered in this work. We can, therefore, not exclude the possibility that this asymmetrical pattern would actually give rise to the more familiar wave pattern in the limit of long times (in which case the asymmetrical pattern would constitute an intermediate asymptote of the system's dynamics). In contrast, we find that the structural properties of the antisymmetric cluster lane patterns do converge, which strongly suggest that these patterns constitute a genuine limit-cycle solution of the underlying Boltzmann equation, Eq. (6). In the following we will, therefore, restrict ourselves to examine antisymmetric cluster lane patterns which, moreover, are more frequently observed than their asymmetric counterparts. In particular, we give a detailed discussion highlighting the basic mechanisms stabilizing these remarkable patterns.
Before we embark on a more thorough discussion on the basic mechanisms stabilizing these nematic ordered patterns, we briefly comment on the observation statis- tics of antisymmetric cluster lane patterns in our numerical solutions. As discussed in detail in Appendix B, the broken symmetry axis in such cluster lane patterns is effectively constrained to one of the two coordinate directions, due to the finite size of the computational grids used in our numerical studies. To obtain a counting statistics, we chose different sets of model parameters inside the DSP regime and used 350 different seeds per set of parameters to prepare independent systems with random initial conditions. We then let each of these systems evolve until a structurally stable, density segregated state had developed. From all 350 patterns, we selected those with a mean polar / nematic order along one of the coordinate axes and computed the observation frequency p of antisymmetric cluster lane patterns among these "parallel patterns" [47]. While the observation frequency of antisymmetric cluster lane patterns is rather low for the smallest systems sizes considered in this work (L = 50; p ≈ 5%), these patterns make up a considerable fraction of all parallel patterns for system sizes L 100, in which case we find p ≈ 25% (we find p ≈ 40% if we count all cluster lane patterns, including asymmetric cluster lanes).
We shall now give a more precise definition of the term "cluster lane" and elucidate the dynamical mechanisms underlying their formation and stabilization. To this end, consider the antisymmetric cluster lane pattern of Figs. 6(b,c) with axis of nematic order along the x-direction. To assess the particle exchange dynamics between the various clusters, we recorded the transverse particle currents g y (x, t), Eq. (3), and computed its lon-gitudinal and time average according to where the time scales t 0 = ∆T = 10 000 λ −1 are chosen such that the computation of the time averaged particle currents is done on the basis of a structurally stable cluster lane pattern, i.e. in the limit of long times. Fig. 6(a) shows j y (y) along with the (instantaneous) particle currents j y (y, t) corresponding to two particular time instances, indicated by the snapshots in Fig. 6(b,c). From inspection of the temporal average of the particle current profile j y (y) [red curve in Fig. 6(a)], we extract the following structural features: The overall functional form of j y (y) resembles a sawtooth shape with large positive slopes (i.e. divergence of j y ) in narrow depletion zones close to the cluster lane boundaries, and a shallow negative slope (i.e. convergence of j y ) throughout the spatially extended cluster zones. The borders y = y i b (with i = 1, 2, . . . labeling the different border lines) delineating adjacent cluster lanes are located at the zeros of the particle current inside the depletion zones, for which j (s) y (y, t) y=y i b = 0 for all t and i (41) holds true in the stationary state (indicated by the superscript s). The sawtooth shape of the time averaged current j y therefore confers the following picture illustrating the balance of particle exchanges across different cluster lanes: The bulk of the clusters, i.e. the cluster zones inside each lane, receive particles at the expense of a narrow depletion zone at the cluster lane borders. Within these depletion zones, large transverse, inwardly directed particle currents, "feeding" the bulk of the clusters, arise due to periodically recurring "cluster grazing"; cf. solid black curve in Fig. 6(a) and snapshot in Fig. 6(b). As time progresses, these particle currents propagate toward the center of the clusters. Concomitantly, rotational diffusion attenuates these propagative modes and gives rise to net diffusive currents at the tails of these modes which entail diffusive spreading of each cluster across the lane boundaries; cf. dashed black curve in Fig. 6(a) and snapshot in Fig. 6(c). In steady state, these diffusive currents balance such that there is no net current across the lane boundaries. When the inwardly directed (attenuated) transverse currents eventually arrive at the cluster centers (which, again, correspond to zeros of the transverse particle currents for all t), the aligning collision rule converts them into longitudinal currents, thus establishing a feedback mechanism to the macroscopic order parameter.
To demonstrate the structural stability of these antisymmetric cluster lane patterns, we further computed the net particle exchanges [δC: boundary delineating the right moving cluster in Figs. 6(b,c)] across the cluster lane boundaries at y 2 b , y 1 b (y 2 b > y 1 b ). Since in the stationary state j (s) y (y i b , t) = 0, we expect the magnitude of j y (y i b , t) to converge to zero in the limit t → ∞. In Fig. 7, the net ratesṅ at which the left and right moving cluster lanes in Figs. 6(b,c) gain and loose particles are plotted as a function of time in the long time regime λ t ∈ [10 000, 20 000], where the cluster lane pattern is virtually stationary. Due to overall particle conservation, the red and black curves in Fig. 7 sum up to zero. The net rate at which both cluster lanes gain and lose particles are extremely small in the displayed long time limit and decay exponentially as t → ∞. This indicates that a balance of particle numbers is being established between both lanes, such that the corresponding cross-lane particle currents cancel in the limit of asymptotically long times. The closeup view in the inset of Fig. 7 reveals that this balancing process is driven by cluster grazing events (spikes) during which the right moving lane gives off particles to the left moving lane, and transverse diffusion (bumps) which leads to a net current of particles into the right moving cluster lane. We stress, however, that even in the stationary state large amounts of particles are transported laterally across the cluster lane boundaries, and only the net exchange of particles vanishes due to a global cancellation of periodically oscillating particle currents which are constantly maintained by cluster grazing and transverse diffusion.
In summary, cluster lane patterns constitute an intriguing, stable limit-cycle solution of the Boltzmann equation. The transport processes entailed by these cluster lane patterns differ from those of the more familiar solitary wave patterns in two important respects: First of all, in contrast to the solitary wave solution, where a macroscopic transport of particles is confined to the longitudinal direction, the antisymmetric cluster lane pattern approaches a stationary state in which macroscopic, periodically oscillating particle currents give rise to an interesting local net transport of particles in lateral directions, which is balanced only on global scales. Second, the global order of this limit-cycle solution is nematic rather than polar. The cluster lane solution of the Boltzmann equation thus highlights a physical mechanism by which a system of self-propelled polar particles and purely polar interactions is capable of eliciting a macroscopic state of global nematic order. Put differently, clusters could be viewed as polar quasi-particles with only nematic interactions, thus rendering the system nematic on macroscopic scales. For a sample movie of a cluster lane solutions refer to the Supplemental Material [42].
Before we conclude, some additional remarks are in order. Our numerical results, presented in this section, provide strong evidence for the actual stability of these novel cluster lane patterns. Moreover, we checked that the emergence of cluster lane patterns is not a singular property of the specific features of the BDG model: As a matter of fact, we observed the formation of cluster lane patterns over a large range of parameters for a class of more general collision rules (cf. Ref. [34]) deviating from the half-angle alignment rule, Eq. (31), employed in the BDG model (data not shown). This suggests that the cluster lane solution is of an even more gen- eral scope within the framework of the Boltzmann equation approach to active polar systems. Nevertheless, it is important to keep in mind that the Boltzmann equation itself rests on a number of specific assumptions, including the absence of multi-particle interactions beyond binary collisions and the neglect of correlation effects among particle states as a consequence of the molecular chaos assumption. To which extent the emergence of cluster lane patterns is hinged on such hidden assumption constitutes an important question to be addressed in order to clarify on the role of such cluster lane patterns for actual active matter systems. We, therefore, hope that future studies employing alternative modeling approaches will shed light on this important question.

F. Coarsening dynamics
Having accessed the frequency of the stationary patterns of phase-segregated clusters and waves, in this section we address the dynamics of the coarsening process for both patterns. For classical phase-separation of fluid mixtures in the droplet-regime (asymmetric quench, [48]), coarsening of structures is driven by Ostwald-ripening [49] and Brownian-coalescence [50]. Since both mechanisms rely on diffusion of either droplet material or droplets, the corresponding growth laws exhibit identical scaling as a function of time, i.e. the characteristic size of structures (clusters or droplets) (t) ∝ t 1/3 [48]. Interestingly, the same coefficient is reported in a model of self-propelled particles exhibiting run-andtumble dynamics on a lattice [51] and a similar (but slightly smaller) coarsening exponent has been found in actively propelled hard spheres [52,53] coinciding with theoretical predictions [54]. Here, we examine the ques-tion whether there is dynamical scaling for the phasesegregated, polar-moving cluster and wave-patterns.
To this end, the characteristic size of a structure (t) must be defined: We scan over the lattice along x while fixing y, and vice versa, and determine the corresponding probability distribution p x/y (l) of connected areas l for which ρ/ρ 0 > 1, i.e. the density is larger than the homogenous density ρ 0 . The characteristic length for scanning over x or y, respectively, is then defined as the first moment of the probability distribution, x/y (t) = dl p x/y (l; t) l.
However, in contrast to a classical phase-separation in fluid mixtures [48] or actively propelled hard spheres [52,53], the emerging phase-segregated patterns discussed in this manuscript are not isotropic. Consequently, the characteristic length has to be determined with respect to the mean polarity, i.e. transversal and longitudinal to it. In the following we analyze the characteristic length transversal to the stationary polarity ⊥ (t), whose long time asymptote gives the lateral correlation length of the pattern, i.e. ⊥ (t → ∞) = ξ ϕ . Speaking in terms of Eq. (43), the transversal characteristic length ⊥ (t) is equal to x/y (t) when the respective pattern moves along the y/x-axis. The corresponding results for wavepatterns and cluster patterns are shown in Fig. 8. Analyzing several random initial conditions, the characteristic transversal length follows roughly ⊥ (t) ∝ t 1 for both, wave and cluster lane pattern morphologies. This linear scaling is compatible with previous results from agentbased simulations [14]. While ⊥ (t) coarsens up to the system size L in case of wave patterns [ Fig. 8(a)], it saturates significantly below L for cluster patterns [ Fig. 8(b)]. Moreover, the pronounced asymmetric fluctuations of the characteristic length in time for the cluster patterns correspond to configurations where anti-parallel moving clusters are positioned directly transversal to each other.
Here, we would like to derive a simplified picture to explain the observed linear scaling for the domain size ⊥ (t) ∝ t 1 . In the following we are led by the coarsening dynamics in classical binary mixtures [48]. Say all polar high-density domains move ballistically with a velocity ∼ 1 (measured in units of v 0 ) and are of characteristic size ∼ ⊥ . Moreover, we assume that the total surface coverage of high-density, polar domains is approximately conserved during the coarsening process, and, especially, during each coalescence event of two polar domains. This implies that the characteristic time-scale between two coalescence events scales according to τ ∼ ⊥ . On average, coalescence causes the radius polar droplets to increase by a factor of √ 2 per coalescence event in 2D, where we assume colliding domains of approximately equal size. We then infer the following scaling picture for relative change of the polar domain size implying that ⊥ (t) ∼ t 1 .
Our analysis for the coarsening process within the phase-segregated regime of the phase diagram [see Fig. 2(c)] is of course limited by the system sizes investigated (here: L = 150). We hope that future studies involving larger system sizes confirm our observation of dynamical scaling in the coarsening dynamics of polarmoving phase-segregated morphologies.

IV. CONCLUSION
In this work, we have developed a general numerical framework to solve the Boltzmann equation for two dimensional systems of self-propelled particles, which we referred to as SNAKE (solving numerically active kinetic equations). Since we have chosen to solve the Boltzmann equation in real space, our algorithm is applicable even far beyond the onset of macroscopic order, so long as the basic assumptions underlying the Boltzmann equation itself (e.g. molecular chaos) are fulfilled. Specifically, our numerical framework lends itself to the implementation of arbitrary boundary geometries, and we discussed how boundary conditions can be formulated directly in terms of microscopic particle-wall interactions.
We applied our algorithm to study an archetypal twodimensional model system of actively propelled particles with binary, polar interactions, as previously proposed by Bertin et al. [22,26], and which we refer to as "BDG model". The BDG model is a kinetic model which is inspired by Vicsek's agent based model [5]. It considers a collection of spherical, self-propelled particles moving at a constant speed v such that each particle's state can be captured by a two-dimensional position vector and an angular variable characterizing the current particle orientation, i.e. direction of its velocity vector. In the absence of noise, particle orientations change in response to binary particle collisions which cause the velocity vectors of both collision partners to align along the average orientation of both particles prior to their collision. In addition, noise effects are assumed to influence both the free motion of particles between subsequent collisions ("selfdiffusion"), as well as the alignment process of the binary collisions themselves.
We showed that our algorithm is fully consistent with previous analytical results on the BDG model itself [26], as well as with the most pertinent results from extensive simulation studies on the Vicsek model [14] at the onset of collective motion. In particular, our numerical results corroborate the first order nature of the phase transition toward collective motion and reproduces hysteresis effects and, therefore, phase coexistence between a spatially homogeneous, isotropic state and the emergence of polar ordered, density segregated patterns; cf. [14,38]. In the framework of our numerical investigations of the BDG model, we have established two novel key insights.
First, we have studied the emergence of macroscopic collective motion, including parameter regions well be-yond the onset of macroscopic order. Our results convey a comprehensive picture, highlighting the core features of the transition between a fully isotropic, disordered state at low densities / high noise values, and a long-range ordered, spatially homogeneous state at high densities / low noise values. We found that this ordering transition bears strong structural similarities to liquid-gas transitions in equilibrium systems. There, if the density of a gas is isothermally increased beyond the binodal, the gas phase becomes metastable and, therefore, susceptible to large magnitude density fluctuations. If present, such large fluctuations lead to nucleation and growth of high density liquid droplets and eventually lead to phase separation into spatially coexisting liquid and gas phases. At still larger densities, the gas phase becomes unstable and spinodal decomposition occurs. In this case, arbitrarily small density fluctuations are amplified and the equilibrium system immediately separates into spatially coexisting low-density gas and high-density liquid phases. The mole fractions of these coexisting phases depend linearly on the system's specific volume such that the liquid phase completely takes over for large enough densities. The same line of arguments applies to the opposite case, where the density of a liquid system is being decreased. Intriguingly, we find the exact same overall behavior in the context of the ordering transition in the BDG model; cf. Fig. 5(b). At constant noise values ("temperature"), we identified binodal (ρ i ) and spinodal (ρ t ) density thresholds, beyond which the isotropic "gas" phase becomes metastable (susceptible to large enough density fluctuations) and unstable, respectively. Beyond the spinodal density threshold, ρ t , arbitrarily small density fluctuations give rise to the formation of polar droplets, whose size grows linearly over time until a macroscopic, density segregated state has developed. Within this regime, well ordered high density waves or clusters ("liquid phase") are moving on a lowdensity, isotropic background ("gas phase"). Similar to the equilibrium case, the surface fractions of both phases linearly depend on the system's overall density and the "gas phase" gets extinct at high enough overall densities. Conversely, any spatially homogeneous polar ordered system at high densities becomes metastable when the overall density falls below the second binodal density threshold, ρ h , and spinodal decomposition occurs upon further reduction of the overall density below the second spinodal density, ρ dh . An analogous discussion, highlighting the similarity between flocking and liquid-gas transitions, can be found in a recent work on a lattice model of active Ising spins [20]. Despite the fact that both, the BDG model and the active Ising spin model of Ref. [20] belong to different symmetry classes, the qualitative features of the flocking transition to spatially homogeneous polar order remain virtually unchanged. This, in turn, strongly suggests a universal character of the flocking transition to homogeneous polar order across different modeling classes and underlying symmetries.
Second, our results indicate the emergence of in-triguing spatial structures after spinodal decomposition, which are genuinely distinct from wave-like patterns of infinite lateral extent. Inside the parameter regime where density segregated patterns are found, we observed a novel type of patterns, we referred to as "cluster lane patterns". These patterns consist of parallel lanes of polar clusters, traveling in opposite directions. We gave a detailed discussion of the physical mechanisms underlying the stability of this limit-cycle solution of the Boltzmann equation in terms of the transverse cross-lane transport of particles. Most importantly, these cluster lane patterns differ in two essential respects from the familiar solitary wave patterns, with which they seem to coexist in the DSP regime of the parameter space: First of all, while the solitary wave solution of the Boltzmann equation reflects the polar symmetry of the constituent particles and their mutual interactions, the cluster lane solution of the Boltzmann equation renders the overall symmetry of the ordered state nematic. Secondly, while the solitary wave solution gives rise to purely diffusive particle currents in the lateral direction which "die out" in the limit of long times, the cluster lane solution is hallmarked by periodically recurring propagative lateral particle transport modes which persist even for asymptotically long times. Regarding their potentially spectacular consequences, we hope that future research will further clarify on the role of these intriguing patterns.
To conclude, the Boltzmann equation provides a convenient platform to implement particle-level descriptions of active model systems, and to assess the ensuing spatiotemporal dynamics on macroscopic scales. While this latter step is prohibitively difficult to achieve in analytical treatments of the Boltzmann equation, our numerical framework allows to study the macroscopic properties of active systems even far beyond the onset of collective order. SNAKE, therefore, provides the basis to study a countless variety of binary collision models with arbitrary interaction symmetries, which can be achieved by appropriately (re)defining the collision and diffusion operators; cf. sections II C 2 and II C 3. More general modeling classes can be constructed by redefining the convection operator, section II C 1, which, for instance, allows for straightforward realization of density dependent motilities [19]. In particular, the kinetic framework discussed in this work is easily extended to a multi species description of active matter [33]. Therefore, SNAKE grants direct access to a wealth of new and exciting questions dealing with the spatiotemporal dynamics of actively propelled systems in conjunction with chemical reactions or in the context of mutually competing species both, in free space and in arbitrary confinement geometries.

ACKNOWLEDGMENTS
We thank Hugues Chaté and Eric Bertin for stimulating discussions. This work was supported by the Deutsche Forschungsgemeinschaft in the framework of the SFB 863 "Forces in Biomolecular Systems" (Project No. B2), and the German Excellence Initiatives via the program "NanoSystems Initiative Munich (NIM)".

Appendix A: Wave patterns revisited
In this appendix, we briefly comment on momentum and velocity fields for wave-like patterns inside the DSP regime, as well as on their bending instability at the DSP → IHP transition.
In section III C we investigated the range of stability of wave-like patterns, as the overall density of the system is gradually decreased or increased. In both cases, our results suggested a first order transition toward a spatially homogeneous isotropic or polar state, respectively. Thereby, the former case, where the system undergoes a phase transition toward a spatially homogeneous, isotropic state upon decreasing the density below the threshold value ρ i , is accompanied by a "reinforcement" of the wave pattern just before the DSP → IHP phase transition occurs; cf. highlighted data points in Fig. 4(b) and corresponding data points in Fig. 10(b). This reinforcement effect can be explained by the finite system size and is illustrated by the sequence of snapshots displayed in Fig. 9. There, the initial state of the system corresponds to a wave pattern of relatively large lateral extent, which is wrapped around the torus (representing the system's domain Ω) twice [ Fig. 9, left panel]. As the overall density is decreased toward the transition threshold, ρ 0 ρ i , this laterally elongated pattern eventually becomes unstable toward a bending instability [ Fig. 9, middle panel]. However, on a periodic domain, this instability does not necessarily destroy the wave pattern, but rather can lead to a reorientation of the wave so as to contract its lateral extension [ Fig. 9, right panel]. This lateral contraction, in turn, leads to a reinforcement of the wave pattern, with both, particle density and momentum across the contracted pattern being increased as compared to the initially elongated wave pattern; cf. Figs. 4(b) and 10(b). Note, however, that this reinforcement mechanism is an artifact of the finite system size and has no analog in infinite systems.
We conclude this appendix by briefly reviewing the characteristics of the momentum and velocity fields for wave patterns inside the DSP regime. To this end, we have recorded the waves' momentum profiles at various overall density ρ 0 inside the DSP regime, Fig. 10(a), and the characteristic particle velocities inside the LD phase (v min ) and HD phase (v max ), together with the system's spatially averaged velocity, v , as a function of ρ 0 , Fig. 10(b). The waves' momentum profiles closely resemble their respective density profiles both, with respect to their asymmetric shape and their growing asymmetry with increasing density. As for the characteristic velocities inside the DSP regime, we observe an approximately linear increase of the spatially averaged velocity v after an apparently discontinuous jump at the DSP → IHP transition. The characteristic velocity inside the HD phase remains virtually constant over the whole span of the DSP regime. In contrast, the characteristic velocity inside the LD phase is vanishingly small only up to densities ρ 0 0.3 ≡ ρ v , in which case the wave patterns can, in fact, be characterized as high density bands moving on an isotropic low-density sea of particles. For higher densities inside the DSP regime (i.e. ρ 0 < ρ h ), wave patterns persist but with increasing longitudinal extensions such that the waves' tails start to infiltrate the LD phase. As a result of this infiltration, the characteristic velocity inside the LD phase is being increased to a density dependent, finite value, although the corresponding characteristic density of the LD phase remains virtually unchanged at a subcritical value, ρ min < ρ t ; cf. Figs. 4(a,b). One possible explanation for this observation of finite velocities inside the LD phase could be as follows: Since ρ min < ρ t , the dynamics of the locally averaged momentum inside the LD phase is such that |g min | (and thus v min ) decays exponentially. Above a certain threshold width of the density waves (i.e. for ρ 0 > ρ v ), however, the time scale between subsequent arrivals of wave maxima becomes comparable to the characteristic decay time scale for local momenta / velocities at ρ = ρ min . Hence, in the stationary limit, the LD phase attains nonvanishing average particle velocities, despite the fact that local conditions in the LD phase are such that polar order tends to be destroyed.

Appendix B: Cluster lane orientations
Here, we will give a brief account of possible cluster lane orientations in systems with periodic boundary conditions. In particular, we will illustrate that the observation rates for cluster lane patterns do strongly depend on cluster lane orientations, and, therefore, explain why only a small number of typical cluster lane orientations is expected to be observed in simulations with finite system sizes.
For the sake of greater clarity, we will consider systems with a 1 : 1 aspect ratio. Further, we will use the onedimensional "midline" of each cluster lane [i.e. the zero of the time averaged particle current inside the cluster zone; cf. Fig. 6(a)] to represent the entire (two-dimensional) lane. Since clusters moving on different lanes (and antiparallel directions) must not collide, each cluster lane on the system's periodic domain Ω must be closed. This can expressed by requiring where ∆ denote the distance along the x-axis between subsequent intersections of the cluster lane's midline with the lines y = 0 and y = L, and where w and l are coprime integers; cf. Fig. 11(a). Since all values of ∆ with ∆ 1 /L = L/∆ 2 are equivalent upon interchanging x-and y-axis, we can assert ∆ ≤ 1 without loss of generality. In Eq. (B1) the integers w and l characterize the topological and geometrical properties of the cluster lane, as is shown in Fig. 11(a). Specifically, w gives the cluster lane's winding number in the equatorial plane of the torus representing the periodic domain Ω; cf. denotes the cluster lane's orientation with respect to the y-axis. According to Eq. (B1), the ratio ∆/L can take on arbitrary rational values. At first sight, this does not seem to imply a severe restriction, since the set of rational numbers is a dense subset of the real numbers. In other words, arbitrary cluster lane orientations α are "infinitely close" to meeting the requirement Eq. (B1). While this is true in the hydrodynamic limit, L → ∞, finite size effects considerably restrict the number of possible cluster lane orientations α w,l . To see this, note that the area s cluster (α) each cluster lane of orientation α covers on Ω is given by (Λ: cluster lane's transversal extension) s cluster (α w,l ) = l · (α w,l ) · Λ.
(B3) Since the total system size is finite, ||Ω|| = L 2 , and since the minimum number of cluster lanes is 2 (one cluster lane alone corresponds to a wave pattern), we have Therefore, noting that the cluster lane width Λ is proportional to the transverse correlation length ξ ϕ of the polar order parameter, and that w ≤ l, we get For the system sizes used in our numerical computations, we found ξ ϕ ∼ L = O(10 2 ). We, therefore, observe cluster lane patterns characterized by 0 ≤ w ≤ l ∼ 1, i.e. w/l ∈ {0, 1}, amounting to cluster lane orientations α w,l ∈ {π/2, π/4}; cf. Figs. 11(b,c). This corresponds to the cluster lane geometries indicated in the top left plot in Fig. 12. Visual inspection of Fig. 12 illustrates the rapidly increasing space requirements for cluster lane patterns with increasing values of w and l. These space requirements, in turn, directly affect the competition between cluster lane and wave patterns during the pattern selection process: Although Eq. (B1) has to be fulfilled for any line, representing the lateral extension of a wave pattern, as well, the resulting space requirements [as captured by Eqs. (B3) and (B4)] are far less restrictive for waves. More precisely, a single wave band can attain any orientation α w,l such that s wave (α w,l ) = l · (α w,l ) · ξ L 2 , where ξ denotes the longitudinal extension of the wave. Since ξ Λ = O(ξ ϕ ), Eq. (B6) is therefore compatible with a much larger set of orientations α w,l than Eqs. (B3) and (B4) (for finite system sizes L). Finite system sizes L, therefore, imply a bias in favor of the formation of wave patterns.
Appendix C: Numerical protocol to measure the high density hysteresis loop In this appendix, we briefly explain the numerical protocol used to obtain the data points shown in Fig. 5 (open symbols). Special care must be taken to properly assess the location of the density scale ρ dh , where spatially homogeneous polar ordered states become linearly unstable toward the emergence of density segregated patterns. At the point where the homogeneous system turns linearly unstable, the development of spatial heterogeneities out of a virtually homogeneous state is extremely slow. As a sensitive means to detect the buildup of spatial heterogeneities, we introduce the max-imum absolute deviation, for the density matrix f α n , where This maximum absolute deviation approaches zero for spatially homogeneous systems and attains finite values [typically of order O(10 −1 )] for states inside the DSP regime.
To estimate the value of ρ dh , we proceed as follows. We start from a spatially homogeneous base state, equilibrated at ρ 0 = 0.415 with ∆ f = O(10 −8 ). We then quench the overall density ρ 0 to successively lower values and let the system evolve for a time span of at least T = 30000 λ −1 [55]. During this time, we record the evolution of the quantity ∆ f . For all data points in Fig. 5 with V hd /||Ω|| = 1, the maximum absolute deviation ∆ f is a monotonically decreasing function of time, and the system approaches a spatially homogeneous state. Conversely, in systems being represented by data points with V hd /||Ω|| < 1, ∆ f increases such that the spatially homogeneous initial states eventually give rise to density segregated patterns.
As mentioned earlier, for density quenches down to final densities ρ 0 ρ dh in the vicinity of ρ dh , the actual buildup of spatial heterogeneities happens on a very large time scale. In actual computations, this time scale turns out to be much slower than the local "equilibration" of the density matrix f α n toward the new density level ρ 0 . We, therefore, do not expect that the results shown in Fig. 5 should depend on whether they are obtained from a "quasi-static" protocol, or a density quench protocol as discussed here.