## Abstract

We describe a numerical investigation of a continuum model of an active nematic, concentrating on the regime of active turbulence. Results are presented for the effect of three parameters, activity, elastic constant and rotational diffusion constant, on the order parameter and flow fields. Defects and distortions in the director field act as sources of vorticity, and thus vorticity is strongly correlated to the director field. In particular, the characteristic length of decay of vorticity and order parameter correlations is controlled by the defect density. By contrast, the decay of velocity correlations is determined by a balance between activity and dissipation. We highlight the role of microscopic flow generation mechanisms in determining the flow patterns and characteristic scales of active turbulence and contrast the behaviour of extensile and contractile active nematics.

## 1. Introduction

Active systems produce their own energy and thus operate out of thermodynamic equilibrium. Examples include suspensions of microtubules and molecular motors, cellular and bacterial suspensions, vibrating granular rods, flocks of birds and schools of fish [1–6]. In all these systems, the motion is generated at the level of individual constituents but energetic, hydrodynamic or intelligent interactions then result in a collective dynamical behaviour. In particular, the spontaneous motion generated by a wide variety of dense active systems is highly chaotic, consisting of fluid jets and swirls (see figure 1). The flow pattern visually resembles that of high *Re* number turbulence and hence this state is often termed ‘active turbulence’ [7]. It is not clear whether the similarities in the turbulent states observed in active systems across a wide range of length scales are superficial or represent underlying universal properties.

Many active systems have nematic symmetry and exhibit orientational order, and concepts from liquid crystal physics have been used to describe them. The equations of nematohydrodynamics have proved useful in predicting the occurrence of hydrodynamic instabilities [8–12] and in simulating the turbulent flow behaviour [13–15]. There is also increasing evidence that topological defects and their dynamics play a role in active turbulence [13,14,16,17]. Recent experiments using microtubule bundles and kinesin molecular motors have demonstrated the presence of topological defects in active systems [1]. Similar defects have also been observed in experiments on living liquid crystals [18] and vibrating granular rods [19] and in the simulations of self-propelling ellipsoidal particles [20].

In this paper, we report extensive simulations of active nematics to investigate and characterize the flow field and the order parameter field, and the coupling between them. We find that two length scales are relevant. For extensile systems, scaling arguments suggest that the length describing the decay of velocity correlations results from a balance between activity and dissipation. Correlation functions of the director and vorticity fields have very similar forms which correspond to a different length scale, the defect spacing. We report different behaviour for extensile and contractile systems, showing the sensitivity of the active state to details of the coupling between the director and flow fields.

In §2, we summarize the theory of nematohydrodynamics and show that it can be extended to the active case by including a simple, additional, dipolar source term in the stress. Section 3 gives details of the numerical approaches used to solve the equations of motion. In §4, we give a qualitative account of the characteristic features of active turbulence in extensile and contractile systems explaining, in particular, how defects are continually created, and subsequently annihilated, to give a steady state with a finite defect density. Section 5 presents extensive quantitative results showing how the magnitude and correlations of velocity, vorticity and the order parameter field depend on model parameters. We also measure the number of defects and their rate of creation and annihilation. In §6, we show that an analysis of the defect density and dynamics can be used to relate the different scalings found in the numerics. We also list open questions that warrant attention in the future.

## 2. Equations of motion

We first describe the hydrodynamic equations of motion of an active nematic system. This may be a dense suspension of microtubule bundles and molecular motors [1], or a dense suspension of bacteria [7], or a mixture of bacteria and liquid crystal molecules [18]. Active systems with nematic symmetry, and density and momentum as conserved variables, can be modelled using the standard equations of nematohydrodynamics, modified to incorporate an active term which produces flows in response to gradients in orientational order [3,8]. The magnitude and direction of the orientational order are described by a second rank tensor **Q** which is symmetric (*Q*_{ij}=*Q*_{ji}) and traceless (*Q*_{ii}=0). For uniaxial nematics **Q**=(*q*/2)(3**nn**−**I**), where *q* is the largest eigenvalue of **Q** and is a measure of nematic degree of order, **n** is the director field and **I** is the identity tensor.

