Hydrodynamics of Active Defects: from order to chaos to defect ordering

TL;DR

This paper develops a hydrodynamic theory for 2D active nematics with many defects, revealing transitions from active turbulence to defect-ordered phases and showing how activity gradients can control defect dynamics.

Contribution

It introduces a comprehensive hydrodynamic model capturing defect-driven flows, chaos, and ordered phases in active nematics, and explores defect manipulation via activity gradients.

Findings

Characterization of defect chaos and turbulence in active nematics.

Discovery of a defect-ordered polar liquid phase at high activity.

Demonstration of activity gradients controlling defect charge sorting.

Abstract

Topological defects play a prominent role in the physics of two-dimensional materials. When driven out of equilibrium in active nematics, disclinations can acquire spontaneous self-propulsion and drive self-sustained flows upon proliferation. Here we construct a general hydrodynamic theory for a two-dimensional active nematic interrupted by a large number of such defects. Our equations describe the flows and spatio-temporal defect chaos characterizing active turbulence, even close to the defect unbinding transition. At high activity, nonequilibrium torques combined with many-body screening cause the active disclinations to spontaneously break rotational symmetry forming a collectively moving defect ordered polar liquid. By recognizing defects as the relevant quasiparticle excitations, we construct a comprehensive phase diagram for two-dimensional active nematics. Using our hydrodynamic approach, we additionally show that activity gradients can act like "electric fields", driving the sorting of topological charge. This demonstrates the versatility of our continuum model and its relevance for quantifying the use of spatially inhomogeneous activity for controlling active flows and for the fabrication of active devices with targeted transport capabilities.

Equations

$$\Psi (\mathbf{r})=S(\mathbf{r})e^{2i\theta (\mathbf{r})},$$

$$n_{\pm }(\mathbf{r},t)$$

$${\mathbf{j}}_{\pm }(\mathbf{r},t)$$

$${\mathbf{v}}_n=\frac{1}{2|\Psi {|}^2}\mathrm{Im}\left({\Psi }^{\ast }\nabla \Psi \right)=\nabla \theta .$$

$$\hat{\mathbf{z}}\cdot \left(\nabla \times {\mathbf{v}}_n\right)=2\pi \rho ,$$

$$\mathbf{p}(\mathbf{r},t)=⟨\sum _i{\mathbf{e}}_i(t)\delta (\mathbf{r}-{\mathbf{r}}_i^{+}(t))⟩.$$

$${\partial }_t{n}_{\pm }+\nabla \cdot {\mathbf{j}}_{\pm }=W_c-W_a,$$

$${\mathbf{j}}_{+}$$

$${\mathbf{j}}_{-}$$

$${\oint }_{\Gamma }\mathrm{d}\theta \equiv {\oint }_{\Gamma }\mathrm{d}\mathbf{s}\cdot {\mathbf{v}}_n=2\pi \sum _{i=1}^N{q}_i.$$

$${\partial }_t\rho +\nabla \cdot {\mathbf{j}}_{\rho }=0,$$

$${\partial }_{t}n+\nabla \cdot {\mathbf{j}}_n=W_c-W_a,$$

$${\mathbf{j}}_{\rho }=\frac{({\mathbf{j}}_{+}-{\mathbf{j}}_{-})}{2}=\frac{v}{2}\mathbf{p}+\nu \kappa n\mathbf{ϵ}\cdot {\mathbf{v}}_n-D_0\nabla \rho ,$$

$${\mathbf{j}}_n=\frac{({\mathbf{j}}_{+}+{\mathbf{j}}_{-})}{2}=\frac{v}{2}\mathbf{p}+\nu \kappa \rho \mathbf{ϵ}\cdot {\mathbf{v}}_n-D_0\nabla n.$$

$${\partial }_t{\mathbf{v}}_n=2\pi \mathbf{ϵ}\cdot {\mathbf{j}}_{\rho }+\kappa \nabla (\nabla \cdot {\mathbf{v}}_n),$$

$${\partial }_t\mathbf{p}-\nu \kappa \hat{\mathbf{z}}\cdot ({\mathbf{v}}_n\times \nabla )\mathbf{p}=-\left[D_R+2\pi \nu \kappa \rho +\beta |{\mathbf{v}}_n{|}^2\right]\mathbf{p}-v_2{n}_{+}\mathbf{ϵ}\cdot {\mathbf{v}}_n-\frac{v_1}{2}\nabla n_{+}-2\kappa (\nabla \cdot {\mathbf{v}}_n)\mathbf{ϵ}\cdot \mathbf{p}+D_0{\nabla }^2\mathbf{p}.$$

$$v_1=v\frac{2T}{K}f(\zeta ),v_2=v\frac{3\nu T}{4K}f(\zeta ),$$

$$\zeta =\frac{v^2\gamma T}{D_R{K}^2},$$

$${\partial }_{t}n$$

$${\partial }_t{\mathbf{v}}_n$$

$${\partial }_t\mathbf{p}$$

$$D_R\left[a_2-a_4|{\mathbf{p}}_0{|}^2\right]{\mathbf{p}}_0=0,$$

$$a_2=\frac{3}{8}\zeta f(\zeta )-1,a_4=\zeta \frac{5K}{32T{n}_0^2}.$$

$${\zeta }_{c}f({\zeta }_c)=\frac{8}{3}⟹{\zeta }_c\approx 0.701,$$

$$|v_{c_1}|=\sqrt{\frac{2\nu \kappa D_R\left(1-\tilde{T}\right)}{\pi \left[1+3\nu \tilde{T}/32\right]}},$$

$$|v_{c_2}|=\sqrt{{\zeta }_c\frac{D_R{K}^2}{\gamma T}}\simeq 0.84\sqrt{\frac{D_R{K}^2}{\gamma T}},$$

$$\left(\frac{T}{K}\right)\frac{\gamma {\ell }_p^2}{K}\gtrsim D_R^{-1}.$$

$$|{\mathbf{p}}_0|\simeq n_0\sqrt{\frac{T}{K}}{\left(\frac{\zeta -{\zeta }_c}{{\zeta }_c}\right)}^{1/2},$$

$${\mathbf{v}}_n^0=\frac{v}{2\nu \kappa n_0}\mathbf{ϵ}\cdot {\mathbf{p}}_0\propto {\left(\frac{\zeta -{\zeta }_c}{{\zeta }_c}\right)}^{1/2},$$

$$\xi =\frac{1}{\sqrt{2\pi \nu n_0}}.$$

$${\partial }_t\delta p_L$$

$${\partial }_t\delta {\mathrm{v}}_n^T$$

$${\partial }_t\delta p_T$$

$${\partial }_t\delta {\mathrm{v}}_n^L$$

$$i{\omega }_{\pm }=\frac{1}{2}\left[\frac{1}{\tau }+D_R\pm \sqrt{{\left(D_R-\frac{1}{\tau }\right)}^2+\frac{3D_R}{2\tau }\zeta f(\zeta )}\right],$$

$$⟨|\delta {\mathbf{p}}_{\mathbf{q}}(t){|}^2⟩\approx n_0\frac{2T}{K}f(\zeta )+\frac{n_0{D}_0(v_2/\nu \kappa {)}^2}{2D_R\left(1+q^2{\xi }^2\right)\left[1+\tau D_R+q^2{\xi }^2\right]}.$$

Full text

I Introduction

In recent years, the framework of active fluids has had significant success in describing the emergent collective motion of bacterial and cellular assemblies wensink2012meso ; saw2017topological ; kawaguchi2017topological ; blanch2018turbulent ; li2019data . Such spontaneous self-organization into complex patterns on large scales is a characteristic feature of active matter, whose constituent units are individually self-driven ramaswamy2010mechanics ; marchetti2013hydrodynamics . Much interest has recently focused on active nematics ramaswamy2003active ; doostmohammadi2018active - collections of head-tail symmetric elongated units that exert active forces on their surroundings and organize in states with apolar orientational order. A rapidly growing list of experimental realizations of active nematics ranges from suspensions of cytoskeletal filaments and associated motor proteins sanchez2012spontaneous ; keber2014topology ; ellis2018curvature ; kumar2018tunable to vibrated granular rods narayan2007long , and colonies of living cells Zhou2014 ; nishiguchi2017long ; duclos2017topological ; saw2017topological ; kawaguchi2017topological ; blanch2018turbulent ; li2019data . In all these systems the interplay of orientational order and self-sustained active flows yields a rich collection of dynamical states, including spontaneous laminar flows voituriez2005spontaneous ; duclos2018spontaneous and spatio-temporal chaos or “active turbulence” with the proliferation of topological defects thampi2013velocity ; giomi2013defect ; shankar2018defect .

