# Spontaneous Self-Constraint in Active Nematic Flows

Louise C. Head  School of Physics and Astronomy, The University of Edinburgh,
Peter Guthrie Tait Road, Edinburgh, EH9 3FD, UK.  Claire Doré  Laboratoire
Gulliver, UMRS 7083, ESPCI Paris, PSL Research University, 75005 Paris, France
Ryan Keogh  School of Physics and Astronomy, The University of Edinburgh,
Peter Guthrie Tait Road, Edinburgh, EH9 3FD, UK.  Lasse Bonn  Niels Bohr
Institute, University of Copenhagen, Blegdamsvej 17, Copenhagen 2100, Denmark
Amin Doostmohammadi  Niels Bohr Institute, University of Copenhagen,
Blegdamsvej 17, Copenhagen 2100, Denmark  Kristian Thijssen  Niels Bohr
Institute, University of Copenhagen, Blegdamsvej 17, Copenhagen 2100, Denmark
Teresa López-León  Laboratoire Gulliver, UMRS 7083, ESPCI Paris, PSL Research
University, 75005 Paris, France  Tyler N. Shendruk  School of Physics and
Astronomy, The University of Edinburgh, Peter Guthrie Tait Road, Edinburgh,
EH9 3FD, UK.  e-mail<EMAIL_ADDRESS>
###### Abstract

Active processes drive and guide biological dynamics across scales—from
subcellular cytoskeletal remodelling, through tissue development in
embryogenesis, to population-level bacterial colonies expansion. In each of
these, biological functionality requires collective flows to occur while self-
organized structures are protected; however, the mechanisms by which active
flows can spontaneously constrain their dynamics to preserve structure have
not previously been explained. By studying collective flows and defect
dynamics in active nematic films, we demonstrate the existence of a self-
constraint—a two-way, spontaneously arising relationship between activity-
driven isosurfaces of flow boundaries and mesoscale nematic structures. Our
results show that self-motile defects are tightly constrained to viscometric
surfaces—contours along which vorticity and strain-rate balance. This in turn
reveals that self-motile defects break mirror symmetry when they move along a
single viscometric surface, in contrast with expectations. This is explained
by an interdependence between viscometric surfaces and bend walls—elongated
narrow kinks in the orientation field. Although we focus on extensile nematic
films, numerical results show the constraint holds whenever activity leads to
motile half-charge defects. This mesoscale cross-field self-constraint offers
a new framework for tackling complex 3D active turbulence, designing dynamic
control into biomimetic materials, and understanding how biological systems
can employ active stress for dynamic self-organization.

##  1 Introduction

Disorderly turbulent flow exists in many classes of fluids, and characterizing
inertial turbulence [1] remains an outstanding challenge due to chaotic flows
across spatial and temporal scales [2]. However, the challenges are compounded
in rheologically complex fluids because couplings between fields introduces
additional nonlinearities. In elastic [2], granular [3], magnetohydrodynamic
[4], quantum [5] and liquid crystaline [6] turbulence, velocity is strongly
coupled to other fields, some of which possess their own topologies with
associated defects [7, 8, 9]. There is growing evidence that many biological
systems spontaneously exhibit turbulence-like disorderly flow states that are
coupled to local orientation (nematic) fields [10, 11]. In these systems,
active turbulence [12] is accompanied by continual creation/annihilation of
topologically protected defects in the orientational field, including self-
propelled $+1/2$ defects in 2D [13] and disclination lines in 3D [14].

To simplify the complexities of active nematic flows, studies often prioritize
one of the fields. Namely, since self-motile defects have been clearly
identified as crucial to active turbulence, they are often viewed as
controlling their own dynamical evolution, with only perturbations—either due
to local deformations or neighboring defects—causing deviations from their
ideal trajectories [15, 16, 17]. From this perspective, hydrodynamic flows
follow directly from the governing defect configuration. The antithetical
approach is to integrate the orientational field dynamics directly into the
hydrodynamic stresses creating defect-free fluid dynamical models [18, 19,
20]. Such velocity-fixated studies focus on long-lived and spatially extended
structures in the velocity/vorticity fields and attempt to identify scaling
laws for the energy/enstrophy spectra [21]. Together, these two approaches
have made substantial progress towards understanding active turbulence, but
the crucial mesoscale bridge between them remains to be developed.

In this article, we reveal that non-linear coupling between flow and
orientational fields in active nematics leads to a strong, two-way,
spontaneous self-constraint. On the one hand, we report that self-motile
topological defects are tightly constrained to specific flow boundaries,
while, on the other hand, these surfaces are driven by mesoscale defect-
associated nematic deformations. Specifically, our results demonstrate that
self-motile $+1/2$ defects are found solely on specific isosurfaces of flow
boundaries identified as viscometric surfaces — contours where vorticity and
strain-rate balance. These surfaces dictate evolution of the self-motile
topological defects, which align with and move along the contours since there
is no deformation of the streamlines along these paths. We uncover how this
causes the defects to break their ideally predicted mirror symmetry and be
classified according to their handedness. However, the spontaneous self-
constraint is a co-dependency and we find that the viscometric surface, in
turn, follows the mesoscopic deformation network of the orientational field.
While we focus on 2D extensile active nematic turbulence, our conclusions hold
whenever activity leads to motile half-integer topological defects, such as in
vortex lattices, flow-tumbling liquid-crystals, contractile active stress,
friction-induced ordering, and 3D structures. The generality of these novel
perspectives may provide greater insight into a more universal understanding
of fluidic systems with non-linear field couplings and topologically protected
states, as well as suggest a mechanism by which biological morphology could be
dynamically protected.

##  2 Executive summary of experiments and simulations

To explore the spontaneous self-constraint between disclinations and flow
structures in active nematics, we employ a combination of experimental and
numerical techniques. Experimentally, we create an active nematic film by
self-assembling labelled filamentous microtubules bundles at oil-water
interface with kinesin motor clusters [22], which act as cross-linkers and
hydrolyse ATP suspended in the aqueous phase, fueling extensile active
stresses in the microtubule film (§ 8.1). Microtubule bundles are imaged via
confocal fluorescence microscopy, which allows both the nematic tensor field
$\bm{\bm{Q}}$ to be inferred through coherence-enhanced diffusion filtering
[23] and the velocity field $\vec{u}$ with an optical flow method.
Numerically, we simulate active nemato-hydrodynamics via a hybrid lattice
Boltzmann approach [24] (§ 8.2). We focus on 2D extensile nematics in steady-
state active turbulence, which give rise to disorderly flow profiles and
continuous defect creation/annihilation (Fig. 1 and Supplementary Video 1).

To explore the interplay between topology and flow, we calculate mesoscopic
structures in each field. Defects are readily apparent (Fig. 1), and their
positions and orientation are found from the nematic field $\bm{\bm{Q}}$ (§
8.3.1). Likewise, active flow fields have mesoscale structure with well-
studied vorticity scales, which are visualized using the
$\mathcal{Q}$-criterion [25]

$\displaystyle\mathcal{Q}$
$\displaystyle=\frac{1}{2}\left(\lVert\bm{\bm{\Omega}}\rVert^{2}-\lVert\bm{\bm{E}}\rVert^{2}\right),$
(1)

where $\bm{\bm{\Omega}}$ and $\bm{\bm{E}}$ are the vorticity and strain-rate
tensors (§ 8.2). Vorticity-dominated regions are identified by
$\mathcal{Q}>0$, while $\mathcal{Q}<0$ in strain-rate-dominated regions (§
8.3.2). This technique has been used to identify vortex structures in
confluent cell layers [26], but any quantitative interplay between
disclination positions, the $\mathcal{Q}$-criterion, and director structure
has yet to be fully understood.

##  3 Plus-half defects are constrained to viscometric surfaces

Figure 1:  Plus-half topological defects reside on isolines of
$\mathcal{Q}=0$. (a) Scanning confocal microscopy of experimental active film
(§ 8.1). Plus-half defects marked as green comet-shaped symbols and minus-half
defects by dark blue trefoil-shaped symbols (§ 8.3.1). Zero-isolines of
$\mathcal{Q}$-criterion are shown as solid lines, colored by the handedness of
the enclosed vortex—red for clockwise and blue for anti-clockwise (§ 8.3.2).
Scale bar: $100\text{\,}\mathrm{\SIUnitSymbolMicro m}$. (b) Numerical
simulations of active turbulence (§ 8.2). Circulation of each vortex colored
representing the vortex strength. Nematic director field plotted as grey line
field. (c) Probability distribution function (PDF) of nearest distances of
$\pm 1/2$ defects to $\mathcal{Q}=0$ boundary (§ 8.3.6). (d) Snapshots of pair
creation event in microtubule-kinesin film. (e) Experimental snapshots of pair
annihilation event. (f) Same snapshot as (b), but vorticity field is colored,
with grey arrows for the velocity field. Zero-vorticity contours are shown as
orange dash-dot lines (§ 8.3.3). Left: Examples where the vorticity is non-
zero (top) and zero (bottom) at the $+1/2$ defect position.

By partitioning flow domains into vorticity-dominated ($\mathcal{Q}>0$) and
strain-rate dominated ($\mathcal{Q}<0$) flow, we confirm the previously
reported qualitative observation that defects tend to be found near the edge
of vortices [25, 27]. However, by imaging $\mathcal{Q}=0$ isolines (§ 8.3.2),
we identify that $+1/2$ defects are always found precisely on the border where
$\mathcal{Q}=0$ (Fig. 1). We refer to the lines of $\mathcal{Q}=0$ as
viscometric surfaces, as these are points where the velocity gradient tensor
is singular and nilpotent indicating simple shear flow with no stretching or
elongation of streamlines (§ 8.3.2). In the experimental system (Fig. 1a;
Supplementary Video 1), roughly one hundred defects are in frame and all of
the $+1/2$ defects are located on viscometric surfaces (solid lines colored by
the handedness of the enclosed vortex) at all times. Viscometric surfaces
enclose vorticity-dominated regions, which exhibit a broad distribution of
sizes, shapes and circulation (Ext. Data Fig. 1). Simulations reveal
$\mathcal{Q}=0$ contours that have an associated $+1/2$ defect tend to enclose
vortices with larger circulation (Fig. 1b; Supplementary Video 2). Both
experiments and simulations demonstrate that, while $+1/2$ defects are
constrained to viscometric surfaces, $-1/2$ defects are not (Fig. 1c). Defect
pair creation/annihilation events occur on $\mathcal{Q}=0$ boundaries (Fig.
1d-e), with the motile $+1/2$ defect pinned to the $\mathcal{Q}=0$ viscometric
surface throughout the process, while the $-1/2$ defects are unbound from
$\mathcal{Q}=0$ lines.

Though previous studies observed that defects reside near the edges of
vortices [25, 27], the implications have not been quantified. Commonly, active
flow states are visualized by their vorticity field (Fig. 1f; Supplementary
Video 3). While these visualizations demonstrate the importance of vorticity
in the characteristic active length and time scales, they are not observed to
strictly constrain the configuration of defects. Identifing edges of vortices
as the point where the vorticity $\omega$ direction inverts (§ 8.3.3), it is
seen that most-but-not-all $+1/2$ defects exist on $|\omega|=0$ (Fig. 1f;
dash-dotted lines). While this agrees with previous statements that defects
tend to sit near the edges of vortices, our experiments and simulations
demonstrate that this is not a strict colocalization constraint, since not
every $+1/2$ defect sits on zero-vorticity contours (Fig. 1f; top inset).
Thus, it is not only where the vorticity is zero but rather where the
vorticity and strain rate balance to produce simple shear.

As a result, in both experimental and simulated results, $+1/2$ nematic
defects reside directly on this flow boundary (Fig. 1c). Likewise, $-1/2$
nematic defects are also most likely to be found on $\mathcal{Q}=0$; however,
this is because they are initially created at $\mathcal{Q}=0$ due to the
constraint on their $+1/2$ partner. In contrast to $+1/2$ defects, $-1/2$
defects are not constrained to remain on specific flow contours and can be
associated with any value of $\mathcal{Q}$, existing in both strain-rate- and
vorticity-dominated regions (Fig. 2a). Thus, they have a wider range in their
distance distribution, falling off exponentially with larger distances (Fig.
1c). Indeed, the ideal picture of a solitary $-1/2$ defect is not located on
$\mathcal{Q}=0$, but instead is encircled by a floweret of six vortices of
alternating handedness (Fig. 2b; § 8.4.1 [28]). As a result, $-1/2$ defects
reside on any sign of $\mathcal{Q}$-criterion (Fig. 2a; Supplementary Video
2). In contrast, the striking coincidence of $+1/2$ defects on $\mathcal{Q}=0$
flow contours suggests an intrinsic relationship between flow structures and
the trajectory of the self-propelled defects.

