|  Flow states of two dimensional active gels driven by external shear
---|---
|  Wan Luo∗ab, Aparna Baskaranc, Robert A. Pelcovitsde, and Thomas R.
Powersabde†
|  Using a minimal hydrodynamic model, we theoretically and computationally
study active gels in straight and annular two-dimensional channels subject to
an externally imposed shear. The gels are isotropic in the absence of
externally- or activity-driven shear, but have nematic order that increases
with shear rate. Using the finite element method, we determine the possible
flow states for a range of activities and shear rates. Linear stability
analysis of an unconfined gel in a straight channel shows that an externally
imposed shear flow can stabilize an extensile fluid that would be unstable to
spontaneous flow in the absence of the shear flow, and destabilize a
contractile fluid that would be stable against spontaneous flow in the absence
of shear flow. These results are in rough agreement with the stability
boundaries between the base shear flow state and the nonlinear flow states
that we find numerically for a confined active gel. For extensile fluids, we
find three kinds of nonlinear flow states in the range of parameters we study:
unidirectional flows, oscillatory flows, and dancing flows. To highlight the
activity-driven spontaneous component of the nonlinear flows, we characterize
these states by the average volumetric flow rate and the wall stress. For
contractile fluids, we only find the linear shear flow and a nonlinear
unidirectional flow in the range of parameters that we studied. For large
magnitudes of the activity, the unidirectional contractile flow develops a
boundary layer. Our analysis of annular channels shows how curvature of the
streamlines in the base flow affects the transitions among flow states.

##

††footnotetext: ∗ Email<EMAIL_ADDRESS>† Email:
<EMAIL_ADDRESS>a School of Engineering, Brown
University, Providence, RI 02912, USA.††footnotetext: b Center for Fluid
Mechanics, Brown University, Providence, RI 02912, USA.††footnotetext: c
Martin Fisher School of Physics, Brandeis University, Waltham, MA 02453
USA.††footnotetext: d Department of Physics, Brown University, Providence, RI
02912, USA.††footnotetext: e Brown Theoretical Physics Center, Brown
University, Providence, RI 02912, USA.

##  1 Introduction

The defining property of an active fluid is that energy is added to the system
at the small length scales of the particles that make up the fluid, instead of
at the large length scales of the bounding walls or inlets of the system.1
Commonly studied examples include cytoplasm 2 or its reconstituted
components,3, 4, 5 collections of swimming microorganisms,6, 7, 8 and model
two-dimensional layers of cells.9 The interplay of the energy injected at
small scales and the interactions among the constituent particles lead to
nonequilibrium collective behavior, including spontaneous coherent flows,10,
11, 12 sustained oscillations,1, 13 active turbulence,14, 15, 16 and two-
dimensional4 or three-dimensional 17, 18 topological defects in active liquid
crystalline fluids. These phenomena suggest that active fluids may be used for
novel microfluidics applications, including fluids that pump themselves or mix
themselves. Since these applications require a degree of control over active
fluids, recent investigations have studied how confinement of active fluids
affects flows and the formation of defects.19, 20, 21, 13, 22 In this paper,
we build on these investigations by studying the flow states of an active gel
in a channel with moving boundaries to see how an imposed shear affects the
possible flow states and the transitions among them.

By ‘active gel’ we mean a model liquid crystal which tends to the isotropic
phase away from boundaries with strong anchoring conditions and in the absence
of shear flow. The motionless, isotropic state of an unbounded two-dimensional
active gel is unstable to spontaneous flow and nematic ordering above a
critical activity.23, 24 Recent numerical calculations have identified the
spontaneous flow states in straight three-dimensional 25, 26 and two-
dimensional channels 25, 13, 27 with stationary walls. In a two dimensional
channel with no-torque anchoring conditions at the walls, the critical
activity for spontaneous flow increases as the channel width decreases.25
Thus, confinement is stabilizing, as has been found in other related
situations.28 For a given value of the activity parameter, new flow states
emerge as the channel width increases, with the flow progressing through
unidirectional, undulating (also known as ‘oscillatory’,13) and dancing flow
states.29, 25 A similar sequence of flow states is found for fixed channel
width and increasing activity.25

Our work is motivated by the experimental observation that imposed shear can
prevent 26 the spontaneous instability of a solution 12 of microtubule bundles
and kinesin motors in the presence of the molecular fuel ATP. Instead of a
motionless state, our base state is the state of simple shear in which the
flow field is given by the solution to the Stokes equation for our straight or
annular channel geometry. Working at fixed channel width, we find that
increasing the activity leads to a sequence of flow states which are
reminiscent of the ones seen in the case of no external shear, but with some
important new elements. For example, the imposed shear rate can be stabilizing
in the same sense that confinement is stabilizing: for an extensile active
gel, we find that the critical activity for the imposed simple shear flow to
develop a spontaneous flow component increases with the imposed shear rate. A
similar result was established using linear stability analysis of a polar
system by Muhuri, Rao, and Ramaswamy.30 Here we give a more systematic
treatment of this problem for the apolar case, revealing that the imposed
shear also leads to oscillatory behavior in the unstable modes. For a
contractile active gel, we find that shear is destabilizing. Earlier work has
also examined the rheology of active nematics and gels, showing that polar
active particles have a nonmonotonic stress-strain relation at high
activity,31 and illuminating the nature of shear banding in apolar active
gels.32 Our work extends these investigations to the case of an annular
channel, illustrating the role of the curvature of the streamlines of the base
flow.

Our paper begins with a minimal hydrodynamic model for active gels. We then
study the linear stability of an active gel in a straight channel subject to a
uniform shear flow imposed by a moving plate. In the stable region, the linear
rheology, orientational order, and the shear stress exerted by the active
fluids on the moving boundary are analytically calculated for the state of
uniform shear. Then we turn to the other flow states using the finite element
method to characterize the flow transitions for the extensile and contractile
fluids. Next, we turn to an annular channel and carry out similar analytical
and numerical studies to assess the effects of the curvature of the
boundaries.

##  2 Minimal hydrodynamic model

We use a well-studied continuum hydrodynamic model for nematic liquid crystals
33, 34 to describe apolar microtubules, adding a term corresponding to non-
equilibrium active forces as was done in the "minimal" model used by Varghese
et al.25 In two dimensions, the orientational order of apolar active matter is
described by a traceless, symmetric tensor—the tensor order parameter that is
used in the theory of nematic liquid
crystals—$Q_{ij}=S(2n_{i}n_{j}-\delta_{ij})$, with $i,j=x,y$. 35 The unit
vector $\mathbf{n}(\mathbf{x})$ is the director at position $\mathbf{x}$ and
the scalar order parameter $S$ represents the degree of alignment. The
equilibrium state of the microtubule bundles is governed by a Landau-Ginzburg
free energy density,

$\displaystyle\mathcal{F}$
$\displaystyle=\frac{K}{2}\partial_{i}Q_{jk}\partial_{i}Q_{jk}+\frac{A}{2}Q_{ij}Q_{ij}+\frac{C}{4}\left(Q_{ij}Q_{ij}\right)^{2},$
(1)

where repeated indices are summed over. The single Frank elastic constant $K$
penalizes gradients of $Q_{ij}$. Since we focus on a low concentration
isotropic phase, $A$ will be positive to guarantee that the minimizing state
is disordered. In two dimensions there is no term cubic in $Q_{ij}$, and the
isotropic-nematic transition is continuous. In the isotropic phase we consider
in this paper, the term proportional to $C$ can be neglected, as was done in
previous studies of two-dimensional and three-dimensional channel flow.25, 26

A minimal hydrodynamic model for incompressible flow in two dimensions is
given by25

$\displaystyle 0$ $\displaystyle=$
$\displaystyle\boldmath{\nabla}\cdot\mathbf{v}$ (2) $\displaystyle 0$
$\displaystyle=$
$\displaystyle-\bm{\nabla}p+\eta\nabla^{2}\mathbf{v}-a\nabla\cdot\textsf{Q}$
(3) $\displaystyle 0$ $\displaystyle=$
$\displaystyle-\nu\big{(}\partial_{t}\textsf{Q}+\textbf{v}\cdot\nabla\textsf{Q}+\textsf{Q}\cdot\Omega-\Omega\cdot\textsf{Q}\bigr{)}-A\textsf{Q}+K\nabla^{2}\mathsf{Q}$
(4) $\displaystyle+$ $\displaystyle 2\lambda\nu\mathsf{E},$

where $\eta$ is the shear viscosity, $\nu$ is the rotational viscosity, $p$ is
pressure,
$(\mathbf{v}\cdot\gradient{\mathsf{Q}})_{ij}=v_{k}\partial_{k}Q_{ij}$,
$\mathsf{E}=(\nabla\mathbf{v}+(\nabla\mathbf{v})^{\rm{T}})/2$ is the strain
rate tensor,
$\mathsf{\Omega}=(\nabla\mathbf{v}-(\nabla\mathbf{v})^{\rm{T}})/2$ [i.e.
$\Omega_{ij}=(\partial_{j}v_{i}-\partial_{i}v_{j})/2$] is the vorticity
tensor, and $a$ is the strength of the activity. A positive value of $a$
corresponds to extensile particles, and a negative value of $a$ corresponds to
contractile particles. The shape parameter $\lambda$ is positive for prolate
particles and negative for oblate particles; $\lambda=1$ corresponds to
needle-like particles. Note that in three dimensions there will be additional
nonlinear terms proportional to $\lambda$ appearing in eqn (4).

We disregard inertial effects because the Reynolds number of the typical
active flows we study is small. In this minimal hydrodynamic model, passive
backflow effects are neglected and the order parameter field $\mathsf{Q}$ only
affects the flow through the active stress $-a\mathsf{Q}$.

Fig. 1: The base state for the flow field and tensor order parameter field
$Q_{ij}$ of an active gel in a straight channel with the bottom wall moving at
a fixed speed $\dot{\gamma}W$. The double-headed arrows correspond to the
director field $\mathbf{n}$ of the extensile apolar active bundles. The tensor
order parameter field is uniform throughout the channel because the flow is
uniform and because we impose Neumann boundary conditions on $Q_{ij}$.

The active time scale which results from the competition between viscosity and
activity is given by $\eta/|a|$. From the dynamical equation for $\mathsf{Q}$,
eqn (4), it is apparent that the relaxation time $\tau$ for distortions away
from the equilibrium isotropic state is $\tau=\nu/A$. Likewise, $\sqrt{K/A}$
is a correlation length for the liquid crystalline order, which we write in
nondimensional form as $\ell=\sqrt{K/A}/W$, where $W$ is the width of the
straight or annular channel. The factor $\lambda\nu$ characterizes the flow
birefringence of a passive ($a=0$) liquid crystal.36 When weak shear
$\dot{\gamma}\ll 1/\tau$ is applied to a nematic liquid crystal in the
isotropic state, the rods align such that $A\mathsf{Q}\approx
2\lambda\nu\mathsf{E}$, which implies that the scalar order parameter is
proportional to the shear rate: $S\propto\dot{\gamma}\tau$.

