Morphology of active deformable 3D droplets

We numerically investigate the morphology and disclination line dynamics of active nematic droplets in three dimensions. Although our model only incorporates the simplest possible form of achiral active stress, active nematic droplets display an unprecedented range of complex morphologies. For extensile activity finger-like protrusions grow at points where disclination lines intersect the droplet surface. For contractile activity, however, the activity field drives cup-shaped droplet invagination, run-and-tumble motion or the formation of surface wrinkles. This diversity of behaviour is explained in terms of an interplay between active anchoring, active flows and the dynamics of the motile dislocation lines. We discuss our findings in the light of biological processes such as morphogenesis, collective cancer invasion and the shape control of biomembranes, suggesting that some biological systems may share the same underlying mechanisms as active nematic droplets.

We numerically investigate the morphology and disclination line dynamics of active nematic droplets in three dimensions. Although our model only incorporates the simplest possible form of achiral active stress, active nematic droplets display an unprecedented range of complex morphologies. For extensile activity finger-like protrusions grow at points where disclination lines intersect the droplet surface. For contractile activity, however, the activity field drives cup-shaped droplet invagination, run-and-tumble motion or the formation of surface wrinkles. This diversity of behaviour is explained in terms of an interplay between active anchoring, active flows and the dynamics of the motile dislocation lines. We discuss our findings in the light of biological processes such as morphogenesis, collective cancer invasion and the shape control of biomembranes, suggesting that some biological systems may share the same underlying mechanisms as active nematic droplets.

I. INTRODUCTION
Active particles use energy from their surroundings to do work. Examples range from eukaryotic cells, bacterial suspensions and motor proteins to active colloids and shaken granular rods [1,2]. Active systems have recently received considerable attention in the community because of their potential in designing mesoscopic engines, as a way of interpreting biological mechanics, and as examples of non-equilibrium statistical physics [3,4]. Here we focus on dense active systems, and in particular the continuum theory of active nematics, which describes active systems with hydrodynamic interactions. The archetypal experimental example is microtubules driven by kinesin motor proteins [5,6]. Other active nematics include myosin-driven actin-microtubule networks [7], swimming bacterial swarms and confluent eukaryotic cells [8][9][10][11].
A key property of active nematics, which distinguishes them from passive liquid crystals, is active turbulence. This is a chaotic flow state characterised by strong vorticity and motile topological defects which are continually created and destroyed. Considerable experimental and theoretical work has been devoted to understanding the properties of active turbulence, and the associated topological defects, in two dimensions. More recently it has proved possible to design an active material that allows imaging of a three-dimensional active nematic and, in particular the associated motile disclination loops and lines [12]. This was achieved by dispersing forcegenerating microtubule bundles in a passive colloidal liquid crystal based on filamentous viruses. The temporal evolution of disclination lines was measured using lightsheet microscopy revealing that the primary topological excitations in bulk active nematics are charge-neutral disclination loops. * liam.ruske@physics.ox.ac.uk Numerical simulations have revealed that flows and morphological dynamics of disclination lines in three dimensional (3D) active nematics are governed by the local director profile surrounding the disclination line and that defect loops in extensile systems are generally formed via the well-known bend instability [13]. Three-dimensional active nematic turbulence and disclination line dynamics in spherical confinement has also been investigated using numerical modeling in both achiral [14] and chiral systems [15,16]. Enforcing strong in-plane surface alignment allowed the formation of defects on the surface and was used to highlight the coupling of surface and bulk topological defect dynamics.
The theories of active materials are increasingly being used to describe biological systems in two dimensions, with examples including biofilm initiation [17], topological defects in cell monolayers [8][9][10] and epithelial expansion [18]. This suggests that in three dimensions there may be relevance to the collective motion of groups of cells, in morphogenesis or to the growth and spread of tumours. Therefore to underpin extension of these approaches to 3D, in this paper we investigate active, self-deforming droplets in three dimensions, describing the interplay between disclination dynamics and droplet morphology. We find that extensile droplets form protrusions at points where disclination lines reach the surface. In contractile droplets, however, dislocations lead to a wrinkled drop surface. Moreover for small, contractile droplets, we find that invagination, to form a cup-like configuration, can be driven by nematic activity alone. We present evidence explaining the reasons behind the different behaviours stressing, in particular, the role of active anchoring, an effective surface alignment resulting from active flows.
We start by presenting the mathematical description of an active nematic and the equations of motion which we solve numerically, and follow this by a summary of 3D active nematics and the classification of disclination lines. The main results are presented in section IV, which is split into four subsections: First we investigate the nematic structure on the surface of an active droplet (A) and how this affects the dynamics of disclination lines in the bulk of spherical droplets (B). We then move on to show how active forces can deform the interface by highlighting mechanisms causing droplet deformations triggered by both extensile (C) and contractile (D) activity. The last section of the paper summarises the key results and points out possible connections to biological systems.

