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

We adapt the Halperin–Mazenko formalism to analyze two-dimensional active nematics coupled to a generic fluid flow. The governing hydrodynamic equations lead to evolution laws for nematic topological defects and their corresponding density fields. We find that ±1/2 defects are propelled by the local fluid flow and by the nematic orientation coupled with the flow shear rate. In the overdamped and compressible limit, we recover the previously obtained active self-propulsion of the +1/2 defects. Non-local hydrodynamic effects are primarily significant for incompressible flows, for which it is not possible to eliminate the fluid velocity in favor of the local defect polarization alone. For the case of two defects with opposite charge, the non-local hydrodynamic interaction is mediated by non-reciprocal pressure-gradient forces. Finally, we derive continuum equations for a defect gas coupled to an arbitrary (compressible or incompressible) fluid flow.

Recent work has focused on formulating an effective description of defects in active nematics as quasiparticles [18][19][20][21][25][26][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 is slaved to the gradient of the active stress. Here and below we use the term 'compressible flow' to refer to flows that do not satisfy incompressibility (hence ∇ · u = 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 processes. 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 reference [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 reference [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 reference [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 references [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 section 2 we briefly summarize the well-established continuum model of active nematic hydrodynamics. In section 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 section 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 section 5 we derive general hydrodynamic equations for a defect gas, applicable to both compressible and incompressible flows. Finally, we summarize and conclude in section 6. The appendices detail analytical computations of the flow field induced by a variety of specific defect configurations.

Continuum model of 2D active nematics
In continuum models, the hydrodynamics of uniaxial active nematics is described by the Q-tensor-the order parameter for the orientational field θ(r, t) and its associated nematic director (line field) n = [cos θ, sin θ]. The head-tail symmetry of the local director (n ≡ −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(2n inj − δ 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 Edwards-Beris equation, which in 2D reads [1,44] The last two terms on the right-hand side (rhs) determine the conjugate field H ij = −δF /δQ ij that governs the relaxational dynamics to minimize the free energy with the isotropic elastic constant K > 0 and the uniform ordered state, S 2 0 = 1. Here γ 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 u with vorticity 2Ω ij = ∂ i u j − ∂ 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 = ∂ i u j + ∂ j u i − δ ij ∂ k u k , and to the flow divergence, ∇ · u. We have neglected higher-order flow alignment terms which are nonlinear in the order parameter. The flow alignment parameters λ and λ 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 λ, directly affects the defect velocity, whereas the coupling to the flow divergence, controlled by λ , 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 (μ), viscous dissipation controlled by the shear viscosity η 0 , potential forces due to pressure, Π, and active stress, σ a . Pressure becomes important for incompressible flows, when it is determined by the incompressibility constraint ∇ · u = 0. The active stress is proportional to the order parameter σ a = α 0 Q, where α 0 is the activity coefficient, such that α 0 < 0 corresponds to extensile activity, and α 0 > 0 to contractile activity. Finally, in equation (3) we neglect elastic liquid crystalline stresses that are higher order in the gradients of the order parameter field compared to the active stress. In general, there are other stresses that contribute to the Stokes flow, such as the elastic stress σ e = λH ij + Q ik H kj − Q jk H ki and the Ericksen stress σ E = −K∂ i Q kl ∂ j Q kl . We neglect them here because they are of higher order in the gradients of the alignment tensor Q compared to the active stress. The full hydrodynamic equations contain three length scales: (i) the nematic correlation length ξ = K/g that controls the minimal length scale of spatial variations of the order parameter (and therefore determines the defect core size), (ii) the viscous screening length v = η/γ, and (iii) the active length a = K/|α 0 | determined by the balance of active and elastic stresses. Another important length scale in simulations and experiments is the system size, L. Extensive simulations of hydrodynamic equations, as well as experiments have established that a also controls the mean separation between defects (see, for instance, reference [2]). In the absence of substrate drag, the uniform active nematics are unstable to flow at small activity where a > L and transition to a state of spatio-temporal chaotic flow with a steady number of topological defects when a L. The presence of substrate friction stabilizes the system at small activity, setting a finite threshold value of activity for the set-up of the spontaneous flow instability. Here we are working deep in the uniform nematic state, where the defect core size ξ a is the smallest length scale. It is only in this regime that one can treat defects as quasiparticles. It corresponds to the low defect density approximation where defects are well-separated from each other, which implies that a L. We also set v = 0, since we neglect viscosity compared to drag. Dimensional analysis reduces the parameter space to two independent dimensionless numbers. We rescale time t = τt by the nematic relaxation time τ = γ/g, and space r = ξr by the nematic coherence length ξ = K/g, the only relevant lenghscale in this formulation. The fluid velocity then scales as u = (ξ/τ )ũ and the fluid pressure asp = p/(μξ 2 /τ ). The relevant parameters for the dimensionless equations then become the rescaled activity α = α 0 /(μξ 2 /τ ), measuring the strength of the active force relative to friction, and the rescaled shear viscosity η = η 0 /(μξ 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 η → 0.
The nematic order field has an elegant complex representation ψ = S e i2θ = Q xx + iQ xy [26,27,45]. The dynamics of the complex ψ-field following from equation (1) is given by a generalized complex Ginzburg-Landau equation with and the incompressibility condition ∂zu + ∂ zū = 0, with z = x + iy and ∂ z = (∂ x − i∂ y )/2. Topological defects appear in the ψ-field as localized regions (of size ξ) where S vanishes and around which the nematic orientation θ winds by π. The method described in the next section allows us to derive the kinematics of the ±1/2 nematic defects from the dynamics of the zeros of the ψ-field.

Defect kinematics
The complex nematic order parameter ψ is always a single-valued function even in the presence of defects. This restricts any phase discontinuity in θ(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 β contained inside C. We consider only energetically stable defects of charge q β = ±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 [19,47,48]. 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 δ(ψ), which picks up the defect positions r β , with β = 1, 2, . . . the particle number index. The field variables map into discrete defect positions {r β }, where the Jacobi determinant D(r, t) = det(∂ψ) is with 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 ψ-field and readily obtained by taking time derivatives on both sides of equation (8) [40], with the result where the conservative current J (ψ) i is given by or in complex form The defect charge is related to the sign of the determinant evaluated at the defect position, q β = 1 2 sgn(D)| r=r β [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 ρ 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 v β =ṙ β : The defect velocity is therefore determined by the field D and its current density: which explicitly relates the defect velocity to derivatives of ψ evaluated at defect positions. In complex form 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 ψ in the vicinity of the defect.

Defect velocity
To relate the defect velocity to nematic distortions and to flow, we need to evaluate gradients of ψ at the defect core (equation (18)). Focusing on the βth defect, the order parameter can be written as where andθ is the smooth phase distortion due to all other defects and to boundary conditions. The core function S(|z|) has the generic asymptotic behaviors: S(|z|) ≈ |z| for |z| 1 and S(|z|) ≈ 1 for |z| 1. Notice that in equation (18), the density current and D-field, which depend on the derivatives of ψ, are evaluated at the defect position where S(|z| = 0) = 0. To leading order in the phase gradients it is then sufficient to consider only the linear profile of the core function near the defect. Including the full-nonlinear profile will introduce additional contributions from the core structure.
From equation (9), the determinant D at the defect position is simply D| z=z β = 2q β . The defect velocity from equation (18) is determined directly from the evolution of the ψ-field. At the defect core, we have ∂ t ψ| z=z β = [4∂ z ∂zψ − u∂ z ψ −ū∂zψ + λ∂zu] z=z β , since all the local terms in ψ vanish. As a result, the defect velocity does not depend on the local fluid vorticity and ordering potential, and is given by The nematic phaseθ β =θ(z β (t),z β (t)|{z γ (t),z γ (t)}), where γ = β labels all the other defects, is evaluated at the position of the tagged defect. To zeroth-order in activity α, the first term on the rhs of each of equation (21) represents the Coulomb interactions with the other defects. This is easily seen by noting that for α = 0 the steady-state nematic phase field can be written as the superposition of the phase induced by each point defect [28],θ where θ 0 is a fixed external phase. In particular, when the system is topologically neutral, θ 0 is equal to the far field orientation of the nematic texture. The first term on the rhs of equation (22) determines the Coulomb pairwise interaction force F β,γ between any two defects β and γ as The other two terms on the rhs of each of equation (21) represent the contribution to the defect velocity from the flow velocity u. The terms linear in u describes advection of both positive and negative by flows at defect cores. The terms proportional to flow gradients arise from the shear flow alignment, and drive defect motion in a direction determined by the flow gradient field relative to the local nematic orientation.

Defect polarity
It has been shown that defect polarity plays an important role in the dynamics of defects in active nematics [20,27,28]. The polarity e + β of the βth +1/2 defect is determined by gradients of the order parameter For a −1/2 defect, ∇ · Q vanishes at the defect core due the defect's trifold symmetry. The polarity of a −1/2 defect can, however, be defined as a triatic [31] or as a vector e − β that points along one of the three equivalent directions that define the defect's orientation [30], defined as where σ 3 is the Pauli matrix, In general the polarity of a defect depends on the orientation and position of all other defects, as well as on θ 0 . Near a +1/2 defect located at the origin and far from all other defects, the order parameter takes the form ψ + (z,z) = (zz) 1/2 z z 1/2 e 2iθ 0 and the polarization e + = 2 e 2iθ 0 is determined entirely by the boundary conditions, as shown in figure 1(a). Similarly, near an isolated −1/2 defect located at the origin, ψ − (z,z) = (zz) 1/2 z z −1/2 e 2iθ 0 and the associated polarization is e − = 2 1/3 e 2iθ 0 /3 (see figure 1(b)). Whereas, for a dipole consisting of a positive defect at z + and a negative defect at z − , we have Denoting by w = z − − z + = |w|e iφ the separation between the two defects, the respective polarizations depend both on the far-field boundary condition and the dipole orientation as given by e + = 2 e i[−φ+π+2θ 0 ] and e − = 2 1/3 e i[φ+2θ 0 ]/3 , as shown in figure 1(c).
We evaluate the dynamics of a given defect assuming thatθ(z,z; {z γ ,z γ }) is quasistatic, meaning that it is constant in time during the motion of the tagged defect or of any of the other defects. The defect polarization is then not an independent variable since its dynamics is determined by the dynamics of {z β (t)}. Nonetheless, it provides a useful effective degree of freedom for describing the dynamics of defects as quasiparticles, as well as the hydrodynamics of the defect gas [20,27].  3 (black arrows, see equation (27)) to the velocity of defects induced by shear flow is shown for specific configuration for a +1/2 (panel (a)) and a −1/2 defect (panel (b)). The orange arrows are the defect polarizations as defined in equations (24) and (25). The background arrows represent the flow velocity gradient.
For a positive defect, the time evolution of the defect polarity follows simply from taking the time derivative of equation (24), with the result [28] where v + β is the defect velocity and F β,γ the pair Coulomb force defined in equation (23). It is also useful to rewrite the equations governing the defect dynamics in real coordinates. Namely, the defect velocity from equation (21) in real coordinates reads as where E σ = σ 3 : E. The coupling to shear flow drives the ±1/2 defects to move in a direction controlled by the angle between their polarity and the flow gradient. For example, when the defect polarization points either along the direction of flow or along the direction of flow gradient, the coupling to shear flow drives the ±1/2 defects to move normal to their polarization, and in opposite directions, as shown in figure 2. Similarly the equation for the polarization of the q β = +1/2 defect takes the real space forṁ The first term on the rhs of equation (28) shows that the polarity of the +1/2 defect rotates at a rate dependent on the degree of alignment between defect velocity and the nematic distortion, as found in [20]. Equivalently, the net torque is given by the cross-product of the Coulomb force with the relative velocity of the tagged positive defect with respect to the rest of defects. A result of this form was recently obtained in reference [28] for the case of compressible flow. On the other hand, equation (28) holds for any flow condition since flow enters implicitly through the defect velocity. In the next section we obtain a set of closed equations for the defect positions and orientation by expressing the Stokes flow in terms of the defect degrees of freedom.

Defect-induced flow
In this section we obtain explicit expressions for the flow induced by specific defect configurations, which, in turn, drives defect dynamics via equation (27). 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 (∇ · 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 reference [18].

Single defect in compressible flow
For compressible flows given only by the active stress u = 2α∂ z ψ, the flow velocity at a positive defect is simply proportional to the defect's polarization e + β , while it vanishes at a negative defect. The flow alignment term in equation (27) can be expressed in terms of the rotated phase gradients using equation (19) for the profile of ψ near the defect core. This yields the defect propulsion velocity obtained in previous work [18-21, 27, 28], with defect dynamics governed by where K λ = 1 + λα/2 is the effective elastic constant (in dimensionless units) incorporating the effect of flow alignment. The dynamics of the polarization can be obtained from equations (28) and (29) aṡ and has precisely the form of the active torque given in reference [20], but with K replaced by K λ . Note that K λ can change sign for extensile activity (α < 0), resulting in a change of sign of the active torque. Such corrections are, however, nonperturbative in activity.

Incompressible flow
For incompressible flows, evaluating the flow velocity requires solving the Poisson equation for the pressure: where the source term acquires contributions from the topological defects and any smooth deformations of the order parameter, such as kink walls. The solution of equation (31) can be written as the convolution of the logarithmic Green's function G(r) = 1/(2π)log(|z|) with the source term: The corresponding incompressible velocity field follows from equation (5): where arises from the pressure gradient. The first term on the rhs of equation (33) is the defect self-propulsion, which has the same form for both compressible and incompressible flows. Solving for the pressure field induced by an arbitrary defect configuration quickly becomes analytically intractable. The calculation can be carried out, however, for simple configurations, such as a single defect far from all others ('isolated' defect) and a defect dipole. These cases are discussed below, with technical details given in appendix A. Isolated defect. We consider a single positive defect located at the origin with a fixed external phase θ 0 . The distortions generated by the defect produce a surrounding flow determined by equation (33), with the induced pressure gradient calculated explicitly in appendix A.1. The velocity from the active stress u = 2α∂ z ψ follows from equation (19). The resulting incompressible flow field at distances |z| larger than the core size is then The first term on the rhs of equation (35) is the same as obtained for compressible flow and remains finite at the defect core, contributing to the defect's propulsive velocity. The two terms in square brackets arise from pressure gradients. The second of these two terms yields a correction to the flow at the core that reduces the defect's propulsion. The net self-propulsion velocity of the positive defects due to both active stress and pressure is then given by  where e + is the defect polarity from equation (24). Pressure gradients drive defect motion in the direction opposite to that determined by active stresses, hence reducing the magnitude of the defect's propulsion velocity as compared to compressible flows (equation (29)). The net propulsion remains, however, in the direction determined by the sign of activity. Finally, the shear flow due to pressure depends on derivatives of the I-function, which vanishes at the positive defect (see appendix A.1), and so has no effect on an isolated positive defect. The relation between the defect's self-propulsion velocity and the defect polarity remains valid for a slowly-varying orientationθ(z,z). The net velocity of an isolated positive defect in an incompressible active nematic is then given by z z e 2iθ 0 , where again the first term is the flow induced by the active stress (which vanishes at the core) and the terms in square brackets arise from pressure gradients. The net velocity of an isolated −1/2 defect is then the same as in compressible nematic and determined solely by the local nematic distortion, On the other hand, when other defects are present, the I-integral is finite at the defect core, as we see next.
Similar expressions for the incompressible flow field around an isolated defect have been derived in reference [26] by solving for the Green's function in the real space.
In figures 3 and 4 we compare the flow field induced by a single defect in a compressible (panel (a)) and an incompressible (panel (b)) active nematic. In the compressible case, the flow is determined entirely by the active stress (u = 2α∂ z ψ). In the limit of zero viscosity considered here, the vorticity field remains the same with or without the incompressibility constraint, since the pressure gradient is curl-free.

Defect dipole
Another configuration for which one can obtain an analytic solution is a pair of oppositely-charged defects separated by w = z − − z + and embedded in an otherwise uniform nematic orientation θ 0 (see figure 5(a)). The details of the evaluation of the integrals arising from pressure gradients are presented in appendix A.2. The flow velocity at the defect cores are given in terms of the relative distance w as The first term on the rhs of equation (40) is the self-propulsion of the +1/2 directed along its polarity e + = 2 e 2iθ 0 w/w. Note that the polarity of the +1/2-defect depends on the orientation of the negative defect through the factor w/w. This dependence survives even when the −1/2 defect is infinitely far away (|w| → ∞). The second term arises from the non-local hydrodynamic interaction between the two defects due to pressure gradients. The velocity of the −1/2 given in equation (41) is determined entirely by hydrodynamic interactions. The different geometric structure of the two defects renders these interactions non-reciprocal, with the velocity of the negative defect five times smaller than that of the positive one. In the limit where one defect is moved to infinity relative to the other, the hydrodynamic interaction terms vanish and we recover zero drift for the isolated −1/2 defect and the self-propulsion velocity for the isolated +1/2 defect. Substituting the explicit expressions for the incompressible flow in terms of the defect positions and polarity in equation (21), we obtain wherer ≡ (r − − r + )/|r − − r + |. The first terms on the rhs of each of equation (42) represent the familiar Coulomb interactions among the two defects. The other terms arising from corrections to active flows induced by pressure gradients are non-reciprocal in nature and yield forces normal to the line joining the two defects. We emphasize that these forces are distinct from non-reciprocal forces obtained in reference [28] from a perturbative analysis of multidefect interactions in a compressible active nematic. These forces are a new result of our work. Finally, notice that there is no need to evolve the defect polarization separately, since it can be reconstructed directly from the defect positions and the uniform nematic phase θ 0 . Figure 5(b) shows the trajectories of the ±1/2 defects determined by equation (42). For comparison we also plot the defect trajectories in a compressible flow field driven by the active stress alone. Incompressibility reduces the transient active torque and the velocity of the positive defect, while slightly increasing the velocity of the negative defect. A generalization to many defects presents many challenges since the pressure field encodes all the non-local interactions between defects through integrals with many more complicated branch cuts and singularities.

Defect hydrodynamics
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 reference [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,ρ andn, their associated current densities,Ĵ (ρ) andĴ (n) , and defect polarization densityp: 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({r β }), 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 u(r) = û(r, t; {r β }) and a coarse-grained nematic texture v n (r) = ∇θ(r, t; {r β }) . 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 where for any vector V = (V x , V y ) we denote the clockwise rotation corresponding to the multiplication with the Levi-Civita tensor ij as V ⊥ = · V = (V y , −V x ). These equations couple defect dynamics to the coarse-grained flow and nematic texture, with f βγ = r βγ /|r βγ | 2 , which r βγ = r β − r γ , which is the interaction force which relates the nematic distortion with defect density as is the interaction force from the defect at position r γ on the defect at position r β . The phase gradient v n contains the singular phase induced by all defects and relates to the macroscopic charge density ρ through Stokes' theorem applied in equation (6), This determines v n up to a potential vector field which, in the quasistatic approximation, satisfies The macroscopic fluid flow u follows the Stokes equation (3), but with the active stress expressed in terms of the density and polarization density of positive defects, Incompressibility additionally requires ∇ · 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, equation (43), using the defect dynamics given in equation (49), and applying the coarse-graining procedure. Neglecting defect creation and annihilation, the equations for the defect density fields ρ and n take the form In the presence of defect annihilation/creation, the rhs of equation (57) needs to be supplemented by the kinetic rates describing these processes. The macroscopic current densities are given by = 2v ⊥ n n(r, t) + uρ(r, t), The evolution of the microscopic polarization density follows from taking the time derivative of equation (47), using equation (49) and expressing the sums over defects in terms of the microscopic densities and their current densities, with the result where U = u + 4v ⊥ n is the net flow including the one generated by nematic distortion. Upon coarse graining, we assume p (r)ρ(r ) ≈ p(r)ρ(r ) and p (r)n(r ) ≈ p(r)n(r ), and obtain the continuum equation ∂ t p + U · ∇p = −2(v n · U)p ⊥ − [∇ · U]p − 2p ⊥ (r) dr f ⊥ (r − r ) · J (ρ) (r ).
The polarization equation does not contain terms controlling growth or decay of the magnitude of p nor any stiffness/diffusive terms because we assumed |e β | = 2 (not normalized) and did not include noise in the defect dynamics. The first term on the rhs of equation (60) controls changes in the direction of polarization and arises from the active torque in the evolution of e β . 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, equation (53) implies that there are also local changes in polar order depending on the charge density field, ρ. Finally, the last term on the rhs of equation (60) is a new nonlocal correction to the torque that arises from the long-range nature of nematic distortions.
It is useful to compare equation (60) for the polarization density to the one obtained in reference [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 reference [27] for compressible fluids. This can be seen by noting that in reference [27] the polarization equation was obtained by eliminating u from the defect dynamics before coarse-graining, using u → αe β in equation (59). If we use this in equation (59), and then coarse grain, using the same moment closure as in reference [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 u but also by the Magnus-type force proportional to v ⊥ n . The fact that only positive defects act as a source of flow in the Stokes equations for 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 v n .

Conclusions
To summarize, we have analyzed the effect of arbitrary fluid flow on the motion of ±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.