##  3 Straight channel: start-up problem and linear stability analysis

Let us begin by reviewing the linear stability analysis of an unbounded two-
dimensional active gel 37. An isotropic ($\mathsf{Q}=0$), motionless
($\mathbf{v}=0$) gel is unstable to shear flow and nematic ordering when the
effective shear viscosity ($\eta_{\mathrm{eff}}\equiv\eta-a\lambda\tau$)
vanishes, which occurs for a critical activity
$a_{\mathrm{c}}=\eta/(\lambda\tau)$.23, 24 The form of the effective shear
viscosity shows that extensile particles tend to reduce the shear viscosity,
whereas contractile particles tend to increase it. In the unstable state of
the unconfined geometry, the pattern of alignment of the bundles follows a
sine wave, appearing like a bent filament, or like the nematic configuration
of bend.35

Next, let us consider an active gel confined to an infinite straight channel
of width $W$ and subject to a steady uniform shear flow
$\mathbf{v}_{0}=\dot{\gamma}(W-y)\hat{\mathbf{x}}$ as shown in Fig. 1. We
assume no-slip boundary conditions on the channel walls for the velocity
field, and Neumann conditions, ($\partial_{i}Q_{jk}=0$) or “zero-torque
conditions" for the order parameter field on the walls. Given the parallel
planar channel walls and zero-torque boundary conditions, the nematic order
parameter is uniform and divergenceless for the imposed uniform shear flow. In
our hydrodynamic model, activity only appears in eqn (3), and thus, when
activity is below the critical value for the instability, the order parameter
field is unaffected by the activity.

Fig. 2: Director angle $\phi$ as a function of time for various shear rates
for the startup problem of the liquid crystal order parameter in the case of
steady simple shear. From top to bottom, the shear rates are
$\dot{\gamma}\tau=5$ (green curve), $\dot{\gamma}\tau=1$ (gold curve), and
$\dot{\gamma}\tau=0.2$ (blue curve).

Before considering the stability of simple shear flow, we solve the startup
problem, assuming an initially stationary isotropic gel with activity below
the critical value (to be deduced below). Since the Reynolds number is assumed
to be small, the flow immediately assumes its steady-state value
$\mathbf{v}_{0}$. But the order parameter field attains its steady-state value
only after a time comparable to the liquid crystal relaxation time $\tau$. 38
Given the boundary conditions on the order parameter, we may assume that
$\mathsf{Q}$ is uniform in space. Since $\mathsf{Q}$ is uniform, the
divergence of the active stress vanishes and the flow remains simple shear as
the order-parameter field evolves. The order parameter equations (eqn (4))
reduce to

$\displaystyle\partial_{t}Q_{xx}$ $\displaystyle=$
$\displaystyle-\frac{1}{\tau}Q_{xx}-\dot{\gamma}Q_{xy}$ (5)
$\displaystyle\partial_{t}Q_{xy}$ $\displaystyle=$
$\displaystyle\dot{\gamma}Q_{xx}-\frac{1}{\tau}Q_{xy}-\lambda\dot{\gamma}.$
(6)

Assuming $\mathsf{Q}(t=0)=0$, we find

$\displaystyle Q_{xx}$ $\displaystyle=$ $\displaystyle
Q^{(0)}_{xx}\left[1-\mathrm{e}^{-t/\tau}\cos\left(\dot{\gamma}t\right)\right]+Q^{(0)}_{xy}\mathrm{e}^{-t/\tau}\sin\left(\dot{\gamma}t\right)$
(7) $\displaystyle Q_{xy}$ $\displaystyle=$ $\displaystyle
Q^{(0)}_{xy}\left[1-\mathrm{e}^{-t/\tau}\cos\left(\dot{\gamma}t\right)\right]-Q^{(0)}_{xx}\mathrm{e}^{-t/\tau}\sin\left(\dot{\gamma}t\right),$
(8)

where the steady-state order parameter tensor $\mathsf{Q_{0}}$ is given by

$\displaystyle Q^{(0)}_{xx}$ $\displaystyle=$
$\displaystyle\frac{\lambda\dot{\gamma}^{2}\tau^{2}}{1+\dot{\gamma}^{2}\tau^{2}},$
(9) $\displaystyle Q^{(0)}_{xy}$ $\displaystyle=$
$\displaystyle-\frac{\lambda\dot{\gamma}\tau}{1+\dot{\gamma}^{2}\tau^{2}}.$
(10)

The order parameter rises to its steady state, with oscillations that become
apparent when the shear rate is greater than the relaxation rate $1/\tau$.
These oscillations are reminiscent of the oscillations observed 39 in the
apparent viscosity during the startup flow of 8CB, a director-tumbling
nematogen.40 In simple shear, the director of a tumbling nematic makes a
complete revolution, like a rod undergoing a Jeffery orbit in shear flow.40 In
our case, as long as $\tau$ is finite, the directors oscillate about their
final steady state. Fig. 2 shows the director angle
$\phi=\arctan[Q_{xy}/(S+Q_{xx})]$ (measured counterclockwise from the
$x$-axis) as a function of time.

The steady-state scalar order parameter and the director angle are given by

$\displaystyle S$ $\displaystyle=$
$\displaystyle\frac{\lambda\dot{\gamma}\tau}{\sqrt{1+\dot{\gamma}^{2}\tau^{2}}}$
(11) $\displaystyle\phi$ $\displaystyle=$
$\displaystyle-\arctan\left(\frac{1}{\sqrt{1+\dot{\gamma}^{2}\tau^{2}}+\dot{\gamma}\tau}\right).$
(12)

Equations (11) and (12) show that in steady state, the flow aligns the nematic
director at a nonzero angle with the horizontal streamlines, with a degree of
order that increases with increasing shear rate. At low shear rates,
$\dot{\gamma}\tau\ll 1$, the bundles are oriented at an angle of $\phi=-\pi/4$
with the streamlines, and the order is weak ($S\ll 1)$. At high shear rates,
the bundles tend to align parallel to the streamlines, and $S\approx\lambda$.
For needle-like particles, with $\lambda\approx 1$, the order is strong in the
limit of high shear rate. The shear stress on the moving plate in the stable
region is

$\sigma_{\rm{W}}=-\eta\dot{\gamma}-aQ^{(0)}_{xy}=\dot{\gamma}\left(-\eta+\frac{a\lambda\tau}{1+\dot{\gamma}^{2}\tau^{2}}\right),\quad
a<a_{c}.$ (13)

From eqn (13), it is easy to see the wall shear stress increases linearly with
activity but the dependence on the imposed shear is not linear when the
activity is below the critical value.

To analyze the stability of the base configuration with flow rate
$\mathbf{v}_{0}$ and the confinement $W$ , we consider a perturbation that is
independent of $x$, the channel axis. ***A more general assumption would be to
suppose the perturbation depends on both $x$ and $y$, but here we forbid
$x$-dependence to simplify the analysis. The more general analysis using
pseudospectral methods will be reported elsewhere. Thus,
$\mathbf{v}=\mathbf{v}_{0}+\mathbf{v}_{1}$ and
$\mathsf{Q}=\mathsf{Q}_{0}+\mathsf{Q}_{1}$ , with the perturbations

$\displaystyle\mathbf{v}_{1}$ $\displaystyle=$ $\displaystyle v_{x}\sin(n\pi
y/W)\exp\left(\beta t\right)\hat{\mathbf{x}},$ (14)
$\displaystyle\mathsf{Q}_{1}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}\mathcal{Q}_{xx}&\mathcal{Q}_{xy}\\\
\mathcal{Q}_{xy}&-\mathcal{Q}_{xx}\end{pmatrix}\cos(n\pi y/W)\exp\left(\beta
t\right),$ (15)

where $v_{x}$, $\mathcal{Q}_{xx}$, and $\mathcal{Q}_{xy}$ are constants, $n$
is a nonzero positive integer, and $\beta$ is the growth rate of the
perturbation. With these assumptions, the $x$ component of the force equation
eqn (3) implies

${v}_{x}=\frac{a\mathcal{Q}_{xy}W}{n\pi\eta}.$ (16)

Using eqn (16) in the linearized equations for $\mathsf{Q}_{1}$ yields

$\displaystyle\beta_{\pm}$ $\displaystyle=$
$\displaystyle-\frac{1}{\tau}\left(1+\frac{\pi^{2}Kn^{2}}{AW^{2}}\right)+\frac{\lambda
a}{2\eta(1+\dot{\gamma}^{2}\tau^{2})}$ (17) $\displaystyle\pm$
$\displaystyle\sqrt{\left[\frac{\lambda
a}{2\eta(1+\dot{\gamma}^{2}\tau^{2})}\right]^{2}-\dot{\gamma}^{2}\left(1+\frac{\lambda
a\tau/\eta}{1+\dot{\gamma}^{2}\tau^{2}}\right)}.$

There are two modes. In the limit of a passive fluid, $a=0$, the modes
collapse to a single mode corresponding to oscillations of the order parameter
as it decays to its equilibrium value given by eqn (11):
$\beta_{\pm}=-\lambda/\tau[1+\pi^{2}K/(AW^{2})]\pm\mathrm{i}\dot{\gamma}$.
Note the similarity between these damped oscillations and the damped
oscillations in the startup problem, eqns (7) and (8). A nonzero activity
makes the two modes distinct. In the limit of zero shear rate, $\beta_{-}$ is
negative and independent of activity even if $a\neq 0$, and corresponds to the
decay of the scalar order parameter of a passive isotropic nematic when it is
perturbed from the isotropic value $S=0$. The other mode corresponds to the
spontaneous flow and ordering of an active isotropic nematic when
$a>a_{\mathrm{c}}=[1+\pi^{2}K/(AW^{2})]\eta/(\lambda\tau)$. Note that the
confining channel walls raise the critical activity above the previously
quoted critical value for unbounded space. The elastic constant $K$ only
enters the growth rate if the channel width is finite.

Fig. 3: Linear stability analysis results for a two-dimensional active gel in
a straight channel of width $W$ subject to a shear flow with rate
$\dot{\gamma}$. The Frank elasticity is small: $K=0.01AW^{2}$. Simple shear
flow is stable against perturbations in the shaded blue region, and the
perturbations are oscillatory in the region between two dashed lines.

In general, the critical activity for instability depends on the shear rate,
and is found by determining when $\real(\beta_{+})=0$ for $n=1$. The modes are
oscillatory when the square root in eqn (17) is imaginary, or when
$a_{-}<a<a_{+}$, where

$a_{\pm}=\frac{2\eta\dot{\gamma}}{\lambda}(1+\dot{\gamma}^{2}\tau^{2})\left(\dot{\gamma}\tau\pm\sqrt{1+\dot{\gamma}^{2}\tau^{2}}\right).$
(18)

When $a_{-}<a<a_{+}$, the critical curve $\mathrm{Re}[\beta_{+}(n=1)]=0$ in
the $\dot{\gamma}$-$a$ plane is given by

