Accuracy of XAS theory for unraveling structural changes of adsorbates: CO on Ni(100)

Surface-sensitive x-ray techniques represent an invaluable tool in surface science. The tunability and intensity of modern synchrotron sources result in high-resolution x-ray spectra, allowing the identification of adsorbed atomic and molecular species, as well as their adsorption sites. The chemical environment of a core-ionized atom can be determined using x-ray photoelectron spectroscopy (XPS), while additional information about the local electronic structure can be obtained from x-ray absorption (XAS) and emission spectra (XES). Detailed knowledge about the formation and population of bonding and antibonding states has been gained for a number of atomic and molecular adsorbates (e.g., C, O, N, CO, and N₂ on transition metal surfaces) relevant for common heterogeneous catalytic processes.^1,2

The advent of free-electron lasers delivering ultrashort x-ray laser pulses offers new opportunities for gaining knowledge about chemical reactions on the sub-picosecond time scale where bonds typically form and break. Reactions on the surface, and/or desorption from it, can be initiated using an ultrashort optical laser pulse, while a probing x-ray pulse with a duration of ∼10 fs gives time-resolved spectra during the evolution of the system, with sub-ps resolution.^3–7 This has led to a new understanding of surface reaction pathways.^4,8–10 However, the interpretation of the x-ray probe spectrum relies on accurate theoretical modeling of spectra from whatever species and their structure that may be present and dynamically evolving on the surface. When moving toward systems relevant for industrial catalysis, the number of intermediates and non-equilibrium structures in the given reaction can be high and reliable simulations of the spectra are crucial.

Due to the complicated electronic response to core-hole creation, calculating XAS/XES of condensed phase species formally requires a many-body formalism, and this is often accomplished by solving the Bethe–Salpeter equation.^11,12 However, because many-body effects are generally negligible for the lower levels of adsorbates on metal surfaces, much simpler density functional theory (DFT) calculations are often used to simulate XAS and XES. While this does not allow the treatment of metastable structures where electronic excitations are important, the short lifetime^13,14 (∼fs) of such excitations of adsorbates on metal surfaces, compared to the time resolution (∼100 fs) in pump-probe experiments, means they have a negligible effect on the spectrum compared to adsorbate vibrational excitation, diffusion, and desorption. DFT is known to capture trends in bonding of adsorbate systems remarkably well, and calculated trends of XAS/XES of adsorbates on a surface have shown a good overall agreement with experiment.^15,16 However, the interpretation of the ultrafast dynamics occurring on surfaces that is observed by x-ray spectroscopy requires a quantitative theoretical understanding of the relationship between the structural changes occurring on the surface and the x-ray spectra. In order to assess how quantitative this can be, we investigate the accuracy of a standard DFT analysis of XPS and XAS, employing a Generalized-Gradient-Approximation (GGA) functional, a plane-wave basis set, and core-hole pseudopotentials. This approach is, therefore, scalable to fairly large, periodic systems including itinerant magnetism within the framework commonly used to study periodic solids. We show that this gives reliable and accurate results for the XAS structure dependence on a system where static experimental data are available for the same adsorbate/metal system at a number of different binding sites that have all been well characterized by a number of different surface science techniques.

We choose the CO/Ni(100) system since the CO molecule binds with very similar adsorption energies in top, hollow, and bridge sites. By changing the CO coverage, the molecule moves to different adsorption sites, and the co-adsorption of hydrogen can be used to give a few different stable, high-coverage phases.¹⁷ Detailed experimental XAS spectra are, therefore, available for these well-characterized systems for comparison with our calculations.^18,19 Without hydrogen, CO at a coverage of 0.5 ML binds on-top and forms a c(2 × 2) structure, while at θ = 0.67 ML coverage, a p(3$2$ × $2$⁠)R45 phase is formed with all CO molecules at bridge sites; the coverage depends on the adsorption temperature and background CO pressure.²⁰ Initially, adsorbing 1 ML of hydrogen at 80 K and then dosing with CO lead again to on-top CO in a c(2 × 2) structure with 0.5 ML saturation coverage.²¹ Annealing this structure to 170 K leads to a c(2$2$ × $2$⁠)R45 phase, with the half CO occupying top sites and the other half hollow sites. Atomic hydrogen is found in all remaining hollow sites.¹⁷ 