##  4 Spontaneous self-constraint violates the ideal, solitary view of defect-
generated flow

Figure 2:  Configurations of defects and viscometric surfaces. (a) Snapshot of
the flow geometry around minus and plus half nematic defects from the
numerical simulations, illustrating three minus half defects in a strain-
dominated region, in a vorticity-dominated region and on $\mathcal{Q}=0$
border. Arrows show instantaneous velocity field, dashed lines director field,
solid lines $\mathcal{Q}=0$ contours, with enclosed area colored by
circulation (red for clockwise and blue for anti-clockwise). Blue trefoil
symbols mark $-1/2$ defects and green comet-shaped symbols represent $+1/2$
defects. (b) Prediction for ideal, solitary $-1/2$ defect (§ 8.4.1). (c)-(e)
Snapshots of plus half defects, for the cases of (c) mirror symmetry, (d)
mirror-symmetry-broken with anti-clockwise vortex, and (e) clockwise. (f) Same
as (b) for solitary $+1/2$ defect. (g) Probability distribution function (PDF)
of alignment angles between the defect orientation and the tangent of the
associated viscometric ($\mathcal{Q}=0$) line (§ 8.3.6). The vertical blue
dash-dot line indicates the ideally expected alignment from (f) (§ 8.4.1) and
the grey dashed line indicates the alignment angle for a point active force (§
8.4.4).

Investigating the relationship between the viscometric ($\mathcal{Q}=0$)
surfaces and the $+1/2$ defect positions reveals that the defect/surface
complexes exist in two distinct configurations: Either a $+1/2$ defect is
positioned at an intersection of viscometric lines (Fig. 2c), or a defect lies
parallel to a single $\mathcal{Q}=0$ viscometric line (Fig. 2d-e). In the
first case (Fig. 2c), a $+1/2$ defect is positioned at a cross-roads between
two viscometric lines, and has two equally strong counter-rotating vortex
regions on either side of the defect axis, while the flow directly in front
and behind is strain-rate dominated. Thus, the defect has a mirror symmetry
along its head-tail axis. This corresponds to the ideal picture of a solitary
defect (Fig. 2f; § 8.4.1 [28]), which is tacitly the expectation of flows
around defects in active turbulence. The combination of these flow geometries
ensures the velocity field flows parallel to the defect orientation and self-
propulsion direction. However, this mirror-symmetric configuration is
transient. Although the ideal model of solitary defects describes how active
forces are a symmetric source of vorticity on both sides of the defects, in
experiments and simulations actual vorticity-dominated regions are not self-
propulsively moving along side self-motile defects as the picture might
suggest.

Instead, the mirror symmetry around the isolines is broken (Fig. 2d-e). Each
of the defects in these states orient parallel to a single viscometric line,
with vorticity dominating on one side and strain rate dominating on the other.
This spontaneous handedness of a subpopulation of defects is not predicted by
the ideal model for an solitary defect; however, defects with this
configuration persist for long durations and are observed to be the majority
(Fig. 1; Supplementary Video 2) because a motile defect can persistently move
along the $\mathcal{Q}=0$ line, circulating around a single neighboring
vortex. As a result, local chirality spontaneously emerges in this globally
achiral system. This spontaneous handedness is necessarily erased when
ensemble averaging fields around defect cores. The average flow profiles of
the two non-symmetric cases are distinct from the mirror-symmetric case when
the subpopulations are separately averaged (Ext. Data Fig. 2).

We quantify the distribution of these defect/viscometric surface complexes via
the distribution between the defect alignment angles between the $+1/2$ defect
axis and the $\mathcal{Q}=0$ boundary tangent (Fig. 2g; § 8.3.6). In
experiments and simulations, the most likely alignment is directly parallel,
quantifying the preference for mirror-symmetry-broken configurations in active
turbulence. This demonstrates that the mirror-symmetry-broken configuration is
actually preferred to the symmetric configuration predicted by the ideal
solitary-defect theory (§ 8.4.1). The simulations also show a secondary peak
around $\theta\approx 20$\mathrm{\SIUnitSymbolDegree}$$, representing the
minority-population in the transient symmetric configuration. The measured
angle is smaller than $\theta\approx 35.26$\mathrm{\SIUnitSymbolDegree}$$, the
ideal-model angle for a solitary defect (§ 8.4.1). The experimental alignment
angle distribution does not show a second peak, which may be due a combination
of the $\mathcal{Q}$-criterion boundaries being less smooth or the mirror-
symmetric case being even shorter-lived in experiments. Disclinations
intermittently transform between the mirror-symmetric and the broken-mirror-
symmetry cases (Supplementary Video 2). The results so far evidence that
specific isolines of flow boundaries constrain the motion of topological
defects but, in the next section, we show that there is a two-way
interdependence and explain by what mechanism defects and viscometeric
contours are pinned.

##  5 Bend walls and viscometric surfaces are interdependent

Figure 3:  Viscometric surfaces generated by bend walls. (a) Same snapshot as
Fig. 1b with colormap corresponding to the scalar-bend bend parameter
$S_{\text{SB}}$ (§ 8.3.5). Narrow lines of strong negative $S_{\text{SB}}$ are
interpreted as bend walls. (b) Probability distribution function (PDF) of
$S_{\text{SB}}$ sampled equally for all points and times, against the
distribution of $S_{\text{SB}}$ only on $\mathcal{Q}=0$ contours. (c)-(e)
$\mathcal{Q}$-criterion and velocity field predicted by Stokesian line force
model of bend walls (§ 8.4.4). Red and blue shading indicate clockwise and
anticlockwise vorticity-dominated regions, while the white regions are strain-
rate dominated. The $\mathcal{Q}=0$ contours are solid lines. The unit-vector
velocity field is shown with grey arrows. The line force is shown as thick
solid black line with the direction of forcing indicated by the purple arrows.
Line force taken as: (c) finite line, (d) full-circle and (e) half-circle.

Though topological defects are often identified as dominant localized-sources
of active forcing (§ 8.4.1), more generally force is generated from any
deformation of the orientation field (§ 8.2). In particular, bend walls play
an essential role in extensile active nematics [29]. Bend walls are narrow
lines of high bend deformation forming sharp kinked lines between nematic
domains, representing a nematic Néel wall [30]. As highly localized
deformations, bend walls generate strong active forcing and, because extensile
nematics are hydrodynamically unstable [31], bend constricts into system-
spanning narrow kink lines before defect pairs unbind along these walls and
the self-motile $+1/2$ defects unzip the bend wall as it advances [32]. After
this initial development of active turbulence, it is then presumed that
significant forces only occur in the immediate vicinity of defects. However,
in practice this is not observed in experiments or simulations. Rather, the
hydrodynamic instability is incessantly constricting bend into sharp kink
walls, as demonstrated by the splay-bend parameter $S_{\text{SB}}$, which
indicates the difference between splay and bend deformations (§ 8.3.5). The
splay-bend parameter traces bend walls through strongly negative
$S_{\text{SB}}$ values (Fig. 3a; Supplementary Video 4), and reveals that bend
walls correlate strongly with the $\mathcal{Q}=0$ viscometric lines. The
distribution of $S_{\text{SB}}$ coinciding with viscometric surfaces is more
skewed towards negative (bend) values than for the system as a whole (Fig.
3b).

###  5.1 Physical Model

The extended concurrence of bend walls and $\mathcal{Q}=0$ lines suggests the
physical mechanism by which defect dynamics and viscometric surfaces become
cross-constrained. We posit that the elongated bend walls are crucial to the
pure shearing flows of the viscometric surfaces. To understand how these
mesoscale nematic and hydrodynamic structures are intertwined, we consider a
series of simple models. These models will demonstrate that the bend walls
linking $\pm 1/2$ defects (Fig. 3a) explain both the viscometric contour and
the handedness of elastically bound $+1/2$ defects, as long as the walls are
(i) finite in length and (ii) follow curved paths.

We have already seen how the ideal analytical solution of a solitary $+1/2$
defect explains the mirror-symmetric case but not the mirror-symmetry-broken
case (Fig. 2). In fact, even modelling solitary defects as simple point-forces
is sufficient to qualitatively predict mirror symmetric defects at
$\mathcal{Q}=0$ intersections (§ 8.4.4). A straightforward model that accounts
for the elongated nature of bend walls might treat them as infinitely long,
perfectly straight structures (§ 8.4.2). Indeed, this is consistent with
$\mathcal{Q}=0$; however, this overly simplified model predicts
$\mathcal{Q}=0$ everywhere, rather than just on the bend wall. This is because
perfectly straight, infinitely long bend walls produce only shear.

To address this shortcoming, bend wall curvature is added to the model by
perturbing the bend walls into infinitely long sinusoidal undulations. This
model numerically demonstrates that curved bend walls produce closed
viscometric surfaces (Ext. Data Fig. 3). Furthermore, this model demonstrates
that active viscometric flows drive additional bend constriction, narrowing
the bend walls into kink lines but not further perturbing their conformation
(§ 8.4.3). This re-emphasizes an essential property: $\mathcal{Q}=0$ lines,
where vorticity and strain-rate balance, is also where velocity gradients
exhibit no stretching or elongation of equidistant streamlines. This is not
only why $\mathcal{Q}=0$ is referred to as viscometric, but also why the bend
walls are not further perturbed by the flows and why flows do not misalign
defects from following $\mathcal{Q}=0$ lines. Thus, this simple model suggests
that for the bend walls to be deformed, they cannot be infinitely long.

These results justifies an even simpler approach of modeling $+1/2$ defects
unzipping narrow finite-length kink walls as idealized lines of active force.
We assume active stresses dominate elastic stresses, and hence, model the
flows using Stokes equation (§ 8.4.4). The bend direction, and therefore
active force, is parallel with the tangent of the bend wall. By solving the
Stokes equation, we calculate the velocity field and $\mathcal{Q}=0$ contours
for the following cases: (i) a finite straight line, (ii) a perfect circle,
and (iii) an arc (Fig. 3c-e). Case (i) assumes a defect is the end of a finite
straight kink wall, which accentuates the acute angle (Fig. 3c; § 8.4.4). This
explains why the observed alignment angle (Fig. 2g) is substantially smaller
than the ideal solitary value. However, this bend wall conformation retains
the mirror symmetry of the flow field — in the limit of an infinitely long
line, the two $\mathcal{Q}=0$ lines collapse into a single line but
simultaneously $\mathcal{Q}$ goes to zero everywhere. Case (ii) instead
considers a perfectly circular bend wall (Fig. 3d) for which a $\mathcal{Q}=0$
contour encircles a vorticity-dominated core, with a strain-rate dominated
exterior. While completely circular bend walls are not observed in experiments
or simulations, this model produces a well-defined vortex, which is the basis
of the current understanding of energy spectra in active turbulence [25]. The
intermediate case (iii) of arced bend walls explains the symmetry-broken
regime in which many of the defects exist. Once defects are present in the
system, bend walls continuously curve (Fig. 3a). To simply model curved bend
walls, we numerically solve an arced force line (Fig. 3e). This model
reproduces two viscometric contours, but they no longer intersect. Rather the
viscometric line on the concave side of the arc continues to coincide with the
bend wall, while the $\mathcal{Q}=0$ contour on the convex side separates,
which is the physical picture observed in active turbulence (Fig. 1).

Hence, the straightforward model of a finite, curved force line explains how
bend walls play the crucial role in setting the interdependence of defect
dynamics and viscometric surfaces: First, bend constriction via the
hydrodynamic instability causes active forces to be localized along narrow
lines, resulting in shearing flows that neither stretch nor elongate
streamlines. Self-motile $+1/2$ defects unzip these narrow bend walls and thus
are elastically constrained to $\mathcal{Q}=0$. Finally, the finite length and
curvature of the bend walls causes the mirror-symmetry-broken case.

##  6 Beyond extensile, flow-aligning 2D active nematic turbulence

