Vortex formation and dynamics of defects in shells of active nematics

We present a hydrodynamic model for a thin spherical shell of active nematic liquid crystal with an arbitrary configuration of defects. The active flows generated by defects in the director lead to the formation of stable vortices, analogous to those seen in confined systems in flat geometries, which generate an effective dynamics for four +1/2 defects that reproduces the tetrahedral to planar oscillations observed in experiments. As the activity is increased and the vortices become stronger, the defects are drawn more tightly into pairs, rotating about antipodal points. We extend this situation to also describe the dynamics of other configurations of defects. For example, two +1 defects are found to attract or repel according to the local geometric character of the director field around them, while additional pairs of opposite charge defects can give rise to flow states containing more than two vortices. Finally, we describe the generic relationship between defects in the orientation and singular points of the flow, and suggest implications for the three-dimensional nature of the flow and deformation in the shape of the shell.


I. INTRODUCTION
Active liquid crystals [1][2][3] (ALCs) have proved successful as a paradigm for living systems on the microscale, providing insight into processes like cell motility [4][5][6] and division [7][8][9], development of cell shapes [10,11], and growth of cell colonies [12]. Certain fundamental motifs have been developed such as the instability of uniformly aligned states, the emergence of spontaneous flows, the creation and self-propulsion of topological defects and the shear-thinning character of extensile gels. When ALCs are confined to a circular geometry a prominent feature is the emergence of stable flow vortices. Confinement gives rise to a single vortex state in dense bacterial suspensions [13,14], active nematic suspensions [15,16] and monolayers of migrating cells [17,18]. Circulatory flows are also characteristic of cytoplasmic streaming [6,16,19,20]. When the system size is increased such vortices become unstable [21] and turbulent flows develop, a prevalent feature in bulk active fluids [22][23][24]. In active systems with high frictional dissipation stable vortices can also arise in the absence of spatial confinement [25][26][27]. Recent experiments by Keber et al. [28] are realisations of a different type of confined geometry, in which the ALC adheres to the surface of a vesicle. Four half-integer defects form in the orientation of the microtubule-based extensile active nematic and are found to be in steady motion, oscillating between tetrahedral and planar configurations. A variety of other states is observed, like defect-associated membrane protrusions, two vortex defects in smaller spherical vesicles and two aster defects in spindle-like vesicles with stiffer microtubules. Here, we develop an active hydrodynamic model for an ALC confined to a spherical shell and show that the dynamics of this system is also characterised by the formation of vortices, which reproduces the defect motion from experiments.
Defects in the director are unavoidable on the sphere [29]. In a typical situation there are four, all of strength +1/2, which are known to self-propel in active liquid crystals [30]. This motivates a minimal description of their motion as a point particle dynamics, and such a model was shown to reproduce the main experimental observations [28]. We extend this to a hydrodynamic model, in the confined geometry of a spherical shell, and show that the dynamics is characterised by the formation of two stable counterrotating vortices, one in each hemisphere, paralleling the vortex formation seen in other types of confinement [13,16,31]. A minimal hydrodynamic model takes the positions of defects to construct a profile for the director over the entire sphere, whose associated active flows advect the defects to yield a selfconsistent dynamics. Tetrahedral to planar oscillations of four +1/2 defects are also obtained with this model. As the activity is increased the two vortices become more pronounced and the pair of +1/2 defects within each are pulled closer together in an effective attraction of like-charge defects. These oscillations appear at a finite threshold of the activity, below which the defects form static configurations of distorted tetrahedra. Linear stability analysis captures the mode of deformation and the threshold for defect motion.
Just as there are defects in the director field there are also vortices and stagnation points in the flow field [29]. There is a one-way relationship that assigns to a defect in the orientation a flow singularity whose winding number depends only on the defect's topological strength. The oscillations of four half-defects are found to be stable against additional half-defect pairs created randomly in larger shells. If the defects are instead induced at specific positions, it is possible to generate more complex, metastable flow vortex configurations. The dynamics of polar configurations with only integer strength defects is similar and we find attraction of pairs of aster-like +1 defects in extensile active nematic shells, but repulsion for vortex-like defects. The speed of defects in the polar case is shown to have different scaling than for nematic shells, in particular the type of motion does not depend on the radius in the former case whereas it does in the latter.

