Quantifying Intuition: Bayesian Approach to Figures of Merit in EXAFS Analysis of Magic Size Clusters

a Analysis of the extended X-ray absorption fine structure (EXAFS) can yield local structural information in magic size clusters even when other structural methods (such as X-ray diffraction) fail, but typically requires an initial guess – an atomistic model. Model comparison is thus one of the most crucial steps in establishing atomic structure of nanoscale systems and relies critically on the corresponding figures of merit (delivered by the data analysis) to make a decision on the most suitable model of atomic arrangements. However, none of the currently used statistical figures of merit take into account such significant factor as parameter correlations. Here we show that ignoring such correlations may result in a selection of an incorrect structural model. We then report on a new metric based on Bayes theorem that addresses this problem. We show that our new metric is superior to the currently used in EXAFS analysis as it reliably yields correct structural models even in cases when other statistical criteria may fail. We then demonstrate the utility of the new figure of merit in comparison of structural models for CdS magic-size clusters using EXAFS data.


Introduction
Establishing the atomic structures of materials is a fundamental step in understanding their mechanical, electronic and optical properties and is essential for material applications.However, recovering the atomic structures of nanomaterials is particularly challenging using standard structural analysis techniques (e.g.X-ray and electron diffraction, Raman scattering, etc.) due to loss of periodicity at atomic level and potential presence of novel metastable atomic arrangements.This is especially true of the recently discovered ultra-small truly mono-disperse nanoparticles -magic-size clusters (MSCs) 1 , 2 , 3 .As a consequence, several advanced structural methods such as X-ray absorption spectroscopy (XAS) and pair distribution function (PDF) analysis have been utilised to investigate atomic structure of MSCs 4 , 5 , 6 .XAS, in particular, has been shown to be sensitive to the atomic arrangements and structural changes in MSCs delivering information about sample stoichiometry and cluster symmetry discriminating between variety of structural models 6 .
This new class of nanoscale systems pushes XAS capabilities to the limit both in terms of the quality of the data required and of analysis methods for the two key parts of the X-ray absorption spectrum: X-ray absorption near edge structure (XANES) and extended X-ray absorption fine structure (EXAFS).The former is a QMUL, Mile End Road, London E1 4NS UK; Tel: ; E-mail: apw813@qmul.ac.uk b Diamond Light Source, Diamond House Harwell Science & Innovation Campus, Didcot OX11 0DE, UK.
† Electronic Supplementary Information (ESI) available: See DOI: 00.0000/00000000.sensitive to the symmetry around the absorbing atom of interest (e.g.Cd in CdS MSCs 6 ) and its oxidation state, while the latter provides information about local coordination numbers, interatomic distances and local atomic dynamics (see Fig. 1).
Analysis of EXAFS data typically involves background subtraction and normalisation followed by theoretical EXAFS calculations for a selected structural model (or a selection of model structures), comparisons of the calculated spectrum with the data and parameter refinement to obtain the best fit and the corresponding structural information 7 8 .Theoretical calculations and subsequent refinement are some of the most crucial steps and require a suitable atomistic model, thus implying some prior knowledge of the atomic structure or having an informed guess (e.g. based on molecular dynamics, DFT calculations or similar material, etc.).Recovering atomic structure of MSCs puts particularly stringent demands on model comparison in EXAFS because local atomic arrangements can be quite similar 9 .When theoretical and experimental spectra are compared, in most EXAFS analysis programs (such as Artemis 10 and Larch 11 ) there are a number of figures of merit (FoMs) available to provide quantitative model evaluation to answer the question of whether the model is a suitable match for the experimental data.However, none of the commonly used FoMs take into account parameter correlation.At the same time, it is well-known 12 , 13 , 14 that correlation can have significant negative consequence on data refinement (i.e.larger errors) and, most importantly, on model verification and selection.
This shortcoming of the EXAFS FoMs has long been recognised as a problem and in the latest development of the Artemis (one of the most commonly used EAXFS analysis package) a heuristic "happiness parameter" is offered to provide in-code indication of the fit quality.This parameter is based on decades of EAXFS analysis experience and includes, with varying weighting, an R-factor (a numerical measure of how well the fit over-plots the data), penalties for parameter correlations, restraints, the number of independent parameters, etc. 15 .While recognised as an important guide during data analysis, being a heuristic parameter, it has no firm basis in statistics and therefore cannot be quoted in publications.
In this article we introduce for the first time in EXAFS analysis an FoM that explicitly includes parameter correlations -the Bayes Factor Integral (BFI) 16 We use EXAFS data for crystalline Ge at low temperature to demonstrate that the BFI is more sensitive than the typical FoMs used for EXAFS analysis to model choice.We then demonstrate that the BFI consistently points towards the correct structure as preferred model.We then use the BFI to compare a selection of models for a material with unknown structure: CdS Magic Sized Clusters (MSCs 311 and 322) 6 5 17 .
With these examples, we introduce the BFI as a numerical metric for quantifying intuition in EXAFS model comparison.

