An Adjoint Method for High-Resolution EPMA Based on the Spherical Harmonics (PN) Model of Electron Transport

Abstract

Electron Probe Microanalysis (EPMA) is a nondestructive technique to determine the chemical composition of material samples in the micro- to nanometer range. Based on intensity measurements of x-radiation, the reconstruction of the material composition and structure poses an inverse problem. The reconstruction methods currently applied are based on models that assume a homogeneous or a layered structure of the studied material. To increase the spatial resolution of reconstruction in EPMA the combination of a more sophisticated reconstruction method, which is based on a model that allows complex material structure, together with multiple measurements with varying beam configurations is required. We present a deterministic k-ratio model that is based on the P_N model, an approximation of the radiative transfer equation for electron transport. Our goal is to approximate a maximum likelihood solution of the inverse problem using gradient-based optimization. We detail the application of the model in the context of algorithmic differentiation, in particular by deriving its continuous adjoint formulation. Algorithmic differentiation provides the flexibility to adapt the reconstruction method to various material parametrizations and thus to regularize and take into account prior knowledge. Through examples, we verify our implementation and demonstrate the flexibility of the reconstruction/differentiation method.

Keywords:

• EPMA
• electron probe microanalysis
• spherical harmonics moment approximation
• PN
• algorithmic differentiation
• adjoint method, inverse modeling

1. Introduction

Electron Probe Microanalysis (EPMA) (Heinrich and Newbury 1991; Reimer 1998) is an imaging technique used for the quantitative analysis of the composition of solid material samples at the micro- to nanometer scale. The sample is excited by a focused beam of electrons which induces multiple relaxation processes inside the sample. In EPMA, the emission of characteristic x-rays is of special focus. If an electron that is induced by the beam strikes a bound electron that occupies an atomic shell of an atom inside the specimen, the bound electron is ejected from its shell and the atom is left with a vacancy. Outer shell electrons fill this vacancy by emitting a quantized x-ray with an energy corresponding to the energy level difference of the excited (initial) and the relaxed (final) level. The energy levels of electron shells are characteristic of a specific atom, hence the energy of the emitted x-ray provides information about the composition of the material sample. Reconstructing the chemical composition from intensities of emitted x-radiation (normalized into k-ratios ${(k_{\alpha })}_{\alpha \in \mathbb{A}},\mathbb{A}:$ measurement setups) forms the inverse problem in EPMA. In Figure 1, we outline the main physical processes that occur during an experiment.

Figure 1. A sketch of the physical processes in EPMA. The sample G is rastered by electron beams (blue) with different positions μ_x. Electrons e^– from the electron beam (blue dotted line) scatter inside the sample and strike bound electrons that leave a vacancy (blue circle). Outer shell electrons (blue discs) fill this vacancy and release an x-ray of characteristic energy (red wobbly line). The x-ray travels through the sample and is counted by a detector. Beam electrons only excite a certain volume of the sample, the interaction volume (gray ellipses) which scales with the beam energy ${\mu }_{ϵ}.$

The crucial ingredient to the inverse problem of reconstruction is the definition of a model $k \circ \mathcal{P}:\mathbb{P}\to {\mathbb{R}}^{\left|\mathbb{A}\right|}$ that maps parameters $p\in \mathbb{P}\subseteq {\mathbb{R}}^{n_p}$ ($\mathbb{P}:$ parameter space, n_p: number of parameters) to k-ratios $(k \circ \mathcal{P})\left(p\right)\in {\mathbb{R}}^{\left|\mathbb{A}\right|}$ ($\left|\mathbb{A}\right|:$ number of measurements). Then the reconstruction can be formalized as the following inverse problem: Find the set of parameters $p^{*}$ such that the modeled k-ratios $(k \circ \mathcal{P})\left(p^{*}\right)$ reproduce the experimental k-ratios $k^{\hspace{0.17em}\exp \hspace{0.17em}}$ as good as possible.

The applications of EPMA are manifold (Llovet et al. 2021; Moy, Fournelle, and von der Handt 2019; Pinard et al. 2013). From geology, and material science to electronics, researchers rely on EPMA to determine the material composition and structure of samples. Every application brings with it different prior knowledge and objectives that must be considered in the reconstruction. We subdivide our model into a k-ratio model $k:(G\to {\mathbb{R}}^{n_e})\to {\mathbb{R}}^{\left|\mathbb{A}\right|},\rho \hspace{0.17em}\mapsto \hspace{0.17em}k[\rho ]$ and a material parametrization $\mathcal{P}:\mathbb{P}\to (G\to {\mathbb{R}}^{n_e}),p\hspace{0.17em}\mapsto \hspace{0.17em}\mathcal{P}(p)=\rho .$ Thereby ${\rho }_i(x)$ is the mass concentration of element i in a compound with n_e constituents at a point $x\in G.$ The sample domain $G\subset {\mathbb{R}}^3$ will be specified later. While the k-ratio model k can handle general materials, the reconstruction can be tailored to specific applications via the material parametrization $\mathcal{P}$ (e.g., using prior knowledge). In particular, our definition includes reconstruction approaches currently applied in EPMA, which are however very restrictive on the material parametrizations $\mathcal{P}.$ ZAF-models assume homogeneity of the material inside the interaction volume of each beam. Effectively, this assumption limits the spatial resolution of EPMA to the size of the interaction volume (Buse and Kearns 2020; Carpenter and Jolliff 2015; Moy and Fournelle 2017). Reconstruction methods based on the integration of $\varphi (\rho z)$-curves currently only allow for samples that are layered in depth (Moy and Fournelle 2020). We aim to leverage a more sophisticated k-ratio model and the combination of k-ratio measurements from various beam positions and energies (resp. measurement setups $\mathbb{A}$) to obtain a high-resolution spatial image of the material composition.

In this article, we present our model for k-ratios (taken from Bünger, Richter, and Torrilhon 2022) with a detail on the dependence on the mass concentrations ρ). The main part is dedicated to the derivation of an efficient and modular approach to compute the derivatives of our model, facilitating its integration into differentiable programming (Dorigo et al. 2023). We derive our differentiation approach based on the Adjoint Mode of Algorithmic Differentiation (Griewank 2003; Naumann 2011) and the Continuous Adjoint Method (Plessix 2006). The derivative forms the crucial part in the implementation of gradient-based optimization methods that find a local minimum of an objective function as a candidate solution of the inverse problem. We show the agreement of our gradient computation with finite differences and show the advantages of the adjoint method in terms of computation time. We also provide reconstruction experiments that compare different parametrizations in 1D and show a 2D reconstruction example.

1.1. The inverse problem of reconstruction

In theory (Stuart 2010; Tarantola 2005), the solution of the inverse problems is defined as knowledge of the posterior uncertainty of the parameters based on the likelihood of the data, which describes the uncertainty of the model and the measurement, and the prior uncertainty of the parameters. So far, we consider a Bayesian inversion of the presented problem as intractable. However, the statistical understanding interprets our approach as finding the maximum likelihood estimate under the assumption of an isotropic Gaussian likelihood. We define the reconstruction problem as the following optimization problem:(1) $$p^{*}=\underset{p\in \mathbb{P}}{\text{arg} \min }\underset{=(J\hspace{0.17em}\circ \hspace{0.17em}k\hspace{0.17em}\circ \hspace{0.17em}\mathcal{P})(p)}{\underbrace{{\sum }_{\alpha \in \mathbb{A}}((k_{\alpha } \circ \mathcal{P})(p)-k_{\alpha }^{\hspace{0.17em}\exp \hspace{0.17em}}{)}^2}}$$(1)

Assuming that the objective function $J\hspace{0.17em}\circ \hspace{0.17em}k\hspace{0.17em}\circ \hspace{0.17em}\mathcal{P}$ is differentiable w.r.t. p, gradient-based optimization methods (Steepest Descent, Conjugate Gradient, L-BFGS) (Nocedal and Wright 2006) can be utilized to iteratively approximate a local minimum $p^{*}\in \mathbb{P}.$ For an efficient application of gradient-based methods, the efficient computation of the gradient is crucial. Additionally, we require our code to remain modular, such that either different parts of the model and, more importantly, the parametrization of the mass concentrations $\mathcal{P}$ is exchangeable. We remark that the computation of maximum likelihood estimates can be ambiguous and ill-posed depending on the number of available measurements and the chosen parametrization. Additional prior information about the material would be necessary that could enter our minimization problem (1)(1) $$p^{*}=\underset{p\in \mathbb{P}}{\text{arg} \min }\underset{=(J\hspace{0.17em}\circ \hspace{0.17em}k\hspace{0.17em}\circ \hspace{0.17em}\mathcal{P})(p)}{\underbrace{{\sum }_{\alpha \in \mathbb{A}}((k_{\alpha } \circ \mathcal{P})(p)-k_{\alpha }^{\hspace{0.17em}\exp \hspace{0.17em}}{)}^2}}$$(1) as an additional regularization term $\mathcal{R}(\rho (·),p).$ Regarding the ill-posedness of the inverse problem, we refer to further research and note that the addition of a differentiable regularization term is possible.

2. The forward model

The characteristic x-ray intensity model used in this article combines a model for electron transport based on the continuous slowing down approximation (CSD) of the linear Boltzmann equation (radiative transfer) with a subsequently applied model for x-ray generation and attenuation. In Bünger (2021); Bünger, Sarna, and Torrilhon (2022); Bünger, Richter, and Torrilhon (2022), the model is derived in more detail. Here we derive and summarize the main equations and address the dependence on mass concentrations ρ in detail.

2.1. K-ratios

In EPMA, x-ray intensity is measured by counting all characteristic x-rays emitted by electron scattering induced by the electron beam. Because of multiple uncertain quantities concerning the experimental device, the x-ray intensity $I_{\alpha }$ is normalized into a k-ratio $k_{\alpha }$ using standard intensities $I_{\alpha }^{\text{std}}$ measured from a known reference sample.(2) $$k_{\alpha }=\frac{I_{\alpha }}{I_{\alpha }^{\text{std}}}$$(2)

Therefore, the normalization eliminates uncertain multiplicative factors that influence the intensity, e.g., the detector efficiency.

An experiment yields multiple k-ratios. Each k-ratio corresponds to a specific x-ray transition and a specific experimental setup, e.g., the beam position, energy or angle. To distinguish between k-ratios, we use $\alpha =({\alpha }_L,{\alpha }_B)\in \mathbb{A}$ to denote a multi-index, that contains all information about the transition α_L and the beam setup α_B.(3) $$\begin{array}{cc}\alpha =\hfill & (\underset{={\alpha }_L}{\underbrace{Z_i\in \{Z_1,Z_2,\dots ,Z_{n_e}\},I\in \{K,L1,\dots ,M1,\dots \},F\in \{L1,\dots ,M1,\dots \},}}\hfill \\ \hfill & \underset{={\alpha }_B}{\underbrace{({\mu }_{x,\text{beam}},{\Sigma }_{x,\text{beam}})\in {\mathbb{R}}^2\times {\mathbb{R}}^{2\times 2},({\mu }_{ϵ,\text{beam}},{\Sigma }_{ϵ,\text{beam}})\in \mathbb{R}\times \mathbb{R},({\mu }_{\Omega ,\text{beam}},{\kappa }_{\Omega ,\text{beam}})\in S^2\times \mathbb{R}}})\hfill \end{array}$$(3)

The individual symbols are the atomic number $Z_i\in {\alpha }_L,$ the initial (excited) $I\in {\alpha }_L$ and the final (relaxed) $F\in {\alpha }_L$ x-ray level of the electron-transition,¹ the mean and variance of the beam’s position $({\mu }_{x,\text{beam}},{\Sigma }_{x,\text{beam}})\in {\alpha }_B,$ of the beam’s energy $({\mu }_{ϵ,\text{beam}},{\Sigma }_{ϵ,\text{beam}})\in {\alpha }_B$ and the mean and concentration parameter of the beam’s direction $({\mu }_{\Omega ,\text{beam}},{\kappa }_{\Omega ,\text{beam}})\in {\alpha }_B.$

2.2. Characteristic X-ray emission model

Based on mass concentrations ρ, we aim for a generic model for k-ratios $k_{\alpha }$ which, according to the usual definition in EPMA, are computed as the ratio of a measured $I_{\alpha }$ and a standard intensity $I_{\alpha }^{\text{std}}.$(4) $$k_{\alpha }\left[\rho \right]=\frac{I_{\alpha }\left[\rho \right]}{I_{\alpha }\left[{\rho }_{\alpha }^{\text{std}}\right]}$$(4)

Both intensities can be computed from the same model, only the underlying mass concentrations, ρ resp. ${\rho }_{\alpha }^{\text{std}},$ changes.

Neglecting multiplicative factors that cancel due to the normalization with standard intensities $I_{\alpha }^{\text{std}},$ we formulate the intensity of characteristic x-rays $I_{\alpha }$ as the integral of the attenuation field ${\mathcal{A}}_{\alpha },$ the number of atoms ${\mathcal{N}}_{\alpha }$ and the generation field of characteristic x-rays ${\mathcal{I}}_{\alpha }$ over a domain $G\subset {\mathbb{R}}^3$ of the sample which is chosen such that G covers the support of ${\mathcal{I}}_{\alpha },$ i.e., ${\mathcal{I}}_{\alpha }(x)=0 \forall x\in {\mathbb{R}}^3\setminus G.$(5) $$I_{\alpha }\left[\rho \right]={\int }_G{\mathcal{A}}_{\alpha }\left[\rho \right](x) {\mathcal{N}}_{\alpha }\left[\rho \right](x) {\mathcal{I}}_{\alpha }\left[\rho \right](x)\text{d}x$$(5)

The attenuation field ${\mathcal{A}}_{\alpha }$ describes the intensity reduction of generated x-rays while traveling through the sample toward the detector. Attenuation is quantified by Beer-Lamberts law where the attenuation coefficient ${\mu }_{\alpha }$ is the weighted sum of mass concentrations ρ_i and mass attenuation coefficients ${\tau }_{\alpha ,i}.$²(6) $${\mathcal{A}}_{\alpha }\left[\rho \right](x)=\hspace{0.17em}\exp \hspace{0.17em}(-{\int }_{l(x,x_d)}{\mu }_{\alpha }(y)\mathrm{}\text{d}y)=\hspace{0.17em}\exp \hspace{0.17em}(-{\int }_{l(x,x_d)}\sum _{i=1}^{n_e}{\tau }_{\alpha ,i}{\rho }_i(y)\mathrm{}\text{d}y)$$(6)

In (6)(6) $${\mathcal{A}}_{\alpha }\left[\rho \right](x)=\hspace{0.17em}\exp \hspace{0.17em}(-{\int }_{l(x,x_d)}{\mu }_{\alpha }(y)\mathrm{}\text{d}y)=\hspace{0.17em}\exp \hspace{0.17em}(-{\int }_{l(x,x_d)}\sum _{i=1}^{n_e}{\tau }_{\alpha ,i}{\rho }_i(y)\mathrm{}\text{d}y)$$(6) the integration domain is the straight line $l(x,x_d)\subset G$ from a point $x\in G$ toward the detector position $x_d\in \partial G,$ such that the line $l(x,x_d)$ lies inside the domain G. Neglecting the part outside the domain G assumes that x-rays travel unattenuated in the vacuum surrounding the material sample.

The number of atoms ${\mathcal{N}}_{\alpha }$ (per unit volume) is given by the division of the mass concentration field ρ_i by the atomic mass $A_{\alpha }$ of the element $Z_i\in \alpha .$(7) $${\mathcal{N}}_{\alpha }\left[\rho \right](x)=\frac{{{\rho }_i(x)|}_{Z_i\in \alpha }}{A_{\alpha }}$$(7)

The third factor that influences the x-ray intensities $I_{\alpha }$ is the field of characteristic x-ray generation.(8) $${\mathcal{I}}_{\alpha }\left[\rho \right](x)={\int }_E{\sigma }_{\alpha }(ϵ){\int }_{S^2}\underset{={\psi }_{\alpha }(x,ϵ,\Omega )}{\underbrace{|v(ϵ)|f(x,ϵ,\Omega )}}\mathrm{}\text{d}\Omega \mathrm{}\text{d}ϵ$$(8)

It contains the electron fluence ${\psi }_{\alpha }(x,ϵ,\Omega )=|v(ϵ)|f(x,ϵ,\Omega )$ ($v(ϵ):$ electron velocity, $f(x,ϵ,\Omega ):$ number density of electrons) of beam electrons at $x\in G$ with energy $ϵ\in E=[{ϵ}_{\text{cut}},{ϵ}_{\max }]$ into direction $\Omega \in S^2$ inside the material probe which is weighted by the x-ray production cross-section ${\sigma }_{\alpha }.$³ The energy interval E is defined by the cutoff energy ${ϵ}_{\text{cut}},$ which is chosen so small that no more ionization can occur, and ${ϵ}_{\max },$ which is chosen sufficiently large to capture all beam electrons.

Describing the electron fluence ${\psi }_{\alpha }$ forms the most difficult part of modeling the emission of characteristic x-rays for inhomogeneous samples. We detail our model for the electron fluence ${\psi }_{\alpha }$ in the following section.

2.3. Deterministic electron transport model

