# Exploring regular and turbulent flow states in active nematic channel flow
via Exact Coherent Structures and their invariant manifolds

Caleb G. Wagner  Mechanical and Materials Engineering, University of Nebraska
- Lincoln, Lincoln, NE 68588  Rumayel H. Pallock  Mechanical and Materials
Engineering, University of Nebraska - Lincoln, Lincoln, NE 68588  Michael M.
Norton  School of Physics and Astronomy, Rochester Institute of Technology,
Rochester, NY 14623  Jae Sung Park  Mechanical and Materials Engineering,
University of Nebraska - Lincoln, Lincoln, NE 68588  Piyush Grover
<EMAIL_ADDRESS>Mechanical and Materials Engineering, University of
Nebraska - Lincoln, Lincoln, NE 68588

###### Abstract

This work is a unified study of stable and unstable steady states of 2D active
nematic channel flow using the framework of Exact Coherent Structures (ECS).
ECS are stationary, periodic, quasiperiodic, or traveling wave solutions of
the governing equations that, together with their invariant manifolds,
organize the dynamics of nonlinear continuum systems. We extend our earlier
work on ECS in the preturbulent regime by performing a comprehensive study of
stable and unstable ECS for a wide range of activity values spanning the
preturbulent and turbulent regimes. In the weakly turbulent regime, we compute
more than 200 unstable ECS that co-exist at a single set of parameters, and
uncover the role of symmetries in organizing the phase space geometry. We
provide conclusive numerical evidence that in the preturbulent regime, generic
trajectories shadow a series of unstable ECS before settling onto an
attractor. Finally, our studies hint at shadowing of quasiperiodic type ECS in
the turbulent regime.

##  I Introduction

Active fluids are viscous suspensions of constituents that consume chemical
energy and convert it into mechanical work by generating stresses on the
microscale, resulting in induced flows [1]. They are governed by fully
nonequilibrium dynamics [2], and their emergent spatiotemporal structures span
multiple scales. In recent years, model active fluids have forged an
unprecedented link between physics and biology [3]. Given their unique nature,
there are several outstanding challenges in characterizing and controlling
active fluids. Among various types of active fluids, _active nematics_ [4, 5]
– suspensions of active, rod-like, and apolar components – are of particular
contemporary interest. This interest is in part motivated by the potential of
exploiting the rich phenomenology of active nematics, including spontaneous
coherent flows, dynamical vortex patterns, and chaotic hydrodynamics (i.e.,
low Reynolds number ‘active/mesoscale turbulence’) [6, 7], for design of smart
materials.

Extensile active nematics are inherently unstable to bend fluctuations [8]. In
the absence of boundaries or other influences, active nematic systems exhibit
chaotic flows driven by topological defects [9, 10, 11, 12]. While the details
vary from model to model, the general picture is that small distortions to a
uniformly aligned active nematic system grow through hydrodynamic feedback,
fueled by active stresses [13]. When these distortions saturate, they create
pairs of $\pm\frac{1}{2}$ defects, named according to their winding number. At
steady state, these defects nucleate, self-propel, and annihilate, all while
creating vortices with a correlation length $L_{\alpha}$. As the activity is
increased, defect density increases and the size of vortices concomitantly
decreases as $L_{\alpha}\sim\sqrt{K/{\alpha}}$, where $K$ is the elastic
constant and $\alpha$ is the strength of the active stress such that
$\sigma_{\alpha}\propto\alpha\mathbf{Q}$ [4], with $\mathbf{Q}$ being the
nematic alignment tensor.

Experimental efforts to modify the chaotic flows of active nematic system have
utilized ‘soft’ anisotropic substrate friction [14, 15], light-activated
activity patterning [16, 17, 18] or confinement [19, 20]. Given the importance
of confined geometries in living active systems, the role of geometric
confinement on dynamics has been explored computationally under a variety of
geometries. A key control parameter in moderating the transition from laminar
to chaotic flows is the domain size relative to $L_{\alpha}$. The ratio of
this intrinsic length scale to imposed geometric length scales, such as
channel width or domain radius, largely determines the spatiotemporal dynamics
of the director and flow fields, and has been extensively explored
experimentally and numerically [21, 22, 23, 24, 4, 20].

In two-dimensions, the general trend [25, 26] is that confining geometries
much smaller than $L_{\alpha}$ result in stationary systems while systems much
larger than $L_{\alpha}$ result in chaotic flows or mesoscale turbulence.
Between these extremes, systems exhibit a variety of geometry-dependent
dynamical states, such as coherent flows accompanied by either stationary or
nontrivial periodic dynamics of the director field. For example, in the
channel geometry, the so-called ‘dancing’ state exists in this regime and is
characterized by braiding orbits of $+\frac{1}{2}$ defects and stationary
$-\frac{1}{2}$ defects that decorate the walls of the channel [22, 24, 27].

From an engineering perspective, there is substantial interest in
understanding how to navigate the large space of spatiotemporal structures in
confined active nematics [28, 29, 30, 31, 32], by spatial, or spatiotemporal
modulation of activity. There are also fundamental unanswered questions
related to active turbulence in confinement [26, 33]: how confined active
fluids become turbulent, how to characterize them, and how to promote or
inhibit turbulence. Some recent theoretical and computational studies have
derived coarse-grained statistical descriptions [6, 34, 35, 36, 37] of active
turbulence, focusing on energy spectra and spatial organization of
statistically steady states, rather than deterministic dynamics. These studies
have revealed that in contrast with inertial turbulence, energy transfer
across scales is absent in active nematic turbulence. Furthermore, there exist
universal scaling regimes in active nematic turbulence. While providing
valuable information about the energetics and scaling properties of fully
developed active turbulence, such studies do not provide much insight into
developing capabilities for prediction or control of active turbulence,
especially in confinement.

An alternative approach to the problem of turbulence is based on dynamical
systems theory. The core premise, as originally applied to inertial
turbulence, goes back to [38, 39]. In this approach, one considers the fluid
to be a deterministic dynamical system evolving in an infinite dimensional
phase space [40]. The governing equations of a given hydrodynamic system can
be written as an abstract first-order ODE

$\dot{X}\mathop{=}F(X),$ (1)

where $X$ denotes the $N$-dimensional state of the system. In computations,
$N$ is a large but finite number corresponding to the degrees of freedom of
the discretized system. The associated flow map is
$f^{t}(X_{0})\mathop{=}X_{0}+\int_{0}^{t}F(X(\tau))d\tau$, where $X_{0}$ is
the initial condition.

The dominant flow patterns of this system are understood in terms of _Exact
Coherent Structures_ (ECS) [41, 42] and the dynamical pathways connecting
them. An ECS is a (stable or unstable) equilibrium, (time-) periodic, relative
(time-) periodic, quasi-periodic or traveling wave solution of Eq. 1. Each
type of ECS can be defined as a solution to a corresponding fixed point
equation (FPE). In phase space, equilbria are 0D manifolds (points), periodic
orbits and travelling waves are 1D manifolds (curves), and relative periodic
orbits are 2D manifolds (surfaces). The FPE for an equilibrium solution
$X_{\text{eq}}$ is $F(X_{\text{eq}})\mathop{=}0$, while any point
$X_{\text{P}}$ on a periodic orbit (PO) satisfies
$f^{T}(X_{\text{P}})\mathop{=}X_{\text{P}}$, where $T$ is the time period.
Similarly, a point $X_{\text{RP}}$ on a relative periodic orbit (RPO)
satisfies $f^{T}(X_{\text{RP}})\mathop{=}\tau_{x}(\ell)X_{\text{RP}}$, where
$\tau_{x}(\ell)$ is an operator that shifts the solution in the streamwise
direction by a distance $\ell$. While stable ECS can be computed by time-
dependent simulations of Eq. 1, the computation of unstable ECS requires the
solution of the above mentioned FPEs.

The stability of each ECS can be determined by linearization about its
trajectory in the phase space. In general, if a $d$-dimensional ECS in an
$N-$dimensional phase space has $n_{u}$ linearly unstable directions forming
the local unstable subspace, and $n_{s}$ linearly stable directions forming
the local stable subspace, then $n_{u}+n_{s}+d=N$. Hence, away from
bifurcation points, each unstable ECS behaves like a (high-dimensional)
saddle. Associated to such an ECS are $(n_{u}+d)-$dimensional unstable
manifolds and $(n_{s}+d)-$dimensional stable manifolds. Mathematically, the
stable and unstable manifolds are sets of initial conditions that converge to
the ECS in forward and backward time, respectively. These invariant manifolds
are dynamical pathways connecting distant regions of phase space.

The steady states that are directly observed in experiments or time-dependent
numerical simulations correspond to various types of _stable_ ECSs. It has
been conjectured that turbulent flows correspond to chaotic trajectories
meandering through the phase space and visiting the neighborhoods of different
ECS (all of which are _unstable_) in a recurring fashion [43, 44, 45].
Therefore, the ECS and their invariant manifolds act as an organizing template
for the complicated spatiotemporal motion of the fluid. This minimal
description of complicated flows is well-suited for devising near-optimal
control inputs using activity patterning. Such an approach has been
successfully applied in several low-dimensional nonlinear dynamical systems
[46, 47, 48], where (near-) optimal control trajectories have been computed
solely based on the knowledge of various equilibria or periodic orbits, and
their stable and unstable manifolds, without solving any formal optimization
problem.

