# A mean-field model for nematic alignment of self-propelled rods

Misha Perepelitsa111These authors contributed equally to this project
222corresponding authors Department of Mathematics, University of Houston, TX
77204  $\rm Ilya\,Timofeyev^{\textrm{a}}$ Department of Mathematics,
University of Houston, TX 77204  $\rm Patrick\,Murphy^{\textrm{a}}$ Department
of Bioengineering, Rice University, Houston, TX 77005  Center for Theoretical
Biological Physics, Rice University, Houston, TX 77005  $\rm
Oleg\,A.\,Igoshin^{\textrm{b}\,,}$ Department of Bioengineering, Rice
University, Houston, TX 77005  Center for Theoretical Biological Physics, Rice
University, Houston, TX 77005  Department of Chemistry, Rice University,
Houston, TX 77005  Department of Biosciences, Rice University, Houston, TX
77005

###### Abstract

Self-propelled rods are a facet of the field of active matter relevant to many
physical systems ranging in scale from shaken granular media and bacterial
alignment to the flocking dynamics of animals. In this paper we develop a
model for nematic alignment of self-propelled rods interacting through binary
collisions. We avoid phenomenological descriptions of rod interaction in favor
of rigorously using a set of microscopic-level rules. Under the assumption
that each collision results in a small change to a rod’s orientation, we
derive the Fokker-Planck equation for the evolution of the kinetic density
function. Using analytical and numerical methods, we study the emergence of
the nematic order from a homogeneous, uniform steady-state of the mean-field
equation. We compare the level of orientational noise needed to destabilize
this nematic order and compare our results to an existing phenomenological
model that does not explicitly account for the physical collisions of rods. We
show the presence of an additional geometric factor in our equations
reflecting a reduced collisions rate between nearly-aligned rods reduces the
level of noise at which nematic order is destroyed, suggesting that alignment
that depends on purely physical collisions is less robust.

##  I Introduction

Interactions between self-propelled rods have been the subject of many studies
on what is known as active matter. These studies were motivated by attempts to
understand interactions of rod-shaped cells [13, 4, 3, 34, 1] and motor-driven
filaments such as microtubules [21, 27, 28]. How the mechanisms of rod
interactions in such systems give rise to self-organization is a fundamental
question that continues to draw interest. For example, mechanical interactions
between the rods result in their local alignment and, under certain
conditions, the emergence of nematic order. These phenomena have been
extensively studied in the spirit of the alignment model of Vicsek et al.
[29], see for example Chaté et al. [12], Peruani et al. [11], Ginelli et al.
[10], Bolley et al. [5], and Degond et al. [22]. Such models typically assume
the collective alignment of nearby rods through mean nematic current or mean
nematic orientation, also known as the director. However, the rules proposed
in many of these studies lack first-principle derivations.

Boltzmann models of binary collisions provide a way to rigorously derive the
form of the collective alignment in scenarios where the rods interact
primarily through physical collisions rather than long range interactions. For
example, binary collision models for polar alignment of self-propelled
particles were studied by Bertin et al. [8]. The authors obtained a Boltzmann-
type equation for rod number density which they reduced to a system of
macroscopic (hydrodynamic) equations for the first two moments of the number
density. The steady-states of the latter system were used to extract
information about the phase transition to an ordered phase. In the context of
non-polar alignment, Hittmeir et al. [26] considered a model for the alignment
of myxobacteria, deriving a Boltzmann-type equation for the cell motion, a
closed system of macroscopic equations, and proving a number of analytical
results about the existence of solutions for these equations. Notably, this
model was based on instantaneous alignment and reversal from collisions
depending on the difference in orientations, and the rules for collisions were
structured to include reversals at the expense of nematic alignment.

Motivated by the mechanics of motion of Myxococcus xanthus [13, 19] and
monolayers of Bacillus subtilis [14, 34, 1], which mainly reorient through
physical-contact interactions, we consider a model for the collective behavior
of self-propelled rods in which nematic re-orientation occurs gradually over a
series of sequential binary collisions. These can be collisions between the
same two rods or with others nearby. Our goal is to derive from microscopic
collision rules a tractable continuous-time equation governing the evolution
of the rods’ spatial distribution, which we can then use to obtain
quantitative information about the impact of density and noise parameters on
self-organization. We assume that each collision results in the re-orientation
of rods by a certain amount, assumed to be small. Thus, the change in a rod’s
orientation over some time is the result of cumulative small changes due to
multiple collisions. The nematic alignment we consider differs from [26] in
this respect.

The paper is organized as follows. In section II, we derive an Enskog-type
equation for the rod number density and, in a certain asymptotic regime, the
mean-field Fokker-Planck equation. Our approach is similar to the treatment
for grazing collisions of molecules in gas dynamics [6]. We introduce an
averaging type agent-based model that, at the macroscopic level, is equivalent
to the binary collision model (BC model). The Fokker-Planck model obtained for
the BC-model turns out to be similar to the heuristic liquid crystal model (LC
model) of Peruani et al. [11]. Our kinetic equation differs from the latter by
a factor in the alignment rate that accounts for a decrease in collisions as
cells become nematically aligned. This naturally leads to the question of
comparing qualitative and quantitative differences between dynamics generated
by these two models. We address this question in section IV where we consider
the stability of a spatially homogeneous steady state with a uniform
distribution of orientations for both alignment models augmented by diffusion
(noise) in rod orientation. Linear stability can be calculated explicitly, as
was done in [11]. In section IV.2 we discuss the nonlinear stability for the
models of interest, using a numerical solver for the nonlinear mean-field
equations. We obtain parameter values for the phase transition to a
nematically ordered phase and compute steady-state orientation distributions
for several representative cases. Our findings show that self-propelled rods
are less ordered in the BC model compared to the LC model due to decreased
chance of collisions between cells with similar nematic orientations.
Therefore, the transition to the disordered phase occurs at a lower level of
rotational noise.

##  II Binary collision model

###  II.1 Derivation of Boltzmann-type equation

We consider a collection of $N$ rods moving on a square domain
$[0,L]\times[0,L]$ with periodic boundary conditions. Rods are rigid and move
along their longer axes, with their midpoint denoted by the coordinates
$\mathbf{x}$. Let $l$ be the length of the rod with corresponding unit
orientation vector $\mathbf{e}$. We denote the rod speed by $\bar{v}.$ In our
model we assume that the collision between a rod with orientation
$\mathbf{e}=\mathbf{e}(\theta)=(\cos\theta,\sin\theta)$, which we call a
$\theta$-rod, with a $\theta_{1}$-rod results in a nematic re-orientation of
both rods, yielding in new angles $\theta^{*},$ and $\theta^{*}_{1}$ according
to the rule

$\begin{split}\theta^{*}&=\theta+\delta\phi(\theta_{1}-\theta),\\\
\theta_{1}^{*}&=\theta_{1}-\delta\phi(\theta_{1}-\theta),\end{split}$ (1)

