Instabilities, motion and deformation of active fluid droplets

We consider two minimal models of active fluid droplets that exhibit complex dynamics including steady motion, deformation, rotation and oscillating motion. First we consider a droplet with a concentration of active contractile matter adsorbed to its boundary. We analytically predict activity driven instabilities in the concentration profile, and compare them to the dynamics we find from simulations. Secondly, we consider a droplet of active polar fluid of constant concentration. In this system we predict, motion and deformation of the droplets in certain activity ranges due to instabilities in the polarisation field. Both these systems show spontaneous transitions to motility and deformation which resemble dynamics of the cell cytoskeleton in animal cells.


I. INTRODUCTION
In animal cells, motility and morphology are strongly coupled and are largely due to the activity of the cell cytoskeleton. Research into these areas is broad and has many applications, from studying metastatic cancer cells to wound healing. In order to mimic aspects of these systems we model, both analytically and numerically, examples of active cytoskeletal material confined to droplets. An active material is defined as driven out-of-equilibrium by the internal energy of its constituent particles [1]. We use the hydrodynamic model of an active polar fluid outlined in [2][3][4] to model the behaviour of such a material at long length and time scales.
Over the past decade there have been a number of calculations of instabilities and non-equilibrium steady states in active liquid crystals; thin or 2D flat films [2,[5][6][7][8][9], thin cortical layers [10][11][12][13], confined in emulsion droplets or vesicles [14][15][16][17][18][19][20][21], and simplified models of animal and plant cells [22,[24][25][26][27][28]. In this paper we model deforming active droplets immersed in a passive fluid using linear perturbation theory. By making justified assumptions, we are able to predict nonequilibrium phase transitions in both of the systems we consider, and predict how the droplet deformation couples to these. These analytical calculations are presented for the three-dimensional case and also repeated for the two-dimensional analogue where we find qualitatively similar results. Numerical simulations use the twodimensional Immersed Boundary method used in [29] and are directly compared to the two-dimensional analytical calculation.
The models presented here are relevant to active systems in vitro (constructed using techniques in [30][31][32]) as well mimicking aspects of cell dynamics. The two cases we consider correspond to two limits of active cytoskeletal behaviour (see figure 1) that represent the minimum degrees of freedom required to observe interesting out- Firstly, we consider an isotropic layer of contractile active material confined to an interface between two fluids, which has physical similarities to the actomyosin cortex in cells. The stresses generated are confined to the plane of the interface giving rise to flows in the surrounding fluid and deformation of the interface itself. Interestingly, diffusion of the active particles through the bulk can result in a change in which mode of the perturbation has lowest critical activity, from a single peak instability driving droplet motion to higher modes which produce symmetric deformation. Furthermore, simulations show that advection through the bulk can stabilise such modes. This suggests that droplets with an active interface could spontaneously deform and possibly divide due to the feedback from the fluid flow.
Secondly, we consider a highly ordered active polar liquid crystal confined inside a fluid droplet. In this case the polarisation gradients direct the internal stresses giving rise to fluid flow. A polar anchoring condition at the in-terface means that the deformation of the droplet and polarisation field are strongly coupled. We find in this case there is a separation of swimming and stationary deforming modes, such that extensile activity destabilises the defect position and results in a swimming drop, whereas a contractile activity stabilises the centred defect position and gives rise to deformations of the interface.

II. ACTIVE FLUID INTERFACE
In this section we consider a fluid droplet coated by active particles on its interface that are isotropically ordered. Such systems have been found to self-organise in in vitro experiments using reconstituted active cytoskeletal material contained in vesicles or droplets [33,34]. These experimental systems are a useful tool for understanding the more complex dynamics of cells. The model in this section makes predictions of interesting active phenomena including symmetry breaking, and droplet deformation, that are relevant to the field of cell mechanics.

A. Model
We consider a fluid droplet described by an interfacial surface Σ separating the contained fluid domain Ω 0 and external fluid domain Ω 1 with viscosities η 0 and η 1 respectively. We define a concentration of active matter c(θ, φ, t) on the interface Σ, which alters the droplet surface tension γ such that γ = γ 0 − ζ c c − Bc 2 /2. γ 0 is the bare surface tension, ζ c is the activity (ζ c < 0 for contractile) and B is a passive repulsion force. This higher order repulsive term represents passive pressure, similar to that in [12], which parametrises the compressibility of the active fluid on the interface. We denote the effective surface tension γ 0 = γ 0 −ζ c c 0 −Bc 2 0 /2, which is the value of γ in the stationary state.
The force density on the droplet interface is then: F = κγn + (∇ s γ)t i , wheren =n(θ, φ, t) is the outward surface normal,t i =t i (θ, φ, t) are the orthogonal surface tangent vectors, κ = ∇ ·n is the local curvature, and ∇ s = (t i · ∇) is the surface gradient. It is useful to define the effective activityζ = ζ c +Bc 0 which defines the scale of the force F for small deviations of the concentration c from c 0 . Thus, the interface has net contractility forζ < 0.
The only forces acting on the system originate at the droplet surface Σ, with position R = R(θ, φ, t)ê r assuming this is single-valued with respect to the angular coordinates (θ,φ). Thus, the resulting force density in the fluid is f ext (r, θ, φ, t) = F δ [r − R(θ, φ, t)]. We ignore inertia taking the low Reynolds' number limit, Re = 0, thus the incompressible fluid flow (∇ · v = 0) is described by Stokes' equation η n ∇ 2 v + f ext − ∇P = 0, where n = 0, 1 denotes the domain Ω 0 or Ω 1 , v = v(r, θ, φ, t) is the fluid velocity, f ext = f ext (r, θ, φ, t) denotes any external force densities and P = P (r, θ, φ, t) is the hydrostatic pressure. We take the limit of a zero-thickness interface and assume flow and stress continuity between the two fluids Ω 0 and Ω 1 . This means the active particles act as an active surfactant, rather than a thin viscous layer (as in [8, 11-13, 27, 28]), which allows us to study the dynamics of deformation in a 3D viscous environment analytically.
The evolution of the surface concentration c with respect to time t is: is the interface flow velocity, D is the diffusion constant for the active particles on Σ, and k on,off are binding and unbinding rates of the particles to the interface. The concentration of unbound particles in the bulk of the drop is denoted ρ = ρ(r, θ, φ, t). Binding occurs at the interface where we denote the concentration of unbound prticles ρ b = ρ(r = R, θ, φ, t). Note that k on has units of velocity, as it contains the adsorption depth parameter. We assume that the active particles are insoluble in the external fluid, and so the evolution of the bulk concentration ρ is given by: with the boundary condition D ρ (n · ∇)ρ = k on ρ − k off c at r = R, to ensure conservation of mass. The parameter D ρ is the bulk diffusion constant of the active particles.
Here we assume that the active particles only generate stresses at the interface, so the bulk concentration acts as a buffer to recycle the surface concentration.

B. Linear Stability Analysis
In this section we present the results of a linear perturbation to the stationary ground state of the droplet. The system is in a stationary (velocity v = 0) steady state when the interface is spherical (fixed radius R = R 0 ) with a homogeneous concentration of active particles (c = c 0 ). Then the bulk concentration is ρ 0 = k off c 0 /k on inside the drop, and the hydrostatic pressure inside is P = P ext + (2γ 0 −ζc 0 )/R 0 where P ext is the stationary state pressure in the external fluid. We perform a linear stability analysis by applying a small perturbation to the variables defined at the interface R and c of the form: are the spherical harmonic functions and δg lm g 0 . To first order, the resulting flow is given by Lamb's solutions for flow around a sphere, which can be expressed as vector spherical harmonics [35]. Solving the Stokes equation with flow and stress continuity conditions at the droplet interface gives expressions for δv (i) lm (as defined in [36] and Supplementary Information appendix A) in terms of δc lm and δR lm . The perturbation on the interface is also coupled to a perturbation of the internal concentration ρ such that We obtain analytical solutions for the stability by assuming a quasistatic solution for δρ (takingρ = 0). This assumption corresponds to a fast relaxation of the bulk concentration ρ compared to the timescale of evolution of the surface concentration c. At linear order, the solution for δρ simply satisfies the diffusion equation with a flux condition at the boundary: This solution enables us to predict the effect of the feedback by diffusion through the bulk analytically. The full solutions to the coupled linear equations are solved exactly with Bessel functions as in [11], however these solutions do not permit an analytical calculation of the stability condition, hence we do not consider them here, but instead compare our approximate analytical solutions directly with the full dynamical simulations.
Finally, we evaluate the coupled system of dynamic equations for the concentration (equation (1) in section II A) and radiusṘ = v b .n (the normal velocity at the interface) to first order in the perturbations. We find instabilities by looking for positive eigenvalues of the stability matrix that relatesċ andṘ to δc and δR to first order in the perturbations (see Supplementary Information appendix A for further details of this calculation). From this analysis we find an instability threshold for the effective activityζ < α I where whereη = (η 0 + η 1 )/2 is the mean viscosity of the internal and external fluid. We see that α I is independent of the effective surface tension γ 0 which shows that the coupled droplet deformation does not contribute to the symmetry breaking threshold. However, the corresponding maximum eigenvalue of the stability matrix does weakly depend on the effective surface tension γ 0 for l > 1. This weak positive relation suggests that the instability should evolve more quickly in large surface tension drops when l > 1. In this linear limit there is no contribution from the advection term in (2) and the second term in (3) (proportional to the binding rates) always increases the threshold. This is because the binding terms allows the concentration on the interface to be recycled by unbinding and diffusing into the bulk of the drop. The stability analysis shows how the droplet will initially deform. This deformation is characterised at short times by the maximally unstable mode l max , which can be found exactly when binding is not included (see figure  2 and Supplementary Information appendix A). At linear order the instability is independent of the spherical harmonic parameter m. Generically, l max predicts that as contractile activity is increased, the more concentration peaks will be initially formed on the droplet surface (figure 2). The total droplet activity scales with droplet size, and so l max is more sensitive to the activity parameterζ in larger droplets. Thus it is easier to observe modes with small l in smaller droplets, where the dynamics are less sensitive to small changes in the activity. Note that only the l = 1 mode (k = 1 in 2D) produces net propulsion of the droplet (i.e. ΣṘn dS = 0), so the first unstable mode corresponds to front-back symmetry breaking of the droplet profile.
As shown in Supplementary Information appendix A, one can approximate the maximally unstable mode l max analytically by solvingṘ = 0 for δR lm . This approximation imposes that R always assumes the steady state shape for a given fixed concentration perturbation δc lm (plotted in figure 2). Physically, this assumes that the shape dynamics are much faster than the concentration dynamics, and so can be taken to be quasistatic. Interestingly, while this assumption does not represent the full coupled dynamics of δc lm and δR lm , it does reproduce the critical activity threshold, and also approximates the mode structure well.
When binding is included (k off = 0) the dispersion relation changes, and as we see from (3) the active threshold is non-linear in l, and hence higher (non-swimming) modes can have lower activity thresholds than the l = 1 (swimming) mode.
Within the assumptions made here, the binding and unbinding dynamics always increase the activity threshold. We see that if the binding is fast k on D ρ , the critical activity takes the same form as the 1D model considered in [11] where the active threshold is always minimal for l = 1 and is proportional to the effective diffusion parameterD = (Dk on +D ρ k off )/k on . However, for fast bulk diffusion, geometrical effects become important. A single peak in the interfacial concentration gives rise to a concentration gradient in the bulk driving diffusion away from it. As the number of peaks on the interface increases the concentration gradients are more localised to the surface, and diffusion has a smaller effect. In this regime, the minimum critical activity can correspond to multi-peak modes (l > 1) when the contribution from bulk diffusion is significant. This is analogous to the findings in [8] for a one-dimensional active fluid consisting of two-components.
The droplet shape instability is enslaved to the concentration (as α I is independent of γ), so we can estimate how the shape will deform due to certain concentration distributions on the interface by solvingṘ = 0 for δR (for l > 1). Plotted in figure 2 is an example of these deformations and the associated flow to linear order. In order to calculate the resulting steady state dynamics we require numerical simulation.

C. Results and Comparison with Simulations
We test these analytical results against the 2D simulations developed in [29]. These use an Immersed Boundary method [37,38] to represent the active interface explicitly as a Lagrangian mesh which is coupled to the Cartesian mesh for the 2D fluid via a numerical Dirac delta function.
Repeating the stability analysis in 2D, we now take perturbations of the form g = g 0 + ∞ k=1 e ikθ . The calculation reveals that surface tension gradients do not deform the drop in 2D (as found in [40]) however the concentration dynamics remain very similar. We compare our predictions in 2D to the results of the Immersed Boundary simulations in figure 3. We run simulations varying the activity, binding rate (taking k off = k on ) and diffusion parameters. At zero binding we observe two steady phases, a stationary state and a steady moving state 3(a) separated by the threshold α I,2D which agrees well with the expected analytical result This moving steady state due to a surface tension gradient is also observed for the the self-propelled droplets studied in [40,41]. The equations of motion we use (see Model section) are similar to those for the self-propelled droplets studied in [40,41] and hence some of the same dynamical behaviour is observed. However, our model predicts new stable states and instabilities corresponding to pure deformation and division as discussed below. This arises due to the advection and diffusion of active particles through the bulk of the drop. Unlike in [40,41] the model here conserves the active particles within the drop making it more relevant to cell cortex dynamics. We next calculate the maximum mode number k max (see Supplementary Information appendix A). In the regime where we predict k max = 2, our simulations show initial formation of 2 peaks in droplet concentration. Without binding, these peaks are unstable and always coalesce to form one (as predicted for a flat active viscous layer in [8]). In this case, the droplet swims persistently and steadily with the concentration peak at its rear. A decomposition of the Fourier modes of this steady state shows that the far field flow is puller like, i.e. its dipole moment is such that it pulls the surrounding fluid inward and pushes it outward along the axis perpendicular to its motion. The activity threshold predicted compares well to that in the simulations for small values of the binding. At larger binding rate, the interior dynamics is not completely diffusion dominated, and the critical activity is underestimated due to the approximation oḟ ρ = 0. As we increase k off andζ we see that eventually the droplet becomes immobile with 2 stable peaks in the concentration (see figure 3). In the intermediate regime the droplet undergoes a 'wandering' motion as the concentration profile oscillates between a single peak and two peaks. Equation (4) predicts a non-trivial k dependence of the active threshold as binding terms become important. For the parameters used in figure 3, this can be seen by the crossing of the lines for the k = 1 and k = 2 modes, meaning that the minimum critical activity is not necessarily for the lowest k mode (k = 1). Note this is very similar to the prediction in 3D in (3).
The simulation results in figure 3 demonstrate that as the binding rate increases, advection of the concentration through the droplet bulk becomes more important. The advection can stabilise the two peaks at diametrically opposite points on the circle, resulting in a stationary droplet. However, we see that in 2D the drop does not deform, as the radial forces from the activity gradients are always cancelled by the hydrostatic pressure P . This is not the case for the full 3D system where we expect concentration gradients to deform the droplets as shown in figure 2. Nonetheless, the 2D simulations show that advection can stabilise the 2 peak configuration, which in 3D would result in symmetric deformation and potentially division of the droplet. Such a 3D simulation is beyond the scope of this work, but would be useful for quantifying the full 3D morphology. Recent work has shown that non-adherent cells exhibit a swimming state similar to the motion described here, and so it would be of interest to test in future work whether the steady state shape in 3D for the model here resembles the 'pear shape' observed in [28,39].

III. ACTIVE POLAR FLUID DROPLET
In this section, we consider a droplet filled with an active polar liquid crystal of constant density everywhere. Realising this system experimentally in droplet systems requires high concentrations of active material so that the polar to isotropic phase transition is localised to the droplet centre. This has been achieved in vitro for microtubule based active nematics but only in thin films thus far [31,32]. In these systems the measured order parameter is approximately constant everywhere except in the vicinity of topological defects. Thus we consider the limit where the active fluid is strongly polarised and restrict the analysis to only the orientational degrees of freedom of the active liquid crystal, and do not consider the density or polarisation magnitude degrees of freedom.

A. Model
We utilise the model of an active polar fluid developed by Kruse et al. in [2][3][4] which has similarities to other continuum models of the cytoskeleton on surfaces (such as [42,43]). We consider the case where the active fluid has strong local ordering and is far from the isotropic phase so that |p| = 1 everywhere (except at defects where p is undefined). This approximation is commonly used to model active and passive liquid crystal systems analytically.
In the Re = 0 limit the total stress in the active polar fluid, σ tot ij = σ visc ij + σ dist ij + σ act ij , has viscous, distortion and active contributions respectively where: The viscous stress is the response to flow assuming a Newtonian fluid. The distortion stress is that of a passive polar liquid crystal due to deviations in filament alignment, where the perpendicular part of the molecular field h i = −(δF/δp j )(δ ij − p i p j ) acts to minimise the free energy functional F = Ω+Σ d 3 rf with respect to p, given |p| = 1. The Ericksen stress, σ e ij = f δ ij − (∂f /(∂(∂ j p n )))(δ ij − p n p k )∂ i p k , is a generalisation of the hydrostatic pressure for complex fluids. Finally, the active stress represents the active dipolar force and thus is second order in p.
The free energy functional F gives the equilibrium properties of the system. Here for simplicity we use the one constant approximation of the Frank free energy: where K is the elastic constant and |p| = 1. Since we are modelling a finite droplet, the surface terms are important. We consider normal anchoring of the filaments to the interface, with surface distortion free energy density f s = W (p ·n − 1) 2 . This form of the surface free energy includes the 'spontaneous splay' term which is allowed in polar liquid crystals [44]. The polarisation flux iṡ where ω ij = (∂ i v j − ∂ j v i )/2 and Γ is the rotational viscosity.

B. Linear Stability Analysis
We contrast the model of an active interface to that of a droplet of active polar fluid of constant density. In this case, rather than the concentration of active particles, the important degree of freedom is the polarisation vector p denoting the average direction of the contractile filaments in the fluid.
We calculate the linear stability of the droplet in the limit of strong anchoring W → ∞ in order to study the effects between the coupling of droplet morphology and polarisation. This equates to the boundary condition p =n at r = R. In the case of weak or no anchoring, instabilities can occur for both extensile (ζ > 0) and contractile (ζ < 0) active polar drops as shown analytically in [45] and in simulations [15]. The condition of fixed polarisation at the interface inhibits certain deformations of the polarisation field at low activities and so the preferred deformation modes are those which can couple to the droplet deformation. This was demonstrated in 2D simulations of active nematic drops in [18]. Here we explain this mechanism analytically in a 3D fluid drop by linear stability analysis. The polar nature of the anchoring produces a "radial hedgehog" topological defect at the droplet centre (or a radial defect with +1 winding number in 2D), giving a simple analytical description of the stationary state. Thus we are able to make analytical predictions about spontaneous symmetry breaking in these systems even in the general 3D case.
Unlike the case of an active interface, the active fluid here fills the drop, and hence active and passive stresses are generated in the bulk. The stationary steady state is given by the polarisation p =r, R = R 0r , and v = 0.
To perform a general linear stability analysis, one would need to consider generic perturbations to both the polarisation field and interface and study the coupled equations for their evolution, this is not analytically tractable in this case. However, we can perform restricted perturbations that we expect to be representative of the dynamics in a particular limit. We consider the case where the polarisation field is enslaved everywhere to the shape of the boundary by the anchoring condition. This corresponds to the limit where bulk instabilities in the droplet are suppressed by its size (i.e. small droplets). In larger droplets, (or equivalently for smaller K) the dynamics of the polarisation field becomes more independent of the anchoring condition, and we expect this approximation to break down.
Due to the symmetry of the stationary state, we first need to consider the special case of the translational mode of perturbation, corresponding to the l = 1 spherical harmonic mode. Without loss of generality we consider a perturbation along the z-direction (m = 0). This mode implies a translation of the hedgehog defect away from the droplet centre. If we assume that the defect has some fixed finite core radius R c then we can treat the liquid crystal as contained between two boundary conditions, one at the defect r = R c and one at the droplet interface r = R 0 − δz cos(θ), where δz is a small displacement of the defect position from the droplet centre along the z-direction. The calculation is done in the reference frame of the defect so that it coincides with the origin of our coordinate system. In the equilibrium case (ζ = 0), we can write a polarisation field to first order that minimises the bulk free energy in (5) by solving h = 0 for these boundary conditions: This method equates the defect to a small colloid with (polar) homeotropic anchoring, and in the strong anchoring case we expect the free energy minimum to corre-spond to the defect being positioned at the droplet centre as we observe in simulations, and is reported in [46,47]. Using the polarisation in equation (7) we can estimate what the bulk free energy increase will be for such a deformation (details in Supplementary Information where = R c0 /R 0 is assumed small in the final approximation of the equation. This ∆F is positive for all , suggesting that the free energy minimum corresponds to the defect being positioned at the droplet centre. Note that this polarisation field is only valid to first order in δz and so higher order terms could affect the form of the quadratic term here. We now introduce a small activity ζ, such that equation (7) remains a valid approximation for the form of the polarisation field, then we see that this gives rise to active forces in the drop. We solve the force balance equations (omitting passive contributions, see Supplementary Information Appendix B) to find the active contribution to the flow. We then integrate to find the active contribution to the velocity of the defect core v c and droplet v drop . The relative velocity of the defect is then: We see that extensile activity (ζ > 0) always results in a relative defect velocity that is in the same direction as the initial defect displacement (alongê z ), as shown by figure  4. This implies that extensile activity will destabilise the defect from the centre and give rise to motion of the droplet as a whole (which to linear order is also alongê z ). Conversely, we expect contractile activity to stabilise the defect at the droplet centre, as the flows resulting from contractile activity (ζ < 0) act to restore the defect back to its stationary position at the droplet centre. Thus, within the assumptions made above, one can predict that the active polar droplet will break translational symmetry spontaneously above some finite activity. This mode of symmetry breaking is independent of surface deformations at linear order, and so its critical activity threshold should not depend on the droplet surface tension. Hence the critical activity threshold will only depend on the increase in the passive free energy (equation (8)), which goes to a finite value in the limit of a point defect and scales as the inverse of the droplet size. In general, the parameter is difficult to define, which is a consequence of the assumption of |p| = 1, which breaks down around the defect. This can be avoided by using a Landau-De Gennes type free energy description for the passive part of the dynamics such that there is an polar-to-nematic phase transition at the centre of the droplet. However, such an approach is not analytically tractable, as it requires solving non-linear partial differential equations for the radial dependence of p. Qualitatively though, the predictions here are consistent with what is observed in the simulations. For perturbation modes l > 1 the flow at the origin will always be zero, and so one can assume that in the strong anchoring limit the defect will remain centred at the origin. We again require an assumption for the r-dependence of the polarisation perturbation. Taking R c0 → 0, we can write a general form as δp ∝ r n for arbitrary n ≥ 0. Importantly, for all n, the active flows always give rise to an instability for ζ < 0 (contractile). Considering only active flows, the maximally unstable perturbation is for n = 0. Thus, below we consider only the results of this mode, which allows us to consider the dynamics in the limit where the filament polarisation at the interface and in the droplet are strongly coupled. However it comes at the cost of reducing the quantitative power of our predictions, and is an important restriction to the dynamics considered. Note, in 2-dimensions, the assumption n = 0 gives rise to an infinite passive contribution to the dynamics (proportional to K) and so we use n = 1, which appears consistent with what is observed in simultions.
In the strong anchoring limit, the polarisation has to match the perturbed interface normal at r = R to first order, such that r(∇Y m l (θ, φ)) . (10) We calculate the resulting flows to first order in δR. Since p is enslaved to the deformation we then only need to consider the radius dynamics given byṘ (for details see Supplementary Information appendix B).
In this strong anchoring limit we find that the droplet is unstable if ζ < α P < 0, i.e. the activity threshold, α P , is always contractile. The threshold α P increases linearly with γ and K. Repeating the linear stability analysis calculation in 2D shows the same qualitative prediction, where this time we take δp ∝ r as this is the leading order contribution allowed. The analytical expressions for the activity threshold are given in Supplementary Information appendix B and a full discussion of the eigenvalues of the general stability matrix (for weak anchoring) can be found in [45].
The result of this analysis is somewhat surprising, in this strong anchoring limit we expect the l = 1 mode to be unstable to extensile activity, whereas the higher modes of deformation are unstable for contractile activity. This suggests that, when our assumptions hold, we should see translational symmetry breaking with the defect moving to the droplet front for an extensile drop and symmetric modes of deformation for a contractile drop (see figure 4). This active threshold scales linearly with K and γ 0 , demonstrating the importance of the coupling of the morphology to the polarisation field. Contrast this to the case of the active interface where the shape does not affect the threshold for a phase transition.
This contractile instability can be understood physically by considering the splay in the drop due to perturbations in the interface curvature. High curvature couples to increased splay which couples to outward flow, further increasing the curvature of the interface and hence the splay. A sketch of this is given in figure 5.

C. Results and Comparison with Simulations
In the 2D simulations (see figure 6) we see symmetry breaking corresponding to the k = 1 mode for extensile activity resulting in a steady motile state, as predicted by the stability analysis. This is characterised by the defect centre moving to the front of the drop and is independent of the boundary deformation (and hence γ 0 ). Due to the extensile nature of the activity this droplet is a pusher, pushing fluid out along its axis of motion and thus elongating parallel to its motion.
Conversely contractile activity stabilises the defect at the droplet centre and we observe a k = 2 mode instability characterised by deformation of the droplet into a 'dumbbell' shape. It is also observed that this phase behaviour breaks down as the value of K/R 2 0 is reduced. In this limit the distortions in the droplet bulk are not strongly coupled to those at the interface and so more complex distortions can occur without significant droplet deformation. Our analytical calculations do not predict this as we assume a form for the r-dependence of the polarisation such that it is strongly coupled to the curvature. This behaviour goes beyond the scope of the analytical work here as this corresponds to a transition to an 'active turbulence' state, as numerically simulated in [29].
Finally, we also observe rotational steady states in the simulations (for extensile activity when using ν = −1.1) which can be characterised exactly by rotationally invariant distortions of the polarisation field [2,3], but these are not predicted for the parameter range used in figure 6.

IV. DISCUSSION
We have used analytical linear stability analysis and numerical simulation to characterise instabilities in active droplets and their resulting non-equilibrium steady states. Recent advances in experimental techniques mean that active gels of cytoskeletal material can be produced in vitro. The predictions of our active interface model could be tested by adsorbing an isotropic actin gel onto the interface of an emulsion drop containing myosin and ATP [33,34]. We predict an activity threshold for spontaneous motion, and a further continuous transition to a stable symmetric state mediated by advection of motors through the droplet bulk. We predict that in 3D this symmetric configuration will be coupled to deformation of the drop, however this cannot be observed in the 2D model.
The active polar drop model we use only predicts some of the dynamics of a real active polar drop system as it ignores the density and ordering magnitude degrees of freedom. However, this model system gives us an insight into the intrinsic instabilities when droplet deformation and filament polarisation direction are strongly coupled. In particular, there is a contractile activity threshold that is linearly dependent on surface tension, above which the droplet spontaneously deforms into a characteristic dumbbell shape. We also see persistent motility in the case of extensile activity such that the droplet acts as a pusher, compared to the puller type motion exhibited in the active isotropic interface model. This is consistent with previous active droplet models that show contractile activity resulting in droplets which are pullers and extensile activity resulting in pushers [13,15,18,20,21]. An interesting future extension of this work would be to consider coupling between both of the active phases studied here within a single drop.
The finite active systems we study improve our understanding of how confinement and deformation affect steady state dynamics. Additionally, we see the importance of feedback, driven by advection through the droplet or the internal orientational order, resulting in more complex dynamics. These results should prove useful in characterising future experiments on in vitro cytoskeletal networks and be useful in developing more complex models of multicomponent active systems in nature. The generic perturbation to the active interface model is of the form: where we assume the perturbations δg l,m 1 for g = c, v (i) , and R. Note that all real perturbations will satisfy the constraint δg lm = (−1) m δg * l,−m . From equation (??) we can calculate the perturbed normal and tangential vectors to first order in δR lm :n Substituting (??) and (??) into the expression for the force on the interface gives: where γ 0 = γ 0 −ζc 0 − Bc 2 0 /2 andζ = ζ + Bc 0 as defined in the text. As this force is located at the interface, we can solve the fluid flow in the bulk with Lamb's general solutions for flow around a sphere, which defines the δv (i) lm and δP lm functions: We can then find the coefficients a  2n · η 0 u − η 1 u ext r=R = F lm + (P − P ext )n r=R (S.14) where u = (∇v + (∇v) T )/2 is the strain rate tensor and the superscript 'ext' denotes a property of the external fluid Ω 1 . Solving these for the flow and pressure from equations (??)-(??) we find: where η = [3η 0 + 4ηl(l + 2)] and η * = (η 1 − 2η 0 + 4l 2η ), andη = (η 0 + η 1 )/2 as in the text. In order to satisfy the binding and unbinding dynamics, we consider the limitρ = 0, which gives δρ lm (r, t) = δc lm k off R (1−l) 0 r l /(D b l +k on R 0 ). We then use these solutions in the dynamic equations (1) andṘ = V .n which can be written in matrix form as: where X = −Dl(l + 1)/R 2 0 − D b k off l/(D b l + k on R 0 ) and Y = l(l + 1)/(η η * ). Clearly, when l = 1, there is no deformation and the eigenvalue λ of T is just equal to T 11 . In general the solution for the eigenvalues is the root of a quadratic. One of the eigenvalues is always zero when γ 0 = 0, and is positive for γ 0 < 0. This instability corresponds to a droplet with a negative effective surface tension. The other solution for λ = 0 corresponds tõ as given in the text. We see that the expression for λ in terms of l is complicated. However, since α I does not depend on the bare surface tension γ 0 the instability threshold does not depend on the deformation. Thus, we can find an approximate expression for λ by solvingṘ = 0 for δR lm , which gives: .
Substituting this for δR lm in equation (1) we find: The resulting equation for ∂λ * /∂l = 0 is a fifth order polynomial in l, this simplifies to a cubic equation in the limit k off → 0. The solution to this cubic equation is the analytical line plotted in figure 2. As shown in equation (??), the activity threshold does not increase linearly with l if k off = 0.
In 2-dimensions, the general solutions for flow around a circle are: where r and θ are plane polar coordinates.
Applying the perturbation to the concentration and solving the boundary conditions of continuous velocity and force balance at the interface we can solve for the constants a On substitution of these into the dynamic equations for c and R, we see that the deformation does not depend on δc k , and hence the droplet shape is always stable as a circle for positive effective surface tension (γ 0 > 0). This means that the only instability comes from the dynamics of c on the surface: leading to the active threshold in equation (6). The maximally unstable mode can again be found by solving for the zeros in the gradient of the dispersion relation ∂λ 2D /∂k = 0. The expressions for k max are the solutions to the cubic equation: When binding is omitted (k off = 0) this reduces to a simple monotonic relation in the activity: to be central. Using the expression for the bulk free energy in equation (5), and substituting in the expression for the polarisation (equation (??)) we get: is the bulk free energy for the unperturbed polarisation field. Thus, the net increase in free energy from this perturbation (F bulk − F 0 ) can be written where = R c /R 0 , as given in equation (8) of the text. This is positive for all R c < R 0 and converges to a fixed value at small . Introducing a small activity, the active force for this polarisation field to leading order in δz: v ≈ ζδz 6R 0 (R 0 − R c0 )η 0 R 2 0 7η 0 + 3η 1 3η 0 + 2η 1 − 3rR 0 + 3r 2 η 0 + η 1 3η 0 + 2η 1 cos(θ)r (S.40) − R 2 0 7η 0 + 3η 1 3η 0 + 2η 1 − 9 2 rR 0 + 6r 2 η 0 + η 1 3η 0 + 2η 1 sin(θ)θ (S.41) The initial stability of the defect position will be determined by whether the defect is advected along its perturbation direction, relative to the velocity of the drop. The active part of this velocity is calculated by the integral over the defect (1/V ) S (v act ·n)rr 2 sin(θ) dθdφ at r = R c , where V is the defect volume enclosed by the surface S, this gives: v c ≈ ζδz η 0 (7 − 9 + 3 2 ) + 3η 1 (1 − ) 2 6η 0 (1 − )(3η 0 + 2η 1 )ẑ (S.42) where = R c0 /R 0 . The velocity of the interface is the result of the same integral over the surface Σ at r = R 0 : v drop ≈ ζδz 6(1 − )(3η 0 + 2η 1 )ẑ (S.43) Finally, the relative velocity of the defect is given by v c − v drop : v c − v drop ≈ ζδz η 0 (2 − ) + η 1 (1 − ) 2η 0 (3η 0 + 2η 1 )ẑ (S. 44) as in the main text in equation (8).
This result can be reproduced in 2-dimensions also. The form of the defect displacement is identical to the 3D case in plane polar coordinates (r, θ) (except the perturbation is along x using this convention). Then the polarisation is: The resulting defect speed and free energy change calculated as before are as follows: (S.47) These predict similar behaviour as the 3D case, again with an extensile ζ > 0 instability, and an increase in total bulk free energy.
In order to find this in 2D we work through exactly the same set of operations in polar coordinates, except that we assume that the polarisation scales as p θ ∝ r, as using r 0 results in infinite passive contributions. The resulting activity threshold is plotted in figure 5 for the simulation parameters using the k = 2 mode.