Figures of Merit in EXAFS analysis
Least-squares fitting (LS, the minimisation of the sums of squares of residuals to optimise a model) is a commonly used method to fit data, to estimate parameters and to make decisions about model selection.In EXAFS LS fitting, reported FoMs in Larch and Artemis are χ 2 , χ 2 ν , R-factor, AIC and BIC.The first is the well-known statistical value characterising the residuals between the model and experimental data.It is a simple statistical measure of how small the fit residuals are: i.e. how closely does the model fit the data, However, it is has been well-established that the number of independent variables (fitting parameters) can significantly influence the value of χ 2 .
The total number of parameters available in EXAFS analysis is limited by the sampling theorem of the Fourier analysis 18 (this is also known as the Nyquist criterion/theorem in EXAFS community).Therefore, the most commonly reported fitting statistic in EXAFS is the so-called reduced chi-squared, χ 2 ν , based on χ 2 but with a modification to include normalisation of χ 2 by degrees of freedom such that once the maximum number of free parameters allowed for the data (N ind ) is reached 18 , it will become negative and provides a clear indication of over-fitting.The R-factor is another variation of the χ 2 criterion with a different normalisation factor.
The AIC and BIC are not found in Artemis, but are used as FoMs in Larch 11 to aid model comparison: both are based on the Likelihood function (rather than χ 2 ) while also including a penalty term for adding parameters to the model (adding fitting parameters to a model -physically meaningful or notnormally increases the likelihood of the model while reducing the probability that the model is correct).
There are a number of problems with the figures of merit described above.They treat all parameters alike, whether physically-meaningful or not.Apart from the number of parameters approaching N ind , there is not much help from these FoMs to tell whether one has a physically meaningful fit.Crucially, none of them include parameter correlations, while it is well documented 12 , 13 , 14 that parameter correlations indicate over-fitting and have significant consequences on the refinement errors and model selection.For example, in EXAFS analysis it is well-established that there exist correlations between fitting parameters even when N ind 18 is not exceeded 19 , 7 , 11 , 20 , 21 .
Both Artemis and Larch do provide functionality to calculate parameter correlations, but these are almost never used to assess the quality of the fit nor to aid model justification or selection.
To compensate for that and to help guide users during the refinement in Artemis there is an inbuilt FoM that does include correlations: the Happiness parameter 15 .However, it cannot be reported in publications since it has no mathematical basis: it is an empirical FoM that can be adjusted between fits to accord with the user's preferences.Hence, there is a need for a FoM rooted in statistics that does include parameter correlations.

