Geometry and topology of turbulence in active nematics

The problem of low Reynolds number turbulence in active nematic fluids is theoretically addressed. Using numerical simulations I demonstrate that an incompressible turbulent flow, in two-dimensional active nematics, consists of an ensemble of vortices whose areas are exponentially distributed within a range of scales. Building on this evidence, I construct a mean-field theory of active turbulence by which several measurable quantities, including the spectral densities and the correlation functions, can be analytically calculated. Because of the profound connection between the flow geometry and the topological properties of the nematic director, the theory sheds light on the mechanisms leading to the proliferation of topological defects in active nematics and provides a number of testable predictions. A hypothesis, inspired by Onsager's statistical hydrodynamics, is finally introduced to account for the equilibrium probability distribution of the vortex sizes.

The paradigm of "active matter" [1][2][3] has had notable successes over the past decade in describing selforganization in a surprisingly broad class of biological and bio-inspired systems: from flocks of starlings [4,5] to robots [6], down to bacterial colonies [7][8][9][10][11], motile colloids [12,13] and the cell cytoskeleton [14][15][16]. Active systems are generic non-equilibrium assemblies of anisotropic components that are able to convert stored or ambient energy into motion. Due to the interplay between internal activity and the interactions between the constituents, these systems exhibit a spectacular variety of collective behaviors which are entirely self-driven and do not require a central control mechanism.
A particularly interesting manifestation of collective behavior in active systems is the emergence of spatiotemporal chaos. In active bio-fluids, such as bacterial suspensions or cytoskeletal mixtures, the chaotic dynamics takes place through the formation of structures, such as jets or swirls, reminiscent of turbulence in Newtonian fluids, in spite of the undisputed predominance of dissipation over inertia at the microscopic scale. Examples of low Reynolds number turbulence in active fluids have been first reported for the case of bacterial suspensions [7][8][9][10][11], where this is believed to have an important impact on nutrient mixing and molecular transport at the microbiological scale.
Recently, a series of remarkable experiments on actomyosin motility assays [17] and suspensions of microtubles bundless and kinesin [18,19] (see Fig. 1A), have unveiled a profound link between the topological structure of the orientationally ordered constituents and the flow dynamics, suggesting that active turbulence could be mediated by unbound pairs of topological defects. According to this picture, the strong distortion associated with a defect, fueled by the active stresses, determines a local shear flow which in turn drives the unbinding of more defect pairs. Similar patterns have been observed in "living liquid crystals" obtained from the combination of swimming bacteria and lyotropic liquid crystals [20].
Whether these examples of active turbulence are different realizations of the same universal mechanism or substantially different forms of spatio-temporal chaos, represents Flow velocity (black streamlines) and vorticity (background color). (C) Schlieren texture constructed from the director field n. The red and blue dots mark respectively the +1/2 and −1/2 disclinations. (D) Clockwise rotating (blue) and counterclockwise rotating (red) vortices, detected by measuring the Okubo-Weiss field as described in the text. a profound and yet unsolved problem.
The current efforts toward understanding active turbulence rely on the use of continuum models, such as the Toner-Tu or Swift-Hohenberg model [10,11] or the equations of active nematodynamics [21][22][23][24][25]. Both these approaches have been shown to be able to account for the occurrence of self-sustained low Reynolds number turbu-lence such as that observed in the experiments on bacteria and cytoskeletal fluids, although a systematic comparison between theory and experiments is still in its infancy. The recent numerical work by Thampi et al. [21][22][23], in particular, has provided a convincing demonstration of the correlation between defects dynamics and turbulence in active nematics. The interplay between defects and turbulence has been further investigated in Ref. [24,25] and, following a different approach, in Ref. [26]. Agentbased simulations have been also recently employed to highlight the interplay between defects and dynamics in granular active nematics [27]. The overlap between these "dry" systems and active liquid crystals remains, however, unclear.
In this article I report an exhaustive numerical study of turbulence in active nematics. As a starting point, I demonstrate that, as for inertial turbulence, low Reynolds number turbulence in active nematics is in fact a multiscale phenomenon characterized by the formation of vortices spanning a range of length scales. Within this active range, the areas of the vortices are exponentially distributed, while their vorticity is approximately constant. Building on these observations, I then formulate a mean field theory of turbulence in active nematics that allows the analytical calculation of several measurable quantities, including the mean kinetic energy and enstrophy, their corresponding spectral densities and the velocity and vorticity correlation functions. The connection between the topological structure of the nematic phase and the geometry of the flow is then elucidated through a quantitative description of the defect statistics.