In our earlier work [49], we reported the co-existence of a chaotic attractor,
multiple regular attractors (i.e., stable ECS), and a large number of unstable
ECS in the phase space, at a single value of activity below the turbulent
transition. This paper extends this previous analysis in several directions.
In section III.1, we present a phase diagram constructed via a comprehensive
study of attractors for a wide range of activity values. We find that above a
critical activity value, there are no regular (i.e., non-chaotic) attractors
in the system, and all trajectories are chaotic, signifying the transition to
turbulence. Moving beyond attactors, we discuss the computation and
classification of various stable and unstable ECSs in preturbulent and
turbulent activity regimes in section III.2. We show that the explicit use of
system symmetries enables the computation and classification of system states.
The symmetry viewpoint [50, 51] also aids in visualization and understanding
of the phase space geometry.

The conjecture [40] that turbulent trajectories shadow various ECSs, guided by
the associated invariant manifolds and heteroclinic connections, has been
confirmed in recent numerical [52] and experimental [53, 45] studies in
certain inertial turbulent systems. In section III.3, we review the tools
introduced in these prior studies [52, 53], and employ them to provide
numerical evidence of shadowing of ECSs and their invariant manifolds in the
active nematic system in the preturbulent regime. Finally, we also discuss our
efforts to confirm shadowing of ECS in the turbulent regime in the same
section.

##  II Model, Symmetries and Computational Framework

Following earlier work [54, 55, 56, 22], we model the active nematic system in
terms of the velocity $\mathbf{u}(\mathbf{r},t)\mathop{=}(u,v)$, and nematic
alignment tensor $\mathbf{Q}(\mathbf{r},t)$. The symmetric and traceless
$\mathbf{Q}$ tensor can be parameterized as
$Q_{ij}=q(n_{i}n_{j}-0.5\delta_{ij})$, where the scalar $q$ and unit vector
$\hat{\mathbf{n}}$ describe the degree and direction of nematic ordering. The
governing equations are:

$\displaystyle\begin{split}&\rho\left(\partial_{t}+\mathbf{u}\cdot\boldsymbol{\nabla}\right)\mathbf{u}=-\boldsymbol{\nabla}p+2\eta(\boldsymbol{\nabla}\cdot\mathbf{E})-\alpha(\boldsymbol{\nabla}\cdot\mathbf{Q}),\\\
&\left(\partial_{t}+\mathbf{u}\cdot\boldsymbol{\nabla}\right)\mathbf{Q}+\mathbf{W}\cdot\mathbf{Q}-\mathbf{Q}\cdot\mathbf{W}=\lambda\mathbf{E}+\gamma^{-1}\,\mathbf{H},\\\
&\qquad\boldsymbol{\nabla}\cdot\mathbf{u}=0,\end{split}$ (2)

where

$\displaystyle E_{ij}$
$\displaystyle=\frac{1}{2}\left(\partial_{i}u_{j}+\partial_{j}u_{i}\right);\qquad
W_{ij}=\frac{1}{2}\left(\partial_{i}u_{j}-\partial_{j}u_{i}\right)$ (3)
$\displaystyle\mathbf{H}$
$\displaystyle=A\left[\mathbf{Q}\mathbin{-}b\mathrm{Tr}\left(\mathbf{Q}^{2}\right)\mathbf{Q}\right]\mathbin{+}K\nabla^{2}\mathbf{Q}.$
(4)

The response of the material to activity is captured in the Navier-Stokes
equation as the active stress $-\alpha\mathbf{Q}$, which models an active
particle as a force dipole [57]. We do not include passive elastic stresses in
this equation, since they are subdominant [49, 23, 58]. In the
$\mathbf{Q}$-equation, the l.h.s. is the convective derivative of $\mathbf{Q}$
with respect to $\mathbf{u}$. In addition to the usual
$(\mathbf{u}\cdot\boldsymbol{\nabla})$ contribution, there is also the
commutator product $\mathbf{W}\cdot\mathbf{Q}-\mathbf{Q}\cdot\mathbf{W}$,
which accounts for rotations of $\mathbf{Q}$ in the convected coordinates. On
the r.h.s, $\lambda\mathbf{E}$ is the flow-alignment term, and $\mathbf{H}$
derives from a free energy functional that penalizes distortions and drives
the creation of nematic order. In the present study, we neglect flow-alignment
to focus on essential aspects of the system [12, 59]. Other aspects of the
nematohydrodynamic model, including generalizations and detailed physical
interpretation, can be found in [57]. The domain is a periodic 2D channel,
parameterized as $(x,y)\in\left[0,L\right]\times\left[0,h\right]$, with $x$
the periodic coordinate. The channel walls impose a no-slip boundary condition
on $\mathbf{u}$ and strong perpendicular anchoring on $\mathbf{Q}$.

For our simulations, we adopt the non-dimensionalization used by Ref. [58],
where the fundamental length and time scales $\ell,\tau$ are defined as

$\ell=\frac{K^{1/2}}{A^{1/2}},\qquad\tau=\frac{\gamma}{A}$ (5)

The corresponding velocity scale $v_{0}$ is

$v_{0}=\frac{\ell}{\tau}=\frac{K^{1/2}}{A^{1/2}}\cdot\frac{A}{\gamma}=\frac{K^{1/2}A^{1/2}}{\gamma}$
(6)

In these units, equations (2) become

$\displaystyle\mathrm{Re}_{\mathrm{n}}\left(\partial_{t}+\mathbf{u}\cdot\boldsymbol{\nabla}\right)\mathbf{u}=-\boldsymbol{\nabla}p+2(\boldsymbol{\nabla}\cdot\mathbf{E})-\frac{\mathrm{R}_{\mathrm{a}}}{\mathrm{Er}}\left(\boldsymbol{\nabla}\cdot\mathbf{Q}\right),$
$\displaystyle\left(\partial_{t}+\mathbf{u}\cdot\boldsymbol{\nabla}\right)\mathbf{Q}+\mathbf{W}\cdot\mathbf{Q}-\mathbf{Q}\cdot\mathbf{W}=\lambda\mathbf{E}+\mathbf{H},$
(7) $\displaystyle\qquad\boldsymbol{\nabla}\cdot\mathbf{u}=0.$
$\displaystyle\mathbf{H}=\mathbf{Q}\mathbin{-}b\mathrm{Tr}\left(\mathbf{Q}^{2}\right)\mathbf{Q}\mathbin{+}\nabla^{2}\mathbf{Q},$
(8)

where

$\displaystyle\mathrm{Re}_{\mathrm{n}}$ $\displaystyle=\frac{\rho
v_{0}\ell}{\eta},$ (9) $\displaystyle\mathrm{Er}$ $\displaystyle=\frac{\eta
v_{0}\ell}{K}=\frac{\eta}{K}\cdot\frac{K}{\gamma}=\frac{\eta}{\gamma},$ (10)
$\displaystyle\mathrm{R}_{\mathrm{a}}$
$\displaystyle=\frac{K}{A}\cdot\frac{\alpha}{K}=\frac{\alpha}{A}.$ (11)

$\mathrm{R}_{\mathrm{a}}$ is the dimensionless activity, $\mathrm{Er}$ is the
Ericksen number, and $\mathrm{Re}_{\mathrm{n}}$ was identified in Ref. [58] as
the microscopic Reynolds number. Consistent with prior work [22, 60], we fix
the parameters $\mathrm{Er=1}\text{ and }\mathrm{Re_{n}}=0.0136$. The channel
dimensions are $L\mathop{=}80$ and $h\mathop{=}20$ in the dimensionless units.

###  II.1 Symmetries

A dynamical system of the form given by Eq. (1) is said be equivariant under a
group action $\sigma$ if the time evolution and the group operation commute,
i.e. for every state $X$ and time $T$, the relationship
$\sigma(f^{T}(X))=f^{T}(\sigma(X))$ is satisfied. The system of Eqs. (2) in a
periodic channel is equivariant under the continuous one-parameter group of
$x$ translations, $\tau_{x}(\ell)$, as well as the $x$ and $y$ reflections,
denoted $\sigma_{x}$ and $\sigma_{y}$, where

$\displaystyle\sigma_{x}[u,v,Q_{11},Q_{12}](x,y)=[u,-v,Q_{11},-Q_{12}](x,h-y),$
$\displaystyle\sigma_{y}[u,v,Q_{11},Q_{12}](x,y)=[-u,v,Q_{11},-Q_{12}](L-x,y),$
$\displaystyle\tau_{x}(\ell)[u,v,Q_{11},Q_{12}](x,y)=[u,v,Q_{11},Q_{12}](x+\ell,y),$
(12)

and $\ell$ is any real number. We will use the notation $T_{k}$ to denote the
discrete $k$-fold translation group with operation $\tau_{x}(\frac{L}{k})$. We
use the following convention while naming various ECS: ‘T$k/m$’ refers to an
ECS with $T_{k}$ symmetry, and the integer $m$ is an identifier used to
differentiate between different ECS with the same symmetry. Suppose an initial
condition $X_{0}$ is symmetric with respect to any of the above group
operations or their compositions, e.g., let $X_{0}=\sigma_{x}(X_{0})$. Then,
the future states $f^{T}(X_{0})$ of the system will retain the same symmetries
[42], and hence, $f^{T}(X_{0})=\sigma_{x}(f^{T}(X_{0}))$ for all $T>0$. In
that case, the state $X_{0}$ and all its future iterates are said to belong to
the $\tau_{x}$ symmetry subspace $Fix(\sigma_{x})\subset\mathbb{R}^{N}$, where
$Fix(\sigma_{x})=\\{X\in\mathbb{R}^{n}|\sigma_{x}(X)=X\\}$. Our discussions in
the following sections will reveal that these symmetries are a powerful tool
for finding and categorizing various ECSs, and analyzing the phase space
geometry.

