A continuum model for nematic alignment of self-propelled particles

A continuum model for a population of self-propelled particles interacting through nematic alignment is derived from an individual-based model. The methodology consists of introducing a hydrodynamic scaling of the corresponding mean-field kinetic equation. The resulting perturbation problem is solved thanks to the concept of generalized collision invariants. It yields a hyperbolic but non-conservative system of equations for the nematic mean direction of the flow and the densities of particles flowing parallel or anti-parallel to this mean direction. Diffusive terms are introduced under a weakly non-local interaction assumption and the diffusion coefficient is proven to be positive. An application to the modeling of myxobacteria is outlined.


Introduction
Systems of large numbers of self-propelled particles are commonly found in nature. In three space dimensions, well-known examples include fish schools, flocks of birds and insects [8,10,39]. In two space dimensions, typical examples are bacteria, such as myxobacteria and Bacillus subtilis [5,33,42,46,48]. However, the interest in such systems is not limited to biological examples. It has been found that the underlying principles also often apply to human walking behavior [30,35] and physical systems, such as granular media [26,37]. The main interest in these systems, both mathematically and from an application's point of view, stems from their emergent properties. Emergence refers to the fact that fairly local interaction rules between the agents like attraction, alignment, repulsion, etc. give rise to large scale collective behavior. Emergence is physically described as a phase transition to an ordered state (see e.g. the review [45]). Mathematical models often rely on a particle (agent) based description [1,11,13,41]. Among these models, the time dependent Vicsek model [44] and variants of it have received much attention owing to its structural simplicity, while convincingly describing a wide range of phenomena. The general idea of the Vicsek model is to describe the agents as particles, moving with constant speed. At each time step each particle adjusts its direction towards the mean direction of the particles in its neighborhood subject to some random perturbation.
Whilst particle models are a good tool to reproduce and study the observed phenomena, the high number of agents involved make simulations computationally costly considering the amount of detail one is interested in. Hydrodynamic (macroscopic) models are typically easier to simulate and analyze and thereby offer a powerful alternative. Therefore a lot of effort has been made to formulate macroscopic equations starting from particle models [3,6,9,12,21,23,28,34]. Mathematically one of the main challenges is the lack of conserved quantities, which is typically a key ingredient in the computation of the macroscopic limit (see again [45]). This lack of conservation is frequent in biological systems, which, as opposed to many physical systems, e.g. don't conserve momentum. One strategy to overcome this is the use of generalized collision invariants (GCI), introduced for the first time in [23] to derive a hydrodynamic model of the Vicsek dynamics (later on referred to as the Self-Organized Hydrodynamic (SOH) model). The passage from individual-based to hydrodynamic models involves the construction of an intermediate kinetic model. For classical physical systems, momentum or energy conservation translates at the kinetic level into the concept of collision invariants. This is generalized to the GCI concept if the underlying system does not satisfy these conditions. The rigorous passage from kinetic to SOH models has been achieved in [36]. The GCI concept has been used in a variety of cases [17,16,20,28] and even applied to a model of economics [22]. Modal analysis of the SOH model can be found in [25]. One important feature of the Vicsek model is the emergence of phase transitions. These are studied in the kinetic of macroscopic framework in [2,6,18,19,43]. Here, we concentrate on the derivation of macroscopic models and postpone the analysis of phase transitions to future works.
In this paper we want to apply this method to nematic alignment of polar, rod-shaped particles moving with constant speed in two space dimensions. Nematic alignment describes the situation where upon meeting, two particles either align if their mutual angle is acute or anti-align if their angle is obtuse. This nematic alignment rule typically results from volume exclusion and can be understood as an inelastic collision between polar particles. It differs from polar alignment of polar particles (which is akin to ferromagnetic interactions in spin systems, as is the case of the Vicsek model) and nematic alignment of apolar particles (active nematics, which is typically observed in liquid crystals [14], in fibrous tissues [15] or in systems of interacting disks [24]). Large collections of polar particles which nematically align have been observed to exhibit interesting macroscopic patterns such as high density traveling bands and their properties differ from the other types of alignment [29,38]. Therefore the need for a hydrodynamic model derived from particle dynamics arises and this work suggests how to derive such a model. In the derivation we account for purely local interactions as well as weakly non-local interactions which leads to an additional diffusion term in the macroscopic limit.
To give an example how the nematic SOH model can be applied to biological systems, we consider the case of myxobacteria. Due to their ability to produce several macroscopic structures, such as waves, spirals and clusters [27,33,46], they have already received a lot of attention in terms of mathematical models [7,31,32,47]. Since this type of bacteria can internally reverse its polarity, the model has to be complemented with a reaction term, describing these reversals. We demonstrate how such a term can be incorporated in a systematic way into the original model.
The rest of the document is structured as follows: In Section 2 we describe the particle model for nematic alignment and the resulting mean-field model. The derived alignment operator Q 0 al is examined in Section 3. Whilst it shares some properties with its polar alignment counterpart, there are also some differences, e.g. the set of equilibria has dimension three (instead of two). Mathematically more unexpected, showing the existence and uniqueness of solutions of the related linear operator, which amounts to proving a Poincaré inequality, is slightly more complex than in the polar case [23]. The section finishes with explicitly characterizing the equilibria and generalized collision invariants of Q 0 al . With these ingredients we derive the macroscopic model in Section 4 and prove its hyperbolicity. In Section 5 we assume weakly non-local interactions and derive the corresponding macroscopic model. Finally in Section 6 we apply both macroscopic models to the case of myxobacteria by additionally describing density dependent reversal at the particle level. In the macroscopic model this leads to a reaction term which structurally hints at interesting dynamics expected to take place.

