A Model for Characteristic X-Ray Emission in Electron Probe Microanalysis based on the Spherical Harmonic (PN) Method for Electron Transport

Classical k-ratio models, e.g. ZAF and phi(rho z), used in electron probe microanalysis (EPMA) assume a homogeneous or multi-layered material structure, which essentially limits the spatial resolution of EPMA to the size of the interaction volume where characteristic x-rays are produced. We present a new model for characteristic x-ray emission that avoids assumptions on the material structure to not restrict the resolution of EPMA a-priori. Our model bases on the spherical harmonic (PN) approximation of the Boltzmann equation for electron transport in continuous slowing down approximation. PN models have a simple structure, are hierarchical in accuracy and well-suited for efficient adjoint-based gradient computation, which makes our model a promising alternative to classical models in terms of improving the resolution of EPMA in the future. We present results of various test cases including a comparison of the PN model to a minimum entropy moment model as well as Monte-Carlo (MC) trajectory sampling, a comparison of PN-based k-ratios to k-ratios obtained with MC, a comparison with experimental data of electron backscattering yields as well as a comparison of PN and Monte-Carlo based on characteristic X-ray generation in a three-dimensional material probe with fine structures.


Motivation
Electron probe microanalysis (EPMA) (Heinrich & Newbury, 1991;Reimer, 1998;Ul-Hamid, 2018;Goldstein et al., 2018) terms a popular method to quantify the distribution of different chemical elements in heterogeneous materials at the micro-to nanometer scale, particularly fine structures and grains in metals and alloys. The EPMA method is an imaging technique that bases on solving an inverse problem (Tarantola, 2005;Aster et al., 2013), as the distribution of different chemical elements is not measured directly but reconstructed from calibrated characteristic x-ray intensities (k-ratios) measured in an electron microscope: Given a vector of experimentally measured k-ratios (k Exp ) for a material sample occupying the volume G, the inverse problem consists of finding the mass concentration fields c Z (x) with x ∈ G of the different chemical elements Z ∈ {1 (= ∧ H), 2 (= ∧ He), 3 (= ∧ Li), 4 (= ∧ Be), . . .} =: Z that best reproduce the experimental k-ratios for some k-ratio model k M : Problem 1 (Inverse Problem). Find c * such that with some appropriate norm · .
It is obvious, that the reconstructed concentration c inherits fundamental characteristics like accuracy or assumptions on the material structure from the forward model k M . Moreover, the forward model substantially influences the complexity and computational cost of solving the inverse problem as it can only be solved iteratively, preferable using efficient gradient based minimization techniques. Hence, the choice of the forward model k M (c) is crucial for both computational cost and accuracy of EPMA. Most models used nowadays are deterministic and very simple in the sense that they either assume that the material is homogeneous (Bastin & Heijligers, 1991;Pouchou & Pichoir, 1991) or are based on approximations of φ(ρz) curves (Bastin et al., 1998;Packwood & Brown, 1981;Philibert, 1963;Love et al., 1984;Pouchou & Pichoir, 1984), i.e. curves that model for the ionization probability as a function of depth expressed in mass thickness. On the one hand, these models and the resulting inverse problem are computationally cheap to solve. On the other hand, due to the strong assumptions on the material (homogeneous or varying only in depth), these models essentially limit the spatial resolution of EPMA to pixels that exceed the size of the interaction volume V ⊂ G where characteristic x-rays are generated due to ionization of atoms by high energy beam electrons. k-ratios for materials with structures that deceed the interaction volume are typically computed by Monte Carlo (MC) simulation of a large number of electron trajectories (Ammann & Karduck, 1990;Gauvin et al., 2006;Joy, 1991;Ritchie, 2005;Salvat, 2015). While MC simulations can exhibit high physical accuracy, they are not well-suited to solve the inverse problem as they are computationally expensive and suffer from statistical noise, which makes the application of gradient based minimization techniques difficult. Over the last decades the spatial resolution of EPMA has improved significantly by instrumental enhancements, e.g. improved focusing optics, of electron microscopes, such that nowadays the spatial resolution EPMA is no longer limited by the beam diameter but the interaction volume. The limitation in the spatial resolution is typically tackled using experimental measures to reduce the interaction volume, e.g. low overvoltage ratio (Newbury, 2002;Richter & Pinard, 2014) or low beam energies (Barkshire et al., 2000;Merlet & Llovet, 2012). We follow another approach to improve the spatial resolution of EPMA. Instead of reducing the interaction volume we aim to resolve the material within the interaction volume by developing forward models k M that a priori make no assumption on the material structure and allow to solve the inverse problem efficiently using gradient based minimization methods in combination with the adjoint state method for gradient computation. In this paper we present a new k-ratio model, which bases on a spherical harmonic P Napproximation of the Boltzmann equation for electron transport in continuous slowing down approximation, i.e. the electron transport is modeled by a set of partial differential equations that follow from a deterministic model order reduction of the Boltzmann equation by taking moments in direction space. Attempts on deterministic modeling of electrons inside solids to predict characteristic quantities, e.g. backscattered electrons and x-ray emission, used in quantitative EPMA have been made earlier (Fathers & Rez, 1979;Bünger et al., 2018). Similar approaches to modeling particle transport by moment equations are also common in other fields, e.g. nuclear engineering (Case & Zweifel, 1967), radiotherapy Küpper, 2016;Duclous et al., 2010;Pichard, 2016;Olbrant et al., 2012), rarefied gas dynamics (Grad, 1949;Levermore, 1996;McDonald & Torrilhon, 2013;Torrilhon & Struchtrup, 2004;Struchtrup & Torrilhon, 2003) and semi-conductors (Anile & Pennisi, 1992;Hauck, 2006). Our k-ratio model is a promising alternative to classical models, as P N moment models exhibit several properties that are desirable in terms of balancing accuracy and computational cost of the material reconstruction (problem 1) on scales smaller than the interaction volume: P N equations have a specific and simple structure that allows an efficient numerical solution, form a family of equations that is hierarchical in accuracy, allow to formulate accurate boundary conditions and are well suited for the adjoint state method of gradient computation.