$a_{\mathrm{1c}}=2\frac{\eta}{\lambda\tau}\left(1+\dot{\gamma}^{2}\tau^{2}\right)\left(1+\pi^{2}\ell^{2}\right),$
(19)

where $\ell$ is the dimensionless correlation length defined in the previous
section. When $a<a_{-}$ or $a>a_{+}$, the growth rate is purely real, and the
critical curve $\beta_{+}(n=1)=0$ is given by

$a_{2\rm{c}}=\frac{\eta}{\tau}\frac{\left(1+\dot{\gamma}^{2}\tau^{2}\right)\left[(1+\pi^{2}\ell^{2})^{2}+\dot{\gamma}^{2}\tau^{2}\right]}{\lambda(1+\pi^{2}\ell^{2}-\dot{\gamma}^{2}\tau^{2})}$
(20)

Note that $a_{2\mathrm{c}}>0$ for $\sqrt{1+\pi^{2}\ell^{2}}>\dot{\gamma}\tau$,
and $a_{2\mathrm{c}}<0$ for $\sqrt{1+\pi^{2}\ell^{2}}<\dot{\gamma}\tau$.

The stability boundaries are plotted in Fig. 3 for the case of $\ell=0$ (zero
Frank elasticity). The region of oscillatory growth rates, $a_{-}<a<a_{+}$, is
the region between the dashed lines. The stable region is the shaded blue
region between the solid blue curves, whereas the unstable regions are the
white regions. Note that the upper stability boundary is given by
$a_{1\mathrm{c}}$ in the oscillatory region, and $a_{2\mathrm{c}}$ in the non-
oscillatory region. The lower stability boundary lies wholly in the non-
oscillatory region, and is therefore given by $a_{2\mathrm{c}}$. Since the
upper stability boundary near $\dot{\gamma}=0$ increases with shear rate, our
results are in agreement with Muhuri et al.,30 who found that shear
counteracts the instability for extensile particles. Surprisingly, we also
find that shear can be destablilizing for contractile active particles if the
magnitude of the activity is large enough.

##  4 Straight channel: nonlinear spontaneous flows

The linear analysis of the previous section predicts that simple shear flow
with uniform nematic order is stable as long as the activity and externally
imposed shear rate lie in the shaded region of Fig. 3. However, there may be
transitions to flow states that are not captured by linear stability analysis,
and furthermore, the linear equations cannot describe the fully-developed flow
states. Thus, we explore the activity-induced flow states and the transitions
between them by numerically solving the full nonlinear equations, eqns
(2)–(4). We use the open source finite element software FEniCS41, 42, 43 to
solve the nonlinear equations, employing a backwards Euler scheme to solve for
the time dependence. We characterize the flow states by the spontaneous
volumetric flow rate as well as the wall shear stress.

The system is initialized with a small value of the nematic order parameter
$S$, appropriate for an isotropic state. For sufficiently small values of the
external shear, the direction of the activity-induced flow for $a>a_{c}$
depends on the configuration of the nematic order. We can achieve positive
flow—flow in the same direction the bottom wall moves—or negative flow—flow
against the direction the bottom wall moves—by imposing appropriate initial
conditions on the directors. These conditions will be described below for the
extensile and contractile cases. The initial director fields also have small
random fluctuations. Because we are neglecting inertial effects, we do not
need to initialize the velocity field, which is determined from eqns (2)-(4).
Instead of attempting to simulate a very long channel, we use periodic
boundary conditions on the left and right boundaries of the channel. The
length $L$ of the channel is chosen to be five times the width $W$; we found
this length to be the longest channel length we could simulate in a reasonable
amount of computing time. We focus on situations in which the width $W$ of the
channel is large compared to the correlation length $\sqrt{K/A}$ of the liquid
crystal. Therefore, our simulations are carried out with a small value of the
Frank elasticity, $K/A=0.01W^{2}$ (i.e. $\ell=0.1$). In our numerical
calculations, $W$ is the unit of length, $\tau$ is the unit of time, and
$\eta/\tau$ is the unit of pressure. We also define the dimensionless activity
$\alpha=a{\lambda}\tau/\eta$, and restrict our simulations to the case of
needle-like particles, $\lambda=1$.

Fig. 4: Numerically determined flow states for extensile fluids with
$\ell=0.1$ and $\lambda=1$ in a straight channel. See the text for an
explanation of how the flow states were determined. The blue shaded region and
the dashed line indicate the linearly stable states and the boundary for
oscillatory modes, respectively, for $\ell=0.1$ (compare with Fig. 3). The
location of the transitions is generally insensitive to whether the initial
splay of the nematic directors converges to the right (as in Fig. 5b) or the
left (as in Fig. 5d). Note that depending on the noise in the initial
conditions, some of states at higher $\alpha$ can either be dancing or
oscillatory-like states; an example is shown in the SI.  Fig. 5: Steady states
of spontaneous flow and nematic order for dimensionless activity $\alpha=1.5$
and dimensionless shear rate $\dot{\gamma}\tau=0.2$ (the bottom wall moves to
the right, and the top wall is stationary). In the velocity field plots
(panels (a) and (c)), color denotes flow speed, and arrows denote the
direction of the flow. In the order parameter field (panels (b) and (d)),
color denotes the scalar order parameter and lines denote the director field.
Note that the leftward spontaneous flow in (c) is not sufficient to reverse
the net flow near the moving wall, but leads to negative net flow near the
stationary wall.

###  4.1 Extensile fluids

For extensile fluids, we find three types of flow states when the activity is
above the critical value $a_{c}$: unidirectional, oscillatory, and dancing.
These states are similar to three of the states found by Samui et al., 13 who
studied an active nematic fluid confined to a channel in the absence of
external shear. These authors also found an active turbulent state at high
activity, which we do not explore here. The unidirectional flow is steady,
consisting of a superposition of spontaneous flow and simple shear flow. The
oscillatory flow is unsteady, with a pattern of flow and order that translates
at a constant velocity along the channel, which makes the spatially-averaged
wall stress constant in time. The dancing flow is truly unsteady, with a
spatially-averaged wall stress that oscillates in time. These states will be
described in more detail below. Fig. 4 shows the phase diagram for flow states
for dimensionless activity in the range $0\leq\alpha\lesssim 2.5$ and shear
rate in the range $0\leq\dot{\gamma}\tau\lesssim 1$. To get positive
spontaneous flow, we imposed initial conditions with the directors converging
to the right, as in Fig. 5b. To get negative spontaneous flow, we imposed
initial conditions with the directors converging to the left, as in Fig. 5d.
We ran each simulation until either all transients died out, or $t=600\tau$,
whichever came first. The final state could either be a steady state or a
state with regular periodic behaviour. Then we classified the states as
follows. The simple shear and unidirectional flow states generally emerge at
times $t<600\tau$. Both states are steady with negligible $y$-component of
velocity, and these two flow states are easily distinguished since simple
shear has the standard linear flow profile $v_{x}=\dot{\gamma}(W-y)$, whereas
unidirectional flow has a spontaneous flow component added to the linear flow.
If there is a nonzero $y$-component of the velocity at the end of the
simulation, we check for oscillations in the average wall stress,
$\bar{\sigma}_{\mathrm{w}}=\int_{0}^{L}\mathrm{d}x\sigma_{xy}(x,y=0)/L$, for
times in the range $550\tau$–$600\tau$. Negligible oscillation in the average
wall stress implies an oscillatory flow state, while non-negligible values
imply a dancing state. Most of the points shown in Fig. 4 reached a steady or
regular periodic state by $t=600\tau$, or came very close to doing so. But a
few cases near transitions between flow states needed much longer to fully
develop.

The limit of stability for the simple shear flow states in the Fig. 4 is the
boundary between the region with black crosses and the regions with other
symbols. We observe that the numerical limit of stability for simple shear
flow matches very well with the prediction of linear stability analysis
(filled blue region), but only for the transition from the simple shear to
unidirectional flow, $\dot{\gamma}\tau\lesssim 0.3$. The disagreement between
the linear stability boundary and the transition from simple shear flow to
oscillatory flow may be due to our neglect of the possibility that the
perturbation could depend on $x$ as well as $y$. In the region of simple shear
flow (black crosses in Fig. 4), our numerical results show that the wall
stress decreases with activity, in agreement with eqn 13. Fig. 6 shows the
numerically computed wall stress, normalized by the passive (viscous) stress.
When the flow state is simple shear, activity reduces the total wall stress in
proportion to the activity, in accord with the general understanding that
extensile particles with activity reduces the effective viscosity. 37

Fig. 6: The spatially-averaged wall stress on the bottom wall of a straight
channel, normalized by the passive wall stress. Because the average stress
depends somewhat on the spatial period of the oscillatory and dancing flow
states, only the results for one representative period are shown. For dancing
flow states, the spatially-averaged wall stress oscillates in time, and the
square symbols show the time-averaged value of the wall shear stress.  Fig. 7:
Spontaneous component of velocity profiles in dimensionless units for
unidirectional flow at different rates of externally imposed shear. For each
value of the shear rate, there are two branches, with the upper branch
corresponding to the positive spontaneous flow, and the lower branch
corresponding to the negative spontaneous flow.

Unidirectional flow. When the externally imposed shear is in the range
$0\leq\dot{\gamma}\tau\lesssim 0.3$, and the dimensionless activity is in a
relatively narrow band near $\alpha\approx 1$ (Fig. 4), activity creates a
steady unidirectional flow along the $x$-axis (Fig. 5). The activity-induced
component spontaneously breaks the left-right symmetry of the channel, with
the actual direction of the active flow component determined not by the
imposed external shear but instead by the initial conditions of the directors,
as described above. Since the total shear rate vanishes at the value of $y$ at
which the flow rate has an extremum, the scalar parameter vanishes at this
same value of $y$ (Fig. 5). Fig. 7 shows the flow profile subtracting off the
imposed shear flow for fixed activity and various values of $\dot{\gamma}$ for
both the left-moving and right-moving spontaneous flows. It indicates that the
spontaneous active component of the flow depends on $\dot{\gamma}$; in other
words, the total flow is not simply a superposition of the passive shear flow
$v_{x}=\dot{\gamma}(W-y)$ and the spontaneous flow at _zero_ externally
imposed shear.

Fig. 8: Dimensionless activity-driven volumetric flow rate versus
dimensionless activity in a straight channel. The symbols denote the flow
states and the colors denote the externally imposed shear rate. The wavelength
is $5W/4$ for all cases with nonzero activity-driven flow rate shown in the
figure.  Fig. 9: Examples of the unsteady oscillatory flow state at
$t=600\tau$. Panels (a), (b), (c) and (d) show the flow patterns and order
parameter fields corresponding to positive and negative spontaneous flows for
small external shear rate. Panels (e) and (f) show a case with larger shear
rate. At this shear rate, we only observe positive spontaneous flow.

To better characterize these flows, we subtract the passive volumetric flow
rate from the total volumetric flow rate to get the dimensionless activity-
induced volumetric flow rate (per unit channel width),

$q^{\mathrm{active}}\equiv\left(\int_{0}^{W}\mathrm{d}yv_{x}-\frac{\dot{\gamma}W^{2}}{2}\right)\frac{\tau}{W^{2}},$
(21)