Our structures are shown in Fig. 1. Since XAS has been recorded for all these surface phases using an identical experimental setup,^18,19 it is an ideal system for comparing to theoretical spectra. We, thus, calculate O and C K-edge XAS spectra for all the mentioned surface structures using the DFT methods described below. For comparison, the adsorption onset was also calculated for a pure ¼ ML CO coverage at different binding sites.

We use the DFT code Quantum ESPRESSO,^22,23 which uses a plane-wave basis for the valence electrons and pseudopotentials to represent the core electrons. Plane waves and pseudopotentials are the standard approach for solid state systems due to the computational efficiency; since we need fairly large supercells for calculating x-ray spectra, this is an obvious advantage.

For treating the core-excited system, we use a frozen-core-hole pseudopotential for the core-ionized atom. The XAS absorption onset equals the binding energy of the core electron relative to the Fermi level,² which for a metallic system is equal to the experimental XPS level. Calculating it requires two single-point DFT calculations: one for the neutral ground state and the other for the core-excited state. Since different pseudopotentials (with and without core hole) are being employed in these two DFT calculations, absolute binding energies cannot be predicted, but the relative shifts due to structural changes can be obtained.

To preserve the charge neutrality of our supercell in the DFT calculations, the core-excited electron is added as an extra valence electron at the Fermi level. Physically, this corresponds to immediate screening by the metal electrons; the electron screening response to the creation of a core hole occurs on a time scale much faster than the core-hole lifetime (6–7 fs for C^24,25). Therefore, using a charged cell to describe adsorbates on metals incorrectly neglects the availability of electrons from the metallic Fermi level, and it has previously been seen that neutral cells indeed do yield better agreement with experiments.¹⁵ 

Instead of introducing an explicit core hole in electronic structure calculations, core-electron binding energy shifts may be obtained using the Z + 1 approximation,^1,26 where the core-excited atom is represented by the next element in the periodic table. This has been found to work well both for gas phase molecules and adsorbates on metal surfaces.²⁷ In our calculations, this amounts to using the pseudopotentials of N and F for C and O, respectively, instead of the explicit core hole.

Calculating XAS requires scattering matrix elements involving unoccupied states far above the Fermi level in addition to the above described computation of the core-electron binding energy. Instead of explicitly calculating the unoccupied band structure, which requires a much larger computational effort than simply converging occupied electronic bands, we calculate the spectrum using a Lanczos recursion Green’s function technique, as implemented in the xspectra code.^28–30 For spectrum calculations, a half-core-hole pseudopotential was used on the excited atom (transition-potential approach¹⁶), which has been shown to give reliable results.^16,31 The spectra were then shifted according to our calculated absorption onsets (Δ-Kohn–Sham method). This gives a better common reference scale for comparing different adsorption sites.³¹ All spectra were finally shifted by a common reference energy, chosen by matching the calculated and experimental c(2 × 2) spectra. By comparison with the experimental c(2 × 2) spectra, a broadening of 0.4 eV FWHM for C and 0.8 eV for O was applied to account for core-hole lifetime, vibrational,¹ and instrumental broadening. The x rays are polarized parallel to the surface, probing the LUMO (2π^*) of the CO molecule, as in the experiments.