Modeling of Characteristic X-ray Emission
In EPMA a material probe is usually modelled by the quantities of interest: Mass density ρ : G → R + and the mass concentration fields c Z : G → [0, 1]. The macroscopic model given by mass density ρ and concentrations c Z can be translated into a mesoscopic description given by the number density functions of the different atoms where A Z is the weight of a Z-atom. It is worth noting that this material model has no notion of correlation between atoms, in particular no notion of correlated arrangement of atoms like lattices or molecules. It is assumed that the material can be modelled isotropic, e.g. oriented atomic structures within the material are sufficiently small and randomly distributed or the orientation has a negligible impact. The emission of characteristic X-rays from the material probe due to excitation by an electron beam is the consequence of various physical processes on the atomic level. The processes are usually described by cross-sections of single free atoms σ Z and the extension to cross sections of compounds is performed by the additivity approximation (Salvat et al., 2007), i.e. the incoherent sum of atomic cross section, which we abbreviate by We also use a compact notation for electron shells and characteristic X-rays: We use tuple (Z, i) and (Z, j) to abbreviate the i th electron shell and characteristic X-ray j of a Z-atom. Our model of characteristic X-ray intensities I (Z,j) can be decomposed into four consecutive parts: 1. Computation of the fluence ψ of beam electrons inside the material 2. Determination of the ionization intensity of Z-atoms I Z 4. Calculating the X-ray absorption A (Z,j) before reaching the detector This results in the model The respective k-ratio model follows by normalization with characteristic x-ray intensities I Std. obtained for a standardized material under the same experimental conditions The critical part in modelling the emission of characteristic X-rays is to determine the electron fluence ψ(x, , Ω) = |v( )| n(x, , Ω), where v( ) denotes the electron velocity as a function of energy and n is the number density of beam electrons at point x ∈ G with energy moving in direction Ω ∈ S 2 . ψ is a fluence in the sense that ψ(x, , Ω)d dΩdAdt is the number of electrons which have an energy within [ , + d ], a direction of travel within the solid angle dΩ around Ω and traverse an area dA with normal Ω positioned around x during a time interval dt.
The intensity of (primary) characteristic X-ray generation induced by an electron fluence ψ is given by where σ ion. (Z,i) is the ionization cross section of a (Z, i) electron shell and w Z,i→j is the probability of (Z, j) X-ray emission after a (Z, i)-shell ionization. It is worth noting, that the electron fluence ψ enters eq. (8) solely by its zeroth angular moment, i.e. integrated over all directions, and hence no information on the angular distribution of the electron fluence is required. Compared to the mean free path of beam electrons the mean free path of the characteristic X-rays inside the material probe is large and interactions with the material are dominated by the photoelectric effect. Hence, characteristic X-rays essentially traverse the material along straight lines until leaving the material or being absorbed by an atom and ionizing it. Taking the X-ray absorption into account, the intensity of primary (Z, j) X-rays reaching a detector of small effective surface ∆A located around a point x far away from the interaction volume (distance L), can be approximated by where Σ abs. j denotes the absorption cross section.
From eqs. (7) to (9) it becomes clear, that modelling and computing ionization, and absorption is straight forward and the challenging part in modelling the emission of characteristic X-rays is to determine the total electron fluence Keeping in mind that the prediction of characteristic X-rays for a given material sample states a forward problem for the inverse problem (problem 1), it is essential that the model of the electron fluence ψ (or total electron fluence ψ (0) ) can be evaluated at low computational cost and is suited for modern gradient-based methods of solving inverse problems. In the following two sections we will motivate the method of moments as a model order reduction technique for the Boltzmann transport equation as well as present and discuss the particular case of P N moment models. Note that the above model misses secondary characteristic X-rays, i.e. characteristic X-rays emitted as a result of atom ionization by primary characteristic X-rays. However, secondary fluorescence can easily be included by augmenting the detector intensities I (Z,j) by F (Z,j) = A (Z,j) • Φ (Z,j) •Ĩ Z , whereĨ Z is the intensity of Z-atom ionization caused by primary X-rays.