shown in Fig. 8. This quantity serves as an order parameter describing the
transitions among the various flow states. Fig. 8 shows that the activity-
driven flow rate has the same magnitude for the left-moving and right-moving
flows, and also that the amplitude of the unidirectional flows increases as
the activity increases.

Examining Fig. 6 for the case of $\dot{\gamma}\tau=0.2$ reveals that reduction
of the normalized wall stress with increasing activity ceases at the onset of
the unidirectional flow, and the normalized wall stress at $y=0$ starts to
increase slightly as activity increases further. The active component of the
wall stress at $y=0$ in the unidirectional flow has the opposite sign compared
to that of the simple shear flow, as can be seen from the opposite orientation
of the directors near the wall $y=0$ in Fig. 5b and Fig. 1. Also, the active
component of the flow changes the sign of the flow gradient near the wall, as
can be seen from Fig. 7. These two effects together lead to the rise in the
normalized wall stress at the onset of unidirectional flow.

Oscillatory flow. Our phase diagram of flow states shows that for
$\dot{\gamma}\tau\lesssim 0.3$, there is a transition with increasing activity
from the unidirectional flow states to two-dimensional oscillatory flows (Fig.
9). When $\dot{\gamma}\tau\gtrsim 0.3$, the simple shear states transition
directly to two-dimensional oscillatory flows as activity increases. Although
the oscillatory flow states are unsteady, with the velocity and order
parameter taking the form of a traveling wave, the flow pattern and order
parameter configuration rigidly translate in the $x$ direction with wave speed
$v_{\mathrm{T}}$. In other words, in the frame moving relative to the channel
walls with speed $v_{\mathrm{T}},$ the streamlines meander in space but are
steady. Likewise, the configuration of the order parameter tensor is steady in
this frame. Because we use periodic boundary conditions, the flow field and
orientational order parameter must have period in $x$ equal to the channel
length $L$. But these fields could also have a shorter period, which must
evenly divide the total channel length. Since we use a channel length $L=5W$,
the possible wavelengths for a periodic configuration are $5W$, $5W/2$,
$5W/3,$ …. Different wavelengths are selected in the dynamical final state
depending on the initial state of the nematic directors, as well as the value
of the activity and the imposed shear. Because it is difficult to determine
the relationship between the random fluctuations imposed on the initial
directors and the wavelength that is finally selected, we did not make a
systematic study of all the possible wavelengths. It is natural to worry that
the steady translation of the flow field and order parameter pattern could be
an artifact of the periodic boundary conditions. In Sec. 5, we study an
annular geometry as a single domain without the need for periodic boundary
conditions. Since we also observe an oscillatory flow state with constant
angular wave speed in that situation, we are confident the constant wave speed
$v_{T}$ we see in the straight channel is not an artifact of the period
boundary conditions.

We measured the volumetric flux for times in the range $t=550$–$600\tau$,
which is when the system is generally in its final dynamically stable state.
In the final state, the volumetric flow rate and wall shear stress of the
oscillatory flows are constant. For small externally imposed shear (e.g.
$\dot{\gamma}\tau\leq 0.2$ in Fig. 4), the spontaneous activity-induced flows
can be either positive or negative, depending on the form of the splay in the
initial conditions for nematic order, as for the unidirectional flows. See
Figs. 9a–d (movies are in the SI). For positive spontaneous flow, the
streamlines undulate, but the externally imposed shear breaks the up-down
symmetry of the waves with respect to the horizontal centerline of the
channel. The velocity at the valleys of the waves is higher than at the peaks.
For negative spontaneous flow, since the activity-induced flow is opposite to
the direction of the externally imposed shear flow, the flow more easily forms
circular streamlines. Thus, for $\dot{\gamma}\neq 0$, the absolute value of
the activity-driven flux of negative spontaneous flows is slightly smaller
than the flux for the positive spontaneous flows, as can be seen by looking
very closely at Fig. 8.

The direction of the spontaneous flow not only determines the shape of the
streamlines, but also determines the direction of translation of the total
flow pattern, including the passive viscous flow. For positive spontaneous
flow, the total flow pattern translates in the $+x$ direction, while for the
negative spontaneous flow case, it translates in the $-x$ direction. The
activity-driven volumetric flow rate is nonzero but generally decreases with
increasing activity as shown in Fig. 8. Fig. 10 shows that the wave speed
$v_{\rm{T}}$ pattern is faster than the flux, and the difference between these
two quantities decreases with the growth of the activity.

We now turn to larger externally imposed shear (e.g. $\dot{\gamma}\tau\geq
0.3$ in Fig. 4). In this case, only the positive spontaneous flow appears; the
symmetry is broken by the flow imposed by the external shear. The activity-
driven volumetric flow rate is zero because the imposed shear rate is large
enough to close the streamlines. Interestingly, our numerical results indicate
that the wave speed is equal to the average volumetric flow rate of simple
shear, $v_{\mathrm{T}}=\dot{\gamma}W/2$.

Fig. 10: Wave speed $v_{T}$ of the flow pattern relative to the average
dimensionless flow rate for the oscillatory flow state with various (small)
externally imposed shear rates.

Dancing flow. At higher activity, the flow field and tensor order parameter
field become unsteady in any frame, and we find states (Fig. 11; movies are in
the SI) analogous to the dancing flows found by Shendruk et al. 29 and Samui
et al. 13 in their study of active nematic flow in a two-dimensional channel.
The volumetric flow rate of dancing flow is still constant with time.
Additionally, in the range we study ($\alpha\leq 2.5)$, when activity is large
enough to dynamically close all streamlines for the part of the flow that is
activity-driven, the total flux is the same as in the passive case. As in the
case of the oscillatory flows, sometimes we find multiple states at the same
values of parameters. For example, noise in the initial conditions may cause
the system to exhibit oscillatory-like states in the region of the phase
diagram where dancing flows are also found.

Given a director configuration
$\hat{\mathbf{n}}=\cos\phi\hat{\mathbf{x}}+\sin\phi\hat{\mathbf{y}}$, we may
define the topological charge inside a closed loop by computing
$\int\mathrm{d}\phi=2\pi m$ around the loop, where $m$ is the charge. Applying
this definition to the configuration in Fig. 11b may be problematic because
the order parameter $S$ vanishes not just in small cores but in extended two-
dimensional regions. If the loop drawn to encircle a potential topological
defect crosses a region where $S$ vanishes, the angle $\phi$ and the
topological charge are ill-defined. Nevertheless, we can simply look at the
director configuration of dancing flow and see that there are parts of the
configuration around the regions of small $S$ near the center of the channel
that closely approximate the director field of $+1/2$ defects. The $+1/2$
defects appear in pairs, and the two defect cores move with undulations of the
flow in opposite directions leading to the pairs exchanging partners with the
pair to the immediate left and immediate right, consistent with the Ceilidh
dance observed by Shendruk 29 and Samui 13.

Fig. 11: Snapshot of a representative unsteady dancing flow state at
$t={599.8}\tau$. In the velocity field plots, colors denote flow speed, and
arrows denote flow direction. In the order parameter field, colors denote the
scalar order parameter, and lines denote directors.  Fig. 12: Flow states for
contractile particles in a straight channel with $\ell=0.1$. As in Fig. 4, the
blue shaded region is linearly stable, and the modes of the linearized
equations are damped but oscillatory above the dashed line. The finite element
results are insensitive to the initial conditions of the director field.  Fig.
13: An example of the velocity and order parameter fields of positive and
negative spontaneous flow for contractile particles. In the velocity field
plots (a) and (c), color denotes the flow speed, and arrows denote flow
direction. In the order parameter field plots (b) and (d), color denotes the
scalar order parameter, and lines denote the directors.

The spatially averaged shear stress imposed by the active flow on the moving
wall also oscillates in time. The average wall shear stress no longer
decreases linearly with activity in the spontaneous flow region.

###  4.2 Contractile fluids

Negative activity corresponds to contractile particles. When the activity is
sufficiently negative and the shear rate is large enough,
$\dot{\gamma}\tau>\sqrt{1+\pi^{2}\ell^{2}}$ , we observe unidirectional flow
states in our finite-element calculations. The stability boundary that we find
in our numerical calculations is consistent with the results of our linear
stability analysis (Fig. 12). As in the extensile case, we get both positive
and negative flows, depending on whether the initial configuration of the
directors bends downward as in Fig. 13b, or upward as in Fig. 13d. After
transients have died out, the active component of the volumetric flow rate is
equal in magnitude for the positive and negative flows, and the amplitude of
the flow rate increases as the magnitude of the activity increases. It is
well-known that contractile elongated particles in a shear flow enhance the
shear viscosity. 37 Thus, the wall stress (normalized by passive stress)
increases linearly with the magnitude of the activity when the flow is simple
shear, according to eqn 13. When the flow transitions to unidirectional flow,
we also find that the normalized wall stress increases linearly with the
magnitude of the activity, however with a slightly smaller absolute value of
slope. The figures showing the dependence of the active component of flow rate
and the dependence of the normalized wall stress on activity are in the SI.

When the magnitude of the activity becomes large, we observe a boundary layer
in the flow velocity. Since we found only steady-state unidirectional flow
states for contractile activity, it is computationally more efficient to
reduce the governing partial differential equations to ordinary differential
equations [see eqns (23–26) below] and solve them using the bvp5c solver of
MATLAB.44 Fig. 14 shows the active component of the flow for the positive and
negative spontaneous flows of contractile gels. When the absolute value of the
activity is large, we observe that the spontaneous component of the flow
approaches simple shear flow, with a boundary layer of dimensionless thickness
$\ell_{\delta}$ near one of the walls, which we define as the displacement
boundary layer thickness 45

$\ell_{\delta}\equiv\frac{\int_{0}^{W}\mathrm{d}y\left(\dot{\gamma}_{0}y-(v_{x}-v_{x}^{\rm{passive}})\right)}{\int_{0}^{W}\mathrm{d}y\dot{\gamma}_{0}y},$
(22)

where $\dot{\gamma}_{0}=\mathrm{d}(v_{x}-v_{x}^{\rm{passive}})/{\mathrm{d}y}$
at $y=0$ for positive spontaneous flow. The boundary layer thickness is the
same for positive and negative spontaneous flow. Fig. 14 shows that the peak
flow speed of the active component is higher and the boundary layer is thinner
for larger magnitudes of the activity. From Fig. 15, we find that
$\ell_{\delta}\propto|\alpha-\alpha_{c}|^{\zeta}$, where $\zeta$ is close to
$-0.5$, but its magnitude increases with $\dot{\gamma}\tau$. This dependence
will be studied in another publication.