II. EQUATIONS OF MOTION
We summarise the continuum equations of motion describing the dynamics of a deformable, nematic droplet in an isotropic fluid background. In the absence of activity the system is described by a free energy F = f dV = f GL + f LC dV. The first contribution, which controls the formation of the nematic droplet, is chosen to take the Ginzburg-Landau form This describes phase separation into two stable phases with concentrations ϕ = ±1 with an interface of width separating the two phases [19,20]. The bending rigidity κ of the interface is related to κ * by κ * = (4 3 /3 √ 2)κ. The final term penalizes gradients in the concentration field and introduces a surface tension σ ∝ K ϕ .
The second contribution to the free energy density is which includes the usual Landau-de Gennes bulk energy of the liquid crystal expanded in terms of the nematic tensor Q which has elements Q ij = 3S 0 /2(n i n j − δ ij /3), where S 0 is the magnitude and n the direction of the nematic order, and a term which penalizes elastic deformations of the director field [21]. A LC sets the bulk energy scale and we use the one-elastic-constant approximation which assigns the same elastic constant K LC to twist, splay and bend deformations. The function η(ϕ), which quantifies the dependence of bulk energy on the the local liquid crystal concentration predicts a first order phase transition between a nematic and an isotropic phase of the liquid crystal at the value η c = 2.7. We define η(φ) := η c + η 1 ϕ, and choose a value of η 1 that ensures that the system forms a nematic droplet (ϕ = +1) which is surrounded by isotropic fluid (ϕ = −1).
Since the total concentration ϕdV is assumed to stay constant, diffusive transport follows Model-B dynamics and the time evolution of ϕ is governed by the following reaction-diffusion equation [22]: Here u is the velocity field and the mobility Γ ϕ quantifies how fast ϕ responds to gradients in the chemical potential µ = δF/δϕ. Unlike the total concentration, the local nematic alignment Q is not a conserved quantity and its time evolution follows modified Model A dynamics [23]: where Γ Q is the rotational diffusivity and H, the molecular field, is given by Rod-like particles can not only be advected by the fluid, but also rotate in response to flow gradients. This behaviour is accounted for by the co-rotational term [23] where D ij = (∂ j u i + ∂ i u j )/2 and Ω ij = (∂ j u i − ∂ i u j )/2 are the symmetric and antisymmetric parts of the velocity gradient tensor W ij = ∂ i u j , respectively. The parameter ξ, determines whether the director aligns with, or tumbles in, a shear flow. Its value depends on the molecular details of the liquid crystal, and it is related to the flow alignment parameter λ defined within Ericksen-Leslie theory by λ = 3S0+4 9S0 ξ. This work mostly focuses on flow tumbling nematics, for which the flow alignment parameter λ < 1.
We assume that the build-up of active flows occurs over time scales which are much longer than any viscoelastic relaxation time of the system. We therefore consider the liquid limit and solve the incompressible Navier-Stokes equations to obtain the flow field u: where the stress tensor Π = Π viscous + Π elastic + Π capillary + Π active . The passive contributions, well known from liquid crystal hydrodynamics [24], are: where ρ is the density, η the viscosity, p the bulk pressure andQ = Q + 1 3 I . In addition to the passive terms, the stress due to the dipolar flow fields produced by the active particles is [25] Π active = −ζQ , where ζ quantifies the magnitude of active stress. For extensile activity, ζ > 0, nematic particles direct the fluid outwards along their elongated direction and inwards along the two perpendicular axes. The flow direction is reversed for contractile activity, ζ < 0. We solved the equations of motion using a hybrid lattice Boltzmann-finite difference method. This involved solving eqs. (3,4) using finite difference methods, and eqs. (7,8) using a lattice Boltzmann algorithm [21,23,26,27]. Simulations were performed on a threedimensional lattice of size 100x100x100 and discrete space and time steps were chosen as unity. We used periodic boundary conditions for the simulation box and the system initially consisted of a spherical, nematic drop of radius R = 30 (ϕ = +1) embedded in an isotropic fluid (ϕ = −1), unless otherwise stated. Initially we let the system relax for t = 500 time-steps such that the droplet interface and the nematic tensor Q reached their equilibrium profile before we switched on activity. We used the following parameter set for all simulations:

III. DISCLINATION LINES IN THREE-DIMENSIONAL ACTIVE NEMATICS
Active stress leads to an instability of the nematic phase ( Fig. 1). This hydrodynamic instability constantly pushes the system out of equilibrium, leading to a chaotic steady state termed active turbulence [28]. 3D active turbulence is characterised by spatiotemporally chaotic flows and the presence of disclination lines which constantly undergo transformation events such as breakup, recombination, nucleation and annihilation [14]. In bulk systems disclination lines typically form closed, chargeneutral loops [12]. However, in the presence of a confining interface, as for nematic droplets, disclination lines can also terminate at the boundary and the dynamics of the resultant defects on the surface is coupled to disclination line dynamics in the bulk by elastic interactions and flows.
Unlike in two-dimensional active nematics where ±1/2 defects carry topological charge and can therefore only FIG. 1. Bend (a) and splay (b) perturbations of the director field n cause active forces which set up active flows which further enhance or stabilise the respective perturbation. Forces in extensile or contractile systems are denoted by red or blue arrows, respectively. It is apparent that extensile (contractile) systems are unstable to bend (splay) deformations.
nucleate or annihilate in pairs of opposite charge, disclination lines in three dimensions can continuously transform from a local −1/2 configuration (in the plane perpendicular to the line) into a +1/2 configuration through an intermediate twist winding as indicated in Fig. 2. As one moves around the core of a disclination in the plane perpendicular to the local disclination line segment, the director field winds around a specific axis, the rotation vector Ω (black arrows) by an angle π. The angle β between Ω and the local line tangent t (yellow arrow) is called the twist angle and can be used to locally characterise the disclination line. For −1/2 (+1/2) wedge-type defects the twist angle corresponds to β = 0(π), while line segments with local twist-type defects are indicated by β = π/2.
Due to activity disclination lines act as self-propelled entities moving through the fluid. Based on a simplified model neglecting elastic interactions, each disclination line segment can be associated with a local self-propulsion velocity which depends on the local director profile [13]. The self-propulsion velocity component perpendicular to the local tangent of the line depends on twist angle β as Thus, unlike in two-dimensional active turbulence where the dynamics is dominated by two species of quasiparticles (±1/2-defects), disclination lines in threedimensional active turbulence act as quasi-particles with an intrinsic degree of freedom, the twist angle 0 ≤ β ≤ π. Line segments with β = 0 are passive while β = π line segments are most active in pushing around the surrounding fluid.

A. Activity leads to preferred director alignment at an interface
In 2D active nematics, the flow induced by active stresses leads to an alignment of the director at an interface with an isotropic phase [29]. This active anchoring of the director field is parallel to the interface for extensile stress and perpendicular for contractile stress. For 2D systems, the active alignment can be significant, but it is not immediately obvious whether the effect will persist in a 3D geometry, or how it will be affected by any defects present on the surface. Therefore, to investigate the effects of active anchoring in 3D, we measure the angle θ between the director field n and the surface normal at the interface of spherical droplets. Since there are more possible configurations for in-plane alignment (θ = π/2) than for perpendicular alignment (θ = 0), we use the distribution of cos θ to quantify alignment effects. A randomly aligned director field results in a uniformly distributed cos θ ∼ U[0, 1], while preferred perpendicular or in-plane alignment leads to bias of the distribution towards 1 or 0, respectively. The results in Fig. 3 show clear evidence that both extensile and contractile activity lead to strong active anchoring on the surface, which weakly depends on the flow-alignment parameter λ (Fig. S1 [30]).
The director orientation also shows signatures of the places where disclination lines end at the surface. Consider first the case of extensile activity. The director lies in plane over most of the surface, although there are localised regions with perpendicular alignment (see Fig. 3(a), Movie S1 [30]). The disclination lines present in the bulk of the droplet create quasi two-dimensional defects at the positions where they terminate at the interface which we shall term surface defects. The distribution of the twist angle β of surface defects has distinct peaks at β = 0 and β = π, showing that most defects on the surface are of wedge type (Fig. 4), corresponding approximately to 2D +1/2 and -1/2 configurations. There are less twist-type surface defects with β ≈ π/2 as these introduce a perpendicularly aligned region in the vicinity of the defect, which is suppressed by active anchoring.
A two-dimensional nematic sheet confined to the surface of a sphere always has at least four +1/2 defects present as the total topological charge is conserved and must add up to two, the Euler-characteristic χ of a sphere. Active anchoring does not confine the director field to the surface everywhere so the topological charge is not strictly conserved.
FIG. 3. Surface alignment of the director field for extensile (a) and contractile activity (b). Surface defects, shown as blue arrows, are connected via disclination lines running through the bulk. The surface alignment is indicated by colour bar, where orange (black) indicates in-plane (perpendicular) director alignment with respect to the interface. Upper right: surface alignment in projection. Lower right: distribution of surface angle cos(θ) over the total surface area By contrast, contractile activity leads to a remarkable stripe pattern in the director alignment on the surface. While most parts of the surface show strong perpendicular (homeotropic) active alignment, this is interspersed with thin stripes of clear in-plane ordering (see Fig. 3(b), Movie S1 [30]). The stripes start and terminate at sur- face defects which are mostly twist-type with β ≈ π/2 (Fig. 4). Wedge-type surface disclinations with β = 0, π are suppressed as they would create a region of in-plane alignment in the vicinity of the defect. Twist-type disclinations on the other hand only introduce a small region of in-plane alignment along one specific direction (Fig. 2) and are therefore favoured. Adjacent twist-type disclinations tend to align with their rotation vectors antiparallel, in a way that minimises the area of unfavourable in-plane alignment, thereby creating the observed stripe pattern (Fig. 5).
Perfect perpendicular alignment at the surface would prevent disclination lines in the bulk terminating there as all disclination lines introduce some degree of in-plane alignment at the surface. Indeed, because topological charge is conserved, a sphere with perfect homeotropic alignment at the interface would force the system to form a +1 defect loop in the bulk [13,14]. However, +1 defect loops in the bulk are associated with large elastic energy of the liquid crystal which usually cannot be overcome by active anchoring.

