X-ray Magnetic and Natural Circular Dichroism from first principles: Calculation of K- and L1-edge spectra

TL;DR

This paper presents a first-principles computational method for calculating X-ray magnetic and natural circular dichroism spectra at K- and L1-edges in solids, incorporating relativistic effects and a new spin-position coupling term.

Contribution

The authors develop an efficient relativistic first-principles approach that includes a novel spin-position operator for accurate XMCD and XNCD spectra calculations, applied to magnetic metals and a non-centrosymmetric compound.

Findings

Electric quadrupole contributions are significant in XMCD spectra.

The spin-position coupling term accounts for 28% of the XMCD amplitude in iron.

The method accurately reproduces angular dependence of XNCD spectra.

Abstract

An efficient first principles approach to calculate X-ray magnetic circular dichroism (XMCD) and X-ray natural circular dichroism (XNCD) is developed and applied in the near edge region at the K-and L1-edges in solids. Computation of circular dichroism requires precise calculations of X-ray absorption spectra (XAS) for circularly polarized light. For the derivation of the XAS cross section, we used a relativistic description of the photon-electron interaction that results in an additional term in the cross-section that couples the electric dipole operator with an operator $\mathbf{\sigma}\cdot (\mathbf{\epsilon} \times \mathbf{r})$ that we name spin-position. The numerical method relies on pseudopotentials, on the gauge including projected augmented wave method and on a collinear spin relativistic description of the electronic structure. We apply the method to the calculations of K-edge XMCD spectra of ferromagnetic iron, cobalt and nickel and of I L1-edge XNCD spectra of $\alpha$-LiIO3, a compound with broken inversion symmetry. For XMCD spectra we find that, even if the electric dipole term is the dominant one, the electric quadrupole term is not negligible (8% in amplitude in the case of iron). The term coupling the electric dipole operator with the spin-position operator is significant (28% in amplitude in the case of iron). We obtain a sum-rule relating this new term to the spin magnetic moment of the p-states. In $\alpha$-LiIO3 we recover the expected angular dependence of the XNCD spectra.

Tables

Table 1. Table 1: Order of magnitude of the operators in Eq. ( 54 ) evaluated at the core state radius r c subscript 𝑟 𝑐 r_{c} compared to the electric dipole operator. The mean radius of core orbitals is deduced from the effective nuclear charge: r c = 3 2 ​ a 0 Z eff subscript 𝑟 𝑐 3 2 subscript 𝑎 0 subscript 𝑍 eff r_{c}=\frac{3}{2}\frac{a_{0}}{Z_{\rm eff}} . Clementi and Raimondi ( 1963 ); Clementi et al. ( 1967 ) B 0 subscript 𝐵 0 B_{0} has been fixed to 2 × 10 4 2 superscript 10 4 2{\times}10^{4} T (1.2eV) which is two order of magnitude larger than the exchange splitting observed for Fe K -edge. The Coulomb potential is V = − Z eff ​ e 4 ​ π ​ ϵ 0 ​ r c 𝑉 subscript 𝑍 eff 𝑒 4 𝜋 subscript italic-ϵ 0 subscript 𝑟 𝑐 V=\frac{-Z_{\rm eff}e}{4\pi\epsilon_{0}r_{c}} and the core state energy E i subscript 𝐸 𝑖 E_{i} is evaluated in a planetary model E i = − Z eff ​ e 2 8 ​ π ​ ϵ 0 ​ r c subscript 𝐸 𝑖 subscript 𝑍 eff superscript 𝑒 2 8 𝜋 subscript italic-ϵ 0 subscript 𝑟 𝑐 E_{i}=\frac{-Z_{\rm eff}e^{2}}{8\pi\epsilon_{0}r_{c}} .

Edge		L1	L2			K
		I	Fe	Gd	Bi	O	Fe
Energy (keV)		5.19	0.72	7.898	15.71	0.53	7.11
$Z_{\mathrm{eff}}$		39.067	22.089	29.8527	39.2335	7.6579	25.381
$T_Q$ (43)	$k{r}_c/2$	$2.7\times {10}^{-2}$	$6.6\times {10}^{-3}$	$5.3\times {10}^{-2}$	$8.1\times {10}^{-2}$	$1.4\times {10}^{-2}$	$5.7\times {10}^{-2}$
$T_{\mathrm{SP}}$ (48)	$\hslash \omega /4m{c}^2$	$2.6\times {10}^{-3}$	$3.5\times {10}^{-4}$	$3.9\times {10}^{-3}$	$7.7\times {10}^{-3}$	$2.6\times {10}^{-4}$	$3.5\times {10}^{-3}$
$T_{A_0}$ (45)	$ek{r}_c{B}_0/4m\omega$	$6.0\times {10}^{-6}$	$1.1\times {10}^{-5}$	$7.8\times {10}^{-6}$	$6.0\times {10}^{-6}$	$3.0\times {10}^{-5}$	$9.2\times {10}^{-6}$
$T^e$ (53)	$\left(E_i-eV\right)/2m{c}^2$	$1.3\times {10}^{-2}$	$4.3\times {10}^{-3}$	$7.9\times {10}^{-3}$	$1.3\times {10}^{-2}$	$5.2\times {10}^{-4}$	$5.7\times {10}^{-3}$

Equations

$${\sigma }^{\mathrm{CD}}=\sigma ({\mathbf{ϵ}}_2)-\sigma ({\mathbf{ϵ}}_1)$$

$$\sigma =4{\pi }^2{\alpha }_0\hslash \omega \sum _f|⟨f|T|i⟩{|}^2\delta (E_f-E_i-\hslash \omega )$$

$$\begin{gathered} H^{\mathrm{FW}}=m{c}^2+\frac{{\mathbf{p}}^2}{2m}+eV-\frac{e\hslash }{2m}\mathbf{\sigma }\cdot {\mathbf{B}}_0 \\ -\frac{e\hslash }{4m^2{c}^2}\mathbf{\sigma }\cdot ({\mathbf{E}}_0\times \mathbf{p})-\frac{e{\hslash }^2}{8m^2{c}^2}\nabla \cdot {\mathbf{E}}_0. \end{gathered}$$

$$T=\mathbf{ϵ}\cdot \mathbf{r}+\frac{i}{2}\mathbf{ϵ}\cdot \mathbf{r}\mathbf{k}\cdot \mathbf{r}+\frac{i\hslash \omega }{4m{c}^2}\mathbf{\sigma }\cdot (\mathbf{ϵ}\times \mathbf{r})$$

$${\sigma }_{\mathrm{D}-\mathrm{D}}=4{\pi }^2{\alpha }_0\hslash \omega \sum _f|⟨f|\mathbf{ϵ}\cdot \mathbf{r}|i⟩{|}^2\delta (E_f-E_i-\hslash \omega ).$$

$$\begin{gathered} {\sigma }_{\mathrm{Q}-\mathrm{Q}}={\pi }^2{\alpha }_0\hslash \omega \\ \sum _f|⟨f|(\mathbf{k}\cdot \mathbf{r})(\mathbf{ϵ}\cdot \mathbf{r})|i⟩{|}^2\delta (E_f-E_i-\hslash \omega ). \end{gathered}$$