Particle Model
As a first step towards a macroscopic model, we start with the description of a timecontinuous particle (or agent-based) model. The model's structure is similar to the timecontinuous version of the Vicsek model presented in [23] and we will discuss similarities and differences. We image the particles to resemble stiff rods of equal length and describe the i-th particle (i ∈ N := {1, . . . , N }) at time t ≥ 0 by its center of mass X i (t) ∈ R 2 and its orientation v(Θ i (t)), where we define as the polar vectors associated with the angle θ in a fixed reference frame and where Θ i (t) is defined modulo 2π. All particles are assumed to move in the direction of their orientation with constant speed, v 0 . In particular this means that they cannot drift sideways (as opposed to the model in [4]). Under these hypotheses, the motion of the particles is supposed to satisfy the following system of stochastic differential equations is the average nematic alignment direction. It is defined as follows: First define a mean 'nematic' current where R > 0 is the interaction radius. Next defineΘ i (t) as an angle of lines (i.e. an angle defined modulo π) such that: We note that the definition of the vector on the left-hand side of (2.3) does not depend on the choice of the reference frame used to measure the angles Θ i asΘ i (t) can be equivalently defined by and we notice that this relation is invariant under a uniform translation of all angles. To obtain a unique representation of the angles, we will use the convention that the particle angles Θ i (t) are taken in [−π, π) and the nematic anglesΘ i (t) are taken in [−π/2, π/2). For simplicity, we will identify [−π, π) with R/2πZ and with the unit circle S 1 . The notation dB i t describes independent Brownian motions with intensity d cos 2 (Θ i −Θ i ) with d > 0. We are aiming for a separation of the particles into two (locally defined) groups, one with cos(Θ i −Θ i ) > 0 and one with cos(Θ i −Θ i ) < 0. Therefore, diffusion must vanish near cos(Θ i −Θ i ) = 0, to avoid particles crossing the boundary between these two angular regions. The factor cos 2 (Θ i −Θ i ) is intended to avoid that the two groups mix. The first term in the velocity evolution describes alignment with the bacteria within the interaction radius R > 0 with an alignment frequency ν Sign(cos(Θ i −Θ i )) sin(Θ i −Θ i ) with ν ≥ 0. Due to the factor Sign(cos(Θ i −Θ i )), alignment occurs parallel or antiparallel to the mean angleΘ i , depending on whether Θ i is closer toΘ i orΘ i + π. Note that Sign(cos(Θ i −Θ i )) sin(Θ i −Θ i ) is a π-periodic function and so, does not depend on the choice ofΘ i orΘ i + π as the representative of the corresponding angle of lines.
Remark 2.1 Also other notions of nematic alignment between several particles are possible. In [40] the nematic mean direction was defined as The main difference between this definition and (2.3) is that the nematic angle defined in (2.3) depends on the particle index i only via the position X i (t), but is independent of the direction of the i-th particle, Θ i . The outcome of (2.4) on the other hand depends on both the position and the direction of the i-th particle. Whilst the definition from [40] might be closer to reality, if the dispersion around the mean angle is small, both types of nematic alignment lead to the same behavior. However, the independence from Θ i makes (2.3) mathematically much easier to handle.

