Antipolar ordering of topological defects in active liquid crystals

ATP-driven microtubule-kinesin bundles can self-assemble into two-dimensional active liquid crystals (ALCs) that exhibit a rich creation and annihilation dynamics of topological defects, reminiscent of particle-pair production processes in quantum systems. This recent discovery has sparked considerable interest but a quantitative theoretical description is still lacking. We present and validate a minimal continuum theory for this new class of active matter systems by generalizing the classical Landau–de Gennes free-energy to account for the experimentally observed spontaneous buckling of motor-driven extensile microtubule bundles. The resulting model agrees with recently published data and predicts a regime of antipolar order. Our analysis implies that ALCs are governed by the same generic ordering principles that determine the non-equilibrium dynamics of dense bacterial suspensions and elastic bilayer materials. Moreover, the theory manifests an energetic analogy with strongly interacting quantum gases. Generally, our results suggest that complex nonequilibrium pattern-formation phenomena might be predictable from a few fundamental symmetry-breaking and scale-selection principles.

Understanding the emergence of such topological super-structures is crucial for the development and control of new materials, as recently demonstrated for colloidal liquid crystals [49][50][51].
We here develop and test a closed continuum theory for dense ALCs by generalizing the higher-order scalar and vector theories of soft elastic materials [52] and bacterial fluids [15,41] to matrix-valued fields.Specifically, we propose a modification of the commonly adopted Landau-de Gennes (LdG) free energy to account for the inherently different microscopic buckling behaviors of passive and active LCs [5].While bending is energetically unfavorable in passive LCs and hence penalized by the LdG energy functional, kinesin-driven ALCs buckle spontaneously even at low concentrations due to the extensile motor action (figure 1(a)).This experimental observation [5] implies that the classical LdG framework is, by construction, ill-suited to describe experiments in which microtubule bundles are sheared relative to each other by motor proteins [5,47,48].The inclusion of active bending effects in the LdG functional yields a tensor version of the Swift-Hohenberg theory [53] of pattern formation.The resulting minimal model has only two dimensionless parameters, thus allowing a detailed comparison with recent experimental data [5,48].
In addition to the traditional Q-tensor formulation, we present an equivalent complex scalar field representation [14,54] that manifests an analogy with a generalized Gross-Pitaevskii theory [55,56] of strongly coupled many-body quantum systems [57][58][59][60].In the case of normal dispersion, the celebrated LCsuperconductor correspondence [54,61] has helped elucidate profound parallels between the smectic phase in passive LCs and the Abrikosov vortex lattices in type-II superconductors [62,63].The results below indicate that a similar analogy may exist between ALCs and Bose-Einstein/Fermi condensates with double-well dispersion [58][59][60], suggesting that ALCs could offer insights into the dynamics of these quantum systems and vice versa.

Experimental conditions
Recent experiments [5,48] show that ATP-driven microtubule-kinesin bundles can self-assemble into a dense quasi-2D ALC layer at a surfactant-supported oil-water interface parallel to a planar solid boundary(figure 1(b)).This 'wet' ALC was found to exhibit local nematic alignment of bundles, persistent annihilation and creation dynamics of topological defects [5], and remarkable nematic order of the defect orientations in thin layers [48].Although a large number of unknown parameters has prevented detailed quantitative comparisons between theory and experiment, several recently proposed multi-order-parameter models of 2D ALC systems [14,20,64,65] were able to reproduce qualitatively selected aspects of these observations, such as defect-pair creation and separation [65].Despite providing some important insights, traditional models often do not account for three relevant details of the experiments [5,48].First, those models typically assume divergence-free 2D fluid flow within the ALC layer, which is a valid approximation for isolated free-standing film experiments [66] but neglects fluid exchange between the 2D interface and bulk in the ALC experiments (figure 1(b)).Indeed, the surfactant-stabilized interface causes the microtubule-kinesin bundles to assemble into a quasi-2D layer, but places no such constraint on the fluid.As is known for classical turbulence [67,68], small-scale energy input can trigger turbulent upward cascades in incompressible 2D flow.Thus, topological defect dynamics in the current standard models may be dominated by artificially enhanced hydrodynamic mixing due to a simplifying 2D incompressibility assumption that is unlikely to hold under realistic experimental conditions [5,48].Second, a relevant yet previously ignored effect is damping from the nearby boundaries, which may promote topological defect ordering.Third, as already mentioned above, the commonly adopted standard LdG free-energy functional does not account for motor-driven spontaneous buckling [5] of microtubule bundles (figure 1(a)), which is one of the key differences between passive and active LCs (Z.Dogic, private communication).To overcome such limitations and achieve a quantitative description of the experiments [5,48], we next construct a closed continuum theory for ALCs described by a nematic tensor field ( ) r Q t, .The theory accounts for the different buckling behaviors of passive and active LCs and builds on a self-consistent hydrodynamic closure condition.

