# Weakly nonlinear analysis of pattern formation in active suspensions

Laurel Ohm1<EMAIL_ADDRESS>Michael J. Shelley2,3  1Department of
Mathematics, Princeton University, Princeton, NJ 08540 2Courant Institute of
Mathematical Sciences, New York University, New York, NY 10012 3Center for
Computational Biology, Flatiron Institute, Simons Foundation, New York, NY
10010

###### Abstract

We consider the Saintillan–Shelley kinetic model of active rodlike particles
in Stokes flow (Saintillan & Shelley, 2008a, b), for which the uniform,
isotropic suspension of pusher particles is known to be unstable in certain
settings. Through weakly nonlinear analysis accompanied by numerical
simulations, we determine exactly how the isotropic steady state loses
stability in different parameter regimes. We study each of the various types
of bifurcations admitted by the system, including both subcritical and
supercritical Hopf and pitchfork bifurcations. Elucidating this system’s
behavior near these bifurcations provides a theoretical means of comparing
this model with other physical systems which transition to turbulence, and
makes predictions about the nature of bifurcations in active suspensions that
can be explored experimentally.

###### keywords:

##  1 Introduction

Inherently nonequilibrium suspensions of active particles abound in biological
and experimental settings (Gompper et al., 2020; Marchetti et al., 2013). For
example, motile bacteria such as _E. coli_ and _Bacillus subtilis_ propel
themselves through their surrounding fluid environment, interacting through
their induced flow fields (Lauga & Powers, 2009; Lushi et al., 2014; Mendelson
et al., 1999; Pedley & Kessler, 1992), while likewise immersed microtubule
bundles slide and extend, driven by ATP-driven molecular motors (DeCamp et
al., 2015; Henkin et al., 2014; Needleman & Dogic, 2017; Opathalage et al.,
2019; Sanchez et al., 2012). These active suspensions are remarkable because,
despite the near lack of inertial effects relative to viscous ones, the
activity of the particles can produce large-scale coherent flows and even so-
called _active_ or _bacterial turbulence_ , characterized by chaotic
fluctuations in particle concentration and fluid velocity (Dombrowski et al.,
2004; Doostmohammadi et al., 2017; Dunkel et al., 2013; Gachelin et al., 2014;
Koch & Subramanian, 2011; Nishiguchi et al., 2017; Peng et al., 2021; Simha &
Ramaswamy, 2002; Sokolov & Aranson, 2012; Sokolov et al., 2007; Stenhammar et
al., 2017; Thampi et al., 2014; Zhang et al., 2010).

Here we consider the kinetic model for a dilute suspension of active elongated
particles developed by Saintillan and Shelley (Saintillan & Shelley, 2008a,
b). We note that a similar model was developed independently by Subramanian &
Koch (2009), and that both models share many similarities with the Doi–Edwards
model for passive polymers (Doi, 1981; Doi & Edwards, 1986). In the dilute
limit, particles only interact with each other hydrodynamically by exerting an
‘active stress’ on the surrounding fluid. Even in this setting, changes in
particle density and activity are known to cause the suspension to transition
from a uniform, isotropic state to more complex states – many of which are
observed experimentally – involving large-scale patterns in particle alignment
and concentration. The kinetic model (Saintillan & Shelley, 2008a, b) thus
presents an opportunity to elucidate some of the fundamental mechanisms behind
the transition to collective dynamics and active turbulence.

To understand this transition, we perform a multiple timescales expansion to
determine exactly how the uniform, isotropic steady state in the 2D kinetic
model loses stability in different parameter regimes. The variety of predicted
behaviors near the onset of instability, which we verify through numerical
simulations, is surprisingly complex. Linear stability analysis alone shows
that, depending on the (fixed) ratio of particle diffusivity to concentration,
the uniform isotropic state can lose stability through either a pitchfork or
Hopf bifurcation. Here the bifurcation parameter is a ratio of the particle
swimming speed to the particle concentration and magnitude of active stress
the particles exert. Our weakly nonlinear analysis shows that both the
pitchfork and Hopf parameter regions can be further subdivided into
subcritical and supercritical regions, again depending on the ratio of
particle diffusivity to concentration.

Numerically, we find hysteresis in the subcritical Hopf region, where far-
from-isotropic quasiperiodic patterns of particle alignment are bistable with
the uniform isotropic state. The patterns in this region are perhaps
precursors to active turbulence. However, the dimensionality of the initial
perturbation to the isotropic steady state makes a difference. If the initial
perturbation is one dimensional, i.e. purely in the $x$ or $y$ direction, then
only a supercritical Hopf bifurcation can occur, and we numerically locate the
stable limit cycle that arises. An example of a 2D limit cycle is also located
numerically within the region in which both 1D and 2D perturbations give rise
to a supercritical Hopf bifurcation. In the supercritical pitchfork setting,
which includes immotile (but active) particles, we identify the stable steady
states emerging just beyond the bifurcation. These steady states resemble the
steady vortex found in Wioland et al. (2013), although we consider periodic
boundary conditions rather than confinement.

A key takeaway is that even this simple kinetic model is capable of capturing
many different types of transitions to collective behavior in an active
suspension. The different bifurcations analyzed here can be compared with
systematic numerical studies of phase transitions in other active suspension
models (Forest et al., 2004a, b; Giomi et al., 2011, 2012; Xiao-Gang et al.,
2014; Yang et al., 2014) and help explain the 1D and 2D patterns – including
limit cycles and other attractors – located numerically in Forest et al.
(2014, 2013, 2015) for a similar version of the kinetic model. The weakly
nonlinear analysis performed here also facilitates comparison with the normal
forms arising in more classical pattern formation processes in fluid
mechanics, especially thermal convection (Crawford & Knobloch, 1991; Cross &
Hohenberg, 1993; Knobloch, 1986; Pomeau, 1986; Schöpf & Zimmermann, 1993;
Swift & Hohenberg, 1977; Swinney & Gollub, 1981), but also other phenomena
arising in complex fluids such as electrohydrodynamic convection in nematic
liquid crystals (Bodenschatz et al., 1988) and the transition from sub- to
supercritical instability in viscoelastic pipe flow (Wan et al., 2021).

The paper begins by introducing the kinetic model (section 2.1) and recapping
the well-studied linear stability analysis (sections 2.2, 2.3) and role of
rotational diffusion (section 2.4). Readers familiar with the model may wish
to skip directly to the outline of results in section 2.5, where the types of
bifurcations are mapped out in greater detail.

##  2 Background

###  2.1 Kinetic model of an active suspension

In the kinetic model (Saintillan & Shelley, 2008a, b), a suspension of $N$
elongated particles in a 2-dimensional periodic box of length $L$ is modeled
as a number density $\Psi(\boldsymbol{x},\boldsymbol{p},t)$ of particles with
center-of-mass position $\boldsymbol{x}$ and orientation $\boldsymbol{p}\in
S^{1}$ at time $t$. Due to conservation of the number of particles, the
density $\Psi$ evolves according to a Smoluchowski equation:

$\p_{t}\Psi=-\bnabla\bcdot(\dot{\boldsymbol{x}}\Psi)-\bnabla_{p}\bcdot(\dot{\boldsymbol{p}}\Psi);\quad\bnabla_{p}:=({\boldsymbol{I}}-\boldsymbol{p}\boldsymbol{p}^{\rm
T})\p_{\boldsymbol{p}}.$ (1)

Here $\bnabla_{p}\bcdot$ denotes the divergence on the unit sphere. The
translational and rotational fluxes are given by

$\displaystyle\dot{\boldsymbol{x}}$
$\displaystyle=V_{0}\boldsymbol{p}+\boldsymbol{u}-D_{T}\bnabla(\log\Psi)$ (2)
$\displaystyle\dot{\boldsymbol{p}}$
$\displaystyle=({\boldsymbol{I}}-\boldsymbol{p}\boldsymbol{p}^{\rm
T})(\bnabla\boldsymbol{u}\,\boldsymbol{p})-D_{R}\bnabla_{p}(\log\Psi).$ (3)

The translational velocity consists of a particle swimming term with speed
$V_{0}$ in direction $\boldsymbol{p}$, particle advection by the surrounding
fluid flow, and translational diffusion. For simplicity we take the
translational diffusion to be isotropic. The rotational velocity depends on a
Jeffery term for the rotation of an elongated particle in Stokes flow
(Jeffery, 1922), written here in the infinitely slender limit, along with
rotational diffusion.

Finally, the surrounding fluid medium satisfies the Stokes equations with
active forcing:

$\displaystyle-\mu\Delta\boldsymbol{u}+\bnabla q$
$\displaystyle=\bnabla\bcdot{\boldsymbol{\Sigma}}^{a},\;\bnabla\bcdot\boldsymbol{u}=0$
(4) $\displaystyle{\boldsymbol{\Sigma}}^{a}$
$\displaystyle=a_{0}\int_{S^{1}}\Psi(\boldsymbol{x},\boldsymbol{p},t)(\boldsymbol{p}\boldsymbol{p}^{\rm
T}-\textstyle\frac{1}{2}{\boldsymbol{I}})\,d\boldsymbol{p}.$ (5)

Here $\boldsymbol{u}(\boldsymbol{x},t)$ and $q(\boldsymbol{x},t)$ are the
fluid velocity and pressure, $\mu$ is the fluid viscosity, and
${\boldsymbol{\Sigma}}^{a}(\boldsymbol{x},t)$ is the trace-free active stress
exerted by the particles on the fluid. The active stress is the orientational
average of the force dipoles exerted by the particles on the fluid, and the
sign of the coefficient $a_{0}$ corresponds to the sign of the dipoles:
$a_{0}>0$ for _puller_ particles, while $a_{0}<0$ for _pushers_. Note that
${\boldsymbol{\Sigma}}^{a}$ vanishes when the particles are in complete
nematic alignment.

###  2.2 Nondimensionalization and quantities of interest

We choose to nondimensionalize equations (1)–(5) over slightly different
characteristic velocity, length, and time scales from those commonly used in
the literature (Ezhilan et al., 2013; Hohenegger & Shelley, 2010; Saintillan &
Shelley, 2008b, a). In particular, letting $N$ denote the number of particles
in the system and $L$ denote the length of the periodic box in which the
particles are suspended, we nondimensionalize the model (1)–(5) according to

$\Psi^{\prime}=\frac{\Psi}{N/L^{2}},\quad\boldsymbol{x}^{\prime}=\frac{\boldsymbol{x}}{L/2\upi},\quad
t^{\prime}=\frac{|a_{0}|N}{L^{2}\mu}t,\quad\boldsymbol{u}^{\prime}=\frac{2\upi\mu
L}{|a_{0}|N}\boldsymbol{u},$

which results in the following system of equations:

$\displaystyle\p_{t^{\prime}}\Psi^{\prime}$
$\displaystyle=-{\rm{div}\,}^{\prime}(\p_{t^{\prime}}\boldsymbol{x}^{\prime}\Psi^{\prime})-\bnabla_{p}(\p_{t^{\prime}}\boldsymbol{p}\Psi^{\prime})$
(6) $\displaystyle\p_{t^{\prime}}\boldsymbol{x}^{\prime}$
$\displaystyle=\beta\boldsymbol{p}+\boldsymbol{u}^{\prime}-D_{T}^{\prime}\bnabla^{\prime}(\log\Psi^{\prime})$
$\displaystyle\p_{t^{\prime}}\boldsymbol{p}$ $\displaystyle=({\bf
I}-\boldsymbol{p}\boldsymbol{p}^{\rm
T})(\bnabla^{\prime}\boldsymbol{u}^{\prime}\,\boldsymbol{p})-D_{R}^{\prime}\bnabla_{p}(\log\Psi^{\prime})$
$\displaystyle-\Delta^{\prime}\boldsymbol{u}^{\prime}+\bnabla^{\prime}q^{\prime}$
$\displaystyle=\pm\int_{S^{1}}(\boldsymbol{p}\boldsymbol{p}^{\rm
T}-\textstyle\frac{1}{2}{\bf
I})\bnabla^{\prime}\Psi^{\prime}(\boldsymbol{x},\boldsymbol{p},t)\,d\boldsymbol{p},\quad{\rm{div}\,}^{\prime}\boldsymbol{u}^{\prime}=0.$

Here we note the presence of three dimensionless parameters: the diffusion
coefficients

$D_{T}^{\prime}=\frac{4\upi^{2}\mu D_{T}}{\left\lvert a_{0}\right\rvert
N}\quad\text{and}\quad D_{R}^{\prime}=\frac{\mu L^{2}D_{R}}{\left\lvert
a_{0}\right\rvert N},$

and a nondimensional ‘swimming speed’

$\beta=\frac{2\upi\mu V_{0}L}{\left\lvert a_{0}\right\rvert N}.$ (7)

We choose this nondimensionalization in order to easily incorporate the
immotile state ($V_{0}=0$) into the analysis. Without swimming, equations
(1)–(5) describe a suspension of _shakers_ (Ezhilan et al., 2013; Stenhammar
et al., 2017), particles which do not swim but still exert an active stress on
the surrounding fluid. We also fix the domain to be the 2-dimensional torus
$\mathbb{T}^{2}:=\mathbb{R}^{2}/(2\upi\mathbb{Z})^{2}$.

The parameter $\beta$ contains more information than just the particle
swimming speed: it is really the ratio of swimming speed to the active stress
magnitude and particle concentration. Note that $\beta$ may be related to the
more familiar nondimensionalization used in Hohenegger & Shelley (2010);
Saintillan & Shelley (2008b) via $\beta=(2\upi\ell)/(|\alpha|L\nu)$, where
$\ell$ is the length of typical swimmer, $\alpha=a_{0}/(V_{0}\mu\ell)$ is a
dimensionless signed active stress coefficient ($\alpha>0$ for pullers and
$\alpha<0$ for pushers), and $\nu=N\ell^{2}/L^{2}$ is the relative volume
concentration of swimmers.

In (6), the active stress coefficient in the Stokes equations is scaled to
unit magnitude but retains the sign of the force dipole exerted by the
particles on the fluid: $+1$ for puller particles and $-1$ for pusher
particles.

Dropping the prime notation, the system (6) may be written more succinctly as

$\displaystyle\p_{t}\Psi$
$\displaystyle=-\beta\boldsymbol{p}\bcdot\bnabla\Psi-\boldsymbol{u}\bcdot\bnabla\Psi-{\rm{div}}_{p}\big{(}({\bf
I}-\boldsymbol{p}\boldsymbol{p}^{\rm
T})(\bnabla\boldsymbol{u}\boldsymbol{p})\Psi\big{)}+D_{T}\Delta\Psi+D_{R}\Delta_{p}\Psi$
(8) $\displaystyle-\Delta\boldsymbol{u}+\bnabla q$
$\displaystyle=\pm\int_{S^{1}}(\boldsymbol{p}\boldsymbol{p}^{\rm
T}-\textstyle\frac{1}{2}{\bf
I})\bnabla\Psi(\boldsymbol{x},\boldsymbol{p},t)\,d\boldsymbol{p},\quad{\rm{div}\,}\boldsymbol{u}=0.$

One way to measure deviations of the swimmer density $\Psi$ from the uniform,
isotropic steady state $\Psi_{0}=1/(2\upi)$ is to consider the relative
entropy

$\mathcal{S}(t)=\int_{\mathbb{T}^{2}}\int_{S^{1}}\frac{\Psi}{\Psi_{0}}\log\bigg{(}\frac{\Psi}{\Psi_{0}}\bigg{)}\,d\boldsymbol{p}d\boldsymbol{x}.$
(9)

Using (8), the entropy can be shown to evolve according to

$\dot{\mathcal{S}}(t)=\mp\frac{4}{\Psi_{0}}\int_{\mathbb{T}^{2}}\left\lvert\mathsfbi{E}\right\rvert^{2}d\boldsymbol{x}-\int_{\mathbb{T}^{2}}\int_{S^{1}}\big{(}D_{T}\left\lvert\bnabla(\log\Psi)\right\rvert^{2}+D_{R}\left\lvert\bnabla_{p}(\log(\Psi))\right\rvert^{2}\big{)}\frac{\Psi}{\Psi_{0}}d\boldsymbol{p}d\boldsymbol{x}$
(10)

where
$\mathsfbi{E}=\frac{1}{2}(\bnabla\boldsymbol{u}+(\bnabla\boldsymbol{u})^{\rm
T})$ (see (Saintillan & Shelley, 2008b)). The first term in (10) arises from
the viscous dissipation of the active stress exerted by the particles and is
negative for pullers and positive for pushers. The two diffusive terms are
negative; hence we expect puller suspensions to always relax to isotropy over
time. For pushers, however, we may expect to see some more interesting
behaviors: indeed, simulations show that patterns and fluctuations in particle
alignment and concentration arise in certain parameter regimes (Hohenegger &
Shelley, 2010; Saintillan & Shelley, 2008b, a, 2013, 2015). We aim to study
the onset of pattern formation in these active pusher suspensions.

It will be useful to first define some system quantities that can be measured
numerically and used to verify analytical predictions. One such quantity is
the active power input, defined for pusher particles by

$\mathcal{P}(t)=\int_{\mathbb{T}^{2}}\int_{S^{1}}\boldsymbol{p}^{\rm
T}\mathsfbi{E}(\boldsymbol{x},t)\boldsymbol{p}\,\Psi(\boldsymbol{x},\boldsymbol{p},t)\,d\boldsymbol{p}d\boldsymbol{x}.$
(11)

The sign is opposite for puller particles. This quantity can be understood as
the perturbative power input due to the interaction of the active particles
with the flow (as opposed to the power input of each individual particle).
Using the Stokes equations in (8), the active power balances the rate of
viscous dissipation in the fluid:

$\mathcal{P}(t)=\int_{\mathbb{T}^{2}}2\left\lvert\mathsfbi{E}(\boldsymbol{x},t)\right\rvert^{2}\,d\boldsymbol{x}.$
(12)

Note that $\mathcal{P}(t)=0$ for the uniform, isotropic steady state, and
tracking the growth of $\mathcal{P}(t)$ (or, equivalently, the rate of viscous
dissipation) serves as a measure of the instability of the uniform state
(Saintillan & Shelley, 2015).

We also define the concentration field

$c(\boldsymbol{x},t)=\int_{S^{1}}\Psi(\boldsymbol{x},\boldsymbol{p},t)\,d\boldsymbol{p}$
(13)

and nematic order tensor

$\mathsfbi{Q}(\boldsymbol{x},t)=\frac{1}{c(\boldsymbol{x},t)}\int_{S^{1}}(\boldsymbol{p}\boldsymbol{p}^{\rm
T}-{\bf I}/2)\Psi(\boldsymbol{x},\boldsymbol{p},t)\,d\boldsymbol{p}.$ (14)

We can then define the scalar-valued nematic order parameter as

$\mathcal{N}(\boldsymbol{x},t)=\text{maximum eigenvalue of
}\mathsfbi{Q}(\boldsymbol{x},t).$ (15)

The nematic order parameter measures the local degree of nematic alignment:
$\mathcal{N}(\boldsymbol{x},t)=0$ when the particle orientations are isotropic
and $\mathcal{N}(\boldsymbol{x},t)=1$ when all surrounding particles are
exactly aligned.

###  2.3 Linear stability: the eigenvalue problem

From here, we restrict our attention to the pusher case in two spatial
dimensions. We begin by recalling the results of the eigenvalue problem
resulting from a linear stability analysis about the uniform, isotropic steady
state $\Psi_{0}=1/(2\upi)$ and $\boldsymbol{u}=0$ (Hohenegger & Shelley, 2010;
Saintillan & Shelley, 2008b; Subramanian et al., 2011). Linearizing (8) about
this state, we obtain

$\displaystyle\p_{t}\Psi$
$\displaystyle=-\beta\boldsymbol{p}\bcdot\bnabla\Psi+2\boldsymbol{p}^{\rm
T}\bnabla\boldsymbol{u}\boldsymbol{p}+D_{T}\Delta\Psi+D_{R}\Delta_{p}\Psi$
(16) $\displaystyle-\Delta\boldsymbol{u}+\bnabla q$
$\displaystyle=-\frac{1}{2\upi}\int_{S^{1}}(\boldsymbol{p}\boldsymbol{p}^{\rm
T}-\textstyle\frac{1}{2}{\bf
I})\bnabla\Psi(\boldsymbol{x},\boldsymbol{p},t)\,d\boldsymbol{p},\quad{\rm{div}\,}\boldsymbol{u}=0.$

We insert the plane wave ansatz $\Psi=\psi(\boldsymbol{k},\boldsymbol{p}){\rm
e}^{{\rm i}\boldsymbol{k}\bcdot\boldsymbol{x}+\sigma t}$,
$\sigma\in\mathbb{C}$, into (16), and choose coordinates such that the
wavevector $\boldsymbol{k}=k\boldsymbol{e}_{x}$ and the particle orientation
$\boldsymbol{p}=\cos\theta\boldsymbol{e}_{x}+\sin\theta\boldsymbol{e}_{y}$,
$0\leq\theta<2\upi$. Defining $\lambda_{k}:=\sigma+D_{T}k^{2}$, we obtain an
eigenvalue problem for $\lambda_{k}$ in the form of an integrodifferential
equation on $S^{1}$:

$\lambda_{k}\psi(k,\theta)=-{\rm
i}k\beta\cos\theta\psi(k,\theta)+\frac{1}{\upi}\cos\theta\sin\theta\int_{0}^{2\upi}\psi(k,\theta)\,\sin\theta\cos\theta\,d\theta+D_{R}\p_{\theta\theta}\psi(k,\theta).$
(17)

We note that while it may be more suggestive to write $\psi$ in terms of
$\sin(2\theta)$ rather than $\sin\theta\cos\theta$, we find that the details
of the weakly nonlinear calculations in the following sections are slightly
simpler in terms of $\cos\theta=\boldsymbol{p}\cdot\boldsymbol{e}_{x}$ and
$\sin\theta=\boldsymbol{p}\cdot\boldsymbol{e}_{y}$ only.

In the absence of rotational diffusion ($D_{R}=0$), we may solve (17) for
$\psi(k,\theta)$ as

$\displaystyle\psi(k,\theta)=\frac{\cos\theta\sin\theta}{\lambda_{k}+{\rm
i}k\beta\cos\theta},$

where $\lambda_{k}$ is such that

$F[\psi]:=\int_{0}^{2\upi}\psi(k,\theta)\,\cos\theta\sin\theta\,d\theta=\upi.$
(18)

In particular, $\lambda_{k}$ satisfies the implicit dispersion relation

$\frac{\lambda_{k}\beta^{2}k^{2}+2\lambda_{k}^{3}-2\lambda_{k}^{2}\sqrt{\beta^{2}k^{2}+\lambda_{k}^{2}}}{\beta^{4}k^{4}}=1.$
(19)

Recalling that $\lambda_{k}=\sigma+D_{T}k^{2}$, we may then numerically solve
for the relationship between $\sigma$ and $\beta$ (see figure 1).

(a)

(b)

Figure 1: (a) Real part (shifted by $D_{T}k^{2}$) and (b) imaginary part of
the growth rate $\sigma(k)$ versus $\beta\left\lvert k\right\rvert$. Figure
(a) can be read as follows: fix a value of $0<D_{T}k^{2}<0.25$. Look at the
corresponding horizontal green line. As $\beta$ is lowered, the green line
intersects the blue curve. This is where the eigenvalue $\sigma(k)$ becomes
unstable.
Different types of bifurcations are possible depending on the value of
$D_{T}k^{2}$: for example, point B corresponds to a purely real eigenvalue
crossing, while the presence of nonzero Im$(\sigma)$ signals a Hopf
bifurcation at points D and E. At point C, indicated by a red dot, two real
eigenvalues meet and become complex. In the case of shakers ($\beta=0$), we
may consider the real eigenvalue crossing at point A along the $y$-axis.

A similar calculation for $\boldsymbol{k}=k\boldsymbol{e}_{y}$ shows that the
eigenvalues $\sigma(k)$ plotted in figure 1 also correspond to a $y$-direction
eigenfunction whose eigenmodes are a 90-degree rotation (in $\theta$) of the
$x$-direction eigenmodes. The solutions of (16) arising from eigenfunctions of
the linearized operator are thus given by

$\Psi(\boldsymbol{x},\theta,t)=c_{x}\psi_{x}(k,\theta){\rm e}^{{\rm
i}kx+\sigma t}+c_{y}\psi_{y}(k,\theta){\rm e}^{{\rm i}ky+\sigma t},\quad
c_{x}^{2}+c_{y}^{2}=1,$ (20)

and all scalar multiples of (20), where

$\displaystyle\psi_{x}(k,\theta)=\frac{\cos\theta\sin\theta}{\sigma+D_{T}k^{2}+{\rm
i}k\beta\cos\theta},\quad\psi_{y}(k,\theta)=\frac{\cos\theta\sin\theta}{\sigma+D_{T}k^{2}+{\rm
i}k\beta\sin\theta}.$

In particular, besides point C, each eigenvalue $\sigma(k)$ has a
4-dimensional eigenspace over $\mathbb{C}$, spanned by the $x$\- and
$y$-direction components with $\pm k$. When $\beta\left\lvert
k\right\rvert<\sqrt{3}/9$, there are two distinct real-valued eigenvalues
$\sigma(k)$ which satisfy the required condition
$\int_{0}^{2\upi}\psi_{x}\,\cos\theta\sin\theta
d\theta=\int_{0}^{2\upi}\psi_{y}\,\cos\theta\sin\theta d\theta=\upi$, both of
which then have a 4-dimensional eigenspace over $\mathbb{R}$.

Note that the dispersion relation (19) is exact only in the absence of
rotational diffusion, but for the purposes of this paper, we will consider
$0<D_{R}\ll 1$. This alters figure 1, but only slightly (see next section 2.4
for details). In return, we do not have to contend with the continuous
spectrum present in the $D_{R}=0$ spectral analysis of Subramanian et al.
(2011). Furthermore, when $D_{R}>0$, the $k=0$ mode is always linearly stable,
as it satisfies the heat equation in $\boldsymbol{p}$: if
$\Psi=\Psi(\boldsymbol{p},t)$ then $\p_{t}\Psi=D_{R}\Delta_{p}\Psi$.

The eigenvalue relation (19) ceases to be valid for Re$(\lambda_{k})\leq 0$;
i.e. when Re$(\sigma(k))\leq-D_{T}k^{2}$ (Hohenegger & Shelley, 2010;
Subramanian et al., 2011). For fixed $0<D_{T}<1/4$, however, the eigenvalue
analysis does capture a sign change in the real part of the growth rate
$\sigma(k)$ as $\beta\left\lvert k\right\rvert$ is varied. For each $k$, this
sign change occurs where the blue curve in figure 1 (corresponding to
$\lambda_{k}$) intersects the green line corresponding to the fixed value of
$D_{T}k^{2}$. Numerical simulations indicate that if the point
$(D_{T}k^{2},\beta\left\lvert k\right\rvert)$ lies to the right of the blue
contour in figure 1(a) for each $k$, then the uniform, isotropic steady state
is stable to small perturbations (see sections 4 and 5). In this region,
particle diffusion and swimming – processes which tend to decrease order among
the particles – are large relative to both the active stress magnitude and the
particle concentration, which favor local alignment. As $\beta$ is decreased,
Re($\sigma(k)$) crosses the green line corresponding to the (fixed) value of
$D_{T}k^{2}$ and the uniform, isotropic state becomes unstable.

Since we are considering the system (6) on a periodic domain and have
nondimensionalized $\boldsymbol{x}$ according to the length of the domain, the
$\left\lvert k\right\rvert=1$ mode is the lowest nontrivial mode in the system
and hence the first to lose stability as $\beta$ decreases. As noted in Koch &
Subramanian (2011), “the finiteness of the domain may act to stabilize the
system.” Indeed, the stabilizing effects of mixing on $\mathbb{T}^{d}$
($d=2,3$) due to swimming are studied in Albritton & Ohm (2022) and play a
role in the patterns observed here.

We aim to understand the onset of pattern formation in (6) by exploring the
many different ways the $\left\lvert k\right\rvert=1$ mode can lose stability
as $\beta$ is decreased. From figure 1, we can see that depending on $D_{T}$,
the type of bifurcation that we expect to see for the $\left\lvert
k\right\rvert=1$ mode will change. We aim to characterize all of the types of
initial bifurcations admitted by the system through a weakly nonlinear
analysis. According to the dispersion relation, if $1/9<D_{T}<1/4$, as $\beta$
decreases, the purely real eigenvalue $\sigma$ will change sign from positive
to negative across the top branch of the blue curve, between points A and C.
For $0<D_{T}<1/9$, however, we expect to see a Hopf bifurcation as $\beta$
crosses the blue curve roughly between points C and E, since for these values
of $\beta$, $\sigma$ has nonzero imaginary part.

As noted, figure 1 does not quite give the exact locations of the bifurcations
that we will consider here, since we still need to consider the effects of
(small) $D_{R}>0$. This will be the subject of the next section.

###  2.4 Role of rotational diffusion

The dispersion relation (19) was obtained in the absence of rotational
diffusion; however, studying pattern formation near the isotropic steady state
will require $D_{R}>0$. When $D_{R}=0$, the system (8) has infinitely many
steady states; in particular, any spatially uniform swimmer distribution
$\Psi=\Psi(\boldsymbol{p})$ is a steady state. The continuum of nearby steady
states obscures the mechanism by which the isotropic state $\Psi=1/(2\upi)$
loses stability; indeed, any function $\Psi(\boldsymbol{p})$ belongs to the
kernel of the linearized operator. When we add in $D_{R}>0$ (along with the
assumption that the total number of swimmers is conserved), this kernel is
eliminated. Thus we aim to determine when the effect of small $D_{R}>0$ can be
considered as a perturbation of the dispersion relation (19) and figure 1.

In particular, given a wavenumber $k$, for small $D_{R}>0$, we wish to
determine when an expansion (in $D_{R}$) of the form

$\sigma=\sigma_{0}+D_{R}\sigma_{1}+O(D_{R}^{2}),\quad\psi=\psi_{0}+D_{R}\psi_{1}+O(D_{R}^{2})$
(21)

is valid for some $(\sigma_{1},\psi_{1})$.

Plugging this expansion into the eigenvalue problem (17) and separating
scales, at $O(1)$ we obtain the $D_{R}=0$ eigenfunctions and eigenvalue
relation

$\displaystyle\psi_{0}(k,\theta)=\frac{\cos\theta\sin\theta}{\sigma_{0}+D_{T}k^{2}+{\rm
i}k\beta\cos\theta}$

where $\sigma_{0}$ is such that

$F[\psi_{0}]=\int_{0}^{2\upi}\psi_{0}\cos\theta\sin\theta\,d\theta=\upi.$ (22)

At $O(D_{R})$ we obtain the expression

$\displaystyle\psi_{1}$
$\displaystyle=\frac{1}{\upi}\psi_{0}\int_{0}^{2\upi}\psi_{1}\cos\theta^{\prime}\sin\theta^{\prime}\,d\theta^{\prime}-\frac{\sigma_{1}\psi_{0}}{\sigma_{0}+D_{T}k^{2}+{\rm
i}k\beta\cos\theta}+\frac{\p_{\theta}^{2}\psi_{0}}{\sigma_{0}+D_{T}k^{2}+{\rm
i}k\beta\cos\theta}.$

Taking $F[\bcdot]$ of both sides and using (22) then yields an integral
expression for $\sigma_{1}$:

$\displaystyle\sigma_{1}\int_{0}^{2\upi}\frac{\cos^{2}\theta\sin^{2}\theta}{(\sigma_{0}+D_{T}k^{2}+{\rm
i}k\beta\cos\theta)^{2}}d\theta$
$\displaystyle=\int_{0}^{2\upi}\frac{\p_{\theta}^{2}\psi_{0}\cos\theta\sin\theta}{\sigma_{0}+D_{T}k^{2}+{\rm
i}k\beta\cos\theta}d\theta.$

As long as Re$(\sigma_{0}+D_{T}k^{2})\neq 0$ (which holds for the $\left\lvert
k\right\rvert=1$ modes as long as $D_{T}\neq 0$), we obtain an expression for
$\sigma_{1}$:

$\displaystyle\sigma_{1}$
$\displaystyle=\frac{\frac{5}{6}\beta^{2}k^{2}}{\beta^{2}k^{2}-3(D_{T}k^{2}+\sigma_{0})^{2}}+\frac{\frac{1}{2}\beta^{2}k^{2}}{\beta^{2}k^{2}+(D_{T}k^{2}+\sigma_{0})^{2}}$
$\displaystyle\qquad+\frac{5(D_{T}k^{2}+\sigma_{0})^{3}}{(\beta^{2}k^{2}-3(D_{T}k^{2}+\sigma_{0})^{2})\sqrt{\beta^{2}k^{2}+(D_{T}k^{2}+\sigma_{0})^{2}}}-\frac{7}{3},$

which is finite as long as $\beta^{2}k^{2}\neq
3(D_{T}k^{2}+\text{Re}(\sigma_{0}))^{2}$. The line

$\sigma_{0}+D_{T}k^{2}=\frac{\beta\left\lvert k\right\rvert}{\sqrt{3}}$ (23)

is plotted along with the real part of the $D_{R}=0$ dispersion relation (19)
in figure 2(a). Away from this line, for small $D_{R}>0$, we may consider the
solutions of the eigenvalue problem (17) as perturbations of the $D_{R}=0$
eigenvalues and eigenfunctions. As we can see, this line (23) passes through
the point C from figure 1 where the two real eigenvalues meet and become two
complex conjugate eigenvalues. A very precise choice of $D_{T}$ and $D_{R}$
should correspond to a codimension 2 bifurcation at point C; however, a
different scaling than (21) with respect to $D_{R}$ is likely needed to study
this point in detail. We will not attempt to study point C in detail here, and
will instead focus on the more generic bifurcations between points A and C and
between points C and E.

(a)

(b)

Figure 2: The (a) real part and (b) imaginary part of the perturbed dispersion
relation $\sigma_{0}+D_{R}\sigma_{1}$ (dotted black curve) is plotted for
$D_{R}=0.001$ on top of the unperturbed relation with $D_{R}=0$ (solid light
blue curve). The perturbative expression fails to be valid along the gray line
$\sigma_{0}+D_{T}k^{2}=\beta\left\lvert k\right\rvert/\sqrt{3}$ plotted in
(a). The inset in figure (a) shows in greater detail the behavior of the
perturbative expression $\sigma_{0}+D_{R}\sigma_{1}$ near the point
$(\sqrt{3}/9,1/9)$ where the gray line intersects the unperturbed expression.
In particular, the perturbed expression blows up as $\beta\left\lvert
k\right\rvert$ approaches $\sqrt{3}/9$.

In figures 2(a) and 2(b) we plot the dispersion relation for
$\sigma_{0}+D_{R}\sigma_{1}$ on top of $\sigma_{0}$ using $D_{R}=0.001$. Away
from the crossing with the line (23), the $\sigma_{0}+D_{R}\sigma_{1}$ curve
aligns very closely with the $\sigma_{0}$ curve. Hereafter, for most numerical
purposes, we will take $D_{R}=0.001$ – small enough to use figure 1 as our
roadmap for determining roughly where in the $(D_{T},\beta)$ parameter space
to look to see different system behaviors. We note that, while the expansion
in $D_{R}$ is not rigorous, it is backed later on by the close agreement of
the predicted behavior near bifurcation points (from amplitude equations) with
the observed behavior from numerical simulations. In particular, it appears
that, away from point ${\bf C}$, the small $D_{R}>0$ picture is largely
captured by this expansion.

###  2.5 Outline of results

The remainder of the paper is devoted to a weakly nonlinear analysis of the
different possible bifurcations apparent in figure 1, which we map out in
greater detail in figure 3. We begin in section 3 by considering the immotile
case $\beta=0$. We examine the real eigenvalue crossing at point A for all
values of $D_{R}$ for which a bifurcation occurs and show that the resulting
pitchfork bifurcation is always supercritical. In section 4 we assume $D_{R}$
is very small and analyze the Hopf bifurcation occuring along the curve
between points C and ${\bf E}$. We show that for initial perturbations to the
uniform isotropic state with both $x$ and $y$ components (i.e. both $c_{x}$,
$c_{y}\neq 0$ in (20); see line labeled 2D in figure 3), the bifurcation
transitions from supercritical to subcritical at roughly point D, but is
always supercritical for initial perturbations with either $c_{y}=0$ or
$c_{x}=0$ (labeled 1D in figure 3). In section 5, we again consider very small
$D_{R}$ and study real eigenvalue crossing occuring along the curve between
points A and C. The pitchfork bifurcation also transitions from supercritical
to subcritical in both the 2D and 1D cases, with the change occuring roughly
at point B in the case of an initial perturbation in both $x$ and $y$, and
just before point C in the case of an $x$-only or $y$-only initial
perturbation.

Figure 3: Diagram of the various types of bifurcations through which the
uniform isotropic steady state loses stability, depending on the location of
the bifurcation value $\beta_{T}$. Here the subscript $T$ is used to reflect
that the value of $\beta_{T}$ depends on the translational diffusion $D_{T}$
through the dispersion relation plotted in figure 1. Note that the letters A –
E correspond to the positions in figure 1(a), which we also repeat here for
clarity. The upper line labeled 2D corresponds to the evolution of initial
perturbations to the uniform isotropic state with both components in both the
$x$ and $y$ direction ($c_{x}$, $c_{y}\neq 0$ in (20)), while the lower line
labeled 1D corresponds to perturbations with $x$ or $y$ component only
($c_{y}=0$ or $c_{x}=0$ in (20)).

Each section is accompanied by numerical simulations verifying the predicted
behaviors near the different bifurcations and illustrating the various states
that arise. The numerics are a pseudo-spectral implementation of equations (6)
with time-stepping via a second order implicit-explicit backward
differentiation scheme.

##  3 Immotile particles: supercritical pitchfork bifurcation

The simplest scenario for studying the onset of pattern formation in the model
(8) is in the case $\beta=0$; i.e. the particles are immotile (or _shakers_
(Ezhilan et al., 2013; Stenhammar et al., 2017)) but still exert a dipolar
force on the surrounding fluid. The uniform, isotropic steady state in a
suspension of immotile particles undergoes a bifurcation to an alignment
instability, indicated by point A in figure 1, which we study in greater
detail. We first show via weakly nonlinear analysis that the resulting
pitchfork bifurcation is always supercritical, and we then numerically explore
examples of the emerging nontrivial steady state.

###  3.1 Weakly nonlinear analysis

In the case of immotile particles, we can explicitly calculate the eigenvalues
and eigenfunctions of the linearized operator when $D_{R}>0$. The (purely
real) eigenvalues $\sigma(k)$ and eigenmodes $\psi_{x}(k,\theta)$,
$\psi_{y}(k,\theta)$ are given by

$\sigma(k)=\frac{1}{4}-D_{T}k^{2}-4D_{R},\quad\psi_{x}(k,\theta)=\psi_{y}(k,\theta)=\cos\theta\sin\theta.$
(24)

From the immotile dispersion relation (24), when $D_{T}>\frac{1}{4}-4D_{R}$,
all eigenvalues $\sigma(k)$ are negative, but as $D_{T}$ is decreased, the
$\left\lvert k\right\rvert=1$ modes are the first to change sign as
$D_{T}=D_{T}^{*}=\frac{1}{4}-4D_{R}$ is crossed. We study the nature of this
bifurcation for different $D_{R}$ via the method of multiple scales. For
$0<\epsilon\ll 1$, we fix $D_{T}=\frac{1}{4}-4D_{R}-\epsilon^{2}$, so the
$\left\lvert k\right\rvert=1$ modes are just barely growing, and define the
slow timescale $\tau=\epsilon^{2}t$. We then assume the following expansions
in $\epsilon$:

$\Psi=\frac{1}{2\upi}(1+\epsilon\Psi_{1}+\epsilon^{2}\Psi_{2}+\epsilon^{3}\Psi_{3}+\cdots),\qquad\boldsymbol{u}=\epsilon\boldsymbol{u}_{1}+\epsilon^{2}\boldsymbol{u}_{2}+\epsilon^{3}\boldsymbol{u}_{3}+\cdots.$
(25)

Plugging each of these expansions into equation (8) with $\beta=0$ and
separating by orders of $\epsilon$, at $O(\epsilon)$ we obtain the $\beta=0$
version of the eigenvalue equation (17), evaluated at the critical value
$D_{T}^{*}=\frac{1}{4}-4D_{R}$ where $\sigma(1)=0$:

$\mathcal{L}[\Psi_{1}]:=-2\boldsymbol{p}^{\rm
T}\bnabla\boldsymbol{u}_{1}\boldsymbol{p}-\bigg{(}\frac{1}{4}-4D_{R}\bigg{)}\Delta\Psi_{1}-D_{R}\Delta_{p}\Psi_{1}=0.$
(26)

This equation is satisfied by the $\left\lvert k\right\rvert=1$ modes of the
eigenfunctions (20), where we recall that the eigenmodes in the immotile case
are given by (24). Thus $\Psi_{1}$ and $\boldsymbol{u}_{1}$ have the form

$\displaystyle\Psi_{1}$
$\displaystyle=\cos\theta\sin\theta\big{(}c_{x}A_{x}(\tau){\rm e}^{{\rm
i}x}+c_{y}A_{y}(\tau){\rm e}^{{\rm i}y}\big{)}+{\rm c.c.},$ (27)
$\displaystyle\boldsymbol{u}_{1}$ $\displaystyle=-\frac{\rm
i}{8}\big{(}c_{x}A_{x}(\tau){\rm e}^{{\rm
i}x}\boldsymbol{e}_{y}+c_{y}A_{y}(\tau){\rm e}^{{\rm
i}y}\boldsymbol{e}_{x}\big{)}+{\rm c.c.}$

Here we have inserted the complex amplitudes $A_{x}(\tau)$, $A_{y}(\tau)$
which depend solely on the slow timescale $\tau$, and for which we aim to find
an equation. Throughout, we use ${\rm c.c.}$ to denote the complex conjugate
of each of the preceding terms.

