Active topological defect absorption by a curvature singularity

We leverage the Born–Oppenheimer approximation to present a general description of topological defects dynamics in p-atic materials on curved surfaces. Focusing on the case of an active nematic, we find that activity induces a geometric contribution to the motility of the +1/2 defect. Moreover, in the case of a cone, the simplest example of a geometry with curvature singularity, we find that the motility depends on the deficit angle of the cone and changes sign when the deficit angle is greater than π, leading to the change in active behavior from contractile (extensile) to extensile (contractile) behavior. Using our analytical framework, we then identify for positively charged defects the basin of attraction to the cone apex and present closed-form predictions for defect trajectories near the apex. The analytical results are quantitatively corroborated against full numerical simulations, with excellent agreement when the capture radius is small compared to the cone size.

In contrast to its passive counterpart, activity renders a comet-shaped +1/2 defect motile, driving it to selfpropel along its axis of symmetry, in a direction dictated by the sign of the active stress: extensile (contractile) active stress drives the defect to move towards the head (tail) of the comet [2].A trefoil-shaped −1/2 defect, on the other hand, due to its three-fold rotational symmetry does not become motile, unless external forces or boundary conditions violate the three-fold symmetry.
Previous theoretical work has described multi-defect active nematics by treating the topological defects as quasiparticles, with elastic distortions of the nematic texture and active flows mediating the effective interactions in phenomenological models [4,19,20].Recently, a systematic formulation based on the Born-Oppenheimer approximation was employed for active nematics [21][22][23] by treating the defect positions as slow degrees of freedom and the nematic texture as a fast degree of freedom that instantaneously responds to the slow motion of the defects.
The discussion thus far has focused on flat surfaces; however, much of biology takes stage on curved surfaces.Specifically, various natural and synthetic developmental processes involve formation of three-dimensional curved structures from cell monolayers.Striking examples include epithelial dome formation [24], embryo gastrulation [25], and in-vitro tissue regeneration [26].In addition, more recently, nematic organization of active entities on flat surfaces has been linked to threedimensionalization processes and in particular to the formation of protrusions in confined myoblast cells [27], from cytoskeleton of animal Hydra during morphogenesis [15], and from the membrane of eukaryotic cells in the form of cellular fingers, known as filopodia [28].
Some progress has been made on curved surfaces in the passive context.Ref. [29] analytically derived the ground state defect configurations on a sphere and torus, and the interaction between liquid crystalline order and curved substrates was studied in [30,31], where it was shown that curvature gives rise to an effective topological charge density.For a cone, this corresponds to negative topological charge concentrated at the apex, and a simple argument was recently presented in [32] for the case of free boundary conditions.Ref. [33] re-derived the induced charge result of Vitelli and Turner [31] and in the context of a cone with tangential boundary conditions used it to determine the ground state defect configuration given a fixed number of topological defects.As a cone, which has a curvature delta function singularity at the apex, is the simplest example of nontrivial curved geometry, we continue to study the dynamics of active nematics on a cone geometry here.
This paper is organized as follows.We begin in Sec.II by formulating a minimal model of a general p-atic texture on an arbitrarily curved surface.Working deep in the ordered limit, and taking advantage of isothermal coordinates (recently introduced in the context of liquid crystals [33,34]), we write down explicitly the quasistatic multi-defect solution in the passive setting on a curved surface.Then in Sec.III we study the dynamics of the texture by working within the Born-Oppenheimer approximation: we assume that the defects move slowly (valid in the limit of low defect density and low activity) and that the nematic texture instantaneously readjusts itself in response.The partial differential equation for the nematic texture is thus reduced to a set of ordinary differential equations for the effective defect positions.Upon introducing activity in Sec.IV, we first derive the hydrodynamic equations of a compressible active nematic film on a curved surface.Working within the Born-Oppenheimer approximation (as in the passive case) yields a number of new results.In particular, applying the framework to the case of active nematics we show that the motility of a +1/2 defect picks up an active geometric contribution, and changes sign as the curvature is increased.For a +1/2 defect on a cone, we identify the basin of attraction to the apex and predict its trajectory.Throughout the paper, we quantitatively check our theoretical predictions against full numerical simulations of active nematics on a cone, explore the effect of boundary conditions, and scrutinize the validity of the Born-Oppenheimer approximation, finding excellent agreements.Sec.V summarizes our main results and comments on applications.Most of the technical details are relegated to Appendices A-D.