###  II.2 Computational Framework

We have developed an open-source computational toolbox titled ‘Exact Coherent
Strucures in Active Matter’ (ECSAct)[61] using the open-source pseudospectral
code Dedalus [62]. Dedalus is a Python-based, MPI-parallelized general purpose
solver for initial-value, boundary-value, and eigenvalue problems on
spectrally representable domains. Nearly arbitrary equations can be
efficiently solved via state-of-the-art algorithms that exploit sparsity,
arithmetic trees, and transforms. Pseudospectral methods are well-suited for
solving nonlinear PDEs in simple geometries [63]. For channel geometries,
Dedalus implements a Fourier basis for the periodic directions and Chebyshev
polynomials for the wall-normal direction.

While Dedalus can natively solve time-dependent problems, it currently lacks
sophisticated solvers that can handle highly nonlinear multi-dimensional FPEs
such as the ones that need to be solved to compute unstable ECSs. To solve
such FPEs, ECSAct uses modified Newton-Raphson algorithms [64]. Two key
ingredients are adaptive ‘hookstep’ step-size selection to improve global
convergence [65], and a matrix-free GMRES [66, 67] algorithm for solving the
linear BVP at each iteration. The matrix-free methods are essential because
they scale efficiently to the large problem dimensions encountered in
hydrodynamic simulations. ECSAct also contains routines for linear stability
analysis of the various types of ECSs, as well as for performing symmetry
operations such as projections into a symmetric subspace, etc. More details on
these algorithms and implementation of ECSAct can be found in the ECSAct
documentation available on the Github repository [61].

All ECS and heteroclinic connections reported here were computed using $128$
Fourier modes and $32$ Chebyshev modes, corresponding to a phase space
dimension of $\sim 4\mathbin{\times}128\mathbin{\times}32\mathop{=}16384$.

##  III Results

###  III.1 Attractors

####  III.1.1 Search strategy

Nonlinear dynamical systems can have multiple attracting states. To discover
an attractor by forward time integration, the initial conditions must be
within the corresponding basin of attraction. However, for multiple
attractors, there is no general method for ensuring that all basins of
attraction have been sampled. Thus, heuristic methods are usually required. In
this section, we explain the heuristics we have used for mapping out
attractors in AN channel flow.

As a starting point, we refer to our previous work [49] on AN channel flow at
a single value of activity. There, we found three coexisting attractors:
three- and four-fold vortex lattice (‘dancing disclinations’) PO states, and a
localized chaotic attractor. Because all three were well-separated from each
other in phase space, a cursory search by forward time integration could
easily have overlooked one by failing to provide initial conditions in each
basin of attraction. Finding all three required a more principled approach, in
which initial conditions were chosen in a set of discrete symmetry subspaces.
The idea is that the dynamics tends to prefer configurations that match the
system’s length scales, corresponding to flow states with exact or approximate
discrete translational symmetry. In a channel with aspect ratio width/height
$\approx 4$, it is plausible that the non-chaotic attractors have 3, 4 or 5
fold translational symmetry: for example, a row of 3, 4 or 5 identical (or
nearly so) vortices. Indeed, this configuration coincides with the stable
vortex lattice PO consistently observed in previous studies. This physical
reasoning concurs with the more general dynamical systems context, where the
role of symmetry subspaces in organizing the phase space geometry is well-
established.

Figure 1: The three equilibria with zero or 1D flow for $\mathrm{R_{a}}=3.4$.

Our strategy, therefore, is to search for attractors on or near symmetry
subspaces that match the system’s preferred length scales. One method is to
choose initial conditions in a symmetry subspace and explicitly enforce that
symmetry under forward time integration. However, doing so would overlook key
features that are only represented in the full phase space. These include
dynamical connections between ECS in different symmetry subspaces, as well as
any attractors that are small perturbations from a symmetric state, e.g., a
vortex lattice with one of the vortices slightly distorted. Therefore, we
perform our searches for attractors in the full phase space.

$n_{x}$ |  $\text{Re}(\lambda)$ |  |  $n_{x}$ |  $\text{Re}(\lambda)$ |  |  $n_{x}$ |  $\text{Re}(\lambda)$
---|---|---|---|---|---|---|---
6 |  0.623 |  |  10 |  0.372 |  |  9 |  0.207
7 |  0.603 |  |  3 |  0.296 |  |  10 |  0.137
5 |  0.595 |  |  11 |  0.254 |  |  12 |  0.12
8 |  0.55 |  |  7 |  0.249 |  |  5 |  0.115
4 |  0.497 |  |  8 |  0.246 |  |  11 |  0.042
9 |  0.471 |  |  6 |  0.208 |  |  2 |  -0.021
Table 1: Features of the leading eigenvectors of the dynamics linearized about
UNI at $\mathrm{R_{a}}=3.4$. The eigenvectors are pure Fourier modes in $x$,
of the form $e^{\pm 2\pi n_{x}/L}g(y)$, where $L=80$ is the channel width.
$\lambda$ is the eigenvalue, so positive (negative) $\text{Re}(\lambda)$ means
the eigenvector represents an unstable (stable) direction.

In practice, the set of equilibria with 1D flow, henceforth referred to as ‘1D
equilibria’, shown in Fig. 1 serve as convenient ‘entry points’ to the
symmetry subspaces. These three equilibria [68, 69, 49] have continuous
translational symmetry, so they belong to every discrete translational
symmetry subspace. For the same reason, the unstable eigenvectors of the
system linearized about these equilibria are of the form $e^{\pm 2\pi
n_{x}/L}g(y)$ and therefore possess discrete translational symmetry. These
symmetries also tend to be the most dynamically relevant ones, as they are
precisely those that grow and persist as a trajectory leaves the equilibrium.
By contrast, the dynamically unfavorable symmetries tend to be represented by
the stable eigenvectors, which are damped out in the neighborhood of the
equilibrium. Table 1 illustrates this idea for UNI at $\mathrm{R_{a}}=3.4$. In
this case, the unstable eigenvectors have wavelengths ranging from one-third
to one-twelfth of the channel width $L$, corresponding to discrete symmetries
$T_{3},T_{4},\ldots,T_{12}$. As demonstrated in the following sections, these
symmetries are indeed important for characterizing typical time-dependent
trajectories, while higher order symmetries, $T_{n>12}$, are inconsequential.

In the end, our search consisted of the following steps:

  1. 1.

For each $\mathrm{R_{a}}\in\\{0.1\,n,n=1\ldots 50\\}$, compute the unstable
eigenvectors of UNI, LAN, and ZF.

  2. 2.

For each unstable eigenvector $\mathbf{v}$ construct the initial condition
$\mathbf{X}_{0}+\epsilon(\mathbf{v}/||\mathbf{v}||)$, where $\mathbf{X}_{0}$
is one of UNI, LAN, or ZF; and $\epsilon=10^{-5}$. From this initial
condition, run a time-dependent simulation for $10^{4}\tau$ time units. Denote
the final state $\mathbf{X}_{\infty}$.

  3. 3.

Run the Newton–Raphson solver with $\mathbf{X}_{\infty}$ as the initial guess.
If the solver converges, compute the stability and see if it is an attractor.

We also performed branch continuation for several ECS, checking whether an
unstable solution branch became stable at higher or lower activity. However,
most attractors were found using the time-dependent method described above.

####  III.1.2 Description and Classification

Figure 2: Snapshots of some of the attractors: The vortical flows (top row)
and chaotic flow (bottom left) are at $\mathrm{R_{a}}=3.0$ and are shown with
the same colorbar scales. The modulated wave (bottom right) is at
$\mathrm{R_{a}}=1.1$; due to the smaller activity, it is shown with different
colorbar scales. The alternating vortex lattice RPO and the modulated wave RPO
both have shift-reflect symmetry $\sigma_{x}T_{2}$, the rolling vortex lattice
RPO has three-fold translation symmetry $T_{3}$, while the chaotic flow does
not have any symmetry.

The non-chaotic attractors can be categorized by the symmetries of Eq. 2. For
$0<\mathrm{R_{a}}<0.7$, one or more of the 1D equilibria (ZF, UNI, LAN) is
stable. The equilibria or POs borne as a result of bifurcations from these
attractors break the continuous translational symmetry, while picking a
discrete translational symmetry. The _stable_ ECS (i.e., the attractors) will
be the ones that pick up symmetries with dynamically favorable length scales,
which in turn depend on the system’s intrinsic length scales. For
$\mathrm{R_{a}}\approx 0.7$, where UNI becomes unstable, an RPO with $T_{4}$
symmetry is born through a homoclinic bifurcation. This RPO cycles between the
neighborhoods of the UNI equilibrium and the PO corresponding to a 4-fold
vortex lattice state. The next attractor to appear (near
$\mathrm{R_{a}}\approx 1.1$) is an RPO invariant with respect to
$\sigma_{x}T_{2}$, which describes a _shift-reflect_ symmetry [70]. Though
this particular attractor (labeled ‘modulated wave’ (MW)) occupies a narrow
range of activity (5th row in Fig. 3), the $\sigma_{x}T_{2}$ symmetry appears
repeatedly in other attractors all the way to $\mathrm{R_{a}}\approx 4.1$. A
snapshot of a MW RPO is shown in Fig. 2 (bottom right).

For the channel dimensions and material constants considered here, the most
robust attractors have a $T_{3}$ symmetry. Although there is a stable PO with
$T_{4}$ symmetry (green bar in Fig. 3), it is stable over a relatively small
range of activity, and even then only weakly so.