The nematohydrodynamic equations can be derived using a Poisson bracket formalism or from linear irreversible thermodynamics with phenomenological models assumed for the relation between fields and currents. The resulting equations governing the evolution of the order parameter **Q** and momentum *ρ***u** are [21,22]
2.1
and
2.2
The generalized nonlinear advection term in equation (2.1) is
where *E*_{ij}=(∂_{i}*u*_{j}+∂_{j}*u*_{i})/2 is the strain rate tensor and *Ω*_{ij}=(∂_{j}*u*_{i}−∂_{i}*u*_{j})/2 is the vorticity tensor. These characterize the deformational and the rotational components of the flow field, respectively. We also define the vorticity, ** ω**=∇×

**u**, which measures the rotation of the fluid elements. It is related to the vorticity tensor by

*ϵ*

_{ijk}

*ω*

_{k}=−2

*Ω*

_{ij}, where

**is the Levi-Civita symbol. We will use vorticity extensively to characterize the flow field in this paper.**

*ϵ*The alignment parameter λ controls the degree of coupling between the orientation field and velocity gradients [12,22]. Mathematically, λ determines the character of the objective time derivative of **Q**. The values λ=1 and λ=−1 correspond to an upper and lower convected derivative and λ=0 corresponds to the corotational derivative which is obtained as an average of the lower and upper convected derivatives. On mapping to the Leslie–Ericksen equation which describes the liquid crystal dynamics in terms of the director field **n** only, (3*q*+4)λ/9*q* is the flow alignment parameter [11] which determines whether the liquid crystal molecules (or here, the active constituents) align or tumble in a shear flow.

In linear approximation, the relaxation of **Q** is driven by the molecular field
2.3
defined as the variational derivative of a free energy functional. **H** is also symmetric and traceless. The phenomenological constant of proportionality, which corresponds physically to the rotational diffusion coefficient, is *Γ*. The Helmholtz free energy,
2.4
has a bulk contribution which is a truncated polynomial expansion in invariants of **Q** [21,22] and a term arising from the gradients in **Q**. This free energy expression is general in that it depends on gradients of the scalar order parameter in addition to the gradients in the director orientation. *A*, *B* and *C* are phenomenological coefficients which are functions of concentration and temperature. Here, we use a single elastic constant (*K*) approximation.

In equation (2.2), the total stress generating the hydrodynamics includes:

(i) The viscous stress, 2.5 where

*μ*is the Newtonian viscosity of the suspension.(ii) The orientational elastic stresses of a passive liquid crystal, 2.6 As a result of this term, the orientational order can affect the flow field, an effect often referred to as the ‘back-flow’. The back-flow can modify the director dynamics, for example accelerating the defect annihilation process [23] and affecting the motion of active defects [16]. The term

*P*=*ρT*−(*K*/2)(∂_{k}*Q*_{ij})^{2}in equation (2.6) is the modified bulk pressure.(iii) The active stress, 2.7 is introduced in [8]. This corresponds to a force dipole density, which is the leading-order contribution of activity to the stress. Gradients in

**Q**thus produce a flow field, in addition to that due to the passive stress terms, which is the source of the hydrodynamic instability in active nematics. The strength of this term is prescribed by the activity coefficient*ζ*, which is +ve for extensile and −ve for contractile systems. For example, molecular motors such as myosin, dynein or kinesin generate contractile stresses as they move on actin filaments or microtubules. As shown recently [1,24], depletent molecules cause formation of microtubule bundles in a suspension and kinesin molecular motors cause the filaments to slide past each other generating extensile stresses in the bundles. Microswimmers also behave effectively as force dipoles in a fluid; for example, bacterium*Escherichia coli*generates extensile stress and alga*Chlamydomonas*generates contractile stress.