II. FORMULATION OF A MINIMAL MODEL IN ISOTHERMAL COORDINATES
Here we review the framework for describing a p-atic texture on a curved surface, following the presentation in Ref. [33].The case of a nematic texture is recovered by simply setting p = 2.

A. Metric
In two dimensions it is always possible to choose local complex coordinates z and z, known as isothermal (or conformal) coordinates, such that the length of the interval squared can be written as [35], where e ϕ is the conformal factor that describes positiondependent isotropic stretching, and g z z = e ϕ /2 is the metric.In Cartesian coordinates, from which follows that ds 2 = e ϕ(x,y) (dx 2 + dy 2 ). ( Conical geometry.The conical geometry is a prime example that we study in this paper.For a cone with half angle β, the metric can be obtained from where χ = 1 − sin β.As such, 2πχ describes the deficit angle of the cone, and for example χ = 0 corresponds to a disk (no fraction missing).Another set of useful coordinates are physical coordinates z, which correspond to unrolling a cone to form a planar disk with a missing sector, are related to z via the coordinates In these coordinates, which is flat everywhere except at the origin.The distances from the apex in these two coordinate systems are related by At the origin, z is not defined and there is a conical singularity with deficit angle 2πχ.See Fig. 1 for a schematic of various coordinate systems for a cone.

B. Tensors and covariant derivatives
Let T denote a traceless real symmetrized rank-p tensor that describes the p-atic order.Then in isothermal coordinates, since T is traceless (contraction of any pair of indices vanishes), T has only two non-zero components T ≡ T z...z and T ≡ T z...z , where here ellipses denote p copies.Also, by reality, T = ( T ) * .
For ease of notation, let ∇ ≡ ∇ z and ∇ ≡ ∇ z denote the covariant derivatives with respect to z and z, respectively.Covariant derivatives of T are quite simple using isothermal coordinates where ∂ ≡ ∂/∂z and ∂ ≡ ∂/∂ z.Note the asymmetry between the first and second equations of Eq.(8a) and Eq.(8b).In Cartesian coordinates, we have and See Table I for the correspondence between isothermal coordinates and Cartesian coordinates.

C. Free energy
The free energy F (in a one Frank constant approximation for p = 2) in terms of expansion in powers of the order parameter T and its gradients can be written as (11) where explicitly Here K, K > 0 are Frank elastic type terms [36,37], and in regions of zero Gaussian curvature (such as any point on a cone other than the apex), the two terms are equivalent by integration by parts.The last term is the Landau-De Gennes type free energy that governs the isotropic-p-atic transition, with controlling the microscopic coherence length and S 0 sets the equilibrium magnitude of the p-atic order.Without loss of generality we set S 0 = 2 p .On using the free energy in Eq. ( 11), the p-atic order parameter T can be determined by minimizing the free energy.Deep in the ordered limit ( 1), we have Writing the order parameter in terms of its amplitude A and phase α, we find T z...z = A z...z e iα = Ae iα , which leads to In other words, the contribution of the potential term to the free energy vanishes, and the free energy simplifies to [38] Upon substitution of Eq. ( 15) into Eq.( 16), the free energy can be written in terms of ϕ (the log of the con-symbol in isothermal coordinates physical interpretation Cartesian representation (z, z) coordinates (x + iy, x − iy) formal factor) and the order parameter's phase α as where we have used Minimizing F (Eq. ( 17)) with respect to the order parameter phase α gives In the presence of a topological defect of charge σ ∈ Z/p, the phase α will wind by 2πpσ.Thus a solution to Eq. ( 19) with a single defect j at z j with charge σ j is (neglecting image charges that allow us to impose various boundary conditions) which results in the energy of a single p-atic topological defect of charge σ j [33] Note that ϕ(z j ) is a purely geometrical contribution that is set by the shape of the curved space on which the patic texture resides.Recall from Eq. ( 4) that, for a cone with half-angle β, we have Having set up the isothermal coordinates and having defined the energy of topological defects for p-atics on arbitrary geometries, we now turn to defect dynamics.Our goal is to establish an effective description of topological defects dynamics through a set of ODEs that describe Langevin-type behavior for the positions of topological defects that behave like charged quasi-particles, where the liquid crystal textures relax instantaneously.We begin with the passive case, and then extend to active defect dynamics.

III. EQUATIONS OF DEFECT DYNAMICS: PASSIVE CASE
We begin by considering the passive case, while keeping the description general for p-atic on an arbitrary geometry with conformal factor e ϕ(z,z) .We obtain p-atic dynamics in the limit of low defect density by the following procedure: we minimize the mean squared deviation of the dynamics on the inertial manifold T 0 from the exact equation of motion.Specifically, we assume that the defects move slowly and that in response the p-atic texture instantaneously readjusts itself.When studying the quantum mechanics of light-weight electrons bonding atoms with much heavier nuclei, this is known as the Born-Oppenheimer approximation [39].Similar approximations were made to study the dynamics of vortices in superfluid helium films driven out of equilibrium by an oscillating substrate [40] and the dislocation-mediated elongation of the peptidoglycan cell walls of bacteria [41].

A. Born-Oppenheimer approximation
In the passive case, and in the absence of active flows, we assume the exact equation of the motion for the defect phase α is controlled by the relaxational dynamics where γ is the rotational diffusion coefficient and the dependence on geometry is manifest through the metric g.We take as our Born-Oppenheimer ansatz [22,23], where z i (t) are time-dependent positions for defects, and find z i (t) by minimizing the mean squared deviation between the dynamics on the inertial manifold and the exact equation of motion (Eq.( 23)), with respect to żi , where α 0 is given by Eq. ( 20).Doing so [22] leads to the following coupled set of simple ODEs for the defect dynamics in isothermal coordinates where − ∂F0 ∂ zi is the Coulombic force computed by differentiating the Coulombic potential (Eq. ( 21)).(Were we to include thermal fluctuations, Langevin noise sources would appear in Eq. ( 26).)The tensors M ij and N ij are collective mobility matrices that describe the effective response of a topological defect to the Coulombic forces, and are geometry dependent where g is the metric.Taken together, Eqs. ( 26)-(27b) along with Eq. ( 20) and Eq. ( 21) for α 0 and F 0 , respectively, describe the dynamics of passive topological defects of charge σ i located at z i for a p-atic texture on a curved geometry that is described by the metric g z z = e ϕ /2.We next apply this framework to first describe defect dynamics on a cone geometry, and will then extend the framework to investigate the effect of activityinduced flows.

B. Evaluation of dynamical equation
We first compute the collective mobility M ii and N ii for a single defect of charge σ i located at z i on a cone, as the simplest example of nontrivial curved geometry.Up to now, our discussion has been general.Here we will continue to be general (for example, write in terms of p and generic defect charge σ i ), but for simplicity we consider the geometry of a cone with cone apex at the origin, described by ϕ = − ln(z z) and the cone radius R a 1/(1−χ) , where a is the defect core size.We assume that the defect is sufficiently far from the cone base so we can ignore any image charges, despite the long-range interactions.The computations of M ii and N ii are presented in Appendix B, and since N ii is subleading to M ii , we shall neglect N ii .
When the defect is sufficiently far from the apex, i.e., and thus żi = −3γ −1 (K + K )χ 1 When the defect is close to the apex, instead we find and thus where δ = ar −χ i ∼ a 1/(1−χ) for r i ∼ δ.Taken together, the equations of defect dynamics for a p-atic texture on a cone with deficit angle 2πχ become ) where r i = |z i |, a ∼ δ 1−χ is the defect core size near the apex, R is the cone radius, γ is the orientational diffusion, and K, K denote orientational elasticities.We comment that in the case of a disk (χ = 0), the mobility reduces to the well-known defect friction [42].

