Enhancement of droplet ejection from molten and liquid plasma-facing surfaces by the electric field of the sheath

Maintaining the stability of a liquid surface in contact with a plasma is of crucial importance in a range of industrial and fusion applications. The most fundamental feature of a plasma-surface interaction, the formation of a highly-charged sheath region, has been neglected from the majority of previous studies of plasma-liquid interactions. This paper considers the effect of the electric field of the sheath on the ejection of micron-scale droplets from bubbles bursting at the liquid surface. A numerical simulation method, based on the ideal electrohydrodynamic model, is introduced and validated against the well-known Taylor cone theory. This model is then used to include the electrical effects of the sheath in simulations of bubble bursting events at a plasma-liquid interface. The results show a significant enhancement in droplet ejection at modest electric fields of between 10% and 20% of the critical field strength required for a solely electrohydrodynamic instability. This finding is in good qualitative agreement with experimental observations and its importance in a wide range of fusion and atmospheric-pressure plasma-liquid interactions is discussed. The inclusion of sheath physics in future studies of plasma-liquid interactions is strongly advocated.

remain an excellent long-term proposition for handling the immense heat fluxes of a fusion plasma [13].
A clear understanding of the processes which cause droplet ejection from liquid metals is required in order to optimise the design of plasma-facing components in magnetic fusion devices. Hassanein and Konkashbaev first mentioned the bursting of bubbles on the surface of liquid metals as a mechanism for droplet ejection from melted surfaces in tokamaks [14,15]. These bubbles might originate from boiling of the liquid or the absorption of gases which can create blisters on the surface of plasma-facing metal components [16]. The basic mechanism is illustrated in figure 4 of [15]: the film covering the surface of the bubble bursts, which can release some very fine droplets [17], and the subsequent collapse of the bubble forms a rising jet of liquid which pinches off, according to the Rayleigh-plateau instability, to form a spray of droplets. This process is also observed, for instance, in the bursting bubbles at the surface of a glass of champagne; the ejected droplets are alleged to 'complement the sensual experience of the taster' [18].
The comprehensive parametric study of Duchemin et al [19] indicates that 10-100 µm bubbles at metallic surfaces eject droplets with velocities of 10-100 m s −1 and radii around one-tenth of the initial bubble radius. Shi, Miloshevsky and Hassanein simulated the specific cases of tungsten and aluminium surface bubbles and found consistency with Duchemin's more general study [20]. The agreement between these predicted droplet velocities and those observed in de Temmerman's tungsten and aluminium melting experiments on Magnum-PSI [21] provides strong evidence that the bubble bursting mechanism is responsible for high-velocity droplet production in tokamaks. However, de Temmerman provides a further intriguing observation: the formation of liquid droplets is dependent on the electric potential, and hence the electric field, of the metal surface. The ejection rate of aluminium droplets increased by a factor of 10 when the surface was at floating potential, rather than being grounded, and further increased by another factor of 10 when the target was biased at −40 V.
These experimental results can be investigated further by including electrical effects in simulations of plasma-liquid interactions. Langmuir's observation of the highly-charged sheath region at a plasma-surface boundary was one of the earliest discoveries in plasma science [22,23]. This sheath contains a large electric field which draws ions out of the plasma in order to give an equal flux of ions and electrons, which move much faster than the ions, to the plasma-facing surface. Any credible theory of plasma-liquid interactions must account for the large electric fields and particle fluxes of sheath region but their role has been broadly neglected until recently [24][25][26]. These studies used a model of cold ions and Boltzmann electrons to account for the charge density and electric fields at a plasma-liquid interface with a planar initial configuration. However, the surface of a liquid in contact with a plasma is unlikely to maintain an idealised planar shape due to effects such as liquid/plasma flow, surface rippling or bubble formation. The present paper builds on these early studies by considering numerical simulations of the effect of the sheath's electric field on droplet ejection from the bubble bursting mechanism where the initial surface geometry contains a spherical void. The enhancement in droplet ejection rates in these simulations is in good qualitative agreement with the findings of de Temmerman's experiment and strongly advocates the inclusion of electric fields in future studies of plasma-liquid interactions.