More details about this model and its application to passive and active systems can be found in [11,21,22,25–31]. The nematic state has been shown to be unstable to flow induced by the activity [8] and the perturbations of the order parameter field that drive the instability have been studied [10,12]. The onset of activity driven flow in a one-dimensional channel was analysed in [9] and the flow patterns generated by the activity were studied numerically [11,31,32]. The rheology of active materials has also been addressed [30,33]. The role of topological defects in active systems has also attracted much attention recently [13,14,16].

## 3. Simulations

The governing equations (2.1) and (2.2) form a coupled system, and we use a hybrid numerical method to solve them [11,30,31]. Instead of directly solving the momentum equation, a lattice Boltzmann method is used to determine the hydrodynamics. For the order parameter evolution equation, we employ a method of lines where all spatial derivatives are calculated using a finite difference approach. The resulting system of ordinary differential equations in time is solved using Euler’s integration method [29,33]. At every time step, the order parameter field obtained from the method of lines is an input to the lattice Boltzmann simulations and then the velocity field obtained from the lattice Boltzmann is used in the next time step of the method of lines for updating **Q**. The computations are performed until a statistical steady state is reached.

For the lattice Boltzmann simulations, we use a *D*3*Q*15 model where the space and time discretizations are chosen as unity. Simulations are performed on a two-dimensional domain of size 400×400 with periodic boundary conditions in all directions. The parameters used are *Γ*=0.34, *A*=0.0, *B*=−0.3, *C*=0.3, *K*=0.02, unless mentioned otherwise. We choose λ=0.7 corresponding to aligning flows [22,12]. All quantities are reported in lattice Boltzmann units and depending upon the system of interest, for example microtubule or bacterial suspensions, appropriate conversion to physical units may be performed [14,27,29].

## 4. Active turbulence in extensile and contractile systems

Here, we give a qualitative description of active turbulence, before presenting and discussing quantitative data for the turbulent state.

Figure 1 shows typical snapshots of fully developed turbulent patterns for a two-dimensional active nematic comparing extensile and contractile stresses and high and low values of the activity. For each set of parameters, the flow field is represented by firstly the streamlines which reveal the circulating patterns and secondly the colour-shaded vorticity field which reveals the direction of rotation. The flow fields are highly disordered in both the extensile and contractile cases, with characteristic swirling patterns. The corresponding director fields are also shown. These also exhibit a high degree of disorder although, for the parameters we use, the corresponding passive liquid crystal is nematic (not shown).

Owing to active stress, flow is generated in response to distortions in the order parameter field. The flow further disturbs the director field thus resulting in active turbulence. As the dominant instability and hence the mechanism for flow generation is bend distortions of the director field for extensile systems and splay distortions for contractile systems [10,12], both the flow and director fields show qualitative differences in their structures.

In the transition to active turbulence, active nematic regions develop *walls*, which are elongated distortions in the director field. Correspondingly, the flow field also exhibits elongated structures of vorticity (figure 1*c*,*e*). We shall call this active turbulence at ‘low activity’. In low activity regime, the walls persist but are continuously advected and deformed due to flow and flow gradients.

By contrast, at ‘high activity’, the walls are strongly deformed by the flow generated by the activity. They are unstable and give rise to the formation of pairs of oppositely charged defects. The defect formation is driven by both elasticity and flow [15]. These topological defects of charge are distinguished by red and blue colouring in the director field in figure 1*b*,*d*,*f*,*h*. As a result of the stronger velocity fields and the defect formation, the vorticity structures in the flow field are also more isotropic (figure 1*a*,*g*) compared with those in the low activity regime. Topological defects in active systems have been observed in a variety of experimental systems and simulation studies [1,18,20].

Details of defect dynamics for an extensile and a contractile system are shown in figure 2. Consider the extensile system first (figure 2*a*). Some of the walls are indicated in the figure. The and defects are highlighted with red and blue colour and their trajectories are also shown. Yellow continuous lines are their past, and magenta dashed lines their future paths. It is apparent that the defects are formed in walls. They preferentially move along the wall releasing the associated bend free energy. During this motion, defects encounter oppositely charged defects and annihilate, thus reinstating nematic regions. All these events are illustrated in figure 2*a*; for example: a creation event is labelled as ‘A’, defects moving along the walls as ‘B’ and a defect annihilation event as ‘C’. A similar snapshot of defect dynamics for a contractile system is shown in figure 2*b*. As the dominant mode of instability for a contractile nematic is splay deformations in the director field, the walls appear different in this case. As can be seen in figure 2*b*, they form the borders of nematic regions of different mean orientations. Again defects are formed in the wall, and preferentially move along the wall. As in the extensile case, they then encounter defects of opposite charge and annihilate re-establishing nematic regions.