C. Global vs. local defect position
It is important to elaborate upon what we mean by the defect position z i .The standard definition of the defect position is determined by the zero of the tensor order parameter T where the phase of T winds around.This gives a local definition of defect position which does not depend on the profile of the texture far away.However, one might instead be interested in a defect position determined by fitting the global texture to a defect profile ansatz.This global definition of course may be different from the local definition.The equations that we have derived have been obtained by fitting the global ansatz and thus the defect dynamics given by Eq. ( 32) depends on the global texture, and it not necessarily the same as the local defect position.
To the extent that our ansatz for the defect profile is good locally, the two notions of local vs global defect positions should give similar results.When the two no-tions differ substantially, the Born-Oppenheimer ansatz breaks down for the local description of the texture.We will check when these two notions agree and disagree.

D. Comparison with simulations
In order to check our analytical formulation, we compare the dynamics of a single topological defect of charge σ calculated from Eq. ( 26) using the mobility M ii from Eq. (30), with full numerical simulations of the passive nematic, i.e. p = 2, texture on a cone.
Although we wrote our equations in terms of the phase α of the order parameter, in our simulations we numerically evolve the full nematic texture T = Q (p = 2) itself explicitly, according to where F is given in Eq. ( 11), S 0 = 2 p , and Q is the nematic tensor, which in isothermal coordinates can be written in terms of its Cartesian components as Without loss of generality, since a cone has zero Gaussian curvature everywhere except at the apex, we set K = 0 in the simulations.We choose periodic boundary conditions on a square box, the simplest boundary conditions to simulate.We expect the choice of boundary conditions to not be important for trajectories and times such that the defect distance r(t) to the apex is much smaller than the size of the box R, i.e. r(t) R. We also checked this expectation by checking that defect trajectories (discussed later in the paper) don't change if we change the boundary conditions to free [32] or tangential boundary conditions [33] at the boundary of a disk.By tangential boundary conditions, we mean that the nematic director is tangential to the boundary, and by free boundary conditions, we mean that the radial component of the gradient of Q vanishes.Explicitly, we impose at the boundary of a square grid of the isothermal coordinate space of Fig. 1(a), where the boundary is given by z = Re iφ and φ is the azimuthal coordinate): where ∇ n is the covariant derivative in the direction normal to the boundary.In our simulations, since defects are initially placed near the apex and far from the boundary, the precise choice of boundary condition should not matter for defect trajectories that remain near the apex.
In our simulations, we solve Eq. ( 33) numerically us-simulation analytics FIG. 2. Initial physical global defect velocity of a +1/2 defect at x/R = −0.08 in the presence of another +1/2 defect already at the origin.The blue points are obtained by fitting the texture at each time to the ansatz and determining zi(t) which minimize the mean squared deviation E in Eq. ( 25).The red curve is fit of Eq. ( 36) with a/( √ K + K ) = 0.5.The y-axis is rescaled by a passive characteristic velocity ṽ0 ing the method of lines [43].The temporal evolution is performed through a predictor-corrector scheme [44] and spatial derivatives are evaluated using five-point stencil central differences.
For χ > 0, a defect of charge σ is attracted to the apex with geometric force proportional to χ(σ − 1 2 σ 2 )/r [33], which is indeed what we observe in our simulations: once a defect of charge σ is absorbed by the apex, then the effective charge at the apex is reduced.Now the net force on any remaining charge or defect on the cone flank is proportional to The force is repulsive for sufficiently small χ, and attractive for sufficiently large χ; in particular, the force vanishes at χ = σ 2 σ− 1 2 σ 2 .For σ = 1/2 (relevant to a p = 2 nematic in the one Frank constant approximation), the force vanishes at χ = 2/3.For a passive +1/2 defect at the apex, upon setting p = 2 and σ i = 1/2, Eq. ( 32) becomes We checked this formula (with K = 0) using simulations (see Fig. 2), where we find that the force does indeed change sign for χ ≈ 2/3.