Figure 4:  The spontaneous self-constraint between motile $+1/2$ defects and
viscometric surfaces holds more generally. (a) Active extensile turbulence
with an alignment parameter in the flow-tumbling regime. (b) Contractile
active turbulence, in an operating regime with rare instances of bound $+1/2$
defect pairs creating effective $+1$ defects, which do not need to be on
$\mathcal{Q}=0$. Same visualization as Fig. 3a. (c) Lanes of alternating
circulation form with anisotropic friction in the flow-aligning regime.
(d)-(e) Dancing defect dynamics in 2D confining channel, with $+1/2$ defects
switching between vortices and at $\mathcal{Q}=0$ intersections in (d) and
aligned along vortex boundaries in (e). (a)-(e) show numerical simulations of
2D active nematics, with colormaps, lines and markers the same as Fig. 1b. (f)
Experimental verification that the spontaneous self-constraint is retained in
a channel confinement. Scale bar: $100\text{\,}\mathrm{\SIUnitSymbolMicro m}$.
(g) Three-dimensional vortex lattice with disclinatines colored by the twist
angle ($\cos\beta$ from § 8.3.1). Yellow indicates local $+1/2$ wedge profile,
green twist-type profile and blue $-1/2$ wedge profile. Grey shading
illustrates vorticity-dominate regions ($\mathcal{Q}>0$), and white where
strain-rate dominates ($\mathcal{Q}<0$).

This study has focused on the spontaneous self-constraint between defects and
flow structures in unconfined 2D extensile nematic turbulence. However, our
evidence suggests this phenomenon is substantially more general, apparently
holding whenever $+1/2$ topological defects are present in active nematics
(Fig. 4).

The conclusions arrived at for flow-aligning active turbulence also hold in
the flow-tumbling regime (Fig. 4a). Once again, defects reside on viscometric
contours, and either travel along a single vortex boundary or transiently
cross an intersection (Supplementary Video 5). The interdependence strongly
resembles the flow-aligning regime, which is consistent with the expectation
that active turbulence in flow-aligning and flow-tumbling behaves similarly.
Futhermore, the conclusions likewise generally hold for contractile activity
(Fig. 4b). Minus half defects are still non-motile and unbound from
viscometric lines. Motile $+1/2$ defects now chase their tail (move in the
opposite direction compared to extensile defects) and unzip splay walls; yet,
they are still most likely to be found self-propulsively moving where
$\mathcal{Q}=0$ (Supplementary Video 6). However in flow-aligning contractile
systems, there is an intriguing exception: In a narrow parameter regime, a
pair of similarly charged $+1/2$ defects can temporarily form impermanent $+1$
complexes [33]. These effectively immotile $+1$ complexes produce a flow
profile that does not necessitate $\mathcal{Q}=0$ and, thus, these effectively
$+1$ compounds do not reside where $\mathcal{Q}=0$, allowing them to briefly
elude the confinement condition. Thus, self-motility of $+1/2$ defects,
actively unzipping deformation walls, is central to the reported
interdependence.

How robust is the interdependence between $+1/2$ defects and $\mathcal{Q}=0$
to extrinsic conditions? Homogeneous friction, for example, can modify the
properties of active turbulence [34, 35] and anisotropic friction can order
flow states [36, 37]. Anisotropic friction creates an easy flow axis, forming
lanes with preferred circulation handedness encircled by $\mathcal{Q}=0$ (Fig.
4c; Supplementary Video 7). The vortex-dominated regions form a lattice of
viscometric contours, with intersection points along the boundary at advective
lanes. Defects still reside on $\mathcal{Q}=0$ and, moreover, now appear to be
predominantly found on intersections. This suggests a distinction between the
low and high friction regimes; defects are more likely to be in the mirror-
symmetry-broken configuration without friction, while these high friction
regimes combat bending curvature/instabilities, which our theoretical model
suggests will make mirror-symmetric defects more likely.

Impermeable, no-slip boundary conditions can also have profound effects on
spontaneous active flow states, creating complex spatio-temporally ordered
dynamics, even in simple geometries, such as 2D channel confinement [27, 38].
The resulting vortex lattice and associated defect dancing is clear in
simulations (Fig. 4(d,e); Supplementary Video 8), but also apparent in
experiments (Fig. 4f; Supplementary Video 9). In both simulations and
experiments, the $+1/2$ defects orient parallel to the viscometric line, when
traversing a vortex (Fig. 4e), and are able to cross to neighboring vortices
by passing through the fleeting $\mathcal{Q}=0$ intersection (Fig. 4d;
Supplementary Video 10). Similarly, in experiments $+1/2$ defects weave
through the center, never departing from $\mathcal{Q}=0$ and transversing to
neighboring vortices through viscometric intersections. These dancing dynamics
beautifully illustrate the transitions between the mirror-symmetric and
mirror-symmetry-broken configurations (Ext. Data Fig. 4). Plus-half defects
orient and travel parallel to the edge of a vortex, until the vortices
instantaneously touch forming a short-lived $\mathcal{Q}=0$ intersection,
allowing the defects to transfer to the next vortex of opposite handedness—all
the while, the defects remain on the $\mathcal{Q}=0$ lines.

In fact, the interdependent spontaneous self-constraint holds in three-
dimensions, where the point defects become disclination lines and viscometric
contours become 2D surfaces. The vortex lattice exists in a 3D duct [39].
Analogously to the two-dimensional case, wedge-type $-1/2$-profile
disclination lines exist at the walls, while disclinations with a wedge-type
$+1/2$ profile dance around the central vortex lattice (Fig. 4g). The
disclinations at the center with wedge-type $+1/2$ profiles continually
maintain contact with the two-dimensional $\mathcal{Q}=0$ isosurface
(Supplementary Video 11), demonstrating that motile defects are constrained to
lie on $\mathcal{Q}=0$ surfaces even in 3D. Altogether, the examples explored
in Fig. 4 indicate that the self-emergent constraint between motile defect
dynamics and viscometric surfaces is not only a property of extensile active
turbulence but is more general, apparently holding whenever motile defects are
present in active nematics.

##  7 Conclusion

In conclusion, we have identified and explained a spontaneously arising
constraint between motile defect dynamics and viscometric surfaces, where
vorticity and strain-rate balance. This work challenges the idea that nematic
defects are solely responsible for their own dynamics — ultimately, neither
defects nor hydrodynamics alone govern the multi-field dynamics that
spontaneously arise in active nematics; rather, they are intrinsically
interdependent. Identifying this spontaneous self-constraint revealed that the
ideal picture of solitary self-motile defects generating a pair of mirror-
symmetric vortices is not strictly true in practice. Instead, motile defects
more often move along a single viscometric surface, where there is no
stretching nor elongation of the fluid, a fact that appears to hold for all
motile defects in active nematics and not only 2D extensile active turbulence.
This shows that defects can be classified into three conformations based on
local handedness with respect to their viscometric surface surroundings. Until
now, this variation, as well as any possible different properties, had been
negelected. Not only do the hydrodynamics not force the plus-half defects off
these contours, but they coincide with the very bend walls that these defects
are unzipping — the field dynamics are codependent. Furthermore, these results
underscore the continual role of bend wall constriction and unzipping in
steady state dynamics, and highlight the centrality of mesoscale structure in
collective dynamics, providing a mechanism to manipulate the degree of orderly
dynamics in a system.

We anticipate that our observation of co-dependence will aid efforts to
understand the highly complicated structure of 3D active turbulence [40], and
complement recently developed Lagrangian descriptions of coherent structures
[41]. The self-constraint identified here could lead to designs for novel
active functionality [42, 43] or biointerfaces [44]. Biological systems,
including active transport in microbial colonies [45], mitotic spindle
assembly [46] tissue responses [47], epithelial reorganisation [48] and
morphogenesis [49], all involve a careful balance between collective motion
and protected-structures for their functionality. Ultimately, spontaneous
constraints in active materials could have far-reaching implications as a
framework for understanding how biological systems employ active stresses for
simultaneous dynamics and restraint.

##  8 Method section

###  8.1 Experiments

####  8.1.1 Materials

Short and stable microtubules (MTs) were polymerized by incubating, at
$37\text{\,}\mathrm{\SIUnitSymbolCelsius}$ for 30 minutes, a mixture
containing $8\text{\,}\mathrm{mg}\text{\,}{\mathrm{mL}}^{-1}$ of recycled
tubulin from bovine brain (Brandeis Materials Research Science and Engineering
Center), $6\text{\,}\mathrm{mM}$ of
guanosine-5’-[($\alpha$,$\beta$)-methyleno]triphosphate (GMPCCP, Jena
Biosciences, NU-405S), and $1\text{\,}\mathrm{mM}$ of DL-dithiothreitol (DTT,
Sigma, 43815) in M2B buffer ($80\text{\,}\mathrm{mM}$ of PIPES [Sigma, P1851],
$1\text{\,}\mathrm{mM}$ of EGTA [Sigma, E3889], $2\text{\,}\mathrm{mM}$ of
MgCl2, pH=6.8). After incubation, the solution was kept at room temperature
for 5 hours, then frozen in liquid nitrogen and stored at
$-80\text{\,}\mathrm{\SIUnitSymbolCelsius}$ for future use. For confocal
fluorescence microscopy, 3% of the tubulin was labelled with Alexa-647, a
bright, far-red–fluorescent dye. Biotinylated kinesin was prepared in the
BioNMR group at the institute of bioengineering of Catalunya (IBEC).
Drosophila melanogaster heavy-chain kinesin-1 K401-BCCP-6His was expressed in
Escherichia coli using the plasmid WC2 from the Gelles Laboratory (Brandeis
University). Biotinylated kinesins were purified with a nickel column and
dialysis against $500\text{\,}\mathrm{mM}$ imidazole aqueous buffer. The
kinesin concentration was estimated at $2.5\text{\,}\mathrm{\SIUnitSymbolMicro
M}$ by absorption spectroscopy, and was finally stored in a 40%w/v sucrose
solution at $-80\text{\,}\mathrm{\SIUnitSymbolCelsius}$ for future use.
Kinesin clusters were obtained by preparing a mixture with
$1\text{\,}\mathrm{\SIUnitSymbolMicro M}$ of biotinylated kinesin dimers and
$31\text{\,}\mathrm{\SIUnitSymbolMicro g}\text{\,}{\mathrm{mL}}^{-1}$ of
tetrameric streptavidin in a M2B buffer supplemented with
$0.2\text{\,}\mathrm{mM}$ of DTT, this stoichiometric ratio corresponding
approximately to 2 kinesin dimers per cluster. This mixture was incubated on
ice for 30 minutes.

The active microtubule/kinesin gel consisted in a M2B preparation with
$1.6\text{\,}\mathrm{mg}\text{\,}{\mathrm{mL}}^{-1}$ of MTs, motor clusters
(the concentration of streptavidin was $8.2\text{\,}\mathrm{\SIUnitSymbolMicro
g}\text{\,}{\mathrm{mL}}^{-1}$), $1.4\text{\,}\mathrm{mM}$ of adenosine
triphosphate (ATP, Sigma, A2383) to fuel the molecular motors and 1.6% w/v of
depleting agent polyethylene glycol (PEG, 20 kDa, Sigma, 9517) to induce the
bundling of the MTs. An ATP-regenerating system (2.8% v/v of pyruvate
kinase/lactate dehydrogenase [PK/LDH, Sigma, P02] and
$26.2\text{\,}\mathrm{mM}$ of phosphoenolpyruvate [PEP; Sigma; P7127]) was
employed to keep a fixed concentration of ATP during the experiments. To
prevent photobleaching and oxidation, some antioxydants were included in the
active gel: $5.4\text{\,}\mathrm{mM}$ of DTT,
$3.3\text{\,}\mathrm{mg}\text{\,}{\mathrm{mL}}^{-1}$ of glucose (Sigma;
G8270), $38\text{\,}\mathrm{\SIUnitSymbolMicro g}\text{\,}{\mathrm{mL}}^{-1}$
of catalase (Sigma, C40),
$0.22\text{\,}\mathrm{mg}\text{\,}{\mathrm{mL}}^{-1}$ of glucose oxidase
(Sigma G2133) and $2.0\text{\,}\mathrm{mM}$ of trolox (Sigma, 238813). The
MgCl2 concentration was raised with 4.7% v/v of Mix solution
($69\text{\,}\mathrm{mM}$ of MgCl2 in M2B). The pH of the active gel was
adjusted to 6.8.

####  8.1.2 Active nematic assembly

