Numerical implementation of a cold-ion, Boltzmann-electron model for nonplanar plasma-surface interactions

Plasma-surface interactions are ubiquitous in the field of plasma science and technology. Much of the physics of these interactions can be captured with a simple model comprising a cold ion fluid and electrons which satisfy the Boltzmann relation. However, this model permits analytical solutions in a very limited number of cases. This paper presents a versatile and robust numerical implementation of the model for arbitrary surface geometries in cartesian and axisymmetric cylindrical coordinates. Specific examples of surfaces with sinusoidal corrugations, trenches, and hemi-ellipsoidal protrusions verify this numerical implementation. The application of the code to problems involving plasma-liquid interactions, plasma etching, and electron emission from the surface is discussed.Plasma-surface interactions are ubiquitous in the field of plasma science and technology. Much of the physics of these interactions can be captured with a simple model comprising a cold ion fluid and electrons which satisfy the Boltzmann relation. However, this model permits analytical solutions in a very limited number of cases. This paper presents a versatile and robust numerical implementation of the model for arbitrary surface geometries in cartesian and axisymmetric cylindrical coordinates. Specific examples of surfaces with sinusoidal corrugations, trenches, and hemi-ellipsoidal protrusions verify this numerical implementation. The application of the code to problems involving plasma-liquid interactions, plasma etching, and electron emission from the surface is discussed.


I. INTRODUCTION
Interactions between plasmas and surfaces form one of the foundations of plasma science with near-ubiquitous occurrence in industrial, technological, fusion, and space plasmas. The widespread use of plasmas in material processing is based on the high ion fluxes and electric fields in the sheath region between the surface and the bulk plasma which allow sputtering, etching, heating, charging, and coating of the surface. 1 Furthermore, the high-energy particles in the sheath are efficient chemical catalysts with applications such as CO 2 conversion, abatement of toxic waste gases, and the formation of functional materials such as nanotubes and nanoparticles. 2 The exposure of surfaces to immense heat fluxes in magnetic fusion devices remains one of the most critical issue facing the development of a fully functioning, power-generating device. 3 Finally, the interaction of objects in space, such as planets, artificial satellites, and dust particles, with sparse interplanetary plasma or solar wind provides phenomena such as the homogenization of cosmic dust, 4 spontaneous self-alignment of dust clouds, 5 lunar swirls, 6 and damage-induced limitations to the lifetime of the spacecraft. 7 This list is, of course, by no means exhaustive but it clearly illustrates the extensive occurrence and economic importance of plasma-surface interactions. Current research topics in the field of plasma-surface interactions include the growth of large-scale graphene sheets, 8 the effect of electron emission on the plasma-sheath structure, 9 and plasma interactions with liquids. 10 This work is motivated, in particular, by the latter topic which finds applications in nanoparticle synthesis, catalysis of chemical reactions, material processing, water treatment, sterilization, and plasma medicine. 10 Existing models of plasma-liquid interfaces have, until very recently, 11 neglected all effects involving charge separation and electric fields in the plasma.
The conventional methods for modelling plasma-surface interactions use particle-based simulations. It would be impractical to simulate the motion and interactions of every electron and ion in the plasma so approximate particle-incell 12 or treecode 13 methods are commonly used. However, these methods are subject to small amounts of random noise and can require lengthy runtimes. An alternative simulation method is provided by kinetic codes which discretize and solve the Vlasov equation in order to directly obtain the electron and ion distribution functions near a surface without particle-based approximations; 14,15 these codes are also impeded by large memory and runtime requirements due to the large number of dimensions of their position-velocity phase space. Lengthy runtimes are a particular issue for modelling surface erosion processes or plasma-liquid interactions where the surface can deform into rapidly changing shapes with complicated geometries. A robust method for modelling a plasma in the vicinity of such surfaces is the first step in a full simulation of a plasma-liquid interface or the erosion of nonplanar surface features. Fluid-based techniques provide an alternative to particles with the advantages of smooth solutions and quick runtimes at the expense of losing information about particle energy distributions and hence kinetic effects. However, they still capture many physical phenomena such as the breakdown of quasineutrality and formation of a sheath. 16,17 Further possible applications of fluidbased approaches to plasma-surface interactions include selfconsistent modelling of dust clouds and cathode arc spots. This paper introduces a numerical scheme for solving the coupled cold-ion fluid equations, the electron Boltzmann relation, and the Poisson's equation in two dimensions in cartesian or axisymmetric cylindrical coordinates.  PLASMAS 25, 043514 (2018) applications of these solutions to problems involving plasmaliquid interactions, plasma etching by ion sputtering, and electron emission from the surface are also discussed.

