# Theory of X-ray absorption spectroscopy: a microscopic Bloch equation
approach for two-dimensional solid states

Dominik Christiansen1 Malte Selig1 Jens Biegert2,3,4 Andreas Knorr1 1Institut
für Theoretische Physik, Nichtlineare Optik und Quantenelektronik, Technische
Universität Berlin, 10623 Berlin, Germany  2ICFO-Institut de Ciencies
Fotoniques, The Barcelona Institute of Science and Technology, 08860
Castelldefels (Barcelona), Spain  3ICREA-Institució Catalana de Recerca i
Estudis Avançats, Barcelona, Spain  4Fritz-Haber Institute of the Max Planck
Society, Berlin, Germany

###### Abstract

We develop a self-consistent Maxwell-Bloch formalism for the interaction of
X-rays with two-dimensional crystalline materials by incorporating the Bloch
theorem and Coulomb many-body interaction. This formalism is illustrated for
graphene, by calculating the polarization-dependent XANES, formulating
expressions for the radiative and Meinter-Auger recombination of core-holes,
and the discussion of microscopic insights into the spectral oscillations of
EXAFS beyond point scattering theory. In particular, the correct inclusion of
lattice periodicity in our evaluation allows us to assign so far uninterpreted
spectral features in the Fourier transformed EXAFS spectrum.

##  I Introduction

At the beginning of the last century X-ray absorption spectroscopy (XAS) has
been the key instrument to acquire knowledge about atomic energy levels
Bohr1913 ; Sommerfeld1916 ; Moseley1913 ; Lundquist1925 ; Coster1927  and
later for the discovery and systematization of rare-earth elements Roehler1987
. The constant development of intense synchrotron radiation sources during
these times led to a steady growth in X-ray experiments Lytle1965 ;
Eisenberger1978 ; Bordwehr1989 ; Sauer2008 ; Young2018 ; Gu2021 . At the end
of the century a connection between the X-ray absorption spectrum and the
local structure of the investigated material has been drawn. This promoted XAS
from a spectroscopic to a structural characterization technique Sayers1970 ;
Sayers1971 : X-ray absorption spectroscopy, probing the core electronic states
of atomic systems, is a highly element-sensitive and environment-specific
spectroscopic technique with applications in atom and molecular physics,
chemistry, biology and material science Eisenberger1978_2 ; Yano2009 ;
Kraus2018 ; Diller2017 ; Hua2019 . For this, a large range of tunable X-ray
synchrotron sources Esarey1993 ; Schoenlein2000 ; Buades2018  sensitive for
most elements in the periodic table are available.

In this contribution, we develop, based on a Bloch function approach in the
Heisenberg equation formalism, a microscopic description of X-ray induced
electronic transitions and ionization for solid state two-dimensional
semiconductors, Dirac semimetals or insulators Novoselov2005 ; Balleste2011 ;
Xu2013 ; Gu2018 . We include nonlinear excitation effects, discuss X-ray
absorption dynamics in solids beyond Fermi’s golden rule, and include
explicitly the lattice periodicity of crystalline solids. In particular, we
shine light on the microscopic origin of the spectral oscillations observed in
absorption and identify quantum-interference peaks, which can not be described
in point scattering theory Buades2018 .

Figure 1: (a) Lattice periodic potential resulting from periodically arranged
atoms in a solid. Sketched are exemplary flat core states, continuous bands
forming higher lying valence bands (VB) and conduction bands (CB) and the
ionization continuum (IC) above the ionization threshold. Transitions into the
conduction bands and the ionization continuum form the XANES and EXAFS,
respectively. (b) Sketch of a typical XAS spectrum with its (i) XANES and (ii)
EXAFS part.

Typical experimental X-ray absorption spectra exhibit two general features:
(i) resonance lines in the absorption coefficient at specific energies below
the ionization threshold and (ii) a continuous absorption coefficient above
the ionization threshold over a large photon energy range (Lindsay1928, ;
Coster1929, ; Kievit1930, ; Norman1986, ; Rehr2000, ; Buades2018, ), including
spectral oscillations. The latter appear only for molecules or crystals but
not for single atoms. Figure 1(a) shows the possible excitations, while Fig.
1(b) sketches the corresponding spectral features.

The first feature (i) is related to the transition energy of most-inner core
shell electrons to unoccupied states in the conduction band (CB) or partly
filled valence bands (VB). This characteristic resonant feature is usually
referred to as absorption edge and the spectroscopic technique exploiting such
transitions is called X-ray absorption near-edge structure (XANES) or near-
edge X-ray absorption fine structure (NEXAFS). In the existing literature, to
describe XANES, typically Fermi’s golden rule for the transition probability
from core bands to unoccupied bands under illumination with a X-ray frequency
$\omega$ is used Stohr1992 ; Chowdhury2012 ; Schnohr2015 :

$\displaystyle\alpha(\omega)=\frac{\omega\pi}{\epsilon_{0}cn}\int\int
d\mathbf{k}~{}d\mathbf{k}^{\prime}~{}|\mathbf{e}\cdot\mathbf{d}_{\mathbf{k},\mathbf{k}^{\prime}}|^{2}\delta(E_{\mathbf{k}^{\prime}}-E_{\mathbf{k}}-\hbar\omega)\;.$
(1)

Here, the transition probability $\alpha(\omega)$ involving initial core
states $\mathbf{k}$ to final conduction band states $\mathbf{k}^{\prime}$ is
determined by the product of the transition dipole moment
$\mathbf{d}_{\mathbf{k},\mathbf{k}^{\prime}}$ projected on the incident X-ray
polarization $\mathbf{e}$ and an energy conserving delta-function during the
transition. The prefactor incorporates the dielectric constant $\epsilon_{0}$,
speed of light in vacuum $c$, and the refractive index $n$ of the material.

The part of the absorption spectrum containing spectral oscillations above the
ionization energy (ii) is named extended X-ray absorption fine structure
(EXAFS) spectrum. It results from the transition of core electrons to the
ionization continuum (IC) above the ionization threshold of the material, cf.
Fig. 1(a). The explanation of the EXAFS spectrum is based on the theory
introduced by Kronig for molecular gasesKronig1931 ; Kronig1932 . This
description explains the oscillations in the spectrum as interference effects
emerging from secondary photoelectron waves emitted by an X-ray absorbing atom
and waves back scattered by neighboring atoms. The corresponding
parametrization of the absorption by Sayers et al. Sayers1970 ; Sayers1971

$\displaystyle\alpha_{k}$
$\displaystyle=S_{0}^{2}\sum_{i}N_{i}\frac{|f_{i}(k)|}{kR_{i}^{2}}\sin\left(2kR_{i}+\phi_{i}\right)e^{-\frac{2R_{i}}{\lambda_{k}}}e^{-2\sigma_{i}^{2}k^{2}}$
(2)

has become standard in describing EXAFS. Equation (2) describes the absorption
coefficient as function of the electronic wave number $k$ of the ionized, free
electron. Here, for X-ray energies higher than the ionization energy $E_{Ion}$
of the material (EXAFS), the excess energy is transfered to the excited
photoelectron in form of kinetic energy. The wave number of the emitted
electrons is related to the X-ray excitation frequency $\omega$ by
$k^{2}=2m_{0}(\hbar\omega-E_{Ion})/\hbar^{2}$ with the free (vacuum) electron
mass $m_{0}$. Here, $\lambda_{k}$ stands for the X-ray wavelength. The
structural parameters are the number of equivalent scatterers $N_{i}$ of type
$i$, the interatomic distances $R_{i}$, and the bond length $\sigma_{i}$. The
Debye-Waller factor accounts for thermally or disorder induced changes in the
bond length. The quantity $|f_{i}(k)|$ describes the back scattering
amplitude, resulting from back scattering of neighboring atoms, and $S_{0}$
describes the reduction factor due to multielectron processes. The exponential
factor $\exp(-2R_{i}/\lambda_{k})$ accounts for the finite lifetime of the
photoelectron, which is only scattered elastically over a short distance.
Finally, $\phi_{i}$ stand for phase shifts due to the initially excited and
the back scattering atom Sayers1971 ; Rehr2000 . The wave functions used to
obtain this formula are not lattice periodic and therefore do not exploit the
Bloch theorem in crystalline matter. While this approach might be appropriate
for molecules, where only a few $R_{i}$ have to be accounted for, a X-ray
absorption theory for solids should take the lattice periodicity into account.

In this manuscript, to address recent advances in pulsed X-ray
spectroscopyPopmintchev2012 ; Teichmann2016 ; Pertot2017  of solids we develop
a dynamical approach to account for the description of ultrafast many-body
solid state phenomena in time resolved X-ray spectroscopy Angela2015 ;
Picon2019 ; Chang2021 : We develop a microscopic, temporal resolved
formulation of the underlying mechanisms of the X-ray absorption in a two-
dimensional solid. For this, we use the method of second quantization to also
have a basis to sequentially include many-body interaction at a later stage of
theory. In particular, we shine light on the origin of peaks in the Fourier
transformed EXAFS spectrum, which are not captured by the point scattering
theory but observable in experiment Buades2018 . Incorporating Maxwell
equations for the X-ray field a self-consistent coupling between the light
field and the microscopic core-transition including radiative damping and
lineshifts can be achieved. Including also many-body interactions the
resulting Maxwell-Bloch formalism provides a powerful tool to investigate core
excitations induced by classical light in solids.

The manuscript is organized as follows: We start in Section II by setting up
the Hamiltonian on the basis of single-particle electronic wave functions and
their interaction matrix elements for the XANES and EXAFS processes. This
section is separated into two parts: the first is using a general Bloch wave
approach suited for periodic atomically-thin solids, while the second part
bridges the Bloch description with a tight binding approximation often used in
electronic structure theory. In Section III we derive the dynamical Bloch
equations for the X-ray induced electronic transitions. Here, we do not only
consider the X-ray induced core excitations but also display the Coulomb-
induced mean-field Hartree-Fock and relaxation channels of the excited core
electrons. In particular, the equations of motion include the spatial
resolution of the X-ray radiation. In Section IV, we solve the wave equation
to couple self-consistently the Maxwell equations to the non-local microscopic
excitations, including nonlinear effects and finally apply in Section V the
developed theory to the exemplary material of graphene.

##  II Solid state electron–X-ray interaction

###  II.1 Hamiltonian

To develop a microscopic theory of electron-X-ray interaction in solids, we
start by deriving the interaction Hamiltonian. A schematic picture of the
electronic structure is sketched in Fig. 1(a). Typically, electrons in
crystalline solids can be separated into two groups: core and valence
electrons. Core electrons occupy filled orbitals and are spatially localized
at the nuclei, both contributing to core-ions, cf. Fig. 1(a). In contrast,
valence electrons are less localized. In addition to the bound states (core,
valence and conduction bands), if ionization processes take place, the
ionization continuum (IC, cf. Fig. 1(a)) for electronic states above the
ionization edge must be taken into account. The corresponding field operators
for electrons can be expanded into a complete set $\\{\alpha\\}$ of (i) core,
(ii) valence, conduction and (iii) ionization continuum states:
$\hat{\Psi}^{(\dagger)}(\mathbf{r},t)=\sum_{\alpha}\Psi^{(*)}_{\alpha}(\mathbf{r})a^{(\dagger)}_{\alpha}(t)$
with quantum number $\alpha$, single particle wave functions
$\Psi_{\alpha}(\mathbf{r})$ and the fermionic annihilation (creation)
operators $a^{(\dagger)}_{\alpha}$. The Hamiltonian contains the single
electron energies $E_{\alpha}$ in the lattice periodic potential, the carrier-
carrier interaction with Coulomb matrix element
$V^{\alpha\beta}_{\alpha^{\prime}\beta^{\prime}}$ and the light-matter
interaction $\Omega_{\alpha\alpha^{\prime}}$ in second quantization:

$\displaystyle H$
$\displaystyle=\sum_{\alpha}E_{\alpha}a^{\dagger}_{\alpha}a^{\mathstrut}_{\alpha}-\hbar\sum_{\alpha,\alpha^{\prime}}\Omega_{\alpha\alpha^{\prime}}~{}a^{\dagger}_{\alpha}a^{\mathstrut}_{\alpha^{\prime}}$
$\displaystyle+\frac{1}{2}\sum_{\alpha,\alpha^{\prime},\beta,\beta^{\prime}}V^{\alpha\beta}_{\alpha^{\prime}\beta^{\prime}}~{}a^{\dagger}_{\alpha}a^{\dagger}_{\beta}a^{\mathstrut}_{\beta^{\prime}}a^{\mathstrut}_{\alpha^{\prime}}\;.$
(3)

The light-matter interaction, described in length gauge discussed in the
supplemental material (SM Sec. I), is given by the Rabi frequency
$\Omega_{\alpha\alpha^{\prime}}$, and reads
$\Omega_{\alpha\alpha^{\prime}}=e_{0}\langle\Psi_{\alpha}|\mathbf{r}\cdot\mathbf{E}(\mathbf{r},t)|\Psi_{\alpha^{\prime}}\rangle/\hbar$,
where $e_{0}$ is the elementary charge. The Coulomb matrix element reads
$V^{\alpha\beta}_{\alpha^{\prime}\beta^{\prime}}=\langle\Psi_{\alpha}\Psi_{\beta}|V(\mathbf{r-r^{\prime}})|\Psi_{\beta^{\prime}}\Psi_{\alpha^{\prime}}\rangle$
with the Coulomb potential $V(\mathbf{r-r^{\prime}})$ in real space discussed
in the supplemental material. The single-particle eigenbasis
$\Psi_{\alpha}(\mathbf{r})$ of the one-particle Hamiltonian and the
corresponding eigenenergies $E_{\alpha}$ are obtained from the Schrödinger
equation

$\displaystyle
H_{0}\Psi_{\alpha}(\mathbf{r})=\left(-\frac{\hbar^{2}\nabla^{2}}{2m}+U(\mathbf{r})\right)\Psi_{\alpha}(\mathbf{r})=E_{\alpha}\Psi_{\alpha}(\mathbf{r})$
(4)

describing electrons in a lattice periodic potential $U(\mathbf{r})$ of the
ions: $U(\mathbf{r}-\mathbf{R}_{\parallel})=U(\mathbf{r})$ with a lattice
vector $\mathbf{R}_{\parallel}$ in the in-plane motion of the atomically thin
two-dimensional structure. The vector $\mathbf{r}$ of the electron is still
three-dimensional to have access to the orbitals extending perpendicular to
the plane. Eigenfunctions of the lattice periodic Hamiltonian have to fulfill
the Bloch condition
$\Psi(\mathbf{r+R_{\parallel}})=\exp(i\mathbf{k_{\parallel}\cdot
R_{\parallel}})\Psi(\mathbf{r})$. Therefore, the eigenvectors $\Psi_{\alpha}$
and eigenenergies $E_{\alpha}$, are classified by two quantum numbers
$\\{\alpha\\}=\\{\mathbf{k}_{\parallel},\lambda\\}$, namely wave vector
$\mathbf{k}_{\parallel}$ in the plane of the two-dimensional material and band
index $\lambda$. The normalized eigenfunctions in the lattice periodic
potential are Bloch waves Bloch1929 ; Wang2018 ; SchaferBuch

$\displaystyle\Psi_{\lambda,\mathbf{k}_{\parallel}}(\mathbf{r})=\frac{1}{\sqrt{Al_{z}}}e^{i\mathbf{k}_{\parallel}\cdot\mathbf{r}_{\parallel}}u_{\lambda,\mathbf{k}_{\parallel}}(\mathbf{r})$
(5)

with lattice periodic Bloch factor
$u_{\lambda,\mathbf{k}_{\parallel}}(\mathbf{r})$, the sample area $A$ and
length of the quantization volume in perpendicular direction to the material
$l_{z}$ (sample volume $Al_{z}$). The lattice periodic function
$u_{\lambda,\mathbf{k}_{\parallel}}(\mathbf{r})$ are orthonormalized on the
in-plane unit cell: $\langle
u_{\lambda,\mathbf{k}_{\parallel}}|u_{\lambda^{\prime},\mathbf{k}_{\parallel}}\rangle_{UC}=V_{UC}\delta_{\lambda,\lambda^{\prime}}$,
with unit cell volume $V_{UC}$. Figure 2 represents a conceptional sketch of a
two-dimensional crystal built by consecutive identical unit cells. Equation
(5) represents solutions for all states $\alpha$ – (i) core, (ii)
valence/conduction band and (iii) ionization continuum state. This accounts
for the translation invariance of the crystal, a property which is also
existing for vanishing coupling among the electrons in the core states of
different atoms. In the following we discuss the properties of the different
wave functions and the approximations made in the following separately:

Figure 2: Sketch of a two-dimensional crystal built by a series of unit cells,
highlighted in green, with an area $A$ and thickness $l_{z}$. The space vector
can be decomposed into a lattice vector $\mathbf{R}_{n\parallel}$ pointing
towards the n-th unit cell and a unit cell local vector $\mathbf{r}_{n}$. The
complete crystal can then be expressed by a sum over the number of repeating
unit cells.

