# Real-space Green’s function approach for intrinsic losses in x-ray spectra

J. J. Kas  Dept. of Physics, Univ. of Washington Seattle, WA 98195  J. J. Rehr
Dept. of Physics, Univ. of Washington Seattle, WA 98195

###### Abstract

Intrinsic inelastic losses in x-ray spectra originate from excitations in an
interacting electron system due to a suddenly created core-hole. These losses
characterize the features observed in x-ray photoemission spectra (XPS), as
well as many-body effects such as satellites and edge-singularities in x-ray
absorption spectra (XAS). However, they are usually neglected in practical
calculations. As shown by Langreth these losses can be treated within linear
response in terms of a cumulant Green’s function in momentum space. Here we
present a complementary ab initio real-space Green’s function (RSGF)
generalization of the Langreth cumulant in terms of the dynamically screened
core-hole interaction $W_{c}(\omega)$ and the independent particle response
function. We find that the cumulant kernel $\beta(\omega)$ is analogous to
XAS, but with the transition operator replaced by the core-hole potential with
monopole selection rules. The behavior reflects the analytic structure of the
loss function, with peaks near the zeros of the dielectric function,
consistent with delocalized quasi-boson excitations. The approach simplifies
when $W_{c}(\omega)$ is localized and spherically symmetric. Illustrative
results and comparisons are presented for the electron gas, sodium, and some
early transition metal compounds.

Green’s function, cumulant, screening

###### pacs:

71.15.-m, 31.10.+z,71.10.-w

##  I Introduction

Intrinsic losses in x-ray spectra are fundamental to the photoabsorption
process.Rehr and Albers (1990) They originate from the dynamic response of the
system to a suddenly created core-hole, leading to dynamic screening by local
fields and inelastic losses. This transient response is responsible for
observable effects in x-ray photoemission spectra (XPS) and x-ray absorption
specta (XAS). These include satellites due to quasi-bosonic excitations such
as plasmons, charge-transfer, and shake-processes, as well as particle-hole
excitations responsible for edge-singularity effects. These features are
signatures of electronic correlation beyond the independent particle
approximation.Zhou et al. (2020) Various theoretical techniques have been
developed for treating these losses, including plasmon models, quasi-boson
approximations, fluctuation potentials, determinantal approaches, dynamical-
mean-field theories, configuration-interaction methods, and coupled-cluster
approaches.Langreth (1970); Bardyszewski and Hedin (1985); Hedin (1999);
Campbell et al. (2002); Fujikawa (2001); Lee et al. (1999); Hedin et al.
(1998); Klevak et al. (2014); Combescot and Nozieres (1971); Tzavala et al.
(2020); Liang et al. (2017); Biermann and van Roekeghem (2016); Bagus et al.
(2022) Recently cumulant Green’s function methods have been developedKas et
al. (2015) based on a real-space real-time (RSRT) generalization of the
Langreth cumulant.Langreth (1970) While the approach gives good results for
the satellites observed in XPS, even for moderately correlated systems such as
transition metal oxides,Kas et al. (2015); Woicik et al. (2020, 2020) it
depends on computationally demanding real-time time-dependent density
functional theory (TDDFT) calculations of the density-density response
function. Thus despite these advances, quantitative calculations remain
challenging, and intrinsic losses are usually neglected in current
calculations of x-ray spectra.

In an effort to facilitate these calculations, we present here an ab initio
real-space Green’s function (RSGF) generalization of the Langreth cumulant
complementary to the RSRT approach, in which the calculations are carried out
using a discrete site-radial coordinate basis. The formalism of the cumulant
kernel $\beta(\omega)$ is analogous to x-ray absorption spectra $\mu(\omega)$,
except that the transition operator is replaced with the core-hole potential
$V_{c}({\bf r})$, and the transitions are between valence and conduction
states with monopole selection rules. The generalized RSGF approach thereby
permits calculations of many-body effects in x-ray spectra in parallel with
RSGF calculations of XAS.Rehr and Albers (1990); Kas et al. (2020) Several
representations of $\beta(\omega)$ are derived, which are useful in the
analysis and comparison with other approximations. For example, we show that
$\beta(\omega)$ can be expressed either in terms of the bare core-hole
potential $V_{c}({\bf r})$ and the full density response function $\chi({\bf
r},{\bf r}^{\prime},\omega)$, or the dynamically screened core hole potential
$W_{c}({\bf r},\omega)$ and the independent particle response function
$\chi^{0}({\bf r},{\bf r}^{\prime},\omega)$. In adddition, we derive the link
between the Langreth cumulant, and the commonly used approximation based on
the $GW$ self-energy.Aryasetiawan et al. (1996); Guzzo et al. (2011); Hedin
(1999) The potential $W_{c}({\bf r},\omega)$ is a key quantity of interest in
this work. However, its real-space behavior does not appear to have been
extensively studied heretofore. This quantity is closely related to the
dynamically screened Coulomb interaction $W({\bf r},{\bf r}^{\prime},\omega)$
used e.g., in Hedin’s $GW$ approximation for the one-electron self-
energy.Hedin (1999) The RSGF approach simplifies when $W_{c}({\bf r},\omega)$
is well localized and spherically symmetric. This leads to a local model for
the cumulant kernel on a 1-$d$ radial grid. The local approach is tested with
calculations for the homogeneous electron gas (HEG), and illustrative results
are presented for nearly-free-electron systems and early $3$d transition metal
compounds. We find that the local model provides a good approximation for
$\beta(\omega)$. The model also accounts for the Anderson edge-singularity
exponent in metals. The behavior of the cumulant kernel reflects the analytic
structure of the loss function, with pronounced peaks near the zeros of the
dielectric function. This structure is consistent with interpretations of
intrinsic excitations in terms of plasmons or charge-transfer excitations.

The remainder of this paper is organized as follows. Sec. II. summarizes the
Langreth cumulant and the RSGF and RSRT generalizations. Sec. III. and IV.
respectively describe the calculation details and results for the HEG and
charge-transfer systems. Finally Sec. V. contains a summary and conclusions.

##  II Theory

###  II.1 Cumulant Green’s function and x-ray spectra

Intrinsic inelastic losses in x-ray spectra including particle-hole, plasmons,
shake-up, etc., are characterized by features in the core-hole spectral
function $A(\omega)=\Sigma_{n}|S_{n}|^{2}\delta(\omega-\varepsilon_{n})$,
where $S_{n}$ is the amplitude for an excitation of energy $\varepsilon_{n}$
due to the creation of the core-hole. Equivalently $A_{c}(\omega)$ is given by
the Fourier transform of the core-hole Green’s function, i.e., the one-
particle Green’s function with a deep core-hole $c$ created at $t=0$
$g_{c}(t)=\langle 0|c_{c}^{\dagger}e^{iHt}c_{c}|0\rangle\theta(t)$,

$A_{c}(\omega)=-\frac{1}{\pi}{\rm Im}\int dt\,e^{-i\omega t}g_{c}(t).$ (1)

Here $H$ is the Hamiltonian of the system while $c^{\dagger}_{c}$ and $c_{c}$
are creation and annihilation operators, respectively. The core-level XPS
photocurrent $J_{k}(\omega)\sim A_{c}(\omega)$ is directly related to the
spectral function, which describes both the asymmetry of the quasi-particle
peak and satellites in the spectra. Intrinsic losses in x-ray absorption
spectra (XAS) $\mu(\omega)$ and related spectra (e.g., EELS) from deep core-
levelsCampbell et al. (2002) can be expressed in terms of a convolution of
$A_{c}(\omega)$ and the single (or quasiparticle) XAS
$\mu_{1}(\omega)$Campbell et al. (2002)

$\mu(\omega)=\int
d\omega^{\prime}\,A_{c}(\omega^{\prime})\mu_{1}(\omega-\omega^{\prime}).$ (2)

This accounts for effects such as satellites and the reduction factor
$S_{0}^{2}$ in the XAS fine structure.Rehr and Albers (1990) Here and below,
unless otherwise noted, we use atomic units $e=\hbar=m=1$ with distances in
Bohr = 0.529 Å and energies in Hartrees = 27.2 eV. A cumulant Green’s
function, which is a pure exponential in the time-domain, is particularly
appropriate for the treatment of intrinsic losses,

$g_{c}(t)=g_{c}^{0}(t)e^{C_{c}(t)},$ (3)

where $g^{0}_{c}(t)=e^{i\epsilon_{c}t}$ is the independent particle Green’s
function for a given core-level $c$, and $C_{c}(t)$ is the cumulant. This
representation naturally separates single (or quasi-particle) and many-
particle aspects of the final-state of the system with a deep core hole
following photoabsorption, where many-body effects are embedded in the
cumulant. It is convenient in the interpretation to represent the cumulant in
Landau form,Landau (1944)

$C_{c}(t)=\int d\omega\,\beta(\omega)\frac{e^{i\omega t}-i\omega
t-1}{\omega^{2}}.$ (4)

where the cumulant kernel $\beta(\omega)$ characterizes the strength of the
excitations at a given excitation energy $\omega$. This representation yields
a normalized spectral function, with a quasi-particle renormalization constant
$Z=\exp(-a)$, where $a=\int d\omega\,\beta(\omega)/\omega^{2}$ is a
dimensionless measure of correlation strength, and $\Delta=\int
d\omega\beta(\omega)/\omega$ is the relaxation energy shift of the core-
level.Hedin (1999) Hedin (1999)

###  II.2 Langreth cumulant

For a deep core-hole coupled to the interacting electron gas Langreth showed
that within linear response, the intrinsic inelastic losses can be treated in
terms of a cumulant Green’s function, with a cumulant kernel in frequency and
momentum space given by

$\beta(\omega)=-\sum_{q}|V_{q}|^{2}S({\bf q},\omega).$ (5)

Here $V_{q}$ is the Fourier transform of the core-hole potential and $S({\bf
q},\omega)$ is the dynamic structure factor, which is directly related to
time-Fourier transform of the density-density correlation function $\chi({\bf
q},\omega)$ and the loss function $L({\bf q},\omega)=-{\rm
Im}\,\epsilon^{-1}({\bf q},\omega)$, i.e.,

$\displaystyle S({\bf q},\omega)$ $\displaystyle\equiv$
$\displaystyle-\frac{1}{\pi}{\rm Im}\,\chi({\bf q},\omega)\theta(\omega)$ (6)
$\displaystyle=$ $\displaystyle-\frac{1}{\pi v_{q}}{\rm
Im}\,\epsilon^{-1}({\bf q},\omega)\theta(\omega),$ $\displaystyle\chi({\bf
q},\omega)$ $\displaystyle=$ $\displaystyle\int dt\,e^{i\omega
t}\langle\rho_{\bf q}(t)\rho_{\bf q}(0)\rangle.$ (7)

