Flow around topological defects in active nematic films

We study the active flow around isolated defects and the self-propulsion velocity of +1/2 defects in an active nematic film with both viscous dissipation (with viscosity η) and frictional damping Γ with a substrate. The interplay between these two dissipation mechanisms is controlled by the hydrodynamic dissipation length ℓd=η/Γ that screens the flows. For an isolated defect, in the absence of screening from other defects, the size of the shear vorticity around the defect is controlled by the system size R. In the presence of friction that leads to a finite value of ℓd, the vorticity field decays to zero on the lengthscales larger than ℓd. We show that the self-propulsion velocity of +1/2 defects grows with R in small systems where R<ℓd, while in the infinite system limit or when R≫ℓd, it approaches a constant value determined by ℓd.


Introduction
Active matter consists of collections of individuals that dissipate energy taken from the environment to generate motion and forces and self-organize into a rich variety of ordered phases. Many active systems exhibit nematic order interrupted by orientational defects and advected by spontaneous flows driven by intrinsic activity of the self-propelled individuals. This behaviour is found in reconstituted systems, such as mixtures of cytoskeletal filaments and motor proteins [1][2][3][4], bacterial suspensions [5,6] and cell sheets [7,8], as well as synthetic systems, like vertically vibrated layers of granular rods [6,9]. A central feature of active nematics is the feedback between active stresses, which distort orientational order and the spontaneous flow generated by such distortions. In hydrodynamic descriptions [6], the active stress σ a ij exerted by elongated active entities on the surrounding fluid is proportional to the nematic order parameter tensor Q ij , namely σ a ij = α 0 Q ij [10,11]. The activity coefficient α 0 embodies the microscale biomolecular processes that convert chemical energy into mechanical forces, and depends on the concentration of active entities, which in general may vary in space and time [12,13]. The sign of α 0 distinguishes between contractile (α 0 > 0) stress generated by 'puller' swimmers, such as the algae Chlamydomonas, versus extensile (α 0 < 0) stress generated by 'pusher' swimmers, e.g. most flagellated bacteria. Its magnitude controls the strength of the active flow. Fluctuations in orientational order yield active stresses and associated flows, which can in turn enhance the orientational distortions. The resulting feedback loop destabilizes the nematic order, driving the system to a state of self-sustained spatio-temporally chaotic flow, with proliferation of topological defects, and termed active turbulence [14,15].
The lowest-energy orientational defects in nematic films have half-integer topological charge and opposite sign. The +1/2 defects have a comet-like shape, while the −1/2 defects have a tri-fold symmetry (figures 1 and 2). Defects strongly disrupt orientational order and induce longrange nematic distortions. In active systems, such distortions generate flows with symmetry and profiles controlled by the defect geometry. The nematic distortion created by a +1/2 defect yields an active flow that is finite at the defect core. A +1/2 defect then rides along with the flow it itself generates, behaving like a motile particle with a non-vanishing self-propulsion velocity v a + , even in the absence of external drive [12,16]. On the other hand, the active backflow generated by a −1/2 defect vanishes at the core due to the defect's threefold symmetry (figure 2).
Flow streamlines (white arrow) around a −1/2 defect for α < 0 obtained from (a) full solution and (b) asymptotic one. The nematic director field is shown in black lines and the background colour map denotes vorticity. To show the structure of the near-field, the vorticity scale is saturated at ±0.2 in (a) and at ±0.1 in (b). (c) Cross section of vorticity obtained from the exact solution (solid blue line) and the asymptotic limit (dotted black line) at x = 0 as a function of y. Note that the vortices changes sign rapidly at origin due to its multi-valued phase (see equation (4.6)). (d) Cross section of the velocity obtained from the exact solution (solid blue line) and the asymptotic limit (dotted black line) at y = 0 as a function of x. (Online version in colour.) −1/2 defects behave like passive particles and have no spontaneous motility in the absence of external driving. A simple estimate demonstrates that v a + is directed along the polar axis of the +1/2 defect and is proportional to the activity α 0 . In an extensile medium +1/2 defects self-propel in the direction of the head of the comet, while in a contractile system they move towards the comet's tail [12,16,17]. The direction of motion of +1/2 defects can then be used as a metric for determining the nature of active stress in the system. Such measurements have, for instance, revealed the surprising dominance of extensile stresses in confluent tissue composed of tightly bound contractile individual cells [18][19][20][21][22].
The flow generated by defects and the resulting propulsive speed of the +1/2 also vary depending on the dissipative processes at play in the system and the role of fluid incompressibility. Specifically, important differences exist between 'dry' systems, where dissipation is dominated by friction Γ with a substrate or an external medium [23,24] and 'wet' systems, where dissipation is mainly controlled by viscosity η, resulting in long-range hydrodynamic effects [14,17,24,25]. In incompressible wet systems, activity is also a source of pressure gradients, which in turn contribute alongside with the nematic distortion to the selfmotility of positive defects. In the limit of viscous dominated flows with no friction with the substrate, the self-propulsion speed scales as |v a + | ∼ (|α 0 |/η) , where is a length scale given by the system size for an isolated defect [17] or by the mean separation between defects, which is, in turn, controlled by the active length scale a = √ K/|α 0 |, with K the nematic stiffness [17]. In overdamped (dry) systems, where viscosity is negligible compared with frictional damping with the substrate, |v a   calculation of the active flows associated with defect configurations and of the propulsive speed of the +1/2 defect that bridges between the two limits is, however, not available. The need for such a calculation is further motivated by recent work that has shown that tuning frictional damping relative to viscous dissipation leads to different dynamical regimes and ordering behaviour of interacting defects [27].
In this paper, we present a detailed calculation of the flow around isolated ±1/2 defects and of the defect's self-propulsion velocity in an incompressible nematic film. We incorporate both viscous dissipation and frictional damping and examine the interplay between the two, as well as the long-range hydrodynamic effects arising from incompressibility. We evaluate the +1/2 self-propulsive speed |v a + | as a function of the hydrodynamic dissipation length d = √ η/Γ , which measures the competition between viscous dissipation and frictional damping. The result is summarized in figure 3. When dissipation is controlled by friction ( d ξ ), one recovers the simple dimensional estimate |v a + | ∼ |α 0 |/(ξΓ ). We show, however, that to obtain this result it is not sufficient to consider the far flow field which diverges near the defect, but one must resolve the full flow field near the defect core. On the other hand, when viscous stresses dominate, the defect propulsive speed depends on the order of limits. If Γ = 0 from the outset, then a simple estimate yields v a x ∼ r due to the long-range nature of defect distortions. This limit, however, corresponds to a 'floating' layer and does not describe experimental situations where the active nematic film is supported by a substrate [25] or in contact with other fluids. It has been argued before that this unbounded growth should be cut off either by the system size or by the defect separation [17]. Our work shows that a finite friction cuts off the large-scale divergence of the defect self-propulsion speed at the scale d , with |v a + | ∼ (|α 0 |/η) d in the limit ζ = d /ξ 1, where viscous dissipation exceeds frictional drag and provides an analytical expression for the defect self-propulsion over all values of friction and viscosity. We find that the structure of the flow field around a defect is also affected by the competition between viscosity and friction. At distances large compared with d , the flow velocity decays in the far-field as ∼ 1/r, due to friction with the substrate [25]. At distances smaller than d , viscous dissipation dominates and smooths out the velocity field near the defect core. Our work is relevant to defects in thin film of microtubule nematics on a substrates, as well as to dense cell layers.
In §2, we describe the hydrodynamic model. In § §3 and 4, we provide analytical derivations of closed expressions for the velocity and pressure fields induced by ±1/2 defects in an infinite system. One implication of the long-range interactions present in active nematics is that there are strong finite-size effects on the single defect flow field. This is discussed in §5, where we compare the analytical predictions with numerical integration of the Stokes equations in a disc of finite radius. Finally, the main results are discussed with concluding remarks in §6.