At $O(\epsilon^{2})$ we obtain the equation

$\mathcal{L}[\Psi_{2}]=-\boldsymbol{u}_{1}\bcdot\bnabla\Psi_{1}-{\rm{div}}_{p}\big{(}({\bf
I}-\boldsymbol{p}\boldsymbol{p}^{\rm
T})(\bnabla\boldsymbol{u}_{1}\boldsymbol{p})\Psi_{1})\big{)}$ (28)

where the operator $\mathcal{L}$ is as defined in (26). Using the expressions
(27) for $\Psi_{1}$ and $\boldsymbol{u}_{1}$, the right hand side of equation
(28) can be calculated explicitly (see equation (70) in Appendix A). Noting
that the right hand side expression contains only exponential terms of the
form ${\rm e}^{\pm 2{\rm i}x}$, ${\rm e}^{\pm 2{\rm i}y}$, and ${\rm
e}^{\pm{\rm i}(x\pm y)}$, with no terms proportional to ${\rm e}^{\pm{\rm
i}x}$ or ${\rm e}^{\pm{\rm i}y}$, equation (28) is solvable without any
additional conditions on the coefficients $A_{x}$ and $A_{y}$. In particular,
due to the form of the right hand side expression, we look for $\Psi_{2}$ and
corresponding $\boldsymbol{u}_{2}$ of the form

$\displaystyle\Psi_{2}$ $\displaystyle=\psi_{2,1}{\rm e}^{{\rm
i}(x+y)}A_{x}A_{y}+\psi_{2,2}{\rm e}^{{\rm
i}(x-y)}A_{x}\overline{A}_{y}+\psi_{2,3}A_{x}^{2}{\rm e}^{2{\rm
i}x}+\psi_{2,4}A_{y}^{2}{\rm e}^{2{\rm i}y}$ (29) $\displaystyle\hskip
85.35826pt+\psi_{2,5}\left\lvert A_{x}\right\rvert^{2}+\psi_{2,6}\left\lvert
A_{y}\right\rvert^{2}+{\rm c.c.},$ $\displaystyle\boldsymbol{u}_{2}$
$\displaystyle=-\frac{\rm i}{8\upi}\int_{0}^{2\upi}\bigg{(}\psi_{2,1}{\rm
e}^{{\rm
i}(x+y)}A_{x}A_{y}(\boldsymbol{e}_{x}-\boldsymbol{e}_{y})+\psi_{2,2}{\rm
e}^{{\rm
i}(x-y)}A_{x}\overline{A}_{y}(\boldsymbol{e}_{x}+\boldsymbol{e}_{y})\bigg{)}(\cos^{2}\theta-\sin^{2}\theta)\,d\theta$
$\displaystyle\qquad-\frac{\rm i}{4\upi}\int_{0}^{2\upi}\bigg{(}\psi_{2,3}{\rm
e}^{2{\rm i}x}A_{x}^{2}\boldsymbol{e}_{y}+\psi_{2,4}{\rm e}^{2{\rm
i}y}A_{y}^{2}\boldsymbol{e}_{x}\bigg{)}\sin\theta\cos\theta\,d\theta+{\rm
c.c.},$

where each $\psi_{2,j}=\psi_{2,j}(\theta)$. Plugging these expressions (29)
into the left hand side of (28) (see Appendix A, equation (71) for the full
expression), after matching exponents with the right hand side, we can
explicitly solve for each $\psi_{2,j}$:

$\displaystyle\psi_{2,1}$
$\displaystyle=-\frac{c_{x}c_{y}}{1+16D_{R}}(\sin^{4}\theta+\cos^{4}\theta)-\frac{c_{x}c_{y}}{2(1-8D_{R})}\sin\theta\cos\theta+\frac{3c_{x}c_{y}}{4(1+16D_{R})}$
(30) $\displaystyle\psi_{2,2}$
$\displaystyle=-\frac{c_{x}c_{y}}{1+16D_{R}}(\sin^{4}\theta+\cos^{4}\theta)+\frac{c_{x}c_{y}}{2(1-8D_{R})}\sin\theta\cos\theta+\frac{3c_{x}c_{y}}{4(1+16D_{R})}$
$\displaystyle\psi_{2,3}$
$\displaystyle=\frac{c_{x}^{2}}{4}\bigg{(}-2\cos^{4}\theta+\frac{3(1-16D_{R})}{2(1-12D_{R})}\cos^{2}\theta+\frac{3D_{R}}{1-12D_{R}}\bigg{)}$
$\displaystyle\psi_{2,4}$
$\displaystyle=\frac{c_{y}^{2}}{4}\bigg{(}-2\sin^{4}\theta+\frac{3(1-16D_{R})}{2(1-12D_{R})}\sin^{2}\theta+\frac{3D_{R}}{1-12D_{R}}\bigg{)}$
$\displaystyle\psi_{2,5}$
$\displaystyle=-\frac{c_{x}^{2}}{16D_{R}}\cos^{4}\theta+\frac{3c_{x}^{2}}{128D_{R}}$
$\displaystyle\psi_{2,6}$
$\displaystyle=-\frac{c_{y}^{2}}{16D_{R}}\sin^{4}\theta+\frac{3c_{y}^{2}}{128D_{R}}.$

Inserting each of the coefficients (30) in the expression (29) for
$\boldsymbol{u}_{2}$, we have that each $\theta$-integral vanishes and
therefore $\boldsymbol{u}_{2}=0$. At $O(\epsilon^{3})$ we thus obtain the
following equation for $\Psi_{3}$:

$\displaystyle\mathcal{L}[\Psi_{3}]$
$\displaystyle=-\p_{\tau}\Psi_{1}-\boldsymbol{u}_{1}\bcdot\bnabla\Psi_{2}-{\rm{div}}_{p}\big{(}({\bf
I}-\boldsymbol{p}\boldsymbol{p}^{\rm
T})(\bnabla\boldsymbol{u}_{1}\boldsymbol{p})\Psi_{2}\big{)}-\Delta\Psi_{1}.$
(31)

Letting $\mathcal{R}(\boldsymbol{x},\theta,\tau)$ denote the right hand side
of (31), we have that $\mathcal{R}$ may be calculated explicitly using (27)
and (29); in particular, $\mathcal{R}$ is of the form

$\displaystyle\mathcal{R}(\boldsymbol{x},\theta,\tau)$
$\displaystyle=R_{x}(\theta,\tau){\rm e}^{{\rm i}x}+R_{y}(\theta,\tau){\rm
e}^{{\rm i}y}+R_{2x^{+}}(\theta,\tau){\rm e}^{{\rm
i}(2x+y)}+R_{2y^{+}}(\theta,\tau){\rm e}^{{\rm i}(x+2y)}$ (32)
$\displaystyle\quad+R_{2x^{-}}(\theta,\tau){\rm e}^{{\rm
i}(2x-y)}+R_{2y^{-}}(\theta,\tau){\rm e}^{{\rm
i}(-x+2y)}+R_{3x}(\theta,\tau){\rm e}^{{\rm i}3x}+R_{3y}(\theta,\tau){\rm
e}^{{\rm i}3y}+{\rm c.c.}$

By the Fredholm alternative, equation (31) admits a solution $\Psi_{3}$ as
long as

$\int_{0}^{2\upi}\int_{\mathbb{T}^{2}}\mathcal{R}(\boldsymbol{x},\theta,\tau)\Phi(\boldsymbol{x},\theta)d\boldsymbol{x}d\theta=0\qquad\text{for
all }\Phi(\boldsymbol{x},\theta)\text{ such that }\mathcal{L}^{*}[\Phi]=0.$
(33)

Since the operator $\mathcal{L}$ defined in (26) is self-adjoint in the
immotile case, we have that any such $\Phi$ has the form
$\Phi=\cos\theta\sin\theta(\alpha_{x}{\rm e}^{{\rm i}x}+\alpha_{y}{\rm
e}^{{\rm i}y})+{\rm c.c.}$ for any $\alpha_{x}^{2}+\alpha_{y}^{2}=1$.

Thus (33) is automatically satisfied for each term of $\mathcal{R}$ except for
$R_{x}(\theta,\tau){\rm e}^{{\rm i}x}+R_{y}(\theta,\tau){\rm e}^{{\rm
i}y}+{\rm c.c}$. The exact form of $R_{x}$ and $R_{y}$ is given in Appendix A,
equation (72).

Since the ratio $\alpha_{x}/\alpha_{y}$ is arbitrary, we need that both

$\displaystyle\int_{0}^{2\upi}R_{x}(\theta,\tau)d\theta=0\quad\text{and}\quad\int_{0}^{2\upi}R_{y}(\theta,\tau)d\theta=0.$

These two conditions together lead to a coupled system of ODEs for the
amplitudes $A_{x}$, $A_{y}$:

$\displaystyle c_{x}\big{(}\quad\p_{\tau}A_{x}$
$\displaystyle=A_{x}+(c_{x}^{2}M_{1}\left\lvert
A_{x}\right\rvert^{2}+c_{y}^{2}M_{2}\left\lvert
A_{y}\right\rvert^{2})A_{x}\quad\big{)}$ (34) $\displaystyle
c_{y}\big{(}\quad\p_{\tau}A_{y}$
$\displaystyle=A_{y}+(c_{y}^{2}M_{1}\left\lvert
A_{y}\right\rvert^{2}+c_{x}^{2}M_{2}\left\lvert
A_{x}\right\rvert^{2})A_{y}\quad)$

where

$M_{1}=-\frac{3(3-28D_{R}-32D_{R}^{2})}{1024D_{R}(1-12D_{R})},\quad
M_{2}=\frac{7-136D_{R}-2432D_{R}^{2}}{1024D_{R}(1-8D_{R})(1+16D_{R})}.$ (35)

If $M_{1}+M_{2}<0$, the system (34) has real, nonzero steady states of the
form

$\displaystyle A_{x}$
$\displaystyle=\pm\frac{1}{c_{x}\sqrt{-(M_{1}+M_{2})}},\quad
A_{y}=\pm\frac{1}{c_{y}\sqrt{-(M_{1}+M_{2})}}.$

We have that

$M_{1}+M_{2}=\frac{-1-104D_{R}+560D_{R}^{2}+9600D_{R}^{3}-6144D_{R}^{4}}{512D_{R}(1-8D_{R})(1-12D_{R})(1+16D_{R})},$
(36)

which, as we can see from figure 4, is indeed negative for all values of
$D_{R}$ for which a bifurcation occurs. Thus for any relevant level of
rotational diffusion, the uniform, isotropic steady state loses stability
through a supercritical pitchfork bifurcation and nontrivial stable steady
states emerge.

Figure 4: In the immotile setting, the coefficients $M_{1}+M_{2}$ (36) and
$M_{1}$ (35) are both negative for all values of $D_{R}$ for which a
bifurcation occurs ($0<D_{R}<1/16$), indicating that a supercritical pitchfork
bifurcation occurs for both 2D ($x$ and $y$) and 1D ($x$-only) initial
perturbations to the uniform isotropic state.

To leading order in
$\epsilon=\sqrt{D_{T}^{*}-D_{T}}=\sqrt{\frac{1}{4}-4D_{R}-D_{T}}$, the stable
steady states which bifurcate from the uniform isotropic state are of the form

$\Psi=\frac{1}{2\upi}\pm\frac{\epsilon}{2\upi}\sqrt{\frac{512D_{R}(1-8D_{R})(1-12D_{R})(1+16D_{R})}{1+104D_{R}-560D_{R}^{2}-9600D_{R}^{3}+6144D_{R}^{4}}}\cos\theta\sin\theta\big{(}{\rm
e}^{{\rm i}x}\pm{\rm e}^{{\rm i}y}\big{)}+{\rm c.c.}$ (37)

Note that as long as the initial perturbation coefficients $c_{x}$ and $c_{y}$
are both nonzero, the form of (37) does not depend on $c_{x}$ or $c_{y}$. If
either $c_{x}=0$ or $c_{y}=0$ in (20), the bifurcating stable steady states
take on a different form. Without loss of generality, we consider $c_{y}=0$.
In this case, the coupled system (34) reduces to the single amplitude equation

$\p_{\tau}A_{x}=A_{x}+M_{1}\left\lvert A_{x}\right\rvert^{2}A_{x}$ (38)

where $M_{1}$ is as in (35). As shown in figure 4, $M_{1}$ is negative for all
$0<D_{R}<1/16$, and therefore the uniform steady state still loses stability
through a supercritical pitchfork bifurcation for all meaningful choices of
$D_{R}$. To leading order in $\epsilon$, the stable steady states which emerge
are of the form

$\Psi(\boldsymbol{x},\theta)=\frac{1}{2\upi}\bigg{(}1\pm\epsilon\sqrt{\frac{1024D_{R}(1-12D_{R})}{3(3-28D_{R}-32D_{R}^{2})}}\cos\theta\sin\theta{\rm
e}^{{\rm i}x}\bigg{)}+{\rm c.c.}$ (39)

Numerical evidence of supercriticality along with simulated examples of the
emerging steady states (37) and (39) are presented in the following section.

###  3.2 Numerics

To study the immotile bifurcation numerically, we begin by checking for
supercriticality. We first fix $D_{R}=0.0125$, so that by (24),
$D_{T}^{*}=0.2$ is the bifurcation value. Taking our initial condition to be a
random, small-magnitude perturbation to the uniform isotropic state in both
$x$ and $y$, we begin by running the simulation with $D_{T}=0.1$ until
$t=500$. Then the value of $D_{T}$ is increased by $0.02$ every $100t$ until
$D_{T}=0.3$. The bifurcation value $D_{T}^{*}=0.2$ is reached at $t=1000$.

Over the course of the simulation we keep track of the $L^{2}$ norm of the
velocity field

$\displaystyle\left\lVert\boldsymbol{u}(\bcdot,t)\right\rVert_{L^{2}(\mathbb{T}^{2})}^{2}=\int_{\mathbb{T}^{2}}\left\lvert\boldsymbol{u}(\boldsymbol{x},t)\right\rvert^{2}\,d\boldsymbol{x},$

which is plotted continuously over time in figure 5(a). We also keep track of
the time-averaged rate of viscous dissipation in the fluid for each value of
$D_{T}$, which we recall from (12) balances the active power input
$\mathcal{P}$. Given a constant value of $D_{T}$ over the time interval
$(t_{1},t_{2})$, we measure the value of

$\overline{\mathcal{P}}=\frac{1}{t_{2}-t_{1}}\int_{t_{1}}^{t_{2}}\mathcal{P}(t)\,dt=\frac{1}{t_{2}-t_{1}}\int_{t_{1}}^{t_{2}}\int_{\mathbb{T}^{2}}2\left\lvert\mathsfbi{E}(\boldsymbol{x},t)\right\rvert^{2}\,d\boldsymbol{x}dt,$
(40)

which we consider as a function of $D_{T}$. In our case, $t_{2}-t_{1}=100$ for
each different value of $D_{T}$. We plot $\overline{\mathcal{P}}$ in figure
5(b) over the course of the simulation for the various values of $D_{T}$. As
expected for a supercritical bifurcation, we see that $\overline{\mathcal{P}}$
smoothly decays to zero as $D_{T}$ increases toward the bifurcation value
$D_{T}^{*}=0.02$, and remains zero after the bifurcation is reached. This
smooth transition from nontrivial dynamics to the uniform, isotropic steady
state as $D_{T}$ is slowly varied from a very unstable value through the
bifurcation value and beyond may be contrasted with the hysteresis seen later
in the subcritical Hopf region for motile particles (Section 4.2).

(a)

(b)

Figure 5: Numerical evidence of supercriticality in the pitchfork bifurcation
for immotile particles ($\beta=0$). Here $D_{R}=0.0125$ is fixed and the
bifurcation occurs at $D_{T}^{*}=0.2$. The simulation is initiated with
$D_{T}=0.1$ until $t=500$; then every $100t$ the value of $D_{T}$ is increased
by $0.02$, so the bifurcation value is reached at $t=1000$. Figure (a) shows
the $L^{2}$ norm of the fluid velocity field over time. Figure (b) shows the
time-averaged viscous dissipation $\overline{\mathcal{P}}$ (40) for the
different values of $D_{T}$. The apparently smooth decrease to zero as
$D_{T}\to D_{T}^{*}$ indicates that the uniform isotropic steady state loses
stability through a supercritical pitchfork bifurcation.

We next consider what the emerging stable steady states actually look like.
Given $D_{R}$, we choose $D_{T}$ such that
$\epsilon^{2}=D_{T}^{*}-D_{T}=0.02$. We initialize the simulation by
perturbing the uniform, isotropic state with a small random-magnitude
perturbation to a random assortment of the five lowest spatial modes in both
$x$ and $y$ with random orientation $\theta$, and run until the dynamics
settle into a steady state. The resulting nematic order parameter and
direction of preferred local nematic alignment are plotted in figures 6(a) and
6(b) for $(D_{R},D_{T})=(0.0025,0.22)$ ($D_{T}^{*}=0.24$) and
$(D_{R},D_{T})=(0.03125,0.105)$ ($D_{T}^{*}=0.125$), respectively. We see that
the higher spatial modes decay over time, and the resulting steady state
consists only of the unstable $\left\lvert k\right\rvert=1$ mode. For
$(D_{R},D_{T})=(0.03125,0.105)$ we also plot the fluid vorticity field

$\omega(\boldsymbol{x},t)=\bnabla\times\boldsymbol{u}(\boldsymbol{x},t)$ (41)

and a sampling of the velocity field $\boldsymbol{u}(\boldsymbol{x},t)$
throughout the domain.

The supplementary video Movie1 shows the system approaching an example of the
steady state shown in figure 6 as the bifurcation is approached from below. In
this case $D_{R}=0.0125$ so $D_{T}^{*}=0.2$. The steady state is reached after
approaching $D_{T}=0.18$ from below.

(a)

(b)

(c)

Figure 6:  Plots of the nematic order parameter
$\mathcal{N}(\boldsymbol{x},t)$ and the direction of local nematic alignment
for (a) $(D_{R},D_{T})=(0.0025,0.22)$ ($D_{T}^{*}=0.24$) and (b)
$(D_{R},D_{T})=(0.03125,0.105)$ ($D_{T}^{*}=0.125$) demonstrate the dependence
of the emerging immotile steady state on $D_{R}$, as predicted by the form of
(37). In both cases, $\epsilon^{2}=D_{T}^{*}-D_{T}=0.02$. Figure (c) shows the
vorticity field $\omega(\boldsymbol{x},t)$ (colors) and velocity field
$\boldsymbol{u}(\boldsymbol{x},t)$ (arrows) for the same steady state pictured
in figure (b). Note the clear extensile flow produced by the aligned dipoles
in the top right and bottom left of the domain.

We also consider the same parameter combinations as in figure 6 for an
$x$-only initial perturbation, and plot the results in figure 7. Given the
form of the calculated steady state (37) for a 2D ($x$ and $y$) perturbation
and (39) for a 1D ($x$-only) perturbation, we expect the particles to display
a slight preference for the (nematic) orientations
$\theta=\frac{\upi}{4}(\equiv\frac{5\upi}{4})$ and
$\theta=\frac{3\upi}{4}(\equiv\frac{7\upi}{4})$ where
$\pm\cos\theta\sin\theta$ is maximized. This preference is clearly visible in
figures 6(a) and 6(b) as well as in figures 7(b) and 7(a).

In addition, in the 2D case, figure 4 suggests that given $\epsilon$ there
should be a $D_{R}\in(0,\frac{1}{16})$ which maximizes the magnitude of the
deviation from isotropy in the emerging stable steady state (37), and that as
$D_{R}\to 0$ the correlation in particle alignment should disappear. We can
see numerical evidence of this $D_{R}$ dependence in the difference in the
magnitude of $\mathcal{N}(\boldsymbol{x},t)$ between figures 6(b) and 6(a).

Similarly, in the case $c_{y}=0$, the $D_{R}$ dependence in (39) suggests that
as $D_{R}$ increases, the deviation from isotropy at a distance $\epsilon^{2}$
from the bifurcation continues to increase. This increase is shown in figures
7(b) and 7(a).

(a)

(b)

(c)

Figure 7: In the case of an $x$-only initial perturbation, a similar
dependence of the emerging steady state on $D_{R}$ can be seen in the two
plots of the nematic order parameter $\mathcal{N}(\boldsymbol{x},t)$ and the
direction of local nematic alignment for (a) $(D_{R},D_{T})=(0.0025,0.22)$
($D_{T}^{*}=0.24$) and (b) $(D_{R},D_{T})=(0.03125,0.105)$
($D_{T}^{*}=0.125$). Again, in both cases,
$\epsilon^{2}=D_{T}^{*}-D_{T}=0.02$. Figure (c) shows the vorticity field
$\omega(\boldsymbol{x},t)$ (colors) and velocity field
$\boldsymbol{u}(\boldsymbol{x},t)$ (arrows) for the $x$-only steady state with
$(D_{R},D_{T})=(0.03125,0.105)$.