As walls tend to be very long, more than one pair of defects may be produced from a single wall, at the same or different times. Defects tend to move along walls and therefore are more likely to encounter defects of opposite charge formed in the same wall and quickly annihilate each other. Some of the defects escape from the walls. The search for a defect of opposite sign is then two dimensional and the defects’ lifetimes are considerably longer.

## 5. Results

### (a) Correlation functions

Vorticity quantifies the local rotation of the fluid elements in a flow field. Many different active systems exhibit continually changing patterns of vorticity [4,7,31,34] (N. S. Rossen, M. H. Jensen & L. B. Oddershede 2013, personal communication), and it was suggested in [7,4] that vorticity might be helpful in characterizing active turbulence. Indeed, vorticity is a fundamental concept in inertial turbulence because regions of high vorticity are characteristic features of turbulent flows. Here, we argue that it is also a key quantity in active turbulence, but for a different reason.

In order to investigate and quantify the structure of the flow and ordering, we plot spatial correlation functions of the velocity, the vorticity and the order parameter field in figure 3*a* and temporal correlation functions of the same quantities in figure 3*b*. To allow an easier comparison, the correlation functions are normalized to unity (i.e. *C*_{vv}(*R*)=〈**v**(*R*,*t*)⋅**v**(0,*t*)〉/〈**v**(0,*t*)^{2}〉, *C*_{ωω}(*R*)=〈** ω**(

*R*,

*t*)⋅

**(0,**

*ω**t*)〉/〈

**(0,**

*ω**t*)

^{2}〉) and the order parameter correlation functions are further scaled to lie between 0 and 1 (i.e. , etc.). There is no scaling of the

*x*-axis in any plot.

For each correlation function, extensile and contractile systems are compared. (Here, and throughout the paper, we use continuous (dashed) lines to depict results for extensile (contractile) systems.) The only difference between the parameters for the extensile and contractile simulations is the sign of the activity coefficient *ζ*. As is clear from the figure, extensile and contractile systems exhibit different characteristic length and time scales of decay for all the quantities measured. This is not surprising considering that the dominant mechanisms for flow generation vary between the two cases.

Note, however, that the director field and vorticity field have a very similar dependence on *R* and *t*, whereas the behaviour of the velocity correlations is quite different. This close association between the order parameter and the vorticity is seen for both the extensile and the contractile simulations. This suggests that vorticity is a fundamental quantity in active turbulence, which is closely related to the microstructure of the system. We shall present more evidence for this assertion in §5*b* when we find that the characteristic lengths associated with the vorticity field and the order parameter field scale in the same fashion with several system parameters.

The connection between vorticity and local order also follows from the governing equations of motion. Generally, active systems such as suspensions of microtubules and molecular motors or bacteria operate at small, often negligible, Reynolds number indicating that these are inertia-less systems and that diffusion of momentum is instantaneous. In this limit, the Navier—Stokes equations considerably simplify to obtain ∂_{j}*Π*_{ij}=0. Taking the curl of this equation and substituting the expressions (2.5)–(2.7) for the viscous, passive and active stresses, gives
5.1
This indicates that gradients in **Q** act as sources of vorticity which diffuse in the system generating vorticity patterns and thus it is reasonable to expect that vorticity and order parameter fields may show a close association in their behaviours.

### (b) Dependence on parameters