II. MODEL
The collisionless cold-ion, Boltzmann-electron model comprises the coupled equations 16 @n @t þ r Á ðnuÞ ¼ 0; (1) which are the ion continuity equation, the ion equation of motion, and Poisson's equation with Boltzmann electrons, respectively, where n and u are the ion density and velocity, / is the electrostatic potential, and n 0 is the density of the quasineutral plasma far away from the surface. It is useful to consider the applicability of these equations for describing real-world plasma-surface interactions before progressing to find their numerical solutions. The origins of the cold-ion, Boltzmann-electron model can be traced back to the plasma theory of Tonks and Langmuir 18 who developed a model based on Boltzmann electrons and ions with no thermal motion; these ions are generated at rest and acquire a velocity given by the conservation of energy as they fall through the electrostatic potential of the plasma. The different initial positions, and hence the initial potential energy, of these ions results in a nonthermal distribution of ion velocities. Bohm later assumed that all of the ions are generated with the same potential energy and therefore all have the same velocity at a given point. 19 The celebrated Bohm criterion, that ions must travel faster than their sound speed in order to break the quasineutrality of the plasma, emerged from this cold-ion theory. A similar model, of Boltzmann electrons and ions with a single-valued velocity, was employed by Allen, Boyd and Reynolds (ABR) to calculate the ion currents to cylindrical and spherical probes. 20 The ions were treated as a collection of individual particles in all of these early studies; it was not until a pair of papers published in 1965, by Woods 21 and Kino and Shaw, 22 that the treatment of the ions as a fluid first originated and Eqs. (1)-(3) appeared in the form given here which is amenable to cases in two or three dimensions.
Experimental measurements have vindicated the use of the cold-ion, Boltzmann-electron model for describing plasma-surface interactions in low pressure plasma discharges which are common in industrial processes. Ion currents collected by cylindrical Langmuir probes with radii a $ k D in low-pressure argon and neon RF plasmas are in excellent agreement with the cold-ion ABR theory 23,24 for temperature ratios at least as high as T e =T i ¼ 0:3. 24 A lowpressure helium plasma begins to depart from the ABR theory for these elevated temperature ratios due to the enhanced thermal velocity of the light helium ions which causes nonnegligible angular motion. 24 Furthermore, the validity of the cold-ion approximation in the limit T i ( T e has been verified by particle-in-cell (PIC) simulations of cylindrical Langmuir probes. 25 Most computational models of warm-ion plasmas near surfaces retain the Boltzmann electron assumption in socalled hybrid PIC codes such as those of Sheridan 26 and Patacchini et al. 27 for warm-ion plasma-surface interactions in planar and spherical geometries, respectively. The validity of the Boltzmann electron assumption requires that most of the electrons are repelled from the surface. In this case, most of the electrons in the vicinity of the surface are returned to the bulk plasma; the repeated repulsion of electrons from the surface allows them to thermalize to a Maxwellian distribution even if the physical size of the plasma discharge is smaller than the electron collision length. The Maxwellian distribution extends all the way to the surface if it is fully reflecting, 27 so the Boltzmann relation is valid for highly negative surfaces which absorb negligible numbers of electrons. The case of a partially absorbing planar surface was considered by Kolobov where /ð0Þ is the potential of the surface located at z ¼ 0.
for various values of /ð0Þ. The error in electron density is shown in Fig. 1 and clearly tends to zero as e/ð0Þ=k B T e ! À1. For a typical surface potential of e/ð0Þ=k B T e ¼ À3, the error never exceeds 2.5% and the Boltzmann approximation is a reasonable assumption. Note, additionally, that Kolobov and Arslanbekov's ion distribution function, Eq. (35) of Ref. 14 tends to a delta function in velocity space in the limit T i ! 0. The peak of the delta function is determined by the conservation of ion energy in exact accordance with the cold-ion model.

III. IMPLEMENTATION
The approximations made by the cold-ion, Boltzmannelectron model have been examined in detail; the implementation of a numerical method for their solution is considered next. The form of the divergence and Laplacian terms depends on whether the equations are being solved in cartesian (x, z) coordinates or axisymmetric cylindrical (r, z) coordinates. The time-dependence provides a useful tool to allow relaxation of the solutions to a steady-state but it also allows the simulation of dynamic processes such as those occurring in plasma-liquid interactions and near the periodically biased electrodes of an RF discharge. The examples considered in this paper focus on the steady-state solutions of the code which are obtained after numerous units of time have elapsed.
Equations (1)-(3) are simplified with dimensionless variables which are specified by upper-case symbols and are related to the dimensional variables by where c s ¼ ðk B T e =m i Þ 1=2 is the ion sound speed and k D ¼ ð 0 k B T e =n 0 e 2 Þ 1=2 is the Debye length. The resulting set of dimensionless equations is where r now indicates differentiation with respect to the normalized spatial coordinates.
Equations (8)-(10) are discretized using a finite difference method on a uniform staggered grid. This grid stores values of the ion density and electrostatic potential at the centre of cells, while ion velocities are stored at cell edges. This arrangement reduces numerical instabilities caused by grid-scale oscillations. The spatial derivatives are taken in the upwind direction, i.e., towards the source of the plasma flow, and time-stepping is performed with an explicit forward Euler method. Full details of the numerical scheme are given in Appendix A.
A layer of ghost cells is placed around the edge of the simulation domain and the variables are carefully chosen in these cells in order to maintain the boundary conditions of each problem. This paper only considers Bohm-like sheaths with two-dimensional features in cartesian or axisymmetric coordinates. The equality form of the Bohm criterion 19 supplies an asymptotic boundary condition, that U ! 0; N ! 1; U ! 0, and V ! À1, which applies as z ! 1. These values cannot be directly applied at the top of the simulation domain which is a finite distance from the surface; however, a series expansion for small values of U is given in Appendix B which does provide a suitable upper boundary condition. The residual error from the series expansion is shown to have a negligible effect on the sheath solutions in Sec. IV A.
The cartesian simulations are assumed to be periodic in the x-direction and to have odd symmetry as shown in Fig. 2. The simulation domain can therefore be taken as half of the period width with reflective boundary conditions applied at the left and right edges: u is zero and the x-derivatives of v, n, and / vanish at these boundaries. These boundary conditions must also be enforced on the axis of a cylindrical simulation, at r ¼ 0, but the outer boundary does not provide any symmetries to simplify the boundary conditions. However, if the surface is nearly planar and the outer edge is far away from any nonplanar features, as is the case in Sec. IV D, then the boundary conditions of Fig. 2 remain suitable. The numerical implementation of these boundary conditions is given in Appendix B.
The surface in contact with the plasma is taken to be perfectly conducting and hence is an equipotential surface. This surface acts as a sink of ions and electrons from the plasma so the ions always flow towards the surface. The upwind scheme for the spatial derivatives therefore means that the ion velocity and density do not need to be specified at the surface. Indeed, determining these values is often the motivation for simulations of plasma-surface interactions. Some of the examples in this paper are for electrically isolated surfaces; the floating potentials of these surfaces are found by equating the electron and ion fluxes to the surface to give where the fraction A s =A 0 is a modification to the equivalent expression for a planar surface. 28 This ratio is simply the fractional increase in area due to nonplanar surface features.
Simulations of plasma-surface interactions with a fixed voltage applied to the surface are performed with the potential of the wall set accordingly. The cold-ion, Boltzmann-electron model is implemented in the Cþþ code Cold Ions with Boltzmann Electron Relation (CIBER) which is freely available at http://www.github.com/ joshholgate/CIBER along with all of the data from the examples given in this paper. The simplicity of the numerical scheme results in a compact code which allows rapid modification and extension of the scheme.

IV. RESULTS
The cold-ion, Boltzmann-electron model, as implemented in the code CIBER, is used to simulate a number of different plasma-surface interactions. These simulations serve several purposes: to verify that the code is producing physically realistic and numerically accurate results, to show the versatility of the numerical scheme and, most importantly, to demonstrate the usefulness of the code in modelling different physical processes at the plasma-surface interface.

A. Planar sheath
The code is tested against the one-dimensional planar sheath solution before proceeding to two-dimensional examples which do not facilitate such simple comparisons with the test cases. Equations (8)-(10) can be solved implicitly, for the planar case, in terms of U as which follow from the conservation of ion energy, ion flux, and total stress through the sheath; the conservation equation for stress is derived in Appendix C. The explicit solutions, in terms of Z, are found by numerically integrating the dimensionless electric field F ¼ ÀdU=dZ with a shooting method starting at the wall in order to provide a test case for planar simulations using CIBER. This provides the classic Bohm sheath solution shown in Fig. 3. The floating potential of a planar surface in a hydrogen plasma, U f ¼ À2:839, is used in this example. The spatial convergence of solutions is tested by simulating a section of a planar sheath with height 20k D and width k D on uniformly spaced grids with varying grid resolutions. The resulting potential profiles along x ¼ 0, which are obtained after 1000 units of time have elapsed, have the known Bohm sheath solution subtracted from them in order to calculate the error in the electric potential. The results are displayed in Fig. 4 and it is evident that the error decreases as the grid resolution increases. The convergence of the solution to a steady state is examined in detail in the next example, specifically in Fig. 7, but it is expected to occur on the timescale taken by an ion to cross the simulation domain. Ions travelling at the Bohm speed will take 20 time units to traverse this 20k D -high domain so 1000 time units is more than sufficient to achieve a steady state.
Although the results show good grid convergence, all of the solutions in Fig. 4 suffer from the same error due to the boundary condition at the top of the simulation domain. This error is investigated using simulations of various heights but with the grid spacing fixed at 0:02k D . The error analysis is performed as before and the results are displayed in Fig. 5. In this case, it is clear that the error in the boundary condition decreases as the simulation height is increased as should be expected. The magnitude of the error at the upper boundary of the simulation domain is also less than the maximum error value, which appears around z % 3k D , when the simulation height is equal to or greater than 20k D . It is finally noted that the results are essentially identical in cartesian and axisymmetric coordinates which is unsurprising as the coordinate parallel to the planar surface is ignorable.

B. Sinusoidal cathode
The plasma sheath of a sinusoidally corrugated cathode has been considered recently in the context of a linear perturbation theory. 11 The motivation of that work was not only to find the linear stability criterion of a conducting liquid surface due to the electric field in the sheath but it also provides a second reference case to compare with the results of CIBER. The simulation is set up in the same way as for the planar sheath in cartesian coordinates but the wall now has a sinusoidal height perturbation with an amplitude of 0:05k D and wavelength of 2k D so that reflective boundary conditions can be applied at x ¼ 0 and x ¼ k D as illustrated by Fig. 2. The area ratio for this surface is A s =A 0 ¼ 1:0015 which gives a hydrogenic floating potential of U ¼ À2:837 which is hardly changed from the planar case. Each simulation uses a different grid resolution and is run for 1000 units of time. Figure 6 displays the profiles of U along a line at z ¼ 0:2k D at the end of each simulation together with the predictions of the linear perturbation theory. In this case, the convergence of solutions to the linear theory is clearly evident without the need for an explicit error calculation. The agreement between the simulations and the linear perturbation theory mutually verifies the results of both approaches.
The errors in the simulation can also be quantified by testing for ion energy conservation. The dimensionless ion energy is calculated as U 2 =2 þ V 2 =2 þ U which should yield a constant value of 0.5 according to the Bohm criterion. The magnitude of the error in energy is calculated in each cell and the maximum error value is displayed in Fig. 7 as a function of simulation time for different grid spacings. The maximum energy error decreases over the course of each simulation as they converge towards the steady solution. The optimal solution is obtained after around 400 units of time beyond which additional timesteps do not improve the energy conservation. The energy error also decreases as the grid resolution increases which provides further evidence for the spatial convergence of the scheme.  The example of a wall with a trench cut into it is considered next. This situation could arise, for example, in the etching of features onto silicon computer chips by ion sputtering. The ever-decreasing size of features on these microchips requires the etching of trenches with large depth-to-width, i.e., aspect, ratios. The study of aspect ratio dependent etching remains an active area of research. 29 The cold-ion, Boltzmann-electron model is employed here to investigate the dependence of plasma parameters on the depth of the trench. The trench is simulated as half of a rectangular wave so that the boundary conditions of Fig. 2 can be used. Industrial plasmas such as those used in semiconductor processing have an externally applied surface potential; therefore a fixed surface potential of -10 is used.
The simulation is initialized in a similar way to the previous two examples but with the surface given by a rectangular wave with a wavelength of 2k D , i.e., a trench width equal to the Debye length, and with a depth between 1 and 32 Debye lengths so the boundary conditions of Fig. 2 remain valid. The region above the top of the trench is kept at a constant width of 20k D and a regular grid spacing of 0:02k D is used. An example of the full two-dimensional electrostatic potential profile around a trench of depth 2k D , after 1000 units of time have passed, is shown in Fig. 8. The variation in the potential occurs mostly above the top of the trench and the potential is relatively constant within the trench.
The effect of trench depth is shown in more detail in Fig. 9 which gives z-profiles of ion density along the trench centres. The profiles are indistinguishable on the scale of this figure and differ only in the value of z at the base of the trench. The depth of the trench appears to have no effect on the sheath above the top of the trench and, even within the trench, the profiles appear identical but with a cutoff depending on where the sheath reaches the bottom of the trench. The density drops off quickly at the trench entrance but levels off further into the deep trenches. These observations illustrate that the reduction in ion density is a key contributing factor to the drop off in etching rates for medium aspect ratio trenches although the interaction of the plasma with neutral particles, which is not included in the present model, has also been shown to play a significant role. 29  The final example of plasma-surface interactions to be considered is of a nearly planar surface with a small hemiellipsoidal protrusion which demonstrates the application of CIBER to problems in cylindrical geometries. The physical motivation for these shapes is for problems involving electron emission from the surface. The emitted current density of electrons from the surface is highly sensitive to the local surface field strength which is itself strongly enhanced by the high curvature of small irregular features on the surface. 30,31 The large emitted electron current from these small irregular features is an important trigger for the formation of unipolar arc spots which are believed to provide a significant source of droplets which contaminate and degrade the plasma of magnetic confinement fusion devices. 32,33 The increase in electric field strength at a surface protrusion can be quantified with field-enhancement factors which give the ratio of the surface field strength at the curved feature to the field strength for a perfectly flat surface under the same conditions. These field-enhancement factors have recently been measured for various surface structures formed by plasma exposure, such as bubbles and fuzz, to be greater than 1000 under vacuum conditions. 34 However, the charge in a plasma sheath can mitigate these high field-enhancement factors. 35,36 The fieldenhancement factors of hemi-ellipsoidal protrusions are calculated in this subsection; these shapes are well-characterized for emission in a vacuum 30 and provide good representations of the elongated fuzz fibres and flattened bubbles which form on plasma-exposed surfaces. 34 Simulations are performed in axisymmetric coordinates by placing semi-ellipsoids of half-width 0:2k D and heights in the range ð0-1Þk D at the centre of the surface. A simulation domain with height 20k D and width k D is used with a grid spacing of 0:01k D . The potential of the wall is kept constant as the floating potential for a planar surface in contact with a fusion-relevant deuterium plasma, U f ¼ À3:185 so that the variation in the electric field is solely due to geometric effects. The resulting electric field strength close to a planar surface with a hemi-ellipsoidal protrusion of aspect, i.e., height-to-width, ratio 1 is displayed in Fig. 10. This ellipsoid is, of course, simply a sphere; aspect ratios greater than one give elongated, needle-like protrusions while aspect ratios less than one give flattened surface features. The field strength at the apex of the sphere is greater than at a planar surface but its enhancement factor, which is 1.91, is significantly less than the vacuum value of 3. 30 The field-enhancement factors for hemi-ellipsoids with aspect ratios in the range 0-5 are illustrated in Fig. 11 and compared with the vacuum theory. The field-enhancement factors are strongly suppressed when the plasma is present. This result is encouraging for the use of classical electron emission formulae when studying plasma-surface interactions because it indicates that surface roughness, which is often poorly characterized, does not play as large a role as in vacuum emission. The plasma will therefore suppress the emitted electron current from the surface in accordance with previous studies of field emission using an electron-only fluid model 35 and thermal emission using PIC simulations. 36

V. SUMMARY
A numerical implementation of the two-dimensional cold-ion, Boltzmann-electron model of a plasma in the vicinity of a nonplanar conducting surface has been presented. This implementation is highly versatile as evidenced by the numerous examples of nonplanar sheaths which have been considered in this paper. These examples were principally chosen in order to illustrate the adaptability and numerical accuracy of the model, but they have also provided physical insights into various aspects of plasma-surface interactions. As such, it is anticipated that the cold-ion, Boltzmann FIG. 11. The electric field strength at the apex of the hemi-ellipsoidal protrusion, as a function of the aspect ratio of the ellipsoid, together with the ratio of this value to the field strength at a planar surface i.e., the fieldenhancement factor. The field-enhancement factors are strongly suppressed by the presence of the plasma. electron model, and its two-dimensional implementation in the code CIBER, will provide a useful tool for future studies of the physics of plasma-surface interactions in topics involving, for instance, plasma-liquid interfaces and surface erosion by the plasma. The simplicity of the code should also allow its extension to include additional physical processes such as electron emission from the surface, collisions and ionization in the plasma, the effect of magnetic fields on the ions, and coupling of the code to moving surfaces. The outlook for fluid-based models of plasma-surface interactions holds significant promise.

ACKNOWLEDGMENTS
This work was supported by the UK's Engineering and Physical Sciences Research Council and the Imperial College President's Ph.D. Scholarship Scheme. Useful discussions with Professor John Allen are gratefully acknowledged. The code and data underlying this work are available at http://www.github.com/joshholgate/CIBER.

APPENDIX A: DISCRETIZATION OF EQUATIONS
This appendix contains details of the finite difference scheme for discretizing the cold-ion, Boltzmann-electron equation on the staggered grid shown in Fig. 12. The ion continuity equation, Eq. (8), is expanded in cartesian coordinates as This is discretized using where H is the uniform spacing of the grid and the fluxes at the cell-edges are given by upwind differences as C iþ1=2;j ¼ maxðU iþ1=2;j ; 0ÞN i;j þ minðU iþ1=2;j ; 0ÞN iþ1;j ; C iÀ1=2;j ¼ maxðU iÀ1=2;j ; 0ÞN iÀ1;j þ minðU iÀ1=2;j ; 0ÞN i;j ; C i;jþ1=2 ¼ maxðV i;jþ1=2 ; 0ÞN i;j þ minðV i;jþ1=2 ; 0ÞN i;jþ1 ; which ensures that the value of the ion density used in the calculation is taken from the direction closest to the source. In axisymmetric cylindrical coordinates, Eq. (8) is and the discretization follows as: where R i ¼ ði À 1=2ÞH, so that R ¼ 0 is located at the boundary between the i ¼ 0 and i ¼ 1 cells, and the upwind fluxes are the same as in the cartesian case. The time derivative is discretized using the explicit forward Euler method, using k as the time step index, as where the derivative on the right-hand-side is taken from either Eq. (A2) or Eq. (A4). The ion momentum equation, Eq. (9), is also discretized using an upwind method with explicit Euler time steps. The X-component of this equation is discretized as Some of the derivatives cannot be taken to the second order near the boundaries where i ¼ 3/2, i ¼ i max À 3=2, j ¼ 1 or j ¼ j max À 1 and, where the upwind scheme demands it, they are replaced by first order derivatives pointing towards the edge of the simulation domain. Boundary values are fixed in cells The derivatives in the ion momentum equation take exactly the same form in cartesian and axisymmetric cylindrical coordinates so the cylindrical version of this equation is obtained by simply replacing X with R. The Z-component of the ion momentum equation follows from a discretization of using the same method as for the X-component.
The evaluation of the new velocity values is clearly dependent on the electrostatic potential. The potential is found from the Poisson-Boltzmann equation, Eq. (10), which is evaluated using overrelaxation according to where x is the relaxation parameter and the index l is incremented until all U i;j values have been found to the desired degree of accuracy. The value of x is set between 1 and 2 to determine the strength of the overrelaxation. Alternatives to overrelaxation, such as fast multigrid methods, are often used in modern codes but these are difficult to implement with arbitrary curved boundaries. Furthermore, the speed bonus from such methods is mitigated because the convergence rate of relaxation methods depends on the quality of the initial estimate of the U i;j values and, by using the solution from the previous timestep, very good initial estimates are readily available. It therefore takes only a couple of relaxation steps for the solutions to reach the desired level of accuracy.
The preceding numerical scheme is robust and stable provided that the Courant-Friedrichs-Lewy (CFL) condition is satisfied. This can be written as where U max is the highest ion velocity at any position in the simulation. The value of dT is chosen adaptively at each timestep to be one quarter of the value given by the marginal CFL condition.

APPENDIX B: SHEATH EDGE EXPANSION AND BOUNDARY-CELL VALUES
The grid values in the boundary layer must be chosen carefully in order to ensure that the necessary boundary conditions are imposed on the simulation as described in Sec. III and illustrated in Fig. 2. This Appendix provides details of the numerical implementation of these conditions in the layer of ghost cells surrounding the main simulation domain. However, the boundary conditions given by the Bohm criterion are applicable in the asymptotic limit z ! 1 so a simple expansion of the planar sheath solution, given by Eqs. (12)- (14), is first developed at U ¼ 0 which allows the application of the Bohm criterion at a large ð) k D Þ but finite distance from the surface.
The series expansion of Eq. (14) at which may be integrated and rearranged to give where C is the constant of integration. The Bohm condition gives U ¼ 0 when Z ! 1 so the constant of integration is zero. The approximate potential at the sheath edge follows as: which is valid provided that U ( 1. This approximate solution is plotted against the exact planar sheath solution in Fig. 13 and is seen to be an excellent approximation far FIG. 13. The variation of electric potential in a planar sheath as given by numerical integration of Eq. (14) and by the small-U approximation, Eq. (B3), together with the difference between the two expressions. Note that the values of z begin at 15k D rather than at the wall. The approximate solution becomes an appropriate boundary condition for small values of U. away from the plasma-facing surface. All of the simulations in this paper use this approximation as their upper boundary condition so any nonplanar perturbations are assumed to be negligible at this surface. The series expansion of Eqs. (12) and (13) gives simple approximations for the ion velocity and density as which are again used as upper boundary conditions for the simulation. The boundary condition at the top of the simulation domain is therefore given by U iþ1=2;j max þ1 ¼ 0 with U i;j max þ1 ; V i;j max þ1=2 , and N i;j max þ1 given by Eqs. (B3)-(B5), respectively, for i in the range ½1; i max .
The boundary conditions at the left and right edges, as illustrated in Fig. 2, require the u component of velocity and the x-derivatives of v, n, and / to vanish. These conditions are enforced by the boundary values U 1=2;j ¼ U i max þ1=2;j ¼ 0; V 0;jþ1=2 ¼ V 1;jþ1=2 ; V i max þ1;jþ1=2 ¼ V i max ;jþ1=2 ; N 0;j ¼ N 1;j ; N i max þ1;j ¼ N i max ;j ; U 0;j ¼ U 1;j , and U i max þ1;j ¼ U i max ;j for j in the range ½1; j max . The lower-edge boundary conditions require only a constant value of U at the interface which forms the equipotential plasma-facing surface; the examples in this paper all use the floating surface potential given by Eq. (11).

APPENDIX C: PLASMA STRESS TENSOR
The plasma stress tensor is a conserved, i.e., divergencefree, object in a steadily flowing, collisionless plasma with no sources or sinks of particles. The divergence-free property of the stress tensor can be derived from Eqs. (1)-(3) as follows. The electron Boltzmann relation is first differentiated to give rn e ¼ en e k B T e r/ ¼ À en e k B T e E; which is multiplied by k B T e and added to the ion equation of motion to give m i n i ðu Á rÞu þ k B T e rn e ¼ eðn i À n e ÞE: The ion continuity equation, r Á ðn u uÞ ¼ 0, allows the first term to be written as a single divergence while Gauss's law allows the right hand side to be rewritten to produce r Á ðm i n i u uÞ þ k B T e rn e ¼ 0 ðr Á EÞE: The electron density term can be expressed as a divergence with the use of the identity matrix. The time-invariant Maxwell equation r Â E ¼ 0 then allows the whole expression to be formulated as r Á ðm i n i u u þ k B T e n e IÞ ¼ 0 ðr Á EÞE À 0 E Â ðr Â EÞ: Finally, the vector identity together with the electron Boltzmann relation, yields the result where the stress tensor is I À m i nu u À n s k B T e exp e/ k B T e I: (C7) The individual elements of this tensor are identified as the electrostatic Maxwell stress, the ion ram pressure, and the isotropic electron thermal pressure. Equation (C6) can be integrated in the planar case in order to produce the algebraic expression for the electric field in terms of the potential in Eq. (14). The divergence-free property of the stress tensor in a steadily flowing plasma indicates that, within a closed volume, all of the electrostatic forces are balanced by the electron pressure and change in momentum of the ions. However, the introduction of an absorbing object or surface into the plasma causes an imbalance in the divergence of the stress tensor which, if integrated over a closed volume, is equal to the net force on the plasma within that volume or, equivalently, the negative of the force on the object. The stress tensor therefore provides a straightforward method for calculating the net force exerted by a plasma on an object. It is also extremely useful in determining the stress balance at a plasma-liquid interface by shrinking the volume under consideration to an infinitesimal region enclosing the interface and, after including the surface tension, viscous stress, and pressure of the liquid, applying the divergence theorem to give the interface jump conditions. 11