Theory
Traditional multi-field models [64,65] aim to describe the 2D nematic phase of a dense ALC suspension by coupling the filament concentration ( ) r c t, and the nematic order tensor ( ) x y , .The nematic order parameter ( ) r S t, is proportional to the larger eigenvalue of Q, and the filaments are oriented along the corresponding eigenvector, or director ( ) n r t, .To construct an alternative closed-form theory for the symmetric traceless 2×2-tensor fieldQ, we start from the generic transport law AB BA , the commutator of two matrices and an effective free energy.Focussing on dense suspensions as realized in the experiments [5, 48], we neglect fluctuations in the microtubule concentration,  º c 0. A derivation of the advection term  • ( ) vQ from the probability conservation laws underlying generic advection-diffusion models is outlined in the supplementary information.It is important, however, that which is typically the case when fluid can enter and leave the interface.Combining LdG theory [69] with Swift-Hohenberg theory [53], we consider the effective non-equilibrium free-energy density (supplementary information) with > a b , 0for the nematic phase.Assuming g 2 can have either sign, ultraviolet stability requires g > 0 4 . For g < 0 2 , F penalizes bending and buckling, as is appropriate for passive LCs and possibly 'dry' shaken nematics [44,45], forcing the system dynamics towards a homogeneous nematic ground-state manifold.By contrast, motor-induced spontaneous buckling [5] (figure 1(a)) of kinesin-driven ALCs demands g > 0 2 , and consequently patterns of characteristic wavelength g g L ~4 2 become energetically favorable, as shown below.
The requirement g > 0 2 for ALCs has an intrinsically microscopic origin, as the ALC assembly consists of microtubules that grow against each other and spontaneously buckle due to the motor-induced extensile shear dynamics of adjacent bundles (figure 1(a)).To explain this important point in more detail, let us recall that passive LCs are modeled using g < 0 2 , as the corresponding term in the free energy penalizes variations in Q and thus inhibits spatial inhomogeneities at damping rate g in Fourier space.Microscopically, the alignment dynamics of two rod-like passive LC molecules roughly corresponds to a time-reversal of the ALC pairinteraction image sequence in figure 1(a), implying that a corresponding ALC system would develop buckling instabilities at growth rate g . Additional empirical support for this theoretical picture comes from a comparison of the experimentally observed length scales in ALC systems: the microtubule-kinesin bundles realized in the ALC experiments [5,48] are approximately 10-30 mm in length (figure 1(b) in [5]), which is on the same scale as both the spontaneous buckling wavelength (figures 1(c) and (d) in [5], and figure 1(a)) and the typical separation distance between defects (bottom panels of figure 3(d) in [5], reproduced in figure 4(a) below).That is, there exists no substantial scale separation between the buckling microscopic constituents and emergent ALC dynamics.
Similar buckling phenomena are generically observed in many systems that are subjected to external or internal stresses, for example in elastic films and sheets [52,70] and in geometrically confined cellular networks [71,72].The ALCs experience an effective compressive stress due to the extensile 'growth' of the filament pairs in a confined geometry, which arises from their motor-induced shear dynamics.Application of such a compressive stress leads to buckling of the network's constituents [71,72].It has been shown that out-of-plane buckling of an elastic sheet due to an effective compressive stress may be quantitatively modeled by a Swift-Hohenberg-type equation with g > 0 2 in the corresponding free energy [52].We expect the same to be true for the in-plane buckling of microtubules confined to a planar interface, and thus analyze here an effective theory that incorporates this spontaneous motor-induced buckling phenomenologically through g > 0