Active nematic layers were assembled at a water/oil interface following two
different protocols. The unconfined active nematic layer (as shown in Fig. 1a)
was prepared in a closed observation chamber, consisting of a bottom
hydrophobic glass slide treated with Aquapel and a hydrophilic cover slip
coated with a polyacrylamide brush, separated by
$120\text{\,}\mathrm{\SIUnitSymbolMicro m}$ of double-sided tape, forming a
rectangular chamber of size
$3\text{\,}\mathrm{mm}$$\times$$22\text{\,}\mathrm{mm}$. The chamber was first
filled with fluorinated oil HFE7500 with 1.8% v/v of 008-FluoroSurfactant (RAN
Biotechnologies), then $10\text{\,}\mathrm{\SIUnitSymbolMicro L}$ of the
active gel was introduced by capillarity in the chamber, replacing most of the
oil but letting a lubricating film of oil at the surface of the bottom glass
slide. The chamber was sealed with UV-curable glue (Norland NOA-81). The
assembly of the microtubule bundles at the water/oil interface was
subsequently sped up by centrifuging the sample 10 minutes at 215$\times$g.
For laterally confined active nematics (as shown in Fig. 4f), we prepared an
open observation chamber by sticking with UV-curable glue circular
polydimethylsiloxane (PDMS) walls ($1\text{\,}\mathrm{cm}$ in diameter) to a
polyacrylamide coated glass slide. $2.5\text{\,}\mathrm{\SIUnitSymbolMicro L}$
of the active gel supplemented with 2.2% w/v of surfactant pluronic F-127
(Sigma, P-2443), was deposited in the well and spread on the hydrophilic glass
surface, then was immediately covered with
$200\text{\,}\mathrm{\SIUnitSymbolMicro L}$ of silicone oil with a dynamic
viscosity of $20\text{\,}\mathrm{mPa}\text{\,}\mathrm{s}$. The assembly of the
microtubules bundles at the water/oil interface occurred on their own within
40 minutes.

####  8.1.3 Micro-printed grid walls

To confine laterally the active nematic to channels, as shown in Fig. 4f, we
used millimeter-sized platforms in photoresist encompassing the rectangular
enclosures. Such grids were introduced in the open observation chamber from
the top and placed at the water/oil interface, consequently trapping the AN in
the rectangular channels [38]. Grids were 3D printed with a two-photon
polymerization printer, a Nanoscribe GT Photonic Professional device, with a
negative-tone photoresist IP-S (Nanoscribe GmbH, Germany) and a $25\times$
objective. The grids were directly printed on silicon substrates without any
preparation to avoid adhesion of the resist to the substrate. After developing
30 minutes in Propylene Glycol Monomethyl Ether Acetate (PGMEA 99.5%, Sigma,
484431) and 5 minutes in isopropanol, a batch polymerization was performed
with UV exposure (5 minutes at 80% of light power). The pattern resolution
achieved was about $50\text{\,}\mathrm{nm}$. The grid thickness was
$100\text{\,}\mathrm{\SIUnitSymbolMicro m}$, to ensure a good resistance
during manipulation.

####  8.1.4 Imaging and image processing

Images of the active nematic layer were acquired with a Nikon eclipse Ti-E
inverted microscope, equipped with a confocal spinning disk CSU-X1 (Yokogawa)
and a motorized stage Nikon Ti-S-EJOY, and a $10\times$ objective (Nikon). The
sample was excited with a laser beam at $647\text{\,}\mathrm{nm}$ with an
exposure between 200-500 $\mathrm{ms}$ depending on the sample fluorescence.
Images were captured with a camera Hamamatsu Orca Flash4.0, and acquired with
the software NiS-Elements. The acquisition frame rate was 2 frame per second.

Image processing, including background correction, was carried out with the
software ImageJ. The nematic tensor $\bm{\bm{Q}}\left(\vec{r},t\right)$ was
computed from the images using custom Matlab scripts [23] that infer local
alignment of the microtubules using coherence-enhanced diffusion filtering.
First, noise in the image is filtered out with a Gaussian blur of standard
deviation $\sigma_{1}$. The pixel-level orientation is obtained by finding the
direction along which the fluorescence intensity is most homogeneous and
smoothed by means of a Gaussian blur of standard deviation $\sigma_{2}$.
Finally, the nematic tensor $\bm{\bm{Q}}$ is constructed by ensemble averaging
the molecular orientation tensor over a small circle of size $\beta$ around
each pixel. The three parameters $\sigma_{1}$, $\sigma_{2}$ and $\beta$ were
manually adjusted by visual inspection of the resulting director field.
Velocity fields $\vec{u}\left(\vec{r},t\right)$ were calculated with the
optical flow method detailed in ref.[50] and associated Matlab codes, using
the ’classic+nl-fastp’ method.

The high resolution images exhibit pixel-level substantial noise, which
hinders the calculation of gradients of the fields, required for the
identification of defects and viscometric surfaces. Therefore, we generate a
coarse-grained field which only takes the value once every chosen number of
pixels. The same resolution was also used for the velocity fields for the
velocity-gradient based analysis. We test the ideal resolution using the
defect tracker (§ 8.3) and choose a value that identifies the most number of
defects correctly (while avoiding overestimating). For both bulk active
turbulence and channel confinement, the chosen resolution was 16.

###  8.2 Simulations

Our active nematic simulations employ a hybrid lattice Boltzmann approach [51]
to reproduce the continuum nematohydrodynamics model [52]. The dynamics of the
orientational order parameter
$\bm{\bm{Q}}\left(\vec{r},t\right)=S\left(\vec{n}\otimes\vec{n}-\bm{\bm{\delta}}/3\right)$
at each position $\vec{r}$ and time $t$ in 3D is described by the Beris-
Edwards transport equation

$\displaystyle\mathcal{D}_{t}\bm{\bm{Q}}$ $\displaystyle=\Gamma\bm{\bm{H}},$
(2)

where the right-hand side is the relaxation to equilibrium through a
relaxation parameter $\Gamma$ and molecular field
$\bm{\bm{H}}=-\left(\frac{\delta
F}{\delta\bm{\bm{Q}}}-\frac{1}{d}\bm{\bm{\delta}}\;\mathrm{tr}\left[\frac{\delta
F}{\delta\bm{\bm{Q}}}\right]\right)$ is the molecular field corresponding to
the nematic free energy $F$. On the left-hand side,
$\mathcal{D}_{t}\bm{\bm{Q}}=\left(D_{t}-R\right)\bm{\bm{Q}}$ is a co-variant
derivative that accounts for material advection by
$D_{t}=\partial_{t}+\vec{u}\cdot\vec{\nabla}$, which couples directly to the
velocity field, and a corotational operator
$R\left\\{\bm{\bm{Q}}\right\\}=\left\\{\lambda\bm{\bm{E}}+\bm{\bm{\Omega}},\tilde{\bm{\bm{Q}}}\right\\}_{+}-2\lambda\left(\bm{\bm{Q}}:\bm{\bm{L}}\right)\tilde{\bm{\bm{Q}}}$.
The co-rotational term couples the orientation to the velocity gradient
$\bm{\bm{L}}=\vec{\nabla}\otimes\vec{u}$ through the vorticity
$\bm{\bm{\Omega}}=\left(\bm{\bm{L}}-{\bm{\bm{L}}}^{T}\right)/2$ and strain-
rate $\bm{\bm{E}}=\left(\bm{\bm{L}}+{\bm{\bm{L}}}^{T}\right)/2$ tensors,
written using a (anti)commutator
$\left\\{\bm{\bm{A}},\bm{\bm{B}}\right\\}_{\pm}=\bm{\bm{A}}\cdot\bm{\bm{B}}\pm\bm{\bm{B}}\cdot\bm{\bm{A}}$,
with a flow-aligning parameter $\lambda$ and
$\tilde{\bm{\bm{Q}}}\equiv\bm{\bm{Q}}+\bm{\bm{\delta}}/d$. The Beris-Edwards
equation (Eq. (2)) is numerically solved by finite difference.

The velocity $\vec{u}\left(\vec{r},t\right)$ evolves according to the
generalized Navier-Stokes equations

$\displaystyle\rho D_{t}\vec{u}$
$\displaystyle=\vec{\nabla}\cdot\bm{\bm{\Pi}}-\bm{\bm{\kappa}}\cdot\vec{u}$
(3)

with density $\rho$ and friction tensor $\bm{\bm{\kappa}}$. The stress
$\bm{\bm{\Pi}}$ consists of a viscous term $2\eta\bm{\bm{E}}$ with viscosity
$\eta$, a pressure term $-p\bm{\bm{\delta}}$, an elastic term
$K\left[2\lambda\left(\bm{\bm{Q}}:\bm{\bm{H}}\right)\tilde{\bm{\bm{Q}}}-\lambda\left\\{\bm{\bm{H}},\tilde{\bm{\bm{Q}}}\right\\}_{+}+\left\\{\bm{\bm{Q}},\bm{\bm{H}}\right\\}_{-}-\vec{\nabla}\left(\bm{\bm{Q}}:\tfrac{\delta
F}{\delta\vec{\nabla}\bm{\bm{Q}}}\right)\right]$, and, most relevantly an
active term,

$\displaystyle\bm{\bm{\Pi}}_{\zeta}$ $\displaystyle=-\zeta\bm{\bm{Q}}.$ (4)

Both the elastic and active terms couple to velocity field to the nematic
field, but the active stress is directly proportional to the nematic field,
with an activity coefficient $\zeta$ that leads to a bend instability when
$\zeta>0$ (or a splay instability when $\zeta<0$), which drives deformation,
wall formation and defect unbinding, leading to the active defect turbulence.
The lattice Boltzmann algorithm simulates Eq. (3) and achieves near
incompressiblity $\vec{\nabla}\cdot\vec{u}\approx 0$.

The free energy includes a bulk component and a surface anchoring term
$F=\int_{V}f\ dV+\int_{A}f_{A}\ dA$ The bulk free energy density is
constructed as a Landau-de Gennes expansion with an one-constant elastic
Oseen-Frank term

$\displaystyle f$
$\displaystyle=A_{0}\left[\left(\frac{1-\gamma/3}{2}\right)\bm{\bm{Q}}:\bm{\bm{Q}}-\frac{\gamma}{3}\bm{\bm{Q}}:\left(\bm{\bm{Q}}\cdot\bm{\bm{Q}}\right)+\frac{\gamma}{4}\left(\bm{\bm{Q}}:\bm{\bm{Q}}\right)^{2}\right]+\frac{K}{2}\left(\vec{\nabla}\otimes\bm{\bm{Q}}\right)\vdots\left(\vec{\nabla}\otimes\bm{\bm{Q}}\right)$
(5)

with $\gamma$ controlling the distance from the isotropic-nematic phase
transition and $K$ the isotropic elastic constant. The surface free energy
density is

$\displaystyle f_{A}$
$\displaystyle=W_{1}\left(\bm{\bm{Q}}-\bm{\bm{Q}}^{\perp}\right):\left(\bm{\bm{Q}}-\bm{\bm{Q}}^{\perp}\right)+W_{2}\left(\bm{\bm{Q}}:\bm{\bm{Q}}-S_{\text{eq}}^{2}\right),$
(6)

where $W_{1,2}$ are anchoring energies, $6S_{\text{eq}}=1+3\sqrt{1-8/3\gamma}$
is the equilibrium scalar order parameter and
$\bm{\bm{Q}}^{\perp}=\bm{\bm{P}}\cdot\left(\bm{\bm{Q}}+S_{\text{eq}}\bm{\bm{\delta}}/3\right)\cdot\bm{\bm{P}}$
for the surface normal projection operator $\bm{\bm{P}}$.

Except where otherwise stated, our 2D extensile nematic simulations use the
parameters (in lattice Boltzmann units) density $\rho=1$, rotational
diffusivity $\Gamma=0.3375$, flow-aligning parameter $\lambda=1$, bulk energy
density scale $A_{0}=1$, distance from the isotropic-nematic transition
$\gamma=3$, dynamic viscosity $\eta=4/3$, Frank elasticity $K=0.05$, friction
$\bm{\bm{\kappa}}=\bm{\bm{0}}$ and extensile activity $\zeta=0.1$. Periodic
boundary conditions are used in a system of $80\times 80$. The velocity and
director fields are randomly initialized. A warmup time of $T_{W}=1\times
10^{4}$ is used before data is collected for $T=2\times 10^{4}$ lattice
Boltzmann times steps.

In Fig. 4 we show that the conclusions arrived at for 2D extensile active
turbulence hold under other conditions. As far as possible, we attempt to keep
our simulation parameters unchanged from the extensile case and only make a
few substitutions.

Flow tumbling (Fig. 4a; Supplementary Video 5)

Substitute $\lambda=0.6$.

Contractile (Fig. 4b; Supplementary Video 6)

Substitute $\zeta=-0.1$ and system size $160\times 160$. We note that these
contractile simulations permit a small non-zero component to the out-of-plane
director and velocity field. This escape into the third dimension is
negligible everywhere except for a small region at the top of Supplementary
Video 6, which persisted in time but did not grow.