(i) As solution of the Schrödinger equation, Eq. (5) can be used to describe
the core state wave functions. However, core states are strongly localized at
the atomic site and it is useful to built a localized wave packet as
superposition of Bloch waves of different $\mathbf{k}_{\parallel}$ to use
Wannier functions Wannier1937 ; Marzari2012

$\displaystyle
w_{\lambda}(\mathbf{r},\mathbf{R}_{n\parallel})=\sum_{\mathbf{k}_{\parallel}}\exp(-i\mathbf{k}_{\parallel}\cdot\mathbf{R}_{n\parallel})\Psi_{\lambda,\mathbf{k}_{\parallel}}(\mathbf{r})/\sqrt{N}$
(6)

as alternative representation to Eq. (5) for Bloch states. The number of unit
cells building the crystal is denoted by $N$. Wannier functions are
essentially a real-space representation of localized orbitals and provide an
extension of the concept of atomic orbitals into solids. They have assigned as
quantum numbers the lattice vectors $\mathbf{R}_{n\parallel}$ of the $n$th
cell, where the orbital is localized, and the band index $\lambda$. Using Eq.
(5) an expression for the Bloch factor
$u_{\lambda,\mathbf{k}_{\parallel}}(\mathbf{r})$ as a function of the Wannier
orbitals $w_{\lambda}(\mathbf{r},\mathbf{R}_{n\parallel})$ reads:
$u_{\lambda,\mathbf{k}_{\parallel}}(\mathbf{r})=\sqrt{V_{UC}}\sum_{\mathbf{R}_{n\parallel}}e^{-i\mathbf{k}_{\parallel}\cdot(\mathbf{r}-\mathbf{R}_{n\parallel})}w_{\lambda}(\mathbf{r},\mathbf{R}_{n\parallel})$.
Since core electronic wave functions of neighbored atoms do barely overlap and
electron hopping between neighboring sites is unlikely, the band dispersion is
small and can be treated as a flat band. This band is built up from all the
core levels of the different atoms of the crystal to account for the
translation invariance. Therefore, irradiating a solid state should not be
seen as light being absorbed by one individual atom with a core-hole that sits
exactly at the individual atom but as a superposition of all core states
represented by a collective flat band excitation.

(ii) For valence and conduction band wave functions, typically, a full band
calculation needs to be performed Kohn1965 ; Hybertsen1984 ; Haastrup2018  and
we will use the full Bloch equation structure Eq. (5) to determine the
interaction matrix element $\Omega_{\alpha\alpha^{\prime}}$ and
$V^{\alpha\beta}_{\alpha^{\prime}\beta^{\prime}}$. These calculation yields
also an expression for the single-particle band structure
$E_{\lambda,\mathbf{k}_{\parallel}}$. At this point, the description of this
work can be connected to ab initio electronic structure theory.

(iii) The wave function of the ionization continuum can, to a good
approximation, be expressed by plane waves in vacuum orthogonalized to core
and band states. The orthogonalization ensures that all states, including
core, valence and conduction band, and ionization continuum form a complete
orthogonal basis. For the ionization continuum, where the plane wave has been
orthogonalized to the Bloch wave functions we exploite the Gram-Schmidt
processHerring1940 ; Callaway1955 ; Woodruff1957

$\displaystyle\Psi_{k_{z},\mathbf{k}_{\parallel}}(\mathbf{r})$
$\displaystyle=\frac{1}{\sqrt{V}}e^{i\mathbf{k\cdot
r}}-\sum_{\lambda}\eta_{\lambda,\mathbf{k}}\Psi_{\lambda,\mathbf{k}_{\parallel}}(\mathbf{r})$
(7)

with orthogonalization coefficient
$\eta_{\lambda,\mathbf{k}}=\int_{-\infty}^{\infty}d^{3}r~{}\exp(i\mathbf{k\cdot
r})\Psi^{*}_{\lambda,\mathbf{k}_{\parallel}}(\mathbf{r})/\sqrt{V}$. The
ionization continuum
$E_{k_{z},\mathbf{k}_{\parallel}}=\hbar^{2}\mathbf{k}^{2}/2m_{0}+E_{Ion}$
constitute the manifold of three-dimensional parabolas with energies starting
at the ionization energy $E_{Ion}$ of the material.

In principle, Eq. (5) (alternatively using Eq. (6)) and Eq. (7) constitute a
complete set of orthonormal basis functions. In the basis of the Bloch
functions Eq. (5), the Hamiltonian Eq. (3) reads

$\displaystyle H$
$\displaystyle=\sum_{\lambda,\mathbf{k}_{\parallel}}E_{\lambda,\mathbf{k}_{\parallel}}~{}a^{\dagger}_{\lambda,\mathbf{k}_{\parallel}}a^{\mathstrut}_{\lambda,\mathbf{k}_{\parallel}}$
$\displaystyle+\frac{1}{2}\sum_{\begin{subarray}{c}\lambda,\nu,\lambda^{\prime},\nu^{\prime}\\\
\mathbf{k}_{\parallel},\mathbf{q}_{\parallel},\mathbf{k}_{\parallel}^{\prime},\mathbf{q}_{\parallel}^{\prime}\end{subarray}}V^{\lambda\nu\nu^{\prime}\lambda^{\prime}}_{\mathbf{k}_{\parallel}\mathbf{q}_{\parallel}\mathbf{q}_{\parallel}^{\prime}\mathbf{k}_{\parallel}^{\prime}}~{}a^{\dagger}_{\lambda,\mathbf{k}_{\parallel}}a^{\dagger}_{\nu,\mathbf{q}_{\parallel}}a^{\mathstrut}_{\nu^{\prime},\mathbf{q}_{\parallel}^{\prime}}a^{\mathstrut}_{\lambda^{\prime},\mathbf{k}_{\parallel}^{\prime}}$
$\displaystyle-\sum_{\lambda,\lambda^{\prime},\mathbf{k}_{\parallel},\mathbf{k}_{\parallel}^{\prime},\mathbf{Q}_{\parallel}}\mathbf{d}^{\lambda\lambda^{\prime}}_{\mathbf{k}_{\parallel},\mathbf{k}_{\parallel}^{\prime}}(\bm{Q}_{\parallel})\cdot\mathbf{E}_{-\bm{Q}_{\parallel}}(z_{0},t)~{}a^{\dagger}_{\lambda,\mathbf{k}_{\parallel}}a^{\mathstrut}_{\lambda^{\prime},\mathbf{k}_{\parallel}^{\prime}}\;,$
(8)

where the $\lambda$ and $\nu$-sums include all three types of states – core,
valence, conduction band and ionization continuum. Note that we use the out-
of-plane wave vector component $k_{z}$ as quantum number for the three-
dimensional states above the ionization threshold: In case that the band index
$\lambda$ describes an unbound state we write $\lambda\rightarrow k_{z}$,
which is conform with the notation used for the wave function Eq. (7).

The Coulomb interaction in Eq. (8) couples two initial
$(\lambda,\mathbf{k}_{\parallel})$ and $(\nu,\mathbf{q}_{\parallel})$ with two
final electron states $(\lambda^{\prime},\mathbf{k}_{\parallel}^{\prime})$ and
$(\nu^{\prime},\mathbf{q}_{\parallel}^{\prime})$. The corresponding matrix
element is defined by

$\displaystyle
V^{\lambda\nu\nu^{\prime}\lambda^{\prime}}_{\mathbf{k}_{\parallel}\mathbf{q}_{\parallel}\mathbf{q}_{\parallel}^{\prime}\mathbf{k}_{\parallel}^{\prime}}=\langle\Psi_{\lambda,\mathbf{k}_{\parallel}}\Psi_{\nu,\mathbf{q}_{\parallel}}\mid
V(\mathbf{r-r^{\prime}})\mid\Psi_{\nu^{\prime},\mathbf{q}_{\parallel}^{\prime}}\Psi_{\lambda^{\prime},\mathbf{k}_{\parallel}^{\prime}}\rangle$
(9)

with the Coulomb potential $V(\mathbf{r-r^{\prime}})$.

In Eq. (8), for the light-matter interaction at the sample position the in-
plane Fourier transform of the X-ray electric field
$\mathbf{E}(\mathbf{r},t)=\sum_{\bm{Q}_{\parallel}}\mathbf{E}_{\bm{Q}_{\parallel}}(z_{0},t)\exp(-i\bm{Q}_{\parallel}\cdot\mathbf{r}_{\parallel})$,
the wave vector of the X-ray radiation $\bm{Q}_{\parallel}$ and $z_{0}$ as
sheet position are introduced. Since we investigate two-dimensional materials
we assume an atomically thin sheet lying at $z_{0}$ that the spatial variation
of the electric field in $z$-direction is negligible. This approximation can
be discussed in more detail: For two-dimensional materials the layer thickness
lies in the range of the lattice constant. For many materials such as
graphene, TMDCs or antimonene the thickness lies between 0.3 nm and 0.6 nm.
Here, we focus on soft X-ray radiation up to an energy of 500 eV corresponding
to a wavelength of 2.5 nm justifying the neglect of the spatial variation of
the electric field perpendicular to the plane. Without the loss of generality,
we set $z_{0}=0$ at the coordinate center. The dipole matrix element reads

$\displaystyle\mathbf{d}^{\lambda\lambda^{\prime}}_{\mathbf{k}_{\parallel},\mathbf{k}_{\parallel}^{\prime}}(\bm{Q}_{\parallel})=e_{0}\langle\Psi_{\lambda,\mathbf{k}_{\parallel}}\mid\mathbf{r}~{}e^{i\bm{Q}_{\parallel}\cdot\mathbf{r}_{\parallel}}\mid\Psi_{\lambda^{\prime},\mathbf{k}_{\parallel}^{\prime}}\rangle\;,$
(10)

which carries also the spatial in-plane component of the X-ray field. To
evaluate the dynamics of X-ray electron interactions, we need to determine the
interaction matrix element
$\mathbf{d}^{\lambda\lambda^{\prime}}_{\mathbf{k}_{\parallel},\mathbf{k}_{\parallel}^{\prime}}(\bm{Q}_{\parallel})$.
While the Hamiltonian Eq. (8) contains all possible electronic transitions, we
discuss at this point the matrix elements for the XANES and EXAFS processes
separately. We start with the XANES process in Sec. II.1.1 and discuss the
EXAFS transition in Sec II.1.2. The corresponding transitions are sketched in
Fig. 1(a).

####  II.1.1 XANES matrix element

In XANES the X-rays have sufficient energy to excite electrons from core band
states to unoccupied excited states below the ionization threshold (CB-
states), cf. Fig. 1(a). This process gives rise to sudden absorption edges in
the XAS spectrum. To calculate the XANES matrix element, we restrict the
quantum number $\lambda$ to bands energetically below the ionization threshold
and use the wave functions Eq. (5) in the form of the Bloch functions for the
definition of the dipole matrix element
$\mathbf{d}^{\lambda\lambda^{\prime}}_{\mathbf{k}_{\parallel},\mathbf{k}_{\parallel}^{\prime}}(\bm{Q}_{\parallel})$,
Eq. (10). This implies a momentum selection rule of
$\mathbf{k}^{\prime}_{\parallel}=\mathbf{k}_{\parallel}-\bm{Q}_{\parallel}+\mathbf{G}_{\parallel}$
for the optical transition, where the momentum is conserved up to a reciprocal
lattice vector $\mathbf{G}_{\parallel}$ and detailed in the SM Sec. II. The
wave vector $\mathbf{k}_{\parallel}$ is a reduced wave vector lying in the
first Brillouin zone. Already soft X-rays have a wave number of tens of
percent of the Brillouin zone. Therefore, when adding $\mathbf{k}_{\parallel}$
and $\bm{Q}_{\parallel}$ the resultant vector can lie outside of the first
Brillouin zone. The reciprocal lattice vector $\mathbf{G}_{\parallel}\neq 0$
accounts for Umklapp processes, refolding the resultant vector to its
equivalent wave vector in the first Brillouin zone MadelungBuch . A detailed
derivation of the XANES dipole matrix operator acting partly as a derivative
on
$a^{\dagger}_{\lambda,\mathbf{k}_{\parallel}}a^{\mathstrut}_{\lambda^{\prime},\mathbf{k}_{\parallel}^{\prime}}$
is provided in the SM. Up to this point, coming from Eq. (8), the X-ray wave
vector is defined by a general Fourier transform over the complete space.
Therefore, to be consistent with the definition of the electronic wave vector
$\mathbf{k}_{\parallel}$ we replace
$\bm{Q}_{\parallel}\rightarrow\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel}$,
where $\mathbf{Q}_{\parallel}$ is now defined within the first Brillouin zone.
We obtain

$\displaystyle\mathbf{X}^{\lambda\lambda^{\prime}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}(\mathbf{G}_{\parallel})$
$\displaystyle=-\frac{ie_{0}}{V_{UC}}\langle\begin{pmatrix}\nabla_{\mathbf{k}_{\parallel}+\mathbf{Q_{\parallel}}}\\\
iz\end{pmatrix}u_{\lambda,\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}\mid
e^{i\mathbf{G}_{\parallel}\cdot\mathbf{r}_{\parallel}}\mid
u_{\lambda^{\prime},\mathbf{k}_{\parallel}}\rangle-
ie_{0}\delta_{\lambda,\lambda^{\prime}}\nabla_{\mathbf{k}_{\parallel}+\mathbf{Q_{\parallel}}}\;.$
(11)

To clarify that the matrix operator Eq. (11) describes only X-ray induced
transitions with initial and final electronic state below the ionization
threshold (XANES) from now on the labeling
$\mathbf{d}^{\lambda\lambda^{\prime}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}(\mathbf{G}_{\parallel})\rightarrow\mathbf{X}^{\lambda\lambda^{\prime}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}(\mathbf{G}_{\parallel})$
is chosen. In Eq. (11), the first term describes interband transitions between
core and unoccupied conduction band states. XANES transitions are labelled K,
L, M, etc. depending on the principle quantum number of the initial core band
$\lambda$, based on the Siegbahn notation Siegbahn1916 ; Jenkins1991 . For
example, K transitions involve 1$s$ electrons, while excitation of 2$s$ and
2$p$ electrons occurs at L edges. The latter, as energetically following
edges, are divided into L1 for 2$s$ electrons and L2, L3 for 2$p$ electrons
(total angular momentum $J=1/2$ and $J=3/2$). The momentum conservation
accounts for a transfer of the in-plane field momentum to the optically
excited electron. The second term with the momentum gradient acting on the
creation operator in Eq. (8), leads to a wave vector gradient, which changes
according to the acceleration theorem with a rate proportional to the X-ray
field Bloch1929  and describes the coupling strength of X-ray excitation to
intraband transitions Golde2008 ; Song2020 . Since X-ray radiation includes
excitation energies from a hundred of eV up to tens of keV, the intraband
dynamics do not significantly change the optical response in a rotating wave
approximation.

Note, that the Bloch factors of the core states in Eq. (11) can also be
expressed as functions of the Wannier orbitals Eq. (6) if this representation
is desired.

####  II.1.2 EXAFS matrix element

After having discussed the XANES dipole matrix element we turn our attention
to the EXAFS matrix element between a Bloch-like core state
$\\{\mathbf{k}_{\parallel},\lambda\\}$, Eq. (5), and ionization continuum (IC)
state $\\{\mathbf{k}_{\parallel},\lambda^{\prime}\equiv k_{z}\\}$, cf. Fig. 1.
Inserting the orthogonalized plane waves Eq. (7) and the Bloch band states Eq.
(5) into the definition of the dipole matrix element Eq. (10) yields formally

$\displaystyle\mathbf{Y}^{\lambda
k_{z}}_{\mathbf{k}_{\parallel},\mathbf{k}_{\parallel}^{\prime}}(\bm{Q}_{\parallel})$
$\displaystyle=e_{0}\langle\Psi_{\lambda,\mathbf{k}_{\parallel}}\mid\mathbf{r}~{}e^{i\bm{Q_{\parallel}\cdot
r_{\parallel}}}\mid\mathbf{k}^{\prime}\rangle$
$\displaystyle+e_{0}\sum_{\lambda^{\prime}}\eta_{\lambda^{\prime},\mathbf{k}^{\prime}}\langle\Psi_{\lambda,\mathbf{k}_{\parallel}}\mid\mathbf{r}~{}e^{i\bm{Q_{\parallel}\cdot
r_{\parallel}}}\mid\Psi_{\lambda^{\prime},\mathbf{k^{\prime}}_{\parallel}}\rangle$
(12)