In the last section, we showed that the decay of the normalized correlation functions of the director field and the vorticity are governed by the same length scale, but that the velocity correlation function decays over a longer length scale. We now show how the correlation functions depend on three of the most important parameters of the system, the activity coefficient *ζ*, the elastic constant *K* and the rotational diffusion constant *Γ*. The results are collected in figure 4. The three columns of the figure show *C*_{vv}(*R*), *C*_{ωω}(*R*) and *C*_{nn}(*R*). (Note that we plot the director correlation function *C*_{nn}(*R*)=〈**n**(*R*,*t*)⋅**n**(0,*t*)〉 instead of that of the order parameter *C*_{QQ}(*R*). Both behave in the same fashion.)

Figure 4*a*,*b* shows how the correlation functions depend on activity *ζ* for extensile and contractile nematics, respectively. Figure 4*c*,*d* shows results for the dependence on the elastic constant *K* and figure 4*e*,*f* those for the dependence on the diffusion constant *Γ*. In each case, the horizontal distance axis is scaled with the relevant parameter to achieve optimal data collapse. The inset in each panel shows the unscaled data. As the dominant hydrodynamic instability in contractile systems is to splay configurations which yield lower velocities than the bend configurations which occur more frequently in extensile nematics [12], we have been able to cover a larger range of activity and elastic coefficients in the contractile systems.

The optimum scalings with respect to different parameters are based on visual observation of the best collapse and the accuracy of the scaling exponent approximately ±0.05. A nearest rational number is always chosen for brevity of discussion.

The resemblance between the second and third columns of figure 4 is immediately apparent, confirming the similar behaviour of the vorticity and the order parameter correlations. As the latter two quantities scale in the same way, we shall just compare the scaling of the velocity and director correlations (columns I and III) in the following. The characteristic length for the director field is for both extensile and contractile systems. The length scale ℓ_{vel} governing decay of velocity correlations is independent of *ζ* for both extensile and contractile nematics [4,14,35]. It scales as ∼*K* for extensile systems but we have not found any obvious scaling behaviour for the contractile ones. None of the correlation functions depend on *Γ* for either extensile or contractile systems.

Note that in both extensile and contractile systems, activity coefficients 0.001 correspond to the low activity regime. In this regime, no topological defects form yielding elongational flow patterns that are visually different from the high activity patterns (figure 1). Thus, the scalings are expected to be different in the low activity regime as is indeed seen in figure 4. Similarly, at sufficiently large values of *Γ*, the order parameter field relaxes fast and very few defects have time to form in the walls. This again results in elongated structures of both director and flow fields similar to those formed in the low activity regime. Owing to this, only a small range of *Γ* could be scanned in contractile systems.

We also plot the root mean square (RMS) value of the velocity and vorticity as a function of *ζ* and *K* in figure 5, for both extensile and contractile systems. The RMS velocity increases with *ζ* and *K*, with extensile systems showing a stronger dependence than contractile systems. The RMS vorticity increases with *ζ* but is insensitive to changes in *K* for both extensile and contractile systems.

