The role of fluid flow in the dynamics of active nematic defects

Nematic order has been widely observed in two-dimensional (2D) realizations of active matter [1, 2], from vertically vibrated rods [3] to mixtures of cytoskeletal filaments and associated motor proteins [4, 5, 6], bacterial suspensions [7, 8, 9], cell sheets [10, 11, 12, 13], and even developing organisms [14]. Nematic order has only quasi-long-range (power-law) falloff in 2D active systems [15, 16, 17], as it does in equilibrium, and is easily destroyed by active stresses that drive spontaneous flows accompanied by the proliferation of topological defects [4]. A distinctive property of active nematics is that the comet-like $+1/2$ disclination becomes motile [4, 18, 19], allowing for an activity-driven defect unbinding transition [20] to a turbulent-like state with chaotic spatio-temporal dynamics. Both experiments and numerical studies based on the solution of continuum models have demonstrated that this complex dynamics can be characterized by focusing either on the statistics of flow vortices or on that of topological defects [2, 21, 22, 23, 24, 6].

Recent work has focused on formulating an effective description of defects in active nematics as quasiparticles [18, 25, 19, 21, 26, 20, 27, 28]. This approach parallels well-established work in equilibrium systems, where defects can be described as a gas of Coulomb charges and are known to drive order-disorder transitions in 2D [29].

An important new ingredient in 2D active nematics is that the active self-propulsion of the $+1/2$ defects requires taking into account the local geometric polarity of the nematic texture near the defect core [18, 19, 30, 31]. The self-propulsion arises from local active flows generated by the distortion of the texture due to the defect itself, and is absent for the $-1/2$ defects due to their threefold symmetry. Shankar et al. [20] have shown that the dynamics of the gas of unbound defects can be mapped onto that of a mixture of self-propelled (the $+1/2$ defects) and passive (the $-1/2$ defects) quasiparticles with Coulomb interactions and aligning torques, provided that the texture gradients generated at the core of one defect by the other defects can be treated as quasistationary. More recently, the inclusion of multi-defect interactions, within the same quasistatic approximation, has revealed new non-central and non-reciprocal forces in multi-defect textures [28]. Important corrections to the effective dynamics also arise from the fact that the phase field induced at the core of a defect by all the others depends on all the defect velocities, as well as on their past history of accelerations, but including such effects remains a formidable challenge. These previous works all treat the case of an active nematic on a substrate, where friction dominates over viscous dissipation and the compressible flow ¹¹1Here and below we use the term ‘compressible flow’ to refer to flows that do not satisfy incompressibility (hence $\bm{\nabla}\cdot\bm{u}\not=0$), although we assume the density to be constant. This could be achieved, for instance, in systems where density is not conserved due to birth and death. is slaved to the gradient of the active stress. Recent numerical works [32] show that friction with the substrate promotes the kink wall formation in the trail of $+1/2$ defects and this may have a persistent memory effect on the collective behavior of many defects [33, 34]. On the other hand, this limit case of a compressible flow offers an analytically tractable regime, where the flow can be readily eliminated in favor of the nematic texture [35, 36] and generates active corrections to the nematic stiffness and new nonlinear terms in the equation for the nematic order parameter [35, 36]. In most experimental realizations, however, the flow is incompressible, which results in long-range hydrodynamic effects neglected in previous work.

In this paper we formulate the dynamics of defects as quasiparticles in the presence of arbitrary flows and apply them to flows governed by the Stokes equation, which balances dissipative forces (friction and viscosity) with pressure and active stress gradients. We find that, as shown in Ref. [26], incompressibility reduces the active propulsive speed of an isolated $+1/2$ defect. We also show, for the first time, that non-local hydrodynamic effects become important for multi-defect configurations and yield a finite mobility for the $-1/2$ defects. We demonstrate this explicitly by calculating the pressure-gradient forces on a defect dipole and showing that incompressibility gives rise to hydrodynamic interactions that renormalize the propulsion of positive defects and mediate non-reciprocal active defect-defect interactions additional to the familiar Coulomb force. These couplings due to pressure gradients are distinct from the perturbative multi-defect interactions obtained in Ref. [28] for a compressible nematic.

By coarse graining the defect equations, we also obtain continuum equations for the defect densities and the polarization density of the $+1/2$ defects coupled to an arbitrary fluid flow. We find that the $+1/2$ defects indeed behave like polar particles in a fluid, though advected and rotated by the flow velocity as well as by the Magnus-type force generated by all the other defects [27]. Our work extends the continuum equations obtained in Ref. [27] for a compressible nematic where flows are slaved to texture deformation. The defect hydrodynamics obtained here applies to both compressible and incompressible flows and provides a useful starting point for describing experimentally relevant situations with spatially inhomogeneous activity, where defects can be trapped and guided by activity gradients [27, 37], as observed for instance in active actomyosin suspensions [38].