Fig. 14: Active component of the velocity profile in dimensionless units for
contractile particles in the unidirectional flow state, for
$\dot{\gamma}\tau=2.0$ and various activities. The upper branches correspond
to positive spontaneous flow, and the lower branches correspond to negative
spontaneous flow. The boundary layer thickness is denoted by $\ell_{\delta}$.
Fig. 15: Log-log plot of dimensionless boundary layer thickness
$\ell_{\delta}$ vs. $\alpha-\alpha_{c}$, for various shear rates (legend).
These results indicate that $\ell_{\delta}\propto|\alpha-\alpha_{c}|^{\zeta}$,
where $\zeta$ gradually changes from $-0.44$ to $-0.56$ with increasing
$\dot{\gamma}\tau$.

###  4.3 Weakly nonlinear analysis for $\dot{\gamma}=0$.

To conclude this section, we turn to a weakly nonlinear analysis of the
spontaneous steady unidirectional flow near the transition from the motionless
isotropic state.46 We continue to assume $\lambda=1$ and only consider the
case of zero shear rate, $\dot{\gamma}=0$, leaving the case of nonzero
$\dot{\gamma}$ for another publication. Assuming that the velocity field,
order parameter tensor, and pressure depend only on the coordinate $y$, the
dimensionless governing equations are

$\displaystyle v_{x}^{\prime\prime}-\alpha Q_{xy}^{\prime}$ $\displaystyle=$
$\displaystyle 0$ (23) $\displaystyle-p^{\prime}+\alpha Q_{xx}^{\prime}$
$\displaystyle=$ $\displaystyle 0$ (24)
$\displaystyle\ell^{2}Q_{xx}^{\prime\prime}-Q_{xx}+v_{x}^{\prime}Q_{xy}$
$\displaystyle=$ $\displaystyle 0$ (25)
$\displaystyle\ell^{2}Q_{xy}^{\prime\prime}-Q_{xy}-v_{x}^{\prime}Q_{xx}+v_{x}^{\prime}$
$\displaystyle=$ $\displaystyle 0$ (26)

with no-slip boundary conditions $v_{x}(0)=v_{x}(1)=0$ and no-torque (Neumann)
boundary conditions $Q_{ij}^{\prime}(0)=Q_{ij}^{\prime}(1)=0$. The prime
denotes a derivative with respect to $y$. We already saw in Sec. 3 that the
motionless, distorted state at zero imposed shear rate is unstable when
$\alpha>\alpha_{\mathrm{c}}$, where $\alpha_{c}=(1+\pi^{2}\ell^{2})$ is the
dimensionless critical activity. Here we study the spontaneous flow and weak
ordering for $\alpha=\alpha_{\mathrm{c}}+\delta\alpha$, with $\delta\alpha>0$.
Assuming the balance $Q_{xx}\approx v_{x}^{\prime}Q_{xy}$ in eqn (25) suggests
that to leading order, $v_{x}=\mathcal{O}(\delta\alpha^{1/2})$,
$Q_{xy}=\mathcal{O}(\delta\alpha^{1/2})$, and
$Q_{xx}=\mathcal{O}(\delta\alpha)$. Thus, we expand in powers of
$\delta\alpha^{1/2}$:

$\displaystyle v_{x}$ $\displaystyle=$
$\displaystyle\delta\alpha^{1/2}v_{x}^{(1)}+\delta\alpha
v_{x}^{(2)}+\delta\alpha^{3/2}v_{x}^{(3)}+\dots$ (27) $\displaystyle Q_{ij}$
$\displaystyle=$ $\displaystyle\delta\alpha^{1/2}Q_{ij}^{(1)}+\delta\alpha
Q_{ij}^{(2)}+\delta\alpha^{3/2}Q_{ij}^{(3)}+\dots.$ (28)

At $\mathcal{O}(\delta\alpha^{1/2})$, we find the steady versions of the
linearized equations we used in Sec. 3 to solve for the growth rate,

$\displaystyle v_{x}^{(1)\prime\prime}-\alpha_{\mathrm{c}}Q^{(1)\prime}_{xy}$
$\displaystyle=$ $\displaystyle 0$ (29)
$\displaystyle\ell^{2}Q^{(1)\prime\prime}_{xx}-Q^{(1)}_{xx}$ $\displaystyle=$
$\displaystyle 0$ (30)
$\displaystyle\ell^{2}Q^{(1)\prime\prime}_{xy}-Q^{(1)}_{xy}+v_{x}^{(1)\prime}$
$\displaystyle=$ $\displaystyle 0.$ (31)

The Neumann boundary conditions on $Q_{ij}$ together with eqn (30) imply that
$Q^{(1)}_{xx}(y)=0$. Integrating eqn (29) yields
$v_{x}^{(1)\prime}-\alpha_{\mathrm{c}}Q^{(1)}_{xy}=\sigma^{(1)}$, where
$\sigma^{(1)}$ is a constant. Eliminating $v_{x}^{(1)}$ from eqn (31) leads to

$\ell^{2}Q_{xy}^{(1)\prime\prime}+(\alpha_{\mathrm{c}}-1)Q^{(1)}_{xy}=-\sigma^{(1)}.$
(32)

To get a solution for $Q^{(1)}_{xy}$ that satisfies the Neumann boundary
conditions, we must have

$\displaystyle Q^{(1)}_{xy}$ $\displaystyle=$ $\displaystyle
c_{1}\cos(\sqrt{\alpha_{\mathrm{c}}-1}y/\ell)-\frac{\sigma^{(1)}}{\alpha_{\mathrm{c}}-1}$
(33) $\displaystyle=$ $\displaystyle c_{1}\cos\pi
y-\frac{\sigma^{(1)}}{\alpha_{\mathrm{c}}-1},$ (34)

Using eqn (29) and the no-slip boundary conditions implies $\sigma^{(1)}=0$
and $v^{(1)}=(c_{1}\alpha_{\mathrm{c}}/\pi)\sin\pi y$. Note that to leading
order, $v$ and $Q_{xy}$ are $\mathcal{O}(\delta\alpha^{1/2})$, but $Q_{xx}$ is
at most $\mathcal{O}(\delta\alpha).$ At the next order, the equations are

$\displaystyle v_{x}^{(2)\prime\prime}-(1+\pi^{2}\ell^{2})Q_{xy}^{(2)\prime}$
$\displaystyle=$ $\displaystyle 0$ (35)
$\displaystyle-\ell^{2}Q_{xy}^{(2)\prime\prime}+Q_{xy}^{(2)}-v_{x}^{(2)\prime}$
$\displaystyle=$ $\displaystyle 0$ (36)
$\displaystyle-\ell^{2}Q_{xx}^{(2)\prime\prime}+Q_{xx}^{(2)}$ $\displaystyle=$
$\displaystyle c_{1}^{2}{(1+\pi^{2}\ell^{2})}\cos^{2}\pi y,$ (37)

with solutions

$\displaystyle Q_{xx}^{(2)}$ $\displaystyle=$ $\displaystyle
c_{1}^{2}\frac{1+\pi^{2}\ell^{2}}{2}\left(1+\frac{\cos 2\pi
y}{1+4\pi^{2}\ell^{2}}\right)$ (38) $\displaystyle Q_{xy}^{(2)}$
$\displaystyle=$ $\displaystyle c_{2}\cos\pi y$ (39) $\displaystyle
v_{x}^{(2)}$ $\displaystyle=$ $\displaystyle
c_{2}\frac{1+\pi^{2}\ell^{2}}{\pi}\sin\pi y,$ (40)

where $c_{2}$ is a constant.

To determine $c_{1}$, we must expand to $\mathcal{O}(\delta\alpha^{3/2})$:

$\displaystyle
v_{x}^{(3)\prime\prime}-{(1+\pi^{2}\ell^{2})}Q^{(3)\prime}_{xy}$
$\displaystyle=$ $\displaystyle-c_{1}\pi\sin(\pi y)$ (41)
$\displaystyle\ell^{2}Q_{xy}^{(3)\prime\prime}-Q_{xy}^{(3)}+v^{(3)\prime}$
$\displaystyle=$ $\displaystyle
c_{1}^{3}C_{0}\left[\left(\frac{3}{2}+4\pi^{2}\ell^{2}\right)\cos\pi y\right.$
(42) $\displaystyle+$ $\displaystyle\left.\frac{1}{2}\cos 3\pi y\right],$

where $C_{0}=(1+\pi^{2}\ell^{2})^{2}/[2(1+4\pi^{2}\ell^{2})]$. Integrating eqn
(41) yields

$v_{x}^{(3)\prime}={(1+\pi^{2}\ell^{2})}Q^{(3)}_{xy}+c_{1}\cos\pi
y+\sigma^{(3)},$ (43)

where the constant $\sigma^{(3)}$ appears in the expansion of the stress,
$\sigma=v_{x}^{\prime}-\alpha
Q_{xy}=\delta\alpha^{1/2}\sigma^{(1)}+\delta\alpha\sigma^{(2)}+\delta\alpha^{3/2}\sigma^{(3)}+\dots$.
The solutions we have already found at lower order imply that
$\sigma^{(1)}=\sigma^{(2)}=0$. The no-slip boundary conditions on
$v_{x}^{(3)}$ also imply that $\sigma^{(3)}=0$. Thus, the stress vanishes not
only at the critical value of the activity, but also as $\alpha$ is increased
above $\alpha_{\mathrm{c}}$. Our numerical computations give the same result
just above the critical activity. Using eqn (43) to eliminate $v^{(3)}$ from
eqn (42) yields

$\ell^{2}Q_{xy}^{(3)\prime\prime}+\pi^{2}\ell^{2}Q_{xy}^{(3)}=C_{1}\cos\pi
y+C_{2}\cos 3\pi y,$ (44)

where

$\displaystyle C_{1}$ $\displaystyle=$
$\displaystyle\frac{c_{1}^{3}(1+\pi^{2}\ell^{2})^{2}(3+8\pi^{2}\ell^{2})}{4(1+4\pi^{2}\ell^{2})}-c_{1}$
(45) $\displaystyle C_{2}$ $\displaystyle=$
$\displaystyle\frac{c_{1}^{3}(1+\pi^{2}\ell^{2})^{2}}{4(1+4\pi^{2}\ell^{2})}.$
(46)

To find $c_{1}$, we use the Fredholm alternative, 47 which implies that the
right-hand side of eqn (44) must be orthogonal to the solution of the
corresponding homogeneous equation. Thus, $C_{1}=0$, and

$\displaystyle Q_{xx}$ $\displaystyle=$
$\displaystyle\frac{2\delta\alpha}{\alpha_{\mathrm{c}}}\frac{(1+4\pi^{2}\ell^{2})}{3+8\pi^{2}\ell^{2}}\left(1+\frac{\cos
2\pi y}{1+4\pi\ell^{2}}\right)+\mathcal{O}(\delta\alpha^{3/2})$ (47)
$\displaystyle Q_{xy}$ $\displaystyle=$
$\displaystyle\pm\frac{2\delta\alpha^{1/2}}{\alpha_{\mathrm{c}}}\sqrt{\frac{(1+4\pi^{2}\ell^{2})}{3+8\pi^{2}\ell^{2}}}\cos\pi
y+\mathcal{O}(\delta\alpha)$ (48) $\displaystyle v_{x}$ $\displaystyle=$
$\displaystyle\pm\frac{2\delta\alpha^{1/2}}{\pi}\sqrt{\frac{(1+4\pi^{2}\ell^{2})}{3+8\pi^{2}\ell^{2}}}\sin\pi
y+\mathcal{O}(\delta\alpha),$ (49)