Hydrodynamic model
We consider a hydrodynamic model of an active nematic that couples flow velocity u(r) to the nematic order parameter Q ij = S(n inj − 1 2 δ ij ), where S quantifies the degree of order and n(r) = (cos θ(r), sin θ (r)) is the orientational director field with head-tail symmetry. In the simplest formulation, we consider that the Q-tensor is a minimizer of the de Gennes-Landau free energy [6] with isotropic elastic constant K > 0 and g the strength of the local ordering potential. The uniform nematic ordered state corresponds to S 2 0 = 1. The flow field satisfies a Stokes equation that balances forces on a fluid element, given by [6] ( where Γ is a friction coefficient per unit area, η is the dynamic viscosity and α 0 is the activity parameter, with dimensions of stress. For simplicity, we neglect the elastic stress as being of higher order in the gradients of Q compared with the active stress and a more important contribution for nematic textures with many defects. Here, we consider the flow field generated by an isolated ±1/2 defect embedded in an otherwise uniform nematic field. In two dimensions, the traceless Q-tensor has two independent components and can be represented equivalently as a complex scalar order parameter ψ = Q xx + iQ xy . The configuration of a defect located at the origin can be written in terms of the ψ-field as ψ(r) = S(r)e 2iθ(r) , where r ≡ |r|. The detailed form of core function S(r) depends on the specific terms retained in the free energy, but it has the important generic asymptotic behaviours that S(r) → 1 for r ξ and S(r) ≈ ar/ξ when r → 0, where ξ = K/g is the coherence length that sets the scale of the defect core and a is a numerical constant O(1). Below we set a = 1, without loss of generality. The coherence length provides an ultraviolet cutoff to separate inner core-solution from outer-core solution. On long distances, the nematic orientation is a potential field that has a branch cut starting at the origin where there is an isolated defect of charge q = ±1/2 and can be written as [24,28] where θ 0 is the uniform background orientation. Without loss of generality, we set θ 0 = 0. We rescale the Stokes equation in units of the nematic relaxation time τ = γ /g (where γ is the inverse of the rotational diffusivity) and the coherent length ξ , such that the dimensionless momentum equation takes the form where F ± = α∇ · Q is the active force generated by a defect. The rescaled activity and pressure are given by The components of the Q tensor for an isolated +1/2 defect are given by Q xx (r) = S(r)(x/r) and Q xy = S(r)(y/r). The active force density then reduces to Similarly, for a negative defect Q xx = S(r)(x/r) and Q xy = −S(r)(y/r), corresponding to an active force density given by 2xy r e y , r 1. (2.6) The solutions for the flow velocity and pressure can be written in terms of the corresponding Green functions as and where u a and u p are the contributions to the flow velocity induced by the active stress and pressure gradients, respectively. Note that the latter also depends (indirectly) on activity. In the limit of no friction, equations (2.7) and (2.8) reduce to equations (3.7) and (3.8) of [17].

Positive nematic defect in an infinite system (a) Defect self-propulsion
The net active flow at the defect core acts as an advective velocity that propels the defect with a velocity v a , which in turn is controlled by both the active stress and pressure gradients. Thus we write v a = u a (0) + u p (0). The flow induced by the active stress at the origin is given by equation (2.7) evaluated at r = 0. The y-component vanishes due to symmetry considerations, and the x-component is given by where ζ = d /ξ , K n (x) are modified Bessel functions and L n (x) modified Struve function. The integral determining the pressure field given by equation (2.8) can be performed by a mapping to complex coordinates (x , y ) → (w,w), (x, y) → (z,z) and then using the substitution to polar coordinates w = r ŵ,ŵ = e iθ . This yields with γ a contour of unit radius centred at origin. The pole atŵ = 0 is always inside the unit disc |ŵ| < 1, whereas the poles atŵ = zr −1 andŵ = r z −1 are inside the unit disc when |z| < r or |z| > r , respectively. The contour integrals are then evaluated using the residue theorem. Consequently, the defect self-propulsion induced by pressure gradient has only an x-component, which counteracts that induced by the active stress, and given by Combining these results, we find that the self-propulsion velocity of an isolated +1/2 defect oriented along the x-axis is v a = v a xê , where v a x has the following scaling form: v a x = αF(ζ ), (3.5) where When ζ 1, we can simplify the expression by expanding in powers of ζ −1 , and, to leading order, we obtain, where γ ≈ 0.577 is the Euler-Mascheroni constant. Similarly, we also take the other limit ζ 1, where the scaling function approaches a constant value. The dependence of the scaling function F on ζ is plotted in figure 4 and its asymptotic scaling at ζ 1 as F ∼ ζ −1 is included as the dotted line. We can discuss the implications of these results better when we use dimensional quantities and write the asymptotic behaviour of the self-propulsion speed as As anticipated from dimensional analysis, v a x ∼ (α 0 /Γ ξ ), in the overdamped limit where dissipation is controlled only by frictional drag [23,24,26]. In the underdamped limit, where the effect of drag is much smaller than viscous dissipation, hydrodynamic lengthscale becomes important in screening the divergence of the self-propulsion speed with system size, such that v a x scales instead as v a x ∼ α 0 / √ ηΓ . In this case, the self-driven motion of +1/2 defect is reduced by both friction and viscosity. As discussed in the introduction, the presence of a finite drag always cuts off the large-scale divergence of the speed of a single defect obtained in a purely viscous two-dimensional layer at the dissipation length d . When the flow equations for a thin nematic film of thickness h on a substrate are derived via a lubrication approximation, the effective friction coefficient relates to the film thickness and the viscosity of the substrate bulk fluid (oil), and scales as Γ ∼η/h 2 [29]. A more detailed calculation relevant to active microtubule suspensions confined between water and oil shows that the bulk viscosity plays an important role as an additional source of dissipation in the nematic layer affecting the individual defect self-propulsion [30], as well as the vortex statistics in the active turbulence regime [31]. Note that [30] shows that the +1/2 defect speed decays algebraically with the bulk oil viscosity (that controls the drag) in the regime where the flow dissipation comes from the viscous dissipation in the nematic layer, consistent with our formulation. When the flow dissipation is dominated by the oil bulk viscosity, there is, however, a logarithmic decay with increasing oil viscosity and, indirectly, drag.

(b) Flow field away from the defect
Outside the core, we treat the defect as a point source. From symmetry considerations, the flow velocity due to σ a is again non-zero only along the x-direction and it is given by The flow velocity associated with pressure gradients is finite also in the y-direction and it is given by To evaluate this integral, we use a complex representation u = u x + iu y and evaluate the resulting contour integrals as shown in appendix A where we express them in terms of complete elliptic integrals of first and second kind. We further use the power series representation of these elliptic integrals, which allows us to write the active fluid velocity as a series expansion in integrals over the zeroth-order modified Bessel function, namely (3.12) The K 0 (x) integrals are computed in appendix B. After some mathematical manipulations the velocity reduces to −(4n + 1)(4k + 3) (n + k + 1) 2 (2n − 1 − 2k) 2 (3.14) and The corresponding vorticity is given by (3.16) Both velocity and vorticity are shown in figure 1.

(i) Asymptotic far-field flow
The flow field greatly simplifies in the far-field r/ζ 1, corresponding to distances much larger than the hydrodynamic dissipation length. Then, the second term in equation (3.12) vanishes due to the exponential decay of the Bessel function. In the first integral, we can replace the upper limit r with ∞ and perform it analytically with the result given as where we have kept the two first terms in the expansion. The slow 1/r-decay term in equation (3.17) is independent of viscosity η and identical to the one derived in Ref. [24] in the frictiondominated regime. Corrections due to viscosity give rise to faster 1/r 3 decay. The corresponding far-field vorticity is (3.18) The far-field solutions are singular at the origin, which is not the case for the full series solution that resolves the near core field. This is demonstrated visually in figure 1c,d where we plot cross sections of the velocity and vorticity profiles for both the full solution and the far-field solution. The form of the expressions makes it natural to scale the position, velocity and vorticity with ζ , ζ /|α| and ζ 2 /|α| respectively. The only free parameter is then the sign of α. Panels (a) and (b) show the flow streamlines and the vorticity field in the background for the full and the far-field solutions, respectively, for an extensile system (α < 0). The velocity magnitude is highest near the defect core and decays as a power law following the far-field asymptote. The velocity streamlines point towards the defect in the right half-plane, and away from the defect in the left half-plane. For positive α, the flow direction is reversed. In an infinite system, the flow streamlines around an isolated defect are not closed. On the other hand, as discussed later, in bounded domains, the system size controls the size of the eddies formed around the defect. For more realistic configurations with many defects, the system size is typically replaced by the mean defect separation. It may be that other intrinsic length scales controlled by elastic stresses are also important in stabilizing finite-size vortices. These effects are left for future investigation.

Negative nematic defect in an infinite system
By similar calculations as in §3, we find that the velocity induced by the active stress at the position of the negative defect vanishes as expected from symmetry consideration. After performing the integral in the complex plane and subsequently integrating over the integrand with the Bessel function, we determine the pressure field induced by the −1/2 defect vanishes inside the defect core and non-zero outside given by and its gradient vanishes at the origin, hence no advective pressure-flow of the negative defect. Thus, an isolated −1/2 defect is stationary in a uniform nematic field, regardless of activity.
(a) Flow field away from the defect The flow field induced by the −1/2 defect can also be expressed analytically as a series expansion of the elliptic integrals as detailed in appendix C, with the resulting expression of the velocity field in the complex representation u − = u − x + iu − y given as The integrals over the Bessel functions are evaluated in appendix B, and the final expression is then given as with the coefficients and (4.5) The corresponding vorticity field as function of the polar coordinates follows as: (4.6)

(b) Asymptotic far-field flow
As with the +1/2 defect, the far-field asymptotic flow is dominated by the leading order terms in the expansion, which can also be computed directly from equation (4.2) in the limit of r/ζ → ∞. The result of this calculation is that As for equation (3.17) the 1/r term here was also obtained in [24]. The vorticity related to this velocity is In this asymptotic approximation, the flow field is singular at the origin. This singularity is however lifted by the higher order terms in the series expansions, so that the exact flow is smooth everywhere. Figure 2 shows the flow streamlines with the vorticity field as the colour map for the asymptotic (in a) and the exact solutions (in b), with the values scaled in the same way as for figure 5. Cross sections of the vorticity and velocity at y = 0 are plotted in (c,d) showing the singular behaviour of the asymptotic approximation at the origin, while it captures very well the far-field behaviour. The plots correspond to an extensile system with α < 0. The sin(3φ) factor in the vorticity divides the plane in six regions where the sign of the vorticity is altered and making it multi-valued at the origin. The size of the velocity is zero at origin as we discussed above. It increases a bit outside before it starts to decay with increasing r following the far-field asymptotic behaviour.
As for the +1/2 defect, the flow streamlines never closed in an infinite system, thus there are no finite size vortices. In the next section, we discuss how the picture changes once the defect is placed in bounded domain.

Isolated defect in a bounded active nematic
The problem of finding the flow field around defects in a bounded domain is challenging to solve analytically. Thus, we resort to numerical solutions of the Stokes flow given by equation (2.4) in a disc of radius R using finite-element methods and homogeneous boundary conditions (zero velocity). In addition, we use the simplification that a single defect is imprinted in an uniform nematic field, while the changes in the nematic orientation induced by confinement are ignored [17].
The Stokes flow equation (2.4) is solved with FEniCS using Taylor-Hood elements, which are quadratic for the velocity and linear for the pressure and vorticity [32,33].  Figure 5 shows the flow streamlines induced by a single +1/2 (a,b) and −1/2 (c,d) defect in a disc of radius R for an extensile system for η = 0 and Γ = 0. The left and right columns correspond to R = 1 and R = 50, respectively. In a bounded system, the vortical flows around each defect span the system size, as also reported in [17] for Γ = 0. However, due to friction with the substrate, the flow decays on length scales larger than d . This is evident by comparing the values in the far-field of vorticity in the left and right columns from figure 6, corresponding to R = 1 (in units of d ) in (a,c) and R = 50 (in units of d ) in (b,d). We note that the centre of a vortex is not fixed at the maximum of the vorticity. This is due to the fact that the ±1/2 defects generate shear flows that localize shear vorticity next to the defect cores. However, unlike curvature vorticity in rotating flows which peaks at the vortex core, shear vorticity is not necessarily an indication of the presence of vortices or their location. In fact, with increasing R, the flow gradients near the defect cores become sharper, the streamlines near the cores are 'stretched' in the radial direction, and the eyes of vortices move further from the origin. In the limit R → ∞, we expect vortices to get stretched out so that flow streamlines close at infinity, and we recover the analytic flow profiles.
In figure 7, we compare the cross-sectional profiles of velocity and vorticity obtained from the numerical solution for a large system, R = 50 ( d ) to the analytical solutions. The analytical solutions is obtained by truncating the summation in the full solution at n = 5000 and k = 500 up to r = 15 and then using the asymptotic solution for r > 15. The plots of velocity in panels (a,c) show that the numerical and analytical solutions agree very well close to the defect cores, but deviate from each other near the boundary. This is due to the imposed boundary conditions of the vanishing velocity field. The vorticity in panels (b,d) agrees well in the entire domain, with a  small boundary effect due to vanishing velocity and vortices spanning the system size. This effect is perhaps more visible for the negative defect and decreases with increasing R.
The self-propulsion speed v a x of the +1/2 defect is also affected by the system size. If Γ = 0 at the outset v a x ∼ R, as noted in [17]. Frictional damping screens out this divergence, yielding the finite value given in equation show the linear scaling with R in the limit of zero friction. We note that viscosity η determines the slope for R dependence in small systems, while friction Γ controls the cross-over to the intrinsic constant speed. Note that the asymptotic constant values of v a x agree very well with the analytical prediction at ζ 1 because in the numerical computations the vortex core is actually set to zero (hydrodynamic regime with S = 1). For comparison, we also show in figure 8b the defect propulsion speed in the absence of friction Γ = 0 from the outset, where the speed increases linearly with the system size. The dotted black line represents the analytical prediction as found in [17].

Conclusion
In summary, we have evaluated the flow field induced by an isolated ±1/2 defect in an incompressible active nematic film on a substrate both for an infinite system and a finite-size disc. While the self-propulsion speed of a +1/2 defect diverges with system size for an isolated film, we show analytically that the presence of finite substrate friction Γ cures this divergence resulting in a finite speed v (a) x ≈ (π/4)(α 0 / √ ηΓ ) = (π/4)(α 0 /η) d that increases with the hydrodynamic dissipation length d . This is also confirmed numerically in a finite disc with R > d . For small discs with R < d , the active speed scales instead linearly with R. Stable shear vortical flows are formed around the defects. In finite systems, the size of the flow vortices is controlled by the dissipation length d , hence spans the whole system if d > R. The eye of the vortices shifts away from the defect core with increasing R. For infinite-size systems, the flow streamlines close at infinity as predicted by the far-field analytical solution. In the same limit, we showed that the absolute value of the velocity decreases as 1/r for distances that are large compared with the dissipation length scale, in agreement with previous studies. The 1/r far-field decay of the flow created by defects may seem surprising as it suggests that a defect acts like a point force. This behaviour arises from the long-range nature of the distortion of the texture created by defects. When other defects are present (as required in the plane to guarantee zero net topological charge), this decay is cut off by the defect separation. In finite domains, it is cut off by the system size. The 1/r decay indicates, however, that a multi-defect approach is needed to describe the defect gas, as attempted in [34,35].
In this work, we have neglected the effect of the elastic stress. An interesting extension would be to study the effects it would have on the flow field, and also considering the effect of having multiple interacting defects.
Data accessibility. This is primarily theoretical work and does not have any experimental data. The computational data and codes for FEniCS are available on GitHub: https://github.com/jonasron/Defect-Flows. Note that all branch pointsẑ = 0, −a and −b are now located on the real axis. We now have to consider what branch cuts we want to use to perform this integral. We consider the integral over the domain as shown in figure 5. Here, we have cut out a area around the branch cut in order to avoid problems. The key hole consists of a circle C 1 with radius around −a, the circle C 2 around the origin and the lines connecting them which is above or below the real line as shown in figure 5. Since there are no poles in the domain between the two contours, the integral of them has to be the same [36]. The contour integral becomes When → 0, the integrals over C 1 and C 2 disappears. The integral above the real line is just above the branch cut and therefore positive, while the one below is negative. We therefore get