Our work relies on describing topological defects and their equation of motion as the zeros of an appropriate nematic order parameter. This method was originally introduced by Halperin [39] and extended by Mazenko to $O(n)$-symmetric order parameters in the context of phase-ordering kinetics [40]. It was shown in Refs. [41, 42, 43] that this is an efficient numerical tool for tracking the positions and velocities of vortices in superfluids and dislocations in crystals. The Halperin-Mazenko (HM) method allows us to derive both the discrete defect dynamics and the dynamics of the corresponding defect density fields from the coupled continuum equations for the flow and the nematic order parameter. For overdamped compressible flows we recover the equations previously derived in [20].

The rest of the paper is structured as follows. In Sec. 2 we briefly summarize the well-established continuum model of active nematic hydrodynamics. In Sec. 3, we present the HM method based on the conservation of topological defects as zeros of the order parameter and relate the defect velocity to the nematic texture. Here we also derive the equation of motion for the polarization of a $+1/2$ defect in the quasistatic-phase approximation. In Sec. 4, we evaluate the defects propulsive velocity arising from activity-induced flows. We show that our method recovers previous results for isolated defects for both compressible [18, 19, 20] and incompressible [26] flows. We also present new results for a defect dipole in an incompressible nematic, where flows driven by pressure gradients yield new contributions to the motility of both the positive and the negative defect, as well as non-reciprocal active forces. In Sec. 5 we derive general hydrodynamic equations for a defect gas, applicable to both compressible and incompressible flows. Finally, we summarize and conclude in Sec. 6. The Appendices detail analytical computations of the flow field induced by a variety of specific defect configurations.

In continuum models, the hydrodynamics of uniaxial active nematics is described by the $Q$-tensor $\textendash$ the order parameter for the orientational field $\theta(\mathbf{r},t)$ and its associated nematic director (line field) $\hat{n}=[\cos\theta,\sin\theta]$. The head-tail symmetry of the local director ($\hat{n}\equiv-\hat{n}$) and the property that $Q$ vanishes in the isotropic phase constrain the $Q$-tensor to be symmetric and traceless. In 2D, $Q_{ij}=S(2\hat{n}_{i}\hat{n}_{j}-\delta_{ij})$ with $S$ the magnitude of the nematic order. $Q$ has two independent components $Q_{xx}=-Q_{yy}$ and $Q_{xy}=Q_{yx}$. The evolution of the $Q$-tensor is given by the the Edwards-Beris equation, which in 2D reads [44, 1]

The last two terms on the right hand side (RHS) represent the relaxation of the order parameter to minimize isotropic distortions, with elastic constant $K>0$, and deviations from the uniform ordered state, $S_{0}^{2}=1$. Here $\gamma$ is the rotational friction coefficient of the nematic. The remaining terms are given by the material time derivative in the presence of the fluid flow velocity $\mathbf{u}$ with vorticity $2\Omega_{ij}=\partial_{i}u_{j}-\partial_{j}u_{i}$ and the tendency of nematics to align with the flow through linear coupling to the symmetric traceless part of the flow strain rate, $2E_{ij}=\partial_{i}u_{j}+\partial_{j}u_{i}-\delta_{ij}\partial_{k}u_{k}$, and to the flow divergence, $\bm{\nabla}\cdot\bm{u}$. We have neglected higher-order flow alignment terms which are nonlinear in the order parameter. The flow alignment parameters $\lambda$ and $\lambda^{\prime}$ are dimensionless but nonuniversal. They are controlled by the degree of nematic order and the microscopic shape of the underlying mesogens. The alignment with shear flow, controlled by $\lambda$, directly affects the defect velocity, whereas the coupling to the flow divergence, controlled by $\lambda^{\prime}$, is linear in the $Q$-tensor, which vanishes at the defect core. As a result this term may be neglected in determining quasistatic defect dynamics.

Flows in active nematics are generally characterized by very low Reynolds numbers and follow the Stokes equation,