Bayes Theorem, Bayes Factor and Bayes Factor Integral
The goal of EXAFS analysis can be described as "to find the best model parameters that fit the data" or, more generally, "to select the best model that fits the data".This lends itself naturally to the Bayesian statistical analysis and the use of Bayes theorem: where P(M|D) represents the conditional probability of the model M, given the data D. In the case of multi-parameter fitting (in- cluding EXAFS) this can be rewritten as (see, for example 22 ): where w is the vector of parameter values.Models can then be compared by, for example, taking a ratio of their conditional probabilities P(w|D, M).This ratio of probabilities of two models (e.g.i and j) represents the odds ratio in favour of one model over the other 23 : where I is prior information we have about the models.Using Eq. 2 it is straightforward to show that: where BF i j is the Bayes factor 23 .The ratio on the right hand side is the prior odds ratio of the two models and throughout this work we consider this to be unity (i.e.no preference of one model over another).Thus, to compare the models we need to compute P(D|M i , I) -probability of the data given the model and prior information.However, expressions for P(D|M i , I) can be rather complicated and for analysis involving many (in general correlated) parameters they include the evaluation of a multi-dimensional integral over the parameter space (the MLI -marginal likelihood integral).Assuming uniform top-hat priors and a Gaussian error distribution for independent identically distributed experimental data points gives (see for example 23 , p. 276): where ∆p i are prior parameter ranges, n is the number of the data points, m is the number of the model parameters p = [p 1 p 2 ...p m ] T and Cov p is the parameter covariance matrix.Thus, although the Bayesian approach has already been demonstrated in application to EXAFS analysis 24 25 , it has not been used to any significant extent, as far as we can tell, on account of its complexity.Indeed, parameter correlation is almost always the case in EXAFS and would normally require evaluation of the multidimensional integral in Eq. 5.That can be addressed by constructing an orthonormal set of the model basis functions 23 (model parameters) so that the new parameters will have no correlation and hence the multidimensional integral can be replaced by the product of multiple single integrals.However, this would require redefining the problem in terms of the new (orthogonal) parameter set, repeating the fit and then back-transforming the new parameters to recover the original ones.
Here we propose a simple alternative FoM that requires only trivial modifications to the statistical procedures already existing in EXAFS analysis.We note that the multidimensional integral on the far right of the Eq. 5 constitutes the volume in the parameter distribution space.We also note that Cov p is symmetric positive definite, hence its diagonalization involves basis rotation.However, the volume (and also det(Cov p )) does not change under rotation of the parameter space required for the transformation to the orthonormal basis.Hence, the following expression normally corresponding to the orthonormal parameter set can be used for calculation of P(D|M i , I) -Bayes Factor Integral (BFI) -for a model with parameter correlations present 16 : where ∆p i are the initial parameter ranges and L max is the likelihood for the model.Thus defined BFI can then be used for model comparison (giving preference for a model with the larger value of BFI) following EXAFS data fitting without the need to redefine the problem in the new orthonormal parameter set.We call this BFI (rather than, for example, MLI) to distinguish from a more common case when the orthogonal parameter set is used to obtain Eq. 6 (and therefore the Cov p is a diagonal matrix).Crucially, the FoM in Eq. 2 naturally incorporates the Occam factor: that accounts for parameter correlations as well as parameter ranges and provides a penalty for a model with significant parameter correlations and/or large initial parameter ranges (parameter uncertainty).However, since values of BFI can vary drastically, a more convenient way of evaluation is through comparison of ln(BFI) of the corresponding models.In such a case model evaluation is reduced to calculation of the ln(BFI) ratios -designated in this paper as ln(BF) -with the following scale 16 26 for ln(BF) values that differentiate between the models: • < 1 -barely worth considering, • 1 − 2 -substantial, • 2 − 5 -strong evidence, • > 5 -decisive.
We proceed below by testing this approach on a reference data set (crystalline Ge, c-Ge) before applying the procedure to the Magic Size clusters of CdS.
3 Results and Discussion