Ideal electrohydrodynamics (EHD)
If the lengthscale of surface deformations is smaller than the Debye length then the electric field close to a conducting surface is dominated by the free charges on the surface. The charge density of electrons and ions in the plasma is insufficient to modify the electric field significantly over these short, sub-Debye distances so the electric field close to the surface can be well approximated as a vacuum field. The fluid equations describing this system are those of a charged liquid in a vacuum; these equations have been referred to as the ideal electrohydrodynamic (EHD) equations elsewhere [27] and are restated here. Previous studies have shown that a complete plasma model, which includes ion and electron fluxes onto the liquid surface, recovers the ideal EHD behaviour in the short-perturbation regime (see [25], figure 3). The fact that the plasma does not appear explicitly in the equations here does not mean that the plasma has completely disappeared; it is still responsible for the accumulated charge on the liquid surface and hence is the original source of the electric field.
The fluid equations of an incompressible liquid, with viscosity η and density ρ, in contact with a vacuum are given by the incompressibility condition and Navier-Stokes equation where p and v denote the liquid pressure and velocity, respectively, and g is the acceleration due to gravity. The electric potential φ in the vacuum region is given by Laplace's equation which produces the electric field as E = −∇φ. The conducting liquid provides an equipotential boundary condition for Laplace's equation. The pressure at the surface of the liquid is determined by the jump condition where γ is the surface tension value and κ is the local curvature of the surface. As noted previously, the omission of the ion ram pressure and electron thermal pressure in this description of plasma-liquid interactions is valid provided that deformations of the liquid surface occur over lengthscales shorter than the plasma Debye length. The ideal EHD equations, equations (1)-(4), have been solved in numerous different situations. Among the most celebrated solutions are the Rayleigh criterion for disruption of charged droplets [28], the static Taylor cone configuration of a liquid surface in an electric field [29], and Melcher's linear analysis of the dispersion and stability of EHD surface waves [30]. The latter result gives the dispersion relation for a wavelike perturbation, with frequency ω and wavenumber k, on the surface of a conducting liquid under the application of a uniform electric field of strength E 0 as The perturbations are unstable when ω 2 < 0 which gives the marginal condition for instability This stability criterion is the same for both plane-wave and cylindrical, i.e. Bessel function, surface perturbations [31]. Note that this stability criterion reduces to that of the standard Rayleigh-Taylor instability; setting E 0 = 0 and reversing the sign of g gives the marginal condition for a Rayleigh-Taylor instability of a liquid layer, with density ρ and surface tension γ, placed above a field-free vacuum region as ρg > γk 2 .

Simulation method
A code to solve the ideal EHD equations, equations (1)-(4), has been developed recently and applied to the dynamics of charged and rotating droplets [27]. This original code is freely available at www.github.com/joshholgate/iEHD. The workings of this simulation method are briefly outlined here; full details can be found in the original publication and source code. A few adjustments to the boundary conditions and initial configuration of the liquid surface allow simulations of the bubble bursting droplet ejection process with the inclusion of electric fields. The simulation is based upon the finite-difference, level-set (FDLS) method which has emerged as a popular scheme for simulating EHD flows over the past decade [27]. The position and motion of the liquid surface is captured by defining an additional scalar field ψ, known as the level-set function [32], which is positive inside the conducting liquid 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 sign across the interface. The level-set function is convected with the velocity field according to and the surface of the liquid moves with it. The level-set function allows straightforward calculation of the normal vector and curvature of the free surface using the formulae n = −∇ψ/|∇ψ| and κ = −∇ ·n. The ideal EHD equations, together with level-set interface tracking equation, are discretised using the finite-difference method on a uniform staggered grid configuration in axisymmetric cylindrical coordinates and solved using the Chorin projection method as outlined in [27]. No-slip boundary conditions are utilised at the points where the liquid comes into contact with the edges of the simulation region.

