Multiplicity of Inertial Self-Similar Conical Shapes in an Electriﬁed Liquid Metal

Liquid-metal ion sources (LMIS) are widely used for ion implantation in semiconductors, focused ion-beam systems for lithographic patterning and now even small in-space thruster units for precision pointing. Above a critical ﬁeld strength, a liquid metal just prior to ion emission forms a conical protrusion, which undergoes accelerated sharpening of the tip due to ﬁeld self-enhancement. Despite decades of interest in this phenomenon, the inﬂuence of inertial eﬀects on the shape and ﬂow ﬁeld prior to emission remain poorly understood. Zubarev [Formation of conic cusps at the surface of liquid metal in electric ﬁeld, JETP Lett. 73 , 544 (2001)] showed that when local Maxwell and capillary forces prevail at an electriﬁed tip of a perfectly conducting liquid in inviscid ﬂow, the asymptotic behavior of the potential ﬁelds and interface shape far from the tip can be described by a one-parameter family of self-similar series solutions describing a dynamic Taylor cone with radially convergent ﬂow toward the tip. The Maxwell and capillary pressure undergo divergent growth in ﬁnite time, characteristic of blowup phenomena in self-focusing singularity ﬂows. Suvorov and Zubarev [Formation of the Taylor cone on the surface of liquid metal in the presence of an electric ﬁeld, J. Phys. D: Appl. Phys. 37 , 289 (2004)] later found solutions that retained inertial eﬀects and showed that the self-focusing, divergent dynamics leading to power-law growth are preserved. In this work, we focus especially on the inﬂuence of inertial eﬀects and extend the analysis to include time-reversal symmetry. We provide an interweaved procedure for calculating the coeﬃcients of the series expansions for the potential ﬁelds and interface shape and use these to deduce a compact relation for bounds on the minimum and maximum liquid apex height achievable. We then develop a boundary integral patching technique, which yields the complete self-similar solution, valid throughout the near-and far-ﬁeld domain. The resulting two-parameter family of numerical solutions reveals a multiplicity of ﬂuid conﬁgurations, which we coin subconical, superconical, and mixed conical, that exhibit not only tip sharpening but tip bulging, tip separation ﬂow, interface stagnation points, and receding shapes reminiscent of recoil after capillary pinchoﬀ. The existence of such multiple conﬁgurations may ultimately help explain experimental observations during LMIS operation involving tip pulsation, droplet emission, liquid recoil and collapse.


I. INTRODUCTION
The strong surface distortion accompanying electrically stressed liquids has fascinated researchers for centuries dating back to experiments in the early 1600s by Gilbert [1], who reported emission of a fine jet of liquid when water was attracted to a highly charged piece of amber, a charged drop is in a state of unstable equilibrium and will undergo ejection of one or more slender jets. Several years later, Larmor [3] set out to quantify the effect of surface charge on the propagation of capillary waves on a deep layer of fluid. Using the unsteady form of the inviscid surface Bernoulli equation and noting that the presence of a surface charge diminishes the surface pressure, he derived the velocity of propagation of capillary waves along an electrified interface. The result indicated that the presence of surface charge is mathematically equivalent to a reduction in the surface tension by 4σ 2 λ where σ is the electric charge surface density and λ the surface undulation wavelength. These two early studies were instrumental in launching the field of electrohydrodynamics, which couples the nonlinear system of equations describing fluid flow with Maxwell's equations.
Over the centuries, physicists have explored a wealth of related phenomena and over time have established key differences in response between perfect conducting, perfectly insulating, and so-called leaky dielectric liquids, the influence of surface kinetic energy on protrusion shapes, conditions required for liquid-jet initiation, liquid-thread breakup and progeny droplet formation, and more. The advent of high-resolution imaging techniques such as highvacuum transmission electron microscopy (HV TEM) is also helping researchers examine in more depth the cusplike formations in liquid-metal alloys prior and during ion emission. A particularly noteworthy feature of the dynamical process became evident during the last two decades through the work of Zubarev and colleagues [4][5][6] who demonstrated by asymptotic analysis that field selfenhancement near the tip of an electrified protrusion progresses through a runaway process controlled not by the external field but the local field at the liquid tip. This work revealed so-called self-similar blowup in finite time leading to accelerated divergence of the local Maxwell (electric), capillary, and kinetic pressure, which presumably culminates in such high electric field strengths that ion emission is triggered.
The complexity of the underlying coupled set of nonlinear equations are not amenable to analytic solution except in certain asymptotic limits. Zubarev determined the asymptotic form of the electric potential, velocity potential, and interface shape applicable to distances far from the tip apex for an inviscid, perfectly conducting liquid and demonstrated power-law divergence in finite time of the surface pressure, a characteristic of blowup phenomena in self-focusing singularity flows. In what follows, we first review in Sec. II key studies of the 20th century, which established principles necessary to the understanding of such phenomena and helped differentiate systems of interest. We provide there estimates showing why liquid metals are well described as perfectly conducting fluids and discuss the validity of the inviscid approximation in examining the dynamic behavior of an electrified liquid metal. In Sec. III B, we outline the original electrohydrodynamic analysis based on the surface Bernoulli equation by Zubarev, which reveals a one-parameter family of asymptotic self-similar solutions for the potential functions and interface shape. We show why these solutions capture only flow configurations for which the surface kinetic energy per unit volume is negligible in comparison to the Maxwell and capillary pressure. In Sec. IV, the analysis is extended to include explicitly the kinetic energy term and augmented to allow for time reversal symmetry for capturing solutions beyond the so-called blowup point. This leads to a two-parameter family of asymptotic self-similar solutions which expand allowable flow configurations to include subconical, superconical, and mixed-conical interface shapes. We provide an interweaved procedure for calculating the coefficients of the series expansions for the velocity potential, electric potential and interface function and use these expressions to deduce a compact relation for bounds on the minimum and maximum liquid height possible. In Sec. V, we provide the complete self-similar solutions, valid throughout the near and far field from the liquid tip, obtained by a boundary integral patching technique, which directly incorporates the asymptotic solutions into the formulation.