II. MODEL
We consider an active nematic in a thin spherical shell of thickness h 0 and inner radius R, with h 0 /R 1. The three-dimensional flow u = (u r , u ⊥ ) in the shell is driven by gradients in the active stresses and can be found as the solution of the generalised Stokes and continuity equations, −∇p + µ∆u + ∇ · σ = 0 and ∇ · u = 0, where p is the pressure and µ the viscosity. The active stress σ a = −σ 0 P P − 1 3 I is extensile throughout this paper, σ 0 > 0, in order to relate with microtubule-based active nematics [22,28], although we comment on the contractile case at the end. If the polarisation P is specified one can solve for the active flow generated by it in a thin film approach [32][33][34][35][36], decribed in Appendix V A. We take the polarisation to be tangential throughout the shell thickness, P = cos(ψ)ê θ + sin(ψ)ê φ , and construct an explicit form from the positions of the defects. This can be done using stereographic projection from the complex plane, z(θ, φ) = R cot(θ/2)e iφ . In the plane a nematic director n = (cos α, sin α) with n def defects with topological strengths m j and positions z j = x j + iy j is given by α = α 0 + j Im (ln(z − z j ) mj ) [37], where the phase α 0 ∈ [0, π) parameterises whether the local geometry of the director around a defect is more splay-like or more bend-like. Finally, stereographic projection of n onto the sphere yields a polarisation field in the spherical shell via Parametrised in this way, P is an exact minimiser of the elastic energy of a nematic on a sphere in the one-elasticconstant approximation [38]. Moreover, it consists only of those defects from which α is constructed explicitly, provided j m j = 2.
In dimensionless variables the tangential component of the flow then has the form with the radial profile f (r) =r 2 2 −r, wherer ∈ [0, 1] is the radial position within the shell in units of h 0 . This solution corresponds to a no-slip inner surface and a vanishing tangential stress on the outer surface. Figure 1 gives an example of the director for four +1/2 defects in a planar configuration and the corresponding active flow given by Eq. (2), which is seen to consist of two counterrotating vortices. This emergence of stable vortices is the germane feature of the active flows on spherical shells.
The director dynamics is dominated by the motion of defects, when the orientational dynamics is rapid. In our approach the director is instantaneously given by the parameterisation above as the defect positions change. The defects are advected by the flow they create and we describe their motion in a point-particle description [30,39]. Each defect moves due to the tangential component of the active flows, given by Eq. (2), and due to standard nematic elasticity. The overdamped dynamical system for the defect positions r k (t) is The resultant dynamics is similar to [28] except that here we obtain the advective flow u def k from a self-consistent hydrodynamics in the spherical shell and generalise to an arbitrary collection of defects. The flow is divergent at the defect locations, therefore we introduce a cut-off and obtain the defect velocity as an average of the flow over a small circle γ k (s) centered at the defect For the defect motion the flow is evaluated at the outer surface, where f (r = 1) = −1/2. The circle γ k (s) = (θ k + ρ cos(s), φ k + ρ sin(s)/ sin(θ k )), s ∈ [0, 2π], has the opening angle ρ, which can be associated with the core size r c of the defect through the relation The core size could be measured for a particular experimental system, for instance as the size of the region around a defect which is devoid of active nematogens. The elastic force F k (t) provides attraction or repulsion of defects depending on their topological strength, with an effective friction coefficient ξ and elastic constant K takes the form This choice of time scale sets the scale of the elastic terms toK = τ K/ξR 2 = 1. This identifies the scaling of τ u def /R as the defining parameter for the defect dynamics, which represents the ratio of active to elastic effects and differs depending on the topological strength of the defect. Equations (6) and (7) are integrated numerically for different defect configurations using a standard Runge-Kutta method.

