On the Stability of Rotating Drops

We consider the equilibrium and stability of rotating axisymmetric fluid drops by appealing to a variational principle that characterizes the equilibria as stationary states of a functional containing surface energy and rotational energy contributions, augmented by a volume constraint. The linear stability of a drop is determined by solving the eigenvalue problem associated with the second variation of the energy functional. We compute equilibria corresponding to both oblate and prolate shapes, as well as toroidal shapes, and track their evolution with rotation rate. The stability results are obtained for two cases: (i) a prescribed rotational rate of the system (“driven drops”), or (ii) a prescribed angular momentum (“isolated drops”). For families of axisymmetric drops instabilities may occur for either axisymmetric or non-axisymmetric perturbations; the latter correspond to bifurcation points where non-axisymmetric shapes are possible. We employ an angle-arc length formulation of the problem which allows the computation of equilibrium shapes that are not single-valued in spherical coordinates. We are able to illustrate the transition from spheroidal drops with a strong indentation on the rotation axis to toroidal drops that do not extend to the rotation axis. Toroidal drops with a large aspect ratio (major radius to minor radius) are subject to azimuthal instabilities with higher mode numbers that are analogous to the Rayleigh instability of a cylindrical interface. Prolate spheroidal shapes occur if a drop of lower density rotates within a denser medium; these drops appear to be linearly stable. This work is motivated by recent investigations of toroidal tissue clusters that are observed to climb conical obstacles after self-assembly [Nurse et al., Journal of Applied Mechanics 79 (2012) 051013].


Introduction
Analyses of the dynamics of a rotating liquid drop held together by surface tension were initiated by Plateau [1]. In his work a liquid drop was immersed in an immiscible liquid which has about the same density but a much smaller viscosity than the drop. The drop was then spun at a controllable rate using a rotating rod that threaded the drop axis. Spinning of the drop produces significant deviations from the spherical equilibrium shape that is obtained for stationary drops. By matching the density of the drop and its surrounding medium, gravitational effects can be minimized in a terrestrial experiment. Assuming rigid body motion and taking into account solely axisymmetric drop configurations, the drops evolve from spherical configurations at zero rotation rate through a family of axisymmetric shapes that progressively flatten at the poles while developing an equatorial bulge. At large enough rotation rates, Plateau observed transient toroidal configurations that tended to break up into smaller drops (see also [2]).
Studies of rotating drops with significant density contrasts have also been performed in microgravity environments [3]. Due to the qualitative similarity between Plateau's rotating liquid masses held together by surface tension and self-gravitating celestial bodies, there have been many theoretical studies of the resulting equilibrium configurations and their stability, including work by Rayleigh [4], Maclaurin [5], Lyttleton [6], Chandrasekhar [7], Ross [8] and Brown and Scriven [9], and Myshkis et al. [10].
Recently there has been interest in the stability of toroidal shapes. Experimentally, macroscopic liquid toroidal droplets [11] and nanoscale toroids of varying sizes have been carefully generated and observed [12,13]. Analysis indicates that the stability of these toroidal shapes is related to the aspect ratio of the http://dx.doi.org/10.6028/jres.120.007 major and minor radii. Both groups report that toroids with a small aspect ratio tend to coalesce to form a single spherical droplet, while thinner toroids, i.e., those with large aspect ratio, mainly break up into smaller droplets.
Experiments on neonatal fibroblast cells that have self assembled into a toroidal cluster about the base of a conical pillar, have shown that the toroidal cluster will actively do work to climb the pillar to become a sphere or will remain at the base of the pillar and break up to form smaller clusters [14]. Subsequent theoretical work on this self assembled system points to the surface energy as the configurational driving force for the climbing motion of the cluster [15]. This suggests that the fate of the cluster is determined by its size. As such, thinner toroidal clusters do not climb the conical pillar and the development of localized deformations along its circumference is considered to be a result of the unstable growth of surface fluctuations.
These studies have revived interest in the stability of a rotating toroidal drop held together by surface tension. In this paper the stationary points of an energy functional are used to determine the equilibrium shapes of spheroidal and toroidal drops. The stability of the drops is then determined by examining the second variation of their energy functionals. The loss of stability of an equilibrium drop can indicate a transition to another equilibrium shape or signal the impending loss of the existence of the equilibrium drop itself.
The energy functional that we employ takes the form of a Lagrangian representing the difference between a drop's potential energy, due to capillary forces, and the kinetic energy of rotation, subject to a volume constraint. This formulation is analogous in many ways to that for the classical problem for the equilibrium and stability of rotating, self-gravitating drops, where the potential energy is then due to a Newtonian gravitational potential [6,16]. In that case there is a class of equilibria given by axisymmetric ellipsoids that take the form of oblate spheroids (shaped like "hamburger buns" with respect to the rotational axis). In our case there also are axisymmetric solutions that resemble oblate ellipsoids. These occur when the density of the drop exceeds that of the surrounding medium, so that due to centripetal acceleration the drop tends to bulge at the equator and flatten at the pole. On the other hand if the density of the surrounding medium exceeds that of the drop a family of solutions that instead resemble prolate ellipsoids ("cigar-shaped" with respect to the rotational axis) are possible. We will refer to these families of solutions as oblate spheroids and prolate spheroids, respectively, with the understanding that in our case these solutions are not literally axisymmetric ellipsoids of revolution. More generally we will refer to our solutions as spheroidal or toroidal depending on their topology.
For the drop undergoing rigid body rotation with a fluid of a different density, the application of an external torque may be necessary to maintain a given rotation rate as the moment of inertia of the drop changes due to variations in the shape of the drop. Here, such a drop is referred to as a driven drop. Alternatively, if there is no applied external torque, the angular momentum of the drop is conserved. If the drop rotates with constant angular momentum, it is termed an isolated drop. The energy functional for the driven drop is a function of the drop's shape and the rate of rotation, while the energy functional of the isolated drop is a function of the drop's shape and its angular momentum. This latter energy functional is formulated by a Legendre transform of the former functional to account for the constant angular momentum constraint, forming what is termed the Routhian in classical mechanics [17].
We first describe the variational formulation of the problem, including the Euler equation governing the equilibrium drops that results from the first variation of the energy functional. In § 3 we discuss the second variation that governs the stability of the equilibrium drops. Our numerical techniques are described in § 4, followed by numerical results in § 5. Some conclusions are given in § 6.