A. Active nematdynamics
Let us consider an incompressible uniaxial active nematic liquid crystal in two spatial dimensions. The twodimensional setting is appropriate to describe the experiments such as that by Sanchez et al. [18], where the microtubule bundles are confined to a water-oil interface forming a dense active nematic monolayer, but also of considerable interest in its own right. Let then ρ and v be the density and velocity of an incompressible nematic fluid. Incompressibility requires ∇·v = 0. Nematic order is described by alignment tensor Q ij = S(n i n j − δ ij /2), with n the director and 0 ≤ S ≤ 1 the nematic order parameter. The tensor Q ij is by construction traceless and symmetric and has only two independent components in two dimensions. The hydrodynamic equations of an active nematic can be constructed from phenomenological arguments [14,28,29], or derived from microscopic mod-els [30,31] in the form: Here D/Dt = ∂ t +v·∇ indicates the material derivative, p is the pressure, η the shear viscosity, λ the flow alignment parameter and γ the rotational viscosity [32]. In Eq.
, are the strain rate and vorticity tensors corresponding to the symmetric and antisymmetric parts of the velocity gradient, while H ij = −δF LdG /δQ ij is the so called molecular tensor, governing the relaxational dynamics of the nematic phase and obtained from the two-dimensional Landau-de Gennes free energy [32]: with K and C material constants. Finally, the stress tensor σ ij = σ e ij +σ a ij is the sum of the elastic stress σ e ij = −λH ij + Q ik H kj − H ik Q kj , due to the entropic elasticity of the nematic phase, and an active contribution σ a ij = αQ ij describing the contractile (α > 0) and extensile (α < 0) stresses exerted by the active particles in the direction of the director field.
Eqs. (1) have been numerically integrated in a square domain of size L with periodic boundary conditions. To render the equations dimensionless, all the variables have been normalized by the typical scales associated with the viscous flow. Distances are then scaled by the system size L, time by the time scale of viscous dissipation τ = ρL 2 /η and stress by the viscous stress scale Σ = η/τ . Finally, low Reynolds number is imposed by setting Dv i /Dt = ∂ t v i in Eq. (1a). The integration is performed by finite differences on a square grid of 256×256 points via a fourth order Runge-Kutta method. To make contact with the recent and ongoing experiments on microtubules suspensions [18,19], I restrict the discussion to the case of extensile systems (α < 0). The contractile case was found to be nearly identical and is briefly described in Appendix D. Unless stated otherwise the parameter values used in the numerical simulations are λ = 0.1, K = 1, γ = 10 and C = 4 × 10 4 , in the previously described units.

