Kohn-Sham density functional theory calculations of non-resonant and resonant x-ray emission spectroscopy

Kohn-Sham density functional theory calculations of non-resonant and resonant x-ray emission spectroscopy Magnus W. D. Hanson-Heine,1 Michael W. George,1,2 and Nicholas A. Besley1,a) 1School of Chemistry, University of Nottingham, Nottingham NG7 2RD, United Kingdom 2Department of Chemical and Environmental Engineering, University of Nottingham Ningbo China, 199 Taiking East Rd., Ningbo 315100, Zhejiang, China


I. INTRODUCTION
The development of advanced synchrotron sources and free-electron lasers has greatly advanced the capability of spectroscopic techniques in the X-ray region. These techniques provide a local probe of structure and are used in a wide range of areas including materials science and biological chemistry. Furthermore, these methods can contribute to the understanding of ultrafast chemical processes through time-resolved measurements. 1-5 X-ray emission spectroscopy (XES) probes the occupied orbitals and provides a complementary technique to X-ray absorption spectroscopy which probes the unoccupied orbitals. In non-resonant X-ray emission spectroscopy, the energy of the incident photon is sufficient to ionise a core electron, and the subsequent relaxation of the core-ionised state occurs with the emission of an X-ray photon. X-ray emission spectroscopy can also be performed where the incident photon is resonant with an excited core-hole state in resonant inelastic X-ray scattering (RIXS), and recent work has studied the RIXS of small molecules in gas-phase and solution, for example, see Refs. 5-15. The rapid advances in experimental measurements have focused attention on the development of computational methods that can provide accurate X-ray absorption and emission spectra. Several approaches have been established for the simulation of X-ray absorption spectroscopy, including static-exchange, transition potential, and Bethe-Salpeter methods. 16-21 X-ray absorption spectra can also be computed using time-dependent density functional theory (TDDFT) 22-26 as a) E-mail: nick.besley@nottingham.ac.uk well as wavefunction based methods that include electron correlation. [27][28][29][30] It is well established that TDDFT with standard exchange-correlation functionals underestimates core excitation energies, and this has been associated with the self-interaction present with approximate exchange functionals. [31][32][33][34][35][36][37][38] More recently it has been demonstrated that accurate core excitation energies can be computed using TDDFT with short-range corrected (SRC) functionals, which have a large fraction of Hartree-Fock (HF) in the short range. 39 However, these calculations can become computationally demanding for large systems, and it has been shown how the cost of these calculations can be reduced. 40 Similarly, a range of computational methods have also been applied to study X-ray emission spectroscopy. It has been shown that accurate non-resonant X-ray emission spectra can be computed within the framework of equation of motion coupled cluster theory including single and double (EOM-CCSD) excitations. This is achieved by using a reference determinant that describes the core-ionised state, and the X-ray emission transitions appear as negative eigenvalues. 41 While this approach can provide accurate emission energies, its applicability is limited by its computational cost and the difficulty in converging a CCSD calculation for a core-hole state. These problems can be overcome through the use of TDDFT; however, standard exchange-correlation functionals lead to an overestimation of the valence to core transition energies. This approach has been applied to study organic and inorganic systems, and functionals have been parameterised to more accurately reproduce experiment. [42][43][44] One potential limitation of this approach is that inaccurate intensities can arise when the transition is described by a mixture of single excitations within the TDDFT calculation. 44 Alternatively X-ray emission spectra can be computed directly from Kohn-Sham DFT using the following: where φ c is a core orbital with energy c and φ is a valence orbital with energy . This and closely related approaches have been used to study the XES of a wide range of systems. [45][46][47][48] Equation (1.2) shows the oscillator strength computed within the dipole approximation; for heavy nuclei it has been observed that the quadrupole contribution to the transition moment can be significant. 49 Simulating XES directly from the Kohn-Sham calculation is computationally less expensive than TDDFT, since the additional TDDFT calculation is not required and the spectra for all core orbitals can be obtained from the same calculation. The simplicity and low computational cost of obtaining X-ray emission spectra directly from a Kohn-Sham DFT calculation are very useful since it allows large systems to be studied and extensive averaging over conformation to be performed. The latter is particularly relevant for the study of the XES of liquids. 6,50,51 X-ray emission spectra computed using standard gradient corrected or hybrid functionals within this methodology do not predict the correct energy scale, and it is necessary to shift the computed spectra to align with experiment. In this work we show that Kohn-Sham calculations of X-ray emission spectra using SRC functionals designed for the prediction of X-ray absorption spectra provide accurate spectra with the correct energy scale. A current challenge for computational methods is the simulation of RIXS, and we explore the extension of Kohn-Sham DFT to describe RIXS spectra.