Simulation validation
The method described above is used to simulate the motion of a liquid surface in an electric field and the results compared with Melcher's linear stability results, equations (5) and (6), in order to verify that the code is producing physicallyrealistic results. The simulation domain is a cylinder of radius 10 µm and height 20 µm which is half-filled by a conducting liquid with the same viscosity, density and surface tension as tungsten: η = 7 mPa s −1 , ρ = 17 600 kg m −3 and γ = 2.5 N m −1 [20]. The Debye length in tokamak plasmas is larger than this radius (appendix A in [32]), so the dimensions of the simulation correspond to the ideal EHD regime described in section 2. A fixed potential difference is applied between the top of the cylinder and the liquid surface which has an initial height perturbation given by a zeroth-order Bessel function with peak-to-trough amplitude 0.5 µm. The radial wavenumber is k = 3.8217 × 10 5 m −1 which corresponds to the first minimum of the Bessel function occurring at r = 10 µ m. The linear perturbation theory predicts, from the dispersion relation in equation (5), that this perturbation grows exponentially at a rate of where gravity is neglected. Growth rates can be extracted from the simulations and compared with this prediction in order to make a quantitative comparison between theory and simulation. The growth of a Taylor cone from this initial setup is shown in figure 1 for an electric field given by 0 E 2 0 /γk = 1.1. The opening angle of the flat slopes of this cone at 46 µm is measured as 46 ± 2 • , in close agreement with Taylor's theoretical value of 49.3 • [29]. Subsequent frames show the extrusion of the cone tip into a fine jet of droplets which replicates the process seen in experimental observations [34]. The exact mechanism of droplet pinch-off is likely to be adversely affected by the grid spacing used in the simulation as the apex of the cone approaches the size of the numerical grid. As such the size of the ejected droplets tends to be determined by the grid resolution, which places a lower limit on the size of the simulated droplets, rather than by any physical effects.
The simulation is run for various electric field strengths as characterised by the factor 0 E 2 0 /γk. The amplitude of the surface perturbation increases exponentially during the early stages of each simulation, as predicted by the linear perturbation theory, so growth rates can be extracted by least-squares fitting exponential curves to these growing amplitudes. The results are shown in figure 2 along with the growth rates predicted by equation (8). The simulation results lie on a line with the formula 0.702( 0 E 2 0 /γk − 1) 1/2 whereas the linear perturbation theory predicts that the factor before the square root should be exactly one. Therefore the growth rates observed in the simulation are 30% lower than predicted the linear perturbation theory; this discrepancy is tentatively assigned to the inclusion of viscosity in the simulations. However the marginal stability limit, 0 E 2 0 = γk, is accurately reproduced by the code which gives a strong indication that the simulations are producing physically-realistic results.
At this point it is natural to wonder whether the basic EHD instability described by these simulations is solely responsible for the ejection of droplets from plasma-facing liquid metal surfaces in tokamaks. The answer is almost certainly no. The critical electric field required to make a tungsten surface become unstable, determined according to equation (8) as 0 E 2 c = γk, is over 100 MV m −1 for the 10 µm radius perturbation shown in figure 1. The electric field in the sheath of a typical tokamak plasma, with electron temperature T e = 50 eV and Debye length λ D = 20 µm (appendix A in [33]), is approximately where the sheath width is on the order of λ D and k B T e /e provides an approximate value for the potential drop across the sheath [23]. This electric field strength is clearly insufficient to initiate a short-wavelength EHD instability of the surface. If the initial perturbation has a larger wavelength than these simulations, and hence a smaller value of k, then the marginal stability criterion predicts a smaller critical electric field strength which could feasibly be experienced in a typical tokamak sheath. However such a perturbation would require surface deformations on a lengthscale larger than the Debye radius and the ideal EHD equations would no longer apply. The linear stability theory of the full system of plasma-liquid fluid equations associated with these large scale instabilities is considered in [24] which shows that ion bombardment can stabilise long-wavelength instabilities against the electric field. On the other hand some technological uses of plasmaliquid interactions involve the application of a large electrical voltage, which far exceeds the floating value, to a liquid cathode. Vyalykh provided the first experimental study of electrically-driven oscillations and instabilities of liquid cathode surfaces in low-pressure glow discharges [35]. Bruggeman later applied the EHD marginal stability condition, equation (6), to argue that water cathodes are always electrically-unstable for typical plasma glow discharge param eters and that this instability is manifest in observations of the plasma splitting into a filamentary structure [36]. The detailed experimental studies of Shirai provide images and analysis of Taylor cone formation and droplet ejection in an atmospheric negative corona discharge [37,38]; plumes of ejected droplets are also reported in the atmospheric-pressure solution-cathode glow discharge of Schwartz, who notes the similar characteristics of these micron-scale jets to those produced by Taylor cones [39]. Finally, the importance of charge accumulation on plasma-facing liquid surfaces is further demonstrated by the recent experiments of Dubinov who observed unusual deformations of a liquid surface in a capillary tube and between thinly-spaced capillary plates when a spark discharge was applied to the surface of an ionic solution [40].