B. The active range
Eqs. (1) contain two important length scales, in addition to the system size L. These are the coherence length of the nematic phase n = K/C and the active length scale a = K/|α|. The former determines how quickly the nematic order parameter drops in the neighborhood of a topological defect and can be taken as a measure of the defect core radius. The quantity a , on the other hand, is the length scale over which active and passive stresses balance, leading to spontaneous elastic distortion and hydrodynamic flow [22,23,25]. As a consequence, a quiescent uniformly oriented configuration becomes unstable to a laminar flowing state once a ∼ L [15,28,29,33]. As a becomes lower than the system size L, the laminar flow is eventually replaced by a turbulent flow. Depending on the values of the various parameters in Eqs. (1), the onset of turbulence can be characterized by the formation of "walls", narrow regions where the director is highly distorted, and their breakup into pairs of ±1/2 disclinations [21,25]. The number of unbound defects increases with activity until saturation when a ≈ n , as the reduction of the nematic order parameter due to the defects compensates the activity increase [25]. Here we overlook on the problem of the onset and focus on the regime where turbulence if fully developed, but still far from saturation: thus n a L. Fig. 1B shows the typical structure of the turbulent flow arising from Eqs. (1) for large (negative) values of the active stress α. The velocity field appears to be decomposed in vortices of various sizes and shapes, while the director field is highly distorted by the presence of several ±1/2 disclination pairs (Fig. 1C).
In order to demonstrate the multiscale structure of the flow, I have measured the distribution of the vortex area. Calling a the area of a vortex, this can be described by a density function n(a), such that dN = da n(a) is the total number of vortices of area in between a and a + da. The area of a vortex can be measured from the numerical data by introducing the so called Okubo-Weiss field [34,35] with ψ the stream function, such that v x = ∂ y ψ and v y = −∂ x ψ. The quantity Q is related with the Lyapunov exponent of a tracer particle advected by the flow: where Q > 0, the distance between two initially close particles will diverge exponentially in time, while for Q < 0 the trajectories will remain close. Thus coherent regions in the flow are defined as regions in which Q < 0 [34]. The Okubo-Weiss field has been widely used for analyzing atmospheric and oceanic circulations and realized as an important quantity to characterize two-dimensional flows [36,37]. To identify the vortex cores, on the other hand, one can calculate the angle the velocity field rotates in one loop around each cell of the computational grid [38]. If the cell contains the core of a vortex, this angle is equal to 2π, regardless of whether the vortex is left or right-handed. The combination of these two criteria allows to formulate the following vortex detection algorithm: 1) from the velocity field v the cores of the vortices are initially located; 2) the area of a vortex is then defined as the area of the region surrounding a vortex core where Q < 0. Fig. 1D shows the vortices detected by this method. (1). The data show a prominent exponential distribution of the form: where a min and a max ∼ L 2 are respectively the minimal and maximal area of an active vortex and a * a suitable scale parameter. By analogy with the inertial range in classic turbulence, hereafter, we will refer to this range as the active range. Here by active vortex we mean a vortex resulting directly from mechanical work performed by the active stresses. Beside active vortices, other vortices might form due to the strong shear in the space between active vortices. These secondary vortices are expected to lie outside the active range, thus where a < a min . The quantities N = amax amin da n(a) and Z = amax amin da exp(−a/a * ) in Eq. (3) represent respectively the total number of active vortices and a normalization constant. In the Conclusions I will speculate about the physical origin of this exponential distribution. The vortices mean vorticity Fig. 2 as a function of a. Unlike n(a), ω v remains roughly constant across the scales and shows some dependence only for large activity values, where the mean vortex size is substantially smaller than the size of the system (see the blue line in the inset of Fig. 2).
As activity is increased, the vortices become smaller and faster as indicated by the dependence of a min , a * and ω v on α. Being a the only length scale associated with activity, intuitively one could expect that a min ≈ a * ∼ 2 a . Analogously, the balance of active and viscous stresses over the scale of a vortex, suggests that ω v ∼ α/η. These expectations are confirmed from the numerical data shown in Fig. 3. The multiscale organization and the exponential distribution of the vortex areas have striking consequences on the overall statistical properties of the flow. From a gross application of the central limit theorem we could expect the velocity components to be Gaussianly distributed. The numerical data shown in Fig. 4A support this expectation. As in classic high Reynolds number turbulence, on the other hand, vorticity and in general any function of the velocity gradients do not obey the Gaussian distribution due to the spatial correlation introduced by the derivatives [39]. In this particular case, the vorticity probability density function (PDF) exhibits a visible deviation from Gaussianity along the tails (see Fig. 4C). Fig. 4C-D show the normalized velocity and vorticity correlation functions: C vv (r) = v(0) · v(r) / |v 2 (0)| and C ωω (r) = ω(0)ω(r) / ω 2 (0) , where the angular brackets · indicate an average over space and time. These quantities have played a central role in the study of active turbulence starting from the experimental work by Sanchez et al. [18]. In this work it was argued that, after rescaling by the mean squared value, the correlation functions no longer depend on activity suggesting that the underling geometrical structure of the flow is due to a passive mechanism, while activity only controls the flow average speed. This scenario found support in the numerical work of Thampi et al. [21], who also provided an elegant interpretation based on the creation and annhilation dynamics of topological defects. The numerical data shown in Fig. 4, combined with what reported in the previous section, demonstrate that the geometrical structure of the flow does in fact depend on activity through the active length scale a . Such a dependence, which is clearly marked by the intersection of the curves in Fig. 4C-D, is however subtle and could have been missed before due to the limited activity range explored. An analytical approximation of the correlation functions C vv (r) and C ωω (r) is given in Appendix C.
To gain further insight about how activity affects the flow, I have measured the average kinetic energy v 2 /2 and enstrophy ω 2 /2 per unit area for varying α values ( Fig. 5A-B). These exhibit respectively a clear linear and quadratic dependence on activity. This scaling properties can be straightforwardly understood from the geometrical picture previously described. As the vorticity is decomposed in a discrete number of vortices having ω v ≈ α/η, the total enstrophy can be expressed as: where N is the total number of vortices and ( · ) = da n(a)( · )/N indicates the vortex ensemble average. Thus: where we have used the fact that a ≈ L 2 /N . Analogously, the total kinetic energy of a single vortex is given by: with the approximation becoming an equality in the case of a circular vortex. Averaging over the vortex ensemble thus gives E tot = 1/(16π) ω 2 v N a 2 , from which: While changing the resolution of the vortex ensemble, activity does not affect the spectral structure of the flow. Fig. 5C-D, show the enstrophy Ω(k) and energy E(k) = Ω(k)/k 2 spectra [39] for various activity values. Analogously to what is observed in bacterial turbulence [10], the spectra are non-monotonic with a peak around k a = 2π/ a dividing the growing regime at small k−values from the decay regime at large k−values. In the latter regime, the data show a clear power-law decay with Ω(k) ∼ k −2 and E(k) ∼ k −4 . In the next section we will illustrate the origin of these exponents in a meanfield framework.