$$\begin{gathered} {\sigma }_{\mathrm{D}-\mathrm{Q}}=-4{\pi }^2{\alpha }_0\hslash \omega \\ \sum _f\Im [⟨f|(\mathbf{k}\cdot \mathbf{r})(\mathbf{ϵ}\cdot \mathbf{r})|i⟩⟨i|{\mathbf{ϵ}}^{\star }\cdot \mathbf{r}|f⟩]\delta (E_f-E_i-\hslash \omega ). \end{gathered}$$

$$\begin{gathered} {\sigma }_{\mathrm{D}-\mathrm{SP}}=-\frac{2{\pi }^2{\alpha }_0{\hslash }^2{\omega }^2}{m{c}^2} \\ \sum _f\Im [⟨f|\mathbf{\sigma }\cdot (\mathbf{ϵ}\times \mathbf{r})|i⟩⟨i|{\mathbf{ϵ}}^{\star }\cdot \mathbf{r}|f⟩]\delta (E_f-E_i-\hslash \omega ). \end{gathered}$$

$$\mathcal{T}=11+\sum _{\mathbf{R},n}(|{\varphi }_{\mathbf{R},n}⟩-|{\tilde{\varphi }}_{\mathbf{R},n}⟩)⟨{\tilde{p}}_{\mathbf{R},n}|.$$

$$\tilde{\mathcal{H}}=E_{\text{kin}}+e{\tilde{V}}^{\text{loc}}(\mathbf{r})+\sum _{\mathbf{R}}e{\tilde{V}}_{\mathbf{R}}^{\text{nl}}+{\tilde{\mathcal{H}}}_{\text{SO}}$$

$${\tilde{\mathcal{H}}}_{\text{SO}}$$

$$=\frac{e\hslash }{4m^2{c}^2}\left(\mathbf{\sigma }\cdot (\nabla {\tilde{V}}^{\text{loc}}(\mathbf{r})\times \mathbf{p})+\sum _{\mathbf{R}}{F}_{\mathbf{R}}^{\text{nl}}\right).$$

$$\begin{gathered} F_{\mathbf{R}}^{\text{nl}}=\sum _{n,m}|{\tilde{p}}_n^{\mathbf{R}}⟩\mathbf{\sigma }\cdot (⟨{\varphi }_{\mathbf{R},n}|\nabla v_{\mathbf{R}}(\mathbf{r})\times \mathbf{p}|{\varphi }_{\mathbf{R},n}⟩ \\ -⟨{\tilde{\varphi }}_{\mathbf{R},n}|\nabla {\tilde{v}}_{\mathbf{R}}^{\text{loc}}(\mathbf{r})\times \mathbf{p}|{\tilde{\varphi }}_{\mathbf{R},n}⟩)⟨{\tilde{p}}_m^{\mathbf{R}}| \end{gathered}$$

$$\begin{gathered} F_{\mathbf{R}}^{\text{nl}}=\sum _{n,m}|{\tilde{p}}_n^{\mathbf{R}}⟩\mathbf{\sigma }\cdot (⟨{\varphi }_{\mathbf{R},n}|\frac{1}{r}\frac{\partial v_{\mathbf{R}}}{\partial r}\mathbf{L}|{\varphi }_{\mathbf{R},n}⟩ \\ -⟨{\tilde{\varphi }}_{\mathbf{R},n}|\frac{1}{r}\frac{\partial {\tilde{v}}_{\mathbf{R}}^{\text{loc}}}{\partial r}\mathbf{L}|{\tilde{\varphi }}_{\mathbf{R},n}⟩)⟨{\tilde{p}}_m^{\mathbf{R}}|. \end{gathered}$$

$$\frac{e\hslash }{4m^2{c}^2}\sum _{n\mathbf{R}m}\mathbf{\sigma }\cdot |{\tilde{p}}_n^{\mathbf{R}}⟩⟨{\varphi }_{\mathbf{R},n}|\frac{1}{r}\frac{\partial v_{\mathbf{R}}}{\partial r}\mathbf{L}|{\varphi }_{\mathbf{R},m}⟩⟨{\tilde{p}}_m^{\mathbf{R}}|.$$

$$\mathbf{\sigma }={\sigma }_z{\mathbf{e}}_z.$$

$$\sigma (\omega )=4{\pi }^2{\alpha }_0\hslash \omega \sum _f|⟨f|O|i⟩{|}^2\delta (E_f-E_i-\hslash \omega )$$

$$\sigma (\omega )=4{\pi }^2{\alpha }_0\hslash \omega \sum _f|⟨\tilde{f}|{\tilde{\phi }}_{{\mathbf{R}}_0}⟩{|}^2\delta (E_f-E_i-\hslash \omega )$$

$$|{\tilde{\phi }}_{{\mathbf{R}}_0}⟩=\sum _n|{\tilde{p}}_n^{{\mathbf{R}}_0}⟩⟨{\varphi }_n^{{\mathbf{R}}_0}|O|i⟩.$$

$$\Im [D{B}^{\star }]=\frac{1}{4}(|D+iB{|}^2-|D-iB{|}^2)$$

$$⟨{\varphi }_{i}s|{\mathbf{ϵ}}^{\star }\cdot \mathbf{r}|{\varphi }_f{s}^{\prime }⟩=⟨{\varphi }_i|{\mathbf{ϵ}}^{\star }\cdot \mathbf{r}|{\varphi }_f⟩{\delta }_{s{s}^{\prime }}.$$

$$⟨{\varphi }_{i}s|\mathbf{\sigma }\cdot (\mathbf{ϵ}\times \mathbf{r})|{\varphi }_{f}s⟩=⟨{\varphi }_i|(\mathbf{ϵ}\times \mathbf{r})|{\varphi }_f⟩\cdot ⟨s|\mathbf{\sigma }|s⟩$$

$$S{P}_{\mathrm{col}}(\mathbf{ϵ})$$

$$={\sigma }_z\frac{4i\pi }{3}r(Y_1^{-1}(\mathbf{ϵ})Y_1^1({\mathbf{u}}_{\mathbf{r}})-Y_1^1(\mathbf{ϵ})Y_1^{-1}({\mathbf{u}}_{\mathbf{r}})).$$

$$S{P}_{\mathrm{col}}({\mathbf{ϵ}}_1)$$

$$S{P}_{\mathrm{col}}({\mathbf{ϵ}}_2)$$

$${\sigma }_{\mathrm{D}-\mathrm{SP}}({\mathbf{ϵ}}_1)$$

$${\sigma }_{\mathrm{D}-\mathrm{SP}}({\mathbf{ϵ}}_2)$$

$$\begin{gathered} {\sigma }_{\mathrm{D}-\mathrm{D}}^s(\mathbf{ϵ})=4{\pi }^2{\alpha }_0\hslash \omega \\ \sum _f|⟨f^s|\mathbf{ϵ}\cdot \mathbf{r}|i^s⟩{|}^2\delta (E_f-E_i-\hslash \omega ) \end{gathered}$$

$$\begin{gathered} {\sigma }_{\mathrm{D}-\mathrm{SP}}(\mathbf{ϵ})=\frac{2{\pi }^2{\hslash }^2{\alpha }_0{\omega }^2}{m{c}^2}\sum _{\nu =-1}^1\Re [Y_1^{-\nu }({\mathbf{ϵ}}^{\star }) \\ \left(Y_1^1(\mathbf{ϵ}){\zeta }_{\mathrm{D}-\mathrm{SP}}^{1\nu }-Y_1^{-1}(\mathbf{ϵ}){\zeta }_{\mathrm{D}-\mathrm{SP}}^{-1\nu }\right)] \end{gathered}$$

