Bayesian sparse modeling of extended x-ray absorption fine structure to determine interstitial oxygen positions in yttrium oxyhydride epitaxial 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_xH_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_xH_y thin films belongs to a face-centered cubic (fcc) lattice¹⁵ and its lattice constant takes 5.34¹⁶ or 5.35 Å.¹⁵ It is also reported that the oxygen coordinates around the Gd atom in gadolinium oxyhydride thin films.¹⁷ However, how the oxygen coordinates around the Y atom in YO_xH_y thin films remains elusive.

Epitaxial thin films of YO_xH_y(111)/CaF₂(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_xH_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,¹⁹ 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 $γαβ(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−E0+ΔE)/ℏ$ as a function of photon energy E, where E₀, Δ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 $γαβ(2)(r,k)$ is written as an oscillating function, A_αβ(r, k)sin Φ_αβ(r, k).¹⁹ The amplitude A_αβ(r, k) and the phase shift Φ_αβ(r, k) are calculated by GNXAS,^20–22 in which a complex Hedin–Lundqvist potential²³ 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χα(ki)|i=1,…,p},$ and the EXAFS analysis becomes a linear regression problem with a given design matrix $Mij=kiwAαβ(rj,ki)sinΦαβ(rj,ki)Δ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₁ regularization is particularly effective in the EXAFS analysis of crystalline solid materials⁷ 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)²⁴ 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.⁹ 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,¹⁰ 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.²⁶ 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,28F(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).¹⁰ For P(x|c), we use a normal distribution $N(0,σxβ)$ 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 $σxβ=rj2(zxβ0)−1/2,$ where $zxβ0$ of each element β will be estimated by BFE minimization.

For P(y|x, c), we assume a multivariate normal distribution $N(M(c◦x),Σ)$ for the superimposed noise in y and also introduce a model for its variance–covariance matrix Σ. Since we analyze ${yi=kiwχα(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)², and the off-diagonal elements are zero. 1/z and σ(k)² 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 $σ(k)=(k/knoise)w(σdata2−1/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¹⁰ 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 Å⁻¹, and the region to obtain σ_data is k = 8.0–10.5 Å⁻¹. 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.²⁹ 

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≔Power(cSpM)$ of c_SpM through the minimization of the BFE, $ĉBI(ΔE,λ)=argminc⊆CF(c;ΔE,λ).$ The exhaustive search is performed under ‖c‖₀ < 20 for $c⊆C$ because of the computational cost. (iii) Finally, we obtain the minimized BFE $[F̂(ΔE)=minλF(ĉBI(ΔE,λ))]$ 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̂(Δ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³χ(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̂BI$ 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.³⁰ 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₀, 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 ith 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₄/d₂ = 5.325/3.750≈1.420, which is very close to $2≈1.414$⁠, characterizes the fcc structure. This result confirms a previous work¹⁵ that the Y atom takes the fcc structure. Additionally, a lattice constant of fcc is also obtained from $2d2$ and d₄, 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₂/d₁, 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₂/d₁ is $8/3≈1.633$⁠. On the other hand, d₂/d₁ for octahedral becomes $2≈1.414$⁠. In our analysis, d₂/d₁≈1.648(=3.750/2.275) is quite close to the tetrahedral one, and we can conclude that the interstitial O sites in YO_xH_y are tetrahedral. A similar result has been obtained in a previous study¹⁷ 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₁ 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_xH_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).