Inversion model for extracting chemically resolved depth profiles across liquid interfaces of various configurations from XPS data: PROPHESY

The PROPHESY framework for obtaining absolute depth profiles from discrete X-ray photoelectron spectra is introduced. It is composed of (i) a model for an X-ray photoelectron spectroscopy experiment accounting for sample geometry, and (ii) an inversion model to reconstruct the concentration profile at the surface of a liquid sample. As a proof of concept, it is applied to simulated data.


Assumption leading to the measurement model
In the work of Ottosson et al. 2010 (Ottosson et al., 2010), the PE signal is defined as a function of electron kinetic energy, but it does not include the dependence of the parameters which is left implicit.However, in this work, we need this important information in order to understand where discrepancies may arise and where approximations may come from.In the work of Dupuy et al. 2021 (Dupuy et al., 2021), the PE signal is modeled with integration over the volume of the sample of the product of the density ρ and the exponential attenuation.The beam profile is left out of the volume integral which implies some (left unknown) assumptions.However, the photon flux is a quantity that varies over the integration volume because the beam profile is not flat (Fedoseenko et al., 2003;Kachel, 2016) Also, the amplitude of the flux is attenuated in the sample.Though, the attenuation of the photon flux can be neglected over the thickness (less than 100 µm) of the sample, e.g. the linear attenuation length of 1 [keV] X-ray in water is around 2.5 mm (Berger, 1998).Therefore, neglecting the attenuation of the photon flux in the sample is a weak assumption for thin samples such as LJ.The variation in space of the differential cross-section density σ χ (ν, K e , θ) in the sample is also left out because the analyzer is sufficiently far from the sample (e.g. a few radii away) that the measurement angle is constant over the sample.
In order to account for the spatial approximations in the photon flux, the alignment parameter α [m −2 ] is introduced as a multiplicative factor (Ozon et al., 2023;Dupuy et al., 2021).The alignment parameter depends on the acquisition setup and should be estimated for every acquisition, e.g. from raw experimental data.
Our presented model accounts for the spatial and spectral variability of the model parameters, i.e. the photoionization differential cross-section density and the photon flux density.Hence, the electron flux J [electron s −1 ] reaching the aperture of the measurement device (not including the device itself) is given by the summation of the contribution of all photon energy over the spectral domain Ω ν , the contribution over the area where the photon beam extends on the sample, and the angle domain Ω θ that cover the aperture of the kinetic analyzer.We write the photo-electron flux as it was modeled in (Ozon et al., 2023) curve M s leading to the surface of the liquid in direction of the analyzer located in P = (x 0 , y 0 , z 0 ) at a distance The emitted photoelectrons are assumed to be traveling in straight lines, and the attenuation length characterizes the losses by elastic or inelastic interaction with the sample.In this formulation, it is implicitly assumed that the scattering properties are uniform across the sample, and that outside the sample (where ρ tot = 0) the photoelectrons do not undergo scattering events.The integration domain Ω V covers the sample and dV is the infinitesimal volume around M .Considering any point M in the liquid and any point P outside of the liquid, the straight line joining M to P can be parameterized as: x(s) y(s) z(s) where the direction angles ω and β are depicted in fig.
(2) and s [m] is the parameter of the curve that represents the signed distance from point M .The angle ω is between the z-axis and the projection M P onto the plane zOx, and β is taken as the angle between the plane zOx and M P .Note that the concentrations ρ and ρ tot in the model are not the same, one represents the density of all species ρ tot (overwhelmingly water for aqueous solutions) and the other ρ is the concentration for a given target orbital in the sample.The spectral integration domain Ω ν covers the support of the spectral density of the light source.The integration domain Ω θ of the emission angle of the PE depends on the sample and the aperture of the kinetic energy analyzer.
In the case of synchrotron light sources, it is possible to assume that the light is monochromatic (Fedoseenko et al., 2003;Kachel, 2016).The monochromaticity is measured relatively to the central frequency value ν k as the ratio of the spread ∆ ν k and ν k .However, the spread ∆ ν k should be compared to the spread of the kinetic energy of the emitted electrons because this is the signal of interest which should not IUCr macros version 2.1.15:2021/03/05 be blurred by the exciting light.We write , and δ the Dirac distribution modeling the spectral density [eV −1 ], then the electron flux is ρ 0 λe(Ke) dτ dθdV (3) In the most general case, the angle θ, and its domain of integration Ω θ , depend on the location in the sample, however, we approximate θ by that of the center of the sample and Ω θ by the apparent angle of the aperture of the analyzer from the center of the sample.It was shown early on that the location of the sample relative to the analyzer (Siegbahn & Siegbahn, 1973) plays an important role in the measured signal.Also, the location in the sample, e.g. the relative angle of the surface and the solid angle of the aperture (Olivieri et al., 2017), affects the amount of signal that can be measured.In the latter, it is shown that the integration domain should be the intersection between the illuminated sample volume with the observation cone.
Rigorously, the signal reaching the aperture of the spectrometer is originating from the volume irradiated by the X-ray source convoluted with the area analyzed by the spectrometer (Guilet et al., 2022).Note that the size of the photon beam also has an effect on the amount of signal coming from the vapor in LJ (Olivieri et al., 2015).
From this assumption, we separate the two integrals Denoting the first integral as α θ σ χ (ν k , K e ), then we have ρ 0 λe(Ke) dτ dV. (5) Finally, we use the alignment parameter α [m −2 ] an average probability density of interaction between the photon beam and the sample (Ozon et al., 2023) to account for approximation of the photon beam profile f by a flat profile with effective photon IUCr macros version 2.1.15:2021/03/05 density αF (ν k ) and the limited angular opening of the analyzer α θ we have For the sake of clarity, the volume integral in (6), which bears the geometry information, will be denoted H(ρ, λ e ).