IV. EQUATIONS OF DEFECT DYNAMICS: ACTIVE CASE
So far we have limited the description to the passive dynamics of topological defects, where the motion is controlled by the relaxational dynamics of the p-atic texture (Eq.( 23)).The effect of activity enters through coupling of the relexational dynamics to the velocity field that is locally generated by the active stresses.In the overdamped limit [45][46][47], the velocity field can be directly calculated from the balance of active force and the frictional damping where µ is the friction coefficient and σ active denotes activity-induced stresses.For a general case of a p-atic, the active stress can comprise of terms proportional to the contractions of p-atic tensor T that are allowed by symmetry [48].The active stress can therefore be considered as a function of the p-atic order parameter tensor T as σ active = σ active [T ].It thus follows from Eq. ( 38) that in the general case, the activity-induced velocity field v can be written in terms of the p-atic order parameter as v = v[T ].While here we keep the description general for active dynamics of any p-atic texture on an arbitrary geometry, in Sec.IV B we present a detailed analyses of the dynamics for active nematics (p = 2) on the cone and compare the theoretical predictions with full numerical simulations of the active nematics.

A. Derivation of dynamical equation
Once the active velocity field is found from Eq. ( 38), we can determine the contribution from the active velocity to the dynamics of the defect phase.To this end, Eq. ( 23) is modified to account for balancing activity-induced flows and relaxational dynamics where D/Dt is the advective derivative.Naively, one would expect However, this is not correct as α is anomalous: indeed, we are not advecting a conventional scalar!Under rotation of coordinates z → e i∆α z, α transforms as α → α + p∆α.
We therefore need to derive the correct form of Dα/Dt.
We start from the advective derivative of a p-atic tensor, and upon substitution of T = e − pϕ 2 +iα (Eq.( 15)), we find that The last term in the above equation arises from the geometric contributions to the advection and the corotation, and cancels out only for p = 2, i.e. nematics.On dividing through by iT we obtain where the velocity field v in our overdamped limit is set by the active stresses (Eq.( 38)).
With Eq. ( 43), we can now obtain active dynamics in the limit of weak activity and low defect density by following the same procedure as in Sec III: we minimize the mean squared deviation of the dynamics on the inertial manifold T 0 from the exact equation of motion Eq. ( 43).This procedure [22] leads to the following coupled set of ODEs for the active defect dynamics and the only difference from the passive case is an important new term which accounts for the contribution of the activityinduced velocity field, where We emphasize that Equations ( 44)-( 46) present a general description of the dynamics of active topological defects for any p-atic texture on an arbitrary geometry characterized by the metric g z z = e ϕ /2.Within this general description, the defect phase α 0 , the free energy F 0 , and the velocity field v are, respectively B. Active stress for a nematic While for a general p-atic, multiple contributions to the active stress are allowed by virtue of the p-fold rotational symmetry of the order parameter [48], the most well-established form of the active stress is the stresslet contribution, which has been extensively observed in biological systems [2,16] and studied theoretically.The stresslet contribution leads to dipolar active forces and is proportional to the nematic (p = 2) tensor where ζ Q is the scalar activity coefficient that characterizes the strength of the stresslet (and corresponding force dipole), and Q zz is the nematic tensor, i.e. a special case of the generic rank p tensor T in Sec.II for p = 2.As such, the nematic tensor in isothermal coordinates Q zz can be expressed in terms of its more familiar Cartesian components through Eq. ( 34).Therefore, using the definition of the active stress in terms of the stresslet contribution (Eq.( 50)), the velocity field in the overdamped limit (Eq.( 38)) can be written in terms of the nematic tensor in isothermal coordinates as From this point onward, for ease of notation, we will suppress z indices.Thus we write the activity-induced velocity field as v = ζ∇Q using this simplified notation.
Active nematic defects in flat geometry.In order to apply the above framework to active nematics, we first review and rederive the well-known motile force for a +1/2 defect in flat space, where the geometric potential ϕ = 0.In complex coordinates z and z, the nematic tensor T = Q for a defect of charge σ i = ±1/2 at the origin takes the form Fig. 3 presents the vector field corresponding to the Consistent with the theoretical [49], numerical [2], and experimental [14] velocity fields for active nematic defects on a flat space, the flow field for a comet-shaped +1/2 defect points along the head-tail axis of the comet.However (unless external stresses or boundary conditions break the 3-fold rotational symmetry), on average any net flow for a −1/2 vanishes.These results are consistent with the intuition that a +1/2 defect, which has an asymmetric comet shape, moves along the comet axis (with a direction depending on the sign of ζ), whereas a −1/2 defect, which has three-fold rotational symmetry, doesn't move by symmetry.
Active nematic defects on a cone.We can now generalize the dynamics of active nematic defects to curved space.As before, v = ζ∇Q, where the sign of ζ dictates contractile vs. extensile activity; in particular, Here we will focus on the case of a cone, in which case where 2πχ is the deficit angle (see Fig. 1).We are as usual working in isothermal coordinates z and z.Upon focusing on the case of a single +1/2 nematic defect, we have , where ϕ = 0 reflects the fact that the metric is no longer flat.
We comment that in the far-field, which is equivalent to taking a defect at the apex, i.e. z i → 0, v z reduces to  i.e. cones formed from disks with exactly half their area removed.Thus v changes sign as χ goes from below 1/2 to above 1/2, and the entire nematic texture moves accordingly along the real axis of our isothermal coordinate system.The vanishing velocity for χ = 1/2 is reminiscent of (but not the same as) the passive case, where the apex charge is completely screened for χ = 1/2.See Fig. 4 for plots of the velocity profile for a +1/2 defect for χ = 0.25 and χ = 0.75.Although the velocity fields generated by activity are similar near the cone apex, they are quite different in the far field.