Deterministic Modeling of Electron Transport
The electron fluence is governed by the steady state Boltzmann transport equation (Cercignani, 2012) in which the left hand side describes the free-flight of electrons and the two terms on the right hand side model elastic and inelastic collisions of electrons with material atoms. The collision operators describe the in-and out-scattering in phase space and are given by As inelastic scattering cross-sections are highly peaked around small energy losses, i.e. electrons most probably loose a significant amount of energy in a sequence of small energy losses, we apply the common continuous slowing down approximation: We split inelastic collisions into energy losses and deflections, model the energy loss as a continuous process governed by the average energy loss per path length (stopping power) and replace the inelastic scattering operator by its continuous slowing down approximation , −ω, µ) dω accounts for angular deflections in inelastic collisions. The motivation for the continuous slowing down approximation is to reduce the coupling in energy space by replacing the energy integral with a differential operator and removing the necessity of very fine discretization in energy space to resolve inelastic scattering cross sections. The continuous slowing down approximation transfers eq. (11) into an evolution equation in energy space, which we will call Boltzmann equation in continuous slowing down approximation (Larsen et al., 1997) where Σ := Σ inel. CSD + Σ el. . Due to the high dimension of the phase space (ψ is a function of six independent variables) and the coupling in direction space through the in-scattering term in the scattering operator Q (eq. (17)), it is computationally very expensive to solve BCSD (eq. (16)) numerically and we will apply the method of moments in order to further reduce the model order.
The BCSD model for the evolution of the electron fluence can be equipped with natural initial and boundary conditions. Initial conditions prescribe the fluence at the maximal energy max ψ(x, max , Ω) = ψ 0 (x, Ω).
Boundary conditions prescribe the incoming half of the angular distribution, i.e. at a boundary point x ∈ ∂G with outward-pointing boundary normal n the boundary condition is of the form with a given ψ in , that prescribes the angular distribution of the electrons entering the material domain G at all energies < max .

The Spherical Harmonic P N -Approximation
The spherical harmonic P N -approximation of the BCSD eq. (16) is based on spherical harmonic functions, which form a complete orthonormal set of functions on the sphere surface. A real spherical harmonic Y k l : S 2 → R of degree l ∈ N 0 and order k (|k| ≤ l) is given by where P |k| l denotes the associated Legendre function of respective degree and order, and θ and ϕ are the polar and azimuthal angle of direction Ω. These functions are well-studied, e.g. (Müller, 1966;Atkinson & Han, 2012).
In the P N -method the electron fluence ψ will be replaced by the moments In order to keep the notation compact, we introduce notation for vector functions that gather spherical harmonic functions, either of a certain degree l or up to a certain degree NῩ and with this the fluence can be approximated up to an order N by Spherical harmonic functions are a very convenient basis for the discretization of the angular distribution in eq. (16). For one, spherical harmonic functions are eigenfunctions of the scattering operator Q (17) which reduces the computational cost of evaluating the discrete scattering operator from a O(N 2 ) operation to a O(N ) operation. Note that the eigenvalues q l depend only on the degree l. Also, the function space spanned by spherical harmonics of degree l is rotational invariant which allows to maintain the rotational invariance of the BCSD in the P N approximation. Furthermore, spherical harmonic functions satisfy a recursion relation of the form where matrices A (i) , having only up to 4 non-zero entries per row, are very sparse. Moreover, a spherical harmonic Y k l is either even or odd in the direction of a Cartesian axes, where we call a function f : It is clear, that two spherical harmonics can couple in the recursion relation eq. (26) only if their even-odd classification differs in direction of the i-th axes and complies in the direction of the two other two coordinate axes. Based on the classification of spherical harmonics into even and odd moments we distinguish the moments in u P N into even and odd moments, which we denote by u e,n P N and u o,n P N . For example, a spherical harmonic Y k l is odd/even in the z-direction if the degree l is odd/even, such that the even/odd moments u e/o,n P N for direction n = e 3 and odd approximation order N > 3 are The even/odd classification of the moments is crucial to the formulation of our P N boundary conditions (Bünger et al., 2020). The coupling structure between even and odd spherical harmonics in the recursion relation (26) can be exploited to derive an efficient higher order discretization of the transport operator using staggered grids and central finite differences, see appendix A.

Derivation of P N Moment Model Equations
The P N approximation of the BCSD (eq. (16)) can be derived by a spectral Galerkin discretization (Canuto et al., 2006) of the fluence ψ in direction space. We replace the fluence ψ in the BCSD eq. (16) by the expansion ψ P N from eq. (24) in spherical harmonics up to an order N and derive equations for the expansion coefficients u P N . This is done by weakly enforcing the BCSD for ψ P N on the test space spanned by the spherical harmonics As indicated by the underbraces in eq. (30), we can use the orthonormality and recursion relation (26) of spherical harmonics to write the P N model as Equation (31) is a linear hyperbolic system of partial differential equations that models the evolution of the spherical harmonic moments (21) in energy space.