Discretization
We estimate the different cross-section densities σm,k χ and the profile concentration ρ from a very limited number of acquisitions.Typically, the number of frequencies ν k used for probing the sample is no more than K = 5.Ideally, we should have access to sufficiently many (K = 20) frequencies so that the collection of attenuation lengths the maximum depth at which we want to reconstruct the concentration profile.

Depth and kinetic energy
To limit the challenge, we do not seek solutions in an infinite space of functions, but rather in a finite subspace.We approximate the sought functions by ρ n e n (r) and σm,k the approximation: where The discretization noise ε m,k is the error due to the approximation of the functions ρ and σm,k χ by their piecewise linear approximations.The matrix element for the k th frequency is: Note that H k ,n does not depend on the specific peak m, rather it is a sample-geometry factor.

Attenuation length
In practice, for the measurement of a single spectrum (collection of PE signals centered around a reference kinetic energy, K k e.g.

2
), the variation in attenuation length λ e is rather small (Thürmer et al., 2013;Ottosson et al., 2010) compared with the central value.For instance, for λ e (K k ) = 2 [nm] the variation is of the order of 10 −2 [nm] over a range of kinetic energy of a few eV.Therefore, we assume that it can be approximated by with λ k = λ e (K k ).The Landau notation O (K e − K k ) 2 stands for the higher order terms, i.e. all the terms of order higher than or equal to 2 in this case.Substituting SI(10) in the model SI(9), we obtain which is in essence the outer product of a discretization matrix over the kinetic energy space and that over the depth space.Therefore, the rank of H k is 1, making it impossible to recover depth information from one spectrum.The coefficient T k c k is the discretization of the kernel function ϕ k onto the basis functions (f ) 1 L , and does not depend on the index because the kernel and basis functions have the same form for all index .From here, eq.SI(8) becomes where ε m,k now also includes the approximation errors described by (ι k ,n ) 1 L .

