Curvature controlled defect dynamics in topological active nematics

We study the spatiotemporal patterns that emerge when an active nematic film is topologically constraint. These topological constraints allow to control the non-equilibrium dynamics of the active system. We consider ellipsoidal shapes for which the resulting defects are 1/2 disclinations and analyze the relation between their location and dynamics and local geometric properties of the ellipsoid. We highlight two dynamic modes: a tunable periodic state that oscillates between two defect configurations on a spherical shape and a tunable rotating state for oblate spheroids. We further demonstrate the relation between defects and high Gaussian curvature and umbilical points and point out limits for a coarse-grained description of defects as self-propelled particles.

Active systems are characterized by constant input of energy, which is converted by autonomous constituents into directed motion, leading to spatiotemporal patterns.These phenomena range from the macroscale, e.g.flocks of birds 1 or schools of fish 2 to the microscale, e.g.bacterial colonies 3 , migrating tissue cells 4 or active nematic films 5 .If such systems are confined on curved surfaces, topological constraints strongly influence the emerging spatiotemporal patterns.Using these topological constraints to guide collective cell behavior might be a key in morphogenesis and active nematic films on surfaces have been proposed as a promising road to engineer synthetic materials that mimic living organisms 6 .However, the complex dynamics of such topological active systems remains wildly unexplored.As in passive systems the mathematical Poincaré-Hopf theorem forces topological defects to be present in the nematic film.On a sphere this leads to an equilibrium defect configuration with four +1/2 disclinations arranged as a tetrahedron [7][8][9] , see Figure 1 The disclinations repel each other and this arrangement maximizes their distance.In active systems unbalanced stresses drive this configuration out of equilibrium.But in contrast to planar active nematics with continuous creation and annihilation of defects 10,11 the creation of additional defect pairs can be suppressed on curved surfaces, which is demonstrated in 6 for an active nematic film of microtubules and molecular motors, encapsulated within a spherical lipid vesicle.This provides an unique way to study the dynamics of the four defects in a controlled manner and led to the discovery of a tunable periodic state that oscillates between the tetrahedral and a planar defect configuration.We confirm this finding by computer simulations, see Figure 1.
Within a coarse-grained model +1/2 disclinations in planar active nematic films can be effectively described by self-propelled particles with a velocity proportional to the activity 5 .
In 6 this relation is extended to spherical nematics.Four self-propelled particles on a sphere also oscillate between the planar and tetrahedral configuration.Both descriptions can be quantitatively linked to each other, but also differences can be pointed out, which become more evident for more general surfaces.For non-constant Gaussian curvature constraints local geometric properties influence the position of the defects and thus can be used to control defect dynamics.We are concerned with a systematic investigation of the impact of such constraints on the emergence of complex patterns and oscillations.