Initial and Boundary Conditions
Analogous to the derivation of P N moment equations, we obtain initial condition for the moments u k l by multiplication of the BCSD initial condition (18) with the weighting functions Y k l and integration over all directions The derivation of boundary conditions is more evolved and different strategies to translate the BCSD boundary condition (19) into boundary conditions for the P N equations exist in the literature (Mark, 1944(Mark, , 1945Marshak, 1947;Egger & Schlottbom, 2012;Schlottbom & Egger, 2015). In (Bünger et al., 2020) we describe the derivation of boundary conditions that are consistent with the characteristic waves of the P N equations, uniquely determine the incoming characteristic waves and assure an energy bound of the solution u P N . The boundary conditions follow from testing the BCSD boundary condition (19) by spherical harmonics that are odd in the direction of the boundary normal. For a point x ∈ ∂G with boundary normal n at an energy < max the boundary conditions take the following form where the boundary matrix B n N is related to the transport matrices A (i) N,N , see (Bünger et al., 2020) for details. The source term g n in is given by half moments of the incoming electron fluence ψ in with odd weighting functions

Comparison to M N Moment Models
It is fair to say that the P N moment model is a "simple" model in comparison to the popular family of entropy-based moment models (M N ), that are obtained by closing moment equations by minimizing the Boltzmann entropy functional which leads to accurate models even for low order moment approximations. At the same time the simple structure of P N equations is the major advantage over M N moment models. Increasing the moment order N corresponds to augmenting the moment system by additional equations of the same simple structure and therefore the computational cost of solving the P N equations primarily depends on the moment order through the number of equations (N + 1) 2 . Increasing the model accuracy by simply increasing the moment order in the P N equations therefore comes at a reasonable additional computational cost. Furthermore, the simple structure of the P N equations can be exploited to build efficient numerical solvers. In contrast, the solution of higher order M N model is computationally very expensive, as analytical solutions of the minimization problem in the closure problem is known only for the first order model M 1 (Mevenkamp, 2013;Pichard, 2016;Bünger et al., 2018) and higher order M N models require the numerical solution of many minimization problems to simulate the evolution of the moments in energy space (Garrett et al., 2015). The simplicity of P N equations also facilitates the analysis and allows to define compatible and accurate boundary conditions, which is essential in the context of EPMA considering that the beam electrons enter the material over a free material surface.
It is well known that P N solutions may exhibit oscillations around discontinuities and large gradients, if the expansion is truncated at too small N and the neglected higher order moments are still significant. In some cases the oscillations can even lead to negative electron concentrations. The fact that P N models do not preserve realizability of moments, i.e. P N models can produce sets of moments that do not follow from a non-negative distribution, is considered a major drawback over M N models. It is worth mentioning, that discretized M N equations not necessarily assure realizability and preserving realizability in the numerical solution of M N equations is part of current research. A common approach to improve the accuracy of spectral approximations around unresolved gradients is to dampen, or filter, higher order terms in the expansion. In many cases where a classical P N solutions suffer from oscillation for too small order N a filtered P N approximation, denoted by F P N , gives good results and hence allows to reduce the computational cost.
A convenient strategy to apply filtering is to simply modify the transport coefficients in the P N equations, e.g. (Küpper, 2016), with a smooth filter function f with f (0) = 1 and parameter α such that f (a) (0) = 0 for a = 1, . . . , α ∈ N in order not to dampen the zeroth order moment but primarily higher order terms of the expansion. The transport coefficients become where the filtering cross-section σ f is an additional tuning parameter. A typical filter function is strictly decreasing on [0, 1], like e.g. the exponential filter (Küpper, 2016) f where M is the machine accuracy. An alternative filtering strategy is to augment the numerical integration scheme by a filtering step on the spherical harmonic expansion at the end of an update. One can show that a filtering step of the form u k l → f ( l N +1 ) σ f ∆ u k l leads to a discretization that in the limit ∆ → 0 is equivalent to the modification (36), which actually led to the filtering strategy of modifying transport coefficients.

Testing Model Reduction by Moments
The idea of the method of moments is to reduce the phase space by replacing the direction Ω as an independent variable by a finite, preferably small, set of moments, in the hope to still be able to accurately model the electron fluence or, more importantly for our case, the total electron fluence ψ (0) . We demonstrate that this hope is indeed well justified by looking at the multiple-scattering distribution of the direction of travel for an electron moving in a homogeneous material. As the diffusion of electrons is a multiple-scattering process dominated by elastic collisions we neglect inelastic collisions for the moment. The probability distribution for the direction of travel Ω after one collision is given by where µ(Ω) =Ω · Ω = cos(θ), with deflection angle θ, measures the deviation from the initial directionΩ andˆ is the energy of the electron. The model reduction by the method of moments can now be motivated by inspecting the probability distribution f n (µ) for the direction of travel Ω after exactly n collisions or alternatively the distribution F s (µ) Cu computed with the ELSEPA code system (Salvat et al., 2005) and a fitted Wentzel scattering cross section (Salvat et al., 2007); And (b) the multiple scattering distribution F (eq. (40)) for the Wentzel cross section at different path length s (measured in mean free path length λ).
after a path length s (Goudsmit & Saunderson, 1940;Salvat et al., 2007) where P l denotes the Legendre polynomial of degree l. We will use these realistic distributions to test the quality of the moment expansion (24). Note that the coefficients in the series expansions of f n and F are determined by moments of the angular distribution f 1 with Legendre polynomials as weighting functions. We know that K 0 = 1, G 0 = 0 and that |K l | < 1 and G l > 0 for l ≥ 1 (P 0 (µ) = 1 and ∀l ≥ 1 : |P l (µ)| < 1, ∀µ ∈ (−1, 1)) and for realistic scattering cross section both series converge, i.e. the coefficients K l approach zero for large enough l. Hence, as the number of collisions n or the path length s increase the angular distribution is more and more dominated by the first terms, which can be modelled by low order moments. In the limit of infinite number of collisions or infinite path length the angular distributions become isotropic: lim n→∞ f n (µ) = lim s→∞ F (s; µ) = (4π) −1 .
In order to confirm that the essential character of a multiple-scattering angular distributions can be captured by a small number of moments we computed the multiple scattering distributions f n and F for an electron traveling in pure copper atˆ = 20keV. To avoid numerical difficulties in the computation of the moments K l for large l, we used analytical moments for a Wentzel scattering cross section (Salvat et al., 2007) that was fitted to a scattering cross section computed with the ELSEPA code system (Salvat et al., 2005). See Figure 1a for a comparison of the two scattering cross section. Figure 1b shows the angular distribution F obtained with the first 10 and first 150 terms in the series expansion (40) at three different path length measured in mean free path length λ = 1/σ tot. . We can observe the relaxation of the initially very peaked distribution towards an isotropic distribution and that the angular distribution is very well captured by only 10 moments already at a moderate path length of 30λ. Even at a path length of 10λ the first 10 moments give a reasonable approximation of the angular distribution. Similar observations can be made for the angular distribution after exactly 5, 15, 25 and 50 collisions plotted in Figures 2a to 2d. Except for a very small number of collision like 5, the angular distribution after 15 and more collisions is fairly well modelled by solely the first 13 moments. Furthermore the graphs confirm that the number of moments required to approximate the angular distribution decreases rapidly with the number of collisions.  This section contains simulation results of five test cases, in which we 1. compare P N solutions of electron transport to M 1 and MC solutions, 2. illustrate the accuracy of our P N boundary conditions, 3. compare our P N k-ratio model with k-ratios obtained using MC, 4. compare electron backscattering yields obtained with P N to experimental data and 5. compare P N to MC based on characteristic X-ray generation intensities in a fictive 3-dimensional sample with fine structures.
Unless stated otherwise we use the following setting. The P N results we obtained using stopping powers from the ICRU Report 77 (Salvat et al., 2007), and elastic scattering cross section generated with ELSEPA code system (Salvat et al., 2005). The presented MC results were obtained from sampling at least 10 7 electron trajectories using the DTSA-II MC code (Ritchie, 2009) with the following sampling strategies for the electron trajectories: Bethe1930Strict (Hans, 1930) with mean ionization potentials, Sternheimer64 (Berger & Seltzer, 1964) for the energy loss of electrons and NISTMottScatteringAngle (Jablonski et al., 2003) for scattering of electrons. In MC the total electron fluence is obtained from a histogram in energy and space of the sampled electron trajectories. Generation and absorption of characteristic X-rays were computed using absolute ionization potential Casnati82 (Casnati et al., 1982), shell transition probabilities Relax (Cullen et al., 1997) and X-ray mass absorption coefficients Chantler2005 (Chantler, 1995(Chantler, , 2003 extracted from the DTSA-II (Ritchie, 2009) MC code. The k-ratios of eq. (5) are obtained using pure materials as standard. Finally, we assumed a Gaussian electron beam in negative z-direction at the material surface (x 3 = 0) centered around x 1 = x 2 = 0 where B is the beam energy, σ D , σ and σ Ω parametrize the width of the beam in space, energy and angle, andΩ = −e 3 is the average direction of travel.

Comparison of Electron Transport Models
Cr Ni B = 10keV e − Figure 3 In this test case we consider a two-phase material consisting of pure chromium and nickel separated by a vertical interface. The material is excited by a narrow 10keV Gaussian electron beam (σ = 0.05keV, σ D = 10nm, σ Ω = 0.1) that is aligned with the interface separating the two material phases, see Figure 3. In Section 6.1 we compare different approximations of the total electron fluence ψ (0) at different energy levels: Three P N solutions of different order (P 5 , P 13 and F P 5 with the exponential filter (37) using α = 6 and σ f = 100nm −1 ), a M 1 solution and a MC solution as reference. The results demonstrate several of the aforementioned differences between the models.
The P 5 solution clearly deviates from the MC results. Especially initially, meaning for energies close to the beam energy, the P 5 solution shows unphysical oscillations with negative electron concentrations in some parts, as the Gaussian beam is too peaked to be resolved well by moments up to order 5 only. The solution of the M 1 model does not suffer from unphysical negative electron concentrations nor oscillations and, considering that the M 1 model contains only moments up to first order, reproduces essential features of the MC solution quiet well, e.g. the area of highest electron density and the volume occupied by electrons. However, the overall shape of the M 1 solution differs from the MC solution at 8.5keV and 7.5keV, in particular: The relaxation of the total electron fluence towards the centre of the interaction volume seems to be too slow such that, in comparison with the MC solution, the shape of the M 1 solution is more peaked close to the lower boundary of the interaction volume. We can increase the accuracy of the P 5 approximation by damping the highest moments, i.e. using the F P 5 approximation, or simply exploiting the hierarchy of P N models and increase the moment order N in the P N approximation. In fact, the unphysical oscillations of the P 5 solution are not present in the F P 5 solution and the predicted total electron fluence is in good agreement with the MC solution. Likewise, the P 13 solution almost looks like a noise-free version of the MC solution. It is worth mentioning that we chose a very narrow Gaussian electron beam to stretch the differences between the moment models. For a broader beam, lets say 50nm, solutions of the low order moment models P 5 and M 1 would already approximate a MC result better. Note, that deviations in the total electron fluence at certain energies are not necessarily representative for deviations in the resulting shell ionization intensities or characteristic X-rays generation intensities, as deviations in the fluence at different energy levels can cancel out in the integral over all energies. Figure 5 As beam electrons enter the material through a polished material surface and possibly leave the material over the same surface back into the vacuum chamber, it is clear that accurate boundary conditions are a crucial ingredient to an overall accurate description of the beam electrons. Especially initially, i.e. for energies close to the beam energy, the electron fluence strongly depends on the translation of the incoming electron distribution modelling the electron beam into boundary conditions for the P N field equations (31).

P N Boundary Condition
Furthermore, when the beam electrons enter the material the distribution in direction space is very peaked around the direction of the electron beam, and the approximation by a spectral method like P N can be challenging in the sense that a large approximation order N might be necessary. In this test case we consider a 2-dimensional material sample of pure copper excited by a Gaussian electron beam of 14keV see Figure 5. We aim to show that the boundary conditions proposed in (Bünger et al., 2020) are not only convenient from a mathematical point of view, but also allow to accurately model a Gaussian electron beam. We address the accuracy of the boundary conditions by comparing the total electron fluence obtained with P 13 to a reference solution obtained with MC. For comparison we have set the stopping power to a constant value of S = −2 · 10 −13 (J/m)/(kg/m 3 ) in the both, the P 13 and MC, computation. Snapshots of the total electron fluence along the z-axis for a narrow Gaussian beam of (σ = 0.1keV, σ Ω = 0.1, σ D = 50nm), are shown in Figure 5. We observe that the P 13 total electron fluence agrees very well with the MC result at the different energies. Indeed, one can hardly distinguish the P 13 solution from the MC solution.

Electron Backscatter Yield
Besides validating PN models by comparison with MC, we also aim to validate the applicability of the P N approximation by comparing with experimental data. This is however not easy: It is not possible to measure the fluence of beam electrons inside a material probe experimentally and therefore we can only try to disproof applicability of an electron fluence model by comparing with experimentally measurable quantities determined or influenced by the electron fluence; Furthermore, experimental data are not guaranteed to be more accurate than computational results obtained with a mathematical model and hence deviations do not necessarily indicate a poor model. In this test case we consider the backscattering (BS) yield η of beam electrons, which is defined as the fraction of beam electrons that re-emerge from the surface of the material sample. Backscattered electrons, just like characteristic X-rays, form a signal used for the acquisition of BSE images supporting EPMA. Given a P N solution for the electron beam fluence (ψ P N ) we obtain the respective P N BS yield from where A ⊂ ∂G denotes the material surface and n the surface normal pointing outside the material.
In Figure 7 we compare BS yields of carbon (C), copper (Cu) and gold (Au) obtained with the P 13 moment model to experimental measurements extracted from a database (Joy, 2008). For moderate and high beam energies the BS yields computed with P 13 are in very good agreement with the experimental measurements for all three elements. The evaluation for low beam energies is difficult as the experimental data show strong deviations from each other. However, the P 13 results lie within the range of the different measurements and thereby have to be more accurate than some of the measurements. It is worth mentioning, that the backscattering yield depends strongly on the scattering cross sections. Hence, deviations from the measurements are not necessarily caused solely by modeling errors introduced by the P N approximation but also inaccuracies in the scattering cross-section.  We now compare k-ratios obtained with our k-ratio model (eq. (5)) using the P 13 moment model for electron transport to k-ratios computed with MC. As test case we consider a thin iron layer of height h ∈ [50nm, 200nm] on top of a silicon substrate (see the schematic drawing in Figure 8) excited by electron beams of different configuration: Two different beam energies ( B ∈ {10keV, 20keV}) and various beam positions. For a fair comparison of the P N and MC simulation we simplify the discretization of the material: We assume that the height of the iron layer varies very little on the scale of the interaction volume, such that it can locally, i.e. for a given beam position, be approximated by a constant height. Note, that the height of the iron layer is of same order as the beam diameter and much smaller than the interaction volume, i.e. the layer can be seen as a fine material structure. Also note, that the material properties change considerably at the interface between the iron and silicon phase.

Layer System
The results for Kα X-rays are plotted in Figure 9. For both beam energies the P 13 and MC k-ratios for iron an silicon match very well over the whole range of layer heights. It seems that the MC solutions deviates spuriously from the very systematic course of the k-ratio graphs obtained with the deterministic P N model. For the low beam energy of 10keV the MC k-ratio for iron exceeds one, which is unrealistic as we obtained k-ratios by normalization with characteristic X-ray intensities of pure materials. We suspect, that the deviations between the P 13 and MC results stem from stochastic noise of the MC method. Fe (10keV) Si (10keV) Fe (20keV) Si (20keV) h

[nm]
Kα k-ratio P13 MC Figure 9: Kα k-ratios obtained for a Fe-layer on a Si-substrate for different layer thicknesses h and two beam energies (10keV and 20keV) using the P 13 moment model and MC.

Characteristic X-ray Generation Intensity in a Small Iron Particle
We now compare the P N approximation of electron transport to MC sampling of individual electron trajectories based on characteristic X-ray generation intensities Φ (Z,·) . The test case we consider is a fictive 3-dimensional cuboid-shaped material probe (G = [−750nm, 750nm]×[−750nm, 750nm]×[−750nm, 0nm]) that contains six small iron particles embedded in a pure copper matrix. The probe is discretized by a Cartesian 66×66×33 material grid consisting of small rectangular material cells with an approximate volume of (22.7nm) 3 . We assume a beam of 20keV electrons (σ D = 100nm, σ E = 0.2keV, σ Ω = 0.1) that is focused onto the point (0, 0, 0) T ∈ ∂G. The material structure is visualized by Figure 10, which shows a P 13 result of the quantity of interest in this test case: The (Fe,Kα) X-ray generation intensity Φ (F e,Kα) caused by beam electrons with energy ≥ 15keV. We will demonstrate a fundamental advantage of deterministic moment models like P N over MC sampling of individual electron trajectories: The absence of stochastic noise in the characteristic X-ray generation intensity Φ (Z,·) inside the material probe. We see two difficulties of noisy MC electron trajectory sampling for EPMA imaging with a high resolution, i.e. when the material grid is very fine. Firstly, to obtain reliable intensities Φ (Z,·) in each material cell the number of sampled trajectories must be high enough to assure that the statistical noise in the material cells is small, i.e. each material cell is crossed by sufficiently many trajectories. The likely hood for an individual trajectory traversing a material cell reduces with the cell size and the distance from the focus point of the electron beam on the material surface. Hence, to reduce statistical noise in small material cells located near the boundary of the interaction volume requires sampling a very large number of electron trajectories. This can make the MC approach computationally too expensive for high-resolution EPMA material imaging, i.e. for solving the inverse problem (problem 1) with fine material grids. Secondly, the statistical noise greatly complicates the use of gradient based minimization methods in solving the inverse problem: Even when the number of electrons is sufficiently large to obtain reliable Figure 10: Normalized Fe-Kα X-ray generation intensity Φ (Fe,Kα) in a 3-dimensional material probe consisting of six small iron particles embedded in a pure copper matrix indicated by the gray background; Φ (Fe,Kα) computed with the P 13 moment model of electron transport.
intensities Φ (Z,·) in each material cell, the computation of reliable gradients using conventional methods like finite-differences, which greatly amplify noise, might not be possible and require sampling an even higher number of electron trajectories.
We demonstrate the before mentioned drawback of noisy or computationally expensive MC trajectory sampling by looking at the characteristic X-ray generation intensity Φ (Fe,Kα) inside the iron particle located around the point (200nm, −250nm, −325nm) ∈ G, i.e. the rightmost iron particle in Figure 10. We computed the intensity Φ (Fe,Kα) by MC sampling of n e = 10 p , p = 3, 4, 5, 6 electron trajectories. The results are visualized in Figures 11a to 11d and the respective computing times obtained on a standard desktop computer with an Intel i7-8700K CPU are given in the captions. Note that the simulation of 10 6 electron trajectories still corresponds to a small sample size, as even at a low beam current of only 1nA more than 10 10 beam electrons penetrate the material per second. From Figure 11a it becomes clear that for the lowest number of only 10 3 trajectories many material cells in the particle are not crossed by any trajectory and therefore would not contribute to detector signals. In Figure 11b we see that when the number of trajectories is increased to 10 4 almost all material cells constituting the particle are crossed by at least one trajectory. However, the intensity Φ (Fe,Kα) is almost dominated by statistical noise as the number of trajectories traversing the material cells is still small. By further increasing the number of trajectories to 10 5 the noise level is sufficiently reduced to guess the underlying intensity Φ (Fe,Kα) , see Figure 11c. Nevertheless, the pollution of the intensity Φ (Fe,Kα) by stochastic noise is still significant. The intensity Φ (Fe,Kα) obtained for 10 6 trajectories, visualized in Figure 11d, indicates that the reduction of the stochastic noise to an acceptable level might require to sample a very high number of electron trajectories with long computing times on standard desktop computers.
We simulated the same test case with our research P N solver implemented in Matlab. Figure 12a shows an intensity Φ (Fe,Kα) obtained with the P N model. The result demonstrates the advantage of a deterministic moment model: Fast computation of completely (a) n e = 10 3 (0min 3s) (b) n e = 10 4 (0min 23s) (c) n e = 10 5 (3min 33s) (d) n e = 10 6 (42min 37s) Figure 11: Contour plots of normalized characteristic X-ray generation intensity Φ (Fe,Kα) in an iron particle computed from MC sampling of n e = 10 p , p = 3, 4, 5, 6 electron trajectories.
smooth intensities Φ (Fe,Kα) . The efficiency of our solver bases on two parts. For one we vectorized computing intense parts of the solver to make use of Matlabs fast matrix operations. Secondly, we use mesh and moment order adaptivity to employ that, in the context of EPMA, the evolution of the electron distribution follows a pattern: An initially very peaked electron distribution diffuses into a volume whose dimensions are usually several times larger than the beam diameter. This pattern can be exploited by starting with a high moment order N on a small and fine computational mesh, and increasing the mesh size while reducing the mesh resolution and moment order N during a simulation. The results shown in Figure 12a were obtained by starting with a high initial approximation order of N = 11 to resolve the peaked distribution at energies close to the beam energy and, while the electrons diffuse into the interaction volume, gradually decreasing the order to N = 6 by dropping higher order moments when their contribution becomes insignificant for the evolution of the total electron fluence.
(a) P 11 -P 6 (0min 29s) Figure 12: Contour plot of normalized characteristic X-ray generation intensity Φ (Fe,Kα) in an iron particle computed with P N moment model.

Conclusion and Outlook
The presented simulation results demonstrate the potential of improving the spatial resolution of EPMA by omitting a-priori assumptions on the material structure using k-ratio models that base on deterministic moment models of electron transport. The P N model of electron transport is particularly promising, as its simple structure allows to implement very efficient numerical solvers and to exploit the hierarchy of the moment system to adjust accuracy at an acceptable increase of computational cost. The comparisons of P N results to reference solutions generated with Monte Carlo showed a good level of accuracy already at moderate moment orders N and the possibility to essentially reproduce Monte Carlo solutions without the statistical noise and reduced computational cost.
In the near future we will widen our comparisons with experimental data. Furthermore, we will implement a gradient-base minimization algorithm to solve the inverse problem of EPMA material imaging. The P N moment model together with the boundary conditions presented in (Bünger et al., 2020) are well suited for the formulation of adjoint equations, which we will exploit by using the adjoint-state method for gradient computation in the minimization algorithm. We plan to study various aspects of the inverse problem, e.g. the sensitivity of different experimental parameters like the beam energy, diameter or displacement on the convexity of the error function k M (c) − k Exp . Our long term goals include: Propagation of uncertainties in experimental measurements into uncertainties in the material image, studying the impact of modelling assumptions like the continuous slowing down approximation and investigation of inevitable uncertainties in modelling X-ray emission introduced by uncertainties in fundamental cross-sections e.g. shell ionization or X-ray absorption cross-section.

Acknowledgements
The authors acknowledge funding of the German Research Foundation (DFG) under grant TO 414/4-1.

A Numerical Method
In this section we present a discretization of the P N moment approximation of BCSD that is inspired by the discretization of P N equations of radiative transfer proposed in (Seibold & Frank, 2014). The discretization exploits the special coupling structure between odd and even moments in P N equations using staggered grids and central finite difference for the discretization in space. Here, we sketch the discretization for the 1-dimensional case as the extension to the 3-dimensional case is straight forward. An alternative discretization using summation-by-parts (SBP) finite difference (FD) operators on staggered grids for the discretization of spatial derivatives (Kreiss & Scherer, 1974;Strand, 1994;Olsson, 1995;Fernández et al., 2014) and simultaneous approximation terms (SAT) (Carpenter et al., 1994) to include boundary conditions can be found in (Bünger et al., 2020). For the discretization of the material domain we choose a Cartesian coordinate system such that the origin sits on the polished surface of the material and the negative z-axis enters the material orthogonal to the surface. When grouping even and odd moments of u P N the 1-dimensional P N equations in z-direction take the form where the vector function u o/e contains the moments of spherical harmonics Y k l that are odd/even in the z-direction, i.e. the degree l is odd/even. The boundary conditions take the form As proposed in (Seibold & Frank, 2014), we exploit that even and odd moment couple solely through their derivatives in the transport operator: We discretize even and odd moments on grids that are staggered to one another (see Figure 13), such that we can approximate spatial derivatives of odd/even moments appearing in the evolution equation of even/odd moments by central finite differences. This discretization is efficient in the sense that central FD provide second order accuracy at the cost of a forward/backward first order finite difference operation. The semi-discrete system for internal grid points is given by where S h o/e and Q h o/e denote the grid functions obtained by restriction of S and Q on the odd/even grid points, u h denotes the numerical approximation of u and the indices i and z z 0 = z L ∆z z i− 1 2 z i z i+ 1 2 z N = z R Figure 13: Illustration of one-dimensional staggered grid: Red circles accommodate even moments (even grid) and blue crosses accommodate odd variables (odd grid). Note that the odd grid has two ghost points lying outside the domain [z L , z R ].
j denote function values at grid nodes, e.g.
and (Q e h e ) j = Q e h e (z j , · ). The semi-discrete evolution equation (47) of the even moments requires the odd moments at the ghost points. We prescribe the odd moments at the ghost points, such that the boundary conditions are satisfied for interpolated odd moments, e.g. at z = z L Note that this discretization is equivalent to using odd grids with, instead of external ghost nodes, additional boundary nodes at which the odd moments are prescribed by evaluating the right hand side in the boundary conditions (45) and using forward/backward finite differences to approximate the spatial derivatives of odd moments at the boundaries points z L and z R , e.g. at z = z L For the numerical integration in energy space we use the Strang splitting (Strang, 1968) strategy as proposed in (Seibold & Frank, 2014): The update operator S ∆ for an energy step ∆ is a composition of even and odd update operators whereby here the notion of odd and even is in terms of the classical parity of spherical harmonics, i.e. a moment u k l is termed odd/even if the degree l is odd/even. Before an update → + ∆ the energy-dependent stopping power S and transport coefficients q l are fixed by evaluation at the mid-energy n + 1 2 ∆ . For the update of odd moments the even moments are fixed. Likewise the odd moments are fixed for the update of the odd moments. For fixed even/odd moments and evolution equations for the odd/even moments turn into ODEs that can be solved analytically, see (Seibold & Frank, 2014) for the analytical solution and update operators S o and S e . Compared to alternative time integration schemes, e.g. iterative Runge-Kutta methods, this splitting scheme is very memory efficient as it doesn't rely on storing intermediate updates. Since the memory consumption depends quasi quadratically on the approximation order N , this is particularly important for P N approximations, where it is not uncommon to use a moderate or even a high moment order N .