This article presents a Bayesian sparse modeling method to analyze extended x-ray absorption fine structure (EXAFS) data with basis functions built on two-body signals. This method does not require any structural model and allows us to evaluate regression coefficients proportional to the radial distribution functions of the respective elements and their errors and is very effective for analysis of EXAFS with weak absorption intensity and severe signal-to-noise ratios. As an application example, we used it to analyze the EXAFS of an yttrium oxyhydride (YO_{x}H_{y}) epitaxial thin film. These EXAFS data show weak absorption intensity and a severe signal-to-noise ratio due to the small amount of x-ray absorption in the thin film sample. However, this approach revealed that the radial distance ratio of the second neighbor yttrium to the first neighbor oxygen coincides with that of a tetrahedral configuration. This result demonstrates that the interstitial oxygen position is tetrahedral in the YO_{x}H_{y} thin film.

The extended x-ray absorption fine structure (EXAFS)^{1,2} is a primary method to study the atomic-scale geometry around an element specified by the energy of the x-ray absorption edge. While the peak-fitting approach and the regularization method^{3,4} have often been used to analyze the EXAFS, deep-learning applications^{5,6} to the EXAFS would afford the key to improve the study of atomic-scale geometry. However, there is another alternative machine learning approach for measured EXAFS data analyses to avoid overfitting that attempts to reproduce even noise. Although efforts for reducing the measurement noise are important, the overfitting becomes a serious problem, especially when the absorption intensity is weak, such as in thin films.

In this article, a unified method of sparse modeling (SpM)^{7–9} and Bayesian inference (BI)^{10–12} is proposed to analyze the EXAFS data.

Our method has the following advantages: (1) only the elemental species information is required for the analysis, and no structural information is needed; (2) the phase shift of the electron wave due to back-scattering is incorporated so that the radial distance can be accurately estimated; (3) a model for *k*-dependent noise (due to *k*^{w} weight) is introduced so that the noise intensity, which varies with *k*, can be properly evaluated; and (4) as a result, even EXAFS signals with a severe signal-to-noise ratio can be analyzed with high noise immunity.

We apply this method onto the EXAFS data measured on an yttrium oxyhydride YO_{x}H_{y} epitaxial thin film in order to determine the atomic-scale structure of the oxygen (O) atom around the yttrium (Y) atom. This material is promising for switchable optical properties such as smart windows and optical sensors^{13,14} because it exhibits a reversible light-induced insulator-to-metal transition. In addition, by fabricating thin films, the photoresponsive character is expected to be enhanced than that of polycrystalline thin films. Recent studies have shown that the Y atom in YO_{x}H_{y} thin films belongs to a face-centered cubic (fcc) lattice^{15} and its lattice constant takes 5.34^{16} or 5.35 Å.^{15} It is also reported that the oxygen coordinates around the Gd atom in gadolinium oxyhydride thin films.^{17} However, how the oxygen coordinates around the Y atom in YO_{x}H_{y} thin films remains elusive.

Epitaxial thin films of YO_{x}H_{y}(111)/CaF_{2}(111) substrates were grown using reactive magnetron sputtering. Details for sample preparation are provided in Ref. 18. Yttrium K-edge EXAFS measurements were performed at room temperature in air at BL14B1 in SPring-8. A Si(311) double-crystal monochromator and a 36-element solid-state detector were used for EXAFS observation with the fluorescence mode. Rh-coated double mirrors were applied for cutting higher-order reflection in incident x rays.

Figure 1 shows the EXAFS data *k*^{w}*χ*(*k*) (*w* = 1, 2, and3) of an YO_{x}H_{y} epitaxial thin film. Abscissa is the wavenumber *k* of photoelectrons. Although the EXAFS data are analyzed with the absorption intensity *χ*(*k*) multiplied by *k*^{w} to emphasize the oscillation at larger *k*, the weight *w* has been empirically selected. As seen in Fig. 1, the signal fluctuations are large in the noise region at larger *k* because of the weak absorption intensity due to a thin film sample and the measurement in the fluorescence mode: in such cases, the EXAFS analysis is difficult with conventional Fourier transform methods. As discussed later, the proposed method makes it possible to determine which order *w* should be used to interpret the data as well as to study the atomic-scale geometry even in such a tough analysis with high noise intensity.

Although the EXAFS oscillation of *χ*(*k*) comes from interferences of photoelectron waves, we consider two-body terms,^{19} which provide a dominant EXAFS oscillation,

