A Novel Reconstruction Method to Increase Spatial Resolution in Electron Probe Microanalysis

: The spatial resolution of electron probe microanalysis (EPMA), a non-destructive method to determine the chemical composition of materials, is currently restricted to a pixel size larger than the volume of interaction between beam electrons and the material, as a result of limitations on the underlying k-ratio model. Using more sophisticated models to predict k-ratios while solving the inverse problem of reconstruction offers a possibility to increase the spatial resolution. Here, a k-ratio model based on the deterministic M1-model in Boltzmann Continuous Slowing-Down approximation (BCSD) will be utilized to present a reconstruction method for EPMA which is implemented as a PDE-constrained optimization problem. Iterative gradient-based optimization techniques are used in combination with the adjoint state method to calculate the gradient in order to solve the optimization problem efﬁciently. The accuracy of the spatial resolution still depends on the number and quality of the measured data, but in contrast to conventional reconstruction methods, an overlapping of the interaction volumes of different measurements is permissible without ambiguous solutions. The combination of k-ratios measured with various electron beam conﬁgurations is necessary for a high resolution. Attempts to reconstruct materials with synthetic data show challenges that occur with small reconstruction pixels, but also indicate the potential to improve the spatial resolution in EPMA using the presented method.


Increasing the Spatial Resolution in EPMA
Electron probe microanalysis (EPMA) [1,2] is a non-destructive method to determine the chemical composition of a material sample that is applied in various fields, for example geophysics, material science, or engineering. The composition is not measured directly, but reconstructed from the intensity of X-rays, which are produced by exciting a material sample with a focused beam of electrons. Atoms inside the sample are ionized by collision with the highly energetic beam electrons and subsequently relax from the excited state by emitting characteristic X-ray radiation. A part of the radiation is absorbed, the other part leaves the sample and is measured by a detector. The characteristic radiation is quantified in the form of k-ratios k = I /I std , where the measured radiation intensity I is normalized by standard intensities I std measured for a known reference sample. A sketch of the physical processes is shown in Figure 1.
We characterize the chemical composition of a material by its mass fraction field c(x). Due to the statistical nature of the k-ratio measurements, the mass fractions c(x), as the solution to the inverse problem of reconstruction, are also probabilistic and all information gained by the experiment is encoded in the statistical moments of c(x) [3]. However, the probabilistic inversion is challenging; thus, this article deals with the minimization of the squared error of k-ratios obtained from an experiment k exp and a forward model k(c(x)). The result can be interpreted as a maximum likelihood estimate under Currently, forward models typically used in reconstruction (ZAF, φ(ρz)) [2,[4][5][6][7] assume either a fully homogeneous or a layered structure interaction volume (the part of the sample in which beam electrons excite atoms). The ZAF model uses correction factors for the atomic number Z(c), the absorption A(c), and the fluorescence F(c) and is implemented as a fixed-point iteration to reconstruct homogeneous mass fractions from k-ratios obtained from a single beam configuration. This idea is extended to layered samples using a correction model based on the depth distribution of X-ray generation φ(ρz). The material is scanned by the electron beam, ensuring that their interaction volumes do not overlap, and each interaction volume is reconstructed independently. In this procedure, the forward models limit the size of the reconstruction pixels to a size that corresponds approximately to the size of the interaction volume.
Previous analysis [8][9][10] of the analytical spatial resolution highlight the need to resolve small structures inside the sample, e.g., the inclusion of spherical droplets in lunar rock which forms due to exposure to energetic particles or small crystals of magnetite embedded in glass. In [8][9][10], the analytical spatial resolution of EPMA is defined as the size of the interaction volume and approaches which physically decrease this volume are described. In contrast, we propose to decouple the spatial resolution and the size of the interaction volume by employing a more sophisticated reconstruction method that harnesses k-ratio measurements from multiple experiments with overlapping interaction volumes.
To analyze the size of the interaction volume of inhomogeneous structures, Monte Carlo models [11][12][13][14][15] are used. Additionally, Monte Carlo models can be used to implement a k-ratio model that allows inhomogeneous material structures by first estimating the electron transport on the basis of Monte Carlo sampling and then approximating the measured X-ray intensity. However, Monte Carlo models suffer from statistical noise which restricts the use of gradient based optimization techniques to solve the inverse problem [16].
In this work, we use a deterministic k-ratio model that bases on the M 1 -model [16][17][18][19][20] in Boltzmann Continuous Slowing-Down approximation (BCSD) and is able to deal with inhomogeneous interaction volumes. The model has previously been validated against Monte Carlo models [16]. Our model additionally allows a fast and precise gradient computation via the adjoint state method, hence efficient gradient based optimization techniques can be applied. The combination of a deterministic model and gradient based optimization methods allows the implementation of a reconstruction method with a higher spatial resolution. Instead of prohibiting overlapping interaction volumes, this method merges sequential k-ratio measurements with multiple beam configurations to accomplish a higher resolution with a pixel size smaller than the size of the interaction volume.

Parameterized Material Description
A sample occupying a domain Ω ⊂ R 3 , which consists of n e chemical elements, is characterized by the mass fractions of all elements c(x; p) : Ω × P → [0, 1] n e at every point x ∈ Ω. The mass fraction describes the ratio of partial m i to total mass m tot in a (infinitesimal small) volume dx around x for one element i = 1, . . . , n e and is therefore constrained by 0 ≤ c(x; p) ≤ 1 and ∑ n e i=1 c i (x; p) = 1. The choice of a global set of parameters p ∈ P ⊆ R n p discretizes the mass fractions and limits the unknown to finite dimensions.
In this work, the mass density ρ(x) is modeled by a linear combination of the mass fractions and the density of the pure elements. Such a model neglects the influence of different molecular structures and assumes that only mass fractions sufficiently describe the material.

k-Ratios as a Function of the Electron Number Density
The measured quantities in EPMA, the radiation intensities of X-rays with wavelengths characteristic to elements occurring in the material, are considered in terms of k-ratios k = I /I std . The measured intensity I gets normalized by the intensity I std measured for a known reference sample, to account for experimental inaccuracies, e.g., the detector efficiency. Furthermore, an experiment is conducted with n s setups (here, we consider different beam energies and beam positions) yielding different measurements. An enumeration of all n k k-ratios allows for collecting all measured k-ratios into k exp ∈ R n k .
We model a k-ratio k j by representing the integral of the zeroth moment of the electron fluence ψ 0 (x, ) over the spatial domain x ∈ Ω and all considered energies ∈ [ 0 , 1 ], weighted by attenuation and emission/fluorescence factors. The presence of element i at a position x is taken into account by the number density of atoms per unit volume N V i = c i (x)ρ(x) /A i , with the atomic mass A i . The ionization cross-section σ emiss i , which combines fluorescence yield and ionization cross-section, describes the fraction of electron-atom collisions leading to X-ray generation. X-ray attenuation is approximated by the Beer-Lambert law with the linear attenuation coefficient being averaged over paths d(x) from each point x to the detector. The attenuation coefficient of a compound affecting an X-ray j is given by µ j (c(y; p)) = ρ(y) ∑ n e m=1 c m (y; p)( µ /ρ) m,j . Emission cross-section and mass attenuation coefficients ( µ /ρ) are interpolated from tabulated data; for details, see [22][23][24][25]. Note that the integrals in Equation (3) model the intensity I j ; hence, the standard intensities I std j can also be calculated this way by inserting the respective mass fractions.

M1-Model of Electron Transport
To model ψ 0 (x, ), we use the M 1 -model [16][17][18][19][20] in BCSD, which describes the first two moments of the electron fluence inside the material. They are defined as (4) using the velocity v( , Θ) of electrons (with magnitude defined by the electron energy and direction Θ) and the stationary number density of electrons n(x, v). Based on the Boltzmann equation in Continuous Slowing-Down approximation, the M 1 -model belongs to a class of models M N which are derived from closing the system of moments by minimization of the Boltzmann entropy. With a closure after the first moment, M 1 , the system takes the form of a nonlinear hyperbolic partial differential equation (PDE) in space and energy, given by is the anisotropy parameter and χ( α ) is the Eddington factor, which implicitly depends on α and is approximated by a rational function. The identity matrix is denoted by I 3 ∈ R 3×3 . Stopping power S and transport coefficient T for a compound are obtained from the mass fractions c, the density ρ, and the respective coefficients of pure elements S i and T i : The analytic forms of S i , T i and the function that approximates χ(||α||) as well as the Jacobian of the flux function F(ψ) are given in Appendix A. For a detailed description of the M 1 -model, we refer to [16][17][18][19][20]. A solution ψ = (ψ (0) , ψ (1) ) T : Ω × [ 0 , 1 ] → R 1+3 of the M 1 -model depends on the material parameters p because stopping power and transport coefficient depend on the mass fractions which are parameterized by p.

Boundary Conditions
Initial 1 and cutoff energy 0 of the computational domain Ω × [ 0 , 1 ] are chosen such that the beam is sufficiently captured beam < 1 and all minimal ionization energies are considered 0 < min i=1...n e edge,i . The edge energy edge,i of an element i is the lowest energy at which ionization still occurs. Initially, no beam electrons are present inside the domain Ω, hence ψ(x, = 1 ) = 0.
The spatial boundary ∂Ω consists of one surface that corresponds to the samples surface S impact where the electron beam impacts. The other boundaries lie inside the sample and are assumed to be sufficiently far from the interaction volume so as not to influence the solution. A schematic of the 2D setup is provided by Figure 2. The electron beam is modeled as Gaussian in space and energy and Dirac in angle, yielding the moments where D beam ∈ R 3 is a unit vector describing the direction of the beam, ϕ the normal distribution, given the respective means, beam and x beam , and variances, σ and Σ beam . This leads to the following spatial boundary condition: which will in the following be referred to as B k (ψ(x, )) = 0 for multiple beam configurations k = 1, . . . n s . For the implementation of the M 1 -model with its boundary conditions using the finite volume library Clawpack [26][27][28], we refer to Section 4.

PDE-Constrained Optimization Problem
Given measured k-ratios k exp , the reconstruction of material parameters can now be written as a PDE-constrained optimization problem. Incorporating the k-ratio model, consisting of the integral equation (Equation (3)) and the M 1 -model G (Equation (5)) with boundary conditions B k , the inverse problem of reconstruction is: This reads as follows. Find the parameters p * such that the objective function H(p), the least-squares error of the residuals between modeled and measured k-ratios, is minimal. Later, we will refer to the residuals as h := k − k exp . Here, the moments of the electron fluence ψ k , which are used to determine the k-ratios k, are subject to the implicit constraint to fulfill the M 1 -model G(ψ k , p) = 0 with its boundary conditions B k (ψ k ) = 0. To obtain a set of k-ratios k, a solution of the PDE-constaints G k is calculated and integrated according to Equation (3).
Note the direct influence of the mass fraction parameters p on the modeled k-ratios k(ψ 1 , . . . , ψ n s , p) via Equation (3) and the indirect influence via the moments of the electron fluence ψ k . Therefore, a perturbation in p will not only change the objective function H due to their direct dependency, but also due to the implied perturbation in ψ k .

Iterative Gradient Based Optimization
To find a solution to the inverse problem, an iterative optimization method of the form p k+1 = p k + ∆p k starting with an initial guess p 0 is used. Many optimization methods calculate the update ∆p based on the gradient of the objective function H(p) and possibly also higher order derivatives-for example, using the negative gradient of the objective in combination with a line search method to find the step length yields the method of steepest descent [29]. The Levenberg-Marquardt [29] algorithm additionally exploits the least-squares form of the objective function and can be interpreted as a combination of a gradient descent approach on the objective function H and a Gauss-Newton method acting on the residuals h. In contrast to the method of steepest descent, which only requires the evaluation of the objective function H and its gradient, the Levenberg-Marquardt method additionally needs the computation of the Jacobian of the residuals J p h(p) with respect to p. In [29,30], these two and several other methods of a similar form are extensively described, with most methods sharing the need to calculate derivatives. Note that such methods only converge to a local optimum and do not guarantee a global minimization.
For the PDE-constrained problem (10), the cost of such methods is dominated by the cost of evaluating the PDE-constraint as part of the objective function and the cost of evaluating its derivatives. It is advantageous to choose an optimization method which requires a low number of evaluations of the objective function and its derivatives, but the efficiency of each evaluation still remains crucial. A naive approach to determine the derivative is a finite difference approximation, which, applied to our problem, requires at least n p + 1 evaluations of the objective function. Each evaluation of the objective function includes solving n s computationally expensive PDEs. Given many parameters n p or experimental setups n s , which are required for high resolution, this approach becomes computationally too expensive. Furthermore, the choice of the spacing for the finite difference approximation is non-trivial. Another option to determine derivatives is the implementation of additional sensitivity equations [31] besides the PDE-constraint. However, the cost also scales with the number of parameters n p in this case, which inhibits their application for our problem.

The Adjoint State Method
The adjoint state method provides an alternative to calculate the gradient. Details on the method can be found in [31][32][33], while here we give the general steps that will later be needed to apply it for problem (10).
The Levenberg-Marquardt method requires the Jacobian of the residuals h, while, for the steepest descent, the gradient of the objective H is needed. Covering both cases, we consider the gradient of a scalar valued functional ∇ p h(ψ, p) with respect to p. The functional h could be one component of the residual h, to gradually calculate the Jacobian for the Levenberg-Marquardt method, or the objective H for the gradient used in steepest descent. Without loss of generality, but, for simpler notation, the discussion here is limited to only one PDE constraint where we consider the state variable ψ, which is implicitly related to the parameters p by G(ψ, p) = 0 (omitting the index k for the beam configurations).
The introduction of a Lagrange multiplier, the adjoint state variable λ, allows for stating the Lagrange function of the optimization problem: Here, ψ, φ = 1 0 Ω ψ(x, ) T φ(x, ) dx d denotes the scalar product in function spaces, ·, · is the scalar product in R n . For G(ψ, p) = 0, one recognizes the equality of L and h, which implies the equality of their gradients ∇ p L = ∇ p h.
Consider the directional derivative dL dq = ∇ p L, q of the Lagrangian L in a direction q. Using ψ q , the direction of change of the state variable corresponding to a small perturba-tion of the parameters, such that G(ψ + τψ q , p + τq) = O(τ 2 ), the directional derivative can be written as (12) Here, δh /δψ denotes the functional derivative, δh /δp the vector of partial derivatives with respect to p, not including the variation due to ψ p , and [ δG /δ·] the corresponding Fréchet derivative operators. Multiple calculations of Equation (12), while iterating over the unit vectors as direction q, accumulate to the gradient ∇ p h, but would include multiple expensive calculations of the direction ψ q . To avoid this calculation, Equation (12) is rearranged to such that the first scalar product in Equation (13) including ψ q vanishes and only the second scalar product remains as the directional derivative of L. Then, the gradient ∇ p h is identified as To calculate the full gradient, one evaluation of the adjoint state Equation (14) is sufficient because it does not depend on the direction q. Its structure is usually similar to the structure of the forward equation G(ψ, p) = 0 but depends on the specific form of the adjoint operator. Being an essential component of the method, the adjoint operator for our problem (Equation (10)) will now be derived under consideration of the form of the PDE-constraint G (Equation (5)).

Adjoint State Method for the M1-Model Constraint
The adjoint state method bases on the identity The adjoint operator [ δG /δψ] * allows for isolating the direction ψ q in the scalar product of Equation (13), which subsequently enables stating the adjoint state Equation (14). To derive the adjoint state equation for the M 1 -model, this identity must be considered, taking G from Equation (5) into account. The Fréchet derivative of G can be found using the Jacobians of the flux functions J ψ F d (ψ) =: A d (ψ). Inserting into Equation (16) and proceeding using integration by parts yields The boundary integrals, which are enforced to be zero, yield boundary conditions for the adjoint state variable λ. At energy 1 , the state variable ψ is specified by the boundary condition of the M 1 -model which does not depend on the parameters p; therefore, ψ q ( 1 ) = 0. At energy 0 , we have to prescribe the adjoint state variable λ( 0 ) = 0. In addition, the boundary term for the spatial boundary ∂Ω with outer normal vector n has to vanish, which yields boundary values for λ on ∂Ω. Their numerical implementation is discussed along with the boundary conditions for the forward model in Section 4.
Then, the adjoint state equation to be solved for the adjoint state variable λ is given by To solve the adjoint equation, the solution ψ of the forward equation G(ψ, p) = 0 has to be known, since A T d (ψ) and δh(ψ,p) /δψ depend on it. The term δh /δψ acts as a source to the adjoint state variable and depends on the choice of the functional h to be differentiated. For the Levenberg-Marquardt algorithm, the calculation of the gradient of each residual is required, thus we choose h = h j = k j (ψ, p) − k exp j . Since in Equation (3) the first component of the state variable ψ (0) enters the integral linearly, the functional derivative of one residuum h j is given by Here, ξ j collects all weighting factors of the integral in Equation (3).
To implement the steepest descent method, the gradient of the objective H is required; therefore, h = H = h 2 2 and, by a similar argument as above, its functional derivative is δH δψ Having the adjoint equation determined by the adjoint operator and the respective source for the adjoint state variable, the calculation of the gradient via the adjoint state method can be summarized by 1.
the solution of the adjoint Equation (18) for λ and 3. the calculation of the gradient ∇ p h(ψ, p) = δ δp (h(ψ, p) + λ, G(ψ, p) ) In comparison to a finite difference approximation with n p + 1 solutions of a PDE, the adjoint state method only requires the solving of the forward (first step) and the adjoint (second step) equation to calculate the gradient. The calculation of the gradient (third step) can usually be implemented efficiently. Consequently, the adjoint state method offers an efficient way to determine the gradient to be used in an iterative optimization method for our problem (10).

Material Parameterization
The subdivision of the material into reconstruction pixels (see Figure 1) motivates the discretization of the mass fractions as piecewise constant functions. Each parameter describes one of the constant mass fractions in one of the reconstruction pixels. Note that, due to the conditions on c (c.f. Section 2), n e − 1 parameters per pixel are sufficient.

Solving the M1-Model Using Clawpack
The form of Equation (5) motivates the application of finite volume methods to approximate its solution. Our solver of the M 1 -model in two spatial dimensions is implemented using the finite volume library Clawpack/pyclaw [27]. A pseudo time variable t( ) = 1 − transforms the energy interval to be compatible with Clawpack's formulation. The term −∂ (Sψ) can then be split into −ψ∂ S and S∂ t ψ, where the first acts like a source term and the second fits to Clawpack's formulation with S being the capacity function. The capacity coefficient and the source term are applied to the solution via the 'before-step' and the 'source' function defined by Clawpack. Furthermore, it is important to choose the grid of finite volume cells at least as fine as the reconstruction pixels and to assure that their cell/pixel boundaries coincide.
To implement a Riemann solver, we, as suggested in [26], approximate the Jacobian A ≈ J ψ F d ( ψ L +ψ R 2 ) at each cell interface using the average of the state variable ψ. (See Appendix A.4 for the analytical form of the Jacobian of the flux function J ψ F 1 (ψ).) Then, A is numerically decomposed in eigenvalues and eigenvectors A = RΛR −1 on the basis of which the wave speeds Λ, the wave strengths R −1 ∆ψ = R −1 (ψ R − ψ L ) and right and left going fluctuations apdq = R max{Λ, 0}R −1 ∆ψ and amdq = R min{Λ, 0}R −1 ∆ψ can be identified.
Spatial boundary conditions on ∂Ω are implemented in Clawpack using ghost cells surrounding the computational domain. At the boundary where the beam enters the domain, the ghost cells are updated using Equation (8). The remaining ghost cells are set to zero, justified by the assumption that the computational domain is large enough such that no electrons reach those boundaries.

Implementation of the k-Ratio Model
Clawpack discretizes the electron fluence into 'snapshots' of the solution at specific energies with piecewise constant cell values in space. The spatial integral for the k-ratios k in Equation (3) reduces to a weighted summation over all cells. The attenuation coefficients for each cell is approximated using the path from its midpoint to the detector, which is split in segments overlapping the reconstruction pixels. For the integral in energy, the trapezoidal rule is applied.

Adaptions to Solve the Adjoint State Equation and Compute the Gradient
Analogous to the M 1 -model, the adjoint Equation (18) is also implemented using Clawpack/pyclaw taking into consideration the reversed direction of computation. Using another pseudo time variable for the adjoint τ( ) = − 0 , the term S∂ λ can be written as S∂ τ λ so that it matches Clawpack's interface. For the wavespeeds, wavestrengths and fluctuations, the same Riemann solver as for the forward problem is employed. It is only modified to calculate the waves of the transposed Jacobians A T d (ψ) instead of the Jacobians A d (ψ). To evaluate the Jacobians, the whole forward solution is stored in memory. Via Clawpack's 'source' function, the source term of the adjoint state variable, depending on the required gradient, is integrated. At energy 0 , we choose λ( 0 ) = 0 as the initial condition, while the ghost cells of the spatial boundary are set to Clawpack's 'extrapolate' setting. The partial derivatives in Equation (15) are computed using the algorithmic differentiation tool jax [34].

Numerical Experiments
The following numerical examples illustrate the application of the presented methods and demonstrate their capabilities and limitations. We use synthetic measurements in our experiments which are computed using the same forward model as used for the inversion. The utilized optimization method does not guarantee convergence to a global minimum, and no constraints are enforced on c. Overcoming those limitations and implementing a more robust reconstruction method is left for future work. Bombarding this sample with an electron beam from above in a negative y-direction (beam energy beam = 12 keV, σ = 0.1 keV, beam position x beam = (500, 0) T nm, σ x = 30 nm), produces the zeroth moment of the electron fluence ψ 0 shown in Figure 3. It displays three snapshots of ψ 0 in energy space as a heat map over the domain Ω. Figure 3 also illustrates the interaction volume, the volume where ψ 0 > 0, which is larger than the reconstruction pixels p 1,... 4 .
While the number of electrons with high energy is highest near the point of electron impact, electrons with lower energy are more widely distributed in the sample. The contribution of each pixel to the k-ratios scales linearly with the zeroth moment of the electron fluence (Equation (3)), hence a measurement carries most information from pixels close to the surface near the point of electron impact. It is impossible to gain information about the chemical composition in pixels where ψ 0 is zero, so only pixels lying sufficiently close to the surface can be reconstructed.

Convergence to a Particular Solution
In order to show convergence of the M 1 -model solver and the k-ratio model to a particular solution, we investigate the difference of a solution ψ 0 calculated with different resolutions to a reference solution ψ 0 * calculated with high resolution (320 × 240, cell size: 3.125 nm × 3.3 nm). Additionally, the k-ratios k based on the variable resolution ψ 0 are compared to the k-ratios k * computed with the high resolution ψ 0 * . Figure 4 shows the convergence to the reference solution with decreasing cell sizes. The different resolutions as well as other model parameters are shown in Table 1.

Synthetic Measurements
To mimic a reconstruction problem, we generate synthetic k-ratios (Cu Kα and Mn Kα) from the six different beam configurations given by the beam positions and parameters in Table 2. This yields a set of 12 k-ratios, which we consider as synthetic data k exp . All other model parameters are taken from the previous example (Table 1).  In order to give an insight into the behavior of the inverse problem of reconstruction, we will evaluate the objective function k(ψ, p) − k exp 2 2 . All k-ratios, synthetic measurements and modeled k-ratios are computed using the same discretization (spatial grid: 80 × 60, energy steps: 100) of the M 1 -model. Both plots in Figure 5 visualize the dependency of the objective function on two of the four parameters describing the mass fractions in the reconstruction pixels. The parameters p 1 and p 2 describe reconstruction pixels that are next to each other while the pixel belonging to p 3 is above that of p 2 , cf. the marked cells in Figure 3. The parameters in all other reconstruction pixels are kept to the values used to calculate the synthetic data. The crosses mark the hidden truth, i.e, the parameter values used to calculate the synthetic data. Both contour plots show a unique minimum with a convex objective function, which, however, cannot be expected for the case where all parameters are simultaneously variable. Nevertheless, Figure 5 gives an intuition into the nature of the problem. The shape of the objective function, which increases less rapidly in the direction in which the parameters cancel each other out, e.g., increasing p a while decreasing p b , shows that the pixels can partially compensate for each other. It is not detectable in which pixel the x-rays were generated, so an increase in the mass fraction of an element in one pixel may be counteracted to some extent by a decrease in its mass fraction in an adjacent pixel. Therefore, a high spatial resolution can only be achieved by combining measurements obtained from experiments with different electron beam positions and energies.
In addition, the greater elongation in the right plot indicates that it is more difficult to resolve the vertical direction than the lateral direction. The intensity of detected X-rays generated deeper inside the material is lower due to the decreased electron energy and the increased absorption, hence the sensitivity of the k-ratios to changes in the mass fractions of deeper pixels decreases. Using different electron beam energies, the depth of pixels which get excited can be varied and better discrimination of vertical pixels can be achieved. However, the reconstruction of vertical pixels is more difficult than the reconstruction of lateral pixels because, in lateral reconstruction, different electron beam positions can be used to excite only certain pixels. In vertical reconstruction, higher lying pixels are always excited as well. In the left-hand plot, two parameters describing pixels lying next to each other are varied, and, in the right-hand plot, two parameters lying below each other are varied, see Figure 3 for the scenario, pixel size, and position. In both cases, all other parameters are fixed to their hidden truth values. The cross marks the hidden truth of the variable parameters.

Reconstruction of Four Parameters
Using the Levenberg-Marquardt optimization algorithm to minimize the objective function of the artificial reconstruction problem yields the optimum illustrated in the left-hand plot of Figure 6, together with the initial values and their evolution during the optimization. Synthetic measurements are computed using the following discretization of the M 1 -model: 160 × 120 spatial grid, 300 energy steps. For the reconstruction, multiple discretizations are used: a 40 × 30 spatial grid with 70 energy steps, a 80 × 60 spatial grid with 100 energy steps and a 160 × 120 spatial grid with 300 energy steps.
The hidden truth values are marked with solid lines, which after around 10 steps are identified by the optimization algorithm. The dotted, dashed, and dashdotted lines show the evolution of the parameters during the optimization for the different discretizations. Using the same resolution as used for the synthetic measurements (160 × 80, dashdotted line) yields a perfect reconstruction. For other resolutions, the reconstructed parameters differ from the hidden truth but converge to similar values. This is due to the fact that the k-ratios also differ when being calculated from the same parameters but with different discretizations. The right-hand plot shows the value of the objective function over the steps of the Levenberg-Marquardt iterations. Again, deviations are visible for different discretizations.
The size of the reconstruction pixel is smaller than the size of the interaction volume (c.f. Figure 3); nevertheless, the iterative reconstruction method is able to reveal the mass fractions. The combination of measurements with different electron beam configurations and overlapping interaction volumes allows the reconstruction with a high spatial resolution.

Reconstruction on Small Vertical Layers
We now consider the reconstruction of a material consisting of iron, nickel, and chromium arranged in vertical layers with a size of 50 nm, which is smaller than the size of the interaction volume. The material is excited by beams with an energy of 17 keV and positions centered on each of the vertical layers. All model parameters and discretizations used in this example are given in Table 3. In Figure 7, the iterative reconstruction of the mass fractions is illustrated. The hidden truth configuration of all layers is shown as the background color, where blue is chromium, red is nickel, and the remaining orange part is iron. The Levenberg-Marquardt iterations for the mass fractions of chromium and nickel in each layer are shown by the blue solid and red dotted lines. The mass fraction of iron can be deduced (c.f. Section 2). From the initial configuration (the leftmost point in each layer), all parameters are reconstructed successfully. The rightmost point in each layer shows the mass fractions after six iterations. The algorithm was able to identify the mass fractions in all vertical layers.

Conclusions and Outlook
The results presented in this paper illustrate the potential of using high-resolution deterministic models (e.g., the M 1 -model in BCSD) for the implementation of a reconstruction method that improves the spatial resolution of EPMA. Although a compromise must always be found between higher resolution and the number of measured values or the computing effort, the size of the reconstruction pixels is no longer limited by the size of the interaction volume and can be selected independently.
The reconstruction can be formalized as a deterministic PDE-constrained optimization problem and can be implemented using iterative gradient-based minimization techniques. Hereby, the computation of the derivative of the objective function forms the major part of the calculation of the next iterate. Hence, the adjoint state method, which efficiently calculates the gradient, yields significant reductions in computing time over e.g., the finite difference method.
Having the possibility to resolve fine structures in materials using EPMA raises the question about the confidence in a reconstruction. By means of the numerical examples, it was illustrated that reconstruction pixels close to the surface are easier to reconstruct than pixels deeper inside the material, and other pixels are not reconstructible. A high number of parameters, for example in a 3D reconstruction, increases the problem of underdetermination and raises the need to include prior knowledge about the material, e.g., in the form of additional regularization. Furthermore, an uncertainty quantification of the reconstruction result based on the uncertainty of the model and the measurements would provide deeper insight. In any case, the noise of real experimental data must be considered since the strength of the noise influences the accuracy of the reconstruction. Reconstruction experiments based on real measurements should be accompanied by a careful parameter study, since the uncertainty for many of the parameters used in this model is unclear. Nevertheless, the examples presented in this paper indicate the potential of using deterministic k-ratio models to improve the spatial resolution of EPMA compared to conventional methods by implementing a more sophisticated reconstruction method.  In the following, we will derive an analytical expression of the Jacobian of the flux function with respect to ψ = (ψ (0) , (ψ (1) ) T ) T = (ψ (0) , ψ (1) x , ψ (1) y , ψ (1) z ) T . Exemplarily, we will only consider the flux function into direction of the first unit vector F 1 (ψ) = F(ψ)e 1 .