Finally, we make a quantitative comparison between the expression (37) and the
numerical solution to the full system (8) with $\beta=0$ as
$\epsilon^{2}=D_{T}^{*}-D_{T}$ is varied. We fix $D_{R}=0.03125$, so
$D_{T}^{*}=0.125$, and consider five different values of $\epsilon^{2}$. For
each $\epsilon$, we allow the system to reach a steady state, and in figure 8,
we plot the difference between the predicted steady state $\Psi_{\rm
p}(x,y,\theta)$ (37) and the computed steady state $\Psi_{\rm c}(x,y,\theta)$
over a 1D slice of $x$-values for fixed $y$ and $\theta$. As $\epsilon$ is
decreased, the pointwise difference between the predicted and computed steady
states decreases like $\epsilon^{2}$, as expected. This is further
corroborated by table 1, which displays the maximum difference $\max_{0\leq
x,y,\theta\leq 2\upi}\left\lvert\Psi_{\rm p}(x,y,\theta)-\Psi_{\rm
c}(x,y,\theta)\right\rvert$.

(a)

(b)

Figure 8: Difference between the predicted steady state $\Psi_{\rm p}$ (37) and the computed steady state for (8) with $\beta=0$ for a fixed value of $y=\upi$ and two different fixed values of $\theta$: (a) $\theta=\upi$ and (b) $\theta=\frac{23\upi}{16}$. As expected, the difference between $\Psi_{\rm p}-\Psi_{\rm c}$ is $O(\epsilon^{2})$.  $\epsilon^{2}$ | 0.02 | 0.01 | 0.005 | 0.0025 | 0.00125
---|---|---|---|---|---
$\left\lVert\Psi_{\rm p}-\Psi_{\rm c}\right\rVert_{L^{\infty}(\mathbb{T}^{2}\times S^{1})}$ | 0.01667 | 0.00833 | 0.00421 | 0.00207 | 0.00107
Table 1: Maximum difference between the predicted steady state $\Psi_{\rm p}$
(37) and the computed steady state for (8) with $\beta=0$ for five different
values of $\epsilon$. The difference scales with $\epsilon^{2}$, as expected.

##  4 Motile particles: Hopf bifurcation

When $\beta>0$, the system (8) can experience a much richer catalogue of
bifurcations, as evidenced by figure 1. We focus first on the Hopf region
(roughly between C and E on figure 1), where the eigenvalue corresponding to
the $\left\lvert k\right\rvert=1$ mode is complex-valued. Fixing a very small
level of rotational diffusion (e.g. $D_{R}=0.001$), we choose a value of
translational diffusion $D_{R}\ll D_{T}<\frac{1}{9}$ such that for some value
of $\beta=\beta_{T}$, we have $\sigma(\pm 1)=\pm{\rm i}b_{T}$; i.e. the real
part of $\sigma(\pm 1)$ vanishes and a Hopf bifurcation occurs. Here the
subscript $T$ is used to denote that the values of $\beta_{T}$ and $b_{T}$
depend on the choice of $D_{T}$ through the implicit expression (19).

We perform a weakly nonlinear analysis of the Hopf bifurcation for the
$\left\lvert k\right\rvert=1$ mode in Section 4.1. We show that for relatively
large $D_{T}$, the bifurcation at $\beta=\beta_{T}$ is supercritical and a
stable limit cycle arises just beyond the bifurcation. However, for small
$D_{T}$, the bifurcation is subcritical for general initial perturbations in
both $x$ and $y$, but supercritical for $x$-only perturbations. We find
numerical evidence of hysteresis in this subcritical region, which we explore
in Section 4.2. In the supercritical region, we numerically locate an example
of the emerging stable limit cycle (see Section 4.4).

###  4.1 Weakly nonlinear analysis

As in the immotile case, we consider $\beta=\beta_{T}-\epsilon^{2}$ for
$\epsilon\ll 1$, so the $\left\lvert k\right\rvert=1$ modes are very slightly
unstable. We again consider a slow timescale $\tau=\epsilon^{2}t$ and expand
$\Psi$ and $\boldsymbol{u}$ in $\epsilon$ as (25). Plugging these expressions
into (8), at $O(\epsilon)$ we obtain the equation

$\mathcal{L}[\Psi_{1}]:=\p_{t}\Psi_{1}+\beta_{T}\boldsymbol{p}\bcdot\bnabla\Psi_{1}-2\boldsymbol{p}^{\rm
T}\bnabla\boldsymbol{u}_{1}\boldsymbol{p}-D_{T}\Delta\Psi_{1}-D_{R}\Delta_{p}\Psi_{1}=0.$
(42)

Note that when $\beta>0$ we no longer have an explicit expression for the
eigenvalues and eigenfunctions of the linearized operator (16) when $D_{R}>0$,
and therefore the following analysis must be performed in the limit of very
small rotational diffusion, which we treat perturbatively using Section 2.4.

Using the form (20) of the eigenvalues and eigenfunctions of the linearized
operator when $D_{R}=0$, and using the perturbative expression for $D_{R}>0$
from Section 2.4, we have that $\Psi_{1}$ and $\boldsymbol{u}_{1}$ are given
up to $O(D_{R})$ by

$\displaystyle\Psi_{1}$ $\displaystyle=c_{x}\psi_{x,1}(\theta)A_{x}(\tau){\rm
e}^{{\rm i}x+{\rm i}b_{T}t}+c_{y}\psi_{y,1}(\theta)A_{y}(\tau){\rm e}^{{\rm
i}y+{\rm i}b_{T}t}+{\rm c.c.}$ (43)
$\displaystyle\quad\psi_{x,1}(\theta)=\frac{\cos\theta\sin\theta}{D_{T}+{\rm
i}b_{T}+{\rm
i}\beta_{T}\cos\theta}+O(D_{R}),\quad\psi_{y,1}(\theta)=\frac{\cos\theta\sin\theta}{D_{T}+{\rm
i}b_{T}+{\rm i}\beta_{T}\sin\theta}+O(D_{R})$
$\displaystyle\boldsymbol{u}_{1}$ $\displaystyle=-\frac{\rm
i}{2}\big{(}c_{x}A_{x}{\rm e}^{{\rm i}x+{\rm
i}b_{T}t}\boldsymbol{e}_{y}+c_{y}A_{y}{\rm e}^{{\rm i}y+{\rm
i}b_{T}t}\boldsymbol{e}_{x}\big{)}+{\rm c.c.}+O(D_{R}).$

Again, $A_{x}(\tau)$ and $A_{y}(\tau)$ are complex-valued amplitudes which
depend only on the slow timescale $\tau$ and for which we wish to obtain an
equation. Here again ${\rm c.c.}$ denotes the complex conjugate of the
preceding terms.

At $O(\epsilon^{2})$ we obtain the equation

$\mathcal{L}[\Psi_{2}]=-\boldsymbol{u}_{1}\bcdot\bnabla\Psi_{1}-{\rm{div}}_{p}\big{(}({\bf
I}-\boldsymbol{p}\boldsymbol{p}^{\rm
T}(\bnabla\boldsymbol{u}_{1}\boldsymbol{p})\Psi_{1})\big{)},$ (44)

where the operator $\mathcal{L}$ is as defined in (42). Using (43), we can
explicitly calculate the expression on the right hand side of (44) up to terms
of $O(D_{R})$ (see Appendix B, equation (73) for the full expression). As in
the immotile case, using the form of the right hand side as a guide, we look
for $(\Psi_{2},\boldsymbol{u}_{2})$ of the form

$\displaystyle\Psi_{2}$ $\displaystyle=\psi_{2,1}{\rm e}^{{\rm i}(x+y)+2{\rm
i}b_{T}t}A_{x}A_{y}+\psi_{2,2}{\rm e}^{{\rm
i}(x-y)}A_{x}\overline{A}_{y}+\psi_{2,3}A_{x}^{2}{\rm e}^{2{\rm i}x+2{\rm
i}b_{T}t}+\psi_{2,4}A_{y}^{2}{\rm e}^{2{\rm i}y+2{\rm i}b_{T}t}$ (45)
$\displaystyle\hskip 85.35826pt+\psi_{2,5}\left\lvert
A_{x}\right\rvert^{2}+\psi_{2,6}\left\lvert A_{y}\right\rvert^{2}+{\rm c.c.},$
$\displaystyle\boldsymbol{u}_{2}$ $\displaystyle=-\frac{\rm
i}{8\upi}\int_{0}^{2\upi}\bigg{(}\psi_{2,1}{\rm e}^{{\rm i}(x+y)+2{\rm
i}b_{T}t}A_{x}A_{y}(\boldsymbol{e}_{x}-\boldsymbol{e}_{y})+\psi_{2,2}{\rm
e}^{{\rm
i}(x-y)}A_{x}\overline{A}_{y}(\boldsymbol{e}_{x}+\boldsymbol{e}_{y})\bigg{)}(\cos^{2}\theta-\sin^{2}\theta)d\theta$
$\displaystyle\qquad-\frac{\rm i}{4\upi}\int_{0}^{2\upi}\bigg{(}\psi_{2,3}{\rm
e}^{2{\rm i}x+2{\rm i}b_{T}t}A_{x}^{2}\boldsymbol{e}_{y}+\psi_{2,4}{\rm
e}^{2{\rm i}y+2{\rm
i}b_{T}t}A_{y}^{2}\boldsymbol{e}_{x}\bigg{)}\sin\theta\cos\theta\,d\theta+{\rm
c.c.},$

where $\psi_{2,j}=\psi_{2,j}(\theta)$. Inserting the ansatz (45) into the left
hand side of (44), we then match exponents with the right hand side and solve
for each of the coefficients $\psi_{2,j}$. Further details are contained in
Appendix B, equations (75) and (76), but we obtain that both $\psi_{2,5}$ and
$\psi_{2,6}$ are $O(D_{R}^{-1})$ for small $D_{R}$, while each of the other
coefficients are $O(1)$ in $D_{R}$ as $D_{R}\to 0$. Thus for sufficiently
small $D_{R}$, the coefficients $\psi_{2,5}$ and $\psi_{2,6}$ dominate the
behavior of $\Psi_{2}$. We may solve for $\psi_{2,5}$ and $\psi_{2,6}$
explicitly up to terms of $O(1)$ in $D_{R}$:

$\displaystyle\psi_{2,5}$
$\displaystyle=\frac{c_{x}^{2}}{D_{R}}\bigg{[}a_{1}\cos\theta+a_{2}(2\cos^{2}\theta-1)+a_{3}\cos^{3}\theta$
(46) $\displaystyle\quad+a_{4}\bigg{(}\log(D_{T}+{\rm i}b_{T}+{\rm
i}\beta_{T}\cos\theta)-\frac{1}{2\upi}\int_{0}^{2\upi}\log(D_{T}+{\rm
i}b_{T}+{\rm i}\beta_{T}\cos\theta)\,d\theta\bigg{)}+O(D_{R})\bigg{]}$
$\displaystyle\psi_{2,6}$
$\displaystyle=\frac{c_{y}^{2}}{D_{R}}\bigg{[}a_{1}\sin\theta+a_{2}(2\sin^{2}\theta-1)+a_{3}\sin^{3}\theta$
$\displaystyle\quad+a_{4}\bigg{(}\log(D_{T}+{\rm i}b_{T}+{\rm
i}\beta_{T}\sin\theta)-\frac{1}{2\upi}\int_{0}^{2\upi}\log(D_{T}+{\rm
i}b_{T}+{\rm i}\beta_{T}\sin\theta)\,d\theta\bigg{)}+O(D_{R})\bigg{]}$

where

$a_{1}=-\frac{{\rm i}(D_{T}+{\rm i}b_{T})^{2}}{2\beta_{T}^{3}},\quad
a_{2}=-\frac{D_{T}+{\rm i}b_{T}}{8\beta_{T}^{2}},\quad a_{3}=\frac{\rm
i}{6\beta_{T}},\quad a_{4}=\frac{(D_{T}+{\rm i}b_{T})^{3}}{2\beta_{T}^{4}}.$
(47)

Note that we must enforce
$\int_{0}^{2\upi}\psi_{2,5}\,d\theta=\int_{0}^{2\upi}\psi_{2,6}\,d\theta=0$ in
order to have
$\int_{\mathbb{T}^{2}}\int_{0}^{2\upi}\Psi_{2}\,d\theta\,d\boldsymbol{x}=0$,
i.e. to enforce that the total concentration of particles in the system is
preserved at $1/(2\upi)$.

We may thus rewrite (45) as

$\displaystyle\Psi_{2}$ $\displaystyle=\psi_{2,5}(\theta)\left\lvert
A_{x}\right\rvert^{2}+\psi_{2,6}(\theta)\left\lvert
A_{y}\right\rvert^{2}+O(D_{R}^{0}),\quad\boldsymbol{u}_{2}=O(D_{R}^{0}).$

At $O(\epsilon^{3})$ we obtain the equation

$\displaystyle\mathcal{L}[\Psi_{3}]$
$\displaystyle=-\p_{\tau}\Psi_{1}+\boldsymbol{p}\bcdot\bnabla\Psi_{1}-\boldsymbol{u}_{1}\bcdot\bnabla\Psi_{2}-\boldsymbol{u}_{2}\bcdot\bnabla\Psi_{1}$
(48) $\displaystyle\quad-{\rm{div}}_{p}\big{(}({\bf
I}-\boldsymbol{p}\boldsymbol{p}^{\rm
T})(\bnabla\boldsymbol{u}_{2}\boldsymbol{p})\Psi_{1})\big{)}-{\rm{div}}_{p}\big{(}({\bf
I}-\boldsymbol{p}\boldsymbol{p}^{\rm
T})(\bnabla\boldsymbol{u}_{1}\boldsymbol{p})\Psi_{2}\big{)}.$

Using (43) and (45), we may calculate the form of the right hand side. First,
we note that

$\displaystyle\p_{\tau}\Psi_{1}$
$\displaystyle=c_{x}\psi_{x,1}(\theta)(\p_{\tau}A_{x}){\rm e}^{{\rm i}x+{\rm
i}b_{T}t}+c_{y}\psi_{y,1}(\theta)(\p_{\tau}A_{y}){\rm e}^{{\rm i}y+{\rm
i}b_{T}t}+{\rm c.c.}$ (49) $\displaystyle\boldsymbol{p}\bcdot\bnabla\Psi_{1}$
$\displaystyle={\rm i}c_{x}\cos\theta\psi_{x,1}(\theta)A_{x}{\rm e}^{{\rm
i}x+{\rm i}b_{T}t}+{\rm i}c_{y}\sin\theta\psi_{y,1}(\theta)A_{y}{\rm e}^{{\rm
i}y+{\rm i}b_{T}t}+{\rm c.c.},$

where $\psi_{x,1}$ and $\psi_{y,1}$ are as in (43). The remaining terms on the
right hand side of (48) may be written

$\displaystyle-\boldsymbol{u}_{1}\bcdot\bnabla\Psi_{2}$
$\displaystyle-\boldsymbol{u}_{2}\bcdot\bnabla\Psi_{1}-{\rm{div}}_{p}\big{(}({\bf
I}-\boldsymbol{p}\boldsymbol{p}^{\rm
T})(\bnabla\boldsymbol{u}_{2}\boldsymbol{p})\Psi_{1})\big{)}-{\rm{div}}_{p}\big{(}({\bf
I}-\boldsymbol{p}\boldsymbol{p}^{\rm
T})(\bnabla\boldsymbol{u}_{1}\boldsymbol{p})\Psi_{2}\big{)}$ (50)
$\displaystyle=\mathcal{R}_{1}(\boldsymbol{x},\theta,t,\tau)+\mathcal{R}_{3}(\boldsymbol{x},\theta,t,\tau)+{\rm
c.c.},$ $\displaystyle\mathcal{R}_{1}(\boldsymbol{x},\theta,t,\tau)$
$\displaystyle=\frac{1}{D_{R}}\bigg{(}R_{xx}(\theta)c_{x}^{3}\left\lvert
A_{x}\right\rvert^{2}A_{x}{\rm e}^{{\rm i}x+{\rm
i}b_{T}t}+R_{xy}(\theta)c_{x}c_{y}^{2}\left\lvert
A_{y}\right\rvert^{2}A_{x}{\rm e}^{{\rm i}x+{\rm i}b_{T}t}$
$\displaystyle\quad+R_{yx}(\theta)c_{y}c_{x}^{2}\left\lvert
A_{x}\right\rvert^{2}A_{y}{\rm e}^{{\rm i}y+{\rm
i}b_{T}t}+R_{yy}(\theta)c_{y}^{3}\left\lvert A_{y}\right\rvert^{2}A_{y}{\rm
e}^{{\rm i}y+{\rm i}b_{T}t}\bigg{)}$
$\displaystyle\mathcal{R}_{3}(\boldsymbol{x},\theta,t,\tau)$
$\displaystyle=R_{2x^{+}}(\theta,\tau){\rm e}^{{\rm i}(2x+y)+{\rm
i}3b_{T}t}+R_{2y^{+}}(\theta,\tau){\rm e}^{{\rm i}(x+2y)+{\rm
i}3b_{T}t}+R_{2x^{-}}(\theta,\tau){\rm e}^{{\rm i}(2x-y)+{\rm i}b_{T}t}$
$\displaystyle\quad+R_{2y^{-}}(\theta,\tau){\rm e}^{{\rm i}(-x+2y)+{\rm
i}b_{T}t}+R_{3x}(\theta,\tau){\rm e}^{{\rm i}3x+{\rm
i}3b_{T}t}+R_{3y}(\theta,\tau){\rm e}^{{\rm i}3y+{\rm i}3b_{T}t}.$

The expressions for $R_{xx}$, $R_{xy}(\theta)$, $R_{yx}(\theta)$, and
$R_{yy}(\theta)$ are written up to $O(D_{R})$ in Appendix B, equation (77). As
in the immotile setting, the exact form of the coefficients in
$\mathcal{R}_{3}(\boldsymbol{x},\theta,t,\tau)$ will not be important due to
the solvability condition for equation (48).

In particular, in order for the $O(\epsilon^{3})$ equation (48) to have a
solution $\Psi_{3}$, the right hand side of (48) must satisfy the same
Fredholm condition (33) as in the immotile case, except now the operator
$\mathcal{L}$ defined in (42) is no longer self-adjoint. The adjoint
$\mathcal{L}^{*}$ is given by

$\displaystyle\mathcal{L}^{*}[\Phi]=-\p_{t}\Phi-\beta_{T}\boldsymbol{p}\bcdot\bnabla\Phi-2\boldsymbol{p}^{\rm
T}\bnabla\boldsymbol{u}\boldsymbol{p}-D_{T}\Delta\Phi-D_{R}\Delta_{p}\Phi,$

and, by the same calculation as in Sections 2.3 and 2.4, any $\Phi$ satisfying
$\mathcal{L}^{*}[\Phi]=0$ has the form

$\displaystyle\Phi=\alpha_{x}$ $\displaystyle\overline{\psi}_{x,1}(\theta){\rm
e}^{{\rm i}x+{\rm i}b_{T}t}+\alpha_{y}\overline{\psi}_{y,1}(\theta){\rm
e}^{iy+{\rm i}b_{T}t}+{\rm c.c.},\qquad\alpha_{x}^{2}+\alpha_{y}^{2}=1$ (51)
$\displaystyle\overline{\psi}_{x,1}$
$\displaystyle=\frac{\cos\theta\sin\theta}{D_{T}-{\rm i}b_{T}-{\rm
i}\beta_{T}\cos\theta}+O(D_{R}),\quad\overline{\psi}_{y,1}=\frac{\cos\theta\sin\theta}{D_{T}-{\rm
i}b_{T}-{\rm i}\beta_{T}\sin\theta}+O(D_{R}).$

Due to the form of $\mathcal{R}_{3}$ (50), we have that

$\displaystyle\int_{\mathbb{T}^{2}}\mathcal{R}_{3}(\boldsymbol{x},\theta,t,\tau)\Phi(\boldsymbol{x},\theta,t)\,d\boldsymbol{x}=0$

for $\Phi$ as in (51), so the Fredholm condition (33) is automatically
satisfied. Thus it remains to ensure that

$\displaystyle\int_{0}^{2\upi}\int_{\mathbb{T}^{2}}\bigg{(}-\p_{\tau}\Psi_{1}+\boldsymbol{p}\bcdot\bnabla\Psi_{1}+\mathcal{R}_{1}\bigg{)}\Phi\,d\boldsymbol{x}d\theta$
$\displaystyle=0$

for any $\Phi$ as in (51). Using the forms of (49), (50), and (51), this leads
to the following coupled system of equations for the amplitudes $A_{x}$ and
$A_{y}$:

$\displaystyle c_{x}\big{(}\quad M_{0}(\p_{\tau}A_{x})$
$\displaystyle=M_{3}A_{x}+\frac{1}{D_{R}}(c_{x}^{2}M_{1}\left\lvert
A_{x}\right\rvert^{2}+c_{y}^{2}M_{2}\left\lvert
A_{y}\right\rvert^{2})A_{x}\quad\big{)}$ (52) $\displaystyle c_{y}\big{(}\quad
M_{0}(\p_{\tau}A_{y})$
$\displaystyle=M_{3}A_{y}+\frac{1}{D_{R}}(c_{y}^{2}M_{1}\left\lvert
A_{y}\right\rvert^{2}+c_{x}^{2}M_{2}\left\lvert
A_{x}\right\rvert^{2})A_{y}\quad)$