where *α* and *β* mean elements of the photo-absorber and scattering atoms, respectively. *n*_{αβ}(*r*) is the coordination number distribution of *β* as a function of the distance *r* from *α*, and $\gamma \alpha \beta (2)(r,k)$ contains the two-body multiple scattering term corresponding to the same geometrical configuration. Here, the wavenumber *k* is given by $k=2m(E\u2212E0+\Delta E)/\u210f$ as a function of photon energy *E*, where *E*_{0}, Δ*E*, and *m* are the absorption edge energy, the energy offset of the theoretical zero for the *k*-space, and electron mass, respectively. In our method, this Δ*E* can be optimized as described below. Each signal $\gamma \alpha \beta (2)(r,k)$ is written as an oscillating function, *A*_{αβ}(*r*, *k*)sin *Φ*_{αβ}(*r*, *k*).^{19} The amplitude *A*_{αβ}(*r*, *k*) and the phase shift *Φ*_{αβ}(*r*, *k*) are calculated by GNXAS,^{20–22} in which a complex Hedin–Lundqvist potential^{23} is used for inelastic loss effects in the framework of the multiple scattering theory.

We convert the integral of *r* in Eq. (1) into a Riemann sum with *N* subintervals whose each length is Δ*r*,

Let **y** be ${yi=kiw\chi \alpha (ki)|i=1,\u2026,p},$ and the EXAFS analysis becomes a linear regression problem with a given design matrix $Mij=kiwA\alpha \beta (rj,ki)sin\Phi \alpha \beta (rj,ki)\Delta r.$ We introduce regression coefficients **x**_{β} = {*x*_{β}(*r*_{j})} for an element *β* and **x** = {**x**_{β}} for all elements in proportion to the coordination numbers.

This linear regression is an ill-posed problem. However, the *L*_{1} regularization is particularly effective in the EXAFS analysis of crystalline solid materials^{7} because neighbor atoms will coordinate at discrete (sparse) distances from the photo-absorbers on the basis of symmetry and chemical structure. Therefore, we employ the least absolute shrinkage and selection operator (LASSO)^{24} method, and a sparse solution **x**_{SpM} is obtained by the formula at some regularization parameter *λ* (>0),

Iesari *et al.* used this LASSO method for the EXAFS signals of Cu foil, Fe foil, and cuprous oxide and determined *λ* by cross-validation.^{9} As a result, they succeeded in obtaining the peak positions of the radial distribution for their materials and structural parameters, such as average distances and Debye–Waller factors, without using any prior information about the structure.

However, the LASSO method would select jointly correlated *x*_{β}(*r*_{j}) when their features in *M* are highly correlated^{24,25} and *λ* is often chosen heuristically. To overcome these problems by BI,^{10} we introduced a binary vector **c** indicating the nonzero components of **x** and applied linear regression **y** = *M*(**c**◦**x**) + ** ϵ**, where

**is the noise superimposed on**

*ϵ***y**and the operator ◦ means the Hadamard product between vectors. In BI, we maximize the posterior probability

*P*(

**c**|

**y**), where it can be expanded as

*P*(

**c**|

**y**) =

*P*(

**y**|

**c**)

*P*(

**c**)/

*P*(

**y**) on the basis of Bayes’ theorem.

^{26}To perform an unbiased analysis, the prior probability

*P*(

**c**) should be uniform. Consequently,

*P*(

**c**|

**y**) is proportional to the likelihood

*P*(

**y**|

**c**), and it is given by the following equation:

We use the Bayes free energy (BFE),^{10,27,28}*F*(**c**) = −ln *P*(**y**|**c**), for model selection of **c**.

To calculate Eq. (4), we introduce appropriate models for the prior probability *P*(**x**|**c**) and the likelihood *P*(**y**|**x**, **c**).^{10} For *P*(**x**|**c**), we use a normal distribution $N(0,\sigma x\beta )$ for the nonzero components in **x**_{β}, where we incorporate the sparsity of **x**_{β} by being the mean to zero. On the other hand, the coordination numbers of nonzero components are expected to increase with $rj2$ because the volume of a spherical shell increases with $rj2$. So we set to $\sigma x\beta =rj2(zx\beta 0)\u22121/2,$ where $zx\beta 0$ of each element *β* will be estimated by BFE minimization.