where the two signs for $v_{x}$ and $Q_{xy}$ correspond to the two different
spontaneous directions of flow, and the corresponding orientation of the
directors. These analytical solutions agree well with our numerical solutions
for the spontaneous unidirectional flow state with activity just above the
critical activity.

##  5 Annular channel: nonlinear spontaneous flows

In our work on the straight channel, we saw that simple shear flow led to a
spatially uniform order parameter field $\mathsf{Q}$ when the activity is less
than a critical value. Uniform $\mathsf{Q}$ leads to zero active force on the
fluid. In contrast, if the shear rate in the flow is spatially nonuniform, the
alignment and degree of ordering of the directors will also be spatially
nonuniform, leading to an active force. This situation arises in the case of
curved boundaries—as in an annular channel—for any nonzero value of the
activity, no matter how small. Previous theoretical studies involving curved
boundaries have focused on the case of motionless walls. For example,
Woodhouse and Goldstein found spontaneous circular flow in a circular
chamber,10 and Norton et al. showed that the nature of topological defects in
the director field is determined by the flow rather than the director
anchoring conditions at the wall of a circular chamber.20

Fig. 16: Examples of spontaneous flow states of an active gel in an annular
channel. In these examples, the inner boundary rotates with dimensionless
speed $\omega\tau=0.2$. The Couette-like 1 and 2 states are steady. In the
oscillatory flow state, the flow pattern and order parameter configuration
rotate at a steady rate. The dancing flow state is unsteady. Videos of the
oscillatory and the dancing flow states can be found in the SI, sections S2.1
and S2.2 respectively. In the velocity field plots (top panels), colors denote
flow speed and arrows denote flow direction. In the order parameter field
plots (bottom panels), colors denote the scalar order parameter and lines
denote the directors.

In this section, we introduce curvature by considering the flow states of a
two-dimensional active gel in the Taylor-Couette geometry of an annular
channel between two concentric circular boundaries of radius $R$ and $R+W$. We
impose external shear by rotating the inner boundary with steady angular
frequency $\omega$, leaving the outer boundary stationary. Stokes flow in this
geometry, known as Couette flow, is given by 48

$v_{\theta}=\frac{\omega R^{2}}{(2R+W)W}\left[\frac{(R+W)^{2}}{r}-r\right],$
(50)

where $r$ is the radial polar coordinate. The second term of eqn (50)
corresponds to rigid body rotation and does not lead to any strain rate, but
the first term leads to a nonuniform strain rate, and thus induces a
nonuniform order parameter field and an active force on the fluid for any
nonzero value of the activity. To study the nonlinear flow states of active
flows in the annular channel, we again employ the finite element software
FEniCS to solve the the full nonlinear equations, eqns (2)-(4). We set
$\ell=0.1$, $\lambda=1$ and $R/W=1$.

###  5.1 Extensile fluids

We begin our discussion of the flow states in the annulus with extensile
active gels, $\alpha>0$. As in the case of the straight channel, we give the
initial director field some splay to induce counterclockwise or clockwise
spontaneous flow, with the flow direction depending on the sense of the splay.
For example, splay with the rods converging as we move counterclockwise around
the annulus (Fig. 16d) leads to counterclockwise active flow (Fig. 16c). For
the activities we used, we find the same kinds of active flow states as in the
straight channel: Couette-like states which have no radial component of flow
and are the analogs of the unidirectional states in the straight channel
(Figs. 16a–d), oscillatory states (Figs. 16e and f), and dancing states (Figs.
16g and h). We run the simulations until $t=600\tau$, and characterize the
flow states as we did in the case of the straight channel (Sec. 4). For the
Couette-like flows, we distinguish two flow states by checking whether the
maximum velocity is at the moving wall or in the interior of the annulus. If
the flow is fastest on the wall, we label it a “Couette-like 1” flow state;
otherwise the label is “Couette-like 2”. If the transverse component of the
velocity $v_{r}$ is nonnegligible, we check whether the torque exerted by the
total flow on the inner boundary oscillates during the time interval
$550\tau$-$600\tau$. If it oscillates, then the state is dancing, otherwise it
is oscillatory. There are a few flow states near transitions that need a
longer time to equilibrate. We also find multiple solutions for particular
values of $\omega\tau$ and $\alpha$ for the oscillatory and dancing flows.
Fig. 17 shows flow transitions in the annular channel in range of
$0\leq\alpha\lesssim 2.5$ and $0\leq\dot{\gamma}\tau\lesssim 1$. The
transition from Couette-like to oscillatory flow states is relatively robust,
with the transition states showing little dependence on the initial
conditions. However, comparing with the case of the straight channel, the
states observed in the transition from oscillatory to dancing flow are more
sensitive to the choice of initial conditions.

Fig. 17:  Flow states for extensile fluids in an annular channel with
$\ell=0.1$ and $R/W=1$. The Couette states at $\omega=0$ are states of zero
flow and zero order.

In the case of a straight channel, our numerical calculations always yielded
the Newtonian simple shear state solution as long as the magnitude of the
activity was small enough. The situation is different for the annular channel:
our numerical calculations only yield the Newtonian Couette flow state
solution (eqn (50)) when the activity vanishes. As emphasized earlier, any
nonzero value of activity leads to active force and an active component of the
flow because the order parameter field is nonuniform for nonzero wall rotation
speed $\omega$. Green, Toner and Vitelli examined a similar phenomena for
active nematics in which a surface of nonvanishing Gaussian curvature
generates a spontaneous flow at arbitrarily low values of the activity
parameter.49 As long as $\omega\tau$ is sufficiently small, the flow profile
varies continuously between the Couette, Couette-like 1, and Couette-like 2
states as the activity increases (Fig. 18). Note that the flow velocity
increases with activity for a given imposed rotation rate, as expected because
extensile activity reduces the effective shear viscosity.37 Also, the change
from the Newtonian Couette flow profile is small as long as the activity is
modest, $\alpha\lesssim 0.9$ (Fig. 18), which we examine in Sec. 6.

Fig. 18:  Velocity profiles in dimensionless units for the extensile case of
Couette ($\alpha=0$), Couette-like 1 ($\alpha=0.3$, $0.6$ and $0.9$) and
Couette-like 2 ($\alpha=1.2$ and $1.5$) flow states, for $\omega\tau=0.2$.

The oscillatory flow in an annular channel (Figs. 16e and f; movies are in SI)
is similar to the oscillatory flow in a straight channel. The flow and order
parameter patterns are steady in a frame that rotates at constant speed, and
the average volumetric flow rate ($\int\mathrm{d}rv_{\theta}/W$) is constant
in time. Since we solve the equations in the annular domain without applying
periodic boundary conditions, and still see steady rotation of the flow
pattern and order parameter pattern, we can be confident that the constant
wave speed we saw in the case of the oscillatory flows in the straight channel
is not an artifact of the periodic boundary condition.

In the dancing state, the flow and order parameter patterns periodically
change in time, similar to the case of the straight channel. Unlike the
straight channel, the volumetric flow rate of the dancing flow state (Figs.
16g and h; movies are in SI) in the annular case is not constant in time. This
time dependence arises because the difference in curvatures of the inner and
outer boundaries of the annulus breaks the reflection symmetry of the
boundaries of the straight channel that relates the dancing flow at the top
wall to the dancing flow at the bottom wall. Also, as in the straight channel,
we observe moving pairs of $+1/2$ defect-like patterns with an exchange of
partners in the annular dancing flow. In the straight channel, the defect
pairs are mirror images of each other (see Fig 11b), but in the annulus, the
different curvatures of the two boundaries spoils this symmetry. Joshi et al.
also found similar oscillatory and dancing flow states for active nematics by
changing the curvatures of the annular channel without external shear. 50

Fig. 19:  Active component of the volumetric flow rate for extensile fluids in
dimensionless units of $v_{\theta}\tau/W$ as a function of activity in the
annular channel. The symbols denote the flow states and the colors denote the
externally imposed dimensionless shear rate. For the dancing flows, the square
symbols show the average value of the oscillation of the volumetric flow rate.

Fig. 19 shows the active component of the average flow rate (defined as before
as the average flow rate of the total flow minus the average flow rate of the
$a=0$ case) for the various flow states we studied in the annular channel. For
the case of zero applied shear ($\omega\tau=0$), there are positive and
negative spontaneous flows when the activity exceeds a critical value. But for
$\omega\tau\neq 0$, the flow rate has no bifurcation: it continuously
increases from zero as the activity increases from zero. Another striking
difference with the straight channel is that for nonzero rotation rates of the
inner curved wall, we only observe positive spontaneous flows (Fig. 19), even
when we attempt to reverse the direction of flow by altering the initial
conditions of the directors. This rectification arises because in the curved
channel, the non-uniform alignment of the directors arising from the applied
shear leads to spontaneous flow with the same rotation sense as the rotating
wall. Furthermore, since the wave translation direction corresponds to the
direction of the spontaneous component of the flow, the oscillatory flow
patterns all translate in the $+\theta$ direction when $\omega\tau\neq 0$.
Another difference from the straight channel case is that the active
contribution to the average flow rate does not disappear in the annular
channel for larger shear rate.

Fig. 20 shows the torque exerted by the active fluid on the inner boundary,
normalized by the wall torque in the passive case. The relation of the wall
torque to the activity is very similar to the relation of the wall stress to
the activity in the straight channel case, i.e. the normalized wall torque
decreases with increasing activity for the Couette-like 1 flow state.

The change in slope in the active-flow rate vs. $\alpha$ curve in Fig. 19, or
the normalized wall torque vs. $\alpha$ in Fig. 20 indicates the transition
from the Couette-like flow state to the oscillatory flow state. As noted
earlier, sometimes our numerical approach finds oscillatory patterns of
different wavelengths for the same values of the parameters, which would
likely result in values of the volumetric flow rate and wall torque different
from those shown in Fig. 19 and Fig. 20. Some of the variation in the
normalized torque in the oscillatory and dancing flow regimes in Fig. 20
arises from abrupt changes in wavelength as $\alpha$ was varied.

We compare the wall torque and wall stress of annular and straight channels in
Fig. 21 to show the effect of curvature on the wall stress as a function of
external shear in the range of $0<\alpha\leq 1$. The normalized wall torque
and wall stress are close to each other for small external shear rate and both
increase with external shear rate, but the increase is larger in the annular
channel, i.e. normalized wall torque is closer in value to the passive case.
Thus, the curvature of the channel reduces the effect of activity on the wall
with increasing external shear.

Fig. 20:  The torque imposed by the active flow of an extensile fluid on the
rotating disk normalized by the passive torque. For the dancing flows, the
square symbols denote the average value of the oscillation of the torque.
Fig. 21:  Comparison of the normalized torque on the inner boundary of the
annulus as a function of the dimensionless frequency of rotation of the inner
disk and the normalized wall stress on the bottom wall of the straight channel
as a function of the dimensionless shear rate.  Fig. 22:  Velocity profiles in
dimensionless units for Couette flow ($\alpha=0$) and the contractile Couette-
like 1 ($\alpha<0$) flows states, for $\omega\tau=2.4$.

