Dry active turbulence in a model for microtubule-motor mixtures

We study the dynamics and phase behaviour of a dry suspension of microtubules and molecular motors. We obtain a set of continuum equations by rigorously coarse graining a microscopic model where motor-induced interactions lead to parallel or antiparallel ordering. Through numerical simulations, we show that this model generically creates either stable stripes, or a never-settling pattern where stripes periodically form, rotate and then split up. We derive a minimal model which displays the same instability as the full model, and clariﬁes the underlying physical mechanism. The necessary ingredients are an extensile ﬂux arising from microtubule sliding and an interfacial torque favouring ordering along density gradients. We argue that our minimal model uniﬁes various previous observations of chaotic behaviour in dry active matter into a general universality class. comprising particles so – corresponding equations is cur-rently research For systems with orientational order active liquid crystals), two important classes of

Here we propose a new universality class of dry active systems defined by the following continuum equations * ∂ t ρ =∇ 2 1 32 ∂ t Q i j = 4 ρ ρ cr − 1 − αQ kl Q kl + κ∇ 2 Q i j where a tensorial field Q i j quantifies the nematic (apolar) ordering of the active particles, and ρ is their density; i, j = {x, y} denote the two-dimensional Cartesian components, and D i j = ∂ i ∂ j − (1/2)δ i j ∂ k ∂ k . We demonstrate that these equations arise naturally from a microscopic kinetic theory of model mixtures of microtubules (MTs) and molecular motors (MMs) 1,11 . While wet incompressible active gels are generically unstable to orientational fluctuations ultimately resulting in "active turbulence" 7,12,13 , we show that compressible dry MT-MM mixtures undergo seemingly similar chaotic dynamics, which we name dry active turbulence. The underlying mechanism is, however, completely different and we use Eqs. (1) and (2) to elucidate it: importantly, in our case concentration inhomogenenities play a central role. We note that similar dynamical patterns were observed, mainly at the level of a kinetic theory, for a very different physical system, a suspension of flocking self-propelled particles with nematic alignment 8,[14][15][16] . We therefore argue that dry active turbulence unifies various previous observations of chaotic behaviour in nematically ordered microtubules and flocking self-propelled particles.
Here, we consider the dynamics of pattern formation in model MT-MM mixtures that have actively been studied as biological and synthetic instances of active matter 1,11,[17][18][19][20] . On the one hand, they incorporate the essential ingredients of the mitotic spindle [21][22][23] , on the other hand, they closely mirror the so-called "hierarchical active matter", which can be self-assembled in the lab from MTs and MMs, in the presence of polyethylene glycol [24][25][26] . Whilst the continuous description of the overdamped active biofilaments can be postulated on symmetry grounds 9,27 , it can also be derived by rigorously coarse-graining a specific underlying microscopic model [28][29][30][31][32][33][34][35] . This avenue -which has also been often used for systems of self-propelled particles 36 -is useful as it allows one to determine the effective parameters of the continuum theory in terms of the geometrical and physical quantities appearing in the microscopic model. Here, we follow this approach and consider a caricature model of a MT-MM mixture, inspired by actin-myosin 37,38 and MT-kinesin/dynein 24,25,[39][40][41][42] systems. While missing many specific biological ingredients, our model instead focuses on the symmetries of the microscopic interactions that, ultimately, define the universality class of Eqs. (1) and (2). Since most of the experiments either studied MT-MM mixtures directly on surfaces 24,[37][38][39][40][41] or observed spontaneous clustering of bulk systems unto interfaces 25 , we consider our model in 2D.
We treat MTs as rigid rods of fixed length l with distinct ends, denoted as "+" and "−", and consider "+"-directed MMs, which are described by their distribution along individual MTs, as detailed below. Unlike in MT motility assays 43  The nature of motor-mediated interactions between filaments depends on the underlying microscopic detail. For instance, it was demonstrated that MMs able to associate with two filaments simultaneously can either cluster MTs, or actively separate them, depending on the initial configuration ( Fig. 1). Importantly, our dynamics captures both possible outcomes -consistent with the current view of most kinesin motors 41,44 , and unlike previous work 34 , which solely focused on the case of polar clustering. Specifically, our interaction rule is the following. If the initial relative angle between rods exceeds some critical value (in our case π/2) then MTs first align in an anti-parallel way and then slide apart. Otherwise, MTs cluster to acquire the same position and orientation (Fig. 1a).
Within our model, MTs are covered by a steady-state static distribution of motors. Motor coverage may either be homogeneous or inhomogeneous 45,46 , and is parametrized by two geometrical quantities: Ξ = m + − m − /(2m − ), and τ 0 (m + and m − denote motor concentrations at the two ends of a MT, Fig. 1b). Hereafter, we refer to Ξ = 0 and Ξ = 0 as the isotropic and anisotropic cases, respectively.
In what follows, we work in two-dimensional Cartesian coordinates. Assuming that motor-induced rearrangements of MTs are fast with respect to diffusion, we treat them as instantaneous collisions. The probability distribution function, P(r, φ ), for a MT to be at a position r with an orientation n = (cos φ , sin φ ), given by the angle φ , obeys the following Boltzmann-like kinetic equation The first two terms on the r.h.s. of Eq. (3) represent contributions from diffusion, where D r is the rotational diffusion coefficient and D i j are components of the translational diffusion tensor 34,47 . The rest of the equation encodes our collision rules between MTs, and includes clustering and sliding, see Figure 1a. The positions of the colliding MTs are given by r 1,2 = r ∓ ξ 2 , while their orientations are defined by the angles φ 1,2 = φ ∓ ω 2 ; ξ and ω parametrise separations between MT centres and their orientations, respectively. The parameter η determines the final relative displacement of MTs after sliding -henceforth we consider η = 1, corresponding to full separation. For needle-like MTs considered here, the collision rates W + 1 , W + 2 , W − only differ from zero when two MTs intersect in 2D; see 34,48 for their explicit dependence on ξ , ω, Ξ, and τ 0 . Also, we note here that while the value π/2 of the critical angle separating clustering and sliding dynamics of MTs is hard-wired into our theory, some idea of phase behaviour corresponding to different values of that angle can be obtained by using η = 1. We performed exploratory simulations with η in the range between 1 and 0.5 and found no qualitative difference with the phase behaviour reported below.
We proceed by applying a rigorous coarse-graining procedure developed in 34 to Eq. (3) to derive a system of mean-field equations for the density of filaments ρ, their mean orientation p i , and the nematic tensorial field Q i j . These variables are defined as the first three moments of P(r, φ ): where i, j = {x, y}, and we introduced dimensionless units 48 . The resulting equations contain a very large number of terms, as is often the case with kinetic theories, and their explicit form is given in 48 . To study the dynamics predicted by this approach, we perform numerical simulations of the hydrodynamic equations and discuss representative results below (see Fig. 2.) Simulations are initialised from an isotropic uniform MT suspension with overall density ρ 0 and a small amount of noise. Without loss of generality, we set τ 0 = 1/2, and vary Ξ and ρ 0 .
A linear stability analysis 48 shows that the uniform isotropic state is linearly unstable towards the emergence of a globallyordered nematic state, when ρ 0 > ρ cr = 6π/(1 + Ξ(1 − τ 0 )). Additionally, for ρ cr < ρ 0 < ρ N , this nematic state is itself unstable. Simulations demonstrate that the latter instability leads to co-existence between high-density, nematically-ordered elongated domains and a low-density isotropic background ( Fig. 2a and Suppl. Movie 1). The outcome of this phase separation at late times depends on the value of the anisotropy parameter, Ξ. For small Ξ, domains coarsen to leave a single static band, whose size scales with the system size (Fig. 2a). Inside the band, MTs are ordered nematically, with residual polar order confined at the interface with the isotropic phase. For large enough Ξ, we instead observe an ever-evolving pattern ( Fig. 2b and Suppl. Movie 2), superficially reminiscent of "active turbulence" 49 in wet active gels. To characterise the properties of this spatiotemporal pattern, which we call dry active turbulence, we plot the time evolution of the domain size, computed via the first moment of the structure factor 48 , and its Fourier transform (Figs. 2f and h respectively). It is apparent that there is a selected lengthscale in the isotropic case, while the dynamics in the anisotropic case appear to be chaotic (as the Fourier transform in Fig. 2h contains all frequencies). Our findings are summarised in the phase diagram in Figure 2d.
The kinetic pathway associated with dry active turbulence becomes apparent in simulations with smaller domains (Fig. 2c,  Suppl. Movie 3). These shows that the self-assembled nematically ordered MT bands undergo a cyclic process where they stretch perpendicular to their long direction, rotate, stretch and split, to reform later on. This process is quasi-periodic in smaller system, but appears to be chaotic in larger ones.
To identify the fundamental mechanism leading to pattern formation in our system, we now search for a minimal model. We define the latter as a set of simple equations, which simultaneously satisfies two conditions. First, it needs to have qualitatively similar dynamics as the full model (Figs. 2a and b): it should retain both a transition between a uniform and a phase separated nematic state, as well as a regime with chaotic dynamics; in small domains, it should exhibit features similar to Figure 2c. Second, we require that the location of the phase boundaries in the mini-mal and full models (Figs.2d and e), is quantitatively similar. As a first step, we exploit the observation that polar order plays a minor role (Fig. 2c), and adiabatically eliminate p i in favour of ∂ i ρ and ∂ j Q i j , keeping only the lowest order terms in spatial gradients (as in a hydrodynamic expansion 50 ). Then, we systematically switch off each term individually in the resulting equations, and compute the phase diagram; the term is only reinstated if its exclusion leads to a substantial change in the phase boundary location.
This procedure yields the minimal model defined in Eqs. (1) and (2), and its phase diagram is presented in Figure 2e. The eight constantsµ, χ, λ , α, κ, ζ , β 1 , β 2 -are not arbitrary, and are explicitly derived in terms of the microscopic parameters ρ 0 , Ξ, τ 0 and η (see SI for details 48 ). Within this set, ζ is the only parameter that can change sign -the others are always positive.
We now discuss the physical meaning of each term in Eqs. (1) and (2). First, µ and λ determine the non-equilibrium chemical potential of our mixture: their main role is to set the values of the densities in the isotropic and nematic phases. Next, α is a non-equilibrium Landau coefficient setting the magnitude of order in the bulk (together with the term 4 (ρ/ρ cr − 1) Q i j ), while κ is the nematic elastic constant. Similar terms are also present in a purely passive Model C 51 describing, for instance, phase separation in passive liquid crystals † . The key qualitative ingredients that produce chaotic behaviour in our model are the terms proportional to χ, ζ , β 1 and β 2 . Among them, χ is an "extensile flux", whose role is similar to that of an extensile stress in active gels 1,7,52 . This term enhances diffusion along the direction of the local nematic order (i.e., the eigenvector of Q i j corresponding to its positive eigenvalue), and decreases it along the perpendicular direction. Second, ζ creates an effective torque at the interface, as the associated term depends on density gradients, which are largest at the interface. When ζ is positive (negative), it tends to orient MTs parallel (perpendicular) to an isotropic-nematic interface. Finally, β 1 and β 2 create modulation of the nematic ordering (i.e., the positive eigenvalue of Q i j ). These terms promote activity-induced disorder, and act similarly to a negative elastic constant 10 . Additionally, they contribute to the interfacial torque at the boundary of a nematic band, where Q kl Q kl drops sharply to zero, following the density field. While similar terms could in principle be derived from an effective free energy, a thermodynamic framework would lead to additional constraints between the values of β 1 and β 2 .
The minimal model is now simple enough for us to dissect the mechanisms underlying pattern formation. The kinetic pathway leading to non-equilibrium phase separation proceeds as follows. Starting from a uniform disordered solution with ρ > ρ cr , MTs rapidly acquire orientational order, through the Landau coupling in Eq. (2). At this point, the extensile active flux, arising from MT sliding, enhances diffusion along the nematic direction, and hinders it perpendicularly. When this effect is strong enough, the perpendicular diffusion becomes effectively negative, causing MT † Here we generically use the term "model C" to describe a set of equations coupling the dynamics of a conserved field to that of a non-conserved one. bundling and the formation of one or more nematically ordered high-density bands (see Fig. 2 and Suppl. Movies 1, 4). Notably, although the phase separation is driven by a non-equilibrium phenomenon (MM activity), the kinetic growth laws resemble canonical Model C phase separation in passive mixtures of liquid crystalline and isotropic fluids 48,51,53 .
Second, when Ξ is sufficiently large, the β 1,2 terms dominate over both the restoring elastic constant κ and the ζ term: the associated torque rotates the MTs at the nematic-isotropic interface, so that they tend to orient perpendicular to the band border. This interfacial alignment conflicts with the direction of the nematic order within the bulk of the band; it couples to the extensile flux to yield locally synchronous rotation (and stretching) of nematic bands as observed in our simulations. This cycle repeats, creating a never-settling pattern, as seen in our simulations in the dry active turbulent regime (Figs. 2 and 4a, and Suppl. Movies 3 and 5). As the sense of the emerging band rotation (clockwise or anticlockwise) is selected by spontaneous symmetry breaking, it may be different in different regions of our simulation domain, yielding a chaotic pattern (Fig. 3b, Suppl. Movies 2 and 6). Measuring the time evolution of the domain size in this regime yields statistically the same results as for the full model (Figs. 2g and h).
There is also a second mechanism that can destabilise a homogeneous nematic state, again dependent on β 1,2 . A linear stability analysis starting from the uniform nematic phase 48 shows that when these terms are large enough, they trigger the development of a modulation in Q i j -in the direction parallel to that of the nematic order, for β 1,2 > 0. This instability is independent of density fluctuations and ultimately fragments the nematic phase into infinitely small microdomains. This pathway to chaos is related to that identified in 9,10 for dry active matter with near-uniform density. However, in our model this instability is only relevant for ρ 0 ρ cr , and for lower ρ 0 is superseded by the turbulent phase separation dynamics discussed above. While our minimal model is the result of a systematic coarsegraining, we can view Eqs. (1) and (2), more generally, as a phenomenological model that contains the lowest terms of the correct tensorial nature 54 . Upon coarse-graining, a microscopic model within the same universality class as the one studied here would, therefore, provide the expressions for the parameters in Eqs. (1) and (2), but would not generate extra terms. Indeed, setting β 1,2 = λ = µ = 0 shows that our equations, in this limit, reduce to the minimal model for flocking of self-propelled particles with nematic order 8,[14][15][16] . We, therefore, propose Eqs. (1) and (2) as a unifying model for dry active systems with nematic order. Recently, similar arguments were used to propose active versions of Models B and H 4,5 in Hohenberg-Halperin classification 51 . We follow this analogy and refer to Eqs. (1) and (2) as active Model C. This model is in a different universality class with respect to active gel theory 1 , which exhibits instabilities in an active nematic fluid with constant density, whereas in our case patterns are always associated with a non-equilibrium phase separation. We want to stress that while previous work reported types of chaotic behaviour similar to the limiting cases of our model, either based on hydrodynamic 8,14 or kinetic theories 15,55 , active model C unifies all this into a general universality class.
Analysis of active Model C with phenomenological coefficients re-enforces our physical interpretation of the instability modes. First, nematic-isotropic phase separation also occurs with ζ = β 1,2 = 0, confirming that this phenomenon relies solely on a nonzero extensile flux, χ = 0 (Fig. 3c, Suppl. Movie 7). Second, setting χ = 0 whilst retaining β 1,2 and ζ only leads to a uniform nematic phase, confirming that χ is necessary for any patterning (Suppl. Movie 8). Third, switching off only ζ leads to chaotic dynamics for a much wider parameter range, including Ξ = 0 (Fig. 3d, Suppl. Movie 9), as now β 1,2 only need to compete with the elastic constant κ. Fourth, switching off only β 1,2 whilst retaining ζ > 0 does not lead to chaotic dynamics as in Fig. 2b and 3a,b, as there is no competition between the orientational order in the bulk and at the interface (see 48 ) . This case however, yields another interesting instability associated with interfacial undulations and an elastic bend deformation in the nematic order (Fig. 3e, Suppl. Movie 10). The ensuing patterns may also be chaotic for sufficiently large ζ (Suppl. Movies 11,12), and are similar to the structures seen experimentally in microtubulekinesin mixtures 25 .
For various values of its parameters, active Model C serves as a catalogue of patterns in dry active systems. As mentioned above, a sub-set of terms in Eqs.(1) and (2) was obtained in models of flocking of self-propelled particles with nematic order 8,14-16 . Within those models, rigorous coarse-graining shows that χ and ζ are both positive, and, accordingly, the generic outcome found by numerical simulations 8,14-16 is non-equilibrium phase separation and chaos through band undulations (as in Fig. 3e). Based on our phenomenological model and interpretation, we also expect dry active turbulence with contractile active flux, χ < 0, and interfacial torques favouring parallel alignment at the interface, as would occur when ζ , or β 1,2 are positive. This scenario may be relevant for pattern formation in MT-MM mixtures where the underlying microscopic collision rules differ from those in Figure 1. Further work is required to identify the criteria for a microscopic model to belong to the universality class of our active Model C.
Discussions with Hugues Chaté are kindly acknowledged. AG acknowledges funding from the Biotechnology and Biological Sciences Research Council of UK (BB/P01190X, BB/P006507). DM acknowledges support from ERC CoG 648050 (THREEDCELL-PHYSICS).

Conflicts of interest
There are no conflicts to declare.