Shapes, stability, and hysteresis of rotating and charged axisymmetric drops in a vacuum

The behavior of rotating and/or charged drops is a classic problem in fluid mechanics with a multitude of industrial applications. Theoretical studies of such liquid drops have also provided important insights into fundamental physical processes across nuclear and astrophysical lengthscales. However, the full nonlinear dynamics of these drops are only just beginning to be uncovered by experiments. These nonlinear effects are manifest in the high sensitivity of the breakup mechanisms to small perturbations of the initial drop shape and in observations of hysteresis in the transition between different drop shape families. This paper investigates the equilibrium shapes and stability of charged and rotating drops in a vacuum with an energy minimization method applied to spheroidal shapes and with numerical simulations using a finite-difference, level-set method. A good working formula for the stability limit of these drops is given by Lmax = 1.15 − 0.59x − 0.56x2, where L is the dimensionless angular momentum and x is the charge fissility parameter. These methods also provide a firm explanation for the hysteresis of rotating and charged drops.


I. INTRODUCTION
Investigations into the shape and stability of liquid drops have a long and illustrious history with early work including Plateau's (1863) experiment on rotating drops 1 and Rayleigh's (1882) derivation of the stability of charged drops. 2 Indeed, the first studies of the stability of liquid surfaces under electrostatic forces can be traced back to before 1600 when Gilbert described 3 (p. 89) "the case of a spherical drop of water standing on a dry surface; for a piece of amber applied to it at a suitable distance pulls the nearest parts out of their position and draws it up into a cone." Remarkably, the mathematical description of these cones dates back to only 50 years ago 4 and their dynamic behavior remains an active area of research. 5 Liquid drop models have subsequently been used to model phenomena across a vast range of lengthscales: from the fission of atomic nuclei 6 to the shapes of astrophysical bodies. 7 An irresistible link between these lengthscales is made by noting the similarity of the electrostatic self-repulsion of a uniformly charged drop to the gravitational attraction of an astrophysical body; the equations describing these objects differ only by a scale factor and a single minus sign. 8 Liquid drop models continue to provide physical insights into, for example, the explosive chemical reactions of alkali metals and water 9 and the biophysics of living cells. 10 The hydrodynamics of drops is also vital for a range of engineering applications such as electrosprays for mass spectrometry, 11 film deposition, 12 inkjet printing, 13 and atomization of fuels. 14 Additional naturally occurring examples include the bursting of charged drops in thunderstorm clouds, a) Electronic mail: j.holgate14@imperial.ac.uk b) Electronic mail: m.coppins@imperial.ac.uk which is believed to be a trigger for lightning strikes, 15 and the shapes of rapidly spinning molten micrometeorites. 16,17 This work is motivated, in particular, by the behavior of drops in plasmas that are either formed deliberately, in processes such as plasma spraying, 18 vapour deposition, 19 and laserinduced water condensation, 20 or occur unintentionally due to, for example, erosion of drops from metal surfaces in plasma arc cutters 21 and magnetic fusion devices. 22 These drops can accumulate significant amounts of charge and rotation from the plasma and are subsequently susceptible to electrostatic 23 and rotational breakup; compelling evidence for the latter process has been seen in forked trajectories in camera images of metal drops in fusion plasmas such as those in Fig. 10 of Ref. 24.
Recent experiments have revealed intriguing behaviors of charged and rotating drops, such as hysteretic shape transitions and hybrid and asymmetric breakup processes, 25 which have not yet been examined by theoretical and computational methods. These experiments, and other previous studies of spinning and charged drops, are reviewed in Sec. II. Theoretical and computational methods for describing these drops are then described in Sec. III, followed by results for the shapes, stability, and hysteresis of charged and rotating drops in Sec. IV.