While a great deal of understanding has been gained from extensive numerical work thampi2013velocity ; giomi2013defect ; shi2013topological ; thampi2014instabilities ; gao2015multiscale ; giomi2015geometry ; hemingway2016correlation ; alert2019universal and experiments guillamat2017taming ; lemma2018statistical ; martinez2019selection , the theoretical analysis of low Reynolds number turbulence in active nematics, beyond deterministic linear instabilities, remains an open problem. A notable exception is Ref. giomi2015geometry , where a phenomenological mean field theory was proposed to describe active nematic turbulence on short scales by using the fact that topological defects actively drive flow. In contrast to equilibrium nematic liquid crystals, disclinations of strength $\pm 1/2$ are spontaneously generated in pairs by activity. Importantly, the comet-shaped $+1/2$ defect acquires self-propulsion giomi2013defect which allows for an activity-driven defect unbinding transition to a turbulent state of spatio-temporal chaos shankar2018defect . Studies of active suspensions of microtubule bundles have also reported a remarkable state where $+1/2$ defects may themselves orientationally order in a nematic fashion on length and time scales larger than the mean free path or lifetime of an individual defect decamp2015orientational . This observation remains the subject of debate in the literature, as numerical simulations of continuum nematodynamic equations and of particle-based models predominantly report only polar defect ordering decamp2015orientational ; putzig2016instabilities ; srivastava2016negative ; patelli2019understanding or defect lattices doostmohammadi2016stabilization , whereas nematic defect order is found to be generally transient and short-lived oza2016antipolar ; srivastava2016negative . The absence of a clear physical picture for the mechanism of defect ordering has hindered a resolution of this debate. Finally, recent experiments have suggested that topological defects may serve a biological function as centers of cell extrusion or accumulation in epithelia saw2017topological ; kawaguchi2017topological , seed mound formation in bacteria yaman2018emergence , or control the morphology of the interface of growing cell colonies doostmohammadi2016defect . The nature and dynamics of active defects have therefore been the subject of intense research efforts in recent times.

In this paper we focus on topological defects and formulate a theory of two dimensional (2D) active nematics on a substrate in terms of the large scale dynamics of an interacting gas of unbound disclinations. Our work yields a complete analytical phase diagram that includes both the defect mediated melting of the active nematic and the defect ordering phase transition as a function of activity and noise. We additionally demonstrate that the defect hydrodynamic equations derived here provide a versatile framework for quantifying the behavior of systems with spatially varying activity.

In the remainder of the Introduction we motivate our approach and summarize the main results of our work.

I.1 Dynamics of Active Nematics

A now well accepted continuum description of active nematics builds on the established hydrodynamics of passive liquid crystals, augmented by nonequilibrium active stresses. This approach typically involves hydrodynamic equations for the fluid flow velocity ${\bf u}$ and the nematic order parameter ${\bf Q}$, a rank-2 symmetric traceless tensor in 2D, whose independent components are conveniently written in terms of a complex field

[TABLE]

with $Q_{xx}=-Q_{yy}=S\cos 2\theta$ and $Q_{xy}=Q_{yx}=S\sin 2\theta$. Here $S$ is the scalar order parameter and $\theta$ is the angle of the nematic director $\hat{{\bf n}}=(\cos\theta,\sin\theta)$.

Numerous numerical studies putzig2016instabilities ; srivastava2016negative ; hemingway2016correlation have shown that this continuum model is capable of reproducing the phenomenology of 2D active nematics on a substrate, including the proliferation of topological defects, as both ${\bf Q}$ and ${\bf u}$ remain well defined even in the presence of defects. On the other hand, the order parameter field and flow velocity evolve in highly complex and nonlinear ways for large activity, impeding any analytical progress. As a consequence, a theory addressing both defect chaos and defect ordering has so far proved difficult.

An alternative and fruitful strategy has been to focus on the topological defects, which both in and out of equilibrium constitute elementary, yet nontrivial, excitations of the homogeneous ordered state. This approach explicitly recognizes defects as the relevant excitations driving complex active flows and aims at developing an effective description of defects as particles, akin to the well established mapping of point-like topological defects in 2D equilibrium systems onto a Coulomb gas chaikin2000principles . In 2D nematic liquid crystals the lowest energy topological defects are strength $\pm 1/2$ point disclinations, corresponding to a distortion of the orientation where the angle $\theta$ becomes multivalued, acquiring a $\pm\pi$ jump as one encircles the respective defect, while the magnitude $S=1$ almost everywhere, vanishing only at the defect core.

The additional challenge in 2D active nematics is that the $+1/2$ disclination becomes motile, with the local geometric polarity of the defect, defined as ${\bf e}_{i}=a\bm{\nabla}\cdot{\bf Q}({\bf r}_{i}^{+})$ shankar2018defect (see Fig. 1), dictating its polarization. The $-1/2$ disclination is not self-propelled by virtue of its three-fold symmetry. As a result, an effective particle model for defects must incorporate the angular dynamics of the $+1/2$ defect polarization. Although the geometric polarity of the $+1/2$ disclination is obviously present even in passive nematics, there it remains a fast mode, rapidly relaxing on a short time scale. In contrast, when active, the dynamical nature of the polarity allows it to be a slow mode with qualitatively new physics.

Recently, working perturbatively in activity, we and collaborators mapped the dynamics of active defects onto that of a mixture of motile ($+1/2$) and passive ($-1/2$) particles with Coulomb-like interaction forces and aligning torques shankar2018defect , putting on firm ground previous purely phenomenological models giomi2013defect ; pismen2013dynamics ; keber2014topology . Using this effective model, we were able to determine the critical activity separating the quasi-ordered two-dimensional nematic from the turbulent-like state of unbound defects. This previous work demonstrates that, like in 2D equilibrium systems, treating defects as quasiparticles affords a new description of 2D ordered active media dual to the more conventional one based on order parameter fields. Far from equilibrium, where analytical progress is often limited, such techniques are very valuable and can provide insight into complex phenomena that may be hard to rationalize otherwise. In the present paper, we have adopted this viewpoint to build a unified theory for the dynamical states of 2D active nematics by considering topological defects as the relevant degrees of freedom that drive large scale flows in the system.

I.2 Results and Outline

A key new result of our work is the derivation of a hydrodynamic theory for active defects that describes defect dynamics on length scales larger than the mean defect separation or the nematic correlation length $\xi$. This is obtained by systematically coarse-graining the effective particle model for the defects derived in Ref. shankar2018defect . We then show that our hydrodynamic theory provides a complete theoretical description of defect organization in active nematics.

Defect Hydrodynamics.

The defect hydrodynamic equations are constructed in the spirit of previous classic approaches used to study the dynamic response of superfluid films ambegaokar1980dynamics , planar magnets and rotating Helium volovik1980hydrodynamics , flux liquids marchetti1990hydrodynamics and the melting of 2D crystals zippelius1980dynamics . They are formulated in terms of continuum fields, given by the number and current densities of the $+1/2$ and $-1/2$ defects defined as

[TABLE]

where ${\bf r}_{i}^{\pm}$ and $\dot{{\bf r}}_{i}^{\pm}$ are the position and velocity of the $i$th $\pm 1/2$ defect, respectively, and the phase gradient or “superfluid velocity” 111Our presentation closely follows Refs. ambegaokar1980dynamics ; zippelius1980dynamics

[TABLE]

Note that, unlike in superfluid He films, where the analogue of ${\bf v}_{n}$ genuinely represents the flow velocity of the condensate, here ${\bf v}_{n}$ is simply the gradient of the angle of the director. In this sense, it captures the distortion of the nematic and is perhaps closer in spirit to the phase gradient defined in other liquid crystalline phase, such as the hexatic zippelius1980dynamics , relevant to the study of 2D melting at equilibrium.

Now, in the presence of defects, the phase $\theta$ of the order parameter field $\Psi$ is multivalued, while ${\bf v}_{n}$ remains smooth and single-valued everywhere away from the defect cores, and provides a useful description of director deformations. Defects are continuously created and annihilated, hence the total number density of defects $n=(n_{+}+n_{-})/2$, is not conserved. On the other hand, since defects are created and annihilated only in pairs, the topological charge density $\rho=(n_{+}-n_{-})/2$ is always conserved and related to ${\bf v}_{n}$ through the important topological constraint

[TABLE]

that allows an analogy with superfluid hydrodynamics ambegaokar1980dynamics and electrostatics, through Gauss’ law Jackson1975 . An important distinguishing feature of the 2D active nematic is that that the geometric polarization of the $+1/2$ disclination provides a new dynamically relevant internal degree of freedom. This requires the addition of a new hydrodynamic field, the defect polarization density, defined as

[TABLE]

The polarization directs the self-propulsion of the $+1/2$ disclinations through self-induced active backflows. In addition, it experiences active torques that reorient the defect in response to the elastic forces from other defects. Both these properties underlie much of the phenomena explored here.

The defect hydrodynamic equations are presented in Sec. II and their derivation is shown in Appendix A.

Isotropic defect chaos.

Upon analyzing the steady states of the defect hydrodynamic equations and their stability (Secs. III and IV), we find two stable states with a finite number of unbound defect. The first is the spatio-temporally chaotic state referred to as “active turbulence” in the literature doostmohammadi2018active . In this state the $+1/2$ defects, although motile, have no preferred direction of polarization. We refer to this state as “isotropic defect chaos”, where ‘isotropic’ refers here to the fact that the gas of $+1/2$ defects has zero mean polarization ($\langle{\bf p}\rangle={\bf 0}$). Our equations validate previous phenomenological scaling hypotheses and provide a well founded theoretical description of defect chaos in the active turbulent regime by extending the mean-field approach of Ref. giomi2015geometry to large scales. Many-body screening crucially impacts the decay of the conserved defect charge density due to the presence of the topological constraint given by Eq. 5, unique to defect hydrodynamics. As the decay rate is governed by the density $n$ of free disclinations, we find that the dominant length scale of correlations in nematic order or velocity are all controlled mainly by the mean defect spacing $\xi\sim 1/\sqrt{n}$, in agreement with previous numerical results giomi2015geometry ; hemingway2016correlation .

Polar defect order.