Variational Formulation
We describe the variational formulation for the equilibrium and stability of rotating spheroidal and toroidal drops. The case of a rigidly rotating system consisting of a drop confined by a co-rotating surrounding medium is considered first, followed by the modifications necessary to treat the rotation of an isolated drop that conserves its angular momentum. An example of the former would be a two-fluid system http://dx.doi.org/10.6028/jres.120.007 in a container that is attached to a platter revolving at a constant rate. An example of the latter would be a isolated drop rotating in vacuum in a microgravity environment.

Forced Rotation at a Prescribed Rate
The equilibrium of a driven axisymmetric drop that undergoes rigid body rotation at a specified angular velocity in tandem with a surrounding co-rotating medium can be described in terms of the stationary points of an energy functional that includes contributions from the kinetic energy and surface energy of the system. For simplicity we begin by formulating the variational principle in a cylindrical coordinate system in which the drop surface has the form = ( ) z f r for 0 1 < < r r r and > 0 z ; we will assume that the drops have a mid-plane of symmetry about =0 z . We later generalize to a body-fitted set of coordinates employing angle/arclength variables that avoids difficulties associated with infinite slopes and is more suitable for the stability determination of distorted drops. For a spheroidal drop 0 =0 r corresponds to the axis of rotation and 1 r is the equatorial radius, with a vanishing slope =/  The effective energy functional is written as where γ is the surface energy of the drop, Ω is the given rotation rate, and is the total surface area of the drop. The effective kinetic energy of the system is 2 /2 Ω  , where . We are assuming there are no gravitational effects. Equilibrium of the drop is then described by requiring the energy functional to be stationary to perturbations f δ in the shape that conserve the volume and satisfy appropriate boundary conditions at 0 = r r and 1 = r r , which leads to the Euler equation for the shape ( ) f r and the Lagrange multiplier P , which is chosen so that the volume constraint is satisfied. The Euler equation is equivalent to the Laplace-Young boundary condition inner outer = K p p γ − at a fluid-fluid interface, where K is the mean curvature of the interface and inner p and outer p are the local pressures on the inside and outside of the drop, respectively. In a hydrodynamic description of the motion based on the Navier-Stokes equations [18], in each phase the pressure satisfies a radial momentum balance 2 =/ r dp dr r Ω for a rigid body motion rΩ in the azimuthal direction, which integrates to 2  There is an analogous expression relating C and the outer radius 1 r where r f tends to negative infinity. The existence of toroidal solutions was proved in 1984 by Gulliver [19]. The spheroidal and toroidal solutions will be described in more detail below when we discuss the stability results in § 5.