where the complex constants $M_{0}$, $M_{1}$, $M_{2}$, and $M_{3}$ are given
by

$\displaystyle M_{0}$
$\displaystyle=\int_{0}^{2\upi}\psi_{x,1}^{2}(\theta)\,d\theta=\int_{0}^{2\upi}\psi_{y,1}^{2}(\theta)\,d\theta$
(53) $\displaystyle M_{1}$
$\displaystyle=\int_{0}^{2\upi}R_{xx}(\theta)\psi_{x,1}(\theta)\,d\theta=\int_{0}^{2\upi}R_{yy}(\theta)\psi_{y,1}(\theta)\,d\theta$
$\displaystyle M_{2}$
$\displaystyle=\int_{0}^{2\upi}R_{xy}(\theta)\psi_{x,1}(\theta)\,d\theta=\int_{0}^{2\upi}R_{yx}(\theta)\psi_{y,1}(\theta)\,d\theta$
$\displaystyle M_{3}$ $\displaystyle={\rm
i}\int_{0}^{2\upi}\psi_{x,1}^{2}(\theta)\cos\theta\,d\theta={\rm
i}\int_{0}^{2\upi}\psi_{y,1}^{2}(\theta)\sin\theta\,d\theta.$

Explicit expressions for each $M_{j}$ depending on $D_{T}$, $b_{T}$, and
$\beta_{T}$ are given in Appendix (B), equation (79). Note that although
$\psi_{x,1}(\theta)\neq\psi_{y,1}(\theta)$, etc., the coefficients $M_{j}$ for
the $x$ and $y$ directions are equal.

In the Hopf setting, each coefficient $M_{j}$ is now complex-valued. Writing
$A_{x}(\tau)=\rho_{x}(\tau){\rm e}^{{\rm i}\varphi_{x}(\tau)}$ and
$A_{y}(\tau)=\rho_{y}(\tau){\rm e}^{{\rm i}\varphi_{y}(\tau)}$, we may
separate (52) into equations for the magnitudes $\rho_{x}$, $\rho_{y}$ and the
phases $\varphi_{x}$, $\varphi_{y}$. We then look for conditions on the
coefficients $M_{j}$ such that (52) admits a nontrivial limit cycle satisfying
$\p_{\tau}\rho_{x}=\p_{\tau}\rho_{y}=0$.

We find that if $M_{0}\neq 0$, Re$\big{(}(M_{1}+M_{2})/M_{0}\big{)}\neq 0$,
and

$\frac{{\rm Re}\big{(}\frac{M_{3}}{M_{0}}\big{)}}{{\rm
Re}\big{(}\frac{M_{1}+M_{2}}{M_{0}}\big{)}}<0,$ (54)

then (52) gives rise to a nontrivial limit cycle with magnitudes $\rho_{x}$
and $\rho_{y}$ given by

$\rho_{x}=\pm\frac{\sqrt{D_{R}}}{c_{x}}\sqrt{-\frac{{\rm
Re}\big{(}\frac{M_{3}}{M_{0}}\big{)}}{{\rm
Re}\big{(}\frac{M_{1}+M_{2}}{M_{0}}\big{)}}},\quad\rho_{y}=\pm\frac{\sqrt{D_{R}}}{c_{y}}\sqrt{-\frac{{\rm
Re}\big{(}\frac{M_{3}}{M_{0}}\big{)}}{{\rm
Re}\big{(}\frac{M_{1}+M_{2}}{M_{0}}\big{)}}},$ (55)

while the phases satisfy

$\p_{\tau}\varphi_{x}=\p_{\tau}\varphi_{y}={\rm
Im}\bigg{(}\frac{M_{3}}{M_{0}}\bigg{)}-{\rm
Im}\bigg{(}\frac{M_{1}+M_{2}}{M_{0}}\bigg{)}\frac{{\rm
Re}\big{(}\frac{M_{3}}{M_{0}}\big{)}}{{\rm
Re}\big{(}\frac{M_{1}+M_{2}}{M_{0}}\big{)}}=:\varphi_{T}.$ (56)

In particular, if (54) holds, the stable limit cycle arising just after the
bifurcation, to leading order in $\epsilon=\sqrt{\beta_{T}-\beta}$, is given
by

$\Psi=\frac{1}{2\upi}\bigg{(}1\pm\epsilon\sqrt{D_{R}}\sqrt{-\frac{{\rm
Re}\big{(}\frac{M_{3}}{M_{0}}\big{)}}{{\rm
Re}\big{(}\frac{M_{1}+M_{2}}{M_{0}}\big{)}}}\bigg{(}\psi_{x,1}(\theta){\rm
e}^{{\rm i}x+{\rm i}(b_{T}+\epsilon^{2}\varphi_{T})t}\pm\psi_{y,1}(\theta){\rm
e}^{{\rm i}y+{\rm i}(b_{T}+\epsilon^{2}\varphi_{T})t}\bigg{)}\bigg{)}.$ (57)

If the condition (54) does not hold, the bifurcation at $\beta=\beta_{T}$ is
subcritical and, while the system behavior for $\beta<\beta_{T}$ is less
predictable, we may expect to see hysteresis if $\beta$ is then increased
beyond $\beta_{T}$. In particular, the system may remain in a stable,
nontrivial state well after the uniform isotropic state has also become
stable.

As in the immotile setting, we also consider the effects of an initial
perturbation in the $x$-direction only. When $c_{y}=0$, we obtain the single
equation

$M_{0}(\p_{\tau}A_{x})=M_{3}A_{x}+M_{1}\left\lvert A_{x}\right\rvert^{2}A_{x}$
(58)

for the $x$-direction amplitude $A_{x}$. In this case a stable limit cycle
emerges beyond the bifurcation if $M_{0}\neq 0$,
Re$\big{(}M_{1}/M_{0}\big{)}\neq 0$, and

$\frac{{\rm Re}\big{(}\frac{M_{3}}{M_{0}}\big{)}}{{\rm
Re}\big{(}\frac{M_{1}}{M_{0}}\big{)}}<0.$ (59)

If these conditions are satisfied, the magnitude and phase of the emerging
limit cycle are given by

$\rho_{x}=\pm\sqrt{D_{R}}\sqrt{-\frac{{\rm
Re}\big{(}\frac{M_{3}}{M_{0}}\big{)}}{{\rm
Re}\big{(}\frac{M_{1}}{M_{0}}\big{)}}},\quad\p_{\tau}\varphi_{x}={\rm
Im}\bigg{(}\frac{M_{3}}{M_{0}}\bigg{)}-{\rm
Im}\bigg{(}\frac{M_{1}}{M_{0}}\bigg{)}\frac{{\rm
Re}\big{(}\frac{M_{3}}{M_{0}}\big{)}}{{\rm
Re}\big{(}\frac{M_{1}}{M_{0}}\big{)}}=:\varphi_{T},$ (60)

and, to leading order in $\epsilon$, the emerging solution after the
bifurcation has the form

$\Psi=\frac{1}{2\upi}\bigg{(}1\pm\epsilon\sqrt{D_{R}}\sqrt{-\frac{{\rm
Re}\big{(}\frac{M_{3}}{M_{0}}\big{)}}{{\rm
Re}\big{(}\frac{M_{1}}{M_{0}}\big{)}}}\psi_{x,1}(\theta){\rm e}^{{\rm i}x+{\rm
i}(b_{T}+\epsilon^{2}\varphi_{T})t}\bigg{)}$ (61)

In figure 9, we plot each of $M_{0}$, Re$\big{(}M_{3}/M_{0}\big{)}$,
Re$\big{(}(M_{1}+M_{2})/M_{0}\big{)}$, and Re$\big{(}M_{1}/M_{0}\big{)}$ using
the perturbed dispersion relation (21) (see figure 2) with $D_{R}=0.001$ for
$\beta_{T}\in[0.2,0.7]$. This range of $\beta_{T}$ essentially covers all
$\beta_{T}$ for which a Hopf bifurcation exists and for which the perturbative
expression in $D_{R}$ is valid.

From figure 9(a), we see that $M_{0}\neq 0$ for all values of $\beta_{T}$ in
the region of interest, so division by $M_{0}$ always makes sense.
Furthermore, from figure 9(b), we see that Re$\big{(}M_{3}/M_{0}\big{)}$ is
always positive. Therefore Re$\big{(}(M_{1}+M_{2})/M_{0}\big{)}$ and
Re$\big{(}M_{1}/M_{0}\big{)}$ will determine the sign of the quantities of
interest in conditions (54) and (59), respectively. Interestingly, figure 9(c)
indicates that Re$\big{(}(M_{1}+M_{2})/M_{0}\big{)}$ changes sign for some
$\beta_{T}\in[0.2,0.7]$. Note that since $M_{1}$ and $M_{2}$ are calculated
only up to $O(D_{R})$, the precise location where
Re$\big{(}(M_{1}+M_{2})/M_{0}\big{)}=0$ cannot be determined since the
location may depend on lower order terms in $D_{R}$. However, for a
sufficiently small bifurcation value $\beta_{T}$ (determined by choosing
$D_{T}$ sufficiently large), we should see a stable limit cycle emerge beyond
the bifurcation, while for $\beta_{T}$ sufficiently large ($D_{T}$
sufficiently small), the bifurcation should be subcritical.

In contrast, Re$\big{(}M_{1}/M_{0}\big{)}<0$ for all $\beta_{T}\in[0.2,0.7]$
(see figure 9(d)), indicating that in the Hopf region, for initial
perturbations in the $x$-direction only, a stable limit cycle should always
arise immediately beyond the bifurcation.

(a)

(b)

(c)

(d)

Figure 9:  The relevant relationships among the coefficients $M_{0}$, $M_{1}$,
$M_{2}$, and $M_{3}$ of the amplitude equations (52) are plotted over the Hopf
bifurcation range $\beta_{T}\in[0.2,0.7]$. In particular, since the curves for
(a) $M_{0}$ and (b) Re$(M_{3}/M_{0})$ are both strictly positive for all such
$\beta_{T}$, the type of Hopf bifurcation is determined by (c)
Re$\big{(}(M_{1}+M_{2})/M_{0}\big{)}$ (for 2D initial perturbations) or (d)
Re$\big{(}M_{1}/M_{0}\big{)}$ (for 1D initial perturbations). In the 2D case,
there is a transition from supercritical to subcritical at some value of
$\beta_{T}\in[0.2,0.7]$, indicated by the vertical dotted line in (c).

In the following sections we explore these predictions numerically.

###  4.2 Subcritical region and bistability

We first consider the subcritical region for initial perturbations in both $x$
and $y$, where we find strong numerical evidence of hysteresis. For the
following simulations, we fix $D_{R}=0.001$ and $D_{T}=0.02$, so the
bifurcation occurs at roughly $\beta_{T}\approx 0.63$. According to figure
9(c), this value of $\beta_{T}$ lies well within the subcritical region for
generic initial perturbations with both $x$ and $y$ components.

Similar to the immotile bifurcation study in figure 5, we take our initial
condition to be a small random perturbation in both $x$ and $y$ to the uniform
isotropic state, and begin by running the simulation with $\beta=0.48$ for
$500t$, allowing the system to move away from the isotropic state. We then
increase $\beta$ by $0.03$ every $100t$ until $\beta=0.93$, well beyond the
bifurcation value of $\beta_{T}\approx 0.63$. We then decrease $\beta$ by
$0.05$ every $100t$ until we reach $\beta=0.53$, again passing through the
bifurcation point. The bifurcation value $\beta_{T}\approx 0.63$ is reached
from below at $t=1000$ and again from above at $t=2600$.

Again we keep track of the time-averaged rate of viscous dissipation
$\overline{\mathcal{P}}$ (40), except now the average is taken over each
constant value of $\beta$ throughout the simulation. We plot
$\overline{\mathcal{P}}$ versus $\beta$ in figure 10(a). In contrast to the
supercritical bifurcation seen in the immotile setting (figure 5), here we can
see clear hysteresis: as $\beta$ increases past the bifurcation at
$\beta_{T}\approx 0.63$, the system remains in a nontrivial state – away from
the uniform, isotropic state – well beyond the bifurcation point (up to about
$\beta=0.81$) before finally dropping down to the uniform, isotropic state
($\overline{\mathcal{P}}=0$). Then as $\beta$ is decreased, the system remains
in the uniform isotropic steady state until the uniform state loses stability
at $\beta=0.63$, after which the system transitions to a nontrivial state
again. This apparent region of bistability between $0.63<\beta<0.81$ is
characteristic of a subcritical bifurcation.

(a)

(b)

(c)

Figure 10: (a) The plot of the time-averaged viscous dissipation
$\overline{\mathcal{P}}$ versus $\beta$ displays bistability in the system
above the subcritical Hopf bifurcation at $\beta_{T}\approx 0.63$. (b)
Hysteresis is also evident in the plot of
$\left\lVert\boldsymbol{u}\right\rVert_{L^{2}(\mathbb{T}^{2})}^{2}$ over the
entire simulation described above. The bifurcation value of $\beta_{T}\approx
0.63$ is reached from below at $t=1000$ and again from above at $t=2600$.
After the bifurcation is passed from below, the system remains in a nontrivial
state well beyond the bifurcation value. (c) An almost periodic structure
appears after the bifurcation value is passed, which persists until
$\beta=0.84$ at $t=1700$. Here we plot
$\left\lVert\boldsymbol{u}\right\rVert_{L^{2}(\mathbb{T}^{2})}^{2}$ from
$t=1300$ to 1500, where $\beta=0.72$ ($t=1300$ to 1400) and $\beta=0.75$
($t=1400$ to 1500).

We also plot the $L^{2}$ norm of the velocity field
$\left\lVert\boldsymbol{u}\right\rVert_{L^{2}(\mathbb{T}^{2})}^{2}$ over the
entire course of the simulation in figure 10(b). Note that fluctuations in the
velocity field persist well past the bifurcation point, which is reached at
$t=1000$ ($\beta=\beta_{T}=0.63$). The fluctuations remain until $t=1700$
(where $\beta$ is increased to 0.84), when they quickly begin to decay. After
$\beta$ reaches a maximum of $\beta=0.93$ from $t=2000$ to 2100, the velocity
field stays motionless as $\beta$ is decreased again. After $\beta=0.63$ is
reached, $\left\lVert\boldsymbol{u}\right\rVert_{L^{2}(\mathbb{T}^{2})}^{2}$
begins to increase again. We note that this hysteretic behavior is replicated
even if we start the above procedure much closer to the bifurcation point –
e.g. at $\beta=0.6$. This indicates that the bistability we are seeing
directly relates to the subcritical nature of the bifurcation.

Looking more closely at the region immediately following the bifurcation at
$t=1000$ (see figure 10(c)), the system develops a very regular, nearly
periodic structure, especially from $t=1300$ to 1500 ($\beta=0.72$ and
$0.75$). This structure persists until $\beta$ is increased to 0.84 at
$t=1700$. Snapshots of the nematic order parameter
$\mathcal{N}(\boldsymbol{x},t)$ and particle concentration field
$c(\boldsymbol{x},t)$ along this upper solution branch at $\beta=0.75$ are
displayed in figure 15. In addition, the supplementary videos Movie2 and
Movie3 show the quasiperiodic nature of the dynamics along this upper branch.

(a)

(b)

Figure 11: Snapshots of (a) the nematic order parameter
$\mathcal{N}(\boldsymbol{x},t)$ and the direction of local nematic alignment,
and (b) the concentration field $c(\boldsymbol{x},t)$ along the nontrivial
upper solution branch which emerges above the subcritical Hopf bifurcation at
$\beta_{T}=0.63$. Here $\beta=0.75$, and figures are taken at successive local
peaks and valleys in the velocity $L^{2}$ norm (every 5 to $5.5t$), starting
with a peak.

The unexpectedly regular structure of the temporal dynamics in figures 10(c)
and 11 prompts a closer look at the nontrivial hysteretic state at
$\beta=0.75$. We simulate this state over a long time and, among other values,
record the velocity $\boldsymbol{u}(\boldsymbol{x},t)$ evaluated at the center
point of the computational domain. The value of the $x$-coordinate $u_{x}$
over 1500 time units is plotted in figure 12(a); the $y$-coordinate $u_{y}$
behaves similarly. The near-perfect periodicity here is striking. We plot the
power spectrum of $u_{x}(t)$ in figure 12(b) and note that, remarkably,
$u_{x}(t)$ decomposes into essentially just two temporal modes: a large mode
at frequency $0.096=1/10.4$ and a small mode at $0.62=1/16.2$. In figure 12(c)
we plot the signal

$s(t)=0.2\bigg{(}\sin\bigg{(}\frac{2\upi(t-2065)}{10.4}\bigg{)}+0.3\sin\bigg{(}\frac{2\upi(t-2065)}{16.2}\bigg{)}\bigg{)}$
(62)

on top of $u_{x}(t)$ for $200t$, where the value $2065$ was chosen to
qualitatively match with $u_{x}$. The overlap is nearly exact. The simplicity
of the signal in Figure 12 is surprising given that these dynamics occur in a
region beyond the predictive scope of the preceding weakly nonlinear analysis,
where we do not necessarily expect such a regular structure.

(a)

(b)

(c)

Figure 12: The surprisingly regular temporal dynamics in the nontrivial
hysteretic state at $\beta=0.75$. Figure (a) displays the near-periodic
dynamics of the $x$-coordinate $u_{x}$ of the velocity field
$\boldsymbol{u}(\boldsymbol{x},t)$ evaluated at the center point of the
computational domain. Figure (b) is the power spectrum of $u_{x}$ over the
time interval plotted in (a). The signal decomposes into just two temporal
modes. In figure (c), we plot $u_{x}(t)$ along with the simple signal $s(t)$
(62) composed of the two modes in (b). The agreement is nearly perfect.

###  4.3 Supercriticality for 1D initial perturbations

While an initial perturbation in both $x$ and $y$ gives rise to a subcritical
bifurcation at $\beta_{T}=0.63$, as predicted by figure 9, an initial
perturbation in only the $x$-direction should result in a supercritical
bifurcation. Indeed, if we perform the same numerical test as in figure 10 but
with an initial perturbation in only the $x$-direction instead, the resulting
relationship between the average viscous dissipation $\overline{\mathcal{P}}$
and $\beta$ (figure 13(a)) is characteristic of a supercritical bifurcation.
In particular, $\overline{\mathcal{P}}$ smoothly decreases to zero as the
bifurcation is approached from below, similar to the behavior seen in the
immotile bifurcation (figure 5). Furthermore, we can numerically locate the
limit cycle which emerges just below the bifurcation value of
$\beta_{T}=0.63$. Snapshots from a single period of this limit cycle are shown
in figure 13, and a few periods of the cycle are documented in the
supplementary video Movie4. The alignment among particles is very weak, but
they display a clear preferred direction which oscillates over time. Note that
the period of the limit cycle corresponds to every other peak of the velocity
$L^{2}$ norm (figure 13(b)); roughly $15$–$15.5t$. This may be compared with
the predicted period of $2\upi/b_{T}=2\upi/0.43\approx 14.6t$.

(a)

(b)

(c)

(d)

Figure 13:  (a) Plot of $\overline{\mathcal{P}}$ versus $\beta$ when
$D_{R}=0.001$ and $D_{T}=0.02$ (so $\beta_{T}\approx 0.63$) for initial
perturbations in the $x$-direction only. The behavior here can be contrasted
with figure 10 for 2D ($x$ and $y$) initial perturbations. (b) Plot of
$\left\lVert\boldsymbol{u}\right\rVert_{L^{2}(\mathbb{T}^{2})}^{2}$ over time
and (c) plot of the $x$-component (blue) and $y$-component (red) of the
velocity field $\boldsymbol{u}$ evaluated at the center point of the
computational domain. Here $\beta=0.6$ is fixed. Both plots show that a stable
limit cycle develops following an $x$-only initial perturbation. (d) Snapshots
of $\mathcal{N}(\boldsymbol{x},t)$ over the period of one limit cycle. The
particles are very weakly aligned here, but their preferred direction
oscillates.

###  4.4 Supercritical region (2D)

While initial perturbations in only the $x$-direction give rise to a
supercritical bifurcation for all $\beta_{T}\in[0.2,0.7]$, generic 2D
perturbations in both $x$ and $y$ should transition from subcritical to
supercritical near $\beta_{T}=0.5$. We fix $D_{R}=0.001$ and $D_{T}=0.075$ so
that $\beta_{T}\approx 0.40$, which according to figure 9(c) lies within the
supercritical region for 2D perturbations. Setting $\beta=0.38$, we simulate
the longtime system dynamics starting with a 2D initial perturbation. After an
initial period of slow growth, a stable limit cycle develops, as shown in
figure 14).

(a)

(b)