For higher activity, our equations predict a second stable state, where the motile $+1/2$ defects themselves order. This occurs via a continuous transition in which the active $+1/2$ disclinations collectively align and condense into a liquid with long-range polar order. At the same time, the underlying nematic develops a periodic modulation of kink walls, as shown in the sketch of Fig. 2. This state, which we refer to as “polar defect order”, has no giant fluctuations in either the defect charge or number density and provides an intriguing realization of a “Malthusian defect flock”. While polar defect order has been reported before in numerical models of active nematics decamp2015orientational ; putzig2016instabilities ; srivastava2016negative ; patelli2019understanding , the mechanism driving it has remained unexplained. Our work identifies a mechanism for polar order as arising from both active self-aligning torques (derived perturbatively in activity in Ref. shankar2018defect ) and many-body screening. Overall, our result is a complete analytical description of 2D active nematics that includes both the defect mediated melting of the active nematic and the defect ordering phase transition as a function of activity and noise. The various states and transitions are summarized in the phase diagram shown in Fig. 3 in terms of parameters of the coarse-grained theory. To make contact with possible experiments, we have also reformulated the phase diagram in Fig. 6 in terms of activity and the liquid crystal stiffness, quantities that are both directly accessible and tunable in experiments. We note here that nematic order of defects as reported in the experiments of Ref. decamp2015orientational cannot appear in our model via a transition from the isotropic defect chaos state. The reasons for this are discussed in Sec. III.2 (see also Fig. 4) and further elaborated in Appendix B using an analogy with electrostatics. In Sec. VI, we expand on future directions and current challenges, including the theoretical possibilities (or lack thereof) for the existence of apolar defect order.

Spatially-varying activity and defect trapping.

Finally, in Sec. V we demonstrate the versatility of our hydrodynamic approach by employing it to describe defect dynamics in systems with spatially inhomogeneous activity. A simple motif we study is an active-passive interface. Due to the self-propelled nature of the $+1/2$ disclinations, they are found to accumulate on the passive (low activity) side of the interface. The consequent charge segregation and local polarization at the interface indicates that an activity gradient can be thought of as a local “electric field” driving charge sorting. An extension of the same phenomenon is also realized at an extensile-contractile interface across which activity changes sign. Such basic principles can be combined with more complicated activity patterns to position and move defects in a programmable fashion. The use of activity gradients to control and guide defect dynamics is an important technique to develop active microfluidic devices with targeted transport capabilities. The spatio-temporal modulation of activity with light schuppler2016boundaries ; ross2018controlling is a very promising approach in functionalizing active matter to engineer new metamaterials. We expect our theoretical results and the proposed hydrodynamic framework to be useful tools in predicting the collective behavior of active defects in both such inhomogeneous backgrounds and complex geometries.

II Active defect hydrodynamics

We consider a 2D active nematic on a substrate with a finite concentration of unbound disclinations. Defects are unbound in pairs to maintain charge neutrality, hence the system contains an equal number of $+1/2$ and $-1/2$ defects. In Ref. shankar2018defect , the full nematodynamic equations for an active nematic on a substrate were recast into an effective model for defects as interacting quasiparticles. While the particle model was explicitly derived only perturbatively for small activity, its general structure is expected to survive for large activity as well. Different models for active defect dynamics of varying complexity have been proposed by others as well cortese2018pair ; tang2019theory but the basic qualitative features remain the same. As the defect equations of motion presented in Ref. shankar2018defect are the easiest to interpret and most amenable to direct coarse-graining, we proceed with them to compute the required transport coefficients. The details of the calculation are given in Appendix A. We shall present the final equations directly here.

The number of $+1/2$ or $-1/2$ defects can change through pair creation or annihilation events, hence the individual number densities $n_{\pm}$ evolve according to

[TABLE]

where $W_{c}$ and $W_{a}$ are the rates of defect creation and annihilation, respectively. Both $W_{c}$ and $W_{a}$ depend on activity and $n_{\pm}$; we refrain from specifying their explicit form which only serves as an input to the hydrodynamic theory. The constitutive relation for the defect currents ${\bf j}_{\pm}$ is (see Appendix A for derivation)

[TABLE]

where $\bm{\epsilon}$ is the 2D Levi-Civita tensor, $\nu=\pi\mu\gamma$ a dimensionless number involving the defect mobility $\mu$ and the rotational viscosity $\gamma$, and $D_{0}=\mu T$ is the bare defect diffusion constant ($T$ is the corresponding effective temperature). Liquid crystal elasticity controls the nematic diffusion constant $\kappa=K/\gamma$, where $K$ is a Frank elastic constant. Activity is encoded in the self-propulsion speed $|v|$ of the $+1/2$ defect. Note that $v$ can be of either sign depending on the nature of active stresses in the medium. Extensile systems have $v<0$ and contractile systems have $v>0$. Even in the passive limit ($v=0$), defects move transverse to the local phase gradient ${\bf v}_{n}$. This response is akin to the Magnus force on vortices lamb1993hydrodynamics or the Peach-Koehler force on dislocations peach1950forces .

Defect hydrodynamics crucially differs from conventional long-wavelength hydrodynamics due to the presence of a topological constraint. For an arbitrary closed curve $\Gamma$ enclosing $N$ unbound disclinations carrying charges $q_{i}=\pm 1/2$ located at positions ${\bf r}_{i}^{\pm}$, the net accumulated director phase is given by the line integral

[TABLE]

Using Stokes’ theorem and Eq. 2, this gives the topological constraint given in Eq. 5. Using Eqs. 7, 8 and 9, we obtain equations for the charge and number densities $\rho$ and $n$ as

[TABLE]

Finally, the dynamics of the phase gradient ${\bf v}_{n}$ and of the defect polarization ${\bf p}$ is derived by coarse-graining (see Appendix A) to give

[TABLE]

In Eq. 15, the defect charge current ${\bf j}_{\rho}$ explicitly breaks the “conservation” of the phase gradient ${\bf v}_{n}$. This, as we shall see later, in turn also causes the charge density $\rho$ to relax on a finite time-scale notwithstanding the local conservation of topological charge (Eq. 11). The second term on the right hand side of Eq. 15 describes the relaxation of smooth director deformations due to liquid crystal elasticity. For the polarization dynamics (Eq. 16) we have adopted a simple isotropic closure for the second moment while neglecting all higher order correlations (see Appendix A). This is sufficient for our purposes of demonstrating defect ordering. The relaxation rate of the polarization is set by the rotational diffusion $D_{R}$ of the $+1/2$ defect and it receives nonlinear passive corrections proportional to $\rho$ and $|{\bf v}_{n}|^{2}$ with $\beta=5\nu^{2}\kappa/8$, as described by the last two terms in square brackets. Note that, although the defect charge density $\rho$ can be locally negative, this doesn’t lead to an instability in the polarization equation at equilibrium, because $\langle\rho\rangle=0$ in a charge neutral system and deviations from this mean value relax on a finite timescale. The convective-like term $\hat{{\bf z}}\cdot({\bf v}_{n}\times\bm{\nabla}){\bf p}$ accounts for the change in the polarization due to the passive motion of $+1/2$ defects and the penultimate term on the right hand side of Eq. 16 is a passive elastic torque that rotates the defect polarization in response to elastic distortions (see also Appendix A.3). Activity enters in two places: in the pressure like term $\sim(v_{1}/2)\bm{\nabla}n_{+}$ and as an active torque $\sim v_{2}n_{+}\,\bm{\epsilon}\cdot{\bf v}_{n}$ that builds up local polar order due to transverse director deformations. The active coefficients $v_{1}$ and $v_{2}$ are defined as

[TABLE]

where $\zeta$ is a nondimensional activity. The positive dimensionless function $f(\zeta)=1+4\zeta+\mathcal{O}(\zeta^{2})$ captures the nonlinear dependence of the torque on activity and can be computed from the moment hierarchy (see Appendix A). Below we shall analyze the steady states of these equations and their stabilty.

III Homogeneous states and phase transitions

We consider a state with a homogeneous density of unbound disclinations. By charge neutrality in the plane and the topological constraint (Eq. 5), we have $\rho=0$ and hence $n_{+}=n_{-}=n$. Setting all the gradient terms to zero, we obtain

[TABLE]

At steady state, we set $n=n_{0}$ such that $W_{c}=W_{a}$, ${\bf v}_{n}={\bf v}_{n}^{0}$ and ${\bf p}={\bf p}_{0}$. From Eq. 20 we have ${\bf v}_{n}^{0}=(v/2\nu\kappa n_{0})\bm{\epsilon}\cdot{\bf p}_{0}$, which corresponds to a vanishing charge current (${\bf j}_{\rho}={\bf 0}$). Eliminating ${\bf v}_{n}^{0}$ from Eq. 21 and after some algebraic manipulations we get

[TABLE]

where

[TABLE]

Note that $a_{2}$ can change sign at high activity ($\zeta\sim v^{2}$) and $a_{4}>0$. For small activity $a_{2}<0$ and the only solution is ${\bf p}_{0}={\bf 0}$ and ${\bf v}_{n}^{0}={\bf 0}$. At large activity $a_{2}>0$ and we obtain a solution with $|{\bf p}_{0}|\neq 0$ and therefore $|{\bf v}_{n}^{0}|\neq 0$, corresponding to a uniformly polarized defect ordered state. The change in sign of $a_{2}$ occurs at a critical value of activity set by

[TABLE]