Hydrodynamic closure
To obtain a closed Q-model, we relate the 2D flow field v to Q through the linearly damped Stokes equation [73,74] where η is the viscosity and the rhs represents active stresses [10,64] with z > 0 for extensile ALCs(figure 1(c)).A pressure term does not appear in equation (3) because the interfacial flow is not assumed to be incompressible and concentration fluctuations are neglected.The ν-term in the force balance (3) has been used to model interfacial damping in other contexts, such as surfactant membranes on a solid substrate [73], and accounts for friction from the nearby no-slip boundary in the Hele-Shaw [74] approximation (figure 1(b)).In the overdamped regime n h L 1 2  , we deduce from equation (3) the closure condition Equation ( 4) is conceptually similar to closure conditions proposed previously for active polar films [75].Importantly, equation (4) predicts divergent interfacial flow, ¹ • v 0, and hence fluid transport perpend- icular to the interface wherever  ¹ Q : 0. Inserting(4) into(1) yields a closed Q-theory in which periodic director patterns corresponding to local minima of the free energy  can become mixed by self-generated interfacial flow.

Complex representation and ALC-quantum analogy
The traditional characterization of 2D nematic order in terms of the symmetric traceless 2×2 matrix field , is redundant, for only two real scalar fields l ( ) r t, and m ( ) r t, are needed to specify the nematic state at each position = ( ) r x y , .To obtain an irreducible representation [14,54] we define the complex position coordinate = + z x y i , velocity field = + ( ) v t z u w , i and complex order parameter , the 2D Laplacian takes the form  = ¶ ¶ 4 where the self-advection operator is given by and the free energy For g < 0 2 and g  0 4 , equation (7) reduces to the energy density of the Gross-Pitaevskii mean-field model [55,56] for weakly interacting boson gases.Historically, this limit case has been crucial [54,61] for elucidating the analogy between the smectic phase of passive LCs and the Abrikosov flux lattice in type-II superconductors [62,63].For g g > , 0 2 4 , equation (7) effectively describes double-well dispersion [57], as recently realized for quasi-momenta in spin-orbit-coupled Bose-Einstein condensates [58,60] and Fermi gases [59].This fact establishes an interesting connection between dense ALCs and strongly coupled quantum systems: when selfadvection is negligible  ( ) D 0 , the fixed point configurations of equation ( 5) coincide with the 'eigenstates' of generalized Gross-Pitaevskii models that incorporate wavelength selection.

Stability analysis
The qualitative model dynamics is not significantly altered for moderate values of κ (Movies S1 and S2), so we neglect the commutator term by setting k = 0 from now on (see supplementary information for k > 0).To understand the properties of equations ( 5) and (6) when self-advection is relevant, we perform a fixed point analysis of the rescaled dimensionless equation (supplementary information) by focussing on the uniform state y = (supplementary information).For subcritical self-advection, , the dominant instability is driven by modes with wavenumber , suggesting the formation of stripe patterns with typical wavelength p g g L » 8 2  4 2 .By contrast, for supercritical advection, > D D c , the most unstable mode propagates perpendicular to the director, , suggesting the possibility of transverse mixing.

Phase diagram
To investigate the nonlinear dynamics of equation (8), we implemented a Fourier pseudospectral algorithm with modified Runge-Kutta time-stepping [76] (Methods) and so evolved the real and imaginary parts of y ( ) t z , in time for periodic boundary conditions in space.A numerically obtained g ( ) D , 2 -phase diagram for random initial conditions confirms the existence of a turbulent nematic phase if active self-advection is sufficiently strong (figures 2(a) and (b); Movie S1).Ordered configurations prevail at low activity (figures 2(a), (c) and (d); Movies S3, S4 and S6).Although the ground-states of the free energy (2) are in general not homogeneous, the critical curve separating the two regimes is in fair agreement with the estimate 2 from linear stability of the homogenous state (white line in figure 2(a)).For subcritical values of the advection parameter D, we observe either defect-free ground-states or long-lived lattice-like states exhibiting ordered defect configurations (figures 2(c) and (d)).Regarding the subsequent comparison between theory and experiment, it is important to note that the lattices are also found in simulations with a large domain (figure 3; Movie S5).These spatially periodic states generally exhibit antipolar long-range ordering of +

