## Abstract

Topological defects are distinctive signatures of liquid crystals. They profoundly affect the viscoelastic behaviour of the fluid by constraining the orientational structure in a way that inevitably requires global changes not achievable with any set of local deformations. In active nematic liquid crystals, topological defects not only dictate the global structure of the director, but also act as local sources of motion, behaving as self-propelled particles. In this article, we present a detailed analytical and numerical study of the mechanics of topological defects in active nematic liquid crystals.

## 1. Introduction

Active liquid crystals are non-equilibrium fluids composed of internally driven elongated units. Examples of active systems that can exhibit liquid crystalline order include mixtures of cytoskeletal filaments and associated motor proteins, bacterial suspensions, the cell cytoskeleton and even non-living analogues, such as monolayers of vibrated granular rods [1]. The key feature that distinguishes active liquid crystals from their well-studied passive counterparts is that they are maintained out of equilibrium not by an external force applied at the system's boundary, such as an imposed shear, but by an energy input on each individual unit. The energy fed into the system at the microscopic scale is then transformed into organized motion at the large scale. This type of ‘reverse energy cascade’ is the hallmark of active systems. In active liquid crystals, the large-scale self-organized flows resulting from activity further couple to orientational order, yielding a very rich behaviour. Novel effects that have been predicted theoretically or observed in simulations and experiments include spontaneous laminar flow [2–4], large density fluctuations [5–7], unusual rheological properties [8–10], excitability [11,12] and low Reynolds number ‘turbulence’ [11–17].

Ordered liquid crystalline phases of active matter can be classified according to their symmetry and to the nature of the forces that the active units exert on the environment. Active particles are often elongated objects with a head and a tail, hence intrinsically polar, such as bacteria or birds. Such systems can order in states with ferromagnetic (polar) order, where all units are on average aligned in a fixed direction. In this case, the ordered state is also a macroscopically moving state. Polar active particles can also order in states with nematic or apolar order, where the particles are aligned along the same axis, but with random head/tail orientation. In this case, the ordered state has apolar or nematic symmetry: if the direction of mean order is denoted by a unit vector ** n**, then the ordered state is invariant under and has zero mean velocity. Some active units are intrinsically apolar, such as vibrated rods [7], melanocytes [18] and some fibroblasts [19], and order in states with nematic symmetry.

In addition to the distinction based on symmetry of the ordered state, active systems can be further classified by the type of stresses or flows they impose on the surrounding medium. These can be contractile, as in actomyosin networks or in migrating cell layers, or extensile, as in suspensions of microtubule bundles or in most bacteria [1]. In the context of swimming microorganisms or artificial swimmers, units that exert extensile forces on the surrounding fluid are known as pushers, whereas those that exert contractile forces are known as pullers. Most bacteria are pushers, whereas the alga *Chlamydomonas* is an example of a puller. The extensile or contractile nature of active stresses affects the stability of ordered states [20].

In this paper, we focus on the rich dynamics of active liquid crystals with nematic symmetry. We consider both extensile [21,22] and contractile [19] systems. Previous work has highlighted the rich dynamics that arises in active nematics from the interplay of activity, orientational order and flow [2–4,11,12,23]. More recently, it was suggested that topological defects play an important role in mediating and driving turbulent-like active flows [14–17,21,24].

Topological defects are inhomogeneous configurations of the order field that provide a distinctive signature of liquid crystalline order and have been extensively studied in passive nematics. For passive nematics, defects may be generated through boundary conditions, externally applied fields, or via sufficiently rapid quenches from the disordered to the ordered state [25–27]. When the constraints are removed, or the system is given time to equilibrate, the defects ultimately annihilate, and the system reaches a homogeneous ordered state that minimizes the free energy. The structure of the topological defects is intimately related to the broken symmetry of the ordered state and effectively provides a ‘fingerprint’ of such a symmetry. When the ordered state has ferromagnetic (polar) symmetry, the lowest energy defect configurations are the charge +1 vortices and asters (monopoles), whereas in states with nematic symmetry charge defects, known as disclinations, are possible [28] and have the lowest energy. The structure of the topological defects therefore provides us an important tool for classifying the broken symmetry of liquid crystalline states.

In active liquid crystals, in contrast to passive ones, defect configurations can occur spontaneously in the bulk and be continuously regenerated by the local energy input, as demonstrated in experiments [18,19,21,22]. While the aster and vortex defects that occur in polar active systems [29,30] have been studied for some time [23,31–33], the properties of defects in active nematics have only recently become the focus of experimental and theoretical attention. Disclinations have been identified in monolayers of vibrated granular rods [7], in active nematic gels assembled *in vitro* from microtubules and kinesins [21], in dense cell monolayers [18,19], and in living liquid crystals obtained by injecting bacteria in chromonic liquid crystals [22]. In bulk suspensions of microtubule bundles, the defects drive spontaneous flows [21]. When confined at an oil/water interface, furthermore, the same suspensions form a two-dimensional active nematic film, with self-sustained flows resembling cytoplasmic streaming and the continuous creation and annihilation of defect pairs [21].