with the notation $\langle\mathbf{r}|\mathbf{k}\rangle=\exp(i\mathbf{k\cdot
r})/\sqrt{V}$ resulting from the unorthogonalized plane wave character of the
final state. The second contribution originates from the orthogonalization
coefficients defined below Eq. (7). To clarify that the dipole matrix element
now describes transitions into the ionization continuum it is denoted by
$\mathbf{d}^{\lambda
k_{z}}_{\mathbf{k}_{\parallel},\mathbf{k}_{\parallel}^{\prime}}(\bm{Q}_{\parallel})\rightarrow\mathbf{Y}^{\lambda
k_{z}}_{\mathbf{k}_{\parallel},\mathbf{k}_{\parallel}^{\prime}}(\bm{Q}_{\parallel})$.
The EXAFS matrix element Eq. (12) can be calculated in similar way as the
XANES matrix element, however the final state wave vector is, in contrast to
the initial state, three-dimensional. After some manipulations, detailed in
the SM Sec. III, and restricting $\bm{Q}_{\parallel}$ to the first Brillouin
zone the EXAFS matrix element is obtained as

$\displaystyle\mathbf{Y}^{\lambda
k_{z}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}(\mathbf{G}_{\parallel})$
$\displaystyle=\sum_{\lambda^{\prime}}\eta_{\lambda^{\prime},\mathbf{k}}\mathbf{X}^{\lambda\lambda^{\prime}(\text{inter})}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}$
$\displaystyle-\frac{ie_{0}}{\sqrt{V_{UC}}}\langle\begin{pmatrix}\nabla_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}\\\
iz\end{pmatrix}u_{\lambda,\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}\mid
e^{i\mathbf{G}_{\parallel}\cdot\mathbf{r}_{\parallel}}e^{ik_{z}z}\rangle\;.$
(13)

Similar to Eq. (11) (XANES), we find an interband term describing the optical
transition between a Bloch and a plane wave state (second line). Note that we
consider for EXAFS only transitions from the material into the ionization
continuum. Therefore the matrix element Eq. (13) contains no intraband
interaction because of the forced orthogonality of the ionization continuum to
all band states. Consequently, the appearing XANES matrix element in Eq. (13),
stemming from the orthogonalization, is restricted to its interband part. In
principle, also an acceleration of the free electrons is possible and a third
contribution to the matrix element
$\mathbf{Y}^{k_{z}k_{z}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}(\mathbf{G}_{\parallel})$
is present, which is proportional to the X-ray field times the momentum
gradient of the electronic creation operator Golde2008 . However, since this
is no EXAFS transition, goes beyond our description, and is not visible in the
absorption spectrum this contribution is neglected.

####  II.1.3 Hamiltonian, discriminating XANES and EXAFS processes

Using the calculated matrix elements for the XANES and EXAFS processes, the
Hamiltonian containing the electron dispersion and the interaction of X-ray
electric field with the crystal electrons reads

$\displaystyle H$
$\displaystyle=\sum_{\lambda,\mathbf{k}_{\parallel}}E_{\lambda,\mathbf{k}_{\parallel}}~{}a^{\dagger}_{\lambda,\mathbf{k}_{\parallel}}a^{\mathstrut}_{\lambda,\mathbf{k}_{\parallel}}+\sum_{\mathbf{k}_{\parallel},k_{z}}E_{k_{z},\mathbf{k}_{\parallel}}~{}a^{\dagger}_{k_{z},\mathbf{k}_{\parallel}}a^{\mathstrut}_{k_{z},\mathbf{k}_{\parallel}}$
$\displaystyle-\sum_{\begin{subarray}{c}\lambda,\lambda^{\prime}\\\
\mathbf{k}_{\parallel},\mathbf{Q}_{\parallel},\mathbf{G}_{\parallel}\end{subarray}}\mathbf{X}^{\lambda\lambda^{\prime}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}(\mathbf{G}_{\parallel})\cdot\mathbf{E}_{-\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel}}(t)~{}a^{\dagger}_{\lambda,\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}a^{\mathstrut}_{\lambda^{\prime},\mathbf{k}_{\parallel}}$
$\displaystyle-\sum_{\begin{subarray}{c}\lambda,k_{z}\\\
\mathbf{k}_{\parallel},\mathbf{Q}_{\parallel},\mathbf{G}_{\parallel}\end{subarray}}\left(\mathbf{Y}^{\lambda
k_{z}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}(\mathbf{G}_{\parallel})\cdot\mathbf{E}_{-\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel}}(t)~{}a^{\dagger}_{\lambda,\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}a^{\mathstrut}_{k_{z},\mathbf{k}_{\parallel}}\right.$
$\displaystyle\left.+\text{H.c.}\right)$
$\displaystyle+\frac{1}{2}\sum_{\begin{subarray}{c}\lambda,\nu,\lambda^{\prime},\nu^{\prime}\\\
\mathbf{k}_{\parallel},\mathbf{q}_{\parallel},\mathbf{k}_{\parallel}^{\prime},\mathbf{q}_{\parallel}^{\prime}\end{subarray}}V^{\lambda\nu\nu^{\prime}\lambda^{\prime}}_{\mathbf{k}_{\parallel}\mathbf{q}_{\parallel}\mathbf{q}_{\parallel}^{\prime}\mathbf{k}_{\parallel}^{\prime}}~{}a^{\dagger}_{\lambda,\mathbf{k}_{\parallel}}a^{\dagger}_{\nu,\mathbf{q}_{\parallel}}a^{\mathstrut}_{\nu^{\prime},\mathbf{q}_{\parallel}^{\prime}}a^{\mathstrut}_{\lambda^{\prime},\mathbf{k}_{\parallel}^{\prime}}\;.$
(14)

The unpertubed single particle energies in a lattice periodic atomic potential
are described by the first line of Eq. (14). The crystal band index $\lambda$
in the first term incorporates all bound bands, starting at the core level
1$s$ state up to the comparable delocalized valence and conduction band
states. The second line describes XANES transitions between two-dimensional
core and conduction band states. The third line describes EXAFS transitions
between two-dimensional initial and three-dimensional final states lying above
the ionization threshold of the crystal. It should be remembered, that $k_{z}$
describes the out-of-plane wave vector and acts as quantum number for the
unbound states in a three-dimensional continuum. The last line describes
Coulomb interaction between the carriers. It couples two initial states
$(\lambda,\mathbf{k}_{\parallel})$, $(\nu,\mathbf{q}_{\parallel})$ and two
final states $(\lambda^{\prime},\mathbf{k}_{\parallel}^{\prime})$,
$(\nu^{\prime},\mathbf{q}_{\parallel}^{\prime})$. Computing the Coulomb matrix
element as is shown in the SM we can also derive momentum selection rules
holding during Coulomb scattering. In detail, we can use
$\sum_{\mathbf{G}_{\parallel}}\delta_{\mathbf{k}_{\parallel},\mathbf{k}_{\parallel}^{\prime}+\mathbf{p}_{\parallel}-\mathbf{G}_{\parallel}}$
and
$\sum_{\mathbf{G}_{\parallel}^{\prime}}\delta_{\mathbf{q}_{\parallel},\mathbf{q}_{\parallel}^{\prime}-\mathbf{p}_{\parallel}-\mathbf{G}_{\parallel}^{\prime}}$,
where $\mathbf{p}_{\parallel}$ describes the momentum transfer during the
carrier-carrier interaction and the reciprocal lattice vectors
$\mathbf{G}_{\parallel}$ and $\mathbf{G}_{\parallel}^{\prime}$ takes into
account for possible Umklapp processes. The Coulomb selection rules can be
inserted in the Hamiltonian Eq. (14) if Coulomb scattering wants to be
evaluated explicitly.

###  II.2 Dipole matrix elements in tight binding approximation

The character of the dipole transitions, Eq. (11) and Eq. (13), including the
lattice geometry, which is of essential importance for X-ray experiments in
the solid state, as well as the orbital band composition are encoded in the
lattice periodic function $u_{\lambda,\mathbf{k}_{\parallel}}(\mathbf{r})$.
Access to the Bloch functions (or Wannier functions) typically requires
computationally expensive methods. To evaluate explicitly the matrix element
and to obtain more analytical insights into X-ray induced electronic
transitions we use the tight binding method Jones1934 ; Slater1954 , in
chemistry usually referred to as Hückel theory Huckel1931 , for core and
conduction band wave function. The method consists of the assumption to
approximate the Wannier functions with atomic orbitals. The Bloch states
$\Psi_{\lambda,\mathbf{k}_{\parallel}}(\mathbf{r})$ are then expanded in terms
of a linear combination of atomic orbitals
$\langle\mathbf{r}|\lambda,\beta,j,0\rangle=\phi^{\lambda}_{\beta,j}(\mathbf{r})$
(orbital of type $j$ on atom $\beta$) of the composing atoms, which are
eigenstates of the single atom Hamiltonian:

$\displaystyle\Psi_{\lambda,\mathbf{k}_{\parallel}}(\mathbf{r})=\frac{1}{\sqrt{N}}\sum_{\beta,j,\mathbf{R}_{\beta\parallel}}C_{\beta
j,\mathbf{k}_{\parallel}}^{\lambda}e^{i\mathbf{k}_{\parallel}\cdot\mathbf{R}_{\beta\parallel}}\phi^{\lambda}_{\beta,j}(\mathbf{r}-\mathbf{R}_{\beta\parallel}).$
(15)

The atomic orbitals spatially decay with a typical constant of
$\zeta=Z/a_{B}$, where $a_{B}$ denotes the hydrogen Bohr radius and $Z$ is the
effective nuclear charge, incorporating electron screening effects. Values for
an effective nuclear charge $Z$ in atoms or ions is provided by the Slater
rule Slater1930 , which can be used as a first estimate also in crystals. For
more accurate values first-principle calculations are necessary. In Eq. (15),
the coordinates of the atoms in the crystal lattice are denoted by
$\mathbf{R}_{\beta\parallel}$ and the tight binding coefficients
$C^{\lambda}_{\beta j,\mathbf{k}_{\parallel}}$ determine the weight of the
different orbitals. For core states the atomic orbitals decay rapidly from the
atomic position that their overlap is negligible small. In the following we
revisit the XANES and EXAFS matrix elements exploiting the wave function Eq.
(15) for core and conduction band wave function and derive the solid state
electron-X-ray Hamiltonian in tight binding approximation.

####  II.2.1 XANES

The XANES matrix element is evaluated by inserting the tight binding wave
function Eq. (15) for initial and final state into the formal definition Eq.
(10). Expanding the integral into a sum over unit cells at the lattice vector
$\mathbf{R}_{\beta\parallel}$ leads to the matrix element

$\displaystyle\mathbf{X}^{\lambda\lambda^{\prime}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}(\mathbf{G}_{\parallel})$
$\displaystyle=-e_{0}\sum_{\alpha,\beta,i,j}\sum_{\bm{\delta}_{\beta\alpha}}C^{*\lambda}_{\beta
j,\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}C^{\lambda^{\prime}}_{\alpha
i,\mathbf{k}_{\parallel}}e^{i\mathbf{k_{\parallel}}\cdot\bm{\delta}_{\beta\alpha}}\langle\lambda,\beta,j,0\mid\mathbf{r}~{}e^{i(\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel})\cdot\mathbf{r}_{\parallel}}\mid\lambda^{\prime},\alpha,i,\bm{\delta}_{\beta\alpha}\rangle$
$\displaystyle-
ie_{0}\sum_{\alpha,\beta,i,j}\sum_{\bm{\delta}_{\beta\alpha}}C^{*\lambda}_{\beta
j,\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}C^{\lambda^{\prime}}_{\alpha
i,\mathbf{k}_{\parallel}}e^{i\mathbf{k}_{\parallel}\cdot\bm{\delta}_{\beta\alpha}}\langle\lambda,\beta,j,0\mid
e^{i(\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel})\cdot\mathbf{r}_{\parallel}}\mid\lambda^{\prime},\alpha,i,\bm{\delta}_{\beta\alpha}\rangle\nabla_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}$
(16)

including the same momentum selection rule as previously, Eq. (11) and using
that $\mathbf{Q}_{\parallel}$ is restricted to the first Brillouin zone. The
projection onto atomic orbital basis yields
$\langle\mathbf{r}|\lambda,\beta,j,\bm{\delta}_{\beta\alpha}\rangle=\phi_{\beta
j}^{\lambda}(\mathbf{r}-\bm{\delta}_{\beta\alpha})$. The electron momentum-
dependence is carried out by the tight binding coefficients and the phase
factor $\exp(i\mathbf{k}_{\parallel}\cdot\bm{\delta}_{\beta\alpha})$, where
$\bm{\delta}_{\beta\alpha}=\mathbf{R}_{\beta\parallel}-\mathbf{R}_{\alpha\parallel}$
stands for the next-neighbor vectors connecting the atoms. The sum over the
sublattices $\alpha$, $\beta$ includes neighboring atoms of arbitrary order.
Since we investigate the interaction of strongly localized core electrons,
their spatial extension decays rapidly even compared to the wavelength of soft
and medium X-rays up to an energy of $\sim 4$-5 keV. Therefore, it is
reasonable to treat the transition integral in dipole approximation, i.e.
perform an expansion of the radiation field in zeroth order
$\exp(i(\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel})\cdot\mathbf{r}_{\parallel})\approx
1$. We have carefully checked, see SM Sec. IV, that the inclusion of the
$\mathbf{Q}_{\parallel}$-dependence to the transition integral does not change
the result at this point. Since the spatial localization of the core electrons
increases with the atomic weight, the dipole approximation becomes better the
heavier the constituting atoms are. Treating the transition integral within
the dipole approximation, the optical selection rules known from atomic
spectroscopy are recovered. However, thanks to the tight binding wave function
of the crystalline solid and the plane wave decomposition in Eq. (16) we
include solid state properties, in particular the lattice periodicity. Because
the core orbitals are generally more localized compared to other states the
optical transitions involving core states are generally weaker than
transitions between valence and conduction band adressed in optical
experiments with visible light. Applying the dipole approximation to the
second line of Eq. (16) the integral turns into an overlap integral
$\langle\lambda,\beta,j,0\mid\lambda^{\prime},\alpha,i,\bm{\delta}_{\beta\alpha}\rangle$.
The overlap is generally small and we may assume that the chosen orbitals are
orthogonal to each other showing that the second line in Eq. (16) describes
intraband transitions similar to Eq. (11).

####  II.2.2 EXAFS

As for the XANES transitions, we evaluate the EXAFS matrix element for tight
binding wave functions. The starting point is the EXAFS dipole matrix element
Eq. (12), now with the initial state electronic wave function
$\Psi_{\lambda,\mathbf{k}_{\parallel}}(\mathbf{r})$ in the tight binding
approach. First we investigate the orthogonalization coefficients defined
below Eq. (7), which appears as a sum involving all bound bands. With the main
$n$, angular $l$ and magnetic $m$ quantum numbers, which determine the atomic
orbitals, the coefficients can be calculated to be
$\eta_{\lambda,\mathbf{k}}=i^{l}\sqrt{(2\zeta)^{2n+1}/(2n)!}J_{n,l}(k)Y_{lm}(\vartheta_{k},\varphi_{k})/2\pi^{2}$
with the spherical harmonics $Y_{lm}$ and the radial function
$J_{n,l}(k)=\sqrt{\pi/2k}(\zeta^{2}+k^{2})^{-(n+1)/2}\Gamma(n-1+l)P^{-l}_{n}[\zeta(\zeta^{2}+l^{2})^{-1/2}]$
with the associated Legendre polynomials $P_{n}^{l}(x)$. The orthogonalization
coefficient to the 1$s$ state is proportional to
$\eta_{1s,\mathbf{k}}\propto(\zeta^{2}+k^{2})^{-2}$. Here,
$\eta_{1s,\mathbf{k}}$ decreases to a value of 10 % at approximately
$k=\sqrt{2}\zeta$. Generally, because the 1$s$-orbital has the smallest extent
in real space, its orthogonalization contribution determines the strength of
the orthogonalization contribution to the EXAFS dipole matrix element.
Further, we see that the orthogonalization coefficient depends on the atomic
number and gains importance with increasing atomic weight. With an exemplary
effective atomic number of $Z=5.7$, stemming from the Slater rule for carbon
Slater1930 , we obtain a wave number of 140 nm-1 corresponding to an energy of
5 eV. We see that the coefficients decrease rapidly with increasing energy
(starting at the ionization threshold) since the orthogonalization contributes
only close to the surface. Since the coefficients $\eta_{\lambda,\mathbf{k}}$
decrease rapidly to zero with increasing energy, we use in the following free
electronic continuum states for simplicity and obtain for the first term of
Eq. (12)

$\displaystyle\mathbf{Y}^{\lambda
k_{z}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}(\mathbf{G}_{\parallel})$
$\displaystyle=\frac{-e_{0}}{\sqrt{V_{UC}}}\sum_{\beta,j}C^{*\lambda}_{\beta
j,\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}$
$\displaystyle\times\langle\lambda,\beta,j,0\mid\mathbf{r}\mid\mathbf{k}+\mathbf{Q}_{\parallel}-\mathbf{G}_{\parallel}\rangle\;.$
(17)

However, the plane wave approximation (17) has to be investigated carefully
for each material independently.

To illustrate the approach, we focus on K-shell transitions, namely
$\lambda=1s$ for the initial state. The calculation of the transition integral
from a 1$s$ core electron to the plane wave state is included in the SM. We
obtain