where $\phi(\theta)$ is a $\pi$–periodic, odd function, positive on the
interval $(0,\pi/2)$ and negative on the interval $(\pi/2,\pi),$ and
$\delta\in(0,1)$ is small numeric parameter that controls the strength of
alignment. We are implicitly assuming here that $\delta$ is the result of
rescaling $\phi$ to order $O(1)$. The function $\phi=\sin(2\theta)$ is a
typical example with a microscopic interpretation that we will use in section
IV. We will assume that a collision between rods occurs when the head of one
rod is in contact with the other (see Figure 1). We will assume that
interactions a predominantly binary and the contribution from simultaneous
interaction of more than two rods is negligible. We also allow for rods to
overlap each other spatially. We will refer to this model of alignment as a
binary collision (BC) model. The probabilistic description of alignment is
based on the one-particle distribution function of rod positions and
orientations, $f,$ with corresponding rod number density
$F(\mathbf{x},\theta,t)=Nf(\mathbf{x},\theta,t)$. As collisions happen between
pairs of particles, we will also need the 2-particle distribution function
$f^{(2)}(\mathbf{x},\theta,\mathbf{x}_{1},\theta_{1},t)$ which stands for the
density of the joint probability distribution for a random pair of cells, and
the corresponding function
$F^{(2)}(\mathbf{x},\theta,\mathbf{x}_{1},\theta_{1},t)=N(N-1)f^{(2)}(\mathbf{x},\theta,\mathbf{x}_{1},\theta_{1},t)$,
which stands for the number of pairs of cells (averaged over an ensemble) when
given two spatial regions and two sets of orientations.

We first must calculate the change in $F=Nf$ due to collisions. The value of
$f$ at a given orientation $\theta$ decreases when a $\theta$-rod interacts
with a $\theta_{1}$-rod ($\theta_{1}\not=\theta$) (Figure 1a) and increases
when a $\theta^{*}$-rod and a $\theta^{*}_{1}$-rod produce a $\theta$-rod
post-collision (Figure 1b). To describe the latter we need to swap the roles
of $\theta$ and $\theta^{*}$ in (1), then use the system to find $\theta^{*}$
in terms of the pre-collision orientation $\theta^{*}_{1}$ as well as the
given angle $\theta$. First,

$\begin{split}\theta&=\theta^{*}+\delta\phi(\theta_{1}^{*}-\theta^{*}).\end{split}$
(2)

It follows using substitution that
$\theta^{*}_{1}-\theta=\theta^{*}_{1}-\theta^{*}-\delta\phi(\theta^{*}_{1}-\theta^{*}).$
Denote the solution $\theta$ of

$\omega=\theta-\delta\phi(\theta)$

by $\psi_{\delta}(\omega)$. Then, we compute

$\partial_{\omega}\psi_{\delta}=(1-\delta\phi^{\prime}(\psi_{\delta}(\omega)))^{-1},$
(3)

which shows that the inverse function $\psi_{\delta}$ exists when $\delta$ is
sufficiently small. Given $\theta$ and an arbitrary pre-collisional
orientation $\theta_{1}^{*},$ we can determine the other pre-collisional
orientation $\theta^{*}$ so that $(\theta^{*},\theta_{1}^{*})$ collision
results in a $\theta$-rod, given by

$\theta^{*}=\theta_{1}^{*}-\psi_{\delta}(\theta_{1}^{*}-\theta).$ (4)

We now need to calculate the change in the number of rods over the time
interval $(t,t+\Delta t)$, denoted by $N\Delta f$. We do this by splitting the
change into decreasing and increasing contributions $\Delta
f=\Delta^{-}f-\Delta^{+}f$, with $\lim_{\Delta t\to
0}(\Delta^{-}f-\Delta^{+}f)/\Delta
t=(\frac{df}{dt})^{+}-(\frac{df}{dt})^{-}=\frac{df}{dt}$, the total or
Lagrangian derivative of $f$. This captures the total change in time of the
given quantity with respect to all variables that $f$ depends on, and can
alternatively be expressed via the chain rule as

$\frac{df}{dt}=\partial_{t}f+\sum_{i=1}^{2}\partial_{x_{i}}f\,\frac{dx_{i}}{dt}+\partial_{\theta}f\,\frac{d\theta}{dt}.$
(5)

In our context we have $\frac{d\theta}{dt}=0$ and
$(\frac{dx_{1}}{dt},\frac{dx_{2}}{dt})=\bar{v}(\cos(\theta),\sin(\theta))$ due
to our assumptions about the self-propelled nature of the rods.

First, we approximate the number of rods within a given spatial region and
orientation range by $Nf(\mathbf{x},\theta,t)d\mathbf{x}\,d\theta$. Figure 1a
shows the geometry involved when calculating possible collisions corresponding
to loss with $P_{1}$ and $P_{2}$ as the parallelograms where collisions can
occur with a given rod. Figure 1b corresponds similarly to gain of cells with
a given orientation with parallelograms $P_{1}$ and $P_{2}$. The magnitude of
decrease in the number of rods with orientation angle in
$(\theta,\theta+d\theta)$ in $\Delta t$ time equals the number of rods in this
range times the number of collisions, per rod, with another rod in
$(\theta_{1},\theta_{1}+d\theta_{1})$. Since the differentials $d\mathbf{x}$
and $d\theta$ ultimately drop out, we do keep track of them for these steps.
The expected number of collisions with rods in this range is given by summing
the expected number of pairs where the second rod is in $P_{1}$ or $P_{2}$.
This is made formal by integrating the pairwise joint density $N(N-1)f_{2}$
over these domains and over $\theta_{1}$.

$N\Delta^{-}f(\mathbf{x},\theta,t){}={}N(N-1)\left(\int_{-\pi}^{\pi}\int_{P_{1}}f^{(2)}(\mathbf{x},\theta,\mathbf{x}_{1},\theta_{1},t)\,d\mathbf{x}_{1}d\theta_{1}{}+{}\int_{-\pi}^{\pi}\int_{P_{2}}f^{(2)}(\mathbf{x},\theta,\mathbf{x}_{1},\theta_{1},t)\,d\mathbf{x}_{1}d\theta_{1}\right),$
(6)

where the parallelograms $P_{1}$ and $P_{2}$ are as in Figure 1.

Figure 1: Alignment geometry for a colliding rod changing from (A) or to (B) a
given orientation $\theta$. A) Parallelogram $P_{1}$ is formed by vectors
$l\mathbf{e}(\theta)$ and $\bar{v}\Delta t\,\mathbf{e}(\theta_{1}).$
Parallelogram $P_{2}$ is formed by vectors $l\mathbf{e}(\theta_{1})$ and
$\bar{v}\Delta t\,\mathbf{e}(\theta).$ A $\theta$–rod (solid arrow) re-orients
to a different orientation (dotted arrow) upon interaction with a
$\theta_{1}$-rod from $P_{1}$ or $P_{2}.$ $\theta_{1}$–rod re-orients by the
same angle in the opposite direction (dotted arrow). B) Same as A) but with a
collision in $P_{1}$ or $P_{2}$ resulting in a $\theta$ rod.

To obtain a closed-form equation for $f$, we assume the equivalent of
molecular chaos, expressed as statistical independence of the pairwise joint
distribution $f^{(2)}$ in terms of the single particle distribution

$f^{(2)}(\mathbf{x},\theta,\mathbf{x}_{1},\theta_{1},t)=f(\mathbf{x},\theta,t)f(\mathbf{x}_{1},\theta_{1},t).$
(7)

This yields

$\Delta^{-}f(\mathbf{x},\theta,t){}={}(N-1)\int_{-\pi}^{\pi}\left(\int_{P_{1}}f(\mathbf{x},\theta,t)f(\mathbf{x}_{1},\theta_{1},t)\,d\mathbf{x}_{1}{}+{}\int_{P_{2}}f(\mathbf{x},\theta,t)f(\mathbf{x}_{1},\theta_{1},t)\,d\mathbf{x}_{1}\right)d\theta_{1}.$
(8)