The response function $\chi({\bf q},\omega)$ can be expressed in terms of the
non-interacting response $\chi^{0}({\bf q},\omega)$ using the relation
$\chi=\chi^{0}+\chi^{0}K\chi$. Here the particle-hole interaction kernel $K$
within TDDFT is given by $K=v+f_{xc}$, where $v=4\pi/q^{2}$ is the bare
Coulomb interaction, and $f_{xc}=\delta v_{xc}[n]/\delta n$ is the TDDFT
kernel; this is obtained from the exchange-correlation potential used in the
definition of the independent particle response function $\chi^{0}(\omega)$.

Although Langreth’s expression for $\beta(\omega)$ in Eq. (5) is elegant for
its simplicity, calculations of the full density response function $\chi({\bf
q},\omega)$ are challenging, being comparable to that for a particle-hole
Green’s function or the Bethe-Salpeter equation (BSE). Moreover, the core-hole
potential $V_{c}({\bf r})$ has a long range Coulomb tail to deal with. On the
other hand, the Thomas-Fermi screening length $r_{0}=0.64\,r_{s}^{1/2}$ is
short-ranged, so one may wonder to what extent a local approximation for the
dynamically screened interaction might be applicable?

To this end, we note that the cumulant kernel can be expressed equivalently in
terms of the dynamically screened core-hole interactions $W_{\bf
q}(\omega)=V_{q}/\epsilon({\bf q},\omega)$ and the independent particle
response function $\chi^{0}({\bf q},\omega)$,

$\displaystyle\beta(\omega)$ $\displaystyle=$
$\displaystyle\frac{1}{\pi}\sum_{q}|V_{q}|^{2}{\rm Im}\,\chi({\bf
q},\omega)\theta(\omega),$ $\displaystyle\equiv$
$\displaystyle\frac{1}{\pi}\sum_{q}|W_{q}(\omega)|^{2}\left[{\rm
Im}\,\chi^{0}({\bf q},\omega)+\right.$ $\displaystyle\ +\left.|\chi^{0}({\bf
q},\omega)|^{2}{\rm Im}f_{xc}\right]\theta(\omega).$ (9)

This equivalence is implicit in Langreth’s analysis of the low energy behavior
of $\beta(\omega)$ for the homogeneous electron gas (HEG) in the random phase
approximation (RPA), where an adiabatic approximation is also valid
$\epsilon({\bf q},\omega)\approx\epsilon({\bf q},0)$. If the exchange
correlation kernel $f_{xc}$ is taken to be real, as in typical implementations
of TDDFT or ignored as in the RPA ($f_{xc}=0$), Eq. (10) reduces to

$\beta(\omega)=\frac{1}{\pi}\sum_{q}|W_{q}(\omega)|^{2}{\rm Im}\,\chi^{0}({\bf
q},\omega)\theta(\omega).$ (10)

This result in terms of the screened-core-hole potential $W_{q}(\omega)$ can
be advantageous for practical calculations in inhomogeneous systems. For
example, in the adiabatic approximation $W_{q}(\omega)\approx W_{q}(0)$, only
a single matrix inversion is needed to obtain $\epsilon^{-1}(\omega=0)$,
rather than an inversion at each frequency.

As a practical alternative to momentum-space, calculations of the Langreth
cumulant for inhomogeneous systems have recently been carried out by
transforming to real-space and real-time (RSRT).Kas et al. (2015) The approach
first calculates the time-evolved density response $\delta\rho({\bf r},t)$
with the Yabana-Bertsch reformulation of TDDFT that builds in a DFT exchange-
correlation kernelYabana and Bertsch (1996)

$\delta\rho({\bf r},t)\equiv\int d^{3}r^{\prime}dt^{\prime}\,\chi({\bf
r},t;{\bf r}^{\prime},t^{\prime})V_{c}({\bf r^{\prime}},t^{\prime}).$ (11)

A Fourier transform then yields the cumulant kernel $\beta(\omega)$

${\beta(\omega)}={\rm Re}\,\frac{\omega}{\pi}\int_{0}^{\infty}dt\,e^{-i\omega
t}\int d^{3}r\,V_{c}({\bf r})\delta\rho({\bf r},t)\,\theta(\omega).$ (12)

This approach has been shown to give good results for a number of systems.Kas
et al. (2015) However, the method requires a computationally demanding long-
time evolution of the density response using a large supercell, with Kohn-Sham
DFT calculations at each time-step. In addition, a post-processing convolution
is needed for XAS calculations.

###  II.3 Real-space Green’s function Cumulant

Our primary goal in this work is to develop an alternative, real-space Green’s
function formulation of the Langreth cumulant and it’s key ingredients. In
particular we aim to investigate the behavior of the cumulant kernel
$\beta(\omega)$ and the dynamically screened core-hole potential $W_{c}({\bf
r}.\omega)$. The method is complementary to the RSRT formulation but is based
on a similar approach for the response function, and can be carried out in
parallel with RSGF calculations of XAS.Rehr and Albers (1990); Kas et al.
(2021) Within the RPA, the RSGF formulation of $\beta(\omega)$ can be derived
from the perhaps better known GW approximation to the cumulant, based on the
core-level self-energy $\Sigma_{c}(\omega)$,Aryasetiawan et al. (1996);
Almbladh and Hedin (1983); Guzzo et al. (2011)

$\beta(\omega)=\frac{1}{\pi}{\rm Im}\,\Sigma_{c}(\epsilon_{c}-\omega),$ (13)

where $\Sigma_{c}(\omega)=\langle c|\Sigma(\omega)|c\rangle$. This matrix
element can be evaluated in real-space using the GW approximation

$\displaystyle{\rm Im}\,\Sigma({\bf r},{\bf r^{\prime}},\omega)=$
$\displaystyle-\sum_{i}^{occ}\psi_{i}({\bf r})\psi_{i}({\bf
r^{\prime}})\times$ $\displaystyle{\rm Im}[W({\bf r},{\bf
r^{\prime}},\epsilon_{i}-\omega)]\theta(\epsilon_{i}-\omega),$ (14)

where $W({\bf r},{\bf r^{\prime}},\omega)$ is the dynamically screened Coulomb
interaction. Within the decoupling approximation,Guzzo et al. (2011) the core-
level wave-function is assumed to have no overlap with any other electrons,
and Eq. (13) becomes

$\displaystyle\beta(\omega)$ $\displaystyle=$ $\displaystyle\frac{1}{\pi}{\rm
Im}\\!\\!\int\\!\\!d^{3}rd^{3}r^{\prime}\rho_{c}({\bf r})W({\bf r},{\bf
r^{\prime}},\omega)\rho_{c}({\bf r}^{\prime})\theta(\omega)$ (15)
$\displaystyle=$ $\displaystyle\frac{1}{\pi}{\rm Im}\\!\\!\int
d^{3}r\,\rho_{c}({\bf r})W_{c}({\bf r},\omega)\theta(\omega),$ (16)

where $W_{c}({\bf r},\omega)=\int d^{3}r^{\prime}\,W({\bf r},{\bf
r^{\prime}},\omega)\rho_{c}({\bf r}^{\prime})$ and $\rho_{c}({\bf
r})=|\psi_{c}({\bf r})|^{2}$. Then noting that ${\rm Im}\,W_{c}={\rm
Im}\,[K\chi V_{c}]$ (indices suppressed for simplicity) and within the RPA
($K=v$), we obtain a real-space generalization of the Langreth cumulant in Eq.
(II.2)

$\beta(\omega)=\frac{1}{\pi}\int d^{3}rd^{3}r^{\prime}\ V_{c}({\bf
r})V_{c}({\bf r^{\prime}})\,{\rm Im}\,\chi({\bf r},{\bf
r^{\prime}},\omega)\theta(\omega).$ (17)

This result is also equivalent to the RSRT expression in Eq. (12).
Alternatively, in analogy with Eq. (10), and again within the RPA,
$\beta(\omega)$ can be expressed in terms of $W_{c}({\bf r},\omega)$ and the
independent particle response function

$\displaystyle\beta(\omega)=$ $\displaystyle\frac{1}{\pi}\int
d^{3}rd^{3}r^{\prime}\,W_{c}^{*}({\bf r},\omega)W_{c}({\bf
r}^{\prime},\omega)\times$ $\displaystyle\times{\rm Im}\,\chi^{0}({\bf r},{\bf
r}^{\prime},\omega)\theta(\omega).$ (18)

While Eqs. (16-18) are formally equivalent within the RPA, they differ if the
interaction kernel $K$ is complex (which is the case for most non-adiabatic
kernels), in which case it is not obvious which approximation is best. Here we
focus on the RPA expressions only, although a generalization to adiabatic
$f_{xc}$ would be relatively simple.

The static limit of Eq. (18) is interesting in itself. At frequencies well
below $\omega_{p}$, the core-hole potential is strongly screened beyond the
screening length $r_{0}$ and nearly static. On expanding $\chi({\bf r},{\bf
r}^{\prime},\omega)$ about $\omega=0$ and keeping only the leading terms, one
obtains the adiabatic approximation

$\beta(\omega)\approx\frac{1}{\pi}\int d^{3}rd^{3}r^{\prime}\,W_{c}^{*}({\bf
r},0)W_{c}({\bf r}^{\prime},0)\,{\rm Im}\,\chi^{0}({\bf r},{\bf
r}^{\prime},\omega).$ (19)

This limiting behavior is similar to the adiabatic TDDFT approximation for
XAS,Ankudinov et al. (2005) but with the replacement of the dipole operator
$d({\bf r})$ by the statically screened core-hole potential $W_{c}({\bf
r},0)$. This potential is also used in calculations of XAS to approximate the
particle-hole interaction, and is similar to the final state rule
approximation for the static core-hole potential.

The transformation in Eq. (16) emphasizes the localization of $\beta(\omega)$,
which only depends on the imaginary part of the screened core-hole potential
${\rm Im}\,W_{c}({\bf r},\omega)$ over the range of the core density
$\rho_{c}({\bf r})$. The above identities also illustrate the connection
between the local and extended behavior of the dynamical response function.
This connection is analogous, e.g., to the origin of fine-structure in XAS,
where back-scattering is responsible for the fine-structure in photoelectron
wave-function at the origin.

Not surprisingly, since both quantities are physically related to dielectric
response, calculations of $\beta(\omega)$ are formally similar to those for
XAS.Prange et al. (2009)

$\mu(\omega)=\frac{4\pi}{V}\int d^{3}rd^{3}r^{\prime}\,d({\bf r})d({\bf
r}^{\prime})\,{\rm Im}\,\chi({\bf r},{\bf r}^{\prime},\omega),$ (20)