Photoionization cross-section density estimation
We assume that the background of the spectra have been estimated and removed from the spectral data, for instance using the package SPANCF (Kukk et al., 2001;Kukk et al., 2005) or any other algorithm for background removal (Baek et al., 2015).We propose an alternative to current C1s peak fitting methods (Kukk et al., 2001;Kukk et al., 2005;Major et al., 2020) that does not rely on parametric peaks to be fitted.is denoised and can be used as input for a regular peak fitting routine, e.g.SPANCF (Kukk et al., 2001;Kukk et al., 2005).The method is written in terms of C1s, but is not limited to this orbital or element.
For each kinetic energy in the k th spectrum, by combining eq.16, eq.22 and eq.SI(12), we have By definition, adding up the contribution of each kinetic energy leads to the total count I k : where the discretization step δ Ke [eV] is assumed constant.The discrete sum over is identified as the Riemann quadrature, hence 1 = . By construction, the approximation error is at most linear in the kinetic energy step, , and from the definition we have δ Ke σ T k and δ Ke < < ∆ Ke .Hence, the total count becomes: Therefore, the spectrum acquisition model can be simplified using the total count I k so that We assume that the discretization error are negligible compared with the noise level, and that the measurement noise terms ε k ∼ N (0, (σ k ) 2 ) are mutually independent.
The total count I k and its variance σ 2 k can then be estimated from the measurements IUCr macros version 2.1.15:2021/03/05 From here, the coefficients of the photoionization cross-section density can be estimated by solving: with I L the L-order identity matrix.We write the more compact form: We assume that all the perturbation elements in ε are Gaussian distributed with covariance matrix Γ I that is block diagonal.The first block is diagonal with entries can not conclude that the maximum of the posterior σk C1s |I k , y is also its expectation.