Mean Field Model
We follow the usual procedure to derive the mean field limit as described for example in [23]. The procedure starts with the definition of the empirical distribution function f N (x, θ, t), with x ∈ R 2 , θ ∈ R (modulo 2π) and t > 0, according to where δ (X i (t),Θ i (t)) (x, θ) is the Dirac delta in (x, θ) space located at (X i (t), Θ i (t)) at time t.
In the absence of noise (d = 0), f N is a deterministic measure that satisfies a first order kinetic equation if we assume that all (X i , V i ) fulfill (2.1)-(2.3). If d = 0, f N is a random measure whose law satisfies a second order kinetic equation. Then as N → ∞, f N approximates a distribution function f (x, θ, t), satisfying the following Kolmogorov-Fokker-Planck type equation: where the collision operatorQ al is given bỹ and whereΘ f is given by We note from (2.7) thatΘ f is an angle of lines, defined modulo π and a representative ofΘ f will be chosen in [−π/2, π/2). We also note that (2.7)-(2.8) can be equivalently written At this point the alignment interaction is non-local which is reflected in the fact that the computations for J f (x, t) involve integrals over a ball centered around x. As a next step we introduce dimensionless variables. Since we are interested in macroscopic phenomena, we use a hydrodynamic scaling. In the following we will do the non-dimensionalization and the hydrodynamic scaling in one step. A more detailed description for a similar equation can be found in [21]. On a microscopic scale we use the reference time t 0 = 1/ν and reference length x 0 = v 0 t 0 . We also introduce a scaled diffusion constant D = d t 0 .
On the macroscopic scale we want to use units t 0 = t 0 ε and x 0 = x 0 ε , which are large compared to the microscopic ones. We introduce the new (dimensionless, macroscopic) variablesx = x .
Note that the scaling of R supposes that the interactions between the particles are local. Another scaling which takes better account of the non-locality of the interactions will be investigated later on (see Section 5). With this scaling and dropping the hats for better readability we find that f = f ε satisfies the problem: where (2.14) Again,θ f is an angle of lines, defined modulo π and a representative ofθ f is chosen in [−π/2, π/2). As above, (2.13)-(2.14) can be equivalently written Dropping the terms of order ε 2 or smaller, we find that f = f ε (x, θ, t) satisfies the perturbation problem: 3 Properties of the alignment operator Q 0 al
We now give a more precise functional setting in which to define the collision operator Q 0 al . For this purpose, we introduce the following: ) be a given angle of lines. The linear operator L θ 0 is defined as