describing force balance between friction with the substrate ($\mu$), viscous dissipation controlled by the shear viscosity $\eta_{0}$, potential forces due to pressure, $\Pi$, and active stress, $\bm{\sigma}^{a}$. Pressure becomes important for incompressible flows, when it is determined by the incompressibility constraint $\bm{\nabla}\cdot\bm{u}=0$. The active stress is proportional to the order parameter $\bm{\sigma}^{a}=\alpha_{0}\bm{Q}$, where $\alpha_{0}$ is the activity coefficient, such that $\alpha_{0}<0$ corresponds to extensile activity, and $\alpha_{0}>0$ to contractile activity. Finally, in Eq. (2) we neglect elastic liquid crystalline stresses that are higher order in the gradients of the order parameter field compared to the active stress.

Dimensional analysis reduces the parameter space to two independent dimensionless numbers. We rescale time $t=\tau\tilde{t}$ by the nematic relaxation time $\tau=\gamma/g$, and space $\mathbf{r}=\xi\tilde{\mathbf{r}}$ by the nematic coherence length $\xi=\sqrt{K/g}$, the intrinsic length scale that determines the defect core size and the minimal length scale of spatial variations of the order parameter. The fluid velocity then scales as $\mathbf{u}=(\xi/\tau)\tilde{\mathbf{u}}$ and the fluid pressure as $\tilde{p}=p/(\mu\xi^{2}/\tau)$. The relevant parameters for the dimensionless equations then become the rescaled activity $\alpha=\alpha_{0}/(\mu\xi^{2}/\tau)$, measuring the strength of the active force relative to friction, and the rescaled shear viscosity $\eta=\eta_{0}/(\mu\xi^{2})$, measuring viscous forces relative to friction. From now on we will drop the tildes on the dimensionless variables and treat the friction-dominated regime $\eta\rightarrow 0$. The general case with an interplay between friction and viscous dissipation will be treated separately elsewhere.

The nematic order field has an elegant complex representation $\psi=Se^{i2\theta}=Q_{xx}+iQ_{xy}$ [45, 26, 27]. The dynamics of the complex $\psi$-field following from Eq. (1) is given by a generalized complex Ginzburg-Landau (CGL) equation

with

and the incompressibility condition $\partial_{\bar{z}}u+\partial_{z}\bar{u}=0$, with $z=x+iy$ and $\partial_{z}=(\partial_{x}-i\partial_{y})/2$.

Topological defects appear in the $\psi$-field as localized regions (of size $\xi$) where $S$ vanishes and around which the nematic orientation $\theta$ winds by $\pi$. The method described in the next section allows us to derive the kinematics of the $\pm 1/2$ nematic defects from the dynamics of the zeros of the $\psi$-field.

The complex nematic order parameter $\psi$ is always a single-valued function even in the presence of defects. This restricts any phase discontinuity in $\theta(\mathbf{r})$ to have a half-integer fundamental winding number. Taking an arbitrary contour $C$ surrounding a set of elementary defects, the net phase jump is

where the sum runs over all defects $\beta$ contained inside $C$. We consider only energetically stable defects of charge $q_{\beta}=\pm 1/2$. Higher charge defects are topologically stable in 2D but are energetically unstable to disassociation into elementary defects, since their intrinsic energy grows as the square of the charge [46].

Defects are accurately located by the zeros of the complex order parameter. Their dynamics is typically derived by solving the hydrodynamic field equations in the comoving frame of the defect together with asymptotic matching of inner- and outer-core solutions under appropriate initial and boundary conditions [47, 48, 19]. The defect velocity and the phase gradients are then related through the mobility coefficient which acquires a logarithmic dependence on velocity. The field solutions are derived perturbatively in phase gradients and defect velocity, and are challenging to generalize when a generic fluid flow is coupled with the evolution of the order parameter. The HM method [40, 49], in which the defect dynamics is identified with the motion of the zeros of the complex order parameter, allows us to directly include the coupling to fluid flow.

The starting point is to represent a configuration of defects by the Dirac delta function $\delta(\psi)$, which picks up the defect positions $\mathbf{r}_{\beta}$, with $\beta=1,2,\cdots$ the particle number index. The field variables map into discrete defect positions $\{\mathbf{r}_{\beta}\}$,

where the Jacobi determinant $D(\bm{r},t)=\det(\partial\psi)$ is

with $\epsilon_{ij}$ the Levi-Civita anti-symmetric tensor. In complex coordinates the determinant is

The field $D$ captures the zeros of the order parameter and vanishes everywhere except at the defect positions. It obeys a conservation law determined by the hydrodynamic equation of the $\psi$-field and readily obtained by taking time derivatives on both sides of Eq. (7) [40], with the result