###  5.2 Contractile fluids

We studied contractile active fluids in a two-dimensional annulus with the
parameters in the range $-16\leq\alpha<0$ and $0<\omega\tau\leq 2.4$. When
$\alpha<0$, we only found Couette-like states with no radial component of the
flow. Since contractile activity is effectively shear thickening,37 the effect
of the activity is always to reduce the flow relative to Newtonian Couette
flow (Fig. 22). As in the extensile case, the direction that the active
component of the flow travels around the annulus is independent of the initial
conditions, but unlike the extensile case, the active component of flow is
negative (against the direction imposed by the externally applied shear). The
magnitude of the negative flow is always less than the magnitude of the
externally imposed Couette flow; therefore, the total flow never reverses. In
this sense, the contractile annular flows are Couette-like 1 states rather
than Couette-like 2 states. In accord with the effective shear-thickenening of
contractile active fluids, the total torque on the inner boundary is always
greater than the hydrodynamic torque in Couette flow (SI Fig. S4).

Fig. 23: Analytical results for flow velocity of the Couette-like 1 state in
an annular channel with weak order for $\ell=0.1$, $\lambda=1$, and various
values of $R/W$ for (a) a contractile fluid with $\alpha=-0.9$ and (b) an
extensile fluid with $\alpha=0.99$. The dashed lines in the panel (a) show the
results for Newtonian flows ($\alpha=0$) for comparison.  Fig. 24: Dependence
of normalized wall torque on $R/W$ for different $\alpha$ from linear analysis
at low shear with $\ell=0.1$ and $\lambda=1$. The colors denote different
activities.

##  6 Annular channel: Linear analysis of curvature at low shear rate

In this section we study the limit in which the flow in the annular channel is
slow enough that the induced order is small, $S\ll 1$. For slow enough flow,
it is valid to neglect the nonlinear terms in eqn (4). On the one hand, this
analysis offers a theoretical explanation of some of the observations in Sec.
5; on the other hand, it gives some insight into the role of the curvature of
the boundaries, which we did not vary in the previous section. For
convenience, here we restate the modified Stokes equation (eqn (3)) in
dimensionless form, along with the dimensionless form of the steady linearized
equation for $\mathsf{Q}$:

$\displaystyle 0$ $\displaystyle=$
$\displaystyle-\bm{\nabla}p+\nabla^{2}\mathbf{v}-{\frac{\alpha}{\lambda}}\nabla\cdot\textsf{Q}$
(51) $\displaystyle 0$ $\displaystyle=$
$\displaystyle-\textsf{Q}+\ell^{2}\nabla^{2}\mathsf{Q}+2\lambda\mathsf{E},$
(52)