The Case of Crystalline Germanium
To verify the utility of BFI we first used the XAS data 27 for Germanium collected at 12K (the x-ray absorption spectrum is shown in Figure 2).The data were selected on the account of their high quality.The structure of Ge at this temperature is well-established and has been verified by previous publications 28 29 30 27 .The data analysis was performed in Larch 11 for background removal, normalisation and the actual fitting of the models to the experimental data since we found Larch to provide the most comprehensive fitting statistics.Figure 1 shows the EXAFS equation used to fit data where N is number of nearest neighbours, R is absorber-scatterer atom distance, S 2 0 is an amplitude reduction factor, σ 2 is Debye-Waller factor, F(k) is photoelectron scattering amplitude, λ (k) is photoelectron mean free path, and φ (k) is the phase shift.The latter three parameters (F(k), λ (k) and φ (k)) are calculated using the FEFF 9 code 31 32 and therefore are not refined.Three structural models have been selected for comparison: (i) the actual structure of crystalline Ge at 12 K known to be a 4coordinated face-centered diamond cubic type 29 33 (Figure ??, Fd3m, Model 1); (ii) 6-coordinated high-pressure Ge phase VI structure (Cmma, Model 2); (iii) 6-coordinated high pressure β -Sn structure of Ge 34 (I4/amd, Models 3a, 3b).For Model 3 two different refinements were carried out: one (Model 3a) was for a single shell of 6 nearest neighbours, while for the other (Model 3b) 3 shells of 2 atoms were used to reflect the actual nearest neighbour configuration in the β -Sn structure.This was also used to gauge the effect on the BFI of increasing the number of model parameters.
To enable a fair model comparison, for each model we only looked at the first peak in the R-space (corresponding to the Ge-Ge bond length of 2.45 Å; the atomic shell structure beyond the first shell is very different in the three selected models) and we used single-scattering paths only (see Figure 2).The data were fitted over the range of 2.00 < k < 22.93 −1 in k-space (see Fig. 2).This ensured that only the first-shell EXAFS were fitted.Parameter ranges are given in Table 1 and are defined as follows.The amplitude reduction factor S 2 0 corrects for inelastic effects in the absorbing atom 35 .This is empirically established to be in the range 0.8-1, and is well-covered by 0.5 range.The shift in the edge position E 0 accounts for errors in experimental calibration and for empirical convention in determination of the absorption edge position 19 7 -the range typically does not exceed 10 eV.Relative change in the nearest-neighbor interatomic distance ∆R is not expected to exceed 10% as the interatomic distances are determined by the covalent radii of elements and the pressuretemperature conditions (as an example, 10% bond length variation is well above that expected on melting or under pressures as high as 10 GPa in a typical semiconductor material such as Ge 30 ).Mean squared relative displacements of atoms due to atomic vibrations σ 2 (and static disorder, if any) accounts for damping effects on χ(k).The initial value can be calculated using e.g.correlated Debye or Einstein approximations 36 37 38 and for c-Ge at 12 K this is around 0.003 Å 2 27 hence the range of ±0.003 is selected to make it positive-definite.The number of nearest neighbours N was set according to the structural models and was not refined.
For Model 1 (zinc-blende structure), one single-scattering single-shell path was used to fit the spectrum.For Model 2 (the high-pressure Cmma structure), the spectrum was fit with 3 single-scattering single-shell paths between (in total) 6 atoms in the first shell to describe the signal.For (β -Sn) Model 3a, one    The summary of the results for the ln(BF) (the difference between ln(BFI) values) are shown in Fig. 3.One can see that Model 1 is favoured over all other models except for 3a (the single path β -Sn fit): the ln(BF) between Model 1 and all other models (except 3a) are found to be > 3 providing strong evidence for Model 1 being the preferred structure.The ln BF value between Model 1 and 3a is 0.67 is slightly in favour of Model 1 but not statistically significant according to the criteria outlined at the end of the previous section.However, the currently available fitting statistics FoMs found in the corresponding tables † favour other models: Model 2 has lowest χ 2 and χ 2 ν , Model 3b has the lowest value of R-factor, while AIC and BIC favour Model 2 over Model 1.This shows that reliance only on the currently used FoMs in EXAFS analysis can lead to an incorrect atomic structure model as the best solution.At the same time, we see that the BFI is able to deliver the correct result in this relatively complicated caseafter all we used a single peak only in the EXAFS FT magnitude in order to differentiate between the models.Having verified the utility of the proposed BFI-based FoM in case of the reference system, in the next section we apply the Bayes approach to analysis of EXAFS of MSCs.