where we have used the leading order expansion of $f(\zeta)=1+4\zeta+\mathcal{O}(\zeta^{2})$. Including higher order terms in the expansion of $f(\zeta)$ will change the numerical value of $\zeta_{c}$, which will, however, remain finite. Altogether we have three distinct homogeneous steady states:

• •

Homogeneous nematic order with no unbound disclinations ($n_{0}=0$, ${\bf v}_{n}^{0}={\bf p}_{0}={\bf 0}$).

• •

Isotropic defect chaos ($n_{0}>0$, ${\bf v}_{n}^{0}={\bf p}_{0}={\bf 0}$).

• •

Defect polar ordered state ($n_{0}>0$, ${\bf v}_{n}^{0},\,{\bf p}_{0}\neq{\bf 0}$).

The three phases along with the intervening phase boundaries discussed below are shown in the phase diagram in Fig. 3, as a function of $v$, $T$ and $D_{R}$, keeping $\nu$ fixed. Changing $\nu$ does not change the global topology of the phase diagram, but affects the phase boundaries only in a quantitative way.

III.1 Defect unbinding transition

The nematic order-disorder transition is mediated by an activity driven unbinding of defect pairs shankar2018defect . The resulting defect-ridden state is isotropic and disordered (${\bf p}_{0}={\bf 0}$) with a finite nematic correlation length $\xi\sim 1/\sqrt{n_{0}}$. Given the disordered motion of the $+1/2$ disclinations, we identify this state with the spatio-temporally chaotic dynamics of “active turbulence”. The activity threshold for defect unbinding was obtained in Ref. shankar2018defect to be

[TABLE]

where $\tilde{T}=T/T_{c}^{\mathrm{eq}}$ is a normalized effective temperature with $T_{c}^{\mathrm{eq}}=\pi K/8$ the equilibrium Kosterlitz-Thouless transition temperature kosterlitz1973ordering ; stein1978kosterlitz . This transition line is marked in blue in Fig. 3.

The location of the defect unbinding transition can be understood by a simple argument first given in Ref. shankar2018defect and repeated here for completeness. In a naïve one-dimensional picture, where the two defects of a neutral pair ($\pm 1/2$) unbind by moving away from each other along a straight line, the self propulsion of the $+1/2$ defect can always overcome the passive Coulomb attraction, resulting in defect unbinding at any activity. On the other hand, rotational diffusion ($D_{R}$) can spoil this process by endowing the $+1/2$ defect with a finite persistence length $\ell_{p}=|v|/D_{R}$. When $\ell_{p}<r_{c}$, with $r_{c}\sim\mu K/|v|$ being the pair separation where the propulsive force $|v|/\mu$ and the attractive Coulomb forces $K/r$ balance, rotational noise disrupts the straight path of the $+1/2$ defect before it can overcome the energy barrier required for unbinding, allowing defect pairs to remain bound. The condition for unbinding due to the activity can then be estimated as $\ell_{p}\sim r_{c}$, which coincides with Eq. 25 for small $T$.

III.2 Defect ordering transition

For $\zeta>\zeta_{c}$ as given by Eq. 24 the isotropic gas of defects spontaneously breaks rotational symmetry by ordering into a polar, collectively moving liquid. Expressing the transition point in terms of the original model parameters, we find that the defect ordering transition occurs at

[TABLE]

which is shown as a red line in Fig. 3.

A simple yet physical way to understand this threshold for defect ordering is as follows. Disregarding the numerical constant $\zeta_{c}$, the condition for defect ordering ($|v|>|v_{c_{2}}|$) can be written as

[TABLE]

The factor $T/K$ accounts for the relative strength of fluctuations $\sim T$ to the defect core energy $\sim K$ chaikin2000principles , and it corresponds to the cost to nucleate a defect pair. Once created, the $+1/2$ disclination self-propels itself away from the $-1/2$ disclination, quasi-deterministically on time scales shorter than the rotational diffusion time $D_{R}^{-1}$. In doing so, it distorts the underlying nematic in its wake, on the scale of a defect persistence length $\ell_{p}$. As $\gamma\ell_{p}^{2}/K$ is the time it takes the underlying nematic to relax distortions on a length scale of $\ell_{p}$, the threshold for defect ordering is a simple balance of time scales.

When the nematic rapidly heals the distortion created by the swarming $+1/2$ disclination, faster than the defect reorientation time ($\gamma\ell_{p}^{2}/K\ll D_{R}^{-1}$), we obtain an isotropic disordered state of defect chaos. In the opposite limit, the underlying nematic responds too slowly and is unable to relax the distortion left by the active defects. This leads to the formation of a long-lived kink wall (also called a $\pi$- or Néel-wall) that terminates at the defect pair. This locally frozen-in distortion feeds back into the defect motion, leading to a build up of polar order through a combination of many-body screening and active torques. Eventually, for strong enough activity, the underlying nematic cannot catch up with the persistent dynamics of the defects, which then condense into a spontaneously flowing defect polarized liquid, i.e., a $+1/2$ defect flock. The active self-aligning torques derived in Ref. shankar2018defect are a crucial ingredient to this mechanism of defect ordering. Reminding ourselves of the physics of the active torque, we note that it arises from the active backflow induced defect motion advecting the defect polarization itself. This results in an effective torque which tries to align the defect orientation to its velocity, with a magnitude controlled directly by activity. While similar torques with short ranged interaction forces have been used to model flocking in both cells szabo2006phase ; henkes2011active and vibrated polar grains weber2013long ; lam2015self , it is only in the presence of many-body screening that the active torques can cause collective motion of defects in active nematics.

Close to the defect ordering transition, for $\zeta>\zeta_{c}$ we have

[TABLE]

with the usual mean-field exponent, although fluctuations are expected to decrease it. As the phase gradient ${\bf v}_{n}=\bm{\nabla}\theta$ is slaved to the polar order, we also have

[TABLE]

with ${\bf p}_{0}\cdot{\bf v}_{n}^{0}=0$. The appearance of polar defect order also spontaneously breaks translational symmetry of the underlying nematic in the direction orthogonal to that of defect order, with $\theta({\bf r})\simeq\theta_{0}+{\bf v}_{n}^{0}\cdot{\bf r}$. While the defect liquid itself has no translational order, the underlying nematic undergoes a concomitant modulational instability and develops a smectic array of splay-bend kink walls, similar to that seen in numerical simulations decamp2015orientational ; putzig2016instabilities ; srivastava2016negative ; patelli2019understanding . A cartoon of the defect and kink wall structure in the defect ordered state is sketched in Fig. 2.

Finally, in light of the apolar nature of an active nematic, the appearance of a directed polar current might seem surprising. The polarization density is an emergent property of nonlinear topological excitations in the system which then permits spontaneous flow due to the absence of detailed balance. On the other hand, experiments in microtubule suspensions have reported apolar (i.e., nematic) ordering of $+1/2$ disclinations decamp2015orientational that may seem a more natural possibility. Nematic defect order is not possible in our model due to the nature of defect interactions, particularly the active torques. To see this, consider a test $+1/2$ disclination placed in the vicinity of a neutral defect pair. The $\pm 1/2$ defect pair instantaneously creates a background phase gradient ${\bf v}_{n}$ orthogonal to the line joining the two defect cores as shown in Fig. 4. When a test $+1/2$ disclination is placed in this background distortion, one immediately sees that the active torque preferentially stabilizes its orientation to be aligned with the $+1/2$ defect in the neutral pair. This is illustrated in Fig. 4 for the case of extensile activity ($v<0$) and is discussed in more detail in Appendix B. The contractile case works analogously. Passive elastic torques $\sim\kappa\bm{\nabla}\cdot{\bf v}_{n}$ cannot change this as they vanish for spatially homogeneous phase gradients, unlike the active torques. Hence active torques always induce an effective polar alignment between $+1/2$ defects, thereby structurally preventing the appearance of nematic defect order. This should be contrasted with other examples of defect order in equilibrium - the Abrikosov vortex lattice in type II superconducting films abrikosov1957magnetic and the twist-grain boundary phase in smectics renn1988abrikosov , both of which result from the local breakdown of Meissner like effects, leading to the penetration of either the magnetic field or twist, respectively.

IV Fluctuations and linear stability

We shall now examine the linear stability of the homogeneous steady states to small spatial fluctuations. In both the isotropic and defect ordered states, fluctuations in the average defect number density ($n$) relax on a short time scale set by the balance of defect pair creation and annihilation and are henceforth neglected. In the presence of unbound defects, the nematic order parameter has a finite correlation length $\xi$ which sets the mean separation between defects. Fixing $n=n_{0}$ we define

[TABLE]

The numerical factor of $2\pi\nu$ is introduced to simplify the notation below.

IV.1 Isotropic defect chaos

Linearizing for small fluctuations about the isotropic steady state, we have $\rho=0+\delta\rho$, ${\bf v}_{n}={\bf 0}+\delta{\bf v}_{n}$, ${\bf p}={\bf 0}+\delta{\bf p}$. Upon Fourier transforming in space ($\Phi_{{\bf q},\omega}=\int\mathrm{d}^{2}r\int\mathrm{d}t\ e^{i\omega t-i{\bf q}\cdot{\bf r}}\Phi({\bf r},t)$), we define the longitudinal and transverse components of $\delta{\bf p}$ as $\delta p_{L}=\hat{{\bf q}}\cdot\delta{\bf p}_{{\bf q}}$ and $\delta p_{T}=\hat{{\bf z}}\cdot\left(\hat{{\bf q}}\times\delta{\bf p}_{{\bf q}}\right)$ respectively. A similar decomposition is performed for $\delta{\bf v}_{n}$ as well. The resulting linearized equations are