A change of variables from Cartesian $\mathbf{x}=(x_{1},x_{2})$ to an
$\mathbf{e}(\theta)$, $\mathbf{e}(\theta_{1})$ coordinate system gives an
extra factor of $|\sin(\theta_{1}-\theta)|$ in the integrand, taking into
account the relative nematic orientations when determining the collision
frequency of two rods with orientations $\theta$ and $\theta_{1}$. This is
analogous to the collision factor derived in [32]. Since the length of
integration in the $\mathbf{e}(\theta)$ direction is proportional to
$\bar{v}\Delta t$, by dividing both sides of (6) by $\Delta t$ and taking the
limit $\Delta t\to 0$, we obtain the magnitude of the negative contribution to
the rate of change of $f(\mathbf{x},\theta,t)$ over all possible $\theta_{1}$
as

$\bigg{(}\frac{df}{dt}\bigg{)}^{-}{}={}\bar{v}(N-1)\int_{-\pi}^{\pi}\int_{-l/2}^{l/2}|\sin(\theta_{1}-\theta)|f(\mathbf{x},\theta,t)\left(f(\mathbf{x}-\frac{l}{2}\mathbf{e}(\theta_{1})+z\mathbf{e}(\theta),\theta_{1},t)\right.\\\
\left.{}+{}f(\mathbf{x}+\frac{l}{2}\mathbf{e}(\theta)+z\mathbf{e}(\theta_{1}),\theta_{1},t)\right)\,dzd\theta_{1},$
(9)

where $z$ is the variable of integration in the
$\mathbf{e}(\theta_{1})$-direction.

Similarly, the increase $N\Delta^{+}f$ can be expressed using the expected
number of collisions for a rod in $(\theta^{*},\theta^{*}+d\theta^{*})$ with a
$(\theta_{1}^{*},\theta_{1}^{*}+d\theta_{1}^{*})$-rod in $\Delta t$ time such
that the result is a rod in $(\theta,\theta+d\theta)$. For this derivation, we
need to proceed carefully, taking into consideration the relationship between
$\theta$, $\theta^{*}$, and $\theta_{1}^{*}$. In particular we need use the
mapping (4) between $\theta$ and $\theta^{*}$ when relating
$f(\mathbf{x},\theta,t)d\mathbf{x}\,d\theta$ to
$(f(\mathbf{x},\theta^{*},t)d\mathbf{x}\,d\theta^{*})(f(\mathbf{x}_{1},\theta_{1}^{*},t)d\theta_{1}^{*})$.
Since the differential $d\mathbf{x}$ again drops out by the end, we do keep
track of it. As in the previous case, there are two types of collisions
depending on whether $\theta_{1}^{*}\in P_{1}$ or $\theta_{1}^{*}\in P_{2}$
(see Figure 1b):

$\Delta^{+}f(\mathbf{x},\theta,t)d\theta{}={}-(N-1)\int_{-\pi}^{\pi}\bigg{(}\int_{P_{1}}f(\mathbf{x},\theta^{*},t)d\theta^{*}f(\mathbf{x}_{1},\theta_{1}^{*},t)\,d\mathbf{x}_{1}\\\
{}+{}\int_{P_{2}}f(\mathbf{x},\theta^{*},t)d\theta^{*}f(\mathbf{x}_{1},\theta_{1}^{*},t)\,d\mathbf{x}_{1}\bigg{)}d\theta_{1}^{*}.$
(10)

Using (4) and (3) at $\omega=\theta_{1}^{*}-\theta$, we find

$d\theta^{*}{}={}\psi^{\prime}_{\delta}(\theta_{1}^{*}-\theta)d\theta=(1-\delta\phi^{\prime}(\psi_{\delta}(\theta_{1}^{*}-\theta)))^{-1}d\theta.$
(11)

Note that $\theta$ and $\theta_{1}^{*}$ are independent, so $d\theta$ in (10)
subsequently drops out. We then obtain the positive rate term in the limit
$\Delta t\to 0$ as

$\bigg{(}\frac{df}{dt}\bigg{)}^{+}{}={}\bar{v}(N-1)\int_{-\pi}^{\pi}\int_{-l/2}^{l/2}\psi^{\prime}_{\delta}(\theta_{1}^{*}-\theta)|\sin(\psi_{\delta}(\theta_{1}^{*}-\theta))|f(\mathbf{x},\theta^{*},t)\times\\\
\bigg{(}f(\mathbf{x}-\frac{l}{2}\mathbf{e}(\theta_{1}^{*})+z\mathbf{e}(\theta^{*}),\theta_{1}^{*},t){}+{}f(\mathbf{x}+\frac{l}{2}\mathbf{e}(\theta^{*})+z\mathbf{e}(\theta_{1}^{*}),\theta_{1}^{*},t)\bigg{)}\,dzd\theta_{1}^{*}.$
(12)

Finally, using that the total derivative of $f$ can be expressed as
$\partial_{t}f+\bar{v}\mathbf{e}\cdot\gradient_{\mathbf{x}}f$, we get the
kinetic equation

$\partial_{t}f{}+{}\bar{v}\mathbf{e}(\theta)\cdot\gradient_{\mathbf{x}}f\\\
{}={}(N-1)\bar{v}\int_{-\pi}^{\pi}\int_{-l/2}^{l/2}\psi^{\prime}_{\delta}(\theta_{1}^{*}-\theta)|\sin(\psi_{\delta}(\theta_{1}^{*}-\theta))|f(\mathbf{x},\theta^{*},t)\left(f(\mathbf{x}-\frac{l}{2}\mathbf{e}(\theta_{1}^{*})+z\mathbf{e}(\theta^{*}),\theta_{1}^{*},t)\right.\\\
\left.{}+{}f(\mathbf{x}+\frac{l}{2}\mathbf{e}(\theta^{*})+z\mathbf{e}(\theta_{1}^{*}),\theta_{1}^{*},t)\right)\,dzd\theta_{1}^{*}\\\
-(N-1)\bar{v}\int_{-\pi}^{\pi}\int_{-l/2}^{l/2}|\sin(\theta_{1}^{*}-\theta)|f(\mathbf{x},\theta,t)\left(f(\mathbf{x}-\frac{l}{2}\mathbf{e}(\theta_{1}^{*})+z\mathbf{e}(\theta),\theta_{1}^{*},t)\right.\\\
\left.{}+{}f(\mathbf{x}+\frac{l}{2}\mathbf{e}(\theta)+z\mathbf{e}(\theta_{1}),\theta_{1}^{*},t)\right)\,dzd\theta_{1}^{*}.$
(13)

Note that we also changed the notation from $\theta_{1}$ to $\theta_{1}^{*}$
in the last integral to match variables of integration. We will take care of
the remaining occurrences of $\theta^{*}$ in the next section’s asymptotics.

###  II.2 Asymptotic expansion and mean-field derivation

Let $T=\frac{L}{\bar{v}}$ be the characteristic time scale associated with
traversing the domain length $L.$ We scale the variables by

$\tilde{\mathbf{x}}=\frac{\mathbf{x}}{L},\quad\tilde{t}{}={}\frac{t}{T},\quad\tilde{f}(\tilde{\mathbf{x}},\theta,\tilde{t}){}={}L^{2}f(\mathbf{x},\theta,t).$
(14)

The kinetic equation scales accordingly and takes the following form (with a
slight abuse of notation, we keep using the old notation $\mathbf{x},t,f$ for
the new variables $\tilde{\mathbf{x}},\tilde{t},\tilde{f}).$