Bulk CdS k-space Fitting
Before proceeding on CdS MSCs we further tested the utility of the new FoM in k-space fitting of the bulk crystalline CdS.As reference data, bulk crystalline CdS EXAFS data at 90K at Cd K-edge were fit in Larch 11 .The first shell (Cd-S scattering paths) in k-space was fit using several different structures respectively: a zinc-blende, wurtzite 39 , NaCl-like 40 and cmcm sturcutres 41 , latter two being high-pressure derived structures.
Parameter Initial value Range ∆R 0 0.1 E0 0 10 S0 2 0.9 0.5 σ 2 0.006 0.006 K-space noise (ε k ) was was evaluated from the signal between 6.50<k<18.30Å and the fit to the EXAFS data was carried out in the region 2.50<k<15.0Å, parameter ranges for the BFI calculations are shown in Table 3.The results of the fit are shown in the Table 4.The BFI-based FoMs support zinc-blende and wurtzite structures significantly over the NaCl-like and cmcm models.This is consistent with our previous results where XPDF analysis of bulk CdS (and of regular CdS quantum dots) has shown CdS to be a mix of zinc-blende and wurtzite structures 5 .

Magic Sized Clusters
Magic Sized Clusters (MSCs) are ultra-small (<3nm) colloidal semiconductor systems 4 .They are materials of interest due to their monodisperse nature 3 that suggests one can deliver atomiclevel control of system size using colloidal synthesis route.Their atomic structure is still under debate as are the methods of their structural verification.One of the key challenges for the latter is the possibility for stable and multiple meta-stable) atomic arrangements that are size-and temperature-dependent 5,6,9 .
The MSCs under investigation in this work are CdS.These MSCs exhibit a sharp UV-vis absorption peak at 311 nm (MSC 311) but when heated to 60 • C 4 this peak shifts to 322 nm (MSC 322) 6 and the shift is accompanied by atomic structure rearrangement as indicated by x-ray pair distribution function (xPDF) and XAS analysis 5,6 .Due to their small size leading to the lack of long-range order, establishing the atomic-level structure of MSCs is challenging 42   found in the bulk can be a starting point if there is no other information (this is frequently the case).
ii) It has been observed 30 that average interatomic distances in small nanoparticles can be reduced compared to their bulk counterparts.This can be interpreted as an effective pressure on these systems.Such compression may results in distortion towards the β -Sn structure 43 , hence it is reasonable to use it as one of the structural models.
iii) It has recently been shown that an InP-like structure 44 provides the best fit to PDF 5 data in CdS MSCs.
All clusters have been cut as spherical regions of appropriate size from the corresponding bulk crystalline structure and were terminated with oxygen (except for InP-like cluster where the structure from our recent work 5,6 was used).The bond length was set to that found in zinc-blende and wurtzite structure at ambient condition (2.55 Å).After fitting with the initial structures they were all relaxed over 500 steps using Avogadro 45 Universal Force Field 46 and fitting process has been repeated.Data analysis has been carried out in Larch 11 and Mathematica 13.0 (for BFI calculations).The initial parameter values and their ranges for all models are given in Table 3 (if the standard errors returned by Larch exceed the half range set for ∆p i terms, the range is corrected to 2 × standard error).
For each model the five highest-importance single-scattering paths were fitted up to the second cumulant (σ 2 ), we then fitted the third cumulant for the highest-importance Cd-O single scattering path of those included in the fit already.While at low temperatures (90 K in our case) 3 rd cumulants for the atoms in the bulk of a nanoparticle are expected not to be significant, this may not be the case for the surface Cd-O coordination shell.The ln(BFI)s for and ln(BF)s between the models were calculated using the same method as in the Ge model comparisons with the fits carried out in the range of 2 − 17 Å −1 for the the coordination shell 1 − 3 Å (see Fig. 4).For the noise value used in EXAFS chi-squared calculations we used the standard deviation of the kspace spectra for the MSC 311 and 322 data at 14.50 − 15 Å −1 and 13 − 14 Å −1 respectively.
The results for ln(BFI) for as-prepared and relaxed model clusters for 311 and 322 MSCs are given in the Tables 5 and 6 below and summarised in Figs.5a and 5b      as the most probable structure; ii) is detecting the difference between the EXAFS data of the two MSCs (as reflected in changing values of FoM and model preferences).This is a very significant finding considering the differences between the EXAFS signals for MSC 311 and 322 are very small (see Figs. 4) and the corresponding χ 2 -values are also quite close (see Supporting Information).
The preference for InP-like structure is certainly in line with previous work 6 5 where it was shown to provide the best model to fit xPDF, EXAFS and XANES data.All this provided very strong support to using the new Bayesian FoM we introduced in this work as a universal metric for model comparison and selection.