[TABLE]

Note that the four fields are only coupled in pairs; $(\delta p_{L},\,v_{n}^{T})$ and $(\delta p_{T},\,v_{n}^{L})$. Here we have used the topological constraint $\hat{{\bf z}}\cdot\left(\bm{\nabla}\times\delta{\bf v}_{n}\right)=2\pi\delta\rho$ and introduced $\tau=(2\pi\nu\kappa n_{0})^{-1}=\xi^{2}/\kappa$ as the finite relaxation time of the phase gradient. The thermal diffusion length scale $\ell^{2}=D_{0}/D_{R}$ is typically microscopic and expected to be of the order of the defect core size ($\ell\sim a$). Although the charge density $\rho$ is locally conserved (Eq. 11), its fluctuations are slaved entirely to $\delta\mathrm{v}_{n}^{T}$ which decays on a finite time scale $\tau$, due to complete screening of the long-ranged Coulomb (passive) interaction between the defects.

Both pairs of coupled modes have the same dispersion relation at $\mathcal{O}(q^{0})$, given by

[TABLE]

differing from $\mathcal{O}(q^{2})$ terms onwards. While $i\omega_{+}>0$ always, $i\omega_{-}$ can go negative for large enough activity ($\zeta$) triggering an instability of the isotropic state. This instability coincides with the defect ordering transition at $\zeta=\zeta_{c}$ (obtained in Eq. 24). One can check that ${\bf q}$-dependent terms do not change this and no further instabilities arise.

As we approach the defect unbinding transition from above, within the defect chaos state, the nematic correlation length diverges ($\xi\to\infty$) and the density of free defects vanishes ($n_{0}\to 0$). Polarization fluctuations relax with a finite rate due to rotational diffusion, even at the unbinding transition ($i\omega_{+}\simeq D_{R}$), while the phase gradient exhibits critical slowing down with $i\omega_{-}\propto\xi^{-2}\to 0$. Of course the region of hydrodynamic validity $q\xi\ll 1$ shrinks rapidly as we approach the defect unbinding transition.

Reinstating the lowest order additive noise as computed from coarse-graining ($2\pi\sqrt{D_{0}n_{0}}\,\bm{\Lambda}_{1}$ in Eq. 15 and $\sqrt{2TD_{R}n_{0}f(\zeta)/K}\,\bm{\Lambda}_{2}$ in Eq. 16, with both $\bm{\Lambda}_{1}$ and $\bm{\Lambda}_{2}$ uncorrelated unit white Gaussian noise), we compute the fluctuation spectra in the isotropic defect gas. Using $q=|{\bf q}|$, the equal time correlator of the polarization at steady state is given by,

[TABLE]

As we are far from the defect ordering threshold ($\zeta\ll\zeta_{c}$), we have only retained the most dominant terms involving $\delta p_{T}$ and neglected $\ell\ll\xi$ for simplicity. In this limit, $\delta p_{L}$ has correlations only on scales much smaller than $\xi$ and mainly contributes an additive constant in Eq. 36. Note the dependence on $n_{0}$ that arises from the variance of the noise in the polarization dynamics. The polarization inherits its dominant spatial correlations from ${\bf v}_{n}$, which survives on the scale of $\xi$. Note that this results in a breakdown of dynamic scaling close to the defect unbinding transition. The charge density fluctuations for ${\bf q}\to{\bf 0}$ are similarly given by

[TABLE]

upon assuming $\zeta\ll\zeta_{c}$, well below the defect ordering threshold. Charge density fluctuations vanish on large scales as a consequence of the unbound defect gas behaving as a screened conducting plasma albeit with an activity enhanced effective “dielectric constant” via a sum rule martin1988sum . While such a dielectric constant relates to an effective elastic constant or stiffness in 2D equilibrium superfluids or XY magnets, the equivalent correspondence in the active case requires care and will be addressed elsewhere shankar2020 .

A common feature of defect chaos is the characteristic decay of the flow velocity correlator and the kinetic energy spectrum. We assume that the defects serve as faithful tracers of the average flow field and write

[TABLE]

where the last approximate equality is obtained by linearizing ${\bf j}_{n}$ for small $\delta{\bf p}$ and $\delta{\bf v}_{n}$ and neglecting higher order gradient terms. This approximation of course breaks down close to the defect unbinding transition when $n_{0}\to 0$.

Two distinguishing properties are generally associated with “active turbulence”. The first is that the typical length scale of the flow depends on the magnitude of activity $|\alpha|$ through the mean defect spacing $\xi$, with $\xi\sim|\alpha|^{-1/2}$ sanchez2012spontaneous ; thampi2013velocity ; giomi2015geometry ; hemingway2016correlation ; lemma2018statistical . The second is the scaling of the kinetic energy spectrum $E({\bf q})=\langle|{\bf u}_{{\bf q}}|^{2}\rangle/2$, with $E(q\to 0)\propto|\alpha|$. In our work $|v|\sim|\alpha|$ (see Appendix A) and $\xi\sim n_{0}^{-1/2}$ away from the defect unbinding threshold. Deep in the regime of defect chaos, the mean density of defects is known to scale with the magnitude of activity giomi2015geometry ; hemingway2016correlation , hence $\xi\sim|v|^{-1/2}\sim|\alpha|^{-1/2}$. To derive this scaling from semi-microscopic considerations requires a full theory of melting, far beyond the scope of this paper, and will be presented elsewhere shankar2020 . As can be seen directly from Eq. 38 and Eq. 36, in our model the characteristic length scale of both flow and defect polarization indeed scales as $\xi$ far from the defect unbinding threshold. For the second point, we can obtain the scaling of the kinetic energy using Eq. 38 and 36, with the result $E(q\to 0)\propto|v|(1+\mathcal{O}(|v|))$. The leading $|v|$ scaling of the average kinetic energy crucially relies on the fact that the polarization fluctuations $\propto n_{0}$ (Eq. 36), a consequence of the central limit theorem. Both these scalings no longer hold near the defect unbinding transition. It is worth emphasizing that our calculation based on a systematic derivation of the defect dynamics is consistent with the phenomenological mean-field picture of Ref. giomi2015geometry . The two approaches are complementary, with the mean-field construction working best on short scales where correlations are neglected, while our approach works best on large scales where hydrodynamic treatments are applicable.

IV.2 Defect polar order

Deep in the defect ordered state, both ${\bf v}_{n}$ and $\rho$ are fast modes and rapidly relax to their steady-state values (${\bf v}_{n}^{0}$ and [math] respectively) on the now short time scale $\tau$. It is instructive (although not essential) to slave both these fields to the polarization. In doing so, we set ${\bf v}_{n}=\pi v\tau\,\bm{\epsilon}\cdot{\bf p}+\mathcal{O}(\bm{\nabla})$ and $\rho=-(v\tau/2)\bm{\nabla}\cdot{\bf p}+\mathcal{O}(\nabla^{2})$ to obtain an effective Toner-Tu like equation for the defect polarization,

[TABLE]

The resulting hydrodynamic coefficients are

[TABLE]

with $a_{2},a_{4}$ given in Eq. 23. Note that, to leading order as $n_{0}\propto|v|$ ($|v|\sim|\alpha|$), we find that $\lambda_{1,2,3}\sim\mathrm{sgn}(v)$ is nonanalytic in activity, and the effective splay elastic constant $D_{1}\sim|v|/K$ is controlled by the well-known active length scale doostmohammadi2018active . As noted earlier in Sec. III.2, for $\zeta>\zeta_{c}\simeq 0.701$, $a_{2}>0$ and the isotropic defect gas spontaneously orders into a polarized liquid. The finite relaxation time for the defect charge density $\rho$ implies that the defect ordered liquid behaves as a Malthusian flock toner2012birth . From Eq. 39, we also see that the important convective nonlinearity $\lambda_{1}{\bf p}\cdot\bm{\nabla}{\bf p}$ naturally appears in our framework and it allows the existence of long-ranged polar order of the $+1/2$ disclinations. The spontaneous breaking of rotational symmetry is accompanied by a simultaneous breaking of translational symmetry in the underlying nematic, characterized by the appearance of periodic kink walls with a nonzero average phase gradient $\langle{\bf v}_{n}\rangle={\bf v}_{n}^{0}$. In Fig. 2, we see that the $+1/2$ defects preferentially move along the kink walls, suggesting a structure akin to the “active smectic-P” state recently reported in Ref. romanczuk2016emergent . This superficial similarity goes no further as neither the defect number nor charge density reflect the necessary periodic modulation. This is an important distinction, since the periodic arrangement of kink walls and associated long-ranged smectic order (unlike all other 2D active smectics adhyapak2013live ; chen2013universality ; romanczuk2016emergent ) is not an independent broken symmetry, but rather the result of the LRO of the polarization. The latter in turn arises as in all Toner-Tu models from the convective nonlinearity in Eq. (39).

Linearizing for small fluctuations about the defect ordered steady state, we only have ${\bf p}={\bf p}_{0}+\delta{\bf p}$ as fluctuations in both ${\bf v}_{n}$ and $\rho$ have already been enslaved. Decomposing $\delta{\bf p}$ along ($\delta p_{||}$) and transverse ($\delta p_{\perp}$) to the polar order, we obtain two coupled modes whose long-wavelength dispersion relation is of the form