RESULTS
For active systems in flat geometries various theoretical descriptions have been proposed, see e.g. 12,13.One of the most studied approaches are Vicsek-like models 14 .They consider particles, which travel at a constant speed to represent self-propulsion, whose direction changes according to interaction rules which comprise explicit alignment and noise.In contrast to equilibrium systems long-range order emerges for two dimensional systems with low noise.We consider an extension of these models which includes excluded volume [15][16][17][18] and classify systems by the head-tail symmetry of their particles in polar or nematic.For active polar particles these models have been formulated on a sphere 19 and on ellipsoidal surfaces 20 .
In these situations a robust rotating-band structure around the waist, with two +1 defects at the poles is found on a sphere.On an ellipsoid the location of the defects is linked to local geometric properties, similar to vortices in surface fluids [21][22][23][24] .The defects are related to the Gaussian curvature and to the umbilical points of the surface (see Materials and Methods for a geometric description).For spheroidal ellipsoids there are two umbilical points, which locate the two +1 defects.This configuration is more stable for prolate spheroids, where the umbilical points are at the points of maximal Gaussian curvature at the poles and less stable for oblate spheroids, where the umbilical points and the maximum in Gaussian curvature are separated.As in the spherical case a rotating-band structure is formed, with possible subbands which counter rotate depending on the initial condition.New dynamical features are found for non-spherical ellipsoids.They have four umbilical points.For lower velocities the defects encircle pairs of umbilical points and for larger velocities the defects are found at the high Gaussian curvature regions between each pair of umbilical points.With this richness in dynamics found for active polar particles on non-constant Gaussian curvature surfaces, we expect similar behavior for active nematic particles and ask up to which complexity of the geometry the dynamics of the four 1/2 disclinations can be effectively described by self-propelled polar particles.
To answer these questions, we first analyze the spherical case in more detail.In addition to the oscillation between the planar and tetrahedral defect configuration on a spherical vesicles and a tunable frequency by the activity and self-propulsion velocity we also track the positions of the defects.Computing the power spectrum from the time series for the average angle < α >= 1 showing the planar and tetrahedral defect configuration within a simulation of 1.000 particles (the four 1/2 disclinations are highlighted, the director field is shown -black lines -and the color coding corresponds to the nematic order parameter P , with minima in the four defects).The results are in excellent agreement with the experimental results in 6 .A video is provided in the SI.
depend on the activity.The same results, but with a small offset and a different slope are obtained for the coarse-grained description by self-propelled particles, see Figure 2. As a consequence for each activity in the nematic film a self-propulsion velocity can be determined in the coarse-grained description, which resamples the frequency of the planar-tetrahedral defect oscillation.Differences between both descriptions are found if we compare the trajectories of the defects and self-propelled particles.Within the considered time interval the 1/2 disclinations are locally confined, each defect only covers part of the vesicle.This is in contrast to the trajectories of the self-propelled particles, which rotate within a band We next consider spheroidal ellipsoids.They are characterized by the aspect ratio a/c and a = b, with a, b and c the length of the major axis.Due to the symmetry all geometric properties can be characterized with respect to the polar axis.As the geometry is topologically equivalent to a sphere we expect for passive systems again a minimal energy configuration with four 1/2 disclinations.They still try to maximize their distance, but are now also influenced by local geometric properties.The 1/2 disclinations tend to accumulate in regions of high Gaussian curvature 25,26 .Computer simulations for thin nematic shells have shown that for prolate ellipsoids pairs of defects are located at opposite ends close to the poles.The defects in each pair arrange at opposite sides of the surface and tend to align perpendicular to the pair at the other pole 25 .As the distance between the defects is no longer maximized, the geometric effect seems to dominate the repulsion in this case.
For oblate ellipsoids the 1/2 disclinations are found near the waist, where the Gaussian curvature is largest.Again two pairs of defects are found, one on each side.They repel each other and are mutually perpendicular to the other pair, leading to an alternating ring of 1/2 disclinations, one above and one below the waist.This behavior seems to be independent of the film thickness 25 , we have confirmed this behavior by our surface model without activity.
For active systems we observe again oscillatory behavior, see Figure 3.For prolate spheroids (a/c < 1) only two 1/2 disclinations are located at the poles, whereas the other two oscillate around the waist.The oscillations are very noisy and can not be tuned by the activity.Even if the distance between the two 1/2 disclinations at the waist is not optimal the average distance between all four defects is larger than in the passive case.While the 1/2 disclinations are still attracted by the high curvature regions at the poles, the active forces push one of the defects away leading to the observed metastable configuration.Within a transition zone (a/c ≈ 1) we observe similar behavior as in the spherical case (a/c = 1) without any defect localization.The behavior changes for oblate spheroids (a/c > 1), where all four 1/2 disclinations are along the waist, maintaining a maximal distance to each other.This behavior is similar to the passive system.However, the defects now oscillate between both sides.The frequency of the alternating oscillations above and below the waist can be extracted for various activities.However, a clear functional dependency on the activity could not be found.If the aspect ratio is further increased the situation changes to pairs of 1/2 disclinations which rotate around the umbilical points at the poles.The defects are no longer located at positions of maximal Gaussian curvature.The high curvature value at the waist creates a distortion of the nematic film, which can be seen from the nematic order parameter.It somehow serves as a barrier for the 1/2 disclinations preventing them from crossing the waist.The rotation is a consequence of the activity and the unfavorable short distance with respect to each other.The frequency of the rotation depends on the activity and can be tuned, see Figure 4. Also the transition to this rotating state depends on the strength of the activity.As stronger the activity as longer it is possible for the defects to cross the barrier at the waist.A tendency to locate the defects away from the high Gaussian curvature waist can also be seen for the passive case.The four different regimes are shown in Figure 5. using the order parameter with N the number of particles, [ts, te] an appropriate time interval and h i (t) the height of the defect i along the polar axis with respect to the waist at time t.We have η = 1 if all defects are at the poles, η = 0 if they are at the waist and η = 0.5 if they are homogeneously distributed along the polar axis.
Within the coarse-grained description by self-propelled polar particles, using the corresponding self-propulsion velocity according to Figure 2, we obtain a qualitatively different behavior.Within the considered parameter regime, the values for η are independent of the self-propulsion velocity.For aspect ratios a/c < 0.5 the particles rotate on closed trajectories, well separated from each other at approximately equal distance along the polar axis.
The transition zone with sphere-like behavior is more extended than for the nematic de- propelled particles are independent of the activity in the corresponding regime to the considered velocities v 0 .From left to right we have (blue) the situation for prolate shapes with location of two defects at the poles, leading to η > 0.5, (green) spherical like shapes with no clear location of the defects, leading to η ≈ 0.5, (yellow) oblate shapes with location of the defects along the waist, leading to η < 0.5, for larger a/c we obtain a phase transition towards the rotating state, with the defects located around the poles, leading to η > 0.5.The transition towards this state depends on the activity.fects.For 0.5 < a/c < 2 a band structure is formed around the waist, which shrinks with increasing aspect ratio.For a/c > 2 all particles are positioned at the waist, rotating in one direction and maintaining their distance.The regime with pairwise rotating defects around the umbilical points could not be found within the coarse-grained model (A more detailed description is given in the SI).Non-spherical ellipsoids, which are characterized by a = b, a = c and b = c, have four umbilical points.They are either prolate-like or oblate-like but in any case have two distinct points of maximal Gaussian curvature.We thus analyze the distance of the four 1/2 disclinations with respect to the umbilical points and the points of maximal Gaussian curvature using the average geodesic distances < DDU > and < DDG >, respectively.Figure 6, which is inspired by 20 shows the distances as a function of the aspect ratios a/b and a/c.Spheroids are also included, the first column shows the previous results for oblate and the diagonal for prolate geometries.Each row in between thus corresponds to a transition from oblate-like to prolate-like geometries.In most cases the 1/2 disclinations are closer to the high Gaussian curvature points than to the umbilical points, with the only exception for oblate-like ellipsoids with a large aspect ration a/c ≥ 4.This leads to the conclusion that 1/2 disclinations tend to be attracted by points of high Gaussian curvature.