where ${d}(\bf r)$ is the (e.g., dipole) transition operator. Analogous
expressions have been derived for TDDFT approximations of atomic
polarizabilities Zangwill and Soven. (1980); Stott and Zaremba (1980) and
optical absorption spectra.Ankudinov et al. (2005) The quasi-particle XAS
$\mu_{1}(\omega)$ in our RSGF calculationsRehr and Albers (1990); Kas et al.
(2021) is obtained by evaluating $\mu(\omega)$ using the final-state rule,
i.e., with a screened core-hole in the final state, which corresponds to the
adiabatic approximation of Eq. (19). Physically the absorption in
$\mu(\omega)$ can be viewed in terms of the damped oscillating electric dipole
moment of the system induced by an external electric dipole potential
oscillating at frequency $\omega$. The kernel $\beta(\omega)$ can be viewed
similarly, except that the dipole potential $d({\bf r})$ is replaced by the
core-hole potential $V_{c}({\bf r})$. Thus a major difference is that the
suddenly turned on core-hole induces a damped oscillating monopole response
field about the absorbing atom.

Formally $W_{c}({\bf r},\omega)$ can be calculated in terms of the inverse
TDDFT dielectric matrix in real-spaceZangwill and Soven. (1980)

$\displaystyle W_{c}({\bf r},\omega)$ $\displaystyle=$ $\displaystyle\int
d^{3}r^{\prime}\,\epsilon^{-1}({\bf r},{\bf r}^{\prime},\omega)V_{c}({\bf
r}^{\prime}),$ (21) $\displaystyle\epsilon({\bf r},{\bf r}^{\prime},\omega)$
$\displaystyle=$ $\displaystyle\delta({\bf r},{\bf
r}^{\prime})\\!-\\!\\!\int\\!d^{3}r^{\prime\prime}K({\bf r},{\bf
r}^{\prime\prime})\chi^{0}({\bf r}^{\prime\prime},{\bf r}^{\prime},\omega),$
(22)

Related expressions for $W_{c}({\bf r},\omega)$ in atoms have also been
reported.Shirley and Martin (1993) Alternatively $W_{c}({\bf r},\omega)$ can
be obtained by iterating the integral equation $W_{c}=V_{c}+K\chi^{0}W_{c}$ to
self-consistency.Zangwill and Soven. (1980) Yet another tack is the use of
fluctuation potentials, e.g., in the quasi-boson approach.Hedin (1999) These
are obtained by diagonalizing the dielectric matrix using an eigenvalue
problem for each frequency $\omega$,

$\int d^{3}r^{\prime}\epsilon({\bf r},{\bf r}^{\prime},\omega)w_{s}({\bf
r}^{\prime},\omega)=\lambda_{s}(\omega)w_{s}({\bf r},\omega),$ (23)

The eigenfunctions $w_{s}(\omega)$ are the same as those for $K\chi^{0}$,
which has eigenvalues $\kappa_{s}(\omega)=1-\lambda_{s}(\omega)$. This
approach is particularly useful near the quasi-bosonic resonances
$\omega_{s}=\omega_{p}$ where ${\rm Re}\,\lambda_{p}(\omega)$ crosses zero and
matrix inversion can be numerically unstable. The corresponding excitation
energies $\omega_{p}$ and eigenfunctions $w_{p}({\bf r},\omega)$ define the
fluctuation potentials.Hedin (1999); Hedin et al. (1998) Close to
$\omega_{p}$, $\lambda_{p}(\omega)$ varies linearly,

$\lambda_{p}(\omega)\approx\lambda^{\prime}_{p}(\omega-\omega_{p})+i\Gamma_{p},$
(24)

where $\lambda^{\prime}_{p}=d[{\rm
Re}\,\lambda_{p}(\omega)]/{d\omega}|_{\omega_{p}}$. This approximation yields
a Lorentzian behavior for the quasi-bosonic peaks in $\beta(\omega)$ of width
$\gamma_{p}=\Gamma_{p}/\lambda_{p}^{\prime}$, which also limits the range of
the density fluctuations and the screened potential $w_{p}({\bf
r},\omega_{p}$). Nevertheless, we have found that an explicit matrix inversion
on a finite real-space basis (see below) usually converges well for all
frequencies due to the finite imaginary part of the dielectric screening
$\epsilon(\omega)$.

Formally, the Lindhard expression for the independent particle response
function on the real-axis can be expressed in terms of the one-particle
Green’s function Prange et al. (2009)

$\displaystyle\chi^{0}({\bf r},{\bf r}^{\prime},\omega)=$ $\displaystyle
2\int^{E_{F}}_{-\infty}dE\,[\rho({\bf r},{\bf r}^{\prime},E)G({\bf
r}^{\prime},{\bf r},E+\omega)$ $\displaystyle+\rho({\bf r^{\prime}},{\bf
r},E)G^{*}({\bf r},{\bf r}^{\prime},E-\omega)],$ (25)
$\displaystyle\frac{-1}{\pi}{\rm Im}\,\chi^{0}({\bf r},{\bf
r}^{\prime},\omega)$ $\displaystyle=2\int_{E_{F}-\omega}^{E_{F}}dE\,\rho({\bf
r},{\bf r}^{\prime},E)\rho({\bf r}^{\prime},{\bf r},E+\omega).$ (26)

Here $\rho({\bf r^{\prime}},{\bf r},E)$ is the spectral density of the one-
particle Green’s function which under the integrals can be replaced with
$(-1/\pi){\rm Im}\,G({\bf r^{\prime}},{\bf r},E)$ (cf. Eq. (33) and (34) of
Ref Zangwill and Soven., 1980). The convolution in ${\rm
Im}\,\chi^{0}(\omega)$ defines a particle-hole spectral function. The
calculation of $W_{c}({\bf r},\omega)$ thus requires both Re and ${\rm
Im}\,\chi^{0}({\bf r},{\bf r}^{\prime},\omega)$, and hence expressions for
both $\rho({\bf r},{\bf r}^{\prime},E)$ and $G({\bf r},{\bf
r}^{\prime},\omega)$. Explicit derivations and algorithms for these functions
are given in Ref. Prange et al., 2009, along with algorithms for calculating
XAS and optical response, based on summations over a finite cluster
surrounding the absorbing atom. Below, we show how $\chi^{0}({\bf r},{\bf
r}^{\prime},\omega)$ and $\chi({\bf r},{\bf r}^{\prime},\omega)$ can be
calculated within a generalized real-space multiple-scattering Green’s
function (RSGF) approach using a discrete site-radial coordinate basis.

###  II.4 Real-Space Multiple Scattering Formalism

The key ingredients in the RSGF calculation of $W_{c}({\bf r},\omega)$ are the
bare response function $\chi^{0}({\bf r},{\bf r}^{\prime},\omega)$ which is
defined in Eq. (23) in terms of the independent particle Green’s functions
$G({\bf r},{\bf r}^{\prime},E)$, and the interaction kernel $K=v+f_{xc}$ which
here is approximated by the RPA $K=v$. The most demanding step is the
inversion of the non-local dielectric matrix $\epsilon({\bf r},{\bf
r}^{\prime},\omega)$ on a 3-$d$ grid, which can be computationally formidable.
In order to simplify the calculation, we employ a generalization of real-space
multiple-scattering (RSMS) formalism with a discrete site-radial coordinate
basis, analogous to that developed for x-ray spectra and optical response.Rehr
and Albers (1990); Prange et al. (2009) In this approach, space is partitioned
into cells $i$ with cell boundaries defined by $\Theta_{i}({\bf r})=1\ (0)$
for points inside (outside) a given cell. Formally the cell boundaries should
be defined by Voronoi or equivalent partitioning. However, for simplicity, the
cells are taken to be Norman spheres (i.e., spheres of charge neutrality).Jr.
(1976); Prange et al. (2009) Thus $\Theta_{i}({\bf r})=\theta(r_{N}-r_{i})$,
where $\theta(r)$ is the unit step function, and $r_{N}$ denotes the Norman
radius. The points within each cell are then represented in spherical polar
coordinates ${\bf r}_{i}=(r_{i},\Omega_{i})$. In addition the radial
coordinates are typically defined on a discretized logarithmic grid, e.g.,
$r_{i,n}=\exp(-x_{0}+n\delta x)$, as in SCF atomic calculations,Desclaux
(1975) while the angular dependence is represented by spherical harmonics
$Y_{L_{i}}(\hat{r}_{i})$, where $L_{i}=(l_{i},m_{i})$. Thus the spatial points
${\bf r}$ are represented by the discrete indices $I=(i,n_{i},L_{i})$, which
have a finite volume $\Delta_{i}=4\pi r_{i}^{2}\delta_{i}$, with
$\delta_{i}\approx r_{i}\delta x$ being the radial grid spacing (e.g.,
$x_{0}=-8.8$ and $\delta x$ = 0.05). This construction conserves charge and
simplifies volume integration, i.e., for a given volume $\int
d^{3}r=\sum_{i}\int^{r_{N}}r_{i}^{2}dr_{i}d\Omega_{i}$. With this
representation, functions of coordinates $F({\bf r})\rightarrow F_{I}$ are
represented by vectors, and $F({\bf r},{\bf r}^{\prime})\rightarrow F_{I,J}$
by matrices in $I$ and $J$ of rank $N_{r}N_{i}(l_{max}+1)^{2}$, where
$l_{max}$ denotes the maximum angular momentum used in the calculation.

This contruction is then implemented within the standard RSMS theory, which
assumes spherical symmetry of the scattering potentials within each cell. The
independent particle Green’s functions are given by,Rehr and Albers (1990)

$\displaystyle G({\bf r},{\bf r^{\prime}},E)=$
$\displaystyle-2k\sum_{iLjL^{\prime}}\left[\delta_{ij}\delta_{LL^{\prime}}R_{iL}({\bf
r}_{i<}{\bf r}_{i>})+\right.$ $\displaystyle\left.R_{iL}({\bf
r_{i}})\tilde{G}_{iLjL^{\prime}}(E)R_{jL^{\prime}}({\bf
r}^{\prime}_{j})\right]\Theta_{i}({\bf r})\Theta_{j}({\bf r^{\prime}})$
$\displaystyle=$ $\displaystyle\sum_{iLjL^{\prime}}\left[G^{c}_{iL}({\bf
r}_{i},{\bf r}^{\prime}_{j},E)\delta_{ij}\delta_{LL^{\prime}}\right.+$
$\displaystyle\left.G^{sc}_{iLjL^{\prime}}({\bf r}_{i},{\bf
r}^{\prime}_{j},E)\right]\Theta_{i}({\bf r})\Theta_{j}({\bf r}^{\prime}).$
(27)

Here $R_{iL}({\bf r})=i^{l}(R_{il}(r)/r)Y_{L}({\bf\hat{r}})$ and $H_{iL}({\bf
r})=i^{l}(H_{l}(r)/r)Y_{L}({\bf\hat{r}})$ are the regular and irregular
solutions to the single site Dirac equation for the atom at site $i$, ${\bf
r_{i}=r-R_{i}}$ is the position relative to the center of the $i^{th}$ cell,
and $\tilde{G}_{iLjL^{\prime}}(E)$ is the multiple scattering matrix.Prange et
al. (2009) The radial wave-functions $R_{L}({\bf
r},E)=i^{l}(R_{l}(r)/r)Y_{L}(\hat{r})$ are scattering-state normalized such
that $R_{l}(r)\rightarrow\sin(kr-l\pi/2+\delta_{l})$ beyond the muffin-tin
radius $r_{mt}$,Ankudinov (1996) where $\delta_{l}$ are partial wave phase
shifts and $k=\sqrt{2E}$ is the photoelectron momentum. Thus,