$$\begin{gathered} {\zeta }_{\mathrm{D}-\mathrm{SP}}^{\lambda \nu }=(-1{)}^{\nu }{\left(\frac{4\pi }{3}\right)}^2\sum _f⟨f|\sum _i({{\sigma }_z}_i{r}_i{Y}_1^{\lambda }({\mathbf{u}}_{{\mathbf{r}}_{\mathbf{i}}}){)}^{\star }|g⟩ \\ ⟨g|\sum _i{r}_i{Y}_1^{\nu }({\mathbf{u}}_{{\mathbf{r}}_{\mathbf{i}}})|f⟩\delta (E_f-E_g-\hslash \omega ). \end{gathered}$$

$$\begin{gathered} ⟨g|\sum _i{r}_i{Y}_1^{\lambda }({\mathbf{u}}_{{\mathbf{r}}_{\mathbf{i}}})|f⟩=\sum _{lm\sigma l_0{m}_0{\sigma }_0^{\prime }}\sqrt{\frac{3(2l+1)}{4\pi (2l_0+1)}} \\ (10l0|l_{0}0)(1\lambda lm|l_0{m}_0)⟨g|a_{l_0{m}_0\sigma }^{†}{a}_{lm\sigma }|f⟩{\mathcal{D}}_{l_0,l} \end{gathered}$$

Full text

I Introduction

A dichroic (“two-coloured” in Greek) material has the property to absorb light differently depending on its polarization. X-ray Circular Dichroism is the difference between X-ray absorption spectra (XAS) obtained from left and right circularly polarized light so it describes the dependence of the absorption cross-section on the state of circularly-polarized light.

In a magnetic sample, the breaking of time-reversal symmetry permits X-ray Magnetic Circular Dichroism (XMCD). It is a powerful tool for the study of the magnetic structure of complex systems as it gives element-specific information. Almost all synchrotron facilities around the world have a beamline dedicated to XMCD. Rogalev and Wilhelm (2015) The existence of well-established magneto-optical sum-rules that allow to obtain the spin and orbital contribution to the magnetic moment directly from the integral of the spectra Thole et al. (1992); Carra et al. (1993); Altarelli (1993) made it an essential technique to study the magnetic properties of matter. These sum-rules are widely and successfully applied at spin-orbit split L2,3-edges of transition metals Vogel et al. (1997); Stöhr (1999); Edmonds et al. (2005); Prado et al. (2013) and M4,5-edges of actinides. Wilhelm et al. (2013) On the other hand, in the absence of spin-orbit splitting of the core state (as for K- or L1-edge), only the orbital magnetization sum-rule Thole et al. (1992); Altarelli (1993) can apply and a quantitative analysis of the spectra is far from being straightforward. Yet, for 3$d$ transition elements, measurements of XMCD at the K-edge is the main way to probe magnetism under pressure and it is a widely used technique despite the interpretation difficulties. Baudelet et al. (2005); Torchio et al. (2011a); Cafun et al. (2013)

X-ray Natural Circular Dichroism (XNCD) occurs in non-centrosymmetric materials (for which the inversion symmetry is not a symmetry of the system). Up to now, it has been less widely used than XMCD but it presents a fundamental interest as it gives access to element specific stereochemical information. Rogalev et al. (2010) In the domain of molecular magnetism, a renewed interest for this technique has recently grown Sessoli et al. (2015) with the emergence of new materials that are both chiral and magnetic. Contrary to optical activity to which a large number of mechanisms contribute,Goulon et al. (1998) XNCD is largely dominated by a single contribution.Rogalev et al. (2010) At L1- and K-edges, XNCD exists only if $p$ and $d$ orbitals are mixed,Natoli et al. (1998) yielding a unique measure of the mixing of even and odd orbitals.

The starting point of our work is a Density Functional Theory (DFT) based pseudopotential method. Using Projector Augmented Wave (PAW) reconstruction, Lanczos algorithm and a continued fraction calculation, Taillefumier et al. (2002); Gougoussis et al. (2009); Giannozzi et al. (2009) it has proved successful for the calculation of absorption (XAS) spectra at the K-edge. Taillefumier et al. (2002); Gougoussis et al. (2009); Cabaret et al. (2010); Bordage et al. (2010) The L1-edge, that corresponds to a 2$s$ core-hole, is expected to have the same behavior. In this paper, we propose the same kind of DFT-based approach for the calculation of XMCD and XNCD spectra in the near-edge (XANES) region.

Several calculations of XMCD at K-edge in the XANES region can be found in the literature. Most of these calculations are based on fully relativistic Ebert et al. (1988); Ebert (1996); Stähler et al. (1993); Gotsis and Strange (1994); Guo (1996); Šipr and Ebert (2005) or semi-relativistic Brouder and Hikam (1991); Brouder et al. (1996) multiple-scattering approaches with muffin-tin potentials even if efforts have been made to go beyond this approximation. Natoli et al. (1986); Joly et al. (2009)

The technique presented in this paper allows the use of a free-shape potential. Relativistic perturbations were taken into account both in the band structureGerstmann et al. (2014) and in the photon matter interaction.Bouldi and Brouder (2016) The method has been implemented within a highly efficient reciprocal space code that allows the modelling of a large range of systems.Giannozzi et al. (2009)

In section II, the terms that enter the absorption cross section up to the electric quadrupole approximation are listed. Section III is dedicated to the presentation of the computational method. Results obtained for XAS and XNCD at L1-edge of iodine in both enantiomers of $\alpha$-LiIO3 and for K-edge XAS and XMCD spectra in 3$d$ ferromagnetic metals are presented in section IV. Finally, in section V, the relativistic operator is examined in detail within the collinear spin approximation. Its corresponding sum-rule is derived and evaluated numerically and an expression that allows for a simple implementation of this term is given.

II Contributions to the cross section

In the case of a fully circularly polarized light with a wave vector $\mathbf{k}$ along $z$, the circular dichroism (XMCD and XNCD) cross section writes:

[TABLE]

where ${\boldsymbol{\epsilon}}_{2}=1/\sqrt{2}(1,i,0)$, ${\boldsymbol{\epsilon}}_{1}={\boldsymbol{\epsilon}}_{2}^{\star}=1/\sqrt{2}(1,-i,0)$ and $\sigma({\boldsymbol{\epsilon}})$ is the X-ray absorption (XAS) cross section of the material. XMCD effect at K-edge of 3$d$ transition elements results at most in an asymmetry in absorption of the order of $10^{-3}$. For this study, it is therefore important to compute the absorption cross section very accurately.

In a monoelectronic semi-relativistic framework the contribution to the XAS cross section from a given core-state of energy $E_{i}$ is given by (see appendix):

[TABLE]

where $\alpha_{0}$ is the fine structure constant, $|i\rangle$ is the two-components wave function that corresponds to the large components of the Dirac wave function of the core state and the sum runs over unoccupied final states with energy $E_{f}$. The wave functions $|f\rangle$ are eigenstates of the time-independent Foldy-Wouthuysen Hamiltonian of the electron in the presence of an electromagnetic field $\mathbf{E}_{0},\mathbf{B}_{0}$:Itzykson and Zuber (1980); Bjorken and Drell (1964)

[TABLE]

Finally, $T$ is the sum of three operator: (i) the electric dipole operator, (ii) the electric quadrupole operator and (iii) a new light-matter interaction term that we named the spin-position operator:

[TABLE]

where ${\boldsymbol{\sigma}}$ is the vector of Pauli matrices.

The absorption cross section expands in six terms among which four terms are significant (see the orders of magnitude in appendix).

The dominant term is the electric dipole-electric dipole (D-D) term:

[TABLE]

It is usually the only term that is taken into account in calculations of XAS and XMCD spectra at the L2,3-edges and sometimes at the K-edge. Rehr and Albers (2000); Natoli et al. (1980); Brouder et al. (1996); Fujikawa and Nagamatsu (2003)

The electric quadrupole-electric quadrupole (Q-Q) term is:

[TABLE]

At the K-edge, it can reach a few percent of $\sigma_{\mathrm{D-D}}$ in amplitude. It contributes mainly to the pre-edge region. It is sometimes included in X-ray absorption calculations. Taillefumier et al. (2002); Bunǎu and Joly (2009)

When neglecting spin-orbit coupling and in the absence of an external magnetic field, it is possible to choose real wave functions. In that case, the D-D and Q-Q terms verify $\sigma({\boldsymbol{\epsilon}})=\sigma({\boldsymbol{\epsilon}}^{*})$, which leads to a zero contribution to circular dichroism. For this reason it is crucial to account for relativistic effects in the wave functions calculation in order to compute XMCD.

On the other hand, the two following terms can give a non-vanishing contribution to the circular dichroism cross section even when wave functions can be chosen real.

The electric dipole-electric quadrupole cross term (D-Q) is:

[TABLE]

If $|i\rangle$ and $|f\rangle$ are parity invariant (i.e. if inversion $\mathbf{r}\rightarrow-\mathbf{r}$ is a symmetry of the system) then $\sigma_{\mathrm{D-Q}}=0$. It is however this term that is responsible for XNCD Natoli et al. (1998) because the electric dipole-magnetic dipole term (that is responsible for optical activity in the optical range) is very small in the X-ray range.

The cross term between the electric dipole and the relativistic operator that we named spin-position (D-SP) is:

[TABLE]

It exists only in magnetic materials. Like the spin-orbit coupling term in the FW Hamiltonian, it arises from the coupling of the small components of the Dirac wave functions. To our knowledge, it has never been evaluated before. We will show in the following that, despite the small prefactor of this term, its contribution to XMCD at the K-edge of 3d metals can account for up to one third of the XMCD intensity near the edge.

III Method

In the framework of the final state rule von Barth and Grossmann (1982) the absorption cross section is obtained from one-electron wave functions. Within the frozen core approximation, the 1$s$ (K-edge) or 2$s$ (L1-edge) unperturbed core states $|i\rangle$ can be determined from an all-electron isolated atom calculation. The stationary final states $|f\rangle$ are calculated self-consistently in the presence of a core hole. Here, they are calculated within a semi-relativistic pseudopotential based DFT and Projector Augmented Wave (PAW) reconstruction framework.Gerstmann et al. (2014) The absorption cross-section is then calculated in a continued fraction scheme using Lanczos algorithm.Taillefumier et al. (2002); Gougoussis et al. (2009)

III.1 Collinear semi-relativistic self-consistent field calculation

Self-consistent field calculations in this study are based on Density Functional Theory (DFT) with a plane-wave basis set, and pseudopotentials as implemented in Quantum ESPRESSO Giannozzi et al. (2009) including the spin-orbit coupling (SOC) term.Gerstmann et al. (2014) Since an accurate implementation of SOC plays a crucial role for the evaluation of XMCD spectra, we briefly describe the underlying approach in the following.

In pseudopotential-based methods the potential near the nuclei is replaced by a fictitious smooth potential. The valence electrons wave functions are replaced by pseudo-wave functions that are exempt from the rapid oscillations near the core. The size of the plane-waves basis set needed to describe the system is therefore considerably lowered which leads to a much better computational efficiency compared to an all-electron approach making possible an ab initio description of large systems with thousands of electrons.

In the PAW formalism, as described by Blöchl,Blöchl (1994) the physical valence wave functions $|{\Psi}\rangle$ can be reconstructed from the pseudo-wave functions $|{\tilde{\Psi}}\rangle$ as they are related through a linear operator ${\cal{T}}$: $|{\Psi}\rangle={\cal{T}}|{\tilde{\Psi}}\rangle$ with

[TABLE]

In our case, the set of all-electron partial waves centered on atomic site $\mathbf{R}$, $|\phi_{\mathbf{R},n}\rangle$, are solutions of the Dirac equation for the isolated atom within a scalar relativistic approximation,Koelling and Harmon (1977) $|\tilde{\phi}_{\mathbf{R},n}\rangle$ are the corresponding pseudo-partial waves and $\langle\tilde{p}_{\mathbf{R},n}|$ form a complete set of projector functions. The operator $\cal{T}$ only acts in augmentation regions enclosing the atoms. Outside the augmentation regions the all-electron and pseudo-wave functions coincide.

The pseudo-Hamiltonian is given by ${\cal{T}}^{\dagger}H^{\mathrm{FW}}{\cal{T}}$:Ceresoli et al. (2010); Gerstmann et al. (2014)

[TABLE]

where $E_{\text{kin}}$ is the kinetic energy as implemented in Quantum ESPRESSO and $\tilde{V}^{\text{loc}}$ and $\tilde{V}_{\mathbf{R}}^{\text{nl}}$ are the local part and the nonlocal part in separable form of the pseudopotentials. $\tilde{{\cal{H}}}_{\text{SO}}$ is the pseudo-Hamiltonian corresponding to the time independent spin-orbit term in the Foldy-Wouthuysen transformed Hamiltonian:Ceresoli et al. (2010)

[TABLE]

The $F_{\mathbf{R}}^{\text{nl}}$ at the atomic site $\mathbf{R}$ are: Pickard and Mauri (2002)

[TABLE]

where $v_{\mathbf{R}}$ and $\tilde{v}_{\mathbf{R}}^{\text{loc}}$ are the atomic all-electron and local channel pseudopotentials respectively. As these potentials are spherical, $F_{\mathbf{R}}^{\text{nl}}$ rewrites:

[TABLE]

The local potential $\tilde{V}^{\text{loc}}(\mathbf{r})=\sum_{\mathbf{R}}\tilde{v}_{\mathbf{R}}^{\text{loc}}(\mathbf{r})$ and the quantity $\frac{1}{r}\frac{\partial\tilde{v}_{\mathbf{R}}^{\text{loc}}}{\partial r}$ decreases in $1/r^{3}$ so that the action of the operator $\tilde{V}^{\text{loc}}(\mathbf{r}){\times}\mathbf{p}$ in the augmentation region is, at first order, the same as the action of $\nabla\tilde{v}_{\mathbf{R}}^{\text{loc}}(\mathbf{r}){\times}\mathbf{p}$. In the PAW framework any pseudo-wave function in the augmentation region can be expanded according to $|\tilde{\Psi}\rangle=\sum_{n}|\tilde{\phi}_{n,\mathbf{R}}\rangle\langle\tilde{p}_{n}^{\mathbf{R}}|\tilde{\Psi}\rangle$. Therefore, the term proportional to $\tilde{v}_{\mathbf{R}}^{\text{loc}}$ and the term proportional to $\tilde{V}^{\text{loc}}(\mathbf{r})$ partially compensate each other so that the dominant contribution arises from the term:

[TABLE]

