Spontaneous self-constraint in active nematic flows

Active processes drive biological dynamics across various scales and include subcellular cytoskeletal remodelling, tissue development in embryogenesis and the population-level expansion of bacterial colonies. In each of these, biological functionality requires collective flows to occur while self-organised structures are protected. However, the mechanisms by which active flows can spontaneously constrain their dynamics to preserve structure are not known. Here, by studying collective flows and defect dynamics in active nematic films, we demonstrate the existence of a self-constraint, namely a two-way, spontaneously arising relationship between activity-driven isosurfaces of flow boundaries and mesoscale nematic structures. We show that self-motile defects are tightly constrained to viscometric surfaces, which are contours along which the vorticity and the strain rate are balanced. This in turn reveals that self-motile defects break mirror symmetry when they move along a single viscometric surface. This is explained by an interdependence between viscometric surfaces and bend walls, which are elongated narrow kinks in the orientation field. These findings indicate that defects cannot be treated as solitary points. Instead, their associated mesoscale deformations are key to the steady-state coupling to hydrodynamic flows. This mesoscale cross-field self-constraint offers a framework for tackling complex three-dimensional active turbulence, designing dynamic control into biomimetic materials and understanding how biological systems can employ active stress for dynamic self-organisation.


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.

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 Q Q Q to be inferred through coherence-enhanced diffusion filtering 23 and the velocity field ⃗ 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 Q Q Q ( § 8.3.1).Likewise, active flow fields have mesoscale structure with well-studied vorticity scales, which are visualized using the Q-criterion 25 where Ω Ω Ω and E E E are the vorticity and strain-rate tensors ( § 8.2).Vorticity-dominated regions are identified by Q > 0, while 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 Q-criterion, and director structure has yet to be fully understood.