where the conservative current $J_{i}^{(\dot{\psi})}$ is given by

or in complex form

The defect charge is related to the sign of the determinant evaluated at the defect position, $q_{\beta}=\frac{1}{2}\,\text{sgn}(D)|_{\bm{r}=\bm{r}_{\beta}}$[49]. This relation allows us to connect the determinant field $D$ with the defect charge density field for a given configuration of defects, namely

The defect charge density $\rho$ also satisfies a conservation law [40, 50]:

where the current density is

Finally, for a given defect configuration, the density current can be also expressed in terms of the velocity of the defects $\bm{v}_{\beta}=\dot{\mathbf{r}}_{\beta}$:

The defect velocity is therefore determined by the field $D$ and its current density:

which explicitly relates the defect velocity to derivatives of $\psi$ evaluated at defect positions. In complex form $v_{\beta}=v_{x,\beta}+iv_{y,\beta}$ is given by

This is the starting point for obtaining an explicit form for the defect velocity in terms of the flow field and nematic distortions, using the profile of $\psi$ in the vicinity of the defect.

In this section we obtain explicit expressions for the flow induced by specific defect configurations, which, in turn, drives defect dynamics via Eq. (3.2). By eliminating the flow in favor of defect degrees of freedom we can identify the advection at the defect core as a defect propulsion velocity. We consider both compressible and incompressible ($\nabla\cdot\mathbf{u}=0$) flows. In the latter case, we calculate for the first time the flow induced by a defect dipole, which was previously studied numerically in Ref. [18].

We have derived ordinary differential equations describing the dynamics of defects as quasiparticles. In the turbulent-like nematic state where defects proliferate, it is useful to formulate a continuum model of the defect gas that describes the dynamics at length scales large compared to the mean defect separation. This was carried out in Ref.  [27] for an overdamped compressible active nematic. Here we generalize to the experimentally important case where the fluid flow is incompressible.

We describe a defect configuration in terms of microscopic density fields: defect charge and number densities, $\hat{\rho}$ and $\hat{n}$, their associated current densities, $\mathbf{\hat{J}}^{(\rho)}$ and $\mathbf{\hat{J}}^{(n)}$, and defect polarization density $\mathbf{\hat{p}}$:

where we use a hat to distinguish microscopic fields from coarse-grained ones obtained by averaging over ensembles of defect configurations. The microscopic charge density satisfies the requisite conservation law.

The coarse-graining method consists of taking ensemble averages over different defect configurations described by a distribution function $P(\{\mathbf{r}_{\beta}\})$,

and factorizing high moments to close the equations.

To obtain the defect hydrodynamic equations, we proceed in two steps. First, we consider a tagged defect moving in the effective mean-field induced by the rest of the defects that is described by a coarse-grained flow field $\bm{u}(\mathbf{r})=\langle\hat{\mathbf{u}}(\mathbf{r},t;\{\mathbf{r}_{\beta}\})\rangle$ and a coarse-grained nematic texture $\bm{v}_{n}(\mathbf{r})=\langle\bm{\nabla}\hat{\tilde{\theta}}(\mathbf{r},t;\{\mathbf{r}_{\beta}\})\rangle$. The dynamics of a defect in this mean field is then governed by the equations (for simplicity we consider the case with no flow alignment $\lambda=0$)

where for any vector $\bm{V}=(V_{x},V_{y})$ we denote the clockwise rotation corresponding to the multiplication with the Levi-Civita tensor $\epsilon_{ij}$ as $\bm{V}^{\perp}=\bm{\epsilon}\cdot\bm{V}=(V_{y},-V_{x})$. These equations couple defect dynamics to the coarse-grained flow and nematic texture, with $\bm{f}_{\beta\gamma}=\mathbf{r}_{\beta\gamma}/|\mathbf{r}_{\beta\gamma}|^{2}$, which $\mathbf{r}_{\beta\gamma}=\mathbf{r}_{\beta}-\mathbf{r}_{\gamma}$, which is the interaction force which relates the nematic distortion with defect density as

is the interaction force from the defect at position $\mathbf{r}_{\gamma}$ on the defect at position $\mathbf{r}_{\beta}$.

The phase gradient $\bm{v}_{n}$ contains the singular phase induced by all defects and relates to the macroscopic charge density $\rho$ through Stokes’ theorem applied in Eq. (5),

This determines $\bm{v}_{n}$ up to a potential vector field which, in the quasistatic approximation, satisfies

The macroscopic fluid flow $\mathbf{u}$ follows the Stokes Eq. (2), but with the active stress expressed in terms of the density and polarization density of positive defects,