$\displaystyle\mathbf{Y}^{1sk_{z}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}(\mathbf{G}_{\parallel})$
$\displaystyle=\frac{-e_{0}}{\sqrt{V_{UC}}}\sum_{\beta,n}C^{*1s}_{\beta
1s,\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}$
$\displaystyle\times\frac{32\sqrt{\pi}i\zeta^{5/2}(\mathbf{k}+\mathbf{Q}_{\parallel}-\mathbf{G}_{\parallel})_{n}}{(\zeta^{2}+(\mathbf{k}+\mathbf{Q}_{\parallel}-\mathbf{G}_{\parallel})^{2})^{3}}\hat{\mathbf{e}}_{n}$
(18)

with the unit vector $\hat{\mathbf{e}}_{n}$ in cartesian coordinates. Equation
(18) consists of the tight-binding coefficients $C^{*1s}_{\beta
1s,\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}$, which carry the lattice
periodicity, and the transition integral (second line), which carries the unit
cell information. The magnitude and width of this form factor peaks are
determined by the spatial electron distribution. Therefore, the form factor is
weakened by the effective inverse Bohr radius $\zeta$, cf. Eq. (18)
denominator. Further, we see that the optical transition into the ionization
continuum is unpolarized in contrast to transitions within the material.
Therefore, the optical transition of internal states into the vacuum are
independent of the incident angle of the light.

Figure 3: Sketch of possible Meitner-Auger transitions included in the
Hamiltonian. (a) Relaxation of a conduction band electron and core-hole
accompanied by an excitation of a core/valence band electron (blue arrows
depending on Pauli blocking) or by an excitation of conduction band electrons
(yellow arrows) previously injected by the X-ray excitation. (b) Recombination
of a core-hole with a core electron accompanied by an interband excitation of
a core/valence electron.

It is interesting to compare X-ray absorption to the complementary technique
X-ray diffraction used to measure the structure of crystals. The difference
between EXAFS and X-ray diffraction lies in the choice of the observable:
While for the EXAFS a transmission or reflection is measured under the same
angle as the X-ray incident angle, for X-ray diffraction also the signal under
a different angle than the incident is measured. In particular, we have
checked that the conventional X-ray diffraction of the ground state is
included in our description for $\omega\rightarrow\infty$. The diffraction is
determined by the Fourier transform of the electron distribution and often
called form factor. As consequence the X-ray propagation involves a reduced
sensitivity to light atoms with stronger localized and only weakly screened
inner shell electronsWang2020 . Moving to electron diffraction the form factor
can be calculated from the X-ray form factor by the Mott-Bethe equation
Mott1929 ; Bethe1930 , which takes additionally to the elastic scattering at
the electron clouds also nucleus scattering into account. Then also crystals
of lighter element can be resolved and measured.

##  III X-ray Bloch equations

The observable describing the X-ray response of the material is the detected
X-ray field (Sec. IV), which can be calculated from the incident field
interfering with the excited dipole density in the sample in reflection or
transmission geometry. To derive the dipole density $\mathbf{P}(\mathbf{r},t)$
as a function of the electric field $\mathbf{E}(\mathbf{r},t)$, we start from
the light-matter interaction Hamiltonian (cf. SM). From Eq. (8) we can
identify the two-dimensional macroscopic polarization density

$\displaystyle\mathbf{P}^{2D}_{\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel}}(t)$
$\displaystyle=\frac{1}{A}\sum_{\lambda_{1}\neq\lambda_{2},\mathbf{k}_{\parallel}}\mathbf{d}^{\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}p^{\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}(t)\>,$
(19)

which determines the X-ray response (Sec. IV) in wave number space. Here, we
explicitly exclude the intraband contribution as discussed in Sec. II.1.
Please note that $\mathbf{Q}_{\parallel}$ is defined within the first
Brillouin zone. For the X-ray dynamics, the relevant quantities occurring in
Eq. (19) are the transition
$p^{\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}(t)=\langle
a^{\dagger}_{\lambda_{1},\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}a^{\mathstrut}_{\lambda_{2},\mathbf{k}_{\parallel}}\rangle(t)$
defined as expectation value for a single electronic transition between the
states $|\lambda_{1},\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}\rangle$ and
$|\lambda_{2},\mathbf{k}_{\parallel}\rangle$, which can be excited if the
corresponding dipole matrix element
$\mathbf{d}^{\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}$
does not vanish. Due to the spatial resolution of the X-ray light and a
possibly non-orthogonal incidence, the non-diagonal character of the
transitions and occupations in momentum space ($\mathbf{Q}_{\parallel}\neq 0$)
have to be included. The electronic transitions can be derived from the
Heisenberg equation of motion using the Hamiltonian Eq. (14). The
corresponding equation of motion for the microscopic transition reads

$\displaystyle
i\hbar\frac{d}{dt}p^{\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}$
$\displaystyle=\left(E^{\lambda_{2}}_{\mathbf{k}_{\parallel}}-E^{\lambda_{1}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}\right)p^{\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}$
$\displaystyle-\sum_{\lambda,\mathbf{Q}_{\parallel}^{\prime},\mathbf{G}_{\parallel}}\mathbf{E}_{-\mathbf{Q}_{\parallel}^{\prime}+\mathbf{G}_{\parallel}}(t)\cdot\left(\mathbf{d}^{\lambda_{2}\lambda}_{\mathbf{k}_{\parallel},\mathbf{k}_{\parallel}-\mathbf{Q}_{\parallel}^{\prime}}(\mathbf{G}_{\parallel})\sigma^{\lambda_{1}\lambda}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}-\mathbf{Q}_{\parallel}^{\prime}}-\mathbf{d}^{\lambda\lambda_{1}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}+\mathbf{Q}_{\parallel}^{\prime},\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}(\mathbf{G}_{\parallel})\sigma^{\lambda\lambda_{2}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}+\mathbf{Q}_{\parallel}^{\prime},\mathbf{k}_{\parallel}}\right)$
$\displaystyle+\sum_{\begin{subarray}{c}\mu,\lambda,\nu\\\
\mathbf{k}_{\parallel}^{\prime},\mathbf{q}_{\parallel},\mathbf{p}_{\parallel}\end{subarray}}\left(V^{\lambda_{2}\mu\lambda\nu}_{\mathbf{k}_{\parallel},\mathbf{k}_{\parallel}^{\prime},\mathbf{q}_{\parallel},\mathbf{p}_{\parallel}}-V^{\lambda_{2}\mu\nu\lambda}_{\mathbf{k}_{\parallel},\mathbf{k}_{\parallel}^{\prime},\mathbf{p}_{\parallel},\mathbf{q}_{\parallel}}\right)\sigma^{\mu\nu}_{\mathbf{k}_{\parallel}^{\prime},\mathbf{p}_{\parallel}}\sigma^{\lambda_{1}\lambda}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{q}_{\parallel}}$
$\displaystyle-\sum_{\begin{subarray}{c}\mu,\lambda,\nu\\\
\mathbf{k}_{\parallel}^{\prime},\mathbf{q}_{\parallel},\mathbf{p}_{\parallel}\end{subarray}}\left(V^{\mu\lambda\nu\lambda_{1}}_{\mathbf{k}_{\parallel}^{\prime},\mathbf{q}_{\parallel},\mathbf{p}_{\parallel},\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}-V^{\lambda\mu\nu\lambda_{1}}_{\mathbf{q}_{\parallel},\mathbf{k}_{\parallel}^{\prime},\mathbf{p}_{\parallel},\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}\right)\sigma^{\lambda\nu}_{\mathbf{q}_{\parallel},\mathbf{p}_{\parallel}}\sigma^{\mu\lambda_{2}}_{\mathbf{k}_{\parallel},\mathbf{k}_{\parallel}}$
$\displaystyle+i\hbar\partial_{t}p^{\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}\big{|}_{\text{coll}}\;.$
(20)

Equation (20) describes the dynamics of an X-ray induced electronic
transitions
$p^{\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}$
within a two-dimensional material defining the response in Eq. (19). The first
term incorporates the free oscillation of the transition with the single-
particle energies of initial and final state. The second line of Eq. (20)
describes the coupling to the X-ray field. For now,
$\sigma^{\lambda\lambda^{\prime}}_{\mathbf{k}_{1},\mathbf{k}_{2}}=\langle
a^{\dagger}_{\lambda,\mathbf{k}_{1}}a^{\mathstrut}_{\lambda^{\prime},\mathbf{k}_{2}}\rangle$
describes a general density matrix element. In case that
$\lambda\neq\lambda^{\prime}$ this expectation value describes a non-local
($\mathbf{k}_{1}\neq\mathbf{k}_{2}$) interband transition
$p^{\lambda\lambda^{\prime}}_{\mathbf{k}_{1},\mathbf{k}_{2}}$. In contrast,
for $\lambda=\lambda^{\prime}$ the quantity stands for a non-local electron
occupation $f^{\lambda}_{\mathbf{k}_{1},\mathbf{k}_{2}}$. The two different
wave number indices indicate the existence of spatial correlations in a
spatially inhomogeneous system. Therefore, Eq. (20) allows for the description
of translational non-invariant, spatially localized, X-ray excitation of the
crystal. The last three lines describe the Coulomb contribution. The many-
particle interaction leads to a coupling to the dynamics of higher-order
expectation values, known as hierarchy problem, which is treated by exploiting
the cluster expansion Axt1994 ; Lindberg1994 ; Fricke1996 . The action of the
Coulomb interaction is divided in Hartree-Fock (third and fourth line) and
collision contributions
$p^{\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}\big{|}_{\text{coll}}$
including many-body interaction and scattering-induced dephasing beyond the
Hartree-Fock level Kochbuch . Here, we see that our approach includes
naturally many-body interaction in a self-consistent way: Depending on the
band index combination of the Coulomb matrix element, the third and fourth
line implicitly include different kinds of Coulomb interaction mechanisms. For
example: i) when all band indices are equal we find a band renormalization.
ii) In case of two $\lambda_{1}$ and two $\lambda_{2}$ indices we can have
electron-core-hole interaction giving rise to core-excitons. iii) With three
occupied and one unoccupied band indices (or vise versa) the third and fourth
line describe Meitner-Auger-type interaction. Meitner-Auger transitions are
recombinations of electron and core-hole accompanied by an energetic elevation
of a second electron. Here, depending on the excitation conditions and band
structure, one has to carefully investigate, which terms contribute to the
dynamics of the X-ray induced transition. In principle, different kinds of
Meitner-Auger transitions are included in the Hamiltonian. Figure 3 shows
exemplary de-excitations of the core-hole. Via Coulomb coupling this interband
transitions can excite either transitions within the core and valence bands
(blue arrows) or excite transitions within the conduction bands (yellow
arrows), cf. Fig. 3(a). In the following we term transitions within the
core/valence bands or conduction/vacuum bands as intersubband transitions to
distinguish from interband transitions, where an electron changes between
occupied core/valence and unoccupied conduction/vacuum bands. Since the X-ray
irradiation can also excite electron occupations, discussed at the end of the
section (cf. Eq. (25)), a core-hole recombination can also excite intersubband
transitions in the conduction band as sketched in Fig. 3(a) via Coulomb
coupling. However, the core-hole can also recombine with a core-electron from
another shell and excite a transition of a core-hole to the conduction bands
or even out from the sample, cf. Fig. 3(b). By inverting the arrows the
complex conjugated processes are sketched, which are also included in the
Hamiltonian and are called impact ionization.

For spectrally sufficient sharp X-ray pulses exciting only the transition
between the bands $\lambda_{1}$ and $\lambda_{2}$ we can obtain more
analytical insights from Eq. (20):

$\displaystyle
i\hbar\frac{d}{dt}p^{\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}$
$\displaystyle=\left(E^{\lambda_{2}}_{\mathbf{k}_{\parallel}}-E^{\lambda_{1}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}\right)p^{\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}$
$\displaystyle-\sum_{\mathbf{Q}_{\parallel}^{\prime},\mathbf{G}_{\parallel}}\mathbf{E}_{-\mathbf{Q}^{\prime}_{\parallel}+\mathbf{G}_{\parallel}}(t)\cdot\left(\mathbf{d}^{\lambda_{2}\lambda_{1}}_{\mathbf{k}_{\parallel},\mathbf{k}_{\parallel}-\mathbf{Q}_{\parallel}^{\prime}}(\mathbf{G}_{\parallel})f^{\lambda_{1}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}-\mathbf{Q}_{\parallel}^{\prime}}-\mathbf{d}^{\lambda_{2}\lambda_{1}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}+\mathbf{Q}_{\parallel}^{\prime},\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}(\mathbf{G}_{\parallel})f^{\lambda_{2}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}+\mathbf{Q}_{\parallel}^{\prime},\mathbf{k}_{\parallel}}\right)$
$\displaystyle+\sum_{\mathbf{k}_{\parallel}^{\prime},\mathbf{q}_{\parallel},\mathbf{p}_{\parallel}}\left(V^{\lambda_{1}\lambda_{1}\lambda_{1}\lambda_{1}}_{\mathbf{k}_{\parallel}^{\prime},\mathbf{q}_{\parallel},\mathbf{p}_{\parallel},\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}f^{\lambda_{1}}_{\mathbf{k}_{\parallel}^{\prime},\mathbf{p}_{\parallel}}p^{\lambda_{1}\lambda_{2}}_{\mathbf{q}_{\parallel},\mathbf{k}_{\parallel}}-V^{\lambda_{2}\lambda_{2}\lambda_{2}\lambda_{2}}_{\mathbf{k}_{\parallel},\mathbf{k}_{\parallel}^{\prime},\mathbf{q}_{\parallel},\mathbf{p}_{\parallel}}f^{\lambda_{2}}_{\mathbf{k}_{\parallel}^{\prime},\mathbf{p}_{\parallel}}p^{\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{q}_{\parallel}}\right)$
$\displaystyle-\sum_{\mathbf{k}_{\parallel}^{\prime},\mathbf{q}_{\parallel},\mathbf{p}_{\parallel}}\left(V^{\lambda_{2}\lambda_{1}\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel},\mathbf{k}_{\parallel}^{\prime},\mathbf{q}_{\parallel},\mathbf{p}_{\parallel}}f^{\lambda_{1}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{q}_{\parallel}}p^{\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel}^{\prime},\mathbf{p}_{\parallel}}-V^{\lambda_{1}\lambda_{2}\lambda_{2}\lambda_{1}}_{\mathbf{k}_{\parallel}^{\prime},\mathbf{q}_{\parallel},\mathbf{p}_{\parallel},\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}f^{\lambda_{2}}_{\mathbf{q}_{\parallel},\mathbf{k}_{\parallel}}p^{\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel}^{\prime},\mathbf{p}_{\parallel}}\right)$
$\displaystyle+\sum_{\mu,\lambda,\mathbf{k}_{\parallel}^{\prime},\mathbf{q}_{\parallel},\mathbf{p}_{\parallel}}\left[V^{\lambda_{2}\mu\lambda\lambda_{1}}_{\mathbf{q}_{\parallel},\mathbf{k}_{\parallel}^{\prime},\mathbf{p}_{\parallel},\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}p^{\lambda_{2}\lambda}_{\mathbf{q}_{\parallel},\mathbf{p}_{\parallel}}p^{\mu\lambda_{2}}_{\mathbf{k}_{\parallel}^{\prime},\mathbf{k}_{\parallel}}-V^{\lambda_{2}\lambda\mu\lambda_{1}}_{\mathbf{k}_{\parallel},\mathbf{q}_{\parallel},\mathbf{k}_{\parallel}^{\prime},\mathbf{p}_{\parallel}}p^{\lambda\lambda_{1}}_{\mathbf{q}_{\parallel},\mathbf{p}_{\parallel}}p^{\lambda_{1}\mu}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}^{\prime}}\right]$
$\displaystyle+\sum_{\mu,\lambda,\mathbf{k}_{\parallel}^{\prime},\mathbf{q}_{\parallel},\mathbf{p}_{\parallel}}\left[V^{\lambda\mu\lambda_{2}\lambda_{1}}_{\mathbf{q}_{\parallel},\mathbf{k}_{\parallel}^{\prime},\mathbf{p}_{\parallel},\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}p^{\lambda\lambda_{2}}_{\mathbf{q}_{\parallel},\mathbf{p}_{\parallel}}p^{\mu\lambda_{2}}_{\mathbf{k}_{\parallel}^{\prime},\mathbf{k}_{\parallel}}-V^{\lambda_{2}\lambda_{1}\mu\lambda}_{\mathbf{k}_{\parallel},\mathbf{q}_{\parallel},\mathbf{k}_{\parallel}^{\prime},\mathbf{p}_{\parallel}}p^{\lambda_{1}\lambda}_{\mathbf{q}_{\parallel},\mathbf{p}_{\parallel}}p^{\lambda_{1}\mu}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}^{\prime}}\right]$
$\displaystyle+i\hbar\partial_{t}p^{\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}\big{|}_{\text{coll}}\;.$
(21)