Plus-half defects are constrained to viscometric surfaces
By partitioning flow domains into vorticity-dominated (Q > 0) and strain-rate dominated (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 Q = 0 isolines ( § 8.3.2),we identify that +1/2 defects are always found precisely on the border where Q = 0 (Fig. 1).We refer to the lines of 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 vorticitydominated regions, which exhibit a broad distribution of sizes, shapes and circulation (Ext.Data Fig. 1).Simulations reveal 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 Q = 0 boundaries (Fig. 1d-e), with the motile +1/2 defect pinned to the Q = 0 viscometric surface throughout the process, while the −1/2 defects are unbound from 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 ω direction inverts ( § 8.3.3), it is seen that most-but-not-all +1/2 defects exist on |ω| = 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 Q = 0; however, this is because they are initially created at 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 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 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 Q-criterion (Fig. 2a; Supplementary Video 2).In contrast, the striking coincidence of +1/2 defects on Q = 0 flow contours suggests an intrinsic relationship between flow structures and the trajectory of the self-propelled defects.

Spontaneous self-constraint violates the ideal, solitary view of defect-generated flow
Investigating the relationship between the viscometric (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 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 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 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 θ ≈ 20 • , representing the minority-population in the transient symmetric configuration.The measured angle is smaller than θ ≈ 35.26 • , 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 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.

Bend walls and viscometric surfaces are interdependent
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 SB , which indicates the difference between splay and bend deformations ( § 8.3.5).The splay-bend parameter traces bend walls through strongly negative S SB values (Fig. 3a; Supplementary Video 4), and reveals that bend walls correlate strongly with the Q = 0 viscometric lines.The distribution of S SB coinciding with viscometric surfaces is more skewed towards negative (bend) values than for the system as a whole (Fig. 3b).

Physical Model
The extended concurrence of bend walls and 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 ±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 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 Q = 0; however, this overly simplified model predicts 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: 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 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 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 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 Q = 0 lines collapse into a single line but simultaneously Q goes to zero everywhere.Case (ii) instead considers a perfectly circular bend wall (Fig. 3d) for which a 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 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 Q = 0. Finally, the finite length and curvature of the bend walls causes the mirror-symmetry-broken case.

Beyond extensile, flow-aligning 2D active nematic turbulence
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 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 Q = 0 and, thus, these effectively +1 compounds do not reside where 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 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 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 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 Q = 0 intersection (Fig. 4d; Supplementary Video 10).Similarly, in experiments +1/2 defects weave through the center, never departing from 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 Q = 0 intersection, allowing the defects to transfer to the next vortex of opposite handedness-all the while, the defects remain on the 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 Q = 0 isosurface (Supplementary Video 11), demonstrating that motile defects are constrained to lie on 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.

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 selfconstraint 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.

Materials
Short and stable microtubules (MTs) were polymerized by incubating, at 37 • C for 30 minutes, a mixture containing 8 mg mL −1 of recycled tubulin from bovine brain (Brandeis Materials Research Science and Engineering Center), 6 mM of guanosine-5'-[(α,β )-methyleno]triphosphate (GMPCCP, Jena Biosciences, NU-405S), and 1 mM of DL-dithiothreitol (DTT, Sigma, 43815) in M2B buffer (80 mM of PIPES [Sigma, P1851], 1 mM of EGTA [Sigma, E3889], 2 mM of MgCl 2 , pH=6.8).After incubation, the solution was kept at room temperature for 5 hours, then frozen in liquid nitrogen and stored at −80 • C 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 mM imidazole aqueous buffer.The kinesin concentration was estimated at 2.5 µM by absorption spectroscopy, and was finally stored in a 40%w/v sucrose solution at −80 • C for future use.Kinesin clusters were obtained by preparing a mixture with 1 µM of biotinylated kinesin dimers and 31 µg mL −1 of tetrameric streptavidin in a M2B buffer supplemented with 0.2 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 mg mL −1 of MTs, motor clusters (the concentration of streptavidin was 8.2 µg mL −1 ), 1.4 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 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 mM of DTT, 3.3 mg mL −1 of glucose (Sigma; G8270), 38 µg mL −1 of catalase (Sigma, C40), 0.22 mg mL −1 of glucose oxidase (Sigma G2133) and 2.0 mM of trolox (Sigma, 238813).The MgCl 2 concentration was raised with 4.7% v/v of Mix solution (69 mM of MgCl 2 in M2B).The pH of the active gel was adjusted to 6.8.

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 µm of double-sided tape, forming a rectangular chamber of size 3 mm×22 mm.The chamber was first filled with fluorinated oil HFE7500 with 1.8% v/v of 008-FluoroSurfactant (RAN Biotechnologies), then 10 µ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×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 cm in diameter) to a polyacrylamide coated glass slide.2.5 µ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 µL of silicone oil with a dynamic viscosity of 20 mPa s.The assembly of the microtubules bundles at the water/oil interface occurred on their own within 40 minutes.

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× 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 nm.The grid thickness was 100 µm, to ensure a good resistance during manipulation.

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× objective (Nikon).The sample was excited with a laser beam at 647 nm with an exposure between 200-500 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 Q Q Q (⃗ r,t) was computed from the images using custom Matlab scripts 23 that infer local alignment of the microtubules using coherenceenhanced diffusion filtering.First, noise in the image is filtered out with a Gaussian blur of standard deviation σ 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 σ 2 .Finally, the nematic tensor Q Q Q is constructed by ensemble averaging the molecular orientation tensor over a small circle of size β around each pixel.The three parameters σ 1 , σ 2 and β were manually adjusted by visual inspection of the resulting director field.Velocity fields ⃗ u (⃗ r,t) were calculated with the optical flow method detailed in ref. 50and 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.

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 Q Q Q (⃗ r,t) = S (⃗ n ⊗⃗ n − δ δ δ /3) at each position⃗ r and time t in 3D is described by the Beris-Edwards transport equation where the right-hand side is the relaxation to equilibrium through a relaxation parameter Γ and molecular field is the molecular field corresponding to the nematic free energy F. On the left-hand side, derivative that accounts for material advection by D t = ∂ t +⃗ u • ⃗ ∇, which couples directly to the velocity field, and a corotational operator The Beris-Edwards equation (Eq.( 2)) is numerically solved by finite difference.
The velocity ⃗ u (⃗ r,t) evolves according to the generalized Navier-Stokes equations with density ρ and friction tensor κ κ κ.The stress Π Π Π consists of a viscous term 2ηE E E with viscosity η, a pressure term −pδ δ δ , an elastic term , and, most relevantly an active term, 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 ζ that leads to a bend instability when ζ > 0 (or a splay instability when ζ < 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 ⃗ ∇ •⃗ u ≈ 0.
The free energy includes a bulk component and a surface anchoring term F = ´V f dV + ´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 with γ controlling the distance from the isotropic-nematic phase transition and K the isotropic elastic constant.The surface free energy density is ) where W 1,2 are anchoring energies, 6S eq = 1+3 1 − 8/3γ is the equilibrium scalar order parameter and Q Q Q ⊥ = P P P• Q Q Q + S eq δ δ δ /3 • P P P for the surface normal projection operator P P P.
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 λ = 0.6.Contractile (Fig. 4b; Supplementary Video 6) Substitute ζ = −0.1 and system size 160×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 κ κ κ = [[0.015,0] , [0, 0]].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 × 25. 3D channel confinement (Fig. 4g; Supplementary Video 11) Same as the 2D channel wall but now in a 3D square duct.

Defect Analysis
Topological disclinations in both experiments and 2D simulations are identified by calculating the topological charge density [53][54][55] 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 with ⟨•⟩ • denoting averaging a closed loop around the disclination 56 .In 3D, defect lines are identified and characterised using the disclination density tensor 53 , D D D, constructed from gradients in the nematic in terms of a non-negative scalar field s, the rotation vector ⃗ R, and tangent vector ⃗ 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 β = ⃗ R • ⃗ T , where β is known as the twist angle.Values of cos β 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 .

Q-criterion
We characterize the principle behavior of flow fields by the Q-criterion, which is defined to be the second invariant of the velocity gradient In the incompressibile limit, tr [L L L] = ⃗ ∇ •⃗ u = 0 and so this definition reduces to Eq. ( 1)), which we repeat here for convenience where the squared norms are ∥A A A∥ 2 = A i j A i j .The Q (⃗ r,t) is shown in Supplementary Video 12 and iso-surfaces of Q = 0 are identified using the skimage python library 59 .