Free Rotation of an Isolated Drop
If an isolated drop is freely rotating rather than being driven by an external torque that provides a constant rotation rate, it is appropriate to formulate the problem in terms of the drop's angular momentum L , which is conserved by the motion. With the explicit form for the energy functional [ , ] f Ω  given in Eq. (1), the angular momentum functional [ , ] L f Ω is given formally by We define the Routhian functional [ , ] f L  (see, e.g., [6,9]) via a Legendre transformation [20] with respect to Ω , [ , ] , where in the right hand side the rotation rate is now regarded as a functional [ , ] = / [ ] f L L f Ω  that is obtained by inversion of the relation (10). This leads to the expression whose first variation, taken at constant L , leads to the same Euler equation (6) as for the driven drop, since we have = L Ω  in each case. Thus the equilibrium states for the driven drops and the isolated drops are the same, although their stability differs, as we describe next.

Second Variation
We next describe the stability of the drops in terms of the second variation of their energy functionals. The equations for the second variation involve the equilibrium shape and its spatial derivatives, and it is convenient to first re-express the unperturbed shape in terms of angle/arclength variables to obtain a more tractable version of the stability equations. http://dx.doi.org/10.6028/jres.120.007

Angle/Arclength Coordinates on the Drop
The axisymmetric equilibrium shapes can be parametrized in terms of their arclength s as = ( ) r R s and = ( ) where the mean curvature is given by = sin( )/ s K r ψ ψ + . We note that the first integral (8) can be written in the form 2 3 sin = . 2 8 In the spheroidal case with =0 C this expression can be used to eliminate the singular term sin( )/R ψ in Eq. (14) to give the alternate expression On the other hand, since for the toroidal drop 0 (0) = > 0 R r this singularity does not arise in that case and Eq. (14) can be used as written. In either case the total volume is given by

Body-Fitted Coordinate System
To compute the second variation we employ a body-fitted coordinate system ( , , ) s w θ where s is arclength, θ is the azimuthal angle about the rotation axis, and w is distance measured along the local outward normal to the drop surface. The mapping from ( , , ) s w θ to cylindrical coordinates ( , , ) r z θ is then given by where the outward normal has components ( , ) = ( sin , cos )   To compute the first and second variations of the energy functional, we formally expand the interface displacement in terms of a small parameter ε , The various functionals can written in terms of ( , , ) , and similarly expanded as

Driven Drop
The first variation of the energy functional for the driven drop is then (1) and requiring (1) =0  for arbitrary variations (1) W recovers the Euler equation (14). Similarly, the second variation (2)  is given by http://dx.doi.org/10.6028/jres.120.007 The latter integrand proportional to (2) W vanishes because of the Euler equation, so that (2)  is a quadratic functional in the variation (1) W . The stability of the system is then determined by the sign of the second variation. This sign is most easily determined by diagonalizing the quadratic form, subject to the constraint The diagonalization of the quadratic form is equivalent to an eigenvalue problem, which leads to the Sturm-Liouville equation where λ is the eigenvalue and µ is a Lagrange multiplier for the volume constraint in Eq. (29), which must also be enforced. Here we have again used the Euler equation to simplify the final expression; note that P is absent from this equation. Expanding a general perturbation in terms of the orthonormal eigenmodes ( ( , ), ) then produces the diagonalized expression Here j a is the expansion coefficient of (1) W with respect to the eigenmode j W , Stability of the driven drop is obtained, viz. (2) > 0  , if j λ is positive for all eigenmodes, and instability occurs if any eigenvalue j λ is negative.