The second line in Eq. (21) describes the coupling to the X-ray field. The
X-ray transition is initiated by the core occupation
$f^{\lambda_{1}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}-\mathbf{Q}_{\parallel}^{\prime}}$
and blocked by the final band occupation
$f^{\lambda_{2}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}+\mathbf{Q}_{\parallel}^{\prime},\mathbf{k}_{\parallel}}$.
The third line describes an energy renormalization, due to intraband Coulomb
interaction, for both bands. They are well-known in the literature as a core-
hole renormalization Oji1998 ; Mizoguchi2000 ; Mauchamp2009  and lead to an
effective energetic blue-shift of the transition energy
$(E^{\lambda_{2}}_{\mathbf{k}_{\parallel}}-E^{\lambda_{1}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}})$
of the electronic transition Eq. (21) for occupied valence states. The fourth
line describes a Coulomb-induced renormalization of the Rabi-frequency
$\mathbf{E}_{\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel}}(t)\cdot\mathbf{d}^{\lambda_{2}\lambda_{1}}_{\mathbf{k}_{\parallel},\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}(\mathbf{G}_{\parallel})$,
which can be interpreted as local field contribution and leads to the
formation of excitons Bechstedt1980 ; Olovsson2009 ; Christiansen2019 . The
crystalline excitons are quasi-particles built up by a Coulomb-correlated wave
number distribution of core-holes and electrons by two occupied/unoccupied
bands, respectively. The equation of motion Eq. (21) implicitly includes the
Wannier equation Christiansen2019 ; Katsch2018  describing the formation of
core-excitons bleached by the occupation difference between electrons and
core-holes. Depending on the bleaching this contribution to the equation of
motion leads to excitonic peaks below the absorption edge. Such bound
electron-core-hole quasi-particle can play a considerable role in the
interpretation of X-ray spectroscopy Hjalmarson1981 ; Olovsson2009 ;
Cocchi2015 ; Begum2021 ; Chang2021 ; Grunes1983 , for example in the case of
metal oxides Grunes1983 ; Biswas2018 ; Geneaux2020 . The interpretation of the
corresponding two lines can be verified in the limit of a spatially
homogeneous system Kirabuch  and are discussed in the SM Sec. V. An
advantageous representation for space-dependent phenomena is the Wigner
representation HessKuhn1996 . The occupations and transitions can be Fourier
transformed with respect to their relative momentum. Performing a gradient
expansion of the Fourier phase factor and going beyond the zeroth order yields
spatially-resolved Bloch equations for occupation and transition in an
inhomogeneous system HessKuhn1996 . The fifth and sixth line include Meitner-
Auger-like terms in the Hartree-Fock approximation. As previously discussed,
Meitner-Auger interaction couples interband transitions between core and
conduction band to transitions within the occupied (or unoccupied) bands. In
line five $\mu$ and $\lambda$ need to correspond either both to occupied or
unoccupied bands. This leads to a product of two (intersubband) transitions
between different conduction or valence bands. Therefore, line five gives rise
to nonlinearities in second order. For the first term of the sixth line the
band indices have to correspond to two occupied bands, while for the second
both need to be unoccupied bands. Also here, the source is quadratic in
intersubband transitions and therefore goes beyond the linear optics limit.
Within a rotating wave approximation only sources on the right hand side of
Eq. (21) contribute, which oscillation energy matches the energy difference of
$\lambda_{1}$ and $\lambda_{2}$ of the left hand side of Eq. (21). Line five
and six include the generation of second harmonics resulting from the
excitation with intense fieldsGarnik2022 .

So far, for the Coulomb contribution no assumption concerning the momentum
selection rules were made. Making use of the momentum conservation law
$V^{\mu\lambda\nu\rho}_{\mathbf{k}_{\parallel},\mathbf{q}_{\parallel},\mathbf{p}_{\parallel},\mathbf{k}^{\prime}_{\parallel}}=\sum_{\mathbf{q}^{\prime}_{\parallel},\mathbf{G}_{\parallel},\mathbf{G}_{\parallel}^{\prime}}V^{\mu\lambda\nu\rho}_{\mathbf{k}_{\parallel},\mathbf{q}_{\parallel},\mathbf{p}_{\parallel},\mathbf{k}^{\prime}_{\parallel}}\delta_{\mathbf{G}_{\parallel},\mathbf{k}_{\parallel}^{\prime}-\mathbf{k}_{\parallel}+\mathbf{q}^{\prime}_{\parallel}}\delta_{\mathbf{G}_{\parallel}^{\prime},\mathbf{p}_{\parallel}-\mathbf{q}_{\parallel}-\mathbf{q}^{\prime}_{\parallel}}$
the last four lines (see SM Sec. VI for a detailed derivation) can be written
in a more convenient form. Higher-order contributions to the Coulomb
interaction, which are included in the collision term, need to be treated on
the same level. The contribution
$\partial_{t}p^{\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}\big{|}_{\text{coll}}$
describes the many-particle scattering beyond the Hartree-Fock interaction.
The collision term contributes to diagonal
$\gamma_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}$
and off-diagonal
$\mathcal{U}^{\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}$
dephasingRossi2002  of the microscopic transition and acts as a dephasing:

$\displaystyle\frac{d}{dt}p^{\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}|_{\text{coll}}$
$\displaystyle=-\gamma_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}p^{\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}+\mathcal{U}^{\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}\;.$
(22)

The diagonal part is determined by the time- and momentum-dependent Coulomb
scattering rates

$\displaystyle\gamma_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}(t)$
$\displaystyle=\Gamma^{\text{in}}_{\lambda_{1},\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}+\Gamma^{\text{out}}_{\lambda_{1},\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}+\Gamma^{\text{in}}_{\lambda_{2},\mathbf{k}_{\parallel}}+\Gamma^{\text{out}}_{\lambda_{2},\mathbf{k}_{\parallel}}.$
(23)

The off-diagonal contribution couples to all coherences in the Brillouin zone
and reads

$\displaystyle\mathcal{U}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}^{\lambda_{1}\lambda_{2}}(t)$
$\displaystyle=\sum_{\mathbf{q},\mathbf{p}}\left(\mathcal{V}^{\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel},\mathbf{q},\mathbf{p}}p^{\lambda_{1}\lambda_{2}}_{\mathbf{q},\mathbf{p}}+\text{c.c}\right).$
(24)

For the carrier relaxation processes the Coulomb interaction is treated up to
second order Born-Markov approximation. The scattering rates in Eq. (23) and
$\mathcal{V}^{\lambda_{1}\lambda_{2}}_{\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{q},\mathbf{p}}$
are specified in the SM Sec. VII. The efficiency of the scattering channels is
determined by the Coulomb matrix element and the occupation probabilities of
the involved states. The scattering rates feature a sum over different band
indices, which include all possible Meitner-Auger relaxation channels, which
fulfill momentum and energy conservation at the same time Winzer2010 ;
Malic2011 . The Meitner-Auger effect Meitner1922 ; Auger1923 ; Grant2004 ;
Matsakis2019  as a non-radiative relaxation mechanism of the core-hole is
characterized by the filling of the inner-shell vacancy accompanied by the
emission of an electron into the unoccupied conduction band states and
possibly even out of the sample into the ionization continuum.

Finally, we present the dynamics in a many band system of the occupations in
band $\lambda_{1}$:

$\displaystyle\frac{d}{dt}f^{\lambda_{1}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}$
$\displaystyle=-\frac{1}{i\hbar}\sum_{\lambda,\mathbf{Q}_{\parallel}^{\prime},\mathbf{G}_{\parallel}}\mathbf{E}_{-\mathbf{Q}_{\parallel}^{\prime}+\mathbf{G}_{\parallel}}(t)\cdot\left(\mathbf{X}^{\lambda_{1}\lambda}_{\mathbf{k}_{\parallel},\mathbf{k}_{\parallel}-\mathbf{Q}_{\parallel}^{\prime}}(\mathbf{G}_{\parallel})p^{\lambda_{1}\lambda}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}-\mathbf{Q}_{\parallel}^{\prime}}-\mathbf{X}^{\lambda\lambda_{1}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}+\mathbf{Q}_{\parallel}^{\prime},\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}(\mathbf{G}_{\parallel})p^{\lambda\lambda_{1}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}+\mathbf{Q}_{\parallel}^{\prime},\mathbf{k}_{\parallel}}\right)$
$\displaystyle-\frac{2}{\hbar}\sum_{\lambda,\mathbf{k}_{\parallel}^{\prime},\mathbf{q}_{\parallel},\mathbf{p}_{\parallel}}\Im\textnormal{m}\left(\left[V^{\lambda_{1}\lambda\lambda\lambda_{1}}_{\mathbf{q}_{\parallel},\mathbf{k}_{\parallel}^{\prime},\mathbf{p}_{\parallel},\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}p^{\lambda\lambda_{1}}_{\mathbf{k}_{\parallel}^{\prime},\mathbf{k}_{\parallel}}-V^{\lambda_{1}\lambda_{1}\lambda\lambda}_{\mathbf{k}_{\parallel},\mathbf{q}_{\parallel},\mathbf{k}_{\parallel}^{\prime},\mathbf{p}_{\parallel}}p^{\lambda_{1}\lambda}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}^{\prime}}\right]p^{\lambda_{1}\lambda}_{\mathbf{q}_{\parallel},\mathbf{p}_{\parallel}}\right)$
$\displaystyle+\sum_{\begin{subarray}{c}\lambda,\mu,\nu\\\
\mathbf{k}_{\parallel}^{\prime},\mathbf{q}_{\parallel},\mathbf{p}_{\parallel}\end{subarray}}\left[V^{\lambda\mu\nu\lambda_{1}}_{\mathbf{q}_{\parallel},\mathbf{k}_{\parallel}^{\prime},\mathbf{p}_{\parallel},\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}\sigma^{\lambda\nu}_{\mathbf{q}_{\parallel},\mathbf{p}_{\parallel}}\sigma^{\mu\lambda_{2}}_{\mathbf{k}_{\parallel}^{\prime},\mathbf{k}_{\parallel}}-V^{\lambda_{1}\lambda\mu\nu}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{q}_{\parallel},\mathbf{k}_{\parallel}^{\prime},\mathbf{p}_{\parallel}}\sigma^{\lambda\nu}_{\mathbf{q}_{\parallel},\mathbf{p}_{\parallel}}\sigma^{\lambda_{1}\mu}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}^{\prime}}\right]$
$\displaystyle+\sum_{\mathbf{k}_{\parallel}^{\prime}}\left[\Gamma^{\text{in}}_{\lambda_{1},\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}^{\prime}}\left(\delta_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}^{\prime}}-f^{\lambda_{1}}_{\mathbf{k}_{\parallel}^{\prime},\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}\right)+\Gamma^{\text{in}}_{\lambda_{1},\mathbf{k}_{\parallel},\mathbf{k}_{\parallel}^{\prime}}\left(\delta_{\mathbf{k}_{\parallel},\mathbf{k}_{\parallel}^{\prime}}-f^{\lambda_{1}}_{\mathbf{k}_{\parallel},\mathbf{k}_{\parallel}^{\prime}}\right)\right]$
$\displaystyle-\sum_{\mathbf{k}_{\parallel}^{\prime}}\left[\Gamma^{\text{out}}_{\lambda_{1},\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}^{\prime}}f^{\lambda_{1}}_{\mathbf{k}_{\parallel}^{\prime},\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}+\Gamma^{\text{out}}_{\lambda_{1},\mathbf{k}_{\parallel},\mathbf{k}_{\parallel}^{\prime}}f^{\lambda_{1}}_{\mathbf{k}_{\parallel},\mathbf{k}_{\parallel}^{\prime}}\right]\;.$
(25)

The first line describes the excitation of a non-equilibrium electron
distribution in the band $\lambda_{1}$. The second line describes the
nonlinear Coulomb sources of the carrier population. The third line includes
Meitner-Auger-type terms. Therefore, the band indices in the Coulomb element
need to correspond to three occupied (unoccupied) and one unoccupied
(occupied) bands. Then line three carries a product of occupation and
interband transition. Note that, again the momentum conservation for the
Coulomb interaction can be inserted. The last two lines expresse the Coulomb
interaction described by a microscopic Boltzmann-like scattering equation. The
scattering rates explicitly include Pauli-blocking terms and are explicitly
given in the SM.

Investigating Eq. (21) we see that all quantities, except for the electric
field, are known. To determine the electric field as an observable self-
consistently, we have to solve the wave equation, which is performed in the
following section.

##  IV Observables and solution of the wave equation

Figure 4: Atomically thin material at $z=z_{0}$ under X-ray radiation
$\mathbf{E}^{0}_{\bm{Q}_{\parallel}}(z,\omega)$ and two-dimensional response
$P^{2D}_{\bm{Q}_{\parallel}}$ of the material. An incident angle of
$\theta=0^{\circ}$ denotes a perpendicular irradiation, while
$\theta=90^{\circ}$ corresponds to an in-plane propagating field.
$\mathbf{E}^{R}_{\bm{Q}_{\parallel}}(z,\omega)$ and
$\mathbf{E}^{T}_{\bm{Q}_{\parallel}}(z,\omega)$ denote the reflected and
transmitted X-ray field, respectively.

Figure 4 sketches the investigated geometry. A freestanding atomically thin
material is irradiated by an external X-ray field
$\mathbf{E}^{0}_{\bm{Q}_{\parallel}}(z,\omega)$ under an angle of incidence
$\theta$. We assume an infinitely thin two-dimensional material located at
$z_{0}$ and use a two-dimensional polarization
$\mathbf{P}^{2D}(\mathbf{r}_{\parallel},t)$ to describe the material response
$\mathbf{P}(\mathbf{r},t)=\mathbf{P}^{2D}({\mathbf{r}_{\parallel}},t)\delta(z-z_{0})$.
The measurable observables are the reflected
$\mathbf{E}^{R}_{\bm{Q}_{\parallel}}(z,\omega)$ and transmitted
$\mathbf{E}^{T}_{\bm{Q}_{\parallel}}(z,\omega)$ fields containing information
on the X-ray material response.

To determine the observables, we solve the Maxwell equations. In the absence
of free charges or currents and a spatially homogeneous dielectric environment
the corresponding wave equation reads

$\displaystyle\nabla^{2}\mathbf{E}(\mathbf{r},t)-\frac{\epsilon}{c^{2}}\frac{\partial^{2}}{\partial
t^{2}}\mathbf{E}(\mathbf{r},t)$
$\displaystyle=\mu_{0}\frac{\partial^{2}}{\partial
t^{2}}\mathbf{P}(\mathbf{r},t)$
$\displaystyle-\frac{1}{\epsilon_{0}\epsilon}\nabla\left(\nabla\cdot\mathbf{P}(\mathbf{r},t)\right)\;.$
(26)

Here, for simplicity, we used a uniform background with constant permittivity
$\epsilon$ representing off resonant electronic transitions and possible
substrate effects. For a free standing layer the dielectric constant reads
$\epsilon=1$. Clearly, $\mathbf{P}^{2D}({\mathbf{r}_{\parallel}},t)$ acts as a
link between the microscopic X-ray Bloch (cf. Eq. (19)) and the wave equations
Eq. (26). To study experimental observables, we need to solve the Bloch
equation Eq. (21) and wave equation Eq. (26) self-consistently.

A formal solution of the wave equation Eq. (26) is given by the Green’s
function. Transforming into momentum and frequency space
$\mathbf{E}_{\bm{Q}_{\parallel}}(z,\omega)=\int
dz^{\prime}~{}\mathcal{G}_{\bm{Q}_{\parallel}}(z-z^{\prime},\omega)\bm{\Pi}_{\bm{Q}_{\parallel}}(z^{\prime},\omega)$
yields an expression for the in-plane and out-of-plane electric field
components as a function of the in-plane radiation wave vector
$\bm{Q}_{\parallel}$, spatial $z$-component and frequency $\omega$ Mills1975 ;
Knorr1999 . $\bm{\Pi}_{\bm{Q}_{\parallel}}(z,\omega)$ denotes the Fourier
transform of the source acting on the right hand side in Eq. (26). The
corresponding Green’s function can be determined to be
$\mathcal{G}_{\bm{Q}_{\parallel}}(z-z^{\prime},\omega)=i\exp(-i\kappa|z-z^{\prime}|)/2\kappa$
with $\kappa^{2}=\epsilon\omega^{2}/c^{2}-|\bm{Q}_{\parallel}|^{2}$. The
electric field, separated in in-plane and out-of-plane components and X-ray
wave vector restricted to the first Brillouin zone, reads

$\displaystyle\mathbf{E}^{\parallel}_{\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel}}(z,\omega)$
$\displaystyle=\frac{i}{2\kappa}\left(\frac{(\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel})\otimes(\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel})}{\epsilon\epsilon_{0}}-\mu_{0}\omega^{2}\right)\mathbf{P}^{2D}_{\parallel\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel}}(\omega)e^{-i\kappa|z-z_{0}|}$
$\displaystyle-\frac{i(\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel})}{2\epsilon\epsilon_{0}}\text{sgn}(z-z_{0})e^{-i\kappa|z-z_{0}|}P^{2D}_{\perp\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel}}(\omega)$
(27) $\displaystyle
E^{\perp}_{\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel}}(z,\omega)$
$\displaystyle=\frac{i}{2}\left(\frac{\kappa}{\epsilon\epsilon_{0}}-\frac{\mu_{0}\omega^{2}}{\kappa}\right)e^{-i\kappa|z-z_{0}|}P^{2D}_{\perp\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel}}(\omega)-\frac{i}{2\epsilon\epsilon_{0}}\text{sgn}(z-z_{0})e^{-i\kappa|z-z_{0}|}(\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel})\cdot\mathbf{P}^{2D}_{\parallel\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel}}(\omega).$
(28)