5)
for any smooth function ϕ defined on R/(2πZ) with values in R.
(ii) The space H θ 0 is defined by It is a Hilbert space when endowed with the inner product and associated norm g 2 where ∂ θ g is defined in the distributional sense. It is a Hilbert space when endowed with the inner product and associated norm g 2 for any ϕ, ψ ∈ V θ 0 . Here, ·, · V θ 0 ,V θ 0 denotes the duality bracket between V θ 0 and its dual. This duality bracket extends the inner product of By abuse of notation, we will simply write (3.6) as follows: Now we investigate the solvability of the equation for a given function ψ. In view of (3.6), this problem can be recast into the following variational formulation: For this variational formulation, we have the following Proposition 3.4 (i) (Homogeneous problem) Suppose ψ = 0. Then, the unique solutions of (3.9), (3.10) are piecewise constant functions, i.e. there exist C ± ∈ R such that (ii) (Nonhomogeneous problem) Let ψ ∈ V θ 0 . Then, there exists a solution to (3.9), (3.10) if and only if the following two solvability conditions hold simultaneously: Under these solvability conditions, there exists a unique solution ϕ to (3.9), (3.10) which satisfies the additional conditions and the general solution in V θ 0 to (3.9), (3.10) is the sum of the unique solution satisfying (3.13) and the general solution (3.11) of the homogeneous problem.
Proof. We easily verify that we can choose θ 0 = 0 and thus we drop the indices θ 0 in the definitions of the operators and spaces.
We now show that the variational formulation (3.14), (3.15) admits a unique solution. For this purpose, we apply Lax-Milgram's Lemma. First, we note thatṼ is a closed subspace of V. Indeed, it is obvious that the linear forms ϕ → ± cos θ>0 ϕ M dθ are continuous. Therefore,Ṽ is a Hilbert space. Then, it is clear that the left-hand side of (3.15) defines a symmetric continuous bilinear form a(ϕ, χ) onṼ and that the right-hand side of (3.15) defines a continuous linear form L(χ) onṼ. To apply Lax-Milgram's Lemma, it remains to prove that a is coercive. This amounts to proving a Poincaré inequality, i.e. that there exists a constant C > 0 such that for all ϕ ∈Ṽ we have: We show that cos θ>0 Obviously, if such a result is true, it is also true if the integration domains are changed to the set {cos θ < 0} in both integrals, and this implies (3.16). We now concentrate on obtaining (3.17). Let ψ ∈Ṽ and let θ 1 and θ 2 be such that cos θ 1 > 0, cos θ 2 > 0. We have: Multiplying by M (θ 1 ), integrating with respect to θ 1 on (−π/2, π/2) and using that Taking the square, multiplying by M (θ 2 ) and integrating with respect to θ 2 on (−π/2, π/2) leads to: and by the Cauchy-Schwarz inequality with respect to the inner product on H, again using that cos θ 1 >0 M (θ 1 ) dθ 1 = 1, we obtain where I 1 , . . . I 4 are the integrals with respect to the sets Ω 1 to Ω 4 with Ω 1 = {| sin ∈ Ω 1 , using Cauchy-Schwarz inequality, we have: and consequently, (3.20) where here and in the following, we denote by C generic constants. Take now (θ 1 , θ 2 ) ∈ Ω 4 . Then, using Cauchy-Schwarz inequality, we can write: Now, we remark that sin θ cos 2 θ , and deduce that Then, we have The first term is bounded on Ω 4 . For the second term we remark that Now, 1/ cos θ is a convex function and θ ∈ [θ 1 , θ 2 ]. Assuming, without loss of generality, that θ 1 < θ 2 , we have: and we deduce that Therefore, M (θ 1 )M (θ 2 )/M (θ) ≤ 1 and since cos θ sin 2 θ is also bounded for θ such that sin θ > 1/2, we deduce that the second term of (3.23) is also bounded. Consequently, for (θ 1 , θ 2 ) ∈ Ω 4 , we get which implies, upon integrating with respect to (θ 1 , θ 2 ), that Now, we take (θ 1 , θ 2 ) ∈ Ω 2 (obviously the case of (θ 1 , θ 2 ) ∈ Ω 3 can be treated similarly). In this case, I 2 is decomposed into I 2 = I 1 2 + I 2 2 where I 1 2 and I 2 2 correspond to integrals over the domains Ω 1 2 = (−π/6, π/6) × (−π/2, −π/6) and Ω 2 2 = (−π/6, π/6) × (π/6, π/2) respectively. We consider I 2 2 and a similar proof can be made for I 1 2 . For (θ 1 , θ 2 ) ∈ Ω 2 2 , we can decompose: and consequently we have: (3.26) Each of the terms will be evaluated separately. For the first term, we proceed as for I 1 and notice that M (θ) and M (θ) | cos θ| are bounded from above and from below away from zero, so that we can estimate it similarly to (3.19) and its contribution to I 2 2 leads to a similar estimate as (3.20). For the second term, we use (3.21) with θ 1 replaced by π/6 and likewise, we find (3.22) still with θ 1 replaced by π/6. We can keep on going through all the steps of the proof and we eventually get an estimate similar to (3.24) with θ 1 replaced by π/6. This shows that the second term of (3.26) leads to an estimate similar to (3.25). Therefore, we get |∂ θ ψ(θ)| 2 M (θ) cos 2 θ dθ. (3.27) Now, inserting (3.20), (3.25), (3.27) and its counterpart for I 3 into (3.18) leads to the Poincaré estimate (3.17) and eventually to (3.16), which ends the proof of the Poincaré estimate.
(iii) We restore the indices θ 0 . Then it is obvious that ψ and ϕ satisfy the conditions (3.12), (3.13) respectively and that they belong to V θ 0 . So, they belong toṼ θ 0 . Furthermore, by the change of variable θ = θ + θ 0 in the variational formulation satisfied by ϕ 0 (i.e. with θ 0 = 0), we easily see that ϕ satisfies the variational formulation (3.9), (3.10). Since the solution of this variational formulation in the spaceṼ θ 0 is unique, ϕ is this solution.

Equilibria
We now define the notion of equilibria and characterize them: We denote by E the set of equilibria of Q 0 al .

Generalized Collision Invariants
We first define the notion of a collision invariant, thanks to the We deliberately leave the function setting vague as this notion will have to be enlarged later on. Clearly, the set of collision invariants is a vector space. At this point, however, we only have a one-dimensional space of collision invariants, spanned by the constant function Ψ(θ) = 1, corresponding to mass conservation. Following the approach in [23], we therefore generalize the concept of a collision invariant and define a generalized collision invariant: Definition 3.9 Given an angle of lines θ 0 ∈ [−π/2, π/2), a function Ψ ∈ V θ 0 is a generalized collision invariant (GCI) associated with θ 0 if and only if Here the equalityθ f = θ 0 is an equality of angles of lines i.e. holds modulo π. The set of GCI associated with θ 0 is a linear space denoted by G θ 0 .
We note that for f ∈ V θ 0 such that Since Ψ being a CI equivalently means that we pass from the definition of a CI to that of a GCI by replacing the orientationθ f of f by an arbitrary orientation θ 0 and requesting that the integral is zero only for those functions f whose mean orientationθ f coincides with θ 0 . Because this means we restrict the set of functions f , we expect to find more Ψ, fulfilling the definition of a GCI as compared to that of a CI. We can reformulate the definition of a GCI in the following way: Lemma 3.10 Given an angle of lines θ 0 ∈ [−π/2, π/2), a function Ψ ∈ V θ 0 is a GCI associated with θ 0 if and only if: Proof. We note that the constraintθ f = θ 0 can be equivalently rephrased as the linear constraint: By a classical duality argument already used in [23], statement (3.31) is equivalent to the existence of a constant a ∈ R such that Using the formal self-adjointness of L θ 0 (see Definition 3.3 (iv)), this is equivalent to which implies that Ψ satisfies (3.32). The converse is obvious.
(ii) The functions χ ± θ 0 and g θ 0 can be written: In the remainder, we will write χ ± and g for χ ± 0 and g 0 . In addition, g satisfies the symmetry conditions: Proof. (i) First, we note that the solution space of Problem (3.32) is the linear space spanned by the solution of the same problem with a = 1. Therefore, we restrict ourselves to a = 1. Now, solving Problem (3.32) with a = 1 amounts to solving the variational formulation (3.9), (3.10) with right-hand side ψ(θ) = sin 2(θ − θ 0 ). We apply Prop. 3.4 (ii). A necessary and sufficient condition for the existence of a solution is that the right-hand side satisfies the solvability condition (3.12), which in the present case, with ψ = sin 2(θ − θ 0 ) reads: This condition is obviously satisfied by symmetry. Therefore, there exists a unique solution g θ 0 of the problem satisfying the cancellation condition (3.13) and all solutions of the problems differ from this one by the addition of a solution of the homogeneous problem, i.e. of the type (3.11). Such solutions are arbitrary linear combinations of χ ± θ 0 . (ii) The statement is obvious for χ ± θ 0 . For g θ 0 , we simply apply Prop. 3.4 (iii). Finally, we remark that sin 2θ satisfies the symmetry conditions (3.34). Therefore, using the uniqueness part of Prop. 3.4 (ii), we deduce that g satisfies the same symmetry conditions.
In the present simple two-dimensional case, we can have an explicit formula for g, given by the Lemma below and depicted in Figure 2.

Derivation
This section is concerned with the limit ε → 0 in Problem (2.15)-(2.16). The first step is the observation that the dominant term Q 0 al f ε = 0 requires the limiting distribution to lie in the kernel of Q 0 al , characterized by Lemma 3.7. To determine the equations satisfied by the macroscopic quantities ρ + (x, t), ρ − (x, t) andθ(x, t), it is necessary to consider the terms of order ε in (2.15). This is done by multiplication with the three independent GCI given in Prop. 3.11 and subsequent integration over θ. The computations are similar to those in [23] and are described in more detail in the Appendix (Section 8.1). Here, we simply state the result: Theorem 4.1 Taking the (formal) limit ε → 0 in (2.15)-(2.16), we obtain where the macroscopic quantities ρ ± (x, t) andθ(x, t) have values in [0, ∞) and [−π/2, π/2) respectively and fulfill the following system of equations: where we recall that v(θ) = (cosθ, sinθ) and v ⊥ (θ) = (− sinθ, cosθ). The coefficients d 1 , d 2 and µ are given by: where g is the GCI defined in Prop. 3.11 and explicitly given by (3.35). For a function ϕ(θ), we denote by ϕ M the following average: where M is Mθ forθ = 0.
The proof can be found in the Appendix (Section 8.1).
As shown in the proof of Lemma 3.7, the macroscopic quantityθ(x, t) is the local mean angle of lines as defined in (2.13) of the equilibrium distributionf . The functions ρ + (x, t) and ρ − (x, t) describe the local densities of particles that point in direction [θ(x, t) − π/2,θ(x, t) + π/2) and [θ(x, t) + π/2,θ(x, t) + 3π/2) respectively. Eqs. (4.37) and (4.38) simply describe number density transport in the direction of v(θ) (for ρ + ) or −v(θ) (for ρ − ). These equations express local conservations of the number densities of the particles moving in the direction of v(θ) or −v(θ) respectively. Since GCI do not express local conservations (unless they are CI), in general one cannot expect that they will lead to conservation laws. Here, Eqs. (4.37) and (4.38) are conservation laws in spite of the fact that they are derived from using the GCI χ ± θ f . This unusual circumstance is due to the particular expression of the collision operator Q 0 al . Note that whenever ρ + = ρ − , it follows that ∂ tθ = 0 and ∂ t (ρ + + ρ − ) = 0. An alternative formulation of System (4.37)-(4.39) can be obtained by defining ρ = ρ + + ρ − and δ = ρ + − ρ − . We get: It is instructive to compare with the Self-Organized Hydrodynamic (SOH) system derived in [23], which, in dimension two, takes the form: where the constantsd 1 ,d 2 andμ depend on κ and the collision invariants of the corresponding alignment operator. This model is based on polar alignment. Here by contrast to System (4.37)-(4.39),θ is an angle of vectors, i.e.θ ∈ R/(2πZ) and the particles are distributed according to a classical von Mises-Fisher distribution (i.e. their angular distribution is proportional to e κ cos(θ−θ) ) and they move (with a certain dispersion described by κ) in the direction ofθ.
The fact the nematic alignment leads to an equilibrium angle of lines as opposed to an angle of vectors is reflected in the existence of two opposing local maxima atθ and θ + π of the GVM Mθ(θ), whilst the classical von Mises-Fisher distribution has only one. The factor cos 2 (θ −θ f ) in front of the diffusion operator in Q 0 al allows for this separation into two groups and causes Mθ(θ) = 0 for θ =θ ± π/2, i.e. in equilibrium there are no particles moving perpendicular to the equilibrium angle of lines. On the other hand, both the von Mises-Fisher distribution and the GVM share their general "bump-like" shape and the convergence to a Dirac delta for κ → ∞. At the level of the macroscopic equations further similarities and differences become evident. In the nematic case both groups of particles have their own conservation law, consequently we have three instead of two macroscopic equations. The general structure remains the same: Transport equations for the density complemented by an equation describing the angle evolution, consisting of a time derivative, a convection term and a pressure term. However for the nematic SOH model both the convection and the pressure term depend on the local density difference δ. In the case where there is only one group of particles, e.g. ρ + = ρ and ρ − = 0 system (4.37)-(4.39) reduces to (4.45)-(4.46) (with different constants), since then the nematic character of the collisions does not play a role.

Hyperbolicity of the Macroscopic Model
An important property of the original SOH system (4.45)-(4.46) is its hyperbolicity. In [28] a variant of the SOH model was examined. It was found that in some regions hyperbolicity is lost, which affects the behavior of the solutions. Therefore, we examine the hyperbolicity of the macroscopic model (4.37)-(4.39). To simplify the analysis we rescale the system with t → d 1 t. In the following we therefore abbreviated 2 = d 2 /d 1 andμ = µ/d 1 . Proof. We use the alternative formulation (4.42)-(4.44), which leads to simpler formulas.
To understand its roots, we use the following definition of a discriminant of a third order polynomial: For p(z) = a 3 z 3 + a 2 z 2 + a 1 z + a 0 we define Next we observe that the discriminant of the characteristic polynomial of A, which we define as ∆ = discr charPoly A , can be interpreted as a quadratic polynomial in X = (δ/ρ) 2 with coefficients dependent on c = cos 2 θ, namely: If ∆(X) > 0 for all c ∈ [0, 1] and all X ∈ [0, 1] we have three distinct, real eigenvalues and the system is (strictly) hyperbolic. The roots of ∆(X) are given by Simple calculations show that for c > 9(μ+d 2 ) 9μ+μd 2 +8d 2 the discriminant ∆(X) has no real roots and is positive for all X ∈ [0, 1]. Otherwise we can distinguish between c > s :=μ µ+d 2 and 0 ≤ c ≤ s. In the first case both roots X 1,2 (c) are negative are therefore not relevant. In the second case we lose hyperbolicity for X > X 2 (c). However a simple analysis shows that X 2 (c) ≥ 1 for all 0 ≤ c ≤ s. In fact X 2 (c) → ∞ for c → 0 + and c → s − and min c∈[0,s] X 2 (c) = 1. Since X ∈ [0, 1] this shows the hyperbolicity of the system.

Weakly non-local scaling
In this section we want to investigate how the model changes, if we allow for a relatively larger interaction radius. With the scaling choice below, this leads to interactions that are no longer local in the limit, but rather weakly non-local. This and other scaling choices for different models of self-propelled particles have been examined in [21]. In the scaling performed in Section 2.2, we usedR = O(ε). If we useR = O( √ ε) instead, namelŷ R = √ ε r, we are led to the following expression for J ε f : andΘ ε f is still given by (2.11). Then, we need to consider higher order terms in the Taylor expansion of J ε f andΘ ε f . These are given by the following:

The last line (5.61) can be split into
We start by rewriting A 1 using integration by parts. This yields where the second equality is obtained by noting that the boundary term vanishes. Note that since dg dθ ≤ 0 on [0, π/2], which can be seen from its definition in (3.35), it also holds that We continue by expanding A 2 1 : which shows that in fact A 2 1 = −A 2 as defined in (5.62). Using this and (5.63) we can conclude Finally, since g ≤ 0 on [0, π/2], we observe that the denominator B as defined in (5.60) (again omitting Z κ ) fulfills which shows the claim and finishes the proof.

Example: Myxobacteria
In this section we want to give an example of how to use and adapt the derived nematic SOH model to a particular biological problem, the movement of myxobacteria. This type of rod-shaped bacteria can glide without the aid of flagella on surfaces. Under starvation conditions, prior to the formation of fruiting bodies, they go through a socalled rippling phase. During this stage traveling waves of high density bands are observed that seemingly pass through each other unaffected. However upon close examination it has been found that a large number of bacteria at the interface of two colliding waves in fact reverse their direction. These reversals can happen spontaneously, but a higher number of head-to-head contacts with other, oppositely moving myxobacteria increases their likelihood [32,46]. Biologically these reversals are quite fascinating, because they happen by the reorganization of the internal movement machinery as opposed to an actual turning process. The interaction of two bacteria that don't reverse can be described by nematic alignment. In terms of other communication between bacteria no diffusible chemotactic signal is known to coordinate their movement. At this point the nematic SOH model can describe the (nematic) alignment of two colliding bacteria, however the reversals still need to be included.

Particle and Mean Field Model Particle Model
We complete the description of the particle model from Subsection 2.1 by assuming that particles reverse their direction of motion with a frequency, dependent on the density of oppositely moving bacteria ρ i ± . So we let Θ i be changed into Θ i + π at Poisson random times with frequency λ(ρ i −Sign cos(Θ i −Θ i ) ), where where Card stands for the cardinal of a set. In other words, we have . In these formula, the function λ(ρ ± ) is a reversal frequency. The way it depends on ρ ± will be discussed later.

Mean Field Model
The above considerations change equation (2.5) to where the collision operatorsQ al is given by (2.6) and and whereΘ f is given by (2.7)-(2.8) and σ f,± given by To obtain local interactions, we assume as in Section 2.2 an interaction radius fulfillingR = O(ε) and choose for the reversal frequency the scalingλ(σf ,± ) = t 0 λ(σ f,± ) = t 0 ε λ(σ f,± ). The local density is then given by and it can be shown that This changes equation (2.15) into where Q 0 al f is given by (2.16) and whereθ f and ρ f,± are given by (2.13) and (6.68) respectively.
Remark 6.1 In (6.69) the contributions from Q 0 rev appear at the same scale as the transport operator, which is a consequence of the scaling choice for the reversal rate λ. Usinĝ λ(σ f,± ) = t 0 λ(σ f,± ) instead, would lead to a different model, in which Q 0 al and Q 0 rev would be of the same order. Consequently also the shape of the equilibria and the macroscopic model would look very different. We leave this path for future work.
The proof consists of some simple additions to the proof of Theorem 4.1 and is given in the Appendix (see Section 8.2).
Macroscopic Model -Weakly non-local Scaling To obtain the hydrodynamic limit of the myxobacteria model for the weakly non-local scaling described in Section 5, we note that in this scaling It holds that We observe that (5.55) is changed into with Q 1 al and Q 0 rev given by (5.54) and (6.70) respectively. The resulting limit is summarized in Proposition 6.2 Taking the (formal) limit ε → 0 in (6.69)-(6.70), we obtain f ε (x, θ, t) →f ρ + (x,t),ρ − (x,t),θ(x,t) (θ), withf ρ + ,ρ − ,θ given by (3.28). The macroscopic quantities ρ ± andθ fulfill where d 1 , d 2 and µ are given by (4.40) and d 3 is given by (5.59). The definition of . M can be found in (4.41).
The proof is a simple combination of the proofs of Theorem 5.3 and Proposition 6.1.
One can observe that the additional reversal term leads to a reaction term on the right-hand side of equations (6.71) and (6.72), with the precise shape depending on the choice of λ(ρ). Motivated by biological findings we expect λ(0) = λ 0 > 0 to account for spontaneous reversals in absence of other bacteria. Since head-to-head contacts increase the reversal frequency, λ(ρ) should be an increasing function of ρ. To demonstrate what type of dynamics could be expected, we give an example, where we choose with λ 0 , λ 1 > 0. For a constant direction,θ(x, t) ≡θ and moving along the characteristics, the corresponding phase portrait is depicted in Figure 3. One can observe that wherever the total density ρ + + ρ − > 2 λ 0 λ 1 one of the two groups dominates the other one, whilst for smaller values of ρ + + ρ − both attain the same number. This is similar to the waves observed in myxobacteria, where at the high density peaks of the waves, almost all bacteria move in the same directions, whilst in between the waves, there can be small numbers of bacteria moving in different directions. A detailed analytical and numerical analysis for this and other choices of λ(ρ) will be examined in a forthcoming work.

Conclusion
In this work we have derived macroscopic equations for self-propelled, rod shaped particles moving in two space dimensions starting at the particle level. The type of interaction considered is nematic alignment of polar particles, which is commonly found in biological and physical systems, but systematically derived hydrodynamic equations for this type of alignment are largely missing from the literature. In order to take the limit in the hydrodynamic scaling, we used the concept of generalized collision invariants. Two types of scalings have been examined: a local scaling and a weakly non-local scaling. The resulting macroscopic equations share some similarities with their polar alignment counterparts, but differ significantly in other aspects. A key difference is the splitting of the particles into two groups moving into opposite directions. Lastly, including density dependent reversals, we have demonstrated how the equations can be used to describe the behavior of myxobacteria. The research in the area of pattern formation and non-equilibrium states of large assembles of particles is very active and more macroscopic models are needed to better understand the phase transitions and underlying causes. Several extensions of this work are possible: to describe a wider range of phenomena it could be generalized to three and more space dimensions. Further at this point the limits examined are on a formal level and lack a rigorous treatment. Extensive numerical experiments need to be conducted to compare the particle model with the macroscopic models and to characterize the different phases and states. To better understand the rippling behavior of myxobacteria different choices of reversal rates can be examined, both numerically and analytically. and, by Definition 3.9, this integral equals zero. Upon division by ε we are left with Letting ε → 0, using (8.76) and the fact thatθ f ε →θ, we are led to ± cos(θ−θ)>0 Note that the integration domain depends on x and t. We calculate which can also be expressed as: We note that in spite of the presence of indicator functions in (3.29),f ρ + ,ρ − ,θ is continuous (and even C ∞ ) across the points θ where cos(θ −θ) = 0, so that there is no delta distribution appearing at these points when we differentiate it. In the following we perform the calculations for the integral over θ such that cos(θ −θ) > 0, the case of cos(θ −θ) < 0 being treated in the same manner. The left-hand side of (8. Thanks to the normalization condition (3.3), we have L 1 = ∂ t ρ + . Using that v(θ) = cos(θ −θ) v(θ) + sin(θ −θ) v ⊥ (θ), (8.80) together with cancellations due to symmetries, L 2 simplifies into where d 1 is defined in (4.40). We easily get L 3 = 0 by symmetry. Finally, using again (8.80), cancellations due to symmetry and integration by parts, we find that Finally, collecting the contributions from above, we get (4.37) and (4.38).

Use of GCI gθ f ε
We start as in the previous section. We interpret the integral of Q 0 al f ε against gθ f ε as and again by Definition 3.9 this integral vanishes. Upon division by ε we get ∂ t f ε + v(θ) · ∇ x f ε gθ f ε dθ = 0.
and finallyθ