Estimating the uncertainty in σk
C1s |I k , y due to the uncertainty in the model, i.e. the variability of the value I k , may be formulated as the covariance of the estimation these quantities can be computed by sampling the model, e.g. 4. Detail of the optimization problem eq. ( 26) The optimization problem defined in eq. ( 26) of the paper relies on two probability densities, i.e. the measurement likelihood P(y|A m , ρ) and the a priori P(A m , ρ).Here we detail the meaning and definition of each term.As stated in sec.3 of the paper, the measurement noise is approximated by a Gaussian distribution because the number of counts is assumed to be greater than 30, making the Gaussian approximation to the Poisson distribution acceptable.We chose this approximation to simplify the implementation of the optimization algorithm, however, it is not a critical assumption and can be modified in a straightforward manner.Here, we have the approximated IUCr macros version 2.1.15:2021/03/05 likelihood Additionally, the noise is assumed to be independent because each measurement is acquired from different experimental setups (e.g.different photon energy) and different samples (the piece of the sample used for acquisition is not the same for different photon energy even though the bulk solution is the same).Therefore, the covariance matrix Γ is diagonal with entries where the variance (σ k ) 2 is defined in eq.SI (19).Note that accounting for the Poisson distributed noise implies changing P(y|A m , ρ) in eq.SI(24), with expected values I k (see eq. SI( 19)).Maximizing the likelihood is an under-determined problem because an infinity of solutions ρ lead to y = A m ρ, hence, we need to constrain the space of possible solutions.The process of constraining the possible solutions is known as regularizing and often relies on a priori (Leong et al., 2023).The second probability density in the optimization problem is the a priori P(A m , ρ), which reflects the knowledge of the state ρ without data, and the uncertainty in the measurement operator.The state ρ and the measurement operator A m are stochastically independent, hence P(A m , ρ) = P(A m )P(ρ).The state a priori P(ρ) represents the probability density of the state ρ.It is interpreted as the plausibility of concentration profiles.This term does not involve nor requires the knowledge of the ground truth.Instead, it represents the properties we expect from a concentration profile.Here, we assume that ρ is not chaotic and is rather smoothly varying with depth (human bias).This is reasonable at the scale/granularity the sample is observed, i.e. averaged over the dimension other than the depth.To represent mathematically this assumption, we resort to the second order difference operator IUCr macros version 2.1.15:2021/03/05 an (N −2)×N matrix which is closely related to the second-order derivative.The choice of difference operator for regularizing optimization problem is ubiquitous and well established in inverse problem (Stolzenburg et al., 2022;Nicholls et al., 2012;Rudin et al., 1992;Twomey, 1963).From here, we write with y D the (N − 2)-vector of expected values of second order differences, and Γ D its covariance matrix.Numerically, we choose y D as the vector whose entries are all 0, which implies that the expected profile is linear (at least piecewise linear).The covariance matrix plays the role of the moderator of the linear-profile assumption.The diagonal of Γ D expresses by how much the profile can deviate from linearity, and the off-diagonal elements represent the correlation between the values of the second-order differences at different depth; it is a control over the smoothness of the second-order differences.Formally, the entries of Γ D are where σ D [m −3 ] controls the amplitude of the second order difference, the correlation length δ D [nm] controls its smoothness, and the ratio N K ensures scalability.The values for σ D and δ D are semi-arbitrary, they are chosen in the ballpark of acceptable amplitude of the second order difference and arbitrarily determined so that the reconstruction is acceptable.The choice for σ D could be automated with a criteria such as the L-curve (Stolzenburg et al., 2022).
For the optimization problem eq. ( 26), the probability density P (A m ) does not play role directly because the optimization is against the concentration profile ρ.However, for the quantification of uncertainty in the reconstruction due to the uncertainty in the model, this probability is central.From the peak area model eq.SI( 14), the IUCr macros version 2.1.15:2021/03/05 measurement probability density can be written as where ) gathers all the multiplicative terms of the measurement model.Here, we are interested only in the uncertainty associated with the attenuation lengths (λ k ) 1 k K .As a working assumption, we state that the uncertainties in the attenuation length and those of τ k and ρ tot are stochastically inde- Further, we assume that τ k and ρ tot are perfectly know, therefore, the probability density of the measurement operator reduces to that of the attenuation lengths We are interested in the consequences of the deviation of the available attenuation lengths values relative to the true values (λ 0 k ) 1 k K , therefore, we choose to write the probability density conditionally to the true values The term P((λ 0 k ) 1 k K ) is an a priori term that is uninformative for this work, so, we focus on the conditional P((λ k ) We consider two possible probability densities for describing the uncertainties: 1) independent errors, and 2) global error.The independent errors represent the errors that the granularity of the current models cannot represent.These error are small variations of the order of a few percents (2.5%).On the other hand, the global error reflects uncertainties that could result from a shift in a fit, or differences between two models.These errors are plausibly of the order of 25%.For the independent error, we write for each attenuation length IUCr macros version 2.1.15:2021/03/05 where τ λ is the relative uncertainty level, and (κ k ) 1 k K are independent uniformly distributed random variables.The global error is also modeled with a relative error term, however, it the same value applies across all attenuation lengths A more refined sampling model may be used for investigating the effect of EAL uncertainty.For instance, the parameters of the following semi-empirical attenuation length model from the work of (Emfietzoglou & Nikjoo, 2007) could be sampled: The parameters A, B and C are fitted from experimental datasets, e.g.IXS-D2, therefore, their accuracy is limited.This formula is an approximation of the attenuation length variation with respect to the kinetic energy and does not capture the fine variations of the EAL.Furthermore, the model parameters depends both on the dataset used for fitting as well as the fitting algorithm.Therefore, the parameters A, B, and C bear uncertainty.Using this model, the attenuation length uncertainty can be represented with where the probability distribution for each parameter may be modeled as a uniform distribution centered on the most likely value of the parameter, e.g.A 0 , B 0 , and C 0 .
The attenuation length error investigated in section 4.2.2 is equivalent to studying the effect of the parameter B in eq.SI(32) with A = 0 and C = 0.
Finally, the data probability P(y) is not necessary to compute the MAP estimate.
Actually, in most practical cases, this probability is intractable.We choose to focus on the noise and write P(y) = P(y|y 0 )P(y 0 ) (34 IUCr macros version 2.1.15:2021/03/05 where P(y|y 0 ) is the noise distribution given the non-noisy observation y 0 For the noise marginalization, we use the above with y 0 = A m ρ GT .

Truncated peak area model
Because the optimization problem eq. ( 26) is not numerically advantageous separate/truncate the known from the unknown values of the concentration profile.The augmented model is where ε D is also a zero-mean Gaussian random vector which bears the meaning of tolerance to deviation from the expected values y D .The covariance matrix Γ D represents the strength of the a priori.
The estimation in this form is unstable because we need to discretize the geometry factor H fairly deep in order to capture most of the signal, e.g.20 nm for a maximum penetration depth of 5 nm.However, the signal is only informative over the first layers of the surface, at most 5 nm.Hence, instead of solving for all the entries in ρ, we will focus a the subset and set the entries for the deeper layers to the bulk concentration ρ B , hence the concentration vector can be written as ρ = [ρ 1 , . . .ρ 1 , ρ t S , ρ B , . . .ρ B ] t .Furthermore, the first component is assumed to be known, ρ 1 = 0.By reorganizing and truncating eq.SI(36), we finally obtain: IUCr macros version 2.1.15:2021/03/05 where the truncated matrices are and The data vector y m S is different in size from the data vector of eq.SI(36) since only the subset of rows mapping ρ S in the difference operator D are retained.Now, the values of the data vector are corrected, the covariance too must be adjusted.For the sake of example, we will only account for the uncertainty in the presumed known concentration values and leave the measurement uncertainty aside.Assuming that both ρ 1 and ρ B are random variables with variance σ 2 B , the covariance matrix Γ m S is then given by: Note that accounting for the model uncertainties in the covariance model is not required since the optimization problem eq. ( 26) assumes the model to be known.

Algorithm for computing the inversion
For solving the optimization problem eq. ( 26) for the formulation described in to iteratively converge to an optimal point.The proximal operators can be understood as projections onto convex sets.In the algorithm ALG2, the formulation alternate between projections in the primal space and dual space.In the reformulation 28, the projection in the dual space makes the iterates evolve toward the unconstrained solution, and the projection in the primal space enforces the positivity constraint.Other formulations are possible, there is no unique solution to solve such problem, but the proposed one is sufficiently efficient and intuitive.
and the proximal operators for the convex conjugate F and for G are where each entry of the vector (x) + is the corresponding entry of x if the entry is positive and 0 otherwise.
Beside the apparent complexity of the concepts deployed for the definition of the proximal operator and the convex conjugate, the formula for the operators (and their implementations) are rather simple.The inverse of I+ σ 2 Γ m S only needs to be computed once, and it can take advantage of the eigen decomposition of Γ m S .

Sampling using Metropolis-Hastings
The goal of sampling is to estimate the mean µ ρ|A m ,y and covariance Γ ρ|A m ,y of the posterior distribution P(ρ|A m , y) from samples (ζ i ) 1 i N sample generated by the algorithm SI6.2.Other quantities defined in section 3 also rely on them.The Metropolis- IUCr macros version 2.1.15:2021/03/05 Hastings algorithm has been described several times, notably by the authors it is named after, Metropolis in 1949 (Metropolis & Ulam, 1949) and Hastings in 1970 (Hastings, 1970).MH has also been refined in some cases to improve its performances (Pereyra et al., 2015).It is a sampling procedure that is based on the Monte Carlo method.The samples proposed by a transition mechanism, i.e. a way to jump from one state to another, are always accepted if they increase the probability distribution, and they are not rejected if they induce a decrease of the density.When the proposed sample ρ prop is less favorable, it is accepted with a probability that decreases with the ratio of probability densities The transition mechanism q M H should be designed similar to the target sampled distribution so that it efficiently samples P(ρ|A m , y).Here we choose a symmetric kernel -q M H (ρ curr |ρ prop ) = q M H (ρ prop |ρ curr ) -defined by a Gaussian distribution, and the proposed sample is generated as: with Γ M H a covariance matrix of a correlated process, i.e. the off diagonal coefficients are significantly non-zero.It is defined by the entries where Require: ζ 0 initial concentration profile, number of iterations N sample , q M H communication mechanism IUCr macros version 2.1.15:2021/03/05 the distribution P(ρ|A m , y), we set the initial state ζ 0 of the collection (ζ i ) 1 i N sample to ρ|A m , y, so that the burn in period is very short.and c) and d) τ λ = 2.5%.The green curves represent the GT.In panels a) and c), the profile reconstruction are plotted in blue (ρ|A m , y), orange (ρ|A m , y 0 ) and red (ρ|A m 0 , y) with their respective variabilities as shaded areas.In panels b) and d), the a posteriori (P(ρ|A m , y)) is represented in blue and the marginals in orange (P(ρ|A m , y 0 )) and red (P(ρ|A m 0 , y)).c) and e), the profile reconstruction are plotted in blue (ρ|A m , y), orange (ρ|A m , y 0 ) and red (ρ|A m 0 , y) with their respective variabilities as shaded areas.In the panels b), d), and f), the a posteriori (P(ρ|A m , y)) in represented in blue and the marginals in orange (P(ρ|A m , y 0 )) and red (P(ρ|A m 0 , y)).c) and e), the profile reconstruction are plotted in blue (ρ|A m , y), orange (ρ|A m , y 0 ) and red (ρ|A m 0 , y) with their respective variabilities as shaded areas.In the panels b), d), and f), the a posteriori (P(ρ|A m , y)) in represented in blue and the marginals in orange (P(ρ|A m , y 0 )) and red (P(ρ|A m 0 , y)).c) and e), the profile reconstruction are plotted in blue (ρ|A m , y), orange (ρ|A m , y 0 ) and red (ρ|A m 0 , y) with their respective variabilities as shaded areas.In the panels b), d), and f), the a posteriori (P(ρ|A m , y)) in represented in blue and the marginals in orange (P(ρ|A m , y 0 )) and red (P(ρ|A m 0 , y)). ).In the panels a) and c), the profile reconstruction are plotted in blue (ρ|A m , y), orange (ρ|A m , y 0 ) and red (ρ|A m 0 , y) with their respective variabilities as shaded areas.In the panels b) and d), the a posteriori (P(ρ|A m , y)) in represented in blue and the marginals in orange (P(ρ|A m , y 0 )) and red (P(ρ|A m 0 , y)). ).In the panels a) and c), the profile reconstruction are plotted in blue (ρ|A m , y), orange (ρ|A m , y 0 ) and red (ρ|A m 0 , y) with their respective variabilities as shaded areas.In the panels b) and d), the a posteriori (P(ρ|A m , y)) in represented in blue and the marginals in orange (P(ρ|A m , y 0 )) and red (P(ρ|A m 0 , y)).