Field-enhanced droplet ejection from bursting bubbles
The linear stability theory and simulations considered so far make the assumption that the surface of the liquid is almost flat before the instability begins. This idealised initial configuration is not always accurate; in particular the bubble bursting mechanism for droplet formation from a liquid surface, as introduced in section 1, begins with the air cavity left at the surface of the liquid after its covering thin film has burst. Previous works have identified this mechanism as a leading contender for droplet ejection from molten surfaces in tokamaks but have not yet accounted for the role of the electric field  in the sheath. The significance of this omission is indicated by the experiments of de Temmerman et al where the ejection of droplets from melted aluminium and tungsten targets in the Magnum-PSI device was attributed to bubble bursting [21].
The ejection rate of aluminium droplets increased by a factor of 10 when the surface was at floating potential, rather than being grounded, and further increased by another factor of 10 when the target was biased at −40 V. The axisymmetric simulation described in the previous section can be used to investigate the role of electric fields in enhancing the bubble bursting mechanism and to explain the findings of these experiments.
The evolution of an aluminium surface with a spherical cavity, of radius 10 µm and with its centre positioned 7.5 µm below the liquid surface, is shown in the top row of figure 3. The liquid region is 4 bubble radii (40 µm) deep and 4 bubble radii from its centre to its outer radial edge. A small potential difference is applied between the liquid and the top of the simulation domain which is located 8 bubble radii (80 µm) above the liquid surface. Gravity is ignored in this simulation. The electric field at the start of the simulation has some enhancement at the corners where the bubble meets the surface while the electric field inside the bubble is almost zero due to the shielding effect of the bubble cap. The bubble collapses rapidly after the simulation is initialised which causes a jet of liquid to shoot upwards but, before it can pinch off into a droplet, surface tension pulls this jet back down. This contradicts the findings of Shi, Miloshevsky and Hassanein who did observe pinching-off of the liquid aluminium column [20]. This difference could be due to their smaller simulation domain, which is 3 bubble radii in depth and 2 bubble radii in radius, which concentrates energy into the upwards motion of the jet.
The second and third rows of figure 3 show the development of the bubble-bursting process with the same initial conditions as before but with a larger potential difference applied between the liquid and the top of the simulation domain. If the surface was perfectly-flat then the electric field generated by this potential difference would still be insufficient to initiate an EHD instability with a wavelength equal to the bubble radius; the ratio of the applied field to the critical field, E c , required for this EHD instability is indicated for each simulation. The electric field strength is significantly enhanced at highly-curved points on the liquid surface as evidenced by the small nodules which are pulled out of the jets in the 1 µs frames. The second row shows the electric field removing a single large droplet from the rising jet of liquid whereas the jet is elongated before being pinched-off in the third row. The resulting column of liquid in the bottom-right frame then becomes unstable to a Rayleigh-plateau instability where the liquid column pinches-off into a chain of droplets. Note that the charge on these ejected droplets will be determined by electron currents within the plasma, with a characteristic charging timescale of 10 −9 s [41], rather than any charge relaxation effects occurring over the much longer bubble-bursting timescale of 10 −6 s. The ejection of liquid droplets from the bubble bursting mechanism is clearly facilitated and enhanced by sub-critical electric fields in accordance with de Temmerman's experimental observations [21]. The highly-curved surfaces created by the bubble bursting mechanism lead to significant enhancements in the electric field strength over that of a planar surface; such geometric effects can be just as important in determining the stability of a plasma-facing surface as the raw strength of the sheath's electric field.
The simulation results discussed here are also important in non-fusion atmospheric plasma discharges. A direct application of this mechanism occurs in plasma welding and cutting devices which use the high temperatures of the plasma to melt the metal surface [42]. As the metal surface boils it will produce surface bubbles which, together with the large electric field of the plasma arc, provide the same initial conditions as shown in figure 3. The removal rate of melted material in the efflux plasma is crucial for decreasing surface roughness and improving the quality of the cut and remains a topic of active research [43].
Furthermore it is noted by Bruggeman that 'plasmas interacting with liquid water can drive instabilities leading to formation of bubbles in solution in addition to affecting the micron-sized bubbles already present' [2]. As such the bubble-bursting mechanism will play a role in the injection of micrometre droplets in atmospheric-pressure low temperature discharges interacting with liquids. However the electric field in these atmospheric plasmas is unlikely to reach the strength required for significant enhancement in droplet ejection rates: a plasma with electron density of 10 19 m −3 and temperature of 1 eV has a field strength of approximately 0.4 MV m −1 , as given by to equation (9), while the critical field strength of an electrohydrodynamic instability, 0 E 2 c = γk, is 20 MV m −1 for a 10 µm water surface perturbation with γ = 0.072 N m −1 . The typical field strength, at 2% of the critical field, falls below the 10%-20% level which is required for enhanced droplet emission in the simulation results of figure 3.