The symmetries loosely correlate with additional flow characteristics. In
particular, the $T_{n}$ symmetries are associated with one of the following:
(a) cycling between a nearly unidirectional flow and a vortex lattice PO (4th
row in Fig. 3); (b) left- or right-drifting vortex lattice RPOs (labeled
‘rolling’ vortex lattice, row 7); or (c) vortex lattice PO (row 8). In cases
(b) and (c), the vortices are co-rotating. A snapshot of a rolling vortex
lattice RPO is shown in Fig. 2 (top right). In contrast, because $\sigma_{x}$
switches the sign on the vorticity, the shift-reflect symmetry
($\sigma_{x}T_{2}$) tends to produce a pair of counter-rotating vortices (row
6). A snapshot of this RPO, labeled ‘alternating’ vortex lattice, is shown in
Fig. 2 (top left).

We did not find any non-chaotic attractors in the activity range
$1.1\leq\mathrm{R_{a}}\leq 1.4$ for the zero flow alignment case considered in
this work. By repeating similar calculations for several non-zero values of
flow-alignment, we verified that generically, at least one non-chaotic
attractor can be found for all values of activities below the turbulent
transition.

Figure 3: Attractors as a function of activity $\mathrm{R_{a}}\equiv\alpha/A$
(Eq. 11). The channel dimensions are $80\ell\times 20\ell$, where $\ell$ is
defined from the material constants of the nematic (Eq. 5). The flow alignment
$\lambda$ is $0$. The horizontal bars group the attractors by symmetries and
other qualitative flow characteristics. A blue color indicates an equilibrium,
orange an RPO, and green a PO. Some of the bars are composed of more than one
ECS; for example, there are two qualitatively similar 3-fold rolling vortices
between $\mathrm{R_{a}}\\!=\\!3.34$ and $\mathrm{R_{a}}\\!=\\!3.86$. Real-
space flow snapshots are shown in Fig. 2, and the full dataset of ECS is given
in the Supplementary material [71].

####  III.1.3 Chaotic attractors

Surprisingly, there is evidence of chaotic attractors for activity as low as
$\mathrm{R_{a}}\approx 1.1$. By contrast, stable ECS exist all the way up to
$\mathrm{R_{a}}\approx 4.3$. In between, the phase space can be partitioned
into basins of attraction for the stable ECS and any chaotic attractors. The
former consists of all points in phase space that converge to a stable ECS as
$t\rightarrow+\infty$, and the latter consists of all points that converge to
a chaotic attractor as $t\rightarrow+\infty$. For visualization purposes, we
project trajectories on to a reduced 3D phase space $(\langle u\rangle,\langle
v^{2}\rangle,\langle Q_{11}\rangle)$ [49], where $\langle.\rangle$ denotes the
instantaneous channel average. Visual inspection suggests that chaotic
trajectories in this phase space projection are localized to a region separate
from the stable ECS, see Fig. 4 for the case of $\mathrm{R_{a}}=3.4$. The
symmetry-reduced distance (Eq. 14) between a chaotic trajectory and the
attractors at $\mathrm{R_{a}}=3.4$ was found to be $O(1)$, providing further
proof of our assertion.

Interestingly, the unstable manifold of the unidirectional equilibrium appears
to play a role in dividing these two regions. In our previous work [49], we
found robust heteroclinic connections between UNI and certain stable ECS, as
well as trajectories connecting UNI to the chaotic attractor. Overall, this
phase space partitioning is reminiscent of the subcritical transition to
turbulence in pipe flow [72, 73, 43]. These studies posit that a codimension-1
manifold separates the basins of attraction of the laminar and the chaotic
states. Such a manifold corresponds to the stable manifold of an unstable
equilibria (or travelling wave) with one unstable direction. The question of
existence and characterization of a similar edge manifold in the phase space
of active nematic channel flow is left for future work.

Figure 4: A chaotic trajectory and the three attractors projected on the
$(\langle U\rangle,\langle V^{2}\rangle)$ (left) and $(\langle
V^{2}\rangle,\langle Q_{xx}\rangle)$ (right) planes, for $\mathrm{R_{a}}=3.4$.

###  III.2 Unstable ECS

In addition to the regular (stable ECS) and chaotic attractors, the phase
space is populated by several unstable ECSs. The unstable ECSs are coherent
steady states that are not _directly_ observed in experiments or time-
dependent numerical simulations. Our interest in computing and characterizing
these states is motivated by their crucial role organizing the global dynamics
in both the non-turbulent and turbulent regime.

In several classes of dynamical systems, the natural dynamical pathways
between a pair of stable steady states are often found to be mediated by
_unstable_ steady states. Consider the recent experimental study of the
dynamics of a ball rolling on a 2D multi-well surface [74]. This conservative
two degree-of-freedom system has a 3D phase space for each fixed value of
total energy. A trajectory can transition between a pair of two wells
(potential energy minima) via a bottleneck that opens up at an energy level
above that of the intermediate saddle point of the energy surface, see Fig. 5.
In addition to this energetic requirement however, the desired transition can
happen only if the trajectory lies inside the 2D cylindrical stable manifold
of the _unstable_ periodic orbit that exists around the saddle point of the
energy surface. If this latter condition is met, the trajectory approaches the
saddle point region while travelling inside the stable manifold, and is
eventually transferred to the second well via the unstable manifold of the
same periodic orbit. If the condition is not met, the trajectory is reflected
back to first the well, even if it satisfies the energetic requirement.

Figure 5: Left: (Top) Potential energy surface of a 2 DOF conservative
dynamical system with two minima and three saddles highlighted. (Bottom) The
phase space near the energy surface saddle can be locally described as a
direct product of a saddle fixed point and a center fixed point on orthogonal
planes. Right: Periodic orbit around the saddle in yellow has 2D stable and
unstable manifolds. A trajectory transiting from the left potential well to
the right potential well travels inside the stable manifold towards the PO,
and inside the unstable manifold while moving away from the PO. Hence, the PO
here acts like an extended saddle point.

The above mentioned example highlights the role of invariant manifolds of
saddle type (unstable) steady states in determining global phase space
transport. Similar results have been obtained in orbital mechanics [46, 75],
ship dynamics [76], beam buckling [77], and solitary wave propagation [78],
etc. This body of work suggests that the unstable ECSs are expected to play a
similarly important role in energy efficient dynamical pathways that transfer
the system between various attractors in the preturbulent regime. The geometry
of phase space can also assist in the discovery of new flow states that may
become stable at different parameter values.

####  III.2.1 Search strategy

We design our search for unstable ECS keeping in mind two groups of objects in
phase space: (1) the unstable manifolds of the 1D equilibra and (2) any
chaotic attractors. By computing ECS on or near these manifolds, we hope to
gain information on nearby time-dependent trajectories. In particular, group 1
is relevant for trajectories that begin near one of the 1D equilibria; such
trajectories, at least initially, will approximately follow the equilibrium’s
unstable manifold.

The Newton-Raphson solver used for solving the FPEs requires good initial
guesses to converge to an unstable ECS. Given our goal of using the ECSs to
characterize time-dependent trajectories, it is natural to take initial
guesses from the latter. For ECSs on or near the unstable manifold of an
equilibrium, we take guesses from trajectories that begin near the equilibrium
and are perturbed along one of its unstable directions. These are actually the
same trajectories defined in section III.1, which we used to find attractors.
The only difference is that now, we are taking guesses from near the beginning
of the trajectory (roughly the first $0<t<1000\tau$) rather than the end. For
ECSs on or near chaotic attractors, we use the same set of trajectories, but
exclude those that end on an attracting ECS. Our initial guesses are then
taken from the asymptotic part of the remaining trajectories.

Depending on the trajectory used to provide initial guesses, the fraction of
solver instances that converge can vary from $0$ to $1$. Among those that
converge, not all give unique ECSs, so it is essential to detect and remove
duplicates. For our case, this is complicated by the fact that, due to
invariance with respect to space and time translations, POs and RPOs are
2-dimensional manifolds in the full phase space. In particular, two separate
solver instances generally converge to different points on the surfaces, and
we need a criterion to detect when two points are on the same surface. To
account for this situation, we use the distance measure described in section
III.3.2, which is $0$ if and only if the two ECSs are equivalent up to space
and time translations. Thus, if the distance between two ECSs is $0$, they
reside on the same ‘ECS surface’, and we remove one of them from our dataset
as a duplicate. For a few values of activity, we also ran the Newton-Raphson
solver in one or more symmetry subspaces, which helped to identity some of the
more unstable ECS. Finally, for many of the ECSs, we performed branch
continuation in $\mathrm{R_{a}}$, which sometimes uncovered ECSs that had been
missed by other searches.

However, none of these searches should be considered exhaustive. In the
presence of chaos, there will generally be an infinite number of ECS with
arbitrarily large periods. Due to the unpredictable convergence behavior of
the Newton-Raphson solver, there is also no guarantee that a given number of
searches will uncover all ECS with a given property, such as periods less than
some value. A more realistic goal is to design a search strategy that is at
least biased toward the most dynamically relevant ECSs, and our method of
searching along time-dependent trajectories indeed serves this purpose.
Nevertheless, there are other biases that may skew our ECS dataset in
different directions. In particular, the Newton-Raphson solver cannot look for
ECSs with arbitrarily large periods because it must integrate over at least
one period. Such ECSs could be important in some settings, such as near a
homoclinic bifurcation. More generally, the dynamical relevance of an ECS is
more closely related to its Floquet exponents than its Floquet multipliers
[52]. The former describe the local divergence of trajectories near a point on
an ECS, whereas the latter quantifies growth or decay of perturbations over an
entire period. By contrast, the Newton-Raphson solver is more likely to
converge to ECSs with small Floquet multipliers. In future work, we will
explore the use of optimization-based methods [79, 80] for finding ECS with
larger periods.