A. Active flow at the defects
In addition to the singularities in the director, the vortices in Fig. 1 (b) contain singularities in the flow field, about which the flow circulates. Such flow singularities are topologically required [29] and can be generated at the locations of defects in the director. A general relationship between defects and flow singularities follows from evaluating (2) on the small circle γ k (s) and expanding in powers of ρ, the angular distance to the k-th defect (see Appendix V A). We make use of the stereographic projection to writẽ The dominant contribution to the flowũ ⊥ near the k-th defect diverges as ∼ 1/ρ and has the winding number Unit strength defects produce a vortex-like (I = 1) singularity in the flow, whose character is sink-or source-like according to whether the defect resembles an aster or a vortex, respectively. When there are two such defects, at antipodal positions, they generate two counterrotating vortices with no other flow singularities. On the other hand, simple stagnation points (I = −1) cannot be created at defect locations. For half-integer defects relation (9) was shown in [30,39] and the flow around a single spiral defect in active polar gels was studied in [41].
In a typical situation the flow singularities at defects are not sufficient to generate a total winding of +2. This is most evident for four half-defects, as seen in Fig. 1 where flow vortices form inbetween the defects, because for m k = 1/2 the flow is non-winding (I = 0). Instead, it is directed along the defect's symmetry axis For these defects, we approximate the advective flow u def in (4) by this well-defined flow direction and the magnitude where U 0 = h 2 0 σ 0 /Rµ is the typical active flow magnitude in the thin film approach (see equation (21) in Appendix V A) and we replaced ρ = r c /R. The speed of +1/2 defects does not depend on the shell radius R, because they generate their own advection locally, where the defining length scales are the core size r c and the shell thickness h 0 . In equations (6) and (7) the scaling of the dimensionless advective term for a +1/2 defect is The next term in the expansion (8), which is O(1), is non-winding only for m k = 1 (see equation (26) in Appendix V A). Therefore, unit strength defects are advected with a flow ∼ U 0 , and the relevant parameter becomes This predicts a different scaling of the defect dynamics in thin polar shells compared to nematic shells. In the former only integer strength defects are present and, notably, the type of motion does not depend on the radius. For all other defect types active advection scales at most as ∼ U 0 r c /R, which makes it negligible compared to the active motion of +1/2 defects. In particular, −1/2 defects can be approximated with u def = 0 in a collection of ±1/2 defects.

B. Four +1/2 defects
In the minimal case of four +1/2 defects, the dynamics is determined by the parameter ν, defined in (12). We increase ν through the activity σ 0 , keeping all other parameters constant, in particular the radius, in order to fix the time scale τ . The phase α 0 also affects how the defects move. In the ranges (0, π/4) and (π/4, π/2) the dynamics is similar and we choose α 0 = π/2 − 0.  For intermediate activity the positions of the four defects periodically pass through tetrahedral and planar configurations, as shown in Fig. 2 (a-c), which is the dynamics found in experiments [28]. The motion is characterised by the formation of two counterrotating flow vortices that separate the defects into two pairs, in which they rotate around each other. This effect becomes more pronounced as the activity is increased, as shown in Fig.  2 (d-f). The separation of defects within each pair decreases significantly with ν. There is also a gradual change in the shape of the trajectories, from square-like to more ellipsoid, such that the tetrahedral configurations are no longer approached and the defects oscillate between two different planar arrangements. As the defects in each pair are drawn closer with increasing activity the dynamics approaches the situation for two antipodal spirals in the director, which in the limit generate a perfectly symmetric flow vortex pair. This behaviour is summarised in Figure 3, where the mean and the minimal angular distances are plotted against ν, the latter reflecting the decreasing separation between defects in each pair.
The total speed of the +1/2 defects, which also includes motion due to elasticity, is dominated by their active speed v 0 given by (11). The frequency of the defect   (36) and (37) (lines). The data suggests that at ν * = 0.7 the eigenvalue λ7 (boxed) becomes positive, which renders the "skewed tetrahedron" linearly unstable. This is qualitatively confirmed by the theoretical prediction, albeit with an overestimated transition point.
oscillations is thus without accounting for the small changes in the orbit shape with increasing ν.
The effective attraction of defects into pairs is mediated by the active flow vortices that form inbetween them, which in turn are controlled by the underlying director. In the tetrahedral configuration the nematic has a characteristic tennis ball texture [40,42], but for generic values of α 0 this texture is skewed, such that each two defects form a separated spiral. The flow vortices accordingly acquire a sink-or source-like component, depending on the tilt in the spiral. As can be seen for instance in Figure 2 (b), the paired defects have a sink-like vortex inbetween them which keeps them together. This active attraction mechanism requires the possibility of radial flows to accomodate this influx, which is guaranteed in the thin film approach.
The choices α 0 = 0, π/2 produce zero tilt in the texture of the initial tetrahedron and the resulting dynamics lacks the contraction of defect trajectories in one of the directions, such that they continue passing through tetrahedra for high activity. Finally, α 0 = π/4 does not have a dynamical regime and defects relax into increasingly tight, but stationary pairs.