11/21
If Q > 0, vorticity dominates over strain-rate, where as if Q < 0, strain-rate dominates.Where Q = 0, the magnitudes of vorticity and strain-rate are equal.In 2D, isolines of Q = 0 must form closed loops, and, in 3D, Q = 0 contours must form closed iso-surfaces enclosing vorticity dominated regions.In 2D, there are only two invariants of L L L, and the first tr [L L L] = ⃗ ∇ •⃗ u = 0 due to incompressibility.Therefore, if Q = 0, all the invariants of L L 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 Q = 0 contours as viscometric surfaces.

Comparing zero-isosurfaces Q with vorticity
Throughout this work, we employ the zero-isolines of 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 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 Q-criterion, we explore contours where contributions from Ω Ω Ω go to zero.In any two-dimensional coordinate basis, Ω Ω Ω only has off diagonal components Ω 21 = −Ω 12 = ω in terms of the vorticity pseudo-vector ⃗ ω.Hence, locations where the rotation-rate tensor vanishes can be simply obtained from a sign change in ω.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 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.

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 Q > 0).After identifying the list of points defining the Q = 0 boundary, vortices are identified as the enclosed regions where 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 Q > 0. We calculate the square of the relative shape anisotropy κ 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 The enclosed circulation is defined as 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 Γ 0 = ω 0 ℓ 2 ζ .The vorticity scale, ω 0 , was identified as the standard deviation of the vorticity distribution.We calculate the active length scale as where the defect number density σ = N +1/2 + N −1/2 / T L 2 for N ±1/2 nematic defects, total timesteps T and system size L.