####  III.2.2 Description and classification

ECS data for $\mathrm{R_{a}}=\\{0.1\,n,n=1\ldots 50\\}$ can be found in the
Supplementary material [71]. For $\mathrm{R_{a}}<0.9$, we documented six
equilibria, including the three 1D equilibria in Fig. 1. This number is much
smaller than the number of ECS in our datasets at higher activity. While one
does indeed expect fewer ECS at small activity, this result may also be
influenced by the smaller number of searches performed in this regime.

Our focus instead was on $\mathrm{R_{a}}>0.9$, where RPOs and higher-
dimensional invariant objects begin to dominate the phase space geometry. For
each value of $\mathrm{R_{a}}$ above this threshold, we identified at least 10
ECS. Most of these are RPOs, though POs, TWs, and EQ are also present, as
discussed below. For the analysis in section III.3 below, more thorough
searches were conducted at a few values of $\mathrm{R_{a}}$. For example, the
$\mathrm{R_{a}}=3.4$ dataset (used in Figs. 8 and 9) contains 56 ECS, and
$\mathrm{R_{a}}=4.0$ and $\mathrm{R_{a}}=4.5$ have 66 and 210, respectively.

Figure 6: Left: An unstable defectless vortex lattice steady state
(equilibrium) with 3-fold translational symmetry $T_{3}$. Right: An unstable
travelling wave with 4-fold translational symmetry $T_{4}$.

The unstable ECS embody a much larger variety of spatiotemporal patterns
compared with the attractors. The discrete translational symmetry $T_{n}$ is
present for $n$ up to 12, and includes two new types of ECS not found among
the attractors. The first are defectless, vortex lattice equilibria found by
the Newton solver only for relatively small activity ($<1.5$). As a shorthand,
we call these states ‘static vortex equilibrium’ (SVE). These equilibria are
not to be confused with stationary vortex (SV) periodic orbits with motile
defects discussed earlier. The second are traveling waves; these occur over a
wider range of activity, but were found to be unstable in every case. The
snapshots of both these types of ECS are shown in Fig. 6. The figure
‘mosaic_R45.png’ in the Supplementary material [71] visually expresses the
phase space diversity with a mosaic of all 210 ECS found at
$\mathrm{R_{a}}=4.5$, all of which are unstable.

####  III.2.3 Distribution in phase space

Figure 7: Trajectory ‘densities’ at $\mathrm{R_{a}}=4.5$ for the 201 unstable
ECS (left) and turbulent attractor (right), projected on the $(\langle
V^{2}\rangle,\langle Q_{xx}\rangle)$ plane. Each figure is generated from
evenly spaced time series of the corresponding trajectories, which are binned
and normalized as a 2D probability density function in $\langle V^{2}\rangle$
and $\langle Q_{xx}\rangle$. Trajectories tend to spend more time in the
yellow regions versus the violet regions.

Finally, it is instructive to see the phase space distribution of various ECSs
in the turbulent regime. In Fig. 7, we compare the probability density of a
long lived turbulent trajectory with that of trajectories on all the computed
ECS. The trajectory for each ECS has duration equal to its respective period.
In this 2D phase space projection, these ECSs tend to be distributed near the
unstable manifolds of the 1D equilibria. In section III.3, we provide evidence
that this conclusion holds in the full phase space as well, i.e., it is not
just an accidental feature of the low-dimensional projection in Fig. 7.

###  III.3 Transient flows

Our interest in unstable flow patterns is driven by their potential to give a
reduced-order description of trajectories that are transient or chaotic. By
_reduced-order_ , we mean theoretical descriptions containing substantially
fewer degrees of freedom than the full, time-dependent hydrodynamic equations.
Such a description is, by necessity, usually incomplete or approximate, but
will at least capture some key features of the dynamics. By _transient
trajectories_ , we mean finite-time trajectories connecting the neighborhood
of an unstable ECS to one of the attractors. For example, one might be
interested in the behavior of a system initialized near an unstable
equilibrium and subsequently allowed to evolve toward one of the attractors in
Fig. 3. As we demonstrate below, there is substantial variety in the pathways
a system might take, so there is no single representative trajectory.
Similarly, there is no finite representative trajectory for a turbulent flow;
in both cases, the complicated nonlinear dynamics makes it difficult to
develop a reduced-order description.

As mentioned in the introduction, the dynamical systems framework [41]
postulates that such a reduced-order description can be built out of unstable
ECS. A recurring motif is a typical turbulent trajectory approaching an
unstable, saddle-like ECS along the latter’s stable manifold, before departing
along the corresponding unstable manifold towards another ECS. There is
mounting experimental [53] and numerical [52, 81, 82, 83] in support of this
idea for the case of weak inertial turbulence in classical Newtonian fluids.
Refs [52, 53] refined the definition of ‘close approach’ to an ECS, taking
into account continuous symmetries and the local geometry in the neighborhood
of the ECS. In this way, they defined ‘shadowing events’ as finite time
intervals during which a turbulent trajectory and an ECS have quantitatively
similar time evolution. By applying this formulation to inertial turbulence,
they were able to experimentally detect shadowing events in turbulent Taylor-
Couette flow [53]. These shadowing events could be used, for instance, as
nodes in a network-based model of trajectory statistics. In a complementary
work, Ref. [81] employed tools from computational topology to construct a
Markov chain where the nodes represent neighborhoods of various ECSs, and the
transition rates between a pair of nodes are governed by the heteroclinic
dynamics between those ECSs. The invariant distribution of the Markov Chain is
a coarse grained representation of the probability density function of a
generic turbulent trajectory.

Although ECS have been found numerically in non-Newtonian fluid turbulence
[84, 85], the usefulness of unstable ECS in such physical settings is still a
hypothesis, potentially requiring simplifications such as confinement to small
domains and precise experimental control. For active matter systems, these
issues have not yet been explored in detail. Here, we provide two pieces of
evidence that certain unstable ECS from section III.2 form a scaffolding of
the unstable manifold of UNI, guiding trajectories _en route_ from UNI to one
of the attractors. Specifically, a typical trajectory from UNI to one of the
attractors experiences one or more shadowing events, which could be used as
basis elements of a reduced-order model, or as access points for control
inputs. Finally, we also report on our attempts to find shadowing in the
turbulent regime.

####  III.3.1 Trajectory densities in the preturbulent regime

Figure 8: Trajectory ‘densities’ along the unstable manifold of UNI in the
preturbulent regime projected on the $(\langle V^{2}\rangle,\langle
Q_{xx}\rangle)$ plane. The dimensionless activity is $\mathrm{R_{a}}=3.4$. The
full unstable manifold is 34-dimensional and impractical to compute directly.
Following the method described in section III.3.1, we instead define an
ensemble of trajectories that maps out a low-dimensional subspace. The density
map (color gradient) is generated from this ensemble by binning evenly spaced
time series over the corresponding interval: $0<t<100$ (top) or $100<t<10^{4}$
(bottom). The apparent overlap of high density regions (yellow) with unstable
ECS suggests the occurrence of shadowing, which we have confirmed for selected
cases (see section III.3.2). Note that there are three stable ECS (or
attractors): T1/9, T3/1, and T3/2, while the rest of the ECS shown are
unstable.

A simple way to study the unstable manifold of an equilibrium like UNI is to
sample it with representative trajectories. These can be generated by
perturbing the equilibrium with a linear combination of unstable eigenvectors,
which span the tangent space of the unstable manifold in the neighborhood of
the equilibrium. Here, we focus on activity $\text{R}_{\text{a}}=3.4$, where
UNI has 34 unstable eigenvectors, and therefore $34-1=33$ coefficients in a
general (infinitesimal) perturbation. However, sampling these coefficients
uniformly will produce trajectories dominated by the most unstable mode.
Hence, to also sample the subdominant directions along the unstable manifold,
we construct perturbations from pairs of eigenvectors
$\\{\mathbf{e}_{n},\mathbf{e}_{m}\\}_{n\neq m}$ only. The initial conditions
$\mathbf{X}_{\text{init}}$ then take the form:

$\mathbf{X}_{\text{init}}=\mathbf{X}_{\text{UNI}}+\epsilon\left(\mathbf{e}_{n}+\alpha\,\mathbf{e}_{m}\right)$
(13)

where $\epsilon=10^{-5}$, $n\neq m$, and $\alpha$ is a random real number
between $-1$ and $1$.

In our case, we ran several hundred simulations with these initial conditions
and let them evolve from $t=0$ to $t=10^{4}\tau$. To synthesize the large
amount of data that was generated, we constructed time-averaged _trajectory
densities_ in the 2D phase space projection $(\langle
Q_{\text{xx}}\rangle,\langle v^{2}\rangle)$. Specifically, we computed a 2D
density map by binning the $(\langle Q_{\text{xx}}\rangle,\langle
v^{2}\rangle)$ time series over a given time interval. Fig. 8 shows the
results for the time intervals $0<t<100\tau$ and $100\tau<t<10^{4}\tau$. For
comparison, the ECS we tabulated at $\mathrm{R_{a}}=3.4$ have periods between
about $2.6\tau$ and $22\tau$. Note that the top panel of 8 has two rectangular
regions cut out from the upper and lower left corners. These regions contain
the UNI and ZF equilibria, respectively, and for early times $t<100\tau$ are
characterized by especially high trajectory density. Hence, to focus on
behavior in other regions of phase space, we exclude these areas from the
overall density map.