C. Evaluation of dynamical equation
With the activity-induced velocity-field discussed, the active contribution to the dynamics of nematic defects can be calculated using Eq. ( 45), where setting p = 2 and σ = +1/2 for active nematics leads to (Appendix C) where a is a short distance cutoff we take to be the spacing between liquid crystal molecules.Note that I vanishes for χ = 1 (corresponding to a long cylindrical shell) and χ = 1/2.We now combine the active velocity contribution, together with the calculations of the collective mobilities on the cone in Sec.III (Eqs.( 28)-( 30)), and the equations of motion (Eq.( 44)) for an active nematic +1/2 topological defect on a cone geometry to obtain our final equation for +1/2 defect dynamics on a cone, In terms of the physical coordinates shown in Fig. 1 we have It is easy to interpret the different contributions to the defect dynamics displayed in Eq. (57) (or equivalently, Eq. ( 59)): The LHS is the mobility, with a logarithmic correction that depends on the distance to the apex, unchanged by the activity.The first term on the RHS is the motile force caused the activity parameter ζ, and the second term on the RHS is the attractive Coulomb force that a +1/2 defect feels from the cone apex with effective topological charge −χ.The geometric corrections to the mobility, active force, and the interaction force are manifest through the dependence on the deficit angle 2πχ of the cone.In particular, the motility term reverses sign for χ > 1/2, consistent with the velocity profile in Eq. (55).In addition, the motility speed is reduced by a factor of 1 − χ.Note that the angle of self-propulsion depends on the position, and, remarkably, for a critical defect position zi , a +1/2 defect can have a stationary radius, which occurs when the RHS vanishes.This is where the Coulomb and motile forces balance each other and (as discussed further below) defines a basin of attraction for the positive defects.
The radius of the basin of attraction (when the magnitude of the RHS vanishes), is Note that the critical radius rc diverges as the activity ζ tends to zero.Interestingly, for finite activity ζ, the defect is always absorbed by the apex in the limit χ → 1, which is the limiting geometry of a cylinder, regarded here as an extremely pointed cylinder.