Recent work by us [24] and others [14–17,34] has begun to systematically examine the effect of activity on the dynamics of disclinations in a nematic liquid crystalline film. We demonstrated that disclinations in active liquid crystals behave like self-propelled particles with an active speed proportional to activity. The direction of active motion is controlled by the extensile or contractile nature of the active stresses. For certain relative orientations of pairs of defects of opposite strength this self-propulsion can overcome the equilibrium repulsive interaction among pairs of opposite-sign defects, allowing for dynamical states with an average sustained concentration of defect–antidefect pairs. In related work, Thampi *et al*. [14] suggested that the mean distance between defects in these turbulent states may be strongly correlated with the correlation length of fluctuations in the flow velocity, and only weakly dependent on activity. In the rest of this paper, we first review the hydrodynamic description of active nematics and the instabilities of the homogeneous ordered state. In §5, we present the results of a systematic numerical study of the various flow regimes induced by activity. Each regime is characterized in terms of flow patterns and defect proliferation. The chaotic regime exhibits a steady number of defects that persist in time. This is made possible by active flows that drive directed motion of the comet-like defects, generating, for certain relative orientations of two opposite sign defects, an effective repulsive interaction between the pair. The mechanisms for this active defect dynamics are analysed in §§3 and 4, where we discuss individual defect dynamics and pair annihilation in active nematics, respectively. We conclude with a brief discussion highlighting open questions.

## 2. Active nematodynamics

### (a) Governing equations

We consider a uniaxial active nematic liquid crystal in two dimensions. The two-dimensional limit is appropriate to describe the experiments by Sanchez *et al*. [21], where the microtubule bundles confined to a water–oil interface form an effectively two-dimensional dense nematic suspension, but also of considerable interest in its own right. The hydrodynamic equations of active nematic liquid crystals have been derived by coarse-graining a semi-microscopic model of cytoskeletal filaments cross-linked by clusters of motor proteins [35]. They can also simply be obtained from the hydrodynamic equations of passive systems by the addition of non-equilibrium stresses and currents owing to activity [1,3,11,12,32,36,37]. We consider here an incompressible suspension where the total density *ρ* of active bundles and solvent is constant. The equations are formulated in terms of the concentration *c* of active units, the flow velocity ** v** of the suspension and the nematic tensor order parameter

*Q*

_{ij}=

*S*(

*n*

_{i}

*n*

_{j}−

*δ*

_{ij}/2), with

**the director field. The alignment tensor**

*n**Q*

_{ij}is traceless and symmetric, and hence has only two independent components in two dimensions. The constraint of constant density

*ρ*requires ∇⋅

**=0. The hydrodynamic equations are given by [11] 2.1a 2.1b 2.1c where**

*v**D*/

*Dt*=∂

_{t}+

**⋅∇ indicates the material derivative,**

*v**D*

_{ij}=

*D*

_{0}

*δ*

_{ij}+

*D*

_{1}

*Q*

_{ij}is the anisotropic diffusion tensor,

*η*the viscosity,

*p*the pressure, λ the nematic alignment parameter and

*γ*the rotational viscosity. Here,

*u*

_{ij}=(∂

_{i}

*v*

_{j}+∂

_{j}

*v*

_{i})/2 and

*ω*

_{ij}=(∂

_{i}

*v*

_{j}−∂

_{j}

*v*

_{i})/2 are the symmetric and antisymmetric parts of the velocity gradient, corresponding to the strain rate and vorticity tensors. The molecular field

*H*

_{ij}=−

*δF*

_{LdG}/

*δQ*

_{ij}embodies the relaxational dynamics of the nematic obtained from the variation of the two-dimensional Landau–de Gennes free energy [28], with 2.2 where

*K*is an elastic constant with dimensions of energy. For the sake of simplicity, we restrict ourselves here to the one-elastic constant approximation to the Frank free energy, i.e. equal bend and splay moduli. The coefficients

*A*and

*C*determine the location of the continuous transition from a homogeneous isotropic state with

*S*=0 to a homogeneous nematic state with a finite value of

*S*given by (where we have used

*tr*

*Q*^{2}=

*S*

^{2}/2). We are interested here in a system where the transition is driven by the concentration of nematogens, as is the case for a fluid of hard rods of length ℓ which exhibits an isotropic-nematic (IN) transition at a concentration

*c*

^{⋆}=3

*π*/2ℓ

^{2}in two dimensions. Noting that

*A*and

*C*have dimensions of energy density (in two dimensions), we choose

*A*=

*K*(

*c*

^{⋆}−

*c*)/2 and

*C*=

*Kc*[12]. This gives in a homogeneous state of density

*c*

_{0}, so that

*S*

_{0}=0 for

*c*

_{0}<

*c*

^{⋆}and

*S*

_{0}≈1 for

*c*

_{0}≫

*c*

^{⋆}. The ratio defines a length scale that corresponds to the equilibrium correlation length of order parameter fluctuations and diverges at the continuous IN transition [26]. Here, we restrict ourselves to mean values of concentration well above

*c*

^{⋆}, where this equilibrium correlation length is microscopic and is of the order of the size ℓ of the nematogens, which will be used as our unit of length in the numerical simulation. Finally, the stress tensor is the sum of the elastic stress due to nematic elasticity, 2.3 where, for the sake of simplicity, we have neglected the Eriksen stress, and an active contribution, given by [1] 2.4 which describes stresses exerted by the active particles. The sign of

*α*

_{2}depends on whether the active particles generate contractile or extensile stresses, with

*α*

_{2}>0 for the contractile case and

*α*

_{2}<0 for extensile systems. Activity yields also a curvature-induced current given by the last term on the right-hand side of equation (2.1a),

*j*^{a}=−

*α*

_{1}

*c*

^{2}∇⋅

**, that drives active units from regions populated by fast moving particles to regions of slow moving particles. The**

*Q**c*

^{2}dependence of the active stress and current is appropriate for systems where activity arises from pair interactions among the filaments via cross-linking motor proteins. The active parameters

*α*

_{1}and

*α*

_{2}will be treated here as purely phenomenological. In microscopic models, they are found to depend on the concentration of active cross-linkers and on the consumption rate of adenosine triphosphate [38,39].

### (b) Linear stability