with the dyadic product of the X-ray wave vectors
$(\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel})\otimes(\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel})$.
To obtain the full field acting in the Bloch equations Eq. (21), the incident
electric field
$\mathbf{E}^{0}_{\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel}}(z,\omega)$ has
to be added to Eq. (27) and Eq. (28) as homogeneous solution. In the limit of
$\mathbf{Q}_{\parallel}\rightarrow 0$, we find the result of a perpendicular
incidence of an optical field on a material sheet, where the field travels as
a plane wave varying only in the propagation direction Graphenbuch

To describe reflection and transmission, we are interested in the left
propagating field in front of the material and the propagating field behind
the sample, respectively. Investigating Fig. 4, for the reflected light for
all $z$ holds $z<z_{0}$. Therefore, the sign function is negative and the
phase factor reads $\exp(i\kappa(z-z_{0}))$. The reflection coefficient is
defined as the ratio of the intensity of the left propagating electric field
in front of the material and incident electric field intensity:
$r_{\mathbf{Q}_{\parallel}}(\omega)=|\mathbf{E}^{L}_{\mathbf{Q}_{\parallel}}(z<z_{0},\omega)|^{2}/|\mathbf{E}_{\mathbf{Q}_{\parallel}}^{0}(z,\omega)|^{2}$.
The transmitted field is located behind the sample, that the sign function
becomes positive and propagates in the direction $z>0$, that the phase factor
is $\exp(-i\kappa(z-z_{0}))$. The transmission coefficient is defined as
$t_{\mathbf{Q}_{\parallel}}(\omega)=|\mathbf{E}^{R}_{\mathbf{Q}_{\parallel}}(z>z_{0},\omega)|^{2}/|\mathbf{E}_{\mathbf{Q}_{\parallel}}^{0}(z,\omega)|^{2}$
using the right propagating electric field behind the sample Knorr1996 ;
Katsch2020 .

Note that, we exploited the two-dimensional nature of the material to access
an analytical solution of Maxwell’s equation for the propagating X-rays.
However, an expansion to three-dimensional bulk materials is possible: In
order to expand the theory to three dimensions, the observable Eq. (19) needs
to be defined in three dimensions including the crystal symmetry, i.e.
$\mathbf{P}_{\mathbf{Q}}(t)$. Subsequently, to define observables the Maxwell
equations need to be solved by applying additional boundary conditions
Ginzburgbuch ; Pollard2009 ; LaBalle2019 . Here, X-ray propagation effects
arise due to the finite sample thickness and the interface needs to be treated
carefully, tractable at least numerically. The occurring propagation effects
come from interactions between electronic resonances mediated by reemitted
photons in an optically thick sample. Propagation induced effects in optically
thick samples could be a pulse breakup or polariton beating in case of strong
light-matter coupling Aaviksoo1991 ; Frohlich1991 ; Mishina1993 . When the
electronic excitations decay radiatively the reemitted photon can propagate
giving rise to partial reflections Manzke1988 ; Stroucken1995 .

###  IV.1 Linear X-ray absorption

We rewrite the electric field as
$\mathbf{E}_{\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel}}(z,\omega)=[K_{\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel}}](z,\omega)\mathbf{P}^{2D}_{\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel}}(\omega)+\mathbf{E}^{0}_{\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel}}(z,\omega)$
with the matrix

$\displaystyle[K_{\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel}}](z,\omega)$
$\displaystyle=\frac{e^{-i\kappa|z-z_{0}|}}{2i}\begin{pmatrix}\frac{\mu_{0}\omega^{2}}{\kappa}-\frac{(Q_{x}+G_{x})^{2}}{\epsilon_{0}\epsilon\kappa}&-\frac{(Q_{x}+G_{x})(Q_{y}+G_{y})}{\epsilon_{0}\epsilon\kappa}&\frac{Q_{x}+G_{x}}{\epsilon_{0}\epsilon}\text{sgn}(z-z_{0})\\\
-\frac{(Q_{y}+G_{y})(Q_{x}+G_{x})}{\epsilon_{0}\epsilon\kappa}&\frac{\mu_{0}\omega^{2}}{\kappa}-\frac{(Q_{y}+G_{y})^{2}}{\epsilon_{0}\epsilon\kappa}&\frac{Q_{y}+Gy}{\epsilon_{0}\epsilon}\text{sgn}(z-z_{0})\\\
\frac{Q_{x}+G_{x}}{\epsilon_{0}\epsilon}\text{sgn}(z-z_{0})&\frac{Q_{y}+G_{y}}{\epsilon_{0}\epsilon}\text{sgn}(z-z_{0})&\frac{\mu_{0}\omega^{2}}{\kappa}-\frac{\kappa}{\epsilon_{0}\epsilon}\end{pmatrix}\;.$
(29)

Together with the definition of the macroscopic polarization Eq. (19) and the
solution of the X-ray Bloch equations in frequency space yields the
identification of the susceptibility
$[\chi_{\mathbf{Q}_{\parallel}}](\omega)$. Here, in anisotropic media the
susceptibility is a second rank tensor since polarization and electric field
may not be necessarily collinear anymore. Using the electric field Eq. (27)
and (28) at the position $z_{0}$ and applying the definition for
$\mathbf{P}^{2D}_{\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel}}(\omega)$
yields a self-consistent description of the electric field of the form

$\displaystyle\mathbf{E}_{\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel}}(z,\omega)$
$\displaystyle=\epsilon_{0}[K_{\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel}}](z,\omega)[\chi_{\mathbf{Q}_{\parallel}}](\omega)\left[1\\!\\!1-\epsilon_{0}\sum_{\mathbf{G}_{\parallel}^{\prime}}[K_{\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel}^{\prime}}](z_{0},\omega)[\chi_{\mathbf{Q}_{\parallel}}](\omega)\right]^{-1}\sum_{\mathbf{G}_{\parallel}^{\prime\prime}}\mathbf{E}^{0}_{\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel}^{\prime\prime}}(z_{0},\omega)+\mathbf{E}^{0}_{\mathbf{Q}_{\parallel}}(z,\omega)\;.$
(30)

Equation (30) describes the formation of the reflected and transmitted field
as driven by the incident electric field
$\mathbf{E}^{0}_{\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel}}(z,\omega)$.
The occurring dielectric susceptibility is determined by the microscopic Bloch
equations, calculated in Sec. (III). For the reflected and transmitted light
the sign function and the phase factor in
$[K_{\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel}}]$ has to be selected as
discussed at the end of the last section. We stress that the relation between
X-ray wave vector and angle of incidence is
$Q_{\parallel}^{2}=\omega^{2}\sin^{2}\theta/c^{2}$ obtained from the linear
light dispersion.

A complementary approach to include propagation effects is by inserting the
field Eq. (27) and (28) in the Bloch equations Eq. (21). Here, besides Coulomb
interaction as recombination channel included in the X-ray Bloch equations, a
second mechanism of dephasing is caused by radiative interaction: Conduction
band electron and the core-hole recombine under the emission of an X-ray
photon. This process is included in the X-ray Bloch equation Eq. (21) via the
self-consistently determined electric field. This allows for a self-consistent
description of the radiative dephasing in two-dimensional materials as a
function of the wave vector. Inserting the emitted electric field Eq. (27) and
Eq. (28) into Eq. (21) yields for the radiative contribution (for simplicity
here given without the incident field):

$\displaystyle\hbar\omega
p^{\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}(\omega)\big{|}_{\text{rad}}$
$\displaystyle=\frac{-i}{2\epsilon\epsilon_{0}}\sum_{\mathbf{k}_{\parallel}^{\prime},\mathbf{G}_{\parallel}}\left[\mathbf{d}^{\parallel\lambda_{2}\lambda_{1}}_{\mathbf{k}_{\parallel},\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}(\mathbf{G}_{\parallel})\left(\frac{1}{\kappa}\left((\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel})\otimes(\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel})-\epsilon\frac{\omega^{2}}{c^{2}}\right)\mathbf{d}^{\parallel\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel}^{\prime}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}^{\prime}}(\mathbf{G}_{\parallel})-(\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel})\right.\right.$
$\displaystyle\left.\left.\times
d^{\perp\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel}^{\prime}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}^{\prime}}(\mathbf{G}_{\parallel})\right)+d^{\perp\lambda_{2}\lambda_{1}}_{\mathbf{k}_{\parallel},\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}(\mathbf{G}_{\parallel})\left(\left(\kappa-\frac{\epsilon\omega^{2}}{\kappa
c^{2}}\right)d^{\perp}_{\mathbf{k}_{\parallel}^{\prime}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}^{\prime}}(\mathbf{G}_{\parallel})-(\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel})\cdot\mathbf{d}^{\parallel\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel}^{\prime}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}^{\prime}}(\mathbf{G}_{\parallel})\right)\right]p^{\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel}^{\prime}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}^{\prime}}(\omega)\;.$
(31)

The radiative interaction couples all X-ray induced transitions
$p^{\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}$
to all others. This way, the Bloch equation fully include the self-consistent
light-matter interaction. To gain more insights in Eq. (31), we determine the
main contribution to the dephasing from the diagonal dephasing
$\mathbf{k}_{\parallel}^{\prime}=\mathbf{k}_{\parallel}$. The diagonal
radiative dephasing reads

$\displaystyle\gamma^{\text{rad}}_{\mathbf{k}_{\parallel},\mathbf{Q}_{\parallel}}$
$\displaystyle=\frac{1}{2\epsilon\epsilon_{0}}\sum_{\mathbf{G}_{\parallel}}\left[-\epsilon\frac{\omega^{2}}{c^{2}}|\mathbf{d}^{\parallel\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}(\mathbf{G}_{\parallel})|^{2}\right.$
$\displaystyle\left.+\frac{1}{\kappa}\mathbf{d}^{\parallel\lambda_{2}\lambda_{1}}_{\mathbf{k}_{\parallel},\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}(\mathbf{G}_{\parallel})(\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel})\otimes(\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel})\mathbf{d}^{\parallel\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}(\mathbf{G}_{\parallel})\right.$
$\displaystyle\left.+(\kappa-\frac{\epsilon\omega^{2}}{\kappa
c^{2}})|d^{\perp\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}(\mathbf{G}_{\parallel})|^{2}\right.$
$\displaystyle\left.-2\Re\textnormal{e}((\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel})\cdot\mathbf{d}^{\parallel\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}(\mathbf{G}_{\parallel})d^{\perp\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel},\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}(\mathbf{G}_{\parallel}))\right]\;.$
(32)

For a vanishing $\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel}$ (normal
incidence) only the first term of Eq. (32) survives. In accordance to a first
order perturbation theory, Eq. (32) resembles Fermi’s golden rule for the
radiative broadening and energy shifts, both, determined by the square of the
dipole matrix element.

To compare our result with the literature, the rate of spontaneous emission
$W_{X}$ between an initial band $\lambda$ and the 1$s$-band for atoms is
calculated in dipole approximation by Fermi’s golden rule HendersonBuch :

$\displaystyle W_{X}\propto\omega^{3}\sum_{\beta,j}|\langle
1s,\beta,1s,0\mid\mathbf{r}\mid\lambda,\beta,j,0\rangle|^{2}$ (33)

where $\omega$ stands for the transition energy between initial and final
state. While we find for the two-dimensional semiconductor a
$\omega$-dependence of the dephasing, cf. Eq. (32) with inserted
$\mathbf{Q}_{\parallel}$, for atoms a $\omega^{3}$-dependence is well known.
The difference originates from the different dimensionality of the emitting
systems. To connect our results to the existing literature we specify to the
atomic case: Here, we can argue that the electronic wave functions have
strongly contributing values only for $|\mathbf{r}|<a_{B}/Z$ and we can
roughly approximate the transition integral to be proportional to $Z^{-1}$.
From Moseley’s law Moseley1913 , being an extension of the Rydberg formula, we
know that the transition energy $\omega$ is proportional to $Z^{2}$.
Consequently, we see that the spontaneous K-shell emission rate is
proportional to $Z^{4}$ KuzmanyBuch . It can be shown that for atoms the
Meitner-Auger yield $W_{A}$ is almost independent of the effective nuclear
charge FeldmannBuch . A semi-empirical expression for the X-ray yield is
introduced in Ref. Fink1966 :

$\displaystyle w_{X}=\frac{W_{X}}{W_{X}+W_{A}}\;,$ (34)

which is proportional to $Z^{4}$. Hence, for atomic systems Meitner-Auger
transitions offer higher sensitivity to core level transitions for low $Z$ and
luminescence is more likely to occur with increasing nuclear charge. The
advantage of the formalism presented here is that for radiative and Meitner-
Auger recombination, it might be possible to extent and establish similar
statements for crystalline solids. For example, the luminescence yield can be
calculated from the time-integrated conduction band electron occupation, which
decays radiatively
$W_{X}=\sum_{\mathbf{k}_{\parallel},\mathbf{Q}_{\parallel}}\gamma^{\text{rad}}_{\mathbf{k}_{\parallel},\mathbf{Q}_{\parallel}}f^{c}_{\mathbf{k}_{\parallel}}$.
In a similar way also the amount of non-radiatively decaying electrons can be
calculated from the microscopic scattering rates expressed by the X-ray Bloch
equations.

##  V Application to graphene

We apply the derived theory to the exemplary material graphene Graphenbuch ;
Reichbuch . Graphene is a monolayer of carbon atoms in a two-dimensional
honeycomb lattice. Carbon exhibits six electrons with an electronic ground
state configuration of 1$s^{2}$2$s^{2}$2$p^{2}$. The 1$s$-electrons form the
core bands, the 2$s$ and 2$p$ electrons are valence electrons and dominate the
energy dispersion around the Fermi level, which is set to 0 eV for undoped
graphene. The 2$s$ and 2$p_{x/y}$-electrons undergo a sp2-hybridization with
one further electron in the 2$p_{z}$-orbital. The in-plane hybridized
electrons form the $\sigma$-bands, which are responsible for the in-plane
covalent bonding between the atoms. The 2$p_{z}$-electrons form the
$\pi$-bands Dresselhausbuch .

Figure 5(a) illustrates the unit cell, spanned by the lattice vectors
$\mathbf{a}_{1}$ and $\mathbf{a}_{2}$, containing two atoms and constituting
the equivalent sublattices A and B, which are rotated with respect to each
other by $\pi/3$. The vectors $\mathbf{R}_{A\parallel}$ and
$\mathbf{R}_{B\parallel}$ denote the position of the corresponding atoms
obtained by linear combination of the lattice vectors. In nearest-neighbor
approximation the constructed vectors $\bm{\delta}_{BA}=\mathbf{R}_{B\parallel
i}-\mathbf{R}_{A\parallel}\equiv\bm{\delta}_{i}$ connect the atoms on
sublattice A with the three surrounding atoms on sublattice B, cf. Fig 5(a).
The six second nearest next-neighbors are located on the same sublattice.

Figure 5: (a) Hexagonal lattice of graphene. Empty circles denote the atoms on
sublattice A and filled circles on sublattice B. The vectors $\bm{\delta}_{i}$
denote the nearest-neighbor vectors connecting atom A to the three surrounding
B atoms. (b) Electronic dispersion above the Fermi level (grey dashed) and for
orientation also the filled $\pi$-band below 0 eV. The light green curves show
the $\pi$-bands and the blue curves show the $\sigma$-bands. Dark green
denotes the degenerate core band and violet the ionization continuum for
$k_{z}=0$.

###  V.1 Band structure and dipole moments of graphene

The energy dispersion and corresponding tight binding coefficients are
obtained from the Schrödinger equation Eq. (4) Dresselhausbuch , cf. Fig. 5
(b) which also includes the corresponding XANES and EXAFS transitions. The
Hamiltonian of graphene is symmetric with respect to the $x$-$y$ plane.
Consequently, the $\pi$ and $\sigma$-bands decouple since the former have an
odd symmetry and the latter even with respect to the reflection. For the
$\pi$-bands we obtain analytically Reich2002

$\displaystyle E_{\pi,\mathbf{k}_{\parallel}}=\pm
t_{\pi}|\xi_{\mathbf{k}_{\parallel}}|$ (35)