The revised Perdew-Burke-Ernzerhof (RPBE) functional was used throughout this study due to its accuracy in determining surface adsorption energies.³² All calculations are spin-polarized (unrestricted Kohn–Sham). For c(2 × 2) phases, we use a four-layer 2 × 2 slab model for initial geometry optimization with 20 Å of vacuum separating the slab from its periodic image. The experimental lattice constant (3.52 Å) is used, and a 4 × 4 × 1 Monkhorst–Pack grid for Brillouin zone integration. For excited-state calculations, we double the cell in the x- and y-directions to 4 × 4 × 4 (correspondingly reducing the k-point grid to 2 × 2 × 1) to reduce the interaction of the core hole with its periodic images. For p(3$2$ × $2$⁠)R45, a 32-atom cell and a 2 × 4 × 1 k-point grid were used for optimization; for c(2$2$ × $2$⁠)R45 23 atoms and 3 × 6 × 1 k-points were used for optimization. For excited states, the supercells of Fig. 1 were used with a k-point grid of 2 × 2 × 1 points. XAS spectra and densities of states were computed using a refined 11 × 11 × 1 k-point grid while keeping the density fixed (non-SCF or Harris calculation). Ultrasoft pseudopotentials³³ were used for atoms without core excitations. For sites where core excitations were considered, ultrasoft half- and full-core-hole pseudopotentials for O, and due to better tested performance, norm-conserving half- and full core-hole pseudopotentials for C were generated using the “atomic” code included with the Quantum Espresso distribution.²² Full- and half-core pseudopotentials are generated by occupying the 1s shell in the single-atom DFT solution by only, respectively, 1 and 1.5 electrons. The constructed 2s-pseudo wavefunctions are nodeless, enabling representation in a plane-wave basis with a reasonable kinetic energy cutoff on the order of 500 eV. The plane-wave and density cutoffs were correspondingly chosen for surface calculations to be 500 eV and 5000 eV, respectively. The geometries were optimized until all forces were less than 0.03 eV/Å with the bottom two layers kept fixed.

The core-level binding energy for ¼ ML coverage was calculated using the GPAW^34,35 code, using PAW:s with an explicit core-hole. The finite-difference mode with a real-space grid spacing of 0.2 Å was used; all other parameters were the same as in the Quantum Espresso calculations described above.

Table I shows the calculated relative binding energies, compared with experimental XPS energy shifts taken from Ref. 18. We use the c(2 × 2) phase as a reference for all other energies. The calculated shifts are consistently within 0.2 eV of the experimental value, which is within the accuracy of typical DFT calculations for this type of surface bond.³² There is no consistent error in the simulations between the different phases, or between the C and O shifts, indicating that the theoretical methods are free of systematic errors regarding the excitation of a certain atom or certain kind of bonding. Co-adsorbing hydrogen represents a significant change in the surface composition and the electronic structure, but the accuracy remains similar. Note that this uncertainty is still fairly large compared to the energy resolution of modern experiments. The good agreement with experimental values indicates that the uncertainty of GGA-level DFT for adsorption energies of CO on transition metals^36,37 does not severely affect our results since the XPS binding energy shifts are significantly larger than the adsorption energy differences [which we find to be at most 70 meV for CO/Ni(100)]. Furthermore, some of the error is expected to be similar in the initial and final states and, thus, cancels out.

The Z + 1 approximation yields remarkably good agreement in the case of core-excited O, where the bond is more ionic and the core-hole effect is largely electrostatic. For C, it also gives accurate results in most cases, but the high degree of hybridization for C^* in the hollow site turns out to be incompletely described in the Z + 1 approximation. This can be clearly seen in the atomic orbital-projected density of states (pDOS), as shown in Fig. 2. While the oxygen levels are shifted in a similar way between the top and hollow sites when using a core-hole pseudopotential and the Z + 1 approximation, the pDOS at the core-ionized carbon atom shows no such simple relation. In particular, the lowest shown (4σ-derived) level has the same energy in the two sites using the Z + 1 approximation, while the more realistic core-hole pseudopotential gives different hybridization in different chemical environments so that the valence electronic structure in the final state is not the same in the Z + 1 vs core-hole case. The electrostatic contribution to the Kohn–Sham effective potential is very similar for the core-hole and Z + 1 approximations and nearly identical at sufficient distance from the nucleus. The exchange-correlation contribution to the potential, however, is more negative near the core in the Z + 1 case due to more core-electronic charge (e.g., we calculate the 2s-state of an isolated N⁺-ion to be ∼2.5 eV lower in energy than for core-hole C). The 2s-states have a non-zero amplitude at the core and are, therefore, affected more strongly than the 2p-states by the lower Kohn–Sham effective potential, which may affect the hybridization of strongly interacting CO in the hollow site.³⁸ This shows that caution may be needed when using the Z + 1 approximation for calculating quantitative shifts in the core-level binding energy in situations where the hybridization of low-lying energy levels takes place.