The trajectory density plot for $0<t<100\tau$ captures the transient
statistics, and regions of high density in that plot are where trajectories
spend the most time within the initial transient period. Many of such high
density regions are visible in the top panel of Fig 8, several of which
coincide with unstable (saddle-type) and stable ECS. Nevertheless, other high-
density contours (particularly in the top panel corresponding to transients)
could not be correlated with ECSs, despite running the Newton-Raphson solver
at points along the underlying trajectories. We conjecture that these are
unstable 2-tori (quasiperiodic orbits), which are a type of ECS that our
current solver is not equipped to compute. The density plot for
$100\tau<t<10^{4}\tau$ shown in the bottom panel of Fig. 8 captures the long
term statistics. As expected, the trajectory density is higher on attracting
ECS (zoomed-in view in the lower right panel), though strong shadowing of
three unstable ECS is also apparent.

####  III.3.2 Shadowing along unstable manifolds in the preturbulent regime

The second piece of evidence, Fig. 9, focuses on two of the trajectories used
to generate Fig. 8. Our aim is to confirm that the _apparent_ shadowing of
unstable ECS (as suggested by the trajectory densities in the reduced phase
space) indeed describes a feature of the full phase space dynamics, and not
just an artifact of the low-dimensional projection in Fig. 8. To do so, it is
necessary to define a measure of distance between two trajectories in phase
space. The natural metric is simply the $L^{2}$ norm of the fields. However,
due to the continuous translational symmetry in the streamwise direction ($x$
coordinate), this metric is inadequate for identifying _dynamical_ similarity
between trajectories. For example, the $L^{2}$ norm may assign a large
distance between two trajectories that evolve nearly identically, but only up
to a constant streamwise translation. To account for this situation, we use
the following time-dependent continuous-symmetry-reduced distance between a
trajectory $X_{1}(t)$ and an ECS $X_{E}(t)$:

$d_{1E}(t)=\min_{\begin{subarray}{c}0\leq l\leq L\\\ 0\leq\tau\leq
T\end{subarray}}\|\tau_{x}(l)(X_{1}(t))-X_{E}(\tau)\|_{2},$ (14)

where $T$ is the time-period of the ECS ($T=0$ if it is an equilibria or a
travelling wave). The top left and bottom left panels in Fig. 9 show the
results of applying this measure to two trajectories initialized from Eq. 13.
The symmetry-reduced distance (color gradient) is computed at a set of fixed
timepoints (horizontal axis) and with respect to several ECSs (vertical axis).
Several shadowing events are clearly visible as horizontal red/black segments,
where the symmetry-reduced distance is small. The corresponding reduced phase
space trajectories are shown on the right. The full time series data of
distances of both trajectories from all the ECS is provided in the
Supplementary material [71].

Figure 9: Examples of shadowing along the unstable manifold of the
unidirectional equilibrium in the preturbulent regime at $\mathrm{R_{a}}=3.4$.
The color gradient plots on the left show the distances (as defined by
equation (14)) between a time-dependent trajectory and a selection of ECS at
the same value of activity. The top left and top right illustrate shadowing of
the unstable ECS T1/5 from approximately $200\tau$ to $250\tau$. The bottom
left and bottom right show two shadowing events: first, the unstable ECS T6/1
from about $30\tau$ to $70\tau$, followed by the unstable ECS
$\sigma_{x}\sigma_{y}(\text{T3/3})$ from $250\tau$ to $600\tau$. Both
simulations eventually approach T1/9, which is an attractor.

####  III.3.3 Testing for shadowing in the turbulent regime

To test whether the turbulent state shadows any ECS, we consider typical
trajectories in the turbulent regime at $\mathrm{R_{a}}=4.5$. As discussed in
Section III.2.2, we computed $210$ distinct unstable ECS at this parameter
value. Each ECS $X_{E}(t)$ has three discrete-symmetry related counterparts,
i.e., $\sigma_{x}(X_{E}(t)),\sigma_{y}(X_{E}(t))$ and
$\sigma_{y}\sigma_{x}(X_{E}(t))$. We computed continuous-symmetry-reduced
distances using Eq. 14 between each of the the 840 ECS, and the 22
trajectories starting on the unstable manifold of UNI. These trajectories are
obtained by picking initial conditions that are small perturbations of UNI
along each of its 22 unstable directions.

Unlike the preturbulent regime, we found no instances of shadowing of any of
the ECS by the turbulent trajectories. Although there are some instances of
close approaches to certain ECS, the evolution of the chaotic trajectory does
not resemble any of the ECS even for times much smaller than the typical time
period of the corresponding ECS [52].

Figure 10: Left: A subset of the distance plot chosen to highlight one of the
closest encounters between a chaotic trajectory and an ECS in the turbulent
regime at $\mathrm{R_{a}}=4.5$. The trajectory visits the neighborhood of an
ECS (labeled T4/1) at $t/\tau\approx 1500$. Right: The trajectory from
$t/\tau=\text{0 to }t/\tau\approx 1800$ in black along with the ECS T4/1.
There is clearly no shadowing of the ECS, and the trajectory appears to be
quasiperiodic. The full time series data of distances of this trajectory from
all the ECS is provided in the Supplementary material [71].

By visual inspection of the 3D phase space plot, we found that some segments
of the chaotic trajectories appear to be quasiperiodic, see the right panel of
Fig. 10 for one example. This leads us to conjecture that these trajectories
might be shadowing quasiperiodic ECS of the system. Recall that the trajectory
density plot shown Fig. 7 of the preturbulent regime was also suggestive of
shadowing of quasiperiodic ECS.

Most prior work on ECS in turbulent fluids has been limited to computing
equilibria, POs and RPOs [86]. Two recent works [87, 88] have highlighted the
importance of quasiperiodic (QPO) ECS in determining the statistics of chaotic
dynamical systems. A quasiperiodic trajectory with n incommensurate
frequencies densely covers a n-tori in the phase space. The authors argue that
since n-tori are structurally unstable for $n>2$, equilibria, POs, RPOs and
2-tori form the complete set of invariant objects for understanding the phase
space transport in chaotic systems. These two works have proposed efficient
algorithms for computing 2-tori, and also demostrated that approximate
statistics of the chaotic forced Kuramoto-Sivashinsky PDE can be captured
using a modified periodic orbit theory adapted to QPOs.

Our set of ECS could also be missing some POs/RPOs that lie in the region
where the local curvature of the chaotic attractor is too high for the Newton-
Rhapson hookstep algorithm to converge. The computation of QPO ECS in the AN
channel flow, as well as further work on shadowing of ECS in turbulent active
nematics will be taken up in a future publication.

##  IV Discussion and Concluding Remarks

We have used the ECS framework to provide a comprehensive study of the stable
and unstable steady states of the active nematic system in a 2D periodic
channel. We found that in the preturbulent regime, multiple regular attractors
co-exist with a chaotic attractor. This implies that a subcritical transition
to chaotic flow is possible in such systems. This also hints at the presence
of an edge manifold in the phase space, separating the regular states from the
chaotic ones.

Our attractor computations confirm that as activity is increased, the general
sequence of stable states i.e., zero flow, 1D flow, vortex lattice, followed
by transition to turbulence, identified in previous studies [22, 89, 26], is
correct. However, visual inspection is not enough to identify dynamic
similarity of various flow states, and phase space analysis aided by symmetry
considerations provide additional insight. In addition to the above mentioned
attractors, our computations revealed the existence of counter-rotating
lattice vortex lattice, traveling wave and modulated wave attractors.

We demonstrated that unstable ECS organize the short term dynamics of
perturbations around the unstable unidirectional steady state in the
preturbulent regime. Typical trajectories lying on the unstable manifold of
the unidirectional equilibrium were found to shadow a subset of the unstable
ECS in the transient phase, and eventually converge onto one of the attractors
in steady state.

In the turbulent regime, which is characterized by the lack of any attractors,
we found over 200 unstable ECS at a single set of parameters. Low dimensional
projections suggest that these ECS exist in the same region of phase space as
the typical turbulent trajectories of the system. However, unlike the
preturbulent regime, there was no evidence of shadowing of the PO/RPO/TW type
ECS by the turbulent trajectories, and we conjecture that these trajectories
are instead shadowing quasiperiodic ECS.

Future work will focus on computing quasiperiodic ECS as well as PO/RPOs of
larger periods using extensions of our existing Newton-Rhapson solver, or
optimization-based methods. It would also be interesting to compare the ECS
discovered in the current work to dynamic modes obtained in active turbulence
using data-driven dynamical systems approaches such as Proper Orthogonal
Decomposition and Dynamic Mode Decomposition [33]. Finally, detection of
stable and unstable ECS in an experimentally realizable active nematic system
such as the 2D annulus [24, 90, 27] will be considered.

###### Acknowledgements.

This material is based upon work supported by the U.S. Department of Energy,
Office of Science, Office of Basic Energy Sciences under award number DE-
SC0022280. This work was completed utilizing the Holland Computing Center of
the University of Nebraska, which receives support from the Nebraska Research
Initiative. We thank Aparna Baskaran, Seth Fraden, Mike Hagan, Chaitanya
Joshi, and Zahra Zarei for helpful discussions.