Anisotropic friction (Fig. 4c; Supplementary Video 7)

We substitute
$\bm{\bm{\kappa}}=\left[\left[0.015,0\right],\left[0,0\right]\right]$.

2D channel confinement (Fig. 4d-e; Supplementary Video 8)

Periodic boundary conditions kept at the ends of the channels but the channel
walls are impermeable, no-slip boundaries with anchoring energy
$W_{1}=W_{2}=0.3$. System dimensions $130\times 25$.

3D channel confinement (Fig. 4g; Supplementary Video 11)

Same as the 2D channel wall but now in a 3D square duct.

Ext. Data Fig. 2

We substitute density $\rho=40$, rotational diffusivity $\Gamma=0.05$, dynamic
viscosity $\eta=1/6$, activity $\zeta=0.05$, system size $256\times 256$,
warmup time $T_{W}=5\times 10^{2}$ and data collection time $T=1.95\times
10^{4}$ lattice Boltzmann time steps.

###  8.3 Analysis

####  8.3.1 Defect Analysis

Topological disclinations in both experiments and 2D simulations are
identified by calculating the topological charge density [53, 54, 55]

$\displaystyle q$
$\displaystyle=\frac{1}{4\pi}(\partial_{x}Q_{x\alpha}\partial_{y}Q_{y\alpha}-\partial_{x}Q_{y\alpha}\partial_{y}Q_{x\alpha}),$
(7)

using summation convention on Greek indices. Points where $q$ is greater than
(less than) $1.5$ ($-1.5$) standard deviations are identified as $k=+1/2$
($-1/2$) defects. Once each 2D defect is identified, its orientation is found
via the angle

$\displaystyle\psi$
$\displaystyle=\frac{k}{k-1}\arctan\left(\frac{\langle\mathrm{sgn}\left[k\right]\partial_{x}Q_{xy}-\partial_{y}Q_{xx}\rangle_{\circ}}{\langle\partial_{x}Q_{xx}+\mathrm{sgn}\left[k\right]\partial_{y}Q_{xy}\rangle_{\circ}}\right),$
(8)

with $\langle\cdot\rangle_{\circ}$ denoting averaging a closed loop around the
disclination [56]. In 3D, defect lines are identified and characterised using
the disclination density tensor [53], $\bm{\bm{D}}$, constructed from
gradients in the nematic $\bm{\bm{Q}}$-tensor. In index notation,
$D_{ij}=\epsilon_{i\mu\nu}\epsilon_{jlk}\partial_{l}Q_{\mu\alpha}\partial_{k}Q_{\nu\alpha}$.
The tensor is interpreted as $\bm{\bm{D}}=s\vec{\mathcal{R}}\otimes\vec{T}$,
in terms of a non-negative scalar field $s$, the rotation vector
$\vec{\mathcal{R}}$, and tangent vector $\vec{T}$ along the disclination. The
values of $s$ approach a maximum at defect cores. We identified the defect
lines as isosurfaces where $s=0.1$. We characterise the local geometry of the
disclination line through the angle between the rotation and tangent vectors
$\cos\beta=\vec{\mathcal{R}}\cdot\vec{T}$, where $\beta$ is known as the twist
angle. Values of $\cos\beta$ smoothly transition between $-1$ for minus-half
wedge profiles, through to $+1$ for plus-half wedge profiles. Visualisations
in 3D used the Mayavi python library [57].

####  8.3.2 $\mathcal{Q}$-criterion

We characterize the principle behavior of flow fields by the
$\mathcal{Q}$-criterion, which is defined to be the second invariant of the
velocity gradient $\bm{\bm{L}}=\vec{\nabla}\otimes\vec{u}$ [58]

$\displaystyle\mathcal{Q}$
$\displaystyle=\frac{1}{2}\left(\left(\mathrm{tr}\left[\bm{\bm{L}}\right]\right)^{2}-\mathrm{tr}\left[\bm{\bm{L}}^{2}\right]\right).$
(9)

In the incompressibile limit,
$\mathrm{tr}\left[\bm{\bm{L}}\right]=\vec{\nabla}\cdot\vec{u}=0$ and so this
definition reduces to Eq. (1)), which we repeat here for convenience

$\displaystyle\mathcal{Q}$
$\displaystyle=\frac{1}{2}\left(\lVert\bm{\bm{\Omega}}\rVert^{2}-\lVert\bm{\bm{E}}\rVert^{2}\right),$
(10)

where the squared norms are $\lVert\bm{\bm{A}}\rVert^{2}=A_{ij}A_{ij}$. The
$\mathcal{Q}\left(\vec{r},t\right)$ is shown in Supplementary Video 12 and
iso-surfaces of $\mathcal{Q}=0$ are identified using the skimage python
library [59].

If $\mathcal{Q}>0$, vorticity dominates over strain-rate, where as if
$\mathcal{Q}<0$, strain-rate dominates. Where $\mathcal{Q}=0$, the magnitudes
of vorticity and strain-rate are equal. In 2D, isolines of $\mathcal{Q}=0$
must form closed loops, and, in 3D, $\mathcal{Q}=0$ contours must form closed
iso-surfaces enclosing vorticity dominated regions. In 2D, there are only two
invariants of $\bm{\bm{L}}$, and the first
$\mathrm{tr}\left[\bm{\bm{L}}\right]=\vec{\nabla}\cdot\vec{u}=0$ due to
incompressibility. Therefore, if $\mathcal{Q}=0$, all the invariants of
$\bm{\bm{L}}$ are zero, which indicates all the eigenvalues are zero. This in
turn means the velocity gradient tensor is singular and nilpotent. Simple
shear flow is an example, as is any viscometric flow in which streamlines are
equi-distance apart because there is no stretching or elongation of the
fluid[60, 61]. Thus, we refer to $\mathcal{Q}=0$ contours as viscometric
surfaces.

####  8.3.3 Comparing zero-isosurfaces $\mathcal{Q}$ with vorticity

Throughout this work, we employ the zero-isolines of $\mathcal{Q}$-criterion
to identify viscometric surfaces as the boundaries between flow structures. It
is a fact that the magnitude of the symmetric and antisymmetric velocity
gradients balance when $\mathcal{Q}=0$ to give finite pure shear. However, it
is also a possibility that the flow has completely vanishing gradients in the
velocity. To compare these two possibilities of the $\mathcal{Q}$-criterion,
we explore contours where contributions from $\bm{\bm{\Omega}}$ go to zero. In
any two-dimensional coordinate basis, $\bm{\bm{\Omega}}$ only has off diagonal
components $\Omega_{21}=-\Omega_{12}=\omega$ in terms of the vorticity pseudo-
vector $\vec{\omega}$. Hence, locations where the rotation-rate tensor
vanishes can be simply obtained from a sign change in $\omega$. Obtaining
contours where the strain rate vanishes is less direct since the tensorial
components vary with coordinate basis as there is no clear handedness to this
object and no pseudovector can be constructed in 2D. At locations where the
vorticity zero-lines and $\mathcal{Q}$-criterion zero-lines coincide, the
strain-rate contours must be zero by incompressibility. Otherwise, these
viscometric lines are where there is finite shear from the balance between
vorticity and strain-rate.

####  8.3.4 Distribution of vortex characteristics

Since viscometric surfaces form the boundary encircling vortex-dominated
domains from strain-rate-dominated regions, we quantify the distribution of
fluid elements composing each vortex domain (where $\mathcal{Q}>0$). After
identifying the list of points defining the $\mathcal{Q}=0$ boundary, vortices
are identified as the enclosed regions where $\mathcal{Q}>0$. A clustering
algorithm is applied that identifies the constituent vortex points. Vortex
area distributions are found as the total number of lattice Boltzmann nodes or
experimental pixels identified inside a viscometric surface with
$\mathcal{Q}>0$. We calculate the square of the relative shape anisotropy
$\kappa$ of each vortex, written in terms of the first ($I_{1}$) and second
($I_{2}$) invariants of the moment of inertia tensor [62], in $d$ dimensions
as

$\displaystyle\kappa^{2}$
$\displaystyle=1-\frac{2d}{d-1}\frac{I_{2}}{I_{1}^{2}}.$ (11)

The enclosed circulation is defined as $\Gamma=\oint_{c}\vec{u}\cdot
d\vec{l}=\oiint\vec{\omega}\cdot d\vec{S}$. We apply the second definition to
simplify the handling of non-zero genus vortices. The circulation distribution
in Fig. 1 is non-dimensionalized to compare the experimental and simulated
distributions by dividing by $\Gamma_{0}=\omega_{0}\ell_{\zeta}^{2}$. The
vorticity scale, $\omega_{0}$, was identified as the standard deviation of the
vorticity distribution. We calculate the active length scale as

$\displaystyle\ell_{\zeta}$ $\displaystyle=\sigma^{-1/2},$ (12)

where the defect number density
$\sigma=\left(N_{+1/2}+N_{-1/2}\right)/\left(TL^{2}\right)$ for $N_{\pm 1/2}$
nematic defects, total timesteps $T$ and system size $L$.

####  8.3.5 Splay-bend parameter

The splay-bend parameter [37, 63]

$\displaystyle S_{\text{SB}}$ $\displaystyle=\partial_{i}\partial_{j}Q_{ij}$
(13)

is used to visualize bend walls as line-like structures with large negative
values [37]. To interpret $S_{\text{SB}}$, consider a uniaxial 2D nematic with
constant scalar order $S$. Under these conditions
$S_{\text{SB}}=S\vec{\nabla}\cdot\left[\vec{n}\left(\vec{\nabla}\cdot\vec{n}\right)-\vec{n}\times\left(\vec{\nabla}\times\vec{n}\right)\right]$.
The first term inside the divergence is the splay and the second term is the
bend, leading to the interpretation that $S_{\text{SB}}$ represents the
difference between the divergence of splay and bend. Because the active force
density is the divergence of the active stress (Eq. (4))

$\displaystyle\vec{f}_{\zeta}$
$\displaystyle=-\zeta\vec{\nabla}\cdot\bm{\bm{Q}},$ (14)

in active nematics the splay-bend parameter can be understood to be
proportional to the divergence of the active force

$\displaystyle S_{\text{SB}}$
$\displaystyle\sim\vec{\nabla}\cdot\vec{f}_{\zeta}.$ (15)

####  8.3.6 Defects and viscometric surfaces

To find the probability distribution of nearest distances between defects and
viscometric surfaces ($\mathcal{Q}=0$ contours), we calculate the distances
for all defects to each point on the boundary and included only the smallest
distance (to avoid introducing cut-off distances). In the case of mirror-
symmetric configurations (Fig. 2c), only the nearest boundary is counted.

To calculate the distribution of alignment angles shown in Fig. 2g, we first
group $+1/2$ defects to $\mathcal{Q}=0$ contours. Each defect is associated to
a boundary if the defect resides within 2 simulation or experimental sites
away from any point of the boundary (experimental sites use the resolution
details explained in § 8.1.4). The closest point $j$ between the defect and
the list of points on the $\mathcal{Q}=0$ boundary is found and is the
reference point for tangent vector to be calculated from. The tangent vector
is then calculated using a forward derivative with the next point in the
boundary list $j+1$ (with the forward direction matching the $+1/2$ defect’s
orientation) as

$\displaystyle\vec{t}$
$\displaystyle=\frac{\vec{r}(j+1)-\vec{r}(j)}{\left|\vec{r}(j+1)-\vec{r}(j)\right|},$
(16)

where $\vec{r}(j)$ is the position vector of the chosen point $j$ on the
$\mathcal{Q}=0$ boundary. We chose the forward derivative over a centered
derivative because the flow structure ahead of the defect is more significant
to the future dynamics of the defect trajectory. The $+1/2$ defect orientation
vector is found by Eq. (8). Each alignment angle is therefore found from the
angle between these two vectors and presented in the angle interval of
$0$\mathrm{\SIUnitSymbolDegree}$$ and $90$\mathrm{\SIUnitSymbolDegree}$$. If
the defect is associated to the mirror-symmetric $\mathcal{Q}=0$ regime, then
the same $+1/2$ defect orientation vector is included for the alignment angle
calculation with both boundaries.