Theory versus experiment
To test our theory systematically against existing experimental data [5,48], we analyze defect-pair dynamics, global defect ordering and defect statistics in the turbulent nematic phase.Spontaneous defect-pair creation and subsequent propagation, as reported in recent ALC experiments [5] and observed in our simulations, are compared in figure 4. In the experimental system [5], a + - ( ) 1 2 -defect pair is created when fracture along incipient crack regions [20] becomes energetically more favorable than buckling.After creation, the + 1 2 -defect moves away rapidly whereas the position of the - 1 2 -defect remains approximately fixed for up to several seconds (figure 4(a)).We note that an asymmetry in the speeds of topological defects has also been observed in passive liquid crystals [77][78][79].Our simulations of the minimal model defined in equation (8) accurately reproduce the details of the experimentally observed dynamics (figure 4(b); Movies S1 and S2).
Another striking and unexplained experimental observation [48] is the emergence of orientational order of + 1 2 -defects in thin ALC layers (figure 5).Using the setup illustrated in figure 1(b), recent experiments [48] demonstrated nematic alignment of + 1 2 -defects in thin ALC layers of thickness h 250 nm (figure 5(a)), whereas thicker ALC layers with m h 1 m showed no substantial orientational order on large scales (figure 5(b)).To investigate whether our theory can account for these phenomena, we tracked defect positions r i (Methods) and defect orientations =    5(c)).In our simulations, this ordering decreases as the effective mixing strength D increases (figure 5(d)), consistent with the experimental results [48] for thicker ALC layers (figure 5(b)).Similar ordering was observed in simulations that incorporated alignment of the nematic field near the horizontal boundaries of the simulation box (supplementary information; figure S5).To quantify the degree of orientational order, we recorded the distances d ij between all + 1 2 -defect pairs (i, j) as well as their relative orientation angles q . The resulting pair-orientation distributions q ( | ) p r , and polar and nematic correlation functions, = á ñ ( ) , are shown in figure 6, with á ñ • r denoting an average over pairs of defects separated by a distance r.The local maxima in the orientation distribution at q = 0 and q p = signal antipolar ordering (figure 6   Lastly, we test the theoretically predicted defect statistics against a separate experimental data set kindly provided by DeCamp et al (private communication).Since our simulations are performed in dimensionless units, there is freedom to choose a characteristic lengthscale l 0 and timescale t 0 .To relate theory and experiments, we determine ( ) l t , 0 0 such that the joint mean speed and mean lifetime of  as the best-match parameters, although nearby parameter values and simulations with k = 1 produce fits of similar quality, corroborating the robustness of the model (figure S4).For - 1 2 -defects, we find adequate agreement between experiment and theory for speed and lifetime probability density functions (PDFs), as evident from figures 7(a) and (b).For + 1 2 -defects, simulation results also agree well with the experimental measurements (figures 7(c) and (d)), but one notices two systematic differences.First, while the peak heights of the PDFs agree within a few percent, experimentally measured speed values for + 1 2 -defects are on average slightly larger than theoretically predicted values (figure 7(c)).Second, simulation data predict a miniscule tailfraction of long-living + 1 2 -defects not detected in the experiment (figure 7(d)).In addition, based on the experimental density estimate of 30 defects mm −2 [48], we find that the defect density at any given time in the 'best-fit' simulation is ~2.3 lower than in the experiments.As discussed below, such deviations can be explained plausibly by specific model assumptions.Taken together, the above results confirm that the minimal model defined by equation (8) provides a satisfactory qualitative and quantitative description of the main experimental results [5,48].

Discussion
Pattern-formation mechanism Equations (2) and( 7) epitomize the idea of 'universality' in spatio-temporal pattern formation, as known from Swift-Hohenberg-type scalar field theories [53,81].The free-energy expressions contain the leading-order terms of generic series expansions in both order-parameter space and Fourier space, consistent with spatial and nematic symmetries.When considering passive systems with a preference for homogenization (g < 0 2 ), it usually suffices to keep only the quadratic gradient terms.By contrast, for pattern forming systems, the coefficient in front of the lowest-order gradient contribution can change sign [52,53], and one must include higher-order derivatives to ensure stability.We hypothesize that the sign change of g 2 is directly related to the motor-induced buckling of microtubule bundles (figure 1(a)), an effect that is not captured by the standard LdG free-energy for passive liquid crystals.In a few select cases, expressions of the form(2) and( 7) can be systematically derived [52,53,82].Generally, one can regard the free-energy expansion (7) as an effective field theory whose phenomenological parameters can be determined from experiments.This approach has proved successful for dense bacterial suspensions [15,41] and now also for ALCs, indeed suggesting some universality in the formation and dynamics of topological defects in active systems.

Nematic defect order
Although g ( ) D , 2 are varied as independent effective parameters in the simulations, they are likely coupled through underlying physical and chemical parameters.For example, it is plausible that a change in ATPconcentration or film thickness would affect bothg 2 and D. The parameter D can also be interpreted as an effective Reynolds number.In our numerical exploration of the g ( ) D , 2 -parameter space, we observe for subcritical advection D either long-lived lattice-like states exhibiting nematically aligned + 1 2 -defects or defectfree ground-states (figures 2 and 3; Movies S3, S4, S5, S6 and S7).Ordered defect configurations correspond to local minima or saddles in the free-energy landscape and have only slightly higher energy than defect-free states (figure 2(c)).When the activity ζ is sufficiently large that advection is marginally supercritical,  D D c , chaotic system trajectories spend a considerable time in the vicinity of these metastable lattice states, which provides a physical basis for the orientational order of defects in thin ALC films [48].For D D c  , the ALC system can access a wider range of high-energy states, leading to increased disorder in the defect dynamics.Although the strongly turbulent regime D D c  requires high time-resolution and is thus difficult to realize in long-time simulations, the inhibition of nematic defect order at larger values of D is evident from the reduced peak heights in figures 6(c) and (d).

Defect statistics
The systematic speed deficit in figure 7(c) likely reflects the overdamped closure condition(4), which suppresses the propagation of hydrodynamic excitations.Since flow in a newly created defect pair generally points from the -1 2 to the + 1 2 -defect, the minimal model(8) can be expected to underestimate the speeds of + 1 2 -defects.This effect could be explored in future experiments through a controlled variation of the thickness and viscosity of the oil film (figure 1(b)).The low-probability tail of long-living + 1 2 -defects in the simulation data (figure 7(d)) may be due to the fact that they can be tracked indefinitely in the simulations but are likely to leave the finite field of view in the experiments.In the future, the minimal theory presented here should be extended systematically by adding physically permissible extra terms [83] to the free-energy, explicitly simulating the full hydrodynamics in equation (3), or incorporating additional terms into equation (1) that account for the interaction between the nematic field and the induced flow [65,83].

Future extensions
The minimal model formulated in equations (1) and (2) can be systematically extended to improve further the quantitative agreement between experiment and theory.For instance, one may append to the right-hand side of (1) an additional fourth-order linear term of the form     denoting the symmetric traceless part of the operator .Such a term only affects the high-wavenumber damping at order k 4 and thus is not expected to alter significantly the results obtained here.We also note that extra terms coupling the nematic field to the induced flow may be added to equation (1), an example being SE, where = is the symmetrized strain rate tensor [25,65].Our above analysis neglected such secondary hydrodynamic effects in the interest of constructing a minimal mathematical theory capable of capturing key experimental observations.Moreover, this simplification may be justified on the physical basis that steric interactions and motor-induced buckling are expected to dominate over flow-alignment effects at high microtubule densities 3 .The effect of microtubule bending may be enhanced by appending a hydrodynamic interaction term proportional to + E [25], since the closure condition , which augments the bending term g -DQ 2 in equation (2).However, the scale separation between the experimentally observed filament buckling wavelength and the flow structures in the isotropic phase at low microtubule concentrations (see figure 1(d) in [5]) suggests that such hydrodynamic effects play a secondary role.A natural next step would be to derive systematically the bending term from a microscopic model of motors and filaments as introduced in [20].Finally, it will be worthwhile to attempt constructing a fully 3D theory for the ALCs and fluid and subsequently project on the quasi-2D interface, although additional assumptions are required then to obtain a closed 2D system of equations.

Conclusions
Recent experimental and theoretical studies showed that fourth-order PDE models for scalar and vector fields provide an accurate quantitative description of surface-pattern formation in soft elastic materials [52] and orientational order in dense bacterial fluids [15,41].Here, we have generalized these ideas to matrix-valued fields describing soft active nematics.The above analysis demonstrates that a generic fourth-order Q-tensor model can shed light on experimental observations in 2D ALCs [48], including the emergence of orientational order of topological defects.Physically, the higher-order generalization(2) becomes necessary because the commonly adopted LdG free-energy, which was designed to describe passive liquid crystals, does not account for the experimentally observed spontaneous buckling of motor-driven ALCs [5].More generally, the fact that three vastly different soft matter systems can be treated quantitatively in terms of structurally similar higherorder PDEs [15,41,52] promises a unified mathematical framework for the description of pattern formation processes in a broad class of complex materials.In addition, the free-energy analogy [54] between dense ALCs and generalized Gross-Pitaevskii models suggests that the self-organization principles [40] of mesoscopic active matter and microscopic quantum systems [57][58][59][60] could be more similar than previously thought.

Numerical solver and defect tracking
To simulate equation (8), we implemented a numerical algorithm that evolves the real and imaginary parts l ( ) r t, and m ( ) r t, of the complex order parameter ψ in time for periodic boundary conditions in space.The algorithm solves equation (8) pseudospectrally in space using = ℓ N 256 or = ℓ N 512 lattice points in each direction and a simulation box of size L = 6π (figures 2 and 4) or L = 18π (figures 3, 5-7).Spectral analysis shows that = ℓ N 256 is generally sufficient to resolve the fine-structure of the numerical solutions(supplementary figure S6).The algorithm steps forward in time using a modified exponential time-differencing fourth-order Runge-Kutta method [76] with time step D - t 2 10 .Simulations were initialized with either a single defect pair or random field configurations l m { ( ) ( )} r r 0, , 0, .Defects are located at the intersections of the zero-contours of l and μ, their positions tracked by implementing James Munkres' variant of the Hungarian assignment algorithm [86] (supplementary information).

Figure 1 . 2 .
Figure 1.(a)Image sequence showing spontaneous buckling of a microtubule bundle (dashed) caused by extensile ATP-driven motor activity, adapted with permission from figure 1(c) in [5].Time interval 11.5 s, scale bar 15 μm.For a passive or contractile bundle, one would expect straightening instead of bending, approximately corresponding to a time-reversal of the depicted sequence.(b)Schematic of the experimental setup reported in [5, 48], not drawn to scale.A thin oil film (thickness m ~3 m) separates a 2D ALC layer (~-0.2 1.0 μm) at the oil-water interface from a solid glass cover.Liquid can be exchanged between the ALC layer and bulk fluid, resulting in compressible 2D interfacial flow that is strongly damped by the nearby no-slip glass boundary and the viscous oil layer.(c)Extensile 2D dipole flow in the interface as predicted by the overdamped closure condition(4) for > D 0 and l l = -( ) Q , 0; 0, with l = -( ) r exp 2 .The central horizontal bar indicates the unit director axis, and background colors the nematic order parameter l S .

Figure 2 . 1 2 1 2
Figure 2. Phase diagram obtained from simulations of equation (8) for one particular set of random initial conditions showing the emergence of turbulent nematic states for supercritical active self-advection.(a)We observe convergence to defect-free stripes (blue, panel (c), Movie S3), long-lived static and oscillatory defect lattice solutions (green, panel (d), Movie S4), oscillatory defect creation and annihilation events (black, Movies S6 and S7), and chaotic dynamics (red, panel (b), Movie S1).The white line indicates the analytical estimate  g g = -+ + ( ) D 2 2 c 2 2 2 g

2 2 .
and the closure condition(4) reduces to y = - ¶ v D Denoting the real and imaginary parts of an operator  by R{ }  and I{ }  , equations (1) and(2) may be equivalently expressed as

1 
corresponds to a nematic order parameter value S = 1and homogeneous director angle θ relative to the x-axis.Considering wave-like perturbations y y = + ˆ( ) • *  t e k r i with | ˆ|  and extensile ALCs with > D 0, one finds that y * is unstable when g > 0 2

1 2 - 1 2
defects(figures 2(d); 3(a)) accompanied by vortex flow lattices (figure3(c)).Numerical free-energy calculations show that defect-free states (figure 2(c)) typically have slightly lower energies than the lattice states (figure2(d)), leaving open the possibility of a very slow decay of the latter.However, regardless of whether such lattice states are extremely long-lived metastable or truly stable states, these simulation results confirm that anti polar ordering of + -defect pairs can persist over experimentally relevant time-scales.

Figure 3 . 1 2 1 2 2 .Figure 4 . 2 . 1 2 1 2
Figure 3. Vortex lattice states with antipolar long-range ordering of nematic defects, see also Movie S5.(a) Long-lived nematic order parameter field with periodically aligned -1 2 -defects (black) and + 1 2 -defects (red).(b) Line integral convolution (LIC) plot of the corresponding director field as a proxy for microtubule-bundle patterns.(c) LIC plot of the corresponding fluid velocity field, colorcoded by normalized vorticity, demonstrates the formation of a vortex flow lattice.Dimensionless simulation parameters are g = = D 0.25, 1.875 2 (a)), which is also reflected in the oscillatory behavior of the polar and nematic correlation functions (figures 6(c) and (d)).The diminished intensity of the local maxima for larger values of D indicates that enhanced hydrodynamic mixing reduces orientational order (figures 6(c) and (d)).