Figure 14: The stable limit cycle which develops just below the Hopf
bifurcation at $\beta_{T}=0.40$. Figure (a) shows the $x$-component (blue) and
$y$-component (red) of the velocity $\boldsymbol{u}$ evaluated at the center
point of the computational domain. Figure (b) shows the $L^{2}$ norm of the
velocity field
$\left\lVert\boldsymbol{u}\right\rVert_{L^{2}(\mathbb{T}^{2})}^{2}$ over time.

Although the peaks in
$\left\lVert\boldsymbol{u}\right\rVert_{L^{2}(\mathbb{T}^{2})}^{2}$ over time
are not perfectly equal in height, they occur at regular intervals ($\sim
6t$), and give rise to a regular, repeated pattern in the particle alignment
and generated velocity field, which can be viewed in the supplementary videos
Movie5 (nematic order parameter) and Movie6 (vorticity field and velocity
direction). Snapshots of a single period of the cycle are plotted in figure
15.

(a)

(b)

Figure 15: Snapshots of (a) the nematic order parameter
$\mathcal{N}(\boldsymbol{x},t)$ and the direction of local nematic alignment
and (b) the fluid vorticity and velocity fields over one period of the limit
cycle. The snapshots alternate between peaks and valleys in
$\left\lVert\boldsymbol{u}\right\rVert_{L^{2}(\mathbb{T}^{2})}^{2}$ in figure
14, starting with a peak. As in figure 6, the extensile flow produced by the
aligned dipoles is clear.

From figure 15 we can see that the period of the limit cycle corresponds to
every 4 peaks in
$\left\lVert\boldsymbol{u}\right\rVert_{L^{2}(\mathbb{T}^{2})}^{2}$, so every
$24t$. This can be compared with the predicted imaginary part of the growth
rate $b_{T}\approx 0.24$, which yields a period of $2\upi/b_{T}\approx 26t$.

The periodic behavior displayed in figure 15 may be compared to the noisy
quasiperiodic behavior along the upper solution branch in the bistable region
for the subcritical Hopf bifurcation in figure 11. The dynamics in the
supercritical case are much more regular over time.

##  5 Motile particles: pitchfork bifurcation

Keeping $D_{R}$ fixed and small, we now fix translational diffusion within the
range $1/9<D_{T}<1/4$ such that one of the (now purely real) eigenvalues
corresponding to the $\left\lvert k\right\rvert=1$ modes crosses zero for some
$\beta_{T}\in(0,\sqrt{3}/9)$. Again, the subscript $T$ is used to denote that
the bifurcation value of $\beta$ depends on the choice of translational
diffusion $D_{T}$.

In this parameter regime, the same weakly nonlinear calculations as in the
Hopf setting may be performed (see Section 4.1), except now $b_{T}=0$. We thus
arrive at the same form of amplitude equations as (52):

$\displaystyle c_{x}\big{(}\quad M_{0}(\p_{\tau}A_{x})$
$\displaystyle=M_{3}A_{x}+\frac{1}{D_{R}}(c_{x}^{2}M_{1}\left\lvert
A_{x}\right\rvert^{2}+c_{y}^{2}M_{2}\left\lvert
A_{y}\right\rvert^{2})A_{x}\quad\big{)}$ (63) $\displaystyle c_{y}\big{(}\quad
M_{0}(\p_{\tau}A_{y})$
$\displaystyle=M_{3}A_{y}+\frac{1}{D_{R}}(c_{y}^{2}M_{1}\left\lvert
A_{y}\right\rvert^{2}+c_{x}^{2}M_{2}\left\lvert
A_{x}\right\rvert^{2})A_{y}\quad),$

but now $M_{0}$, $M_{1}$, $M_{2}$, and $M_{3}$ given in Appendix B, equation
(79) are all real-valued for each $(D_{T},\beta_{T})$ in the region of
interest. Since the coefficients $M_{j}$ are real-valued, similar to the
immotile setting, we may look for conditions under which (63) admits a
nontrivial steady state solution $\p_{\tau}A_{x}=\p_{\tau}A_{y}=0$. If the
coefficients $M_{j}$ satisfy $M_{0}\neq 0$, $M_{1}+M_{2}\neq 0$, and

$\frac{M_{3}}{M_{1}+M_{2}}<0$ (64)

for initial perturbations in both $x$ and $y$, or $M_{0}\neq 0$, $M_{1}\neq
0$, and

$\frac{M_{3}}{M_{1}}<0$ (65)

for initial perturbations in $x$ only, then (63) admits nontrivial steady
states of the form

$A_{x}=\pm\frac{\sqrt{D_{R}}}{c_{x}}\sqrt{-\frac{M_{3}}{M_{1}+M_{2}}},\quad
A_{y}=\pm\frac{\sqrt{D_{R}}}{c_{y}}\sqrt{-\frac{M_{3}}{M_{1}+M_{2}}}$ (66)

for 2D ($x$ and $y$) initial perturbations or

$A_{x}=\pm\sqrt{D_{R}}\sqrt{-\frac{M_{3}}{M_{1}}}$ (67)

for $x$-only perturbations. Then, to leading order in
$\epsilon=\sqrt{\beta_{T}-\beta}$, a stable steady state emerges after the
real eigenvalue crossing of the form

$\Psi=\frac{1}{2\upi}\bigg{(}1\pm\epsilon\sqrt{D_{R}}\sqrt{-\frac{M_{3}}{M_{1}+M_{2}}}\big{(}\psi_{x,1}(\theta){\rm
e}^{{\rm i}x}\pm\psi_{y,1}(\theta){\rm e}^{{\rm i}y}\big{)}\bigg{)}+{\rm
c.c.}$ (68)

for initial perturbations in both $x$ and $y$ and

$\Psi=\frac{1}{2\upi}\bigg{(}1\pm\epsilon\sqrt{D_{R}}\sqrt{-\frac{M_{3}}{M_{1}}}\psi_{x,1}(\theta){\rm
e}^{{\rm i}x}\bigg{)}+{\rm c.c.}$ (69)

for initial perturbations in the $x$-direction only. If the conditions (64)
and (65) do not hold for 2D and 1D perturbations, respectively, then the
bifurcation is subcritical and the system behavior beyond the bifurcation
value is less predictable.

As in the Hopf setting, using the expressions from Appendix B, equation (79),
we plot $M_{0}$, $M_{3}$, $M_{1}+M_{2}$, and $M_{1}$ over the desired range of
$(D_{T},\beta_{T})$ in figure 16. Here we again use the perturbed dispersion
relation (21) of Section 2.4 (see figure 2) with $D_{R}=0.001$. We find that
both $M_{0}$ and $M_{3}$ are always positive within our region of interest,
although $M_{0}$ appears to approach 0 as $\beta_{T}\to\sqrt{3}/9\approx
0.192$ while $M_{3}\to 0$ as $\beta_{T}\to 0$. Thus, as in the Hopf case, the
existence of a new stable steady state following the real eigenvalue crossing
is determined by the sign of $M_{1}+M_{2}$ (for 2D initial perturbations in
$x$ and $y$) or the sign of $M_{1}$ (for $x$-only perturbations). From figures
16(b) and 16(c), we see that both $M_{1}+M_{2}$ and $M_{1}$ are negative for
small $\beta_{T}$, indicating the emergence of a nontrivial stable steady
state after the bifurcation, but both $M_{1}+M_{2}$ and $M_{1}$ are positive
for larger $\beta_{T}$, indicating that the bifurcation type switches to
subcritical somewhere in the interval $0<\beta_{T}<\sqrt{3}/9$.

(a)

(b)

(c)

Figure 16: Plots of the coefficients (a) $M_{0}$ and $M_{3}$, (b)
$M_{1}+M_{2}$, and (c) $M_{1}$ given by the expressions (53) in the region
$0<\beta_{T}<\sqrt{3}/9$ (or $1/9<D_{T}<1/4$), where each $M_{j}$ is real-
valued. The vertical dotted lines in (b) and (c) indicate the value of
$\beta_{T}$ where the pitchfork bifurcation transitions from supercritical to
subcritical for 2D and 1D initial perturbations, respectively.

We again explore the supercritical and subcritical regions numerically in the
following sections.

###  5.1 Supercritical region

When $\beta_{T}$ is small, we expect dynamics near the bifurcation to look
very similar to the immotile case, where the isotropic steady state loses
stability to a supercritical pitchfork bifurcation. We fix $D_{R}=0.001$ and
$D_{T}=0.23$, so the bifurcation occurs at roughly $\beta_{T}\approx 0.09$.
According to figure 16, this value of $\beta_{T}$ lies within the
supercritical region for both 2D ($x$ and $y$) and 1D ($x$-only) initial
perturbations.

We begin with a small random perturbation to the uniform, isotropic state in
both $x$ and $y$. We initialize the simulation with $\beta=0$ until $t=1500$
and then increase $\beta$ by 0.02 every $100t$ until $\beta=0.2$. We plot the
average viscous dissipation $\overline{\mathcal{P}}$ (40) versus each
different value of $\beta$ in figure 17(a). Again, the relationship between
$\overline{\mathcal{P}}$ and $\beta$ supports the expectation of
supercriticality. We note that the viscous dissipation
$\overline{\mathcal{P}}$ behaves qualitatively the same as figure 17(a) using
an $x$-only initial perturbation.

The emergent stable steady state for $\beta<\beta_{T}$ is plotted in figure
17(b) for a 2D initial perturbation and in figure 17(c) for an initial
$x$-only perturbation. Not surprisingly, in both cases the steady state
essentially looks like the stable states which arise following the immotile
bifurcation (figure 6). The calculation in the immotile case helps to explain
the very weak alignment seen in the emerging steady state, since the
rotational diffusion $D_{R}$ is very small in this setting.

(a)

(b)

(c)

Figure 17: (a) Plot of $\overline{\mathcal{P}}$ versus $\beta$ using
$D_{R}=0.001$, $D_{T}=0.23$ (so $\beta_{T}\approx 0.09$), and a 2D ($x$ and
$y$) initial perturbation. The bifurcation is supercritical, and we plot the
nematic order parameter $\mathcal{N}(\boldsymbol{x},t)$ and preferred local
alignment direction for the stable state which emerges just beyond the
bifurcation in the case of (b) a 2D initial perturbation and (c) an $x$-only
initial perturbation.

###  5.2 Subcritical region

We next fix $D_{R}=0.001$ and $D_{T}=0.13$, so $\beta_{T}\approx 0.19$ is the
bifurcation value. According to figure 16, both 2D ($x$ and $y$) and 1D
($x$-only) initial perturbations to the isotropic state give rise to a
subcritical bifurcation at this value of $\beta_{T}$.

Using a small, random 2D perturbation to the uniform, isotropic state as our
initial condition, we simulate the model dynamics until $t=3000$ using
$\beta=0.17$, just on the unstable side of $\beta_{T}$. As predicted by the
amplitude equation coefficients in figure 16, and in contrast to the
supercritical setting, the system does not settle into a stable nontrivial
steady state. Rather, the system undergoes a relatively quick spike in
activity before settling into what appears to be a stable limit cycle, as
shown in figure 18 as well as the supplementary video Movie7).

(a)

(b)

(c)

Figure 18: For an initial perturbation in $x$ and $y$, after a quick spike in
activity, the system settles into a limit cycle following the subcritical
pitchfork bifurcation. Figure (a) shows
$\left\lVert\boldsymbol{u}\right\rVert_{L^{2}(\mathbb{T}^{2})}^{2}$ over time.
Figure (b) shows the $x$-component (blue) and $y$-component (red) of
$\boldsymbol{u}$ evaluated at the center of the domain over time. Figure (c)
displays snapshots of the nematic order parameter
$\mathcal{N}(\boldsymbol{x},t)$ and direction of local nematic alignment at
successive peaks and valleys in
$\left\lVert\boldsymbol{u}\right\rVert_{L^{2}(\mathbb{T}^{2})}^{2}$ over one
period of the cycle, starting with a peak.

The fast initial spike in activity may be contrasted with the supercritical
Hopf bifurcation, where the oscillations grow slowly in amplitude over a much
longer time before reaching the near-periodic dynamics shown in figure 14.
Since the isotropic state is predicted to lose stability through a subcritical
pitchfork bifurcation for this parameter combination, we may be seeing the
results of a second bifurcation. Unlike the subcritical Hopf setting, this
behavior does not persist beyond the bifurcation value – the system appears to
quickly return to the isotropic state when $\beta$ is increased above
$\beta_{T}$.

Similar behavior is observed for $x$-only initial perturbations, as shown in
figure 19. In particular, using the same parameter values as in the 2D
setting, we see a quick initial spike in activity which then settles into what
appears to nearly be a limit cycle (the amplitude here is very slightly
decreasing over time). Again, this near-periodic behavior may be the result of
a second bifurcation beyond the subcritical pitchfork. In both the 2D ($x$ and
$y$) and 1D ($x$-only) cases, however, we find no evidence of any type of
hysteresis like what is seen in the subcritical Hopf setting.

(a)

(b)

(c)

Figure 19: An initial $x$-only perturbation also results in a limit cycle
following the subcritical pitchfork bifurcation. Figure (a) shows
$\left\lVert\boldsymbol{u}\right\rVert_{L^{2}(\mathbb{T}^{2})}^{2}$ over time.
Figure (b) shows the $x$-component (blue) and $y$-component (red) of
$\boldsymbol{u}$ evaluated at the center of the domain over time. Figure (c)
shows the nematic order parameter $\mathcal{N}(\boldsymbol{x},t)$ and
direction of local nematic alignment at the $L^{2}$ norm peak and valley,
respectively. The system oscillates slowly between the two states pictured.

##  6 Discussion

We have determined exactly how the uniform isotropic steady state loses
stability in the Saintillan–Shelley kinetic model of a dilute suspension of
active rodlike particles. Our weakly nonlinear analysis reveals a surprisingly
complex array of possible bifurcations that the system may undergo as the
nondimensional swimming speed (ratio of swimming speed to concentration and
active stress magnitude) is decreased below the instability threshold. The
type of bifurcation depends on the particle diffusivity, and passes from
supercritical pitchfork (relatively high diffusion or immotile particles) to
subcritical pitchfork to supercritical Hopf to subcritical Hopf as the
diffusivity is decreased. The analysis is supported by numerical examples of
each different bifurcation in the model. The numerical examples include
uncovering surprising structures in the hysteretic solution that is bistable
with the uniform isotropic state in the subcritical Hopf region, as well as
locating different 1D and 2D limit cycles and steady states which emerge
following the supercritical Hopf and pitchfork bifurcations.

The bifurcation analysis presented here provides mathematical insight into the
onset of collective particle motion and active turbulence. It would be
interesting to see if the full range of different transitions to collective
motion can be realized experimentally. The patterns observed here may depend
strongly on the aspect ratio of the periodic domain, as evidenced especially
by the difference in behavior for 1D versus 2D initial perturbations in the
Hopf region (section 4.2). In particular, a long, thin 2D domain may give rise
to patterns more similar to the 1D case. Furthermore, while a similar linear
stability picture to section 2.3 holds in 3D (Hohenegger & Shelley, 2010), the
effects of an additional spatial and orientational dimension on sub- versus
supercriticality remain unclear. These questions both warrant further study.

From a modeling perspective, it would be useful to perform a similar
bifurcation analysis for different closure models derived from the kinetic
theory to see to what extent different closures can capture the complexity of
the transition to collective behavior seen in the full kinetic model. It also
may be interesting to analyze how the inclusion of steric interactions among
particles (Gao et al., 2017) affects the nature of the possible types of
bifurcations to collective motion. From a mathematical perspective, it would
be interesting to prove that the uniform isotropic steady state is indeed
stable above the bifurcation at $\beta=\beta_{T}$ and that this stability
boundary is sharp. Our numerical tests indicate that the eigenvalues
calculated in Section 2.3 provide a fairly full picture of the stability
boundaries, but it would be reassuring to verify this rigorously, particularly
in the absence of translational or rotational diffusion. See forthcoming work
(Albritton & Ohm, 2022) addressing the stabilizing effects of the swimming
term in the Saintillan–Shelley model.

[Supplementary data]Supplementary movies are available online.

[Acknowledgements]We thank Scott Weady for supplying a base code for numerics.

[Funding]LO acknowledges support from NSF postdoctoral fellowship DMS-2001959.
MJS acknowledges support from NSF grants DMR-2004469 and DMR-1420073 (NYU-
MRSEC).

[Declaration of interests] The authors report no conflict of interest.

##  Appendix A Immotile bifurcation

Using the expressions for $\Psi_{1}$ and $\boldsymbol{u}_{1}$ (27), we can
calculate the right hand side of equation (28), given by

$\displaystyle-\boldsymbol{u}_{1}\bcdot\bnabla\Psi_{1}$
$\displaystyle-{\rm{div}}_{p}\big{(}({\bf I}-\boldsymbol{p}\boldsymbol{p}^{\rm
T})(\bnabla\boldsymbol{u}_{1}\boldsymbol{p})\Psi_{1})\big{)}$ (70)
$\displaystyle=G_{1}(\theta){\rm e}^{{\rm i}(x+y)}A_{x}A_{y}+G_{2}(\theta){\rm
e}^{{\rm i}(x-y)}A_{x}\overline{A}_{y}$
$\displaystyle\qquad+G_{3}(\theta)({\rm e}^{2{\rm i}x}A_{x}^{2}+\left\lvert
A_{x}\right\rvert^{2})+G_{4}(\theta)({\rm e}^{2{\rm i}y}A_{y}^{2}+\left\lvert
A_{y}\right\rvert^{2})+{\rm c.c.},$

where

$\displaystyle G_{1}$
$\displaystyle=-\frac{c_{x}c_{y}}{8}\big{(}\sin^{4}\theta+\cos^{4}\theta-6\cos^{2}\theta\sin^{2}\theta+2\cos\theta\sin\theta\big{)},$
$\displaystyle G_{2}$
$\displaystyle=-\frac{c_{x}c_{y}}{8}\big{(}\sin^{4}\theta+\cos^{4}\theta-6\cos^{2}\theta\sin^{2}\theta-2\cos\theta\sin\theta\big{)},$
$\displaystyle G_{3}$
$\displaystyle=-\frac{c_{x}^{2}}{8}\big{(}\cos^{4}\theta-3\cos^{2}\theta\sin^{2}\theta\big{)},\quad
G_{4}=-\frac{c_{y}^{2}}{8}\big{(}\sin^{4}\theta-3\sin^{2}\theta\cos^{2}\theta\big{)}.$

The form of (70) leads us to consider the ansatz (29) for $\Psi_{2}$ and
$\boldsymbol{u}_{2}$. Plugging this ansatz into the operator $\mathcal{L}$,
the left hand side of the $O(\epsilon^{2})$ equation (28) takes the form

$\displaystyle\mathcal{L}[\Psi_{2}]$
$\displaystyle=-\frac{1}{4\upi}(\cos^{2}\theta-\sin^{2}\theta)\int_{0}^{2\upi}\bigg{(}\psi_{2,1}{\rm
e}^{{\rm i}(x+y)}A_{x}A_{y}+\psi_{2,2}{\rm e}^{{\rm
i}(x-y)}A_{x}\overline{A}_{y}\bigg{)}(\cos^{2}\theta-\sin^{2}\theta)d\theta$
(71)
$\displaystyle\quad-\frac{1}{\upi}\cos\theta\sin\theta\int_{0}^{2\upi}\bigg{(}\psi_{2,3}{\rm
e}^{2{\rm i}x}A_{x}^{2}+\psi_{2,4}{\rm e}^{2{\rm
i}y}A_{y}^{2}\bigg{)}\sin\theta\cos\theta d\theta$
$\displaystyle\quad+\bigg{(}\frac{1}{4}-4D_{R}\bigg{)}\big{(}2\psi_{2,1}{\rm
e}^{{\rm i}(x+y)}A_{x}A_{y}+2\psi_{2,2}{\rm e}^{{\rm
i}(x-y)}A_{x}\overline{A}_{y}+4\psi_{2,3}{\rm e}^{2{\rm
i}x}A_{x}^{2}+4\psi_{2,4}{\rm e}^{2{\rm i}y}A_{y}^{2}\big{)}$
$\displaystyle\quad-D_{R}\bigg{(}(\p_{\theta\theta}\psi_{2,1}){\rm e}^{{\rm
i}(x+y)}A_{x}A_{y}+(\p_{\theta\theta}\psi_{2,2}){\rm e}^{{\rm
i}(x-y)}A_{x}\overline{A}_{y}+(\p_{\theta\theta}\psi_{2,3}){\rm e}^{2{\rm
i}x}A_{x}^{2}$ $\displaystyle\hskip
56.9055pt+(\p_{\theta\theta}\psi_{2,4}){\rm e}^{2{\rm
i}y}A_{y}^{2}+(\p_{\theta\theta}\psi_{2,5})\left\lvert
A_{x}\right\rvert^{2}+(\p_{\theta\theta}\psi_{2,6})\left\lvert
A_{y}\right\rvert^{2}\bigg{)}+{\rm c.c.}$

Matching the exponents in (70) to (71), we obtain a series of six independent
integrodifferential ODEs which may each be solved to yield the coefficients
$\psi_{2,j}$ listed in (30).