D. Mean-field theory
The spectral structure of turbulence in twodimensional active nematics as well as short-scale velocity and vorticity correlation can be satisfactorily described within a mean-field approximation. This approach was introduced by Benzi et al. [34,40] to account for the emergence of self-similar coherent structures in two-dimensional decaying turbulence and can be extended to the non self-similiar case discussed here. Let us consider a two-dimensional flow whose vorticity field can be decomposed in a discrete number of vortices of radius R i and vorticity ω i (r) = ω v,i f (r/R i ), with r the distance from the vortex center and ω v,i a constant.
where r i is the position of the i−th vortex center. The power-spectrum of the function ω(r) can then be expressed as: (6) whereF (kR) = 1/(2π) ∞ 0 dξ ξf (ξ)J 0 (κR ξ), with J 0 a Bessel function of the first kind, is a dimensionless vortex structure factor (see Appendix A). Now, if we neglect the spatial correlation between the vortices only the diagonal terms in the sum survive upon averaging. Then: where we have replaced the summation with an integral over the vortex population. Finally, the enstrophy spectral density, can be calculated from the vorticity power spectrum as: Ω(k) = 4π 3 k |ω(k)| 2 (see Appendix B). Now, the distribution function n(R) can be obtained from Eq. (3) upon setting a = πR 2 , so that n(R) = |da/dR|n(a). Note that assuming the vortices to have circular shape is not generally correct as it is not guaranteed that the ansatz used to parametrize the vorticity field would hold in general. Nonetheless, based on the existence of a single characteristic length scale a , one could hope that both approximations would affect the accuracy of the calculation only through irrelevant prefactors. Next, using the fact that ω v does not depend on R and replacing n(R) into Eq. (7) yields, after some algebraic manipulations: where we have set R * = a * /π ∼ a and κ = kR * . Now, consistently with the previous assumption about the shape of a vortex, we can choose f (r/R) = 1 for r/R ≤ 1 and f (r/R) = 0 otherwise, then the structure factor can be easily calculated in the formF (ξ) = J 1 (ξ)/(2πξ) (see Appendix A). Using this in Eq. (8) and extending for simplicity the integration to the whole positive real axis, yields: where I 0 and I 1 are modified Bessel function of the first kind [41] and C = π 2 ω 2 v N R * 5 /(2Z) is a quantity independent on κ. The spectral structure of the turbulent flow is encoded in the asymptotic behavior of the function in Eq. (9). For κ 1, Ω(κ) ∼ κ −2 (see Appendix B) in perfect agreement with the numerical data. The energy spectral density can be calculated straightforwardly from Ω(κ), this yields E(κ) ∼ κ −4 . For κ ≈ 0, on the other hand, Eq. (9) yields Ω(κ) ∼ κ and E(κ) ∼ κ −1 . These predictions are very difficult to compare with the numerical data as they refer to the narrow range of the spectrum (i.e. k min < k < 10k min , with k min = 2π/L) preceding the crossover region.
From the spectra Ω(k) and E(k), the correlation functions C vv (r) and C ωω (r) can be easily determined (see Appendix C). Due to the mean-field approximation, however, the accuracy of this calculation is limited to the range 0 < r < R * , where the spatial correlation between vortices is negligible. The flow is obtained from an analytical solution of the incompressible Stokes equation in the presence of a body force f ± = ∇ · (αQ ± ) and Q ± the nematic tensor associated with a ±1/2 disclination [25].