In this study, we consider collinear spin along $z$ and only the $z$ Pauli matrix is considered (diagonal spin-orbit coupling approximation):

[TABLE]

In XMCD experiments a magnetic field is usually applied parallel to the beam,Rogalev et al. (2006) which justifies to consider the quantization axis parallel to $\mathbf{k}$.

This semi-relativistic approach, that includes spin-orbit coupling in a two-component approach, is computationally less expensive than a fully relativistic one. It has been shown to reproduce the fully relativistic band structure. Gerstmann et al. (2014) For heavy atoms, the formula can be generalized by substituting $\nabla\tilde{V}^{\text{loc}}$ and $\frac{\partial\tilde{v}_{\mathbf{R}}^{\text{loc}}}{\partial r}$ with reduced gradients, resulting in a ZORA-type of Hamiltonian.Gerstmann et al. (2014)

In this study, the calculations have been performed using Troullier-Martins norm-conserving pseudopotentials and are based on the generalized gradient approximation (GGA) with PBE functionals.Perdew et al. (1996) The charge density is evaluated self-consistently in the presence of a core hole which is described by removing a 1s or 2s electron in the pseudopotential of the absorbing atom. A large unit cell (supercell) must be built to minimize the interactions between periodically reproduced core-holes and the k-points grid can be reduced accordingly.

III.2 Cross-section calculation

We implemented XMCD and XNCD in the Xspectra codeGougoussis et al. (2009) of Quantum ESPRESSOGiannozzi et al. (2009) distribution. The first results of this implementation for the terms D-D and Q-Q can be found in Ref. Gougoussis, 2009.

In the PAW formalism it has been shown Taillefumier et al. (2002); Gougoussis et al. (2009) that the contribution of the operator $O$ to the absorption cross-section,

[TABLE]

can be rewritten, as the initial wave function is localized around the absorbing atoms $\mathbf{R}_{0}$,

[TABLE]

with

[TABLE]

This sum involves in principle an infinite number of projectors but experience demonstrated that two or three linearly independent projectors are in general sufficient in order to achieve the convergence of the D-D term at the K-edge in the XANES region.Bunǎu and Calandra (2013)

The determination of all empty states in Eq. (17) would require a lot of computing resources and, as a consequence, would limit the size of the manageable supercell. To increase the efficiency of the method, the cross section is evaluated as developed in Ref. Taillefumier et al., 2002 and Gougoussis et al., 2009 via the Green’s function using Lanczos algorithmLanczos (1952) which avoids the heavy workload of a large matrix inversion. The cross terms D-SP and D-Q are not in the form of Eq. (16) but they can be determined from two calculations of this type using the relationship:

[TABLE]

where $B$ is either the electric quadrupole or the spin-position operator and $D$ is the electric dipole operator. For the term D-SP within the diagonal spin-orbit coupling approximation, we have checked that this approach yields the same result as the computational time sparing calculation from the D-D spin-polarized contributions presented in section V (Eq. (26)).

The calculated spectra are broadened with a Lorentzian function. Furthermore, the occupied states, that do not contribute to the absorption cross section, are cut according to the method described in paragraph III-B of Ref. Brouder et al., 1996.

For the selected examples below, the different contributions to the cross sections for left- and right-circularly polarized light $\sigma({\boldsymbol{\epsilon}}_{2})$ and $\sigma({\boldsymbol{\epsilon}}_{1})$ were computed accurately in order to obtain circular dichroism.

IV Applications

IV.1 Technical details

For LiIO3, the experimental structure is used:Svensson et al. (1983) the $\Delta$ enantiomer of $\alpha$-LiIO3 belongs to the hexagonal space group $P6_{3}$ with lattice parameter $a=5.48~{}\mathrm{\SIUnitSymbolAngstrom}$ and $c=5.17~{}\mathrm{\SIUnitSymbolAngstrom}$. The atomic positions Goulon et al. (1998) are Li 2(a) (0,0,0.076), I 2(b) (1/3,2/3,0) and O 6(c) (0.247,0.342,0.838). The $\Lambda$ enantiomer is the mirror image of the $\Delta$ one (see Fig. 1) and it belongs to the same space group. A $2{\times}2{\times}2$ supercell (80 atoms) is used so that the smallest distance between a core-hole and its periodic image is $10.344~{}\mathrm{\SIUnitSymbolAngstrom}$. Gamma-centered k-points grids $3{\times}3{\times}3$ for the self-consistent charge density calculation and $9{\times}9{\times}9$ for the spectra calculation are used. A constant Lorentzian broadening, with full width at half maximum set to the core-hole lifetime broadening 3.46 eV,Fuggle and Inglesfield (1992) is applied. As XNCD is a structural effect and not a magnetic effect, the calculation is not spin-polarized.

The XMCD calculations for the 3$d$ ferromagnetic metals are carried out by using the following experimental lattice parameters: $a=2.87~{}\mathrm{\SIUnitSymbolAngstrom}$ for bcc Fe, $a=3.52~{}\mathrm{\SIUnitSymbolAngstrom}$ for fcc Ni and $a=2.51~{}\mathrm{\SIUnitSymbolAngstrom}$ and $c=4.07~{}\mathrm{\SIUnitSymbolAngstrom}$ for hcp Co. The number of atoms per supercell is 64 atoms for Fe and Ni and 96 atoms for Co, so the smallest distance between the periodically repeated core-holes is $9.84~{}\mathrm{\SIUnitSymbolAngstrom}$ in Fe, $9.97~{}\mathrm{\SIUnitSymbolAngstrom}$ in Ni and $10.03~{}\mathrm{\SIUnitSymbolAngstrom}$ in Co. A Methfessel-Paxton cold smearing of 0.14 eV (0.01 Ry) and a centered $2{\times}2{\times}2$ k-points grid are used for the self-consistent charge density calculation. The spectra calculation is performed with a $6{\times}6{\times}6$ grid for Fe and Co and a $8{\times}8{\times}8$ grid for Ni. These calculations are performed with collinear spins along the easy axis of the crystals, that is to say, [001] for bcc Fe and hcp Co and [111] for fcc NiO’Handley (1999) and the wave vector $\mathbf{k}$ is set along the same axis.

The spectra are convolved with a Lorentzian broadening function to simulate the effect of the finite lifetime of the core-hole (constant in energy) and of the inelastic scattering of the photoelectron (additional energy-dependent broadening) for which we use the curves published by Müller et al.Müller et al. (1982)

Experimental and calculated spectra are normalized such that the edge jump is equal to 1.

During the calculation of the spectra the origin of energy is set to the Fermi energy of the material $E_{F}$. For the spectra to be compared with experiment, a rigid shift in energy is applied to the calculated spectra to make the maxima of the calculated XAS correspond to the maxima of the experimental spectra. The same shift is applied to the XMCD spectra. In the plots, the origin of energy $E_{0}$ is therefore the one chosen in the publications from which the experimental spectra are extracted.

IV.2 XNCD at the L1-edge of I in $\alpha$-LiIO3

Natural circular dichroism in the inorganic non-centrosymmetric lithium iodate (LiIO3) crystal have been measured in 1998 Goulon et al. (1998) and it has been attributed to the interference of electric dipole and electric quadrupole transitions.Goulon et al. (1998); Natoli et al. (1998) Previous calculations Natoli et al. (1998); Goulon et al. (1998); Ankudinov and Rehr (2000) were indeed able to reproduce the overall peak positions and intensities in this framework. The agreement is however not entirely satisfactory for the absorption spectra. These discrepancies have been attributed to the use of muffin-tin potentials.Ankudinov and Rehr (2000)

