Projector augmented wave calculation of x-ray absorption spectra at the L2,3 edges

We develop a technique based on density functional theory and the projector augmented wave method in order to obtain the x-ray absorption cross section at a general edge, both in the electric dipole and quadrupole approximations. The method is a generalization of Taillefumier et al., PRB 66, 195107 (2002). We apply the method to the calculation of the Cu L2,3 edges in fcc copper and cuprite (Cu2O), and to the S L2,3 edges in molybdenite (2H-MoS2). The role of core-hole effects, modeled in a supercell approach, as well as the decomposition of the spectrum into different angular momentum channels are studied in detail. In copper we find that the best agreement with experimental data is obtained when core-hole effects are neglected. On the contrary, core-hole effects need to be included both in Cu2O and 2H-MoS2. Finally we show that a non-negligible component of S L2,3 edges in 2H-MoS2 involves transition to states with s character at all energy scales. The inclusion of this angular momentum channel is mandatory to correctly describe the angular dependence of the measured spectra. We believe that transitions to s character states are quantitatively significant at the L2,3 edges of third row elements from Al to Ar.


I. INTRODUCTION
X-ray absorption near edges structure (XANES) spectroscopy is a very powerful tool for condensed matter studies. Its chemical and spatial selectivity allow solving crystallographic structures and investigating the properties of materials, such as magnetism or orbital hybridization. Unambiguous interpretation of x-ray spectra requires theoretical insight, in order to determine the origin of spectral features and to understand the effect of the core-hole on the measured spectrum.
Several theoretical approaches, either in real space 2,3 or based on periodic boundary conditions, [4][5][6][7][8] were proven successful to describe XANES spectra in molecules and solids at the density functional theory (DFT) level. All-electron based implementations, 5 the most accurate within DFT, suffer from significantly larger computational costs than those required by pseudopotentialbased methods. [6][7][8] In this respect, the advantage of pseudopotential-based methods is that they allow to tackle larger systems, e.g. composed of several hundreds of atoms. In pseudopotential-based methods the atomic core-potential is replaced by a fictitious one, smoother and less singular, hence requiring a smaller representation in terms of plane waves and consequently less computation time. As a consequence, in order to describe properties that depend on core-states, like XANES, one needs to reconstruct the all-electron wavefunction of the empty valence states. This is achieved in the framework of the projector augmented wave (PAW) method. 9 The PAW approach for the calculation of x-ray absorption spectra has been implemented in some freely available (for K-edge only) 6 first-principles codes, commercial code distributions 7,10 and in other non-distributed software. 8,11,12 The PAW method has been proven suc-cessful in interpreting experimental data 8,11,[13][14][15][16][17][18][19][20][21] , thanks to its ability to perform the spectral analysis of the calculated XANES. Most of the applications of the existing methods involve K-edge calculations, where usually core-hole is reliably modeled in a supercell approach and electronic excitations are well described by DFT.
PAW studies at L 2,3 edges are not numerous. 7,22,23 One of the main reasons is that, contrary to K edges, many-electron effects are common at L 2,3 edges. For instance, in the case of L 2,3 edges of transition elements and M 4,5 edges of rare earths, excitations cannot no longer be treated at the DFT level, due to the breakdown of the single-particle picture. 24 Consequently, more sophisticated methods based on time dependent DFT, 25,26 multichannel multiple scattering theory 27,28 or the Bethe-Salpeter equation 29 need to be employed. However, none of the aforementioned first principles methods can rigorously deal with open shells. In this particular case, the ligand field multiplet theory 30 (LFM) can successfully interpret the spectra, albeit not from first principles.
When the states probed by the spectroscopy are delocalized -K and L 1 edges in general (with the exception of very light materials), L 2,3 and M 2,3 edges of elements belonging to the 4d and 5d series -x-ray absorption spectra can usually be interpreted in the DFT framework. In these cases, the final states, mainly of p (K and L 1 ) or d character (L 2,3 and M 2,3 edges) show limited correlation effects and, as a consequence, the structures in a XANES spectrum can be associated to the Kohn-Sham states.
In this work we focus on the description of L 2,3 edges within DFT and PAW frameworks. The previous work in Ref . 1 is generalized to the case of a general edge, for a spin-orbit split core-level. In a first step, we validate our method by calculating Cu L 2,3 edges of metallic Cu. We study the convergence of the spectra with respect to the number of projectors considered in the reconstruction step, i.e. the completeness of the projector basis. It is generally agreed 1,6 that at K edges two linearly independent projectors per angular momentum channel are sufficient as to insure an accurate reconstruction of the 1s state. Albeit a pre-requisite for PAW applications, other edges, and L 2,3 in particular, have not benefit from similar studies.
Next, we switch to the more challenging cases of the Cu L 2,3 edges in cuprite (Cu 2 O) and S L 2,3 edges in molybdenite (2H-MoS 2 ). For both compounds, we perform a detailed analysis of the electronic origin of the spectral features and we assess the effect of the 2p core. Concerning Cu 2 O, we investigate how the empty states probed by the spectroscopy are influenced by the inclusion of the Hubbard U term. In the case of the hexagonal 2H-MoS 2 we interpret the effect of the polarization of the beam, i.e. the x-ray natural linear dichroism (XNLD), by decomposing the spectra into various angular moment channels.
The article is organized as follows. In section II we recall formal aspects of the x-ray cross section calculation in the PAW formalism. In section III we give the technical details of the calculations. Next, in section IV, we present results for the three aforementioned compounds. Finally, in the appendix we include the full analytical expression of the electric dipole (E1) and electric quadrupole (E2) matrix elements for a general edge.