D. Defect trajectories
Having derived the equation of motion for the nematic defects in Eq. (57) (or equivalently, Eq. ( 59)), we then calculate defect trajectories.To summarize, the dynamical equations for a +1/2 defect in isothermal coordinates z and physical coordinates z, read, respectively, where for a particular defect i we define Dividing Eq. ( 61) by its complex conjugate results in For χ = 0, we recover the known result of dz = dx + idy = dz, i.e., a defect with dy = 0 moving parallel to the real axis, which we assume is aligned with the "comettail" of the +1/2 defect.We now determine dynamical trajectories to leading order in χ.From Eq. (64), B = O(χ), leading to Upon integrating after cross-multiplying (and ignoring the O(χ) corrections to the above equation), we have where z 0 is the initial position of the defect.Redefining our coordinates as w = A B z and rearranging leads to In terms of w = w R + iw I = re iφ , we arrive at where each set of initial conditions (w I (t = 0), φ(t = 0)) gives a different trajectory.Upon remembering that we can compute defect trajectories, up to corrections of O(χ 2 ).Fig. 5 illustrates the contours of defect trajectories around the cone apex, showing that the defect can get attracted to or repelled from the apex depending on the original position.

simulation analytics
FIG. 6.Initial physical global defect velocity due to motility for a +1/2 defect initially near the apex (r(t = 0) R).Blue points are obtained by fitting the texture at each time to the ansatz and determining zi(t) which minimize the mean squared deviation E in Eq. (25).The red curve is fit of Eq. ( 75) from motile term with δ/R = 0.04 and the size of region of optimization is R /R = 0.08.The y-axis is rescaled by the characteristic active velocity ṽ0 = ζ/ R where R is the physical radius of the base of the cone.