Incompressibility additionally requires $\bm{\nabla}\cdot\mathbf{u}=0$.

The continuum equations for the defect densities and polarization density can then be obtained by taking the time derivative of the microscopic fields, Eqs. (42), using the defect dynamics given in Eqs. (48), and applying the coarse-graining procedure. Neglecting defect creation and annihilation, the equations for the defect density fields $\rho$ and $n$ take the form

In the presence of defect annihilation/creation, the RHS of Eq. (56) needs to be supplemented by the kinetic rates describing these processes. The macroscopic current densities are given by

The evolution of the microscopic polarization density follows from taking the time derivative of Eq. (46), using Eqs. (48) and expressing the sums over defects in terms of the microscopic densities and their current densities, with the result

where $\bm{U}=\mathbf{u}+4\bm{v}_{n}^{\perp}$ is the net flow including the one generated by nematic distortion. Upon coarse graining, we assume $\langle\hat{\bm{p}}(\mathbf{r})\hat{\rho}(\mathbf{r}^{\prime})\rangle\approx{\bm{p}}(\mathbf{r})\rho(\mathbf{r}^{\prime})$ and $\langle\hat{\bm{p}}(\mathbf{r})\hat{n}(\mathbf{r}^{\prime})\rangle\approx{\bm{p}}(\mathbf{r})n(\mathbf{r}^{\prime})$, and obtain the continuum equation

The polarization equation does not contain terms controlling growth or decay of the magnitude of $\bm{p}$ nor any stiffness/diffusive terms because we assumed $|\mathbf{e}_{\beta}|=2$ (not normalized) and did not include noise in the defect dynamics. The first term on the RHS of Eq. (5) controls changes in the direction of polarization and arises from the active torque in the evolution of $\mathbf{e}_{\beta}$. The next term only appears when the flow is compressible and it describes the fact that density variations associated with compressible flow can provide sources and sinks of polarization. In addition, Eq. (52) implies that there are also local changes in polar order depending on the charge density field, $\rho$. Finally, the last term on the RHS of Eq. (5) is a new nonlocal correction to the torque that arises from the long-range nature of nematic distortions.

It is useful to compare Eq. (5) for the polarization density to the one obtained in Ref. [27] for a compressible fluid. As mentioned above, here we do not have terms describing polarization growth/relaxation. The kinematic terms describing coupling of polarization to flow do reduce to those obtained in Ref. [27] for compressible fluids. This can be seen by noting that in Ref. [27] the polarization equation was obtained by eliminating $\mathbf{u}$ from the defect dynamics before coarse-graining, using $\mathbf{u}\rightarrow\alpha{\bm{e}}_{\beta}$ in Eq. (LABEL:eq:phat). If we use this in Eq. (LABEL:eq:phat), and then coarse grain, using the same moment closure as in Ref. [27], we recover the same results.

The continuum equations obtained here apply generally to describe the coupling of defect dynamics to any flow field. They resemble the equations for self-propelled polar particles in a fluid, with the difference that the $+1/2$ defects are advected and rotated not only by the fluid flow velocity $\mathbf{u}$ but also by the Magnus-type force proportional to $\bm{v}_{n}^{\perp}$. The fact that only positive defects act as a source of flow in the Stokes equations for $\mathbf{u}$ is consistent with the fact that only the $+1/2$ defects behave like active quasiparticles. The negative defects are advected and rotated by flows, like passive particles, but also contribute to the phase deformations described by $\bm{v}_{n}$.

To summarize, we have analyzed the effect of arbitrary fluid flow on the motion of $\pm 1/2$ nematic defects and derived both equations describing the dynamics of defects as quasipartciles and hydrodynamic equations describing the dynamics of the defect gas on length scales large compared to the mean defect separation. By adapting the Halperin-Method formalism to a nematic order parameter, we demonstrate that it provides an effective alternative to perturbative methods of matched asymptotic expansions for deriving defect kinematics, especially for incompressible flows. We reproduce previous results for friction-dominated, compressible flows and find new effects for incompressible systems, where pressure gradients incorporated through the incompressibility constraint yield long-range forces among any two defect pair. We specifically examine the case of a defect dipole embedded in a uniform far-field texture, where analytical calculations are possible. We show that in this case the flows arising from pressure gradients reduce the self-propulsion of the $+1/2$ defects and introduce non-reciprocal forces acting on both defects. A future challenge is to extend our work to describe collective behavior and pattern formation of many interacting defects.