II. X RAY ABSORPTION CROSS SECTION
In the single particle approximation, the x-ray absorption cross section reads: 31 where ψ i is the one-particle wavefunction of the core level and ψ f is the all-electron one-particle (empty) final state. The corresponding electronic energies are E i and E f , ω being the photon energy and α is the fine structure constant. The transition operator is taken as: O =D +Q whereD =ε · r and Q = i 2 (ε · r)(k · r) are the electric dipole and quadrupole operators, respectively. The quantities,ε and k are the polarization vector and the wavevector of the incident x-ray beam, respectively, and r is the position coordinate of the electron. In this work we only deal with the pure dipolar absorption term (E1E1). The quadrupolar contribution (E2E2) is generally negligible at the L 2,3 edges of transition elements, but may become important in rare-earths compounds. E1E2 crossterms appear in optically active materials.
In the next paragraphs we briefly introduce the PAW method applied to the calculation of the x-ray absorption cross section. We closely follow Refs. 1 and 6.
In a pseudopotential calculation the pseudowavefunction of the crystal |ψ f is obtained at the end of the self-consistent field run. In order to get the all-electron wavefunctions |ψ f needed in Eq. 1 all electron reconstruction 9 needs to be performed. As shown in Refs. 1, 6 and 32 the cross section may be written as: In Eq. 3, |φ R0,p (|φ R0,p ) are the all-electron (pseudo) partial waves centered on the absorbing atom. The quantities p R0,p | are a complete set of projector functions, labeled by the index p. The following conditions are satisfied:φ The quantity Ω R is the so-called augmentation (or core) region centered on atomic sites R. In our case, the allelectron partial waves are chosen to be the solutions of the Schrödinger equation for the isolated atom. It is worthwhile recalling that Eq. 2 holds under the assumption that the initial state is localized on the absorbing atom, so that the overlap between the stateÔ|ψ i and |φ R can be neglected if R = R 0 . Furthermore, while Eq. 3 is in principle valid for a complete set of projectors, it converges after a few terms.
In the appendix section A we derive the expression of the transition matrix element φ R0,p Ô ψ i for a general edge.
The single particle interpretation of x-ray absorption is no longer pertinent in the presence of a strong interaction between the photoelectron and the core hole. For the E1E1 (electric dipole) channel, this is the case when the spin-orbit splitting of the initial state is small and the adjacent edges overlap. A notable example are the L 2,3 edges of transition elements 26 or M 4,5 edges of rare earths, no longer obeying the statistical branching ratio. 24 At K edges, neglecting the interaction between the photoelectron and the core hole leads to the overestimation of the energy position of the E2E2 (electric quadrupole) spectral structure. 17 Often, the angular dependence of the E2E2 structure is inappropriately described in the single particle picture. 15