Attenuation length sampling domain
f k (ν, M ) [photon m −2 eV −1 ] is the photon spatial and spectral density for the central frequency ν k , and ρ is the concentration expressed in number of molecules per unit volume [m −3 ] and represents the concentration of an orbital, e.g.χ =C1s.The concentration ρ tot [m −3 ] is the molecular concentration of all the species in the sample, e.g.ρ tot = ρ water + ρ SDS for an aqueous solution with SDS.The distance τmax 0 ρtot(Ms(τ )) ρ 0 dτ is the summation of the relative concentration ρtot ρ 0 along the path of the emitted electrons from the emission point M = (x M , y M , z M ) along the parametric IUCr macros version 2.1.15:2021/03/05 ) with (e n ) 1 n N and (f ) 1 L two basis of the linear interpolation function subspace.The coefficient (ρ n ) 1 n N and (σ m,k ) 1 L are the values of the functions ρ and σm,k χ evaluated at the discretization nodes of their respective domain (r n ) 1 n N and (K k e ) 1 L .Using these approximations in the peak model eq.(23), we obtain IUCr macros version 2.1.15:2021/03/05 The second block is the variance of ι c should reflect the Riemann integration error amplitude.The last block corresponding to the terms (ε D i ) 1 i L−2 bears the smoothness strength information.The operator D is the second order difference operator in dimension L, which regularizes the inversion by seeking somewhat smooth solutions.Additionally, the total photoionization cross-section is a positive number and so is its density function σ C1s .The positivity constraint can be enforced in the optimization algorithm used for seeking the solution of the optimization problem σk C1s |I k , y ∈ arg min as VMLM-B (Thiébaut, 2002) can solve the optimization problem SI(22).The a posteriori probability density P(σ k C1s |I k , y) is modelled as the product of three Gaussian functions, one for the PE measurements (likelihood) and two for the a priori (smoothness and integral value).Despite being formed as the product of Gaussian function, the resulting a posteriori is not a Gaussian function.Therefore, one IUCr macros version 2.1.15:2021/03/05