II. COMPUTATIONAL DETAILS
A key feature of this work is the application of SRC exchange-correlation functionals. In these functionals the electron repulsion operator is partitioned according to 39 1 r 12 ≡ C SHF erfc(µ SR r 12 ) r 12 − C SHF erfc(µ SR r 12 ) r 12 (2.1) The first and third terms of Equation (2.1) are treated with HF exchange and the remaining terms with DFT exchange leading to the following functional: and × erfc(µ SR r 12 ) r 12 ψ iσ (r 2 )ψ jσ (r 2 )dr 1 dr 2 , (2.4) respectively. The long and short range DFT exchange energies are computed from modifying the exchange energy 52 to give and and for the Becke functional used here 56 with α and β constants and x is the reduced density gradient (x = |∇ρ σ |/ρ 4/3 ). The functional has four parameters, C SHF , C LHF , µ SR , and µ LR , which determine the amount of HF exchange in the short and long ranges. These parameters were optimised to reproduce a set of core-excitation energies. 39 Two sets of parameters were optimised, the first for excitations at the K-edge of first row nuclei and the second for excitations from second row nuclei, and these functionals are denoted SRC1r1 and SRC1r2, respectively. 53 A further closely related functional was also reported, 39 (2.13) If µ SR µ LR , this functional is distinct from the SRC1 functional. This functional has been similarly parameterised for first and second row nuclei to give SRC2r1 and SRC2r2 functionals. The exchange functional of Becke was used 56 for treating the DFT contribution to the exchange energy with the correlation term given by the following combination of LYP 57 and VWN 58 functionals: (2.14) RIXS spectra are simulated by using a reference determinant that describes the core-excited state, in contrast to the core-ionised state for non-resonant XES, and the transition energies and oscillator strengths are computed according to Equations (1.1) and (1.2). In the calculations presented here, the core-ionised and core-excited states are maintained during the self-consistent field (SCF) process using an overlap criterion called the maximum overlap method (MOM). 54,55 This is a two-step procedure where a set of orbitals are generated through a calculation of the ground state, and then the occupancies of the orbitals are specified to describe the excited state for a subsequent calculation. The MOM procedure is applied in the second calculation and should maintain the orbital occupancies for the specified excited state. For the core-ionised and low-lying core-excited states studied here, this approach successfully prevents the collapse of the core-hole during the SCF procedure. However, for higher lying core-excited states, it is possible that a collapse to lower energy state can occur during the SCF process. It is also possible to perform the calculations within a frozen orbital approach, wherein the core-excited state is generated by populating the ground state orbitals accordingly, and no relaxation of the orbitals (i.e., no SCF calculation) is performed for the core-excited state. The frozen orbital approximation greatly simplifies the calculations; however, if changes in the intensity of the emission bands between different core-excited initial states are to be described, it is necessary to allow for orbital relaxation in the intermediate state.
Non-resonant XES and RIXS spectra computed with the SRC functionals are compared with spectra computed with the B3LYP functional, 62,63 since this functional is a representative of standard hybrid functionals. Basis sets of three different sizes are used. The smallest uses 6-31G* for all atom types; the second, medium sized basis set uses the 6-311G** basis set for all atoms except chromium where the Ahlrichs VTZ basis set is used. The largest basis set uses the cc-pwCVTZ basis set for all atoms except hydrogen, for which the cc-pVTZ basis set is used.
Relativistic effects result in the energy of the core orbitals being lowered relative to the valence orbitals. For the heavier nuclei, this energy change is significant and the computed transition energies for the heavier nuclei have been corrected to account for this effect. The magnitude of the shift is computed from the lowering in energy of the 1s orbital between relativistic and non-relativistic HF calculations with the Ahlrichs VTZ basis set, where the relativistic energy was computed with the Douglas-Kroll-Hess Hamiltonian 59 implemented in MOLPRO. 60 This gives energy corrections to the transition energies of +0.6 eV for fluorine in fluorobenzene, +8.9 eV for chlorine in CF 3 Cl, and +33.9 eV for chromium in the metal complexes. For the lighter nuclei, no correction for relativistic effects has been applied. The evaluation of the line intensities according to Equation (1.2) has been implemented in a development version of the Q-Chem 53 package. Spectra have also been computed using TDDFT within the Tamm-Dancoff approximation. 61 To compute core excitations efficiently, these calculations are performed in a restricted single excitation space that includes only excitations from the relevant core orbital(s). 22 This approach has been used to study a wide range of systems and further details can be found elsewhere. 26 Spectra were generated by convoluting with Lorentzian functions with a width of 1.0 eV. A range of different molecules and K-edges have been studied for which experimental data are available, and these are illustrated in Figure 1. All of the molecules have a singlet ground state, except Cp 2 Cr which has a triplet ground state. Figure 2 illustrates the dependence of the Kohn-Sham DFT computed spectra for fluorobenzene (fluorine K-edge) and Cr(CO) 6 (chromium K-edge) on the quality of the basis set. There is little variation in the computed line shapes between the three basis sets used. However, the spectra for the smallest basis set, 6-31G*, are shifted to higher energy by 0.5 eV and 1.5 eV for fluorobenzene and Cr(CO) 6 , respectively. The shift in energy arises from a combination of a lowering in the orbital energies of the valence orbitals and an increase in energy of the core orbital with the larger basis sets relative to the smaller 6-31G* basis set. For the medium and large basis sets, there is no significant difference in the computed transition energies, and the medium sized basis set (6-311G** and Ahlrichs VTZ for Cr) is used for subsequent calculations presented in this work. The sensitivity of the computed transition energies with respect to the size of the basis set is considerably smaller than for TDDFT calculations of XES. 44 Spectra for nitrobenzene, phenol, and fluorobenzene computed from Kohn-Sham DFT and TDDFT are compared in Figure 3. For these calculations the same functional and basis set (B3LYP/6-311G**) have been used. The overall spectral profiles predicted by the two methods are similar, demonstrating some degree of consistency between the two approaches. However, closer inspection of the spectra does reveal some noticeable differences. This is most evident for nitrobenzene where TDDFT predicts two bands with moderate intensity at higher energy than the most intense band while these bands are not evident in the Kohn-Sham DFT spectrum. Also for phenol, there is a small difference in the shape of the most intense band. Factors such as molecular dynamics can affect the spectra (see later); this combined with the relatively low resolution of experimental spectra means that it is not possible to reach a definitive conclusion to which calculation is most accurate. Although it is reassuring that there is reasonable agreement between the two approaches. Table I gives values for the energy of the most intense peak in the non-resonant X-ray emission spectra computed The variation in the computed spectral profile for the three functionals is illustrated in Figure 4 for four of the molecules. For each molecule, the spectra have been shifted to align with the SRC1 spectrum. The plot shows that the computed spectral profile is much less sensitive than the transition energies to the exchange-correlation functional. For three of the molecules there is no significant difference in the computed spectra. The only exception is for phenol where B3LYP predicts greater intensity for the band at 527 eV, which lies to the high energy side of the most intense band. Figure 5 shows a comparison of the SRC1 computed spectra with experiment. All of the spectral profiles are in excellent agreement with experiment. The largest discrepancy is in the low energy bands for nitrobenzene where the calculations underestimate the experimental intensity. For the majority of the molecules, the computed spectra are closely aligned with experiment. For three molecules, methanol, phenol, and CF 3 Cl, the computed spectra are marginally too high in energy. One of the main benefits of calculations of X-ray emission spectra is that it is possible to assign the peaks in the spectra to molecular orbitals. This is illustrated in Figure 6 using the fluorobenzene fluorine K-edge spectrum as an example. This highlights that fact that the most intense peaks are associated with transitions from orbitals with a large p orbital-like component on the fluorine atom.