Splay-bend parameter
The splay-bend parameter 37, 63 is used to visualize bend walls as line-like structures with large negative values 37 .To interpret S SB , consider a uniaxial 2D nematic with constant scalar order S.Under these conditions The first term inside the divergence is the splay and the second term is the bend, leading to the interpretation that S 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)) in active nematics the splay-bend parameter can be understood to be proportional to the divergence of the active force 12/21

Defects and viscometric surfaces
To find the probability distribution of nearest distances between defects and viscometric surfaces (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 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 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 where⃗ r( j) is the position vector of the chosen point j on the 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 • and 90 • .If the defect is associated to the mirror-symmetric 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 ⟨Q⟩ ∆r 2 as the average Q-criterion in a small region ∆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 ⟨Q⟩ ∆r 2 three different defect populations can be identified: one group where ⟨Q⟩ ∆r 2 is close to zero (with mirror symmetry), and two where ⟨Q⟩ ∆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 ⟨Q⟩ ∆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 ∆r 2 .Ext.Data Fig. 2a shows the ensemble-averaged case of mirror-symmetric defects at the 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 Q > 0 vortex to the right of the defect head and Q < 0 to the left and vice-versa in Ext.Data Fig. 2c.In both cases, the defect aligns parallel with the 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.

Stokesian solitary defect model
The Q-criterion and flow field around +1/2 and −1/2 topological defects is shown in Fig. 2b and f.The Q-criterion field is determined from the velocity field solutions ⃗ u ± for a single solitary ±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 (kφ ) x + sin (kφ ) ŷ 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 ⃗ u = ⃗ 0) in the absence of friction and with negligible elastic stresses, leaving the Stokes equation for the velocity field, as a response to the active force density ⃗ f ζ (Eq.( 14)).For an extensile plus half defect (k = +1/2), the active force is in the direction of the defect − 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 with the two-dimensional Oseen tensor 28 13/21 in a system bound by the length scale L. Substituting in the active force for the fixed, solitary defects gives the velocity fields where ζ is the activity, η 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 L L L = ⃗ ∇ ⊗⃗ u and the Q-criterion field is calculated as the second invariant of L L L. These calculations and corresponding plots were performed using Mathematica 64 .For the +1/2 defect case, the Q-criterion takes the form, In this construction, the +1/2 defect is aligned with a self-propulsion direction in the (negative) x direction, with an angle θ = sin −1 3 −1/2 ≈ 35.26 • to the 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.

Stokesian straight bend-wall model
Since our experimental and numerical results indicate that bend walls are coincident with 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 ⃗ n = cos θ x + sin θ ŷ, where the angle θ is for a modulation length scale λ .This generates the force density field ⃗ f ζ = ζ π λ sin 2y λ (cos 2θ x + sin 2θ ŷ).Since the system has translational symmetry in x, only derivatives with respect to y exist for the pressure and velocity.Through the incompressibility equation, ∇ •⃗ u = 0, the absense of x−gradients requires ∂⃗ u y ∂ y = 0.This leads to only one non-zero contribution to the velocity gradient tensor, ∂⃗ u x ∂ y , which is obtained through integrating the x−momentum equation in Eq. ( 17).The velocity gradient tensor is found as, Since this is upper diagonal at all points (representing simple horizontal shear), it is a viscometric flow and Q = 0 everywhere.Thus, a straight, infinitely long series of kink walls is insufficient to explain the colocation of bend walls and Q = 0.