Fig. 1 .
Fig.1.Examples of the photoionization cross-section density estimations from the PE spectra as described in sec.SI3 for a sample probed with five photon energies.For each photon energy, the estimate and the uncertainty are depicted as solid lines and shaded areas respectively.
sec.SI5, we turn to a primal-dual algorithm described in the paper Chambolle & Pock (2011) because of its convergence properties, simplicity for implementing the positivity constraint and the simplicity of implementation.The algorithm described IUCr macros version 2.1.15:2021/03/05 in algorithm SI6.1 -dubbed ALG2 in the paper Chambolle & Pock (2011) -relies on a primal-dual reformulation of the optimization problem and uses the proximal operator In alg.SI6.1, L A is the norm of the operator A, following the predication for the acceleration, we set σ 0 = 1 τ 0 L 2 A and γ = 2 ρ−ρ 0 τ 0 with the distance ρ − ρ 0 being approximated by a rough upper bound of the true (unknown) distance.The upper bound of the distance ρ − ρ 0 is approximated by the distance between two extreme IUCr macros version 2.1.15:2021/03/05 cases with constant concentrations 0 and ρ B , hence we set it to √ N ρ B .The acceleration parameter should depend on the relative strength between the likelihood and the a priori, but we have fixed that value.It does not seem to be strongly dependent within the tested range.In practice, the symmetric positive definite matrix Γ w is taken as the identity matrix and the relative tolerance in the data space r y and in the primal space r x are both set to an arbitrarily small value, e.g. 10 −3 .Using the notations from the paper by Chambolle and Pock 2011 (Chambolle & Pock, 2011), the functions are and δ M H = 5 for N = 100, and the dimension of ρ S is N S − N b − 1 ∈ {10, 14} depending on the profile.This communication mechanism ensures that the proposed state deviates by a small difference σ M H and that the difference is somewhat smooth δ M H . Algorithm SI6.2 is formulated in terms of probability density, but for numerical reasons, it is here implemented in terms of the logarithm of the distributions.For sampling IUCr macros version 2.1.15:2021/03/05 Algorithm SI6.2: Statement of Metropolis-Hasting for sampling the distribution P(ρ|A m , y) using a symmetric transition mechanism.