C. Linear stability of static configuration
The defects only move above a critical ν * ≈ 0.7 (Fig.  3) and we describe this transition in a linear stability analysis. The system is initialised with the four defects at the vertices of a tetrahedron with θ (0) i = (β, π − β, β, π − β) and φ (0) i = (0, π/2, π, 3π/2), with β = arctan √ 2 . It is evident from simulations that for activities below the threshold the defects settle into an increasingly distorted tetrahedron, which can be described by the coordinates with small deviations δθ and δφ as shown in Fig. 4 (a). Using this ansatz we find analytical solutions for the deviations at linear order (see Appendix V C). The deformation of the tetrahedron is a superposition of two modes -twisting around and stretching along the z-axis, illustrated in the inset of Fig. 4 (a).
At the critical activity the skewed tetrahedron becomes linearly unstable, as seen from the spectrum of the dynamical matrix ∇g (see Appendix V C for definition) shown in Figure 4 (b). The simulation data suggests that one eigenvalue changes sign at ν * , while all others stay non-positive, indicating that the skewed tetrahedron is stable below ν * . The three vanishing eigenvalues correspond to rigid body rotations. Calculating the eigensystem using the analytical solutions for the deviations above the threshold allows to characterise this instability, albeit with an overestimated transition activity. The eigenvalue λ 7 becomes positive at ν ≈ 1.0, which marks the linear instability of the skewed tetrahedron towards a deformation that strongly increases the twist and slightly reverses the stretching. This can be seen from the corresponding eigenvector, which is of the form (a, −a, a, −a, −b, b, −b, b) with b a > 0, and is exactly the dynamics found in simulations at the beginning of the periodic orbits (see Fig. 3 (b) and (c)). (a) Distance of the two +1 defects over time for different local director geometry, controlled byα0, which varies in steps of π/16. Here, the motion of defects is due to active advection only, withK = 0. (b) Perfect asters (α0 = 0) or vortex defects (α0 = π/2) move along geodesics, but in general the defects move on outward or inward spiralling trajectories. Two perfect spirals (α0 = π/4) move along a circular path, without changing their distance. (c) Two spiral-like defects (α0 = 5π/16) are attracted to each other by the flow vortex that forms inbetween them. (d) When elasticity is included, attraction of +1 defects is found only for ν (1) above a threshold and for large enough tiltα0.