Isolated Drop
For the isolated drop, the expansion of the Routhian proceeds in a similar fashion, leading to the second variation which differs from the expression for (2)  by the latter positive term proportional to ( is the unperturbed moment of inertia. In this case diagonalizing (2)  leads to an integro-differential eigenvalue problem, http://dx.doi.org/10.6028/jres.120.007 The isolated drop is stable if the eigenvalues j λ are positive for all eigenfunctions of Eq. (34). Since the coefficients of both Eqs. (30) and (34) are independent of the azimuthal angle θ , normal modes of the form (1) (1) ( , ) = ( ) cos W s w s n θ θ are solutions, which reduces Eq. (30) to an ordinary differential equation (ODE) for the corresponding eigenmodes (1) ( ) w s . The non-axisymmetric modes with 0 n ≠ then automatically satisfy the volume constraint since the integrals over θ vanish, and the associated Lagrange multiplier µ can be taken to vanish. The volume constraint must still be applied for the axisymmetric modes with =0 n . Similarly Eq. (34) becomes an integro-ordinary-differential equation that must be solved along with a volume constraint for axisymmetric modes. The equation for nonaxisymmetric modes with 0 n ≠ also reduces to an unconstrained ODE, and the stability problem for nonaxisymmetric modes are identical for driven and isolated drops.
The normal mode solutions can generally be divided into families that are even or odd about a midplane of symmetry at =0 z . As discussed by Brown and Scriven [9], for the related cases of rotating drops that are held together by self-gravitation it is known that the drops have reflective symmetry about their equatorial plane [21]. Brown and Scriven therefore confined their finite element computations to solutions with even symmetry about =0 z . We have computed normal modes with either even or odd symmetry, and have found instabilities only for modes with reflective symmetry about =0 z . There are also neutrallystable modes with odd symmetry about the equatorial plane that correspond to simple energy-preserving translations along the z -axis. With the exception of these translation modes, the modes with odd symmetry about =0 z are found to be stable, and so we will confine our discussion to the even modes. The numerical procedures used to solve these corresponding eigenvalue problems are summarized next.

Numerical Techniques
The eigenvalue problems for determining the stability of the rotating drops are intractable analytically except in special cases, and we have resorted to numerical techniques for their solution. We have used two complementary approaches. Firstly, a finite difference discretization of the Sturm-Liouville equations can be used to produce a matrix eigenvalue problem, which produces N approximate eigenvalues for a system using N mesh points. Secondly, we have used an ODE solver in tandem with a shooting procedure to compute individual eigenmodes. The latter procedure is quite accurate provided adequate starting values are available to estimate the eigenvalues; we have used the matrix approximation to furnish the necessary initial guesses. Since the coefficients of the ODE's involve the shape of the unperturbed drop, the Euler equations are also solved numerically to provide the appropriate values at the mesh points of the matrix formulation or at the internal integration steps of the ODE solver. Some details are provided in Appendix 7.2 and 7.3.

Numerical Results
We first consider the case of heavier drops rotating in a lighter medium, followed by the case of drops that are lighter than their surroundings. http://dx.doi.org/10.6028/jres.120.007

Base States for
> 0 r ∆ The evolution of the axisymmetric drop shapes with > 0 r ∆ as the rate of rotation Ω is varied has been described by a number of authors [7,9]; some examples are illustrated in Fig. 2. Here we have defined the dimensionless rotation rate * Ω , the moment of inertia *  , and the angular momentum * L via where 0 R is the radius of the sphere with equivalent volume, . For small rates of rotation the drops are nearly spherical, and as the rate of rotation increases the drops develop an equatorial bulge while flattening at the poles. The continual decrease in polar radius eventually produces dimpling of the surface at the pole and the drop becomes non-convex. The family of spheroidal drops terminates at a point in parameter space where the polar radius 0 Z of the drop vanishes and the drop pinches off at the poles. There is also a nearby family of toroidal drops which originate near this point in parameter space; in this case, the inner radius 0 r of the torus tends to zero as the "hole" of the toroid closes up. The pinching of the spheroidal drops and the "healing" of the toroidal hole are illustrated in the sequences shown in Fig. 2. The relation between angular rotation rate and angular momentum for the spheroidal and toroidal drop families is shown in Fig. 3. The angular momentum L of the spheroidal drops initially increases with rotation rate, but the rotation rate eventually decreases as the angular momentum continues to increase (see Fig. 3). As the inner radius of the toroids increases, the angular momentum of the drops initially decreases, then reverses and increases steadily as the cross section of the toroids becomes more and more circular. We note that the spheroidal and toroidal families do not merge with one another, although their curves are very close at their respective terminal points in Fig. 3

Linear Stability of the Oblate Spheroids
The stability of rotating oblate spheroidal drops has been considered previously by a number of authors, including Chandrasekhar [7], Brown and Scriven [9], and Heine [22]; the latter two papers include the computation of non-axisymmetric solutions that bifurcate from the axisymmetric family at specific rotation rates. To validate our numerical procedure, we have reproduced these bifurcation points, which correspond to conditions of marginal stability ( =0 λ ) where the energy functional ceases to be a minimum relative to non-axisymmetric perturbations of a given mode number n . Results are shown in Fig. 4, where in the upper plot the bifurcation points for perturbations with mode numbers = 2,3, 4 n and 5 are indicated on the curve of * Ω versus * L . The numerical values agree with those given by Brown and Scriven for = 2,3, n and 4 to three decimal places; they were unable to compute the =5 n bifurcation because of their use of spherical coordinates, which preclude the computation of highly-dimpled shapes that are nonconvex relative to the origin. The results for =2 n also agree with those given by Heine to five decimal places; Heine also describes a bifurcation to an =6 n perturbation on the extension of the solution curves to self-intersecting drops for * > 1.1194 L . In the middle plot of Fig. 4  As discussed by Brown and Scriven, the stability of the driven drops ( * = Ω constant) and isolated drops ( * = L constant) to non-axisymmetric perturbations ( > 0 n ) are identical, since the perturbed moment of inertia vanishes for non-axisymmetric perturbations. For the case of axisymmetric perturbations ( =0 n ), the stability results do differ, as shown in the lower plot in Fig. 4. The driven drop is unstable to an axisymmetric disturbance at the limit point [23] of the solution branch where * Ω reaches its maximum value, with * = 0.7540 Ω , * = 0.7291 L . The isolated drop is stable with positive values of λ that, for each value of * L , are larger than those for the driven drop, as expected from Eq. (33). No limit point with respect to * L occurs on the solution branch, so there is no analogous axisymmetric instability for the isolated drop.

Linear Stability of the Toroids
The linear stability of the toroidal family of rotating drops is shown in Fig. 5. The upper plot gives the parametric relation between the rotation rate * Ω (solid curve) and the angular momentum * L (dashed curve) of the base state as functions of the dimensionless inner radius 0 r . As the inner radius tends to zero, the rotation rate decreases (over a short interval), and the angular momentum increases. The opposite is true as the inner radius becomes large, and there is a maximum value of * Ω and a minimum value of * L as 0 r varies from zero to infinity. The middle plot shows the lowest eigenvalues for =1 n to =5 n for 0 0 < < 2.5 r . For small values of 0 r toroidal drops are unstable for all five mode numbers, but as 0 r increases the stability of these modes increase and reach maxima near 0 = 0.5 r , where only the first three modes are unstable. With further increases of inner radius, the modes are all destabilized and the toroidal drop is unstable to higher and higher mode numbers; the trends indicated for the lowest five modes are also observed for higher mode numbers n and larger inner radii 0 r . The lower plot in Fig. 5 shows the lowest two eigenvalues for axisymmetric modes ( = 0) n for the case of driven toroidal drops (solid curves) and isolated drops (dashed curves). As expected from the previous discussion of Eq. (33), the isolated drops are more stable than the driven drops in each case, although for large values of 0 r the eigenvalues become nearly identical. For small values of 0 r the difference is more pronounced. The lowest eigenvalues both become very large and negative as 0 r tends to zero, but the lowest mode for the driven drop remains slightly unstable for large 0 r , whereas the isolated drop becomes stable near 0 = 0.5 r and then deceases in magnitude for large 0 r . The second lowest modes are both stabilized with increasing 0 r , although the driven drop is initially unstable for small 0 r . There are two axisymmetric, neutrally-stable modes ( =0 λ ). For the driven drop the neutral mode corresponds to a limit point on the solution branch in which * L is regarded as a function of * Ω . For the isolated drop the neutral mode corresponds to a limit point on the solution branch in which * Ω is regarded as a function of * L . In the top plot in Fig. 5 these points corresponds to extremal values of * Ω and * L regarded as functions of 0 r . These results indicate that the family of rotating driven toroidal drops is entirely unstable, both to axisymmetric and non-axisymmetric disturbances, with an increasing number of non-axisymmetric instabilities with increasing 0 r . The same is true for non-axisymmetric disturbances to the isolated drop, although in that case axisymmetric perturbations are stable for large enough values of 0 r . The geometry of the unstable axisymmetric modes is shown in Fig. 6 for a driven drop with 0 = 0.2 r , where the perturbed shapes for the lowest two modes are shown superimposed upon the base state (solid dots). The lowest mode (solid curve) represents a distortion of the shape that occurs predominantly at small radii, leaving the outer portion of the drop unaffected. Loosely speaking, this mode represents an instability driven by a change in the major radius of the torus. The second lowest mode (dashed curve) represents a perturbation that changes the ellipticity of the cross-section, with distortions at the inner and outer radii that are accompanied by a distortion of opposite sign at intermediate radii that preserves the net volume.

Rayleigh Instability Analogy
With increasing values of angular momentum * L drops are subject to an increasing number of nonaxisymmetric instabilities; only the most dangerous modes in the first five families ( = 1, 2,3, 4,5 n ) of instabilities are shown in the middle plot of Fig. 5. An example of a non-axisymmetric instability is shown in Fig. 7, corresponding to the neutral instability for =5 n that occurs near 0 = 1.2 r in Fig. 5. For toroidal drops with a large major radius, an interpretation of these high-wavenumber modes is possible in terms of a classical surface-tension-driven instability. For example for large values of * L the cross-section of the drops become more and more circular, and the drops increasingly resemble a circular torus. For large values of the effective major radius of the drop, the non-axisymmetric instabilities are then analogous to the capillary-driven Rayleigh instabilities [24]  ) the effective centrifugal force at the equator is inward, and the drops are elongated at the poles rather than the equator; we designate the resulting shapes as prolate spheroids. A dimensionless rotation rate P Ω for the prolate solutions is then defined as We consider only the case of driven drops. An analytic solution in this case was derived by Rosenthal [33] and Princen [34] in terms of incomplete elliptic integrals and is summarized in § 7.1. We have also implemented the previously-described numerical procedure for the base state in this case as well in order to facilitate the stability calculations.
In Fig. 8 we show the evolution of the prolate spheroidal shapes as the rotation rate P Ω is increased. For =0 P Ω the equilibrium is a spherical drop, and with increasing P Ω the equilibria tend to approximately cylindrical shapes that are terminated by roughly spherical caps. The equatorial radius E r decreases monotonically and the polar radius P z increases monotonically with increasing rotation rate, consistent with the imposed constraint of equal volumes for the family. Some numerical results are given in Table 2. For large rotation rates approximate expressions for the equatorial radius ( A E r ) and the polar radius ( tip A z ) can be obtained from an asymptotic evaluation of the elliptic integrals (see § 7.1); the corresponding results are also given in Table 2. Numerical calculations for the linear stability of the prolate drops are shown in Fig. 9. For both axisymmetric perturbations (lower plot) and non-axisymmetric perturbations (middle plot) the drops are found to be stable ( 0 λ > ). The stationary drop with =0 P Ω is again neutrally stable ( =0 λ ) to an =1 n mode that represents a lateral translation of the drop. For very large rotation rates the most dangerous axisymmetric mode is becoming decreasingly stable; the other modes are apparently increasingly stable for increasing rotation rates. http://dx.doi.org/10.6028/jres.120.007 As the rotation rate increases the drops become quite elongated with cylindrical mid-sections; it is therefore interesting to consider the possibility of a Rayleigh instability to axisymmetric perturbations with suitable wavelengths. We note that rotation about the cylindrical axis is known to stabilize the Rayleigh instability of an infinite cylinder if < 0 r ∆ . For example, Gillis and Kaufman [25] show that the rotating cylinder is stable if where C R is the cylinder radius, and k and n are the axial and azimuthal wavenumbers. Axisymmetric modes ( =0 n ) are stable for disturbances with 2 2 where we have expressed the eigenmode as (