Conclusions
In this work we introduced the Bayesian-based statistical metric, the Bayes Factor Integral, for model comparison in EXAFS analysis.We showed for the first time that the new FoM provides a superior tool for model comparison in EXAFS by quantifying the intuitions about the parameter ranges and correlations through the Occam factor.We tested the new Bayesian FoM against reference EXAFS data for Ge and demonstrated that the it is superior to the FoMs typically used in EXAFS analysis and reliably predicts the correct structure.In the process we showed that ignoring  model parameter correlations may result in a selection of an incorrect structural model.We then applied the new FoM for model comparison in analysis of MSCs where we demonstrated that it is sensitive to the differences between EXAFS signals of MSCs 311 and 322 and can point to the most probable structural models for these systems.So far we utilised identical parameter ranges for all models of interest in all our tests (except for the range correction when the standard error ≥ 1/2∆p i ).Imposing stricter model-specific parameter ranges in the calculation of model BFI values should provide better results in model selection.Such model-specific parameter ranges can be obtained, for example, from moleculardynamics or ab-initio simulations.
We note that the current approach is so far limited by the requirement of providing initial guess structures, while the ultimate goal of structural analysis of MSCs (and of non-periodic systems in general) is producing structural models for materials with unknown structures.We believe that newly developed FoM can be the crucial aid on that pathway, for example by combining analysis described in this paper with machine learning/AI methods.

Fig. 1
Fig. 1 EXAFS analysis is the study and interpretation of the fluctuations in the post-edge X-ray absorption spectrum.The fluctuations in the signal (in the purple highlight) are the result of interference of the outgoing photoelectron wave with the portion scattered by the neighboring atoms.The EXAFS equation used for modelling these oscillations.χ(k) is related to the plotted absorption µ(E) by the transform: k = 2me h2 (hω − E 0 ), E 0 being absorption edge energy.
(a) Ge spectrum in E-space.(b) Ge spectrum in k-space.(c) Ge spectrum in R-space.
and in this work we examine the sensitivity of our new EXAFS BFI-based FoM to the structural model selection.To this end we compared 4 models as candidates of possible structures of MSCs 311 and 322: i) Cd 40 S 19 with zinc-blende structure, ii) Cd 40 S 20 with wurtzite structure, iii) Cd 33 S 32 β -Sn-like structure and ii) Cd 37 S 20 with an InP-like structure.The rationale for the model selection is as follows: i) Bulk CdS can possess wurtzite or zinc-blende structure (while regular CdS quantum dots are known to exhibit both characters 5 ).Hence, when constructing a model for an unknown atomic structure of MSCs, the atomic in the case of MSCs 311 InP and Zinc-Blende models are almost equally favoured over the β -Sn and Wurtzite structures.The results for ln(BFI) of 322 MSCs show that InP-like structure is significantly favored over all other models.Thus, our findings here indicate that the new FoM: i) is pointing to the InP-like model (a) R-space EXAFS of the MSCs.(b) k 2 -space EXAFS of the MSCs.

Fig. 5
Fig. 5 Charts of the results for each MSC.

Table 2
Models and Ln(BF) values for the different Ge models.

Table 3
Parameter values and ranges for all MSCs.

Table 4 Ln
(BFI) values for each model in the bulk CdS EXAFS fitting (1st shell).
. These results show that