II. EXISTING WORK
The introduction of dimensionless parameters facilitates a straightforward comparison between existing theories and experiments and between drops on different lengthscales. The conventional definition of the dimensionless angular velocity, as introduced by Chandrasekhar, 26 is where ω is the dimensional angular velocity and ρ, r, and γ are the density, radius, and surface tension of the drop, respectively. The radius is defined as that of a sphere with the same volume as the drop. The dimensionless angular momentum is defined in terms of the dimensional momentum L and the volume of the drop V as using the convention of García-Garrido, Fontelos, and Kindelán 27 and Liao and Hill. 25 Note that this differs from the definition used by others such as Brown and Scriven. 28,29 The dimensionless parameter for charge, which is often referred to as the fissility parameter, is where q is the dimensional charge.
The first comprehensive study of the breakup of uncharged rotating drops is the finite element analysis of Brown and Scriven, 28,29 which has since been extended by Heine. 30 These calculations show that slowly spinning drops assume an axisymmetric flattened shape that is aligned with the axis of rotation. Such shapes are evident in, for example, the slight equatorial bulge of the Earth and can be well approximated by oblate spheroids. When the angular velocity of the drop exceeds a critical value Ω * , or equivalently the angular momentum exceeds L * , the drop begins to elongate in a direction perpendicular to the axis of rotation. The values of Ω * and L * have been determined analytically 26 and numerically 29 as 0.559 and 0.628, respectively, and confirmed by experiments. 25 The drop transitions through triaxial ellipsoidal shapes when Ω * is reached before assuming a shape that is well approximated by an elongated prolate spheroid. The moment of inertia of the drop is increased in the prolate configuration, so the angular velocity will decrease even as the angular momentum is increased. As such the angular velocity of stable drops does not exceed Ω * in numerical models, 29,30 but a slight overshoot has been observed experimentally. 25,31 If additional angular momentum continues to be supplied to the drop, then it will assume a two-lobed dumbbell shape that eventually splits into two large drops. Often these are joined by a thin thread of liquid that subsequently bursts into several smaller drops. This breakup process occurs at L max , which has been determined numerically 29 as 1.08 and experimentally 25 as 1.15.
The first study of the stability of charged drops is due to Lord Rayleigh 2 who predicted that a drop becomes unstable at the critical value x c = 1. The drop remains spherical until this value is reached at which point it elongates and forms a cone at each end of the drop. Jets of highly charged droplets are emitted from the tips of these cones that carry away the excess charge on the drop. This process, and the corresponding value of x c , was not confirmed experimentally until relatively recently. 32,33 Rayleigh arrived at this limit through a linear stability analysis, but Taylor predicted the identical result using a spheroidal approximation. 4 Taylor also calculated the opening angle of the cones formed during the breakup process; these are now known as Taylor cones.
The effect of rotation on the breakup of charged drops was first studied by Natarajan and Brown 34 who extended Rayleigh's analysis to include small deformations of the equilibrium spheroid shape caused by slow rotation. The first experimental study of charged and rotating drops appears to be that of Rhim, Chung, and Elleman 35 who obtained some limited quantitative data. García-Garrido, Fontelos, and Kindelán provided the first comprehensive theoretical study of these drops using a small-L perturbation theory, a spheroidal approximation, and numerical simulations. 27 Their methods are in excellent mutual agreement and reproduce the charge-and rotation-only results of previous studies, but their final bifurcation diagram, Fig. 9, does not cover the full range of angular momenta up to L max , and there appears to be some discrepancy between the values of L * in this bifurcation diagram and their earlier Fig. 8. Liao and Hill have very recently published experimental results for the shapes and stability limits of charged and rotating drops in a diamagnetic levitator for a full range of angular momenta and charges up to the Rayleigh limit. 25 One intriguing feature of experimental studies of rotating drops, both with and without charge, 25,31 is that the point where the drop transitions from an oblate to a prolate shape, or from a prolate to an oblate shape, depends on whether the drop is spinning up, i.e., its angular momentum is increasing, or spinning down. Drops that are spinning up will retain their oblate shape for angular velocities above Ω * , while drops that are spinning down retain their prolate shape at angular velocities below Ω * . The dependence of the shape of the drop on its history is a clear example of hysteresis that has not yet been examined theoretically or computationally.
The mathematical and numerical methods that are required to probe the behavior of charged and rotating drops are described in Sec. III, before moving on to the results generated by these methods. The results for the spinning-up/down hysteresis in Sec. IV D are the first of their kind and are of particular interest.