To investigate the mirror-symmetry breaking, we sort the configurations based
on the three classes of defect/flow configurations (Fig. 2c-e). For each
$+1/2$ defect, we calculate $\left\langle\mathcal{Q}\right\rangle_{\Delta
r^{2}}$ as the average $\mathcal{Q}$-criterion in a small region $\Delta
r^{2}$ set to be a rectangle of seven by six lattice units on either side of
the mirror-symmetry axis of the defect director field (i.e. the defect tail).
From $\left\langle\mathcal{Q}\right\rangle_{\Delta r^{2}}$ three different
defect populations can be identified: one group where
$\left\langle\mathcal{Q}\right\rangle_{\Delta r^{2}}$ is close to zero (with
mirror symmetry), and two where $\left\langle\mathcal{Q}\right\rangle_{\Delta
r^{2}}$ has a large positive or negative value (the mirror-symmetry broken
cases). To ensemble average over these three regimes, defects are binned into
lowest $10\%$ of $\left\langle\mathcal{Q}\right\rangle_{\Delta r^{2}}$, middle
$80\%$ and highest $10\%$. This gives Ext. Data Fig. 2b, Ext. Data Fig. 2a and
Ext. Data Fig. 2c respectively. The results are found to be robust to varying
the exact binned regions $\Delta r^{2}$. Ext. Data Fig. 2a shows the ensemble-
averaged case of mirror-symmetric defects at the $\mathcal{Q}=0$ intersection,
with opposite-handed vortices on either side. Ext. Data Fig. 2b shows the
broken-mirror-symmetry state for defects with clockwise-handed $\mathcal{Q}>0$
vortex to the right of the defect head and $\mathcal{Q}<0$ to the left and
vice-versa in Ext. Data Fig. 2c. In both cases, the defect aligns parallel
with the $\mathcal{Q}=0$ tangent. In addition, the local velocity field in
Ext. Data Fig. 2(b-c) on either side of the defect has the same curvature
direction as the vortex handedness, while the mirror-symmetric regime has an
inversion to the focal point.

###  8.4 Modeling

####  8.4.1 Stokesian solitary defect model

The $\mathcal{Q}$-criterion and flow field around $+1/2$ and $-1/2$
topological defects is shown in Fig. 2b and f. The $\mathcal{Q}$-criterion
field is determined from the velocity field solutions $\vec{u}_{\pm}$ for a
single solitary $\pm 1/2$ defects first derived by Giomi, et al. [28]. To find
the ideal flow field structure around solitary defects of topological charge
$k$, they fixed the director field
$n_{k}=\cos\left(k\phi\right)\hat{x}+\sin\left(k\phi\right)\hat{y}$ and so
Beris-Edwards equation (Eq. (2)) does not evolve. They then solve the
generalized Navier-Stokes equations (Eq. (3)) for the steady-state
($D_{t}\vec{u}=\vec{0}$) in the absence of friction and with negligible
elastic stresses, leaving the Stokes equation for the velocity field,

$\displaystyle\eta\nabla^{2}\vec{u}-\vec{\nabla}p+\vec{f}_{\zeta}$
$\displaystyle=0$ (17)

as a response to the active force density $\vec{f}_{\zeta}$ (Eq. (14)). For an
extensile plus half defect ($k=+1/2$), the active force
$\vec{f}_{\zeta}=-\zeta\hat{x}/2r$ is in the direction of the defect
$-\hat{x}$ (bend-side) and decays inversely with distance $r$ from the defect
core [28]. This rapid decay of the force is the basis for the simple view of
defects as point sources of forcing. While a bulk active nematic cannot exert
a global net force, the conceit of this ideal Stokesian model is that a single
solitary defect exists breaking topological charge neutrality, allowing a net
force. The solution to Eq. (17) can be written as

$\displaystyle\vec{u}\left(\vec{r}\right)$ $\displaystyle=\int
dA^{\prime}\bm{\bm{G}}\left(\vec{r},\vec{r}^{\prime}\right)\cdot\vec{f}_{\zeta}\left(\vec{r}^{\prime}\right)$
(18)

with the two-dimensional Oseen tensor [28]

$\displaystyle\bm{\bm{G}}\left(\vec{r},\vec{r}^{\prime}\right)$
$\displaystyle=\frac{1}{4\pi\eta}\left[\left(\log\frac{L}{\left|\vec{r}-\vec{r}^{\prime}\right|}-1\right)\bm{\bm{\delta}}+\frac{(\vec{r}-\vec{r}^{\prime})\otimes(\vec{r}-\vec{r}^{\prime})}{\left|\vec{r}-\vec{r}^{\prime}\right|^{2}}\right]$
(19)

in a system bound by the length scale $L$. Substituting in the active force
for the fixed, solitary defects gives the velocity fields

$\displaystyle\vec{u}_{+}(r,\phi)$
$\displaystyle=-\frac{\zeta}{12\eta}\left[\left\\{3(R-r)+r\cos
2\phi\right\\}\hat{x}+r\sin 2\phi\hat{y}\right]$ (20)
$\displaystyle\vec{u}_{-}(r,\phi)$ $\displaystyle=-\frac{\zeta r}{12\eta
R}\left[\left\\{\frac{3}{4}(r-R)\cos 2\phi-\frac{R}{5}\cos
4\phi\right\\}\hat{x}-\left\\{(\frac{3}{4}r-R)\sin 2\phi+\frac{R}{5}\sin
4\phi\right\\}\hat{y}\right],$ (21)

where $\zeta$ is the activity, $\eta$ is the viscosity and $R$ is an
integration length scale due to the logarithmic nature of hydrodynamic
interactions in 2D (§ 8.4.4). These are shown in Fig. 2b and f. Taking
derivatives of Eqs. (20)-(21) produces the velocity gradient tensor
$\bm{\bm{L}}=\vec{\nabla}\otimes\vec{u}$ and the $\mathcal{Q}$-criterion field
is calculated as the second invariant of $\bm{\bm{L}}$. These calculations and
corresponding plots were performed using Mathematica [64]. For the $+1/2$
defect case, the $\mathcal{Q}$-criterion takes the form,

$\displaystyle\mathcal{Q}_{+}=\left(\frac{\eta\zeta}{6}\right)^{2}\left(3\sin^{2}\phi-1\right).$
(22)

In this construction, the $+1/2$ defect is aligned with a self-propulsion
direction in the (negative) $\hat{x}$ direction, with an angle
$\theta=\sin^{-1}\left(3^{-1/2}\right)\approx
35.26$\mathrm{\SIUnitSymbolDegree}$$ to the $\mathcal{Q_{+}}=0$ line. This
angle is compared in Fig. 2g against the observed distribution of alignment
angles in simulations and experiments, as a benchmark for an isolated $+1/2$
defect. The solution for this angle is predicted independent of the activity,
viscosity and overall system size.

####  8.4.2 Stokesian straight bend-wall model

Since our experimental and numerical results indicate that bend walls are
coincident with $\mathcal{Q}=0$ and crucial to the spontaneous self-
constraint, we explore a series of models for the active flows generated by
bend walls to identify which features of the bend walls are vital to this
interdependence. Bend walls are known to arise from the hydrodynamic bend
instability [65] leading to arch-like bends of the director field, reminiscent
of Néel walls [30]. First consider bend walls modelled as infinitely long,
perfectly straight bands of alternating bend and splay given by
$\vec{n}=\cos\theta\hat{x}+\sin\theta\hat{y}$, where the angle $\theta$ is

$\displaystyle\theta(y;\lambda)$
$\displaystyle=\pi\cos^{2}\left(\frac{y}{\lambda}\right),$ (23)

for a modulation length scale $\lambda$. This generates the force density
field
$\vec{f}_{\zeta}=\zeta\frac{\pi}{\lambda}\sin\frac{2y}{\lambda}\left(\cos
2\theta\hat{x}+\sin 2\theta\hat{y}\right)$. Since the system has translational
symmetry in $x$, only derivatives with respect to $y$ exist for the pressure
and velocity. Through the incompressibility equation, $\nabla\cdot\vec{u}=0$,
the absense of $x-$gradients requires $\frac{\partial\vec{u}_{y}}{\partial
y}=0$. This leads to only one non-zero contribution to the velocity gradient
tensor, $\frac{\partial\vec{u}_{x}}{\partial y}$, which is obtained through
integrating the $x-$momentum equation in Eq. (17). The velocity gradient
tensor is found as,

$\displaystyle\bm{\bm{L}}$
$\displaystyle=-\zeta\frac{\sin\left(\pi\cos\left[2y/\lambda\right]\right)}{2\eta}\begin{pmatrix}0&1\\\
0&0\end{pmatrix}.$ (24)

Since this is upper diagonal at all points (representing simple horizontal
shear), it is a viscometric flow and $\mathcal{Q}=0$ everywhere. Thus, a
straight, infinitely long series of kink walls is insufficient to explain the
colocation of bend walls and $\mathcal{Q}=0$.

####  8.4.3 Wavy bend wall model

Next we consider if bend wall curvature explains their coincidence with
$\mathcal{Q}=0$. The bend walls from § 8.4.2 are perturbed into infinitely
long sinusoidal waves with director field once again of the form
$\vec{n}=\cos\theta\hat{x}+\sin\theta\hat{y}$, where the angle $\theta$ but
with angle

$\displaystyle\theta\left(x,y;\lambda\right)$
$\displaystyle=\pi\cos^{2}\left(\frac{\pi x}{4\lambda}\right)+\frac{\pi
y}{\lambda}.$ (25)

This director field generates a periodic sequence of bend and splay walls that
modulate sinusoidally about slices of $\hat{y}$ (Ext. Data Fig. 3a). This
initial director field is held fixed and the velocity field is numerically
evolved to steady state according to Eq. (3). This simple model of bend walls
is sufficient to generate vortices enclosed by viscometric contours (Ext. Data
Fig. 3b). The vortices reside on the inside regions of the bend and splay
walls with greatest curvature, with alternating handedness generated through
the coupling between the generated active force and the walls curvature
orientation. Thus, bend wall curvature is essential to the existence of
viscometric surfaces in active nematics.

To further explore the nonlinear coupling between the modulated director field
and the active flow field, we then hold the resulting initial steady-state
velocity field fixed, but allow the director field to evolve according to Eq.
(2). The bend walls proceed to constrict towards a thin kink line, while the
strong divergence in the splay wall breaks open, leaving minimised gradients
in splay around the bend walls (Ext. Data Fig. 3c). This showcases the active
instability that drives growth in bend-type deformations — even with an
iterative progression of the two fields, the bend walls drive flows and the
flows in turn concentrates and enhances bend [65]. Since the centerline of the
bend walls coincides with $\mathcal{Q}=0$, the backbones themselves are
largely undeformed, while the nearby nematic field is rotated towards a
parallel alignment with the bend wall centerline, localising the gradients in
bend. Iterating again, demonstrates that the vortices begin to merge and
elongate. The centerline of the bend wall still resides on $\mathcal{Q}=0$

Overall, this model demonstrates that modulation in bend walls is able to
split distinct $\mathcal{Q}$ regions without needing defects. The geometry of
the flow generated by the bend has $\mathcal{Q}=0$ on the center of the bend
wall, which gives bend wall a protection against deforming — the active flows
can constrict but not cleave the bend walls. The shape of the bend wall can
only be perturbed by the creation of defects.

####  8.4.4 Stokesian line force model of bend walls

The hydrodynamic instability constricts infinitely long bend into narrow kink
walls and their curvature facilitates closed $\mathcal{Q}=0$ contours (§
8.4.3). However, after the initial onset of active turbulence, finitely long
bend walls continue to arise via constriction and be unzipped by $+1/2$
defects (Supplementary Video 4) [32, 15]. The active force density is non-
negligible along the bend walls and so we propose a simplified model of finite
narrow kink walls being unzipped by $+1/2$ defects modelled as an idealized
line of active force. As before (§ 8.4.1), elastic stresses are neglected and
the Stokes equation (Eq. (17)) is solved for the velocity field and
$\mathcal{Q}=0$ contours. This is done for four idealized bend wall
conformations: (i) a point, (ii) a straight finite line, (iii) a circle and an
arc. For all extended bend wall conformations, the Stokes equation (Eq. (17))
is numerically solved. For a straight finite line a $2000\times 1000$ grid is
used and $1000\times 1000$ for the circle and arc. The force lines are
constructed by arranging between $500$ to $1000$ point sources in a linear
formation, with constant forcing magnitude orientated parallel to the line
tangent. Each point on the discretised line contributes to the overall
velocity field through the Oseen tensor (Eq. (19)).

#####  (i) Point force

In the far-field limit away from a solitary bend wall, it is effectively a
point force. From the Oseen tensor (Eq. (19)), crudely modelling a bend walls
as a point force $f_{\zeta}\hat{x}$ at the origin produces a 2D Stokeslet