Conclusions
Maintaining the stability of plasma-facing liquid surfaces is of crucial importance in a wide range of industrial applications and in magnetic confinement fusion devices where metallic components can suffer melt damage due to the high fusion heat loads. The melting and erosion of metal surfaces is also critical for the operation of plasma cutting and welding devices. The plasma region immediately adjacent to the liquid surface, known as the sheath, is highly-charged and exerts a considerable electric stress on the liquid. The electric field of the sheath has only recently been included in models of plasma-liquid interactions. However these electrical effects have been considered in isolation so far; this paper investigates their role in enhancing droplet ejection rates from the bursting of bubbles on the surface of the liquid.
The dynamics of micrometre-scale deformations of the plasma-liquid interface, such as those involved in the bubble bursting mechanism, can be resolved using the ideal electrohydrodynamic equations which are introduced in section 2. A code to solve these equations using the finite-difference, level-set method is briefly described in section 3 and validated against the linear stability theory for small perturbations of a nearly-flat liquid surface in section 4. Full simulations of the bubble bursting process, both with and without an electric field, are then described in section 5 and the results are displayed in figure 3. The number of droplets ejected from the surface increases with the strength of the electric field in accordance with experimental observations of the bubble bursting mechanism [21]. The electric field does not need to be particularly large to cause this enhancement; an electric field between 10% and 20% of the critical field strength required for a solely electrohydrodynamic instability leads to a noticeable enhancement in droplet ejection.
It is clear that the electric field strength of the sheath alone does not entirely determine the stability of a plasma-facing surface; the initial shape of the surface, and any geometric enhancement of the electric field, is vital in determining its behaviour. The example of field-enhanced droplet ejection from the bubble bursting mechanism demonstrates the importance of these geometric effects and explains the increase in droplet emission in experiments where the surface has a slight potential bias [21]. The inclusion of electrical effects in other known droplet-forming mechanisms, such as Kelvin-Helmholtz instabilities [12,44], should certainly be investigated further in the future.