[TABLE]

where $q_{||}={\bf q}\cdot\hat{{\bf p}}_{0}$. Both amplitude ($\delta p_{||}$) and orientational ($\delta p_{\perp}$) fluctuations propagate along the direction of polar order, though with different speeds,

[TABLE]

having used the magnitude of polarization $|{\bf p}_{0}|=p_{0}$. Depending on the activity ($v$), these drift speeds can be of either sign (see Eqs. 40, 41). The relaxation rates for the two modes upto $\mathcal{O}(q^{2})$ are given as

[TABLE]

where $\varphi$ is the angle made by the wave-vector ${\bf q}$ with the direction of polar order ($\hat{{\bf q}}\cdot\hat{{\bf p}}_{0}=\cos\varphi$). As expected, in the defect ordered state, the amplitude mode decays with a finite relaxation rate even as ${\bf q}\to{\bf 0}$, while orientational fluctuations remain soft. Being a broken symmetry variable, this is the only true hydrodynamic mode here. From our derivation of transport coefficients (Eqs. 40, 41), we note that $\lambda_{2}\lambda_{3}<0$ for small $v$, changing sign for larger activity. This is true for both extensile and contractile systems. Hence, while the amplitude mode is always stable ($\Gamma_{||}({\bf q})>0$) 222Though, from Eq. 44, it might seem that $\delta p_{||}$ can go unstable at high activity when $\lambda_{2}\lambda_{3}>0$, this is an artifact of the low ${\bf q}$ expansion. In fact, for $\varphi=\pi/2$, the complete dispersion relation for the amplitude mode is $i\omega_{||}=[\tilde{a}_{4}+(2D_{0}+D_{1})q^{2}+\sqrt{(\tilde{a}_{4}-D_{1}q^{2})^{2}-4q^{2}\tilde{a}_{4}\Lambda}]/2$, where $\tilde{a}_{4}=2a_{4}D_{R}p_{0}^{2}$ and $\Lambda=\lambda_{2}\lambda_{3}/a_{4}D_{R}$. For either sign of $\Lambda$, $\mathrm{Re}(i\omega_{||})>0$ for all $q$, confirming that the only possible instability is the orientational one obtained in Eq. 46 the orientational mode $\delta p_{\perp}$ can develop a long-wavelength splay instability (at $\varphi=\pi/2$) if $\Gamma_{\perp}({\bf q})<0$. Rewriting this condition in terms of our original control parameters, this corresponds to

[TABLE]

for a splay instability to occur. The detailed dependence of $n_{0}$ on system parameters is sensitive to microscopic details, so this instability threshold is model dependent. As $\zeta\propto v^{2}$ and $n_{0}\propto|v|$ along with $f(\zeta)$ monotonic and positive, for large enough activity, we expect the above condition to not be satisfied and thereby allow a stable polar ordered phase of defects. In this paper we only consider the stable situation and shall not discuss the splay instability further.

Putting back noise as before, we can compute the correlation function of the polarization fluctuations. The convective term is a relevant nonlinearity in 2D that dramatically modifies the scaling of the autocorrelation, and is well known for being responsible in stabilizing long-ranged polar order. Using the exact exponents calculated in Ref. toner2012birth , we have at steady-state

[TABLE]

As before we have used $r_{||}=\hat{{\bf p}}_{0}\cdot{\bf r}$ and $r_{\perp}=\hat{{\bf z}}\cdot\left(\hat{{\bf p}}_{0}\times{\bf r}\right)$, and $G(x,y)$ is a scaling function that depends on model parameters in a complicated way. Both defect charge and number density exhibit normal Poissonian fluctuations, as expected of a Malthusian flock. The defect ordered liquid continually turns over due to spontaneous pair creation and annihilation of defects, with polar order persisting for infinitely longer than the finite lifetime of individual disclinations. The conspicuous absence of both giant number fluctuations ramaswamy2003active ; narayan2007long ; marchetti2013hydrodynamics and motility-induced phase separation phenomenology tailleur2008statistical ; fily2012athermal ; cates2015motility for defects is at first glance surprising. The local conservation of topological charge, along with the self-propulsion of $+1/2$ disclinations might lead one to naïvely expect large fluctuations in the charge density upon defect ordering. In addition, the repulsive interaction between motile $+1/2$ disclinations might also raise the question of the possibility of phase separation of a defect liquid. None of these scenarios are realized. The fundamental reason is the nature of the screened Coulomb interaction between defects mediated by Frank elasticity, which renders the charge density a non-hydrodynamic field with a finite relaxation time.

This concludes the analysis of the phase diagram predicted by our model of defect hydrodynamics. Untill now, all the phases considered had no charge separation and the average defect charge current always vanished ($\langle{\bf j}_{\rho}\rangle={\bf 0}$). In 2D passive systems such as superfluid films, a charge current can be generated by the application of an external electric field. While applying a similar external field in active nematics might be tricky, we will show below that spatially varying activity can locally act as an “electric field” causing a local sorting of defects based on their topological charge.

V Inhomogeneous activity

We now demonstrate that our hydrodynamic model provides a useful framework for describing situations with spatially inhomogeneous profiles of activity. Activity gradients provide additional nonequilibrium driving forces that can be used to generate and control spatial patterns at will. Recent advances in engineering optical control of biomolecular or catalytic activity have emerged as a venerable platform for dynamically creating large scale reconfigurable patterns and structures in diverse systems ranging from gels of biofilament-motor complexes schuppler2016boundaries ; ross2018controlling to suspensions of self-propelled colloids palacci2013living and bacteria arlt2018painting ; frangipane2018dynamic . Such techniques, when used in active nematics, can control the chaotic dynamics and provide novel ways to precisely sculpt flow and tune material properties through the spatial organization and positioning of topological defects. This is an important pre-requisite for developing programmable active metamaterials whose local organization dictates its global response and transport.

To include spatially inhomogeneous activity in our defect hydrodynamics, we make use of the fact that activity enters through the defect self-propulsion $v$ and simply replace $v\to v({\bf r})$ (including it under gradients when appropriate). While modulating the microscopic active processes will in general also affect the elasticity and viscosity of the nematic, we shall not consider such modifications for simplicity. There are two distinct phenomena that result when activity is spatially varying in an active nematic. First, activity controls the motility of $+1/2$ disclinations, which behave as active particles that are known to aggregate where they move slowly schnitzer1993theory ; tailleur2008statistical ; cates2015motility . This well understood phenomenon survives even in the presence of interactions and has been used to design self-assembled rectification devices stenhammar2016light and particle traps magiera2015trapping ; sharma2017brownian ; grauer2018spontaneous . The second phenomenon is a distinct property of active nematics, which is that higher activity generates more defect pairs. Hence, low activity leads to low defect motility and a consequent accumulation of $+1/2$ defects, but it also decreases the total defect density. As the two competing effects don’t act symmetrically on both charge defects, we have the possibility of sorting defects by charge in the presence of an activity gradient.

Working for simplicity in a one-dimensional (1D) setting, we consider activity to only vary in the $x$-direction ($v(x)$). We will take the maximum value of activity to correspond to states deep in the regime of defect chaos and to never exceed the defect ordering threshold. As in Sec. IV.1, we assume that the average defect density relaxes to its steady state value on a short time scale 333This neglects changes in $n$ due to a defect flux $\bm{\nabla}\cdot{\bf j}_{n}$, but it can be accounted for easily., and all other fields are small with weak gradients, permitting a linearized analysis. This is done simply for analytical progress - alternatively the equations could be solved numerically. The dependence of $W_{c}$ and $W_{a}$ on $v(x)$ causes the average defect number density to also have a spatial profile $n=n_{0}(x)\propto|v(x)|$. Due to the 1D setup, only ${\bf p}=p(x)\hat{{\bf x}}$ and ${\bf v}_{n}=\mathrm{v}_{n}(x)\hat{{\bf y}}$ are nonvanishing, and at steady state, once again ${\bf j}_{\rho}={\bf 0}$. The linearized steady state equations to leading order in gradients are

[TABLE]

Note that $v_{1}(x)$ and $v_{2}(x)$ vary in space through their dependence on $v(x)$ (see Eq. 17). Eliminating $p(x)$, we obtain a finite $\mathrm{v}_{n}(x)$ controlled by an activity gradient. The topological constraint gives $\rho(x)=\partial_{x}\mathrm{v}_{n}(x)/(2\pi)$, which in analogy with Gauss’ law demonstrates that a transverse phase gradient, here set up by a gradient in activity, acts as a local “electric field”. So we have

[TABLE]

One can check that $\rho(x)$ does not depend on the sign of $v(x)$ whereas $p(x)$ does, as expected. The denominators in the expressions above don’t vanish as we are below the defect ordering transition. The defect charge and number density for two different activity profiles are plotted in Fig. 5, along with a schematic showing the polarization of defects at the active-passive interface. Hence an active-passive interface is both charge and orientationally polarized. As seen in Fig. 5, there is an excess of $+1/2$ disclinations on the passive side of such an interface and charge balanced excess of $-1/2$ disclinations on the active side. This can be understood by recalling that the $+1/2$ defects being motile tend to accumulate where they move slowly. Similarly, an active region flanked by passive regions on either side develops a net negative topological charge concentrated near the interfaces 444A related phenomenon is the spontaneous charging of isotropic tactoids in an active nematic genkin2018spontaneous . While the overall scale of $\rho(x)$ involves the defect persistence length and nematic elasticity, the length scale over which charges separate is directly governed by the width of the interface. One can similarly analyze the defect charge distribution setup by a spatial activity pattern that switches from extensile ($v<0$) to contractile ($v>0$). At the interface the activity is forced to vanish, leading to an accumulation of $+1/2$ defects bordered on either side by a compensating layer of $-1/2$ defects. The charge and polarization distribution continue to be given by Eqs. 51 and 50, respectively, independent of any spatial variations in the sign of activity.