D. Two +1 defects
ALCs can develop unit strength defects in their orientation [41,43] and the confinement to a spherical shell provides a setup where two such defects could be topologically stabilised. We study the motion of two +1 defects in the limit of very strong activity, settingK = 0 and using the active time scale T = R 2 µ/h 2 0 σ 0 . The value of α 0 that is required to generate a particular director geometry at the defects depends on their position. This ambiguity can be removed by setting α 0 = − arg(z 1 −z 2 )+α 0 , where the additional constant is found by imposing an aster-geometry forα 0 = 0 for both defects irrespective of their position. Now,α 0 ∈ [0, π/2] produces increasingly tilted spirals, up to two pure vortex defects forα 0 = π/2. We find that two +1 defects are either attracted to or repelled from each other by active advection, depending on the local director geometry, as shown in Fig. 5 (a). Two defects that are aster-like (0 ≤α 0 < π/4) experience repulsion and relax into an antipodal configuration. Since the asters create sinks in the flow, a source-like flow vortex forms in between them, pushing them apart. Vortex-like defects (π/4 <α 0 ≤ π/2) show the converse effect and are drawn towards each other. In this idealised setting without elasticity, they merge into a +2 boojum, with a local flow structure of a +3 singularity accompanied by a stagnation point at the antipodal point. Two perfect spiral defects (α 0 = π/4) keep a constant distance, rotating around each other on a circular path.
During this motion the defects are typically spiralling inward or outward, as shown in Fig. 5 (b). Only in the two limiting cases do the defects move along their connecting geodesic. Figure 5 (c) shows the active flow with the additional sink-like vortex inbetween the defects (α 0 = 5π/16), that draws the defects inward on a spiralling trajectory. The defects' trajectory rotates in a direction opposite to the rotation of their local flow vortices.
When elastic repulsion is included, withK = 1, the defects relax into the antipodal configuration for allα 0 for activities up to ν (1) ≈ 12. Above this threshold active attraction overbalances the elastic repulsion for large enoughα 0 , as shown in Fig. 5 (d). Interestingly, in such cases the defects again collapse into a very tight pair. In an experimental system, fluctuations in the tilt of a spiral around the limiting value ofα 0 might lead to oscillations between the antipodal and the collapsed configurations.

E. Many-defect states
When the activity |σ 0 | or the shell size R are increased additional ±1/2 defect pairs may be created on top of the four +1/2 defects, as the system approaches the onset of active turbulence [44]. To study such situations we increase the radius as R = γR 0 , with γ > 1, keeping all other parameters fixed. This changes the elastic time scale to τ = γ 2 τ 0 and the activity-to-elasticity ratio to ν = γν 0 . The reference values correspond to parameters in Section III. A for the regime of tetrahedral-planar oscillations, for instance ν 0 = 1. We consider a system with four defects in this oscillatory state and inject one ±1/2 pair at a random position. Figure 6 (a) shows how the dynamics reacts to this perturbation. One of the +1/2 defects very quickly annihilates with the −1/2 and the remaining four defects resume the oscillation, usually in a different pairing. Similarly, when all defects are placed at random positions the annihilation events happen rapidly, leaving the minimal four-defect state in the oscillating regime. The same is found for more than one additional pair of half-defects in the system. This indicates that the oscillatory state is stable, as long as additional defect pairs occur as fluctuations and are not produced constantly. On the other hand, by inducing additional ±1/2 defects at specific locations more complex flow vortex configurations may be constructed, in which elastic forces and active flows are balanced. The simplest many-defect configuration that is metastable has n pair = 2 additional defect pairs, shown in Fig. 6 (b,c). The six +1/2 defects are allocated to three flow vortices arranged equidistantly around the equator, with another three vortices rotating in the opposite direction inbetween them. The three-fold symmetry of this flow field is guided by the flow singularities at the −1/2 defects, which have I = −2 and are located at the poles. This configuration is transient and reduces to the four-defect state due to coalescence of oppositely charged defects. In Figure 6 (d) the times to the first annihilation event for this metastable configuration and for n pair = 2 with random initial defect positions are compared. For the vortex configuration this time increases considerably with γ. This opens an interesting direction of tuning specific many-defect states before the onset of active turbulence by exploiting the topologically required singularities in both flow and director. The metastability of such configurations could be aided by an advantageous manipulation of the shell shape, for instance by trapping positive defects in regions of higher curvature [45].