DISCUSSIONS
In 6 it was shown that in a confined active system, a dense suspension of microtubules and molecular motors on the surface of a spherical lipid vesicle, cyclic oscillations between defect configurations can be observed.They result from topological constraints and the coupling between velocity fields and defect-defect interactions.These findings may push forward the design of systems that harness the ability of nanoscale active matter to transform chemical energy into mechanical work.On non-spherical surfaces defects are known to be strongly influenced by local geometric properties.The induced geometric interaction can lead to locating of defects, which is established for vortices in surface fluids 21,22,24 and vortices, sources and sinks in polar systems 20,27 .For strong variations in geometric properties it has even be found computationally that lower energy minima in passive systems can be formed by creating additional defects 23,27 for surface fluids and surface polar particles, respectively.Our work extends the understanding of the delicate relations between topology, geometry and defect dynamics on non-spherical shapes for the system considered in 6 .We are concerned with ellipsoidal surfaces and identify crucial geometric features which influence collective motion patterns in active nematic films.We have shown that 1/2 disclinations are related to both, maxima in the Gaussian curvature and umbilical points of the surface.On prolate spheroids maxima in Gaussian curvature and umbilical points coincide, they are located at the two poles and attract the 1/2 disclinations.However, the repulsive defect-defect interaction allows only two of the defects to be located at the poles, the other two try to maximize their distance and are located around the waist, where they oscillate.Spherical like shapes lead to similar behavior as observed on a sphere, with no distinguished location of the defects and an oscillation between a tetrahedral and planar defect configuration.For oblate spheroids all 1/2 disclinations are located at the waist, the region of high Gaussian curvature.They again maximize their distance and oscillate.With increasing aspect ratio a/c the situation changes.The defects can no longer cross the waist, where the high Gaussian curvature leads to a distortion of the nematic order.As a consequence pairs of 1/2 disclinations rotate around the umbilical points.The frequency of the rotation depends on the activity and can be tuned.This found rotating state is an other step towards a controllable transformation of chemical into mechanical energy in nanoscale active matter and asks for experimental validation.The results for non-spheroidal ellipsoids confirm these findings, even if the separation of the different states is not as distinct as in Figure 5.A smooth transition of the dynamics between prolate-like and oblate-like shapes is identified in Figure, with a clear tendency of the 1/2 disclinations to locate at points of maximal Gaussian curvature.Only for extrem values of a/c and almost spheroidal shapes the situation changes and the rotating state around the umbilical points could be identified.
We further demonstrate that the proposed coarse-grained description of 1/2 disclinations in active nematic matter by self-propelled particles fails if geometric properties come into play.Already on spherical shapes the trajectories of the defects and the self-propelled particles differ significantly and on spheroidal ellipsoids both descriptions not even qualitatively agree.
In summary we explored the complex interaction of topology, geometry and defect dynamics in nematic films on ellipsoidal surfaces and demonstrated how topological constraints and geometric properties can be used to control the collective behavior in nanoscale active matter.The non-linear coupling between non-constant Gaussian curvature and defect-defect interactions leads to tunable spatiotemporal patterns.Among these findings is a stable rotating state on strongly oblate-like ellipsoids, which suggests an other pathway towards a controllable generation of mechanical work in nanoscale active matter.The richness of physics observed in our work will further increase if the underlying shape is deformable.
First experimental results of such an interplay between activity-driven defect motion and deformability of the vesicle are already shown in 6 and discussed in 28 .However, for theoretical descriptions of these phenomena new methods will be required.