$G({\bf r,\bf
r^{\prime}},E)=\sum_{iLjL^{\prime}}Y_{L}^{*}(\hat{r}_{i})G_{iLjL^{\prime}}({r}_{i},{r}^{\prime}_{j})Y_{L^{\prime}}(\hat{r^{\prime}_{j}}),$
(28)

where the cell functions $\Theta_{i}$ and $\Theta_{j}$ have been absorbed into
the definition of $G_{iLjL^{\prime}}(r_{i},r^{\prime}_{j},E)$. On the real
axis, the density matrix ${\rho={\rm Im}\,G}$ can be expressed similarly,

$\displaystyle\rho({\bf r},{\bf r^{\prime}})$
$\displaystyle=\sum_{iLjL^{\prime}}R_{iL}({\bf
r_{i}})\tilde{\rho}_{iLjL^{\prime}}(E)R_{jL^{\prime}}({\bf
r^{\prime}_{j}})\Theta_{i}({\bf r})\Theta_{j}({\bf r^{\prime}}),$
$\displaystyle=\sum_{iLjL^{\prime}}Y_{L}(\hat{r})\rho_{iLjL^{\prime}}(r_{i},r^{\prime}_{j},E)Y_{L^{\prime}}(\hat{r^{\prime}}).$
(29)

Additional details of the method and algorithms for calculations of the
quantities involved in Eq. (27) can be found elsewhere.Rehr and Albers (1990);
Prange et al. (2009) With the above representation of the Green’s function and
density matrix, the bare response function $\chi^{0}(\omega)$ can be expressed
similarly in a site, radial-coordinate, and angular momentum basis,

$\displaystyle\chi^{0}({\bf r},{\bf
r}^{\prime},\omega)=\sum_{iLjL^{\prime}}Y_{L}(\hat{r}_{i})\chi^{0}_{iLjL^{\prime}}({r}_{i},{r}^{\prime}_{j},\omega)Y_{L^{\prime}}({\hat{r}}^{\prime}_{j}).$
(30)

Likewise, we expand the interaction kernel $K$ in spherical harmonics about
each cell center, and for simplicity given the near spherical symmetry of the
core-hole, keep only the spherical terms $L=L^{\prime}=0$ . If we ignore the
TDDFT contribution $f_{xc}=\partial V_{xc}(\rho)/\partial\rho$, as in the RPA,
the kernel becomes

$\displaystyle K({\bf r},{\bf r}^{\prime})$
$\displaystyle\approx\sum_{ij}Y_{00}(\hat{r}_{i})K_{i0j0}(r_{i},r^{\prime}_{j})Y_{00}(\hat{r}^{\prime}_{j})\Theta_{i}({\bf
r})\Theta_{j}({\bf r^{\prime}}),$ $\displaystyle
K_{i0j0}(r_{i},r^{\prime}_{j})=\left[\frac{4\pi}{r^{i}_{>}}\delta_{ij}+\frac{4\pi}{R_{ij}}(1-\delta_{ij})\right],$
(31)

where $R_{ij}=1/|{\bf R}_{i}-{\bf R}_{j}|$, and
$r_{i>}=\max(r_{i},r^{\prime}_{i})$.

As a further simplification, we assume here that the deep core-hole density
and potential, and hence the screened core-hole potential are spherically
symmetric about the central site,

$\displaystyle\rho_{c}({\bf r})$
$\displaystyle=\rho_{c0}(r)Y_{00}({\bf\hat{r}}),$ $\displaystyle V_{c}({\bf
r})$ $\displaystyle=V_{c0}(r)Y_{00}({\bf\hat{r}}),$ $\displaystyle W_{c}({\bf
r},\omega)$ $\displaystyle=w_{c}(r,\omega)Y_{00}({\bf\hat{r}}).$ (32)

Then $W_{c}({\bf r},\omega)$ as well as $\beta(\omega)$ depend only on the
spherical $L=L^{\prime}=0$ components of the bare response function. The
matrix inversion of ${\bf{\bm{\epsilon=1-K\chi_{0}}}}$ can then be expressed
in terms of a matrix inverse in the site and radial coordinates

$\bm{\chi}=\bm{\chi^{0}}[{\bf 1-K}\bm{\chi^{0}}]^{-1}.$ (33)

From Eq. (24), the matrix elements $\bm{\chi^{0}}_{IJ}$ in the discrete basis
$I=(i,r_{i})$ and $J=(j,r^{\prime}_{j})$ are

$\displaystyle\bm{\chi^{0}}_{I,J}(\omega)=\frac{1}{2\pi}\int_{-\infty}^{E_{F}}\sum_{LL^{\prime}}\Delta_{i}^{1/2}\Delta_{j}^{1/2}\rho_{iLjL^{\prime}}(r_{i},r^{\prime}_{j},E)\times$
$\displaystyle\
\left[G_{jL^{\prime}iL}(r^{\prime}_{j},r_{i},E+\omega)+G^{*}_{jL^{\prime}iL}(r^{\prime}_{j},r_{i},E-\omega)\right].$
(34)

Examining the above equations and Eq. (27), we see that site off-diagonal
terms $i\neq j$ contribute to the full response function in several ways: i)
off-diagonal terms in $\chi^{0}$ corresponding to the free propagation of a
particle-hole from one site to another; ii) off-diagonal terms in the
interaction kernel $K$, corresponding to the coupling of a particle-hole state
at site $i$ with another at site $j$; iii) a combination of these two,
corresponding to the propagation of a particle-hole state from one site to an
intermediate site, then scattering into another particle-hole state at a third
site; and iv) the Green’s function itself has contributions due to scattering
of the single particle (photoelectron or valence hole) from neighboring sites.

To summarize, the expressions for $\beta(\omega)$ in the discrete site-radial
coordinate basis are

$\displaystyle\beta(\omega)$ $\displaystyle={\rm Im}\left[\bm{\rho_{c0}}{\bf
W_{c0}}\right]\theta(\omega),$ $\displaystyle={\rm Im}[{\bf V_{c0}\chi
V_{c0}}]\theta(\omega),$ $\displaystyle={\rm Im}\left[{\bf
W_{c0}^{*}}\bm{\chi^{0}}{\bf W_{c0}}\right]\theta(\omega),$ (35)

where the vectors are given by

$\displaystyle\bm{\rho_{c0}}$ $\displaystyle=\bm{\rho_{cI}}\delta_{i0},$
$\displaystyle{\bf V_{c0}}$ $\displaystyle={\bf K}\bm{\rho_{c0}},$
$\displaystyle{\bf W_{c0}}$ $\displaystyle=({\bf 1-K}\bm{\chi^{0}})^{-1}{\bf
V_{c0}}.$ (36)

The last line in Eq. (35) can also be recast in terms of matrix elements of
the screened core-hole potential

$W_{iL}(E,\omega)=\int^{r_{N}}drR_{il}(r_{i},E)w_{c}(r,\omega)R_{il}(r_{i},E+\omega).$
(37)

This reveals a striking resemblance to the RSMS theory of XAS,Rehr and Albers
(1990) or more precisely, that of optical response,Prange et al. (2009) the
difference being that the dipole transition operator ${d(r)}$ is replaced by
the dynamically screened core-hole potential $w_{c}(r,\omega)$,

$\displaystyle\beta(\omega)=\int_{E_{F}}^{E_{F}+\omega}$ $\displaystyle
dE\sum_{iLjL^{\prime}}W_{iL}^{*}(E,\omega)\rho_{iLjL^{\prime}}(E)\times$
$\displaystyle\rho_{iLjL^{\prime}}(E+\omega)W_{jL^{\prime}}(E,\omega).$ (38)

###  II.5 Local RSGF cumulant approximation

The first line in Eq. (35) shows that the screened core-hole potential is only
needed locally about the absorbing atom, since it is multiplied by the density
of the core-orbital. For example, for the $1s$ state of Na, the required range
is only $\sim 1/11$ Bohr. This suggests that it may be possible, at least in
some systems, to treat the problem locally, neglecting site off-diagonal terms
in $K$ and $\chi^{0}$, while in other systems it may be sufficient to
approximate the contribution from off diagonal terms. To do this, the Green’s
functions in $\chi^{0}$ are taken to be those of an atom embedded in an
electron gas at the interstitial density $r_{s}$, for points $r,r^{\prime}$
outside the Norman radius of the absorbing atom, while the full Green’s
function is used within the Norman radius. Thus outside the absorbing cell,
the Green’s function is approximated simply in terms of phase-shifted
spherical Bessel functions centered about the absorbing site. The dielectric
matrix $\bm{\epsilon}=\bf{1-K}\bm{\chi^{0}}$ is then inverted on a single
radial grid, which greatly simplifies the calculations. In this case the
matrix equations Eq. (33) and (36) still apply, but are limited to $i=0$, and
the radius $R_{\rm max}$ defining the central cell is treated as a convergence
parameter. The equation for $\beta(\omega)$ then becomes

$\beta(\omega)=\frac{1}{\pi}\int_{0}^{R_{\rm
max}}dr\,\rho_{c}^{0}(r)w_{c}^{0}(r,\omega).$ (39)

If we further assume that the Green’s function has the structure of that of a
spherical system, i.e., $G_{L,L^{\prime}}=G_{l}\,\delta_{LL^{\prime}}$, one
can express $\beta(\omega)$ in terms of a joint density of states, [cf. Eq.
(32) of Prange et al.],Prange et al. (2009)

$\displaystyle\beta(\omega)=2\sum_{l}(2l+1)$
$\displaystyle\int_{E_{F}-\omega}^{E_{F}}dE\,|W_{l}(E,\omega)|^{2}\rho_{l}(E)\times$
$\displaystyle\rho_{l}(E+\omega)\theta(\omega),$ (40)

where $\rho_{l}(E)$ is the angular-momentum projected density of states (LDOS)
for a given $l$ calculated at the central atom, and the factor of 2 accounts
for spin degeneracy.

With spherical symmetry of the core-hole density, calculations of
$\beta(\omega)$ reduce to a form with only radial coordinates. The spherical
part of the independent particle response function
$\chi^{0}({r},{r}^{\prime},\omega)|_{l=0}$, Prange et al. (2009) is

$\displaystyle\chi^{0}_{l=0}({r},{r}^{\prime},\omega)$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}\int_{-\infty}^{E_{F}}dE\,\sum_{l}(2l+1)\rho_{l}({r},{r}^{\prime},E)\times$
$\displaystyle\times\,$
$\displaystyle[G_{l}({r},{r}^{\prime},E+\omega)+G_{l}^{*}({r},{r}^{\prime},E-\omega)].$
(41)