with the positive solution describing the valence band and the negative
solution the conduction band formed of 2$p_{z}$-electrons ($t_{\pi}<0$). The
shape of the band structure depends on the lattice symmetry only, described by
the nearest-neighbor form factor
$\xi_{\mathbf{k}_{\parallel}}=\sum_{i}\exp(i\mathbf{k}_{\parallel}\cdot\bm{\delta}_{i})$.
Here, for the analytical treatment we neglected contributions from the overlap
of electronic orbitals, which is generally small and lead to an asymmetry
between the valence and conduction $\pi$-band. Since we include the full form
factor $\xi_{\mathbf{k}_{\parallel}}$ without Taylor expansion trigonal
warping is included Rakyta2010 ; Rose2013 . The electronic hopping between
2$p_{z}$-orbitals from different lattice sites is denoted by $t_{\pi}$. The
band exhibits Dirac cone-like minima at the K-points, saddle points at the
M-points and a maximum at the $\Gamma$-point, all included in Eq. (35). As
discussed in Sec. II.1, the core dispersion is treated as a flat band
$E^{i}=-283\,\mathrm{eV}$ Susi2015  and the ionization continuum as a manifold
of three-dimensional parabolas. Because of the neglected hopping for the core
states between different sublattices the bonding and anti-bonding core state
are degenerate. The tight binding coefficients for $\pi$ and $1s$-state read
$C_{A2p_{z}/1s,\mathbf{k}_{\parallel}}^{\pi/1s}=\mp\xi_{\mathbf{k}_{\parallel}}/\sqrt{2}|\xi_{\mathbf{k}_{\parallel}}|$
and $C^{\pi/1s}_{B2p_{z}/1s,\mathbf{k}_{\parallel}}=1/\sqrt{2}$. In
theoretical approaches focusing on the response with respect to optical
frequencies in the visible range the orbital composition of graphene is
restricted to the 2$p_{z}$-orbital, which governs the optical properties of
graphene. The reason for this approximation that only the K and M-points of
the $\pi$-bands lie energetically in the optical range of the Fermi surface.
However, in X-ray experiments the excitation energy is tuned over several
hundreds of eV. Therefore, also states energetically further away from the
Fermi level have to be taken into account. The $\sigma$-bands are obtained by
numerical diagonalization of the Schrödinger equation. A spaghetti plot of the
used band structure is displayed in Fig. 5(b). Note, that for the numerical
evaluation we focus on an analytical 10-band tight-binding model as proof of
principle and as illustration of the approach. However, to fully describe
linear X-ray absorption experiments over several hundreds of eV an accurate
band structure beyond five orbitals, as it appears in Eq. (19), from ab initio
electronic structure theory, as provided by density function theory (DFT)
calculations are necessary. In this paper we use the 5-orbital system as a toy
model only. In the following, we will always refer with $\pi$ and $\sigma$ to
the unoccupied conduction band states.

With the tight-binding coefficients
$C^{\lambda}_{j\beta,\mathbf{k}_{\parallel}}$ we can calculate the XANES
dipole matrix element Eq. (16) for the core-conduction band transitions. Since
we are interested in excited interband transitions the intraband contribution
is neglected. Restricting the next-neighbor sum to the first nearest-neighbors
we obtain in dipole approximation for the K-edge dipole matrix elements

$\displaystyle\mathbf{X}^{1s\lambda}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}$
$\displaystyle=-e_{0}\sum_{\alpha,\beta,i}C^{*1s}_{\beta
1s,\mathbf{k_{\parallel}+Q_{\parallel}}}C^{\lambda}_{\alpha
i,\mathbf{k}_{\parallel}}$ $\displaystyle\times\left(\langle
1s,\beta,1s,0\mid\mathbf{r}\mid\lambda,\alpha,i,0\rangle~{}\delta_{\alpha,\beta}\right.$
$\displaystyle\left.+e^{i\mathbf{k_{\parallel}}\cdot\bm{\delta}_{\beta\alpha}}\langle
1s,\beta,1s,0\mid\mathbf{r}\mid\lambda,\alpha,i,\bm{\delta}_{\beta\alpha}\rangle\right)\;.$
(36)

The dipole matrix element exhibits an onsite (first term) and an offsite
contribution (second term). The band $\lambda=\pi$ consists just of the
$i=2p_{z}$ orbital. The $\lambda=\sigma$-band is built by $i=2s,2p_{x},2p_{y}$
orbitals. The offsite transition integral between sublattices A and B can be
calculated analytically by transforming the integral to prolate spheroidal
coordinates helping to handle the two-center nature of the integral. The
calculation can be found in the SM Sec. VIII. Here, we briefly present the
result of the calulations: The $z$-polarized $\pi$-transition depends only on
the absolute value of the next-neighbor vector, which is
$\delta_{i}=0.14\,\mathrm{nm}$ for the nearest-neighbors. The corresponding
integral has a value of

$\displaystyle\langle
1s,A,1s,0\mid\mathbf{r}\mid\pi,B,2p_{z},0.14\rangle=\begin{pmatrix}0\\\ 0\\\
0.14\end{pmatrix}\text{pm}\;.$ (37)

In contrast the transition into the $\sigma$-bands are in-plane polarized
because of the mirror symmetry of the graphene plane. The in-plane transitions
are differently weighted for each next-neighbor vector depending on its
angular orientation $\vartheta_{i}$ to the $y$-axis, cf. Fig. 5(a). We find

$\displaystyle\langle
1s,A,1s,0\mid\mathbf{r}\mid\sigma,B,2p_{x},\bm{\delta}_{i}\rangle=$
$\displaystyle\begin{pmatrix}\cos^{2}\vartheta_{i}&\sin^{2}\vartheta_{i}&0\\\
\sin\vartheta_{i}\cos\vartheta_{i}&-\sin\vartheta_{i}\cos\vartheta_{i}&0\\\
0&0&1\end{pmatrix}\begin{pmatrix}-0.12\\\ 0.14\\\ 0\end{pmatrix}\text{pm}$
$\displaystyle\langle
1s,A,1s,0\mid\mathbf{r}\mid\sigma,B,2p_{y},\bm{\delta}_{i}\rangle=$
$\displaystyle\begin{pmatrix}\sin\vartheta_{i}\cos\vartheta_{i}&-\sin\vartheta_{i}\cos\vartheta_{i}&0\\\
\sin^{2}\vartheta_{i}&\cos^{2}\vartheta_{i}&0\\\
0&0&1\end{pmatrix}\begin{pmatrix}-0.12\\\ 0.14\\\ 0\end{pmatrix}\text{pm}\;.$
(38)

In the $\delta\rightarrow 0$ limit the prolate spheroidal coordinates evolve
into the spherical coordinates. The onsite transitions for the $2p_{x}$,
$2p_{y}$ and $2p_{z}$-integrals are purely $x$, $y$ and $z$-polarized,
respectively, and have the value of 4.00 pm. Finally, the 1$s$ to 2$s$
transition does not contribute for vanishing connection vector because of
there even symmetry and the uneven parity of the dipole vector. For the
offsite transition we obtain

$\displaystyle\langle
1s,A,1s,0\mid\mathbf{r}\mid\sigma,B,2s,\bm{\delta}_{i}\rangle=-0.5\begin{pmatrix}\cos\vartheta_{i}\\\
\sin\vartheta_{i}\\\ 0\end{pmatrix}\text{pm}\;.$ (39)

The result is negative because of the node of the $2s$-orbital. Obviously,
offsite transitions $(\delta\neq 0)$ are smaller than onsite transition
$(\delta=0)$ integrals because of the small overlap of the core orbital with
other orbitals. The transition integral of the EXAFS matrix element for
1$s$-electrons has already been calculated in Eq. (18).

Next, we evaluate the X-ray Bloch equations. For a weak X-ray excitation
density we assume that the occupations are unchanged by the optics
corresponding spatially homogeneous system
($f^{\lambda}_{\mathbf{k}_{1},\mathbf{k}_{2}}\rightarrow
f^{\lambda}_{\mathbf{k}_{1},\mathbf{k}_{2}}\delta_{\mathbf{k}_{1},\mathbf{k}_{2}}$)
and set the diagonal occupation in initial and final band to one and zero,
respectively. Then, in linear optics, the solution of the microscopic Bloch
equation Eq. (21) reads, with dipole approximated matrix element, in frequency
space

$\displaystyle
p^{\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}(\omega)$
$\displaystyle=\frac{-\sum_{\mathbf{G}_{\parallel}}\mathbf{X}^{\lambda_{2}\lambda_{1}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}\cdot\mathbf{E}_{\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel}}(\omega)}{\hbar\omega-E^{\lambda_{2}}_{\mathbf{k}_{\parallel}}+E^{\lambda_{1}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}+i\gamma}$
(40) $\displaystyle p^{\lambda
k_{z}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}(\omega)$
$\displaystyle=\frac{-\sum_{\mathbf{G}_{\parallel}}\mathbf{Y}^{k_{z}\lambda}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}\cdot\mathbf{E}_{\mathbf{Q}_{\parallel}+\mathbf{G}_{\parallel}}(\omega)}{\hbar\omega-E^{k_{z}}_{\mathbf{k}_{\parallel}}+E^{\lambda}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}+i\gamma}$
(41)

for the XANES and EXAFS transitions, respectively. The former equation
describes absorption between two-dimensional bands within the material (1$s$
to conduction band) and the latter yields the absorption of 1$s$ to vacuum
transitions. Coulomb-induced renormalizations of the unexcited graphene band
gap are assumed to be included in the single-particle flat band core binding
energy $E_{1s}$ Susi2015 .As a first attack, and as experiments suggest, core-
excitonic effects of graphene are neglected in Eq. (40). This allows to
illustrate the strength of including the Bloch theorem in crystalline solids
via a first analytical approach Eq. (40) and (41). Additionally, we introduced
a phenomenological dephasing constant $\gamma$ attributing to the finite
lifetime of the electronic transition and leading to a broadening of the
absorption line. The radiative lifetime for graphene is of a few picoseconds
obtained from Eq. (32). Different studies report a core-hole lifetime in
graphite and carbon-based molecules of around 10 fs Harada2004 ;
Schlachter2004 ; Mucke2015 . However, the formalism allows for a
straightforward calculation of the linewidth due to Coulomb interaction (cf.
Eq. (22)) or phonon scattering Selig2016 ; Christiansen2017 ; Bernardi2021 .
Together with the definition of the two-dimensional macroscopic polarization,
being proportional to the electric field, we can identify the dyadic
susceptibility

Figure 6: (a) XANES of graphene for different angles of incidence. The first
peak stems from absorption into the $\pi$-bands, while the three succeeding
peaks are transitions into the $\sigma$-bands. (b) The $\pi$-peak as extracted
from the numerics shows a $\sin^{2}\theta$-behavior. The deviation at higher
angles stems from the fact that we assumed an infinite thin graphene layer.

$\displaystyle[\chi_{\mathbf{Q}_{\parallel}}](\omega)$
$\displaystyle=\sum_{\lambda_{1},\lambda_{2},\mathbf{k}_{\parallel}}\frac{\mathbf{X}^{\lambda_{1}\lambda_{2}}_{\mathbf{k}_{\parallel},\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}\otimes\mathbf{X}^{\lambda_{2}\lambda_{1}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}}{\hbar\omega-E^{\lambda_{2}}_{\mathbf{k}_{\parallel}}+E^{\lambda_{1}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}+i\gamma}$
$\displaystyle+\sum_{\lambda,\mathbf{k}_{\parallel},k_{z}}\frac{\mathbf{Y}^{\lambda
k_{z}}_{\mathbf{k}_{\parallel},\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}\otimes\mathbf{Y}^{k_{z}\lambda}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}}{\hbar\omega-E^{k_{z}}_{\mathbf{k}_{\parallel}}+E^{\lambda}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}+i\gamma}$
(42)

with the dyadic product of the dipole matrix elements yielding a $3\times
3$-matrix. The XANES dipole matrix element
$\mathbf{X}^{\lambda_{2}\lambda_{1}}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}$
is determined by Eq. (36) and the EXAFS matrix element $\mathbf{Y}^{\lambda
k_{z}}_{\mathbf{k}_{\parallel},\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}$
by Eq. (18). The X-ray wave number and the excitation frequency are related by
the relation $Q^{2}_{\parallel}=\omega^{2}\sin\theta/c^{2}$, where $\theta$
denotes the incident angle of the X-ray radiation and $c$ the speed of light
in the surrounding medium (cp. Sec. (IV)). Consequently, we perform no further
approximation and include all $\mathbf{Q}_{\parallel}$ over the full range of
the absorption spectrum. With the susceptibility all quantities are known and
the response of the material to the weak excitation is investigated by
calculating the absorption. Therefore, we introduce with the dipole
approximated dipole matrix element the absorption coefficient as
$\alpha_{\mathbf{Q}_{\parallel}}(\omega)=1-r_{\mathbf{Q}_{\parallel}}(\omega)-t_{\mathbf{Q}_{\parallel}}(\omega)$
defined by reflection $r_{\mathbf{Q}_{\parallel}}$ and transmission
$t_{\mathbf{Q}_{\parallel}}$. The reflection is defined by the reflected
intensity of the left propagating field in front of the graphene sheet
$r_{\mathbf{Q}_{\parallel}}(\omega)=I^{L}_{\mathbf{Q}_{\parallel}}(\omega)/I^{0}_{\mathbf{Q}_{\parallel}}(\omega)$
and the transmission by the intensity of the right propagating electric field
behind the graphene sheet
$t_{\mathbf{Q}_{\parallel}}(\omega)=I^{R}_{\mathbf{Q}_{\parallel}}(\omega)/I^{0}_{\mathbf{Q}_{\parallel}}(\omega)$.

###  V.2 X-ray spectrum of graphene

Due to the non-isotropic character of graphene, in the susceptibility all
possible components of the dipole matrix element are coupled, which leads to a
rather difficult expression for the absorption. We start by investigating the
XANES spectrum, i.e. the susceptibility in Eq. (42) is restricted to the first
line and calculating the absorption with Eq.(30).

Figure 7: Comparison between the true absorption coefficient, Fermi’s golden
rule and experiment. Due to a weak X-ray-matter interaction both calculated
absorption curves are almost identical. They reproduce the first four
resonances of the experiment but miss the plateau at higher energies since
these electronic states are not included in the numerical evaluation. The
experimental curve is adjusted in height to the calculated absorption
spectrum.

Figure 6(a) displays the calculated XANES spectrum of graphene for different
angles of incidence $\theta$. The absorption starts at 283 eV corresponding to
the binding energy of the core electrons relative to the Fermi surface of
graphene. The first peak at 286 eV rises with increasing angle and stems from
the transition of 1$s$-electrons into the $\pi$-band close to the M-point,
which exhibits a van Hove singularity due to a saddle point. The following
three peaks at 293 eV, 298 eV, and 303 eV stem from transitions into the three
$\sigma$-bands and are decreasing with angle of incidence. All four peaks are
observed in experiments for graphite Rosenberg1986 ; Buades2018  and graphene
Pacile2008 ; Papagno2009 ; Xu2015  and matches with the here calculated
energetic positions and spectral structure. However, when comparing Fig. 6(a)
to experiments Rosenberg1986 ; Buades2018 ; Pacile2008 ; Papagno2009 ; Xu2015
a clear difference can be noticed (cp. green line in Fig. 7): After the second
peak a plateau in the absorption spectrum is observed. This can be traced back
to the point that we limited our approach to orbitals up to the 2$p$-orbitals
(10 band approximation). Including more orbitals Willis1974  leads to a
continuum of close lying bands and should form an almost continuous
absorption. At this point, we want to point out, that the strength of the
developed Bloch equation approach lies in the description of self-consistent
many-body interaction, as introduced in Sec. III, and the description of
nonlinear optics occurring in ultrafast pump-probe type experiments as
attosecond transient absorption. The use of the tight binding approach to
calculate the band structure as input (as done here for illustration) is not
always sufficient (illustrated later on in Fig. 7) but a connection of the
Maxwell-X-ray Bloch formalism developed here and more advanced DFT methods, as
$\Delta$SCF-DFT Liang2017 ; Klein2021 , for single-particle energies and
coupling elements should be applied to accurately predict X-ray-matter
interaction and nonlinear effects. In our first attack here, the tight binding
method is just used for illustration.

Figure 6(b) shows the peak height of the $\pi$-peak as a function of the angle
of incidence. We obtain that the absorption follows a
$\sin^{2}\theta$-behavior before deviating from this trend at about 50 ∘: This
polarization-dependent absorption can be explained by the symmetry of the
orbital composition of the final band as was already done in different
theoretical and experimental works Chowdhury2012 ; Rosenberg1986 ;
Schiessling2003 ; HemrajBenny2006 . Because of the out-of-plane character of
the 2$p_{z}$-orbital the dipole transition for the $\pi$-peak is
$z$-polarized, which greatly simplify the susceptibility tensor since only one
component remains nonzero. Consequently, we can find an analytical expression
for the absorption into the $\pi$-band. Without the loss of generality we
assume for the in-plane incidence angle zero degree, such that $Q_{y}=0$. Then
we obtain for small $\mathbf{Q}_{\parallel}$ for the absorption