## References

  * Ramaswamy [2019] S. Ramaswamy, Active fluids, Nature Reviews Physics 1, 640 (2019).
  * Kurzthaler  _et al._ [2023] C. Kurzthaler, L. Gentile, and H. A. Stone,  _Out-of-equilibrium Soft Matter: Active Fluids_ (Royal Society of Chemistry, 2023).
  * Balasubramaniam  _et al._ [2022] L. Balasubramaniam, R.-M. Mège, and B. Ladoux, Active nematics across scales from cytoskeleton organization to tissue morphogenesis, Current Opinion in Genetics & Development 73, 101897 (2022).
  * Doostmohammadi  _et al._ [2018] A. Doostmohammadi, J. Ignés-Mullol, J. M. Yeomans, and F. Sagués, Active nematics, Nat. Commun. 9, 1 (2018).
  * Mestre [2022] F. S. Mestre,  _Colloidal Active Matter: Concepts, Experimental Realizations, and Models_ (CRC Press, 2022).
  * Alert  _et al._ [2022] R. Alert, J. Casademunt, and J.-F. Joanny, Active turbulence, Annual Review of Condensed Matter Physics 13 (2022).
  * Martínez-Prat  _et al._ [2019] B. Martínez-Prat, J. Ignés-Mullol, J. Casademunt, and F. Sagués, Selection mechanism at the onset of active turbulence, Nature physics 15, 362 (2019).
  * Aditi Simha and Ramaswamy [2002] R. Aditi Simha and S. Ramaswamy, Hydrodynamic fluctuations and instabilities in ordered suspensions of self-propelled particles, Phys. Rev. Lett. 89, 058101 (2002).
  * Sanchez  _et al._ [2012] T. Sanchez, D. T. N. Chen, S. J. DeCamp, M. Heymann, and Z. Dogic, Spontaneous motion in hierarchically assembled active matter, Nature 491, 431 (2012).
  * Thampi  _et al._ [2014a] S. P. Thampi, R. Golestanian, and J. M. Yeomans, Instabilities and topological defects in active nematics, Europhysics Letters 105, 18001 (2014a).
  * DeCamp  _et al._ [2015] S. J. DeCamp, G. S. Redner, A. Baskaran, M. F. Hagan, and Z. Dogic, Orientational order of motile defects in active nematics, Nature materials 14, 1110 (2015).
  * Shankar  _et al._ [2018] S. Shankar, S. Ramaswamy, M. C. Marchetti, and M. J. Bowick, Defect unbinding in active nematics, Phys. Rev. Lett. 121, 108002 (2018).
  * Baskaran and Marchetti [2009] A. Baskaran and M. C. Marchetti, Statistical mechanics and hydrodynamics of bacterial suspensions, Proceedings of the National Academy of Sciences 106, 15567 (2009).
  * Guillamat  _et al._ [2016] P. Guillamat, J. Ignés-Mullol, and F. Sagués, Control of active liquid crystals with a magnetic field, PNAS 113, 5498 (2016).
  * Guillamat  _et al._ [2017] P. Guillamat, J. Ignés-Mullol, and F. Sagués, Taming active turbulence with patterned soft interfaces, Nat. Commun. 8, 1 (2017).
  * Ross  _et al._ [2019] T. D. Ross, H. J. Lee, Z. Qu, R. A. Banks, R. Phillips, and M. Thomson, Controlling organization and forces in active matter through optically defined boundaries, Nature 572, 224 (2019).
  * Zhang  _et al._ [2021] R. Zhang, S. A. Redford, P. V. Ruijgrok, N. Kumar, A. Mozaffari, S. Zemsky, A. R. Dinner, V. Vitelli, Z. Bryant, M. L. Gardel, _et al._ , Spatiotemporal control of liquid crystal structure and dynamics through activity patterning, Nature materials 20, 875 (2021).
  * Zarei  _et al._ [2023] Z. Zarei, J. Berezney, A. Hensley, L. Lemma, N. Senbil, Z. Dogic, and S. Fraden, Light-activated microtubule-based 2D active nematic, arXiv preprint arXiv:2303.02945 (2023).
  * 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, Reconfigurable flows and defect landscape of confined active nematics, Communications Physics 2, 1 (2019).
  * Opathalage  _et al._ [2019] A. Opathalage, M. M. Norton, M. P. Juniper, B. Langeslay, S. A. Aghvami, S. Fraden, and Z. Dogic, Self-organized dynamics and the transition to turbulence of confined active nematics, PNAS 116, 4788 (2019).
  * Giomi  _et al._ [2012] L. Giomi, L. Mahadevan, B. Chakraborty, and M. Hagan, Banding, excitability and chaos in active nematic suspensions, Nonlinearity 25, 2245 (2012).
  * Shendruk  _et al._ [2017] T. N. Shendruk, A. Doostmohammadi, K. Thijssen, and J. M. Yeomans, Dancing disclinations in confined active nematics, Soft Matter 13, 3853 (2017).
  * Norton  _et al._ [2018] M. M. Norton, A. Baskaran, A. Opathalage, B. Langeslay, S. Fraden, A. Baskaran, and M. F. Hagan, Insensitivity of active nematic liquid crystal dynamics to topological constraints, Phys. Rev. E 97, 012702 (2018).
  * Chen  _et al._ [2018] S. Chen, P. Gao, and T. Gao, Dynamics and structure of an apolar active suspension in an annulus, J. Fluid. Mech. 835, 393 (2018).
  * Samui  _et al._ [2021] A. Samui, J. M. Yeomans, and S. P. Thampi, Flow transitions and length scales of a channel-confined active nematic, Soft Matter 17, 10640 (2021).
  * Thampi [2022] S. P. Thampi, Channel confined active nematics, Current Opinion in Colloid & Interface Science , 101613 (2022).
  * Joshi  _et al._ [2023] C. Joshi, Z. Zarei, M. Norton, S. Fraden, A. Baskaran, and M. Hagan, From disks to channels: Dynamics of active nematics confined to an annulus, arXiv preprint arXiv:2304.04895 (2023).
  * Norton  _et al._ [2020] M. M. Norton, P. Grover, M. F. Hagan, and S. Fraden, Optimal control of active nematics, Phys. Rev. Lett. 125, 178005 (2020).
  * Shankar  _et al._ [2022a] S. Shankar, V. Raju, and L. Mahadevan, Optimal transport and control of active drops, Proceedings of the National Academy of Sciences 119, e2121985119 (2022a).
  * Shankar  _et al._ [2022b] S. Shankar, L. V. Scharrer, M. J. Bowick, and M. C. Marchetti, Spatiotemporal control of active topological defects, arXiv preprint arXiv:2212.00666 (2022b).
  * Mozaffari  _et al._ [2021] A. Mozaffari, R. Zhang, N. Atzin, and J. J. De Pablo, Defect spirograph: Dynamical behavior of defects in spatially patterned active nematics, Physical Review Letters 126, 227801 (2021).
  * Falk  _et al._ [2021] M. J. Falk, V. Alizadehyazdi, H. Jaeger, and A. Murugan, Learning to control active matter, Physical Review Research 3, 033291 (2021).
  * Henshaw  _et al._ [2023] R. J. Henshaw, O. G. Martin, and J. S. Guasto, Dynamic mode structure of active turbulence, Physical Review Fluids 8, 023101 (2023).
  * Alert  _et al._ [2020] R. Alert, J.-F. Joanny, and J. Casademunt, Universal scaling of active nematic turbulence, Nature Physics 16, 682 (2020).
  * Giomi [2015] L. Giomi, Geometry and topology of turbulence in active nematics, Phys. Rev. X 5, 031003 (2015).
  * Linkmann  _et al._ [2019] M. Linkmann, G. Boffetta, M. C. Marchetti, and B. Eckhardt, Phase transition to large scale coherent structures in two-dimensional active matter turbulence, Phys. Rev. Lett. 122, 214503 (2019).
  * Linkmann  _et al._ [2020] M. Linkmann, M. C. Marchetti, G. Boffetta, and B. Eckhardt, Condensate formation and multiscale dynamics in two-dimensional active suspensions, Phys. Rev. E 101, 022609 (2020).
  * Hopf [1948] E. Hopf, A mathematical example displaying features of turbulence, Commun. Pure Appl. Math. 1, 303 (1948).
  * Ruelle and Takens [1971] D. Ruelle and F. Takens, On the nature of turbulence, Les rencontres physiciens-mathématiciens de Strasbourg-RCP25 12, 1 (1971).
  * Cvitanović [2013] P. Cvitanović, Recurrent flows: the clockwork behind turbulence, J. Fluid. Mech. 726, 1 (2013).
  * Graham and Floryan [2020] M. D. Graham and D. Floryan, Exact coherent states and the nonlinear dynamics of wall-bounded turbulent flows, Annu. Rev. Fluid Mech. 53 (2020).
  * Cvitanovic  _et al._ [2005] P. Cvitanovic, R. Artuso, R. Mainieri, G. Tanner, G. Vattay, N. Whelan, and A. Wirzba, Chaos: classical and quantum, ChaosBook. org (Niels Bohr Institute, Copenhagen 2005) 69, 25 (2005).
  * Park and Graham [2015] J. S. Park and M. D. Graham, Exact coherent states and connections to turbulent dynamics in minimal channel flow, J. Fluid. Mech. 782, 430 (2015).
  * Budanur  _et al._ [2019] N. B. Budanur, A. S. Dogra, and B. Hof, Geometry of transient chaos in streamwise-localized pipe flow turbulence, Phys. Rev. Fluids 4, 102401 (2019).
  * Suri  _et al._ [2020] B. Suri, L. Kageorge, R. O. Grigoriev, and M. F. Schatz, Capturing turbulent dynamics and statistics in experiments with unstable periodic orbits, Phys. Rev. Lett. 125, 064501 (2020).
  * Koon  _et al._ [2000] W. S. Koon, M. W. Lo, J. E. Marsden, and S. D. Ross, Heteroclinic connections between periodic orbits and resonance transitions in celestial mechanics, Chaos 10, 427 (2000).
  * Sakamoto [2013] N. Sakamoto, Case studies on the application of the stable manifold approach for nonlinear optimal control design, Automatica 49, 568 (2013).
  * Flaßkamp  _et al._ [2012] K. Flaßkamp, S. Ober-Blöbaum, and M. Kobilarov, Solving optimal control problems by exploiting inherent dynamical systems structures, Journal of Nonlinear Science 22, 599 (2012).
  * Wagner  _et al._ [2022] C. G. Wagner, M. M. Norton, J. S. Park, and P. Grover, Exact coherent structures and phase space geometry of preturbulent 2D active nematic channel flow, Physical Review Letters 128, 028003 (2022).
  * Golubitsky and Stewart [2003] M. Golubitsky and I. Stewart,  _The symmetry perspective: from equilibrium to chaos in phase space and physical space_ , Vol. 200 (Springer Science & Business Media, 2003).
  * Hoyle and Hoyle [2006] R. Hoyle and R. B. Hoyle,  _Pattern formation: an introduction to methods_ (Cambridge University Press, 2006).
  * Krygier  _et al._ [2021] M. C. Krygier, J. L. Pughe-Sanford, and R. O. Grigoriev, Exact coherent structures and shadowing in turbulent taylor–couette flow, Journal of Fluid Mechanics 923 (2021).
  * Crowley  _et al._ [2022] C. J. Crowley, J. L. Pughe-Sanford, W. Toler, M. C. Krygier, R. O. Grigoriev, and M. F. Schatz, Turbulence tracks recurrent solutions, Proceedings of the National Academy of Sciences 119, e2120665119 (2022).
  * Marenduzzo  _et al._ [2007] D. Marenduzzo, E. Orlandini, M. E. Cates, and J. M. Yeomans, Steady-state hydrodynamic instabilities of active liquid crystals: Hybrid Lattice Boltzmann simulations, Phys. Rev. E 76, 10.1103/physreve.76.031921 (2007).
  * Thampi  _et al._ [2014b] S. P. Thampi, R. Golestanian, and J. M. Yeomans, Vorticity, defects and correlations in active turbulence, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 372, 20130366 (2014b).
  * Giomi  _et al._ [2014] L. Giomi, M. J. Bowick, P. Mishra, R. Sknepnek, and M. Cristina Marchetti, Defect dynamics in active nematics, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 372, 20130365 (2014).
  * Shaebani  _et al._ [2020] M. R. Shaebani, A. Wysocki, R. G. Winkler, G. Gompper, and H. Rieger, Computational models for active matter, Nature Reviews Physics 2, 181 (2020).
  * Koch and Wilczek [2021] C.-M. Koch and M. Wilczek, Role of advective inertia in active nematic turbulence, Phys. Rev. Lett. 127, 268005 (2021).
  * Blow  _et al._ [2017] M. L. Blow, M. Aqil, B. Liebchen, and D. Marenduzzo, Motility of active nematic films driven by “active anchoring”, Soft matter 13, 6137 (2017).
  * Doostmohammadi  _et al._ [2017] A. Doostmohammadi, T. N. Shendruk, K. Thijssen, and J. M. Yeomans, Onset of meso-scale turbulence in active nematics, Nature communications 8, 15326 (2017).
  * Wagner and Grover [2022] C. Wagner and P. Grover, https://github.com/DynamicalSystemsLab-UNL/ECSAct (2022).
  * Burns  _et al._ [2020] K. J. Burns, G. M. Vasil, J. S. Oishi, D. Lecoanet, and B. P. Brown, Dedalus: A flexible framework for numerical simulations with spectral methods, Phys. Rev. Res. 2, 023068 (2020).
  * Boyd [2001] J. P. Boyd,  _Chebyshev and Fourier Spectral Methods_ , 2nd ed. (Springer-Verlag, 2001).
  * Viswanath [2007] D. Viswanath, Recurrent motions within plane Couette turbulence, J. Fluid. Mech. 580, 339 (2007).
  * Dennis and Schnabel [1996] J. E. Dennis and R. B. Schnabel,  _Numerical Methods for Unconstrained Optimization and Nonlinear Equations_ (Society for Industrial and Applied Mathematics, 1996).
  * Saad and Schultz [1986] Y. Saad and M. H. Schultz, GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Comput. 7, 856 (1986).
  * Chandler and Kerswell [2013] G. J. Chandler and R. R. Kerswell, Invariant recurrent solutions embedded in a turbulent two-dimensional Kolmogorov flow, J. Fluid. Mech. 722, 554 (2013).
  * Duclos  _et al._ [2018] G. Duclos, C. Blanch-Mercader, V. Yashunsky, G. Salbreux, J.-F. Joanny, J. Prost, and P. Silberzan, Spontaneous shear flow in confined cellular nematics, Nature physics 14, 728 (2018).
  * Walton  _et al._ [2020] J. Walton, G. McKay, M. Grinfeld, and N. J. Mottram, Pressure-driven changes to spontaneous flow in active nematic liquid crystals, Eur. Phys. J. E 43, 10.1140/epje/i2020-11973-8 (2020).
  * Gibson  _et al._ [2008] J. Gibson, J. Halcrow, and P. Cvitanovic, Visualizing the geometry of state space in plane Couette flow, J. Fluid. Mech. 611, 107 (2008).
  * [71] Supplementary material is available at this link, https://osf.io/9cp54/.
  * Pringle  _et al._ [2012] C. C. Pringle, A. P. Willis, and R. R. Kerswell, Minimal seeds for shear flow turbulence: using nonlinear transient growth to touch the edge of chaos, Journal of Fluid Mechanics 702, 415 (2012).
  * Duguet  _et al._ [2013] Y. Duguet, A. Monokrousos, L. Brandt, and D. S. Henningson, Minimal transition thresholds in plane couette flow, Physics of Fluids 25, 084103 (2013).
  * Ross  _et al._ [2018] S. D. Ross, A. E. BozorgMagham, S. Naik, and L. N. Virgin, Experimental validation of phase space conduits of transition between potential wells, Physical Review E 98, 052214 (2018).
  * Grover and Ross [2009] P. Grover and S. Ross, Designing trajectories in a planet-moon environment using the controlled Keplerian map, Journal of guidance, control, and dynamics 32, 437 (2009).
  * Naik and Ross [2017] S. Naik and S. D. Ross, Geometry of escaping dynamics in nonlinear ship motion, Communications in Nonlinear Science and Numerical Simulation 47, 48 (2017).
  * Zhong  _et al._ [2018] J. Zhong, L. N. Virgin, and S. D. Ross, A tube dynamics perspective governing stability transitions: an example based on snap-through buckling, International Journal of Mechanical Sciences 149, 413 (2018).
  * Mohammed and Grover [2022] M. A. Mohammed and P. Grover, Phase space analysis of nonlinear wave propagation in a bistable mechanical metamaterial with a defect, Physical Review E 106, 054204 (2022).
  * Farazmand [2016] M. Farazmand, An adjoint-based approach for finding invariant solutions of navier–stokes equations, Journal of Fluid Mechanics 795, 278 (2016).
  * Azimi  _et al._ [2022] S. Azimi, O. Ashtari, and T. M. Schneider, Constructing periodic orbits of high-dimensional chaotic systems by an adjoint-based variational method, Physical Review E 105, 014217 (2022).
  * Yalnız  _et al._ [2021] G. Yalnız, B. Hof, and N. B. Budanur, Coarse graining the state space of a turbulent flow using periodic orbits, Physical Review Letters 126, 244502 (2021).
  * Page  _et al._ [2021] J. Page, M. P. Brenner, and R. R. Kerswell, Revealing the state space of turbulence using machine learning, Physical Review Fluids 6, 034402 (2021).
  * Page  _et al._ [2022] J. Page, P. Norgaard, M. P. Brenner, and R. R. Kerswell, Recurrent flow patterns as a basis for turbulence: predicting statistics from structures, arXiv preprint arXiv:2212.11886 (2022).
  * Page  _et al._ [2020] J. Page, Y. Dubief, and R. R. Kerswell, Exact traveling wave solutions in viscoelastic channel flow, Phys. Rev. Lett. 125, 154501 (2020).
  * Dubief  _et al._ [2020] Y. Dubief, J. Page, R. R. Kerswell, V. E. Terrapon, and V. Steinberg, A first coherent structure in elasto-inertial turbulence, arXiv preprint arXiv:2006.06770 (2020).
  * Suri  _et al._ [2019] B. Suri, R. K. Pallantla, M. F. Schatz, and R. O. Grigoriev, Heteroclinic and homoclinic connections in a Kolmogorov-like flow, Physical Review E 100, 013112 (2019).
  * Parker and Schneider [2022] J. P. Parker and T. M. Schneider, Invariant tori in dissipative hyperchaos, Chaos: An Interdisciplinary Journal of Nonlinear Science 32, 113102 (2022).
  * Parker  _et al._ [2023] J. P. Parker, O. Ashtari, and T. M. Schneider, Predicting chaotic statistics with unstable invariant tori, arXiv preprint arXiv:2301.10626 (2023).
  * Doostmohammadi and Yeomans [2019] A. Doostmohammadi and J. M. Yeomans, Coherent motion of dense active matter, The European Physical Journal Special Topics 227, 2401 (2019).
  * Zarei  _et al._ [2020] Z. Zarei, C. Joshi, M. Norton, M. Hagan, and S. Fraden, Role of annulus confinement on 2D active nematic behaviour, Bulletin of the American Physical Society 65 (2020).