$\displaystyle\vec{u}\left(\vec{r}\right)$
$\displaystyle=\frac{f_{\zeta}}{4\pi\eta}\left[\left(\log\frac{L}{\left|\vec{r}\right|}-1\right)\hat{x}+\frac{x\vec{r}}{\left|\vec{r}\right|^{2}}\right].$
(26)

Taking the gradient results in the $\mathcal{Q}$-criterion

$\displaystyle\mathcal{Q}$
$\displaystyle=-\left(\frac{f_{\zeta}}{4\pi\eta\left|\vec{r}\right|}\right)^{2}\cos{2\theta},$
(27)

which predicts an intersection of viscometric surfaces at an angle
$\theta=45^{\circ}$.

#####  (ii) Straight finite line

In the vicinity of a finite straight-kink wall, the flow field geometry
differs from the point-source-like far-field approximation. A horizontal line
of constant active force density is constructed using 1000 point sources
arranged across 1000 horizontal lattice points, oriented towards the positive
$\hat{x}$ direction. The velocity field and zero isolines of the
$\mathcal{Q}$-criterion (Eq. (1)) are shown in Fig. 3c. The
$\mathcal{Q}$-criterion retains the mirror symmetry observed for the solitary
defect or far-field point source, but now differs by constraining the
alignment angle $\theta$ towards more acute values. Smaller defect alignment
angles in the proximity of bend walls is more consistent with the peak in Fig.
2g, compared with the angle predicted by the solitary defect model (§ 8.4.1).

#####  (iii) Circle and arc

To provide intuition of the flow geometry generated by bend walls with
curvature, we next consider a circular arc force centered on the origin. The
arc is parameterised in terms of $s$ that can wrap out an angle in a
prescribed interval with endpoints lying between $0$ and $2\pi$. Positions on
the arc are $\vec{r}^{\prime}(s)=R\left(\cos s\hat{x}+\sin s\hat{y}\right)$,
where $R$ is the radius. As before, the force contour is assumed to have
constant magnitude $f_{\zeta}$ and orient azimuthally. By choosing the forces
to follow the arc in a clock-wise sense, the force is
$\vec{f}_{\zeta}^{\prime}(s)=f_{\zeta}\left(\sin s\hat{x}-\cos
s\hat{y}\right)$. For the arc, numerical integration is performed using 500
point sources distributed between $s=-\pi/2$ and $s=\pi/2$, while the circle
uses 1000 point sources between $s=0$ and $s=2\pi$. The results for the
velocity field and $\mathcal{Q}$-criterion are shown in Fig. 3d (circle) and e
(arc). For both cases, the $\mathcal{Q}$-criterion goes to zero when
$\left|\vec{r}\right|=R$. This result (Fig. 3e) shows that this simplified
model can capture the correlation between bend walls and $\mathcal{Q}=0$.
Through comparison with the straight finite line case, curvature facilitates
the transformation from two vortices to one at the source of the active
forcing. Therefore, lines of active force generated by narrow, bend walls, can
capture the essential geometry of the flow field if they are (i) finite and
(ii) curved.

##  9 Extended data figures

Extended Data Figure 1:  Probability distribution functions (PDF) of vortex
properties, shown for experiments and simulations. (a) Shape anisotropy
parameter $\kappa^{2}$ (Eq. (11)). The limiting value $\kappa=0$ represents a
circle and $\kappa=1$ indicates a line. (b) Vortex area non-dimensionalized by
the active length scale $\ell_{\zeta}$. (c) Total circulation within a
$\mathcal{Q}=0$ closed contour (§ 8.3.4). Extended Data Figure 2:  Time and
ensemble averaged director fields and $\mathcal{Q}$-criterion around $+1/2$
defects sorted by handedness. (a) Mirror-symmetric state showing
$\mathcal{Q}$-criterion structure around $+1/2$ defect. The $\mathcal{Q}=0$
contours join towards a point, separating two vortices of opposite handedness.
Strain rate dominated flow regions are located in front and behind the defect.
(b) Broken-symmetry state with a defect following a clockwise vortex. (c) Same
as (b) but following an anti-clockwise vortex. Director field shown as black
lines, contours where $\mathcal{Q}=0$ are shown as solid lines and colormap
corresponds to vorticity field.  Extended Data Figure 3:  Infinitely long,
repeating wavy bend wall model (§ 8.4.3). (a) Initial splay-bend parameter
($S_{\text{SB}}$) field with repeating model bend walls in purple, with same
colorbar as Fig. 3a. (b) Resulting $\mathcal{Q}$ of steady-state active flows
due to fixed $S_{\text{SB}}$ from (a). (c) Resulting steady-state
$S_{\text{SB}}$ due to fixed flows from (b). (d) Resulting steady-state flows
due to fixed $S_{\text{SB}}$ from (c). Director field is shown as black lines,
the velocity field is shown as grey arrows and contours where $\mathcal{Q}=0$
are shown as solid lines.  Extended Data Figure 4:  Numerical simulation of
defect dynamics and circulation injection in two-dimensional dance defect
state. (a) Defect dynamics and circulation time evolution. The pink shaded
region identifies the time periods where at least one defect resides on two
vortices of alternate handedness simultaneously. (top) Average absolute +1/2
defect displacement, $y$, from the channel centre, scaled by the channel
height $L_{y}$. Shaded region represents the standard deviation of the defect
height. (middle) Ensemble average of the circulation, $\Gamma$, contained
inside the central clockwise (red) and anti-clockwise (blue) vortices. Shaded
region represents the standard deviation for each of the clockwise (red
shading) and anti-clockwise (blue shading) vortices. (bottom) Total number of
+1/2 defects, $N_{+1/2}$, associated with the central clockwise (red) and
anti-clockwise (blue) vortices. (b) Reconnecting bend-wall structure supports
the reshaping of viscometric $\mathcal{Q}=0$ surfaces. Colormap identifies the
line-like nematic bend deformations through the strongly-negative splay-bend
parameter $S_{SB}$ (§ 8.3.5). Zero-isolines of the $\mathcal{Q}$-criterion are
shown as solid lines, colored in red for clockwise and blue for anticlockwise.
Nematic $+1/2$ defects are marked by green comet-shaped symbols and minus-half
defects by dark blue-trefoil shaped symbols (§ 8.3.1). (c) Same as (b) for a
later frame in Supplementary Video 10 when $+1/2$ defects are in their mirror-
symetry broken conformation.

##  10 Supplementary Videos

  1. Supplementary Video 1

Fluorescence microscopy video of a bundled microtubule/kinesin network at the
oil-water interface exhibiting active turbulence. Solid lines indicate the
viscometric surfaces ($\mathcal{Q}=0$). Topological $\pm 1/2$ defects marked
by green comet-like and blue trefoil-like symbols. Figure 1a from this video,
and scale is the same.

  2. Supplementary Video 2

Numerically simulated, two-dimensional active turbulence for extensile and
flow-aligning nematic (§ 8.2). Colormap shows the circulation inside
vorticity-dominated regions ($\mathcal{Q}>0$). Strain-rate-dominated regions
($\mathcal{Q}<0$) uncolored. Solid lines indicate the viscometric surfaces
($\mathcal{Q}=0$). Topological $\pm 1/2$ defects marked by green comet-like
and blue trefoil-like symbols. Figure 1b from this video.

  3. Supplementary Video 3

Same as Supplementary Video 2 with colormap showing the vorticity field.
Viscometric surfaces and $\pm 1/2$ defects illustrated as in Supplementary
Video 2. Lines of zero-vorticity contours are shown as orange-dash-dot lines.
Figure 1f from this video.

  4. Supplementary Video 4

Same as Supplementary Video 2 with colormap showing the splay-bend parameter,
$S_{\text{SB}}$, field. Viscometric surfaces and $\pm 1/2$ defects illustrated
as in Supplementary Video 2. Figure 3a from this video.

  5. Supplementary Video 5

Numerically simulated, two-dimensional active turbulence in the flow tumbling
regime (§ 8.2). Circulation, viscometric surfaces and $\pm 1/2$ defects
illustrated as in Supplementary Video 2. Figure 4a from this video.

  6. Supplementary Video 6

Numerically simulated, two-dimensional active turbulence with contractile
activity (§ 8.2). Splay-bend parameter, viscometric surfaces and $\pm 1/2$
defects illustrated as in Supplementary Video 7. Figure 4b from this video.

  7. Supplementary Video 7

Numerically simulated, two-dimensional active turbulence with anisotropic
friction (§ 8.2). Circulation, viscometric surfaces and $\pm 1/2$ defects
illustrated as in Supplementary Video 2. Figure 4c from this video.

  8. Supplementary Video 8

Numerically simulated, two-dimensional dancing defects and vortex lattice in
channel with impermeable, no-slip, homogeneous-anchored walls (§ 8.2).
Circulation, viscometric surfaces and $\pm 1/2$ defects illustrated as in
Supplementary Video 2. Figure 4d-e from this video.

  9. Supplementary Video 9

Experimental microscopy of bundled microtubule/kinesin network confined
between micro-printed grid walls (§ 8.1.3). Fluorescence microscopy,
viscometric surfaces and $\pm 1/2$ defects illustrated as in Supplementary
Video 1. Figure 4f from this video, and scale is the same.

  10. Supplementary Video 10

Same as Supplementary Video 8 with colormap showing the splay-bend parameter,
$S_{\text{SB}}$, field. Splay-bend parameter, viscometric surfaces and $\pm
1/2$ defects illustrated as in Supplementary Video 7. Ext. Data Fig. 4b-c from
this video.

  11. Supplementary Video 11

Numerically simulated, three-dimensional dancing defects and vortext lattice
in square duct with impermeable, no-slip, homogeneous-anchored walls (§ 8.2).
Grey shading demarcates vorticity-dominate ($\mathcal{Q}>0$) and strain-rate
dominate ($\mathcal{Q}<0$) space. Disclinations lines visualized as thin tubes
colored $\cos\beta$ (§ 8.3.1), where yellow indicates $\cos\beta\approx+1$
representing local $+1/2$ wedge profiles, green $\cos\beta\approx 0$ for
twist-type profiles and blue $\cos\beta\approx-1$ for $-1/2$ wedge profiles.
Figure 4g from this video.

  12. Supplementary Video 12

Same as Supplementary Video 2 with colormap showing the
$\mathcal{Q}$-criterion. Viscometric surfaces and $\pm 1/2$ defects
illustrated as in Supplementary Video 2.