III. TECHNICAL DETAILS
The calculation of a general edge was implemented in the XSpectra code 6 belonging to the Quantum Espresso distribution. 33 The method was applied to the calculation of x-ray absorption spectra of fcc Cu (L 2,3 edges of Cu), Cu 2 O (L 2,3 edges of Cu) and 2H-MoS 2 (L 2,3 edges of S).
For Cu, Mo and O, we used ultrasoft pseudopotentials 37 including semicore states for Cu and Mo. We use a Troullier-Martin 38 norm-conserving pseudopotential for S. Three, respectively two projectors for the l = 2 channel, as well as one projector for the l = 0 channel were needed for the proper reconstruction of core states for the absorbing Cu and S atoms.
Cu and Cu 2 O are described within the generalized gradient approximation (GGA) in the Perdew-Burke-Ernzerhof (PBE) parametrization. 39 In the case of 2H-MoS 2 we used both the local density approximation (LDA) via the Perdew-Zunger (PZ) functional 40 and PBE, without any detectable difference.
The kinetic energy cutoff of the wavefunction expansion is chosen as 30 Ry in Cu, 60 Ry in Cu 2 O, and 50 Ry in 2H-MoS 2 . The charge density cutoff is set ten times larger then the kinetic energy one.
In the case of metallic Cu, integration over electronmomentum is performed by using a uniform 20 × 20 × 20 grid. In the case of insulating Cu 2 O and 2H-MoS 2 we use 20 × 20 × 20 and 10 × 10 × 5 electron-momentum grids, respectively.
The core-hole is described within a static approximation. A pseudopotential with a core-hole is generated for the absorbing atom. Furthermore, a supercell is used to minimize the interaction between the absorbing atom and its periodic images. We use a 3 × 3 × 3 supercell for Cu, a 2 × 2 × 2 one for Cu 2 O and finally a 2 × 2 × 1 supercell for 2H-MoS 2 .
In the XANES calculations, electronic integration was performed on a 30 × 30 × 30 grid for Cu and on a 20 × 20 × 10 grid for 2H-MoS 2 . In the case of Cu 2 O we use the same grid as in the self-consistent calculation. These grids refer to the unit cell when no core-hole is involved. For supercells the grids are rescaled, due to the smaller Brillouin zone.
The spin-orbit coupling on the final (valence) states probed by the spectroscopy is neglected. Although full relativistic calculations are needed to describe heavy transition elements of the late 4d and the 5d series, our approximation is reasonable for the applications presented in this paper. In particular, the spin-orbit coupling acting on the 3d states of transition elements subjected to crystal field is negligible. 41,42 However, the description of the spin-orbit coupling for the core states is mandatory in order to distinguish between related edges as the L 2 (transitions from 2p 1/2 ) and L 3 (2p 3/2 ), M 2 and M 3 or M 4 and M 5 . In the appendix we show the transition matrix elements calculated accordingly.
The radial parts of the 2p 1/2 and 2p 3/2 core wavefunc- The raw spectrum is convoluted with a lorentzian of acrtangent-like, energy-dependent width. EF is the cutoff energy, Γ hole is the spectral width at low energy, above EF and Γmax is the value at high (infinite) energy. The inflection point is situated at Ectr.
tions are taken as identical, a very good approximation in the applications treated in the present paper, as it can be explicitly verified from a DFT relativistic all electron calculation on the isolated atom. We account for the spin-orbit splitting of the 2p states by shifting the calculated spectra by the value of the energy separation between the L 3 and the L 2 edges. This value can be either taken from a DFT relativistic all-electron calculation on the isolated atom or from experiments. The two choices are equivalent when the core-hole can be treated at a DFT level. On the contrary, strong Coulomb interaction between the 2p hole and the valence d electrons may lead to an energy separation between the L 3 and the L 2 which is different from the value of the core level spinorbit splitting. 29 In the absence of spin-orbit coupling on the d states, the L 2 and L 3 edges have similar structures, while the intensities are related by a factor of 2 (the statistical branching ratio), according to the manifold of the initial states. The occupied states are eliminated from the XANES spectra by introducing a cutoff, which appears implicitly in Eq. 1 in the summation on the empty states f . The cutoff is chosen at the calculated Fermi level or at the lowest unoccupied band, depending on the metallic character of the compound. In the presence of a core-hole, the self-consistent calculation is performed on a charged cluster. This means that the energy cutoff is calculated by assuming that the photoelectron does not participate to the screening of the core-hole.
The XANES spectrum being calculated from stationary Kohn-Sham states, a broadening must be applied in order to account for the finite life-time of the electronic levels involved in the transition: The Fermi level or the lowest unoccupied band energy is labeled E F . Subsequently, the raw spectrum σ(ω) is convoluted with a lorentzian having an energy dependent width Γ f : The first term Γ hole is the core-level width, taken as energy-independent since screening is treated statically. In practice we use the standard, edge-dependent, tabulated values. 43 The second term γ(Ω) is the spectral width due to the final state. In general, the spectral structures at the L 2 edge are broader than the ones at L 3 due to the presence of an extra desexcitation channel, the Coster-Krönig effect. 44 Hence, distinct convolution widths are needed at the two edges in order to reproduce the spectral shape. Furthermore, to account for inelastic scattering events (e.g. plasmons), a Lorentzian with a more sophisticated energy dependence is sometimes preferred to a step-like function. We choose the convolution with an arctangent-like function, an empirical model close to the Seah-Dench formalism: 45 where e = (Ω − E F )/(E ctr − E F ). The energy dependent broadening is bound between Γ hole and Γ max , the values at low and respectively, high (infinite positive) energies. The inflection point of the arctangent is situated at E ctr . Figure 1 illustrates the generic form of Γ f (Ω) and the significance of the above parameters. In this work, the following values of Γ hole , Γ max and E ctr − E