GTGTGT
Fig. 2. Reconstruction of concentration profile for three different simulated experimental acquisition setups: a) and b) 5 attenuation lengths over the range [1.62, 1.95] nm (N 5 ), d)and d) 5 attenuation lengths over the range [1.28, 5.5] nm (W 5 ), and e) and f) 10 attenuation lengths over the range [1.28, 5.5] nm (W 10 ).The panels a), c), and e) show the estimates and the different variability, with respect to the measurement noise in orange (Γ ρ|A m ,y 0 ), and with respect to the measurement model error in red (Γ ρ|A m 0 y ).The panels b) d) and f) show the conditional posterior probability P(ρ|A m , y) (blue), and the marginals P(ρ|A m , y 0 ) (orange) and P(ρ|A m 0 , y) (red).

GTGTGTGTGT
Fig. 3. Reconstruction of concentration profile for three different simulated experimental acquisition setups: a) and b) 5 attenuation lengths over the range [1.62, 1.95] nm (N 5 ), d)and d) 5 attenuation lengths over the range [1.28, 5.5] nm (W 5 ), and e) and f) 10 attenuation lengths over the range [1.28, 5.5] nm (W 10 ).The panels a),c) and e) show the estimates and the different variability, with respect to the measurement noise in orange (Γ ρ|A m ,y 0 ), and with respect to the measurement model error in red (Γ ρ|A m 0 y ).The panels b) d) and f) show the conditional posterior probability P(ρ|A m , y) (blue), and the marginals P(ρ|A m , y 0 ) (orange) and P(ρ|A m 0 , y) (red).