As in our numerical calculations, we use the width $W$ of the channel as the
unit of length. Since we seek to study the Couette-like flow state, we assume
$\mathbf{v}=v_{\theta}(r)\hat{\bm{\theta}}$. Note that this flow is
incompressible. We also suppose that $p$, $Q_{rr}$, and $Q_{r\theta}$ are
functions of radius only. With these assumptions, the $rr$ component of eqn
(52) is homogeneous, which together with the Neumann boundary conditions
$\partial_{r}Q_{rr}=0$ at $r=R/W$ and $r=R/W+1$ implies $Q_{rr}=0$. Since the
radial component of the modified Stokes equation(eqn (51) with $Q_{rr}=0$
implies that the pressure gradient vanishes, we take $p=0$.

To solve for the velocity and order parameter fields, we take the divergence
of eqn (52) and combine with eqn (51) with $p=0$ to find

$\nabla^{2}\left(\nabla^{2}-\frac{1}{\xi^{2}}\right)\mathbf{v}=0,$ (53)

where $\xi^{2}={\ell^{2}}/(1-\alpha)$. To focus our attention on the Couette-
like states only, we restrict our analysis to $\alpha<1$ in this section.
Thus,

$\displaystyle v_{\theta}=c_{1}r+c_{2}/r+c_{3}I_{1}(r/\xi)+c_{4}K_{1}(r/\xi),$
(54)

where the $c_{i}$ are constants to be determined, and $I_{1}(x)$ and
$K_{1}(x)$ are modified Bessel functions. Inserting the velocity field eqn
(54) into the $r\theta$ component of eqn (52),

$0=\ell^{2}\left(Q_{r\theta}^{\prime\prime}+\frac{1}{r}Q_{r\theta}^{\prime}-\frac{4}{r^{2}}Q_{r\theta}\right)-Q_{r\theta}+\lambda\left(v_{\theta}^{\prime}-\frac{v_{\theta}}{r}\right),$
(55)

yields

$\displaystyle\ell^{2}\left(Q_{r\theta}^{\prime\prime}+\frac{1}{r}Q_{r\theta}^{\prime}-\frac{4}{r^{2}}Q_{r\theta}\right)-Q_{r\theta}$
(56) $\displaystyle=$
$\displaystyle\frac{\lambda}{\xi}\left[\frac{2c_{2}\xi}{r^{2}}-c_{3}I_{2}(r/\xi)+c_{4}K_{2}(r/\xi)\right],$

which has general solution

$\displaystyle Q_{r\theta}$ $\displaystyle=$ $\displaystyle
c_{5}I_{2}(r/\ell)+c_{6}K_{2}(r/\ell)-2c_{2}\lambda/r^{2}$ (57)
$\displaystyle-$ $\displaystyle
c_{3}\frac{\lambda\xi}{\ell^{2}-\xi^{2}}I_{2}(r/\xi)+c_{4}\frac{\lambda\xi}{\ell^{2}-\xi^{2}}K_{2}(r/\xi).$

Inserting this solution into the modified Stokes equation [eqn (51)] shows
that $c_{5}=c_{6}=0$. The rest of the integration constants are determined by
the no-slip boundary conditions on the (dimensionless) velocity,
$v_{\theta}(R/W)=\omega\tau R/W$ and $v_{\theta}(R/W+1)=0$, and the Neumann
boundary conditions on the order parameter field
$Q_{rr}^{\prime}(R/W)=Q_{rr}^{\prime}(R/W+1)=0$. The complete formulas are too
complicated to display, but we plot the velocity in Fig. 23 for various ratios
of $R/W$ for a representative contractile case (top panel) and extensile case
(bottom panel). In both cases, the flow velocity approaches a linear profile
as $R/W$ becomes large, as expected, since in that limit the curvature of the
annulus becomes unimportant, and the flow approaches simple shear flow. For
the contractile case, Fig. 23a, the velocity profile is close to the Newtonian
result, with the agreement between the two cases getting better as $R/W$
increases. For the extensile case, the velocity curves for different values of
$R/W$ get closer to each other as $\alpha$ increases, becoming very close to
the linear profile around $\alpha=0.885$. Above this value of activity, the
order of the curves reverses, with the linear curve lying below all the other
curves. When $\alpha$ gets very close to unity and $R/W$ is small, the maximum
velocity is not at the wall, i.e. the flow continuously changes from the
Couette-like 1 state to the Couette-like 2 state [Fig. 23b]. Fig. 24 shows the
total torque $M=2\pi R^{2}\sigma_{r\theta}$ on the circle $r=R$ as a function
of $R/W$. Note that the limit of a straight channel is almost obtained when
$R$ becomes comparable to $W$. The torque for a contractile fluid is higher
than the passive value since contractile fluids effectively increase the shear
viscosity. Likewise, the torque for an extensile fluid is less than the
passive value since extensile fluids are shear thinning. The torque approaches
the passive value when $R\ll W$. Note that since we use $W$ as the unit of
length, the limit $R\ll W$ corresponds to making the inner cylinder of
vanishing thickness. When $R<\ell$, the term $\ell^{2}\nabla^{2}\mathsf{Q}$
dominates eqn (52), and therefore $\mathsf{Q}\rightarrow 0$. In this limit,
the active force vanishes, and flow is Couette flow.

It is informative to find the velocity and the order parameter field in the
limit $R\gg W$, where the curvature of the annulus is small. Rather than
taking the limit of the formulas used to make Figs. 23 and 24, it is simplest
to solve the equations directly using regular perturbation theory in powers of
$W/R$. Reinstating the dimensions and writing $r=R+y$, we find

$\displaystyle v_{\theta}$ $\displaystyle=$ $\displaystyle\omega
R\left(1-\frac{y}{W}\right)-\frac{\omega y}{2}\left(1-\frac{y}{W}\right)$ (58)
$\displaystyle+$ $\displaystyle\frac{2\alpha\ell^{2}\omega
W}{1-\alpha}\left[\frac{1-\cosh\left[{(1-2y/W)/}{\xi}\right]}{\cosh[{1}/(2\xi)]}\right]$
$\displaystyle Q_{r\theta}$ $\displaystyle=$
$\displaystyle-\frac{\lambda\omega\tau R}{W}$ (59) $\displaystyle+$
$\displaystyle\frac{\lambda\omega\tau}{2}\left[4y/W-3+4\xi\frac{\sinh\left[{(1-2y/W)/}{\xi}\right]}{\cosh[{1}/(2\xi)]}\right].$

The first terms of eqns (58) and (59) correspond to the velocity and order
parameter field, respectively, of a straight channel with an infinitesimal
imposed shear rate $\dot{\gamma}=\omega R/W$. The remaining terms are the
corrections due to the nonzero curvature of the annular channel. Unlike our
weakly analysis of the active flow in the straight channel (Sec. 4.3), which
had spontaneous flow in either direction, here we see that the component of
flow driven by the activity has a definite sign, and is the same direction as
externally imposed flow for extensile fluids.

##  7 Summary

We investigated the stability and flow states of the active gel confined in a
channel subject to a external shear. An externally imposed shear flow can
stabilize an extensile fluid that would be unstable to spontaneous flow when
there is no external shear flow, and destabilize a contractile fluid that
would be stable against spontaneous flow when there is no external shear flow.
In accordance with previous simulations 13, 25 carried out in the absence of
external shear, we find three kinds of nonlinear flow states in the range of
parameters we study: unidirectional flows, oscillatory flows, and dancing
flows for extensile fluids. The unidirectional flow observed in the straight
channel can have a spontaneous active component which is either positive—in
the same direction as the moving wall—or negative—in the opposite direction of
the moving wall. The oscillatory flow states also have two possible directions
for the spontaneous active component when the externally imposed shear rate is
small. For greater imposed shear rates, the spontaneous flow direction will be
the same as the moving wall. For contractile gels, we only observe
unidirectional flow states in the range of parameters that we studied. These
unidirectional flows can have positive or negative spontaneous active
components. In the analysis of the the wall stress caused by the active flow
on the moving boundary, the extensile flow helps the motion of the moving
boundary, while the contractile flow resists the motion. Moreover, the
external shear flow can weaken this effect of activity on the motion.

Our analysis of the curvature shows there are three main differences between
the flows states for the straight channel and the annular channel. First in
the annular channel, there is no critical activity for the system to be stable
against the spontaneous flow given a nonzero external shear. Second, we only
observe one direction of spontaneous flow: positive for extensile gels, but
negative for contractile gels. Last, the average volumetric flow rate of the
annular case oscillates with time for the dancing flow state, while it is
steady in the straight channel. Also, we find increasing the curvature of the
streamlines weakens the dependence of the wall stress on activity.

Our work suggests several directions for future study. An obvious extension is
to work in three dimensions, allowing both the directors and velocity vectors
to point out of the plane and vary in both directions across a channel. Also,
it would be natural to study the effect of aligning flows induced by a
pressure gradient rather than a moving wall, since Poiseuille-like flow may be
easier to study experimentally.

## Acknowledgements

This work was supported in part by the National Science Foundation through
Grant Nos. MRSEC DMR-2011846, CBET-2227361, and PHY-1748958. We are grateful
to Jesse Ault, Kenny Breuer, Guillaume Duclos, Hamid Karani, Jasper Chen,
Alexander Morozov, and Pranay Sampat for helpful discussions. We also thank
the Center for Computation and Visualization (CCV) at Brown university for use
of high performance computing facilities.

## Notes and references

  * Marchetti _et al._ 2013 M. C. Marchetti, J.-F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao and R. A. Simha, _Rev. Mod. Phys._ , 2013, 85, 1143.
  * Goldstein _et al._ 2008 R. E. Goldstein, I. Tuval and J.-W. van de Meent, _Proc. Natl. Acad. Sci., USA_ , 2008, 105, 3663.
  * Nedelec _et al._ 1997 F. J. Nedelec, T. Surrey, A. C. Maggs and S. Leibler, _Nature_ , 1997, 389, 6648.
  * Sanchez _et al._ 2012 T. Sanchez, D. T. Chen, S. J. DeCamp, M. Heymann and Z. Dogic, _Nature_ , 2012, 491, 431–434.
  * Alvarado _et al._ 2017 J. Alvarado, M. Sheinman, A. Sharma, F. C. MacKintosh and G. H. Koenderink, _Soft Matter_ , 2017, 13, 5624–5644.
  * Riedel _et al._ 2005 I. Riedel, K. Kruse and J. Howard, _Science_ , 2005, 309, 300\.
  * Koch and Subramanian 2011 D. L. Koch and G. Subramanian, _Annual Review of Fluid Mechanics_ , 2011, 43, 637–659.
  * Saintillan and Shelley 2013 D. Saintillan and M. J. Shelley, _Comptes Rendus Physique_ , 2013, 14, 497–517.
  * Duclos _et al._ 2016 G. Duclos, C. Erlenkämper, J.-F. Joanny and P. Siberzan, _Nat. Phys._ , 2016, 13, 58.
  * Woodhouse and Goldstein 2012 F. G. Woodhouse and R. E. Goldstein, _Phys. Rev. Lett._ , 2012, 109, 168105.
  * Lushi _et al._ 2014 E. Lushi, H. Wioland and R. E. Goldstein, _Proc. Nat. Acad. Sci. USA_ , 2014, 111, 9733–9738.
  * Wu _et al._ 2017 K.-T. Wu, J. B. Hishamunda, D. T. Chen, S. J. DeCamp, Y.-W. Chang, A. Fernández-Nieves, S. Fraden and Z. Dogic, _Science_ , 2017, 355, eaal1979.
  * Samui _et al._ 2021 A. Samui, J. M. Yeomans and S. P. Thampi, _Soft Matter_ , 2021, 17, 10640–10648.
  * Wensink _et al._ 2012 H. H. Wensink, J. Dunkel, S. Heidenreich, K. Drescher, R. E. Goldstein, H. Löwen and J. M. Yeomans, _Proc. Nat. Acad. Sci. U.S.A._ , 2012, 109, 14308–14313.
  * Dombrowski _et al._ 2004 C. Dombrowski, L. Cisneros, S. Chatkaew, R. E. Goldstein and J. O. Kessler, _Phys. Rev. Lett._ , 2004, 93, 098103.
  * Dunkel _et al._ 2013 J. Dunkel, S. Heidenreich, K. Drescher, H. H. Wensink, M. Bär and R. E. Goldstein, _Phys. Rev. Lett._ , 2013, 110, 228102.
  * Simha and Ramaswamy 2002 R. A. Simha and S. Ramaswamy, _Phys. Rev. Lett._ , 2002, 89, 058101\.
  * Čopar _et al._ 2019 S. Čopar, J. Aplinc, Ž. Kos, S. Žumer and M. Ravnik, _Phys. Rev. X._ , 2019, 9, 031051.
  * Araújo _et al._ 2023 N. A. M. Araújo, L. M. C. Janssen, T. Barois, G. Boffetta, I. Cohen, A. Corbetta, O. Dauchot, M. Dijkstra, W. M. Durham, A. Dussutour, S. Garnier, H. Gelderblom, R. Golestanian, L. Isa, G. H. Koenderink, H. Löwen, R. Metzler, M. Polin, C. P. Royall, A. Šarić, A. Sengupta, C. Sykes, V. Trianni, I. Tuval, N. Vogel, J. M. Yeomans, I. Zuriguel, A. Marin and G. Volpe, _Soft Matter_ , 2023, 19, 1695–1704.
  * Norton _et al._ 2018 M. M. Norton, A. Baskaran, A. Opathalage, B. Langeslay, S. Fraden, A. Baskaran and M. F. Hagan, _Phys. Rev. E_ , 2018, 97, 012702.
  * Opathalage _et al._ 2019 A. Opathalage, M. M. Norton, M. P. Juniper, B. Langeslay, S. A. Aghvami, S. Fraden and Z. Dogic, _Proc. Natl. Acad. Sci. U.S.A._ , 2019, 116, 4788–4797.
  * Hardoüin _et al._ 2019 J. Hardoüin, R. Hughes, A. Doostmohammadi, J. Laurent, T. Lopez-Leon, J. M. Yeomans, J. Ignés-Mullol and F. Sagués, _Commun. Phys._ , 2019, 2, 1–9.
  * Soni _et al._ 2018 H. Soni, R. A. Pelcovits and T. R. Powers, _Phys. Rev. Lett._ , 2018, 121, 178002.
  * Santhosh _et al._ 2020 S. Santhosh, M. R. Nejad, A. Doostmohammadi, J. M. Yeomans and S. P. Thampi, _J. Stat. Phys._ , 2020, 180, 699.
  * Varghese _et al._ 2020 M. Varghese, A. Baskaran, M. F. Hagan and A. Baskaran, _Phys. Rev. Lett._ , 2020, 125, 268003.
  * Chandrakar _et al._ 2020 P. Chandrakar, M. Varghese, S. A. Aghvami, A. Baskaran, Z. Dogic and G. Duclos, _Phy. Rev. Lett._ , 2020, 125, 257801.
  * Chandragiri _et al._ 2020 S. Chandragiri, A. Doostmohammadi, J. M. Yeomans and S. P. Thampi, _Phys. Rev. Lett._ , 2020, 125, 148002.
  * Wioland _et al._ 2016 H. Wioland, E. Lushi and R. E. Goldstein, _New J. Phys._ , 2016, 18, 075002.
  * Shendruk _et al._ 2017 T. N. Shendruk, A. Doostmohammadi, K. Thijssen and J. M. Yeomans, _Soft Matter_ , 2017, 13, 3853–3862.
  * Muhuri _et al._ 2007 S. Muhuri, M. Rao and S. Ramaswamy, _EPL_ , 2007, 78, 48002\.
  * Giomi _et al._ 2010 L. Giomi, T. B. Liverpool and M. C. Marchetti, _Phys. Rev. E_ , 2010, 81, 051908.
  * Fielding _et al._ 2011 S. M. Fielding, D. Marenduzzo and M. E. Cates, _Phys. Rev. E_ , 2011, 83, 041910.
  * Olmsted and Goldbart 1992 P. D. Olmsted and P. M. Goldbart, _Phys. Rev. A._ , 1992, 46, 4966\.
  * Tóth _et al._ 2002 G. Tóth, C. Denniston and J. M. Yeomans, _Phys. Rev. Lett._ , 2002, 88, 105504.
  * de Gennes and Prost 1995 P. G. de Gennes and J. Prost, _The Physics of Liquid Crystals_ , Oxford University Press, Oxford, 2nd edn, 1995.
  * Gennes 1969 P. G. D. Gennes, _Phys. Lett. A_ , 1969, 30A, 454–455.
  * Hatwalne _et al._ 2004 Y. Hatwalne, S. Ramaswamy, M. Rao and R. A. Simha, _Phys. Rev. Lett._ , 2004, 92, 118101.
  * Krieger _et al._ 2015 M. S. Krieger, M. A. Dias and T. R. Powers, _Eur. Phys. J. E_ , 2015, 38, 94.
  * Gu _et al._ 1993 D.-F. Gu, A. M. Jamieson and S. Q. Wang, _J. Rheol._ , 1993, 37, 985\.
  * Larson 1999 R. G. Larson, _The Structure and Rheology of Complex Fluids_ , Oxford University Press, New York, 1999.
  * Logg and Wells 2010 A. Logg and G. N. Wells, _ACM Transactions on Mathematical Software_ , 2010, 37, 1.
  * A. Logg and Hake 2012 G. N. W. A. Logg and J. Hake, _Automated Solution of Differential Equations by the Finite Element Method_ , Springer, 2012, vol. 84, ch. 10.
  * Alnaes _et al._ 2014 M. S. Alnaes, A. Logg, K. Olgaard, M. E. Rognes and G. N. Well, _ACM Transactions on Mathematical Software_ , 2014, 40, 1.
  * Kierzenka and Shampine 2008 J. Kierzenka and L. F. Shampine, _J. Numer. Anal. Ind. Appl. Math._ , 2008, 3, 27–41.
  * Kundu and Cohen 2008 P. K. Kundu and I. M. Cohen, _Fluid Mechanics, 4th ed._ , Elsevier, Amsterdam, 2008.
  * Ohm and Shelley 2022 L. Ohm and M. J. Shelley, _J. Fluid Mech._ , 2022, 942, A53–1.
  * Stakgold and Holst 2011 I. Stakgold and M. J. Holst, _Green’s functions and boundary value problems, 3rd ed._ , Wiley, Hoboken, N. J., 2011.
  * Landau and Lifshitz 1987 L. D. Landau and E. M. Lifshitz, _Fluid mechanics_ , Butterworth-Heinemann, Ltd., Oxford, 2nd edn, 1987.
  * Green _et al._ 2017 R. Green, J. Toner and V. Vitelli, _Physical Review Fluids_ , 2017, 2, 104201.
  * Joshi _et al._ 2023 C. Joshi, Z. Zarei, M. M. Norton, S. Fraden, A. Baskaran and M. F. Hagan, _arXiv preprint arXiv:2304.04895_ , 2023.