Dynamics of active nematic defects on the surface of a sphere
DDynamics of active nematic defects on the surface of a sphere
Yi-Heng Zhang ( 张 一 恒 ), ∗ Markus Deserno, and Zhan-Chun Tu ( 涂 展 春 ) Department of Physics, Beijing Normal University, Beijing 100875, China. Department of Physics, Carnegie Mellon University,5000 Forbes Avenue, Pittsburgh, PA 15213, USA (Dated: July 21, 2020)A nematic liquid crystal confined to the surface of a sphere exhibits topological defects of totalcharge +2 due to the topological constraint. In equilibrium, the nematic field forms four +1 / de-fects, located at the corners of a regular tetrahedron inscribed within the sphere, since this minimizesthe Frank elastic energy. If additionally the individual nematogens exhibit self-driven directionalmotion, the resulting active system creates large-scale flow that drives it out of equilibrium. Inparticular, the defects now follow complex dynamic trajectories which, depending on the strength ofthe active forcing, can be periodic (for weak forcing) or chaotic (for strong forcing). In this paper wederive an effective particle theory for this system, in which the topological defects are the degrees offreedom, whose exact equations of motion we subsequently determine. Numerical solutions of theseequations confirm previously observed characteristics of their dynamics and clarify the role playedby the time dependence of their global rotation. We also show that Onsager’s variational principleoffers an exceptionally transparent way to derive these dynamical equations, and we explain thedefect mobility at the hydrodynamics level. I. INTRODUCTION
Topological defects show up in a surprising variety ofordered systems, and whenever they do give rise to fasci-nating emergent physics. They have been observed in thestudy of liquid crystals [1–4], crystalline solids [5–7], cellassemblies [8–10], superfluid vortices [11–13], magneticskyrmions [14–22], and cosmology [23–25]. Since defectscannot continuously vanish (they typically only pair-create or -annihilate), they constitute long-lived markersof the field that forms them, and hence their dynamicscan provide deep insights into the long-scale time evo-lution of such systems [26]. In recent experiments, Ke-ber et al. [27] have fabricated a spherical nematic byconfining kinesin-microtubule bundles onto the surfaceof a spherical lipid vesicle, and adding ATP to the ki-nesin motors renders this system active, i.e. , self-driven.The ATP-induced activity can drive the system far froma static equilibrium state, leading to a novel defect dy-namics in which the active flow competes with dissipa-tive relaxation processes, as well as elastic forces thatarise from the nematic field itself and the curvature ofthe substrate. Significant theoretical understanding hasbeen gained how such active driving forces affect the de-fect dynamics in a planar nematic liquid crystals [28–36];in contrast, a comprehensive theoretical framework forthe topologically distinct case of spherical confinementlags noticeably behind.The motion of defects can be regarded as a particularway in which the active nematic field collectively mani-fests itself: the field moves, and the defects “ride” withit. Hence, in order to understand this motion, it is nat-ural to start from the traditional nematic hydrodynamic ∗ Corresponding author; [email protected] equations [37–41], amended by a suitably chosen activestress [42–46]. In this way, one arrives at a continuumtheory that describes the evolution of the nematic orderand its underlying flow velocity, and one may therebypredict how the defects are transported along. We alsoemphasize that the curvature of the substrate can playan important role in the active nematic hydrodynam-ics, not only the intrinsic curvature of a sphere, but alsothe extrinsic curvature, like the curvature of a cylindri-cal surface [47, 48]. While conceptually fairly clean, theresulting partial differential equations do not afford illu-minating analytical solutions, and therefore a substantialamount of work has been dedicated to numerically solv-ing them, which has indeed provided much novel insightinto the collective behavior of spherical active nematicsystems [49–51].However, the topological nature of this problemstrongly suggests an alternative viewpoint. Recall thatdefects cannot continuously vanish or arise. Their topo-logical discreteness “protects” their existence, or, inphysics parlance, it equips them with a conservation law .More precisely: since a defect’s topological index canonly in finite steps (say, integral, or half-integral), defectsrequire other defects to change their overall number—meaning that they typically create and annihilate inpairs. This, of course, equips them with a property wecommonly associate with elementary particles (such aselectrons), which can be created or destroyed only inparticle-antiparticle pairs. Their longevity suggests thatit should be expedient to formulate an effective theorythat conceives of these discrete particles as the essen-tial degrees of freedom, which in turn interact by forcefields that emerge as “leftovers” of the originally stronglycoupled field. In our context, this perspective suggeststhat we view topological defects in a nematic as particles,whose dynamics is given by the force balance betweenthe effective friction and the elastic interaction. Such a r X i v : . [ c ond - m a t . s o f t ] J u l an effective theory (in soft matter language more usu-ally called “coarse-grained”) has indeed been successfullyused to describe the defect dynamics in planar nematics[52–66]. Since the Euler characteristic of a (simply con-nected) plane is zero, the total sum of all indices of defectsis zero, so that the number of defects equals that of the“antidefects” (meaning that there are exactly as many +1 / defects as − / defects). Hence, the characteris-tic physical processes are the creation and annihilationof defects in pairs out of and into an otherwise struc-tureless (“trivial”) vacuum ground state—a scenario thatis indeed excellently described by the effective particleapproach [28, 29, 33, 36, 67].Keber et al. have taken the first steps and developeda minimal model of defects as effective particles movingon a spherical surface [27]. Following their lead, variousimportant studies of the dynamics of active defects incurved surfaces have been proposed. In a recent publi-cation Khoromskaia and Alexander [68] refined the Ke-ber model by calculating the active flow via the Stokesequation of the active nematic. Although their theoreti-cal model explained some of the experimentally observedphenomena, it does not yet provide a detailed connectionbetween the coarse-grained dynamics of defects and theactive nematic hydrodynamic equations. Furthermore,nematic defects have a broken rotational symmetry, andhence their effective description must include not justtheir position as the sole degree of freedom, but also avectorial orientation [27, 69, 70]. This is especially ob-vious for the +1 / defects, which self-propel in a well-defined direction. Moreover, broken rotational symmetryimplies that these oriented particles are not just subjectto effective isotropic forces; rather, forces will dependon mutual orientation, and there will also be effective torques —all of which will affect the resulting dynamics[27, 35, 69–71]. Khoromskaia and Alexander have indeedexplored the impact of this orientational degree of free-dom on the active behavior of these systems [68], butthe resulting dynamics, and especially the role of effec-tive elastic torques, is still not understood well. In arecent paper [72], Brown has accounted for the orienta-tion dynamics of defects on the sphere, and partly ex-plained the experimental and simulated trajectories ofdefects. More generally, there is significant interest tounderstand how confinement impacts active nematody-namics in cases other than the surface of a sphere, es-pecially for different topologies and nontrivial Gaussiancurvature, such as a torus, and some important progressalong those lines has been made [73–75].In the present paper, we aim to develop a more de-tailed effective field theory for defects in active nematicsthat are confined to the surface of a sphere, based onactive nematic hydrodynamic. The main philosophy isno different from the planar case, and hence we expectthe essential scale separation to work just as well. How-ever, there is a very significant difference that rendersthis case more difficult, but also more interesting: theEuler characteristic of a spherical surface is +2 , and so we have broken “matter-antimatter symmetry”. In otherwords, the vacuum is not empty . Moreover, since it iswell known that in an active nematic the +1 / defectsare self propelled [33], this also implies that the vacuumis never at rest . Hence, particle pairs are not created intoan otherwise quiescent vacuum but are immediately ad-vected with the preexisting ground-state flow—which wehence need to understand first. We follow the ideas pro-posed by Tang and Selinger [76] to advance the theory ofa dry active nematic liquid crystal confined onto a spher-ical substrate. Specifically, we show how to make use ofan elegant variational principle due to Onsager [76–81]in order to transition from the infinite-dimensional ne-matic field theory to the finite-dimensional effective fieldtheory of oriented defects. As a particular result, we de-rive the anisotropic mobility coefficient matrix of defects.We then show that our effective theory fully reproducesthe complex periodic trajectory of the four +1 / defects,and how it depends on system size and the strength ofthe active forcing, as reported in earlier numerical workand experiments [27, 50, 51]. Our theory includes the fulldynamics of a defect’s orientation, which permits us toidentify its importance for the resulting particle motions.We have organized the content of our paper as fol-lows. In Sec. II, we outline the hydrodynamic modelwe employ in our discussion, followed by the theoreti-cal approach we use to separate the dynamical variablesof the defects from the temporal evolution of the activenematic field. In Sec. III, we look at some of the pre-dictions of our theory, including the well-known periodicoscillation of the defect configuration under weak forc-ing, and the chaotic trajectories under stronger drivingconditions. Sec. IV summarizes our conclusions and listsa number of limitations of our theory. II. MODELS AND METHODSA. Minimal model
The nematic liquid crystal is locally characterized bytwo variables: order and flow. We can describe the orderby the nematic tensor Q ab = q (cid:18) n a n b − g ab (cid:19) , (1)in which q is the magnitude of the nematic tensor repre-senting its average alignment in a small region, and g ab is the (inverse) metric tensor; the unit vector n = n a e a denotes the local nematic director, whose components n a = (cos ψ, sin ψ ) refer to the spherical orthonormal ba-sis { e θ , e φ } [3] and are fully specified by a single scalarfunction ψ ( θ, φ ) of position. (As usual, a repeated upperand lower index implies a tensor contraction.) The (in-compressible) flow is characterized by the velocity field V a . These two variables { Q ab , V a } satisfy the hydrody-namic equations of the active nematic, with a constantdensity ρ on the curved surface [41, 75, 82–85]: ˚ Q ab = S ab − γ δF LdG δQ ab , (2a) ρ ˚ V a = α ∇ b Q ba + ∇ b Π ab − ζV a . (2b)Here, the ring-operator (˚) ≡ D / D t ≡ ∂/∂ t + V c ∇ c de-notes the covariant material derivative , γ is the rotationalviscosity, and F LdG is the Landau-de Gennes free energyof the system [3], which contains the homogeneous phaseenergy F p , as has also recently been derived by dimen-sional reduction from a full three-dimensional Landau-de Gennes functional [86]: F p = (cid:90) d S (cid:18) A Q ab Q ab + B (cid:0) Q ab Q ab (cid:1) (cid:19) ; (3)and, by assuming (for simplicity) the one-constant-ap-proximation, the Frank energy F e caused by spatial dis-tortions: F e = K (cid:90) d S ( ∇ c Q ab )( ∇ c Q ab ) . (4)The tensor S ab in Eq. (2a) represents the couplingbetween the director field, the (symmetric) strain ratetensor A ab = ∇ [ a V b ] and the (antisymmetric) vorticity ω ab = ∇ ( a V b ) , with the parameter ξ reflecting the flowaligning of the nematic field [83, 87–91]: S ab = ξqA ab + Q ac ω cb − ω ac Q cb . (5)The activity enters into the system through the activestress α ∇ b Q ba in the covariant Navier-Stokes equation(2b), where the parameter α controls the strength of theactivity. The system is termed “contractile” if α > , and“extensile” if α < [43, 91, 92]. The passive stress is de-noted by Π ab . The substrate friction ζV a arises from thedamping force between the active nematic and the spher-ical substrate beneath it. If we neglect the inertia termin Eq. (2b) and assume the substrate friction is muchlarger than the other remaining dissipation terms, theonly contribution to the velocity field will be the activeflow [35, 71]: V a = αζ ∇ b Q ba . (6) B. Coarse grained dynamics of +1/2 defects
In our particular (two-dimensional) case, the index orcharge of a vector field’s defect may be defined as follows:take a closed path around an isolated defect and followthe orientation of the vector field as you traverse thatpath once, then the index is the number s of times whichthis vector field rotates before you arrive back at thestarting point (meaning that it acquires a “phase shift”of πs along this loop). We give a schematic of a +1 / defect in Fig. 1. It is clear to see that if we start from FIG. 1. Vector field with a +1 / defect at the center anda discontinuity line along the negative real axis due to the“phase shift” of π . any point at the negative real axis and trace the vectordirection around the original point clockwise, the vectoracquires a “phase shift” of π . The continuum model ofnematic liquid crystals on a spherical surface can be lo-cally described by a vector field on a spherical surface [3].Because the Euler characteristic of a spherical surface is +2 , the existence of defects is inevitable: summing theindices of all defects must yield +2 . The only differenceis that the apolar nature of a liquid crystal’s director field (meaning, it is already invariant under a ◦ in-plane ro-tation) permits half-integral indices for the defect points.In the present work we will restrict our discussion to thecase of four +1 / defects, which is the ground state ofthe nematic director field on the spherical surface.The motion of defects is a collective effect of the ne-matic field, so the dynamic information of the defects iscontained in the hydrodynamic equations we have justgiven. In order to extract the variables we care about,which are the position and the orientation of each de-fect, and obtain an effective field theory of the defects,we will now introduce the framework of Onsager’s varia-tional principle as an alternative description of the activenematic hydrodynamics [47, 77–81].The so-called Rayleighian corresponding to the Beris-Edwards equation (2a) is [93]: R = dd t F LdG + γ (cid:90) d S (cid:16) ˚ Q ab − S ab (cid:17) (cid:16) ˚ Q ab − S ab (cid:17) . (7)Its major use lies in the fact that the Beris-Edwards equa-tion follows from minimizing the Rayleighian (7) withrespect to ∂ t Q ab .Now, if we assume that the timescale of the active flow R tan θ θ θ FIG. 2. The stereographic mapping of the spherical surfaceonto the complex plane. is much slower than the characteristic timescale by whichthe director field n relaxes, it is reasonable to assumethat n stays close to its equilibrium configuration duringits active motion [33, 68]. By using Riemann’s stereo-graphic projection z ( θ, φ ) = 2 R tan θ e i φ , which maps thespherical coordinate ( θ, φ ) of the sphere with radius R tothe complex plane (such that the north pole becomes theorigin and the south pole maps to complex infinity; seeFig. 2), the equilibrium configuration of the director fieldwith four +1 / topological defects at specified positionswas given by [56, 68] ψ ( z ) = Ω − φ + 12 (cid:88) k =1 ψ k = Ω − φ + 12 (cid:88) k =1 Im ln( z − z k ) , (8)in which Im picks out the imaginary part of a complexnumber, z k = 2 R tan θ k e i φ k is the position of the k -thdefect, and ψ k represents the director field excited by asingle +1 / defect at z k . This is nothing but the solutionof the Euler-Lagrange equation of the Frank energy, effi-ciently expressed in complex notation. For the purposeof simplicity, we will take the equilibrium solution (8)as the configuration of the dynamic director field, thusomitting higher order corrections from the velocity of de-fects. This is in accord with our assumption of a smallactive flow.If we look at the director field near the j -th defect, anduse Ψ j to denote the angle between the symmetry axisof the defect and e θ , then Ψ j = 2Ω − φ j + (cid:88) k (cid:54) = j arg ( z j − z k ) . (9) According to Refs. [69, 70], this expression of Ψ j repre-sents the optimal relative orientation of defects. It im-plies an extra constraint for the time evolution of the ori-entation of defects. We will illustrate later that this con-straint is consistent with our assumption of the timescaleseparation. This also shows that Ω describes a global ro-tation of all the defects; we will demonstrate in the nextsection that it has a very important effect on their dy-namics.As a consequence of the small active flow assumption,the scalar order q varies only weakly outside the defects.We hence assume that q is constant away from the defectbut melts to zero inside a small core region, because thedirector would have to assume every orientation at thecenter. We will not try to describe the precise way inwhich the order vanishes towards the center, and it willindeed not be important. Instead, we will assume thatany integral area we are concerned about must excludethe defect core, and that outside the defects q is constant(we will assume q = 0 . in our subsequent numericalexamples to permit easier comparison with Ref. [50]).And since q is now constant over the entire domain ofinterest, the phase energy F p does not contribute to thedynamics. We can therefore restrict the free energy F LdG to the Frank energy contribution F e .With this ansatz, the time evolution of the directorfield is identical to the evolution of the position and thedirection of defects, which means ∂ t ψ = ˙Ω + ˙ θ k ∂ψ k ∂θ k + ˙ φ k ∂ψ k ∂φ k = ˙Ω + 12 (cid:88) k =1 (cid:16) i z − z k ˙ z k + c . c . (cid:17) (10)and the dot ( ˙ ) ≡ d / d t means a total derivative withrespect to time. Now, according to Onsager’s variationalprinciple, the minimization of the Rayleighian (7) withrespect to the dynamic variables of concern (here: veloc-ities of positions and global orientations of defects), willdirectly give us the dynamics of the defects. Hence, theirequations of motion can be succinctly written as ∂∂ ˙ θ k R = 0 , (11a) ∂∂ ˙ φ k R = 0 , (11b) ∂∂ ˙Ω R = 0 . (11c)By doing this, the infinite numbers of degrees of freedomin the nematic field are now reduced to nine degree offreedom for the four defects, which is our main aim inthis work.Now we can rewrite the Rayleighian (7) as R = dd t F e + γq (cid:90) d S ( ∂ t ψ ) − γαq ζ (cid:90) d S ∇ b ( ∇ a ψ − A a ) Q ab ∂ t ψ , (12)in which A a = e θ · ∂ a e φ is the spin connection. Detailsof the derivation are presented in Appendix A.The Frank energy, in turn, is given by [56, 58, 59, 64] F e = − πKq (cid:88) j (cid:54) = k ln (1 − cos β jk ) + const . , (13)where β jk is the angular distance between defects j and k , whose cosine can be expressed via cos β jk = cos θ j cos θ k + sin θ j sin θ k cos( φ j − φ k ) . (14)After making use of these, we can now write the dynamicequations (11a–11c) as M jk ˙ θ k + N jk ˙ φ k + Θ j ˙Ω − T j = − q γ ∂ θ j F e , (15a) N kj ˙ θ k + P jk ˙ φ k + Φ j ˙Ω − S j = − q γ ∂ φ j F e , (15b) Θ k ˙ θ k + Φ k ˙ φ k + 4 πR ˙Ω − L = 0 , (15c)where we introduced the following abbreviations for theeight different integrals that emerge in the process: M jk = (cid:90) d S ∂ θ j ψ j ∂ θ k ψ k , (16a) N jk = (cid:90) d S ∂ θ j ψ j ∂ φ k ψ k , (16b) P jk = (cid:90) d S ∂ φ j ψ j ∂ φ k ψ k , (16c) Θ j = (cid:90) d S ∂ θ j ψ j , (16d) Φ j = (cid:90) d S ∂ φ j ψ j , (16e) L = αζ (cid:90) d S ∇ b ( ∇ a ψ − A a ) Q ab , (16f) T j = αζ (cid:90) d S ∇ b ( ∇ a ψ − A a ) Q ab ∂ θ j ψ j , (16g) S j = αζ (cid:90) d S ∇ b ( ∇ a ψ − A a ) Q ab ∂ φ j ψ j . (16h)Notice that in all cases the integral area excludes a defectcore of size r k = r (1 + | z k | / R ) , which is the projec-tion image of the k -th defect core radius r at the com-plex plane [56]. Complete analytical expressions for allof these integrals, and some of the technically tedious de-tails for how to obtain them, are given in Appendices B,C, and D.If we introduce the characteristic timescale τ = γR / K and define the corresponding dimensionless time ˜ t = t/τ , then the dimensionless dynamic equations of thedefects of an active nematic on a sphere can be succinctlywritten as m jk d θ k d˜ t + n jk d φ k d˜ t − t j = − ∂ θ j f , (17a) n kj d θ k d˜ t + p jk d φ k d˜ t − s j = − ∂ φ j f , (17b) dΩd˜ t = − Φ k πR d φ k d˜ t , (17c)with the corresponding scaled variables m jk = M jk R , n jk = N jk R , p jk = P jk R − Φ j Φ k πR ,t j = γ K T j , s j = γ K S j , f = F e q K .
III. RESULTS
We have established the dynamics of defects in thespherical active nematic and expressed the mobility co-efficient matrices (16a–16c) in terms of active nematichydrodynamics parameters. Now we would like to com-pare our model with previous phenomenological theoriesand numerical studies.The analytical expressions (D7) and (D8) for the inte-grals T j and S j , worked out in Appendix D, suggest thatthese are actually the components of the mean activeback flow around the j -th defect—in agreement with thecalculations given by Khoromskaia and Alexander [68].We hence notice that the diagonal elements of the mo-bility matrices in Eq. (17a) and (17b) are m jj = π (cid:16) Rr − (cid:17) , (18a) p jj = π (cid:16) Rr − (cid:17) sin θ j . (18b)If the radius R of the sphere is very large compared tothe core radius r of the defects, and defects do not ap-proach each other too closely during their motion, thenondiagonal elements are negligible. In this situation,Eq. (17a) and (17b) are identical to the approximationsproposed in Ref. [68], and therefore the physical pictureis similar here. The motion of defects reflects the com-petition between the elastic and the active stresses: thevelocity of defects arises from the reorientation of the di-rector field due to the elastic stresses, and the advectionof the director field is due to the active back flow.Apart from the defect positions , we have also obtainedthe dynamics of their orientation in Eq. (17c). Thedefinition of the defect orientation and its passive dy-namics have been discussed thoroughly and adequatelyin Refs. [69, 70]. These authors discovered that the di-rector field can have an extra distortion because of therelative rotation of defects. Accordingly, our form of thedirector field in Eq. (8) is chosen such that the relativeorientation of the defects is set to their optimal configu-ration, which means it does not have any elastic interac-tion that arises from extra distortions of the director fielddue to a relative rotation of defects. In fact, accordingto Refs. [69, 70], Eq. (8) implies a special choice of theboundary condition of the director field near defects, sothat the relative orientation of defects is locked to a cer-tain angle. The same is implied in Eq. (17c), which doesnot contain the Frank elastic torque. We believe this isa reasonable result in the light of our assumption thatthe characteristic timescale of director field relaxation isfaster than that of defect motion. The activity enters thesystem via the back flow, according to Eq. (6), and thedirection of the flow is always along the symmetry axisof the defect; hence, the flow does not impact the defect transversely , and the elastic locking of defects maintainsthe relative orientation. As a consequence, the relaxationdynamics of the relative orientation of defects is purelypassive and very fast, and so we expect it to only enterat higher order in perturbation theory.In order to illustrate the predictions of our theory, wewill now present numerical results of the dynamic equa-tions (17a)–(17c). Here we chose R = 32 µ m , r = 2 µ m , γ/ K = 0 .
013 ms / µ m and the corresponding character-istic time scale τ b = γR / K = 13 .
312 ms as our refer-ence baseline. The choice of the initial condition of thesystem we discuss is θ i = π , φ i = π × ( i − , Ω = − . , (19)for i ∈ { , , , } . This describes four defects evenlydistributed in the equatorial plane, with their direc-tion deviating by . ◦ ( i = 1 , or . ◦ ( i = 2 , from the north-south direction. We used the Pythonmodule scipy.linalg to numerically solve the dynamicequations (17a–17c) via the fourth-order Runge-Kuttamethod, using a discretization timestep of ∆ t = 0 . τ b .Due to the huge reduction in the number of degrees offreedom (we now deal with nine first-order coupled differ-ential equations, instead of a partial differential equationin dimensions), these calculations are vastly moreefficient than if we had actually solved for the entire ne-matic field. For instance it took us only about
10 min on a laptop with Intel Core i7-6820HQ CPU at .
70 GHz to complete the τ b = 20 000 ∆ t trajectory shown inFig. 3(b).Our numerical results show that if the activity is be-low a certain threshold, the active forces are unable tocontinually overcome the elastic forces. After a brief dy-namic transient, during which the defects re-position intoan approximately tetrahedral pattern ( (cid:104) β (cid:105) > ∼ . ◦ ), thesystem approaches a steady state, in which despite ongo-ing flow the location of defects remains stationary—seeFig. 3(a).Above this threshold, our theory gives an excellent de-scription of the periodic but highly nonlinear ratchet-liketime evolution of the average angular distance (cid:104) β jk (cid:105) jk be-tween defects, as shown in Fig. 3(b), whose characteristicpattern has been previously discovered in both numericaland experimental work [27, 50, 51]. The oscillation is be-tween . ◦ and ◦ , corresponding to a transition be-tween the tetrahedral and the planar configuration. Wealso find that the dynamics of the defect orientation hasa significant influence on their trajectory. Recall that in the early phenomenological theory by Khoromskaia andAlexander [68] the global orientation of defects, Ω , wasfixed; but the authors understood the limitation of thisassumption and cautioned that Ω could be important incodetermining the motion of defects. As a result of thislimitation, their theory predicted that defects approacheach other in pairs. However, earlier simulation work byZhang et al. [50] had already shown that defects do notmove in pairs around their center of mass (in the regimeof low activity), but follow a complex trajectory similarto what had been reported in even earlier experimentsby Keber et al. [27]. We now find the same “zig-zag”type of complex motion in our theory, as can be seen inFig. 4. An illustration for the time evolution of such atrajectory can also be seen in the Supplemental Movie[94]. We believe that the dynamics of the defect orien-tation Ω , described in Eq. (17c), is responsible for thisqualitative change. Indeed, if we artificially fix the valueof Ω , instead of having it evolve via Eq. (17c), we observethe qualitatively different defect trajectories described byKhoromskaia and Alexander in Ref. [68].In both experiment [27] and the numerical work [50]it was found that the activity tunes the frequency of thedefects’ oscillations between the tetrahedral and the pla-nar configuration. We have reproduced this phenomenonin our theory. The numerical result of the dynamic equa-tions (17a)–(17c) reveals that the frequency linearly in-creases with increasing activity, as shown in Fig. 5(a).More remarkably, the size of the spherical substrate alsoaffects the dynamics of defects. If we fix the activity ofthe system at α/ζ = − . µ m / ms , but successively in-crease the sphere radius R , the frequency of the defects’oscillatory trajectory initially increases, until it attains amaximum, beyond which it again declines, as shown inFig. 5(b). To understand this, we notice that the coef-ficients of the mean active back flow components, (D7)and (D8), are both proportional to R/r , so increasing R is equivalent to linearly enhancing a defects’ veloc-ity. Meanwhile, the definition of characteristic timescale τ = γR / K suggests that a bigger R will quadratically slow the frequency of the defects’ oscillation in real time,because defects have larger and larger distances to tra-verse on their periodic trajectories as R increases. Even-tually, the competition of these two effects leads to aspecial substrate radius at which the frequency of thedefects’ oscillation attains its maximum, as shown inFig. 5(b).The coupling of the defect orientation and the activemotion causes a rather diverse set of patterns of motion.We observed the trajectory of defects becoming chaoticwhen the activity or the radius of the sphere substrateis in a certain range. For low activity, the trajectory re-mains closed and the dynamics displays a periodic tem-poral evolution, as shown in Figs. 6(a) and 6(d). Butwhen we increase the activity even further, our simula-tions show that the dynamics will gradually enter a newregime in which the trajectories no longer close; in fact,the numerical evidence strongly suggests that the system time (units of b ) d i s t a n ce ( d e g ree ) (a) < > time (units of b ) d i s t a n ce ( d e g ree ) (b) < > FIG. 3. Time evolution of the average angular distance (cid:104) β jk (cid:105) jk between defects at relatively low activity; recall that the planarconfiguration corresponds to (cid:104) β jk (cid:105) jk = 120 ◦ , while the tetrahedral configuration corresponds to (cid:104) β jk (cid:105) jk = 109 . ◦ . (a) For thesub-threshold activity α/ζ = − . µ m / ms , the defects approach stationary positions after a brief transient, reflecting thebalance between the elastic interaction of defects and the active flow. (b) For the activity α/ζ = − . µ m / ms , the defectsperiodically oscillate between the tetrahedral and the planar configuration. (a) (b) (c) FIG. 4. Zigzag defect trajectories for
R/r = 16 and α/ζ = − . µ m / ms . The four defects trajectories in (a) decomposeinto two pairs of defects trajectories in (b) and (c). becomes ergodic, based on the trajectory of defects inboth real space and phase space [see Figs. 6(b) and 6(e)],supporting the notion that at sufficiently high activitythe dynamics of defects will be chaotic. And just like thechaotic transition when increasing the activity, the defecttrajectory can be open for large sphere radius, as shownin Figs. 6(c) and 6(f). The chaotic transition we found in our theory is con-sistent with what Zhang et al. [50] observed in theirsignificantly more challenging (Lattice-Boltzmann based)simulations for the evolution of the entire nematic field—which is potentially more prone to numerical instabilities.It is hence reassuring to see that the same transition tochaos is observed in our much simpler system of coupledordinary differential equations. IV. CONCLUSIONS AND OUTLOOK
Using Onsager’s variational approach to nonequilib-rium (thermo-)dynamics, we have derived an effectivetheory for the defect dynamics in an active nematic con-fined to the surface of a sphere. Within our formalism thestate of the system is fully determined by the positionsand orientations of a small number of defects, whose dy-namics follows from solving a set of coupled first-orderdifferential equations. This is both conceptually moredirect and numerically more efficient than the alterna- / ( m /ms ) / T ( un i t s o f b ) (a)
40 50 60 70
R ( m ) / T ( un i t s o f b ) (b) FIG. 5. (a) Oscillation frequency of periodic defect trajec-tories [such as those in Figs. 6(a) and 6(d)] as a functionof scaled activity − α/ζ for given sphere size R = 32 µ m .(b) Same oscillation frequency, but now as a function ofthe substrate sphere radius R for a given scaled activity of α/ζ = − . µ m / ms . tive of working with a partial differential equation forthe nematic field, which lives in (2 + 1) dimensions anddescribes either the full order parameter tensor or at leastthe nematic director.We focused on the fundamental question how to ar-rive at such an effective description, and hence we havemade a few approximations along the way in order toavoid cluttering the path with incidental technical diffi-culties. Our theory is therefore no equivalent rewritingof the original problem but rather a first-order approx-imation for the exact time evolution in terms of conve-nient smallness parameters. Nevertheless, it offers severalquantitative and qualitative insights, and it efficientlyreproduces previous numerical and experimental results.For instance, in the weakly driven regime we confirmedthe existence of a threshold activity below which defectmotion is arrested, and we characterized the shape of thedefect trajectories past this threshold. In the high ac-tivity regime we observed transition into chaos. And asan example of a conceptual insight, we have elucidatedthe importance of promoting the global rotation Ω to adynamic variable in order to reproduce the type of tra-jectories observed in experiment.Our theoretical development can be generalized in anumber of ways. Most straightforwardly, one could relaxthe one-constant-approximation and investigate how dif-fering moduli for splay and bend affect configuration anddynamics of an active nematic field with nontrivial topol-ogy [49]. More challenging, in order to understand thesometimes surprisingly large fluctuations of the nematicdensity observed in both experiments [27] and numeri-cal studies [51], one must include density contributionsto the Rayleighian (7) and add the expression that fol-lows from the associated variation to the set of dynamicequations.More interesting, but also much more challenging,would be steps to complete our formalism into a truefluctuating effective field theory, in which particle num-ber is not conserved; in other words, to incorporate pair creation and annihilation events. These occur at suffi-ciently large activity or sphere radius and have been nu-merically observed [51]. The difficulty is that annihilat-ing a +1 / against a − / defect, or creating such a pairout of an initially smooth background field, ultimatelyrests on the local physics our effective theory has inte-grated out. It might seem that at least the annihilation event is relatively easy to account for, as one could simplyeliminate a particle-antiparticle pair as soon as they ap-proach closer than some critical distance; but formulatingconditions on the local smooth field that trigger pair cre-ation is nontrivial. Moreover, since the rates of creationand annihilation will surely satisfy joint thermodynamicconstraints, these two processes are not independent andhence cannot be treated in isolation.Even more challenging questions arise once we lookat actual experimental realizations of active nematics,such as the lipid vesicles covered by microtubule filamentsthat are rendered active by the addition of kinesin mo-tors [27]. The coupling between the fluid-elastic curva-ture energetics of the vesicle and the elasticity underlyingthe nematic liquid crystal impacts the morphology of thevesicle, which need not remain a perfect sphere. In fact,protrusions growing from the location of the defects havebeen observed [27], and the vesicle may deform into aspindle-like structure with two +1 defects at the spin-dle poles [27]. Metselaar et al. [95] have recently pro-posed a continuum model for such a deformable activenematic shell. Our own framework could be extendedto include this as well, by incorporating Helfrich’s cur-vature elastic energy [96] for the vesicle shape and thedissipation associated with changes of this shape into theRayleighian (7). However, this will require some signif-icant new formalism, for instance replacing the simpleRiemann mapping of a sphere to a more general surfaceparametrization. But since throughout our developmentthe description of the geometry is fully covariant, we sus-pect that the framework of Onsager’s variational princi-ple will continue to provide a convenient starting pointfor developing an effective theory. ACKNOWLEDGMENTS
We are grateful for valuable and inspiring discussionswith Masao Doi and Tanniemola Liverpool. We also ac-knowledge financial support from the National NaturalScience Foundation of China (Grants No. 11675017 andNo. 11975050) and partial support from the National Sci-ence Foundation of the United States (Grant No. CHE1764257).
Appendix A: The Rayleighian
According to Onsager’s variational principle, onlyterms that explicitly contain ∂ t Q ab will contribute to thedynamic equations of defects. So the Rayleighian can be (a) x v x (d) (b) x v x (e) (c) x v x (f) FIG. 6. (a–c) Trajectory of an individual defect in three dimensions. (d–f) Projection of the phase space trajectory of the sameindividual defect into the X - V X plane. Specifically: (a, d) Defect trajectories for small sphere radius R/r = 16 and low activity γα/ Kζ = − . . (b, e) Defect trajectories for small sphere radius R/r = 16 and high activity γα/ Kζ = − . . (c, f)Defect trajectories for large sphere radius R/r = 75 and low activity γα/ Kζ = − . . written as R = dd t F + γ (cid:90) d S (cid:0) ∂ t Q ab ∂ t Q ab +2 ∂ t Q ab V c ∇ c Q ab − ∂ t Q ab S ab (cid:1) . (A1)If we substitute the active flow V a = αζ ∇ b Q ba , we obtain ∂ t Q ab S ab = ξq∂ t Q ab ∇ [ a V b ] + ∂ t Q ab Q ac ∇ ( c V b ) − ∂ t Q ab ∇ ( a V c ) Q cb = αξqζ ∂ t Q ab ∇ a ∇ d Q db + αζ ∂ t Q ab Q ac ∇ ( c ∇ d Q db ) − αζ ∂ t Q ab ∇ ( a ∇ d Q dc ) Q cb = αζ ∂ t Q ab Q ac ∇ ( c ∇ d Q db ) − αζ ∂ t Q ab ∇ ( a ∇ d Q dc ) Q cb + αξq ζ ∆ ψ∂ t ψ (A2)The equilibrium configuration of the director field ensuresthat the last term in Eqn. (A2) vanishes, so the flow aligning parameter ξ does not contribute. Finally, all thecoupling terms read γ (cid:90) d S (cid:0) ∂ t Q ab ∂ t Q ab + 2 ∂ t Q ab V c ∇ c Q ab − ∂ t Q ab S ab (cid:1) = γq (cid:90) d S ( ∂ t ψ ) − γαq ζ (cid:90) d S ∇ b ( ∇ a ψ − A a ) Q ab ∂ t ψ. (A3) Appendix B: Riemann sphere representation
Before calculating integrals (16a)-(16h), we must ex-press each term in the Riemann sphere representation byusing the stereographic projection z ( θ, φ ) = 2 R tan θ e i φ shown in the main text. This leads to (cid:90) d S = (cid:90) d z d¯ z (cid:0) z ¯ z R (cid:1) (B1) ∂ θ k ψ k = i4 R (cid:16) z k ¯ z k R (cid:17) (cid:114) z k ¯ z k z − z k + c . c ., (B2)0 ∂ φ k ψ k = − z k z − z k + c . c ., (B3) ∂ φ ψ k = i4 zz − z k − i4 z ( z − z k ) + c . c ., (B4) ∂ θ ∂ φ ψ k = R (cid:16) z ¯ z R (cid:17) (cid:114) z ¯ z z − z k − R (cid:16) z ¯ z R (cid:17) (cid:114) z ¯ z z ( z − z k ) + c . c ., (B5) ∂ θ ψ k = − i8 R (cid:16) z ¯ z R (cid:17) R √ z ¯ z R + z ¯ z (cid:114) z ¯ z z − z k + i4 R (cid:16) z ¯ z R (cid:17) z ¯ z z − z k ) + c . c . . (B6)We can calculate the area integrals (16a-16h) based onthe following theorem:For a given function G in an area S with the boundary ∂S in the complex plane z = x + i y , (cid:90) S ∂ z G ( z, ¯ z ) d z d¯ z = (cid:90) S ∂ z G ( z, ¯ z ) d x d y = (cid:90) S ( ∂ x − i ∂ y ) G d x d y . (B7)Based on Green’s theorem, if the first derivative of G iscontinuous in S , then we have (cid:90) S ∂ z G ( z, ¯ z ) d z d¯ z = (cid:90) S ( ∂ x − i ∂ y ) G d x d y = (cid:73) ∂S ( G d y + i G d x )= i (cid:73) ∂S G d¯ z. (B8)Then we can use (cid:90) S ∂ z G ( z, ¯ z ) d z d¯ z = i (cid:73) ∂S G d¯ z (B9)to transform an area integral into a boundary integral.And if G is not analytic, the integral will depend on theshape of the contour of the boundary, instead of just thetopology. Appendix C: Mobility matrices
In order to obtain the result of integral in (16a-16e),there are two basic integrals we need to calculate first: (cid:90) d S z − z j ) ( z − z k ) , (C1) (cid:90) d S z − ¯ z j ) ( z − z k ) . (C2)According to Eq. (B1)-(B3), the mobility matrices (16a)-(16e) can be considered as linear combinations of theintegrals (C1) and (C2) with their complex conjugates. The boundaries ∂S are those peripheries around each de-fects based on the integral area we require. First, we useEq. (B9) to transform the area integrals into line inte-grals at each boundaries. Then we expand the integrandnear each boundary, and only keep the first order basedon the assumption R (cid:29) r . The detailed procedure isdemonstrated below.
1. The calculation of the integral (C1) a. Diagonal element: k = j (cid:90) d z d¯ z (cid:0) z ¯ z R (cid:1) ( z − z j ) =i (cid:73) C d z R z (4 R + z ¯ z ) ( z − z j ) =i (cid:73) C + C j d z − R z (4 R + z ¯ z ) ( z − z j ) =i (cid:73) C d z − R z ( z − z j ) + i (cid:73) C j d z − R z (4 R + z ¯ z j ) ( z − z j ) = 4 πR ¯ z j (4 R + z j ¯ z j ) (C3)In order to fulfill the condition required by Eq. (B9), weneed to exclude the origin of the complex plane so thatthe integrand is continuous over the domain of integra-tion. Then the boundary C should contain the boundary C around the origin, and the boundary C j around thedefect at z j . And with performing the same procedure,the nondiagonal element can be calculated as follows. b. Nondiagonal element: k (cid:54) = j (cid:90) d z d¯ z (cid:0) z ¯ z R (cid:1) ( z − z j ) ( z − z k )=i (cid:73) C d z R z (4 R + z ¯ z ) ( z − z j ) ( z − z k )=i (cid:73) C + C j + C k d z − R z (4 R + z ¯ z ) ( z − z j ) ( z − z k )=i (cid:73) C d z − R z ( z − z j ) ( z − z k )+ i (cid:73) C j d z − R z (4 R + z ¯ z j ) ( z − z j ) ( z − z k )+ i (cid:73) C k d z − R z (4 R + z ¯ z k ) ( z − z j ) ( z − z k )= 16 πR (¯ z k − ¯ z j ) − πR ¯ z j ¯ z k ( z k − z j )( z j − z k ) (4 R + z j ¯ z j ) (4 R + z k ¯ z k ) (C4)1
2. The calculation of the integral (C2) a. Diagonal element: k = j (cid:90) d z d¯ z (cid:0) z ¯ z R (cid:1) (¯ z − ¯ z j ) ( z − z j )= − i (cid:73) C d z (cid:20) R (4 R + z ¯ z ) (4 R + z ¯ z j ) ( z − z j )+ 8 R (4 R + z ¯ z j ) ( z − z j ) ln z ¯ z − z ¯ z j R + z ¯ z (cid:35) = − i (cid:73) C d z (cid:20) R (4 R + z ¯ z ) (4 R + z ¯ z j ) ( z − z j )+ 8 R (4 R + z ¯ z j ) ( z − z j ) ln ( z − z j ) (¯ z − ¯ z j )4 R + z ¯ z (cid:35) − i (cid:73) C d z R (4 R + z ¯ z j ) ( z − z j ) ln zz − z j = − i (cid:73) C d z (cid:20) R (4 R + z ¯ z ) (4 R + z ¯ z j ) ( z − z j )+ 8 R (4 R + z ¯ z j ) ( z − z j ) ln ( z − z j ) (¯ z − ¯ z j )4 R + z ¯ z (cid:35) =i (cid:73) C d z (cid:34) R z ¯ z j (4 R + z ¯ z j ) ( z − z j )+ 8 R (4 R + z ¯ z j ) ( z − z j ) (cid:18) ln z j ¯ z j − z ¯ z j R + 1 (cid:19)(cid:35) + i (cid:73) C j d z R (4 R + z ¯ z j ) ( z − z j ) (cid:32) ln r j R + z ¯ z j + 1 (cid:33) + i (cid:73) C ∗ j d z (cid:34) − R z j (4 R + z ¯ z j ) ( z − z j ) + 8 R (4 R + z ¯ z j ) ( z − z j ) ln 4 R + z j ¯ z j R (cid:35) = 4 πR z j ¯ z j (4 R + z j ¯ z j ) − πR (4 R + z j ¯ z j ) (cid:18) Rr j R + z j ¯ z j (cid:19) (C5)The integrand indicates that there are three singularpoints: z = 0 , z = z j , and z = − R / ¯ z j , correspondingto boundaries of C , C j and C ∗ j separately. In addition,the first two are also branch points because of the log-arithmic function. In order to use Eq. (B9), we againneed the integrand to be continuous along the boundary.We therefore add a branch cut at both of the boundary C around the branch point z = 0 , and the boundary C j around the branch point z = z j . Once we take theboundary C = C + C j + C ∗ j as shown in Fig.(7), the lineintegral of the sixth line is zero, according to the residuetheorem. i . z j − 4 R ¯ z j FIG. 7. Integral contour of (C5). b. Non-diagonal element: k (cid:54) = j (cid:90) d z d¯ z (cid:0) z ¯ z R (cid:1) (¯ z − ¯ z j ) ( z − z k )= − i (cid:73) C d z (cid:20) R (4 R + z ¯ z ) (4 R + z ¯ z j ) ( z − z k )+ 8 R (4 R + z ¯ z j ) ( z − z k ) ln z ¯ z − z ¯ z j R + z ¯ z (cid:35) =i (cid:73) C d z (cid:34) R z ¯ z j (4 R + z ¯ z j ) ( z − z k )+ 8 R (4 R + z ¯ z j ) ( z − z k ) (cid:18) ln z j ¯ z j − z ¯ z j R + 1 (cid:19)(cid:35) + i (cid:73) C j d z R (4 R + z ¯ z j ) ( z − z k ) (cid:32) ln r j R + z ¯ z j + 1 (cid:33) + i (cid:73) C k d z (cid:20) R (4 R + z ¯ z j ) (4 R + z ¯ z k ) ( z − z k )+ 8 R (4 R + z ¯ z j ) ( z − z k ) ln ( z − z j ) (¯ z k − ¯ z j )4 R + z ¯ z k (cid:35) + i (cid:73) C ∗ j d z (cid:20) − R z j (4 R + z ¯ z j ) ( z − z j ) ( z − z k )+ 8 R (4 R + z ¯ z j ) ( z − z k ) ln 4 R + z j ¯ z j R (cid:35) = − πR (4 R + z k ¯ z j ) (cid:34) (cid:0) R + z k ¯ z j (cid:1) (cid:0) R − z j ¯ z j z k ¯ z k (cid:1) (4 R + z j ¯ z j ) (4 R + z k ¯ z k )+4 R ln 4 R ( z j − z k ) (¯ z j − ¯ z k )(4 R + z j ¯ z j ) (4 R + z k ¯ z k ) (cid:21) (C6)2Notice that the result seems to be divergent when z k = − R / ¯ z j . This is because we treat z k and − R / ¯ z j astwo different singular points when performing the residuetheorem in the calculation. If we instead consider thecase of z k = − R / ¯ z j , we have: (cid:90) d z d¯ z (cid:0) z ¯ z R (cid:1) (¯ z − ¯ z j ) ( z − z k )= (cid:90) d z d¯ z (cid:0) z ¯ z R (cid:1) (¯ z − ¯ z j ) (cid:16) z + R ¯ z j (cid:17) = − i (cid:73) C d z (cid:34) R ¯ z j (4 R + z ¯ z ) (4 R + z ¯ z j ) + 8 R ¯ z j (4 R + z ¯ z j ) ln z ¯ z − z ¯ z j R + z ¯ z (cid:35) =i (cid:73) C d z (cid:34) R z ¯ z j (4 R + z ¯ z j ) + 8 R ¯ z j (4 R + z ¯ z j ) (cid:18) ln z j ¯ z j − z ¯ z j R + 1 (cid:19)(cid:35) + i (cid:73) C j d z R ¯ z j (4 R + z ¯ z j ) (cid:32) ln r j R + z ¯ z j + 1 (cid:33) + i (cid:73) C ∗ j d z (cid:34) − R z j ¯ z j (4 R + z ¯ z j ) ( z − z j )+ 8 R ¯ z j (4 R + z ¯ z j ) ln 4 R + z j ¯ z j R (cid:35) = − πR z j ¯ z j (4 R + z j ¯ z j ) . (C7)We can check that the result in (C7) indeed coincideswith (C6) in the limit in which z k approaches − R / ¯ z j .In fact, the singular point − R / ¯ z j is the antipode of thedefect z j on the sphere, and the singularity does not haveany physical meaning.
3. Results of mobility matrices
Substituting the corresponding coefficients accordingto Eqs. (B2) and (B3), the elements of the mobility ma-trices can be constructed as follows:
Diagonal element: P jj = πR | z j | (4 R + | z j | ) − πR | z j | (4 R + | z j | ) (cid:18) Rr j R + | z j | (cid:19) , (C8) M jj = − πR (cid:18) Rr j R + | z j | (cid:19) . (C9) Nondiagonal element: P jk = − πR z k ¯ z j R + z k ¯ z j ) (cid:34) (cid:0) R + z k ¯ z j (cid:1) (cid:0) R − | z j | | z k | (cid:1) (4 R + | z j | ) (4 R + | z k | )+4 R ln 4 R | z j − z k | (4 R + | z j | ) (4 R + | z k | ) (cid:21) + πR z j z k z j − z k ) (cid:18) ¯ z k R + | z k | − ¯ z j R + | z j | (cid:19) + c . c . . (C10) M jk = 4 πR (¯ z j − ¯ z k ) + π ¯ z k ¯ z j ( z k − z j )64 ( z j − z k ) (cid:114) z j z k ¯ z j ¯ z k − π (cid:0) R − | z j | | z k | (cid:1)
64 (4 R + z k ¯ z j ) (cid:115) ¯ z j z k z j ¯ z k − πR (cid:0) R + | z j | (cid:1) (cid:0) R + | z k | (cid:1)
16 (4 R + z k ¯ z j ) (cid:115) ¯ z j z k z j ¯ z k ln 4 R | z j − z k | (4 R + | z j | ) (4 R + | z k | ) + c . c . , (C11) N jk = − i πRz k (cid:114) z j ¯ z j (cid:20) R (¯ z k − ¯ z j )4 ( z j − z k ) (4 R + | z k | )+ ¯ z j ¯ z k
16 (4 R + | z k | ) (cid:21) + i πR ¯ z k (cid:0) R − | z j | | z k | (cid:1)
16 (4 R + z j ¯ z k ) (4 R + | z k | ) (cid:114) z j ¯ z j + i πR ¯ z k (cid:0) R + | z j | (cid:1) R + z j ¯ z k ) (cid:114) z j ¯ z j ln 4 R | z j − z k | (4 R + | z j | ) (4 R + | z k | ) + c . c . . (C12)It is also easy to calculate the vector Φ j and Θ j : Φ j = (cid:90) d S∂ φ j ψ j = − (cid:90) d z d¯ z R z j (4 R + z ¯ z ) ( z − z j ) + c . c . = − i (cid:73) C d z R z j z (4 R + z ¯ z )( z − z j )= i (cid:73) C d z R z j z − zz j + i (cid:73) C j d z R z j z ( z − z j )(4 R + z ¯ z j ) + c . c . = 2 πR | z j | R + | z j | , (C13) Θ j = (cid:90) d S∂ θ j ψ j = (cid:90) d z d¯ z i R (cid:0) R + z j ¯ z j (cid:1) R + z ¯ z ) ( z − z j ) (cid:114) z j ¯ z j + c . c . = − (cid:73) C d z R (cid:0) R + z j ¯ z j (cid:1) z (4 R + z ¯ z )( z − z j ) (cid:114) z j ¯ z j = (cid:73) C d z R (cid:0) R + z j ¯ z j (cid:1) z ( z − z j ) (cid:114) z j ¯ z j + (cid:73) C j d z R (cid:0) R + z j ¯ z j (cid:1) z ( z − z j ) (4 R + z ¯ z j ) (cid:114) z j ¯ z j + c . c . = − i π R | z j | + i π R | z j | = 0 . (C14) Appendix D: The calculation of T j , R j , and L Here we calculate T j as an example to illustrate theprocedure: T j = − q α R ζ (cid:90) d S (cid:2) sin 2 ψ (cid:0) θ − θ csc θ∂ φ ψ − θ∂ θ ∂ φ ψ ) − cos 2 ψ (cid:0) ∂ θ ψ − csc θ∂ φ ψ − cot θ∂ θ ψ (cid:1)(cid:3) ∂ θ j ψ j . (D1)Under the Riemann sphere representation, we have sin 2 ψ = − i2 exp (2iΩ) ¯ zz (cid:89) p =1 (cid:115) z − z p ¯ z − ¯ z p + c . c . , (D2) cos 2 ψ = 12 exp (2iΩ) ¯ zz (cid:89) p =1 (cid:115) z − z p ¯ z − ¯ z p + c . c . . (D3)It seems very hard to find a global function G , so that T j = i (cid:72) d¯ z G . Our strategy is to expand Eq. (D1) neareach of the defects, and subsequently do the procedurewe performed above. According to Eqs. (B1)-(B6), wecan expand the integrand of T j in terms of z − z j in thefollowing form: T j = i (cid:73) C j d¯ z (cid:90) d z (cid:115) z − z j ¯ z − ¯ z j (cid:88) q =1 q (cid:88) m =0 (cid:34) z − z j ) m (¯ z − ¯ z j ) q − m ∞ (cid:88) l =0 A m,q,l ( z − z j ) l (cid:35) + i (cid:88) k (cid:54) = j (cid:73) C k d¯ z (cid:90) d z (cid:114) z − z k ¯ z − ¯ z k (cid:88) n =0 (cid:34)(cid:18) z − z k ) n + 1(¯ z − ¯ z k ) n (cid:19) ∞ (cid:88) l =0 B n,l ( z − z k ) l (cid:35) + c . c . = i (cid:73) C j d¯ z (cid:88) q =1 q (cid:88) m =0 ∞ (cid:88) l =0 A m,q,l r − m +2 lj − m + 2 l (¯ z − ¯ z j ) m − l − q − + i (cid:88) k (cid:54) = j (cid:73) C k d¯ z (cid:88) n =0 ∞ (cid:88) l =0 B n,l (cid:34) r l − nk l − n (¯ z − ¯ z k ) n − l − + 2 r lk l (¯ z − ¯ z k ) − n − l − (cid:35) + c . c . . (D4)We can also expand it in terms of ¯ z − ¯ z j : T j = − i (cid:73) C j d z (cid:90) d¯ z (cid:115) z − z j ¯ z − ¯ z j (cid:88) q =1 q (cid:88) m =0 (cid:34) z − z j ) m (¯ z − ¯ z j ) q − m ∞ (cid:88) l =0 C m,q,l (¯ z − ¯ z j ) l (cid:35) − i (cid:88) k (cid:54) = j (cid:73) C k d z (cid:90) d¯ z (cid:114) z − z k ¯ z − ¯ z k (cid:88) n =0 (cid:34)(cid:18) z − z k ) n + 1(¯ z − ¯ z k ) n (cid:19) ∞ (cid:88) l =0 D n,l (¯ z − ¯ z k ) l (cid:35) + c . c . = − i (cid:73) C j d z (cid:88) q =1 q (cid:88) m =0 ∞ (cid:88) l =0 C m,q,l r − q +2 m +2 lj − q + 2 m + 2 l ( z − z j ) q − m − l − i (cid:88) k (cid:54) = j (cid:73) C k d z (cid:88) n =0 ∞ (cid:88) l =0 D n,l (cid:34) r l +1 k l + 1 ( z − z k ) − n − l + 2 r l − nk l − n ( z − z k ) n − l (cid:35) + c . c . , (D5)which represents the order change of the area integral. A , B , C and D are the expansion coefficients, which can beobtained by expanding the integrand near the boundariesaround each defects. And we should notice that A m,q, = C m,q, , B n, = D n, , (D6)This relationship means Eqs. (D4) and (D5) should givethe same result, which reflects the order free of the in-tegral with respect to z and ¯ z . It is easy to check thatthe leading term that contributes to T j is − A , , /r j or − C , , /r j .Finally, the expression of T j is T j = Re πq α Rζr j exp (2iΩ) (cid:0) R + | z j | (cid:1) (cid:115) ¯ z j z j (cid:89) l (cid:54) = j (cid:115) z j − z l ¯ z j − ¯ z l = πq αR ζr cos ( w j + 2Ω − φ j ) . (D7)We can perform the same procedure to the calculation of S j and L : S j = Re − i πq α ¯ z j ζr j exp (2iΩ) (cid:89) l (cid:54) = j (cid:115) z j − z l ¯ z j − ¯ z l = πq αR ζr sin θ j sin ( w j + 2Ω − φ j ) , (D8) L =0 . (D9)And w j is the summation of the phase angle between the j -th defect and others: w j = (cid:88) l (cid:54) = j arg ( z j − z l ) . (D10)4 [1] S. Chandrasekhar and G. S. Ranganath, The structureand energetics of defects in liquid crystals, Adv. Phys. , 507 (1986).[2] M. V. Kurik and O. D. Lavrentovich, Defects in liq-uid crystals: homotopy theory and experimental studies,Sov. Phys. Usp. , 196 (1988).[3] P. G. de Gennes and J. Prost, The Physics of LiquidCrystals (Clarendon Press, Oxford, 1993).[4] H. Kikuchi, M. Yokota, Y. Hisakado, H. Yang, and T. Ka-jiyama, Polymer-stabilized liquid crystal blue phases, Na-ture Materials , 64 (2002).[5] G. I. Taylor, The mechanism of plastic deformation ofcrystals. part i.–theoretical, Proc. Royal Soc. Lond. Ser.A , 362 (1934).[6] J. Friedel, Dislocations (Pergamon Press, Oxford, 1964).[7] F. Nabarro,
Theory of Crystal Dislocations (ClarendonPress, Oxford, 1967).[8] G. Duclos, C. Erlenkämper, J.-F. Joanny, and P. Sil-berzan, Topological defects in confined populations ofspindle-shaped cells, Nature Physics , 58 (2016).[9] T. B. Saw, A. Doostmohammadi, V. Nier, L. Kocgozlu,S. Thampi, Y. Toyama, P. Marcq, C. T. Lim, J. M. Yeo-mans, and B. Ladoux, Topological defects in epitheliagovern cell death and extrusion, Nature , 212 (2017).[10] K. Kawaguchi, R. Kageyama, and M. Sano, Topologicaldefects control collective dynamics in neural progenitorcell cultures, Nature , 327 (2017).[11] D. V. Osborne, The rotation of liquid helium ii, Proc.Phys. Soc., Sec. A , 909 (1950).[12] J. Wilks and D. Betts, An Introduction to Liquid Helium (Clarendon Press, Oxford, 1987).[13] D. Vollhardt and P. Wolfle,
The Phases of Helium 3 (Tay-lor and Francis, New York, 1990).[14] N. Bogdanov and D. A. Yablonskii, Thermodynamicallystable “vortices” in magnetically ordered crystals. themixed state of magnets., Sov. Phys. JETP , 101 (1989).[15] A. Bogdanov and A. Hubert, Thermodynamically stablemagnetic vortex states in magnetic crystals, J. Magn.Magn. Mater. , 255 (1994).[16] U. K. Rößler, A. N. Bogdanov, and C. Pfleiderer, Sponta-neous skyrmion ground states in magnetic metals, Nature , 797 (2006).[17] B. Binz, A. Vishwanath, and V. Aji, Theory of the helicalspin crystal: A candidate for the partially ordered stateof mnsi, Phys. Rev. Lett. , 207202 (2006).[18] S. Tewari, D. Belitz, and T. R. Kirkpatrick, Blue quan-tum fog: Chiral condensation in quantum helimagnets,Phys. Rev. Lett. , 047207 (2006).[19] S. Mühlbauer, B. Binz, F. Jonietz, C. Pfleiderer,A. Rosch, A. Neubauer, R. Georgii, and P. Böni,Skyrmion lattice in a chiral magnet, Science , 915(2009).[20] A. Neubauer, C. Pfleiderer, B. Binz, A. Rosch, R. Ritz,P. G. Niklowitz, and P. Böni, Topological hall effect inthe a phase of mnsi, Phys. Rev. Lett. , 186602 (2009).[21] C. Pappas, E. Lelièvre-Berna, P. Falus, P. M. Bentley,E. Moskvin, S. Grigoriev, P. Fouquet, and B. Farago,Chiral paramagnetic skyrmion-like phase in mnsi, Phys.Rev. Lett. , 197202 (2009).[22] N. Nagaosa and Y. Tokura, Topological properties anddynamics of magnetic skyrmions, Nature Nanotech. , 899 (2013).[23] N. Turok, Global texture as the origin of cosmic struc-ture, Phys. Rev. Lett. , 2625 (1989).[24] A. Pargellis, N. Turok, and B. Yurke, Monopole-antimonopole annihilation in a nematic liquid crystal,Phys. Rev. Lett. , 1570 (1991).[25] I. Chuang, R. Durrer, N. Turok, and B. Yurke, Cosmol-ogy in the laboratory: Defect dynamics in liquid crystals,Science , 1336 (1991).[26] A. J. Bray, Theory of phase-ordering kinetics, Advancesin Physics , 357 (1994).[27] F. C. Keber, E. Loiseau, T. Sanchez, S. J. DeCamp,L. Giomi, M. J. Bowick, M. C. Marchetti, Z. Dogic, andA. R. Bausch, Topology and dynamics of active nematicvesicles, Science , 1135 (2014).[28] L. Giomi, M. J. Bowick, X. Ma, and M. C. Marchetti,Defect annihilation and proliferation in active nematics,Phys. Rev. Lett. , 228101 (2013).[29] L. M. Pismen, Dynamics of defects in an active nematiclayer, Phys. Rev. E , 050502 (2013).[30] S. P. Thampi, R. Golestanian, and J. M. Yeomans, Ve-locity correlations in an active nematic, Phys. Rev. Lett. , 118101 (2013).[31] S. P. Thampi, R. Golestanian, and J. M. Yeomans, In-stabilities and topological defects in active nematics, Eu-rophys. Lett. , 18001 (2014).[32] S. P. Thampi, R. Golestanian, and J. M. Yeomans, Vor-ticity, defects and correlations in active turbulence, Phi-los. Trans. R. Soc. A , 20130366 (2014).[33] L. Giomi, M. J. Bowick, P. Mishra, R. Sknepnek, andM. Cristina Marchetti, Defect dynamics in active nemat-ics, Philos. Trans. R. Soc. A , 20130365 (2014).[34] T. Gao, R. Blackwell, M. A. Glaser, M. D. Betterton, andM. J. Shelley, Multiscale polar theory of microtubule andmotor-protein assemblies, Phys. Rev. Lett. , 048101(2015).[35] S. Shankar, S. Ramaswamy, M. C. Marchetti, and M. J.Bowick, Defect unbinding in active nematics, Phys. Rev.Lett. , 108002 (2018).[36] D. Cortese, J. Eggers, and T. B. Liverpool, Pair creation,motion, and annihilation of topological defects in two-dimensional nematic liquid crystals, Phys. Rev. E ,022704 (2018).[37] J. L. Ericksen, Anisotropic fluids, Arch. Ration. Mech.Anal. , 231 (1959).[38] J. L. Ericksen, Conservation laws for liquid crystals,Trans. Soc. Rheol. , 23 (1961).[39] F. M. Leslie, Some constitutive equations for anisotropicfluids, Quart. J. Mech. Appl. Math. , 357 (1966).[40] F. M. Leslie, Some constitutive equations for liquid crys-tals, Arch. Ration. Mech. Anal. , 265 (1968).[41] A. Beris and B. Edwards, Thermodynamics of FlowingSystems (Oxford University Press, Oxford, 1994).[42] R. Aditi Simha and S. Ramaswamy, Hydrodynamic fluc-tuations and instabilities in ordered suspensions of self-propelled particles, Phys. Rev. Lett. , 058101 (2002).[43] M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B.Liverpool, J. Prost, M. Rao, and R. A. Simha, Hydrody-namics of soft active matter, Rev. Mod. Phys. , 1143(2013).[44] J. Prost, F. Jülicher, and J.-F. Joanny, Active gel physics, Nature Physics , 111 (2015).[45] S. Ramaswamy, Active matter, Journal of Statistical Me-chanics: Theory and Experiment , 054002 (2017).[46] A. Doostmohammadi, J. Ignés-Mullol, J. M. Yeomans,and F. Sagués, Active nematics, Nature Comm. , 3246(2018).[47] G. Napoli and L. Vergori, Hydrodynamic theory for ne-matic shells: The interplay among curvature, flow, andalignment, Phys. Rev. E , 020701 (2016).[48] G. Napoli and S. Turzi, Spontaneous helical flows in ac-tive nematics lying on a cylindrical surface, Phys. Rev.E , 022701 (2020).[49] H. Shin, M. J. Bowick, and X. Xing, Topological de-fects in spherical nematics, Phys. Rev. Lett. , 037802(2008).[50] R. Zhang, Y. Zhou, M. Rahimi, and J. J. de Pablo, Dy-namic structure of active nematic shells, Nature Comm. , 13483 (2016).[51] S. Henkes, M. C. Marchetti, and R. Sknepnek, Dynamicalpatterns in nematic active matter on a sphere, Phys. Rev.E , 042605 (2018).[52] C. M. Dafermos, Disinclination in liquid crystals, Quart.J. Mech. Appl. Math. , 49 (1970).[53] H. Imura and K. Okano, Friction coefficient for a movingdisinclination in a nematic liquid crystal, Phys. Lett. A , 403 (1973).[54] L. M. Pismen and J. D. Rodriguez, Mobility of singular-ities in the dissipative ginzburg-landau equation, Phys.Rev. A , 2471 (1990).[55] G. Ryskin and M. Kremenetsky, Drag force on a linedefect moving through an otherwise undisturbed field:Disclination line in a nematic liquid crystal, Phys. Rev.Lett. , 1574 (1991).[56] T.C. Lubensky and Jacques Prost, Orientational orderand vesicle shape, J. Phys. II (France) , 371 (1992).[57] C. Denniston, Disclination dynamics in nematic liquidcrystals, Phys. Rev. B , 6272 (1996).[58] J.-M. Park and T. C. Lubensky, Topological defectson fluctuating surfaces: General properties and theKosterlitz-Thouless transition, Phys. Rev. E , 2648(1996).[59] D. R. Nelson, Toward a tetravalent chemistry of colloids,Nano Lett. , 1125 (2002).[60] D. Svenšek and S. Žumer, Hydrodynamics of pair-annihilating disclination lines in nematic liquid crystals,Phys. Rev. E , 021712 (2002).[61] E. I. Kats, V. V. Lebedev, and S. V. Malinin, Disclinationmotion in liquid crystalline films, J. Exp. Theor. Phys. , 714 (2002).[62] G. Tóth, C. Denniston, and J. M. Yeomans, Hydrody-namics of topological defects in nematic liquid crystals,Phys. Rev. Lett. , 105504 (2002).[63] A. M. Sonnet, Viscous forces on nematic defects, Cont.Mech. Thermodyn. , 287 (2005).[64] V. Vitelli and D. R. Nelson, Nematic textures in sphericalshells, Phys. Rev. E , 021711 (2006).[65] A. M. Sonnet and E. G. Virga, Flow and reorientationin the dynamics of nematic defects, Liq. Cryst. , 1185(2009).[66] L. Radzihovsky, Anomalous energetics and dynamics ofmoving vortices, Phys. Rev. Lett. , 247801 (2015).[67] R. Zhang, N. Kumar, J. L. Ross, M. L. Gardel, and J. J.de Pablo, Interplay of structure, elasticity, and dynamicsin actin-based nematic materials, Proc. Nat. Acad. Sci. (USA) , E124 (2018).[68] D. Khoromskaia and G. P. Alexander, Vortex formationand dynamics of defects in active nematic shells, New J.Phys. , 103043 (2017).[69] A. J. Vromans and L. Giomi, Orientational properties ofnematic disclinations, Soft Matter , 6490 (2016).[70] X. Tang and J. V. Selinger, Orientation of topologicaldefects in 2d nematic liquid crystals, Soft Matter , 5481(2017).[71] S. Shankar and M. C. Marchetti, Hydrodynamics of ac-tive defects: From order to chaos to defect ordering,Phys. Rev. X , 041047 (2019).[72] A. T. Brown, A theoretical phase diagram for an ac-tive nematic on a spherical surface, Soft Matter , 4682(2020).[73] M. Bowick, D. R. Nelson, and A. Travesset, Curvature-induced defect unbinding in toroidal geometries, Phys.Rev. E , 041102 (2004).[74] P. W. Ellis, D. J. Pearce, Y.-W. Chang, G. Goldsztein,L. Giomi, and A. Fernandez-Nieves, Curvature-induceddefect unbinding and dynamics in active nematic toroids,Nature Physics , 85 (2018).[75] D. J. G. Pearce, P. W. Ellis, A. Fernandez-Nieves,and L. Giomi, Geometrical control of active turbulencein curved topographies, Phys. Rev. Lett. , 168002(2019).[76] X. Tang and J. V. Selinger, Theory of defect motion in2d passive and active nematic liquid crystals, Soft Matter , 587 (2019).[77] G. Vertogen, The equations of motion for nematics, Z.Naturforsch. A , 1273 (1983).[78] G. Vertogen and W. H. de Jeu, Thermotropic liquid crys-tals, fundamentals (Springer-Verlag, Berlin, 1989).[79] A. M. Sonnet and E. G. Virga, Dynamics of dissipativeordered fluids, Phys. Rev. E , 031705 (2001).[80] M. Doi, Onsager’s variational principle in soft matter, J.Phys.: Cond. Mat. , 284118 (2011).[81] A. M. Sonnet and E. G. Virga, Dissipative Ordered Flu-ids: Theories for Liquid Crystals (Springer Science &Business Media, 2012).[82] B. J. Edwards, A. N. Beris, and M. Grmela, General-ized constitutive equation for polymeric liquid crystalspart 1. model formulation using the hamiltonian (poissonbracket) formulation, J. Non-Newton Fluid Mech. , 51(1990).[83] P. D. Olmsted and P. M. Goldbart, Isotropic-nematictransition in shear flow: State selection, coexistence,phase transitions, and critical behavior, Phys. Rev. A , 4966 (1992).[84] P. D. Olmsted and C.-Y. D. Lu, Coexistence and phaseseparation in sheared complex fluids, Phys. Rev. E ,R55 (1997).[85] P. D. Olmsted and C.-Y. D. Lu, Phase separation of rigid-rod suspensions in shear flow, Phys. Rev. E , 4397(1999).[86] G. Napoli and L. Vergori, Surface free energies for ne-matic shells, Phys. Rev. E , 061701 (2012).[87] K. Kruse, J. F. Joanny, F. Jülicher, J. Prost, and K. Seki-moto, Asters, vortices, and rotating spirals in active gelsof polar filaments, Phys. Rev. Lett. , 078101 (2004).[88] D. Marenduzzo, E. Orlandini, M. E. Cates, and J. M.Yeomans, Steady-state hydrodynamic instabilities of ac-tive liquid crystals: Hybrid lattice boltzmann simula-tions, Phys. Rev. E , 031921 (2007). [89] T. B. Liverpool and M. C. Marchetti, Hydrodynamicsand rheology of active polar filaments, in Cell Motility ,edited by P. Lenz (Springer, New York, 2008) pp. 177–206.[90] A. W. C. Lau and T. C. Lubensky, Fluctuating hydrody-namics and microrheology of a dilute suspension of swim-ming bacteria, Phys. Rev. E , 011917 (2009).[91] L. Giomi, L. Mahadevan, B. Chakraborty, and M. F. Ha-gan, Banding, excitability and chaos in active nematicsuspensions, Nonlin. , 2245 (2012).[92] S. Ramaswamy, The mechanics and statistics of activematter, Annu. Rev. Cond. Mat. Phys. , 323 (2010).[93] I. W. Stewart, The Static and Dynamic Continuum The-ory of Liquid Crystals (Taylor & Francis, 2004). [94] See Supplementary Movie at http://link.aps.org/supplemental/10.1103/PhysRevE.102.012607 for themovie illustrating the time evolution of a periodic tra-jectory followed by four +1 / defects in an active ne-matic on a sphere, for the parameters R/r = 16 and γα/ Kζ = − . . The arrow represents the orienta-tion of the defect. And we speed up the motion in themovie.[95] L. Metselaar, J. M. Yeomans, and A. Doostmohammadi,Topology and morphology of self-deforming active shells,Phys. Rev. Lett. , 208001 (2019).[96] W. Helfrich, Elastic properties of lipid bilayers: the-ory and possible experiments, Z. Naturforsch. C28