$\partial_{t}f{}+{}\mathbf{e}(\theta)\cdot\gradient_{\mathbf{x}}f\\\
{}={}(N-1)\int_{-\pi}^{\pi}\int_{-l/(2L)}^{l/(2L)}\psi^{\prime}_{\delta}(\theta_{1}^{*}-\theta)|\sin(\psi_{\delta}(\theta_{1}^{*}-\theta))|f(\mathbf{x},\theta^{*},t)\left(f(\mathbf{x}-\frac{l}{2L}\mathbf{e}(\theta_{1}^{*})+z\mathbf{e}(\theta^{*}),\theta_{1}^{*},t)\right.\\\
\left.{}+{}f(\mathbf{x}+\frac{l}{2L}\mathbf{e}(\theta^{*})+z\mathbf{e}(\theta_{1}^{*}),\theta_{1}^{*},t)\right)\,dzd\theta_{1}^{*}\\\
-(N-1)\int_{-\pi}^{\pi}\int_{-l/(2L)}^{l/(2L)}|\sin(\theta_{1}^{*}-\theta)|f(\mathbf{x},\theta,t)\left(f(\mathbf{x}-\frac{l}{2L}\mathbf{e}(\theta_{1}^{*})+z\mathbf{e}(\theta),\theta_{1}^{*},t)\right.\\\
\left.{}+{}f(\mathbf{x}+\frac{l}{2L}\mathbf{e}(\theta)+z\mathbf{e}(\theta_{1}),\theta_{1}^{*},t)\right)\,dzd\theta_{1}^{*}.$
(15)

Assuming small $\delta$, we introduce the series approximations

$\displaystyle\theta^{*}{}$ $\displaystyle=$
$\displaystyle{}\theta-\phi(\theta_{1}^{*}-\theta)\delta+O(\delta^{2}),$ (16)
$\displaystyle\psi^{\prime}_{\delta}(\theta_{1}^{*}-\theta){}$
$\displaystyle=$
$\displaystyle{}1+\phi^{\prime}(\theta_{1}^{*}-\theta)\delta+O(\delta^{2}),$
(17) $\displaystyle\psi_{\delta}(\theta_{1}^{*}-\theta){}$ $\displaystyle=$
$\displaystyle{}(\theta_{1}^{*}-\theta)+\phi(\theta_{1}^{*}-\theta)\delta+O(\delta^{2}),$
(18) $\displaystyle f(\mathbf{x},\theta^{*},t)-f(\mathbf{x},\theta,t){}$
$\displaystyle=$
$\displaystyle{}\partial_{\theta}f(\mathbf{x},\theta,t)(\theta^{*}-\theta)+O((\theta^{*}-\theta)^{2})$
$\displaystyle=$
$\displaystyle{}-\phi(\theta_{1}^{*}-\theta)\partial_{\theta}f(\mathbf{x},\theta,t)\delta+O(\delta^{2}).$

Additionally, for any continuous function $G(\theta)$, we can expand

$\int_{-\pi}^{\pi}|\sin(\psi_{\delta}(\theta_{1}^{*}-\theta))|G(\theta_{1}^{*})\,d\theta_{1}^{*}\\\
{}={}\int_{-\pi}^{\pi}\left(|\sin(\theta_{1}^{*}-\theta)|+{\rm
sign}\,(\sin(\theta_{1}^{*}-\theta))\cos(\theta_{1}^{*}-\theta)\phi(\theta_{1}^{*}-\theta)\delta\right)G(\theta_{1}^{*})\,d\theta_{1}^{*}\\\
{}+{}O(\delta^{2})\sup|G(\theta)|.$ (20)

Finally, we can expand $f$ in powers of $l/L$ to get

$f(\mathbf{x}-\frac{l}{2L}\mathbf{e}(\theta_{1}^{*})+z\mathbf{e}(\theta^{*}),\theta_{1}^{*},t)=f(\mathbf{x},\theta_{1}^{*},t)+O(lL^{-1}).$
(21)

Substituting all these expansions into the kinetic equation (15), we derive
the asymptotic relation

$\partial_{t}f{}+{}\mathbf{e}(\theta)\cdot\gradient_{\mathbf{x}}f{}={}-\frac{(N-1)l\delta}{L}\partial_{\theta}\left(f(\mathbf{x},\theta,t)\int_{-\pi}^{\pi}|\sin(\theta_{1}^{*}-\theta)|\phi(\theta_{1}^{*}-\theta)f(\mathbf{x},\theta_{1}^{*},t)\,d\theta_{1}^{*}\right)\\\
{}+{}O\left(NlL^{-1}\delta^{2}\right){}+{}O\left(Nl^{2}L^{-2}\delta\right).$
(22)

To finalize the model, we postulate two hypotheses. The first is that the
strength of alignment per collision, $\delta,$ is small and that
$Nl\delta{L^{-1}}$ is non-vanishing and finite. Since
$Nl\delta{L^{-1}}=(NL^{-2})(Ll)\delta,$ which can be interpreted as the
average number of collisions when traveling the domain length multiplied by
the maximum change per collision, this condition implies that many small
changes in a rod’s orientation accumulate to a finite magnitude. In this case
$\delta$ is inversely proportional to the number of rods in a narrow band of
width $l$. The second hypothesis is a macroscopic limit assumption:
$l{L}^{-1}$ is small and the number of rods, $N,$ is large. Under these
hypotheses, the leading order approximation of (22) is a mean-field, Fokker-
Planck equation

$\partial_{t}f{}+{}\mathbf{e}(\theta)\cdot\gradient_{\mathbf{x}}f{}+{}\kappa\partial_{\theta}\left(f(\mathbf{x},\theta,t)\int_{-\pi}^{\pi}|\sin(\theta_{1}-\theta)|\phi(\theta_{1}-\theta)f(\mathbf{x},\theta_{1},t)\,d\theta_{1}\right){}={}0,$
(23)

where

$\kappa=\frac{Nl\delta}{L},$ (24)

$\delta,\,l/L\ll 1$ and $N\gg 1.$ Note that unlike what one might expect,
$\kappa$ is not invariant to parameter changes that keep the mean density
$N/L^{2}$ fixed. This is due to the time rescaling being proportional to $L$.
If we double $L$ and keep the mean density fixed, $\kappa$ doubles as well,
but this is a reflection of one unit of time, and the average number of
collisions a rod would undergo in that time interval, also doubling. However,
our phase transition results in the next sections are independent of any
scaling that keep the mean density fixed.

In the derivation of the kinetic equation we assumed that $\delta$ is small
while $\frac{Nl\delta}{L}$ is finite. We would like to describe two model
scenarios in which these assumptions can be justified. The first situation is
a statistical black box approach. Suppose that one has a limited information
about the precise mechanism of alignment of self-propelled rods, except that
a: alignment is nematic and collisions are mostly binary; and b: the change of
angles in collisions is proportional to the angle between interacting cells.
Suppose that statistical analysis of an experiment reveals that there is some
characteristic length $L$ such that when the average angle between cell
orientations is $\theta_{0},$ the changes of cell orientation over distance
$L$ are proportional to $\theta_{0}.$ Estimating the average number of
collisions over distance $L$ as $\frac{Nl}{L},$ one can formulate an effective
law that each binary interaction changes a cell angle by the amount
$\delta\theta_{0},$ so that

$\frac{Nl}{L}\delta\theta_{0}\approx\theta_{0},$