A. Spheroidal energy minimization
Spheroidal approximations have been used to estimate the shape and stability of liquid drops for a wide range of problems involving gravitational, 36 rotational, 37 electric, 4 and even nuclear 38 forces. The assumption that the drop has a spheroidal shape often leads to accurate results that are in excellent agreement with numerical simulations of arbitrarily shaped drops. The simplicity of the spheroidal analysis also provides insights into the physical processes involved in, for example, the hysteresis of the prolate-to-oblate transitions that are not immediately evident from full simulations. The theoretical foundations of the spheroidal approximation are briefly reviewed in this section. An energy minimisation procedure, similar to those of Cheng and Chaddock 39 and Mestel, 13 is used here rather than the pressure balance method of Taylor 4 and Basaran and Scriven. 40 The generalized coordinates of a spheroidal mass of liquid must be determined in order to perform an appropriate minimization of the spheroid's energy. These coordinates can be identified from the kinetic energy of the spheroid; the more general case of an ellipsoid with three distinct time-varying axes is briefly considered here. If the liquid drop rotates like a rigid body, i.e., all fluid elements within the drop have the same angular velocity, then it can be studied in a frame of reference co-rotating with the drop. The liquid is assumed to be incompressible and irrotational in this frame such that its velocity potential, χ, satisfies Laplace's equation. The velocity of the fluid is given by the gradient of this potential as v = −∇ χ. A simple derivation of χ is given by Lamb 37 (p. 147) who attributes it to an earlier study by Bjerknes. Inserting the equation of the ellipsoidal surface, into the surface condition Df /Dt = 0 yields which is satisfied by Laplace's equation is solved by this potential wheṅ which is simply the condition that the volume of the ellipsoid, V = 4πabc/3, is conserved. The kinetic energy follows as and the generalized coordinates are identified accordingly as The equilibrium shape of the liquid can now be determined by energy minimization according to Ref. 37 (p. 713), where U and T 0 are the potential energy and the rotational kinetic energy of the ellipsoid. Note that equilibria being determined are static in the frame of reference co-rotating with the drop, so dynamic effects, such as the viscosity or the nonrotational kinetic energy in Eq. (8), do not need to be included in this energy minimization procedure. The condition for the resulting equilibrium shape to be stable is that the energy is a minimum, so Each condition gives a set of equations that must be solved simultaneously. However, the volume conservation constraint reduces each set of equations to a single condition. If the unique axis of a spheroid is given by a, as illustrated in Fig. 1, then its volume conservation constraint is ab 2 = r 3 , where r is the radius of a sphere with equal volume to the spheroid. Minimizing the energy of the spheroid with respect to b, and applying the chain rule, gives FIG. 1. The unique axis of both oblate (left) and prolate (right) spheroids is taken as a in this work. The spheroid energy is minimized with respect to α = a/r, where r 3 = ab 2 , in order to determine its equilibrium shape and stability.
The second derivatives are analyzed in the same manner to give The first-order derivative is zero when the spheroid is in equilibrium, so the second-order derivatives with respect to a and b must have the same sign. It is therefore sufficient to perform the energy minimization procedure in terms of either a or b. The minimization procedure followed in this paper uses the elongation α = a/r, which simplifies the algebra of the energy terms. Therefore the equations to be solved for stable spheroidal drops are where expressions for U and T 0 are given as follows.
For prolate spheroids, the surface energy γS can be written in terms of the elongation where division by 2πr 2 γ renders the energy dimensionless. The equivalent expression for oblate spheroids is The energy minimization procedure is valid provided that the liquid drop rotates like a rigid body, i.e., all fluid elements within the drop have the same angular velocity. This configuration can be attained experimentally by slowly spinning the drop to allow the viscosity of the drop to damp out any shear stresses between fluid elements rotating with different angular velocities. The rotational kinetic energy, Iω 2 /2, can be expressed in terms of the elongation and the dimensionless angular velocity as T (p) 0 for prolate spheroids and as for oblate spheroids with Ω defined in Eq. (1). The electrostatic energy due to a charge q on a spheroid is given by q 2 /2C, where expressions for the capacitance C can be found in the work of Landau and Lifshitz. 41 This energy can be written as for prolate spheroids and as for oblate spheroids with x defined in Eq. (3). Note that, unlike a rotating drop, there is no preferred direction for the elongation or squashing of a charged drop other than its initial shape perturbation.

B. Finite-difference, level-set (FDLS) simulations
The method described in Subsection III A provides simple and accurate results when the drop shape is well approximated by that of a spheroid, but a complimentary method must be developed, which accounts for significantly non-spherical shapes. The finite-difference, level-set method (FDLS) has emerged as a popular scheme for simulating electrohydrodynamic flows over the past decade. [42][43][44][45][46][47] The popularity of these methods is attributable to their simple implementation and natural handling of complex configurations and merging or breakup of the interface. This section gives an overview of the method and its implementation in the code iEHD. This code utilizes cylindrical (r, z) coordinates so that axisymmetry is imposed on the simulations. The code and the results given in this paper are freely available at http://www.github.com/joshholgate/iEHD.