MATERIALS AND METHODS
We consider a more general approach than the Viscek-like models confined on a sphere or an ellipsoids in 19,20 .

A. Equations of motion
We consider N active particles of mass m i = 1, which are constrained to move on a surface algebraically described by g(q) = 0, with particle positions q = (q 1 , . . ., q N ).Newton's equations of motion (EOM) with holomonic constraint g(q) read: with forces F = (F 1 , . . ., F N ) and velocities v = (v 1 , . . ., v N ).λ = (λ 1 , . . ., λ N ) are the Lagrange multipliers and G(q) = ∇ q g(q) is the Jacobian of g(q).The force F i can be written as: where γ is the translational friction coefficient, F ac i the active force acting on particle i and F ij the pair-interaction force between particle i and particle j.Additionally every particle has an internal degree of freedom, its orientation n i .Denoting by ω i the angular velocity we have the following EOM for the orientational dynamics: where γ a is the rotational friction coefficient and T i (q, n) is the torque acting on particle i, with n = (n 1 , . . ., n N ).Depending on the specific form for the active force F ac i , the pairinteraction force F ij , the torque T i (q, n) and the holomonic constraint g(q) we will be able to describe polar and nematic active systems on various surfaces.

B. Active polar particles
For active polar particles on a sphere of radius R we would specify 19 F ac i = v 0 n i with a constant self-propulsion velocity v 0 , F ij = k(2σ − q g ij ) q i −q j q ij for q g ij < 2σ and F ij = 0 otherwise, a short-range repulsion between spheres of radius σ, with elastic constant k, euclidian distance q ij = |q i − q j | and geodesic distance q g ij = |q i − q j | g .Parallel orientations between neighbouring particles are favored and therefore we use the aligning torque T i (q, n) = −J j∈U (i) (n i × n j ), with J > 0 the strength and U (i) the first shell of neighbors of particle i, identified as all the particles within a cutoff radius of 2.4σ from r i .The holomonic constraint for a sphere of radius R reads g(q i ) = q 2 i,1 + q 2 i,2 + q 2 i,3 − R 2 , with q i = (q i,1 , q i,2 , q i,3 ) ∈ R 3 .This approach can be used to reproduce the results in 19 in which the overdamped limit, the euclidian distance instead of the geodesic distance and an additional noise term are considered.

C. Active nematic particles
To describe active nematic particles we use the tensor order parameter Q j αβ = n j α n j β − δ αβ /3 , where the upper index corresponds to the particles and the lower indices represent the componenets x, y, z.The active force does not distinguish 'head from tail' and it thus has the form: The torque reflects the fact that both parallel and anti-parallel configurations are favored.
It has the form: The pair-interaction force F ij and the holomonic constraint g(q) are the same as in the active polar particles case.The simulation parameters for this case are (J, k, σ, γ, γ a , v 0 ) = (10, 3, 2, 0.1, 2.5, 1.1) unless otherwise specified.