B. Disclination line dynamics in spherical droplets
We now focus on disclination line dynamics inside spherical droplets which are essentially undeformed by active forces. In the regime of active turbulence, disclination lines constantly form and annihilate. Previous investigations of active turbulence in bulk systems found that the dominant excitations of three-dimensional active nematics are charge-neutral disclination loops which undergo complex dynamics and recombination events [12]. In droplets however disclination lines do not need to form loops as they can also exist as growing or shrinking half circles or lines which terminate at the surface at positions [x 1 , x 2 ]. The distance between endpoints (surface defects) of disclination lines D = |x 1 − x 2 | first scales linearly with the total length L of lines before plateauing at a finite value as the lines increase in length. For disclination lines which are short compared to the droplet radius (L < R), the mean endpoint separation scales as D ≈ (2/π)L, showing that the lines mainly nucleate or annihilate at the droplet's surface as half-circles. The endpoints of very long disclination lines (L R) are randomly distributed over the droplet's surface so the mean separation of endpoints converges towards D/R ≈ 4/3 ( Fig. S2 [30]). The local properties of disclination lines can be classified by the twist angle β which varies continuously along the line. We find that the director profiles of disclination lines close to the surface are heavily influenced by active anchoring which favours wedge-(twist)type disclination lines for extensile (contractile) activity (Fig. 6, Movie S1 [30]). This is in contrast to the bulk of the drop where extensile (contractile) activity favours twist-(wedge)-type disclinations. The distributions of β at the surface and in the bulk are compared for extensile and contractile droplets in Fig. 6(c) and Fig. 6(d) respectively.
Behind the apparent disorder of three-dimensional active turbulence we identify an active length-scale ζ ∝ K LC /ζ which controls the average disclination line density and coincides with the active length scale governing the decay of the vorticity correlations in twodimensional active turbulence [3]. Fig. 7 shows the area density of disclination lines, defined as the total length L tot of all disclination lines divided by the droplet volume V, as a function of K LC and ζ, confirming that it scales with the active length-scale as −2 ζ . Moreover, measuring the mean variation of twist angle β along disclination lines shows an exponential convergence towards a plateau with a characteristic length scale ξ β (Fig. S3  [30]). We find that in the turbulent regime, the correlation length ξ β of twist angle β scales with the active length-scale as ξ β ∝ As a consequence of active anchoring, for extensile systems disclinations close to the surface tend to be wedge-type whereas those in the bulk tend to be twist disclinations as is evident from the distribution of β.
(d) For contractile systems disclinations close to the surface slightly tend to be twist-type whereas those in the bulk tend to be wedge disclinations.
three-dimensional space, hence correlations do not decay linearly along lines. Based on these observations, we define the dimensionless activity number A = R ζ/K LC as the ratio of droplet size R to active length-scale ζ . The mean curvatureκ of disclination lines increases with activity number roughly asκ ∼ A 0.6 .
Due to activity disclination lines move through the fluid with an associated transverse self-propulsion velocity v SP ⊥ (eqn. 13). In addition, they are passively advected by the surrounding flow and thus the mean transverse velocity of disclination line segments is ex- where v rms denotes the root-mean-square velocity of the fluid inside the droplet. Surprisingly this simplified model of uncoupled, self-propelled line segments closely matches the velocity profile we observe in turbulent droplets for large β (Fig. 8(a)). As disclination lines consist of line segments with varying self-propulsion velocities, more active line segments are elastically coupled with more passive segments. Passive line segments are thus pulled around by elastic interactions in addition to passive advection by the flow, which is the reason for the deviation of v ⊥ from the theoretical prediction for small β. We have argued that, as a consequence of active anchoring, the distribution of β varies between disclination lines near the surface and those deep into the bulk (Fig. 8(b), solid lines). This is mirrored in the average self-propulsion velocity (Fig. 8(b), diamonds). In extensile systems the wedge-type defects close to the surface on average move faster than the twist-type line segments in the bulk. Similarly, contractile activity causes line segments in the bulk, which are mostly wedge-type, to move faster than surface defects. The relative position is measured as the ratio of Γ, which is the path length between each point and the nearest end-point of the disclination line, and L 1/2 , which is half the total line length.