Wavy bend wall model
Next we consider if bend wall curvature explains their coincidence with Q = 0.The bend walls from § 8.4.2 are perturbed into infinitely long sinusoidal waves with director field once again of the form ⃗ n = cos θ x + sin θ ŷ, where the angle θ but with angle This director field generates a periodic sequence of bend and splay walls that modulate sinusoidally about slices of ŷ (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 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 Q = 0 Overall, this model demonstrates that modulation in bend walls is able to split distinct Q regions without needing defects.The geometry of the flow generated by the bend has 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.

Stokesian line force model of bend walls
The hydrodynamic instability constricts infinitely long bend into narrow kink walls and their curvature facilitates closed 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) 15,32 .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 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 × 1000 grid is used and 1000 × 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 ζ x at the origin produces a 2D Stokeslet Taking the gradient results in the Q-criterion which predicts an intersection of viscometric surfaces at an angle θ = 45 • .
(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 x direction.The velocity field and zero isolines of the Q-criterion (Eq.( 1)) are shown in Fig. 3c.The Q-criterion retains the mirror symmetry observed for the solitary defect or far-field point source, but now differs by constraining the alignment angle θ 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π.Positions on the arc are⃗ r ′ (s) = R (cos s x + sin s ŷ), where R is the radius.As before, the force contour is assumed to have constant magnitude f ζ and orient azimuthally.By choosing the forces to follow the arc in a clock-wise sense, the force is ⃗ f ′ ζ (s) = f ζ (sin s x − cos s ŷ).For the arc, numerical integration is performed using 500 point sources distributed between s = −π/2 and s = π/2, while the circle uses 1000 point sources between s = 0 and s = 2π.The results for the velocity field and Q-criterion are shown in Fig. 3d (circle) and e (arc).For both cases, the Q-criterion goes to zero when |⃗ r| = R.This result (Fig. 3e) shows that this simplified model can capture the correlation between bend walls and 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.

Figure 1 .
Figure 1.Plus-half topological defects reside on isolines of 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 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 µ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 ±1/2 defects to 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.

Figure 2 .
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 Q = 0 border.Arrows show instantaneous velocity field, dashed lines director field, solid lines 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 (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).

Figure 3 .
Figure 3. Viscometric surfaces generated by bend walls.(a) Same snapshot as Fig. 1b with colormap corresponding to the scalar-bend bend parameter S SB ( § 8.3.5).Narrow lines of strong negative S SB are interpreted as bend walls.(b) Probability distribution function (PDF) of S SB sampled equally for all points and times, against the distribution of S SB only on Q = 0 contours.(c)-(e) 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 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.

Figure 4 .
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 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 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 µm.(g) Three-dimensional vortex lattice with disclinatines colored by the twist angle (cos β 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 (Q > 0), and white where strain-rate dominates (Q < 0).
co-rotational term couples the orientation to the velocity gradient L L L = ⃗ ∇ ⊗⃗ u through the vorticity Ω Ω Ω = L L L − L L L T /2 and strain-rate E E E = L L L + L L L T /2 tensors, written using a (anti)commutator {A A A, B B B} ± = A A A • B B B ± B B B • A A A, with a flow-aligning parameter λ and Q

15 / 21 Extended Data Figure 1 .Extended Data Figure 3 .Extended Data Figure 4 .
Probability distribution functions (PDF) of vortex properties, shown for experiments and simulations.(a) Shape anisotropy parameter κ 2 (Eq.(11)).The limiting value κ = 0 represents a circle and κ = 1 indicates a line.(b) Vortex area non-dimensionalized by the active length scale ℓ ζ .(c) Total circulation within a Q = 0 closed contour ( § 8.3.4).Extended Data Figure 2. Time and ensemble averaged director fields and Q-criterion around +1/2 defects sorted by handedness.(a) Mirror-symmetric state showing Q-criterion structure around +1/2 defect.The 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 Q = 0 are shown as solid lines and colormap corresponds to vorticity field.Infinitely long, repeating wavy bend wall model ( § 8.4.3).(a) Initial splay-bend parameter (S SB ) field with repeating model bend walls in purple, with same colorbar as Fig. 3a.(b) Resulting Q of steady-state active flows due to fixed S SB from (a).(c) Resulting steady-state S SB due to fixed flows from (b).(d) Resulting steady-state flows due to fixed S SB from (c).Director field is shown as black lines, the velocity field is shown as grey arrows and contours where Q = 0 are shown as solid lines.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, Γ, 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 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 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.