Discussion
We have computed solutions for axisymmetric equilibrium shapes of spheroidal and toroidal drops or bubbles that correspond to extrema of an energy functional containing surface energy and rotational energy contributions, subject to a volume constraint. Examination of the second variation of the energy functional then determines whether the drops are stable, representing energy minima, or instead represent unstable saddle points or energy maxima. An alternate approach is to determine the linear stability of equilibrium shapes by solving the hydrodynamic equations of motion as given by Newton's law, which provides a dynamical growth rate for normal mode solutions. For example, Pairam and Fernández-Nieves ( [11], see also [26]) are able to interpret their experimental observations of the breakup of toroidal drops by comparing with the theoretical results of Tomotika [27] for the fastest-growing instability of a cylindrical thread of viscous liquid surrounded by another viscous fluid. A number of other authors have discussed the dynamic instability of toroidal drops based on approximate base states that are assumed to have circular cross sections [28,29], or have observed or simulated the temporal evolution of arbitrary (non-equilibrium) toroidal shapes [13,26,30,31]. Our approach focuses on the accurate computation of bifurcation points for self-consistent equilibrium shapes. As discussed by Brown and Scriven [9], the role of viscosity in determining the linear stability of rotating drops by solving the hydrodynamic governing equations can lead to subtle distinctions between "ordinary stability" and "secular instability," wherein an equilibrium that is stable according to the inviscid equations of motion is destabilized by the inclusion of viscous effects [6,32]. The stability results that we compute based on energy minimization correspond to the viscous case in this context.