GTGTGTGTGTGTGT
Fig.6.Profile reconstruction in the case W 10 for two levels of attenuation length uncertainty: a) and b) τ λ = 1%, and c) and d) τ λ = 2.5%.The green curves represent the GT.In panels a) and c), the profile reconstruction are plotted in blue (ρ|A m , y), orange (ρ|A m , y 0 ) and red (ρ|A m 0 , y) with their respective variabilities as shaded areas.In panels b) and d), the a posteriori (P(ρ|A m , y)) is represented in blue, and the marginals in orange (P(ρ|A m , y 0 )) and red (P(ρ|A m 0 , y)).

GTGT
Fig. 9. Profile reconstruction the case 10 for three levels of global attenuation length uncertainty: a) and b) τ λ = 10%, c) and d) τ λ = 20%, and e) and f) τ λ = 30%.The green curves represents the GT.In the panels a), c) and e), the profile reconstruction are plotted in blue (ρ|A m , y), orange (ρ|A m , y 0 ) and red (ρ|A m 0 , y) with their respective variabilities as shaded areas.In the panels b), d), and f), the a posteriori (P(ρ|A m , y)) in represented in blue and the marginals in orange (P(ρ|A m , y 0 )) and red (P(ρ|A m 0 , y)).

GTGT
Fig. 10.Profile reconstruction in the case W 10 for three levels of global length uncertainty: a) and b) τ λ = 10%, c) and d) τ λ = 20%, and e) and f) τ λ = 30%.The green curves represents the GT.In the panels a),c) and e), the profile reconstruction are plotted in blue (ρ|A m , y), orange (ρ|A m , y 0 ) and red (ρ|A m 0 , y) with their respective variabilities as shaded areas.In the panels b), d), and f), the a posteriori (P(ρ|A m , y)) in represented in blue and the marginals in orange (P(ρ|A m , y 0 )) and red (P(ρ|A m 0 , y)).

GTGT
Fig. 12. Reconstruction of concentration profile for two levels of acquisition noise: panels a and b) very low (σ k = 0.01, SNR = E[I m k ] 2 2 k ∈ × 10 3 , 16 × 10 6 ]), and panels c) and d) very high (σ k = 0.5, SNR ∈ [100, 6400]).In the panels a) and c), the profile reconstruction are plotted in blue (ρ|A m , y), orange (ρ|A m , y 0 ) and red (ρ|A m 0 , y) with their respective variabilities as shaded areas.In the panels b) and d), the a posteriori (P(ρ|A m , y)) in represented in blue and the marginals in orange (P(ρ|A m , y 0 )) and red (P(ρ|A m 0 , y)).
This is a model-free method whose output is the probability density σk C1s and the area I k , i.e. the total PE count from the signal of interest.The smooth output σk