Ideal electrohydrodynamic (EHD) equations
The ideal EHD equations for a conducting liquid with viscosity η in contact with a vacuum are given by the incompressibility condition and Navier-Stokes equation, in the liquid region and by the Maxwell equations for the electric field, in the vacuum region. Equation (21) requires that the density ρ is uniform in the liquid region, which removes the need for an energy equation. The last term in Eq. (22) gives the centrifugal force due to rotation of the drop, at angular velocity ω, around an axis in theω direction; the strength of this force depends on the perpendicular distance r of each fluid element from the rotation axis. Equation (24) is valid provided that there are no externally applied magnetic fields or electrical currents flowing, and it allows the electric field to be written as the gradient of the electric potential, E = −∇φ, so that Eq. (23) becomes the Laplace equation, The electric potential takes a constant value at the surface of the conducting liquid and, at the outer edge of the simulation, is determined by a truncated multipole expansion, where q is the charge on the drop, R = (r 2 + z 2 ) 1/2 is the radial distance in spherical coordinates, and the quadrupole moments are estimated by those of a spheroid with the same charge and polar-to-equatorial axis ratio as the drop (Ref. 41, p. 28), The axisymmetric drops simulated in this paper have no dipole (φ ∼ 1/R 2 ) or hexapole (φ ∼ 1/R 4 ) moments, so the leadingorder omitted term is the octupole (φ ∼ 1/R 5 ) term. Taking the gradient of Eq. (26) gives the electric field which is enforced as a Neumann boundary condition at the upper and lower edges of the cylindrical domain and a similar condition, with r ↔ z, at the outer radial surface. The shortest distance from these boundaries to the surface of the drop is three times the radius of the drop. The greatest inaccuracy in these boundary conditions arises from assuming spheroid-like quadrupole moments. This assumption is valid for moderately deformed drops during the early development of the Rayleigh breakup process and, accordingly, it was found that extending the simulation domain did not modify the results for the stability of the drops. The electric field is aligned with the outward normal of the equipotential liquid surface. The liquid pressure at the surface, as found from the normal component of the stress continuity condition [Ref. 48,Eq. (3.9)], is simply given by where γ is the surface tension value and κ is the local curvature of the surface; note that κ is negative for a convex surface. Equation (30) provides a boundary condition for the Navier-Stokes equation that is evolved through time in order to simulate the fluid motion under electric, viscous, surface tension, and gravity forces. An illustration of the ideal EHD configuration is shown schematically in Fig. 2. The E 2 s term means that the EHD behavior is identical for positive or negative drop charges with the same magnitude.
The ideal EHD equations, Eqs. charged and rotating droplet. The simplifications assumed by these equations require some elucidation. First, replacing the vacuum region with a surrounding gas should not modify the EHD behavior of the drop in any tangible way; the low density of a gas compared to the liquid renders it "dynamically inactive," 5 and its permittivity is essentially that of the vacuum. However, when the drop is in a vacuum, evaporation could significantly modify the evolution of the drop. Indeed evaporation is often used experimentally to reduce the size of a drop, while keeping its charge fixed, until it reaches the Rayleigh limit and becomes unstable. 32,33 The disruption of the drop proceeds rapidly once this stability limit is reached. An estimate for the hydrodynamic time scale for deformations of a 1-cm water drop can be found from the speed of sound in water as ∼10 −5 s. Experiments of similarly sized water drops at atmospheric pressure require several hundreds of seconds for complete evaporation. 49 The evaporation rate scales linearly with pressure according to the Hertz-Knudsen relation, 50 so a vacuum pressure below 10 −7 atm is required before the evaporational volume losses become important on the time scales taken for the hydrodynamic deformation of the drop. Finally the omission of electrical breakdown and quantum emission effects can also become important due to strong enhancement of the electric field when the surface becomes highly deformed; these effects are discussed in Sec. IV A.

Level-set interface tracking
The position and motion of the liquid interface are captured using the level-set method. 51 This involves defining an additional scalar field ψ, known as the level-set function, which is positive inside the conducting fluid and negative outside it; the interface between the two regions is then given by the contour ψ = 0. A popular choice for ψ is the signed distance function, as used here to specify the initial surface configuration, which gives the closest distance from the interface at any point and which changes the sign across the interface. The level-set function is advected with the velocity field according to ∂ψ ∂t and the surface of the fluid moves with it. The advection of ψ requires that the velocity field v is given everywhere, but, so far, the velocity has only been defined in the liquid region. A fictional velocity, commonly referred to as the extension velocity, 52 is therefore introduced in the vacuum region that is used solely in Eq. (31) in order to avoid bunching and stretching of ψ in the vacuum region. The fictional velocity does not provide any additional mechanical stresses on the interface. The code described here constructs this fictional velocity by assuming that it satisfies the incompressibility condition and, at the liquid-vacuum interface, takes the same value as the adjacent liquid. Note that this choice of fictional velocity preserves the smoothness of ψ but does not ensure that it remains a signed distance function. Reinitialization algorithms, which are sometimes used to maintain the smoothness of ψ, are therefore not required. The level-set function allows straightforward calculation of the normal vector and curvature of the free surface using the formulaen = −∇ψ/|∇ψ| and κ = −∇ ·n. The disadvantage of the level-set method, compared to alternatives such as the volume-of-fluid method, is that the volume of the fluid tends not to be conserved; however, for the cases tested so far, the code described here conserves more than 99.9% of its initial volume. This better-than-expected volume conservation is most likely to be due to the choice of a divergence-free fictional velocity and the avoidance of any reinitialization algorithms.