A. Cu L2,3 edges in copper
To validate our implementation we first calculate L 2,3 edges for fcc Cu. In fcc Cu, the 3d states of the absorbing atom are almost entirely occupied (atomic Cu is in nominal 3d 10 ). XANES at L 2,3 edges probes the first empty d states just above the main 3d occupied band. These states contribute very weakly to the overall DOS and in a narrow energy region. For this reason, and in agreement with the result in Ref. 47, we find that an extremely dense k-point sampling (electron-momentum grid) is needed to converge the x-ray spectra (see the technical details section).
In Fig. 2 the calculated x-ray absorption (XAS) spectra are compared to experimental data. 46 The effect of the core-hole at the Cu L 2,3 edges in fcc Cu was thoroughly tested. We find that the best agreement between theory and experiment is achieved by neglecting core-hole effects. The global intensity mismatch at the L 2 edge (the calculated intensities are less than the measured one) is most probably due to the background subtraction in the experimental spectrum. The inclusion of the full corehole (not shown) shifts the weight of the d states at several eV below the Fermi level and significantly decreases their contribution above the Fermi level. Even with half core-hole, we find that the agreement with the experiment is worse than in the case without core-hole, with the first peak at the L 3 being less intense than the second one (see Fig. 2).
Our analysis is in partial agreement with the previ- ous interpretation of XANES data on fcc Cu in Ref. 48, based on the Korringa-Kohn-Rostoker method of band structure calculation. The only notable difference between the two calculations is that Ref. 48 produces full relativistic results. In both calculations, the core-hole was not included explicitly. While the positions of the spectral structures agree in the two calculations, it is no longer the case for the intensity of the leading L 3 peak. More precisely, the ratio between the leading peak and the second L 3 structure in the spectrum of Ref. 48 is 1.5 larger than in our calculations. The reason for this disagreement is very likely related to the procedure adopted to eliminate the occupied states from the calculation. We find that the intensity of the first peak in the spectrum strongly depends on this factor, which is often a problem in metallic systems, due to the absence of the gap in the electronic spectrum. As such, the disagreement between the two calculations is not due to fundamental methodological differences in the calculation.
Another difference with Ref. 48 concerns the estimate of the magnitude of transitions from 2p to s states. In Ref. 48 it was estimated to be 5% of the total intensity, whereas they are completely negligible according to our results. The discrepancy may come from treatment of relativistic effects.
Note that the analysis of electron energy loss near edge structure spectroscopy (ELNES) data, 47,49 carried out with a similar formalism to the one used in XANES, leads to the conclusion that a partial half core-hole screening is necessary to explain the experimental data. This apparent disagreement is mostly due to the differences in the experimental XANES and ELNES spectra, and not to our calculation, which actually agrees fairly well with the ones featured in Refs. 47 and 49. Indeed, the first peak L 3 peak in ELNES spectra 47,49 is smaller than the second one, opposite to XANES experimental data 46,48 (see, for instance, Fig. 2).
The fact that we achieve a very good description of the experimental data by completely neglecting the corehole should be interpreted carefully. While one might be tempted to say that in metallic systems the core-hole is almost entirely screened, we recall that at K edges this is not the case. 6,50 It has already been pointed out 51 that it is unlikely to predict whether a XANES calculation with core-hole will agree better to the experiment than the one without hole, due to the inadequacy of the description of the core-hole in the single particle approximation. Nevertheless, one can say that the 2p hole is less profound than the 1s one and consequently less screened. Therefore partial screening is more likely to be needed for the description of L 2,3 edges than at K edges. Fig. 3 features the effect of the convolution procedure on the shape of L 2,3 edge of Cu (see Eq. 7). The inclusion of an energy dependent convolution improves the agreement with the experiment at the L 2 edge. Indeed, a Lorentzian of constant width (Γ hole ) overestimates the resolution at high energies, whereas the arctangent-like energy-dependent Lorentzian can reproduce the correct broadening.
We recall that the sum in Eq. 3 formally involves an infinite number of terms. However, experience demonstrated that at K edges, two linearly independent projectors for the l = 1 channel are enough in order to converge the E1E1 calculated XANES within the first ≈ 50 eV above the edge. 6 The description of the E2E2 contribution sometimes requires three projectors for the l = 2 channel, particularly when the energy of the first projector is set at the one of the semicore states. It is then not obvious that two l = 2 projectors are sufficient for calculations of L 2,3 edges, where the sampled states have mostly d character. Fig. 4 shows the convergence of the L 2 spectrum with respect to the number of projectors with d character considered. In particular, we set the first projector at the energy of the DFT all-electron 3d states of atomic Cu in neutral configuration. The second and third projector are at 2.5 Ry and 10.0 Ry higher energies. This choice leads to linearly independent projectors. As shown in Fig. 4, the use of a single projector underestimate the XANES cross section already at energies just above the main L 2 edge. The use of additional projectors increases the calculated L 2 XANES at higher energies. Three projectors are by far sufficient to describe the full L 2 spectrum 50 eV above the edge.
In Ref. 7 the authors recall that Eq. 2 is not necessarily exact for the 2p core level. Indeed, should the 2p wavefunction no longer be confined within the augmentation region, Eq. 2 is to be replaced with its more general form (see, for instance, Eq. 5 in Ref. 1). However, we find that in the case of Cu both 1s and the 2p core wavefunctions are well localized within the augmentation region, hence the reconstruction according to Eq. 2 is undoubtedly correct, provided that (i) a sufficient Having validated our implementation on fcc Cu, we switch to the more complex case of Cu 2 O. Cu 2 O is a non magnetic insulator where copper is nominally in a Cu +1 , 3d 10 valence state. Density functional theory (DFT) succeeds in capturing the insulating ground state of Cu 2 O, albeit with a smaller gap than the one detected in experiments. Our PBE calculations predict a 0.36 eV direct gap, in agreement with previous calculations, [52][53][54] to be compared to the much larger experimental value, ranging from 2.02 to 2.4 ± 0.3 eV (see Ref. 55

and references therein).
A way to improve the description of the gap is to include the Hubbard U in the calculations, via the DFT+U scheme. While this has been achieved in Ref. 56, the question whether it improves the description of the empty states seen in the x-ray absorption spectroscopy is still open. While the gap increase is very small, the empty states are substantially shifted to higher energy after inclusion of the U term (Fig. 5). This will significantly affect the shape of the calculated XANES spectra, which can be safely interpreted within the single particle picture, given the (nominally) closed shell configuration. Furthermore, by direct comparison with the experiment one can assess whether the DFT+U approximation is appropriate to the description of Cu 2 O.
In order to see the qualitative effect of the Hubbard term, we choose the value of U = 8 eV that is (i) a realistic value for Cu based oxides 6,19 and (ii) allows a substantial shift of spectral features to judge the effect of U in this system. Cu 2 O crystallizes in a cubic structure with four equivalent Cu atoms per unit cell. To calculate the absorption signal, one should consider the contribution of all Cu atoms (one calculation for each absorber). However, given that the four Cu atoms in the unit cell are equivalent, the problem can be simplified by making use of symmetry considerations. Let the atomic E1E1 absorption tensor d defined by with a = b real. The remainder d (j) with j = 2, 3, 4 can be obtained by: where R ij is the 3 × 3 cartesian rotation matrix associated to the symmetry operation relating the prototypical atom i to its equivalent j. In our case the ( 1 4 1 4 1 4 ) Cu atom is related by C 3 rotations to the other Cu atoms. The crystal absorption tensor D can be obtained as the average over all the atomic tensors: In matrix form: Note that D is isotropic, which was to be expected for the absorption tensor of a crystal described by a cubic spacegroup. On the contrary, the atomic absorption tensors  (100) direction. Similar to the case of Cu, we performed XANES calculations with three distinct screening schemes: with corehole, with half core-hole and without core-hole. The results are shown, with or without the inclusion of the Hubbard U , in Figs. 7 (U = 8 eV) and 6 (U = 0 eV) respectively. Calculations were normalized to the edge-jump at high energies and aligned to match the experimental peak at 934 eV. Similarly to the case of Cu, we find that the inclusion of a full core-hole (not shown) shifts the d states completely below the Fermi level and hence yields unphysical results. In Fig. 7 we show that the inclusion of the Hubbard U term in the description of the d states of Cu in Cu 2 O gives a very poor agreement with the experimental data, independently of the screening, contrary to what has been observed in other oxides. 17 We thus bring supplemental evidence to the finding in Ref. 56 according to which GGA+U does not contain the proper physics to describe Cu 2 O. The inclusion of a half core-hole is enough in order to explain the measured XANES at the L 2,3 edges of Cu. In Fig. 6, we show that calculations without core-hole at U = 0 eV overestimate the L 3 structure at 931 eV. On the contrary, a good agreement with the experiments is achieved with a partial screening (half-hole). Indeed, the white line at the L 3 edge is well reproduced in terms of shape and intensity, as well as the positions of the 934, 936 and 943 eV structures. The only remaining source of disagreement with the experiment is related to the underestimation of the background between the two edges.
A deeper analysis of the spectrum reveals that, as expected, the transitions occur to the states formed by the hybridization between the e g of Cu and the p states of O. Although the s states of Cu and O lie in the same energy range as the empty 3d states, they are not seen at the Cu L 2,3 edges. Indeed we find that transitions to states with s character yield a negligible contribution to the spectrum. The 2H polytype of MoS 2 (2H-MoS 2 ) is a layered compound which crystallizes in a hexagonal structure. Each layer is composed of a MoS 2 unit and there are two layers per cell (6 atoms). Layers are bound together by weak Van der Waals interactions.
While isotropic in the case of cubic compounds as Cu and Cu 2 O, the E1E1 x-ray absorption signal depends on the orientation of the polarization in the case of 2H-MoS 2 . We calculated the two independent XANES spectra corresponding to incident light polarized along the principal axis directions, in the plane and parallel to the c axis, at the L 2,3 edges of S. Once again, two distinct sets of calculations have been performed, with and without core-hole.
The comparison with experimental data 58 for the two orientations of the incident beam polarization is shown in Fig. 8. The experimental spectrum atε || c was deduced from the two orientations published in Ref. 58. Calculations including the core-hole are in better agreement with the experiment than the ones without hole, at both orientations. The agreement with the experiment is good, with the position and spectral width of peaks being reproduced at least at the main edge. We note nevertheless that the energy position of the pre-edge peak is slightly overestimated in the theoretical spectra by 1 to 1.5 eV. Second, the intensity of the pre-edge peak is underestimated by the calculation. We believe that both disagreements are a consequence of the inter-edge mixing, 26 due to the fact the S L 2,3 edges are very close in energy (1.1 eV). Third, the intensity of the spectral features above 180 eV is significantly overestimated by the calculations, at both orientations, albeit the energy dependent convolution. We have currently no explanation for this disagreement.
In order to reproduce the quantitative XNLD at the S L 2,3 edges in 2H-MoS 2 , effects of both 2p → s (∆l = −1) and 2p → d (∆l = +1) transitions need to be included in the evaluation of the cross section (Fig. 9). Fig. 9 shows the decomposition of the L 2,3 absorption into the two angular momenta channels, for the two polarization directions. Contrary to the common assumption that at L 2,3 edges transitions to s states are negligible, Fig. 9 proves they have a significant weight in a wide range of energies, albeit lesser than the one of the 2p → d transitions. This behaviour is a consequence of the significant E1 matrix element between the core 2p and valence s states. We expect it to be a general feature of L 2,3 edges of elements belonging to the third period (Al to Ar), as suggested from weighted s and d DOS calculations compared with Si L 2,3 edges in ELNES/XANES spectra of silicon 59 and silicates 60 , as well as from calculations of ELNES at the L 2,3 edges of Ga in GaN. 22 In particular, transitions to s states have been reported previously for XANES at the S L 3 edge in In 2 S 3 , 61 based on a analysis of DOS.
While the difference between the two sets of curves featured in Fig. 9 is due to transitions to s states, one can see that the aforementioned difference depends on the polarization direction albeit the spherical character of the s states. This apparent inconsistency can be understood by noting that, while the ∆l = −1 transitions are indeed isotropic (Fig. 9), the total spectrum also contains cross (with respect to the selection rule), non-negligible ds terms which do have an angular dependence. From Figs. 8 and 9 it is clear that the latter are an essential ingredient for describing the XNLD: in the absence of the cross term, the angular dependence of the 173 and 176 eV structures would nearly disappear.
We draw the attention on a fundamental difference between K and higher rank edges. At K edges, the intensity of XANES structures is directly proportional to the corresponding l projection of the local DOS seen in the spectroscopy -the p-DOS (d-DOS) in the case of E1E1 (E2E2) transitions. This is no longer the case at L 2,3 edges, in the presence of d -s cross terms, when the spectra can no longer be interpreted in terms of pure l-DOS. Note that in the isotropic limit, the d-s interference term is strictly zero. 62 In this sense, the methodology of interpreting polarized XANES spectra in non-cubic samples by confronting them to weighted DOS is not justified.
In this respect, a challenge for the calculation is to disentangle a particular aspect of the spectroscopy, i.e. the interference between the ∆l = 1 and ∆l = −1 channels, from the electronic properties of the material itself. One can recover information about the l-DOS from the spectroscopy calculation and group theory. To begin with, in the following argument we consider for simplicity that (i) there is no spin-orbit coupling on the 2p state and (ii) the local point symmetry at the absorbing S site is described by an abelian group, with real characters.
Defining ψ f | = s| + d|, and |ψ i = |p i where s| and d| are the s and d components of the empty valence states (p being forbidden by the E1 selection rule) and |p i with i = x, y, z is the core state, the E1 matrix element in Eq. 1 is: We focus on the last term in Eq. 10 that represents the d -s interference. It is worthwhile to recall that a term of the form: with α = x, y, z is zero under the above assumptions unless (i) α = i and (ii) d and s belong to the same irreducible representation. This can be understood by noting that, p i and r i being in the same (one-dimensional) irreducible representation γ, γ×γ is the totally symmetric irreducible representation A 1 . The term in (11)  paragraph, can mix. One can easily see that in theε || c direction the cross term is d z 2 -s, while for theε ⊥ c orientation the cross term is due to d x 2 −y 2 -s and d xy -s contributions. This conclusions remains valid if one considers the spin-orbit coupling on the 2p level. Under these circumstances, while the d yz -s and d xz -s contributions are not explicitly forbidden by symmetry, their contribution cancels after summation on the initial states belonging to the 2p 1/2 or 2p 1/2 manifold.
The difference between the d -s interference in the two directions (the linear dichroism) is related to the anisotropy of the d orbitals. More specifically, the linear dichroism of the d -s channel probes the difference between the d z 2 -s term, on one hand, and the joint con- tribution of d x 2 −y 2 -s and d xy -s, on the other hand. This conclusion is equally valid had we considered the true C 3v symmetry. We have already seen in Fig. 9 that the main dichroic effects at 173 and 176 eV are due to the angular dependence of the cross terms. In Fig. 10 we continue the analysis by comparing the calculated XNLD (bottom panel), taken as the difference between theε || c andε ⊥ c spectra, to the cross d -s term in the two orientations (top panel) and its linear dichroism (middle panel). The three insets sharing the same scale, it becomes clear that the XNLD can be safely interpreted in terms of d -s interference. Surprisingly, the pure ∆l = 1 channel has no clear signature in the XNLD, except around 183 eV. While we have shown that the linear dichroism of the cross terms measures the d states anisotropy in 2H-MoS 2 , according to Fig. 10 (middle and bottom panels) this is also the case for the XNLD spectra. One can therefore assign the dichroic feature at 173 eV to d z 2 -s interference, while the feature at 176 eV is attributed to d x 2 −y 2 -s and d xy -s mixing. In spite of the peculiar interference effect, XNLD still measures in-plane / out-of-plane d anisotropy.
On the basis of the angular dependence, the authors of Ref. 58 assign the pre-edge peak and the shoulder at 170 eV to transitions to s states. The decomposition in Fig.  9 contradicts this conclusion: while there exists a contribution of the s states, it is definitely not dominating.
A careful analysis of the DOS (not shown) shows that at the energies corresponding to the pre-edge peak the s states of S hybridize with the d orbitals of Mo. Altogether this proves the importance of calculations in order to accurately interpret the spectroscopic data.
We recall that the interpretation of L 2,3 x-ray absorption spectra in terms of single particle s and d orbitals is not exact in the presence of inter-edge mixing. However, we do not expect it to dramatically affect the main edge, where the probed d states are rather delocalized. This, together with the fact that the agreement with the experimental data is satisfactory, re-confirms the choice of our methodology.
The d -s cross term has also been depicted at the L 2,3 edges of Ag in magnetic multilayers, 63 albeit not in the XANES spectra but in the x-ray magnetic circular dichroism (XMCD). This suggests that a correct description of the cross term is equally important for magnetic spectroscopy. The spectroscopic interpretation would be even richer in the magnetic case, due to the exchange splitting within the 2p 1/2 and 2p 3/2 manifold.