The hydrodynamic equations of an active nematic have two homogeneous stationary solutions, with *c*=*c*_{0} and ** v**=

**0**. For

*c*

_{0}<

*c*

^{⋆}, the homogeneous state is disordered with

*S*

_{0}=0. For

*c*

_{0}>

*c*

^{⋆}, the homogeneous solution is an ordered nematic state with . The linear stability of the ordered state has been studied in detail for the case of a

*L*×

*L*periodic domain [11,12]. It is found that above a critical activity the homogeneous state is unstable to a laminar flowing state. This instability corresponds to the spontaneous flow instability well-studied in a channel geometry [2,4,40]. The critical activity value associated with the instability of the homogeneous state is given by [12] 2.5 For the system is subjected to spontaneous splay deformations characterized by the instability of the first transverse mode. For , on the other hand, the instability is determined by the first longitudinal mode corresponding to bending deformations. Note that the sign of depends both on the sign of the fraction and the sign of the term 1∓λ in the denominator. Thus flow-aligning nematics (|λ|>1) are unstable to splay under the effect of an extensile active stress () and to bending under the effect of a contractile active stress (). Vice versa, flow-tumbling nematics (−1<λ<1) are unstable to bending under the effect of an extensile active stress () and to splay under the effect of a contractile active stress () (figure 1). In this paper, we focus exclusively on flow-tumbling systems and discuss both the cases of extensile and contractile stresses.

### (c) Dimensionless units and numerical methods

To render equations (2.1) dimensionless, we scale distances by the length of the active nematogens ℓ, set by the critical concentration *c*^{⋆}, stresses by the elastic stress of the nematic phase *σ*=*K*/ℓ^{2} and time by *τ*_{p}=*γ*ℓ^{2}/*K* representing the ratio between viscous and elastic stress. In these dimensionless units, we take *α*_{1}=|*α*_{2}|/2 and we introduce
2.6
as the fundamental measure of active stress. This serves as the control parameter throughout this work.

In the numerical calculations presented in §§4 and 5 the incompressibility condition ∇⋅** v**=0 is implemented by expressing equation (2.1)b in terms of vorticity and stream function, by writing
2.7
so that the incompressibility condition is automatically satisfied, whereas the vorticity field is given by
2.8
Equation (2.1)b is then cast in terms of vorticity in the form
2.9
where ∇×

**=∂**

*f*_{x}

*f*

_{y}−∂

_{y}

*f*

_{x}is the two-dimensional curl and we have neglected the inertial convective term to reproduce the typical low Reynolds number situation of cytoskeletal fluids (assuming the typical length of a filament ℓ∼1 μm, its typical velocity

*v*∼1 μm s and the same density and viscosity of water at room temperature,

*ρ*=10

^{3}kg m

^{−3}and

*η*=10

^{−3}Pa s, one finds

*Re*=10

^{−6}).

Finally, the numerical integration is performed via finite differences on a 256×256 square grid. The time integration was performed via a fourth order Runge–Kutta method with time step Δ*t*=10^{−3}. Except where mentioned otherwise, the numerical calculations described in this section use the parameter values *D*_{0}=*D*_{1}=1, λ=0.1, *c*_{0}=2*c*^{⋆} (corresponding to *S*_{0}=0.707) and *L*=20.

## 3. Dynamics of an isolated disclination

Topological defects are spatially inhomogeneous configurations of the director field that cannot be transformed continuously into a uniform state. In equilibrium, defects can occur upon quench from the disordered into the ordered phase or upon application of external electric or magnetic fields. Geometry or suitable boundary conditions can also be used to generate and maintain defects in the system. For instance, in a spherical nematic droplet, with boundary condition such that the director is normal to the droplet surface, the equilibrium configuration consists of a radial director field with a point defect at the centre of the droplet. In the absence of such constraints or external fields and given enough time to equilibrate defects of opposite sign always annihilate and the system settles into a uniform equilibrium state [41,42].

In two dimensions, defects are point-like. The strength of a disclination depends on how much the director field rotates around the defect core in one loop. In two dimensions, this can be expressed in terms of a single scalar field *θ* representing the angle formed by the director with the horizontal axis of a Cartesian frame. This gives
3.1
where the integral is calculated along an arbitrary contour enclosing the defect. The integer *k* is called strength of the defect and is analogous to the winding number of vortex defects in polar systems. In two-dimensional uniaxial nematics, the lowest energy defect configurations consists of half-strength disclinations with . In the presence of an isolated defect located at the origin, a solution *θ*=*θ*_{d} that minimizes the energy (2.2) and satisfies the constraint (3.1) is given by
3.2
with *ϕ* the usual polar angle. The corresponding energy is given by
3.3
where *R* is the size of the system and *a* is the core radius, defined as the radius of the region in the immediate proximity of the defect where the order parameter drops from its equilibrium value to zero. The quantity *ϵ*_{c} is the contribution to the total energy owing to the isotropic defect core.

Point defects in nematic liquid crystals have several particle-like features. Like charged particles, opposite-sign defects attract and same-sign defects repel ([28] and §4). It is well established that the dynamics of an isolated disclination that evolves according to equations (2.1) can be cast in the form of an overdamped equation of motion, given by
3.4
where ** r** is the defect position,

**is the net force acting on the defect, owing to interaction with other defects or externally imposed perturbations and**

*F***the local flow velocity at the position of the defect, which includes both external and self-generated flows. Finally,**

*v**ζ*is an effective drag coefficient proportional to the rotational viscosity

*γ*in equation (2.1)c and possibly space-dependent.

In the absence of fluid flow, an isolated disclination moves only in response to an externally imposed distortion and relaxes to the minimal energy texture *θ*_{d} given in equation (3.2). Following Denniston [43], one can then set *θ*=*θ*_{d}+*θ*_{ext}, where *θ*_{ext} expresses the departure from the optimal defective configuration *θ*_{d}, and calculate the energy variation with respect to a small virtual displacement of the defect core. For small deformations, this gives ** F**=−2

*πkK*∇

_{⊥}

*θ*

_{ext}, with ∇

_{⊥}=(−∂

_{y},∂

_{x}), and 3.5 where

*K*

_{1}and

*I*

_{1}are Bessel functions, and is the Ericksen number at the length scale of the defect core. The second equality in equation (3.5) implies

*a*≈0; for finite core radius, it introduces a dependence of the effective friction

*ζ*on the defect velocity. We refer the reader to references [43–45] for details.

The hydrodynamic coupling between the local orientation of the director and the flow gives rise to a self-generated flow, known as *backflow*, that in turn advects the defect core. When the dynamics of the flow is much faster than the orientational dynamics of the director, and in the absence of external forces, one can neglect ** F** in equation (3.4) and calculate the flow velocity

**from the Stokes equation 3.6 where**

*v***=**

*f*

*f*^{a}+

*f*^{e}=∇⋅(

*σ*^{a}+

*σ*^{r}) is the force arising from the active and elastic stresses acting in the system. In the case of isolated defects, this scenario is generally realistic for

*η*/

