Simulated dynamics of a plasma-sheath-liquid interface

The discovery of a highly-charged sheath region at the boundary between a plasma and a surface is one of the earliest and most important discoveries in plasma science. However sheath physics has almost always been omitted from studies of the dynamics of plasma-facing liquid surfaces which are rapidly assuming a pivotal role in numerous industrial and fusion applications. This paper presents full simulations of the plasma-sheath-liquid interface and finds good agreement with theoretical stability limits and experimental observations of cone formation and pulsed droplet ejection. Consideration of sheath physics is strongly encouraged in all future studies of plasma–liquid interactions.

Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence.
Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
( ) E e n n e k T exp . 3 B e 0 0 A linearisation procedure was performed in [17] in order to make these equations tractable. An alternative approach is to compute solutions directly with numerical simulations such as the C++ finite-difference code 'Cold Ions with Boltzmann Electron Relation' (CIBER) which is freely-available at http://github.com/ joshholgate/CIBER and is described in detail by [20]. This code uses axisymmetric (r, z) coordinates but is otherwise unrestricted to any particular geometry of the conducting surface; as such it can be used to generate solutions for the plasma sheath profile near a liquid surface with arbitrary deformations. The boundary conditions on the plasma sheath region are given by the Bohm condition far away from the liquid surface, implemented according to the specification in [20], appendix B, while the conducting liquid surface is kept at a fixed potential of f w .
The dynamics of the incompressible liquid are determined by where v, p, ρ and η are the liquid velocity, pressure, density and viscosity, respectively, and g is the acceleration due to gravity. The iEHD code, which solves equations (4) and (5) for a conducting liquid with a vacuum electric field applied to its surface, is available at http://github.com/joshholgate/iEHD and is described in detail by [21]. The original publication was concerned with charged droplets but the code has since been applied to a cylindrical container of liquid with a surface electric field [22]. This code uses the level-set method to track the motion of the liquid surface and can naturally handle any breakup of the surface [23]. The iEHD code also uses a finitedifference method on a uniform grid in axisymmetric (r, z) coordinates and can handle arbitrary deformations of the liquid surface. These similarities make it relatively straightforward to couple CIBER and iEHD together into a single simulation of a plasma-sheath-liquid interface. However a suitable pressure jump condition is required in order to combine the two codes. This is provided by [17]  where the s subscripts denote values at the surface between the plasma and the liquid andn is the surface unit normal which points into the plasma region. The terms on the right-hand side of this expression are the ion ram pressure, electron thermal pressure, electric Maxwell stress and the Young-Laplace capillary pressure with surface tension γ and curvature k = - ·n. Note that this is the same as the pressure jump condition in the original iEHD code except for the inclusion of the electron and ion terms. The combined simulation solves equations (4) and (5) for the liquid under the pressure gradients generated by surface tension, electric fields, ion ram pressure and electron thermal pressure at the liquid surface according to equation (6). The surface of the liquid evolves on a long timescale compared to the time taken for rapidlymoving electrons and ions to adjust to the new surface shape. These different temporal scales are accounted for by taking numerous steps in the solution of equations (1)-(3) to find a steady sheath profile at each time step of the liquid simulation. The level-set method is still used to determine surface tension forces and to track the motion of the plasma liquid interface. The combined code has been named, with no little lack of imagination, as CIBER-iEHD and is freely available at http://github.com/joshholgate/CIBER-iEHD together with all of the data generated for this paper.
The full plasma-sheath-liquid simulations are initially configured with a liquid occupying the bottom of a cylindrical domain with a depth of ten Debye lengths, 10λ D . A small zeroth-order Bessel function perturbation is applied to the liquid surface with its maximum at r=0 and its first minimum at the outer radius of the simulation domain. Varying the outer radius of the domain therefore determines the wavenumber k, and hence the dimensionless quantity kλ D , for the initial perturbations. The liquid is taken to have a water-like density, viscosity and surface tension of ρ=1000 kgm −3 , η=0.89 mPAs −1 and γ=0.072 Nm −1 .
The plasma sheath region occupies the upper portion of the computational domain to a height of 30λ D above the liquid surface. The ion mass of this plasma can be shown to have no effect on the dynamics of the system by adopting a dimensionless normalisation of the plasma equations as in [20]. However the plasma density and electron temperature are important in determining the behaviour of the plasma-liquid interface; the plasma Bond number provides a dimensionless measure of the plasma pressure to the capillary pressure of the liquid [17] and is set prior to running each simulation. A fixed potential drop, f w is applied across the sheath region. Increasing the width of the sheath region beyond 30λ D was found to lead to prohibitively long runtimes; this limit on the sheath width means that attention is restricted to potential drops of f -< e k T 1000 . The surface rapidly forms a cone-like protrusion but with a rounded top and a much shallower opening angle than the α=49.3°predicted by the static Taylor cone theory for a vacuum-liquid interface [24]. Furthermore the liquid is not stationary but circulates steadily in response to the pressure gradients generated by the plasma stresses at the liquid surface.
The stable surface shown in figure 1 is on the cusp of becoming unstable; a slight increase in the potential drop of the sheath to f =e k T 32 w B e leads to the formation of a jet of droplets as shown in figure 2. The simulation method imposes that these droplets are at the same potential as the liquid surface which is certainly true at the instant of their formation but which might change as they interact with the plasma and attain their floating potential. Once these droplets are ejected they appear to shield the remaining protrusion from the electric field in the sheath which leads to the collapse of this protrusion. The duration of the entire process is around 10μs for the water-like liquid parameters used in this simulation.
The critical wall potential f c , which causes a transition from stable protrusions to unstable droplet formation, is identified for a range of initial perturbation wavelengths and plasma Bond numbers of 0.01, 0.1 and 1. The results are plotted as points in figure 3 together with solid lines showing the linear stability thresholds of [17]. Additional dashed lines represent the ideal electrohydrodynamic theory for the stability of plasmavacuum interfaces [17,25].
Several conclusions can be made from figure 3. First, the relative importance of the electric stress, which drives the surface instability, to the stabilising effect of the electron and ion pressures can be inferred from the differences between the field-only theory (solid lines) and the full plasma-sheath-liquid system calculations (solid lines and points). In the short-scale, sub-Debye-length regime (kλ D 1) the inclusion of electrons and ions makes no difference to the stability threshold which indicates that particle pressure plays an insignificant role in the stability of such surfaces. However when the wavelength of the instability is on the Debye scale and larger (kλ D 1) the full sheath results have a higher stability threshold that the field-only theory; in this case the magnitude of the electron thermal pressure and, in particular, the ion ram pressure becomes comparable to that of the electric stress. Second, the difference between the linear perturbation theory (solid lines) and numerical simulations (points), both of which model the full plasma-sheath-liquid system, reveals that nonlinear effects lower the stability threshold of the surface: the linear theory flattens at normalised potential values of 12, 21 and 72, for plasma Bond numbers of 1, 0.1 and 0.01 respectively, whereas the equivalent values for the simulations are 2, 8 and 50. There are several possible causes of the quantitative discrepancy between theory and simulation. For example the linear theory in [17] is developed for plane-wave perturbations whereas the simulations use axisymmetric cylindrical geometry. The marginal stability condition is identical for plane-wave and cylindrical perturbations of the ideal electrohydrodynamic liquid-vacuum interface [25,26] but this is no longer true when the electron and ion effects are included. The simulations also require a finite-sized perturbation to initiate the instability, rather than the infinitesimal surface displacement of the linear theory, which could remove the system from a local equilibrium configuration. Most significant, however, is probably the difference in , showing the formation of a conelike protrusion, ejection of a burst of droplets, and subsequent collapse of the protrusion. The duration of the entire process is around 10μs for the water-like liquid parameters used in this simulation. Figure 3. Comparison of the critical wall potential f c , which causes a transition from stable to unstable surfaces, given by the ideal electrohydrodynamic perturbation theory [17,25] (dashed lines), plasma-sheath-liquid perturbation theory [17] (solid lines) and numerical simulations of the full plasma-sheath-liquid system (points). equilibrium shape. The linear theory took a perfectly-flat surface as the unperturbed equilibrium while the simulations can also handle conical equilibrium shapes such as that of figure 1. This new equilibrium shape could become unstable to linear or nonlinear perturbations, both of which are included in the numerical simulation, at a lower threshold than the flat surface.
The recent experiments of Shirai et al have imaged the formation of corona discharges on the surface of a Taylor cone which are accompanied by the ejection of fine droplets from its tip [27,28]. The formation of droplet plumes is also reported by Schwartz et al [29]. It is difficult to make exact comparisons between Shirai's experiments and the simulations reported here but some important similarities can be identified.
Shirai's most recent experiment used a nozzle of diameter 1.5mm filled with an ionic solution as a negative electrode and placed it beneath an upper plate electrode. A DC potential difference of up to 15kV was applied between these two electrodes in order to initiate a plasma in the atmospheric air above the liquid cathode. Various compounds were dissolved in the water: sodium dodecyl sulphate to lower its surface tension and sodium chloride to increase its conductivity. When the applied voltage reached 4.52kV the liquid surface deformed into a Taylor cone [28] and produced a fine jet of droplets. Assuming typical electron temperatures and densities of 1eV and 10 19 m −3 , for a humid air discharge [2], gives a normalised sheath potential drop of f »e k T 4500 w B e and the Debye length as 2.4μm. Furthermore the plasma Bond number is roughly » -Bo 10 P 4 when the lowering of surface tension due to the surfactant is accounted for. These values correspond well to the predicted theoretical stable-to-unstable transition [17].
The cones formed in Shirai's experiments also lead to the formation of droplets and detection of current pulses at the upper plate electrode. These current pulses, which are displayed in [28], figure 4(a), are accompanied by a sharp increase in the light intensity from the discharge. The numerical simulation in figure 2 provides one possible explanation for these current pulses; the ejection of a burst of droplets followed by the collapse of the protrusion over a timescale of around 10μs which roughly matches the duration of the experimental current bursts. A slow recovery of the protrusion back into a cone-like shape over the following 100μs would then account for the gaps between the experimental current pulses.
Two further experimental verifications of the electrohydrodynamic stability theory can also be anticipated; both have already been performed but currently lack the numerical data to make a direct comparison with figure 3. The first setup is described by Allen and Magistrelli [30] where a mercury cathode was employed to supply high electrical currents to a low-pressure helium discarge. A protective umbrella was included in the apparatus to shield the discharge from jets of liquid emanating from the unstable liquid cathode. A second, more unexpected, electrohydrodynamic instability was also observed where a unipolar arc formed on the pool of condensed mercury which accumulated in the U-bend of a cold trap. Vyalykh et al provide a further experimental study of electrically-driven oscillations and instabilities of ionic liquid cathode surfaces in lowpressure DC air discharges [31]. These instabilities are recorded photographically but without quantitative data for the applied voltages which initiate the disruption of the liquid surface. Both of these texts describe their experiments in detail and, given that electrohydrodynamic instabilities were observed in both cases, appear to offer guaranteed prospects for interesting results. This paper has presented full simulations of a plasma-liquid interface which fully account for the strong electric fields and particle flows in the sheath region which forms adjacent to the liquid surface. Sheath effects, which have not been included in previous simulations of plasma-liquid interactions, are found to be crucial in determining the shape and stability of the surface. The electric fields in the plasma sheath can deform the surface into a cone-like protrusion which ejects bursts of droplets from its apex in accordance with experimental observations of technological plasma-liquid interactions. However the ions from the plasma, which are accelerated through the electric field of the sheath, provide a countering force on the liquid surface which can suppress the onset of electrohydrodynamic instabilities.
The simulations presented in this paper are concerned with the solely sheath-driven dynamics of a plasmaliquid interface; lower stability thresholds and increased droplet ejection will be discovered when sheath effects are included in other theories of plasma-liquid surface instabilities and splashing. Sheath physics, particularly the accumulation of surface charge and the generation of large electric fields, must be fully accounted for when considering the shape and stability of any plasma-facing liquid surface.