While both the ¼ ML bridge and hollow sites give similar core-level binding energies as p(3$2$ × $2$⁠) and c(2$2$ × $2$⁠), respectively (except for hollow C Z + 1, for the reasons discussed above), for the top site, the presence of hydrogen gives a significant energy shift compared to the pure c(2 × 2) reference. Disregarding the hydrogen would, thus, underestimate the energy shift between the top and hollow sites, which occurs in the actual experimental c(2$2$ × $2$⁠) co-adsorbed structure. This illustrates the sensitivity of core binding energies to the surface structure and co-adsorbates, and the power of the core-hole pseudopotential DFT approach in discriminating such different structures and identifying the likely observed experimental structure.

Figure 3 shows the calculated XAS at the C and O K-edges, respectively, using the calculated XPS shifts, and compared to the experimental data of Ref. 19. As also seen experimentally, the main resonance narrows when co-adsorbing H, consistent with a weakened chemisorption bond. The double-peak structure of the c(2$2$ × $2$⁠)R45 phase, where the CO molecules occupy two different sites, is also reproduced. However, when using the calculated binding energies to shift the spectrum, the double peak is not as clearly seen in the calculated C K-edge spectrum as in the experiment. This stems from the fact that, as seen in Table I, the experimental C^* shift between top and hollow is 1.1 eV, while the calculated one is only 0.79. In the Z + 1 approximation, the shift is just 0.41 eV, and the double-peak structure would not be resolved at all.

Ultrafast optical pump and x-ray probe experiments offer unique capabilities for gaining insight into reactions on surfaces; however, correct identification of surface species often requires accurate calculated spectra. We have shown that DFT using core-hole pseudopotentials in a plane-wave code gives quite exact core-level binding energy shifts (within 0.2 eV) and x-ray absorption spectra for the experimentally well-studied CO/Ni(100) system. Thus, calculated XPS shifts and x-ray spectrum simulations are powerful tools to identify binding sites from experimental spectra as we have benchmarked here for static experiments on CO/Ni(100). Even if DFT fails to predict small chemisorption energy differences between sites, the XPS shifts are typically larger, which we expect to be useful in guiding the analysis of ultrafast experiments with a dynamic evolution between different binding sites. Resorting to much more expensive approaches to account for explicit many-body effects can likely be avoided as a first level of approximation.

Our shifts are accurate enough to distinguish between different binding sites and reproduce the shape and location of the XAS K-edge. Compared to experiments, the main uncertainty comes from the calculated absorption onset, while the shape of the spectral peaks is very well reproduced. The limitations that turn up in our c(2$2$ × $2$⁠)R45 phase must be kept in mind when comparing different XAS peaks coming from the same atomic species at different sites. Since the uncertainty of about 0.2 eV is independent for the two sites, this method can reliably identify peaks more than ≈0.4 eV apart. This also illustrates the value of any additional experimental information—the O K-edge spectra show a very clear distinction between the peaks, leading to an unambiguous identification of the correct structure. We further see that the Z + 1 approximation, while intuitively appealing and in most cases quantitatively accurate, must be used with caution in situations with the strong covalent bonding and hybridization of low-lying energy levels, where the core-excited state needs to be accurately represented.

This research was supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Chemical Sciences, Geosciences, and Biosciences Division, Catalysis Science Program to the Ultrafast Catalysis FWP 100435, at SLAC National Accelerator Laboratory under Contract Grant No. DE-AC02-76SF00515. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract Grant No. DE-AC02-05CH11231. G.L.S.R. and L.G.M.P. acknowledge support from the Knut and Alice Wallenberg Foundation through Grant No. KAW-2016.0042. The GPAW calculations used resources provided by the Swedish National Infrastructure for Computing (SNIC) at the HPC2N and PDC centers.

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