The dynamics of the beam electrons are described by a linear Boltzmann equation, an integro-differential equation describing the evolution of electron number density in phase-space (position-velocity space), whereby the physical processes governing the evolution of the electron number density are free-flight phases interrupted by elastic and inelastic collisions with atoms of the material probe. Since the timescale of the physical processes is very small in comparison to the duration of the experimental measurements of x-ray intensities, it is valid to neglect the time dependence and assume a steady-state distribution of beam electrons immediately after switching on the electron beam. Furthermore, as inelastic scattering cross-sections are highly peaked around small energy losses, i.e., electrons most probably lose a significant amount of energy in a sequence of small energy losses, it is common in EPMA to model the energy loss of projectile electrons as a continuous deceleration along the trajectory. These simplifications lead to modeling the beam electrons by the Boltzmann equation of electron transport in continuous slowing down (CSD, resp. BCSD) approximation – an evolution equation in energy space given by (Larsen et al. 1997)(9) $$-{\partial }_{ϵ}(S[\rho ](x,ϵ){\psi }_{\alpha }(x,ϵ,\Omega ))+\Omega ·{\nabla }_x{\psi }_{\alpha }(x,ϵ,\Omega )=(\mathcal{Q}[\rho ]{\psi }_{\alpha })(x,ϵ,\Omega ),$$(9) where the energy loss is governed by the stopping power S and the scattering operator $\mathcal{Q}$ accounts for elastic and inelastic collisions (CSD-average). Both, the stopping power S and the scattering operator $\mathcal{Q}$ are derived from cross-sections for compounds, which are defined by means of the additivity approximation. We use the cross-sections from Salvat et al. (2007) that are also used in Bünger (2021); Olbrant (2012). Assuming a well-defined atomic mass, we can write both as the sum of atomic quantities S_i and ${\mathcal{Q}}_i$⁴ weighted by mass concentrations ρ_i.(10) $$S\left[\rho \right](x,ϵ)=\sum _{i=1}^{n_e}{\rho }_i(x)S_i(ϵ)$$(10) (11) $$(\mathcal{Q}[\rho ]{\psi }_{\alpha })(x,ϵ,\Omega )=\sum _{i=1}^{n_e}{\rho }_i(x)({\mathcal{Q}}_i{\psi }_{\alpha })(x,ϵ,\Omega )$$(11)

The BCSD model describes the evolution of the electron fluence ${\psi }_{\alpha }$ in energy space and consequently needs to be equipped with initial and boundary conditions. As the evolution is prescribed from high to low energies an initial condition defines the electron fluence at the maximal energy ${ϵ}_{\max }$(12) $${\psi }_{\alpha }(x,ϵ={ϵ}_{\max },\Omega )={\psi }_{\alpha ,0}(x,\Omega )=0.$$(12)

We will initialize the electron fluence at an energy ${ϵ}_{\max }$ sufficiently larger than the mean beam energy ${\mu }_{ϵ,\text{beam}}$ by ${\psi }_{\alpha ,0}(x,\Omega )=0$ and model the electron beam by boundary conditions. The electron transport model can be equipped with Dirichlet-type boundary conditions by prescribing the incoming half of the electron fluence at the boundary(13) $${\psi }_{\alpha }(x,ϵ,\Omega )={\psi }_{\alpha ,\text{in}}(x,ϵ,\Omega )\hspace{1em}\forall \Omega \in S^2:n·\Omega <0$$(13) where n denotes the outward-pointing normal-vector at $x\in \partial G$ and ${\psi }_{\alpha ,\text{in}}$ is a given distribution of incoming particles. We model electron beams by a product of Gaussian distributions in space and energy together with a von-Mises-Fisher distribution in direction. An electron beam that is aligned with the outward-pointing normal n of the polished material surface and focused onto the point ${\mu }_{x,\text{beam}}\in \partial G$ is modeled by(14) $${\psi }_{\alpha ,\text{in}}(x,ϵ,\Omega )={\mathfrak{N}}_{\alpha ,\mathrm{x}}(x){\mathfrak{N}}_{\alpha ,ϵ}(ϵ){\mathfrak{F}}_{\alpha ,\Omega }(\Omega ),$$(14) where ${\mathfrak{N}}_{\alpha ,x}(x)$ is the probability density function of a Gaussian with mean ${\mu }_{x,\text{beam}}\in \alpha$ and covariance ${\Sigma }_{x,\text{beam}}\in \alpha .$ Respectively the distribution in energy ${\mathfrak{N}}_{\alpha ,ϵ}(ϵ)$ (Gaussian with mean ${\mu }_{ϵ,\text{beam}}\in \alpha$ and variance ${\Sigma }_{ϵ,\text{beam}}\in \alpha$) and the distribution in direction ${\mathfrak{F}}_{\alpha ,\Omega }(\Omega )$ (von-Mises-Fisher with mean ${\mu }_{\Omega ,\text{beam}}\in \alpha$ and concentration ${\kappa }_{\Omega ,\text{beam}}\in \alpha$) is defined.

2.4. Spherical harmonic approximation

A numerical solution of the electron transport model in Equation (9)(9) $$-{\partial }_{ϵ}(S[\rho ](x,ϵ){\psi }_{\alpha }(x,ϵ,\Omega ))+\Omega ·{\nabla }_x{\psi }_{\alpha }(x,ϵ,\Omega )=(\mathcal{Q}[\rho ]{\psi }_{\alpha })(x,ϵ,\Omega ),$$(9) is expensive due to the high dimension of phase space of ${\psi }_{\alpha }(x,ϵ,\Omega )$ and complexity of the scattering operator $\mathcal{Q}.$ Therefore, we employ the method of moments to derive a modal approximation of the linear Boltzmann equation. We reduce the phase space by replacing the direction Ω as an independent variable by a small set of moments to model the electron fluence ψ. Here we consider the spherical harmonic (P_N) moment approximation (Case and Zweifel 1967; Mark 1944; Mark 1945; Seibold and Frank 2014), i.e., we approximate the electron fluence ψ by its series expansion in spherical harmonics $Y_l^k:S^2\to \mathbb{R}$ truncated after some finite degree $l\le N\in \mathbb{N}.$(15) $${\psi }_{\alpha }(x,ϵ,\Omega )\approx {\psi }_{\alpha }^{P_N}(x,ϵ,\Omega )=\sum _{l=0}^N\sum _{k=-l}^l{u}_{\alpha }^{l,k}(x,ϵ)Y_l^k(\Omega )$$(15) where the expansion coefficients are the spherical harmonic moments of the angular distribution(16) $$u_{\alpha }^{l,k}(x,ϵ)={\int }_{S^2}{Y}_l^k(\Omega ){\psi }_{\alpha }(x,ϵ,\Omega )\mathrm{}\text{d}\Omega \hspace{1em}\forall l\le N,k\le \left|l\right|,$$(16) and we model the fluence through the system of evolution equations of its moments $u_{\alpha }^{l,k}$ which we refer to as the P_N equations. The P_N equations are obtained by weakly enforcing Equation (9)(9) $$-{\partial }_{ϵ}(S[\rho ](x,ϵ){\psi }_{\alpha }(x,ϵ,\Omega ))+\Omega ·{\nabla }_x{\psi }_{\alpha }(x,ϵ,\Omega )=(\mathcal{Q}[\rho ]{\psi }_{\alpha })(x,ϵ,\Omega ),$$(9) for ${\psi }_{\alpha }^{P_N}$ on the test space spanned by spherical harmonics up to degree N. By gathering the coefficients $u_{\alpha }^{l,k}$ and ansatz functions $Y_l^k$ into vector functions $u_{\alpha }:G\times E\to {\mathbb{R}}^{{(N+1)}^2}$ and ${\Upsilon }_N:S^2\to {\mathbb{R}}^{{(N+1)}^2}$ such that ${\psi }_{\alpha }^{P_N}={({\Upsilon }_N)}^T{u}_{\alpha },$ the equation for $u_{\alpha }$ can be written as $F_{\alpha }(u_{\alpha },\rho )=0.$(17) $$\underset{=F_{\alpha }(u_{\alpha },\rho )}{\underbrace{-{\partial }_{ϵ}(S[\rho ](x,ϵ)u_{\alpha }(x,ϵ))+\sum _{d=1}^3{A}_{(d)}{\partial }_{x_d}{u}_{\alpha }(x,ϵ)-Q\left[\rho \right](x,ϵ)u_{\alpha }(x,ϵ)}}=0\hspace{1em}\forall x\in G,ϵ\in E$$(17)

With the following definitions for $A_{(d)}$ and $Q[\rho ].$(18) $$A_{(d)}={\int }_{S^2}{\Omega }_d{\Upsilon }_N{({\Upsilon }_N)}^T\mathrm{}\text{d}\Omega$$(18) (19) $$Q\left[\rho \right](x,ϵ)={\int }_{S^2}{\Upsilon }_N\mathcal{Q}\left[\rho \right]{({\Upsilon }_N)}^T\mathrm{}\text{d}\Omega =\sum _{i=1}^{n_e}{\rho }_i(x)Q_i(ϵ)$$(19)

Using the method of moments, we reduced the phase space by two dimensions at the cost of replacing a scalar equation with a system of ${(N+1)}^2$ equations. Additionally, $A_{(d)}$ and Q exhibit structures that can be exploited to develop efficient numerical solvers. $Q[\rho ](x,ϵ)$ is diagonal as spherical harmonics are eigenfunctions of the collision operator $\mathcal{Q}[\rho ]$ such that the computational cost of the scattering operation grows linearly in the number of moments.

The transport matrices $A_{(d)}$ inherit a sparsity pattern from the recursion relation of spherical harmonics. When ordering spherical harmonics that are odd (even) in Cartesian direction d into a vector function ${\Upsilon }_N^{[e,d]}:S^2\to {\mathbb{R}}^{(N+1)(N+2)/2}$ (resp. ${\Upsilon }_N^{[o,d]}:S^2\to {\mathbb{R}}^{N(N+1)/2}$) and choosing ${\Upsilon }_N={({({\Upsilon }_N^{[e,d]})}^T,{({\Upsilon }_N^{[o,d]})}^T)}^T$ the matrix $A_{(d)}$ is of the following form.(20) $$A_{(d)}=\left(\begin{array}{cc}0& {\widehat{A}}_{(d)}^T\\ {\widehat{A}}_{(d)}& 0\end{array}\right)$$(20) (21) $${\widehat{A}}_{(d)}={\int }_{S^2}{\Omega }_d{\Upsilon }_N^{[o,d]}{\left({\Upsilon }_N^{[e,d]}\right)}^T\mathrm{}\text{d}\Omega$$(21)

Hence, the even (odd) moments $u_{\alpha }^{[e,d]}$ couple with the odd (even) moments $u_{\alpha }^{[o,d]}$ only through their spatial derivatives in direction d, which can be employed by central finite difference schemes on staggered grids (Seibold and Frank 2014).

We employ the boundary conditions introduced in Bünger, Sarna, and Torrilhon (2022). The boundary conditions follow from taking half moments of the Boltzmann boundary condition (13)(13) $${\psi }_{\alpha }(x,ϵ,\Omega )={\psi }_{\alpha ,\text{in}}(x,ϵ,\Omega )\hspace{1em}\forall \Omega \in S^2:n·\Omega <0$$(13) with odd (in direction $d\in \{1,2,3\}$) spherical harmonics as weighting functions and a stabilization step that is similar to the truncation of the series expansion in (15)(15) $${\psi }_{\alpha }(x,ϵ,\Omega )\approx {\psi }_{\alpha }^{P_N}(x,ϵ,\Omega )=\sum _{l=0}^N\sum _{k=-l}^l{u}_{\alpha }^{l,k}(x,ϵ)Y_l^k(\Omega )$$(15) . For a given boundary point $x\in \partial G$ with outward pointing boundary normal e_d we formulate boundary conditions for the moment vector $u_{\alpha }={({(u_{\alpha }^{[e,d]})}^T,{(u_{\alpha }^{[o,d]})}^T)}^T.$ The boundary conditions are given by(22) $$u_{\alpha }^{[o,d]}=L_{(d)}{\widehat{A}}_{(d)}{u}_{\alpha }^{[e,d]}+g_{\alpha }$$(22) (23) $$L_{(d)}=2{\int }_{e_d^T\Omega \ge 0}\frac{1}{{\Omega }_d}{\Upsilon }_N^{[o,d]}{\left({\Upsilon }_N^{[o,d]}\right)}^T\mathrm{}\text{d}\Omega \text{ and }{g}_{\alpha }={\int }_{e_d^T\Omega \le 0}{\psi }_{\alpha ,\text{in}}{\Upsilon }_N^{[o,d]}\mathrm{}\text{d}\Omega .$$(23)

These boundary conditions are compatible with the characteristics of the P_N equations (17)(17) $$\underset{=F_{\alpha }(u_{\alpha },\rho )}{\underbrace{-{\partial }_{ϵ}(S[\rho ](x,ϵ)u_{\alpha }(x,ϵ))+\sum _{d=1}^3{A}_{(d)}{\partial }_{x_d}{u}_{\alpha }(x,ϵ)-Q\left[\rho \right](x,ϵ)u_{\alpha }(x,ϵ)}}=0\hspace{1em}\forall x\in G,ϵ\in E$$(17) and allow us to assure energy stability, see Bünger, Sarna, and Torrilhon (2022) for details. From (12)(12) $${\psi }_{\alpha }(x,ϵ={ϵ}_{\max },\Omega )={\psi }_{\alpha ,0}(x,\Omega )=0.$$(12) we derive the following initial condition for the moments $u_{\alpha }.$(24) $$u_{\alpha }(x,ϵ={ϵ}_{\max })={\int }_{S^2}{\Upsilon }_N(\Omega ){\psi }_{\alpha ,0}(x,\Omega )\mathrm{}\text{d}\Omega =0$$(24)

The P_N equations (17)(17) $$\underset{=F_{\alpha }(u_{\alpha },\rho )}{\underbrace{-{\partial }_{ϵ}(S[\rho ](x,ϵ)u_{\alpha }(x,ϵ))+\sum _{d=1}^3{A}_{(d)}{\partial }_{x_d}{u}_{\alpha }(x,ϵ)-Q\left[\rho \right](x,ϵ)u_{\alpha }(x,ϵ)}}=0\hspace{1em}\forall x\in G,ϵ\in E$$(17) describe the evolution of all ${(N+1)}^2$ moments $u_{\alpha }$ of ${\psi }_{\alpha },$ whereas the intensity model (8)(8) $${\mathcal{I}}_{\alpha }\left[\rho \right](x)={\int }_E{\sigma }_{\alpha }(ϵ){\int }_{S^2}\underset{={\psi }_{\alpha }(x,ϵ,\Omega )}{\underbrace{|v(ϵ)|f(x,ϵ,\Omega )}}\mathrm{}\text{d}\Omega \mathrm{}\text{d}ϵ$$(8) only requires the angular average of ${\psi }_{\alpha }$ because we assumed the isotropic emission of x-rays. The spherical harmonic $Y_0^0(\Omega )$ is constant, so we identify the moment $u_{\alpha }^{0,0}\propto {\int }_{S^2}{\psi }_{\alpha }(x,ϵ,\Omega )\mathrm{}\text{d}\Omega$ as part of our x-ray emission model and simplify the integral (8)(8) $${\mathcal{I}}_{\alpha }\left[\rho \right](x)={\int }_E{\sigma }_{\alpha }(ϵ){\int }_{S^2}\underset{={\psi }_{\alpha }(x,ϵ,\Omega )}{\underbrace{|v(ϵ)|f(x,ϵ,\Omega )}}\mathrm{}\text{d}\Omega \mathrm{}\text{d}ϵ$$(8) to the following.(25) $${\mathcal{I}}_{\alpha }(x)={\int }_E{\sigma }_{\alpha }(ϵ){\int }_{S^2}{\psi }_{\alpha }(x,ϵ,\Omega )\mathrm{}\text{d}\Omega \mathrm{}\text{d}ϵ\propto {\int }_E{\sigma }_{\alpha }(ϵ)u_{\alpha }^{0,0}(x,ϵ)\mathrm{}\text{d}ϵ$$(25)

2.5. Mass concentration model

Typically, the quantities of interest in EPMA are weight fractions ${\omega }_i(x)=\frac{M_i(\mathrm{}\text{d}x)}{M(\mathrm{}\text{d}x)}$ that form a partition of unity ${\sum }_{i=1}^{n_e}{\omega }_i(x)=1.$ In order to determine mass concentrations $\rho (x)$ from weight fractions $\omega (x),$ additional assumptions are necessary. We neglect the influence of different molecular or crystal structures and assume that the compound $\{{\omega }_1,\dots ,{\omega }_{n_e}\}$ is formed in a volume-preserving manner. Then the total volume is $V(\mathrm{}\text{d}x)={\sum }_{i=1}^{n_e}{V}_i(\mathrm{}\text{d}x)={\sum }_{i=1}^{n_e}\frac{M_i(\mathrm{}\text{d}x)}{{\rho }_i^{\text{pure}}}={\sum }_{i=1}^{n_e}\frac{{\omega }_i(x)M(\mathrm{}\text{d}x)}{{\rho }_i^{\text{pure}}}.$ We derive the following relation between mass fractions and mass concentrations.(26) $${\rho }_i(x)={\omega }_i(x)\underset{:={\rho }_{\text{tot}}(x)}{\underbrace{{\left(\sum _{j=1}^{n_e}\frac{{\omega }_j(x)}{{\rho }_j^{\text{pure}}}\right)}^{-1}}}$$(26)