E. Validity of the Born-Oppenheimer
approximation for χ ≥ 1/2 It is important to note that for χ ≥ 1/2, since the motility for the global defect texture reverses sign, then it is as if the global texture is moving in the opposite direction of the local defect position.In other words, the Born-Oppenheimer ansatz in the near-field of the defect is not a good approximation for χ ≥ 1/2.It is possible that the Born-Oppenheimer approach could nevertheless be useful for understanding global properties of the texture, far away from the defect.

F. Comparison with simulations
When activity is included, as a check on our analytic results, we numerically evolve the full nematic texture Q according to where ∂ t Q is given in Eq. ( 33) and v = ζ∇Q.
Sign of global defect velocity.Simulations indeed confirm that the local defect velocity never changes direc- tion as χ is varied, as predicted by dynamical equations such as Eq.(57) (or equivalently, Eq. ( 59)).From simulations, we can also extract the global defect position by fitting the defect ansatz to the global texture data.When the defect is near the apex, i.e. r i ∼ δ, the global defect velocity due to motility (i.e., ignoring the Coulomb term) is obtained by replacing the mobility in Eq. (57) with Eq. ( 30), leading to the global defect velocity, where we used δ = ar χ i in the above equation.The global velocity does indeed change sign as χ increases from below 1/2 to above 1/2.See Fig. 6 for comparison of simulations with theoretical prediction of global defect velocity due to the motility term of Eq. (75).The fact that the motility flips sign indicates that the behavior of extensile and contractile activity flips at χ = 1/2, i.e., when the integral of the total Gaussian curvature associated with the cone apex is π.
Defect trajectories.For small χ (χ = 0.1) and small activity (ζ/[γ −1 (K + K )] = 2 × 10 −2 ), we used simulations to evolve the nematic texture consisting of a defect initially near the apex.As before, we extract the global defect position by fitting the defect ansatz to the global texture data and plot in Fig. 7. Eq. (72) predicts that all of the slopes are the same for global defect positions.Comparison to numerical simulations confirms this prediction for defects close to the cone apex (Fig. 7; blue, green, orange, and red lines), with the slope m ≈ 0.3R.On setting B/A = m and solving for a, we find a/( √ K + K ) = 0.67.Simulation results further show that the prediction for the slope becomes worse as the initial distance of the defect to the apex increases (Fig. 7; brown and grey lines), which is consistent with the fact that as time passes we expect the ansatz to become less accurate.

V. DISCUSSION
With the help of the Born-Oppenheimer approximation, which assumes liquid crystal textures can relax instantaneously to defect positions, we have presented a description of topological defects dynamics in compressible active nematic materials deep in the ordered phase on curved surfaces with a focus on the behavior near the apex of cones.We find that activity induces an active geometric contribution to the motility of the global texture of the +1/2 defect which changes sign as the curvature increases, a result we expect to hold more generally.Furthermore, in the case of a cone, we provide a closed form prediction for the radius of the basin of attraction around the apex, and present analytical description of the defect trajectories near the apex.Sufficiently far from the cone base, the analytical results agree well with full numerical simulations.
We would like to emphasize that although our analysis focused on the interaction between an active nematic defect and the cone apex, the formulation presented here is general and can be applied to geometries with more general conformal factors and p-atic textures.Furthermore, since we assumed that the defect was near the apex and far away from the boundary of the cone, we could safely ignore the effect of boundary conditions.The framework can be extended to study defect interactions with the boundaries through introducing image charges, which is the focus of our upcoming study; preliminary investigations have revealed a rich phase diagram of allowed dynamical states, which exhibits not only single flank defect orbits and two flank defect orbits, but also transitions between them via defect absorption, defect emission, and defect pair creation via activity at the apex [50].We expect significantly different results for defect dynamics closer to the boundary, with, say, tangential as opposed to free boundary conditions.
In a related direction, it would also be worth exploring defect configurations on cones with both tangential and free boundary conditions at finite temperatures.With increasing temperatures, entropic effects might cause the cone apex to emit defects or induce defect unbinding.Finally, we note that the formulation presented here could be a useful tool for studying topological defects interacting on non-trivially curved surfaces that biologically active entities experience in various setups.One striking example is the interaction of topological defects with dynamic protrusions that are formed from cellular membranes [28], cell monolayers [27], or the cell cytoskeleton [15].Such protrusions are often characterized by non-trivial, dynamically changing, surface curvature, and even sharp tips, which according to our analyses can result in dynamic geometric contributions to the motility of active topological defects.We find in upcoming work [50] that activity can catalyze sharp, lightning-rodlike conical geometries to create and emit defects.
Furthermore, there is currently intense research in designing in-vitro surfaces for tissue regeneration [26].Topological defects are also finding increasing applications in governing biological functions such as cell death and removal, and cell differentiation.In this vein, our predictive framework can potentially find applications in designing surface geometries that allow for prescribed locations of topological defects.
Upon shifting z → z + z i , Thus, It is interesting that for all χ, the motility is reduced, and in fact reverses sign for χ > 1/2.It is also interesting to note that for χ = 1/2 and χ = 1, I 1 vanishes, in which case a +1/2 defect is not motile.
Combining everything, We now calculate the time-dependence of the defect position for a defect trajectory along the x-axis.We have which is integrated to give trajectory x(t) which satisfies where Li 2 is the dilogarithm.