Finite difference discretization
A finite difference method is chosen to discretize the ideal EHD equations because it is easily adaptable to different liquid configurations and because it will allow a simple extension of the code to describe plasma-liquid interactions in a future piece of work. The discretization of the Navier-Stokes equation is essentially given by the Chorin projection method on a uniform staggered grid, as described in Ref. 53. Provisional values of the liquid velocity at the next time step are calculated using the forward Euler method but with the ∇p terms excluded. These provisional velocities are then substituted into the incompressibility condition to produce a Poisson equation for the pressure at the next step. This is solved using a relaxation method, and the provisional velocities are then updated according to the calculated pressure values at the next step. This method is described as semi-implicit since the velocities are updated explicitly, while the pressures are updated implicitly. The temporal discretization is first order, but the code described here uses second-order upwind differences for the convective term (v · ∇)v in order to give second order accuracy to the spatial discretization. The convective term near the liquid-vacuum interface is calculated using the fictional velocity values in the vacuum region where necessary.
The Poisson equation for the liquid pressures requires the free-surface boundary values to be calculated in the vacuum gridpoints that are adjacent to the surface. The electric field in the vacuum region must therefore be calculated prior to evaluating the pressure Poisson equation. The Laplace equation for potential is solved using relaxation with the appropriate boundary conditions for the problem. No special considerations are made for the exact position of the liquid-vacuum interface, and the φ(ψ = 0) = 0 condition is simply applied to all gridpoints with ψ > 0. However, the calculation of the electric field strength accounts for the orientation of the surface by taking an off-centred second-order finite difference gradient of φ in the direction of the normal to the surface. The normal and curvature of the surface are calculated by taking second-order derivatives in the opposite direction to the normal so that values of ψ are taken from either side of the interface.
The final step is to update the level-set function before repeating this procedure at the next time step. This requires the fictional velocities in the vacuum region to be calculated. This is done by assuming the fictional velocity to be irrotational as well as incompressible to produce a Laplace equation for velocity potential. The velocity potential is found at the centre of gridpoints so that the velocities, which are placed at grid boundaries in the staggered grid configuration, may simply be taken from the difference in velocity potential of neighbouring cells. The fictional velocity is assumed to satisfy the same outer boundary conditions as the real liquid velocity, while the boundary values at the liquid-vacuum interface are given by the liquid velocity at that point. Finally the level-set advection equation is evaluated using the second-order upwind finite difference scheme.

IV. RESULTS
This section presents the results of the spheroid energy minimization and FDLS methods, as described in Sec. III, which are used to calculate the shape and stability of charged and/or rotating drops. The hysteresis of the oblate-to-prolate shape transition is then examined.

A. Stability of charged, nonrotating drops
The first problem to be solved is the disruption of a charged, nonrotating drop. The dimensionless potential energy comprising Eqs. (15) and (19) for prolate spheroids with α > 1, or Eqs. (16) and (20) for oblate spheroids with α < 1, is displayed in Fig. 3 for a range of values of x. The stable and unstable equilibria are easily identified as the minima and maxima, respectively, of the potential energy in accordance with the energy minimization method. It is straightforward to extract the locations of the extreme points from Fig. 3 and plot them in terms of x and α to produce the bifurcation diagram shown in Fig. 4. Two bifurcation points are immediately evident with the transcritical bifurcation at x = 1 corresponding exactly to Rayleigh's stability limit. The other bifurcation point, for significantly elongated spheroids with α = 1.75, is of little relevance to real liquid drops that will become unstable and highly nonspheroidal at such large deformations. The transcritical nature of the Rayleigh-limit bifurcation has been identified previously by Tsamopoulos, Akylas, and Brown, 54 and, subsequently, the family of oblate spheroid shapes at x > 1 was shown to be unstable to nonaxisymmetic perturbations by Natarajan and Brown. 34 Tsamopoulos noted that the existence of the unstable prolate spheroid shape family means that spherical drops charged below the Rayleigh limit are unstable to finite-amplitude perturbations, which helps one to explain the difficulties encountered in the experimental verification of the Rayleigh limit. 32,33 The nonlinear dynamics of the disruption of charged, nonrotating drops can be resolved by numerical simulations using the FDLS method described in Sec. III B. Figure 5 gives an example of the disruption of a drop charged to just above the Rayleigh limit with x = 1.02. The simulation parameters correspond to a water drop of diameter 1 cm. The drop is initiated with a slight prolate elongation of α = 1.02 in order to facilitate the instability. The symmetrical formation of two Taylor cones is observed, in accordance with expectations, followed by the ejection of jets of subdroplets.
It is worthwhile to pause at this juncture and consider whether the FDLS simulations are producing physically realistic results. Recent simulations of thin liquid films in vacuum electric fields have established that jets of droplets cannot be ejected from Taylor cones when the liquid involved is a perfect conductor. 5 This conclusion is supported by additional simulations of drops charged above the Rayleigh limit 55 and uncharged drops in applied electric fields. 56 Prior to these FIG. 5. Simulation results for the symmetric disruption of an overcharged nonrotating drop with charge x = 1.02 via the formation of two Taylor cones. Jets of droplets form at the tips of the cones in the final frame. The size of these droplets depends on the resolution of the numerical grid, which supports a previous study that concluded that they are not physical features of perfectly conducting liquid drops. 5,55,56 simulation studies, Zubarev had shown theoretically that a perfectly conducting Taylor cone is always stable when the ideal EHD equations are taken to fully describe the physics of the drop. 57 These previous studies suggest that the ejection of droplets in the FDLS simulations cannot be a physical EHD effect but is instead an artifact of the finite mesh size of the simulation. The size of these subdroplets decreases as the mesh spacing is reduced but does not disappear entirely even on a high resolution grid with spacings of r/200. Further evidence for the aphysical properties of these droplets is apparent from the discrepancy in the number of ejected droplets in the final frame of Fig. 5; this asymmetry has no physical origin. The ejection of aphysical droplets does not invalidate the results of the FDLS simulations but serves to illustrate their limitations for studying sub-gridscale phenomena. However, the results for the overall shape and stability of the drops should not be affected by these unrealistic small-scale inaccuracies that only occur at the apices of the cones.
The fact that perfectly conducting drops cannot disintegrate by ejecting jets of subdroplets prompts the question of what might instead happen at the tips of the drops. The high curvature of these tips causes a significant enhancement of the electric field from 2.5 MV m −1 for a spherical water drop of diameter 1 cm at the Rayleigh limit to a peak of 26 MV m −1 immediately before the ejection of droplets. This geometric enhancement of the electric field is highly localized, as shown by the electric field profile in Fig. 6, but can lead to physical effects beyond the ideal EHD model at the very tips of the drop. Evidence for electrical breakdown of the surrounding gas is observed in the experiments of Macky 58 62 Such applications use an electrically biased capillary tube to deliver the liquid metals, but the emission process is identical for isolated liquid drops. 63 Once the excess charge of an isolated drop has been emitted, either as electrons, ions, or charged droplets, then it will relax back to an equilibrium spherical shape.