V. CONCLUSIONS
In this work we have developed a technique based on DFT and on the PAW formalism in order to obtain the x-ray absorption cross section at a general edge, in both E1 and E2 approximations. We applied the method to the calculation of Cu L 2,3 edges in fcc Cu and Cu 2 O, as well as to S L 2,3 edges in 2H-MoS 2 . In metallic Cu we find a good agreement with experimental data without having to include a core-hole. On the contrary corehole effects are relevant in insulating Cu 2 O and 2H-MoS 2 , where their inclusion is essential in order to obtain a good description of the experimental data.
In the case of Cu L 2,3 edges in Cu 2 O we equally study the effect of the Hubbard U term on the x-ray absorption spectra. In contrast to results on other oxides, we find that inclusion of U worsens the agreement with experimental data.
Finally in the case of S L 2,3 edges in 2H-MoS 2 , we show that transitions to s states yield a non-negligible contribution to the S L 2,3 spectra in 2H-MoS 2 , in a wide range of energies. We believe that this is a general property of L 2,3 edges of third row elements from Al to Ar. We equally point out that when cross d -s terms are significant, the L 2,3 spectra are no longer to be interpreted in terms of projections of the DOS. The mixing is essential as to interpret the XNLD at S L 2,3 edges in 2H-MoS 2 . In such case it is utmost important to disentangle the two channels, which can only be achieved via first principles calculations.
We have proven that single-particle approaches can be very useful for the interpretation of L 2,3 edges when the inter-edge mixing is negligible. In this case ab initio XANES calculation can assign spectral structures to various angular moment channels, and explore their possible interference.