As we saw in §4 that topological defects appear to be important in controlling active turbulence, we next present numerical results for the number of defects and their rate of formation and destruction. As we find only a weak dependence on *Γ* in these and previous measurements, we exclude this parameter from further discussion. In figure 6*a*,*b*, we plot the number of defects as a function of activity and elasticity for extensile and contractile systems finding
5.2
in both cases. At steady state, the rate of creation and the rate of annihilation of defects are equal. We report this rate in figure 6*c*,*d* giving a scaling with activity and elasticity as *r*∼*ζ*^{2}/*K*.

## 6. Discussion

We now discuss the scaling behaviour of the different dynamical variables and suggest a framework in which they can be connected by considering the motion of defects. We will first consider extensile nematics and argue that ℓ_{vel}, the length scale of the velocity field, is controlled by a balance between activity and dissipation, whereas ℓ_{n}, the length scale of the order and the vorticity, is set by the defect spacing.

We consider the rate of energy input into the system by the active stress and assume that this balances the viscous dissipation over a domain of size ℓ_{vel}. Using equation (2.2), we find a scaling for the magnitude of the velocity
6.1
The results presented in figure 6 show that the rate of creation of defects in active turbulence is ∼*ζ*^{2}/*K*. We estimate the rate of annihilation of defect pairs of opposite charge using simple kinetic theory arguments. If *σ* is the scattering cross section for collisional events leading to annihilation and *n* is the number density of defects, then 1/*nσ* is the mean free path and the rate of annihilation ∼*σvn*^{2} where we have made the implicit assumption that the defect and the fluid velocities scale in the same way. At steady state, the rate of generation and the rate of annihilation of defects are equal and thus we obtain
6.2
assuming that *σ*, *μ* and *q* are independent of *K* and *ζ*. The data in figure 6*a*,*b* indicate that . Therefore, we expect
6.3
in agreement with the numerical results in figure 4*a*,*c*. Substituting the expression (6.3) back into equation (6.1) gives *v*∼*ζK*, in agreement with the results for the RMS velocity reported in figure 5*a*.

(In [14], we made the assumptions that the rate of creation of defects scaled as *r*∼*ζ*/*K* and that the number of defects *n*∼1/*K*. The additional numerical results here suggest that a better fit is *r*∼*ζ*^{2}/*K* and . The extra factors of *ζ* drop out to give the same scaling in equation (6.3).)

A second characteristic length in the system is that which controls the decay of both the order parameter field and the vorticity . Our results show that the number of defects ∼*ζ*^{1/2}/*K*∼ℓ^{−2}_{n}, suggesting that the decay of both the director field and the vorticity is controlled by the distance between defects, and hence that the defects are acting as sources of vorticity. The magnitude of the vorticity, however, would be expected to scale as a derivative of the velocity, ∼*v*/ℓ_{vel}∼*ζ*, as confirmed by the data in figure 5*b*. It is interesting to note that the length scale governing the hydrodynamic instability appears to be unrelated to the scales ℓ_{vel} and ℓ_{n} measured in the fully developed turbulent state [15].

Thus a simple scaling picture agrees well with the numerical results for extensile nematics. However, it may be noted that walls in these systems may also act as sources of vorticity. The role of walls becomes increasingly dominant as we reduce the activity in the system, for example the low activity regime (as in figure 1*c*,*d*) where the defects are not formed and the walls alone determine the entire dynamics. Scaling laws also break down at these low activities as can be observed in figure 6*a*. Moreover, the scaling picture presented above does not hold for contractile systems where, for example, the velocity–velocity correlation function does not show any simple scaling with *K*. This may be because flow in contractile systems is generated by both topological defects which comprise predominantly bend distortions and by the splay deformations which arise from hydrodynamic instabilities. It may also be due to the shear thickening known to occur in contractile systems which will introduce a dependence of *μ* on activity.

Our numerical results suggest many avenues for further research. For example, we lack a theory of why defects are created at a rate *ζ*^{2}/*K* or for the cross section for defect annihilation. Moreover, the mechanisms whereby defects and distortions create vorticity in active systems remain to be better understood. It will also be interesting to obtain similar numerical results for other systems which show active turbulence, for example those with polar symmetry or damped hydrodynamics, to identify which properties of the turbulent state are universal. The extent to which equations (2.1)–(2.7) are appropriate to describe experiments on microtubules [1], which are long elastic filaments, remains an open question. Also, the behaviour of these active systems in three dimensions, where the structure of walls and defects may be very different, remains to be explored.

In conclusion, we have performed numerical simulations of active turbulence in extensile and contractile nematics. We find that the same length scale is associated with correlations of the order parameter field and of the vorticity field. This suggests that defects and distortions are acting as sources of vorticity, in agreement with visual inspection of the simulation results, and provides evidence that vorticity is an important quantity in describing active turbulence [7].

## Funding statement

This work was supported by the ERC Advanced grant MiCE.

## Acknowledgements

The authors thank Z. Dogic, D. Pushkin and D. Chen for helpful discussions. The authors would also like to acknowledge the use of the Advanced Research Computing (ARC) in carrying out this work.

## Footnotes

One contribution of 10 to a Theme Issue ‘New trends in active liquid crystals: mechanics, dynamics and applications’.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.