B. Stability of uncharged, rotating drops
The disruption of uncharged, rotating drops is considered next. The effective potential energy, U − T 0 , is calculated using Eqs. (15) and (17) for prolate spheroids and Eqs. (16) and (18) for oblate spheroids. These energies are plotted in Fig. 7 for spheroids with a range of angular velocities and elongations. The transition between prolate and oblate shapes at α = 1 is no longer smooth because the axis of rotation flips from the unique axis of length a for the oblate spheroids to one of the degenerate axes of length b in the prolate case. The stable equilibria are, once again, determined by the locations of the minima, while the maxima are identified as unstable equilibria.
The bifurcation diagram for uncharged, rotating spheroids is displayed in Fig. 8. This diagram is easily produced from the α and Ω values of the energy minima and maxima. A single saddle-node bifurcation is identified from the bending-back of the curve at Ω = 0.53, which is close to the accepted value of the oblate-to-prolate transition point Ω * = 0.559. This critical point is variously referred to as a turning point or limit point value in bifurcation theory, but the term bendback value is used here. Stable prolate spheroid shapes cannot exist for values of Ω above this transition point. However, the stable oblate shape family extends to values of Ω beyond this transition point; this is in contradiction to the assertion of Brown and Scriven that "all axisymmetric shapes that spin more rapidly than Ω II (this paper's Ω * ) prove to be unstable to a two-lobed perturbation" 29 (p. 345). Drops that are having their angular momentum gradually increased can therefore remain in an FIG. 7. The potential energy of rotating and uncharged spheroidal liquid drops as a function of elongation and dimensionless angular velocity. The oblate spheroid rotates around its unique axis of length a, while the prolate spheroid rotates around one of its degenerate axes of length b; this flipping of the rotation axis leads to non-smooth solutions at α = 1. The prolate shape family loses stability at the saddle-node bifurcation indicated by the cross for Ω * = 0.53.
FIG. 8. Bifurcation diagram for rotating and uncharged spheroidal liquid drops. The prolate spheroid shape family transitions from stable to unstable equilibria at the saddle-node bifurcation point with Ω * and α = 1.4. The oblate shape family remains stable for Ω > Ω * . A pitchfork bifurcation is also identified at Ω = 0. oblate configuration above the Ω * transition point in accordance with experimental observations; 25,31 conversely drops that are having their angular momentum steadily decreased can assume stable prolate shapes below Ω * . The resulting hysteresis of the oblate-to-prolate transition is discussed in Sec. IV D.
The prolate spheroid shape family with α > 1.4 is unstable according to energy minimization arguments. However, real liquid drops do not immediately become unstable at large elongations and instead transition into two-lobed capsule and dumbbell shapes that cannot be captured by the spheroid approximation. Instead, these shapes can be resolved using the FDLS method. Figure 9 shows the simulated change in elongation and angular velocity for a drop that is subjected to a small, time-invariant torque. This procedure replicates commonly FIG. 9. The elongation and dimensionless angular velocity of an uncharged simulated drop that is spun up by the application of a small, time-invariant torque. The transition from prolate to capsule and dumbbell shapes is illustrated by the drop shapes at points A-D as displayed in Fig. 10. The transition to these capsule and dumbbell shapes corresponds to the loss of stability of prolate spheroids in the spheroid model as indicated by the gray cross. used experimental methods. 25,31,35 The parameters for a water drop of initial diameter 1 cm are used, and the z-axis of the simulation is aligned perpendicular to the rotation axis so that prolate and dumbbell shapes are captured. Figure 9 clearly exhibits the bendback of the curve, but the maximum angular velocity has a slightly lower value, Ω * = 0.48, than expected from the spheroid model and from previous theoretical and experimental studies. These discrepancies are likely to be caused by the axisymmetric and spheroidal constraints on the models, whereas real drops can transition through triaxial shapes. However, the dimensionless angular momentum of the drop at this point is 0.618, which is in excellent agreement with the value of Brown and Scriven 29 of L * = 0.628.
The angular momentum is steadily increased beyond the bendback value, and the resulting drop shapes with angular momenta L = {0.6, 0.8, 1.1, 1.146} are shown in Fig. 10. These correspond to the points labeled A-D, respectively, in Fig. 9. The formation of symmetrical dumbbell shapes is reproduced by the code. The drop eventually breaks up when the angular momentum reaches L = 1.148, which is slightly larger than the numerically calculated value of Brown and Scriven 29 of L max = 1.08 but in excellent agreement with Liao and Hill's experimentally determined value 25 of L max = 1.15.