Here we have used the symmetry and reality of $\rho_{l}(r,r^{\prime},E)$ and
$G_{l}(r,r^{\prime},E)$. In order to stabilize this expression, we split the
two terms of the energy integral as follows

$\displaystyle\chi^{0}_{l=0}(r,r^{\prime},\omega)$ $\displaystyle=$
$\displaystyle\frac{1}{2\pi}\sum_{l}(2l+1)[\chi^{0}_{l}(I)+\chi^{0}_{l}(II)],$
(42) $\displaystyle\chi^{0}_{l}(I)$ $\displaystyle=$
$\displaystyle\int_{E_{F}-\omega}^{E_{F}}dE\
\rho_{l}(r,r^{\prime},E)G_{l}(r,r^{\prime},E+\omega),$
$\displaystyle\chi^{0}_{l}(II)$ $\displaystyle=$
$\displaystyle-\frac{1}{\pi}{\rm
Im}\,\int^{E_{F}}dE\,G_{l}(r,r^{\prime},E)G_{l}(r,r^{\prime},E-\omega).$

The second of the above integrals can be performed in the complex energy plane
where the Green’s function is smoother. Note also that these expressions
implicitly include fine-structure in $G_{l}(r,r^{\prime},\omega)$ from
multiple-scattering from atoms beyond the central atom. Finally the screened
core hole potential $w_{c}(r,\omega)$ is defined by a the radial integral,

$w_{c}(r,\omega)=\int
dr^{\prime}\,[1-K\chi^{0}(\omega)]^{-1}_{l=0}(r,r^{\prime})v_{c}(r^{\prime}),$
(43)

which can be calculated directly by matrix inversion in the radial coordinate
basis for the single site $i=0$. We use this local RSGF approximation in all
of our results shown below.

Eq. (40) is also consistent with the Anderson edge-singularity exponent
$\alpha$ for XAS. For metallic systems, $\beta(\omega)\approx\alpha\ \omega$
is roughly linear in frequency for small $\omega$, reflecting the behavior of
the joint density of states near the Fermi energy. The coefficient $\alpha$
can be determined from the zero frequency limit, with matrix elements given by
their values at the Fermi level. Since the screened core-hole potential is
well localized and spherically symmetric near $\omega=0$,

$\alpha=2\sum_{l}(2l+1)|W_{l}(E_{F},\omega=0)|^{2}\rho_{l}(E_{F})^{2}.$ (44)

This result is consistent with that derived by Anderson, Nozìeres and De
DominicisNozìeres and De Dominicis (1969)
$\alpha=2\Sigma_{l}(2l+1){(\delta_{l}/\pi)^{2}}$. This can be verified by
noting that $\rho_{l}(E)=dN_{l}/dE$, and $W_{l}(E_{F})=\Delta E_{l}$ is the
change in energy of levels $l$ due to the screened core-hole potential
$W(r,0)$, so $W_{l}(E_{F})\rho_{l}(E_{F})=\Delta N_{l}=\delta_{l}/\pi$, i.e.,
the screening charge in level $l$, as in the Friedel sum rule. This singular
behavior shows up as an asymmetry in the main peak of the XPS. The XAS has an
additional contribution from the Mahan edge singularity exponent
$\alpha_{l}=-2\delta_{l}/\pi$, where the photoelectron has local angular
momentum $l$.Rehr et al. (2020) Even in insulators, the low energy background
terms tend to grow linearly beyond the gap, leading to an asymmetry in the
quasi-particle peak, though without a true singularity.

##  III Details of the Calculations

The RSGF formalism used here has been implemented within the FEFF10 code,Kas
et al. (2020) with the RPA approximation for the TDDFT kernel $K=1/|{\bf
r}-{\bf r}^{\prime}|$. The one-particle Green’s functions are calculated using
self-consistent-field (SCF) potentials, full-multiple-scattering,Rehr and
Albers (1990); Kas et al. (2021) and unless noted otherwise, the final-state-
rule approximation with a screened core hole in the final state. The RSGF
calculations of $w_{c}(r,\omega)$ depend on the radius $R_{max}$ beyond which
the response functions are truncated, the maximum angular momentum $l_{max}$,
the radius $r_{fms}$ of the cluster of atoms used in the full-multiple-
scattering calculations, and a small imaginary part $\gamma$ added to the
energy of the Green’s functions; the value of $\gamma$ is set to $0.1$ eV for
all calculations except for the electron gas where it is set to $0.01$ eV.
$R_{max}$ is treated as a convergence parameter for metallic systems and set
to the Norman radius for insulators. By default, we use a dense logarithmic
grid e.g., $r_{i,n}=\exp(-8.8+0.05n)$, as in the Dirac-Fock atomic
calculations in FEFF10.Desclaux (1975) Typically this amounts to about 200
points per cell for each frequency, which is computationally manageable,
though a sparser grid may be adequate in some cases. For the calculations of
the cumulant spectral function, broadening by a Voigt function with both
Gaussian and Lorentzian half-widths of 0.25 eV is used for all systems.
Finally, for calculations of the 2p XPS, several additional parameters were
used. First, a parameter is introduced to account for the spin-orbit splitting
between the $2p_{1/2}$ and $2p_{3/2}$ contributions; second, a Gaussian
broadening with width $\Gamma_{exp}$ was applied to the spectrum to account
for experimental broadening, and separate Lorentzian widths
$\Gamma_{2p_{1/2}}$ and $\Gamma_{2p_{3/2}}$ to account the different for core-
hole lifetimes; third, a Shirley backgroundShirley (1972) was calculated from,
and then added to the theoretical spectrum with a scaling parameter $B$ to
match the low energy limit of the experimental data. Below we present results
only for the local RSGF spherical approximation, which is found to be semi-
quantitative in comparison to experiment for the systems treated here. A more
complete treatment including longer-ranged contributions with the generalized
RSMS formalism is reserved for the future.

##  IV Results

###  IV.1 Homogeneous electron gas

As a first illustration and proof of principle, we present results for the
homogeneous electron gas (HEG). An approximate treatment of the cumulant for
this system was described by Langreth,Langreth (1970) and in more detail by
Lundqvist and Hedin Lundqvist (1967); Hedin (1999) for the plasmon-pole (PP)
and RPA approximations. For a deep core-level, the core charge density can be
approximated by a point charge, so that the excitation spectrum
$\beta(\omega)$ is given by a one-dimensional integral in momentum space,

$\beta(\omega)=\frac{2}{\pi^{2}}\int dq\,{\rm
Im}\,[\epsilon^{-1}(q,\omega)]\theta(\omega).$ (45)

This can be computed numerically within the RPA, or analytically with the
plasmon pole model.sup  The results are shown in Fig. 1, where
$\beta(\omega)/\omega$ calculated with RSGF at several truncation radii
$r_{max}$ are compared to the numerically exact $q-$space results and with the
single plasmon-pole model.Lundqvist (1967) Consistent with the $q$-space
formulation, the behavior of $\beta(\omega)$ exhibits a linear part at low
energy due to the particle-hole continuum, where $-{\rm Im}\,\epsilon({\bf
q},\omega)^{-1}\approx{\rm Im}\,\epsilon(q,\omega)/|\epsilon(q,0)|^{2}$ and
$\epsilon(q,0)$ is finite, as well as a sharply peaked structure above the
plasmon onset $\omega_{p}\approx 6.0$ eV corresponding to the zeros of
$\epsilon(q,\omega)$. The adiabatic approximation is also shown, which
compares well at low frequencies.

Figure 1: Cumulant kernel $\beta(\omega)/\omega$ for the HEG at $r_{s}=4$
calculated using the RSGF approach compared to that based on the $q-$space
formula in Eq. (45) and the plasmon-pole model.Lundqvist (1967); Hedin (1999)
RSMS calculations are shown for $R_{max}=4r_{N}$ and $R_{max}=64\,r_{N}$,
where the Norman radius $r_{N}=1.4$ Å was chosen as that for bcc sodium. Note
that the low frequency, nearly constant behavior of $\beta(\omega)/\omega$,
gives the edge singularity index $\alpha$. The dashed purple curve shows the
adiabatic result which gives a good approximation for $\alpha$.

Note that the RSGF formulation (red) performs well at all energies except near
the plasmon peak, where it loses some spectral weight. This is likely due to
the truncations of $R_{max}=64\,r_{N}$ and $l_{max}=25$. The calculations also
agree well with the PP modelHedin (1999) above the plasmon frequency $\sim 6$
(eV); however the PP approximation completely ignores the low energy
contribution.sup  This low frequency tail in the HEG is dominated by the
particle-hole contribution to the dielectric function and is linear in
$\omega$. This behavior accounts for the Anderson edge singularity in x-ray
spectra and an asymmetric quasi-particle line-shape in XPS
$J_{k}(\omega)\sim\omega^{\alpha-1}$. For the HEG for $r_{s}=4$, we obtain
$\alpha=0.24$, which compares well to the experimental value $0.21\pm 0.015$
for the 1s XPS of sodium, a nearly free electron system with $r_{s}\approx 4$.
Citrin et al. (1977) We also analyzed the extent to which the calculation of
$\beta(\omega)$ can be approximated locally. The green curve in Fig. 1 shows
the result from RSGF calculations with $R_{max}=4r_{N}$ in the radial
arguments of $\chi^{0}(r,r^{\prime},\omega)$. While the linear portion is
reproduced very well, the position of the plasmon peak is overestimated by
$\sim 2$ (eV) or $30\%$. Although difficult to see, the high energy tail is
poorly represented and only corrected with much larger $r_{max}=64r_{N}$. The
RSGF spectrum approaches the $q-$space result with increasing $R_{max}$, but
does not converge until $R_{max}=64r_{N}\sim 177$ Bohr. This difference in
convergence in different energy ranges reflects the fact that the sharp peak
in $\beta(\omega)$ is largely defined by the sharp plasmon peak in the loss
function, which is dominated by momentum transfer $q$ near zero. On the other
hand, the linear behavior near $\omega=0$ and the long tail at high frequency
are due to the particle-hole continuum and the dispersion of the plasmon
respectively. These depend on higher values of momentum transfer which results
in faster convergence with $R_{max}$. Fig. 2 shows the screened core-hole
potential $w_{c}(r,\omega)$ (also for $r_{s}=4$) scaled for convenience by the
bare potential $V_{c}(r)$ vs $r$ for frequencies from $0$ eV to well above the
$\omega_{p}$. The RSGF result (solid) are compared with the results of the
plasmon pole (dot dashes) and those from the RSRT formulation (dashes).

Figure 2: Screened core-hole potential from the RSGF approach (solid) relative
to the bare core potential ${\rm Re}\,w_{c}(r,\omega)/V_{c}(r)$ (top) and
${\rm Im}\,w_{c}(r,\omega)/V_{c}(r)$ vs $r$ for the HEG at $r_{s}=4$ for
$\omega=0$, $\omega=8$ eV $\approx\omega_{p}$ , and $\omega=16{\rm
eV}>>\omega_{p}$. Note that the behavior becomes strongly delocalized,
oscillatory, and lossy near $\omega_{p}$. For comparison results from the RSRT
approach (dashes) and the plasmon pole model (dot dashes) are also shown.