Figure 6 . 1 2 1 2
Figure 6.Increasing activity and film thickness decreases antipolar ordering in simulations.(a), (b) Maxima of the numerically obtained local pair orientation PDFs q ( | ) p r ij signal antipolar local ordering of + 1 2 -defects as they are separated by the typical defectlattice spacing.The defect distance r is specified in units of the mean nearest-neighbor distance r 0 between + 1 2 -defects.(c), (d)Polar P (r) and nematic N(r) correlation functions for D=1.5 (red) and D=3 (blue).Increasing the effective hydrodynamic couplingD leads to stronger mixing and hence decreases nematic order, which is corroborated by the nematic correlation length being ~40% shorter for D=3 than for D=1.5 (panel (d)).This is also reflected by the diminished intensity of the maxima in panel (b) relative to panel(a).The simulation parameters correspond to those given in figures 5(c) and (d).

, 2 -
s 1 and t = ¯52.8s.After fixing these global scales, we can compare details of the speed and lifetime distributions (figure7).To this end, we first locate the 'best-fit' simulation parameters in the g ( ) D parameter space explored in the phase diagram (figure2(a)).This procedure identifies g

Figure 7 . 1 2− 1 . 1 2
Figure 7. Quantitative comparison of defect statistics between predictions of equation (8) and experimental data [48], using the parameter estimation procedure described in the text.For -1 2 -defects, both (a) speed distribution and (b) lifetime distribution agree well.(c) For + 1 2 -defects, experimentally measured speed values are slightly larger, as our model assumes a strongly overdamped limit.(d)Simulations with periodic boundary conditions (Movie S8) predict a low-probability tail of large lifetimes which is not visible in the experiment, likely due to its restricted field of view or additional noise.Dimensionless simulation parameters D=1.75 and g = 12