A. Important early studies
In the century that followed the work of Rayleigh and Larmor, significant advances in understanding became possible through the beautiful and insightful experiments by Zeleny from 1914Zeleny from -1920. His work produced not only stunning images but key information about the critical field strengths needed to initiate ejection of a fine charged spray, which he obtained by an elegant measurement technique based on hydrostatic force balance [7][8][9][10]. He developed a number of visualization methods for capturing the distortion and emission process in small hemispherical droplets protruding from the end of a fine capillary when exposed to field strengths of the order of a MV/m. Using liquids such as water, alcohol, dilute solutions of hydrochloric acid or carbonic acid, and capillaries of different radii made of glass or metal, he was able to capture a variety of phenomena. His images revealed that above a critical field strength, the apex of a hemispherical droplet deformed into a cusplike shape and spontaneous emission of a fine liquid thread, which fragmented into a spray of minute droplets. He also observed tip pulsation, thread whipping motion, liquid recoil and collapse. He commented that behavior in such liquids-now called leaky dielectric liquids (i.e., liquids with finite conductivity)-stood in contrast to his studies with liquid mercury showing rather stable current emission presumably from an ion current. Many experimentalists using aqueous droplets and even aqueous soap solutions [11][12][13][14][15] confirmed similar wide-ranging behavior.
Zeleny also derived a stability criterion for the critical surface charge or critical field strength of an isolated ellipsoidal droplet with constant surface charge. This criterion was based on the fact that the first derivative of the sum of the normal pressures acting on the interface with respect to a coordinate specifying the liquid shape must vanish at an unstable point on the surface. His result reduced exactly to the Rayleigh criterion for a spherical liquid droplet. In a comprehensive paper written in 1935, he concluded that the magnitude of the field strength for emission was in line with his prediction. However, measurement techniques at the time lacked the resolution required to determine whether the field strength at the point of discharge is the same as that required to initiate protrusion growth, an issue still under investigation today. While most of these early studies involved aqueous ionic liquids, researchers like Beams and Quarles [16,17] began focusing on liquid metals such as liquid mercury to better understand differences in behavior. They attempted to measure more accurately the critical field strength required for ion emission by varying the externally applied voltage while trying to present large-scale distortion or movement of the liquid. To do so, they applied impulsive fields of the order of 10 7 V/m of short duration times 10 −7 -10 −6 s. These studies helped catalog estimates of the work function required for electron emission but could not shed light on the distorted shapes or field values generated just prior to ion emission. (When liquid metals are positively biased, they emit ions and when negatively biased generate electrons.) By 1935, experimentalists had documented a wealth of phenomena accompanying the distortion of electrified stressed liquids subject to high field strengths. Tonks [18], a plasma physicist working at General Electric Company, was most intrigued by the fact that surface distortion and rupture leading to charge emission in liquid metals required far smaller values of the applied field strength than equally smooth solid metal surfaces. After careful reexamination of key studies by Zeleny and others [7][8][9][10][11][12][13][14][15], Tonks quickly set about developing a dynamic model based on the accelerated pressure imbalance of a small hemispherical protrusion in an otherwise planar, perfectly conducting liquid subject to a critically large uniform electric field. Tonk's analysis yielded a key relation between the amplitude of the initial liquid distortion to the so-called rupture time and field strength, which for liquid mercury he estimated to be 53 kV/cm. He predicted that the linear dimensions of an evolving protrusion will therefore vary inversely as the square of the field strength and the time to rupture will vary inversely as the cube of the field strength. From this, he estimated that a hemispherical bump of liquid mercury with an initial maximum height of 0.4 μm and a radius of 9 μm subject to an initial uniform field of 10 8 V/m will sharpen in time and undergo surface rupture in just 5.4 μs.
Realizing that a rigorous analytic solution to such a complex electrohydrodynamic (EHD) problem was too formidable a task, Tonks instead relied on an insightful approximate treatment of an electrified perfectly conducting liquid that revealed critical aspects of the distortion runaway process at late time. Tonks focused on the dynamics of a single liquid protrusion and developed an equation of motion based on conservation of momentum relating the local acceleration rate at the liquid apex to the degree of eccentricity of the sharpening interface. The analysis revealed that beyond a critical field strength, the fluid pressure at the conical apex cannot maintain an equilibrium state since the Maxwell pressure increases as the square of the apex height while stabilizing forces due to capillary or gravitational pressure at most scale inversely or linearly with apex height, respectively. His model established then that any small advance of an electrified liquid tip will continue to grow in length and narrow in width without bound until the point of emission is reached, due to accelerated growth in the unbalanced pressure at the conical tip.
He realized, correctly so, that the accelerating imbalance causes a runaway process from field self-enhancement driven by the shrinking radius of curvature of the electrified tip. This insight led him to the realization that the growth of the liquid tip should proceed in self-similar fashion, which he deftly portrayed through hand-drawn sketches in Fig. 8 of Ref. [18] showing self-similar growth of a small liquid bump in time. Shortly after publication of this work, Frenkel [19] offered a different perspective on the problem based on a hydrodynamic and not mechanical approach. He assumed the existence of an initially long and flat liquid layer able to sustain multiple periodic protrusions and from the inviscid surface Bernoulli equation derived a dispersion relation for the speed of electrified capillary waves. The wave speed was found to vanish at the exact same critical field strength in Tonks' analysis required for accelerated growth of a bump. The fact that the condition for the instability of an electrified liquid surface reduces to the vanishing of the wave speed seemed at first surprising. Frenkel explained that this condition corresponds to the field strength at which the wave velocity and wavelength both become imaginary, which in the analysis then reduces to a wave amplitude that increases in time. The well-known instability describing periodic distortion of the surface of a perfectly conducting liquid above a critical field strength is now known as the Larmor-Tonks-Frenkel instability.

B. Hydrostatic Taylor cone and subsequent studies by Miscovsky and co-workers
In 1964, Taylor [20] decided to revisit Zeleny's criterion [8] for the condition required for the stability of an isolated charged ellipsoidal droplet. Taylor somehow misinterpreted Zeleny's analysis and mistakenly believed that 044001-3 the derivation required the gauge pressure to vanish at instability. He then argued that a gauge pressure dependent on position was required for instability, which in fact Zeleny had emphasized in his work. Taylor was aware from experiments of his own and others that spray emission tended to ensue immediately after protrusion growth without having to increase the applied field strength. He also knew from images obtained just prior or during liquid emission that an electrified tip tends to adopt a conical shape from whose apex a jet was discharged. He was intrigued by this particular geometric feature and wondered what geometric shape would an infinite mass of liquid in hydrostatic equilibrium have to adopt in order for the negative Maxwell pressure to cancel the positive capillary pressure everywhere on its surface. He found that such a perfect balance was indeed possible for the shape given by a static liquid cone with interior half-angle θ T ∼ = 49.2923 • , now known as the Taylor angle. This solution required however a curved counter electrode forming a coordinate surface with radius r = r o [P 1/2 (cos θ)] −2 , where r o denotes the central distance to the liquid surface and P 1/2 is the Legendre function of the first kind with index 1/2. Taylor did not remark on the fact that the Maxwell and capillary pressure diverge to infinity at the conical singularity point, although he surely realized that such divergence rendered a state of unstable equilibrium. Unfortunately, this simple exercise in electrostatics led to considerable subsequent confusion in the field and the incorrect notion that an electrified conical tip represents a state of static equilibrium. Researchers unfamiliar with Tonks' work tended therefore to neglect the influence of inertial forces during the late stages of tip sharpening. This misconception, which sometimes appears in the modern literature as well, is based on neglect of the dynamic response of the liquid to the electric stress.
By the early 1980s, researchers began using variational analysis to determine whether certain liquid and counter electrode shapes might allow stable equilibrium configurations. In 1983, Miscovsky and co-workers [21] demonstrated that a static mass of perfectly conducting fluid in the shape of a Taylor cone, or a cone with any angle for that matter, is inconsistent with a variational formulation of the equilibrium problem. They traced the problem to two issues with Taylor's analysis: neglect of the fluid pressure and truncation of the Legendre expansion for the electrostatic potential to the first term (P 1/2 ). Using a similar variational approach, Chung et al. [22] later demonstrated that even an ideal cuspidal shape is incompatible with the variational formulation. In another study [23], they reported yet another problem with Taylor's analysis, namely that the field-strength criterion for static equilibrium he derived, if also considered the tipping point for instability, would imply the existence of a global and not a local condition. Therefore, such an instability should occur simultaneously across the entire surface of the cone, an assumption at odds with experimental observations showing emission from a single point. To resolve these and other issues with the Taylor analysis, Chung et al. [23] formulated the problem in terms of the time-dependent inviscid surface Bernoulli equation coupled to the requisite potential field equations and boundary conditions, which gives rise to a set of nonlinear equations not generally amenable to analytic solution. They did not provide any details of the analysis or solution until a couple of years later [24].
During this period, they devoted their efforts to formulating stability criteria based both on an electrohydrostatic (EHS) approach and an EHD analysis. The EHS analysis [25], which assumed static equilibrium over the entire surface of a fluid whose interface is described by coordinate surfaces of a cone, cusp, or a pointed hyperboloid of revolution. Using the Zeleny criterion for instability, based on a vanishing first variation in the sum of the normal stresses, they showed that the variation does not vanish but in fact diverges at the apices of such shapes. The conclusion from that study was that a static Taylor cone is an unstable configuration whose apex is likely to disintegrate spontaneously. They used the word "disintegrate" in reference to Taylor's original description of droplet distortion accompanied by the formation of a pointed end, which discharges small droplets or ions [20].
Their EHD analysis relied on a time-dependent, unbalanced rate of change in the net fluid pressure to drive fluid flow. Their first attempt at a rigorous calculation [23] valid to first order was restricted to long wavelength disturbances acting on a static Taylor cone. The results showed that the tip of a Taylor cone should rapidly deform into a concave shape with a rounded apex and progressive tip sharpening under perturbation from an external electric field. In subsequent work [24,26], they appealed to the full Navier-Stokes equation and carried out a linear stability analysis to derive an expression for the amplification of surface capillary waves in a viscous fluid valid for all wavelength disturbances.
In the years that followed, Driesel and co-workers [27][28][29] were able to successfully modify a 1-MeV HV TEM to allow in situ observation of electrified tip shapes as a function of increasing ion-emission current in various liquid metals and alloys. These high-resolution images revealed details never before seen within just a few microns of the apex just prior and subsequent to ion emission. Aware of Tonks' earlier hypothesis of runaway growth at an electrified tip and armed with intuition gained from these images, a researcher in the Soviet Union by the name of Zubarev [4] made a key contribution to the field. His work established that tip field self-enhancement is in fact driven by a self-similar runaway process leading to blowup in finite time. We shall outline the details of his analysis in Sec. III but first address two assumptions of the model.