As expected physically, there is little screening at small $r<<r_{s}$ for all
frequencies due to the slow build-up of the induced density response. For
$r>r_{s}$ the zero frequency curve (red) exhibits extremely efficient Yukawa-
like screening, while the behavior at high frequencies is weakly screened,
with $w_{c}(r,\omega)$ remaining finite at large $r$. Near the plasma-
frequency $\omega_{p}\approx 6$ eV $w_{c}(r,\omega)$ oscillates sinusoidally,
corresponding to delocalized charge density fluctuations with substantial
contributions to ${\rm Im\,}w(r,\omega)$. The oscillatory behavior matches
reasonably well with that of the plasmon pole model, although there is more
damping and larger screening in the high frequency curves calculated with
RSGF, and the oscillations are slightly phase shifted. Some of these
differences are expected due to differences in the plasmon dispersion and the
lack of particle-hole continuum states in the plasmon pole model. The majority
of the damping seen in the RSGF results comes from the truncation of the real-
space grid, although this damping does not affect the calculated values close
to the origin, which determine $\beta(\omega)$.

###  IV.2 Nearly free-electron system

As an example of the RSGF approach that can be compared with the HEG, we
present calculations for bcc Na, a nearly free electron (NFE) system with
$r_{s}\approx 4$. This is illustrated in Fig. 3 for a range of $r_{fms}$ and
for large $R_{max}$. The behavior of $\beta(\omega)$ is clearly linear below
$\omega_{p}$ and comparable to that of the HEG at high energies as well.
Moreover, the behavior near the plasmon peak also has long range oscillatory
contributions. For such NFE systems, scattering from near-near neighbors
appears to be less important than a large $R_{max}$ alone, as the results
converge reasonably well by setting $r_{fms}=0$. The oscillatory behavior near
the peak reflects the discrete nature of plasmon excitations within a sphere
of finite $R_{max}$.