Appendix A: The E1 and E2 transition matrix elements
For completeness we provide the full analytical expressions of the transition matrix element φ Rp |Ô|ψ i , at a general edge, as it was implemented in the XSpectra module of the Quantum Espresso distribution.
By considering the spin-orbit coupling, the initial state (the core-level) |ψ i is written as: where n i , l i , m i , s = 1/2 and m s are the usual notations for the quantum numbers associated to the initial (core) state in the uncoupled basis and l i m i , sm s |jm j the Clebsch-Gordan coefficients, whose tabulated values can be found in Ref. 64. The functions R nilims (r) and Y m l (r) are the radial and angular (complex spherical harmonics) parts, respectively, of the atomic wavefunction. χ s ms is the m s component of the spinor of spin s. The total momentum quantum numbers j = l i −s, ...l i +s and m j = −j, ..., j define the core-level in the l · s coupled basis. The Clebsch-Gordan coefficients are zero unless the following holds: The all-electron partial waves φ Rp (r) are chosen as solutions of the Schrödinger equation of the isolated atom. As such they are written as where p labels the radial wavefunction at a given energy and n, l, m, s, σ are angular and spin quantum numbers associated to the partial wave.
The electric dipole transition operator can be expanded in the spherical harmonics basis as D =ε · r = 4π 3 r where Ω R is the augmentation region. The electric quadrupole transition operator can be written as: so that the E2 transition matrix element is The generalized Gaunt coefficients are defined as: (Ω) Y m4 l4 (Ω) (A11) and the following relation holds: Note that to calculate the cross section in Eq. 2, the contribution from each initial state ψ i from the 2p 1/2 (2p 3/2 ) must be evaluated independently. Initial states with distinct jm j quantum numbers do not cross in the single particle picture.