*γ*≪1 and even in a system containing multiple defects this purely advective dynamics continues to hold as long as the defects are sufficiently far apart (see §4). The general case in which both

**and**

*v***are non-zero was discussed by Kats**

*F**et al*. [46].

In the remainder of this section, we consider the regime in which *η*/*γ*≪1 and calculate the backflow owing to the stresses arising in the presence of isolated disclinations. Let us then consider a disclination located at the origin of a circular domain of size *R*. The domain might represent either the entire system or, more realistically, the defect-free portion of the system surrounding a given central defect. We refer to this as the *range* of a defect. Because of the linearity of the Stokes equation, the solution of equation (3.6) can be written as ** v**=

*v*^{0}+

*v*^{a}+

*v*^{e}, where

*v*^{0}is the solution of the homogeneous Stokes equation, whereas

*v*^{a}and

*v*^{e}are the flows produced by the active and elastic force, respectively. The solution can be expressed as the convolution of the two-dimensional Oseen tensor with the force per unit area 3.7 where

*G*

_{ij}is the two-dimensional Oseen tensor [47], given by 3.8 with a length scale adjusted to obtain the desired behaviour at the boundary. Taking , with , and assuming uniform concentration and nematic-order parameter outside the defect core, the body force owing to activity can be calculated straightforwardly as where, for the sake of simplicity, we have assumed

*S*

_{0}=1 outside the core. Using this in equation (3.7), and using (3.8), yields, after some algebraic manipulation, 3.9a and 3.9b A plot of these flow fields is shown in figure 2. Setting in equation (3.4), we find that active disclinations self-propel at constant speed along their symmetry axis ( in this setting), and their equation of motion can be written as 3.10 where

*v*

_{0}=

*αR*/(4

*η*). On the other hand, , and so disclinations are not propelled by the active backflow, but rather move solely under the effect of the elastic force produced by other defects. It is crucial to note that the self-propulsion speed

*v*

_{0}scales linearly with the active stress

*α*and so disclinations in contractile (

*α*>0) and extensile (

*α*<0) active nematic suspensions self-propel in opposite directions.

Equations (3.9) can be complemented with various kinds of boundary conditions, which, however, have no effect on the qualitative features of the solution and, more importantly, on the dynamics of the defects. To illustrate this point we consider a slippery interface, such that *v*_{r}(*R*,*ϕ*)=0 and , where *ξ* is the coefficient associated with the frictional force exerted by the interface on the fluid. If the domain of equation (3.6) is interpreted as a container, then *ξ* is the actual frictional coefficient of the container wall. On the other hand, if the domain is interpreted as the *range* of a defect, then *ξ*∼*ηh*, where *h* is the thickness of the boundary layer between the ranges of neighbouring defects. No-slip boundary conditions can be recovered in the limit .

A solution *v*^{0} of the homogeneous Stokes equation enforcing the boundary conditions for the active backflow induced by the disclination can be found from the following biharmonic stream function
3.11
with *a*_{1} and *b*_{1} constants. The corresponding velocity field is given by
3.12
Then, setting and , with , and solving for *a*_{1} and *b*_{1} yields
The associated self-propulsion speed *v*_{0} in equation (3.10) then becomes
3.13
or *v*_{0}=*αR*/(12*η*) in the no-slip limit. Thus, as anticipated, incorporating the effect of the boundary only changes the speed of the defects without altering their dynamics. Similarly, in the case of a disclination, we can consider the biharmonic stream function
3.14
whose associated velocity field is given by
3.15
Setting and and solving for *a*_{3} and *b*_{3} yields
Clearly, this does not change the symmetry of the active backflow, and thus negative-charge disclinations are stationary.

