Density Estimation Techniques for Multiscale Coupling of Kinetic Models of the Plasma Material Interface
DDensity Estimation Techniques for Multiscale Coupling of KineticModels of the Plasma Material Interface
Shane Keniley, Davide Curreli ∗ April 7, 2018
Abstract
In this work we analyze two classes of Density-Estimation techniques which can be used to consis-tently couple different kinetic models of the plasma-material interface, intended as the region of plasmaimmediately interacting with the first surface layers of a material wall. In particular, we handle thegeneral problem of interfacing a continuum multi-species Vlasov-Poisson-BGK plasma model to discretesurface erosion models. The continuum model solves for the energy-angle distributions of the particlesstriking the surface, which are then driving the surface response. A modification to the classical Binary-Collision Approximation (BCA) method is here utilized as a prototype discrete model of the surface, toprovide boundary conditions and impurity distributions representative of the material behavior duringplasma irradiation. The numerical tests revealed the superior convergence properties of Kernel DensityEstimation methods over Gaussian Mixture Models, with Epanechnikov-KDEs being up to two orders ofmagnitude faster than Gaussian-KDEs. The methodology here presented allows a self-consistent treat-ment of the plasma-material interface in magnetic fusion devices, including both the near-surface plasma(plasma sheath and presheath) in magnetized conditions, and surface effects such as sputtering, back-scattering, and ion implantation. The same coupling techniques can also be utilized for other discretematerial models such as Molecular Dynamics.
The near-wall region of a magnetized plasma is a far-from-equilibrium system characterized by a multitudeof coupled physical phenomena. Particles sputtered from the surface will interact with the plasma and eitherenter the plasma bulk or be redeposited, where they can both induce further sputtering events and change theconstitution of the surface layers. A computational model of Plasma-Material Interactions (PMI) must be ∗ SK: [email protected], DC: [email protected],
University of Illinois at Urbana Champaign. a r X i v : . [ phy s i c s . c o m p - ph ] A p r ble to incorporate a kinetic description of the plasma sheath and presheath, surface response, and impuritytransport in order to capture these processes. However, modeling PMI is complicated by the disparate time-and length-scales involved: relevant plasma processes occur over millimeters of length and microseconds oftime, while surface interactions take place on the order of an angstrom and picosecond timescales. Accuratelysimulating PMI is thus not only a matter of implementing the necessary physics, but also developing thetechniques for integrating physical models which necessarily operate at different scales.Multiscale Modeling approaches this issue by using a separate model to describe each physical process,and later coupling the processes together using an opportune numerical strategy. A summary of multiscaleand multiphysics research was compiled by Groen et al. [1]. This paper is aimed at investigating a couplingmethodology based on Density-Estimation (DE) techniques, and apply such techniques to the problem ofconsistently passing continuum plasma distributions to a discrete material model and back. We consider thegeneral problem of dealing with fully-kinetic distributions for both the plasma and the material species, thatis typically referred as the “kinetic-to-kinetic” coupling problem, as opposed to the “kinetic-to-fluid” couplingproblem. For the plasma species we adopted a continuum multi-species electrostatic Vlasov-Poisson solverwith a Bhatnagar-Gross-Krook (BGK) collision operator; such a model provides a fully-kinetic description ofall relevant species in the partially-ionized conditions normally encountered close to a material surface. Forthe material we used a discrete surface model based on the Binary Collision Approximation (BCA) includingdynamic surface composition and surface morphology, such as those implemented in the Fractal-TRIDYNcode [2]. The continuum solver accurately describes the plasma sheath and presheath in front of the materialsurface, including full-orbit effects due to finite ion Larmor radius. The same model also provides the IonEnergy-Angle Distributions (IEAD) of the plasma particles striking the wall, which are discretized and usedas an input to the surface model. The surface model in turn provides discrete energy-angle distributionsof sputtered impurities, together with the distributions of reflected and implanted particles. The statisticalvalidity of the coupling has been analyzed to ensure that the passage from discrete to continuum is performedby guaranteeing continuity of the conservation laws at the interface between the two models.For near-surface plasma problems involving plasma sheaths, a continuum discretization of the Vlasov-Poisson-BGK problem is generally preferable rather than discrete approaches such as Particle-in-Cells (PIC).PIC methods utilize a statistical sampling of the phase-space by tracking a finite number of macroparticles, N p . PIC methods have been previously adopted for near-wall simulations [3, 4, 5], but they tend to under-sample the spatial regions closest to the wall where the plasma density drops of more than an order ofmagnitude with respect to the upstream density. Furthermore, PIC codes are inherently subject to numerical2oise, which can obscure features of the near-wall plasma, such as the steep gradients occurring across themagnetic presheath and inside the Debye sheath. PIC noise can be controlled in a number of ways, suchas increasing the number of computational particles N p , but it typically diminishes only by a factor of (cid:112) N p . In contrast, continuum kinetic solvers directly discretize the Vlasov-Poisson-BGK problem in phase-space with a mesh. They are more computationally expensive than equivalent PIC models, but avoidthe numerical noise problem of PIC codes. The method allows higher resolution even in regions of largegradients such as the plasma sheath. Eulerian schemes applied to the problem have been studied withsemi-Lagrangian methods [6, 7, 8, 9], discontinuous Galerkin schemes [10], and finite difference methods[11]. In this work we will adopt a high-order Cheng-Knorr Strang splitting of the Vlasov operator. Theerosion process itself is traditionally modeled using either Molecular Dynamics (MD) or a BCA approach.The advantage of BCA over MD is mainly computational, since the ion-material interaction is reduced froma multi-body problem to a Monte Carlo sequence of binary collisions, allowing faster computations andstraightforward parallelization. Furthermore, the sputtering yields computed via BCA have been shown tohave good agreement with experimental measurements even at low-energy in regimes close to the sputteringthreshold [12]. For a review on atomic collisions in amorphous targets with BCA methods, see [13].Thanks to the coupling methodology developed in this paper, a self-consistent kinetic model of the near-surface plasma (plasma sheath and presheath) including material response has been obtained. The modelallows kinetic analyses of the near-surface plasma behavior, and numerical predictions of the erosion, re-flection, and implantation rates during PMI activity. In the following sections, first we present the modelequations of both the plasma and the material in order to establish the necessary exchange of informa-tion between the two regions, and then we describe the method of dynamically coupling the multi-speciesBoltzmann-Poisson plasma treated with finite differences to the BCA model describing the plasma-facingsurface (Section 2). Tests to ensure the statistical validity of the coupling method are presented in section3, including an example of the dynamically-coupled model of a plasma sheath. Consider a volume dV = d x d v of the phase space ( x , v ), where x = ( x , x , x ) is the position vector and v = ( v , v , v ) is the velocity vector, populated with N s plasma particles of species s in which the densityof particles is sufficiently large that in an infinitesimally small region the condition N s d x d v (cid:29) t may be described by thedistribution function f s ( x , v , t ). The time evolution of such a distribution function over the phase space( x , v ) is governed by the Boltzmann-Poisson integro-differential problem [14]: ∂f s ∂t + v · ∇ x f s + q s m s ( −∇ x φ + v × B ) · ∇ v f s = (cid:18) ∂f s ∂t (cid:19) coll (1) (cid:15) ∇ x φ = q e n exp (cid:18) q e φKT e (cid:19) − (cid:88) s q s (cid:90) f s d v (2)where m s and q s are the mass and the electric charge of all massive species s (ions and neutrals), φ the electricpotential, and B the magnetic field. The collision operator at the right-hand side of Eq. 1 accounts for theFokker-Planck and atomic modifications to the distribution function f s due to collisional processes. In firstapproximation the collision operator within the plasma region can be treated as a Bhatnagar-Gross-Krook(BGK) operator [15], (cid:18) ∂f s ∂t (cid:19) coll = ν bgk ( f s, − f s ( x , v , t )) (3)allowing explicit handling of elastic scattering, electron-neutral collisions, and charge-exchange between ionsand neutrals, where each process is described by its relaxation frequency ν bgk . Collisions in the materialregion are treated with a dedicated binary-collision operator accounting for ion-matter interaction, describedin Section 2.2.The Vlasov-Poisson-BGK problem of Eqs. (1)-(2)-(3) is numerically solved in one spatial dimensionand three velocity dimensions (1D3V) by means of a finite-volume discretization of the phase-space with aCheng-Knorr [8, 16] high-order Strang splitting of the Vlasov operator [17, 18]:1. ∂ t f s + v x ∂ x f s = 0 x-advection (∆ t/ ∇ φ = − ρ/(cid:15) Poisson solver3. ∂ t f s + q s m s ( E x + v y B z − v z B y ) ∂ v x f s = 0 v x -advection (∆ t )4. ∂ t f s + q s m s ( v z B x − v x B z ) ∂ v y f s = 0 v y -advection (∆ t )5. ∂ t f s + q s m s ( v x B y − v y B x ) ∂ v z f s = 0 v z -advection (∆ t )6. ∂ t f s + v x ∂ x f s = 0 x-advection (∆ t/ ∂ t f s = ν bgk ( f − f ( x, v )) Collision operator4perator splitting reduces the Vlasov-Poisson-BGK problem of Eqs. (1)-(2)-(3) from a multidimensionalintegro-differential problem to a sequence of problems of lower dimensionality, each one with uniquely-definednature (elliptic, hyperbolic, quadrature). High-order upwind schemes are used for the hyperbolic sections ofthe problem, finite-difference discretization for the elliptic portion, and high-order Simpson schemes for thequadratures.A key feature desired for plasma sheath simulations is that the plasma always flows from the bulk towardthe boundary. Upwind-biased numerical schemes exploit this directionality, achieving the highest accuracyper gridpoint in a finite difference stencil [19]. First order upwind-biased methods are stable and convergent,but are too dissipative without a highly-refined grid. This work will instead focus on high-order upwindmethods. In particular, since the rate at which error is decreased with increasing resolution is faster by anorder of magnitude for a high-order method, the fourth-order upwind scheme will be used for our simulations.The third-order method has a less strict stability condition which would allow for larger timesteps to be used[20], but the error control of the higher order scheme was considered more beneficial than including a largertimestep. High order time approximations are also necessary for the time-pushing in order to achieve thesame accuracy of the grid upwinding, the most common example of which are the set of Runge-Kuttaalgorithms [21]. The present work utilizes a fourth-order explicit Runge-Kutta (RK4) method.The simulation of a near-wall plasma requires a grid size small enough to resolve the plasma sheath, whichrequires a resolution on the order of the Debye length (∆ x = O ( λ De ) ). However, since using such a highresolution across the entire domain is computationally expensive, a nonuniform grid that is highly refinedin only a small region is desirable. In this work we adopt a hyperbolic-tangent mesh refinement suitable tothe plasma sheath problem, providing fine resolution on one end of the domain for the plasma sheath andcoarse resolution on the other for the plasma bulk, with a gradual transition between the two regions [22]: g ( s ) = ∆ x + ∆ x − ∆ x α ( s − s x is the refined grid step near the wall, ∆ x is the grid step in the plasma bulk, α is the steepnessof the transition region, and s is the location of the transition. The adoption of a non-uniform spatial grid,Eq. 4 requires a modification of the finite difference schemes used for the discretization of the operators. Theapproximate (discrete) solution u i = u i ( x ) of a generic function of the spatial coordinate ˆ u = ˆ u ( x ) may bedefined on a general numerical grid with nonuniform grid spacing with an N g -order polynomial, where N g is the order of accuracy of the resulting numerical scheme. As an example, a second-order upwind-biased5iscretization at a point i is represented by the following set of polynomials:ˆ u = ax + bx + c (5) u i = c (6) u i − = a ( − h i − ) + b ( − h i − ) + c (7) u i − = a ( − h i − ) + b ( − h i − ) + c (8)with undefined coefficients a , b, and c . The polynomial fit defines a linear system of equations. An expressionfor the coefficients a , b , and c may be found by inverting the resulting matrix, and the numerical discretizationat any point is then derived by taking the derivative of ˆ u : d ˆ udx = 2 ax + bd ˆ udx = 2 a Thus the second-order upwind-biased discretization of the first derivative is defined as the coefficient b .For a uniform grid spacing, it is simple to verify that the parametric fit method reproduces the traditionalsecond-order upwind scheme. The same treatment is used to derive nonuniform spacing for the third- andfourth-order upwind schemes, along with the second- and fourth-order central difference schemes for Eq. 2. The trajectory of each individual ion leaving the plasma and entering into a material surface of particle density N follows a random walk characterized by: (1) large deflections of the incoming particle at each encounterwith the nuclei, and (2) continuous energy losses between one nuclear encounter and the following due toinelastic electronic interactions. Here we consider the problem of evaluating the change to the distributions f s needed for surface erosion and ion implantation considering the following Fokker-Planck problem: ∇ · [ V f s − ∇ ( D f s )] = S (9)The problem of Eq. 9 describes the drift and diffusion of the plasma and material particles traveling throughthe material under a source S determined by the plasma sheath. Eq. 9 can be recast [23] as a StochasticDifferential Equation (SDE), and resolved with a Lagrangian approach by reconstructing the distributions6rom a large number of uncorrelated particle histories experiencing a Wiener process. It can be shown thatthe Fokker-Planck problem of Eq. 9 is equivalent to the following Ito SDE: d x ( t ) = V ( x , t ) dt + D ( x , t ) dW ( t ) (10)where x is the the particle position vector, V ( x , t ) the drift vector, D ( x , t ) the diffusion tensor, and dW ( t )is the differential element of a multi-dimensional Wiener process. When Eq. 10 is projected on a referenceframe (ˆ e , ˆ e , ˆ e ) having the third axis ˆ e parallel to the surface normal, ˆ e // ˆ n , the following Ito equationis obtained (written by components): d x ( t ) = v i ˆ e i dt + (cid:112) D i ˆ e i dW i , i = 1 , .., V = ( v , v , v ); and the three orthonormal components of the diffusion tensor are D i = eig ( DD T ) =( D , D , D ). A discrete first-order approximation to the Ito problem of Eq. 11 is then readily obtained fromthe Euler-Maruyama scheme through the Markov chain x k = x ( t k ) at discrete time intervals t k = k ∆ t : x k +1 = x k + v ki ˆ e k ∆ t + (cid:112) D i ˆ e i ∆ W i , i = 1 , .., W i ( i = 1 , .., W i = W i ( t k + ∆ t ) − W i ( t k ) (13)are either uniformly distributed random variables on a circle (for angular coordinates), or normally dis-tributed random variables (for translations) with zero expected value and variance ∆ t .In Binary Collision Approximation codes such as TRIM [24], TRIDYN [25], SDtrimSP [26], and Fractal-TRIDYN [2], the following simplifying assumptions are made to Eq. 12, v ki ∆ t = λ − s ( i = 1 , , , (14) (cid:15) k +1 = (cid:15) k − ∆ (cid:15) kloss , (15)In Eq. 14 the total path covered by the particle between two subsequent nuclear encounters is given by thedifference of the local mean free path λ = N − / and the distance of asymptotic deflection, s = R sin( θ/ R is the distance of closest approach and θ is the deflection angle of the ion trajectory expressed inthe Center-of-Mass (CoM) reference frame. The distance of closest approach R depends on the details ofthe ion-nucleus interaction potential V , which can be written as the product of a Coulomb-like potential ofconstant V o = q e Z Z / (4 π(cid:15) ), and a nuclear screening function φ ( ξ ) expressed as a series of Yukawa-likeexponentials (typically 3-4 terms), V = V o R φ ( ξ ) , φ ( ξ ) = (cid:88) i c i e d i ξ , where ξ = R/a. (16)The numerical value of R is then obtained by means of a Newton-Raphson scheme from the solution of thefollowing non-linear scalar equation (obtained from conservation of energy and angular momentum), (cid:15) r p + aRφ ( ξ ) − aR(cid:15) r ξ = 0 , (17)where p = r o / ( λπN ) is the collisional radius (or impact parameter), r o ∈ (0 ,
1] a uniformly-distributedrandom number, (cid:15) r = ( a/V o ) (cid:15) CoM the reduced energy, (cid:15)
CoM = (cid:15)m / ( m + m ) the CoM energy, a = a (9 π / / ( Z / + Z / ) − / the atomic screening distance, and a the Bohr radius, a = ¯ h/m e cα . Thedeflection angle θ is calculated in the CoM reference frame by solving the geometry of the binary collisionevent between the ion and the nucleus; the value of θ in the laboratory frame is then given by:tan Ψ = sin θm /m + cos θ (18)The diffusion tensor D = D (ˆ u k | Ψ , Φ) is obtained from the collision operator of binary collisions havinginteraction potential as in Eq. 16, where ˆ u k = (cos α, cos β, cos γ ) is the vector of the direction angles at thecurrent iteration k , the angle Ψ is the the deflection angle in the laboratory frame (Eq. 18), and Φ = 2 πr o (azimuthal deflection angle) is a uniformly distributed random variable on a circle. Additional diffusionprocesses, such as temperature effects, are typically neglected in the BCA approach. Including them in thescheme of Eq. 12 is a straightforward extension of the current treatment, and will not be considered in thepresent work.The integral energy loss (cid:15) loss between one collision and the following (Eq. 15) is given by the sum of theenergy transferred to the nucleus (cid:15) n and the inelastic electronic losses (both local (cid:15) loc and non-local (cid:15) nl ),∆ (cid:15) loss = (cid:15) n + (cid:15) loc + (cid:15) nl (19)8here the nuclear energy transfer is provided by a simple balance of the energy-momentum of the twocolliding particles, (cid:15) n = 4 (cid:15)m m / ( m + m ) sin ( θ/ (cid:15) nl = ( λ − s ) N S el and (cid:15) loc = d exp( d ξ ) S el / (2 πa ) are both functions of Lindhard’s [27] semi-classical electronic stopping power S el ∝ √ (cid:15) .The distributions of interest f s can be reconstructed from a large number of Markov chains ( N mc ),typically N mc ∼ O (10 − ), such as those in Eq. 12 with discrete displacements described by Eqs. 14–15. It is relevant to highlight that the individual Markov chains retain a purely Lagrangian character. Acomputational mesh is only required to collect the cumulative contribution of each individual trajectory andreconstruct the distributions f s . This is traditionally accomplished by means of a track length estimator.For a mesh made of c = 1 , ..., N c cells, an estimate of the distribution f s within each cell is obtained fromthe cumulative of the time τ jc spent by a particle j visiting the cell c : f s,c = 1 V c N mc (cid:88) j =1 w j τ jc , ( c = 1 , ..., N c ) (20)where V c is the volume of the c -esim cell and w j the weight of the particle. The method just described in thissection offers a linearized solution to the Ito problem of Eq. 11, and enables the treatment of the followingplasma-surface interaction processes (schematically depicted in Figure 1 for the case of a helium plasmaimpacting on an inhomogeneous tungsten wall): (1) material sputtering, (2) self-sputtering, (3) sputteringof implanted gas, (4) reflection, (5) implantation, and (6) prompt reemission. (a) Atomic Sputtering (b) Self Sputtering (c) Gas Sputtering(d) Reflection (e) Implantation (f) Prompt Reemission Figure 1: Illustrations of the plasma-surface interaction processes modeled by the method described inSec. 2.2. 9 .3 Plasma-Wall Coupling
In this section we describe a general method to couple the continuum plasma model described in Sec. 2.1to the discrete surface model reported in Sec. 2.2. The possibility to dynamically couple the two modelsand simulate the plasma-wall system “on-the-fly” strictly depends on the efficiency of the numerical schemeadopted for coupling. In this section we discuss the model equations used to merge the two physical models,with particular care to numerical efficiency. The coupling methodology requires: (1) a geometrical definitionof the surface shape; (2) a consistent strategy to convert particle distributions to continuum distributions;(3) consistency criteria to be satisfied in order to enforce the conservation laws at the interface.In the most general case, the shape of a material surface is a three-dimensional surface of topologicaldimension 2 ≤ D <
3, where D = 2 corresponds to an ideal flat surface, and intermediate dimensionsbetween two and three are for rough surfaces of fractal dimension D [ cite ]. The geometrical location of asurface is defined by an implicit function of space G ( x ) = 0. The simplest case of a flat surface correspondsto the function x − L = 0, with L the location of the wall and x a scalar coordinate perpendicular to thewall. More complex surfaces (e.g. fractals) can be defined using the same strategy as well [2]. The function G also provides a logical condition of surface crossing, useful to determine whether a particle or a plasmafluid element is inside or outside the wall, G ( x ) > G ( x ) < N mc individual Markovchains of Eq. 12, G ( x k ) >
0, in order to determine whether a particle traveling inside the wall has crossedthe surface and is thus released into the plasma. This logical test provides a number j = 1 , .., N ex of particlesreleased by the surface, where in general N ex ≤ N mc . The state vector of the particles leaving the surface attime t k is q j ( t k ) = (cid:0) x j ( t k ) , v j ( t k ) (cid:1) , where the particle velocity at time t k is given by v j = ˆ u j (2 (cid:15) j /M j ) / .The fractions of sputtered particles Y s , backscattered particles B s and implanted particles I s are givenrespectively by (per each species s ): Y s = (cid:88) p N ex,p N mc , B s = N ex,s N mc , I s = N mc − N ex,s N mc (23)where p is an index running over all secondary particles produced by the species s . The particle ensemble10 j ( j = 1 , .., N ex ) is converted to a continuum distribution by means of the Density Estimation operator K = K ( q j ), which provides an estimate ˆ f of the continuum distribution function f using information fromdiscrete particle vectors. For particles leaving the surface at time t k at those locations x specified by G ( x ) = 0, the reconstructed velocity distribution at the interface is provided by an operator K = K ( v j )which is a function of the particle velocity only,ˆ f s ( x , v , t k ) = 1 N ex ∆ v N ex (cid:88) j K (cid:18) v − v j ∆ v (cid:19) , at G ( x ) = 0 (24)where ∆ v is the window width (bandwidth) along the velocity coordinates, and v is a continuum velocitycoordinate. The choice of the density estimators K and of the optimal bandwidth ∆ v will be discussed indetail in Sec. 2.5. Note that the reconstructed distribution ˆ f s is by definition normalized to one, (cid:82) ˆ f s d v = 1,so that in order to restore physical units it must be multiplied by the density of particles backscattered bythe surface G ( x ) = 0, ∆ f s ( x , v , t k ) = ∆ n s ˆ f s ( x , v , t k ) , for G ( x ) = 0 , v · ˆ n > n is the surface normal, and the density ∆ n s is given by∆ n s ( x ) = (cid:88) s (cid:48) B s (cid:48) (cid:90) ∆ f s (cid:48) ( x , v , t k ) d v , for G ( x ) = 0 , v · ˆ n < s (cid:48) is an index that runs over all species contributing to backscattering of species s (sputtering plusreflection), including self-interactions of species s with itself. Here the quantity ∆ f s (cid:48) denotes the fraction ofparticles entering the wall between time t k and time t k +1 = t k + ∆ t ,∆ f s (cid:48) = (cid:90) t k +1 t k f s (cid:48) ( x , v , t ) dt, for G ( x ) = 0 , v · ˆ n < v · ˆ n , which is positive for particles leaving the wall, and negativefor particles entering the wall on the surface G ( x ) = 0, obeying the convention v · ˆ n > v · ˆ n < .4 Continuum-to-Particle: Rejection Sampling The conversion from a continuous distribution f s ( x , v , t ) to a particle distribution q j ( j = 1 , .., N mc ) wasobtained through a simple rejection sampling technique [28]. Tests of best utilization of different instrumentdistributions (uniform, normal) suggested that a multivariate normal distribution can achieve maximumutilization in the range of 60-70% by means of a covariance scaling factor comprised between 1 and 2. We have considered two classes of Density Estimation techniques for the calculation of the distribution ˆ f fromthe particle ensemble q j : (1) the parametric Gaussian Mixture Model (GMM), and (2) the nonparametricKernel Density Estimation (KDE). The two classes are representative of the two modern approaches todensity estimation, which may be split into: parametric , which assumes the Probability Density Function(PDF) may be fitted by some parametric family of functions, and non-parametric , which can be applied toan arbitrary data set without having any prior knowledge of the shape or the behavior of the underlyingdistribution function. Here we briefly describe the two techniques and the selection criteria for optimaldensity estimate.The Gaussian Mixture Model is parametrized by the component weights w g , the component means µ g ,and the variance σ g , where the index g runs over the number of Gaussian functions g = 1 , ..., N gmm . Thecomponent weights are constructed such that (cid:80) N gmm g =1 w g = 1. For example in one dimension the expressionof ˆ f is given byˆ f ( x ) = N gmm (cid:88) g =1 w j N ( x | µ g , σ g ) , N ( x | µ g , σ g ) = 1 √ π σ g exp (cid:18) − ( x − µ g ) σ g (cid:19) (1D) (30)The Gaussian mixture model has been paired with an expectation maximization algorithm [29, 30] to estimatethe values of µ g and σ g maximizing the probability of the particle data q j . The optimal number of Gaussian N gmm was decided using the Akaike information criterion [31].In Kernel Density Estimation (KDE), the shape of the kernel is variable, but the most widely used kernelsare PDFs themselves and thus have the desirable property that their integral is equal to one, (cid:82) ∞−∞ K ( y ) dy = 1.With a given kernel K ( y ), the distribution ˆ f can be constructed as [32]ˆ f ( x ) = 1 N ex ∆ x N ex (cid:88) j =1 K (cid:18) x − x j ∆ x (cid:19) , (31)12here ∆ x is the bandwidth. Options of the kernel K ( y ) relevant to our problem are: K ( y ) = 1(2 π ) d/ exp (cid:18) − y T y (cid:19) Gaussian (32) K ( y ) = 12 c − d ( d + 2)(1 − y T y ) Epanechnikov (33) K ( y ) = , for | y | < c d is the volume of a unitary sphere in d dimensions ( c = 2, c = 2 π , c = 4 π/ A ( K ) as A ( K ) = (cid:18) d + 1 (cid:19) d +4 Gaussian (35) A ( K ) = (cid:18) d ( d + 2)( d + 4)(2 √ π ) d (2 d + 1) c d (cid:19) d +4 Epanechnikov (36)The choice of the kernel should also be made considering such traits as differentiability and computationalefficiency [32]; those aspects will be discussed in detail in Sec. 3.
The transfer of information from the continuum plasma code to the discrete material code must be done insuch a way that the conservation laws are satisfied at the interface between the plasma and the material.The interface itself does not have to artificially produce mass, momentum, and energy, or in other words,the numerical method employed for the conversion from continuum to discrete (both forward and backward)does not have to produce artificial moments of the distribution functions up to a given order k . Consideringthe velocity moments up to third order at the surface, n s = (cid:90) f s d v , Γ s = (cid:90) v f s d v , (37) P s = (cid:90) m s vv f s d v , Q s = (cid:90) m s v v f s d v , at G ( x ) = 0 (38)13onsistency is ensured when the artificial moments generated by the conversion of the true distributions f s into their density estimates ˆ f s tend toward zero:∆ n s = (cid:90) ( ˆ f s − f s ) d v → , at G ( x ) = 0 (39)∆Γ s = (cid:90) v ( ˆ f s − f s ) d v → , at G ( x ) = 0 (40)∆ P s = (cid:90) m s vv ( ˆ f s − f s ) d v → , at G ( x ) = 0 (41)∆ Q s = (cid:90) m s v v ( ˆ f s − f s ) d v → , at G ( x ) = 0 (42)In order to enforce such consistency at the interface, the common factor ˆ f s − f s appearing in Eqs. 39-42must be numerically convergent. We considered the Mean Integrated Squared Error ( M ISE ) as a measureof optimal density estimate for such convergence,
M ISE ( ˆ f ) = E (cid:90) ( ˆ f ( x ) − f ( x )) d x (43) M ISE is simply the expected value of the integrated squared error. As the
M ISE measure is not adimension-free quantity, and as such it is not suitable to compare results of different dimensionality, weconstructed a dimensionless quantity
ISE ( ˆ f ) = (cid:82) ( ˆ f ( x ) − f ( x )) d x (cid:82) ( f ( x )) d x (44)by introducing the scaling factor (cid:82) f d x as originally suggested by Epanechnikov [33]. However, the errormeasure of Eqs. 43-44 can be used only when the true distribution f ( x ) is known, such as in analyticalverifications, making the ISE impossible to calculate for all cases where analytical solutions are not available.More generally for our case, a purely data-based estimation of M ISE and
ISE is desired, and may becalculated via bootstrapping [34]. In the present work, the following procedure has been used for thecalculation of the bootstrapped
M ISE (Faraway and Jhun [35]):1. Obtain the N ex samples of particles leaving the surface, q , q , ..., q N ex
2. Estimate the density function of the samples, ˆ f , using one of the kernels K ( y )3. Draw b bootstrap samples ( q ∗ , q ∗ , ... q ∗ b ) from ( q , q , ... q N ex ) with replacement4. Generate b bootstrap estimates ˆ f ∗ b from ( q ∗ , q ∗ , ... q ∗ b )14. Define the bootstrapped ISE , denoted (cid:100)
ISE , as: (cid:100)
ISE = (cid:80) bl =1 (cid:82) ( ˆ f ∗ l ( x ) − ˆ f ( x )) d x b (cid:82) ˆ f ( x ) d x (45)where the integral in the denominator is the normalization factor.If the density estimates used to construct ˆ f and ˆ f ∗ B are convergent, the estimate (cid:100) ISE tends toward the valueof the true
ISE , that is (cid:100)
ISE → ISE , and an estimate of the error may therefore be calculated without anyknowledge of the true distribution from which the original samples were drawn. The convergence propertiesof the (cid:100)
ISE have been numerically characterized in the next section as a function of the sample size of discreteparticles. The computational time for each numerical test is also reported for a direct comparison of thecomputational efficiency of the Density Estimation techniques adopted in this work.
For the dynamically coupled model intended for this work, density estimation must be performed frequentlyand on-the-fly, in the most demanding case at each time step of the Boltzmann solver. The possibility ofconstructing distributions on-the-fly by using particle data strictly depends on the convergence rate and onthe computational cost of the density estimation algorithm; a better convergence rate means that a smallernumber of Markov chains is necessary to reconstruct the distributions. Since this issue is of significantpractical importance, it will be addressed and discussed here in detail. In this section we report the numericaltests aimed at characterizing the convergence behavior and the timing of the coupling scheme of Sec. 2.Two types of tests were performed: (1) first, analytical tests of verification using samples drawn froma known Probability Distribution Function (PDF), so that the mean integrated square error (
M ISE ) ofeach estimate could be calculated analytically; (2) second, the estimation techniques were applied to dis-tributions of sputtered particles obtained from the surface erosion model of Sec. 2.2. Since in this secondcase the distributions of sputtered particles have no analytic solution, the (cid:100)
ISE was estimated through thebootstrapping algorithm described in Sec. 2.6. 15 .1.1 Analytical tests
The analytical tests were performed using samples drawn from 1D, 2D, 3D unbiased unimodal GaussianPDF distributions, which were reconstructed on a uniform grid of N d = 60 d gridpoints ( d = 1 , , N mc . The most evident trend in theconvergence rate plot (Fig. 2) is the “curse of dimensionality” of a probability distribution of dimension d .As expected [36], increasing the dimensionality of the problem both increases the amount of error in thesystem and decreases the rate of convergence. As can be seen from the same figure, increasing the samplesize from 10 to 10 in the one-dimensional Gaussian-KDE case decreases the error from 20% to 0.008%;in three dimensions, the same density estimate only decreases the (cid:92) M ISE from 30% to 0.3%. The GMMmodel has a less drastic decrease in convergence rate. However, while noticeably better with large samplesizes, the GMM performs worse in all three dimensions at a low sample size than the KDE methods; it onlyperforms better than the KDE when the sample sizes are greater than ∼
15, 30, and 100 in one, two, andthree dimensions, respectively. The computational time for the density estimates (Fig. 3) increases withboth the number of particles N mc and the dimensionality of the problem, with approximately an order ofmagnitude increase per dimension. Notably, the cost of the GMM model only marginally increases withdimensionality. At the largest sample size of 10 , the Gaussian-KDE is two orders of magnitude slower thanboth the Epanechnikov-KDE and the GMM, mainly due to the presence of the exponential factor, causingeach call to the Gaussian-KDE function to be more computationally expensive than the Epanechnikov-KDEfunction. Additional tests were performed to characterize the convergence of the bootstrapped (cid:100) ISE to theanalytical
ISE as a function of the number of samples. The tests revealed deviations smaller than 1% for asample size N mc >
30. With 50 samples, the accuracy falls below 0.5%. Similar tests revealed that approxi-matively 50 bootstrapped samples were sufficient for an estimate of the distribution of sputtered particles,as expected from a “well-behaved” smooth and unimodal underlying distribution.
Tests on distributions of sputtered particles (Fig. 4) were performed using data samples from the surfaceerosion model of Sec. 2.2, calculated with the F-TRIDYN code [2]. The sputtering simulations were run overa range of parameters of typical interest in nuclear fusion, considering a helium plasma incident on a tungstensurface at energies comprised between (cid:15) = 10 − θ = 0 ◦ − ◦ (ranging16etween normal incidence and grazing incidence). Sample sizes ranging from N mc = 10 to 10 were drawnfrom the particle data for the calculation of the bootstrapped (cid:100) ISE . Density estimates were performed alsoin this case using both the GMM (Eq. 30) and the KDE (Eqs. 31–33) methods. The distributions werereconstructed on a three-dimensional velocity grid of N d =3 = 60 grid points. Figure 5 (left) shows theresults of the numerical tests for the case of He impacting on W at energy (cid:15) = 1000 eV and angle θ = 60 ◦ ).The figure reports the bootstrapped (cid:100) ISE for the three density estimation techniques under analysis. Thetests show that the mean integrated square error of the GMM with small sample sizes is much larger than theKDE methods. The larger error of the GMM case can be attributed to the presence of both a tail in the dataand skewness of the distribution along the direction perpendicular to the surface, both of which representsignificant deviation from Gaussian distributions. Both KDE methods perform better in this case since theydo not rely on a parametric fit and can resolve these features with a smaller particle sample relative to theGMM. The plot also shows the convergence properties of the different techniques. Figure 6 summarizes alarge set of simulations at different energies and angles for two wall materials, beryllium and tungsten. Thefigure reports the number of samples required to obtain an (cid:100)
ISE of 10% for the estimation of beryllium andtungsten sputtered distributions. Incident angles between 0-80 degrees and energies of 150 eV, 500 eV, and1000 eV were considered for both materials. While GMM techniques require 500-1000 particle samples, thetwo KDE techniques have superior convergence, allowing to obtain an integrated (cid:100)
ISE error of the order of10% using a small sample of 40-80 sputtered particles. For such small particle samples of the order of 10 ,the computational time required for density estimation largely depends on the type of technique adopted(see Fig. 3), with the Epanechnikov-KDE being a best choice for the reconstruction of 3D distributions. The numerical tests revealed the superior convergence properties of Kernel Density Estimation methods overGaussian Mixture Models. KDE methods allow the reconstruction of the distributions with a smaller numberof particle samples. Three dimensional distributions of sputtered particles could be reconstructed with abootstrapped integrated square error (cid:100)
ISE smaller than 0.5% using a sample of just 50 particles. However,the two KDE methods differ considerably in computational time, with the Epanechnikov-KDE being up totwo orders of magnitude faster than the Gaussian-KDE. Such a difference is less dramatic in the range ofinterest of small sample sizes (approximately a factor of 2), but still showing that the Epanechnikov KernelDensity Estimator is the most appropriate method for on-the-fly reconstruction of distributions of sputteredparticles. 17 .2 Sheath formation in front of a wall releasing impurities
The methodology presented in the previous sections has been applied to the study of microscopic erosion ofmaterial walls exposed to a high-density magnetized plasmas. Here we report a numerical example, showinga 1D3V Boltzmann simulation of plasma sheath formation and depletion in front of a grounded wall in astrongly-magnetized environment. A magnetic field B almost parallel to the wall was chosen, with magneticangle ψ = 86 ◦ with respect to the normal to the surface. The simulation was for a helium plasma on aberyllium wall, and it was including four Boltzmann species ( s = 4): He + (He ions), He-I (He neutrals), Be + (Be ions), Be-I (Be neutrals). The plasma is allowed to propagate within a one-dimensional domain of sizeL = 1 mm in physical space and within ± v th in the velocity space, where v th is the thermal velocity of eachspecies, v th = ( k b T s /M s ) , with T s and M s are the temperature and mass of each species respectively. Theplasma was initialized with Maxwellian distributions for He + and He-I, and no Be species (Be + and Be-I)were initially present.Figure 7 shows a snapshot of the four phase spaces of each species during the plasma sheath formation.Several striking features have been observed. The helium ions (Fig. 7.a) are supersonically accelerated towardthe surface, and are responsible for the release of beryllium impurities in neutral state (Fig. 7.d). The Beneutrals are able to cross the Debye sheath and the magnetic presheath mostly remaining in neutral state.They ionize more frequently as they pass from the sheath region to the plasma bulk, turning into berylliumions (Fig. 7.c). Once ionized, the beryllium impurities become electrically and magnetically driven by theplasma; a fraction of the Be+ population flows back toward the wall (redeposition) and the remaining fractionflows toward the plasma bulk contributing to plasma contamination. Both Be + and neutrals are symmetricin the V thy and V thz dimensions, although a slight acceleration in v z of the Be ions due to the magneticfield is visible. One of the most striking physics features is related to the effect of magnetization. Themagnetization of the plasma causes an E × B acceleration along the v z direction, which is in turn responsibleof the modification of the energy-angle distribution functions at the wall. Finally, the velocity of the particlesincident on the walls is smaller than the classical value of 2 . V th , suggesting that the electrostatic portionof the sheath (Debye sheath) might decrease in size as the magnetic angle increases. In this work we have numerically characterized two classes of density estimation techniques useful for multi-scale coupling of kinetic plasma-material interaction models, namely the Gaussian Mixture Model technique18nd Kernel Density Estimation techniques. First, we have constructed a fully-kinetic model of the plasmamaterial interface, including: (1) a continuum 1D3V plasma solver of the multi-species Vlasov-Poisson-BGKproblem, and (2) a discrete ion-matter-interaction model based on Binary Collision Approximation. Then,we have constructed a coupling technique based on Density Estimation, and defined a consistent strategyto quantify the error introduced by such coupling, based on bootstrapping of a normalized Mean IntegratedSquare Error.From the numerical tests we have found that Kernel Density Estimation (KDE) techniques have superiorconvergence properties compared to Gaussian Mixture models, and are thus preferable for the treatmentof distribution functions exhibiting significant deviations from a Maxwellian, such as those encounteredin sputtering problems. Among the KDE methods, we have observed that Epanechnikov’s KDE methodoffers a good compromise between convergence rate and computational time, and it is thus the method ofpreference for distribution reconstruction. Finally, we have applied the kinetic model to a case of interest forfusion applications, simulating a multi-species strongly-magnetized plasma facing a material wall releasingimpurities.
This work was funded by the US Department of Energy through Grant No. de-sc0004736. Numerical testsutilizing density estimation techniques were performed in Python 2.7 with the scikit-learn library [37].
References [1] Derek Groen, Stefan J. Zasada, and Peter V. Coveney. Survey of Multiscale and Multiphysics Applica-tions and Communities. arXiv:1208.6444 [physics] , August 2012. arXiv: 1208.6444.[2] Jon Drobny, Alyssa Hayes, Davide Curreli, and David N. Ruzic. F-TRIDYN: A Binary CollisionApproximation code for simulating ion interactions with rough surfaces.
Journal of Nuclear Materials ,494:278–283, October 2017.[3] R. Dejarnac, M. Komm, J. Stockel, and R. Panek. Measurement of plasma flows into tile gaps.
Journalof Nuclear Materials , 382(1):31–34, November 2008.[4] Charles K. Birdsall and A. Bruce Langdon.
Plasma physics via computer simulation . McGraw-Hill,New York, 1985. 195] Rinat Khaziev and Davide Curreli. Ion energy-angle distribution functions at the plasma-materialinterface in oblique magnetic fields.
Physics of Plasmas , 22(4):043503, April 2015.[6] Gregory J. Parker and Nicholas G. Hitchon. Convected Scheme Simulations of the Electron DistributionFunction.
Jpn. J. Appl. Phys. , 36:4799–4807, July 1997.[7] Eric Sonnendrucker, Jean Roche, Pierre Bertrand, and Alain Ghizzo. The Semi-Lagrangian Method forthe Numerical Resolution of the Vlasov Equation.
J. Comp. Phys. , 149:201–220, 1999.[8] Chio-Zong Cheng and Georg Knorr. The integration of the Vlasov equation in configuration space.
Journal of Computational Physics , 22(3):330–351, 1976.[9] Jing-Mei Qiu and Andrew Christlieb. A conservative high order semi-Lagrangian WENO method forthe Vlasov equation.
Journal of Computational Physics , 229(4):1130–1149, 2010.[10] Petr Cagas, Ammar Hakim, James Juno, and Bhuvana Srinivasan. Continuum Kinetic and Multi-FluidSimulations of Classical Sheaths.
Physics of Plasmas , 24(2):022118, February 2017. arXiv: 1610.06529.[11] Francis Filbet and Chang Yang. Mixed semi-Lagrangian/finite difference methods for plasma simula-tions. arXiv preprint arXiv:1409.8519 , 2014.[12] Yasunori Yamamura and Yoshiyuki Mizuno.
Low-energy sputterings with Monte Carlo program ACAT .Institute of Plasma Physics, Nayoga University, Nagoya, Japan, May 1985.[13] Jon Drobny and Davide Curreli. F-tridyn simulations of tungsten self-sputtering and applications tocoupling plasma and material codes.
Computational Materials Science , 149:301 – 306, 2018.[14] Ludwig Boltzmann. Further Studies on the Thermal Equilibrium of Gas Molecules (Reprint). In
TheKinetic Theory of Gases , volume Volume 1 of
History of Modern Physical Sciences , pages 262–349. Im-perial College Press, July 2003. DOI: 10.1142/9781848161337 0015 DOI: 10.1142/9781848161337 0015.[15] Prabhu Lal Bhatnagar, Eugene P. Gross, and Max Krook. A model for collision processes in gases. I.Small amplitude processes in charged and neutral one-component systems.
Physical review , 94(3):511,1954.[16] G.I. Marchuck. Splitting and Alternating Direction Methods. In
Variational Methods for the NumericalSolution of Nonlinear Elliptic Problems , CBMS-NSF Regional Conference Series in Applied Mathemat-ics, pages 203–461. 2015. 2017] Gilbert Strang. On the construction and comparison of difference schemes.
SIAM Journal on NumericalAnalysis , 5(3):506–517, 1968.[18] Lukas Einkemmer and Alexander Ostermann. Convergence analysis of a discontinuous Galerkin/Strangsplitting approximation for the Vlasov–Poisson equations.
SIAM Journal on Numerical Analysis ,52(2):757–778, January 2014. arXiv: 1211.2353.[19] B. van Leer. Upwind and High-Resolution Methods for Compressible Flow: From Donor Cell toResidual-Distribution Schemes.
Commun. Comput. Phys. , 1(2):192–206, April 2006.[20] Ge-Cheng Zha and Chakradhar Lingamgunta. On the Accuracy of Runge-Kutta Methods for UnsteadyLinear Wave Equation.
American Institute of Aeronautics and Astronautics , 2003.[21] A. Jameson, W. Schmidt, and E. Turkel. Numerical Simulation of the Euler Equations by Finite Volumemethods Using Runge-Kutta Time Stepping Schemes. Palo Alto, California, 1981.[22] S. Devaux and G. Manfredi. Vlasov simulations of plasma-wall interactions in a magnetized and weaklycollisional plasma.
Physics of Plasmas , 13(8):083504, August 2006.[23] Peter E. Kloeden and Eckhard Platen.
Numerical Solution of Stochastic Differential Equations . Stochas-tic Modelling and Applied Probability. Springer-Verlag, Berlin Heidelberg, 1992.[24] J.F. Ziegler and J.P. Biersack.
The Stopping and Range of Ions in Matter, Treatise on Heavy IonScience , 6:95, 1985. cited By 1.[25] W. Mller and W. Eckstein. Tridyn a trim simulation code including dynamic composition changes.
Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials andAtoms , 2(1):814 – 818, 1984.[26] A. Mutzke and R. Schneider. Sdtrimsp-2d: Simulation of particles bombarding on a two dimensionaltarget version 1.0. ipp report 12/4 (garching, max-planck-institute for plasmaphysics, 2009.[27] J. Lindhard, V. Nielsen, M. Scharff, and P.V. Thomsen.
Mat. Fys. Medd., Dan. Vid. Selsk. , 33(10),1963.[28] Christian Robert and George Casella.
Monte Carlo Statistical Methods . Springer-Verlag, New York,2004. 2129] A.P. Dempster, N.M. Laird, and D.B. Rubin. Maximum Likelihood from Incomplete Data via the EMAlgorithm.
Journal of the Royal Statistical Society , 39(1):1–38, 1977.[30] Jonathan Q. Li and Andrew R. Barron. Mixture density estimation. In
Advances in neural informationprocessing systems , pages 279–285, 2000.[31] H. Akaike.
IEEE Transactions on Automatic Control , 19:716–723, 1974.[32] Bernard W. Silverman.
Density estimation for statistics and data analysis , volume 26. CRC press, 1986.[33] V.A. Epanechnikov.
Theory of Probability and its Applications , 14:153–158, 1969.[34] Bradley Efron. The Jackknife, the Bootstrap, and Other Resampling Plans. Technical Report 163,Stanford University, Stanford, California, December 1980.[35] Julian J. Faraway and Myoungshic Jhun. Bootstrap Choice of Bandwidth for Density Estimation.
Journal of the American Statistical Association , 85(412), December 1990.[36] Thomas Nagler and Claudia Czado. Evading the curse of dimensionality in nonparametric densityestimation with simplified vine copulas.
Journal of Multivariate Analysis , 151:69–89, October 2016.arXiv: 1503.03305.[37] Fabian Pedregosa, Gal Varoquaux, Alexandre Gramfort, Vincent Michel, Bertrand Thirion, OlivierGrisel, Mathieu Blondel, Peter Prettenhofer, Ron Weiss, Vincent Dubourg, and others. Scikit-learn:Machine learning in Python.
Journal of Machine Learning Research , 12(Oct):2825–2830, 2011.22igure 2: Integrated squared error for the density estimates as a function of sample size. Samples weredrawn and estimates were fitted in one through three dimensions. The solid lines are the one-dimensionalcalculations, the dashed lines are two-dimensional, and the dotted lines are three-dimensional.Figure 3: Computational time for the density estimates as a function of sample size. Samples were drawnfrom the standard Gaussian distribution in one (a), two (b), and three (b) dimensions.23igure 4: Integrated Squared Error compared to the bootstrapped Mean Integrated Squared Error ( (cid:92)
M ISE )in three dimensions. The dotted line is the (cid:92)
M ISE estimate, and the solid lines are the
ISE shown in Fig. 2.Figure 5: L2-norm of the Integrated Square Error ( (cid:100)
ISE ) of the continuum distributions reconstructed fromdiscrete data of sputtered particles, calculated using three different Density Estimation techniques. The x axis reports the number of sputtered particles, ranging from 10 to 10 .24igure 6: Comparison of the number of samples required to obtain an (cid:100) ISE of 10% for estimated berylliumand tungsten sputtered distributions. Incident angles between 0-80 degrees and energies of 150 eV, 500 eV,and 1000 eV were considered for both materials. 25 a) He Ions (b) He Neutrals(c) Be Ions (d) Be Neutrals
Figure 7: Phase space of He ions (a), He neutrals (b), Be ions (c), and Be neutrals (d) in a 1D3V magnetizedplasma after 30 ns. In each, the top figure is the X − V x plane, the middle is the X − V y plane, and the lastis the X − V zz