which leads to the smallness of $\delta$ when $Nl/L\gg 1,$ i.e. the number of
collisions per length $L$ is large. The other situation, perhaps less
applicable to the motion of myxobacteria but still of sufficient interest, is
the case when the time of the contact of cells in collisions is short and the
torque generated during the collision is finite. For example, bacteria could
be able to crawl over one another without fully aligning. In this case,
$\delta$ can be taken proportional to the time of interaction, and $L$ to be
the distance over which changes in orientations are of finite magnitude.

###  II.3 Alignment to mean nematic orientation

In what follows, we replace the factor $N-1$ in (24) with $N$ since both have
the same asymptotic behavior for $N\gg 1.$ Following [11], we choose the
alignment function to be $\phi(\theta)=\sin(2\theta)$. This functional form is
commonly used and is related to frictional force balance, commonly found in
environments at low Reynolds’ number that bacterial cells experience. It can
be derived in different contexts, such as bacterial alignment in a colony due
to growth [2]. In our case it comes from the force a colliding rod receives
perpendicular to its orientation, proportional to $\sin(\theta)\cos(\theta)$,
from opposing frictional forces generated by the second rod.

With this choice, one can introduce a local mean nematic orientation and (23)
can be interpreted as a continuous alignment to such mean orientation. The
mean orientation is, in general, is specific for each rod, i.e. it depends on
$\theta$ through the rate of collisions with other rods. Indeed, we can write

$\int_{-\pi}^{\pi}|\sin(\theta_{1}-\theta)|\phi(\theta_{1}-\theta)f(\mathbf{x},\theta_{1},t)\,d\theta_{1}{}={}{\rm
Im}\,\int_{-\pi}^{\pi}|\sin(\theta_{1}-\theta)|e^{2i(\theta_{1}-\theta)}f(\mathbf{x},\theta_{1},t)\,d\theta_{1}\\\
{}={}{\rm
Im}\,\left[e^{-2i\theta}\int_{-\pi}^{\pi}|\sin(\theta_{1}-\theta)|e^{2i\theta_{1}}f(\mathbf{x},\theta_{1},t)\,d\theta_{1}\right]{}={}|{\bf
J}|\sin(2(\Theta-\theta)),$ (25)

where

${\bf
J}(\mathbf{x},\theta,t)=\int_{-\pi}^{\pi}|\sin(\theta_{1}-\theta)|e^{2i\theta_{1}}f(\mathbf{x},\theta_{1},t)\,d\theta_{1},$
(26)

is a mean nematic current, with polar angle $2\Theta(\mathbf{x},\theta,t)$
defined by the condition:

$e^{2i\Theta}=|{\bf
J}|^{-1}\int_{-\pi}^{\pi}|\sin(\theta_{1}-\theta)|e^{2i\theta_{1}}f(\mathbf{x},\theta_{1},t)\,d\theta_{1},$
(27)

or equivalently as an angle such that

$\int_{-\pi}^{\pi}|\sin(\theta-\theta_{1})|\sin(2(\theta_{1}-\Theta))f(\mathbf{x},\theta_{1},t)\,d\theta_{1}{}={}0.$
(28)

The kinetic equation in this case is

$\partial_{t}f{}+{}\mathbf{e}(\theta)\cdot\gradient_{\mathbf{x}}f{}+{}\kappa\partial_{\theta}\left(f|{\bf
J}|\sin(2(\Theta-\theta))\right){}={}0.$ (29)

The discussion of similar models based on the mean nematic current can be
found in Chaté et al. [12], Ginelli et al. [10] and Degond et al. [22].

##  III Agent-based models of alignment

In this section we describe a simplified agent-based model that can be used in
Monte Carlo simulations to obtain suitable approximate solutions of (23). In
this model we do not need to trace pairwise interactions of rods, but instead
use a more computationally efficient method of accumulating alignments over a
neighborhood of a rod. While this is the approach often take for
phenomenological models of alignment, we use the form of the mean nematic
current we derived from microscopic rules in previous sections to help ensure
consistency. Consider an agent-based model of $N$ rods moving on square domain
with side length $L$ with periodic boundary conditions. Let the positions of
rod centers $\\{\mathbf{x}_{i}\\}$ and orientation angles $\\{\theta_{i}\\}$
evolve according to the ODEs

$\begin{split}\frac{d\mathbf{x}_{i}}{dt}&=\bar{v}{\bf e}(\theta_{i}),\\\
\frac{d\theta_{i}}{dt}&=\gamma_{a}\sum_{|\mathbf{x}_{i}-\mathbf{x}_{j}|\leq
R}|\sin(\theta_{j}-\theta_{i})|\sin(2(\theta_{j}-\theta_{i}))+\sqrt{2D}\xi_{i}(t),\quad
i=1..N,\end{split}$ (30)

where $\gamma_{a}>0$ is the strength of alignment, $R$ is the interaction
radius, $D$ is the rotational diffusivity (diffusion coefficient), and
$\\{\xi_{i}(t)\\}_{i=1..N}$ is the vector of independent whites noises. The
above model is closely related to a model of nematic alignment of self-
propelled particles by Peruani et al. [11] that was built on the geometry of
alignment of liquid crystals. In particular it assumes only angular diffusion
so that rods move with constant velocity. We will refer to the model of
Peruani et al. [11] as the LC-model of alignment. The difference between the
two models is a factor $|\sin(\theta_{j}-\theta_{i})|$ in the rate of
alignment.

We would like to determine ranges of values of $\gamma_{a}$ and $R$ such that
the agent-based model (30) corresponds to BC-model from the last section. To
this end we compare (23) to the the mean-field equation describing the
dynamics of the probability density distribution function $f$ for model (30).
Following [11], the latter can be expressed as a Fokker-Plank equation

$\partial_{t}f{}+{}\bar{v}\mathbf{e}(\theta)\cdot\gradient_{\mathbf{x}}f\\\
{}+{}\gamma_{a}N\partial_{\theta}\left(V_{R}(\theta)f(\mathbf{x},\theta,t)\right){}={}D\partial_{\theta}^{2}f,$
(31)

with angular transport velocity

$V_{R}(\theta){}={}\int_{|\mathbf{x}_{1}-\mathbf{x}|<R}\int_{-\pi}^{\pi}|\sin(\theta_{1}-\theta)|\sin(2(\theta_{1}-\theta))f(\mathbf{x}_{1},\theta_{1},t)\,d\theta_{1}d\mathbf{x}_{1},$

found by rewriting the function underlying rate of change in angle due to
collisions $\sum_{|\mathbf{x}_{i}-\mathbf{x}_{j}|\leq
R}|\sin(\theta_{j}-\theta_{i})|\sin(2(\theta_{j}-\theta_{i}))$ in (30) as a
local integral with the number density $F=Nf$, then factoring out $N$. $D>0$
is the rotational diffusion coefficient from the rotational noise in the agent
based model.

Finally, rescaling $x$, $t$, and $f$ with (14) yields

$\partial_{t}f{}+{}\mathbf{e}(\theta)\cdot\gradient_{\mathbf{x}}f\\\
{}+{}\frac{\gamma_{a}NL}{\bar{v}}\partial_{\theta}\left(V_{R/L}(\theta)f(\mathbf{x},\theta,t)\right){}={}\frac{DL}{\bar{v}}\partial_{\theta}^{2}f.$
(32)

The additional local integration over a region
$|\mathbf{x}_{1}-\mathbf{x}|<\frac{R}{L}$ in $V_{R/L}$ differs from our model
(23), but they can be matched with a simple series expansion for
$f(\mathbf{x},\theta,t)$. Assuming $\frac{R}{L}\ll 1$, we can expand
$f(\mathbf{x}_{1},\theta,t)$ in a power series $f(\mathbf{x},\theta,t)+O(R/L)$
around $\mathbf{x}$ analogously to section II, yielding to first order