In summary, half-strength disclinations in active nematic liquid crystals can be described as self-propelled particles with overdamped dynamics governed by an equation of the form (3.4). In the absence of external forces, or forces owing to interactions with other defects, the disclination core is advected at constant speed by a self-generated backflow, provided the dynamics of the flow is faster than the relaxational dynamics of the director (i.e. *η*/*γ*≪1). As a result, positive disclinations travel along their symmetry axis at speed *v*_{0}∼*αR*/*η*, whereas negative disclinations are stationary and move only as a consequence of their interaction with other disclinations. The direction of motion is controlled by the sign of the activity *α*, which in turns depends on whether the system is contractile (*α*>0) or extensile (*α*<0). Thus, the comet-like disclination travels in the direction of its ‘tail’ in contractile systems and in the direction of its ‘head’ in extensile systems. We note that active curvature currents in the concentration equation (2.1a) controlled by the parameter *α*_{1} have a similar effect, as noted by Narayan *et al*. in a system of vibrated rods [7] and recently investigated by Shi & Ma [48] through extensive numerical simulations. Such curvature-driven currents control the dynamics in systems with no momentum conservation, but are very small in the regime discussed here, where active effects are driven almost entirely by the active stress given in equation (2.4). This is because the curvature-driven currents controlled by the active diffusivity *D*_{α}=*α*_{1}*c* in equation (2.1a) are negligible compared with diffusive currents controlled by *D*_{0}. To see this, we take the active stress *α* to be of order of the force per unit area exerted by motor proteins *α*∼*pN* μm^{−2}∼1 Pa (note that we use three-dimensional quantities for our estimates here). The active diffusivity *α*_{1} also arises from active stresses (see reference [49]) and can be estimated as *α*_{1}*c*∼*α*/*nζ*, where *n* is the number density of the solvent and *ζ*∼*η*ℓ is the friction between the microtubule bundles and the solvent. Using the Stokes–Einstein relation to estimate *D*_{0}=*k*_{B}*T*/6*πη*ℓ, we find *D*_{α}/*D*_{0}=*α*/*nk*_{B}*T*∼10^{−5}, where we have taken *nk*_{B}*T*=1 atm=10^{5} *Pa* as the hydrostatic pressure of the fluid. Although the effective *D*_{0} is likely to be increased by both high density and activity, we expect active diffusion currents to be negligible compared with passive ones. This is supported by the experiments that are in a regime where the density is very high and shows little spatial variations. Our numerics are indeed performed in this regime and show very little variations of the density even at the largest activity.

Special attention should be devoted to the fact that the self-propulsion speed *v*_{0} depends linearly on the range *R* of a defect, this being defined as the defect-free portion of the system surrounding a defect (and possibly coinciding with the entire system). Although two-dimensional hydrodynamics is known to be plagued with anomalies, such as the Stokes paradox [50], this behaviour does not result solely from the two-dimensionality of the problem. Point defects in three dimensions would also yield a similar behaviour. To see this, we note that the deviatoric part of the Oseen tensor scales like *r*^{2−D}, with *D* the space dimension. On the other hand, the active force *f*^{a}=*α*∇⋅** Q** always scales such as

*r*

^{−1}. The backflow velocity in

*D*dimensions thus scales like , regardless of the space dimension. In three-dimensions disclinations are, however, line defects. In this case, denoting by

*ξ*

_{d}the persistence length of the disclinations, i.e. the length scale over which these line defects can be treated as straight lines, the length

*R*controlling the flow velocity induced by a defect would scale as , with

*L*the system size and

*a*the core radius. These results could also be obtained on the basis of dimensional analysis by noting that

*v*^{a}is always proportional to

*α*/

*η*. Because

*α*has dimensions of stress and

*η*has dimensions of stress over time, the resulting velocity must also be proportional to a length scale. This length scale was first noted experimentally by Sanchez

*et al.*[21] and investigated numerically by Thampi

*et al*. [16,14]. Its nature, however, remains elusive (see §5 for further discussion on this matter).

The growth of the flow field generated by a defect at large distances can also be cut off by a frictional force *f*_{s}=−*ζ*_{s}** v** as may arise from the fact that the nematic film is confined at an oil/water interface [51]. Such a frictional interaction with the subphase removes energy from the flow at the length scale , thus controlling the decay of the velocity field. Finally, the limit where the friction dominates viscous forces corresponding to a no-slip Hele–Shaw geometry has been discussed in detail by Pismen [34]. In this case, the flow generated by a single disclination is found to decay as ∼

*r*

^{−3}at large distances from the defect.

## 4. Annihilation dynamics of defect pairs

Here, we discuss the annihilation of a pair of oppositely charged disclinations. The study of the annihilation dynamics of defect–antidefect pairs is a mature topic in the liquid crystals field and has been subjected to numerous investigations [43,44,52–56]. In the simplest setting [28], one considers a pair of disclinations located at *r*_{±}=(*x*_{±},0) and separated by a distance Δ=*x*_{+}−*x*_{−} (figure 3). The energy of the pair is given by
4.1
Each defect experiences an elastic force of the form and thus, in the absence of backflow, equation (3.4) can be cast in the form
4.2
where *κ*=2*πk*^{2}*K*/*ζ*. This yields
4.3
so that the distance between annihilating defects decreases as a square-root, , with *t*_{a} the annihilation time. More precise calculations have shown that the effective friction is itself a function of the defect separation [44,45], , although this does not imply substantial changes in the overall picture. This simple model predicts that the defect and antidefect approach each other along symmetric trajectories and annihilate at Δ(0)/2 in a time *t*_{a}=Δ^{2}(0)/4*κ*. The backflow produced by the balance of elastic and viscous stresses [52,56], as well as the anchoring conditions at the boundary [53], can produce an asymmetry in the trajectories of the annihilating defects or even suppress annihilation when the defects are initially far from each other or the anchoring is sufficiently strong.

To understand how activity changes the simple annihilation dynamics described so far, we have integrated numerically equations (2.1) for an initial configuration of uniform concentration and zero flow velocity, with two disclinations of charge symmetrically located with respect to the centre of the box along the *x*-axis. Figure 3 shows a snapshot of the order parameter and flow field shortly after the beginning of the relaxation for both a contractile and extensile system, with *α*=±0.8 in the units defined in §2c.

In contractile systems active backflow yields a net speed-up of the defect towards its antidefect for the annihilation geometry shown in figure 3*a*. In extensile systems, with *α*<0, backflow drives the defect to move towards its head, away from its partner in the configuration of figure 3*b*, acting like an effectively *repulsive* interaction (figure 3*b*). If the initial positions of the defects are exchanged, the behaviour is reversed. The effective attraction or repulsion between oppositively charged active defects is thus dictated by both the contractile or extensile nature of the active stresses, which determines the direction of the backflow, and the relative orientation of the defects, as summarized in figure 4. This effect has been observed in experiments with extensile microtubules and kinesin assemblies [21] and can be understood on the basis of the hydrodynamic approach embodied in equations (2.1). In figure 5, we have reproduced from reference [21] a sequence of snapshots showing a pair of disclinations moving apart from each other together with the same behaviour observed in our simulations.

Figure 6*a* and *b* shows the trajectories of the active defects, with the red and blue line representing the and disclination, respectively. The tracks end when the cores of the two defects merge. For small activity and small values of the rotational viscosity *γ*, the trajectories resemble those obtained for passive systems [52,56]. At large values of activity, however, the asymmetry in defect dynamics becomes more pronounced, and when the activity dominates over orientational relaxation, the disclination moves independently along its symmetry axis with a speed *v*_{0}∼*αR*/*η* (see §3) whose direction is dictated by the sign of *α*. This behaviour is clearly visible in figure 6*c*, showing the defect separation Δ(*t*) as a function of time. For *γ* sufficiently large, the trajectories are characterized by two regimes. For large separation, the dynamics is dominated by the active backflow, and thus and Δ(*t*)∝−*αt*. Once the defects are about to annihilate, the attractive force takes over, and the defects behave as in the passive case with .

This behaviour can be understood straightforwardly from the basic concepts of active defect dynamics discussed in §3. Each defect in the pair travels in space according to equation (3.4), with ** v** given by and

*v*_{±}given in equations (3.9) plus a suitable homogeneous solution of the Stokes equation that enforces the periodic boundary conditions. Next, we retain only the active contribution to the backflow and replace the flow profiles by their constant values at the core of the defect, with and

*v*_{−}(

*x*

_{−})=

**0**. This yields the following simple equation for the pair separation 4.4 Equation (4.4) explicitly captures the two regimes shown in figure 6

*c*and described earlier. The solution takes the form 4.5 The pair annihilation time

*t*

_{a}is determined by Δ(

*t*

_{a})=0 and is given by 4.6 For passive systems (

*α*=0), this reduces to . This predicts that the annihilation time, normalized to its value in passive systems, , depends only on Δ(0)

*v*

_{0}/2

*κ*∼

*αγ*. Figure 6

*d*shows a fit of the annihilation times extracted from the numerics to this simple formula. The model qualitatively captures the numerical behaviour. Note also that defects created a finite distance Δ(0) apart will always separate provided the activity and hence the magnitude of the self-propulsion velocity exceeds a critical value .

## 5. Defect proliferation

Pair annihilation is the fundamental mechanism behind defect coarsening in nematic liquid crystals [41]. Once this mechanism is suppressed by activity, as described in §4, the coarsening dynamics is replaced by a new steady state in which pairs of disclinations are continuously produced and annihilated at constant rate. The chaotic dynamics following from continuous defect proliferation and annihilation results in turbulent flow [14–17,24]. In this section, we describe the onset of chaos and the proliferation of defects that are observed from numerical solutions of equations (2.1) upon varying the activity parameter *α* and the rotational viscosity *γ*. As initial configurations, we take a homogeneous state with the director field aligned along the *x*-axis and subjected to a small random perturbation in density and orientation. The equations were then integrated from *t*=0 to *t*=2×10^{3}*τ* (see §c for a description of units).

Figure 7 summarizes the various regimes obtained by exploring the (*α*,*γ*)−plane (i) a homogeneous, quiescent ordered state (H); (ii) a periodic flow marked by the emergence of relaxation oscillations (O); (iii) a non-periodic oscillatory flow characterized by the formation of ‘walls’ in the nematic phase and the unzipping of these walls through the unbinding of defect pairs (W); and (iv) a chaotic or ‘turbulent’ state associated with a constant defect density (T). The latter three regimes are described in more detail in the following.

### (a) Relaxation oscillations

As described in references [11,12], relaxation oscillations occur in active nematics when |*α*| exceeds as the result of the competition of two time scales: the relaxation *τ*_{p}=*γ*ℓ^{2}/*K* of the nematic structure and the time scale *τ*_{a}=*η*/|*α*| that controls the rate at which active stresses are injected in the system. When *τ*_{p}<*τ*_{a}, the microstructure can relax to accommodate the active forcing and the ordered state is stable. This is a quiescent state, with uniform order parameter. Conversely, when *τ*_{p}>*τ*_{a}, the relaxation of the nematic structure lags behind the injection of active stresses, yielding various dynamical states with spatially and temporally inhomogeneous order parameter.

In this regime, the dynamics consists of a sequence of almost stationary passive periods separated by active ‘bursts’ in which the director switches abruptly between two orthogonal orientations (figure 8*b*–*d*). During passive periods, the particle concentration and the nematic-order parameter are nearly uniform across the system, there is no appreciable flow, and the director field is either parallel or perpendicular to the *x*-direction (owing to the initial conditions). Eventually, this configuration breaks down and the director field rotates by 90^{°}. The rotation of the director field is initially localized along narrow-extended regions, generating flowing bands similar to those obtained in active nematic films [2] (figure 8*c*). The temporary distortion of the director field as well as the formation of the bands is accompanied by the onset of flow along the longitudinal direction of the bands, with neighbouring bands flowing in opposite directions. The flow terminates after the director field rotates and a uniform orientation is restored. The process then repeats.

Depending on whether the active stress fuelling the oscillatory dynamics is contractile or extensile, the rotation of the director occurs through an intermediate splay or bending deformation. During bursts, the nematic-order parameter, otherwise equal to its equilibrium value, drops significantly (figure 8*a*, black line). Without this transient melting the distortions of the director field required for a burst are unfavourable for any level of activity.

The frequency of the oscillation is proportional to *k*^{2}*α*, where *k*=2*π*/*L* is the wave number of the longest-wavelength mode to go unstable [11,12]. In spite of the strong elastic deformation and the dramatic drop in the order parameter, this regime contains no unbound defects.

### (b) Wall formation and unzipping

For larger values of activity, the bend and splay deformations of the director at the band boundaries become large enough to drive creation of defect pairs, as shown in figure 8*e*–*g*. The alignment of the bands again oscillates between two orthogonal directions (*x*- and *y*-axis for the initial condition used here), but the switching takes place through an intermediate more complex configuration with lozenge-shaped ordered regions. Defect pairs then unbind and glide along the narrow regions separating two bands where elastic deformations and shear flows are largest. These regions of high distortion and shear rate are commonly referred to as ‘walls’ in the liquid crystals literature [57]. Movies displaying the dynamics of wall unzipping and defect creation for both extensile and contractile systems are included as the electronic supplementary material.

The formation of walls and the ‘unzipping’ of defects along the walls by the creation of pairs of disclinations has being discussed in detail by Thampi *et al*. [15] (see also reference [16]). Defect unbinding along the walls relaxes both the excess elastic energy and the high shear stresses present in these regions, where the nematic order parameter *S* is driven near zero. For this reason, wall formation and unzipping occurs often at the boundary between pairs of vortices of opposite circulation (figure 9). In passive nematic liquid crystals, the strong bending deformation that leads to the formation of walls, and precedes defect unbinding, requires an external action, such as an applied electric field or an externally imposed shear stress [58]. In active nematics, on the other hand, a similar complex spatio-temporal dynamics occurs spontaneously, driven by the *local* injection of active stresses, which are in turn balanced by spontaneous distortions and flows as described in §2*b*. Wall formation and unzipping also combines here with the oscillatory dynamics described in §5*a* to give rise to the periodic creation and annihilation of topological defects that marks the transition to the turbulent regime.

Oscillating band structures of the type observed here are found in passive nematic fluids under externally applied shear flows and are precursors to rheochaos [59,60]. They have been predicted theoretically [61] and observed experimentally [62] in suspensions of worm-like micelles, where the spatio-temporal dynamics is directly correlated to shear banding [63] and to stress/shear fluctuations at the shear band interface. In passive systems it has been argued that the route to rheochaos depends on whether the stress or the strain rate is controlled during the experiment. It would be interesting to analyse in more detail the route to chaos in this case.

### (c) Turbulence

At even higher activity, the bands/walls structure begins to bend and fold and the dynamics becomes chaotic, resembling that of a ‘turbulent’ fluid, as displayed in figure 8*h*–*j* and in the movies included as the electronic supplementary material. The system reaches a dynamical steady state where defect pairs are continuously created and annihilated, but their mean number remains on average constant in time. Owing to topological charge conservation, at any time the system contains an equal number of positive and negative defects. Unlike in equilibrium systems, however, in active nematics, it is possible for opposite-charge defects of certain orientations to repel instead of attracting (see §4). This allows the formation of a defective steady state, with self-sustained flows and a constant mean number of defects.

The defective dynamics observed in the system is similar to that obtained in a passive nematic subject to an externally imposed shear or to electrohydrodynamic instabilities [64]. In active nematics, however, defects are generated in the absence of any externally imposed forces or constraints, as a result of the spontaneous distortion induced by the active stresses. After formation, the defects are convected by the swirling flow and interact with one another through distortional elasticity as well as hydrodynamically through the modifications the defects themselves induce on the flow field.

The sequence of flow regimes observed here is reminiscent of the evolution of the dynamics of sheared tumbling nematic polymers with increasing shear rate, known as the Ericksen number cascade [65]. In a nematic film sheared at a rate between two plates separated by a distance *h*, the Ericksen number *Er* provides a dimensionless measure of the magnitude of viscous torques () relative to elastic torques arising from spatial gradients in the average molecular orientation (∼*K*/*γh*^{2}), with . This suggest the definition of an *active* Ericksen number *Er*_{α} where the active stress *α* takes the place of viscous stresses . The active Ericksen number is then defined as *Er*_{α}=*αγL*^{2}/(*ηK*). With this definition, *Er*_{α} can also be interpreted as the ratio of the time scale *τ*_{p}=*γL*^{2}/*K* for the relaxation of an orientational deformation of the nematic order on the scale of the entire system to the time *τ*_{a}=*η*/*α* controlling the injection of active stresses.