C. Stability of charged, rotating drops
The behavior of drops with both charge and rotation is now considered. This is easily achieved within the spheroid model by including all of the surface, rotational, and electrostatic energies in the expression for the total energy, which is then minimized. The resulting bifurcation diagram for prolate spheroids is displayed for various values of x in Fig. 11 with Ω providing the bifurcation parameter. The results are unsurprising; increasing the charge on the drop reduces the value of Ω at the bendback point where the prolate spheroids lose stability. The decrease in this value of Ω is modest for x < 0.5, but it rapidly approaches zero as x tends to the critical Rayleigh charge.
The FDLS simulations allow the behavior of the drops to be probed beyond the point where prolate shapes become unstable. The simulations are initiated in the same manner as the water drops in Subsection IV B but with a given amount of charge placed on the drop at the start of the simulation. This procedure replicates the experimental method of Liao and Hill. 25 The results are displayed in Fig. 12, which shows the change in angular velocity and elongation as the drops are The results are plotted as dashed lines in Fig. 13, and excellent quantitative agreement between the two methods is observed. These results are also in excellent quantitative agreement with the experimental results of Liao and Hill, 25 which are displayed in the right panel of Fig. 13. The experimental results exhibit a small increase in the transition value of L above that of the spheroid model for x > 0.8, which is attributable to the confining force of the diamagnetic trap; this force is not included in the spheroid model or FDLS simulations.
FIG. 12. The elongation and angular velocity of simulated drops that are charged and slowly spun up through the application of a small constant torque. The behavior resembles that of the spheroid model in Fig. 11. However, the points above the bendback now correspond to stable capsule and dumbbell shapes rather than unstable prolate spheroids.  Figure 13 also shows the values of L and x that cause the drop to become unstable according to the FDLS simulations. These stability limits are, again, in excellent agreement with the experimental results. A least-squares quadratic fit to the upper curve of the left panel of Fig. 13 is L max = 1.15 − 0.59x − 0.56x 2 , which provides a good working formula for the stability limits of drops for engineering applications. However, caution is advised in the use of this formula for highly charged drops with x approaching 1 that are unstable against perturbations that remove the drop from its local energy minimum, as discussed in Sec. IV A; these highly charged drops are rarely observed in experiments despite being below Rayleigh's linear stability limit.