IV. DISCUSSION
In contact with a passive fluid on the outside the active shell may swim due to its self-generated surface flows. In a squirmer approach the swimming velocity of a spherical shell can be calculated as the surface integral of the slip velocity, U (t) = − 1 4πR 2 S u ⊥ (t)dS, with a similar expression for the angular velocity Ω(t) [46]. Evaluating the integrals exactly, we find that a single +2 defect as well as two +1 defects do not generate translational or rotational motion of the shell, irrespective of their position. The same results from numerical integra-tion of the slip velocity for four half-defects along their symmetric trajectories. There is a number of ways in which this symmetry could be broken in order to achieve self-propulsion and rotation, including distortions of the vortices by noise, changes in vesicle shape [47], and interactions with surfaces or other active vesicles [22]. One promising direction is a controlled asymmetric modification of the local director geometry around two +1 defects [48].
Our results can be extended to contractile active fluids by changing the sign of the activity σ 0 . The reversed sign of the flow exchanges the role of splay-like and bend-like distortions in the orientation. The direction of motion of half-integer defects is reversed, but the tetrahedralplanar oscillations and the formation of vortices -with opposite rotation sense -is unchanged. Similarly, the type of active interaction between +1 defects is reversed. The thin film approach used here allows for a non-zero radial flow component, which in general is present in the examples considered and enables the repulsion or attraction of defects due to active flows. The radial component is small, ∼ O(h 0 /R), compared to tangential flows and will result in a dynamic deviation from the spherical shape that complements the defect motion [48]. Stationary shell shapes, e.g. for two asters, should locally resemble profiles found for flat droplets of active nematics with defects [35,36]. We have taken a one elastic constant approximation for simplicity, but in systems of elongated filaments one can expect anisotropy. If this is sufficiently large it could lead to qualitative changes of the dynamics and is certainly an extension worth pursuing. Our work establishes the formation of vortices under confinement as a generic feature also for ALCs on spherical surfaces. It would be interesting to extend to other topologies, for instance tori with additional handles [49].
We thank Carl Whitfield, Daniel Pearce and Julia Yeomans for enligtening discussions on various aspects of this work. This work was supported by the UK EPSRC through Grant No. A.MACX.0002.

A. Flow in a thin active nematic shell
In a thin spherical film of active nematic the generalised Stokes equation given in the main text can be expanded in the small parameter ε = h 0 /R. Using the typical magnitude U 0 = R/T of active flows, with the time scale T , the dimensionless velocity components areũ θ = u θ /U 0 ,ũ φ = u φ /U 0 , andũr = ur/εU 0 . The radial coordinate becomes and the corresponding partial derivatives ∂rf = 1 h 0 ∂rf , for some function f (r). In the dimensionless form, for instance the θcomponent of the Stokes equation becomes and there is a similar expression for the φ-component. The active stress tensor gradient is Therefore, in order for the activity to drive the tangential flows the corresponding prefactor in (19) has to scale as ∼ O(1), leading to the dimensionless prefactorσ and a similar relation for the pressure. Analogous to planar thin films [32,36], the leading order part of the r-component of the Stokes equation yields a constant pressure. With the boundary conditions ∂r sũφ |r s=1 = 0 andũ φ |r s =0 = 0 , the solution (2) for the tangential flow components is obtained.
For a fixedr, for instancer = 1 used for the defect dynamics, the tangential flow can be written in a complex representation making use of the stereographic projection z(θ, φ) = R cot(θ/2)e iφ . In this way, the projection point is the north pole and the plane crosses the sphere along the equator. The complex flow is given bỹ This expression is well-defined through the stereographic projection of the tangential flow onto the plane, which is given in the complex form as u = ux + iuy and relates toũ as We evaluate (23) on the projection of the small circle γ k (s) introduced in the main text, which has the form with s ∈ [0, 2π], provided the circle does not enclose one of the poles on the sphere. An expansion of (23) in powers of ρ reads where the functions h 1 and h 2 only depend on the defect positions and other constants. Integrating this expression over s yields a non-winding contribution at the order O(1/ρ) for m k = 1/2 and at the order O(1) for m k = 1, as discussed in the main text. B. Point-particle-like dynamics of defects The free energy of a nematic on a sphere can be phrased in terms of the defects' pairwise interaction energies and self-energies [38,40,42], which are constant in our model, The angular distance between defect i and j is given by cos β ij = cos θ i cos θ j + sin θ i sin θ j cos(φ i − φ j ) .
For all defect configurations the dynamical systems are integrated using the ordinary differential equation solver ode23s provided by the software MATLAB 2016a, with relative and absolute accuracies set to 10 −6 τ .