$\partial_{t}f{}+{}\mathbf{e}(\theta)\cdot\gradient_{\mathbf{x}}f\\\
{}+{}\frac{\gamma_{a}NL}{\bar{v}}\frac{\pi
R^{2}}{L^{2}}\partial_{\theta}\left(f(\mathbf{x},\theta,t)\int_{-\pi}^{\pi}|\sin(\theta_{1}-\theta)|\sin(2(\theta_{1}-\theta))f_{0}(\mathbf{x},\theta_{1},t)\,d\theta_{1}\right)=\frac{DL}{\bar{v}}\partial_{\theta}^{2}f.$
(33)

Note that in terms of parameters, (23) corresponds to model (33) when
$\frac{R}{L}\ll 1$ in such a way that $NR^{2}L^{-2}$ remains constant and

$\kappa{}={}\pi\frac{\gamma_{a}N}{\bar{v}}\frac{R^{2}}{L},$ (34)

from which we find that the rate of alignment $\gamma_{a}$ and interaction
radius $R$ are functionally dependent:

$\gamma_{a}R^{2}{}={}\frac{l\bar{v}\delta}{\pi}.$ (35)

##  IV Stability of a homogeneous uniform steady state

###  IV.1 Linear stability

Consider the mean-field equation (23) with added diffusion in orientation:

$\partial_{t}f{}+{}\mathbf{e}(\theta)\cdot\gradient_{\mathbf{x}}f{}+{}\kappa\partial_{\theta}\left(f(\mathbf{x},\theta,t)\int_{-\pi}^{\pi}|\sin(\theta_{1}-\theta)|\phi(\theta_{1}-\theta)f(\mathbf{x},\theta_{1},t)\,d\theta_{1}\right)\\\
{}={}\mu\partial^{2}_{\theta}f,$ (36)

where $\mu>0$ is a non-dimensional diffusion coefficient. Following [11], we
consider the problem of stability in the space of spatially homogeneous
solutions $f(\mathbf{x},\theta,t)=g(\theta,t).$ The BC-model’s equation for
$g$ is

$\partial_{t}g{}+{}\kappa\partial_{\theta}\left(g\int_{-\pi}^{\pi}|\sin(\theta_{1}-\theta)|\sin(2(\theta_{1}-\theta))g(\theta_{1},t)\,d\theta_{1}\right){}={}\mu\partial_{\theta}^{2}g.$
(37)

The function $g_{0}=1/(2\pi)$ is a steady-state solution of (37). By
introducing the perturbation from steady state $g=g_{0}+\epsilon g_{1}$ and by
linearizing equation (37) around $g_{0},$ we obtain

$\partial_{t}g_{1}{}+{}\frac{1}{2\pi}\partial_{\theta}\left(\int_{-\pi}^{\pi}|\sin(\theta_{1}-\theta)|\sin(2(\theta_{1}-\theta))g_{1}(\theta_{1},t)\,d\theta_{1}\right){}={}\eta\partial_{\theta}^{2}g_{1},$
(38)

with the rescaled viscosity coefficient $\eta=\frac{\mu}{\kappa}.$ Equation
(38) admits solutions in the form

$g_{1}=e^{\lambda t}e^{in\theta},\quad n\in\mathbb{Z},$ (39)

provided that

$\lambda{}={}\frac{\alpha_{n}}{\pi}-\eta n^{2},$ (40)

where ${\rm\bf i}\alpha_{n}$ is the eigenvalue corresponding to the
eigenfunction $e^{in\theta}$ of the operator
$I[g]{}={}\int_{-\pi}^{\pi}|\sin(\theta_{1}-\theta)|\sin(2(\theta_{1}-\theta))g(\theta_{1},t)\,d\theta_{1}.$
For all odd $n,$ $\alpha_{n}$=0, and for even $n$,

$\alpha_{n}{}={}2n\left(\frac{1}{(n-1)(n+1)}-\frac{1}{(n-3)(n+3)}\right).$
(41)

Among these values, only $\alpha_{2}=32/15$ is positive. Thus, when
$\eta<\eta_{c}=\frac{8}{15\pi},$ perturbations of the constant steady state
$g_{0}$ will grow in the direction of $\sin(2(\theta-\theta_{0})),$ for some
$\theta_{0}.$ For comparison, see [11], we recall that for the LC-model the
critical value is $\eta_{c}=0.25.$ Equation (37) can be used to show that
asymptotically, for $\eta$ close to $\eta_{c},$ the steady-state amplitude of
perturbations is of order $(\eta_{c}-\eta)^{1/2},$ i.e.

$g(\theta)\approx
g_{0}(\theta)+c_{0}(\eta_{c}-\eta)^{\frac{1}{2}}\sin(2(\theta-\theta_{0})).$
(42)

The stability condition $\eta<\eta_{c}$ can be written in terms of parameters
of the agent-based model as a balance between the diffusivity, the density,
and the strength of alignment:

$D>\frac{8}{15}\frac{N}{L^{2}}(\gamma_{a}R^{2}).$ (43)

In the linear case the instabilities grow without bound. The nonlinear
equation however will take the perturbations to a well defined steady state,
which we compute by solving (37) numerically in section IV.2.

###  IV.2 Nonlinear stability

####  IV.2.1 Numerical scheme

Equation (37) can be schematically written as

$\partial_{t}g+\partial_{\theta}\left(\mathbf{e}g\right)=\eta\partial_{\theta}^{2}g,$
(44)

where the velocity $\mathbf{e}\equiv\mathbf{e}(\theta,g,t)$ is computed using
a nonlocal integral operator. We use an operator-splitting method to
numerically integrate equations of his type. Thus, the advection and the
diffusion parts of this equation are integrated using different
discretizations. The diffusive part is integrated numerically using a standard
central finite-difference scheme. We use the local Lax-Friedrichs method from
LeVeque [24] for the advection part. In particular, if we consider the wave
equation written in flux form $\partial_{t}g+\partial_{\theta}F=0$ (where
$F\equiv\mathbf{e}g$), then we use following numerical method to advance the
solution during one time-step

$g_{k}^{j+1}=g_{k}^{j}-\frac{\Delta t}{\Delta\theta}[F_{k+1/2}-F_{k-1/2}]$
(45)

where $g_{k}^{j}=g(k\Delta\theta,j\Delta t)$ and numerical fluxes are given by

$F_{k+1/2}=\frac{1}{2}(v_{k}g_{k}^{j}+v_{k+1}g_{k+1}^{j})-\lambda_{i+1/2}(g_{k+1}^{j}-g_{k}^{j})$
(46)

with $v_{k}=v(k\Delta\theta,j\Delta t)$ and

$\lambda_{i+1/2}=\max(|v_{k}|,|v_{k+1}|)/2.$ (47)

We use $N=400$ points for the discretization and $\Delta t=2.5\times 10^{-3}$
in all simulations. We also verified our simulations with $N=800$ and $\Delta
t=6.25\times 10^{-4}$ for several selected simulations (several different
initial conditions) and did not observe any significant difference between
simulations with low and high resolutions. Thus, we can conclude that
simulations with $N=400$ and $\Delta t=2.5\times 10^{-3}$ are well-resolved
and the numerical diffusion due to the Lax-Friedrichs discretization of the
advection part is not significant.