We also need to keep track of the components of $R_{x}(\theta,\tau){\rm
e}^{{\rm i}x}$ and $R_{y}(\theta,\tau){\rm e}^{{\rm i}y}$ appearing in the
expression $\mathcal{R}$ (32) for the right hand side of the $O(\epsilon^{3})$
equation (31).

Each term of $\mathcal{R}$ contributes the following to
$R_{x}(\theta,\tau){\rm e}^{{\rm i}x}$ and $R_{y}(\theta,\tau){\rm e}^{{\rm
i}y}$:

$\displaystyle\p_{\tau}\Psi_{1}:$
$\displaystyle\quad\cos\theta\sin\theta\big{(}c_{x}(\p_{\tau}A_{x}){\rm
e}^{{\rm i}x}+c_{y}(\p_{\tau}A_{y}){\rm e}^{{\rm i}y}\big{)}+{\rm c.c.}$ (72)
$\displaystyle\Delta\Psi_{1}:$
$\displaystyle\quad-\cos\theta\sin\theta\big{(}c_{x}A_{x}{\rm e}^{{\rm
i}x}+c_{y}A_{y}{\rm e}^{{\rm i}y}\big{)}+{\rm c.c.}$
$\displaystyle\boldsymbol{u}_{1}\bcdot\bnabla\Psi_{2}:$
$\displaystyle\quad-\frac{\sin\theta\cos\theta}{64D_{R}}\bigg{(}c_{x}c_{y}^{2}{\rm
e}^{{\rm i}x}A_{x}\left\lvert A_{y}\right\rvert^{2}+c_{x}^{2}c_{y}{\rm
e}^{{\rm i}y}\left\lvert A_{x}\right\rvert^{2}A_{y}\bigg{)}+{\rm c.c.}$
$\displaystyle{\rm{div}}_{p}\big{(}({\bf I}-\boldsymbol{p}\boldsymbol{p}^{\rm
T})(\bnabla\boldsymbol{u}_{1}\boldsymbol{p})\Psi_{2}\big{)}:$
$\displaystyle\bigg{(}r_{1}\cos\theta\sin\theta(4\sin^{2}\theta-6\sin^{4}\theta-1)+r_{2}\cos^{3}\theta\sin\theta(2-3\cos^{2}\theta)$
$\displaystyle-(r_{3}+r_{4})\sin\theta\cos\theta\bigg{)}\bigg{(}c_{x}c_{y}^{2}{\rm
e}^{{\rm i}x}A_{x}\left\lvert A_{y}\right\rvert^{2}+c_{x}^{2}c_{y}{\rm
e}^{{\rm i}y}\left\lvert A_{x}\right\rvert^{2}A_{y}\bigg{)}$
$\displaystyle-\cos\theta\sin\theta\bigg{(}3r_{5}\cos^{4}\theta+2r_{6}\cos^{2}\theta+r_{7}\bigg{)}\bigg{(}c_{x}^{3}{\rm
e}^{{\rm i}x}\left\lvert A_{x}\right\rvert^{2}A_{x}+c_{y}^{3}{\rm e}^{{\rm
i}y}\left\lvert A_{y}\right\rvert^{2}A_{y}\bigg{)}+{\rm c.c.}$

where

$\displaystyle r_{1}$ $\displaystyle=-\frac{1}{2(1+16D_{R})},\quad
r_{2}=-\frac{1}{64D_{R}},\quad r_{3}=\frac{3}{8(1+16D_{R})},\quad
r_{4}=\frac{3}{512D_{R}},$ $\displaystyle r_{5}$
$\displaystyle=-\frac{1+8D_{R}}{64D_{R}},\quad
r_{6}=\frac{3(1-16D_{R})}{32(1-12D_{R})},\quad
r_{7}=3\frac{32D_{R}^{2}-12D_{R}+1}{512D_{R}(1-12D_{R})}.$

The Fredholm condition (33) then requires that

$\displaystyle\alpha_{x}F[R_{x}(\theta,\tau)]{\rm e}^{{\rm
i}x}+\alpha_{y}F[R_{y}(\theta,\tau)]{\rm e}^{{\rm i}y}+{\rm
c.c.}=0,\qquad\alpha_{x}^{2}+\alpha_{y}^{2}=1,$

where again $F[\bcdot]=\int_{0}^{2\upi}\bcdot\,\cos\theta\sin\theta\,d\theta$.
In particular, we need $F[R_{x}]=F[R_{y}]=0$. The contribution of each term in
(72) to $F[R_{x}]$ is given by

$\displaystyle\p_{\tau}\Psi_{1}:\quad\frac{\upi}{4}c_{x}(\p_{\tau}A_{x}),\qquad\Delta\Psi_{1}:\quad-\frac{\upi}{4}c_{x}A_{x},\qquad\boldsymbol{u}_{1}\bcdot\bnabla\Psi_{2}:\quad\frac{\upi}{4}\frac{c_{x}c_{y}^{2}}{8(1-8D_{R})}\left\lvert
A_{y}\right\rvert^{2}A_{x},$ $\displaystyle{\rm{div}}_{p}\big{(}({\bf
I}-\boldsymbol{p}\boldsymbol{p}^{\rm
T})(\bnabla\boldsymbol{u}_{1}\boldsymbol{p})\Psi_{2}\big{)}:\quad\frac{\upi}{4}\bigg{(}-\frac{c_{x}c_{y}^{2}(7+48D_{R})}{1024D_{R}(1+16D_{R})}\left\lvert
A_{y}\right\rvert^{2}+c_{x}^{3}\frac{3(3-28D_{R}-32D_{R}^{2})}{1024D_{R}(1-12D_{R})}\left\lvert
A_{x}\right\rvert^{2}\bigg{)}A_{x}.$

Similarly, the contribution of each term in (72) to $F[R_{y}]$ is given by

$\displaystyle\p_{\tau}\Psi_{1}:\quad\frac{\upi}{4}c_{y}(\p_{\tau}A_{y}),\qquad\Delta\Psi_{1}:\quad-\frac{\upi}{4}c_{y}A_{y},\qquad\boldsymbol{u}_{1}\bcdot\bnabla\Psi_{2}:\quad\frac{\upi}{4}\frac{c_{x}^{2}c_{y}}{8(1-8D_{R})}\left\lvert
A_{x}\right\rvert^{2}A_{y},$ $\displaystyle{\rm{div}}_{p}\big{(}({\bf
I}-\boldsymbol{p}\boldsymbol{p}^{\rm
T})(\bnabla\boldsymbol{u}_{1}\boldsymbol{p})\Psi_{2}\big{)}:\quad\frac{\upi}{4}\bigg{(}-\frac{c_{x}^{2}c_{y}(7+48D_{R})}{1024D_{R}(1+16D_{R})}\left\lvert
A_{x}\right\rvert^{2}+c_{y}^{3}\frac{3(3-28D_{R}-32D_{R}^{2})}{1024D_{R}(1-12D_{R})}\left\lvert
A_{y}\right\rvert^{2}\bigg{)}A_{y}.$

##  Appendix B Hopf bifurcation

As in the immotile case, we may use the expressions (43) for $\Psi_{1}$ and
$\boldsymbol{u}_{1}$ to compute the right hand side of equation (44) up to
$O(D_{R})$:

$\displaystyle-\boldsymbol{u}_{1}\bcdot\bnabla\Psi_{1}$
$\displaystyle-{\rm{div}}_{p}\big{(}({\bf I}-\boldsymbol{p}\boldsymbol{p}^{\rm
T})(\bnabla\boldsymbol{u}_{1}\boldsymbol{p})\Psi_{1})\big{)}$ (73)
$\displaystyle=H_{1}(\theta){\rm e}^{{\rm i}(x+y)+2{\rm
i}b_{T}t}A_{x}A_{y}+H_{2}(\theta){\rm e}^{{\rm i}(x-y)}A_{x}\overline{A}_{y}$
$\displaystyle\qquad+H_{3}(\theta)({\rm e}^{2{\rm i}x+2{\rm
i}b_{T}t}A_{x}^{2}+\left\lvert A_{x}\right\rvert^{2})+H_{4}(\theta)({\rm
e}^{2{\rm i}y+2{\rm i}b_{T}t}A_{y}^{2}+\left\lvert A_{y}\right\rvert^{2})+{\rm
c.c.},$

where

$\displaystyle H_{1}$
$\displaystyle=-\frac{c_{x}c_{y}}{2}\bigg{[}\frac{-3\cos^{2}\theta\sin^{2}\theta+\sin^{4}\theta+\cos\theta\sin\theta}{D_{T}+{\rm
i}b_{T}+{\rm i}\beta_{T}\cos\theta}-{\rm
i}\beta_{T}\frac{\sin^{4}\theta\cos\theta}{(D_{T}+{\rm i}b_{T}+{\rm
i}\beta_{T}\cos\theta)^{2}}$
$\displaystyle\qquad+\frac{-3\cos^{2}\theta\sin^{2}\theta+\cos^{4}\theta+\cos\theta\sin\theta}{D_{T}+{\rm
i}b_{T}+{\rm i}\beta_{T}\sin\theta}-{\rm
i}\beta_{T}\frac{\cos^{4}\theta\sin\theta}{(D_{T}+{\rm i}b_{T}+{\rm
i}\beta_{T}\sin\theta)^{2}}\bigg{]}+O(D_{R})$ $\displaystyle H_{2}$
$\displaystyle=-\frac{c_{x}c_{y}}{2}\bigg{[}\frac{-3\cos^{2}\theta\sin^{2}\theta+\sin^{4}\theta-\cos\theta\sin\theta}{D_{T}+{\rm
i}b_{T}+{\rm i}\beta_{T}\cos\theta}-{\rm
i}\beta_{T}\frac{\sin^{4}\theta\cos\theta}{(D_{T}+{\rm i}b_{T}+{\rm
i}\beta_{T}\cos\theta)^{2}}$
$\displaystyle\qquad+\frac{-3\cos^{2}\theta\sin^{2}\theta+\cos^{4}\theta-\cos\theta\sin\theta}{D_{T}-{\rm
i}b_{T}-{\rm i}\beta_{T}\sin\theta}+{\rm
i}\beta_{T}\frac{\cos^{4}\theta\sin\theta}{(D_{T}-{\rm i}b_{T}-{\rm
i}\beta_{T}\sin\theta)^{2}}\bigg{]}+O(D_{R})$ $\displaystyle H_{3}$
$\displaystyle=-\frac{c_{x}^{2}}{2}\bigg{[}\frac{-3\cos^{2}\theta\sin^{2}\theta+\cos^{4}\theta}{D_{T}+{\rm
i}b_{T}+{\rm i}\beta_{T}\cos\theta}-{\rm
i}\beta_{T}\frac{-\cos^{3}\theta\sin^{2}\theta}{(D_{T}+{\rm i}b_{T}+{\rm
i}\beta_{T}\cos\theta)^{2}}\bigg{]}+O(D_{R})$ $\displaystyle H_{4}$
$\displaystyle=-\frac{c_{y}^{2}}{2}\bigg{[}\frac{-3\cos^{2}\theta\sin^{2}\theta+\sin^{4}\theta}{D_{T}+{\rm
i}b_{T}+{\rm i}\beta_{T}\sin\theta}-{\rm
i}\beta_{T}\frac{-\sin^{3}\theta\cos^{2}\theta}{(D_{T}+{\rm i}b_{T}+{\rm
i}\beta_{T}\sin\theta)^{2}}\bigg{]}+O(D_{R}).$

Again, as in the immotile case, the form of the right hand side expression
(73) leads us to consider the ansatz (45) for $(\Psi_{2},\boldsymbol{u}_{2})$.
Using this ansatz, the left hand side of the $O(\epsilon^{2})$ equation (44)
takes the following form:

$\displaystyle\mathcal{L}[\Psi_{2}]$ $\displaystyle=2{\rm
i}b_{T}\bigg{(}\psi_{2,1}{\rm e}^{{\rm i}(x+y)+2{\rm
i}b_{T}t}A_{x}A_{y}+\psi_{2,3}{\rm e}^{2{\rm i}x+2{\rm
i}b_{T}t}A_{x}^{2}+\psi_{2,4}{\rm e}^{2{\rm i}y+2{\rm
i}b_{T}t}A_{y}^{2}\bigg{)}$ (74) $\displaystyle\quad+{\rm
i}\beta_{T}\bigg{(}(\cos\theta+\sin\theta)\psi_{2,1}{\rm e}^{{\rm
i}(x+y)+2{\rm i}b_{T}t}A_{x}A_{y}+(\cos\theta-\sin\theta)\psi_{2,2}{\rm
e}^{{\rm i}(x-y)}A_{x}\overline{A}_{y}$ $\displaystyle\hskip
56.9055pt+2\cos\theta\psi_{2,3}{\rm e}^{2{\rm i}x+2{\rm
i}b_{T}t}A_{x}^{2}+2\sin\theta\psi_{2,4}{\rm e}^{2{\rm i}y+2{\rm
i}b_{T}t}A_{y}^{2}\bigg{)}$
$\displaystyle\quad-\frac{1}{4\upi}(\cos^{2}\theta-\sin^{2}\theta)\int_{0}^{2\upi}\bigg{(}\psi_{2,1}{\rm
e}^{{\rm i}(x+y)+2{\rm i}b_{T}t}A_{x}A_{y}+\psi_{2,2}{\rm e}^{{\rm
i}(x-y)}A_{x}\overline{A}_{y}\bigg{)}(\cos^{2}\theta-\sin^{2}\theta)d\theta$
$\displaystyle\quad-\frac{1}{\upi}\cos\theta\sin\theta\int_{0}^{2\upi}\bigg{(}\psi_{2,3}{\rm
e}^{2{\rm i}x+2{\rm i}b_{T}t}A_{x}^{2}+\psi_{2,4}{\rm e}^{2{\rm i}y+2{\rm
i}b_{T}t}A_{y}^{2}\bigg{)}\sin\theta\cos\theta d\theta$
$\displaystyle\quad+D_{T}\bigg{(}2\psi_{2,1}{\rm e}^{{\rm i}(x+y)+2{\rm
i}b_{T}t}A_{x}A_{y}+2\psi_{2,2}{\rm e}^{{\rm
i}(x-y)}A_{x}\overline{A}_{y}+4\psi_{2,3}{\rm e}^{2{\rm i}x+2{\rm
i}b_{T}t}A_{x}^{2}+4\psi_{2,4}{\rm e}^{2{\rm i}y+2{\rm
i}b_{T}t}A_{y}^{2}\bigg{)}$ $\displaystyle\quad-
D_{R}\bigg{(}(\p_{\theta\theta}\psi_{2,1}){\rm e}^{{\rm i}(x+y)+2{\rm
i}b_{T}t}A_{x}A_{y}+(\p_{\theta\theta}\psi_{2,2}){\rm e}^{{\rm
i}(x-y)}A_{x}\overline{A}_{y}+(\p_{\theta\theta}\psi_{2,3}){\rm e}^{2{\rm
i}x+2{\rm i}b_{T}t}A_{x}^{2}$ $\displaystyle\hskip
56.9055pt+(\p_{\theta\theta}\psi_{2,4}){\rm e}^{2{\rm i}y+2{\rm
i}b_{T}t}A_{y}^{2}+(\p_{\theta\theta}\psi_{2,5})\left\lvert
A_{x}\right\rvert^{2}+(\p_{\theta\theta}\psi_{2,6})\left\lvert
A_{y}\right\rvert^{2}\bigg{)}+{\rm c.c.}$

Equating exponents from the right hand side (73) with the left hand side (74),
we obtain six independent integrodifferential ODEs for the coefficients
$\psi_{2,j}(\theta)$.

For each of $\psi_{2,j}$, $j=1,2,3,4$, it can be shown that (small) $D_{R}>0$
results in a small perturbation of the $D_{R}=0$ solution to the following set
of equations:

$\displaystyle\big{(}2D_{T}+2{\rm i}b_{T}+{\rm
i}\beta_{T}(\cos\theta+\sin\theta)\big{)}\psi_{2,1}-\frac{\cos^{2}\theta-\sin^{2}\theta}{4\upi}\int_{0}^{2\upi}\psi_{2,1}(\cos^{2}\theta-\sin^{2}\theta)d\theta-
D_{R}\p_{\theta\theta}\psi_{2,1}$ (75)
$\displaystyle\qquad=-\frac{c_{x}c_{y}}{2}\bigg{[}\frac{-3\cos^{2}\theta\sin^{2}\theta+\sin^{4}\theta+\cos\theta\sin\theta}{D_{T}+{\rm
i}b_{T}+{\rm i}\beta_{T}\cos\theta}-{\rm
i}\beta_{T}\frac{\sin^{4}\theta\cos\theta}{(D_{T}+{\rm i}b_{T}+{\rm
i}\beta_{T}\cos\theta)^{2}}$ $\displaystyle\hskip
56.9055pt+\frac{-3\cos^{2}\theta\sin^{2}\theta+\cos^{4}\theta+\cos\theta\sin\theta}{D_{T}+{\rm
i}b_{T}+{\rm i}\beta_{T}\sin\theta}-{\rm
i}\beta_{T}\frac{\cos^{4}\theta\sin\theta}{(D_{T}+{\rm i}b_{T}+{\rm
i}\beta_{T}\sin\theta)^{2}}\bigg{]}+O(D_{R}),$
$\displaystyle\big{(}2D_{T}+{\rm
i}\beta_{T}(\cos\theta-\sin\theta)\big{)}\psi_{2,2}-\frac{\cos^{2}\theta-\sin^{2}\theta}{4\upi}\int_{0}^{2\upi}\psi_{2,2}(\cos^{2}\theta-\sin^{2}\theta)d\theta-
D_{R}\p_{\theta\theta}\psi_{2,2}$
$\displaystyle\qquad=-\frac{c_{x}c_{y}}{2}\bigg{[}\frac{-3\cos^{2}\theta\sin^{2}\theta+\sin^{4}\theta+\cos\theta\sin\theta}{D_{T}+{\rm
i}b_{T}+{\rm i}\beta_{T}\cos\theta}-{\rm
i}\beta_{T}\frac{\sin^{4}\theta\cos\theta}{(D_{T}+{\rm i}b_{T}+{\rm
i}\beta_{T}\cos\theta)^{2}}$ $\displaystyle\hskip
56.9055pt+\frac{-3\cos^{2}\theta\sin^{2}\theta+\cos^{4}\theta+\cos\theta\sin\theta}{D_{T}+{\rm
i}b_{T}+{\rm i}\beta_{T}\sin\theta}-{\rm
i}\beta_{T}\frac{\cos^{4}\theta\sin\theta}{(D_{T}+{\rm i}b_{T}+{\rm
i}\beta_{T}\sin\theta)^{2}}\bigg{]}+O(D_{R}),$
$\displaystyle\big{(}4D_{T}+2{\rm i}b_{T}+2{\rm
i}\beta_{T}\cos\theta\big{)}\psi_{2,3}-\frac{1}{\upi}\cos\theta\sin\theta\int_{0}^{2\upi}\psi_{2,3}\sin\theta\cos\theta\,d\theta-
D_{R}\p_{\theta\theta}\psi_{2,3}$
$\displaystyle\qquad=-\frac{c_{x}^{2}}{2}\bigg{[}\frac{-3\cos^{2}\theta\sin^{2}\theta+\cos^{4}\theta}{D_{T}+{\rm
i}b_{T}+{\rm i}\beta_{T}\cos\theta}-{\rm
i}\beta_{T}\frac{-\cos^{3}\theta\sin^{2}\theta}{(D_{T}+{\rm i}b_{T}+{\rm
i}\beta_{T}\cos\theta)^{2}}\bigg{]}+O(D_{R}),$
$\displaystyle\big{(}4D_{T}+2{\rm i}b_{T}+2{\rm
i}\beta_{T}\sin\theta\big{)}\psi_{2,4}-\frac{1}{\upi}\cos\theta\sin\theta\int_{0}^{2\upi}\psi_{2,4}\sin\theta\cos\theta\,d\theta-
D_{R}\p_{\theta\theta}\psi_{2,4}$
$\displaystyle\qquad=-\frac{c_{y}^{2}}{2}\bigg{[}\frac{-3\cos^{2}\theta\sin^{2}\theta+\sin^{4}\theta}{D_{T}+{\rm
i}b_{T}+{\rm i}\beta_{T}\sin\theta}-{\rm
i}\beta_{T}\frac{-\sin^{3}\theta\cos^{2}\theta}{(D_{T}+{\rm i}b_{T}+{\rm
i}\beta_{T}\sin\theta)^{2}}\bigg{]}+O(D_{R}).$

In particular, following a similar approach to Section 2.4, it may be shown
that the solutions to each of the above four equations is bounded independent
of $D_{R}$ as $D_{R}\to 0$ for all values of the triple
$(D_{T},b_{T},\beta_{T})$ of interest.

In contrast, the remaining two coefficients $\psi_{2,5}$ and $\psi_{2,6}$ blow
up like $1/D_{R}$ as $D_{R}\to 0$:

$\displaystyle\p_{\theta\theta}\psi_{2,5}$
$\displaystyle=\frac{c_{x}^{2}}{2D_{R}}\bigg{[}\frac{-3\cos^{2}\theta\sin^{2}\theta+\cos^{4}\theta}{D_{T}+{\rm
i}b_{T}+{\rm i}\beta_{T}\cos\theta}+{\rm
i}\beta_{T}\frac{\cos^{3}\theta\sin^{2}\theta}{(D_{T}+{\rm i}b_{T}+{\rm
i}\beta_{T}\cos\theta)^{2}}+O(D_{R})\bigg{]}$ (76)
$\displaystyle\p_{\theta\theta}\psi_{2,6}$
$\displaystyle=\frac{c_{y}^{2}}{2D_{R}}\bigg{[}\frac{-3\cos^{2}\theta\sin^{2}\theta+\sin^{4}\theta}{D_{T}+{\rm
i}b_{T}+{\rm i}\beta_{T}\sin\theta}+{\rm
i}\beta_{T}\frac{\sin^{3}\theta\cos^{2}\theta}{(D_{T}+{\rm i}b_{T}+{\rm
i}\beta_{T}\sin\theta)^{2}}+O(D_{R})\bigg{]}.$

In particular, for sufficiently small $D_{R}>0$, the behavior of $\psi_{2,5}$
and $\psi_{2,6}$ dominate over each of the other terms which are quadratic in
the amplitudes $A_{x}$ and $A_{y}$. We may then integrate the equations (76)
in $\theta$ to obtain the expressions for $\psi_{2,5}$ and $\psi_{2,6}$ given
in equation (46).

Moving on to the $O(\epsilon^{3})$ equation (48), we note that out of all the
cubic-in-$A$ terms in the right hand side expression (50), each term is $O(1)$
in $D_{R}$ except for $-{\rm{div}}_{p}\big{(}({\bf
I}-\boldsymbol{p}\boldsymbol{p}^{\rm
T})(\bnabla\boldsymbol{u}_{1}\boldsymbol{p})\Psi_{2}\big{)}$, which is
$O(D_{R}^{-1})$. For sufficiently small $D_{R}$, this term determines the
behavior of the cubic-in-$A$ terms and gives rise to the following expressions
for $R_{xx}(\theta)$, $R_{xy}(\theta)$, $R_{yx}(\theta)$, and
$R_{yy}(\theta)$:

$\displaystyle R_{xx}$
$\displaystyle=\frac{3}{2}a_{1}\cos^{2}\theta\sin\theta+4a_{2}\cos^{3}\theta\sin\theta-
a_{2}\cos\theta\sin\theta+\frac{5}{2}a_{3}\cos^{4}\theta\sin\theta$ (77)
$\displaystyle\quad+\frac{a_{4}}{2}\frac{{\rm
i}\beta_{T}\sin\theta\cos^{2}\theta}{D_{T}+{\rm i}b_{T}+{\rm
i}\beta_{T}\cos\theta}+a_{4}\cos\theta\sin\theta\,L(\cos\theta)+O(D_{R})$
$\displaystyle R_{xy}$
$\displaystyle=\frac{3}{2}a_{1}\sin^{2}\theta\cos\theta-\frac{a_{1}}{2}\cos\theta-4a_{2}\sin\theta\cos^{3}\theta+a_{2}\sin\theta\cos\theta+\frac{5}{2}a_{3}\cos\theta\sin^{4}\theta$
$\displaystyle\quad-\frac{3}{2}a_{3}\cos\theta\sin^{2}\theta-\frac{a_{4}}{2}\frac{{\rm
i}\beta_{T}\cos^{3}\theta}{D_{T}+{\rm i}b_{T}+{\rm
i}\beta_{T}\sin\theta}+a_{4}\sin\theta\cos\theta\,L(\sin\theta)+O(D_{R})$
$\displaystyle R_{yx}$
$\displaystyle=\frac{3}{2}a_{1}\cos^{2}\theta\sin\theta-\frac{a_{1}}{2}\sin\theta-4a_{2}\cos\theta\sin^{3}\theta+a_{2}\cos\theta\sin\theta+\frac{5}{2}a_{3}\cos^{4}\theta\sin\theta$
$\displaystyle\quad-\frac{3}{2}a_{3}\cos^{2}\theta\sin\theta-\frac{a_{4}}{2}\frac{{\rm
i}\beta_{T}\sin^{3}\theta}{D_{T}+{\rm i}b_{T}+{\rm
i}\beta_{T}\cos\theta}+a_{4}\cos\theta\sin\theta\,L(\cos\theta)+O(D_{R})$
$\displaystyle R_{yy}$
$\displaystyle=\frac{3}{2}a_{1}\sin^{2}\theta\cos\theta+4a_{2}\sin^{3}\theta\cos\theta-
a_{2}\sin\theta\cos\theta+\frac{5}{2}a_{3}\sin^{4}\theta\cos\theta$
$\displaystyle\quad+\frac{a_{4}}{2}\frac{{\rm
i}\beta_{T}\cos\theta\sin^{2}\theta}{D_{T}+{\rm i}b_{T}+{\rm
i}\beta_{T}\sin\theta}+a_{4}\sin\theta\cos\theta\,L(\sin\theta)+O(D_{R}).$

Here the constants $a_{j}$ are given in (47), and the notation $L(\cdot)$ is
used to denote

$L(z):=\log(D_{T}+{\rm i}b_{T}+{\rm
i}\beta_{T}z)-\frac{1}{2\upi}\int_{0}^{2\upi}\log(D_{T}+{\rm i}b_{T}+{\rm
i}\beta_{T}z)\,d\theta,\quad z=\sin\theta,\cos\theta.$ (78)

The solvability condition for equation (48) then requires integrating the
expressions in (77) against the eigenmodes $\psi_{x,1}(\theta)$ or
$\psi_{y,1}(\theta)$ to obtain the coefficients $M_{j}$ defined in (53). Most
of the resulting $\theta$-integrals in (53) may be evaluated analytically,
leading to the following expressions for $M_{0}$, $M_{1}$, $M_{2}$, $M_{3}$:

$\displaystyle M_{0}$
$\displaystyle=-\frac{\upi}{\beta_{T}^{4}}\bigg{(}6(D_{T}+{\rm
i}b_{T})^{2}+\beta_{T}^{2}-(6(D_{T}+{\rm
i}b_{T})^{2}+4\beta_{T}^{2})\frac{(D_{T}+{\rm
i}b_{T})}{\sqrt{\beta_{T}^{2}+(D_{T}+{\rm i}b_{T})^{2}}}\bigg{)}$ (79)
$\displaystyle M_{1}$
$\displaystyle=\frac{\upi}{96\beta_{T}^{8}}\bigg{(}5\beta_{T}^{6}+232\beta_{T}^{2}(D_{T}+{\rm
i}b_{T})^{4}+512(D_{T}+{\rm i}b_{T})^{6}-28\beta_{T}^{4}(D_{T}+{\rm
i}b_{T})^{2}$ $\displaystyle\quad-\frac{8(D_{T}+{\rm
i}b_{T})^{3}}{\sqrt{\beta_{T}^{2}+(D_{T}+{\rm
i}b_{T})^{2}}}\big{(}61\beta_{T}^{2}(D_{T}+{\rm i}b_{T})^{2}+64(D_{T}+{\rm
i}b_{T})^{4}+4\beta_{T}^{4}\big{)}\bigg{)}$
$\displaystyle\quad+\frac{(D_{T}+{\rm
i}b_{T})^{3}}{2\beta_{T}^{4}}L_{1}(D_{T},b_{T},\beta_{T})+O(D_{R})$
$\displaystyle M_{2}$ $\displaystyle=\frac{\upi(D_{T}+{\rm
i}b_{T})^{3}}{4\beta_{T}^{8}}\bigg{(}\sqrt{\beta_{T}^{2}+(D_{T}+{\rm
i}b_{T})^{2}}\bigg{(}4(D_{T}+{\rm
i}b_{T})^{2}+\beta_{T}^{2}+\frac{2\beta_{T}^{2}(D_{T}+{\rm
i}b_{T})^{2}}{\beta_{T}^{2}+2(D_{T}+{\rm i}b_{T})^{2}}$
$\displaystyle\quad-4(D_{T}+{\rm i}b_{T})^{3}-4\beta_{T}^{2}(D_{T}+{\rm
i}b_{T})\bigg{)}\bigg{)}+\frac{(D_{T}+{\rm
i}b_{T})^{3}}{2\beta_{T}^{4}}L_{2}(D_{T},b_{T},\beta_{T})+O(D_{R})$
$\displaystyle M_{3}$ $\displaystyle=\frac{2\upi(D_{T}+{\rm
i}b_{T})}{\beta_{T}^{5}}\bigg{(}4\big{(}D_{T}+{\rm
i}b_{T})^{2}+\beta_{T}^{2}-\frac{(4(D_{T}+{\rm
i}b_{T})^{2}+3\beta_{T}^{2}\big{)}(D_{T}+{\rm
i}b_{T})}{\sqrt{\beta_{T}^{2}+(D_{T}+{\rm i}b_{T})^{2}}}\bigg{)},$

where

$\displaystyle L_{1}(D_{T},b_{T},\beta_{T})$
$\displaystyle=\int_{0}^{2\upi}\frac{\cos^{2}\theta\sin^{2}\theta}{D_{T}+{\rm
i}b_{T}+{\rm
i}\beta_{T}\cos\theta}L(\cos\theta)\,d\theta=\int_{0}^{2\upi}\frac{\cos^{2}\theta\sin^{2}\theta}{D_{T}+{\rm
i}b_{T}+{\rm i}\beta_{T}\sin\theta}L(\sin\theta)\,d\theta$ $\displaystyle
L_{2}(D_{T},b_{T},\beta_{T})$
$\displaystyle=\int_{0}^{2\upi}\frac{\sin^{2}\theta\cos^{2}\theta}{D_{T}+{\rm
i}b_{T}+{\rm
i}\beta_{T}\cos\theta}L(\sin\theta)\,d\theta=\int_{0}^{2\upi}\frac{\sin^{2}\theta\cos^{2}\theta}{D_{T}+{\rm
i}b_{T}+{\rm i}\beta_{T}\sin\theta}L(\cos\theta)\,d\theta$

for $L(\cdot)$ as in (78).

## References

  * Albritton & Ohm (2022) Albritton, Dallas & Ohm, Laurel 2022 On the stabilizing effect of swimming in an active suspension. arXiv preprint .
  * Bodenschatz et al. (1988) Bodenschatz, Eberhard, Zimmermann, Walter & Kramer, Lorenz 1988 On electrically driven pattern-forming instabilities in planar nematics. Journal de Physique 49 (11), 1875–1899.
  * Crawford & Knobloch (1991) Crawford, John David & Knobloch, Edgar 1991 Symmetry and symmetry-breaking bifurcations in fluid dynamics. Annual Review of Fluid Mechanics 23 (1), 341–387.
  * Cross & Hohenberg (1993) Cross, Mark C & Hohenberg, Pierre C 1993 Pattern formation outside of equilibrium. Reviews of modern physics 65 (3), 851.
  * DeCamp et al. (2015) DeCamp, Stephen J, Redner, Gabriel S, Baskaran, Aparna, Hagan, Michael F & Dogic, Zvonimir 2015 Orientational order of motile defects in active nematics. Nature materials 14 (11), 1110–1115.
  * Doi (1981) Doi, Masao 1981 Molecular dynamics and rheological properties of concentrated solutions of rodlike polymers in isotropic and liquid crystalline phases. J. Polym. Sci. B 19 (2), 229–243.
  * Doi & Edwards (1986) Doi, Masao & Edwards, Samuel Frederick 1986 The Theory of Polymer Dynamics, , vol. 73. Oxford University Press.
  * Dombrowski et al. (2004) Dombrowski, Christopher, Cisneros, Luis, Chatkaew, Sunita, Goldstein, Raymond E & Kessler, John O 2004 Self-concentration and large-scale coherence in bacterial dynamics. Physical review letters 93 (9), 098103.
  * Doostmohammadi et al. (2017) Doostmohammadi, Amin, Shendruk, Tyler N, Thijssen, Kristian & Yeomans, Julia M 2017 Onset of meso-scale turbulence in active nematics. Nature communications 8 (1), 1–7.
  * Dunkel et al. (2013) Dunkel, Jörn, Heidenreich, Sebastian, Drescher, Knut, Wensink, Henricus H, Bär, Markus & Goldstein, Raymond E 2013 Fluid dynamics of bacterial turbulence. Physical review letters 110 (22), 228102.
  * Ezhilan et al. (2013) Ezhilan, Barath, Shelley, Michael J & Saintillan, David 2013 Instabilities and nonlinear dynamics of concentrated active suspensions. Physics of Fluids 25 (7), 070607.
  * Forest et al. (2014) Forest, M Gregory, Phuworawong, Panon, Wang, Qi & Zhou, Ruhai 2014 Rheological signatures in limit cycle behaviour of dilute, active, polar liquid crystalline polymers in steady shear. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 372 (2029), 20130362.
  * Forest et al. (2004a) Forest, M Gregory, Wang, Qi & Zhou, Ruhai 2004a The flow-phase diagram of doi-hess theory for sheared nematic polymers ii: finite shear rates. Rheologica Acta 44 (1), 80–93.
  * Forest et al. (2004b) Forest, M Gregory, Wang, Qi & Zhou, Ruhai 2004b The weak shear kinetic phase diagram for nematic polymers. Rheologica acta 43 (1), 17–37.
  * Forest et al. (2013) Forest, M Gregory, Wang, Qi & Zhou, Ruhai 2013 Kinetic theory and simulations of active polar liquid crystalline polymers. Soft Matter 9 (21), 5207–5222.
  * Forest et al. (2015) Forest, M Gregory, Wang, Qi & Zhou, Ruhai 2015 Kinetic attractor phase diagrams of active nematic suspensions: the dilute regime. Soft Matter 11 (32), 6393–6402.
  * Gachelin et al. (2014) Gachelin, J, Rousselet, A, Lindner, A & Clement, E 2014 Collective motion in an active suspension of escherichia coli bacteria. New Journal of Physics 16 (2), 025003.
  * Gao et al. (2017) Gao, Tong, Betterton, Meredith D, Jhang, An-Sheng & Shelley, Michael J 2017 Analytical structure, dynamics, and coarse graining of a kinetic model of an active fluid. Physical Review Fluids 2 (9), 093302.
  * Giomi et al. (2011) Giomi, Luca, Mahadevan, Lakshminarayanan, Chakraborty, Bulbul & Hagan, MF 2011 Excitable patterns in active nematics. Physical review letters 106 (21), 218101.
  * Giomi et al. (2012) Giomi, L, Mahadevan, L, Chakraborty, B & Hagan, MF 2012 Banding, excitability and chaos in active nematic suspensions. Nonlinearity 25 (8), 2245.
  * Gompper et al. (2020) Gompper, Gerhard, Winkler, Roland G, Speck, Thomas, Solon, Alexandre, Nardini, Cesare, Peruani, Fernando, Löwen, Hartmut, Golestanian, Ramin, Kaupp, U Benjamin, Alvarez, Luis & others 2020 The 2020 motile active matter roadmap. Journal of Physics: Condensed Matter 32 (19), 193001\.
  * Henkin et al. (2014) Henkin, Gil, DeCamp, Stephen J, Chen, Daniel TN, Sanchez, Tim & Dogic, Zvonimir 2014 Tunable dynamics of microtubule-based active isotropic gels. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 372 (2029), 20140142.
  * Hohenegger & Shelley (2010) Hohenegger, Christel & Shelley, Michael J 2010 Stability of active suspensions. Physical Review E 81 (4), 046311\.
  * Jeffery (1922) Jeffery, George Barker 1922 The motion of ellipsoidal particles immersed in a viscous fluid. Proceedings of the Royal Society of London. Series A, Containing papers of a mathematical and physical character 102 (715), 161–179.
  * Knobloch (1986) Knobloch, E 1986 Oscillatory convection in binary mixtures. Physical Review A 34 (2), 1538.
  * Koch & Subramanian (2011) Koch, Donald L & Subramanian, Ganesh 2011 Collective hydrodynamics of swimming microorganisms: living fluids. Annual Review of Fluid Mechanics 43, 637–659.
  * Lauga & Powers (2009) Lauga, Eric & Powers, Thomas R 2009 The hydrodynamics of swimming microorganisms. Reports on Progress in Physics 72 (9), 096601.
  * Lushi et al. (2014) Lushi, Enkeleida, Wioland, Hugo & Goldstein, Raymond E 2014 Fluid flows created by swimming bacteria drive self-organization in confined suspensions. Proc. Natl. Acad. Sci. 111 (27), 9733–9738.
  * Marchetti et al. (2013) Marchetti, M Cristina, Joanny, Jean-François, Ramaswamy, Sriram, Liverpool, Tanniemola B, Prost, Jacques, Rao, Madan & Simha, R Aditi 2013 Hydrodynamics of soft active matter. Reviews of Modern Physics 85 (3), 1143.
  * Mendelson et al. (1999) Mendelson, Neil H, Bourque, Adrienne, Wilkening, Kathryn, Anderson, Kevin R & Watkins, Joseph C 1999 Organized cell swimming motions in Bacillus subtilis colonies: patterns of short-lived whirls and jets. J. Bacteriol. 181 (2), 600–609.
  * Needleman & Dogic (2017) Needleman, Daniel & Dogic, Zvonimir 2017 Active matter at the interface between materials science and cell biology. Nature reviews materials 2 (9), 1–14.
  * Nishiguchi et al. (2017) Nishiguchi, Daiki, Nagai, Ken H, Chaté, Hugues & Sano, Masaki 2017 Long-range nematic order and anomalous fluctuations in suspensions of swimming filamentous bacteria. Physical Review E 95 (2), 020601.
  * Opathalage et al. (2019) Opathalage, Achini, Norton, Michael M, Juniper, Michael PN, Langeslay, Blake, Aghvami, S Ali, Fraden, Seth & Dogic, Zvonimir 2019 Self-organized dynamics and the transition to turbulence of confined active nematics. Proceedings of the National Academy of Sciences 116 (11), 4788–4797.
  * Pedley & Kessler (1992) Pedley, TJ & Kessler, John O 1992 Hydrodynamic phenomena in suspensions of swimming microorganisms. Annual Review of Fluid Mechanics 24 (1), 313–358.
  * Peng et al. (2021) Peng, Yi, Liu, Zhengyang & Cheng, Xiang 2021 Imaging the emergence of bacterial turbulence: Phase diagram and transition kinetics. Science Advances 7 (17), eabd1240.
  * Pomeau (1986) Pomeau, Yves 1986 Front motion, metastability and subcritical bifurcations in hydrodynamics. Physica D: Nonlinear Phenomena 23 (1-3), 3–11.
  * Saintillan & Shelley (2008a) Saintillan, David & Shelley, Michael J 2008a Instabilities and pattern formation in active particle suspensions: kinetic theory and continuum simulations. Physical Review Letters 100 (17), 178103.
  * Saintillan & Shelley (2008b) Saintillan, David & Shelley, Michael J 2008b Instabilities, pattern formation, and mixing in active suspensions. Physics of Fluids 20 (12), 123304.
  * Saintillan & Shelley (2013) Saintillan, David & Shelley, Michael J 2013 Active suspensions and their nonlinear models. Comptes Rendus Physique 14 (6), 497–517.
  * Saintillan & Shelley (2015) Saintillan, David & Shelley, Michael J 2015 Theory of active suspensions. In Complex Fluids in biological systems, pp. 319–355. Springer.
  * Sanchez et al. (2012) Sanchez, Tim, Chen, Daniel TN, DeCamp, Stephen J, Heymann, Michael & Dogic, Zvonimir 2012 Spontaneous motion in hierarchically assembled active matter. Nature 491 (7424), 431–434.
  * Schöpf & Zimmermann (1993) Schöpf, Wolfgang & Zimmermann, Walter 1993 Convection in binary fluids: Amplitude equations, codimension-2 bifurcation, and thermal fluctuations. Physical Review E 47 (3), 1739.
  * Simha & Ramaswamy (2002) Simha, R Aditi & Ramaswamy, Sriram 2002 Hydrodynamic fluctuations and instabilities in ordered suspensions of self-propelled particles. Physical review letters 89 (5), 058101.
  * Sokolov & Aranson (2012) Sokolov, Andrey & Aranson, Igor S 2012 Physical properties of collective motion in suspensions of bacteria. Physical review letters 109 (24), 248109.
  * Sokolov et al. (2007) Sokolov, Andrey, Aranson, Igor S, Kessler, John O & Goldstein, Raymond E 2007 Concentration dependence of the collective dynamics of swimming bacteria. Physical review letters 98 (15), 158102.
  * Stenhammar et al. (2017) Stenhammar, Joakim, Nardini, Cesare, Nash, Rupert W., Marenduzzo, Davide & Morozov, Alexander 2017 Role of correlations in the collective behavior of microswimmer suspensions. Phys. Rev. Lett. 119, 028005.
  * Subramanian & Koch (2009) Subramanian, Ganesh & Koch, Donald L 2009 Critical bacterial concentration for the onset of collective swimming. Journal of Fluid Mechanics 632, 359–400.