$\displaystyle\alpha_{\mathbf{Q}_{\parallel}}(\omega)=\frac{\Delta_{\mathbf{Q}_{\parallel}}(\omega)\Im\textnormal{m}\left(\chi^{33}_{\mathbf{Q}_{\parallel}}(\omega)\right)-\frac{Q_{x}^{2}}{2\epsilon}|\chi^{33}_{\mathbf{Q}_{\parallel}}(\omega)|^{2}}{|1-\frac{i}{2}\Delta_{\mathbf{Q}_{\parallel}}(\omega)\chi^{33}_{\mathbf{Q}_{\parallel}}(\omega)|^{2}}$
(43)

with
$\Delta_{\mathbf{Q}_{\parallel}}(\omega)=\frac{\omega^{2}}{c^{2}\kappa}-\frac{\kappa}{\epsilon}$
and $\chi^{33}_{\mathbf{Q}_{\parallel}}(\omega)$ as the $zz$-entry of the
susceptibility tensor. The coefficient
$\Delta_{\mathbf{Q}_{\parallel}}(\omega)$ stems from the matrix
$[K_{\mathbf{Q}_{\parallel}}](\omega)$ following from the Maxwell equations.
In particular, the denominator arises from the self-consistent treatment of
Maxwell and Bloch equations and includes for instance the radiative coupling
in the sample. The in-plane light wave vector $|\mathbf{Q}_{\parallel}|^{2}$,
which is orthogonal to the X-ray polarization vector of the incident light, is
proportional to $\sin^{2}\theta$. When inserting the definition for
$\kappa=\sqrt{\epsilon\omega^{2}/c^{2}-|\mathbf{Q}_{\parallel}|^{2}}$ and
express the X-ray wave vector as function of the angle of incidence we obtain
for the prefactor
$\Delta_{\mathbf{Q}_{\parallel}}=\omega^{2}\sin^{2}\theta/\epsilon
c^{2}\sqrt{\epsilon-\sin^{2}\theta}$. Now, we explicitly observe the
$\sin^{2}$-dependence of the true absorption coefficient for the
$\pi$-transition. For a perpendicular irradiation of the sample
$\theta=0^{\circ}$ we see immediately that
$\Delta_{\mathbf{Q}_{\parallel}}(\omega)$ and $Q_{x}$ vanish and consequently
also the absorption vanishes as expected from the $z$-polarized transition.
The derived absorption formula diverges for an incident angle of
$\theta=90^{\circ}$ in the case of $\epsilon=1$, explaining the observed
deviation in Fig. 6(b). That for $\epsilon>1$ the singularity does not
explicitly appear anymore, suggests that the origin lies in the assumption of
an infinitely thin layer in vacuum. The divergence can be lifted by starting
with a three layer model – supstrate, layer, substrate – all with finite
thickness and solve the Maxwell equation for each region with corresponding
continuity conditions Chewbuch . In our case, because of the weak X-ray matter
interaction, the additional denominator in Eq. (43) including the radiative
coupling plays only a minor role. We elaborate on this in the following.

In Fig. 7, we plot the 20∘-result from Fig. 6(a) and a XANES measurement on
graphene taken from Ref. Papagno2009  (16∘). As discussed previously, we can
see that the four resonances, which are explicitly calculated, reproduces well
the experiment, but that the plateau is missed due to the fail of the
restricted tight binding method (cp. above) applied in our approach. To
evaluate the influence of X-ray propagation effects, included in the
denominator unequal to one in Eq. (43), we plot the result obtained from
Fermi’s golden rule, which neglects propagation effects. To obtain Fermi’s
golden rule, we start from the X-ray Bloch equation Eq. (40) and carry out the
limit $\gamma\rightarrow 0$. The intensity of a light field propagating though
a medium is damped by the imaginary part of the complex permittivity Kochbuch
. Therefore, we can approximate the absorption by
$\alpha(\omega)=\omega\Im\textnormal{m}(\chi(\omega))/cn=\omega\Im\textnormal{m}(P(\omega)/E(\omega))/\epsilon_{0}cn$
Kochbuch , where we projected the polarization density on the polarization
vector and introduced the material refractive index $n$. With the macroscopic
polarization density we obtain an expression for the absorption coefficient
for mall $\mathbf{Q}_{\parallel}$ in terms of Fermi’s golden rule

$\displaystyle\alpha_{\mathbf{Q}_{\parallel}}(\omega)$
$\displaystyle=\frac{\omega\pi}{\epsilon_{0}cn}\frac{A}{4\pi^{2}}\sum_{\lambda}\int
d^{2}k_{\parallel}~{}|X_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel},\mathbf{k}_{\parallel}}^{1s\lambda}|^{2}$
$\displaystyle\times\delta\left(E^{\lambda}_{\mathbf{k}_{\parallel}}-E^{1s}_{\mathbf{k}_{\parallel}+\mathbf{Q}_{\parallel}}-\hbar\omega\right)\;,$
(44)

which we can fairly compare to the self-consistent result. In Eq. (44) we
transformed the $\mathbf{k}_{\parallel}$-sum into an integral with area $A$.
To have a fair comparison with our result, the delta function is
phenomenologically broadened by the same rate $\gamma$ as the X-ray Bloch
equations. Figure 7 shows the result of Fermi’s golden rule Eq. (44), which
well reproduces the self-consistent result: We see that the magnitude of the
results are comparable. The reason lies in the weak X-ray-matter interaction.
Differences between Fermi’s golden rule and the true absorption can be
expected, when non-markovian effects such as quantum kinetics of electron-
phonon interactions (if not included in Fermi’s golden rule) dominate the
lineshape or strong radiative interaction occurs in a way that the
perturbation approach of Eq. (44) breaks down. The latter can for example
occur for two-dimensional materials involving more layers enhancing the light-
matter interaction: For example, taking ten two-dimensional sheets can lead to
an increase of the susceptibility by a factor of 10 Graphenbuch . Increasing
the susceptibility leads to a deviation of Fermi’s golden rule and the true
absorptionStroucken1996 . In this case, the true absorption becomes suppressed
and accumulates an increased radiative broadening, cf. Fig. S2 of the SM. Also
the oscillator strength can be distributed differently, due to lineshifts.
Other scenarios, where Fermi’s golden rule breaks down are X-ray excitations
of materials in a cavity or for materials, which have an intrinsically
stronger light-matter interaction. For the latter nonlinear many-body effects
can show up differently in transmission and reflection geometry Katsch2020_2
and reflection is not included in Fermi’s golden rule.

Figure 8: X-ray absorption spectrum of graphene for $30^{\circ}$ irradiation.
We obtain a dominating XANES contribution stemming from transition between
core and conduction bands. Also, the weaker transitions into the ionization
continuum, displayed in the inset, are obtained. The absorption into the
three-dimensional continuum (grey dashed displays the joint density of states)
is modulated by oscillations with maxima at different energies.

After investigating the XANES contribution, we compute the full X-ray
absorption spectrum by including also the EXAFS part, described by the second
line of Eq. (42). The ionization threshold is set artificially to the maximum
of the energetically highest lying $\sigma$-band. Figure 8 shows the full
absorption spectrum ($\theta=30^{\circ}$) and for better visibility the EXAFS
part as inset, where the absorption involving the three-dimensional electronic
continuum of vacuum states as final states can be observed. The absorption
curve is modulated with oscillations. For the particular situation plotted
here we recognize three maxima. Following the pioneering work by Sayers et al.
it is widely recognized that these wiggles could be used to obtain
quantitative information about the spatial structure of the atomic lattice the
X-rays are interacting with. The current interpretation in literature is that
these oscillations stem from a modulation of the absorption cross section due
to interference of the X-ray waves between neighboring atoms. In our approach,
which is specifically developed for crystalline solids, includes the full
lattice symmetry by the sum over reciprocal lattice vectors, we calculate the
absorption using the susceptibility determined by the second line of Eq. (42),
which naturally includes the Bloch theorem. For a constant EXAFS dipole matrix
element the second line of Eq. (42) summed over the electronic wave vectors
describes a square-root like absorption into the three-dimensional continuum
states of the vacuum electrons, cf. the grey dashed line of Fig. 8 (inset).
However, this simplified absorption line is modulated by the dyadic product of
the EXAFS matrix elements, cf. Eq. (42). The occurring form factor of the
EXAFS matrix element is abbreviated $f_{n}(\mathbf{k})=8\pi
i\zeta^{3/2}\mathbf{k}_{n}\hat{\mathbf{e}}_{n}/(\zeta^{2}+(\mathbf{k})^{2})^{3}$,
cf. Eq. (18). Analytically we obtain for the squared dipole matrix element Eq.
(18) of the entry $nm$ in the susceptibility tensor Eq. (42)

$\displaystyle|(C^{*1s}_{A1s,\mathbf{k}_{\parallel}}+C^{*1s}_{B1s,\mathbf{k}_{\parallel}})|^{2}f_{n}(\mathbf{k})f_{m}(\mathbf{k})=$
$\displaystyle\left(1+\frac{1}{|\xi_{\mathbf{k}_{\parallel}}|}\sum_{i}\cos(\mathbf{k}_{\parallel}\cdot\bm{\delta}_{i})\right)f_{n}(\mathbf{k})f_{m}(\mathbf{k})$
(45)

where we inserted the definitions for the tight binding coefficients, shown
below Eq. (35). Similar to the work of Sayers, in Eq. (45), bracket right hand
side, we clearly find an oscillating behavior in $\mathbf{k}_{\parallel}$ with
the lattice vector frequency $\bm{\delta}_{i}$ comparable to Eq. (2). Due to
the small nuclear number, the broadening of the form factor $f_{n}$ is that
strong that $f_{n}$ is almost constant, $f_{n}(\mathbf{k})\equiv f_{n}$.
Therefore, the reciprocal lattice geometry is not resolved and its influence
on $\alpha_{\mathbf{Q}_{\parallel}}(\omega)$ is negligible. Consequently, Eq.
(45) is determined solely by the squared tight binding coefficients
$(1+\sum_{i}\cos(\mathbf{k}_{\parallel}\cdot\bm{\delta}_{i})/|\xi_{\mathbf{k}_{\parallel}}|)$,
which are responsible for the cosine-like oscillations in Fig. 8. From a
physical point of view, the tight binding coefficients describe the orbital
contribution of a specific lattice site to the band composition, cf. Eq. (15).
Due to the lattice periodicity, the tight binding coefficients are a function
of the wave vector $\mathbf{k}_{\parallel}$, which is determined by the
lattice geometry. We can therefore understand the absorption, summing over all
wave vectors (Eq. (42)) and proportional to Eq. (45), as showing the quantum
interference between orbital electronic wave functions, which can be
constructive or destructive depending on the sign of the tight binding
coefficients.

Summarizing this discussion, we can conclude that the EXAFS oscillations stem
from an interference of standing electronic wave functions of the graphene
sublattices $A$ and $B$, which are summed up in Eq. (45). The interference is
a spatially stationary property of the intrinsic electronic wave function, and
therefore it is independent of the X-ray excitation.

Figure 9: (a) Fourier amplitude of the EXAFS spectrum. The blue curve includes
only the first nearest-neighbor (nn). It peaks at the first nearest-neighbor
distance and its higher harmonics. The green curve includes also the second
and third next-neighbors. The peaks stem from the first three neighbors
(black) and interferences among themselves (yellow). (b) Comparison of Bloch
theory, experiment (adapted from Ref. Buades2018 ) and point scattering theory
(Eq. (47)).

Since for EXAFS, the X-ray light has an excitation energy higher than the
threshold energy, the absorption coefficient
$\alpha_{\mathbf{Q}_{\parallel}}(\omega)$ can be expressed as function of the
wave number of free electrons, i.e. with dispersion $\omega(k)=\hbar
k^{2}/2m_{0}+E_{Ion}/\hbar$ as final states. To discuss the spatial
interference in close analogy to experiments, the absorption coefficient is
Fourier transformed $\int
d^{2}k~{}k^{m}\alpha_{\mathbf{Q}_{\parallel}(\mathbf{k}_{\parallel})}(\mathbf{k}_{\parallel})\exp(i\mathbf{k}_{\parallel}\cdot\mathbf{r}_{\parallel})$
to real space with respect to the wave number. To increase the visibility and
access more conveniently the oscillations at higher $k$, the absorption
spectrum as function of wave number is multiplied by $k^{m}$ with typically
$m\in[1,2,3]$ and Fourier transform subsequently Ertl2008 ; Kas2016 ;
Britz2020 . In the following we discuss the Fourier transform of the
absorption spectrum evaluated by our description, which involves the Bloch
theorem for solid states:

The blue curve in Fig. 9(a) displays the Fourier amplitude of the $k$ ($m=1$)
weighted and background corrected EXAFS spectrum Ravel2005  in first neighbor
approximation. The effective radial distribution function clearly displays a
peak at the distance of the 1${}^{\text{st}}$ neighbor. The second peak at
0.28 nm lying at twice of the 1${}^{\text{st}}$ nearest-neighbor distance
originates from the performed discrete Fourier transformation on a finite
grid.

So far we have considered only the nearest-neighbor hoppings for the
calculation of the matrix elements. As a result we obtain a peak at the first
neighbor distance. However, calculating the matrix elements beyond the
nearest-neighbor approximation should add additional features to the Fourier
amplitude of the EXAFS, cf. Fig. 9(a) green curve. To go beyond the nearest-
neighbor hopping we introduce the second
$\xi^{(2)}_{\mathbf{k}_{\parallel}}=\eta_{2}\sum_{i=4}^{9}\exp(i\mathbf{k}_{\parallel}\cdot\bm{\delta}_{i})$
and third
$\xi^{(3)}_{\mathbf{k}_{\parallel}}=\eta_{3}\sum_{i=10}^{12}\exp(i\mathbf{k}_{\parallel}\cdot\bm{\delta}_{i})$
next-neighbor form factors. For the free parameters of second and third next-
neighbor hoppings we reasonably choose $\eta_{2}=0.05$ and $\eta_{3}=0.01$,
respectively: The ratios of $\eta_{2}$ and $\eta_{3}$ with respect to the
hopping parameter to the nearest-neighbor are chosen such that they coincide
with the relative deviation of the hoppings between 2$p_{z}$-orbitals with
increasing order of neighbors Kundu2017 ; Tran2017 . The peak heights are
uncertain under this assumption, but the peak positions – which is the most
important in our study – are unaffected by that. The tight binding
coefficients up to the third neighbor hopping read

$\displaystyle
C^{1s}_{A1s,\mathbf{k}_{\parallel}}=\frac{\xi_{\mathbf{k}_{\parallel}}+\xi^{(3)}_{\mathbf{k}_{\parallel}}}{\sqrt{2}\left(|\xi_{\mathbf{k}_{\parallel}}+\xi^{(3)}_{\mathbf{k}_{\parallel}}|+\xi^{(2)}_{\mathbf{k}_{\parallel}}\right)},\quad\text{and}\quad
C^{1s}_{B1s}=\frac{1}{\sqrt{2}}$ (46)

and enter the matrix element, described by the first line of Eq. (45), now
including 12 neighbors instead of three.

The green line in Fig. 9(a) displays the effective radial distribution
function including hoppings up to the third neighbor coupling in the EXAFS
matrix element. In contrast to the calculation with first neighbor hopping
(blue curve), we find additional peaks at the second and third neighbor
distance of 0.24 nm and 0.28 nm, respectively. Interestingly, since the
modulus square of the dipole matrix element is observed, also peaks appear,
which do not correspond to next neighbor vectors, but do correspond to
interferences of different next neighbor vectors. Exemplary, we obtain a peak
at 0.38 nm corresponding to the sum of the first and second neighbor distance.
Since 0.38 nm also correspond to the $4^{\text{th}}$ neighbor hopping in
graphene, this peak could also be interpreted as the fourth neighbor in
experiments. At 0.05 nm we can resolve a peak, which arises from the
difference of third and second neighbor. So far, this peak has been explained
as phase shift stemming from a difference between measured and geometrical
interatomic distances and required a theoretical or experimental correction
Rehr2000 . In contrast, in our approach, we show that the peak can be
interpreted as quantum interference between electronic Bloch wave functions of
first and second neighbor. The slight sideband of the first nearest-neighbor
peak at 0.1 nm stems from the difference of the second and first neighbor.
Further, we observe a peak at 0.42 nm, which can be understood as interference
of the $1^{\text{st}}$ with $3^{\text{rd}}$ neighbor.

To manifest our interpretation, in Fig. 9(b) we provide a direct comparison of
our full computational result up to the third nearest-neighbor to experiment,
adapted from Ref. Buades2018 , and additionally display the outcome obtained
from the point scattering theory Sayers1971 . In point scattering theory the
oscillatory part of the EXAFS is solely described by the matrix element
squared corresponding to Eq. (2). The Fourier transform yields the structure-
related function Sayers1971

$\displaystyle
S(r)=\frac{1}{2}\sum_{i}\frac{N_{i}}{\sigma_{i}R_{i}^{2}}e^{-\gamma
R_{i}}e^{-2(r-R_{i})^{2}/\sigma_{i}^{2}}$ (47)

displayed as dashed line in Fig. 9(b). The $\gamma$-factor accounts for the
photoelectron scattering range and leads to an overall decrease with
increasing distance. The result is a sum of Gaussian functions lying at the
next-neighbor distances $R_{i}$.