####  IV.2.2 Numerical simulations

Figure 2: Evolution of $g(\theta,t)$ to a steady state in ustable regime. The
result of numerical simulation of equation (37) at times $t=0,\;10,\;13,\;25$
starting from random initial data. By time $t=25$ the numerical solution
reaches the steady state.

We used the above numercal scheme to obtain solutions of equation (37) and
study it’s dependence on the diffusivity parameter $\eta=\frac{\mu}{\kappa}.$
For initial data we selected $10^{4}$ random orientations from the uniform
distribution on the interval $[-\pi,\pi].$ The PDF was binned on the mesh in
the $\theta$ variable, with bin width $\Delta\theta.$ Figure 2 shows the
result of one such simulation in the unstable regime, when the solution
converges to a non-uniform steady state (the plot corresponding to $t=25$). We
preformed the comparison of numerical solutions of (37) with the numerical
solutions of the liquid crystal model. Figure 3 shows the difference in steady
states obtained from both models starting from the same initial data for two
different levels of noise. For lower noise levels, the BC solution reaches a
less ordered state due to the presence of the factor
$|\sin(\theta_{1}-\theta)|$ in the rate alignment, which decreases the
alignment strength among rods. Increasing the noise to the regime where the
binary collision model is completely disordered but the liquid crystal model
is not, we see the latter is perturbed from the uniform state in the direction
of the mode $\sin(2(\theta-\theta_{0}))$ for some phase $\theta_{0}$.

Figure 3: Steady states according to binary collision (BC) model and liquid
crystal (LC) model. Figure shows plots of function $g(\theta,t)$ at time $t$
when the solutions reach steady states, according to the numerical simulation
of equation (37) and the corresponding equation for the liquid crystal model.
The solutions start with the same random initial data. Left: $\eta=0.05$,
Right: $\eta=0.2$.  Figure 4: Nematic order parameter $S$ as function of
diffusivity $\eta.$ Figure shows nematic order $S$ at the steady states,
computed numerically, averaged over 20 different simulations with random
initial data, when $\eta$ ranges from 0.02 to 0.3 with an increment of 0.02.
Dotted lines mark the critical diffusivites for the corresponding linear
models.

Figure 4 shows the phase transition in the nematic order parameter $S(t)=|{\bf
S}(t)|,$ where

${\bf S}(t){}={}\int_{-\pi}^{\pi}e^{2{\rm\bf i}\theta}g(\theta,t)\,d\theta,$
(48)

is a function of diffusivity coefficient $\eta=\frac{\mu}{\kappa}$ for both
the BC and the LC models. In these figures, $S$ was computed when the solution
reached the steady state, and the average was taken over 1000 distinct
solutions with the initial data corresponding to $10^{4}$ rods with
orientations randomly selected uniformly in the interval $[-\pi,\pi].$ In the
ordered regime, the nematic order is stronger in LC model for the same value
of $\eta$, while for both models the nematic order decreases at approximately
the same rate at the phase transition.

##  V Discussion

We have developed an approach to modeling interactions of self-propelled rods
from a set of microscopic collision rules. Using a Boltzmann-like framework
[6], we derived a mean-field model under the assumption that binary collisions
are the dominant type of interaction, and that a collision result in a small
change in a rod’s orientation. We show for a common choice of the alignment
function that the resulting Fokker-Planck equation can be matched to commonly-
used mean-field models assuming alignment of rods in a local neighborhood
[11], albeit with a modified alignment function. These models have an
additional second-order term from added rotational diffusion. Restricting to
spatially homogeneous solutions, we calculate both analytically and
numerically when the constant solution with a uniform distribution in
orientation space becomes unstable, resulting in the emergence of nematic
order.

Our approach differs from standard treatments for this type of phenomenon [11,
8, 10, 12]. We avoid as much as possible postulating phenomenological rules
for how rods align, instead using first principles to derive the form of the
mean nematic current, which is a purely local quantity spatially. In
particular, our nematic current has an additional factor of
$|\sin(\theta_{2}-\theta_{1})|$ that takes into account the geometry of
collisions. In particular, it accounts for a decrease in the collision rate
between rods with similar nematic orientations $\theta_{1}$ and $\theta_{2}$,
as those rods are closer to being parallel with each other and less likely to
collide. Furthermore, we show that this additional factor does not change the
functional form of the phase transition measured via the nematic order
parameter $S$, but does notably raise the densities and lower rotational noise
levels at which the transition occurs. The physical interpretation is that a
decrease in collisions between rods with similar nematic orientation will
decrease the total rate of alignment in the population and make it more
susceptible to noise. As such, a higher density or lower rotational diffusion
is needed to produce the phase transition to an ordered state.

As a final note on our modeling approach, our treatment of $\delta$ as a small
parameter from rescaling $\phi$ is not needed under conditions where rods are
locally aligned. The function $\phi$ will naturally be small in such
situations assuming the change in orientation on collision is small for small
differences in angles, as a large discontinuity at $\phi(0)$ would imply rods
actively repel each other when nearly aligned. This suggests our model can be
used to analyze situations where the change is angle per collision is not
necessarily small, but rods are locally aligned into streams or clusters.

The alignment of active particles can be generated by two general mechanisms.
The first is interactions at a distance. These can be through hydrodynamic
interactions, e.g. alignment through long range interactions via forces in a
medium [25, 17, 9], or through more general flocking interactions [31, 30].
The second is short-ranged physical interactions, such as the type we
considered in this paper. They depend on the physical geometry of the
particles under consideration, arise from physical collisions, and often
appear when modeling granular media [20, 16] or biological phenomena. These
phenomena range in scale, but are generally found on either the macromolecular
scale when considering interaction of biological filaments [21, 27, 28] or
cellular populations [13, 7, 3, 18, 2, 1]. At larger scales, organisms’
sensory organs tend to result in longer-ranged interactions.

The effect the presence our geometric factor has on the influence of level of
noise has some interesting implications for active nematics. Since it lowers
the threshold at which nematic order is destroyed, this suggesting that
systems that depend on purely physical collisions to align are less robust to
effects that introduce noise into the rods’ orientations. In the context of
cellular interactions, this could influence emergent behaviors such as
swarming [1, 14] or aggregation [15, 7, 23] that often depend on some manner
of cellular alignment. The noise itself could come externally such as from
heterogeneity or external forces in the medium the cells are interacting with,
or from internal noise coming from the biochemistry of the cell. In either
case, having an alternative alignment mechanism that acts at at longer range,
such as through hydrodynamic forces or by actively shaping the environment
through the use of an extracellular matrix or biofilm [33, 3, 18], provides a
greater level of robustness for when alignment is desired.

##  VI Acknowledgments

Research is supported by National Science Foundation Division of Mathematical
Sciences awards 1903270 (to M.P. and I.T.) and 1903275 (to OAI and PM). We
also would like to thank prof. D. Kuzmin (TU Dortmund) for helpful discussions
on the numerical methods for the kinetic model presented in this paper.