FIG. 1 .
FIG. 1. Schematic of coordinate systems for a cone: (a) Isothermal coordinates z = re iφ .(b) Physical coordinates z = re i φ, corresponding to an unrolled cone with a missing sector, with angle fraction χ.(c) A diagram of the cone in 3D with cone half angle β, where sin β = 1 − χ.
55) which reveals the slow falloff of the activity induced velocity field and the special nature of cones with χ = 1/2,

FIG. 4 .
FIG.4.Velocity profile produced by a contractile (ζ > 0) +1/2 defect on a cone in the isothermal coordinates of Fig.1for a +1/2 defect at (1, 0) for χ = 0.25 (top row) and χ = 0.75 (bottom row).The black dot represents the apex and the green dot represents the +1/2 defect.Note that the velocity profile near the apex is similar for both χ = 0.25 and χ = 0.75 (left column), but far away (right column, zoomed out) the velocity profile flips sign for χ = 0.75.The colormap represents velocity magnitude normalized by its maximum value, thus ranging from 0 (dark blue) to 1 (dark red).

FIG. 5 .
FIG.5.Streamlines of defect trajectories for small χ = 1/4.Top row: streamlines in isothermal coordinates z, where (a): small χ analytic trajectories, given by Eq. (72) and (b): numerical trajectories obtained by numerically integrating the exact ODEs, given by Eq. (68).(c): trajectories in physical coordinates wrapped around a 3D cone.(d): trajectories on an unrolled cone in z coordinates where black arrows indicate identification of the two black edges.In all of the plots, the red dot denotes the cone apex, and the green point is where the motile and Coulomb forces on the +1/2 defect balance each other.Any point within a distance given by the radius of the red circle given by rc in Eq. (60) gets attracted to the apex, and the magenta curve is an incoming separatrix.We expect that trajectories outside of the apex domain of attraction go down the cone flanks and are eventually influenced by the boundary conditions at the cone base.Other parameters: A = B = a = 1.

7 .
Left plot: various defect trajectories using isothermal coordinates obtained from simulations and our analytical approximation.Colored markers are the local defect positions determined by simulation (each one corresponding to a different initial position), where the black asterisk markers desimulations that used free boundary conditions for comparison.curves are contours (Eq.(72)) from the analytic approximation of the global defect position.Trajectories on the left side are attracted to the cone apex, while those on the flanks run away toward the base, where they will eventually be influenced by the boundary conditions at the cone base.The cone apex is at the origin (denoted by the black dot).The right plot contains linear best fits for each trajectory.Simulation parameters: χ = 0.1 and ζ/[γ −1 (K + K )] = 2 × 10 −2 , and a/( √ K + K ) = 0.67.Radius of the disk in z coordinates is R = 128, but since this size is large enough, it does not affect the trajectories near the apex.

TABLE I .
Definitions and correspondence between isothermal coordinates and Cartesian coordinates.