The approach presented in this work, that does not rely on the muffin-tin approximation, was applied to compute the XAS and XNCD spectra for $\alpha$-LiIO3. The absorption is dominated by the electric dipole-electric dipole term (D-D) as shown in Fig.2 b). The XNCD spectra, on the other hand, is entirely due to the cross term electric dipole-electric quadrupole (D-Q).

As illustrated by Fig.2 a), both the calculated XAS and XNCD spectra at the I L1-edge are in good agreement with experiment. However, the amplitude of the calculated XNCD is $4{\times}10^{-2}$ compared to the edge jump while the amplitude of the experimental spectra from Ref. Goulon et al., 1998 is $6.5{\times}10^{-2}$. Such an underestimation was also observed in Ref. Natoli et al., 1998 within a multiple-scattering approach.

From Fig.2 a) bottom, it becomes obvious that the XNCD spectra for both enantiomers are opposite. Indeed, it has the same effect for XNCD to change an enantiomer for the other ($\Delta\leftrightarrow\Lambda$) as for XMCD to change the sign of the magnetic field ($\mathbf{B}\leftrightarrow-\mathbf{B}$).

The angular dependence of the calculated XNCD spectra is depicted in Fig. 2 c) and its amplitude is plotted in inset as a function of $\theta$, the angle between $\mathbf{k}$ and the $\mathbf{c}$-axis of the crystal. This amplitude varies as $3\cos^{2}\theta-1$ so it is maximal in the case $\mathbf{k}$ is parallel to the $\mathbf{c}$ axis. This dependence is consistent with the formula derived in Ref. Natoli et al., 1998 for point group $C_{6}$ (point group of the space group of the crystal). Note that, as ${\boldsymbol{\epsilon}}$ is kept perpendicular to $\mathbf{k}$ and $C_{6}$ is a dichroic point group,Brouder (1990) the XAS spectra also present an angular dependence. It does not prevent a comparison of the amplitude of the XNCD spectra because the edge jump remains unchanged.

IV.3 XMCD at the K-edge of 3d transition metals

XMCD was recorded for the first time at the Fe K-edge in magnetized Fe in 1987. Schütz et al. (1987) Ever since, a large number of calculations for the electric dipole term of the XMCD spectra on Fe K-edge in bcc Fe in the XANES region have been reported, for example in Ref. Gotsis and Strange, 1994; Igarashi and Hirai, 1994; Brouder et al., 1996; Guo, 1996; Ebert, 1996; Fujikawa and Nagamatsu, 2003; Šipr and Ebert, 2005; Dixit and Alouani, 2016. Calculations of XMCD at the K-edge in fcc Ni and hcp Co are fewer Igarashi and Hirai (1994, 1996); Guo (1996); Torchio et al. (2011b) and are not really conclusive.

These calculations have been performed with various methods, often within the electric dipole approximation and with muffin-tin potentials. Here, we present the calculation of the three terms (D-D, Q-Q and D-SP) that are likely to contribute to the XMCD cross-section at the K-edge of ferromagnetic 3d transition metals showing the relevance of the D-SP term.

The contribution of the D-SP term to the absorption cross-section is not shown here because it is negligible. On the other hand its contribution to the XMCD spectra (Fig. 3) is significant: it reaches 28% of the D-D term in amplitude. This can be understood considering the sum-rules that are made explicit in the next section: in the XMCD cross-section, the D-SP term probes the spin polarization of the $p$ states whereas the D-D term probes their orbital polarization. In Ref. Igarashi and Hirai, 1996 the $4p$ orbital magnetic moment in Co, Fe and Ni is evaluated to a few $10^{-4}\mu_{B}$ (Fe: 5$\times 10^{-4}\mu_{B}$, Co: 16$\times 10^{-4}\mu_{B}$, Ni: 6$\times 10^{-4}\mu_{B}$) and in Ref. Chen et al., 1995 the $4p$ spin magnetic moment in Fe and Co is evaluated to several $10^{-2}\mu_{B}$ (Fe: 5$\times 10^{-2}\mu_{B}$, Co:6$\times 10^{-2}\mu_{B}$) in the opposite direction. This difference in order of magnitude of both quantities compensate for the smallness of prefactor ($\hbar\omega/4mc^{2}$) of the D-SP term (see Table 1 in appendix).

To check possible numerical problems, we also performed the calculations using the FDMNES codeBunǎu and Joly (2009) where, for this purpose, the D-SP term was introduced in the same way. This code follows Wood and Boring Wood and Boring (1978) to eliminate the small component and obtain a couple of Schrödinger-like equations, including the spin-orbit effect, closely akin but improving the Pauli equation. Despite the very different approach (no pseudo-potential, calculation in real space and no diagonal spin-orbit coupling approximation), we found very similar results for both the shape and relative amplitude of the D-SP contribution.

The agreement with the experimental spectra is fair as illustrated Fig. 5. As usual in independent-particle calculations, the energy axis is slightly compressed Materlik et al. (1983); Grunes (1983); Lengeler and Zeller (1984) due to the energy dependence of the real part of the self-energyMustre de Leon et al. (1991) for which corrections to the calculated spectra could be applied.Kas et al. (2009) Alternatively, the position of the calculated peaks could be improved by phenomenological rescaling.Mustre de Leon et al. (1991); Materlik et al. (1983)

For Fe, the main peaks of the experimental XMCD are reproduced by the calculation. As in calculations by others,Šipr and Ebert (2005) the positive peak at 10 eV is overestimated probably due to the approximate description of the exchange-correlation energy. Indeed, the comparison between the spectra calculated with PBE or LDA functionals (Fig. 4) shows that this peak would be even more enhanced with LDA.

For Ni and Co, a main negative peak is present near the main rising edge in the calculation as in the experiment but the satellite peaks that appear in the calculation are difficult to link to the experiment.

In these calculations, the polarization rate of the light is taken to be 100 % and a single crystal with full $3d$ spin polarization is considered. In Fe, Ni and Co, saturation is reached with usual experimentally applied magnetic field and the anisotropy is quite weak so that the rate of circular polarization of the light $P_{c}$ is expected to account for most of the discrepancy in amplitude between the calculated and the experimental XMCD spectra. The data for Fe and Co were recorded in a 5 T magnet at 5 K and within a setup that reaches 90 % circular polarization rate.Facility (accessed November 30, 2016) The correction on the amplitude of the calculated spectra to fit the experimental condition should therefore be of order 0.9. Here, it is approximately 0.6 in the case of Fe and 1.0 in the case of Co. The data for Ni were recorded at ambient temperature in a 0.7 T magnet within a dispersive setup with a diamond quarter-wave plate for which we can infer that $P_{c}\approx 0.7$.Giles et al. (1994) However, no correction on the amplitude of the calculated spectra is needed to make it correspond to the amplitude of the experimental spectra. So, whereas our calculation overestimates the amplitude of the XMCD spectra in the case of Fe, it underestimates it in the case of Ni.

V Contribution of the D-SP term to XMCD: the case of collinear spins

V.1 The SP operator