## References

  * [1] Be’er A, Ilkanaiv B, Gross R, Kearns D B, Heidenreich S, Bär M, and Ariel G.  A phase diagram for bacterial swarming.  Nat Comm Physics, 3(66), 2020.
  * [2] D. AU Dell’Arciprete, M. L. Blow, A. T. Brown, F. D. C. Farrell, J. S. Lintuvuori, A. F. McVey, D. Marenduzzo, and W. C. K. Poon.  A growing bacterial colony in two dimensions as an active nematic.  Nat Commun, 9(4190), 2018.
  * [3] Rajesh Balagam and Oleg A. Igoshin.  Mechanism for Collective Cell Alignment in Myxococcus xanthus Bacteria.  PLOS Computational Biology, 11(8):e1004474, August 2015.
  * [4] Rajesh Balagam, Douglas B. Litwin, Fabian Czerwinski, Mingzhai Sun, Heidi B. Kaplan, Joshua W. Shaevitz, and Oleg A. Igoshin.  Myxococcus xanthus Gliding Motors Are Elastically Coupled to the Substrate as Predicted by the Focal Adhesion Model of Gliding Motility.  PLOS Computational Biology, 10(5):e1003619, May 2014.
  * [5] François Bolley, José A. Cañizo, and José A. Carrillo.  Mean-field limit for the stochastic Vicsek model.  Applied Mathematics Letters, 25(3):339–343, March 2012.
  * [6] Cercignani C.  The Boltzmann Equation and Its Applications.  Springer-Verlag, New York, 1988.
  * [7] Patrick D. Curtis, Rion G. Taylor, Roy D. Welch, and Lawrence J. Shimkets.  Spatial Organization of Myxococcus xanthus during Fruiting Body Formation.  J Bacteriol, 189(24):9126–9130, December 2007.
  * [8] Bertin E, Droz M, and Grégoire G.  Boltzmann and hydrodynamic description for self-propelled particles.  Phys Rev E Stat Nonlin Soft Matter Phys, 74, 2006.
  * [9] J. Elgeti, R. G. Winkler, and G. Gompper.  Physics of microswimmers, single particle motion and collective behavior: A review.  Rep. Prog. Phys., 78(5):056601, April 2015.
  * [10] Ginelli F, Peruani F, Bär M, and Chaté H.  Large-scale collective properties of self-propelled rods.  Phys Rev Lett, 104, 2010.
  * [11] Peruani F, Deutsch A, and Bär M.  A mean-field theory for self-propelled particles interacting by velocity alignment mechanisms.  Eur. Phys. J. Spec. Top., 157:111–122, 2008.
  * [12] Chaté H, Ginelli F, and Montagne R.  Simple model for active nematics: quasi-long-range order and giant fluctuations.  Phys Rev Lett, 96:180602, 2006.
  * [13] D Kaiser.  Coupling cell movement to multicellular development in myxobacteria.  Nat Rev Microbiol, 1(1):45–54, October 2003.
  * [14] Daniel B. Kearns and Richard Losick.  Swarming motility in undomesticated bacillus subtilis.  Molecular Microbiology, 49(3):581–590, 2003.
  * [15] Maria A. Kiskowski, Yi Jiang, and Mark S. Alber.  Role of streams in myxobacteria aggregate formation.  Phys. Biol., 1(3):173–183, October 2004.
  * [16] Arshad Kudrolli, Geoffroy Lumay, Dmitri Volfson, and Lev S. Tsimring.  Swarming and Swirling in Self-Propelled Polar Granular Rods.  Phys. Rev. Lett., 100(5):058001, February 2008.
  * [17] Nitin Kumar, Harsh Soni, Sriram Ramaswamy, and A. K. Sood.  Flocking at a distance in active granular matter.  Nat Commun, 5(1):4688, September 2014.
  * [18] Xuefei Li, Rajesh Balagam, Ting-Fang He, Peter P. Lee, Oleg A. Igoshin, and Herbert Levine.  On the mechanism of long-range orientational order of fibroblasts.  Proceedings of the National Academy of Sciences, 114(34):8974–8979, August 2017.
  * [19] José Muñoz-Dorado, Francisco J Marcos-Torres, Elena García-Bravo, Aurelio Moraleda-Muñoz, and Juana Pérez.  Myxobacteria: Moving, killing, feeding, and surviving together.  Front Microbiol, 7:781, 2016.
  * [20] Vijay Narayan, Sriram Ramaswamy, and Narayanan Menon.  Long-Lived Giant Number Fluctuations in a Swarming Granular Nematic.  Science, 317(5834):105–108, July 2007.
  * [21] F. J. Ndlec, T. Surrey, A. C. Maggs, and S. Leibler.  Self-organization of microtubules and motors.  Nature, 389(6648):305–308, September 1997.
  * [22] Degond P, Manhart A, and Yu H.  A continuum model for nematic alignment of self-propelled particles.  Discrete Contin. Dyn. Syst. Ser. B, 22:1295–1327, 2017.
  * [23] Fernando Peruani, Jörn Starruß, Vladimir Jakovljevic, Lotte Søgaard-Andersen, Andreas Deutsch, and Markus Bär.  Collective Motion and Nonequilibrium Cluster Formation in Colonies of Gliding Bacteria.  Phys. Rev. Lett., 108(9):098102, February 2012.
  * [24] LeVeque R.  Finite Volume Methods for Hyperbolic Problems.  Cambridge University Press, 2002.
  * [25] Ingmar H. Riedel, Karsten Kruse, and Jonathon Howard.  A Self-Organized Vortex Array of Hydrodynamically Entrained Sperm Cells.  Science, 309(5732):300–303, July 2005.
  * [26] Hittmeir S, Kanzler L, Manhart A, and Schmeiser Ch.  Kinetic modelling of colonies of myxobacteria.  Kinetic and Related Models, 14:1–24, 2021.
  * [27] Volker Schaller, Christoph Weber, Christine Semmrich, Erwin Frey, and Andreas R. Bausch.  Polar patterns of driven filaments.  Nature, 467(7311):73–77, September 2010.
  * [28] Yutaka Sumino, Ken H. Nagai, Yuji Shitaka, Dan Tanaka, Kenichi Yoshikawa, Hugues Chaté, and Kazuhiro Oiwa.  Large-scale vortex lattice emerging from collectively moving microtubules.  Nature, 483(7390):448–452, March 2012.
  * [29] Vicsek T, Czirók A, Ben-Jacob E, Inon Cohen I, and Shochet O.  Novel type of phase transition in a system of self-driven particles.  Phys. Rev. Lett., 75, 1995.
  * [30] John Toner and Yuhai Tu.  Flocks, herds, and schools: A quantitative theory of flocking.  Phys. Rev. E, 58(4):4828–4858, October 1998.
  * [31] Tamás Vicsek, András Czirók, Eshel Ben-Jacob, Inon Cohen, and Ofer Shochet.  Novel Type of Phase Transition in a System of Self-Driven Particles.  Phys. Rev. Lett., 75(6):1226–1229, August 1995.
  * [32] Christoph A Weber, Florian Thüroff, and Erwin Frey.  Role of particle conservation in self-propelled particle systems.  New Journal of Physics, 15(4):045014, apr 2013.
  * [33] Charles Wolgemuth, Egbert Hoiczyk, Dale Kaiser, and George Oster.  How Myxobacteria Glide.  Current Biology, 12(5):369–377, March 2002.
  * [34] H. P. Zhang, Avraham Be’er, E.-L. Florin, and Harry L. Swinney.  Collective motion and density fluctuations in bacterial colonies.  Proceedings of the National Academy of Sciences, 107(31):13626–13630, 2010.

Figure 5: Rate of decay of the nematic order parameter $S$ as function of
diffusivity $\eta$ for $\eta<\eta_{c}$. Red dashed line corresponds to the fit
$\log(S)\sim 0.5\log(\eta_{c}-\eta)$. Note that
$\log(\eta_{c}-\eta)\approx-4.5$ corresponds to $\eta=1.6$.