C. Tip shapes: Perfect conductors or dielectrics versus leaky dielectrics
Recent numerical simulations have revealed in detail why the tip shapes of electrified liquids modeled as perfectly conducting or insulating (i.e., perfect dielectrics) differ significantly from those in leaky dielectric liquids characterized by a finite conductivity [30]. The difference in flow behavior and tip shape is directly related to the difference in time scale for charge relaxation [31], τ charge = o liq /σ liq and fluid flow, τ flow = ( Here, liq and σ liq denote the liquid dielectric constant and conductivity, o ≈ 8.854 × 10 −12 C 2 /N m 2 is the vacuum permittivity constant, γ and ρ denote the liquid surface tension and density, and E o denotes some characteristic electric field strength. For a typical liquid metal like indium [32] used in LMIS systems, γ = 0.556 N/m, ρ = 7.03 × 10 3 kg/m 3 , σ liq = 3.2 × 10 6 S/m (Ref. [33]), and liq ≈ 40 (Ref. [34]). For typical initial operating voltages of about 4 × 10 7 V/m, the ratio τ charge /τ flow ∼ O (10 −20 ). A liquid metal is therefore well modeled as a perfectly conducting liquid in which charge relaxation is instantaneous in comparison to any flow time scale. By contrast, the charge relaxation time for a semi-insulating liquid like corn oil is about one second [31].
Recent simulations and experiments [30,35,36] have shed further light on this distinction and helped elucidate which types of liquids allow and which forbid formation of progeny droplets from an electrified tip. Liquids with finite conductivity allow surface charge transport and redistribution on time scales comparable to the flow, which in turn generates surface tangential stresses leading to EHD tip streaming [30]. Tip streaming describes a process whereby a fluid tip elongates into a very slender filament prone to capillary breakup and droplet detachment. By contrast, perfectly conducting or insulating liquids, which can only sustain electric fields oriented in the direction normal to the interface, can only develop a sharpened tip from whose apex an ion beam is discharged. In the absence of any surface oxidation effects or surface contamination, molten metals cannot support tip streaming and therefore no progeny droplet formation resulting from breakup of a liquid thread. Shown in Fig. 1(a) is a high-resolution image of the tip of a AuGe liquid alloy droplet undergoing ion emimssion obtained in situ by HV TEM. Shown in Fig. 1(b) is a recent scanning electron micrograph of the surface of molten indium after exposure to an extremely large field gradient of the order of 100 MV/m and subsequent solidification. The latter formations are currently of great concern due to their suspected role in causing breakdown of accelerators such as the 30-GHz Compact Linear Collider at CERN. Researchers are actively seeking alternative fabrication and conditioning processes better suited to extreme environments involving very high accelerating gradients such as in rf structures.

D. Validity of the inviscid approximation for liquid metals
All the analysis in this current work is based on the inviscid form of the Bernoulli equation for quantifying the dynamic behavior of liquid metals prior to ion emission. One may ask whether the inviscid approximation is well justified in this case. The following back of the envelope estimate offers good support for this approach. As is well known, an initially quiescent liquid can only generate vorticity through boundary motion [38]. For the problem at hand, vorticity generated both by interface curvature and acceleration diffuses into the bulk liquid through viscous stresses initially confined within a viscous boundary layer whose thickness is of the order of μ 2 /ργ , where μ is the 044001-5 shear viscosity. For liquid metals [32], the viscous boundary layer is estimated to be extremely small, of the order of ten nanometers. Also, conical formations are observed to form almost immediately after exposure to a large electric field, with ion emission occurring within a few nanoseconds to microseconds after startup. This allows little time for surface vorticity to penetrate any significant distance within the bulk. Given these estimates then, it is anticipated that the EHD process is well modeled by inviscid flow.

III. ZUBAREV PREDICTION OF SELF-SIMILAR CONICAL GROWTH
In this section we outline and then assess the implications of the original analysis by Zubarev [4], a leading physicist in the Nonlinear Dynamics Group at the Institute of Electrophysics of the Russian Academy of Sciences. Unlike the many previous studies that viewed the tipgrowth process as an instability, he instead focused on formulation of the late time dynamical behavior of an electrified tip in a liquid metal just prior to ion emission. His main interest was whether the liquid apex undergoes selfsimilar growth leading to a runaway process and eventual blowup in finite time.
Any self-similar process requires that local and not global conditions prevail in the region of interest and that the process allow self-replicating behavior at smaller and smaller scales until there is breakdown due to critical conditions and intervention of another physical mechanism. Zubarev therefore replaced the usual external field uniformity condition in the electrodynamic formulation of the problem with a decay condition specifying vanishing field strength far from the liquid tip. He reasoned correctly that the local electric field strength in the apical region of a liquid tip is expectedly to exceed the external value, rapidly and appreciably, if local conditions prevail. This insight allowed him to analyze fluid motion near the apex without reference to any particular electrode geometry. Scaling of the governing EHD equations under dilation revealed a set of self-similar transformations amenable to asymptotic analysis. He demonstrated that in the laboratory frame, the asymptotic solutions exhibit finite-time blowup wherein the Maxwell and capillary pressure undergo divergence to infinity, a process directly related to the diminishing radius of curvature of the liquid tip. We review Zubarev's original approach in detail below [39].

A. Notational symbols
To keep track of the notation in this work, we note that vector and tensor quantities are indicated by bold-face variables and partial differentiation by subscripts. Dimensional variables are designated by lowercase Roman or Greek letters decorated by a tilde, e.g., surface velocity potentialψ( r,t). Dimensionless variables are designated by uppercase Roman or Greek letters, e.g., (R, T) and dimensionless self-similar variables by lower-case Roman or Greek letters, e.g., ψ(r, z). The exception to these rules is the electric field distribution denoted by E in dimensional form and E in dimensionless form. Unit normal vectors on boundaries point outward from the volume of interest. In this study, it is assumed that the system manifests axisymmetry about the central vertical axis. As such, it proves convenient to introduce both cylindrical coordinates (r, z) and spherical coordinates (r, θ), depending on the form of the expressions desired. Projected lengths along the central axis of symmetry are therefore given by z = r cos θ , which can be a negative value. In switching between coordinate systems as needed, we caution that the symbol r denotes the radial coordinate in spherical coordinates and r denotes the radial coordinate in cylindrical coordinates.

B. Asymptotic one-parameter family of inertialess self-similar solutions
Here we review and expand on Zubarev's analysis for an electrically stressed axisymmetric protrusion emanating from a perfectly conducting fluid in vacuo subject to incompressible, inviscid and irrotational flow. In a perfectly conducting liquid with no net charge, all mobile charges reside on the liquid interface and rearrange instantaneously (in comparison to the time scale for fluid motion) to maintain the liquid mass at constant electric potential ψ. The interior liquid domain therefore defines a Gaussian volume devoid of an electric field. Consequently, the electric field at the surface of the liquid can only sustain a normal component, Eñ, which we everywhere simply designate by E. Assuming all boundaries except that of the moving liquid are held stationary, the electric potential distribution in the vacuum domain varies in time only in response to the instantaneous location and shape of the moving liquid surface. (For an isolated charged liquid mass, the electric potential distribution depends only on the shape of the liquid.) The difference in the corresponding electric field distribution across the liquid-vacuum interface then gives rise to a jump in the Maxwell stress tensor E E T − | E| 2 I/2, which causes a jump in the electrostatic pressure [40] given by − o E 2 /2. The negative sign reflects the fact that the net interfacial electrical stress acts to pull liquid toward the vacuum region. Setting the gauge pressure in the vacuum to be zero, the total pressure acting on the fluid interfacep int is where the mean curvature in dimensional units is defined by H = −( ∇ ·ñ)/2. We note that this expression for the interfacial pressure can also be obtained from the first variation in shape for a liquid volume whose energy is exclusively governed by electrostatic and capillary forces [41]. The electric field distribution E in the vacuum domain and along the moving interfacez =h(r,z,t) is obtained from the gradient of the electric potential ψ(r,z,t) and must satisfy ∇ · E(r,z,t) = 0 ( 3 ) Likewise, for an ideal liquid subject to incompressible ( ∇ · u = 0) and irrotational ( ∇ × u = 0) flow, the velocity field u in the liquid domain and along the moving interface is obtained from the gradient of the velocity potential ψ(r,z,t) and must satisfy The unsteady form of the inviscid Bernoulli equation evaluated at the free surfacez =h(r,t) is prescribed by It is assumed that gravitational forces are negligible in comparison to all other forces in the problem. The underbrace terms describe the various contributions to the surface pressure driving the rate change in the velocity potential. We use the term "kinetic pressure" to signify the kinetic energy per unit volume acting on the moving interface, which derives from the term in the Euler equation representing advective acceleration of the fluid in the laboratory frame.
Zubarev introduced the following scalings for nondimensionalization based on a characteristic value of the field strength E o , according to which the dimensionless equations are and the dimensionless boundary conditions are The last condition is the usual kinematic boundary condition for free surface flows, which requires that the material points on the moving interface move with the local flow speed of the interface. We note that despite the fact that the original system describes a liquid subject to applied external electric field, Zubarev nonetheless invoked boundary conditions specifying vanishing velocity and electric potential in the far field, in line with our earlier discussion at the start of Sec. III describing essential features of any self-similar process. The underlying assumption for this system is that the electric field in the region of high interface curvature at the liquid tip rapidly amplifies from accelerated field selfenhancement, thereby quickly exceeding the magnitude of an external field. By comparison then, the external field 044001-7 far from the tip becomes negligible. This same reasoning is applied to the velocity potential since the flow speed caused by field enhancement also quickly exceeds that of the far field, which by comparison becomes negligible.
Zubarev noted that since the governing equations and boundary conditions are invariant to dilation under the transformations ( , ) → (α , α ), (R, Z, H ) → (α 2 R, α 2 Z, α 2 H ) and T → α 3 T, the system of equations support similarity solutions of the form where The dimensionless blowup time T C defines the asymptotic time at which the Maxwell and capillary pressure diverge to infinity at the liquid tip due to field self-enhancement. Identical scalings as these are known to occur in problems involving capillary-inertial pinchoff in inviscid flow [42][43][44] and have been verified in experimental studies of liquid curvature collapse in a viscous liquid [45]. Finitetime singularities such as these typically represent local divergences in the amplitude or gradient of a physical observable in finite time and are found in many physical systems, not just hydrodynamic ones. The scaling in Eq. (14) also implies that in the laboratory frame, the interface curvature of the liquid tip will shrink rapidly as lim τ →0 τ 2/3 → 0, a process leading to accelerated tip sharpening. The governing equations in the self-similar frame become subject to the boundary conditions where the origin (r = 0, z = 0) of the coordinate systems is made to coincide with the blowup point. Zubarev succeeded in finding a self-consistent set of asymptotic self-similar solutions to this system of equations for the velocity potential ψ, electric field potential φ, and interface height h represented here in cylindrical coordinates [46], according to which where P 1/2 (cos θ) is the Legendre polynomial of the first kind of order 1/2 and θ = tan −1 (r/z) denotes the polar angle in spherical coordinates. We note that as r → ∞, Eq. (21g) reduces to the condition that h r = 0, which defines a surface of constant slope. In addition, for the leading-order term in the electric potential to satisfy a null equipotential condition on the liquid surface, P 1/2 (cos θ 0 ) must vanish, which requires cos θ 0 ≈ −0.6522 or θ 0 ≈ 130.7077 • . Equivalently, this defines a liquid cone with an interior half-angle θ T = 49.2923 • , the classic Taylor angle. Zubarev highlighted the fact that in the laboratory frame this shape describes a dynamic cone, which he remarked was a dynamic analog to the classic static Taylor cone [20]. The first few coefficients of the asymptotic series are found to be Zubarev therefore uncovered a one-parameter family of asymptotic self-similar solutions whose projected length along the central axis, namely h(r) = r cot θ 0 + c 2 /r 5 where c 2 > 0, describes an interface lying above the conical envelope with interior half-angle θ T . As r = (r 2 + z 2 ) 1/2 → ∞, the leading behavior in spherical coordinates is found to be The asymptotic velocity potential ψ therefore depends solely on r and the streamlines converge radially toward the liquid apex, as illustrated in Fig. 3(b). As a check on the flow direction, we note that along the central axis, ψ(r = 0, z) = a 0 /z where a 0 > 0 and z < 0 since the z-coordinate values lie below the blowup point at z = 0. Therefore the fluid velocity (∂ψ/∂z) r=0 = +a 0 /z 2 > 0. This confirms radial flow directed toward the apex of the conical envelope in a radially symmetric manner, as depicted in Fig. 3(b).

C. Physical implication of Zubarev's original solution
A key feature of this one-parameter family of solutions becomes evident in examining the magnitude of terms in the surface Bernoulli equation, here expressed in spherical coordinates: where κ denotes the mean curvature. Substitution of Eqs. (32)-(34) into Eq. (35) and noting that 2κ = cot θ 0 /r, reveals that the kinetic pressure on the interface (i.e., second term on the left side) scales as r −4 , while all the remaining terms scale as r −1 . The surface kinetic energy is therefore a negligible contribution to the flow and to leading order, the surface velocity potential is dominated by the Maxwell and capillary force. In a subsequent publication, Zubarev [5] extended his original analysis to the case of a perfect dielectric liquid and determined that the interior angle of the asymptotic conical interface, which then depends on the dielectric constant, is always smaller than θ T .
At the conclusion of that work, Zubarev noted that the full set of coupled equations for a perfectly conducting or perfect dielectric liquid admit a more general solution but there was no more said. Subsequently, he proposed a more general solution [6] but provided only an outline. The real focus of that work was on numerical solution of the full Navier-Stokes equation for Reynolds number of the order of 100 to see if the blowup behavior he had predicted was supported by numerical simulations. Indeed, finitetime blowup behavior was evident even with inclusion of viscous effects although the power-law exponents seemed to veer somewhat from the scalings in Eqs. (14)- (16).

IV. GENERAL FORMULATION LEADING TO ASYMPTOTIC INERTIAL SELF-SIMILAR SOLUTIONS
In the next section, we detail the approach leading to a more general solution for the case of a perfectly conducting liquid and extend the analysis first proposed [6] to include time-reversal symmetry. The fact that the kinetic pressure to leading order is of the same magnitude as the Maxwell and capillary pressure, when coupled with time-reversal symmetry, introduces a two-parameter family of solutions, which exhibit accelerating, decelerating, and mixed-mode flow configurations. As shown below, the general solution imbues the velocity potential with a leading term that scales as r 1/2 P 1/2 (− cos θ), which introduces both radial and angular dependence to the resulting interface shapes. This stands in contrast to the restricted solution discussed above, which scales simply as r −1 .
Shown in Fig. 2 are illustrations of the geometry of interest and notational symbols used for description of the laboratory and self-similar frames of reference (in nondimensional form). We introduce here some notation, which proves convenient for discussion of the boundary integral technique introduced in Sec. V.
In the laboratory frame, we define a Cartesian level set function F(X, T) = Z − (X , Y, T) = 0, where comprises the set of points on the interface separating the vacuum from liquid domain: The outwardly pointing normal vector is then defined as N = ∇F/|∇F|. The usual kinematic condition applied to a moving boundary, which enforces mass and momentum conservation, requires that all material points on the moving boundary be convected with the local surface velocity U . The level set function F(X, T) must therefore satisfy the relation where D/DT = ∂/∂T + U · ∇ and U denotes the surface flow velocity. This constraint assumes that there is 044001-9 Laboratory frame no phase transformation on the moving interface. Accordingly, the kinematic boundary condition simplifies to The general unsteady form of the inviscid Bernoulli equation for the velocity potential , which includes the term representing the fluid kinetic energy per unit volume |∇ | 2 /2, valid everywhere within the fluid volume liq , is given by [38] ∂ ∂T where P denotes the fluid pressure, which at the moving interface assumes the value Since the flow is assumed to be inviscid, the velocity potential therefore satisfies time-reversal symmetry. Accordingly, we introduce a signed time variable τ to allow description of pre-and postsingularity flows: where T C denotes the dimensionless blowup time signaling the divergence of the pressures acting at the singular apex of the asymptotic conical envelope. For τ = T C − T > 0, the fluid advances toward the singularity from below, i.e., a presingularity event. For τ = T − T C > 0, the fluid retracts toward the singularity from above, i.e., a postsingularity event. We introduce a self-similar coordinate vector where X C denotes the intersection of the central vertical axis with the tangent to the conical envelope defined by the envelope defined by the asymptotic Taylor cone, and a dilated time variable given by The corresponding velocity potential, electric potential, and level set function undergo transformation to where the surface boundary γ separating the semi-infinite liquid and vacuum domains, ω liq and ω vac , respectively, is defined by where f (χ , t) = z − h(r, t). As τ → 0, these field variables diverge to infinity. Any residual time dependence t in the transformed fields ψ, φ, and f indicates time dependency over and above the explicit scalings in τ . Because of the isotropic scaling of the spatial coordinates in Eq. (43), the Laplace equation for the velocity and electric potential remain unchanged: (49) where the operator ∇ is understood to act on the selfsimilar coordinate χ. Since the fluid is assumed to be a perfectly conducting liquid, φ = 0 on γ . The transformed surface Bernoulli equation and kinematic condition, respectively, become where κ is the local mean curvature of γ and n = ∇f /|∇f | is the unit normal vector pointing away from the liquid mass. The pressure within the liquid mass is given by In what follows, we restrict attention to field solutions, which are stationary in the self-similar frame of reference such thatψ =φ =ḟ = 0 where the dot denotes the derivative ∂/∂t. The surface Bernoulli equation, recast in spherical coordinates, is where the electric and velocity potential throughout the liquid and vacuum domain must satisfy the relations These boundary conditions for solving this system of equations are chosen to be Decay lim Decay lim with identical far-field conditions as before. We note that as r → ∞, Eq. (59g) reduces simply to the condition r = 0 defining a surface of constant slope. In analogy to the derivation in Sec. III B, this condition when augmented by the equipotential condition restricts the fluid envelope to that of a dynamic Taylor cone with exterior polar angle θ 0 = π − θ T . There are three regimes of interest embedded in the full set of self-similar equations and boundary conditions specified above. The first limiting case, described next in Sec. IV A, represents the classic static solution derived by Taylor [20] illustrated in Fig. 3(a).That unique solution defines the surface shape for which the capillary and Maxwell pressure everywhere cancel, although the conical apex is a singular point. The second regime, which is derived in Sec. III, describes the one-parameter family of solutions obtained by Zubarev [4] and illustrated in Fig. 3(b) in which the surface kinetic pressure is negligible in comparison to the Maxwell and capillary pressure. The third limit, derived below in Sec. IV B, requires that the kinetic pressure in the self-similar form of the surface Bernoulli equation also contribute to leading order.
Here, we explicitly allow for configurations which obey time-reversal symmetry as well. The resulting asymptotic solutions considerably expands the number of fluid configurations possible for an electrified protrusion in a perfectly conducting liquid.

A. Limit of classic Taylor cone solution
In the laboratory frame, the classic Taylor cone solution describes the shape of a stationary mass of electrified fluid acted upon solely by capillary and Maxwell forces. Throughout the volume ω liq , the velocity potential is maintained at ψ = 0. The general harmonic solution [47] for the electric potential φ is then where P ν (cos θ) is the Legendre function of the first kind of order ν, P ν (cos θ) = P −ν−1 (cos θ) and ν denotes a real nonintegral number. The function P ν (cos θ) has a logarithmic singularity at θ = π , which is excluded since φ is confined to the vacuum domain 0 ≤ θ ≤ (r). For φ to be finite at the origin, ν must be positive. The requirement that the liquid represent an equipotential mass requires P ν (cos θ) = 0. Substitution of ψ = 0 (as well as its derivatives) into Eq. (53) reduces the surface Bernoulli equation to the two competing terms on the right-hand side, which reflect the static balance between the capillary and Maxwell pressure. Since the mean curvature κ scales as 1/r and the Maxwell pressure scales as |∇ψ| 2 ∼ r 2(ν−1) , the only allowable solution is ν = 1/2. For P 1 /2, the only zero in the range 0 < θ < π is given by cos θ 0 = −0.6522 and so θ 0 ∼ = 130.7077 and π − θ 0 ∼ = 49.2923 • = θ T . To leading order, this static solution is given by This solution plotted in Fig. 3(a) contains no free parameters. Of note, this solution is not uniformly valid since the electric field strength ∂φ/∂r diverges to infinity at r = 0. The static Taylor cone is therefore at best a metastable state.

B. Asymptotic two-parameter family of inertial self-similar solutions
Solutions which are stationary in the self-similar frame nonetheless allow liquid configurations in the laboratory frame, which are time dependent and can undergo runaway growth. This section focuses then on asymptotic self-similar solutions in which the surface kinetic pressure, capillary pressure, and Maxwell pressure in Eq. (53) all contribute to leading order in r.
Inspection of Eq. (53) reveals this is possible when φ ∼ r 1/2 F(θ ) and ψ ∼ r 1/2 G(θ ). Substitution of these relations in Eq. (53) leads to cancelation of the first term on the left side since 2rψ r = ψ, leaving only terms proportional to r such that 1 r This same scaling ψ ∼ r 1/2 G(θ ) also occurs in studies of capillary-inertial pinchoff [42]. The leading-order solutions are then given by where terms proportional to Q 1/2 (Legendre function of the second kind of order 1/2) are eliminated to prevent divergence of φ along the axis θ = 0 and divergence of ψ along the axis θ = π . The term P 1/2 (− cos θ), which has a logarithmic singularity at θ = 0, is nondivergent since the velocity potential is confined to the liquid domain 0 < (r) ≤ θ ≤ π . Likewise, P 1/2 (cos θ), which has a logarithmic singularity at θ = π , is nondivergent since the electric potential is confined to the vacuum domain 0 ≤ θ ≤ (r). Applying the equipotential potential condition in Eq. (59a) and the kinematic condition in Eq. (59g) yields the recursion relation In contrast to the inertia-free solution in Sec. III, the leading-order term in the velocity potential no longer scales purely as r −1 but instead as r 1/2 P 1/2 (− cos θ), which imbues the field with an angular dependence, as shown in Fig. 3(c). The inertial contribution to the flow causes potential field lines that are no longer purely radial and no longer aligned with the pressure gradient. The complete asymptotic series valid in the limit r 1, here denoted by the subscript ∞, is given by where P −ν−1 (x) = P +ν (x). The series coefficients are given below. In Eq. (71), we reintroduce the cylindrical radial coordinate r = r sin θ so as to absorb the sine function into the coefficients c k , which simplifies later expressions.
The procedure for determining the coefficients is a nontrivial exercise since a k , b k , and c k are coupled together by three nonlinear equations, which must be evaluated along the actual interface h = c 0 r + ∞ k=1 c k h k and not simply the Taylor cone surface h = r cot θ 0 . The interweaved procedure for obtaining these coefficients in illustrated schematically in Fig. 4.These coefficients are computed term by term using the symbolic manipulation software Mathematica [48]. The notation P ν (·) denotes differentiation of the Legendre function with respect to the argument cos θ . The first few terms in the velocity potential are a 0 = Free parameter, a 2 (c 0 , . . . , c 3 , a 1 The first few coefficients for the electric potential are − 8c 3 csc 4 θ 0 + c 3 1 csc 2 θ 0 + 1 + 12c 1 c 2 cot θ 0 csc 2 θ 0 8 csc 5/2 θ 0 P 3 (cos θ 0 ) The first few coefficients for the liquid interface height are Evaluation of the term c 1 in Eq. (74) gives c 1 = 0.368 43a 0 -the coefficients c 1 and a 0 are therefore simply related by a positive constant. When the fluid interface lies below the Taylor cone envelope specified by the k = 0 term in Eq. (71), then c 1 < 0 and therefore a 0 < 0. The general recursion relations for a k , b k , and c k for k > 1 are uniquely determined by the leading-order coefficients a 0 , b 0 , and c 0 . Moreover, since the coefficient c 0 = cot θ 0 , the solutions are specified by only two independent parameters, namely (a 0 , b 0 ), or equivalently, (b 0 , c 1 ), as for example the coefficient a 2 , where The change of variable − cos θ 0 → x in Eq. (74) also establishes that the coefficient c 2 = 0 since

044001-14
The last equality derives from the fact that the differential equation is none other than the Legendre equation for P 1/2 (x). The fact that c 2 = 0 places special relevance on the leading-order correction to the asymptotic conical envelope given by (h ∞ − r cot θ 0 ) = c 1 / √ r.
Once the coefficients for a k , b k , and c k are computed, the series solution for the fluid pressure can be obtained from the stationary form of Eq. (52) where The pressure field, which is subject to time-reversal symmetry, now includes inertial effects, which introduce the dependence on θ.

C. Estimation of liquid apex height from geometric considerations
Maximum and minimum bounds on the liquid height along the central axis of symmetry can be obtained from geometric considerations. For convenience, we resort in this part of the analysis to cylindrical coordinates. To begin, the kinematic condition in Eq. (21g) (valid even in the inertial case), when evaluated about the central axis of symmetry reduces simply to The fluid velocity at r = 0 is therefore controlled by the sign of h(r = 0), which denotes the projection of the interface function onto the axis of symmetry. For example, when T < T C and the fluid tip is advancing from below toward the apex of the conical envelope defined by the Taylor angle, indicated by the red dot in Fig. 6(a), then h(r = 0) < 0 and ψ z | apex > 0. Likewise when T > T C and the fluid tip is advancing toward the conical apex from above, then h(r = 0) > 0 and ψ z | apex < 0. Separated flow is also possible in which the fluid in the liquid tip is moving upward but the nearby bulk fluid is moving downward or vice versa, similar to capillary-inertial recoil in liquid filaments [44,49]. In his original study, Zubarev [4] developed an elegant geometric argument for obtaining a bound on the vertical distance between the liquid tip and the apex of the conical envelope. He noted that since the velocity potential is a harmonic function, the net flux induced by ∇ψ must vanish upon integration over a closed bounded surface. What simplified the analysis is that the correction to the leadingorder term in the expansion for h − r cot θ 0 ∼ O(r −5 ), a term which could be ignored. For the case in which inertial forces are significant, h(r) − c 0 r = c 1 r −1/2 + · · · , which introduces some difficulties; however, we show next that useful and compact relation can still be obtained by invoking some minor additional assumptions.
We refer to the illustrations and variable definitions in Fig. 6 to develop estimates and bounds on the volumes and distances shown. The interstitial volume ω T is computed from the relation where χ = (r, z) denotes coordinate points. No liquid can traverse the axis of symmetry and therefore that contribution vanishes. The negative sign in the last term reflects the fact that along the boundary γ , n T = −n. Along the conical boundary γ T , it is also the case that χ · n T = 0 and therefore the corresponding integral in Eq. represents the small volume exterior to γ , which is enclosed within the gray rectangular region with horizontal length r =â and vertical length z = z max . The variableâ denotes the radial coordinate corresponding to the turning radius of γ D . The variable z max is related to the minimum elevation of γ at r = 0 as described in the text.
vanishes. The integral over γ * also simply reduces to (2π/3) r 2 * [c 0 r * − z * ]. Evaluation of the integral along γ requires an analytic expression for the boundary function. The stationary form of Eq. (51), along with the observation that ∇ 2 ψ = 0 within ω liq and application of the divergence theorem yields 2 3 γ n · χrdγ = 1 2 γ n · ∇ψ2π rdγ where the quantity γ ∪ γ liq denotes the union of the boundaries enclosing the volume ω liq . Recall too that r = r sin θ. In the limit r → ∞, the velocity potential on γ liq approaches its asymptotic value ψ ∞ given by Eq. (69), which when substituted into Eq. (80) yields S liq n · ∇ψr 2 * sin θdθ The last step is obtained by noting that a relation obtained directly from integration of the Legendre differential equation. Substitution of Eqs. (80) and (81) into Eq. (79) results in the series solution ω T = a 1 π(1 + cos θ * ) − a 0 2π sin 2 θ * 3 r 3/2 * P 1/2 (− cos θ * ) where the domain truncation radius r * = r * sin θ. We now expand each term in powers of r in the limit where the domain truncation radius r * (or r) 1 such that θ * → θ 0 and c 2 = 0 from Eq. (76). After some straightforward algebra, we find Since all terms remain bounded as r * → ∞, the final result is given by The choice a 0 = 0, which from Eq. (74) implies that c 1 = 0, reduces to the result by Zubarev [4] for inertialess systems.
A more intuitive geometric interpretation of the relation in Eq. (85) can be derived by introducing the curved boundary γ D illustrated in Fig. 6(b) where We define a corresponding volume ω D , which encloses the region between the liquid-vacuum interface γ and the leading-order asymptote γ D such that Numerical evaluation of this expression yields lim r * →∞ ω D ≈ 1.093 a 1 − 0.118 a 2 0 . Equation (87), which places an additional constraint on a 1 and c 1 , or equivalently, a 0 and b 0 , establishes that arbitrary combinations of a 0 and b 0 are not all admissible.
A bound on the distance between the liquid top and the apex of the conical envelope can now be derived from Eq. (87). For simplicity, we restrict attention to the case c 1 < 0 and recall that c 0 = cot θ 0 < 0 as well. Consider a 044001-17 boundary γ everywhere bounded below by γ D , according to which c 0 r + c 1 r −1/2 ≤ h(r) ≤ c 0 r for r ≥ 0 and c 1 < 0. (88) Referring to Fig. 6(c), let the radial coordinate r =â define the turning point radius of γ D where the curve attains its maximum amplitude such that which then yields the coordinate valueŝ It follows immediately from evaluation of Eq. (88) on the interval 0 ≤ r ≤â that the maximum elevation in h must occur above the maximum elevation of γ D atâ, and its minimum elevation below z max such that

V. COMPLETE NUMERICAL SOLUTION BY BOUNDARY INTEGRAL PATCHING TECHNIQUE
The asymptotic solutions for the velocity potential, electric potential, and interface shape derived in Sec. IV are valid only at large distances far from the liquid tip. The complete solution to the self-similar system of nonlinear equations describing electrified flow requires a fully numerical approach. In what follows, we outline the boundary integral patching technique [43] used for this purpose, in which the asymptotic solutions derived in the previous section are used to construct the far-field boundary conditions. The interested reader may wish to consult the original study by Lennon and co-workers [50], which lay the ground work for formulating numerical solutions to problems involving nonlinear, free surface, axisymmetric flows, as well as more recent studies by Lister and co-workers [42][43][44], who examined self-similar capillary pinchoff of an inviscid fluid. The textbook by Pozrikides [51] is also an excellent source on the ins and outs of boundary integral techniques. Fig. 7 are illustrations of the geometry and boundary conditions used to obtain the numerical solutions. According to Green's theorem, since the velocity potential ψ and electric potential φ satisfy a harmonic equation, they can be represented by boundary integrals of the form

Shown in
where the axisymmetric Green's function g(χ , χ) is defined symbolically by where ϑ is the cylindrical azimuthal angle, The function g(χ , χ), which represents the potential strength at χ from a ring source at χ , cannot be expressed in terms of simple elementary functions. However, it can be expressed in terms of the complete elliptic integrals of the first and second kind [50]. For numerical integration, the boundary curves γ and γ liq enclosing the fluid volume ω liq , and γ vac and γ enclosing the vacuum volume γ vac , are discretized by a finite number of straight-line segments where β (expressed in radians) represents the interior angle between two adjacent segments. Along smooth portions of these boundaries β = π while at the domain truncation point β = π/2. While in most boundary value problems involving semi-infinite domains, contributions from boundary elements in the far field can be neglected since they vanish identically, this is not the case here. For example, when evaluated along the boundary segment γ liq , the contribution which does not vanish but increases with r * . Similarly for the contribution to the electric potential from the integral containing g and ∂φ/∂n. It is therefore critical that the contributions from the far-field boundaries γ liq and γ vac be included in the computations in order to obtain accurate solutions. The boundaries γ , γ liq , and γ vac are discretized by a set of knots interpolated by quintic splines. Quadratic Lagrange basis functions along the spline arc length are used to approximate ψ, ∂ψ/∂n, φ, and ∂φ/∂n using the Lagrange nodal values along the boundary. For integrals involving Green's function g or its derivative that did not contain a singular point within the boundary element and therefore behaved regularly, integration proceeded by Gauss-Legendre quadratures. When the singular point occurred at either end of a boundary element, integration relied instead on logarithmic-weighted quadratures since the normal derivative of g has a logarithmic singularity.
The illustrations in Fig. 7 indicate the geometry and boundary conditions used in the numerical computations. The boundary conditions applied along the curves γ , γ liq , and γ vac are specified as follows. The velocity potential ψ is required to satisfy a Neumann condition on γ given by the steady form of the kinematic condition in Eq. (51) and a Dirichlet condition on γ liq prescribed by the asymptotic series in Eq. (69) up through k = 4. The electric potential φ on γ is required to satisfy the Dirichlet equipotential condition in Eq. (59a) and a Neumann condition ∂φ/∂n on γ vac obtained from Eq. (70) up through k = 4. The domain truncation radius is situated very far from the liquid tip-simulations were carried out with |χ * | = 50, 75, 100. Although we did not conduct a convergence study, the results shown using |χ * | = 50 reproduced the same values obtained with larger radii.
Functions and their first derivatives evaluated on boundary elements adjacent to a corner can exhibit discontinuous behavior-as such, element distributions sharing a corner node are no longer appropriate [51]. The corner point χ * joining the boundaries γ and γ liq or γ and γ vac in Fig. 7 is a boundary corner where β = π/2. The discontinuity stems from the different boundary conditions on the linear segments joined by that point. In evaluating the corner point where β = π/2, we adopt the numerical approach described in Ref. [52] in which an extra constraint is imposed based on a local approximation (finite difference) to the field variables ψ and φ and their derivatives just prior and subsequent to the point.
The numerical scheme proceeds as follows. Given a parameter pair (a 0 , b 0 ) specified by Eqs. (69) and (70), a corresponding trial liquid interface function and its first two derivatives are obtained from Eq. (71) and evaluated at the domain truncation point χ * . The matrix equations resulting from evaluation of the discretized integral equations for ψ and φ in Eqs. (94) and (95), subject to the boundary values set by the asymptotic forms of ψ and φ and their first derivatives up through the k = 4 terms, are solved by QR decomposition to obtain the nodal values for the velocity potential ψ and the electric field strength ∂φ/∂n on γ . A Newton-Raphson method is then applied iteratively to adjust parametrization of the boundary γ until ψ and φ satisfy the time-independent form of the surface Bernoulli equation in Eq. (50). The Jacobian matrix at each Newton step is numerically approximated by perturbations to γ along the boundary normals.

A. Results showing subconical, superconical, and mixed-conical flow
Plotted side by side in Fig. 8(a) in cylindrical coordinates are results from numerical solutions for two values of a 0 (leading-order coefficient in the velocity potential 044001-19 ψ) and the same value of b 0 (leading-order coefficient in the electric potential φ). The left image in Fig. 8(a) shows subconical flow, where the liquid surface lies below the Taylor cone envelope and where the liquid in the tip and bulk region advance toward the conical apex in unison, as sketched on the right in Fig. 5(a). As discussed earlier, inertial effects generate a velocity field with angular dependence-in this example, the closer the coordinate is to the central axis, the more vertically aligned is the flow. Closer inspection of the pressure contours just below the liquid surface reveals very small interface oscillations due to capillary effects. Keller and Miksis [53] first identified similar oscillations in simulations of two-dimensional fluid wedges recoiling in vacuum. They reported how initial retraction of the wedge tip generates capillary-inertial waves that propagate away from the vertex with diminishing magnitude. Sierou and Lister [44] confirmed similar phenomena in boundary integral simulations of postcapillary pinchoff in a fluid thread. We leave it to future work to establish whether the wavelength of surface oscillations associated with electrified capillary-inertial waves obeys the same scaling with distance and time as the capillaryinertial waves originally reported by Keller and Miksis. We surmise this will be the case based on the fact that in this current study, the capillary, Maxwell and kinetic pressure in the surface Bernoulli equation scale similarly to leading order.
The right image in Fig. 8(a) depicts a velocity potential describing the original solution by Zubarev [4] for inertialess flow. Here, the asymptotic velocity potential scales as 1/r and not r 1/2 (i.e., a 0 = 0 and a 1 = 0.23). Unlike the ballistic flow to the left, the majority of the flow near the tip is now dominated by the radial dependence of the asymptotic potential, which generates isobaric curves that are mostly spherically symmetric. The velocity and interface field manifest superconical flow, where the interface lies above the conical envelope and the fluid in the tip and nearby bulk recede from the conical apex in unison, as sketched in Fig. 5(b). The interface again exhibits small-scale oscillations due to electrified capillary waves. Given negligible inertial effects, the velocity field near the interface seems to be more strongly influenced by the surface stagnation points, which punctuate regions of flow reversal.
Plotted side by side in Fig. 8(b) in cylindrical coordinates are results for equal values of a 0 and two different values of b 0 . The flow profile in the left image corresponds to a mixed-conical configuration, where the liquid near r = 0 advances toward the conical apex from below while the nearby bulk fluid recedes from that point, as sketched in Fig. 5(c). The choice b 0 = 0 enforce no electric field-the flow profile therefore depicts liquid motion resulting solely from capillary and inertial forces. This field is reminiscent of behavior subsequent to capillary pinchoff after strong droplet elongation as reported by Sierou and Lister [44]. The right image in Fig. 8(b) corresponds to a superconical configuration with tip and interior fluid retraction toward the conical apex, as shown in Fig. 5(b). The change in sign of a 0 from that in Fig. 8(a) generates a different tip shape that protrudes beyond the conical envelope, indicative of a postsingularity event.
The numerical solutions yield a multiplicity of selfsimilar flow configurations depending on the input parameter values (a 0 , b 0 ). Since b 0 controls the magnitude of the leading-order term in the asymptotic electric field, larger values of b 0 generate more conical-like tips. And while in the self-similar frame the region about the liquid maximum at r = 0 resembles a spherical cap, we note that in the laboratory frame, the corresponding tip will appear quite conical since the scaling of the spatial coordinates in Eq. (43) dictates that the radial and vertical extent of this region must scale as τ 2/3 . The cap therefore shrinks rapidly as τ → ∞ while maintaining the same asymptotic slope given by Taylor angle. This can be seen by comparing the interface shape in Fig. 9(a)plotted in self-similar coordinates, which matches exactly a case study by Burton and Taborek [36], with their snapshots in real time plotted in the laboratory frame shown in the inset image of Fig.  1(b) in Ref. [36]. What appears like a rounded cap in the self-similar frame corresponds to a fairly sharp conical tip in the laboratory frame.

B. Comparison with Burton and Taborek [36] simulations
For decades, researchers have been interested in quantifying the influence of charge transport on the shape and emission of progeny drops during Coulombic fission. For a perfectly conducting mass like a liquid metal, the electrical conductivity is assumed infinite and any excess charge on the surface redistributes instantaneously in order to maintain equipotential conditions. This is in line with the estimate in Sec. II C. The charge density and electric field therefore depend purely on geometry and not the advection of charge. In a seminal paper in 2011, Burton and Taborek [36] simulated the deformation accompanying Coulombic fission of an isolated inviscid droplet of perfectly conducting liquid of density ρ liq of sufficiently high surface charge density embedded within an exterior fluid of lower of zero density. Their numerical simulations revealed formation of self-sharpening tips at opposite ends of a droplet in which the apical values of the surface charge density and interface curvature underwent strong divergence in finite time. As expected, progeny droplets are only observed with liquids of finite conductivity and not perfectly conducting liquids. Significantly, they uncovered robust power-law growth spanning an incredible 12 decades in time in which the apex curvature scaled as 0.604τ −2/3 and the surface charge density scaled as 0.925τ −1/3 .
We compare our results to their simulations (not shown here) by magnifying the inset image of Fig. 1(b) in Ref. [36] showing multiple snapshots in time of a liquid tip evolving into a dynamic cone. The boundary data points for 20 such snapshots are extracted and transformed to self-similar coordinates. Although no length scale accompanied the plot, we determine the correct overall isotropic scale by conducting a least-squares fit between their transformed data sets and our results for coefficient values a 0 = −2.37 and b 0 = 1.757. This comparison yields a value of h = −0.608 and |n · ∇φ| r=0 = 0.922. Shown in Fig. 9(a) are the self-similar data points extracted from Burton and Taborek [36] superposed on the interface shape from our boundary integral patching technique. The agreement is excellent. Plotted in Fig. 9(b) are results extracted from our boundary integral simulations showing the behavior of the interface height, surface velocity potential, and surface electric field strength with increasing radial distance r, confirming the asymptotic behavior in Eqs. (69)-(71).

VI. CONCLUSION
In this work, we examine in detail the dynamic evolution of an axisymmetric protrusion in an electrically stressed liquid metal. Estimates provided indicate why a liquid metal is well approximated as a perfectly conducting inviscid fluid and so the electrohydrodynamic analysis takes full advantage of the Bernoulli equation and kinematic boundary condition for specifying interface evolution. Based on previous work by Zubarev and co-workers, as well as interface dynamics known to occur in capillaryinertial systems just prior or subsequent to pinchoff, the 044001-21 analysis predicates that protrusion growth is governed by a self-similar blowup process in which the kinetic, capillary, and Maxwell pressure contribute equally to leading order. The complexity of the coupled set of equations and boundary conditions preclude analytic solution except in the far field. Interestingly, these competitive forces allow asymptotic solutions to be found which in the far field to leading order adopts the shape of a Taylor cone with an interior angle given by the classic Taylor angle. This shape, however, when transformed back to the laboratory frame, defines a dynamic cone and not the conventional static Taylor cone. The corresponding asymptotic solutions for the velocity potential, electric field potential and interface function represent a two-parameter family of self-similar solutions prone to blowup in finite time. By invoking timereversal symmetry inherent to inviscid flow, we illustrate liquid formations describing subconical, superconical, and mixed-conical flow in which the interface defining the liquid-tip region appears below, above, or in a mixed state with regard to the Taylor conical envelope.
To obtain the complete solution, valid throughout the near-and far-field domain, we implement a boundary integral patching technique. In contrast to conventional boundary integral calculations on semi-infinite domains in which boundary contributions typically decay to zero and are neglected, the far-field boundary contributions in the vacuum and liquid domain are non-negligible. The formulation of the patching technique therefore directly incorporates information from the far field through the series expansions describing the asymptotic solutions. The accuracy of this approach is confirmed by direct comparison to a numerical study by Burton and Taborek [36] for a system involving Coulombic fission of droplets, which shows excellent agreement.
The results of our numerical simulations highlight the crucial influence of inertial forces in the liquid apical region as it undergoes rapid acceleration or deceleration toward or away from the blowup point. Different parameter choices for the two-parameter family of asymptotic solutions reveal a multiplicity of fluid configurations that exhibit not only apical sharpening but tip bulging, tip separation flow, clusters of interface stagnation points from capillary wave phenomena, and receding interface profiles reminiscent of recoil after capillary pinchoff. The resulting interface shapes confirm that the local interior half-angle within the liquid tip varies with position and can be larger or smaller than the classic Taylor angle defining the asymptotic envelope.
Experimental observations during actual LMIS operation often report phenomena such as tip pulsation, droplet emission, liquid recoil and collapse. We wonder if some of these processes may be related to progressive amplification of the small interface oscillations observed in our numerical solutions. This line of inquiry requires a separate stability analysis, which can be conducted now that the self-similar base state solutions are known. Such an analysis will help further establish whether the electric capillary wave trains observed can trigger nonlinear instability or whether the asymptotic constraint set by the Taylor conical envelope represses such behavior.

ACKNOWLEDGMENTS
S.M.T. acknowledges valuable discussions with Dr. Colleen Marrese-Reading and Dr. James E. Polk at the NASA Jet Propulsion Laboratory about liquid metal ion sources, and the assistance of Dr. Peter Thompson with the LIS2T computing cluster used in this work. S.M.T. also wishes to acknowledge the three journal referees for a close reading of the manuscript and helpful suggestions provided.