As a first possibility, we model weight fractions ${\omega }_i^{pc}(x)$ as a piecewise-constant function. We subdivide the spatial domain G into $n_k\in \mathbb{N}$ disjoint subdomains G_k and introduce the indicator function ${\mathbb{1}}_{G_k}(x).$ To guarantee the partition of unity we only use $n_p=n_k\times (n_e-1)$ parameters $p_{i,k}^{pc}$ and define the parametrization(27) $${\omega }_i^{pc}(x)=\{\begin{array}{cc}\sum _{k=1}^{n_k}{p}_{i,k}^{pc}{\mathbb{1}}_{G_k}(x)\hfill & i\in \{1,\dots n_e-1\}\hfill \\ 1-\sum _{k=1}^{n_k}\sum _{j=1}^{n_e-1}{p}_{j,k}^{pc}{\mathbb{1}}_{G_k}(x)\hfill & i=n_e\hfill \end{array}$$(27) with the additional constraints, that $p_{i,k}^{pc}\ge 0$ and ${\sum }_{j=1}^{n_e-1}{p}_{j,k}^{pc}\le 1.$ Together with (26)(26) $${\rho }_i(x)={\omega }_i(x)\underset{:={\rho }_{\text{tot}}(x)}{\underbrace{{\left(\sum _{j=1}^{n_e}\frac{{\omega }_j(x)}{{\rho }_j^{\text{pure}}}\right)}^{-1}}}$$(26) , (27)(27) $${\omega }_i^{pc}(x)=\{\begin{array}{cc}\sum _{k=1}^{n_k}{p}_{i,k}^{pc}{\mathbb{1}}_{G_k}(x)\hfill & i\in \{1,\dots n_e-1\}\hfill \\ 1-\sum _{k=1}^{n_k}\sum _{j=1}^{n_e-1}{p}_{j,k}^{pc}{\mathbb{1}}_{G_k}(x)\hfill & i=n_e\hfill \end{array}$$(27) forms the piecewise-constant material parametrization ${\mathcal{P}}_{pc}:\mathbb{P}\to (G\to {\mathbb{R}}^{n_e}).$ In the examples, we will replace the weight fraction parametrization ${\omega }^{pc}(x)$ to investigate other parametrizations $\mathcal{P}(p).$

3. Differentiation of the forward model

As described in the introduction, we consider the inverse problem as a minimization problem and apply iterative gradient-based optimization methods to find a candidate solution $p^{*}.$(28) $$p^{*}=\underset{p\in \mathbb{P}}{\text{arg} \min }(J\hspace{0.17em}\circ \hspace{0.17em}k\hspace{0.17em}\circ \hspace{0.17em}\mathcal{P})(p)$$(28)

Gradient-based optimization methods require the gradient of the objective function $\frac{dJ}{dp}.$ An efficient implementation of the gradient of our model, that simultaneously remains extensible and modular (particularly regarding the parametrization) is a challenge. While automatic algorithmic differentiation tools have the potential to transform any numerical code, the current versions of AD tools in julia (Zygote.jl, ReverseDiff.jl) cannot be applied to our model due to memory limitations and lack of performance.

In the following, we present a method, that adopts the structural implementation of derivative code of numerical code from algorithmic differentiation (AD). We focus on the adjoint mode of AD, since we are mainly interested in the gradient of the scalar objective function. The advantage of the adjoint mode is that it scales with the number of outputs, whereas with numerical differentiation (finite differences) or the tangent mode of AD, the complexity scales with the number of inputs. Also, the analytical derivation of our method provides a clear understanding that facilitates its memory-efficient implementation.

The core idea of AD is to base the derivative computation on local function transformations that are automatically composed. This also achieves modularity, i.e., if a particular function is modified, only the corresponding transformation needs to be adjusted.

3.1. Adjoint algorithmic differentiation

Algorithmic differentiation (AD) (Burghardt et al. 2022; Hüser 2022; Innes 2018; Kreyszig 1991; Naumann 2011) defines a framework to compute derivatives of numerical program code. In contrast to symbolic or numerical differentiation (finite differences), algorithmic differentiation is without approximation error and modular, so the code remains extensible. Also, the automatic differentiation of numerical code is possible. We review the basic definitions of AD with a (unusual) focus on high/infinite dimensional spaces. Alongside a motivating example, we describe the systematic composition into code that is able to compute derivatives.

Consider a differentiable function $f:U\subseteq X\to Y,x\hspace{0.17em}\mapsto \hspace{0.17em}f(x)=y,$ where X and Y are Hilbert spaces and U is the open domain of f. Given an input $x\in U,$ we define the tangent operator $\frac{\partial f(x)}{\partial x}[·]:X\to Y$ as the linear bounded operator that satisfies(29) $$\underset{h\to 0}{\lim }\frac{f(x+h\dot{x})-f(x)}{h}=\frac{\partial f(x)}{\partial x}\left[\dot{x}\right]\hspace{1em}\forall \dot{x}\in X.$$(29)

Additionally, given a tangent $\dot{x}\in X$ of the input, we define the tangent $\dot{y}$ of the output by(30) $$\dot{y}=\frac{\partial f(x)}{\partial x}\left[\dot{x}\right].$$(30)

Given the tangent operator $\frac{\partial f(x)}{\partial x}[·],$ the corresponding adjoint operator ${\frac{\partial f(x)}{\partial x}}^{*}[·]:Y\to X$ is defined by(31) $${〈\overline{y},\frac{\partial f(x)}{\partial x}\left[\dot{x}\right]〉}_Y={〈{\frac{\partial f(x)}{\partial x}}^{*}\left[\overline{y}\right],\dot{x}〉}_X\hspace{1em}\forall \dot{x}\in X,\overline{y}\in Y$$(31) where ${〈·,·〉}_X$ and ${〈·,·〉}_Y$ are the respective inner products of X and Y. Similarly to the tangent, given an adjoint $\overline{y}\in Y$ of the output, we define the adjoint $\overline{x}$ of the input by(32) $$\overline{x}={\frac{\partial f(x)}{\partial x}}^{*}\left[\overline{y}\right].$$(32)

Consequently, the identity ${〈\overline{y},\dot{y}〉}_Y={〈\overline{x},\dot{x}〉}_X$ holds for all $\dot{x}\in X$ and $\overline{y}\in Y.$

In real, finite dimensional spaces (e.g., $X={\mathbb{R}}^n,Y={\mathbb{R}}^m$), the tangent operator can be represented by the Jacobian $Df(x)\in {\mathbb{R}}^{m\times n}$ and the adjoint operator by its transpose $Df{(x)}^T\in {\mathbb{R}}^{n\times m}.$ For examples in infinite dimensional (function) spaces, we refer to Sec. 3.2.

The core idea of AD is to subdivide a function into individual subfunctions (with individual tangent and adjoint operators) and intermediate variables (with intermediate tangents and adjoints). As an example, consider the following decomposition of $f:U\to Y,x\hspace{0.17em}\mapsto \hspace{0.17em}f(x)=y$ into two subfunctions $g:U\to V$ and $h:V\times U\to Y$ where we use an intermediate variable $v\in V$ (V is another Hilbert space).(33) $$\begin{array}{cc}v\hfill & =g(x)\hfill \\ y\hfill & =h(v,x)\hfill \end{array}$$(33)

Given $x\in U,v\in V$ and $y\in Y$ as well as an input tangent $\dot{x}\in X,$ the tangent $\dot{y}$ of is given by(34) $$\dot{y}=\frac{\partial f(x)}{\partial x}\left[\dot{x}\right]=\frac{\partial h(v,x)}{\partial v}[\underset{=\dot{v}}{\underbrace{\frac{\partial g(x)}{\partial x}\left[\dot{x}\right]}}]+\frac{\partial h(v,x)}{\partial x}\left[\dot{x}\right]$$(34) where we employ the chain rule to split the tangent operator of f into a combination of tangent operators of g and h. Additionally, we identify the tangent $\dot{v}$ of the intermediate variable.

To find the corresponding adjoint operator, we successively use adjoint operators of the subfunctions g and h to isolate $\dot{x}$ in the inner product identity(35) $$\begin{array}{cc}{〈\overline{y},\dot{y}〉}_Y={〈\overline{y},\frac{\partial h}{\partial v}[\dot{v}]+\frac{\partial h}{\partial x}[\dot{x}]〉}_Y\hfill & ={〈{\frac{\partial h}{\partial v}}^{*}[\overline{y}],\underset{=\dot{v}}{\underbrace{\frac{\partial g}{\partial x}[\dot{x}]}}〉}_V+{〈{\frac{\partial h}{\partial x}}^{*}[\overline{y}],\dot{x}〉}_X\hfill \\ \hfill & ={〈{\frac{\partial g}{\partial x}}^{*}[\underset{=\overline{v}}{\underbrace{{\frac{\partial h}{\partial v}}^{*}[\overline{y}]}}]+{\frac{\partial h}{\partial x}}^{*}[\overline{y}],\dot{x}〉}_X={〈\overline{x},\dot{x}〉}_X.\hfill \end{array}$$(35)

To summarize, given the adjoint $\overline{y}\in Y,$ the adjoints $\overline{v}$ and $\overline{x}$ are(36) $$\overline{x}={\frac{\partial g(x)}{\partial x}}^{*}[\underset{=\overline{v}}{\underbrace{{\frac{\partial h(v,x)}{\partial v}}^{*}\left[\overline{y}\right]}}]+{\frac{\partial h(v,x)}{\partial x}}^{*}\left[\overline{y}\right].$$(36)

Based on the chain rule, the composition of individual tangent and adjoint operators of a larger decomposition into an algorithm to compute the derivative is very systematic and can be implemented automatically. We describe the methodology using an arbitrary decomposition of a function $f:x\hspace{0.17em}\mapsto \hspace{0.17em}f(x)=y$ with n inputs and m outputs into subfunctions ${\varphi }_{i=1,\dots I}.$(37) $$\begin{array}{cc}\hfill & {(v_0,v_{-1},\dots ,v_{-(n-1)})}^T=x\hfill \\ \hfill & v_i={\varphi }_i({(v_j)}_{j\prec i})\hspace{1em}i=1,\dots ,I\hfill \\ \hfill & y={(v_I,v_{I-1},\dots v_{I-m})}^T.\hfill \end{array}$$(37)

With ${(v_j)}_{j\prec i}$ we denote all predecessors of v_i, i.e., all intermediate variables v_j that v_i directly depends on, and with ${(v_j)}_{j\succ i}$ we denote all successors of v_i. Tangent mode AD is the successive computation of tangents(38) $${\dot{v}}_i=\sum _{k\prec i}\frac{\partial {\varphi }_i({(v_j)}_{j\prec i})}{\partial v_k}\left[{\dot{v}}_k\right].$$(38)

Given an input tangent $\dot{x},$ a tangent ${\dot{v}}_i$ represents the directional derivative of v_i into direction $\dot{x}.$ Note that the tangent mode follows the sequence of operations in the evaluation of the primal. In adjoint mode AD, the sequence of operations is reversed. To compute an adjoint ${\overline{v}}_i,$ we consider the adjoints ${({\overline{v}}_k)}_{k\succ i}$ of all intermediate variables, that depend on v_i.(39) $${\overline{v}}_i=\sum _{k\succ i}{\frac{\partial {\varphi }_k({(v_j)}_{j\prec k})}{\partial v_i}}^{*}\left[{\overline{v}}_k\right].$$(39)

Usually, in AD literature the single assignment code is considered, where all intermediate variables $v_i\in \mathbb{R}$ are scalar, and the tangent operators (resp. multiplication with a scalar) is self-adjoint. Then the ${·}^{*}$ can be neglected. For real multivariate intermediate variables the adjoint is the transpose ${·}^{*}={·}^T.$

We briefly comment on some aspects of AD:

• Computational Complexity: Implementations of tangent and adjoint mode AD differ in algorithmic complexity. The Jacobian of a function $f:{\mathbb{R}}^n\to {\mathbb{R}}^m$ is obtained using tangent mode AD by $j=1,\dots ,n$ successive executions of (38)(38) $${\dot{v}}_i=\sum _{k\prec i}\frac{\partial {\varphi }_i({(v_j)}_{j\prec i})}{\partial v_k}\left[{\dot{v}}_k\right].$$(38) where the input tangents are seeded with unit vectors $\dot{x}=e_j.$ Elements of the Jacobian are then retrieved from ${(Df(x))}_{ij}={{\dot{y}}_i|}_{\dot{x}=e_j}.$ In contrast to the tangent mode, using adjoint mode AD the Jacobian is obtained by $i=1,\dots ,m$ successive applications of (39)(39) $${\overline{v}}_i=\sum _{k\succ i}{\frac{\partial {\varphi }_k({(v_j)}_{j\prec k})}{\partial v_i}}^{*}\left[{\overline{v}}_k\right].$$(39) where the output adjoints are seeded with unit vector $\overline{y}=e_i.$ Elements of the Jacobian are then retrieved from ${(Df(x))}_{ij}={{\overline{x}}_j|}_{\overline{y}=e_i}.$ Hence, to compute the full Jacobian, an implementation of tangent mode scales in the number of input variables while the adjoint mode scales in the number of output variables.

• Modularity and Extensibility: An implementation of either tangent or adjoint mode AD is modular and extensible. It consists of a systematic composition of individual tangent and adjoint operators. For example, changing one subfunction in (37)(37) $$\begin{array}{cc}\hfill & {(v_0,v_{-1},\dots ,v_{-(n-1)})}^T=x\hfill \\ \hfill & v_i={\varphi }_i({(v_j)}_{j\prec i})\hspace{1em}i=1,\dots ,I\hfill \\ \hfill & y={(v_I,v_{I-1},\dots v_{I-m})}^T.\hfill \end{array}$$(37) only affects its corresponding tangent and adjoint operator in (38)(38) $${\dot{v}}_i=\sum _{k\prec i}\frac{\partial {\varphi }_i({(v_j)}_{j\prec i})}{\partial v_k}\left[{\dot{v}}_k\right].$$(38) resp. (39)(39) $${\overline{v}}_i=\sum _{k\succ i}{\frac{\partial {\varphi }_k({(v_j)}_{j\prec k})}{\partial v_i}}^{*}\left[{\overline{v}}_k\right].$$(39) . Implementations of the other operators remain.

• Computational Graph: Both modes of AD can be visualized using a computational graph consisting of vertices for each variable (input, output, intermediate) and directed edges for direct dependence of variables. Assume a labeling of the edges with their corresponding tangent/adjoint operator applied to the tail(tangent)/head(adjoint) of the edge. The tangent is defined as the sum of all incoming edge tangent operators, while the adjoint is defined as the sum of all outgoing edge adjoint operators (c.f. Equations (38)(38) $${\dot{v}}_i=\sum _{k\prec i}\frac{\partial {\varphi }_i({(v_j)}_{j\prec i})}{\partial v_k}\left[{\dot{v}}_k\right].$$(38) and (39)(39) $${\overline{v}}_i=\sum _{k\succ i}{\frac{\partial {\varphi }_k({(v_j)}_{j\prec k})}{\partial v_i}}^{*}\left[{\overline{v}}_k\right].$$(39) ).

• Checkpointing – Adjoint Memory Consumption: While the tangents can be implemented alongside the execution of the primal function and intermediate primal variables (x, v) are available, they must be kept in memory for the computation of adjoints. Storing all intermediate variables is expensive. Checkpointing remedies the high memory consumption by recomputing certain intermediate variables while computing the adjoint operators. In the example (36)(36) $$\overline{x}={\frac{\partial g(x)}{\partial x}}^{*}[\underset{=\overline{v}}{\underbrace{{\frac{\partial h(v,x)}{\partial v}}^{*}\left[\overline{y}\right]}}]+{\frac{\partial h(v,x)}{\partial x}}^{*}\left[\overline{y}\right].$$(36) , the primal $v=g(x)$ would be recomputed while evaluating the adjoint of h.

Our model includes the special case, that the mapping from $\rho (x)$ to $u_{\alpha }(x,ϵ)$ is given implicitly by the PDE in Equation (17)(17) $$\underset{=F_{\alpha }(u_{\alpha },\rho )}{\underbrace{-{\partial }_{ϵ}(S[\rho ](x,ϵ)u_{\alpha }(x,ϵ))+\sum _{d=1}^3{A}_{(d)}{\partial }_{x_d}{u}_{\alpha }(x,ϵ)-Q\left[\rho \right](x,ϵ)u_{\alpha }(x,ϵ)}}=0\hspace{1em}\forall x\in G,ϵ\in E$$(17) . We state the two basic steps of calculating adjoints in implicit relations and refer to the Appendix A for a more detailed derivation. Consider the implicit definition $F(y=g(x),x)=0$ (with differentiable F(y, x)), where we assume that it uniquely defines the differentiable mapping $y=g(x).$ We introduce an additional intermediate adjoint variable $\overline{F}.$ Then, given the adjoint $\overline{y}$ of the output, the adjoint $\overline{x}$ of the input can be computed from(40) $${\frac{\partial F}{\partial y}}^{*}\left[\overline{F}\right]=\overline{y}\hspace{1em}\text{ and }\hspace{1em}\overline{x}=-{\frac{\partial F}{\partial x}}^{*}\left[\overline{F}\right]$$(40) without explicit definition of $\overline{x}={\frac{\partial g(x)}{\partial x}}^{*}[\overline{y}].$

3.2. Differentiation method for our model

In the following, we will successively define adjoint operators for the individual parts of the model described in Sec. 2. All derivations follow the same scheme: First, we define the tangent operator using the directional derivative and second we derive the adjoint operator from the inner product identity. We will elaborate on complex derivations and only state the adjoint operator for simple ones. In this paper, the only relevant derivatives are with respect to the mass concentrations ρ, hence all quantities that do not depend on ρ are treated as constants and their adjoints are omitted.

Compared to the forward computation, adjoints are computed in reversed order. We follow this behavior while deriving the adjoint operators for the individual dependencies. For brevity of notation we omit the accumulation of adjoints (tangents) and respectively only consider the operators that describe the direct propagation. The accumulation becomes apparent in the high-level computational graph in Figure 2. If multiple edges leave a node, the adjoints of all outgoing edges have to be summed up.

• $[{\overline{k}}_{\alpha }={\frac{\partial J}{\partial k_{\alpha }}}^{*}\overline{J}]:$ The final operator (1)(1) $$p^{*}=\underset{p\in \mathbb{P}}{\text{arg} \min }\underset{=(J\hspace{0.17em}\circ \hspace{0.17em}k\hspace{0.17em}\circ \hspace{0.17em}\mathcal{P})(p)}{\underbrace{{\sum }_{\alpha \in \mathbb{A}}((k_{\alpha } \circ \mathcal{P})(p)-k_{\alpha }^{\hspace{0.17em}\exp \hspace{0.17em}}{)}^2}}$$(1) of our model computes the squared error J from a set of simulated k-ratios $k_{\alpha },\alpha \in \{{\alpha }_1,\dots \}.$ Given the adjoint $\overline{J},$ which for the computation of derivatives $\frac{\partial J}{\partial ·}$ must be chosen $\overline{J}=1,$ the definition of the adjoint operator of the squared error from Equation (1)(1) $$p^{*}=\underset{p\in \mathbb{P}}{\text{arg} \min }\underset{=(J\hspace{0.17em}\circ \hspace{0.17em}k\hspace{0.17em}\circ \hspace{0.17em}\mathcal{P})(p)}{\underbrace{{\sum }_{\alpha \in \mathbb{A}}((k_{\alpha } \circ \mathcal{P})(p)-k_{\alpha }^{\hspace{0.17em}\exp \hspace{0.17em}}{)}^2}}$$(1) is trivial.(41) $${\overline{k}}_{\alpha }=2(k_{\alpha }-k_{\alpha }^{\hspace{0.17em}\exp \hspace{0.17em}})\overline{J}$$(41)

• $[{\overline{I}}_{\alpha }={\frac{\partial k_{\alpha }}{\partial I_{\alpha }}}^{*}{\overline{k}}_{\alpha }]:$ Also the adjoint operator for normalization with standard intensities in Equation (4)(4) $$k_{\alpha }\left[\rho \right]=\frac{I_{\alpha }\left[\rho \right]}{I_{\alpha }\left[{\rho }_{\alpha }^{\text{std}}\right]}$$(4) is derived by simple means.(42) $${\overline{I}}_{\alpha }=\frac{{\overline{k}}_{\alpha }}{I_{\alpha }^{\text{std}}}$$(42)

• $[\{{\overline{\mathcal{A}}}_{\alpha },{\overline{\mathcal{N}}}_{\alpha },{\overline{\mathcal{I}}}_{\alpha }\}={\frac{\partial I_{\alpha }}{\partial \{{\mathcal{A}}_{\alpha },{\mathcal{N}}_{\alpha },{\mathcal{I}}_{\alpha }\}}}^{*}{\overline{I}}_{\alpha }]:$ We elaborate on the adjoint operators derived from Equation (5)(5) $$I_{\alpha }\left[\rho \right]={\int }_G{\mathcal{A}}_{\alpha }\left[\rho \right](x) {\mathcal{N}}_{\alpha }\left[\rho \right](x) {\mathcal{I}}_{\alpha }\left[\rho \right](x)\text{d}x$$(5) . Linearity of the integral in (5)(5) $$I_{\alpha }\left[\rho \right]={\int }_G{\mathcal{A}}_{\alpha }\left[\rho \right](x) {\mathcal{N}}_{\alpha }\left[\rho \right](x) {\mathcal{I}}_{\alpha }\left[\rho \right](x)\text{d}x$$(5) directly yields the tangent ${\dot{I}}_{\alpha }.$ Using the standard definition of the inner products in $\mathbb{R},{〈a,b〉}_{\mathbb{R}}=ab$ and in $L^2(G),{〈a,b〉}_{L^2(G)}={\int }_G{a}^T(x)b(x)\mathrm{}\text{d}x$ we derive adjoints ${\overline{\mathcal{A}}}_{\alpha },{\overline{\mathcal{N}}}_{\alpha }$ and ${\overline{\mathcal{I}}}_{\alpha }.$ We detail ${\overline{\mathcal{A}}}_{\alpha }.$(43) $$\begin{array}{cc}{〈{\overline{I}}_{\alpha },{\dot{I}}_{\alpha }〉}_{\mathbb{R}}\hfill & ={〈{\overline{I}}_{\alpha },{\int }_G{\mathcal{N}}_{\alpha }(x){\mathcal{I}}_{\alpha }(x){\dot{\mathcal{A}}}_{\alpha }(x)\mathrm{}\text{d}x〉}_{\mathbb{R}}\hfill \\ \hfill & ={\int }_G{\mathcal{N}}_{\alpha }(x){\mathcal{I}}_{\alpha }(x){\overline{I}}_{\alpha }{\dot{\mathcal{A}}}_{\alpha }(x)\mathrm{}\text{d}x\hfill \\ \hfill & ={〈\underset{:={\overline{\mathcal{A}}}_{\alpha }(·)}{\underbrace{{\mathcal{N}}_{\alpha }(·){\mathcal{I}}_{\alpha }(·){\overline{I}}_{\alpha }}},{\dot{\mathcal{A}}}_{\alpha }(·)〉}_{L^2(G)}\hfill \end{array}$$(43)

  Accordingly, the adjoints ${\overline{\mathcal{N}}}_{\alpha }(·)$ and ${\overline{\mathcal{I}}}_{\alpha }(·)$ are derived.(44) $${\overline{\mathcal{N}}}_{\alpha }(x)={\mathcal{A}}_{\alpha }(x){\mathcal{I}}_{\alpha }(x){\overline{I}}_{\alpha }\hspace{1em}\text{and}\hspace{1em}{\overline{\mathcal{I}}}_{\alpha }(x)={\mathcal{N}}_{\alpha }(x){\mathcal{A}}_{\alpha }(x){\overline{I}}_{\alpha }$$(44)

• $[\overline{\rho }={\frac{\partial {\mathcal{N}}_{\alpha }}{\partial \rho }}^{*}{\overline{\mathcal{N}}}_{\alpha }]:$ From Equation (7)(7) $${\mathcal{N}}_{\alpha }\left[\rho \right](x)=\frac{{{\rho }_i(x)|}_{Z_i\in \alpha }}{A_{\alpha }}$$(7) , we find the adjoint operator.(45) $${{\overline{\rho }}_i(x)|}_{Z_i\in \alpha }=\frac{{\overline{\mathcal{N}}}_{\alpha }(x)}{A_{\alpha }}$$(45)

• $[\overline{\rho }={\frac{\partial {\mathcal{A}}_{\alpha }}{\partial \rho }}^{*}{\overline{\mathcal{A}}}_{\alpha }]:$ Deriving the adjoint from the attenuation operator given in (6)(6) $${\mathcal{A}}_{\alpha }\left[\rho \right](x)=\hspace{0.17em}\exp \hspace{0.17em}(-{\int }_{l(x,x_d)}{\mu }_{\alpha }(y)\mathrm{}\text{d}y)=\hspace{0.17em}\exp \hspace{0.17em}(-{\int }_{l(x,x_d)}\sum _{i=1}^{n_e}{\tau }_{\alpha ,i}{\rho }_i(y)\mathrm{}\text{d}y)$$(6) requires more analysis. The tangent of the absorption operator ${\dot{\mathcal{A}}}_{\alpha }[\rho ](x)$ is given by the directional derivative into direction $\dot{\rho }(x).$(46) $${\dot{\mathcal{A}}}_{\alpha }\left[\rho \right](x)=-{\mathcal{A}}_{\alpha }\left[\rho \right](x){\int }_{l(x,x_d)}\sum _{i=1}^{n_e}{\tau }_{\alpha ,i}{\dot{\rho }}_i(y)\mathrm{}\text{d}y$$(46)

  By inserting the tangent ${\dot{\mathcal{A}}}_{\alpha }$ into the definition of the adjoint (31)(31) $${〈\overline{y},\frac{\partial f(x)}{\partial x}\left[\dot{x}\right]〉}_Y={〈{\frac{\partial f(x)}{\partial x}}^{*}\left[\overline{y}\right],\dot{x}〉}_X\hspace{1em}\forall \dot{x}\in X,\overline{y}\in Y$$(31) , we derive the adjoint of the attenuation operator. In the following derivation, we rewrite the line integral from (46)(46) $${\dot{\mathcal{A}}}_{\alpha }\left[\rho \right](x)=-{\mathcal{A}}_{\alpha }\left[\rho \right](x){\int }_{l(x,x_d)}\sum _{i=1}^{n_e}{\tau }_{\alpha ,i}{\dot{\rho }}_i(y)\mathrm{}\text{d}y$$(46) as an integral over the whole space G by introducing $\delta (d(x,x_d,y))$ to denote the Dirac, that indicates the line segment $l(x,x_d).$ So $d(x,x_d,y)$ denotes the distance of y to the line segment $l(x,x_d)$ between x and x_d.(47) $$\begin{array}{cc}{〈{\overline{\mathcal{A}}}_{\alpha },{\dot{\mathcal{A}}}_{\alpha }〉}_{L^2(G)}\hfill & ={〈{\overline{\mathcal{A}}}_{\alpha }(·),-{\mathcal{A}}_{\alpha }\left[\rho \right](·){\int }_{l(·,x_d)}\sum _{i=1}^{n_e}{\tau }_{\alpha ,i}{\dot{\rho }}_i(y)\mathrm{}\text{d}y〉}_{L^2(G)}\hfill \\ \hfill & ={\int }_G-{\overline{\mathcal{A}}}_{\alpha }(x){\mathcal{A}}_{\alpha }\left[\rho \right](x){\int }_{l(x,x_d)}\sum _{i=1}^{n_e}{\tau }_{\alpha ,i}{\dot{\rho }}_i(y)\mathrm{}\text{d}y\mathrm{}\text{d}x\hfill \\ \hfill & ={\int }_G{\int }_G\sum _{i=1}^{n_e}-{\overline{\mathcal{A}}}_{\alpha }(x){\mathcal{A}}_{\alpha }\left[\rho \right](x){\tau }_{\alpha ,i}{\dot{\rho }}_i(y)\delta (d(x,x_d,y))\mathrm{}\text{d}y\mathrm{}\text{d}x\hfill \\ \hfill & ={\int }_G\sum _{i=1}^{n_e}{\int }_G-{\overline{\mathcal{A}}}_{\alpha }(x){\mathcal{A}}_{\alpha }\left[\rho \right](x){\tau }_{\alpha ,i}\delta (d(x,x_d,y))\mathrm{}\text{d}x{\dot{\rho }}_i(y)\mathrm{}\text{d}y\hfill \\ \hfill & ={\int }_G\sum _{i=1}^{n_e}{\int }_{l(y,x_{d^{*}})}-{\overline{\mathcal{A}}}_{\alpha }(x){\mathcal{A}}_{\alpha }\left[\rho \right](x){\tau }_{\alpha ,i}\mathrm{}\text{d}x{\dot{\rho }}_i(y)\mathrm{}\text{d}y\hfill \\ \hfill & ={〈{\int }_{l(·,x_{d^{*}})}-{\overline{\mathcal{A}}}_{\alpha }(x){\mathcal{A}}_{\alpha }\left[\rho \right](x)\left(\begin{array}{c}{\tau }_{\alpha ,1}\\ ⋮\\ {\tau }_{\alpha ,n_e}\end{array}\right)\mathrm{}\text{d}x,\dot{\rho }(·)〉}_{L^2(G)}\hfill \\ \hfill & ={〈\overline{\rho },\dot{\rho }〉}_{L^2(G)}\hfill \end{array}$$(47)

  By swapping the order of integration $\mathrm{}\text{d}y\mathrm{}\text{d}x$ to $\mathrm{}\text{d}x\mathrm{}\text{d}y,$ we are able to separate ${\dot{\rho }}_i,$ but the interpretation of $\delta (d(x,x_d,y))$ changes. For a given y, $\delta (d(x,x_d,y))$ now indicates the points x of all lines $l(x,x_d)$ that contain y. In other words: $\delta (d(x,x_d,y))$ indicates the line $l(y,x_{d^{*}})$ from y to the point reflection $x_{d^{*}}\in \partial G$ of $x_d\in \partial G$ at $x\in G.$ Exemplarily, we illustrate $l(x,x_d),l(x,x_{d^{*}})$ and $d(x,x_d,y)$ in Figure 3. From (47)(47) $$\begin{array}{cc}{〈{\overline{\mathcal{A}}}_{\alpha },{\dot{\mathcal{A}}}_{\alpha }〉}_{L^2(G)}\hfill & ={〈{\overline{\mathcal{A}}}_{\alpha }(·),-{\mathcal{A}}_{\alpha }\left[\rho \right](·){\int }_{l(·,x_d)}\sum _{i=1}^{n_e}{\tau }_{\alpha ,i}{\dot{\rho }}_i(y)\mathrm{}\text{d}y〉}_{L^2(G)}\hfill \\ \hfill & ={\int }_G-{\overline{\mathcal{A}}}_{\alpha }(x){\mathcal{A}}_{\alpha }\left[\rho \right](x){\int }_{l(x,x_d)}\sum _{i=1}^{n_e}{\tau }_{\alpha ,i}{\dot{\rho }}_i(y)\mathrm{}\text{d}y\mathrm{}\text{d}x\hfill \\ \hfill & ={\int }_G{\int }_G\sum _{i=1}^{n_e}-{\overline{\mathcal{A}}}_{\alpha }(x){\mathcal{A}}_{\alpha }\left[\rho \right](x){\tau }_{\alpha ,i}{\dot{\rho }}_i(y)\delta (d(x,x_d,y))\mathrm{}\text{d}y\mathrm{}\text{d}x\hfill \\ \hfill & ={\int }_G\sum _{i=1}^{n_e}{\int }_G-{\overline{\mathcal{A}}}_{\alpha }(x){\mathcal{A}}_{\alpha }\left[\rho \right](x){\tau }_{\alpha ,i}\delta (d(x,x_d,y))\mathrm{}\text{d}x{\dot{\rho }}_i(y)\mathrm{}\text{d}y\hfill \\ \hfill & ={\int }_G\sum _{i=1}^{n_e}{\int }_{l(y,x_{d^{*}})}-{\overline{\mathcal{A}}}_{\alpha }(x){\mathcal{A}}_{\alpha }\left[\rho \right](x){\tau }_{\alpha ,i}\mathrm{}\text{d}x{\dot{\rho }}_i(y)\mathrm{}\text{d}y\hfill \\ \hfill & ={〈{\int }_{l(·,x_{d^{*}})}-{\overline{\mathcal{A}}}_{\alpha }(x){\mathcal{A}}_{\alpha }\left[\rho \right](x)\left(\begin{array}{c}{\tau }_{\alpha ,1}\\ ⋮\\ {\tau }_{\alpha ,n_e}\end{array}\right)\mathrm{}\text{d}x,\dot{\rho }(·)〉}_{L^2(G)}\hfill \\ \hfill & ={〈\overline{\rho },\dot{\rho }〉}_{L^2(G)}\hfill \end{array}$$(47) , we identify the adjoint operator of (6)(6) $${\mathcal{A}}_{\alpha }\left[\rho \right](x)=\hspace{0.17em}\exp \hspace{0.17em}(-{\int }_{l(x,x_d)}{\mu }_{\alpha }(y)\mathrm{}\text{d}y)=\hspace{0.17em}\exp \hspace{0.17em}(-{\int }_{l(x,x_d)}\sum _{i=1}^{n_e}{\tau }_{\alpha ,i}{\rho }_i(y)\mathrm{}\text{d}y)$$(6) .(48) $${\overline{\rho }}_i(x)={\tau }_{\alpha ,i}{\int }_{l(x,x_{d^{*}})}-{\overline{\mathcal{A}}}_{\alpha }(y){\mathcal{A}}_{\alpha }\left[\rho \right](y)\mathrm{}\text{d}y$$(48) Where the integration domain $l(x,x_{d^{*}})$ is the line between x and the reflection $x_{d^{*}}$ of x_d at x.

• $[{\overline{u}}_{\alpha }={\frac{\partial {\mathcal{I}}_{\alpha }}{\partial u_{\alpha }}}^{*}{\overline{\mathcal{I}}}_{\alpha }]:$ In Equation (25)(25) $${\mathcal{I}}_{\alpha }(x)={\int }_E{\sigma }_{\alpha }(ϵ){\int }_{S^2}{\psi }_{\alpha }(x,ϵ,\Omega )\mathrm{}\text{d}\Omega \mathrm{}\text{d}ϵ\propto {\int }_E{\sigma }_{\alpha }(ϵ)u_{\alpha }^{0,0}(x,ϵ)\mathrm{}\text{d}ϵ$$(25) , the only relevant moment for the ionization field ${\mathcal{I}}_{\alpha }$ is $u_{\alpha }^{(0,0)},$ hence the adjoint ${\overline{u}}_{\alpha }.$(49) $${\overline{u}}_{\alpha }^{l,k}(x,ϵ)=\{\begin{array}{cc}{\sigma }_{\alpha }(ϵ){\overline{\mathcal{I}}}_{\alpha }(x)\hfill & l,k=0\hfill \\ 0\hfill & \text{else}\hfill \end{array}$$(49)

• $[\overline{\rho }={\frac{\partial u_{\alpha }}{\partial \rho }}^{*}{\overline{u}}_{\alpha }]:$ The relation between ρ and $u_{\alpha }$ is implicitly given by the P_N equations and the boundary conditions, i.e., by the implicit PDE operator $F_{\alpha }(u_{\alpha },\rho )=0$ (17)(17) $$\underset{=F_{\alpha }(u_{\alpha },\rho )}{\underbrace{-{\partial }_{ϵ}(S[\rho ](x,ϵ)u_{\alpha }(x,ϵ))+\sum _{d=1}^3{A}_{(d)}{\partial }_{x_d}{u}_{\alpha }(x,ϵ)-Q\left[\rho \right](x,ϵ)u_{\alpha }(x,ϵ)}}=0\hspace{1em}\forall x\in G,ϵ\in E$$(17) . We describe the derivation of the adjoint operators (40)(40) $${\frac{\partial F}{\partial y}}^{*}\left[\overline{F}\right]=\overline{y}\hspace{1em}\text{ and }\hspace{1em}\overline{x}=-{\frac{\partial F}{\partial x}}^{*}\left[\overline{F}\right]$$(40) for the P_N equations. For a detailed derivation of (40)(40) $${\frac{\partial F}{\partial y}}^{*}\left[\overline{F}\right]=\overline{y}\hspace{1em}\text{ and }\hspace{1em}\overline{x}=-{\frac{\partial F}{\partial x}}^{*}\left[\overline{F}\right]$$(40) we refer to the Appendix A. Let us introduce the auxiliary adjoint ${\overline{F}}_{\alpha }$ and derive the operators ${\frac{\partial F_{\alpha }}{\partial u_{\alpha }}}^{*}[·]$ and ${\frac{\partial F_{\alpha }}{\partial \rho }}^{*}[·]$ in the following. Successively, both operators form the adjoint ${\frac{\partial u_{\alpha }}{\partial \rho }}^{*}[·].$ In particular, we elaborate on the identities(50) $${〈{\overline{F}}_{\alpha },\frac{\partial F}{\partial u_{\alpha }}{\dot{u}}_{\alpha }〉}_{L^2(G\times E)}={〈\underset{={\overline{u}}_{\alpha }}{\underbrace{{\frac{\partial F}{\partial u_{\alpha }}}^{*}{\overline{F}}_{\alpha }}},{\dot{u}}_{\alpha }〉}_{L^2(G\times E)}\hspace{1em}\text{and}\hspace{1em}-{〈{\overline{F}}_{\alpha },\frac{\partial F}{\partial \rho }\dot{\rho }〉}_{L^2(G\times E)}={〈\underset{=\overline{\rho }}{\underbrace{-{\frac{\partial F}{\partial \rho }}^{*}{\overline{F}}_{\alpha }}},\dot{\rho }〉}_{L^2(G\times E)}.$$(50)

• $[{\frac{\partial F_{\alpha }}{\partial u_{\alpha }}}^{*}{\overline{F}}_{\alpha }={\overline{u}}_{\alpha }]:$ The adjoint operator ${\frac{\partial F_{\alpha }}{\partial u_{\alpha }}}^{*}$ applied to ${\overline{F}}_{\alpha }$ describes the continuous adjoint equation.(51) $$\underset{={\frac{\partial F_{\alpha }}{\partial u_{\alpha }}}^{*}{\overline{F}}_{\alpha }}{\underbrace{S\left[\rho \right](x,ϵ){\partial }_{ϵ}{\overline{F}}_{\alpha }(x,ϵ)-\sum _{d=1}^3{A}_{(d)}{\partial }_{x_i}{\overline{F}}_{\alpha }(x,ϵ)-Q\left[\rho \right](x,ϵ){\overline{F}}_{\alpha }(x,ϵ)}}={\overline{u}}_{\alpha }(x,ϵ)$$(51) In analogy to the reversed application of adjoint operators, the adjoint equation describes the evolution of the adjoint ${\overline{F}}_{\alpha }$ in reversed direction (w.r.t. the energy ϵ). Initial and boundary conditions follow from the derivation.(52) $${\overline{F}}_{\alpha }(x,{ϵ}_{\min })=0$$(52) (53) $${\overline{F}}_{\alpha }^{[o,d]}=-L_{(d)}{\widehat{A}}_{(d)}{\overline{F}}_{\alpha }^{[e,n]}$$(53)

  Derivation Linearity of $F_{\alpha }$ in $u_{\alpha }$ directly yields the tangent $\frac{\partial F_{\alpha }}{\partial u_{\alpha }}{\dot{u}}_{\alpha }.$ Using integration by parts and symmetry of $A_{(d)}$ and Q, we obtain(54) $$\begin{array}{l}〈{\overline{F}}_{\alpha },\frac{\partial F_{\alpha }}{\partial u_{\alpha }}{\dot{u}}_{\alpha }〉={\int }_E{\int }_G-{({\partial }_{ϵ}(S{\dot{u}}_{\alpha }))}^T{\overline{F}}_{\alpha }+\sum _{d=1}^3{(A_{(d)}{\partial }_{x_d}{\dot{u}}_{\alpha })}^T{\overline{F}}_{\alpha }-{(Q{\dot{u}}_{\alpha })}^T{\overline{F}}_{\alpha } \mathrm{}\text{d}x\mathrm{\hspace{0.17em}}\text{d}ϵ\hfill \\ \begin{array}{cc}\hfill & ={\int }_E{\int }_G-{\partial }_{ϵ}(S{\dot{u}}_{\alpha }^T{\overline{F}}_{\alpha })+{\dot{u}}_{\alpha }^T(S{\partial }_{ϵ}{\overline{F}}_{\alpha })\hfill \\ \hfill & \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\hspace{1em}+\sum _{d=1}^3{\partial }_{x_d}({\dot{u}}_{\alpha }^T{A}_{(d)}{\overline{F}}_{\alpha })-{\dot{u}}_{\alpha }^T(A_{(d)}{\partial }_{x_d}{\overline{F}}_{\alpha })-{\dot{u}}_{\alpha }^T(Q{\overline{F}}_{\alpha }) \mathrm{}\text{d}x\mathrm{\hspace{0.17em}}\text{d}ϵ\hfill \end{array}\hfill & \hfill \\ \begin{array}{cc}\hfill & ={\int }_E{\int }_G{\dot{u}}_{\alpha }^T(S{\partial }_{ϵ}{\overline{F}}_{\alpha })-\sum _{d=1}^3{\dot{u}}_{\alpha }^T(A_{(d)}{\partial }_{x_d}{\overline{F}}_{\alpha })-{\dot{u}}_{\alpha }^T(Q{\overline{F}}_{\alpha }) \mathrm{}\text{d}x\mathrm{\hspace{0.17em}}\text{d}ϵ\hfill \\ \hfill & \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\hspace{1em}+\underset{=0}{\underbrace{{\int }_E{\int }_G-({\partial }_{ϵ}(S{\dot{u}}_{\alpha }^T{\overline{F}}_{\alpha }))+\sum _{d=1}^3{\partial }_{x_d}({\dot{u}}_{\alpha }^T{A}_{(d)}{\overline{F}}_{\alpha })\mathrm{}\text{d}x\mathrm{\hspace{0.17em}}\text{d}ϵ}}\hfill \end{array}\hfill & \hfill \\ \hfill & =〈{\frac{\partial F_{\alpha }}{\partial u_{\alpha }}}^{*}{\overline{F}}_{\alpha },{\dot{u}}_{\alpha }〉=〈{\overline{u}}_{\alpha },{\dot{u}}_{\alpha }〉\hfill \end{array}$$(54)

  The last two terms in (54)(54) $$\begin{array}{l}〈{\overline{F}}_{\alpha },\frac{\partial F_{\alpha }}{\partial u_{\alpha }}{\dot{u}}_{\alpha }〉={\int }_E{\int }_G-{({\partial }_{ϵ}(S{\dot{u}}_{\alpha }))}^T{\overline{F}}_{\alpha }+\sum _{d=1}^3{(A_{(d)}{\partial }_{x_d}{\dot{u}}_{\alpha })}^T{\overline{F}}_{\alpha }-{(Q{\dot{u}}_{\alpha })}^T{\overline{F}}_{\alpha } \mathrm{}\text{d}x\mathrm{\hspace{0.17em}}\text{d}ϵ\hfill \\ \begin{array}{cc}\hfill & ={\int }_E{\int }_G-{\partial }_{ϵ}(S{\dot{u}}_{\alpha }^T{\overline{F}}_{\alpha })+{\dot{u}}_{\alpha }^T(S{\partial }_{ϵ}{\overline{F}}_{\alpha })\hfill \\ \hfill & \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\hspace{1em}+\sum _{d=1}^3{\partial }_{x_d}({\dot{u}}_{\alpha }^T{A}_{(d)}{\overline{F}}_{\alpha })-{\dot{u}}_{\alpha }^T(A_{(d)}{\partial }_{x_d}{\overline{F}}_{\alpha })-{\dot{u}}_{\alpha }^T(Q{\overline{F}}_{\alpha }) \mathrm{}\text{d}x\mathrm{\hspace{0.17em}}\text{d}ϵ\hfill \end{array}\hfill & \hfill \\ \begin{array}{cc}\hfill & ={\int }_E{\int }_G{\dot{u}}_{\alpha }^T(S{\partial }_{ϵ}{\overline{F}}_{\alpha })-\sum _{d=1}^3{\dot{u}}_{\alpha }^T(A_{(d)}{\partial }_{x_d}{\overline{F}}_{\alpha })-{\dot{u}}_{\alpha }^T(Q{\overline{F}}_{\alpha }) \mathrm{}\text{d}x\mathrm{\hspace{0.17em}}\text{d}ϵ\hfill \\ \hfill & \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\hspace{1em}+\underset{=0}{\underbrace{{\int }_E{\int }_G-({\partial }_{ϵ}(S{\dot{u}}_{\alpha }^T{\overline{F}}_{\alpha }))+\sum _{d=1}^3{\partial }_{x_d}({\dot{u}}_{\alpha }^T{A}_{(d)}{\overline{F}}_{\alpha })\mathrm{}\text{d}x\mathrm{\hspace{0.17em}}\text{d}ϵ}}\hfill \end{array}\hfill & \hfill \\ \hfill & =〈{\frac{\partial F_{\alpha }}{\partial u_{\alpha }}}^{*}{\overline{F}}_{\alpha },{\dot{u}}_{\alpha }〉=〈{\overline{u}}_{\alpha },{\dot{u}}_{\alpha }〉\hfill \end{array}$$(54) vanish due to the initial and boundary conditions imposed to ${\dot{u}}_{\alpha }$ and ${\overline{F}}_{\alpha }.$ The former term is zero because of the conditions ${\dot{u}}_{\alpha }(·,{ϵ}_{\max })=0$ and ${\overline{F}}_{\alpha }(·,{ϵ}_{\min })=0.$(55) $${\int }_E{\int }_G({\partial }_{ϵ}(S{\dot{u}}_{\alpha }^T{\overline{F}}_{\alpha }))\mathrm{}\text{d}x\mathrm{\hspace{0.17em}}\text{d}ϵ={\int }_G\underset{=0\mathrm{ }\forall {\dot{u}}_{\alpha }}{\underbrace{(S{\dot{u}}_{\alpha }^T{\overline{F}}_{\alpha }){|}_{ϵ={ϵ}_{\max }}}}-\underset{=0\mathrm{ }\forall {\overline{F}}_{\alpha }}{\underbrace{(S{\dot{u}}_{\alpha }^T{\overline{F}}_{\alpha }){|}_{ϵ={ϵ}_{\min }}}}\mathrm{}\text{d}x=0$$(55)

  The latter term in (54)(54) $$\begin{array}{l}〈{\overline{F}}_{\alpha },\frac{\partial F_{\alpha }}{\partial u_{\alpha }}{\dot{u}}_{\alpha }〉={\int }_E{\int }_G-{({\partial }_{ϵ}(S{\dot{u}}_{\alpha }))}^T{\overline{F}}_{\alpha }+\sum _{d=1}^3{(A_{(d)}{\partial }_{x_d}{\dot{u}}_{\alpha })}^T{\overline{F}}_{\alpha }-{(Q{\dot{u}}_{\alpha })}^T{\overline{F}}_{\alpha } \mathrm{}\text{d}x\mathrm{\hspace{0.17em}}\text{d}ϵ\hfill \\ \begin{array}{cc}\hfill & ={\int }_E{\int }_G-{\partial }_{ϵ}(S{\dot{u}}_{\alpha }^T{\overline{F}}_{\alpha })+{\dot{u}}_{\alpha }^T(S{\partial }_{ϵ}{\overline{F}}_{\alpha })\hfill \\ \hfill & \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\hspace{1em}+\sum _{d=1}^3{\partial }_{x_d}({\dot{u}}_{\alpha }^T{A}_{(d)}{\overline{F}}_{\alpha })-{\dot{u}}_{\alpha }^T(A_{(d)}{\partial }_{x_d}{\overline{F}}_{\alpha })-{\dot{u}}_{\alpha }^T(Q{\overline{F}}_{\alpha }) \mathrm{}\text{d}x\mathrm{\hspace{0.17em}}\text{d}ϵ\hfill \end{array}\hfill & \hfill \\ \begin{array}{cc}\hfill & ={\int }_E{\int }_G{\dot{u}}_{\alpha }^T(S{\partial }_{ϵ}{\overline{F}}_{\alpha })-\sum _{d=1}^3{\dot{u}}_{\alpha }^T(A_{(d)}{\partial }_{x_d}{\overline{F}}_{\alpha })-{\dot{u}}_{\alpha }^T(Q{\overline{F}}_{\alpha }) \mathrm{}\text{d}x\mathrm{\hspace{0.17em}}\text{d}ϵ\hfill \\ \hfill & \mathrm{\hspace{1em}}\mathrm{\hspace{1em}}\hspace{1em}+\underset{=0}{\underbrace{{\int }_E{\int }_G-({\partial }_{ϵ}(S{\dot{u}}_{\alpha }^T{\overline{F}}_{\alpha }))+\sum _{d=1}^3{\partial }_{x_d}({\dot{u}}_{\alpha }^T{A}_{(d)}{\overline{F}}_{\alpha })\mathrm{}\text{d}x\mathrm{\hspace{0.17em}}\text{d}ϵ}}\hfill \end{array}\hfill & \hfill \\ \hfill & =〈{\frac{\partial F_{\alpha }}{\partial u_{\alpha }}}^{*}{\overline{F}}_{\alpha },{\dot{u}}_{\alpha }〉=〈{\overline{u}}_{\alpha },{\dot{u}}_{\alpha }〉\hfill \end{array}$$(54) vanishes due to the choice of the boundary conditions for ${\dot{u}}_{\alpha }$ and ${\overline{F}}_{\alpha }.$ With the divergence theorem ($n={(n_1,n_2,n_3)}^T$ the boundary normal)(56) $${\int }_G\sum _{d=1}^3{\partial }_{x_d}({\dot{u}}_{\alpha }^T{A}_{(d)}{\overline{F}}_{\alpha })dx={\int }_{\partial G}\sum _{d=1}^3({\dot{u}}_{\alpha }^T{A}_{(d)}{\overline{F}}_{\alpha })n_{d}dS,$$(56) and the boundary conditions ${\dot{u}}_e{}^{[o,d]}=L_{(d)}{\widehat{A}}_{(d)}{\dot{u}}_{\alpha }{}^{[e,d]}$ (${\dot{g}}_{\alpha }=0$) and ${\overline{F}}_{\alpha }^{[o,d]}=-L_{(d)}{\widehat{A}}_{(d)}{\overline{F}}_{\alpha }^{[e,d]}$ the integrand in (56)(56) $${\int }_G\sum _{d=1}^3{\partial }_{x_d}({\dot{u}}_{\alpha }^T{A}_{(d)}{\overline{F}}_{\alpha })dx={\int }_{\partial G}\sum _{d=1}^3({\dot{u}}_{\alpha }^T{A}_{(d)}{\overline{F}}_{\alpha })n_{d}dS,$$(56) is zero because(57) $$\begin{array}{l}{\dot{u}}_{\alpha }^T{A}_{(d)}{\overline{F}}_{\alpha }={\left(\begin{array}{c}{\dot{u}}_{\alpha }^{[e,d]}\\ {\dot{u}}_{\alpha }^{[o,d]}\end{array}\right)}^T\left(\begin{array}{cc}0& {\widehat{A}}_{(d)}^T\\ {\widehat{A}}_{(d)}& 0\end{array}\right)\left(\begin{array}{c}{\overline{F}}_{\alpha }^{[e,d]}\\ {\overline{F}}_{\alpha }^{[o,d]}\end{array}\right)\hfill \\ ={({\dot{u}}_{\alpha }^{[e,d]})}^T{\widehat{A}}_{(d)}^T{\overline{F}}_{\alpha }^{[o,d]}+{({\dot{u}}_{\alpha }^{[o,d]})}^T{\widehat{A}}_{(d)}{\overline{F}}_{\alpha }^{[e,d]}\hfill \\ =-{({\dot{u}}_{\alpha }^{[e,d]})}^T{\widehat{A}}_{(d)}^T(L_{(d)}{\widehat{A}}_{(d)}{\overline{F}}_{\alpha }^{[e,d]})+{(L_{(d)}{\widehat{A}}_{(d)}{\dot{u}}_{\alpha }{}^{[e,d]})}^T{\widehat{A}}_{(d)}{\overline{F}}_{\alpha }^{[e,d]}\hfill \\ ={({\dot{u}}_{\alpha }^{[e,d]})}^T\underset{=0}{\underbrace{(-{\widehat{A}}_{(d)}^T{L}_{(d)}{\widehat{A}}_{(d)}+{\widehat{A}}_{(d)}^T{L}_{(d)}^T{\widehat{A}}_{(d)})}}{\overline{F}}_{\alpha }^{[e,d]}=0.\hfill \end{array}$$(57)

  Using that $L_{(d)}=L_{(d)}^T,$ defined in (23)(23) $$L_{(d)}=2{\int }_{e_d^T\Omega \ge 0}\frac{1}{{\Omega }_d}{\Upsilon }_N^{[o,d]}{\left({\Upsilon }_N^{[o,d]}\right)}^T\mathrm{}\text{d}\Omega \text{ and }{g}_{\alpha }={\int }_{e_d^T\Omega \le 0}{\psi }_{\alpha ,\text{in}}{\Upsilon }_N^{[o,d]}\mathrm{}\text{d}\Omega .$$(23) , is symmetric.

• $[\overline{\rho }=-{\frac{\partial F_{\alpha }}{\partial \rho }}^{*}{\overline{F}}_{\alpha }]:$ The tangent $\frac{\partial F_{\alpha }}{\partial \rho }\dot{\rho }$ follows directly from linearity of $F_{\alpha }$ (and S and Q) in ρ. For the adjoint we find(58) $$\begin{array}{cc}-〈{\overline{F}}_{\alpha },\frac{\partial F_{\alpha }}{\partial \rho }\dot{\rho }〉\hfill & =-{〈{\overline{F}}_{\alpha },-{\partial }_{ϵ}(\sum _{i=1}^{n_e}{S}_i{\dot{\rho }}_i{u}_{\alpha })-\sum _{i=1}^{n_e}{Q}_i{\dot{\rho }}_i{u}_{\alpha }〉}_{L^2(G\times E)}\hfill \\ \hfill & =-{\displaystyle {\int }_G}\sum _{i=1}^{n_e}{\dot{\rho }}_i{\int }_E\sum _{l\le N,\left|k\right|\le l}{\overline{F}}_{\alpha }^{l,k}(-{\partial }_{ϵ}(S_i{u}_{\alpha }^{l,k})-Q_i{u}_{\alpha }^{l,k})\mathrm{}\text{d}ϵ\mathrm{}\text{d}x\hfill \\ \hfill & ={〈{\int }_E\sum _{l\le N,\left|k\right|\le l}{\overline{F}}_{\alpha }^{l,k}\left({\partial }_{ϵ}\right(\left(\begin{array}{c}{S}_1\\ \dots \\ S_{n_e}\end{array}\right)u_{\alpha }^{l,k})+\left(\begin{array}{c}{Q}_1\\ \dots \\ Q_{n_e}\end{array}\right)u_{\alpha }^{l,k})\mathrm{}\text{d}ϵ,\dot{\rho }〉}_{L^2(G)}\hfill \\ \hfill & ={〈\overline{\rho },\dot{\rho }〉}_{L^2(G)},\hfill \end{array}$$(58) and identify(59) $${\overline{\rho }}_i(x)={\int }_E{\overline{F}}_{\alpha }^T(x,ϵ)({\partial }_{ϵ}(S_i(ϵ)u_{\alpha }(x,ϵ))+Q_i(ϵ)u_{\alpha }(x,ϵ))\mathrm{}\text{d}ϵ.$$(59)

• $[{\overline{\omega }}^{pc}={\frac{\partial \rho }{\partial \omega }}^{*}\overline{\rho }]:$ The adjoint of the mass concentration model (26)(26) $${\rho }_i(x)={\omega }_i(x)\underset{:={\rho }_{\text{tot}}(x)}{\underbrace{{\left(\sum _{j=1}^{n_e}\frac{{\omega }_j(x)}{{\rho }_j^{\text{pure}}}\right)}^{-1}}}$$(26) can again be derived by simple means.(60) $${\overline{\omega }}_i^{pc}(x)={\rho }_{\text{tot}}(x)\overline{{\rho }_i}(x)-\frac{{\rho }_{\text{tot}}^2(x)}{{\rho }_i^{\text{pure}}}\sum _{k=1}^{n_e}{\omega }_k(x){\overline{\rho }}_k(x)$$(60)

• $[\overline{p}={\frac{\partial {\omega }^{pc}}{\partial p}}^{*}{\overline{\omega }}^{pc}]:$ The indicator functions in (27)(27) $${\omega }_i^{pc}(x)=\{\begin{array}{cc}\sum _{k=1}^{n_k}{p}_{i,k}^{pc}{\mathbb{1}}_{G_k}(x)\hfill & i\in \{1,\dots n_e-1\}\hfill \\ 1-\sum _{k=1}^{n_k}\sum _{j=1}^{n_e-1}{p}_{j,k}^{pc}{\mathbb{1}}_{G_k}(x)\hfill & i=n_e\hfill \end{array}$$(27) lead to a partitioning of the integration domain G into G_k for the adjoints ${\overline{p}}_{i,k}.$(61) $${\overline{p}}_{i,k}={\int }_{G_k}{\overline{\omega }}_i^{pc}(x)-{\overline{\omega }}_{n_e}^{pc}(x)\mathrm{}\text{d}x$$(61)

Figure 2. The abstract computational graph of our model. Straight edges (bottom to top) correspond to the primal function evaluation. Bend edges (top to bottom) correspond to the adjoint evaluation and are labeled with the respective adjoint operators. Dashed boxes highlight the individual model evaluation for every experimental setup $\alpha \in \mathbb{A}.$ Note that multiple adjoint edges entering a node correspond to the accumulation of the respective adjoint operators. This is not explicitly addressed in the derivation in section 3.2.

Figure 3. The absorption path $l(x,x_d)=\left\{y\right|d(x,x_d,y)=0\}$ from x toward the detector position x_d (solid line) and the path used for the adjoint $l(x,x_{d^{*}})=\left\{y\right|d(y,x_d,x)=0\}$ (dotted line). G is the computational domain.

This completes the derivation of adjoint formulations for all parts of the model presented in Sec. 2. Note that for nonlinear operators, the primary result is required to calculate the adjoint, and therefore must be saved during the forward execution or recalculated. This is analogous to discrete adjoint tools. Furthermore, all derived operators are modular and only need to be assembled according to the computational graph. When replacing parts of the model, only the corresponding adjoint operator needs to be adjusted. In the examples in Sections 4.3 and 4.4, we will use this modularity and exchange the parametrizations $\mathcal{P}.$

3.3. Computational effort

The motivation for deriving and implementing the adjoint formulation of our model is that of the efficiency advantage over finite differences and the tangent mode of AD. Our implementation inherits the runtime complexity of the adjoint mode, as the runtime of a Jacobian calculation scales in the number of output variables of the primal. The price to pay for faster execution is the higher memory requirement since all primal variables have to be stored during forward execution or recalculated later.

Leveraging the convenience of automated AD tools to derive the discrete adjoint code of our implementation was not possible with current tools in julia. However, our approach somewhat mimics the development approach that tools, e.g., Zygote.jl or ChainRulesCore.jl take. Instead of implementing only tangent and adjoint versions of low-level operations, they aim for efficient implementations of high-level methods. We adopt this approach and derive an adjoint version in the continuous framework and implement discretized versions of the operators in code.

The high-level view offers further insight, and one more possibility to reduce the computational effort for the gradient. Note that the naive adjoint implementation requires $\left|\mathbb{A}\right|$ “solves” of the adjoint PDE (same amount as the forward). We separate the multi-index into $\alpha =({\alpha }_L,{\alpha }_B)\in {\mathbb{A}}_L\times {\mathbb{A}}_B$ where α_L describes the k-ratio line (material, transition) and α_B the beam setup (position, energy). We revisit Equations (49) $${\overline{u}}_{\alpha }^{l,k}(x,ϵ)=\{\begin{array}{cc}{\sigma }_{\alpha }(ϵ){\overline{\mathcal{I}}}_{\alpha }(x)\hfill & l,k=0\hfill \\ 0\hfill & \text{else}\hfill \end{array}$$(49) (44)(44) $${\overline{\mathcal{N}}}_{\alpha }(x)={\mathcal{A}}_{\alpha }(x){\mathcal{I}}_{\alpha }(x){\overline{I}}_{\alpha }\hspace{1em}\text{and}\hspace{1em}{\overline{\mathcal{I}}}_{\alpha }(x)={\mathcal{N}}_{\alpha }(x){\mathcal{A}}_{\alpha }(x){\overline{I}}_{\alpha }$$(44) (49)(49) $${\overline{u}}_{\alpha }^{l,k}(x,ϵ)=\{\begin{array}{cc}{\sigma }_{\alpha }(ϵ){\overline{\mathcal{I}}}_{\alpha }(x)\hfill & l,k=0\hfill \\ 0\hfill & \text{else}\hfill \end{array}$$(49) and (51)(51) $$\underset{={\frac{\partial F_{\alpha }}{\partial u_{\alpha }}}^{*}{\overline{F}}_{\alpha }}{\underbrace{S\left[\rho \right](x,ϵ){\partial }_{ϵ}{\overline{F}}_{\alpha }(x,ϵ)-\sum _{d=1}^3{A}_{(d)}{\partial }_{x_i}{\overline{F}}_{\alpha }(x,ϵ)-Q\left[\rho \right](x,ϵ){\overline{F}}_{\alpha }(x,ϵ)}}={\overline{u}}_{\alpha }(x,ϵ)$$(51) . We realize that the only dependence of all equations from α_B is due to the scalar ${\overline{I}}_{\alpha }.$ Instead of computing ${\overline{\mathcal{I}}}_{\alpha },{\overline{u}}_{\alpha }$ and ${\overline{F}}_{\alpha },$ we consider ${\tilde{\mathcal{I}}}_{{\alpha }_L}=\frac{{\overline{\mathcal{I}}}_{\alpha }}{{\overline{I}}_{\alpha }},{\tilde{u}}_{{\alpha }_L}=\frac{{\overline{u}}_{\alpha }}{{\overline{I}}_{\alpha }},{\tilde{F}}_{{\alpha }_L}=\frac{{\overline{F}}_{\alpha }}{{\overline{I}}_{\alpha }}$ and eliminate thereby the dependence from the beam setup α_B. The operators for ${\tilde{\mathcal{I}}}_{{\alpha }_L},{\tilde{u}}_{{\alpha }_L}$ and ${\tilde{F}}_{{\alpha }_L}$ do not depend on α_B. The accumulation of the adjoint ${\overline{\rho }}_i$ (cf. Equation (59)(59) $${\overline{\rho }}_i(x)={\int }_E{\overline{F}}_{\alpha }^T(x,ϵ)({\partial }_{ϵ}(S_i(ϵ)u_{\alpha }(x,ϵ))+Q_i(ϵ)u_{\alpha }(x,ϵ))\mathrm{}\text{d}ϵ.$$(59) , where we sum over all incoming edges in the computational graph 2), can then be rewritten to(62) $$\begin{array}{cc}{\overline{\rho }}_i(x)\hfill & =\sum _{\alpha \in {\mathbb{A}}_L\times {\mathbb{A}}_B}{\int }_E{\overline{F}}_{\alpha }^T(x,ϵ)(-{\partial }_{ϵ}(S_i(ϵ)u_{\alpha }(x,ϵ))-Q_i(ϵ)u_{\alpha }(x,ϵ))\mathrm{}\text{d}ϵ\hfill \\ \hfill & =\sum _{{\alpha }_L\in {\mathbb{A}}_L}{\tilde{F}}_{{\alpha }_L}^T(x,ϵ){\int }_E\sum _{{\alpha }_B\in {\mathbb{A}}_B}{\overline{I}}_{\alpha }(-{\partial }_{ϵ}(S_i(ϵ)u_{\alpha }(x,ϵ))-Q_i(ϵ)u_{\alpha }(x,ϵ))\mathrm{}\text{d}ϵ.\hfill \end{array}$$(62)

Instead of solving the adjoint PDE for every ${\overline{F}}_{\alpha }$ it is sufficient to compute ${\tilde{F}}_{{\alpha }_L}$ for every α_L. The number of expensive adjoint PDE solves that are required is reduced from $\left|\mathbb{A}\right|$ to $\left|{\mathbb{A}}_L\right|.$

4. Numerical examples

Preliminary results for the examples presented in section 4.3 and 4.4 can be found in (Achuda et al. 2023).

4.1. Implementation and discretization

Conceptually the discretization of the P_N-equation (17)(17) $$\underset{=F_{\alpha }(u_{\alpha },\rho )}{\underbrace{-{\partial }_{ϵ}(S[\rho ](x,ϵ)u_{\alpha }(x,ϵ))+\sum _{d=1}^3{A}_{(d)}{\partial }_{x_d}{u}_{\alpha }(x,ϵ)-Q\left[\rho \right](x,ϵ)u_{\alpha }(x,ϵ)}}=0\hspace{1em}\forall x\in G,ϵ\in E$$(17) and the adjoint P_N-equation (51)(51) $$\underset{={\frac{\partial F_{\alpha }}{\partial u_{\alpha }}}^{*}{\overline{F}}_{\alpha }}{\underbrace{S\left[\rho \right](x,ϵ){\partial }_{ϵ}{\overline{F}}_{\alpha }(x,ϵ)-\sum _{d=1}^3{A}_{(d)}{\partial }_{x_i}{\overline{F}}_{\alpha }(x,ϵ)-Q\left[\rho \right](x,ϵ){\overline{F}}_{\alpha }(x,ϵ)}}={\overline{u}}_{\alpha }(x,ϵ)$$(51) are the most challenging. For the numerical computation of the moments $u_{\alpha }$ and ${\overline{F}}_{\alpha }$ we implemented a 3D staggered grid finite-difference method based on the method StaRMAP presented in (Seibold and Frank 2014). StaRMAP is specifically designed to exploit the structure of P_N-equations. Individual sets of moments $u_{\alpha }$ are discretized on the staggered grids in such a way, that sparsity of the transport matrices $A_{(d)}$ decouples the system and enables the use of second-order central differences that approximate the spatial derivative on the respective displaced grid. Second-order integration in energy is achieved using a splitting of the moments and half steps for each part of the solution.

The boundary conditions in (22)(22) $$u_{\alpha }^{[o,d]}=L_{(d)}{\widehat{A}}_{(d)}{u}_{\alpha }^{[e,d]}+g_{\alpha }$$(22) and (53)(53) $${\overline{F}}_{\alpha }^{[o,d]}=-L_{(d)}{\widehat{A}}_{(d)}{\overline{F}}_{\alpha }^{[e,n]}$$(53) are used to compute moments on ghost nodes, located outside the computational domain G. Values of the moments for ghost nodes are set, such that the boundary condition holds for interpolated moments on the boundary of the domain $\partial G.$ For details see Bünger (2021).

For integrations in energy and space as well as the line integral for the absorption we use the trapezoidal rule.

4.2. Comparison: sensitivities of a 1D material with finite differences

We consider a binary material consisting of copper Cu and nickel Ni and compare sensitivities $\frac{\partial I_{\alpha }}{\partial p}$ computed using our implementation (resp. $\overline{p}$ where ${\overline{I}}_{\alpha }$ is seeded with unit vectors) with sensitivities computed by central finite differences. By means of this comparison of our implementation to finite differences we address several questions at the same time:

• By introducing a second parametrization, we exemplarily show the modularity of our implementation.

• We provide insight into the nature of the model. The sensitivity $\frac{\partial I_{\alpha }}{\partial p}$ is the linear approximation of intensities with respect to model parameters, hence a large sensitivity relates to a strong dependency and a small sensitivity to a weak dependency. This is particularly interesting for the inverse problem because it also unveils possible difficulties in the reconstruction.

• Agreement with finite differences validates our implementation.

• We demonstrate the claimed superiority of our method compared to finite differences in terms of computation time.

We introduce an alternative parametrization, where ${\Lambda }_k(x),k\in \{1,\dots n_k\}$ are triangular functions covering the domain G, in particular ${\sum }_{k=1}^{n_k}{\Lambda }_k(x)=1 \forall x\in G.$ The partition of unity ${\sum }_{i=1}^{n_e}{\omega }_i(x)=1$ is (again) guaranteed by $p_{i,k}^{\text{pl}}\ge 0$ and ${\sum }_{j=1}^{n_e-1}{p}_{j,k}^{\text{pl}}\le 1.$ We refer to the parametrization as piecewise-linear.(63) $${\omega }_i^{\text{pl}}(x)=\{\begin{array}{cc}\sum _{k=1}^{n_k}{p}_{i,k}^{\text{pl}}{\Lambda }_k(x)\hfill & i\in \{1,\dots n_e-1\}\hfill \\ 1-\sum _{k=1}^{n_k}\sum _{j=1}^{n_e-1}{p}_{j,k}^{\text{pl}}{\Lambda }_k(x)\hfill & i=n_e\hfill \end{array}$$(63)

The adjoint of the piecewise-linear parametrization ${\overline{p}}^{\text{pl}}$ is given by(64) $${\overline{p}}_{i,k}^{\text{pl}}={\int }_G{\Lambda }_k(x)({\overline{\omega }}_i^{\text{pl}}(x)-{\overline{\omega }}_{n_e}^{\text{pl}}(x))\mathrm{}\text{d}x.$$(64)

Modularity is set out by the fact, that only the parametrization and its adjoint are implemented, the rest of the code remains unchanged.

Now both parametrizations, the piecewise-constant parametrization (with n_k = 20 parameters) and the piecewise-linear parametrization (also with n_k = 20 parameters) are employed to compute sensitivities $\frac{\partial I_{\alpha }}{\partial p^{\text{pc}}}$ resp. $\frac{\partial I_{\alpha }}{\partial p^{\text{pl}}}.$ The physical setup and the settings of the P_N approximation for this example can be found in Table 1. Parameter values for the parametrizations are sampled independently from a uniform distribution in $[0,1].$

Table 1. Model settings used for the comparison of 1D sensitivities.

Domain G	$(-500 \text{nm},0 \text{nm})$ nz = 100
Energies $[{ϵ}_{\text{cut}},{ϵ}_{\min }]$	$(0.5 \text{keV},16 \text{keV})$
Approximation order PN	N = 9
Beam ${\mu }_{ϵ,\text{beam}},{\Sigma }_{ϵ,\text{beam}}$	$15 \text{keV},0.1 \text{keV}$
Elements	Cu and Ni, ${\omega }_i(x)=0.5$

We visualize the sensitivities by plotting ${\omega }_1^{\text{pc}}(x){|}_{p=\frac{\partial I_{\alpha }}{\partial p^{\text{pc}}}}$ and ${\omega }_1^{\text{pl}}(x){|}_{p=\frac{\partial I_{\alpha }}{\partial p^{\text{pl}}}}$ in Figure 4. Using the parametrizations $\omega (x)$ is gentle abuse but since the parameters in ${\omega }^{\text{pc}}$ and ${\omega }^{\text{pl}}$ are related to space points, an interpretation of the plot is possible. Increasing a parameter, increases the mass concentration of Cu and simultaneously decreases the mass concentration of Ni, hence the positive values for Cu x-rays and the negative value for Ni x-rays. In depth the sensitivity decreases because with greater depth the generation of x-rays decreases and the absorption increases.

Figure 4. A visualization of the sensitivities computed using our method and finite differences. We abuse the parametrization function of each parametrization and plot ${\omega }_1(x){|}_{p=\frac{\partial I_{\alpha }}{\partial p}}.$ The curves computed using our method (solid lines) coincide with the curves computed using finite differences (black dashed lines). A quantitative comparison of sensitivity values is given in Table 2.

Figure 4 also compares sensitivities computed by our method (colored, solid lines) with the sensitivities computed by central finite differences (black, dashed lines). The curves coincide. Additionally, we compare relative errors (65) $$\underset{k}{\max }\left|\frac{1}{\frac{\partial I_{\alpha }}{\partial p_k}}\left(\frac{\partial I_{\alpha }}{\partial p_k}-\frac{I_{\alpha }(p+h{e}_k)-I_{\alpha }(p-h{e}_k)}{2h}\right)\right|$$(65) of sensitivities in Table 2. Derivatives computed by algorithmic differentiation are exact (AD computes the exact derivative of the numerical approximation up to machine precision) (Naumann 2011). Our method as well as finite differences only provide approximations of the derivative. For our method, the underlying approximation of the adjoint state variable introduces errors, for finite differences the approximation error originates from the perturbation h.

Table 2. Relative errors in the sensitivities/derivatives of intensities for multiple x-rays α_L and parametrization methods $\omega (x).$ Sensitivities/derivatives are computed using our implementation and using a central finite difference approximation of order 2.

αL $\omega (x)$	$(Cu,K-L_2)$	$(Cu,K-L_2)$	$(Cu,K-L_2)$	$(Cu,K-L_2)$
Piecewise-constant	0.0015	0.0005	0.0013	0.0003
Piecewise-linear	0.0030	0.0007	0.0008	0.0011

In Table 3, the superiority of our method in terms of computation time compared to finite differences becomes apparent. We compare the runtime of a central finite difference method (2nd order) with our implementation for two different numbers of parameters (n_p = 20 and n_p = 40). The runtime of the forward model is $\approx 0.5s,$ and the computation of the derivative with finite differences clearly scales with the number of parameters n_p. The runtime of our implementation scales in the number of outputs (here $\left|\mathbb{A}\right|=4$) and clearly does not scale with n_p. Note the minor overhead because of the backward execution.

Table 3. Runtimes of the computation of sensitivities/derivatives using our implementation and a central finite difference approximation of order 2 (FiniteDiff.jl).

Method $\omega (x)$	AD (np = 20)	FiniteDiff (np = 20)	AD (np = 40)	FiniteDiff (np = 40)
Piecewise-constant	$2.29\mathrm{s}$	$18.64\mathrm{s}$	$2.33\mathrm{s}$	$39.29\mathrm{s}$
Piecewise-linear	$2.34\mathrm{s}$	$18.92\mathrm{s}$	2.43s	$38.80\mathrm{s}$

Table 4. Model settings used for the reconstruction of the 1D interfaces.

Domain $G$	$(-500 \text{nm},0 \text{nm})$ nz = 100
Energies $[{ϵ}_{\text{cut}},{ϵ}_{\min }]$	$(7.0 \text{keV},15.5 \text{keV})$
Approximation order PN	N = 9
Beam $ϵ,{\sigma }_{ϵ}$	$\{9,10.5,12,13.5,15\}\text{keV},0.1 \text{keV}$
Elements/x-rays	$(Ni,K-L_2),(Ni,K-L_3),(Fe,K-L_2)$ and $(Fe,K-L_3)$

4.3. 1D reconstruction of a sharp and a diffusive interface

We present the reconstruction of a coated material consisting of iron Fe and nickel Ni with a sharp (discontinuous) and a diffusive (continuous) interface between the layers.

• We show the necessity of choosing a suitable parametrization.

• We compare the previously defined parametrizations ((27)(27) $${\omega }_i^{pc}(x)=\{\begin{array}{cc}\sum _{k=1}^{n_k}{p}_{i,k}^{pc}{\mathbb{1}}_{G_k}(x)\hfill & i\in \{1,\dots n_e-1\}\hfill \\ 1-\sum _{k=1}^{n_k}\sum _{j=1}^{n_e-1}{p}_{j,k}^{pc}{\mathbb{1}}_{G_k}(x)\hfill & i=n_e\hfill \end{array}$$(27) and (63)(63) $${\omega }_i^{\text{pl}}(x)=\{\begin{array}{cc}\sum _{k=1}^{n_k}{p}_{i,k}^{\text{pl}}{\Lambda }_k(x)\hfill & i\in \{1,\dots n_e-1\}\hfill \\ 1-\sum _{k=1}^{n_k}\sum _{j=1}^{n_e-1}{p}_{j,k}^{\text{pl}}{\Lambda }_k(x)\hfill & i=n_e\hfill \end{array}$$(63) ) to an additional non-linear parametrization (66)

We introduce another parametrization that is based on the non-linear structure of a neural network (with x as the input and ${\omega }^{\text{nl}}$ as the output. The activation function for one hidden layer is $\tanh ()$ and $\text{sig}()$ for the output layer to enforce $0\le {\omega }_i^{\text{nl}}\le 1$). The parametrization uses n_p = 10 parameters, which we call $p^{\text{nl}}={(a_1,a_2,a_3,b_1,b_2,b_3,{\widehat{a}}_1,{\widehat{a}}_2,{\widehat{a}}_3,\widehat{b})}^T.$(66a) $$v_i(x)=\tanh (\underset{w_i(x)}{\underbrace{a_{i}x+b_i}})\hspace{1em}\forall i=1,2,3$$(66a) (66b) $${\omega }_1^{\text{nl}}(x)=\text{sig}(\underset{=\widehat{w}(x)}{\underbrace{\sum _{i=1}^3{\widehat{a}}_i{v}_i(x)+\widehat{b}}})$$(66b) (66c) $${\omega }_2^{\text{nl}}(x)=1-{\omega }_1^{\text{nl}}(x)$$(66c)

The adjoints for the non-linear parametrization are given in the following.(67a) $${\widehat{\overline{a}}}_i={\int }_G\text{sig}\prime (\widehat{w}(x))v_i(x)({\overline{\omega }}_1^{\text{nl}}(x)-{\overline{\omega }}_2^{\text{nl}}(x))\mathrm{}\text{d}x$$(67a) (67b) $$\widehat{\overline{b}}={\int }_G\text{sig}\prime (\widehat{w}(x))({\overline{\omega }}_1^{\text{nl}}(x)-{\overline{\omega }}_2^{\text{nl}}(x))\mathrm{}\text{d}x$$(67b) (67c) $${\overline{v}}_i(x)=\text{sig}\prime (\widehat{w}(x)){\widehat{a}}_i({\overline{\omega }}_1^{\text{nl}}(x)-{\overline{\omega }}_2^{\text{nl}}(x))$$(67c) (67d) $${\overline{a}}_i={\int }_G\tanh \prime (w_i(x))x{\overline{v}}_i(x)\mathrm{}\text{d}x$$(67d) (67e) $${\overline{b}}_i={\int }_G\tanh \prime (w_i(x)){\overline{v}}_i(x)\mathrm{}\text{d}x$$(67e)

For a fair comparison, the number of parameters in the other parametrizations is also chosen to be 10. For the binary material the scalar ${\omega }_1(x)$ completely describes the material, because ${\omega }_2(x)=1-{\omega }_1(x).$

The artificial measurements for the reconstruction are k-ratios from multiple different beam energies between $9 \text{keV}$ and $15 \text{keV}.$ To compute the measurements we assume the reference material to be given. Additional P_N model settings and the considered x-rays are tabulated in Table 4. We choose the BFGS method implemented in Optim.jl (Mogensen and Riseth 2018) as optimization method (with additional box-constraints for the piecewise-constant and the piecewise-linear parametrization, and no constraints for the non-linear parametrization).

In Figure 5, the reconstructed density of the material with a sharp interface is visualized using the different parametrizations. From initial configurations (piecewise-constant, piecewise-linear: 0.5 Fe and 0.5 Ni; non-linear: random) the parametrizations iterate toward the reference density (black line) during the optimization. We visualize the initial, first, fifth iteration and the 100th iteration. The inability of the piecewise-linear parametrization to represent discontinuities becomes apparent. However, the piecewise-constant parametrization only approximates perfectly because one of the subdomain interfaces matches the interface of the reference material. On the other hand, the non-linear parametrization is flexible enough to identify the location of the interface and to approximate the discontinuity.

Figure 5. The reconstructed density of a material with a sharp interface between an iron Fe layer covering an Ni substrate. We use the piecewise-constant, the piecewise-linear and the non-linear parametrization for the reconstruction. For the piecewise-constant and the piecewise-linear parametrization, the geometry (interfaces) are visualized by gray vertical lines. The black line shows the reference density. All parametrizations converge. The piecewise-constant parametrization has a clear advantage because the interface aligns with an interface of the parametrization. The piecewise-linear parametrization cannot approximate the discontinuous interface, hence it converges to presented smoothed representation. The non-linear parametrization is flexible enough to identify the location of the interface and also approximates the discontinuity using a strong gradient.

In Figure 6, the reconstructed density of the material with a diffusive interface is visualized. Except for the reference material, exactly the same settings as for the reconstruction of the sharp interface are used. Again, the parametrizations iterate toward the reference density (black line). For the diffusive interface, none of the parametrizations can reconstruct the interface perfectly, but the piecewise-linear and the non-linear parametrizations perform better than the piecewise-constant parametrization.

Figure 6. The reconstructed density of a material with a diffusive interface between an iron Fe layer covering an Ni substrate. We use the piecewise-constant, the piecewise-linear and the non-linear parametrization for the reconstruction. For the piecewise-constant and the piecewise-linear parametrization, the geometry (interfaces) are visualized by gray vertical lines. The black line shows the reference density. All parametrizations converge. None of the parametrizations can represent the interface perfectly. The piecewise-linear and the non-linear parametrization perform better for this example, because of their flexibility to approximate continuous functions.

In Figure 7, we visualize the normalized error, the value of the objective function $\left|\right|k(p)-k^{\hspace{0.17em}\exp \hspace{0.17em}}|{|}^2,$ the error of the mass concentrations $\left|\right|\rho (x)-{\rho }^{\text{true}}|{|}^2$ and the error of additional k-ratio measurements $\left|\right|k_{\text{+}}(p)-k_{\text{+}}^{\hspace{0.17em}\exp \hspace{0.17em}}|{|}^2$ (with beam energies $11 \text{keV}$ and $14 \text{keV}$) during the iterations of the optimization. The kinks in the error of the objective function are an artifact of the Optim.jl implementation. Their default implementation is based on the combination of a line search method (Hager and Zhang 2005) and the BFGS algorithm presented in (Nocedal and Wright 2006).

Figure 7. The value of the objective function and additional reconstruction quality measures during the iterations of the optimization for all parametrizations. The objective function $\left|\right|k(p)-k^{\hspace{0.17em}\exp \hspace{0.17em}}|{|}^2$ is visualized as a solid line. Also, the error in the mass concentrations $\left|\right|\rho (x)-{\rho }^{\text{true}}(x)|{|}^2$ (dashed line) and the error of k-ratios $\left|\right|k_{\text{+}}(p)-k_{\text{+}}^{\hspace{0.17em}\exp \hspace{0.17em}}|{|}^2$ (dotted line) from additional beam energies ($11\text{keV}$ and $14\text{keV}$) is plotted.

Differences in the performance of the different parametrizations are clearly visible. The error in the mass concentrations for the piecewise-constant parametrization for the sharp interface is by far the smallest, due to the possible perfect reconstruction. Only the non-linear parametrization performs similarly well for both examples.

Comparing the error in mass concentrations for piecewise-constant and piecewise-linear parametrization for the sharp and the diffusive interface, they vary. A correlation between the error in mass concentrations and the error of additional k-ratios would be beneficial but is hard to obtain visually. Ideally, both would behave similarly, because then we could propose the error in additional k-ratios as a suitable reconstruction quality measure. However, the evaluation of reconstruction quality measures is beyond the scope of this work and is deferred to further research.

4.4. 2D reconstruction of an ellipsoidal inclusion

We describe the use case of a reconstruction of an ellipsoidal copper Cu inclusion in an iron Fe substrate from k-ratios obtained from a $12 \text{keV}$ electron beam and multiple beam positions. In Figure 8, we plot k-ratios obtained from an artificial line-scan of the inclusion. Only the ten k-ratios retrieved from five beam positions ${\mu }_y=\{-50,-25,0,25,50\} \text{nm},$ which are marked by black crosses, are used for the reconstruction. To compute the k-ratios, the reference material is assumed to be given. We visualize the total density ${\rho }_{\text{tot}}$ of the reference material in Figure 9 and additionally show the ionization distributions of $(Cu,K-L_2)$ and $(Fe,K-L_2)$ that define the size of the interaction volume. The ellipsoidal Cu inclusion is clearly smaller than the interaction volumes. Additional settings of the k-ratio model are shown in Table 5.

Figure 8. The k-ratio profile of the elliptical Cu inclusion in the Fe substrate computed using a $12 \text{keV}$ electron beam. Dashed lines show the k-ratio profiles of $(Cu,K-L_2)$ (blue) and $(Fe,K-L_2)$ (orange) with a high beam resolution. Black markers show the k-ratios which are used for the reconstruction. While not directly visible, the shape and height of the k-ratio curves encode information about the shape and location of the inclusion and the structure of the interface between inclusion and substrate.

Figure 9. The total density and two ionization distribution fields of the reference material. The ionization distribution curves determine the size of the interaction volume. The size of the ellipsoidal Cu inclusion is smaller than the interaction volume.

Table 5. Model settings used for the 2D reconstruction of the ellipsoidal Cu inclusion.

Domain G	$(-300 \text{nm},0 \text{nm})\times (-150 \text{nm},150 \text{nm})$ $(n_x,n_y)=(80,80)$
Energies $[{ϵ}_{\text{cut}},{ϵ}_{\min }]$	$(7.0 \text{keV},12.5 \text{keV})$
Approximation order PN	N = 9
Beam ${\mu }_{ϵ},{\sigma }_{ϵ}$	$12 \text{keV},0.2 \text{keV}$
Beam ${\mu }_x,{\sigma }_x$	$(0,\{-50,-25,0,25,50\}) \text{nm},20 \text{nm}$
Beam ${\mu }_{\Omega },$ κ	$[-1,0],$ 10
Elements/x-rays	$(Cu,K-L_2),(Cu,K-L_3),(Fe,K-L_2)$ and $(Fe,K-L_3)$
Standard materials	$(Cu,K-L_{*}):(1.0Cu,0.0Fe),(Fe,K-L_{*}):(0.0Cu,1.0Fe)$

The profile shows an increase in the $(Cu,K-L_2)$ k-ratio for a beam position close to $\approx 25 \text{nm}.$ We assume an elliptical Cu inclusion which we parametrize by $p={({\mu }_1,{\mu }_2,a,b,r,\widehat{a},\widehat{b})}^T,$ where μ₁ and μ₂ specify position, a and b specify scaling factors of the principal axes, r specifies rotation and $\widehat{a}$ and $\widehat{b}$ specify the inclusion interface. The weight fractions $\omega (x)$ are given by the following.(68a) $$v(x)=\frac{{((x_1-{\mu }_1)\hspace{0.17em}\cos \hspace{0.17em}(r)+(x_2-{\mu }_2)\hspace{0.17em}\sin \hspace{0.17em}(r))}^2}{a^2}+\frac{{((x_1-{\mu }_1)\hspace{0.17em}\sin \hspace{0.17em}(r)-(x_2-{\mu }_2)\hspace{0.17em}\cos \hspace{0.17em}(r))}^2}{b^2}$$(68a) (68b) $${\omega }_1(x)=\text{sig}(\underset{\widehat{w}(x)}{\underbrace{\widehat{a}v(x)+\widehat{b}}})$$(68b) (68c) $${\omega }_2(x)=1-{\omega }_1(x)$$(68c)

And their adjoint(69a) $$\widehat{\overline{a}}={\int }_G\text{sig}\prime (\widehat{w}(x))v(x)({\overline{\omega }}_1(x)-{\overline{\omega }}_2(x))\mathrm{}\text{d}x$$(69a) (69b) $$\widehat{\overline{b}}={\int }_G\text{sig}\prime (\widehat{w}(x))({\overline{\omega }}_1(x)-{\overline{\omega }}_2(x))\mathrm{}\text{d}x$$(69b) (69c) $$\overline{v}(x)=\text{sig}\prime (\widehat{w}(x))\widehat{a}({\overline{\omega }}_1(x)-{\overline{\omega }}_2(x))$$(69c) (69d) $$\{{\overline{\mu }}_1,{\overline{\mu }}_2,\overline{a},\overline{b},\overline{r}\}={\int }_G\frac{\partial v(x)}{\partial \{{\mu }_1,{\mu }_2,a,b,r\}}\overline{v}(x)\mathrm{}\text{d}x$$(69d)

For the reconstruction, we use the L-BFGS algorithm implemented in Optim.jl (Mogensen and Riseth 2018). The total density of the material ${\rho }_{\text{tot}}$ during the iteration steps is visualized in Figure 10. From the initial guess (educated guess: higher copper concentration at $x 30 \text{nm},$ see Figure 10(a)) the reconstruction converges to the reference material (see Figure 9(a)). Alongside the total density of the reconstructed material, we visualize the k-ratio profiles of $(Cu,K-L_2)$ and $(Fe,K-L_2)$ in Figure 11. The shapes of the curves quickly coincide with the measured k-ratios (black crosses). At iteration 30 the measured k-ratios are already well approximated, the total density, however, is not. For sufficient reconstruction of the ellipsoidal structure, the error in the k-ratios (the objective function) must be further reduced. In Figure 12, we show the value of the objective function $\left|\right|k(p)-k^{\hspace{0.17em}\exp \hspace{0.17em}}|{|}^2$ and the error in mass concentrations $\left|\right|\rho (p)-{\rho }^{\text{true}}|{|}^2$ during the iterations of the optimization. For the first 50 iterations, the value of the objective function decreases, whereas the error in the mass concentrations remains close to the initial error. Only after the first 50 iterations, the error in the mass concentrations also decreases. In Figure 10, this effect is also observed, where only after iteration 90 the shape of the reference ellipse starts to form.

Figure 10. The reconstructed density ${\rho }_{\text{tot}}$ of the Fe material with a Cu inclusion during the iterations of the L-BFGS optimization.

Figure 11. The k-ratio profile of $(Cu,K-L_2)$ and $(Fe,K-L_2)$ during the iterations of the reconstruction. Black crosses mark the measured k-ratios which are considered in the objective function $\left|\right|k(p)-k^{\hspace{0.17em}\exp \hspace{0.17em}}|{|}^2$ (cf. 8). The visualized iterations are the same as in 10. Note that the k-ratios quickly converge to the measured values, although the shape of the ellipse does not yet coincide with the shape of the ellipse in the true material (9).

Figure 12. Errors of the objective function $\left|\right|k(p)-k^{\hspace{0.17em}\exp \hspace{0.17em}}|{|}^2$ and the error of the mass concentrations $\left|\right|\rho (p)-{\rho }^{\text{true}}|{|}^2$ during the iterations of the L-BFGS optimization. Although the objective function significantly decreases for the first 50 iterations, the error of the mass concentrations remains close to the initial error. The same behavior is noticed in 10 and 11.

5. Conclusion

In this paper, we derive an adjoint approach to compute derivatives for a deterministic k-ratio model to be applied in reconstruction in EPMA. The k-ratio model is based on the P_N-model, a moment expansion of the continuous slowing down approximation of the linear Boltzmann equation that is pioneered for EPMA in Bünger, Richter, and Torrilhon (2022). We formulate the model with a focus on the inverse problem and describe the application of the model within an adjoint algorithmic differentiation framework. This work thus represents a building block in the development of a high-resolution reconstruction method that uses gradient-based optimization techniques. By replacing material parametrizations, experiments can be realized that reconstruct different material parameters and take into account different prior knowledge.

We conclude the paper with a validation of our differentiation method and reconstruction experiments. The experiments show the influence of the material parametrization in the reconstruction and illustrate our idea of an application of our method. The combination of measurements with prior knowledge makes high-resolution reconstruction problems tractable with a reasonable amount of measurement data. At the same time, we demonstrate the extensibility of our implementation to various material structures.

High-resolution imaging, i.e., the reconstruction of material structures that are smaller than the interaction volume, is possible with sufficient data or sufficient material assumptions. But there will always be a balance between the accuracy of the measurement and the model to the accuracy of the reconstruction. Especially for high resolution, the analysis of this relation is key. To achieve a reliable reconstruction result, the quantification and propagation of measurement and model uncertainties are necessary.

Data/code availability

The source code accompanying this article is made available via GitHub: https://github.com/tam724/pnepma.

Disclosure statement

The authors declare no conflict of interest.

Additional information

Funding

The authors gratefully acknowledge funding from the German Research Foundation (DFG) under grant number TO 414/7-1 as well as additional funding contributed under the excellence strategy of the German federal and state governments.

Notes

1 In IUPAC notation (Jenkins et al. 1991), an x-ray level corresponds to an electron configuration that is typically described by the removal of an electron from the configuration of neutral ground state (e.g., K, L₂). A characteristic x-ray is emitted by transitioning from an initial to a final x-ray level, namely the x-ray transition (e.g., $K-L_2$).

2 In the literature the mass attenuation coefficient is commonly denoted using ${(\mu /\rho )}_{\alpha ,i}.$ It depends on the medium i and the x-ray wavelength α. For brevity of notation, we write ${\tau }_{\alpha ,i}.$

3 For x-ray transitions α_L where the initial state is K, the x-ray production cross-section can be replaced by the ionization cross-section of the K shell (k-ratio normalization eliminates the fluorescence yield). For other transitions, e.g. L or M as the final state, the x-ray production cross-section should include the ionization cross-section of all lower shells weighted by transition probabilities.

4 For brevity of notation we divide both quantities by the respective atomic masses.

Appendix A.

Adjoint continuous implicit differentiation

Consider the implicit definition: Given $x\in X,$ find $y\in Y$ such that $F(y=g(x),x)=0.$ We refer to the implicit function with $y=g(x)$ and assume that it is unique at x. Given a (tangent) direction $\dot{x},$ we realize that the directional derivative of F in direction $\dot{x}$ does not change. Employing the chain rule we find the following. (70) $${\frac{\partial }{\partial h}F(g(x+h\dot{x}),x+h\dot{x})|}_{h=0}=\frac{\partial F(y,x)}{\partial y}\left[\frac{\partial g(x)}{\partial x}\left[\dot{x}\right]\right]+\frac{\partial F(y,x)}{\partial x}\left[\dot{x}\right]=0$$(70)

We consider the previous statement in a weak sense, meaning we test the constraint with $\overline{F}\in T,$ where T is an appropriate space.(71) $${〈\overline{F},\frac{\partial F(y,x)}{\partial y}\left[\frac{\partial g(x)}{\partial x}\left[\dot{x}\right]\right]+\frac{\partial F(y,x)}{\partial x}\left[\dot{x}\right]〉}_T=0\hspace{1em}\forall \overline{F}\in T$$(71)

We identify the tangent of y with $\dot{y}=\frac{\partial g}{\partial x}\left[\dot{x}\right]$ and use adjoint operators ${\frac{\partial F}{\partial y}}^{*}[·]$ and ${\frac{\partial F}{\partial x}}^{*}[·]$ of the tangent operators $\frac{\partial F}{\partial y}[·]$ and $\frac{\partial F}{\partial x}[·]$ to derive the following.(72) $$\begin{array}{c}{〈\overline{F},\frac{\partial F(y,x)}{\partial y}\left[\dot{y}\right]〉}_T=-{〈\overline{F},\frac{\partial F(y,x)}{\partial x}\left[\dot{x}\right]〉}_T\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}}\forall \dot{x}\in X\hfill \\ \iff {〈\underset{=\overline{y}}{\underbrace{{\frac{\partial F(y,x)}{\partial y}}^{*}\left[\overline{F}\right]}},\dot{y}〉}_Y\hfill & ={〈\underset{=\overline{x}}{\underbrace{-{\frac{\partial F(y,x)}{\partial x}}^{*}\left[\overline{F}\right]}},\dot{x}〉}_X\hspace{1em}\forall \dot{x}\in X\hfill \end{array}$$(72)

Comparing (72)(72) $$\begin{array}{c}{〈\overline{F},\frac{\partial F(y,x)}{\partial y}\left[\dot{y}\right]〉}_T=-{〈\overline{F},\frac{\partial F(y,x)}{\partial x}\left[\dot{x}\right]〉}_T\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}}\forall \dot{x}\in X\hfill \\ \iff {〈\underset{=\overline{y}}{\underbrace{{\frac{\partial F(y,x)}{\partial y}}^{*}\left[\overline{F}\right]}},\dot{y}〉}_Y\hfill & ={〈\underset{=\overline{x}}{\underbrace{-{\frac{\partial F(y,x)}{\partial x}}^{*}\left[\overline{F}\right]}},\dot{x}〉}_X\hspace{1em}\forall \dot{x}\in X\hfill \end{array}$$(72) with the definition of the adjoint (31)(31) $${〈\overline{y},\frac{\partial f(x)}{\partial x}\left[\dot{x}\right]〉}_Y={〈{\frac{\partial f(x)}{\partial x}}^{*}\left[\overline{y}\right],\dot{x}〉}_X\hspace{1em}\forall \dot{x}\in X,\overline{y}\in Y$$(31) ${〈\overline{y},\dot{y}〉}_Y={〈\overline{x},\dot{x}〉}_X$ we derive the two steps for the implicit adjoint. Given $\overline{y}\in Y,$ if we can find $\overline{F}\in T$ such that (73a)(73a) $${\frac{\partial F(y,x)}{\partial y}}^{*}\left[\overline{F}\right]=\overline{y}$$(73a) holds, the required adjoint $\overline{x}$ is given by (73b)(73b) $$\overline{x}=-{\frac{\partial F(y,x)}{\partial x}}^{*}\left[\overline{F}\right].$$(73b) .(73a) $${\frac{\partial F(y,x)}{\partial y}}^{*}\left[\overline{F}\right]=\overline{y}$$(73a) (73b) $$\overline{x}=-{\frac{\partial F(y,x)}{\partial x}}^{*}\left[\overline{F}\right].$$(73b)

Example

We consider the Poisson-equation in 2D with homogeneous (Dirichlet) boundary conditions. Given a domain $\Omega \subset {\mathbb{R}}^2$ with $\xi \in \Omega$ and $x:\Omega \to \mathbb{R}$ a suitable right-hand side. The primal problem is to find $y:\Omega \to \mathbb{R}$ with $y(\xi )=0 \forall \xi \in \partial \Omega$ such that(74) $$F(y,x)=\Delta y(·)-x(·)=0.$$(74)

To derive the adjoint, we first analyze the left-hand side of (72)(72) $$\begin{array}{c}{〈\overline{F},\frac{\partial F(y,x)}{\partial y}\left[\dot{y}\right]〉}_T=-{〈\overline{F},\frac{\partial F(y,x)}{\partial x}\left[\dot{x}\right]〉}_T\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}}\forall \dot{x}\in X\hfill \\ \iff {〈\underset{=\overline{y}}{\underbrace{{\frac{\partial F(y,x)}{\partial y}}^{*}\left[\overline{F}\right]}},\dot{y}〉}_Y\hfill & ={〈\underset{=\overline{x}}{\underbrace{-{\frac{\partial F(y,x)}{\partial x}}^{*}\left[\overline{F}\right]}},\dot{x}〉}_X\hspace{1em}\forall \dot{x}\in X\hfill \end{array}$$(72) . Using integration by parts twice, we derive the adjoint operator ${\frac{\partial F}{\partial y}}^{*}[·].$(75) $$〈\overline{F},\frac{\partial F(y,x)}{\partial y}\left[\dot{y}\right]〉=〈\overline{F},\Delta \dot{y}〉=\underset{=0}{\underbrace{{\int }_{\partial \Omega }\overline{F}(\xi )\nabla \dot{y}(\xi )n-\dot{y}(\xi )\nabla \overline{F}(\xi )n\mathrm{}\text{d}S}}+〈\underset{={\frac{\partial F(y,x)}{\partial y}}^{*}\left[\overline{F}\right]}{\underbrace{\Delta \overline{F}}},\dot{y}〉$$(75)

 The boundary integral in (75)(75) $$〈\overline{F},\frac{\partial F(y,x)}{\partial y}\left[\dot{y}\right]〉=〈\overline{F},\Delta \dot{y}〉=\underset{=0}{\underbrace{{\int }_{\partial \Omega }\overline{F}(\xi )\nabla \dot{y}(\xi )n-\dot{y}(\xi )\nabla \overline{F}(\xi )n\mathrm{}\text{d}S}}+〈\underset{={\frac{\partial F(y,x)}{\partial y}}^{*}\left[\overline{F}\right]}{\underbrace{\Delta \overline{F}}},\dot{y}〉$$(75) yields boundary conditions for the adjoint $\overline{F}.$ Since $\dot{y}(\xi )=0 \forall \xi \in \partial \Omega ,$ we assign $\overline{F}(\xi )=0 \forall \xi \in \partial \Omega .$ From the right-hand side of (72)(72) $$\begin{array}{c}{〈\overline{F},\frac{\partial F(y,x)}{\partial y}\left[\dot{y}\right]〉}_T=-{〈\overline{F},\frac{\partial F(y,x)}{\partial x}\left[\dot{x}\right]〉}_T\mathrm{\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}}\forall \dot{x}\in X\hfill \\ \iff {〈\underset{=\overline{y}}{\underbrace{{\frac{\partial F(y,x)}{\partial y}}^{*}\left[\overline{F}\right]}},\dot{y}〉}_Y\hfill & ={〈\underset{=\overline{x}}{\underbrace{-{\frac{\partial F(y,x)}{\partial x}}^{*}\left[\overline{F}\right]}},\dot{x}〉}_X\hspace{1em}\forall \dot{x}\in X\hfill \end{array}$$(72) we trivially find(76) $${\frac{\partial F(y,x)}{\partial x}}^{*}\left[\overline{F}\right]=\overline{F}.$$(76)

Computing the adjoint $\overline{x}$ therefore consists of the following steps. Given $\overline{y}(·)$ find $\overline{F}(·)$ with $\overline{F}(\xi )=0 \forall \xi \in \partial \Omega$ such that (77a)(77a) $$\Delta \overline{F}=\overline{y}$$(77a) holds. Then, the required adjoint $\overline{x}(·)$ is given by (77b)(77b) $$\overline{x}=\overline{F}$$(77b) (trivial for this example).(77a) $$\Delta \overline{F}=\overline{y}$$(77a) (77b) $$\overline{x}=\overline{F}$$(77b)