Figure 3: RSGF calculation of the cumulant kernel $\beta(\omega)$ for Na (top)
and TiO2 (bottom, with varying FMS cluster size $r_{fms}$. Note that
scattering from the nearest-neighbor atoms has a drastic effect in TiO2 but is
much smaller for Na. Both spectra show minimal differences with theinclusion
of atoms beyond the nearest-neighbors.

###  IV.3 Transition metal compounds

Next we consider several early transition metal compounds, which are typically
moderately correlated. We first present results for TiO2 (rutile), which was
previously studied with the RSRT approach.Kas et al. (2015) The strong peak in
$\beta(\omega)$ for TiO2 is shown in Figs. 3 (bottom) and 4 (top). Unlike the
behavior of Na, the spectrum is largely independent of $R_{max}$, while
$r_{fms}$ has a drastic effect on the position and intensity of the main peak
in $\beta(\omega)$ as seen in Fig. 3. The independence of the spectrum with
respect to $R_{max}$ suggests a localized nature of the excitations in TiO2.
The main peak is neither directly related to peaks in $\chi^{0}(\omega)$ nor
in ${\rm Im}\,\epsilon(\omega)$. Instead the peak is consistent with zero-
crossings of ${\rm Re}\,\epsilon(\omega)$. For TiO2 we have found that there
is a single dominant eigenvalue of the dielectric matrix at $\omega_{p}$,
although the lack of more than one crossing may be due to the limited spacial
extent of our calculations. Consequently a fluctuation potential treatment
[cf. Eq. (22) and (23)] is appropriate. This PP model yields a Lorentzian
shape for $\beta(\omega)$ and ${\rm Im}\,w_{c}(r,\omega)$, as seen in Fig. 4
(top, green) as well as a finite range of $w_{c}(r,\omega)$ varying inversely
with $\gamma_{p}$.

Figure 4: Cumulant kernel $\beta(\omega)$ for TiO2 (Top) calculated with the
RSGF approach (red), and for comparison the RSRT (blue), the fluctuation
potential approach with a single pole (PP), and the adiabatic limit (dashed).
The lower figure for ScF3 (Bottom) compares results calculated with ground
state electronic structure (GS), and with the final state rule (FSR). The
core-hole shifts the peak to lower energy by $\sim 0.8$ eV, in better
agreement with experiment.

The cumulant kernel $\beta(\omega)$ calculated with Eq. (40) and
$R_{max}=r_{N}$ is shown in Fig. 4 (top), along with a comparison to the RSRT
method.

Figure 5: Spectral functions $A(\omega)$for Na (Top) compared with the HEG
with different $r_{fms}$. The middle figure compares the RSGF and RSRTKas et
al. (2015) approaches for TiO2, and the bottom figure compares the spectral
function of ScF3 calculated with (FSR) and without (GS) the final state rule
core-hole.

In contrast to the HEG, the adiabatic approximation (top, purple) does not
reflect the spectrum of $\beta(\omega)$, even at low frequencies, and is in
error by $\sim 50\%$ by $\approx 2.5$ eV, which is about where the first
excitation feature occurs. Instead, almost all structure in $\beta(\omega)$,
or at least that above $\approx 7$ eV, reflects that of $|W_{c}(\omega)|^{2}$.
Below $7$ eV, the spectrum has little structure. The prominant peak at about
15 eV stems from the single pole of $\epsilon^{-1}(\omega)$ where
$\lambda_{p}(\omega)$ and hence ${\rm Re}\,\epsilon(\omega)$ crosses zero. The
low energy structure is similar for both RSGF and RSRT methods, and the
agreement of the peak near 15 eV is semi-quantitative. Finally, at higher
energies there is appreciable difference between the RSGF and RSRT results,
which could be due to the limited basis set and lack of continuum states in
the RSRT results, which rely on local basis functions, or the local
approximation of the RGSF method. Note however, that the structure above 15 eV
does not contribute significantly to the spectral function, due to the factor
of $1/\omega^{2}$ in the cumulant $C_{c}(t)$ in Eq. (4). Results for ScF3
(bottom) are similar, with a large peak at about $13$ to $14$ eV depending on
whether the calculation was performed with (FSR), or without (GS) a core-hole.
The inclusion of the core-hole red-shifts the spectrum, and increases the main
peak intensity, similar to the effects seen in x-ray absorption near-edge
spectra. This excitonic shift improves the agreement with experiment, which
shows a peak at about $12$ eV (see Fig. 8). The spectral functions, which are
closely related to XPS, are shown in Fig. 5 for Na, TiO2 and ScF3, calculated
using Eq. (1). The spectral function of Na is compared with that of the
homogeneous electron gas, and is shown for several values of $r_{fms}$, while
that of TiO2 (middle) is compared to the result from RSRT, and shows good
agreement. Finally, the spectral function of ScF3 is shown as calculated with
(FSR) and without (GS) a final-state-rule core-hole. In these curves, we see a
large main peak at high energy, corresponding to the quasiparticle peak, and
satellite peaks at lower energy, which reflect the behavior of
$\beta(\omega)/\omega^{2}$. The quasiparticle peak in the spectral function of
Na shows an appreciable asymmetry, corresponding to the edge singularity, and
originating from the linear behavior in $\beta(\omega)$ at low frequency.

Figure 6: Screening fraction
$f(r,\omega)=[w_{c}(r,\omega)-V_{c}(r,\omega)]/V_{c}(r)$ relating the
dynamically screened $w_{c}(r,\omega)$ and bare $V_{c}(r)\approx 1/r$ core-
hole potentials for the Ti core-hole in TiO2 at a few selected energies. Note
that the screening is local and strong at low (black), weak at high (red), and
becomes lossy and delocalized near $\omega_{CT}$ with roughly a Lorentzian
line-shape. The behavior at large $r$ is only semi-quantitative due to the
truncation of $r_{max}$.

The dynamical screening fraction of the core-hole potential which here is
defined as $f(r,\omega)=[w_{c}(r,\omega)-V_{c}(r)]/V_{c}(r)$, calculated from
Eq. (43) with the local RSGF cumulant is shown for Rutile TiO2 in Fig. 6 for a
few selected energies. The bare core-hole potential $V_{c}(r)$ is that for the
$1s$ core-level of Ti and is nearly Coulombic for a point charge, i.e.,
$V_{c}(r)\approx 1/r$ beyond the Ti $1s$ radius $r_{1s}\approx 1/22$ Bohr. The
screening fraction $f(r,\omega)$ is small near $r=0$ (i.e., little screening),
and fairly strong $f(\omega)\approx-1$ beyond the Thomas-Fermi screening
length $r_{0}$. However, near the peak at $\omega_{p}\approx 14.8$ eV, abrupt
changes occur, with the system becoming over-screened ($f(\omega)$ reaching
about $-2$ at $14$ eV) or under-screened ($f(\omega)<0.25$) above $15.4$ eV.
The imaginary part of $W_{c}$ largely reflects the structure in
$\beta(\omega)$, starting at $0$ for $\omega=0$, and rising to a peak at
$\omega_{p}\approx 14.8$ eV, then rapidly subsiding, although it is still
appreciable at high energy. Note that the behavior of $W_{c}(r,0)$ is
metallic, quickly going to zero beyond the Thomas-Fermi screening radius,
which is likely due to both the finite cluster size in the calculation of the
Green’s functions, which produce only a pseudo-gap in the material, and the
neglect of off-diagonal elements in the interaction kernel $K$ and response
function $\chi^{0}$. However, this unphysical behavior does not seem to affect
$\beta(\omega)$ appreciably. We also note that the satellites in the TiO2
spectral function correspond to the $\approx 13$ eV peak seen in experimental
XPS, which have been interpreted both as charge-transfer excitations,Kas et
al. (2015); Woicik et al. (2020, 2020); Bocquet et al. (1996) as well as
plasmons.Vast et al. (2002) Although distinction between plasmon vs charge
transfer excitations may be a question of semantics,Miyakawa (1968); Županović
et al. (1997); Gunn and Inkson (1979) the behavior of the excitation in TiO2
and the other transition metal compounds is rather different from the plasmons
found in FEG-like materials, where the lattice has little effect on the
plasmon energy. The apparently localized nature of the excitations in TiO2 is
consistent with the definition of charge transfer excitations. Nevertheless,
for the systems investigated here, the excitations at $\omega_{p}$ corresponds
to a zero crossing of the dielectric matrix.

Finally, we compare calculated and experimental metal $2p$ XPS of TiO2, ScF3,
and ScCl3 in Figs. 7 and 8, which show the spectra as a function of energy
relative to the $2p_{3/2}$ main peak. The main $2p_{3/2}$ peak is centered at
$E=0$ eV, and the first peak below that is the $2p_{1/2}$ main peak, seen at
$\approx-6$ eV in TiO2, and $\approx-5$ eV in the Sc halides. Below the
$2p_{1/2}$ peaks, two satellites can be seen, corresponding to a many-body
excitation beyond the main $2p_{3/2}$ or $2p_{1/2}$ peak, lowering the
photoelectron energy by the energy of the excitation $\Delta E$.

Figure 7: Comparison of the RSGF Ti 2p spectral function with the RSRT
approach and XPS experiment.de Boer et al. (1984) Figure 8: Comparison of the
Sc 2p XPS of ScF3 (red) and ScCl3 (blue) calculated with the RSGF approach and
with experiment (dots). The experimental data was digitally reproduced from
Ref. de Boer et al., 1984

Experimentally, the satellite energy is the largest in TiO2 at $\Delta
E\approx 13.4$ eV, followed by ScF3 with $\Delta E\approx 12.3$, then ScCl3 at
$\Delta E\approx 10$ eV. This trend is reproduced by the calculations,
although the satellite energies are slightly overestimated for TiO2 and ScF3,
and slightly underestimated for ScCl3, with $\Delta E\approx 14.1$ eV, $13.4$
eV, and $9.5$ eV respectively. The intensities of the satellites are also
slightly underestimated by the local RSMS theory as well as that of the RSRT
(shown for TiO2). The satellite energies are inversely correlated with the
metal ligand bond-length in the material, with the bond-lengths at $1.95$ Å
for TiO2, $2.03$ Å for ScF3, and $2.51$ Å for ScCl3. This is reminiscent of
XAS, where an expansion of the bond contracts the oscillatory fine-structure
in the spectrum. This likely explains the trend between ScF3 and ScCl3,
although the differing oxidation state of Ti in TiO2 may also play a role in
the satellite energy. Overall, the agreement between the local RSMS theory and
experiment is remarkable given the level of approximations made.

##  V Summary and Conclusions

We have developed an ab initio real-space Green’s function approach for
calculations of intrinsic inelastic losses in x-ray spectra. These losses show
up e.g., as satellites in XPS and reduced fine-structure amplitudes in XAS,
but are often neglected in conventional calculations. The approach is based on
the cumulant Green’s function formalism and a generalization of the Langreth
cumulant to inhomogeneous systems analogous to that in the RSRT, TDDFT
approach. The formalism of the cumulant kernel $\beta(\omega)$ is analogous to
x-ray and optical absorption spectra $\mu(\omega)$, except that the (dipole)
transition operator is replaced by the core-hole potential $V_{c}({\bf r})$,
with monopole selection rules. Within the RPA, three equivalent expressions
for $\beta(\omega)$ are derived, and the link to the GW self-energy is
discussed. A generalized real-space multiple scattering (RSMS) formalism with
a discretized site-coordinate basis is introduced to carry out the
calculations with standard linear algebra. The computational effort further
simplifies with spherical symmetry of the screened core-hole potential and
response function. The calculations can then be carried out on a one-
dimensional radial grid, similar to that in conventional atomic calculations
and the RSGF approach for XAS. Solid-state effects from the extended system
are included in terms of RSGF calculations of back scattering contributions
from atoms beyond the absorption site. The behavior of both $W_{c}({\bf
r},\omega)$ and $\beta(\omega)$ reflect the analytic structure of the inverse
dielectric matrix and the loss function, where the peaks in $\beta(\omega)$
arise from the zeros of ${\rm Re}\,\epsilon(\omega)$. This behavior is
consistent with delocalized quasi-boson excitations, which have been
interpreted both as plasmons and charge-transfer excitations. Moreover, the
low energy $\omega\rightarrow 0$ limit, $\beta(\omega)\approx\alpha\omega$ in
metallic systems, is consistent with the Anderson edge singularity exponent.
Our local RSGF cumulant approximation yields results in semi-quantitative
agreement with the complementary RSRT approach, as well as with experimental
XPS, which validates the local approximations. The discrepancy between RSGF
and RSRT results is likely due to the neglect of the long range behavior in
the bare response function and in the interaction kernel. This behavior is
often not well captured by the truncation of the real-space calculations to
small $R_{max}$. Although this is one of the main limitations of the local
model, corrections can be added using the more general RSMS formalism with a
larger site-radial coordinate basis. Nevertheless, in transition metal oxides
like TiO2, ScF3, and ScCl3, the convergence is found to be fairly rapid,
yielding reasonable accuracy compared with the RSRT approach. In conclusion,
the RSGF approach illustrates the nature and localization of the dynamically
screened core-hole $W_{c}(r,\omega)$ and the cumulant kernel $\beta(\omega)$.
The approach permits ab initio treatments of the spectral function
$A_{c}(\omega)$ and the quasi-particle XAS $\mu_{1}(\omega)$ in the
convolution in Eq. (2), which can be calculated in parallel with the same RSGF
formalism. The approach thereby permits a unified treatment of XPS and XAS
that builds in intrinsic losses. The method can also be applied ex post facto
to include intrinsic losses in other calculations of x-ray spectra. Many
extensions are possible. For example, corrections to the dielectric matrix
from next neighbors may be desirable to improve convergence with respect to
$R_{max}$.Prange et al. (2009); Ankudinov et al. (2005) Although we focused on
result within the RPA, generalization to a real TDDFT interaction kernel
$f_{xc}$ poses no formal or numerical difficulties. Our generalized RSMS
formulation opens the possibility of RSMS approaches to optical spectra beyond
the RPA, e.g., within the TDDFT or the Bethe-Salpeter equation formalism.
Finally, contributions from extrinsic losses and interference terms can be
relevant.Campbell et al. (2002); Cudazzo and Reining (2020); Kas et al. (2015)
Further developments along these lines are reserved for the future.

Acknowledgments – We thank T. Fujikawa, L. Reining, E. Shirley and J. Woicik
for comments, and S. Adkins and S. Thompson for encouragement. This work is
supported by the Theory Institute for Materials and Energy Spectroscopies
(TIMES) at SLAC, which is funded by the U.S. DOE, Office of Basic Energy
Sciences, Division of Materials Sciences and Engineering, under contract DE
AC02-76SF0051.

## References

  * Rehr and Albers (1990) J. J. Rehr and R. C. Albers, Phys. Rev. B 41, 8139 (1990).
  * Zhou et al. (2020) J. S. Zhou, L. Reining, A. Nicolaou, A. Bendounan, K. Ruotsalainen, M. Vanzini, J. Kas, J. Rehr, M. Muntwiler, V. N. Strocov, et al., Proc. Nat. Acad. Sci. 117, 28596 (2020).
  * Langreth (1970) D. Langreth, Phys. Rev. B 1, 471 (1970).
  * Bardyszewski and Hedin (1985) W. Bardyszewski and L. Hedin, Physica Scripta 32, 439 (1985).
  * Hedin (1999) L. Hedin, Journal of Physics: Condensed Matter 11, R489 (1999).
  * Campbell et al. (2002) L. Campbell, L. Hedin, J. J. Rehr, and W. Bardyszewski, Phys. Rev. B 65, 064107 (2002).
  * Fujikawa (2001) T. Fujikawa, J. Synchrotron Rad. 8, 76 (2001).
  * Lee et al. (1999) J. D. Lee, O. Gunnarsson, and L. Hedin, Phys. Rev. B 60, 8034 (1999).
  * Hedin et al. (1998) L. Hedin, J. Michiels, and J. Inglesfield, Phys. Rev. B 58, 15565 (1998).
  * Klevak et al. (2014) E. Klevak, J. J. Kas, and J. J. Rehr, Phys. Rev. B 89, 085123 (2014).
  * Combescot and Nozieres (1971) M. Combescot and P. Nozieres, J. de Physique 32, 913 (1971).
  * Tzavala et al. (2020) M. Tzavala, J. J. Kas, L. Reining, and J. J. Rehr, Phys. Rev. Research 2, 033147 (2020).
  * Liang et al. (2017) Y. Liang, J. Vinson, S. Pemmaraju, W. S. Drisdell, E. L. Shirley, and D. Prendergast, Phys. Rev. Lett. 118, 096402 (2017).
  * Biermann and van Roekeghem (2016) S. Biermann and A. van Roekeghem, J. Electron Spectr. Relat. Phen. 208, 17 (2016).
  * Bagus et al. (2022) P. S. Bagus, C. J. Nelin, C. Brundle, B. V. Crist, N. Lahiri, and K. M. Rosso, Phys. Chem. Chem. Phys. 24, 4562 (2022).
  * Kas et al. (2015) J. J. Kas, F. D. Vila, J. J. Rehr, and S. A. Chambers, Phys. Rev. B 91, 121112 (2015).
  * Woicik et al. (2020) J. Woicik, C. Weiland, C. Jaye, D. Fischer, A. Rumaiz, E. Shirley, J. Kas, and J. Rehr, Phys. Rev. B 101, 245119 (2020).
  * Woicik et al. (2020) J. C. Woicik, C. Weiland, A. K. Rumaiz, M. T. Brumbach, J. M. Ablett, E. L. Shirley, J. J. Kas, and J. J. Rehr, Phys. Rev. B 101, 245105 (2020).
  * Kas et al. (2020) J. Kas, F. Vila, and J. Rehr (2020).
  * Aryasetiawan et al. (1996) F. Aryasetiawan, L. Hedin, and K. Karlsson, Phys. Rev. Lett. 77, 2268 (1996).
  * Guzzo et al. (2011) M. Guzzo, G. Lani, F. Sottile, P. Romaniello, M. Gatti, J. J. Kas, J. J. Rehr, M. G. Silly, F. Sirotti, and L. Reining, Phys. Rev. Lett. 107, 166401 (2011).
  * Landau (1944) L. D. Landau, J. Phys. 8, 201 (1944).
  * Yabana and Bertsch (1996) K. Yabana and G. F. Bertsch, Phys. Rev. B 54, 4484 (1996).
  * Kas et al. (2021) J. Kas, F. Vila, C. Pemmaraju, T. Tan, and J. Rehr, J. Synchrotron Rad. 28 (2021).
  * Almbladh and Hedin (1983) C.-O. Almbladh and L. Hedin, 1b, Hrsg. E. Koch (North Holland, Amsterdam 1983) S pp. 607–904 (1983).
  * Ankudinov et al. (2005) A. L. Ankudinov, Y. Takimoto, and J. J. Rehr, Phys. Rev. B 71, 165110 (2005).
  * Prange et al. (2009) M. P. Prange, J. J. Rehr, G. Rivas, J. J. Kas, and J. W. Lawson, Phys. Rev. B 80, 155110 (2009).
  * Zangwill and Soven. (1980) A. Zangwill and P. Soven., Phys. Rev. A 21, 1561 (1980).
  * Stott and Zaremba (1980) M. J. Stott and E. Zaremba, Phys. Rev. A 21, 12 (1980).
  * Shirley and Martin (1993) E. L. Shirley and R. M. Martin, Phys. Rev. B 47, 15404 (1993).
  * Jr. (1976) J. G. N. Jr., Molecular Physics 31, 1191 (1976).
  * Desclaux (1975) J. P. Desclaux, Comp. Phys. Comm. 9, 31 (1975).
  * Ankudinov (1996) A. L. Ankudinov, Ph.D. thesis, University of Washington (1996).
  * Nozìeres and De Dominicis (1969) P. Nozìeres and C. T. De Dominicis, Phys. Rev. 178, 1097 (1969).
  * Rehr et al. (2020) J. J. Rehr, F. D. Vila, J. J. Kas, N. Y. Hirshberg, K. Kowalski, and B. Peng, J. Chem. Phys. 152, 174113 (2020).
  * Shirley (1972) D. A. Shirley, Phys. Rev. B 5, 4709 (1972).
  * Lundqvist (1967) B. I. Lundqvist, Phys. kondens. Materie 6, 206 (1967).
  * (38) See Supplemental Material at [URL will be inserted by publisher] for derivation of the screened potentias within the plasmon pole model.
  * Citrin et al. (1977) P. H. Citrin, G. K. Wertheim, and Y. Baer, Phys. Rev. B 16, 4256 (1977).
  * Bocquet et al. (1996) A. E. Bocquet, T. Mizokawa, K. Morikawa, A. Fujimori, S. R. Barman, K. Maiti, D. D. Sarma, Y. Tokura, and M. Onoda, Phys. Rev. B 53, 1161 (1996).
  * Vast et al. (2002) N. Vast, L. Reining, V. Olevano, P. Schattschneider, and B. Jouffrey, Phys. Rev. Lett. 88, 037601 (2002).
  * Miyakawa (1968) T. Miyakawa, J. Phys. Soc. Jpn. 24, 768 (1968).
  * Županović et al. (1997) P. Županović, A. Bjeliš, and S. Barišić, Z. Phys. B 101, 397–404 (1997).
  * Gunn and Inkson (1979) J. M. F. Gunn and J. C. Inkson, J. Phys. C: Solid State Phys. 12, 1049 (1979).
  * de Boer et al. (1984) D. K. G. de Boer, C. Haas, and G. A. Sawatzky, Phys. Rev. B 29, 4401 (1984).
  * Cudazzo and Reining (2020) P. Cudazzo and L. Reining, Phys. Rev. Research 2, 012032 (2020).

##  VI Supplementary Material: Plasmon-pole model

$\epsilon^{-1}(q,\omega)=1+\frac{\omega_{p}^{2}}{(\omega+i\delta)^{2}-\omega_{q}^{2}}.$
(46)

$W(q,\omega)=\frac{4\pi}{q^{2}}\epsilon^{-1}(q,\omega).$ (47)

$\displaystyle W(r,\omega)$ $\displaystyle=\int\
\frac{d^{3}q}{(2\pi)^{3}}\frac{4\pi}{q^{2}}\epsilon^{-1}(q,\omega)e^{i\bf{q}\cdot\bf{r}}$
$\displaystyle=\frac{2}{\pi}\int_{0}^{\infty}dq\
\epsilon^{-1}(q,\omega)\frac{\sin(qr)}{qr}$ (48)

Expanding $\sin(qr)=[\exp(iqr)-\exp(-iqr)]/2i$, and noting that
$\epsilon^{-1}$ is an even function of $q$ gives,

$W(r,\omega)=\frac{1}{i\pi}\int_{-\infty}^{\infty}\epsilon^{-1}(q,\omega)\frac{e^{iqr}}{qr}$
(49)

The first term in $W$ coming from the $1$ in $\epsilon^{-1}$ gives the
unscreened Coulomb potential $v(r)=1/r$. The second term $\tilde{W}=W-v$ is
given by,

$\tilde{W}(r,\omega)=\frac{1}{i\pi}\int_{-\infty}^{\infty}dq\
\frac{e^{iqr}}{qr}\frac{\omega_{p}^{2}}{(\omega-\omega_{q}+i\delta)(\omega+\omega_{q}+i\delta)}.$
(50)

For simplicity, we take the dispersion relation
$\omega_{q}=\omega_{p}+q^{2}/2$, which gives an integrand with five poles at
$q=0$, $q=\pm\sqrt{2(\omega-\omega_{p}+i\delta)}$, and $q=\pm
i\sqrt{2(\omega+\omega_{p})}$. The appearance of $\exp(iqr)$ in the requires
that the contour be closed in the upper half plane, so that only the poles
with positive imaginary part, as well as the pole on the real axis contribute.
Assuming positive frequency only, we have poles at $q=0$,
$q_{1}=\sqrt{2(\omega-\omega_{p}+i\delta)}$, and
$q_{2}=i\sqrt{2(\omega+\omega_{p})}$. The $q=0$ pole gives a frequency
dependent contribution $\tilde{W}_{0}$ proportional to the unscreened Coulomb
interaction,

$\tilde{W}_{0}(r,\omega)=\frac{1}{r}\frac{\omega_{p}^{2}}{(\omega+i\delta)^{2}-\omega_{p}^{2}}.$
(51)

In the zero frequency limit, this term cancels the $v(r)$ in $W$, which is
necessary if the screened Coulomb interaction is going to fall off faster than
$1/r$ at large $r$. The pole at $q_{1}$ gives a contribution $\tilde{W}_{1}$
which is real and exponentially decaying for $\omega<\omega_{p}$, and complex
and oscillatory for $\omega>\omega_{p}$, i.e.,

$\tilde{W}_{1}(q,\omega)=-\frac{\omega_{p}^{2}}{\omega
q_{1}}\frac{e^{iq_{1}r}}{q_{1}r}.$ (52)

The above term is also the only term that is complex, so that one can
immediately find the quasiboson excitation spectrum by taking the limit of the
imaginary part as $r\to 0$,

$\displaystyle\beta(\omega)$ $\displaystyle=-\frac{1}{\pi}\lim_{r\to
0}\frac{\omega_{p}^{2}}{\omega
q_{1}}\frac{\sin(q_{1}r)}{q_{1}r}\theta(\omega-\omega_{p})$
$\displaystyle=\frac{1}{\pi}\frac{\omega_{p}^{2}}{\omega
q_{1}}\theta(\omega-\omega_{p}),$ (53)

which matches the result found previously for the same model.Hedin (1999) Note
that the $q=q_{1}$ contribution becomes unphysically large at low frequency,
and is singular at zero frequency, and this is precisely where the $q=q_{2}$
term is important. This gives a contribution $\tilde{W}_{2}$,

$\tilde{W}_{2}(r,\omega)=\frac{\omega_{p}^{2}}{\omega
q_{2}}\frac{e^{iq_{2}r}}{q_{2}r}$ (54)

At low frequency $\omega<\omega_{p}$, $q_{1}$ becomes imaginary, so that
$\tilde{W}_{1}$ and $\tilde{W}_{2}$ tend to cancel,

$\displaystyle\tilde{W}_{1}(r,\omega)$ $\displaystyle+\tilde{W}_{2}(r,\omega)$
$\displaystyle=\frac{\omega_{p}^{2}}{\omega
r}\left[\frac{e^{-\kappa_{1}r}}{\kappa_{1}^{2}}-\frac{e^{-\kappa_{2}r}}{\kappa_{2}^{2}}\right],$

where $\kappa_{1}=\sqrt{2(\omega_{p}-\omega)}$, and
$\kappa_{2}=\sqrt{2(\omega_{p}+\omega)}$, $\kappa_{0}=\sqrt{2\omega_{p}}$,
$\delta\kappa=\omega/\kappa_{0}$. Taking the zero frequency limit of the above
gives the zero frequency limit of $W$, since $v-W_{0}=0$ at zero frequency,

$\lim_{\omega\to
0}W(r,\omega)==e^{-\kappa_{0}r}\left[\frac{\kappa_{0}}{2}+\frac{1}{r}\right].$
(55)

Thus in its entirety, the screened potential is,

$\displaystyle W(r,\omega)$
$\displaystyle=\frac{1}{r}\left[1+\frac{\omega_{p}^{2}}{(\omega+i\delta)^{2}-\omega_{p}^{2}}\right]$
(56)
$\displaystyle-\frac{\omega_{p}^{2}}{\omega}\frac{e^{-\kappa_{2}r}}{\kappa_{2}^{2}r}-\frac{\omega_{p}^{2}}{\omega}\frac{e^{iq_{1}r}}{q_{1}^{2}r}.$
(57)

Finally, we can take real and imaginary parts of this function to diagnose its
behavior at the plasmon frequency.

$\displaystyle{\rm Re}[W(r,\omega)]$
$\displaystyle=\frac{1}{r}\left\\{1+\frac{\omega_{p}^{2}}{(\omega+\omega_{p})(\omega-\omega_{p})}\right.$
$\displaystyle-\frac{\omega_{p}^{2}}{2\omega}\left[\frac{e^{-\kappa_{2}r}}{\omega+\omega_{p}}-\frac{e^{-\kappa_{1}r}}{\omega_{p}-\omega}\theta(\omega_{p}-\omega)\right.$
$\displaystyle\left.\left.+\frac{\cos(q_{1}r)}{\omega-\omega_{p}}\theta(\omega-\omega_{p})\right]\right\\},$
(58) $\displaystyle{\rm Im}[W(r,\omega)]$
$\displaystyle=\frac{\omega_{p}^{2}}{r}\left[-\frac{1}{\omega+\omega_{p}}\delta(\omega-\omega_{p})\right.$
$\displaystyle\left.+\frac{\cos(q_{1}r)}{2\omega}\delta(\omega-\omega_{p})-\frac{\sin(q_{1}r)}{\omega
q_{1}^{2}}\theta(\omega-\omega_{p})\right]$
$\displaystyle=-\frac{\omega_{p}^{2}}{\omega
q_{1}}\frac{\sin(q_{1}r)}{q_{1}r}\theta(\omega-\omega_{p}).$ (59)

Reiterating, $W(r,\omega)$ is real for frequencies less than the plasmon
frequency $\omega<\omega_{p}$. At very low frequencies, the potential falls
off exponentially, although a Yukawa term is also present, with a growing
Coulombic term as frequency increases. At frequencies above the plasmon
frequency, the real part also includes an oscillatory term
$\sim\cos(q_{1}r)/q_{1}r$, and an imaginary part appears, which oscillates as
$\sim\sin(q_{1}r)/q_{1}r$. Near $\omega=\omega_{p}$ is also interesting to
analyze. At the plasmon frequency, the second term in the real part cancels
with the last two terms, giving ${\rm Re}W\sim 1/r$ at large r, while the
imaginary part is singular, with a behavior ${\rm Im}W\sim
1/\sqrt{\omega-\omega_{p}}\theta(\omega-\omega_{p})$. This means that
somewhere near the plasmon frequency, the oscillations in $W$ become very
large, and at high frequency die down, similar to what is seen in Fig. 2.
Finally, in the high frequency limit, all terms except the bare Coulomb
potential scale as $\sim 1/\omega^{2}$, so that we essentially retain the bare
Coulomb potential, as expected.