Figure 10 shows a plot of the *area fraction* occupied by defects as a function of the active Ericksen number *Er*_{α} for both extensile and contractile systems. The area fraction is defined as the relative area occupied by the core of the defects, or *Nπa*^{2}/*L*^{2}, with *N* the number of defects. The core radius *a* resulting from the hydrodynamic equations (2.1) is approximatively given by the size of the boundary layer between the position of a defect, where *S*=0, and surrounding space, where *S*=*S*_{0}. This is related to the coefficient *A* in the Landau–de Gennes free energy (2.2)
5.1
For larger activity, the reduction of the nematic order parameter owing to the unbound defects compensates the activity increase by effectively reducing the injected local stress *σ*^{a}≈*αS*.

For the same magnitude of activity, extensile systems contain a larger mean number of defects than contractile ones. In both cases, defect pairs first unbind within the walls, which are regions of large bend and splay deformations in extensile and contractile systems, respectively. As the system evolves towards the regime of chaotic dynamics, the walls begin to deform largely via bend deformations in both systems [15,16]. This asymmetry could be because the severe splay deformations localized at the walls lead to a more drastic reduction of the nematic order parameter. Thus contractile flow-tumbling nematics, whose spontaneous distortion involves mostly splay deformations, are effectively less active than flow-tumbling extensile systems.

## 6. Discussion and conclusion

We have presented a detailed analytical and numerical study of the mechanics of topological defects in active nematic liquid crystals. Topological defects are distinctive signatures of liquid crystals and profoundly affect their viscoelastic behaviour by constraining the orientational structure of the fluid in a way that inevitably requires system-wide (global) changes not achievable with any set of local deformations. In ordered states of both passive and active nematics, the topological defects are fingerprints of the broken symmetry in the ordered state. In particular, the presence of strength defects clearly reveals the nematic nature of the orientational order, in contrast to systems with polar (ferromagnetic) symmetry where the lowest energy defects allowed have strength ±1. Active liquid crystals have the additional feature that defects act as local sources of motion, behaving as self-propelled particle-like objects (see §3). The direction of motion of the strength defects provides, furthermore, a clear signature of the extensile or contractile nature of the active stresses, as the comet-like positive defects are advected towards their head in extensile systems and towards their tail in contractile ones. In passive liquid crystals, defect dynamics is always transient, as oppositely charged defects attract and eventually annihilate. In active nematics, on the other hand, the interplay between active and viscous stresses, modulated by the director geometry induced by the defects, enriches the spectrum of defect–defect interactions by allowing for an effective repulsion between defects and anti-defects (see §4). For highly active systems, this mechanism can arrest the process of coarsening, leading to a state where unbound pairs of defects are continuously created and annihilated, but with their mean density constant in time (see §5). The chaotic dynamics originating from the continuous defect proliferation and annihilation results in spontaneous low Reynolds number turbulence, akin to the so-called director turbulence seen in sheared polymer and micellar nematics [65].

Several open questions remain concerning the defect dynamics of active nematics. The defect area fraction shown in figure 10 exhibits a crossover from growth at low activity to saturation at high activity, but an understanding of the length scales that control this behaviour is still lacking. Both our work and work by Thampi *et al*. [16] suggest that the mean separation between defects in the turbulent regime coincides with the typical vortex size, but more work is needed to elucidate the behaviour of this length scale with activity over a wide range of parameters. While defect generation in sheared nematics has been explained in terms of a simple rate equation that balances creation and annihilation [65], a similar simple model for active defects is still missing. On the basis of numerics, Thampi *et al*. [16], have suggested that the defect creation rate should scale like the square of the activity but no simple argument is available to understand this counterintuitive result. Finally, an important open question is the different behaviour of extensile and contractile systems apparent from figure 10. The turbulent state of contractile systems is much less defective than that of extensile ones, suggesting that for equal magnitude of the activity *α*, contractile nematics are effectively ‘less active’ than extensile ones.

## Funding statement

L.G. was supported by SISSA Mathlab. The work at Syracuse was supported by the National Science Foundation through awards DMR-1305184 (MCM, P.M.) and DGE-1068780 (M.C.M.) and by funds from the Syracuse Soft Matter Program. M.C.M. also acknowledges support from the Simons Foundation. Finally, M.C.M. and M.J.B. thank the KITP at the University of California, Santa Barbara, for support during the preparation of the manuscript.

## Acknowledgements

We thank Zvonimir Dogić, Suzanne Fielding, Jean-Francois Joanny, Oleg Lavrentovich and Tim Sanchez for several illuminating discussions.

## Appendix A. Active backflow of disclinations

Finding the backflow produced by the positive disclination reduces to the calculation of the following integrals
A 1a
and
A 1b
To calculate the first integral, we can make use of the logarithmic expansion
A 2
where and . The integral over the angle can be immediately carried out using the orthogonality of trigonometric functions: . The remaining radial integral can be then straightforwardly calculated
A 3
Thus,
To calculate the integral *I*_{2}, we can make use of the fact that
A 4
Thus, one can write
A 5
The first integral is calculated by differentiating *I*_{1}
A 6
The second integral in equation (A 5) can be calculated again with the help of the logarithmic expansion (A 2)
A 7
where we have used again the orthogonality of trigonometric functions
A 8
Thus, taking the derivative and combining with equation (A 9) yields
A 9
Adding together *I*_{1} and *I*_{2} and switching to polar coordinates gives
A 10a
and
A 10b
Setting and expressing everything in polar coordinates one finally obtains equation (3.9a).

## Appendix B. Active backflow of disclinations

The calculation of the active backflow associated with a negative disclination reduces to the calculation of the following integrals
B 1
and
B 2
The first integral can be calculated with the help of the logarithmic expansion (A 2) as well as the orthogonality condition (A 8). This yields
B 3
To calculate *I*_{4}, we can use again equation (A 4). Thus
B 4
The first and third integral in equation (B 4) are respectively the opposite of the *y*- and *x*-component of *I*_{3}. The remaining two integrals, can be calculated using equation (A 2) and (A 8). This gives
B 5
where and . Combining *I*_{3} and *I*_{4} and switching to polar coordinates finally gives equation (3.9c).

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