III. RESULTS AND DISCUSSION
Now we consider the simulation of RIXS spectra where the initial excitation is resonant with a core-excited state. Table II shows computed emission line energies for two low lying initial core-excited states of water and ammonia. We note that for the 1s → 2e initial transition, the degeneracy of the 1e emission lines is lifted and the values in the table represent an average of the corresponding two lines. For these two molecules, gas-phase RIXS data have been reported. 7,14 For the SRC1 exchange-correlation functional, which was the most accurate for non-resonant X-ray emission spectra, the computed transition energies are too low. If there is no relaxation of the orbitals in the core-excited state in the frozen orbital approximation, the predicted transition energies are too high. Interestingly, the experimental value lies approximately half way between the relaxed and frozen orbital SRC1 values. This relationship is likely to be associated with the success of the half-core approximation that is commonly used in the simulation of X-ray absorption spectroscopy. 69 The B3LYP functional in conjunction with the frozen orbital approximation leads to transition energies that are much too high. However, when the orbital relaxation is allowed in the core-excited state, the predicted excitation energies are in good agreement with experiment. The findings regarding which is the more accurate functional are reversed for RIXS compared with non-resonant XES. One of the key differences between the calculation of FIG. 6. Molecular orbitals associated with the observed bands in the fluorine K-edge non-resonant X-ray emission spectrum of fluorobenzene. XES and RIXS is that the core orbital energy in Equation (1.1) corresponds to an occupied orbital for XES and an unoccupied orbital for RIXS. The SRC functionals were parameterised in the context of TDDFT calculations of X-ray absorption spectroscopy, and in this case the core orbitals are occupied. Examination of the unoccupied core orbital energies in the RIXS calculations indicates that the SRC functionals predict an energy that is too high, resulting in an underestimation of the transition energies. This is perhaps not surprising since it is outside the scope for which the functionals were parameterised; it is however somewhat surprising that B3LYP performs so well. Closer inspection shows that the detailed spectral shifts between the two different initial states are not reproduced accurately. For example, in the case of water, experiment shows a small increase in the energy of the 1b 1 emission line for the 1s → 2b 2 initial transition compared with the 1s → 4a 1 initial transition. 7 However, the calculations suggest a shift of opposite sign. This prediction is consistent for both functionals with or without the frozen orbital approximation. One important factor that is neglected in the calculations presented here is the nuclear dynamics. This can be critical for these systems since some of the core-excited states are dissociative. It is possible that nuclear dynamics may be important in reproducing some of these effects, and this is explored in more detail later. Previous theoretical studies of the RIXS of these systems have used the restricted active space self-consistent field method with multiconfigurational perturbation theory to simulate the states. 14 This is considerably a more detailed and extensive theoretical treatment than the one used here. The major advantage of a simple Kohn-Sham DFT treatment is that it adds negligible computational cost, and extensive averaging over molecular conformation from an ab initio molecular dynamics simulation can be performed. These calculations would not be expected to match the accuracy of correlated multi-determinant based methods; however, as shown here they do describe the spectral features observed in experiment reasonably well.
The application of Kohn-Sham DFT to describe RIXS spectra is illustrated in Figure 7 which shows computed B3LYP/6-311G** and experimental RIXS spectra at the carbon K-edge for methanol following excitation to the 3sa and 3pa orbitals. The experimental data are adapted from Ref. 15 and the figure shows spectra computed based upon the ground state structure and also when averaging over conformations taken from an ab initio molecular dynamics simulation in the core-excited (intermediate) state. The MOM procedure which has been shown describe the structure of excited states well, 70 and the molecular dynamics trajectory is propagated with the MOM procedure used to maintain the trajectory on the core-excited state. This approach neglects the possibility of surface crossing occurring in the excited state trajectory. The figure also illustrates the molecular orbitals associated with the observed bands following the 1s → 3sa excitation. Here the orbital labels refer to the orbitals of the ground state.
The computed spectra based upon the ground state structure reproduce the four distinct bands observed in experiment with approximately the correct energy. One significant discrepancy with experiment is that the calculation predicts a shift to lower energy for the 1s → 3pa spectrum while experiment shows no significant shift between the two states. The origin of this shift in the calculations is the energy of the unoccupied core 1s orbital, which is lower in the state arising from the 1s → 3sa excitation. In the experiment there is also a reduction in the intensity associated with the 5a band following the 1s → 3pa excitation. The calculations show a small reduction in this intensity, but the change in intensity is considerably less than observed in experiment. The ab initio molecular dynamics simulations show that both core-excited states are dissociative leading to CH + 3 and OH . The spectra shown represent an average over 25 structures taken at equal intervals from the first 6.5 fs of the dynamics in the core-excited state. In this time period, dissociation has not occurred. The resulting spectra agree more closely with the experimental measurements. In particular, the shift in energy between the spectra for the two different intermediate states is no longer observed. However, the difference in intensity for the 5a band between the two states is not reproduced by the calculations. This may be a consequence of deficiencies in the Kohn-Sham DFT description of these states but may also be related to accurately describing the core-hole lifetimes of the different intermediate states.

IV. CONCLUSIONS
The calculation of non-resonant and resonant (RIXS) X-ray emission spectroscopy based upon Kohn-Sham DFT has been studied. In this approach the transition energy is computed as the difference in the orbital energies and the intensity is proportional to | φ c |μ|φ | 2 . For non-resonant X-ray emission spectroscopy, short-range corrected functionals that have been designed for the prediction of X-ray absorption energies with TDDFT provide accurate spectra with the correct energy scale. This approach can be extended to RIXS through the use of a reference determinant that describes a core-excited state. For this spectroscopy, the standard B3LYP hybrid functional provides accurate spectra with transition energies close to experiment. However, subtle changes between spectra for different intermediate states are not reproduced accurately based upon the single structure and it is necessary to average over conformations from an ab initio molecular dynamics simulation for the core-excited intermediate state. Overall, the results show that surprisingly accurate spectra can be simulated based solely on a Kohn-Sham DFT calculation. This provides a foundation for the study of large systems and the incorporation of nuclear dynamics via averaging over molecular conformation.