In this section, we study the spin-position operator $SP({\boldsymbol{\epsilon}})={\boldsymbol{\sigma}}\cdot({\boldsymbol{\epsilon}}\times\mathbf{r})$. We consider collinear spins along $z$ with independent spin channels. The spin part of the wave functions $|s\rangle$ can either be the spin up spinor $\begin{pmatrix}1\\ 0\end{pmatrix}$, or the spin down spinor $\begin{pmatrix}0\\ 1\end{pmatrix}$.

The D-SP term is the cross term between the electric dipole and the spin-position operator. Spin does not appear in the electric dipole operator, so it is diagonal in spin:

[TABLE]

This imposes $s^{\prime}=s$. On the other hand, the vector of Pauli matrices ${\boldsymbol{\sigma}}$ appears explicitly in the Spin-Position operator:

[TABLE]

As $\langle s|\sigma_{x}|s\rangle=\langle s|\sigma_{y}|s\rangle=0$, we can exclude a priori the terms that are proportional to $\sigma_{x}$ and $\sigma_{y}$ in the Spin-Position operator. In that case the spin position operator rewrites:

[TABLE]

Its selection rules are almost the same as for the electric dipole Sébilleau (2006) one: $\Delta l=\pm 1$, $\Delta m=\pm 1$.

As $Y_{1}^{-1}({\boldsymbol{\epsilon}}_{1})=0$, $Y_{1}^{-1}({\boldsymbol{\epsilon}}_{2})=\sqrt{3/4\pi}$, $Y_{1}^{0}({\boldsymbol{\epsilon}}_{1})=Y_{1}^{0}({\boldsymbol{\epsilon}}_{2})=0$, $Y_{1}^{1}({\boldsymbol{\epsilon}}_{1})=-\sqrt{3/4\pi}$ and $Y_{1}^{1}({\boldsymbol{\epsilon}}_{2})=0$,

[TABLE]

Hence,

[TABLE]

with

[TABLE]

where $s=\uparrow$ or $\downarrow$. Therefore, in the diagonal spin-orbit coupling collinear spins case, the D-SP term can be computed from the D-D cross section for the up and down spin channels.

V.2 Sum-rule at K-edge

A sum-rule is a formula in which the integral of the circular dichroism spectra due to a given term of the cross-section is expressed as a function of ground state expectation value of some operator. The sum-rules at L2,3-edges are well established Wu et al. (1993); Chen et al. (1995) and are widely used to extract quantitative magnetic ground state properties. Their derivation is based on several approximations among which the fact that the radial integrals are spin and energy independent.Altarelli (1998) At K-edge the sum-rule for the electric dipole-electric dipole term Ankudinov and Rehr (1995); Carra et al. (1993); Altarelli (1993); Igarashi and Hirai (1996) relates the integral of the XMCD spectra to the orbital magnetic moment of occupied $p$ states that is proportional to ${<}L_{z}{>}_{p}$. This sum rule is however almost impossible to apply in practice because the upper limit of the integral is not well defined and, in the case of 3$d$ transition elements, the 4$p$ states are almost unoccupied so ${<}L_{z}{>}_{p}$ is very small and has a minor impact on the magnetic moment of the material. Deriving a similar sum-rule for the D-SP term is nevertheless very useful to understand why, despite its very small prefactor, this term is so large in XMCD. We derive it following the method of Thole et al. Thole et al. (1992); Altarelli (1993) with many body wave functions and operators assuming all spins collinear and within the diagonal spin-orbit coupling approximation.

In a many body framework, using the expression for $SP$ in terms of spherical harmonics Eq (22):

[TABLE]

with

[TABLE]

In a second quantized form with $l$,$m$ and $\sigma$ the usual quantum numbers:Natoli et al. (1998)

[TABLE]

where ${\cal{D}}_{l_{0},l}=\int\mathrm{d}r\,r^{3}R^{\star}_{l_{0}}(r)R_{l}(r)$ is assumed - as usual in sum-rules derivations - to be spin-independent. The experimental procedure enables to obtain the signal corresponding to a specific $l_{0}$. At K-edge $l_{0}=0$ and $m_{0}=0$ so that,

[TABLE]

where ${\cal{D}}={\cal{D}}_{0,1}$.

Similarly, as $\langle\sigma_{0}^{\prime}|\sigma_{z}|\sigma^{\prime}\rangle=\sigma^{\prime}\delta_{\sigma_{0}^{\prime},\sigma^{\prime}}$

[TABLE]

Using the completeness relation $\int\mathrm{d}E\,\sum_{f}|f\rangle\langle f|\delta(E_{f}-E_{g}-E)=1\!\!1-|g\rangle\langle g|$, as the core shell is full and under the assumption that the radial integral ${\cal{D}}$ does not depend on energy:

[TABLE]

The combination of Eq. (33) and Eq. (28) leads to:

[TABLE]

The difference between the two integrals yields the XMCD sum rule for the D-SP term:

[TABLE]

with the operator

[TABLE]

corresponding to a partial spin polarization of the occupied $p$ states.

If one considers the derivative of this sum-rule, we see that the electric dipole - spin position (D-SP) circular dichroism signal probes the spin polarization of the empty $p$ states. Fig. 6 illustrates the correspondence between both quantities. This proves the validity of the D-SP sum-rule. Unfortunately, this sum rule can not be applied directly on experimental spectra, mainly because of the superposition of the D-D contribution to the D-SP contribution.

VI Conclusion

We have developed an efficient computational approach to determine accurate XMCD and XNCD spectra.

The main result is that the contribution from the relativistic term D-SP in the transition operator is significant in XMCD spectra despite being negligible in XAS. This importance is explained by the fact that this term probes the spin of the $p$ states that is two orders of magnitude larger than its orbital counterpart.

For XNCD, the calculated spectra are in good agreement with experiment and the angular dependence corresponds to the expected one.

A big advantage of the method employed in this paper to perform XMCD and XNCD calculations is its wide adaptability that opens opportunities for applications to several kinds of systems such as strongly-correlated materials or molecules absorbed on functionalized surfaces. The same method could be apply to compute X-ray magneto-chiral dichroism (XM$\chi$D) that has been observed in magnetized chiral systems Sessoli et al. (2015). The features of XM$\chi$D differ from the one of XMCD and XNCD making it a promising probe of the interplay between chirality and magnetism.

Acknowledgement

This work was supported by French state funds managed by the ANR within the Investissements d’Avenir programme under Reference No. ANR-11-IDEX-0004-02, and more specifically within the framework of the Cluster of Excellence MATISSE led by Sorbonne Universités. We are grateful to Delphine Cabaret for very interesting and constructive feedback on this work. We also thank François Baudelet and Lucie Nataf for providing reference spectra for the energy scaling. UG and NJV acknowledge support by DFG (FOR 1405). The numerical calculations have been performed using HPC resources from the Paderborn Center for Parallel Computing (PC2) and from GENCI-IDRIS (Grant i2016096863).

Appendix A Semi-relativistic transformation of the relativistic cross section

We start from the expression for the cross section in a relativistic framework Bouldi and Brouder (2016) and we adapt it to the specific need of our numerical method that is the determination of large components of the Dirac wave function for the core state and of Foldy-Wouthuysen (FW) wave functions for the valence states.

A.1 Relativistic cross section

The contribution to the X-ray absorption (XAS) cross section from a given four-components Dirac core-state $|\Psi_{i}\rangle$ of energy $E_{i}$ is given by:Bouldi and Brouder (2016)

[TABLE]

where the sum runs over unoccupied final states $|\Psi_{f}\rangle$ with energy $E_{f}$, $\alpha_{0}$ is the fine structure constant and $T_{D}$ is the transition operator defined as:

[TABLE]

where the polarization vector ${\boldsymbol{\epsilon}}$, the wave vector $\mathbf{k}$ and the energy $\hbar\omega$ describe the incident electromagnetic wave, $\mathbf{r}$ is the position operator and $\boldsymbol{\alpha}=(\alpha_{x},\alpha_{y},\alpha_{z})$ is the vector of Dirac matrices.

Here, as in our numerical calculations, a one-electron scheme is used. In a many-body framework the formula for the cross-section would be the same but with N-electrons wave functions and many-body operators that write as sums over electrons.

In Ref. Bouldi and Brouder, 2016, the transformation into a two-component representation for $|\Psi_{i}\rangle$ and $|\Psi_{f}\rangle$ was performed by applying a time-independent Foldy-Wouthuysen transformation (FW) at order $c^{-2}$. The FW transformation of $\Psi_{l}$ is obtained by applying a unitary operator: $\psi_{l}^{\mathrm{FW}}=U_{\mathrm{FW}}\Psi_{l}$ with,Eriksen (1958)

[TABLE]

where $\beta$ the standard Dirac matrix. In this expression, ${\cal{O}}$ is the odd operator entering the Dirac Hamiltonian: ${\cal{H}}^{D}=\beta mc^{2}+{\cal{O}}+{\cal{E}}$ where ${\cal{E}}$ is even. It is defined as ${\cal{O}}=c\boldsymbol{\alpha}\cdot(\mathbf{p}-e\mathbf{A}_{0})$ where $\mathbf{p}$ is the momentum operator and $\mathbf{A}_{0}$ is the static vector potential.

Only the large components of $\psi_{l}^{\mathrm{FW}}$, denoted $\phi_{l}^{\mathrm{FW}}$, are non zero up to order $c^{-2}$. The cross-section can be written as a function of the large components of $\psi_{i}^{\mathrm{FW}}$ and $\psi_{f}^{\mathrm{FW}}$:Bouldi and Brouder (2016)

[TABLE]

The operator $T_{\mathrm{FW}}$ is the projection on upper components of $U_{\mathrm{FW}}T_{D}U_{\mathrm{FW}}^{\dagger}$:

[TABLE]

where

[TABLE]

and

[TABLE]

are the standard electric dipole and electric quadrupole operators.

The magnetic dipole operator $T_{\mathrm{MD}}$ writes:

[TABLE]

where $\mathbf{L}=\mathbf{r}\times\mathbf{p}$ and ${\boldsymbol{\sigma}}$ is the vector of Pauli matrices. It is proportional to the total magnetic moment operator $(\hbar{\boldsymbol{\sigma}}+\mathbf{L})=(2\mathbf{S}+\mathbf{L})$ where $\mathbf{S}$ is the spin operator. $T_{\mathrm{MD}}$ is also present in common non-relativistic derivations. Brouder (1990); Di Matteo et al. (2005) Its selection rules are $l_{i}=l_{f}$ and $n_{i}=n_{f}$,Brouder (1990) so it vanishes in the X-ray energy range because the states involved in the transition have different principal quantum numbers.

The correction to this term due to the static vector potential $\mathbf{A}_{0}$ is:

[TABLE]

The last term in Eq. (41) is only present when relativistic effects are included in the calculation of the transition operator:

[TABLE]

A similar term was already found in Ref. Gougoussis, 2009, but derived from a semi-relativistic Hamiltonian and this approach presents a conflict with time-dependent perturbation theory.Bouldi and Brouder (2016) It can be rewritten noticing that, in the non-relativistic limit, $|\phi_{i}^{\mathrm{FW}}\rangle$ and $|\phi_{f}^{\mathrm{FW}}\rangle$ are eigenstates of:

[TABLE]

where $\mathbf{B}_{0}$ is the static external magnetic field. This Hamiltonian obeys $\mathbf{p}-e\mathbf{A}_{0}=(m/i\hbar)[\mathbf{r},H_{0}^{0}]$ so that,

[TABLE]

We name ${\boldsymbol{\sigma}}\cdot({\boldsymbol{\epsilon}}\times\mathbf{r})$, the spin-position operator and define the associated transition operator:

[TABLE]

For technical reasons, in the present paper we consider a different situation than in Ref. Bouldi and Brouder, 2016: we use FW wave function for the final states and large components of the Dirac wave function for the initial (core) state. This difference in treatment is linked to the fact that the core wave function is determined from a relativistic atomic code whereas the unocupied states are calculated with a semi-relativistic condensed-matter code.

A.2 Rewriting the cross section with large components of the Dirac wave function for the core state

We note $\phi_{i}$ and $\chi_{i}$ the large and small components of $\Psi_{i}$. The order of magnitude of the ratio between small and large components is $v/c$ where $v$ is the velocity of the particle.Strange (1998) Up to order $c^{-1}$, the small component writes: Strange (1998); van Lenthe et al. (1996)

[TABLE]

Only the second term in $U_{\mathrm{FW}}$ Eq. (39) couples the small and the large components. From Eq. (49) and (39), the large component of the FW transformed wave function can be expressed as a function of the large components of the Dirac wave function up to order $c^{-2}$,

[TABLE]

$[{\cal{O}}^{2}]_{p}$ is the projection of ${\cal{O}}^{2}$ on large components: Strange (1998)

[TABLE]

and ${\cal{O}}_{p}=c{\boldsymbol{\sigma}}.(\mathbf{p}-e\mathbf{A}_{0})$ is the projection of $\beta{\cal{O}}$ on the upper right components.

The identity $c{\cal{O}}_{p}{\boldsymbol{\sigma}}\cdot(\mathbf{p}-e\mathbf{A}_{0})=[{\cal{O}}^{2}]_{p}$ leads to:

[TABLE]

From this relation, the cross section of Eq. (40) can be adapted to the case that we consider here.

In Ref. Bouldi and Brouder, 2016 the expansion was made to order $1/c^{2}$ for the dipole contribution and to order $kr$ for multipole contributions. At the same order,

[TABLE]

where $T_{\mathrm{FW}}^{\prime}$ is:

[TABLE]

There is one extra terms in the cross section compared to $T_{\mathrm{FW}}$ that is related to the use of large components of the Dirac wave function instead of Foldy-Wouthuysen wave function for the core state:

[TABLE]

We show in the next subsection that it is negligible for the core states considered in this work.

As the magnetic dipole term is negligible in the X-ray range, $T_{\mathrm{FW}}$ thus contains four operators (see Eq. (41) and the subsequent comments) so $T_{\mathrm{FW}}^{\prime}$ writes,

[TABLE]

A.3 Order of magnitude of the operators

As the core wave function is very localized, we obtain an idea of the relative order of magnitude of the operators in Eq. (54) by evaluating them at the radius corresponding to the core state. In Table 1 these evaluations are given compared to the electric dipole operator.

When expanding the square modulus of the matrix elements in Eq. (52), we keep the terms with contributions higher than $10^{-3}$ compared to the dominant electric dipole term. We also neglect the term $T^{e}$: as $V(\mathbf{r})$ is almost spherical at the core state radius, it concerns transitions to the same orbitals as the electric dipole term. It does not contain a spin operator so that, even in XMCD, it only yields a negligible correction to the electric dipole contribution. Therefore, we are left with the four significant terms D-D, Q-Q, D-Q and D-SP discussed in section II.