For *P*(**y**|**x**, **c**), we assume a multivariate normal distribution $N(M(c\u25e6x),\Sigma )$ for the superimposed noise in **y** and also introduce a model for its variance–covariance matrix Σ. Since we analyze ${yi=kiw\chi \alpha (ki)}$, *y*_{i} is considered to include uniform *k*-independent noise and *k*-dependent noise in proportion to *k*^{w}, where the variance of the uniform noise and the standard deviation of the *k*-dependent noise are 1/*z* and *σ*(*k*), respectively. The diagonal elements of Σ are the sum of 1/*z* and *σ*(*k*)^{2}, and the off-diagonal elements are zero. 1/*z* and *σ*(*k*)^{2} are formulated as follows: First, the fluctuation in the noise region (*k* ≥ *k*_{noise}) of Fig. 1 is considered to come from noises because the EXAFS oscillation almost decays, and the root mean square deviation of the fluctuation is obtained as *σ*_{data}. Subsequently, *σ*(*k*) is given by $\sigma (k)=(k/knoise)w(\sigma data2\u22121/z)1/2$ for *k* < *k*_{noise}. We also estimate 1/*z* by minimizing BFE.

By the modeling of *P*(**x**|**c**) and *P*(**y**|**x**, **c**), we can obtain the posterior probability and BFE analytically. The regression coefficient **x** is optimized by the maximum *a posteriori* (MAP) estimation^{10} under a non-negative constraint. As the posterior probability *P*(**x**|**y**, **c**) is the Gaussian-type function, the reliability of **x** can also be estimated from the diagonal elements of the variance–covariance matrix of *P*(**x**|**y**, **c**).

The EXAFS oscillation shown in the analysis region of Fig. 1 is chosen for analysis because some x-ray absorption near edge structures would remain in the small *k* region. The analysis region is *k* = 3.0–8.0 Å^{−1}, and the region to obtain *σ*_{data} is *k* = 8.0–10.5 Å^{−1}. For radial distances *r*_{j}, we use a linear gird of *r*_{j} = 2.000–9.975 Å with the step Δ*r* = 0.025 Å.

In BI, the number of combinations for the binary vector **c** that represents the nonzero components of **x** is enormous and their exhaustive search is impractical. So we first use Eq. (3) to get a sparse solution **x**_{SpM}, *i.e.*, **c**_{SpM}. The sparsity control parameter *λ* is set to be in the range of *λ*_{max}–*λ*_{max}/100 in a geometric sequence, where all **x**_{SpM} are zero for *λ* ≥ *λ*_{max}.^{29}

Since the design matrix *M* contains the parameter Δ*E* to be optimized, we obtain the BFE as a function of Δ*E* for each *w* as follows: (i) We obtain **x**_{SpM} (**c**_{SpM}) by solving Eq. (3) at given Δ*E* and *λ* and calculate the MAP solution **x**_{BI} for the obtained **c**_{SpM}. (ii) We perform exhaustive search in the power set $C\u2254Power(cSpM)$ of **c**_{SpM} through the minimization of the BFE, $c\u0302BI(\Delta E,\lambda )=argminc\u2286CF(c;\Delta E,\lambda ).$ The exhaustive search is performed under ‖**c**‖_{0} < 20 for $c\u2286C$ because of the computational cost. (iii) Finally, we obtain the minimized BFE $[F\u0302(\Delta E)=min\lambda F(c\u0302BI(\Delta E,\lambda ))]$ for each Δ*E*. In these procedures, the exhaustive search is necessary because the basis functions are not strictly orthogonal to each other and there is a concern about excess extraction of the correlated components.

Figure 2 shows $F\u0302(\Delta E)$ for *w* = 1, 2, and3 in the range of Δ*E* = 0.0–10.0 eV with 0.1 eV step. The information criterion minimizing the BFE selects Δ*E* = 1.0 eV and *w* = 1. Although the same Δ*E* is selected for *w* = 2, Δ*E* for *w* = 3 differs from others. This is because the noise intensity is overestimated for a large exponent for high *k*-weighting. In Fig. 1, the amplitude of *k*^{3}*χ*(*k*) in the analysis region is smaller than that in the noise region.

The red curve in Fig. 3(a) is the regression curve for the selected model $x\u0302BI$ with Δ*E* = 1.0 eV and *w* = 1, and it is found that the red curve explains appropriately the measured data indicated by open circles without overfitting. In the EXAFS analysis, we have to pay attention to the contributions from three- or more-many body multiple scattering effects.^{30} However, as seen from the scattered open circles in Fig. 3(a), it is considered that such contributions are not important in the case of severe signal-to-noise ratios because of a thin film sample. On the other hand, the blue dashed curve is drawn only by components of the two nearest neighbor atoms (O: 2.275 Å, Y: 3.750 Å) and presents a good overview of the EXAFS oscillation. In our method, the elements can be treated separately as *β* in Eq. (1), and the radial distribution of each element can be obtained separately as shown in Figs. 3(b) and 3(c), which are for Y and O atoms, respectively. Error bars show the reliability estimated from the variance–covariance matrix of *P*(**x**|**y**, **c**). Since we do not consider the effect of the reduction factor *S*_{0}, the regression coefficients will be underestimated. Nevertheless, the coefficients take large values. This is because that the thermal effect is excluded by the exhaustive search. As the EXAFS measurements were performed at room temperature, each peak in Figs. 3(b) and 3(c) should be distributed with a finite width near the equilibrium atomic position. However, it is very difficult to estimate the width of the radial distribution function from the data with a severe signal-to-noise ratio. Because the exhaustive search excludes the correlated coefficients and the contributions of the thermal effect due to a trade-off between excluding correlated coefficients and considering the thermal effect, the coefficients stand only at the equilibrium position despite the finite temperature and become large.

In the used basis functions, the phase shift due to photoelectron wave scatterings is incorporated. Therefore, the distance *d*_{i} for the *i*th neighbor atom can be evaluated correctly, and we can estimate the atomic-scale three dimensional structure based on the distance ratio between the neighbor atoms. When we focus on Y atoms in Fig. 3(b), *d*_{4}/*d*_{2} = 5.325/3.750≈1.420, which is very close to $2\u22481.414$, characterizes the fcc structure. This result confirms a previous work^{15} that the Y atom takes the fcc structure. Additionally, a lattice constant of fcc is also obtained from $2d2$ and *d*_{4}, and their weighted-mean value with the weights of the inverse squares of the error bars in Fig. 3(b) is 5.306 ± 0.025 Å. It is slightly smaller than the previous studies^{15,16} (5.34–5.35 Å), and this difference might come from the epitaxial growth.

The coordination structure of O atoms around the Y atom can be determined from the ratio *d*_{2}/*d*_{1}, where there are two candidates for the interstitial O-sites: tetrahedral and octahedral sites. If the interstitial sites form a regular tetrahedron, the ratio *d*_{2}/*d*_{1} is $8/3\u22481.633$. On the other hand, *d*_{2}/*d*_{1} for octahedral becomes $2\u22481.414$. In our analysis, *d*_{2}/*d*_{1}≈1.648(=3.750/2.275) is quite close to the tetrahedral one, and we can conclude that the interstitial O sites in YO_{x}H_{y} are tetrahedral. A similar result has been obtained in a previous study^{17} for the EXAFS of gadolinium oxyhydride thin films, in which the conventional Fourier transform was used with a structure model predicted by the x-ray diffraction. In contrast, the advantage of our method is obvious since it is based on the accurate atomic distances compensated with the phase shift and structure models are not required.

In summary, we proposed a Bayesian SpM for analysis of EXAFS data with GNXAS basis functions incorporating two-body multiple scatterings. The model selection for the sparse solution of the regression coefficients proportional to the radial distribution function was performed by BFE minimization through the exhaustive search within the nonzero components obtained by *L*_{1} regularization. Using this method, we succeeded in determining that the sublattice of the Y atom is fcc and the interstitial O position is tetrahedral in the YO_{x}H_{y} epitaxial thin film.

The synchrotron radiation experiments were performed using a JAEA experimental station at JAEA (QST) beamline BL14B1, SPring-8, with the approval of the Japan Synchrotron Radiation Research Institute (JASRI) (Proposal No. 2020A3648) and the JAEA Advanced Characterization Nanotechnology Platform under the remit of “Nanotechnology Platform” of the Ministry of Education, Culture, Sports, Science, and Technology (MEXT), Japan (Grant No. JPMXP09A20AE0006). This work was supported by JST CREST (Grant Nos. JPMJCR1861, JPMJCR1761, and JPMJCR21O1), Japan; JST PRESTO (Grant Nos. JPMJPR17N2 and JPMJPR17N6), Japan; JSPS KAKENHI (Grant Nos. JP19H02596, JP19H04689, JP18H05513, and JP18H05514), Japan; and AGC, Inc. (Grant No. KC31AGC08p).

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.