C. Extensile activity triggers the formation of finger-like protrusions
We now consider lower values of the surfaces tension such that the active flows are strong enough to deform the droplet, which leads to the formation of defect-mediated, finger-like protrusions for extensile activity. Unlike in stiff droplets, the active flow field around disclination lines close to the surface can push the interface outwards creating a bulge along the self-propulsion direction of a disclination line ( Fig. 9(a)). As the disclination line continues to move outwards towards the interface, the local bulge extends and forms a thin protrusion. Due to the reduced separation of interfaces compared to the spherical droplet, the properties of disclination lines inside protrusions are dominated by active anchoring which favours in-plane surface alignment. This causes disclination lines to span the width of the protrusion as a straight line with most line segments resembling wedge-type +1/2 defect profiles (β ≈ π), thereby enhancing the directed selfpropulsion of the disclination line even further. The +1/2 defect profile aligns along the protrusion axis as this configuration satisfies the in-plane anchoring condition everywhere at the surface ( Fig. 9(b,d)). This quasi-two dimensional movement of +1/2 defect lines has also been observed in active nematic films confined to a channel below a critical wall separation [31]. Wedge-type −1/2 defect profiles (β ≈ 0) are not observed in protrusions as they are passive and lack the self-propulsion necessary to form a bulge in the first place.
The straight +1/2 defect line eventually reaches the tip of the growing protrusion and moves out of the nematic droplet, leaving behind a homogeneous director field which is aligned along the protrusion axis ( Fig. 9(c,e)). Thereby an area of perpendicular surface alignment is introduced at the end of protrusions, which are points of large negative mean curvature. The surface alignment cos(θ) is therefore correlated to the local mean curvature of the interface (Fig. S4 [30]). In the absence of disclination lines mediating active forces, the aligned protrusions slowly retract due to surface tension and bending energy. The constant formation and retraction of droplet protrusions is also shown in Movie S2 [30].
To further quantify the mechanism of protrusion formation we measured several properties of disclination lines as a function of radial position x r ,with the droplet's centre of mass being the reference point of a spherical coordinate system. For approximately ellipsoidal droplet shapes, x r can be used as a proxy to divide the initially spherical droplet of radius R into a bulk domain (x r < R) and a protrusion domain (x r > R). Disclination lines in the bulk (x r < R) of soft droplets undergo chaotic movement ( v r ≈ 0) and are mostly twist type. By contrast, in protrusions (x r > R) disclination lines mostly consist of +1/2 line segments (≈ 70%) and show persistent selfpropulsion along the radial protrusion axis ( v r 0, see Fig. 10(a)). Disclination lines in protrusions are nearly straight and their total length L is limited by the width of the protrusion (L(x r ) ≈ const for x r > R, see Fig. 10(b)).
Droplets of sizes much larger than the active lengthscale R ζ are strongly deformed and protrusions do not always grow along the radial axis ( Fig. S5(a) [30]), thereby rendering the spherical approximation unsuitable. Still it can be observed that +1/2 defect lines are much more frequent in soft droplets with protrusions than in spherical droplets without protrusions (Fig. S5(b) FIG. 9. Formation of finger-like protrusions by motile disclination lines shown by three snapshots at times t1 < t2 < t3. Panels (a)-(c) show the three-dimensional droplet shape with colour coding showing the local surface alignment of the director field. Disinclination line segments are coloured according to the twist angle β with the same colour coding as in Fig. 6. (a) Disclination line with varying twist angle β moves towards the interface. Active flow pushes the interface outwards creating a bulge along the self-propulsion direction of the disclination line. (b) The small protrusion width combined with in-plane alignment at the surface stabilises the disclination line into an almost straight configuration with β ≈ π. (c) Disclination line moves out of the droplet leaving behind a defect-free protrusion with an aligned director field. Panels (d),(e) show a schematic diagram of the director field and active forces. The red arrow in panel (d) denotes the selfpropulsion direction of the β ≈ π disclination lines. As shown in panel (c), disclinations leave behind an area of perpendicular surface alignment (dark regions) at the end of protrusions which slowly retract due to surface tension and bending rigidity. Snapshots created from simulation with droplet size R = 15, see also Movie S2 [30]. [30]).
The morphology of extensile droplets is determined by two control parameters: The activity number A = R ζ/K LC controls the density of disclination lines in- whose length is limited by the diameter of the protrusion. The radial component of the self propulsion velocity vr was obtained by taking the droplet's centre of mass as the reference point of a spherical coordinate system. The fraction of +1/2 disclination lines was defined as the fraction of disclination line segments with twist angle 3/4π ≤ β ≤ π and the persistence of disclination lines as the ratio of endpoint distance over total line length.
side the bulk of droplets while the ratio of elastic constant to surface tension Ψ = K LC /K ϕ quantifies the energetic cost associated with nematic deformations in the bulk compared to deformations of the interface (assuming surface tension dominates the bending stiffness). The morphology of droplets as a function of A and Ψ is shown in Fig. 11. For very stiff interfaces (Ψ 1), droplets are nearly spherical and host disclination lines if the activity A is sufficiently large (blue diamonds, purple triangles). Soft interfaces (Ψ ≥ 1) however allow the formation of finger-like protrusions as soon as disclination lines are formed inside droplets (yellow circles). If active forces become much larger than the passive restoring force of the interface, active flows tear apart the droplet which breaks up into smaller parts (orange squares). The strength of protrusion formation can be quantified by the gyrification index, which is a function of A and Ψ (Fig. S8(a) [30]).