## References

  * [1] Lukasz Klotz, Grégoire Lemoult, Kerstin Avila, and Björn Hof.  Phase transition to turbulence in spatially extended shear flows.  Physical Review Letters, 128:014502, 2022.
  * [2] Karthik Duraisamy, Gianluca Iaccarino, and Heng Xiao.  Turbulence modeling in the age of data.  Annual Review of Fluid Mechanics, 51(1):357–377, 2019.
  * [3] Luca Brandt and Filippo Coletti.  Particle-laden turbulence: Progress and perspectives.  Annual Review of Fluid Mechanics, 54(1):159–189, 2022.
  * [4] Chuanfei Dong, Liang Wang, Yi-Min Huang, Luca Comisso, Timothy A. Sandstrom, and Amitava Bhattacharjee.  Reconnection-driven energy cascade in magnetohydrodynamic turbulence.  Science Advances, 8(49):eabn7627, 2022.
  * [5] Ladislav Skrbek, David Schmoranzer, Šimon Midlik, and Katepalli R. Sreenivasan.  Phenomenology of quantum turbulence in superfluid helium.  Proceedings of the National Academy of Sciences, 118(16):e2018406118, 2021.
  * [6] Kazumasa A. Takeuchi, Masafumi Kuroda, Hugues Chaté, and Masaki Sano.  Directed percolation criticality in turbulent liquid crystals.  Physical Review Letters, 99:234503, 2007.
  * [7] David N Hosking and Alexander A Schekochihin.  Reconnection-controlled decay of magnetohydrodynamic turbulence and the role of invariants.  Physical Review X, 11:041005, 2021.
  * [8] Juan Ignacio Polanco, Nicolás P Müller, and Giorgio Krstulovic.  Vortex clustering, polarisation and circulation intermittency in classical and quantum turbulence.  Nature Communications, 12(1):7090, 2021.
  * [9] Maruša Mur, Žiga Kos, Miha Ravnik, and Igor Muševič.  Continuous generation of topological defects in a passively driven nematic liquid crystal.  Nature Communications, 13(1):6855, 2022.
  * [10] Lakshmi Balasubramaniam, René-Marc Mège, and Benoît Ladoux.  Active nematics across scales from cytoskeleton organization to tissue morphogenesis.  Current Opinion in Genetics & Development, 73:101897, 2022.
  * [11] Suraj Shankar, Anton Souslov, Mark J Bowick, M Cristina Marchetti, and Vincenzo Vitelli.  Topological active matter.  Nature Reviews Physics, 4(6):380–398, 2022.
  * [12] Henricus H Wensink, Jörn Dunkel, Sebastian Heidenreich, Knut Drescher, Raymond E Goldstein, Hartmut Löwen, and Julia M Yeomans.  Meso-scale turbulence in living fluids.  Proceedings of the National Academy of Sciences, 109(36):14308–14313, 2012.
  * [13] Tim Sanchez, Daniel TN Chen, Stephen J DeCamp, Michael Heymann, and Zvonimir Dogic.  Spontaneous motion in hierarchically assembled active matter.  Nature, 491(7424):431–434, 2012.
  * [14] Guillaume Duclos, Raymond Adkins, Debarghya Banerjee, Matthew S. E. Peterson, Minu Varghese, Itamar Kolvin, Arvind Baskaran, Robert A. Pelcovits, Thomas R. Powers, Aparna Baskaran, Federico Toschi, Michael F. Hagan, Sebastian J. Streichan, Vincenzo Vitelli, Daniel A. Beller, and Zvonimir Dogic.  Topological structure and dynamics of three-dimensional active nematics.  Science, 367(6482):1120–1124, 2020.
  * [15] Suraj Shankar and M Cristina Marchetti.  Hydrodynamics of active defects: From order to chaos to defect ordering.  Physical Review X, 9(4):041047, 2019.
  * [16] Luiza Angheluta, Zhitao Chen, M Cristina Marchetti, and Mark J Bowick.  The role of fluid flow in the dynamics of active nematic defects.  New Journal of Physics, 23(3):033009, 2021.
  * [17] Jonas Rønning, Cristina M. Marchetti, Mark J. Bowick, and Luiza Angheluta.  Flow around topological defects in active nematic films.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 478(2257):20210879, 2022.
  * [18] Vasil Bratanov, Frank Jenko, and Erwin Frey.  New class of turbulence in active fluids.  Proceedings of the National Academy of Sciences, 112(49):15048–15053, 2015.
  * [19] Ricard Alert, Jean-François Joanny, and Jaume Casademunt.  Universal scaling of active nematic turbulence.  Nature Physics, 16(6):682–688, 2020.
  * [20] Siddhartha Mukherjee, Rahul K Singh, Martin James, and Samriddhi Sankar Ray.  Intermittency, fluctuations and maximal chaos in an emergent universal state of active turbulence.  Nature Physics, pages 1–7, 2023.
  * [21] Ricard Alert, Jaume Casademunt, and Jean-François Joanny.  Active turbulence.  Annual Review of Condensed Matter Physics, 13(1):143–170, 2022\.
  * [22] Jerôme Hardoüin, Claire Doré, Justine Laurent, Teresa Lopez-Leon, Jordi Ignés-Mullol, and Francesc Sagués.  Active boundary layers in confined active nematics.  Nature Communications, 13(1):6675, 2022.
  * [23] Perry W Ellis, Daniel JG Pearce, Ya-Wen Chang, Guillermo Goldsztein, Luca Giomi, and Alberto Fernandez-Nieves.  Curvature-induced defect unbinding and dynamics in active nematic toroids.  Nature Physics, 14(1):85–90, 2018.
  * [24] C. Rorai, F. Toschi, and I. Pagonabarraga.  Coexistence of active and hydrodynamic turbulence in two-dimensional active nematics.  Physical Review Letters, 129:218001, 2022.
  * [25] Luca Giomi.  Geometry and topology of turbulence in active nematics.  Physical Review X, 5:031003, 2015.
  * [26] Carles Blanch-Mercader, Victor Yashunsky, Simon Garcia, Guillaume Duclos, Luca Giomi, and Pascal Silberzan.  Turbulent dynamics of epithelial cell cultures.  Physical Review Letters, 120(20):208101, 2018.
  * [27] Tyler N Shendruk, Amin Doostmohammadi, Kristian Thijssen, and Julia M Yeomans.  Dancing disclinations in confined active nematics.  Soft Matter, 13(21):3853–3862, 2017.
  * [28] Luca Giomi, Mark J Bowick, Prashant Mishra, Rastko Sknepnek, and M Cristina Marchetti.  Defect dynamics in active nematics.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 372(2029):20130365, 2014.
  * [29] Berta Martínez-Prat, Jordi Ignés-Mullol, Jaume Casademunt, and Francesc Sagués.  Selection mechanism at the onset of active turbulence.  Nature Physics, 15(4):362–366, 2019.
  * [30] Aurelio Patelli, Ilyas Djafer-Cherif, Igor S. Aranson, Eric Bertin, and Hugues Chaté.  Understanding dense active nematics from microscopic models.  Physical Review Letters, 123:258001, 2019.
  * [31] R. Aditi Simha and Sriram Ramaswamy.  Hydrodynamic fluctuations and instabilities in ordered suspensions of self-propelled particles.  Physical Review Letters, 89:058101, 2002.
  * [32] Sumesh P Thampi, Ramin Golestanian, and Julia M Yeomans.  Instabilities and topological defects in active nematics.  Europhysics Letters, 105(1):18001, 2014.
  * [33] Kristian Thijssen and Amin Doostmohammadi.  Binding self-propelled topological defects in active turbulence.  Physical Review R, 2(4):042008, 2020.
  * [34] Kristian Thijssen, Dimitrius A Khaladj, S Ali Aghvami, Mohamed Amine Gharbi, Seth Fraden, Julia M Yeomans, Linda S Hirst, and Tyler N Shendruk.  Submersed micropatterned structures control active nematic flow, topology, and concentration.  Proceedings of the National Academy of Sciences, 118(38):e2106038118, 2021.
  * [35] Berta Martínez-Prat, Ricard Alert, Fanlong Meng, Jordi Ignés-Mullol, Jean-Fran çois Joanny, Jaume Casademunt, Ramin Golestanian, and Francesc Sagués.  Scaling regimes of active turbulence with external dissipation.  Physical Review X, 11:031065, 2021.
  * [36] Pau Guillamat, Jordi Ignés-Mullol, and Francesc Sagués.  Taming active turbulence with patterned soft interfaces.  Nature Communications, 8(1):564, 2017.
  * [37] Kristian Thijssen, Luuk Metselaar, Julia M Yeomans, and Amin Doostmohammadi.  Active nematics with anisotropic friction: The decisive role of the flow aligning parameter.  Soft Matter, 16(8):2065–2074, 2020.
  * [38] Jérôme Hardoüin, Rian Hughes, Amin Doostmohammadi, Justine Laurent, Teresa Lopez-Leon, Julia M Yeomans, Jordi Ignés-Mullol, and Francesc Sagués.  Reconfigurable flows and defect landscape of confined active nematics.  Communications Physics, 2(1):121, 2019.
  * [39] Ryan R Keogh, Santhan Chandragiri, Benjamin Loewe, Tapio Ala-Nissila, Sumesh P Thampi, and Tyler N Shendruk.  Helical flow states in active nematics.  Physical Review E, 106(1):L012602, 2022.
  * [40] Nika Kralj, Miha Ravnik, and Žiga Kos.  Defect line coarsening and refinement in active nematics.  Physical Review Letters, 130:128101, 2023.
  * [41] Mattia Serra, Linnea Lemma, Luca Giomi, Zvonimir Dogic, and Lakshminarayanan Mahadevan.  Defect-mediated dynamics of coherent structures in active nematics.  Nature Physics, pages 1–7, 2023.
  * [42] Rui Zhang, Ali Mozaffari, and Juan J de Pablo.  Autonomous materials systems from active liquid crystals.  Nature Reviews Materials, 6(5):437–453, 2021.
  * [43] Rui Zhang, Ali Mozaffari, and Juan J. de Pablo.  Logic operations with active topological defects.  Science Advances, 8(8):eabg9060, 2022.
  * [44] De-Qing Zhang, Peng-Cheng Chen, Zhong-Yi Li, Rui Zhang, and Bo Li.  Topological defect-mediated morphodynamics of active–active interfaces.  Proceedings of the National Academy of Sciences, 119(50):e2122494119, 2022.
  * [45] Jayabrata Dhar, Anh LP Thai, Arkajyoti Ghoshal, Luca Giomi, and Anupam Sengupta.  Self-regulation of phenotypic noise synchronizes emergent organization and active transport in confluent microbial environments.  Nature Physics, 18(8):945–951, 2022.
  * [46] Daniel Needleman and Zvonimir Dogic.  Active matter at the interface between materials science and cell biology.  Nature Reviews Materials, 2(9):1–14, 2017.
  * [47] Wang Xi, Thuan Beng Saw, Delphine Delacour, Chwee Teck Lim, and Benoit Ladoux.  Material approaches to active tissue mechanics.  Nature Reviews Materials, 4(1):23–44, 2019.
  * [48] Silke Henkes, Kaja Kostanjevec, J Martin Collinson, Rastko Sknepnek, and Eric Bertin.  Dense active matter model of motion patterns in confluent cell monolayers.  Nature Communications, 11(1):1405, 2020.
  * [49] Farzan Vafa and L. Mahadevan.  Active nematic defects and epithelial morphogenesis.  Physical Review Letters, 129:098102, 2022.
  * [50] Deqing Sun, Stefan Roth, and Michael J. Black.  Secrets of optical flow estimation and their principles.  In 2010 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, pages 2432–2439, 2010.
  * [51] Enzo Orlandini, Michael E Cates, Davide Marenduzzo, Luca Tubiana, and Julia M Yeomans.  Hydrodynamic of active liquid crystals: A hybrid lattice boltzmann approach.  Molecular Crystals and Liquid Crystals, 494(1):293–308, 2008.
  * [52] Sumesh P Thampi and Julia M Yeomans.  Active turbulence in active nematics.  The European Physical Journal Special Topics, 225:651–662, 2016\.
  * [53] Cody D. Schimming and Jorge Viñals.  Singularity identification for the characterization of topology, geometry, and motion of nematic disclination lines.  Soft Matter, 18:2234–2244, 2022.
  * [54] Matthew L. Blow, Sumesh P. Thampi, and Julia M. Yeomans.  Biphasic, lyotropic, active nematics.  Physical Review Letters, 113:248303, 2014.
  * [55] Dario Dell’Arciprete, Matthew L Blow, Aidan T Brown, Fred DC Farrell, Juho S Lintuvuori, Alexander F McVey, Davide Marenduzzo, and Wilson CK Poon.  A growing bacterial colony in two dimensions as an active nematic.  Nature Communications, 9(1):4190, 2018.
  * [56] Arthur J Vromans and Luca Giomi.  Orientational properties of nematic disclinations.  Soft Matter, 12(30):6490–6495, 2016.
  * [57] Prabhu Ramachandran and Gaël Varoquaux.  Mayavi: 3D visualization of scientific data.  Computing in Science & Engineering, 13(2):40–51, 2011.
  * [58] Jinhee Jeong and Fazle Hussain.  On the identification of a vortex.  Journal of Fluid Mechanics, 285:69–94, 1995.
  * [59] Stefan Van der Walt, Johannes L Schönberger, Juan Nunez-Iglesias, François Boulogne, Joshua D Warner, Neil Yager, Emmanuelle Gouillart, and Tony Yu.  scikit-image: Image processing in python.  PeerJ, 2:e453, 2014.
  * [60] Raja R. Huilgol.  A class of motions with constant stretch history.  Quarterly of Applied Mathematics, 29(1):1–15, 1971.
  * [61] Roney L. Thompson and Paulo R. Souza Mendes.  Persistence of straining and flow classification.  International Journal of Engineering Science, 43(1):79–105, 2005\.
  * [62] Doros N. Theodorou and Ulrich W. Suter.  Shape of unperturbed linear polymers: Polypropylene.  Macromolecules, 18(6):1206–1214, 1985.
  * [63] Simon Čopar, Tine Porenta, and Slobodan Žumer.  Visualisation methods for complex nematic fields.  Liquid Crystals, 40(12):1759–1768, 2013.
  * [64] Wolfram Research, Inc.  Mathematica, Version 13.2.  Champaign, IL, 2022.
  * [65] Sriram Ramaswamy and Madan Rao.  Active-filament hydrodynamics: Instabilities, boundary conditions and rheology.  New Journal of Physics, 9(11):423, 2007.