VI Discussion

Topological defects play a foundational role in characterizing ordered media, being fingerprints of broken symmetry. In active nematics, they acquire additional dynamical character due to the breakdown of detailed balance. By emphasizing the dominant role of defects as drivers of flow, we have developed a detailed hydrodynamic theory of active defects to capture the various dynamical states of a noisy 2D active nematic on a substrate. This approach allows us to analytically treat both active turbulence, and defect ordering at higher activity. Our results on spatiotemporal defect chaos in the active turbulent regime are consistent with previous numerical work giomi2015geometry ; hemingway2016correlation , and provide a tractable starting point to address the large scale flow signatures of a strongly interacting defect gas. At high enough activity, torques acting on the $+1/2$ disclinations become strong enough to collectively align the moving defects into a spontaneously flowing liquid. Our analysis identifies a definite physical mechanism that drives the defect ordering transition and explains the underlying reason for polar defect order. Finally, extending our treatment to handle spatially inhomogeneous activity, we demonstrate that activity gradients can act as “electric fields” (see also Appendix B) that can be used to corral defects and segregate them. Understanding situations where activity can be spatially or temporally manipulated is the first step to controlling and patterning structure, along with facilitating targeted transport in active matter.

The phase transitions demarcating the different states of defect organization are depicted in Fig. 3 through a phase diagram constructed in terms of the parameters of our theoretical model. To make contact with possible experimental realizations of our predictions, it is useful to cast the phase boundaries in terms of variables that are easier to control in experiments. Both translational noise as captured by our effective temperature $T$ and rotational noise captured by $D_{R}$ are dominantly of nonthermal origin and a priori unknown. Furthermore, they will generally be controlled by active processes, which renders them dependent on the system’s activity. To construct a phase diagram in terms of two independent and in principle theoretically accessible axes, we assume a generic dependence of the effective temperature and the rotational diffusion on activity as $T=T_{0}+T_{1}v^{2}$ and $D_{R}=D_{R}^{0}+D_{R}^{1}v^{2}$, with $T_{0}$ and $D_{R}^{0}$ as passive contributions, and construct a phase diagram in terms of the nematic stiffness $K$ and activity $|v|$. The resulting phase diagram is shown in Fig. 6. In biofilament-motor complex suspensions sanchez2012spontaneous ; kumar2018tunable , the passive elasticity of the nematic, is primarily controlled by filament length and density, whereas activity can be tuned by changing ATP and possibly motor concentration. Of course activity will also affect the nematic elasticity kumar2018tunable , but this can be considered a higher order effect. A different realization of an active nematic involves perfusing a biocompatible molecular liquid crystal with bacteria Zhou2014 . Here the elasticity of the liquid crystal medium is known and well-controlled, while the activity can be tuned by regulating the bacterial concentration. We look forward to experiments that will direct tests of our results.

Our work opens up several important future directions and leaves open challenges to be resolved. Arguably the most obvious one is the lack of nematic defect order in our model. As discussed in Sec. III.2 and sketched in Fig. 4, the absence of nematic defect order is not the result of the coarse-graining scheme or closure approximations used. It instead is a basic feature of the structure of the dynamics of defect as quasi-particles that we use as our starting point. These equations were derived perturbatively in Ref. shankar2018defect and the resulting active torques only favor an effective polar alignment of $+1/2$ disclinations in an unbound defect gas. One may then ask whether a torque that favors nematic alignment of the $+1/2$ defects may be obtained at higher order in activity. We believe the following arguments demonstrate that this is not possible in the context of a dry nematic. To understand this, we first stress that the mechanism for defect ordering described here relies on two ingredients: the balance of elastic forces and propulsive fluxes to ensure ${\bf j}_{\rho}\approx{\bf 0}$ and prevent charge segregation, and the active torques that favor local buildup of polar order in the presence of an elastic distortion. Any active torque with nematic symmetry that were to arise in a non-perturbative treatment will have to involve an alignment interaction between the $+1/2$ defect polarization and the local phase gradient (${\bf v}_{n}$) as these are the only vector fields in the problem. Such a torque cannot, however, support pure nematic defect order on large scales because many-body screening will force ${\bf v}_{n}$ to decay to zero, making it impossible to balance a finite elastic distortion against an active charge current involving the nematic order parameter of defects. This seems to exclude the possibility of purely nematic defect ordering in an active nematic on a substrate. One possible way out is to include hydrodynamic interactions and momentum conserving flow. Genuine fluid flow can provide nonlocal alignment interactions in addition to new terms controlling charge transport. Solvent flow is expected to be important, at least partly, in experiments involving micotubule-kinesin suspensions, where nematic ordering of defects has been previously reported decamp2015orientational , and is also known to stabilize defect lattices in simulations doostmohammadi2016stabilization . Extending our framework to include viscous flow and account for nematic defect ordering is a significant challenge for the future.

Another direction for future research is to use our framework to quantify design principles for engineering active microfluidic devices. Controlling flows and defect organization through confinement keber2014topology ; ellis2018curvature ; guillamat2017taming ; opathalage2019self invites future investigations into the nontrivial role of boundaries and curvature in active matter. Our work on defect segregation through activity gradients poses new avenues for exploration in the context of spatio-temporal control of activity. Furthermore, in the spirit of metamaterial design, controlling flow and patterning structures in active matter is essential to use active devices for applications. In this regard, marrying optimal transport with active fluid hydrodynamics is the next step towards achieving this goal.

Acknowledgements.

We thank Mark Bowick and Sriram Ramaswamy for illuminating discussions. We also acknowledge useful discussions with Juan de Pablo, Rui Zhang and Margaret Gardel. This work was supported by the National Science Foundation through the award DMR-1609208 (M.C.M and S.S.), at KITP under Grant No. PHY-1748958 and at the Aspen Center for Physics under grant PHY-1607611. S.S. would like to thank both the KITP and the Aspen Center for Physics for their support and hospitality during the completion of some of this work.

Appendix A Derivation of active defect hydrodynamics

We briefly recapitulate the order parameter based description of active nematics to both make contact with conventional modeling strategies and also remind ourselves of the various parameters in the model. The continuum nematodynamic equations are just as for the order parameter ${\bf Q}$ of a passive 2D liquid crystal beris1994thermodynamics , given by

[TABLE]

with ${\bf u}$ the flow velocity, $2\Omega_{\mu\nu}=\nabla_{\mu}u_{\nu}-\nabla_{\nu}u_{\mu}$ the vorticity tensor and $\gamma$ the rotational viscosity. The flow alignment term involves the coupling $\lambda$ and the deviatoric strain rate tensor $(\nabla_{\mu}u_{\nu})_{ST}=\nabla_{\mu}u_{\nu}+\nabla_{\nu}u_{\mu}-\delta_{\mu\nu}(\bm{\nabla}\cdot{\bf u})$. The molecular field ${\bf H}=b(1-S^{2}){\bf Q}+K\nabla^{2}{\bf Q}$ with $b>0$ controls the mean-field isotropic to nematic transition at equilibrium and $K$ is the Frank elastic constant.

In the presence of a frictional substrate and activity, we supplement this with a force balance equation for the flow velocity, involving active stresses, $-\Gamma{\bf u}+\bm{\nabla}\cdot\bm{\sigma}^{a}={\bf 0}$. We have assumed that friction with the substrate $\Gamma$ screens the flow on large scales and have retained only the dominant active stress, $\bm{\sigma}^{a}=\alpha{\bf Q}$, neglecting all passive elastic contributions. The activity $\alpha$ corresponds to the average force dipole exerted by the microscopic apolar active units on the surrounding fluid ramaswamy2010mechanics ; marchetti2013hydrodynamics , with $\alpha>0$ for contractile systems and $\alpha<0$ for extensile systems. For an isolated $+1/2$ disclination positioned at ${\bf r}_{i}^{+}$, the active flow generated at the core of the defect is ${\bf u}({\bf r}_{i}^{+})=v{\bf e}_{i}$, where $v=\alpha/(a\Gamma)$ relates the defect motility to the activity and ${\bf e}_{i}=a\bm{\nabla}\cdot{\bf Q}({\bf r}_{i}^{+})$ is the defect orientation as mentioned in the main text. Below we shall derive the hydrodynamic description of active defects.

A.1 Phase gradient dynamics

The topological constraint (Eq. 5) relates the phase gradient ${\bf v}_{n}$ to the defect charge density through $\hat{{\bf z}}\cdot\left(\bm{\nabla}\times{\bf v}_{n}\right)=2\pi\rho$. As defects are created or annihilated only in pairs, the charge density $\rho$ is locally conserved and we have a continuity equation for the dynamics of $\rho$, as noted in Eq. 11

[TABLE]

Note that this expression is valid at the fluctuating level in the Stratanovich convention. As the phase gradient ${\bf v}_{n}$ is single valued, we can commute derivatives on it. Writing $2\pi\rho=\hat{{\bf z}}\cdot\left(\bm{\nabla}\times{\bf v}_{n}\right)$ and commuting the time and space derivatives, Eq. 53 gives