D. Spinup/spindown hysteresis
A striking conclusion from experimental studies of rotating drops is that the value of Ω at the transition between oblate and prolate shapes, and vice versa, depends on whether the drop spinning up, i.e., it is having its angular momentum FIG. 14. Experimental evidence for hysteresis as reproduced with permission from L. Liao and R. J. A. Hill, "Shapes and fissility of highly charged and rapidly rotating levitated liquid drops," Phys. Rev. Lett. 119, 114501 (2017). Copyright 2017 Author(s), licensed under a Creative Commons Attribution 4.0 License. Note that R * is defined in terms of the longest axis of the drop; it is identical to this paper's definition of α for elongated drops with α ≥ 1 but not for those with α < 1. The variable Ω * corresponds exactly to Ω as defined in this paper.
increased, or spinning down. 25,31 An uncharged drop that is spinning up reaches an angular velocity of Ω ≈ 0.66 ± 0.02 in Fig. 3 of Liao and Hill, reproduced here in Fig. 14, which is significantly higher than the theoretical maximum stable velocity of Ω * = 0.559. This hysteretic behavior is easily explained using energy arguments within the spheroid model. Figures 7  and 8 show that oblate spheroids with angular velocities above Ω * are stable against small shape perturbations but unstable to finite-sized perturbations that allow the spheroid to escape from its local energy minimum. If Fig. 8 is extended to higher velocities, then the oblate spheroid shape family terminates with a saddle-node bifurcation at {Ω = 0.968, α = 0.020}, but real drops are expected to lose stability and assume a prolate shape well before this point is reached. The prediction of the stable oblate family with angular velocities Ω > Ω * is, as previously discussed, in disagreement with the calculations of Brown and Scriven 29 but is strongly supported by the experimental evidence.
The FDLS simulations provide evidence for an additional hysteretic transition between prolate and dumbbell shapes. Figure 15 shows data for the same drops as   Fig. 14 is observed once the difference between R * and α for oblate drops with α < 1 has been accounted for. at some point after they have become dumbbell-shaped. The curves followed subsequently by these drops undershoot the spinning-up bendback value of Ω. Evidence for this undershoot is also seen in experiments. 25,31 The simulations are constrained to axisymmetric shapes, with the axis of symmetry aligned perpendicular to the rotation axis, so they cannot model oblate spheroids rotating about their minor axis; the hysteresis must therefore be in the transition between prolate and dumbbell shapes. The stability of the dumbbell shape family at angular velocities below Ω * could be investigated with an energy minimization procedure applied to dumbbell-like shapes composed of two intersecting spheres and would make for an interesting future study.
A sketch of the full hysteresis loop is displayed in Fig. 16, which incorporates the transitions between oblate and elongated shapes and between dumbbell and prolate shapes. Drops that are spun up from rest assume oblate shapes with increasing deformation, and hence decreasing values of α, as their angular momentum is increased. Oblate spheroids appear to be preferred to prolate spheroids due to the deeper potential well on the oblate side of Fig. 7. The drop then transforms into elongated dumbbell shapes via a triaxial shape family. These triaxial shapes cannot be modeled using the axisymmetric methods of this paper, and the dashed line in Fig. 16 is intended to indicate the expected behavior only. Once the drop has assumed a dumbbell shape, it is spun down, through prolate shapes, back to a static sphere. This sketch fully accounts for the experimental observations in Fig. 14 once the difference between R * and α for oblate drops with α < 1 has been taken into account.

V. SUMMARY
The determination of the shape and stability of charged and rotating drops is an important problem in both fundamental and applied physics; highly charged and rapidly rotating drops are particularly common in plasma applications. This paper has described the development of a spheroid energy minimization method and complementary simulations using the finite-difference, level-set (FDLS) method. The resulting calculations are in excellent agreement with each other and with earlier theoretical, computational, and experimental studies. A good working formula for the dimensionless angular momentum required for the disruption of a charged drop has been proposed accordingly as L max = 1.15 − 0.59x − 0.56x 2 . This stability limit is appropriate for use in numerous engineering applications. Drops just below the Rayleigh limit x c = 1, although stable against very small perturbations, are unstable to finite-sized perturbations that allow the drop to escape its local energy minimum; these drops are rarely observed in experiments.
The transitions between the different stable shape families of the drops have also been considered in detail. The energy minimization procedure shows that the oblate spheroid shape family, contrary to earlier calculations, 29 remains stable for angular velocities above the transition point Ω * . Oblate drops that are spinning up can, therefore, retain a stable oblate shape beyond the transition value; conversely, prolate drops that are spinning down can remain prolate below the transition point. The dependence of the drop shapes on their history is a clear example of hysteresis and is strongly supported by experimental observations. 25,31 The FDLS simulations provide evidence for a further hysteretic transition between prolate spheroids and dumbbell shapes.
In summary, the stability of charged and rotating drops has been considered using a spheroid energy minimization method and FDLS simulations. This work provides the first theoretical study of the hysteresis of these drops that strongly supports the findings of recent experiments. The results are directly applicable to numerous industrial applications and, given the historical successes of liquid drop models, should continue to provide important physical insights into fundamental processes from nuclear to astrophysical scales.

ACKNOWLEDGMENTS
This work was supported by the UK's Engineering and Physical Sciences Research Council, the Imperial College President's Ph.D. Scholarship Scheme, and an Aerosol Society C. N. Davies Award. The code and data underlying this work are available at http://www.github.com/joshholgate/iEHD.