E. Topological structure of active turbulence
As I mentioned in the introduction, the geometry of the flow field is strictly connected with the topological structure of the nematic phase [18,[21][22][23][24]. As it has been stressed in [25], the configuration of the nematic director in the neighborhood of a ±1/2 disclination determines the local vortex structure (Fig. 6). Topological defects serve then as a template for the turbulent flow, which in turn advects the defects themselves leading to chaotic mixing. To provide a quantitative description of the defect chaotic dynamics I have measured the meansquared displacement (MSD) of ±1/2 disclinations as a function of time (Fig. 7A). For both positively and negatively charged defects this shows a substantially diffusive behavior, with a slight super-diffusive trend in the short time dynamics of +1/2 disclinations, due to the self-propulsion provided by the self-induced dipolar flow ( Fig. 6A and Ref. [25]).
The total number of topological defects N d is evidently proportional to the number of vortices. Thus N d ∼ L 2 / a ∼ α, consistently with what found numerically (see Fig. 7B). As already mentioned, this linear growth in the defect population tends to saturate when the active length scale a approaches the defects core radius, proportional to n (not shown here) [25]. Analogously, the defect mean free path (Fig. 7C) is Fig. 7D, shows the defect creating and annihilation rates ν c and ν a versus activity (in units of the viscous time scale τ ). For large activity values these exhibit a quadratic dependence on α: ν c ≈ ν a ∼ α 2 . This behavior can be understood by noticing that topological defects moves predominantly along the edge of the vortices at ap- proximatively constant angular velocity ω v ≈ α/η. During this circulation, they might approach an oppositely charged defect an annihilate. The elastic, Coulomb-like, attraction between oppositely charged defects takes over only when these have become very close to each other, thus the typical time scale of annihilation, t a , is predominantly dictated by the active circulation: t a ∼ 1/ω v . From this we might expect that ν a ≈ N d /t a ∼ α 2 . Once turbulence reaches a steady state, the defects creation and annihilation balance, hence ν c ≈ ν a .

II. DISCUSSION AND CONCLUSIONS
In this article I have reported a thorough numerical and analytical investigation of low Reynolds number turbulence in two-dimensional active nematics. Spectacular experimental realizations of this system are found in cytoskeletal fluids of microtubules and kinesin at the wateroil interface [18] or incapsulated in a lipid vesicle [19]. For large enough activity values (corresponding to high concentrations of motors or adenosine-triphosaphate in cytoskeletal fluids), these systems are known to develop a chaotic spontaneous flow reminiscent of turbulence in viscous fluids (see Fig. 1A).
Here I have demonstrated that, as for inertial turbulence, low Reynolds number turbulence in active fluids is in fact a multiscale phenomenon characterized by the appearance of vortices spanning a range of length scales. Within this active range the vortex areas follow the exponential distribution. This peculiar geometrical structure of the flow leaves a strong signature on all the relevant physical observables. The mean kinetic energy, for instance, scales linearly with activity (and not quadratically as one could have naively expected from a comparison with the laminar case, where v ∼ αL/η) because the vortices becomes smaller as activity is increased. Furthermore, the enstrophy and energy spectra scale as k −2 and k −4 rispectively, thus in net contrast with twodimensional inertial turbulence [39]. The statistics of the vortices, finally, completely determines that of the defects (and vice-versa) making possible the formulation of various scaling relations amenable to experimental scrutiny.
While some questions have been answered in this work, others remain open. How energy and enstrophy flow across the scales? In two-dimensional inertial turbulence it is well known that enstrophy flows toward the small scale where it is eventually dissipated, while energy flows toward the large scale where it is either dissipated by frictional interactions with the wall or condensed in large coherent structures [39]. In complex fluids, on the other hand, kinetic energy can be converted into elastic energy and be dissipated or stored via mechanisms that do not require cascading. While the existence of an active range does imply that of energy and enstrophy flux across scales, the organization of such a flux remains unknown.
Another important question, which I deliberately saved for the end, concerns the origin of the exponential distribution of the vortex areas. Perhaps the most natural explanation appeals to the interpretation of the vortex population as an equilibrium ensemble, subject to the laws of statistical mechanics. This reasoning is not new in two-dimensional turbulence, but goes back to the pioneering works of Onsager [42] and of Joyce and Montgomery [43,44] (see also Refs. [45,46] for a review). To clarify this concept let us consider a system of N active vortices having the same absolute vorticity and let n i be the number of vortices of area a i , so that i n i = N . A microscopic configuration is then characterized by a set of occupancy numbers {n i } ∞ i=1 and, as the vortices are indistinguishable, there are W = N !/ i n i ! different ways to realize the same microscopic configuration. In the limit of large N , corresponding to fully developed active turbulence, W ∼ e S , where S = − da n(a) log n(a) is an analog of the microcanonical entropy for the vortex ensemble. Since the vortices all have the same vorticity, a macroscopic state can be arguably identified by the their total number N = da n(a) and area A = da n(a) a, with a = A/N . The most probable (N, A)−macrostate is that maximizing the entropy S for fixed N and A, hence n(a) = N/a exp(−a/a). Finally, setting a ∼ 2 a yields an expression equivalent to Eq. (3).
This hypothesis is inevitably naive and yet incredibly fascinating in suggesting an unexpected connection between the simplest and the most complex forms of matter. The dimensionless vortex structure factorF (kR), introduced in the mean-field calculation, is defined from the Fourier transform of the vorticity profile function f (r/R), where r is the distance from the vortex center. Thus: where J 0 is a zeroth-order Bessel function of the first kind. Now, the simplest approximation of the function f is evidently: representing a circular vortex of radius R. Replacing this into (A1) yields: The dimensionless functionF (kR) is then defined from f (k) = R 2F (kR), hence: Appendix B: Asymptotic behavior of the spectral densities The expression of the enstrophy spectrum, as obtained within the mean-field approximation, is given by: Now, small κ limit can be easily determined by considering that, for κ 1, I 0 (κ 2 /2) ≈ exp(−κ 2 /2) ≈ 1 and I 1 (κ 2 /2) ≈ 0. Therefore: To calculate the large κ limit we can use the following asymptotic expansion of the modified Bessel function: From this we obtain: (3), can be obtained by setting R * = a * /π, this gives R * = 0.054 rmax slightly larger than the value obtained from a fit of the correlation function. This slight discrepancy is presumably due to the inevitable systematic error in the calculation of the vortex area from the Okubo-Weiss field as well as the circular approximation of the vortex shape.
The exponential term exactly cancels that in Eq. (B1) resulting in a simple power-law behavior: The asymptotic behavior of the energy spectrum follows directly from this by virtue of the relation Ω(k) = k 2 E(k).
of Eqs. (1) and their mean-field approximations given in Eqs. (C3) and (C6). For small distances, the agreement is remarkable. This, however, breaks down for r R * , where the spatial correlation between neighboring vortices (which is neglected in the mean-field framework) becomes crucial. For instance, C ωω (r) becomes negative when r is larger than the average vortex diameter, due to the fact that a given central vortex is surrounded by vortices of opposite vorticity (see the red dots in Fig. 8). This feature is clearly absent in the mean field calculation and the resulting correlation function decays without sign changes (solid black line in Fig. 8).