Elliptic Integral Formulation for Prolate Spheroids
For spheroidal solutions, the explicit first integral of the Euler equation (8) has a vanishing constant of integration. The oblate spheroidal solution, for > 0 r ∆ , was obtained in terms of elliptic integrals by Chandrasekhar [7], and for the prolate spheroidal solution, for < 0 r ∆ , by Rosenthal [33] and Princen [34]. Here we summarize the solution for the prolate spheroid to obtain the asymptotic solution used in Table 2.
Evaluating Eq. (8)  Pr r where Σ is a dimensionless rotation rate used by Chandrasekhar [7], who then goes on to obtain solutions for > 0 Σ which correspond to the oblate spheroid. Here we consider the range 1/2 < < 0 − Σ , corresponding to solutions for the prolate spheroid. The dimensionless rotation rates Σ and P Ω are related by 3 Scaling by the equatorial radius, with 1 =/ r r r , the first integral can now be expressed as 3 In the integral above, the argument of the radical can be simplified as 2 ( For 1/2 < < 0 − Σ , the roots are all real with 1 < < b a . Therefore, Eq. (A4) can be rewritten as The volume of the prolate spheroid also has an analytical expression given by , then we may also write < , > = T 1 w h w as a conventional matrix inner product. The constraint equation is then treated by introducing an explicit projection onto the subspace < , > =0 1 w , defined by the matrix = , where  is the identity matrix and the outer product T 1h is a rank-one matrix.  satisfies =0 , and is symmetric with respect to the inner product: < , The diagonalization of (2) E on the subspace < , > =0 1 w is achieved by computing the N normal modes j w that satisfy the conventional symmetric eigenvalue problem = . We note that with initial conditions (0)   For the isolated non-axisymmetric drop, the integro-differential equation (34) is converted to an ODE by replacing the integral term by an auxiliary variable during the shooting procedure. We set 2  Since the leading order of the additional term is found to be 4 ( ) O s , retaining only one term is sufficient for this additional series, with Here B is the integral term Eq. (A29) from the integro-differential Sturm-Liouville equation and the coefficients 2 w and 4 w remain as in Eq. (A39) and Eq. (A40) with =0 n .