D. Contractile activity triggers droplet invagination and surface wrinkles
We now consider contractile activity which causes droplet invagination or the formation of comb-shaped droplet deformations creating wrinkle patterns on the surface. Unlike protrusion formation in extensile systems, droplet deformations in contractile systems originate from smooth director field deformations and are not mediated via disclination lines.
In contractile systems active anchoring favours normal surface alignment. Every spherical or topologically equivalent surface with normal surface alignment everywhere enforces the formation of at least one +1 disclination-loop in the bulk due to topological constraints. These loops are usually associated with a large elastic energy cost due to strong deformations of the director field and thus are only observed in drops if strong normal anchoring is enforced by additional terms in the free energy. Otherwise, in the absence of disclination lines, the director field forms a ring with in-plane surface alignment encircling the droplet (Fig. 12(a)) to maximize the area of perpendicular surface alignment favoured by active anchoring while avoiding the formation of a +1 defect-loop in the bulk. The ring of in-plane surface alignment is associated with nematic bend deformations in the bulk; thus we will refer to this structure as a bendring from now on (Fig. 12(d)). Contractile activity causes the bend-ring to push outwards, thereby deforming the droplet to an oblate ellipsoid until passive forces arising from surface tension and membrane rigidity counterbalance the active force.
If active forces are small compared to the passive restoring forces arising from membrane and director field deformations, the oblate bend-ring configuration is stable. Perturbations of ring shape and position are associated with increasing elastic deformations of the membrane and the director field, thereby creating an energy barrier which stabilises the configuration (Fig. S6 [30]). The exact free energy profile depends on model parameters such as the surface tension K φ , membrane rigidity κ, nematic elastic constant K LC and bulk properties.
If active forces are sufficiently strong to overcome the passive restoring force, however, the bend-ring at the equator is unstable. It contracts while moving towards one of the poles and its motion results in invagination leading to a cup-shaped configuration of the droplet.
This instability arises from the active flow set up by the director configuration of the bend-ring which pushes outwards along the deformation axis ( Fig. 1(a)). Any small deviations from a perfect ellipsoidal droplet shape cause a component of the active flow to point towards one of the poles. The resulting droplet deformation rotates the axis of the bend-ring further towards the pole thereby in turn further increasing the active force component in the direction of the pole (Fig. 12(b,e)).
As the contracting bend-ring approaches the pole, it creates an area of splay deformation in its centre due to active anchoring enforcing perpendicular surface alignment ( Fig. 12(e)). This splay deformation sets up an active flow pointing inwards ( Fig. 1(b)). The outward pointing forces of the bend-ring together with the inward pointing forces in its centre combine to drive droplet invagination ( Fig. 12(c,f)). If the magnitude of active stress ζ is not reduced after complete invagination, this process eventually leads to droplet break-up as the passive restoring forces resulting from membrane rigidity and surface tension are not sufficient to compensate the active forces. However, if the system is active only for a certain amount of time and ζ is reduced after the initial invagination, the cup-shaped drop can be stabilised, analogous to the natural course of morphogenetic events, which are controlled by biochemical signals and begin and end at a predefined time [32].
Alternatively, if the interface is too stiff for active forces to initiate drop invagination, droplets perform an active run-and-tumble motion (Fig. 13, Movie S3 [30]). Instead of a large cavity, only a small dip is formed at the centre of the bend-ring where the interface resists further deformations and the splay configuration of the director field leads to flows which push the droplet forwards. Eventually, the director deformations relax via the nucleation of a β = π half-loop and the self-propulsive disclination line moves through the droplet until it reaches the FIG. 12. Initiation of droplet invagination by contractile activity shown by three snapshots at times t1 < t2 < t3. Panels (a)-(c) show the three-dimensional droplet shape with colour coding indicating the local surface alignment of the director field. Panels (d)-(f) are schematic diagrams of the director field n, showing the active forces (red arrows) and the position of the bend-ring (red, dotted circle). (d) The ring of in-plane alignment at the equator causes nematic bend deformations in the bulk. Contractile activity causes these bend deformations to push outwards (red arrows), thus deforming the droplet to an oblate shape. (e) While the ring of in-plane alignment moves towards one of the poles, the bend deformations in the bulk point perpendicular to the surface which causes the active forces to develop a component pointing towards the pole. This creates a dent in the centre of the ring. (f) Due to the perpendicular surface alignment caused by active anchoring, the dent in the centre of the ring leads to a director splay deformation in the bulk. Contractile activity causes this splay deformation to push inwards, thus increasing the depth of the dent. Snapshots created from simulation with elastic constant KLC = 0.4. opposite interface. This re-orientates the bend-ring and the process repeats along a random direction.
In contractile droplets which are larger than the active length-scale R ζ , many disclination lines are present in the bulk and the surface shows characteristic stripes of in-plane director alignment connecting the endpoints of disclination lines (Fig 3). The stripes on the surface are associated with nematic bend deformations in the bulk which induce active forces pushing outwards. This creates comb-shaped deformations along FIG. 13. Droplets perform an run-and-tumble motion if the interface is too stiff for active forces to cause complete droplet invagination. Instead, a bend stripe at the equator deforms the droplet to an oblate shape (a), which then shrinks, forming a dip at the centre of the bend-ring where a splay deformation in the bulk pushes the droplet forward (b). Then the elastic deformations relax via the nucleation of a β = π halfloop (c). The self-propulsive disclination line moves through the droplet and eventually reaches the interface. Thereby the bend-ring is re-oriented along a random direction and the process repeats (d). Snapshots created from simulation with droplet size R = 15, see also Movie S3 [30].
the stripes of in-plane alignment, resulting in a surface wrinkle pattern which we will term active wrinkling ( Fig  14, Movie S4 [30]). In addition the droplet deformation causes the perpendicularly aligned surface areas between in-plane stripes to form splay deformations in the bulk, which cause inward pushing, active forces that create dimples and valleys, further enhancing the wrinkle pattern. Along the in-plane ridges the mean curvature is negative while in the centre of holes it is positive. Surface alignment cos(θ) is therefore correlated to the local mean curvature of the interface with in-plane surface alignment (cos(θ) ≈ 0) associated with points of negative mean curvature and perpendicular surface alignment (cos(θ) ≈ 1) predominating at points of positive mean curvature (Fig. S7 [30]).
As for the extensile case, the morphology of contractile droplets is determined by the activity number A = R ζ/K LC and ratio of elastic constant to surface tension Ψ = K LC /K ϕ (Fig. 15). Very stiff interfaces (Ψ 1) cause droplets to be nearly spherical and host disclination lines if activity A is sufficiently large (blue  (Fig. 3). Bend and splay deformations of the director field near the surface create comb-shaped deformations along the stripes of in-plane alignment, resulting in a wrinkle pattern. In between stripes of in-plane alignment, dimples and valleys are driven into the drop as shown in snapshot (c) which depicts the local radius on a semi-transparent drop. The dynamical structure of active surface wrinkles is shown in Movie S4 [30], together with a 3D view of the surface structure.
diamonds, purple triangles). If interfaces are sufficiently soft (Ψ ≈ 1), droplets can be either be static, oblate ellipsoids without disclination lines (blue diamonds), perform an run-and-tumble motion (yellow circles) or form active surface wrinkles (green stars). For very soft interfaces (Ψ 1), droplets are static ellipsoids at low activity and full droplet invagination takes place at larger activity, leading the droplet to eventually break-up (orange squares). The degree of invagination and wrinkling can also be quantified by a gyrification index, which varies with A and Ψ (Fig. S8(b) [30]).