[TABLE]

where $\bm{\epsilon}$ is the 2D Levi-Civita tensor and $\Pi$ is a smooth scalar function free of singularities. Away from defects ${\bf j}_{\rho}={\bf 0}$, and a smooth director fluctuation then obeys

[TABLE]

on large scales, where $f_{\theta}$ is unit white Gaussian noise and $T$ is the noise strength that functions as an effective temperature. All nonlinearities affecting the smooth phase fluctuations in the ordered 2D active nematic are known to be perturbatively irrelevant on large enough length scales mishra2010dynamic ; shankar2018low . Comparing Eq. 55 and Eq. 54, we set $\Pi=(K/\gamma)\bm{\nabla}\cdot{\bf v}_{n}+\sqrt{2T/\gamma}f_{\theta}$. Putting it all together, we then have

[TABLE]

Upon setting $\kappa=K/\gamma$ and neglecting the noise in a mean-field description, we obtain Eq. 15 in the main text.

A.2 Current constitutive equation

The motion of a $\pm 1/2$ disclination in an active nematic was derived in Ref. shankar2018defect to be

[TABLE]

where $\mu\propto\gamma^{-1}$ is the defect mobility and $\bm{\xi}_{i}(t)$ is unit white noise. We emphasize that $v\propto\alpha$ can be of either sign, with extensile and contractile systems propelling the $+1/2$ defect in opposite directions giomi2013defect , while the $-1/2$ defect remains diffusive and non-motile. The particular form of the Magnus-like force in Eq. 57 was chosen to recover the passive Coulomb interaction between bound defects. In the absence of free unbound defects, the phase gradient at any point due to bound defect pairs with charge $q_{i}=\pm 1/2$ is

[TABLE]

where ${\bf r}_{i}$ is the position of the $i$th defect and $a$ is the defect core size that provides a microscopic cutoff. Using Eq. 58 in Eq. 57, one can check that the correct form of the passive elastic force is obtained when considering just bound defect pairs. In the presence of unbound defects, Eq. 58 is no longer applicable and we instead have to use Eq. 56 to obtain the phase gradient.

Using the defect equations of motion (Eq. 57) in the definition of the fluctuating current (${\bf j}_{\pm}=\sum_{i}\dot{{\bf r}}^{\pm}_{i}\delta({\bf r}-{\bf r}_{i}^{\pm})$), we obtain

[TABLE]

where $\bm{\xi}_{\pm}({\bf r},t)$ are unit space-time white noises. For spatially varying $n_{\pm}$, the multiplicative noise must be interpreted in Ito style. Writing $D_{0}=\mu T$, replacing $\pi\mu K$ by $\nu\kappa$, and neglecting noise at the mean-field level, we obtain Eqs. 8 and 9. The corresponding fluctuating expressions for ${\bf j}_{\rho}$ and ${\bf j}_{n}$ can then be trivially obtained.

A.3 Polarization dynamics

The orientational dynamics of ${\bf e}_{i}$ including active torques was derived in Ref. shankar2018defect to be

[TABLE]

This is written here in a form that is appropriate for treating unbound defects. Unpackaging the various terms, the top line of Eq. 61 involves both passive and active terms that relax or enhance $|{\bf e}_{i}|$, which is a measure of local nematic distortion. The second line of Eq. 61 includes orientational torques that leave $|{\bf e}_{i}|$ fixed. The first of these is the important active torque that causes the $+1/2$ defect orientation to self-align with the local phase gradient. The second is a passive elastic torque that causes the $+1/2$ defect to reorient in a splayed phase gradient. This can be easily obtained by noting that, away from defects, an external phase gradient causes the nematic director to precess at a rate given by $\bm{\nabla}\cdot{\bf v}_{n}$ (from Eq. 55). Following the analysis of Ref. shankar2018defect , this generates the required elastic torque in Eq. 61. As can be seen, this elastic torque is subdominant to the active one, as it only arises for $\bm{\nabla}\cdot{\bf v}_{n}\neq 0$. Similar elastic torques have been obtained by different means previously as well keber2014topology ; vromans2016orientational ; tang2017orientation . The final two terms in the last line of Eq. 61 are noise, with $D_{R}$ as the rotational diffusion constant, $\eta_{i}(t)$ as unit white Gaussian noise, and the longitudinal noise $\bm{\nu}_{i}(t)$ is Gaussian with zero mean and correlations

[TABLE]

The form of the noise correlator is obtained by using essentially a fluctuation-dissipation like relation in the limit $v=0$. For simplicity we use the same effective noise strength $\sim T$ here as other choices don’t change the results in any qualitative way (also see Ref. shankar2018defect ). Coarse-graining Eq. 61, we obtain a fluctuating hydrodynamic equation for the defect polarization density ${\bf p}$,

[TABLE]

This equation is not closed as it involves the second order tensorial moment ${\bf M}({\bf r},t)=\sum_{i}{\bf e}_{i}(t){\bf e}_{i}(t)\delta({\bf r}-{\bf r}_{i}^{+}(t))$. The Gaussian noise $\bm{\Lambda}$ is spatio-temporally unit white. We have also neglected subdominant terms in the noise involving $n_{+}|{\bf v}_{n}|^{2}$ and the traceless part of ${\bf M}$, along with $\mathcal{O}(\nabla^{2})$ contributions. These do not affect the results presented.

To close the moment heirarchy, we use a Ginzburg-Landau ansatz taking both ${\bf p}$ and ${\bf v}_{n}$ to be equally small and slowly varying. In addition, we set ${\bf M}=(m/2){\bf 1}$ and disregard all higher anisotropic moments of the orientation that decay on short time scales $\sim D_{R}^{-1}$, neglecting any nematic ordering of the $+1/2$ defects for now. As the active torques only generate alignment between ${\bf p}$ and ${\bf v}_{n}$, within our model it is impossible to obtain nematic order of defects directy from the isotropic state. The last thing left to do is then determine $\mathrm{tr}({\bf M})=m$ in terms of the variables retained. At lowest order, without noise, we find

[TABLE]

where $\zeta$ is a nondimensional activity. The positive dimensionless function $f(\zeta)=1+4\zeta+\mathcal{O}(\zeta^{2})$ involves the leading $\zeta$ correction arising from eliminating the third order moment [$\langle\sum_{i}{\bf e}_{i}|{\bf e}_{i}|^{2}\delta({\bf r}-{\bf r}_{i}^{+})\rangle\simeq-v\pi\mu\gamma(8T^{2}/D_{R}K^{2})n_{+}\,\bm{\epsilon}\cdot{\bf v}_{n}$]. This completes the derivation of defect hydrodynamics.

Appendix B Defect alignment and electric field analogy

In this Appendix we provide further details on the alignment torques acting on the $+1/2$ disclinations, with an intuitive interpretation in terms of an electrostatic analogy. The full equation for the angular dynamics of the $+1/2$ defect orientation is given in Appendix A.3. To develop the electrostatic analogy, we only consider the angular dynamics of the $+1/2$ defect and neglect noise and passive elastic torques. Writing ${\bf e}_{i}=|{\bf e}_{i}|(\cos\psi_{i},\sin\psi_{i})$ for an individual defect polarization, in the presence of a finite phase distortion ${\bf v}_{n}$ the active torques generate an angular velocity given by

[TABLE]

For the last equality, we have recast the active alignment torque as deriving from an effective alignment “energy” $\mathcal{H}=v\sum_{i}\hat{{\bf z}}\cdot({\bf e}_{i}\times{\bf v}_{n})$. This is simply a formal rewriting of the dynamical torque with no real energetic underpinning. Now, from the topological charge constraint $\hat{{\bf z}}\cdot(\bm{\nabla}\times{\bf v}_{n})=2\pi\rho$ (Eq. 5), we identify an “electric field” ${\bf E}=\bm{\epsilon}\cdot{\bf v}_{n}$. This allows us to rewrite the charge constraint in the form of Gauss’ law,

[TABLE]

This identification is the same as in the Maxwell analogy for superfluid dynamics ambegaokar1980dynamics . As discussed in Sec. V, an activity gradient can enforce a locally nonvanishing ${\bf v}_{n}$, which in turn acts as an “electric field” through the above identification.

We can then rewrite the effective alignment interaction $\mathcal{H}$ in terms of the electric field ${\bf E}$ as

[TABLE]

It is then evident that the active alignment of a neutral defect pair takes exactly the form of the alignment of an electric dipole with dipole moment $-v{\bf e}_{i}$ in an electric field ${\bf E}$. This allows us to easily and intuitively interpret the consequences of the active torque. As shown in the schematic in Fig. 4 and in more detail in Fig. 7, a neutral pair of defects generates an elastic distortion or phase gradient ${\bf v}_{n}$ plotted as blue lines in Fig. 7a. Conversely, one can reinterpret the same picture in terms of ${\bf E}$, as shown in Fig. 7b. Here the “electric field” lines are shown in purple and are orthogonal to the local phase gradient. As stated in Eq. 68, the active torque can effectively be seen as akin to a dipole alignment interaction. In the extensile case ($v<0$), the effective dipole moment is $|v|{\bf e}_{i}$ and preferentially aligns with the local electric field. A similar argument works for the contractile case as well with the signs flipped. In both cases, the effective alignment between the $+1/2$ defects is polar. In short, Fig. 7 provides a different interpretation of the alignment induced by active torques (as sketched in Fig. 4) in terms of an electrostatic analogy.