D. Coarse-grained defect description
In the coarse-grained defects description by active polar particles 5 and 6 the elastic energy between defects is E ∼ log(q g ij ), where q g ij is the geodesic distance between the defects.The pair-interaction force is therefore F ij = k q g ij q i −q j q ij , which is no longer short-ranged.Defects align anti-parallel to each other and the restoring torque strength is 6 T i = J j∈U (i) cot( , where θ ij is the angle between n i and n j .The vector form for the torque can be written in terms of the orientations as: Finally the defects are treated as self-propelled particles and the active force is F ac i = v 0 n i .The simulation parameters for this case are (J, k, γ, γ a , v 0 ) = (3, 4, 0.1, 2.5, 0.11) unless otherwise specified.

E. Geometric properties
Besides a sphere we consider two classes of ellipsoidal surfaces: For spheroidal ellipsoids two of these values are equal.The algebraic description reads g(q i ) = i,1 In figure 7 we show three different ellipsoids, where umbilical points are highlighted and the color coding corresponds to the Gaussian curvature K.

F. Numerical methods
Eqs. 1 have been numerically solved using RATTLE discretization 29 .The equation for the orientational dynamic eqs. 3 have been first solved unconstrained with the torque T i projected on to the normal plane of the surface at point r i .Afterwards the orientation n i has been projected on to the tangent plane of the surface at point r i and the angular velocity ω i takes the direction of the normal to the surface at point r i .
We fix the number of particles N = 1000 and the volume fraction φ 1.The surface area for the sphere is equal to A = 4πR 2 , with R = 31.6.The ellipsoid parameters a, b, c have been chosen such that the surface area is equivalent to the surface area of the sphere and the aspect ratio is respected.
The nematic order parameter is defined as where the sum is over the nearest neighbors and w ij = q −1 ij .
Defects are calculated as the local center of mass for regions where the local order parameter P i is smaller then 0.45 (some corrections were required for regions of high Gaussian curvature, due to strong distortion of the director field).
The simulation code is implemented in C++, using the GeographicLib library 30 for the calculation of the geodesic distances.However, for non-spheroidal ellipsoids the euclidian distance has been used.This approximation can be justified by the short-range interactions.
Data have been analyzed using Python, Ovito 31 and Paraview.

FIG. 1 .
FIG. 1. Defect oscillations: a) Top: Kymograph showing the time evolution of the angles α ij , which denote the angle between the radii from the center of the sphere to each of the defect pairs.Bottom: Oscillation of the average angle < α >.The blue and the green line correspond to the planar (< α >= 120 • ) and tetrahedral (< α >= 109, 5 • ) defect configuration.b) Snapshots

FIG. 3 .
FIG. 3. Defect localization on spheroids: a) Snapshot showing the defect configuration within a simulation of 1.000 particles on a prolate spheroid with a/c = 0.25 (the four 1/2 disclinations are highlighted, the director field -black lines -is shown and the color coding corresponds to the nematic order parameter P , with minima in the defects).In addition the trajectories of the four 1/2 disclinations are shown (each color corresponds to one defect).The height h i for each defect with respect to the waist is also shown as a function of time.b) same as a) for a oblate spheroid with a/c = 2.The oscillations of the four defects have the same frequency and alternate with respect to each other.Videos for a) and b) are provided in the SI.

FIG. 4 .
FIG. 4. Defect rotations: a) Snapshots from above and below showing the defect configurtion within a simulation with 1.000 particles on an oblate spheroid with a/c = 6 (the four 1/2 disclinations are highlighted, the director field -black lines -is shown and the color coding corresponds to the nematic order parameter P , with minima in the defects).b) Oscillations of the angle measuring the rotation around the umbilical points (top and bottom) and c) frequency of the oscillation as a function of activity for two different aspect ratios.A video for case a/c = 6 is provided in the SI.
FIG. 5. Phase diagram: Phase diagram for patterns and oscillations on spheroidal ellipsoids for 1/2 disclinations and self-propelled particles.The results for the coarse-grained description by self-

FIG. 6 .
FIG. 6. Relation to geometric properties: Average geodesic distance of 1/2 disclinations to the umbilical points < DDU > (left) and to the points of maximal Gaussian curvature < DDG > (right) for non-spheroidal and spheroidal (first column -oblate and diagonal -prolate) ellipsoids of different aspect ratio.Only for the extreme case of a/c = 4,6 and a/b = 1.1 the disclinations are closer to the umbilical points.Also in these cases a rotating state as in Figure 4 can be observed, which however is not as regular, see SI.In all other situations the disclinations are closer to the points of maximal Gaussian curvature.A video for case a/b = 1.1 and a/c = 6 is provided in the SI.

FIG. 7 .
FIG. 7. Geometric features: Example of an ellipsoid with major axis a) a/b = 1 and a/c = 4, b) a/b = 1.25 and a/c = 4 and c) a/b = 1 and a/c = 0.25 .Umbilical points are shown as points and the Gaussian curvature K is color coded.

c 2 − 2 T
1 = 0.An umbilic point is a point where the maximum and minimum curvatures coincide.At an umbilical point, the surface is "locally spherical".These points are found at±a a 2 − b 2 a 2 − c 2 , 0, ±c b 2 − c 2 a 2 − c