V. DISCUSSION
In this paper we have investigated the morphology and defect dynamics in active deformable droplets in three dimensions, considering both extensile activity, as present in experimental systems such as MT-kinesin mixtures [33], human bronchial epithelial cells [34] and Madine-Darby canine kidney (MDCK) cells [8], as well as contractile activity, which is found in systems such as mouse fibroblast cells [10] or actomyosin gels [35]. We have shown that active stresses cause in-plane (extensile) or perpendicular (contractile) surface-alignment of the director field at the interface of droplets. Unlike thermodynamic surface anchoring, which involves an anchoring free energy, active anchoring is a hydrodynamic effect due to the flows driven by gradients of the nematic ordering Q at the interface. Active anchoring has important consequences for the dynamics of disclination lines close to the interface, such as the preferential formation of wedge-like surface disclinations in extensile systems or the creation of in-plane aligned stripes connecting twist-like surface disclination in contractile systems. This in turn triggers the formation of finger-like protrusions in soft, extensile droplets or the creation of surface wrinkles and droplet invagination in contractile systems. Furthermore, we applied the ideas of 2D active turbulence, such as defect velocities and the existence of an active length-scale, to three dimensions. We identified an active length-scale ζ of three-dimensional active turbulence, which controls both the density of disclination lines and the correlation length of the twist angle β along disclination lines in the bulk of nematic droplets.
Recent work has interpreted several 2D biological systems including bacteria biofilms [17,36], epithelial tissue [8], spindle-shaped cell monolayers [10] and actin fibers in regenerating Hydra [37] in terms of the theories of active nematics. It is interesting to ask whether similar ideas may prove useful in 3D to describe biological systems in which nematic constituents, such as protein filaments, eukaryotic cells or bacteria, are organised as solid, 3D structures with a confining interface. There are several examples where cells migrate collectively as a cohesive group, thereby maintaining supracellular properties such as collective force generation and tissue-scale hydrodynamic flow [38][39][40]. For example, in some modes of collective cell movement implicated in cancer invasion, a blunt bud-like tip consisting of multiple cells that variably change position and lack well-defined leader cells protrudes from the tumour [41,42]. Hence, the natural appearance of finger-like protrusions in active nematic droplets could serve as a new approach to explain collective cell invasion of tumors.
In eukaryotic cells, membrane shapes under mechanical stress are mostly controlled by the mechanics of the cortical actin cytoskeleton underlying the cell membrane which produces contractile stresses. It has been observed that membrane shapes mainly depend on the actin thickness, where thin shells show a cup-shaped deformation and thick shells produce membrane wrinkles [43]. Our work provides a possible link between contractile activity and the emergence of cell-surface ruffles, circular dorsal ruffles (CDRs) and caveolae [44][45][46] at the surface of cells. Surface invagination which plays a role in active transport through cell membranes might also be related to local contractile forces. During macropinocytosis (fluid endocytosis), extracellular fluid is brought into the cell through an invagination of the cell membrane forming a small vesicle inside the cell. Vesicles form from cell-surface ruffles that close first into open cups (ruffle closure) and then into intracellular vesicles (cup closure) [47,48], which is reminiscent of the dimples that are driven into contractile drops. Some morphologies observed in contractile drops, such as partial invagination and run-and-tumble motion are known to also occur in active polar droplets [49,50]. In contrast to nematic systems, which are characterised by a headless director field n ↔ −n, polar systems are described by a polarisation vector p. As a result, polar models do not exhibit ±1/2 topological defects, which have been shown to contribute to the dynamics and flows in cell cultures and bacterial biofilms [8,10,17,34,36].
Up to now, most of the research on three-dimensional active nematics has been focused on the topology and dynamics of closed disclination loops in the bulk or disclination lines in fixed geometries [12][13][14][15][16]. Previous investigations of two-dimensional active nematics which are confined to the surface of a 3D shell explored the connection between 2D topological defects and shell morphology [51]. With this work we expanded the current understanding of three-dimensional active nematics, considering the interplay between motile disclination lines, deformations of the nematic director field and the droplet interface.
Future directions include the consideration of more complex geometries, in particular thick, nematic shells or the addition of an anchoring free energy enforcing a given surface alignment. The motivation behind this is that most biological systems during embryonic development form sheets or shells rather than solid structures. Another aspect when considering cell aggregates such as embryos or tumors is that they are often composed of different cell types with different mechanical properties and cellular responses [52]. It would be interesting to extend our model to two or more distinct, nematic fluid phases with different mechanical properties and activity coefficients, thereby incorporating cellular heterogeneity within cell aggregates. Another important aspect is that cells are usually surrounded by a three-dimensional dense network of macromolecules, the extracellular matrix, which provides structural support and is known to strongly affect cell migration behavior [53,54]. The investigation of how a viscoelastic medium or a dense polymer solution affects active droplets therefore appears to be a promising direction.