# Model dielectric functions for fluctuation potential calculations in
electron gas: a critical assessment

Aditi Mandal1,2 ID , Sylvain Tricot2 ID , Rakesh Choubisa1 ID and Didier
Sébilleau2 ID<EMAIL_ADDRESS>
###### Abstract

In this article, we report a critical assessment of dielectric function
calculations in electron gas through the comparison of different modelling
methods. This work is motivated by the fact that the dielectric function is a
key quantity in the multiple scattering description of plasmon features in
various electron-based spectroscopies. Starting from the standard random phase
approximation (RPA) expression, we move on to correlation-augmented RPA, then
damped RPA models. Finally, we study the reconstruction of the dielectric
function from its moments, using the Nevanlinna and memory function
approaches. We find the memory function method to be the most effective, being
highly flexible and customizable.

##  I Introduction

Plasmon effects occur in many electron-driven spectroscopies and are the
signature of the response of the system to the sudden appearance of an extra
charge (a traveling electron, a hole left behind, etc). In core-level
photoemission for instance, plasmons appear as separate peaks at energies
below the core peak. The energy differences with the core peaks are multiple
of the plasmon energy $n\hbar\omega_{p}$. Here, $n$ is the number of plasmon
losses suffered by the photoelectron before escaping the material under study.
Values of $n$ at least up to 6 have been observed by Barman and coworkers
Biswas  _et al._ (2003) in Aluminium. More recent results using 6 keV Hard
X-ray photoelectron spectroscopy (HAXPES) seem to exhibit up to 14 plasmon
peaks Barman (2021). For surface-sensitive photoemission, the surface plasmon
peak can be distinguished from the bulk one as they appear at different
energies.

The information embedded into plasmon peaks has not been much studied so far
in photoemission. Back in 1990, Osterwalder and coworkers showed that they
exhibited photoelectron diffraction-like features just as their core-level
peak Osterwalder  _et al._ (1990). More recently, David, Godet and coworkers
proposed to use these plasmon peaks to extract from their energy distribution
some information on the system’s dielectric function da Silva Santana  _et
al._ (2018); Godet  _et al._ (2020). They termed this new spectroscopy
Photoemission Electron Energy Loss Spectroscopy (PEELS). Similarly, plasmon
structures originating from valence band electrons have also been studied by
Guzzo and coworkers Guzzo  _et al._ (2011, 2012). Therefore, although not used
that much yet as system’s information provider, plasmon structures seem a
promising tool to extract some new information from spectroscopies, and in our
case, from photoemission.

From the theoretical point of view, Lars Hedin and coworkers introduced two
approaches in order to model plasmon features in photoemission and X-ray
absorption: the GW + cumulant expansion and the quasiboson model Hamiltonian
method Guzzo  _et al._ (2011, 2012); Vigil-Fowler  _et al._ (2016). The two
methods have been shown to be formerly equivalent by Vigil-Fowler et al Vigil-
Fowler  _et al._ (2016). The GW + cumulant expansion has been favoured by
Reining and coworkers Guzzo  _et al._ (2011, 2012) to model plasmon structures
originating from valence bands. They showed in particular that neither the DFT
nor the GW were able to model the plasmon. On the other side, Takashi Fujikawa
and co-workers have used Hedin’s quasiboson model in order to incorporate it
into the multiple scattering description of spectroscopies Fujikawa and Arai
(2002, 2005); Kazama  _et al._ (2014). Following Fujikawa’s approach, it can
be shown that, using reasonable approximations, the first plasmon peak cross-
section can be written as the product of the core peak cross section times a
loss function, which is very convenient as standard multiple scattering codes
such as MsSpec Sébilleau  _et al._ (2011) compute the core peak cross-section.
Therefore, the implementation of the plasmon peak into a core-level code only
involves the computation of this loss function. Fujikawa and coworkers have
shown that this function can be expressed as Kazama  _et al._ (2014)

$\int f_{c}(\mathbf{r})V^{i}(\mathbf{r})d\mathbf{r}$ (1)

where $f_{c}(\mathbf{r})$ is a function involving the core hole wave function
and the escaping photoelectron wave function. $V^{i}(\mathbf{r})$ is the so-
called fluctuation potential that describes the excitation of a plasmon (by
either the core hole or the photoelectron). The only unknown in this approach
is the fluctuation potential. Both Hedin and Fujikawa have used for
$V^{i}(\mathbf{r})$ the analytical expressions derived either by Inglesfield
Inglesfield (1983) or by Bechstedt Hedin  _et al._ (1998); Bechstedt  _et al._
(1983) . These descriptions are based on a 3D metallic-type material with a
well identified surface, and for photoelectrons of low kinetic energies,
originating from the vicinity of the surface. Therefore, they completely
exclude lower dimensional systems (quantum dots, quantum wires, quantum wells,
etc) or semiconductors, heterostructures and Dirac systems such as graphene
and related layers or bilayers. Moreover, they are not suited to HAXPES
experiments where the surface is basically overlooked by the escaping
electron.

In order to treat more diverse type of materials and spectroscopies involving
more energetic electrons, we propose an alternative method to the problem of
the description of the fluctuation potential. For this, we go back to the
definition that was given by Hedin and coworkers. In reference Hedin  _et al._
(1998) they define this potential as

$V^{\mathbf{q}}(\mathbf{r})=\left|\frac{V_{c}(\mathbf{q})}{\displaystyle{\frac{\partial\varepsilon(\mathbf{q},\omega)}{\partial\omega}|_{\omega=\omega_{\mathbf{q}}}}}\right|^{1/2}e^{i\mathbf{q}.\mathbf{r}}$
(2)

where $V_{c}(\mathbf{q})$ is the Fourier transform of the Coulomb potential
and $\varepsilon(\mathbf{q},\omega)$ is the dielectric function of the system.
The derivation is taken along the plasmon dispersion $\omega_{\mathbf{q}}$
which can be obtained by solving $Re[\varepsilon(\mathbf{q},\omega)]=0$.

Through the use of the dielectric function, this definition offers us a
convenient, flexible and general way to compute the fluctuation potential
needed to evaluate the cross-section of plasmon peaks. Moreover, being based
on a proper description of the dielectric function, it makes a direct
connection with the PEELS method developed by Godet, David and coworkers da
Silva Santana  _et al._ (2018); Godet  _et al._ (2020), thereby reinforcing
their idea that this spectroscopy could be used ultimately to map somehow the
bulk or surface dielectric functions.

It is not the purpose of this article to implement the calculation of the
plasmon peak. This will be done in a forthcoming study. Rather here, we
dedicate our work to reviewing and documenting different models of dielectric
functions in order to assess their sensitivity to various parameters. In this
view, we remark that some spectroscopies such as photoemission do not resolve
the momentum of the plasmon and therefore will need to integrate over it,
while others such as Electron Energy Loss Spectroscopy (EELS) are plasmon
momentum-sensitive, at least for small values of $\mathbf{q}$.

Ideally, we could obtain the dielectric function from an ab-initio electronic
structure calculation. Unfortunately, most standard ab-initio codes only
compute the dielectric function at $\mathbf{q}=0$, as they essentially use it
for optical spectroscopies. Non zero $\mathbf{q}$ values are generally not
implemented, and when they are, the calculation uses a lot of computing
resources. In order to bypass this problem, we have developed a computer code
called MsSpec-DFM (Dielectric function module), which will be published as a
separate module of MsSpec Sébilleau  _et al._ (2011). This new code can
compute different model dielectric functions for a very large range of
materials. The purpose of this article is to compare different approaches to
analytically model a dielectric function and try to ascertain the most
accurate. All the methods we use are based either on the homogeneous gas
approach or on a Fermi liquid one. This means that the electron system that
responds to the sudden appearance of either a core hole or a photoelectron has
no structure whatsoever. Later, in a forthcoming work, we will try to assess
the effect of the band structure (non homogeneous distribution of
$\mathbf{q}$) and of the crystal structure (interaction with phonons) on the
plasmon description. Work is currently in progress using the Questaal LMTO
code a (63) to model a band structure-sensitive plasmon dispersion Jackson
_et al._ , and preliminary results are very promising.

In section II, we recall the basics of the Hedin-Fujikawa quasiboson multiple
scattering description of the plasmon photoemission peak. Subsection A
introduces the fluctuation potential and subsection B the essentials on
dielectric function formulation. Section III is devoted to the description of
various model dielectric functions encompassing: (i) the simple plasmon pole
approximation, (ii) the RPA approximation, (iii) the correlation-augmented RPA
approximations. As these models do not conserve the number of particles, we
consider in section IV two approximations that ensure this conservation: (i)
The Mermin approximation and (ii) the Hu-O’Connell approximation. More refined
methods that also conserve the momentum and the energy such as the Atwal-
Ashcroft Atwal and Ashcroft (2002) are not tested here because they can be
easily implemented into the family of dielectric functions described in the
section V. Section V deals with a more unusual pathway to build up a
dielectric function. Indeed, as we mentioned before, most standard model
dielectric functions do not conserve the number of particles. In addition, the
RPA does not contain correlation effects and they have to be added externally
through local field corrections (LFC). It can be shown however that the
$\omega^{1}$ moment of the inverse of the dielectric function is nothing else
than an expression of the conservation of the number of particles. Likewise
the $\omega^{3}$ moment can be shown to reflect the electron pair correlations
Pathak and Vashishta (1973) and the $\omega^{5}$ moment the electron triplet
correlations De Raedt and De Raedt (1978). Consequently, if we can reconstruct
the dielectric function, or its inverse, from its moments, we will will have
”built-in” both the conservation of the number of particles and the
correlations. We document here two methods based on this approach: (i) the
Nevanlinna function method and (ii) the memory function method. We show that
these two methods do improve upon the correlation-augmented RPA methods or the
damping methods (Mermin, Hu-O’Connell). In addition, complex energies can be
used that incorporate plasmon damping. Although in this article we will limit
ourselves to the scalar memory function approach, we note that it can be
augmented to a matrix version that will incorporate the conservation of the
momentum and the conservation of the energy. As the memory function is based
on the way a system relaxes, other features, such as the timescales of
different plasmon decay modes, can also be added into the memory matrix method
which makes it the most flexible and general way to model ”analytically” a
dielectric function. Furthermore, it can be shown that most, if not all, other
methods can be shown as particular cases. A comparative discussion is proposed
in section VI.

Note that all throughout this article, we consider the case of Aluminium for
our calculations. This particular choice comes from the fact that Aluminium is
a much studied system usually considered as a test case, or toy model.

##  II Theoretical background

Following Hedin’s quasiboson model Hamiltonian method Bardyszewski and Hedin
(1985); Hedin  _et al._ (1998); Hedin (1999) we add the coupling to the
plasmon field as a perturbation and treat it within MS theory. The quasiboson
model Hamiltonian in our case can be written as,

$\begin{array}[]{ll}H=\sum\limits_{q}\hbar\omega_{q}b_{q}^{\dagger}b_{q}\,+\,\sum\limits_{\mathbf{k}}\epsilon_{\mathbf{k}}c_{\mathbf{k}}^{\dagger}c_{\mathbf{k}}\\\
+\,\sum\limits_{q\mathbf{k}\mathbf{k^{\prime}}}\Bigg{[}V^{q}_{\mathbf{k}\mathbf{k^{\prime}}}b_{q}^{\dagger}+(V^{q}_{\mathbf{k}\mathbf{k^{\prime}}})^{*}b_{q}\Bigg{]}c_{\mathbf{k}}^{\dagger}c_{\mathbf{k^{\prime}}}-\sum\limits_{q}V^{q}_{cc}(b_{q}+b_{q}^{\dagger})\end{array}$
(3)

where the first term describes the free boson field, the second term
corresponds to the ejected photoelectron, the third term corresponds to the
interaction between the photoelectron and the boson field. The last term
represents the core hole-boson coupling. In this approach, the interaction
between an external charge and the electrons in the system is described
through the fluctuation potential $V^{\mathbf{q}}$ given by Eq. 2 where
$\omega_{\mathbf{q}}$ is the frequency of a plasmon of potential energy
$\hbar\omega_{p}$ and momentum $\hbar\mathbf{q}$. The momentum vector of a
plasmon is given by $\mathbf{q}=\mathbf{k}_{in}-\mathbf{k}_{sc}$ where
$\mathbf{k}_{in}$ is the momentum of the electron before the plasmon loss and
$\mathbf{k}_{sc}$ is the the momentum of the electron after the plasmon loss.

In their respective works, Hedin and coworkers and Fujikawa and coworkers
relied on two analytical fluctuation potentials developed separately by
Bechstedt Bechstedt  _et al._ (1983) and Inglesfield Wikborg and Inglesfield
(1977), which we will describe later for the sake of completeness.

In the Hedin-Fujikawa formalism, the intensity of the first plasmon peak is
given by Kazama  _et al._ (2014),

$I^{1}(\mathbf{k},\hbar\omega_{\mathbf{q}},\hbar\omega_{i})=I^{0}(\mathbf{k},\hbar\omega_{\mathbf{q}})\frac{\alpha(\hbar\omega_{i})}{\hbar\omega_{i}},$
(4)

where $\mathbf{k}$ is the momentum of the detected electron. The ”no-loss”
core-peak cross-section is Kazama  _et al._ (2014)

$\displaystyle
I^{0}(\mathbf{k},\hbar\omega_{\mathbf{q}})=\,4\pi^{2}\alpha_{FS}\,\frac{\hbar}{m^{2}\omega_{\mathbf{q}}}\mid<\tilde{\phi}_{\mathbf{k}}\mid\Delta\mid\,c>^{0}\mid^{2}$
$\displaystyle\,\times\exp\big{[}-\int_{0}^{+\infty}\frac{\alpha(\epsilon)}{\epsilon}d\epsilon\big{]}$
(5)

The term before the multiplication sign is the usual core-level cross-section
which is computed by standard MS codes such as MsSpec.

Within the Hedin-Fujikawa formalism, the loss function
$\alpha(\epsilon)/\epsilon$ can be expressed as:

$\frac{\alpha(\epsilon)}{\epsilon}=\sum_{i}\mid\int
f_{c}(\mathbf{r})V^{i}(\mathbf{r})d\mathbf{r}\mid^{2}\delta(\epsilon-\hbar\omega_{i})$
(6)

where $f_{c}(\mathbf{r})$ is a well-defined core-state related function and
$V^{i}(\mathbf{r})$ is the fluctuation potential corresponding to the
excitation of a plasmon of energy $\hbar\omega_{i}$. So, for quantitative
modelling of plasmon features in spectroscopies it is imperative to have a
good fluctuation potential as it is the only unknown in (6). Both Hedin and
Fujikawa have used for $V^{i}(\mathbf{r})$ the analytical expressions derived
either by Inglesfield Inglesfield (1983) or by Bechstedt Hedin  _et al._
(1998); Bechstedt  _et al._ (1983). These descriptions are based on a 3D
metallic-type material with a well identified surface, and for photoelectrons
of low kinetic energies, originating from the vicinity of the surface.
Therefore, they completely exclude lower dimensional systems (quantum dots,
quantum wires, quantum wells…) or semiconductors, heterostructures and Dirac
systems such as graphene and related layers or bilayers. Moreover, they are
not suited to HAXPES experiments where the surface is basically overlooked by
the escaping electron. In order to treat more diverse type of materials and
spectroscopies involving more energetic electrons, we propose an alternative
method to the problem of the description of the fluctuation potential.

###  II.1 Fluctuation potentials

The fluctuation potential describes the coupling between the electron and the
bosons. It is not the purpose of this article to go into the details of the
analytical derivation of the ones that can be found in the literature. From
equation (2), we see that it can be factorized as

$V^{\mathbf{q}}(\mathbf{r})=V^{\mathbf{q}}\,e^{i\mathbf{q}.\mathbf{r}}$ (7)

The different models of analytical fluctuation potentials available in the
literature have been obtained within linear response theory. In addition, in
the case of Inglesfield and Bechstedt, the choice of a semi-infinite electron
gas implies that the exponential in (7) reduces essentially to the form
$\exp\left[\mathbf{q}_{\parallel}\;z\right]$.

####  II.1.1 Plasmon-pole

The plasmon-pole fluctuation potential was originally derived by Lundqvist
[19]. Its value is

$\left\\{\begin{array}[]{lll}V^{\mathbf{q}}(z)&=&V^{\mathbf{q}}_{\mathrm{PP}}\;e^{i\mathbf{q}.\mathbf{r}}\\\
&&\\\
V^{\mathbf{q}}_{\mathrm{PP}}&=&\sqrt{V_{c}(\mathbf{q})\dfrac{\omega_{p}^{2}}{\omega_{\mathbf{q}}}\dfrac{1}{2V}}\end{array}\right.$

where $V^{\mathbf{q}}_{\mathrm{PP}}$ is the electron plasmon coupling constant
and $V$ is the system volume. For the surface case, we replace the 3D Fourier
transformation of the bare Coulomb potential $V_{c}(\mathbf{q})$ by the 2D
transformation and take into account that $\mathbf{q}$ is parallel to the
surface.

####  II.1.2 Inglesfield

Inglesfield introduced a more realistic analytical potential, using the fact
that the bulk modes in a semi-infinite system are standing waves, i.e., phase-
shifted cosines modified at the surface Wikborg and Inglesfield (1977). The
Inglesfield bulk plasmon potential is given as:

$\begin{array}[]{ll}V^{\mathbf{q}}(z)=&\displaystyle{V^{\mathbf{q}}_{\mathrm{IN}}[\cos(-q_{z}z+\phi_{\mathbf{q}})-e^{q_{\parallel}z}\cos\phi_{\mathbf{q}}]\theta(z)},\,\text{inside},\\\
\\\ \,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,=&0,\,\,\text{outside}.\end{array}$ (8)

where

$\begin{array}[]{lll}\displaystyle{V^{\mathbf{q}}_{\mathrm{IN}}}&=&\displaystyle{\sqrt{\frac{\omega_{p}^{2}}{\omega_{\mathbf{q}}}\frac{1}{V}V_{c}(\mathbf{q})}},\\\
&&\\\
\displaystyle{\phi_{\mathbf{q}}}&=&\displaystyle{\tan^{-1}(\frac{q_{\parallel}}{q_{z}})}\end{array}$
(9)

Similarly, Inglesfield surface plasmon potential is given as:

$\begin{array}[]{lll}\displaystyle{V^{\mathbf{q}}}(z)&=&\displaystyle{V^{\mathbf{q}}_{\mathrm{IN}}\;e^{q_{\parallel}z}},\\\
&&\\\
\displaystyle{V^{\mathbf{q}}_{\mathrm{IN}}}&=&\displaystyle{\sqrt{V_{c}^{2D}\frac{\omega_{p}}{A\sqrt{8}}}}\end{array}$
(10)

The Coulomb potential Fourier transforms being

$V_{c}=\,V_{c}^{2D}\,=\dfrac{e^{2}}{2\epsilon_{0}\sqrt{q^{2}+k_{s}^{2}}}$ (11)

where $k_{s}$ is a screening momentum.

####  II.1.3 Bechstedt

Bechstedt and co-workers derived an expression for the screened potential $W$.
It was recast in terms of the fluctuation potential by Hedin and coworkers
Hedin  _et al._ (1998):

$\begin{array}[]{lll}V^{\mathbf{q}}(z)=&\displaystyle{N_{b}\,[e^{q_{\parallel}z}\theta(z)+\\{(2+C_{1}+C_{3})\,\cos(q_{z}z)}\\\
\\\ &\displaystyle{-C_{2}\,\sin(q_{z}z)-(1+C_{1})e^{q_{\parallel}z}-}\\\ \\\
&\displaystyle{C_{3}\exp(\sqrt{\omega_{p}+\omega_{\mathbf{q}}+q_{\parallel}^{2}}z)\\}\,\theta(-z)]}\end{array}$
(12)

where the coefficients are given by

$\begin{array}[]{lll}C_{1}=\dfrac{\omega_{p}^{2}}{\omega_{\mathbf{q}}^{2}-\omega_{p}^{2}},\\\
C_{2}=\dfrac{q_{\parallel}\omega_{p}^{2}}{2\omega_{\mathbf{q}}(\omega_{\mathbf{q}}-\omega_{p})q_{z}},\\\
C_{3}=\dfrac{q_{\parallel}\omega_{p}^{2}}{2\omega_{\mathbf{q}}(\omega_{\mathbf{q}}-\omega_{p})q_{z}\sqrt{\omega_{p}+\omega_{\mathbf{q}}+q_{\parallel}^{2}}}\end{array}$
(13)

In the surface case,

$V^{\mathbf{q}}(z)=N_{s}\left[e^{-q_{\parallel}z}\theta(z)+a(q_{\parallel},z,\omega)\theta(-z)\right]$
(14)

where $a(q_{\parallel},z,\omega)$ is related to the bulk dielectric function
through Bechstedt  _et al._ (1983)

$a(q_{\parallel},z,\omega)=\frac{2q_{\parallel}}{\pi}\int_{0}^{\infty}\frac{\cos(q_{z}z)}{q^{2}\varepsilon(\mathbf{q},\omega)}dq_{z}$
(15)

where $\varepsilon(\mathbf{q},\omega)$ can be computed within the simple
plasmon-pole model.

These three potentials are the standard analytical fluctuation potentials
available in literature. Their space variations are represented in Fig.1, in
the Aluminium case. However, due to the fact that they are surface-related
(and derived for metals), they cannot be used to describe spectra derived for
reduced-symmetry systems such as quantum dots, quantum wells, graphene and
other 2D materials, Dirac materials, etc. In addition, they are not suited to
high-energy spectroscopies such as HAXPES where the surface can be safely
ignored.

In order to overcome the limitations of these potentials, we will go back to
Hedin’s definition of the fluctuation potential (equation (2)), which
expresses this potential as a function of the dielectric function. This will
allow us a much more flexible approach where the true dimensionality and
structure of the material can be properly taken into account.

Figure 1: $r$-variations of the standard analytical fluctuation potentials.

###  II.2 Dielectric function: General Background

The dielectric function of a system describes the response of this system to
an external perturbation. It is determined by the properties of the system and
its interaction with the perturbing object. Several related quantities can be
found in the literature, together with their relationship with various
spectroscopies. for instance, it is well-known that the cross-section of the
EELS is related to the loss function. Therefore, in the following subsection,
we will introduce the different quantities of interest for spectroscopies and
that are related to the dielectric function.

We start by decomposing the dielectric function into its real and imaginary
parts

$\varepsilon(\mathbf{q},\omega)=\varepsilon_{1}(\mathbf{q},\omega)+i\varepsilon_{2}(\mathbf{q},\omega)$
(16)

Then, the loss function $L(\mathbf{q},\omega)$ is related to dielectric
function of the solid through,

$L(\mathbf{q},\omega)\,\propto\,\Im\left[\varepsilon^{-1}(\mathbf{q},\omega)\right]$
(17)

Similarly, the dynamical structure factor can be expressed as Götze and Wölfle
(1972)

$S(\mathbf{q},\omega)=\Im\left[\frac{-1}{\varepsilon(\mathbf{q},\omega)}\right]=\,\frac{\varepsilon_{2}(\mathbf{q},\omega)}{\mid\varepsilon(\mathbf{q},\omega)\mid^{2}}$
(18)

It describes the spectrum of excitations in the systems as a function of the
momentum transfer $\mathbf{q}$ and the energy transfer.

Likewise, the susceptibility, or density-density response function, can be
defined as

$\chi(\mathbf{q},\omega)=\dfrac{1}{V_{c}(q)}\left[\dfrac{1}{\varepsilon(\mathbf{q},\omega)}-1\right]$
(19)

With these tools, we can access many different ways to model the dielectric
function, and compute cross-sections.

####  II.2.1 The simplest case: Plasmon pole

The plasmon pole dielectric function describes the response of the system
entirely in terms of collective modes. The dielectric function is just the
analytic continuation of a simple pole and is given by Aulbur  _et al._ (2000)

$\varepsilon(\mathbf{q},\omega)=1-\frac{\omega_{p}^{2}}{\omega^{2}+\omega_{p}^{2}-\omega_{q}^{2}}$
(20)

Its real part and imaginary part are represented in figure 2 for the case of
Aluminium.

The plasmon dispersion band corresponds to the upper band of
$\Re[\varepsilon(\mathbf{q},\omega)]=0$.

##  III RPA and beyond

The plasmon pole approximation’s main drawback, as can be seen from the left-
hand part of figure 2, is that it does not incorporate any damping mechanism.
Therefore, we have a plasmon (the upper white band in the real part of
$\varepsilon(\mathbf{q},\omega)$) which, once created, never decays. This is
clearly not physical. Historically, the first expression of the dielectric
function to include a damping mechanism is the Random Phase approximation
(RPA) originally derived by Lindhard Lindhard  _et al._ (1954). This
approximation is based on the description of the delocalized electrons system
as a _homogeneous and non-interacting electron gas_. The real and imaginary
parts of the RPA dielectric function are represented in the middle figure of
figure 2. Now, we have a clear damping inside the so-called Landau region
which is delimited by the two parabola $\omega_{+}$ and $\omega_{-}$. This
damping comes from the decay of the plasmon into an electron-hole pair. There
is however no damping mechanism built-in outside the Landau region, which
means that there, the plasmon will ”live” forever. In addition, being built on
the assumption that the electron gas is non-interacting, the RPA does not
allow for correlation effects.

Correlation effects can be introduced into the RPA through the so-called
_local field corrections_ (LFC) $G(\mathbf{q},\omega)$. These corrections are
related to the exchange and correlation term of the Density Functional Theory
(DFT). For a given correction term, the correlation-augmented dielectric
function can be written as

$\varepsilon(\mathbf{q},\omega)=\frac{V_{c}(\mathbf{q})\Pi^{RPA}(\mathbf{q},\omega)}{1+V_{c}(\mathbf{q})G(\mathbf{q},\omega)\Pi^{RPA}(\mathbf{q},\omega)}$
(21)

where $G(\mathbf{q},\omega)$ is the dynamical local field correction. The RPA
dielectric function is given by

$\varepsilon^{RPA}(\mathbf{q},\omega)=1-V_{c}(\mathbf{q})\Pi^{RPA}(\mathbf{q},\omega)$
(22)

where $\Pi^{RPA}(\mathbf{q},\omega)$ is the RPA polarization.

Figure 2: Real (above) and imaginary (below) part of dielectric function
$\epsilon(\mathbf{q},\omega)$ for (left) plasmon pole, (center) RPA, (right)
RPA + UTI1 respectively.

Very few models of _dynamical_ local field corrections exist in the
literature, while a lot of attention has been devoted to _static_ corrections
$G(\mathbf{q})$. Therefore, we will restrict here to the latter.

In order to investigate the effect of such corrections on the behaviour of the
dielectric function, we consider here three different types of static
corrections, the Hubbard model (HUBB), which takes only into account exchange
effect, the Pathak-Vashista correction (PVHF) Pathak and Vashishta (1973), and
the Utsumi-Ichimaru (UTI1) one Utsumi and Ichimaru (1980). They can be
expressed respectively as

$\left\\{\begin{array}[]{lll}G(\mathbf{q})_{\mathrm{HUBB}}&=&\displaystyle{\frac{1}{2}\frac{x^{2}}{1+x^{2}}}\\\
&&\\\
G(\mathbf{q})_{\mathrm{PVHF}}&=&\displaystyle{\frac{1}{\omega_{p}^{2}}\;J(q)}\\\
&&\\\
G(\mathbf{q})_{\mathrm{UTI1}}&=&\displaystyle{\frac{3(4-x^{2})\,(28+5x^{2})}{16x}}\times\log\left|\dfrac{2+q}{2-q}\right|\Bigg{]}\end{array}\right.$
(23)

where we have used the notation $x=q/k_{F}$.

Here, the Kugler function $J(q)$ is given by Kugler (1975)

$J(q)=\,\frac{e^{2}}{m\pi}\;\int_{0}^{+\infty}\,k^{2}\;\left[S(k)-1\right]\;J(k,q)\;dk$
(24)

where $J(k,q)$ is given by

$J(k,q)=\,\frac{5}{6}-\frac{k^{2}}{2q^{2}}+\frac{q}{4k}\left(\frac{k^{2}}{q^{2}}-1\right)^{2}\;\ln\left|\frac{k+q}{k-q}\right|$
(25)

$S(q)$ is the static structure factor. Here, we approximate it by its Hartree-
Fock value

$\begin{array}[]{llll}S_{HF}(q)&=&\displaystyle{\frac{3}{4}\;\frac{q}{k_{F}}-\frac{1}{16}\;\left(\frac{q}{k_{F}}\right)^{3}}&\text{for
}q<2k_{F}\\\ &&&\\\ &=&1&\text{for }q>2k_{F}\end{array}$ (26)

Figure 3: Comparison of collective excitations dispersion lines for RPA and
some RPA + local field corrections. The plasmon dispersion corresponds to the
upper band.

The right-hand plots of Fig.2 corresponds to the real part (up) and imaginary
part (down) of the UTI1-enhanced RPA dielectric function. We see clearly a
change in the dispersion of the collective excitations (white band), with
respect to the RPA case. The comparison between this dispersion is even
clearer in figure 3 where we give a 2D representation of the collective
excitations bands for the RPA, RPA + Hubbard and RPA + Pathak-Vashista static
local field corrections. Nevertheless, even if we now incorporate correlation
effects into the dielectric function, we see clearly that, despite some minor
changes in the imaginary part of the dielectric function (see Fig.2), we are
still unable to properly describe plasmon damping outside the Landau regime.
In principle, we could go to dynamic LFC, but in the following, we will rather
explore completely other ways to describe the dielectric function with built-
in damping.

##  IV Damping-based dielectric functions

RPA alone has a certain number of limitations: (i) it fails to conserve the
number of particles, (ii) it does not contain correlation effects (they have
to be added externally through local field corrections) and (iii) it does not
incorporate plasmon damping outside the Landau regime. In order to cure (i),
Mermin Mermin (1970) extended the Lindhard dielectric function in the
relaxation-time approximation, where essentially the collisions relax the
electronic density matrix not to its uniform equilibrium value, but to a local
equilibrium density matrix. The Mermin Nersisyan and Das (2004) dielectric
function thus can be written as:

$\varepsilon(\mathbf{q},\omega)\,=\,1+\frac{\displaystyle{\left(1+\frac{i}{\omega\tau}\right)}\;\left[\varepsilon^{0}(\mathbf{q},\omega+i/\tau)-1\right]}{1+\displaystyle{\frac{i}{\omega\tau}\;\left[\frac{\varepsilon^{0}(\mathbf{q},\omega+i/\tau)-1}{\varepsilon^{0}(\mathbf{q},0)-1}\right]}}$
(27)

where $\varepsilon^{0}(\mathbf{q},\omega)$ is the RPA dielectric function and
$\bar{\omega}=\omega+i/\tau$ is a complex frequency incorporating damping
through the relaxation time $\tau$. Mermin’s interest here is in the long
wavelength limit behaviour of the electron gas where the focus is on obtaining
a modified Lindhard dielectric function which reduces to the correct classical
behaviour in the $q\longrightarrow 0$ limit.

Following the same idea, Hu and O’ Connell Hu and O’Connell (1989) generalized
the Lindhard dielectric function to include fluctuation effects arising from
electron-electron and electron-impurity interactions. They also studied
Friedel oscillations as an example of the application of their generalized
version of the Lindhard function and observed that these oscillations are
damped due to the inclusion of fluctuation effects. Due to the complexity of
the Hu-O’Connell equations, which are expressed in terms of the diffusion
coefficient $D$, we do not reproduce them here but refer instead to their
article Hu and O’Connell (1989). As displayed in Fig. 4, we now observe a
damping in the imaginary part of both the Mermin and the Hu-O’Connell
dielectric function in the non Landau region. This damping can be regarded as
the signature of a built-in lifetime for the plasmon. However, being based
upon the standard RPA model, they lack a proper description of correlations.
Refer to 4 for the real and the imaginary part of both the Mermin and the
Hu-O’Connell dielectric function.

Figure 4: Real and imaginary part of the (left) Mermin dielectric function
with a relaxation time of $0.5$ fs, (right) Hu-O’Connell dielectric function
with a diffusion coefficient $D=2.910^{-5}$ m2 s-1

These two approaches are encouraging, as we do need a proper damping of the
plasmon in order to correctly describe the relaxation of the system. This
damping is clearly absent in the standard plasmon pole and RPA methods, even
when the latter is augmented with various _static_ local field corrections in
order to incorporate electron-electron correlations.

##  V Dielectric function: Alternative approach

An alternative family of dielectric functions can be obtained through a
reconstruction from the first few moments. Two independent approaches can be
found in the literature, the Nevanlinna approach Arkhipov  _et al._ (2007) and
the memory function approach Jindal  _et al._ (1977). The main advantage of
these two approaches is that conservation of the number of particles and
correlation are built-in. We note that for larger values of $\omega$,
$\Im\left[\varepsilon^{-1}(\mathbf{q},\omega)\right]$ tends quickly to zero,
so that we have the expansion

$\varepsilon^{-1}(\mathbf{q},\omega)\,=\,1+\sum_{n=1}^{+\infty}\,\frac{\left\langle{\omega^{2n-1}}\right\rangle_{L}}{\omega^{2n}}$
(28)

As mentioned in the introduction, it is well-known that the RPA dielectric
function does not satisfy the compressibility sum rule and the frequency
moment sum rules Pines (2018). This realization is the starting point of the
two above-mentioned approaches which we outline now.

###  V.1 The Nevanlinna function method

The Nevanlinna approach is based on the moments of the loss function. The
standard loss function described in Eq. 17 is related to _Nevanlinna loss
function_

$\mathcal{L}(\mathbf{q},\omega)\,=\,\frac{L(\mathbf{q},\omega)}{\omega}$ (29)

The corresponding moments are defined as Arkhipov  _et al._ (2007)

$C_{n}\;=\,\frac{1}{\pi}\;\int_{-\infty}^{+\infty}\,\omega^{n-1}\;\Im\left[\frac{-1}{\varepsilon(\mathbf{q},\omega)}\right]\;d\omega\,=\,\frac{1}{\pi}\;\left\langle{\omega^{n-1}}\right\rangle_{\mathcal{L}}$
(30)

A dimensional analysis shows that moments $C_{2}(q)$ has the dimension of a
squared frequency, and likewise for every ratio $C_{n+2}(q)/C_{n}(q)$.
Therefore, we introduce the notation

$\omega_{n}^{2}(q)\,=\,\frac{C_{2n}(q)}{C_{2n-2}(q)}$ (31)

The $\omega_{n}(q)$ have the dimension of a frequency and are the
characteristic frequencies in the Nevanlinna function approach. The same
characteristic frequencies will appear in the memory function method. The
connection between the different types of moments can be made using, for the
$T=0$ K case,

$\begin{array}[]{lcrl}L\left(\mathbf{q},\omega\right)\,=\,\Im\left[\displaystyle{\frac{-1}{\varepsilon(\mathbf{q},\omega)}}\right]&=\,V_{C}(q)\;\displaystyle{\frac{\pi\bar{n}}{\hbar}}\;S(\mathbf{q},\omega)\\\
&&&\\\ &=\,-V_{C}(q)\;\Im\left[\chi(\mathbf{q},\omega)\right]\end{array}$ (32)

where we find

$V_{C}(q)\;\frac{\pi\bar{n}}{\hbar}\,=\,\frac{\pi}{2}\;\frac{\omega_{p}^{2}}{\omega_{\mathbf{q}}}$
(33)

implying

$\left\langle{\omega^{n}}\right\rangle_{S}\,=\,\frac{2}{\pi}\;\frac{\omega_{\mathbf{q}}}{\omega_{p}^{2}}\;\left\langle{\omega^{n}}\right\rangle\,=\,-\frac{\hbar}{\bar{n}\pi}\;\left\langle\omega^{n}\right\rangle\chi$
(34)

with $\bar{n}$ being the number density. Now, the frequency moments of the
loss function can be written as Arkhipov  _et al._ (2013a),

$\left\langle\omega^{2n-1}\right\rangle_{L}\,=\,2\;\int_{0}^{+\infty}\,\omega^{2n-1}\;\Im\left[\frac{-1}{\varepsilon(\mathbf{q},\omega)}\right]\;d\omega$
(35)

So, then we have

$\begin{array}[]{llll}\left\langle\omega^{-1}\right\rangle_{L}&=&\pi\;\left[1-\displaystyle{\frac{1}{\epsilon(q)}}\right]\text{
compressibility sum rule}\\\ &&&\\\
\left\langle\omega\right\rangle_{L}&=&\pi\;\omega_{p}^{2}\;\;\;f\text{-sum
rule}\\\ &&&\\\
\left\langle\omega^{3}\right\rangle_{L}&=&\pi\;\displaystyle{\omega_{p}^{2}\;\left[\omega_{\mathbf{q}}^{2}+4\;\omega_{\mathbf{q}}\frac{\left\langle{t}\right\rangle}{\hbar}+\omega_{p}^{2}+J(q)\right]}\\\
&\end{array}$ (36)

$J(q)$ is the Kugler function defined in equation (24) and
$\left\langle{t}\right\rangle$ is the average kinetic energy per electron.

$n=1$ ($f$-sum rule) ensures the conservation of the number of particles and
$n=3$ ensures a proper account of two-body correlation effects. The structure
factor $S(\mathbf{q},\omega)$ for a non-zero temperature system can be
obtained from the loss function through the fluctuation-dissipation theorem

$\Im\left[\frac{-1}{\varepsilon(\mathbf{q},\omega)}\right]\,=\,\frac{\pi}{\hbar}\;V(\mathbf{q})\;\left[1-\exp\left(-\displaystyle{\frac{\hbar\omega}{k_{B}T}}\right)\right]\;S(\mathbf{q},\omega)$
(37)

For low energies, or large temperatures, this can be approximated by,

$\Im\left[\frac{-1}{\varepsilon(\mathbf{q},\omega)}\right]\,\approx\,\frac{\pi}{k_{B}T}\;V(\mathbf{q})\;\omega\;S(\mathbf{q},\omega)$
(38)

The Nevalinna method is mathematically involved. It relies on the mathematical
solution of the so-called _non-canonical solution of the Hamburger moment
problem_. It is not the purpose of the present article to go into the details
of the method. We refer the readers interested by this approach to the review
article by Tkachenko Tkachenko (2018). The main point is that it involves a
reconstruction of the dielectric function in terms of characteristic
frequencies $\omega_{n}$ that are functions of the first moments of the loss
function. Building up on the mathematical theorems, one arrives at the
equation

$\varepsilon^{-1}(\mathbf{q},\bar{\omega})\,=\,1+\frac{\omega_{p}^{2}\;\left(\bar{\omega}+Q_{2}(\mathbf{q},\bar{\omega})\right)}{\bar{\omega}\;(\bar{\omega}^{2}-\omega_{2}^{2})+Q_{2}(\mathbf{q},\bar{\omega})\;(\bar{\omega}^{2}-\omega_{1}^{2})}\,\Im\left[\bar{\omega}\right]\geq
0$ (39)

for the 3-moment expression. Here, $\omega_{p}$ is the plasmon frequency,
$Q_{2}(\mathbf{q},\omega)$ is the _unknown_ Nevalinna function and
$\bar{\omega}$ is the complex frequency whose imaginary part represents the
plasmon damping. We note that $Q_{2}(\mathbf{q},\bar{\omega})$ has the
dimension of a frequency.

Figure 5: Real and imaginary part of the 3-moment Nevanlinna dielectric
function within the static approximation STA3.

In an electron liquid, the Nevanlinna function plays the role of the
_dynamical_ local field correction $G(\mathbf{q},\omega)$. More precisely, we
have Dubovtsev (2019)

$G({\mathbf{q}},\bar{\omega})\,=\,1+\frac{1}{\left[\varepsilon_{\mathrm{RPA}}(\mathbf{q},\bar{\omega})-1\right]}+\frac{\bar{\omega}^{2}}{\omega_{p}^{2}}-\frac{\bar{\omega}\omega_{2}^{2}+\omega_{1}^{2}\;Q_{2}(\mathbf{q},\bar{\omega})}{\omega_{p}^{2}\left(\bar{\omega}+Q_{2}(\mathbf{q},\bar{\omega})\right)}$
(40)

Nevalinna functions must fulfill a number of mathematical properties including
having a Riesz-Herglotz representation Arkhipov  _et al._ (2013b); Tkachenko
(2018). Therefore, finding relevant functions is a complex mathematical
problem that has led to several papers in the literature Tkachenko (2018);
Dubovtsev (2019); Arkhipov  _et al._ (2014); Adamjan  _et al._ (1989).

Here, we make the choice Adamjan  _et al._ (1989)

$Q(q,\bar{\omega})\,=\,i\;\frac{\pi}{2}\;\varepsilon(q)\;\left(\varepsilon(q)-1\right)\;(qa_{0})\;\omega_{\mathbf{q}}\;\left[\frac{\omega_{2}^{2}}{\omega_{1}^{2}}-1\right]$
(41)

which is well suited to our needs.

From Fig.5, we observe that the correlations are properly taken into account
(there is a dramatic change into the plasmon dispersion with respect to RPA-
based methods) and in the imaginary part of the Nevanlinna dielectric function
we see some damping of the plasmon outside the Landau region, especially along
the plasmon dispersion line.

###  V.2 The memory function method

The derivation of the memory function is quite involved so that we describe
here only the different steps necessary in order to arrive to the final
result. For a detailed derivation, we refer the reader to some review articles
Yoshida and Takeno (1989); Umberto Balucani and Tognetti (2003); Kumari and
Singh (2020). This approach was pioneered by Zwanzig Zwanzig (1961) and Mori
Mori (1965) who built the memory function framework upon Kubo’s non-
equilibrium statistical physics. The starting point is to divide the
variables’ space into two sub-spaces, that of slowly moving ones, the so-
called conservative quantities, that will affect the macroscopic experimental
signal, and that of fast-moving ones (non-conservative), which are not
expected to impact the experimental signals. The memory function method
includes the effect of the fast variables into the dynamics of the slow
variables. In order to do this, we work with time-correlation functions or
response functions. The connection to the dielectric function is made through
the realization that (i) the dynamical structure factor $S(\mathbf{q},\omega)$
is the time-Fourier transform of the _classical_ density auto-correlation
function, (ii) the density-density response function
$\chi_{nn}(\mathbf{q},\omega)$ is related to the dielectric function.

The next step is to show that within this variable space partitioning, auto-
correlation functions satisfy an intregro-differential generalized Langevin
equation (GLE)

$\frac{dC_{A}(t)}{dt}\,=\,-\int_{0}^{t}\;K_{A}(t-s)C_{A}(s)\;ds\quad\text{with}\quad
C_{A}(0)\,=\,1$ (42)

where $C_{A}(t)$ is the auto-correlation function of interest, and the unknown
function $K_{A}(t)$ is the so-called memory function which keeps track of what
happened to the system before the present time $t$. The response function
$\chi_{AA}(t)$ is related to the auto-correlation function through

$\chi_{AA}(t)=\frac{dC_{A}(t)}{dt}\;\left(A|A\right)$

where $(A,B)$ is the so-called _Kubo scalar product_.

The usual Langevin equation of Brownian motion corresponds to the Markovian
choice $K_{A}(t)=1/\tau\,\delta(t)$ where at $t$, the dynamics of the system
does not depend on the previous states of the system. The integral term
containing $K_{A}(t)$ describes the influence of the fast moving variables,
presumably out of reach to the experiment, on the (conservative) variable $A$
of interest.

It can be demonstrated that the memory function is also an auto-correlation
function, and as such it satisfies a hierarchy of similar GLE

$\frac{d{K}_{n}(t)}{dt}\,=\,-\int_{0}^{t}\;{K}_{n+1}(t-s)\;{K}_{n}(s)\;ds$
(43)

where we take $C_{A}(t)=K_{0}(t)$ and $K_{A}(t)=K_{1}(t)$.

The next step is to take the Laplace transform of the integro-differential GLE
in order to change it into an easily solvable equation. This leads to the
continued fraction expansion

$\bar{C}_{A}(z)\,=\,\frac{K_{0}(0)}{z+\displaystyle{\frac{K_{1}(0)}{z+\displaystyle{\frac{K_{2}(0)}{z+\displaystyle{\frac{K_{3}(0)}{\ddots}}}}}}}$
(44)

Truncating to order 3 (i.e. expressing the result in terms of 3 moments) and
going to the response function representation, gives

$\bar{\chi}_{A}^{(3)}(z)\,=\,\displaystyle{-\frac{\Delta_{0}^{2}\Delta_{1}^{2}\left[z+\bar{K}_{3}(z)\right]}{z^{3}+z^{2}\bar{K}_{3}(z)+z\left(\Delta_{1}^{2}+\Delta_{2}^{2}\right)+\Delta_{1}^{2}\;\bar{K}_{3}(z)}}$
(45)

where the $\Delta_{n}$ are function of the moments
$\left\langle\omega^{2n-1}\right\rangle_{\chi}$ and $K_{n}(0)=\Delta_{n}^{2}$.
Here, $\bar{K}_{3}(z)$ is the Laplace transform of the memory function
$K_{3}(t)$.

As the moments for $n$ = 0,1,2 are well-known (and to some extent, making use
of some approximations, $n=3$), the only unknown is the memory function
$\bar{K}_{3}(z)$.

Then, we have, taking $A$ as the (constant) particle density $\bar{n}$

$\varepsilon^{-1}(\mathbf{q},\omega)\,=\,1+\frac{\omega_{p}^{2}\;\left[\omega+i\bar{K}_{3}(-i\omega)\right]}{\omega\;\left[\omega^{2}-\omega_{2}^{2}\right]+i\bar{K}_{3}(-i\omega)\;\left[\omega^{2}-\omega_{1}^{2}\right]}$
(46)

We note that under the exchange $Q_{2}(\omega)\leadsto
i\bar{K}_{3}(-i\omega)$, this equation transforms into the Nevanlinna
equation. The difference is that we have now a physically motivated function
$K_{3}$ for which many approximations exist and that is related to the
relaxation of the system, and hence in the various decay mechanisms, which
could, in principle, be studied experimentally.

Figure 6: Real and imaginary part of the 3-moment memory function dielectric
function within the Cole-Cole approximation (COCO) for a value of 0.5 fs.

By definition, through the use of $\left\langle\omega^{1}\right\rangle_{\chi}$
and $\left\langle\omega^{3}\right\rangle_{\chi}$, Eq. 46 contains the
conservation of the number of particles and a proper treatment of pair
correlation effects.

If we replace the auto-correlation function $A$ in Eq. 44 by a vector
containing the density, the momentum and the energy, the GLE transforms into a
matrix equation Mountain (1976); Hansen and McDonald (2005). Applying the same
approach, we can derive an expression of the dielectric function that
conserves the three previous quantities. This makes the memory function/memory
matrix approach a very powerful and flexible tool to build up a physically
meaningful dielectric function. Moreover, in case of need, two or more memory
functions describing relaxations on different timescales can be used within
Götze’s mode-coupling formalism Reichman and Charbonneau (2005); Janssen
(2018); Götze and Sjögren (1995). This allows to incorporate different decay
channels into the modelling of the dielectric function.

Here we consider the Cole-Cole (COCO) type memory function expression Khamzin
_et al._ (2012) for the complex dielectric function describing the relaxation
process occurring in the system. This _classic_ empirical model can be
expressed as Khamzin  _et al._ (2012)

$K_{CC}(t)=\tau^{2}\dfrac{1}{\Gamma(\alpha-1)}(\frac{t}{\tau})^{\alpha-2}$
(47)

with the power-law exponent $\alpha$ being $0<\alpha\leq 1$. Here, $\Gamma$ is
the Gamma function. In the following, we choose $\alpha=0.5$ and a relaxation
time of $\tau=0.5$ fs, to be consistent with the results of the Mermin
dielectric function.

Figure 7: The Cole-Cole (COCO) memory function in the time domain.

In Fig. 7 we observe the behaviour of memory function in the Cole-Cole
approximation at different time-scale. The corresponding dielectric function
is represented in Fig. 6. Here, we observe that the imaginary part of the
dielectric function also shows some damping of the plasmon outside the Landau
region.

Figure 8: Comparison of the collective excitations dispersion bands for
different modelling of the dielectric function.

Figure 9: Respectively (a) RPA, (b) Mermin, (c) Nevanlinna (STA3) and (d)
memory function (COCO) structure factors.

##  VI Discussion

We have explicited in the previous sections different ways to analytically
compute model dielectric functions, starting from the simple plasmon pole and
RPA, and ending up with more involved approaches involving a reconstruction in
terms of the first moments. We recall that our aim in this work is to find a
simple (if possible) and flexible method that would describe plasmon
dispersion as accurately as possible. Indeed, our ultimate goal being to
calculate plasmons fluctuation potentials in the most general possible case,
we need for this to follow Hedin’s definition (2)

$V^{\mathbf{q}}(\mathbf{r})=\left|\frac{V_{c}(\mathbf{q})}{\displaystyle{\frac{\partial\varepsilon(\mathbf{q},\omega)}{\partial\omega}}|_{\omega=\omega_{\mathbf{q}}}}\right|^{1/2}e^{i\mathbf{q}.\mathbf{r}}$

For low values of $\mathbf{q}$, the fluctuation potential will be dominated by
the Coulomb part. But when $\mathbf{q}$ increases (in photoemission, for
instance, we will have to integrate over $\mathbf{q}$ as the plasmon momentum
is not detected), the features of the first derivative of the dielectric
function, taken along the plasmon dispersion will become more and more
important. From this point of view, we see that two features of the dielectric
function could become important: 1) the shape of the plasmon dispersion and 2)
a correct description of the plasmon damping. The former corresponds to the
upper band of the zeros of the real part of $\varepsilon(\mathbf{q},\omega)$
while the latter is embedded into the imaginary part.

All the calculations were done for a value of the dimensionless Wigner-Seitz
radius $r_{s}/a_{0}$=2.079, which corresponds to Aluminium ($r_{s}$= 1.01 Å).

From figure 2, we see that the plasmon pole does not exhibit any damping and
that the plasmon never decays. This means that in principle, we should have to
integrate $q$ to infinity. For RPA and correlation-augmented RPA, the
integration over $q$ will be limited to a $q_{max}$ above which the plasmon
completely decays into an electron-hole pair. This $q_{max}$ is very close to
the intersection of the plasmon dispersion with the Landau region, which means
that in practice no damping will be taken into account as clear from the
imaginary part plots.

The next type of dielectric functions we have considered are the number of
electron-conserving Mermin and Hu-O’Connell ones (figure 4). In the former, we
used a realistic value of the relaxation time $\tau\,=\,0.5$ fs. We note that
this value is about one order of magnitude larger that the approximation often
made in the literature $\tau\,\approx\,1/\omega_{p}$ which gives here
$\tau\,\approx\,0.04$ fs. This relaxation time approach involves the use of a
complex frequency $\bar{\omega}=\omega+i/\tau$ which automatically introduces
damping outside the Landau continuum. The related diffusion coefficient method
(Hu-O’Connell) also gives such a damping. From this point of view, they are
very appealing methods. Unfortunately, being based on the RPA, they lack a
proper description of correlation effects. In the case of Mermin’s approach,
this could be remedied for by replacing in equation (27) the RPA
$\varepsilon^{0}(\mathbf{q},\omega+i/\tau)$ by a correlation-augmented one
such as RPA + UTI1.

The final approach we have documented is that based on the 3-moment
reconstruction of the dielectric function. This method makes use of the
moments $\left\langle\omega^{-1}\right\rangle$,
$\left\langle\omega^{1}\right\rangle$ and
$\left\langle\omega^{3}\right\rangle$ of the loss function. Here,
$\left\langle\omega^{1}\right\rangle$ ensures the conservation of the number
of electrons while $\left\langle\omega^{3}\right\rangle$ incorporates a proper
treatment of pair correlations. As can be seen from figure 5 and figure 6, the
more exact treatment of the correlation changes considerably the shape of the
plasmon dispersion with respect to the RPA one (dotted line) or even with
respect to RPA + UTI1 (figure 2). In addition, this type of approach does give
already some damping of the plasmon outside the Landau regime. But it is
noteworthy that both calculations have been done _without explicit damping_ ,
i.e. using a real frequency $\omega$ in equation (39) and (46), in contrast to
the Mermin calculation that was incorporating an explicit damping in the
frequency. Therefore, these two methods have even more flexibility that we
have used so far.

A quantity related to the dielectric function that can be measured
experimentally is the dynamical structure factor $S(\mathbf{q},\omega)$.
Lemell et al Lemell  _et al._ (2015) for instance, present in their Fig.2
results for the dynamical structure factor derived from optical data of Mg.
There, we seen clearly the plasmon, starting from $q=0$. In Fig. 9, we present
the structure factors for Aluminium computed from our dielectric function
using the (a) RPA, (b) Mermin, (c) Nevalinna-STA3 and (d) memory function-
COCO. These plots show that RPA and RPA based dielectric functions should be
ruled out as they are not able to describe properly the plasmon for low-$q$
values, while the other methods give a much better representation that does
look like the experimental results of Lemell et al Lemell  _et al._ (2015).

Now we come to the more general discussion of the choice of a suitable
dielectric function model. We have not explored all the possibilities of the
models, but we can already outline some important points and perspectives. The
Mermin approach is limited by being RPA-based, although it could be in
principle expanded through the use of LFC. In addition, following the scheme
devised by Götze Götze (1978) for the computation of the susceptibility, it
can be demonstrated that Mermin’s dielectric function is a particular case of
Götze’s scheme when neglecting correlations and using the Markovian form for
the memory function $K(\omega)=1/\tau$.

Nevanlinna functions are hard to find, but as mentioned in the memory function
subsection, they are mathematically related to the memory functions through
$Q_{2}(\omega)\leadsto i\bar{K}_{3}(-i\omega)$. Moreover, following equation
(40), we see that Nevalinna functions and hence memory functions are strongly
related to _dynamical_ LFC.

In practice, these considerations mean that most, if not all, the model
dielectric functions can be viewed as particular cases of the memory function
approach. Another advantage of this method is that it is highly flexible and
customizable. We have shown here only the basic features, limiting ourselves
to a simple memory function and to a scalar generalized Langevin equation,
using the constant density $\bar{n}$ as the only variable of interest. But, as
mentioned in the introduction, by using a vector composed of 1) the number
density $n_{\mathbf{k}}(t)$, 2) the longitudinal current density
$j^{\ell}_{\mathbf{k}}(t)$ and 3) the energy density $e_{\mathbf{k}}(t)$, we
can build a matrix GLE containing a $3\times 3$ memory matrix Mountain (1976);
Hansen and McDonald (2005). Following the same scheme as presented here, we
can obtain the expression of a dielectric function that has embedded: (i) the
conservation of the number of particle, (ii) the conservation of the momentum,
(iii) the conservation of the energy, (iv) pair correlations. Moreover, we
believe that the Atwal-Ashcroft approach Atwal and Ashcroft (2002) will be
found to be a particular case of the $3\times 3$ memory matrix method.

If need be, this method can be further augmented by building on 4 moments, as
it is known that $\left\langle\omega^{5}\right\rangle$ incorporates three-body
correlations De Raedt and De Raedt (1978). The exact value of this term is
very complicated, but reasonable approximations exist that allow to compute
it. Work is in progress to incorporate it into the MsSpec-DFM computer code.

The memory function approach we have used here relies on a single memory
function, which means that it includes the relaxation of the system on a
single timescale. In classical systems, such as the relaxation of a fluid,
molecular dynamics seems to favour two-relaxation-time laws Snook (2006). The
memory matrix method we have outlined in the previous paragraph makes use of
several memory functions. But even in the simpler method discussed here, we
can accommodate several relaxation processes operating on different
timescales, which in our case would describe different decay channels of the
plasmon. This is the so-called mode-coupling framework Reichman and
Charbonneau (2005); Janssen (2018); Götze and Sjögren (1995) developed by
Götze.

##  VII Conclusion

In this work, we have tested different electron gas dielectric functions in
order to ultimately model fluctuation potentials. These fluctuations
potentials are the key quantity in the multiple scattering description of
plasmon features in spectroscopies such as photoemission or EELS.

We have presented three families of dielectric functions, RPA-based, damped
RPA-based and 3-moment reconstructed. We have shown that the simple RPA-based
ones are not suited to our needs. Furthermore, as the damped type of
dielectric functions can be shown to be a particular case of Götze’s memory
function scheme, we have come to the conclusion that the memory function
approach, to which the Nevalinna one is strongly related, has all the features
needed for a precise and accurate modelling. In addition, it can be improved
by passing to the matrix approach and by incorporating the different plasmon
decay channels and their timescale through the mode-coupling framework.

This approach is still based on the homogeneous electron gas model. Therefore,
no band structure and no crystal structure is involved. We are currently
investigating the influence of the band structure by performing ab initio
calculations of $\varepsilon(\mathbf{q},\omega)$ with the Questaal code a
(63); Jackson  _et al._ . Preliminary results point towards a minor influence,
at least in the case of Aluminium which has a simple bandstructure. Another
line of research which we are currently pursuing in order to further improve
our description is to couple this electron dielectric function to a phonon
dielectric function. We are hopeful all these improvements will lead to a
fast, efficient and flexible approach to the modelling of the dielectric
function necessary to incorporate plasmons features in to the multiple
scattering description of spectroscopies.

The MsSpec-DFM code, built during this work, will be published soon as a
separate module of the MsSpec code Sébilleau  _et al._ (2011). In addition to
3D dielectric functions, as exposed here, it will also contain the modelling
of other dimensionalities, including graphene-type and multilayer structures.

## Acknowledgements

A. M. is indebted to Rennes Métropole for providing her with two 6-month
grants as a visiting scientist.

## References

  * Biswas  _et al._ (2003) C. Biswas, A. Shukla, S. Banik, V. Ahire, and S. Barman, Physical Review B 67, 165416 (2003).
  * Barman (2021) S. Barman, “private communication,” (2021).
  * Osterwalder  _et al._ (1990) J. Osterwalder, T. Greber, S. Hüfner, and L. Schlapbach, Physical Review B 41, 12495 (1990).
  * da Silva Santana  _et al._ (2018) V. M. da Silva Santana, D. David, J. S. de Almeida, and C. Godet, Brazilian Journal of Physics 48, 215 (2018).
  * Godet  _et al._ (2020) C. Godet, D. David, V. M. da Silva Santana, J. S. de Almeida, and D. Sébilleau,  _Recent Advances in Thin Films_ , pp 181-210 (Springer, Singapore, 2020) Chap. Photoelectron Energy Loss Spectroscopy: A Versatile Tool for Material Science.
  * 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, Physical Review Letters 107, 166401 (2011).
  * Guzzo  _et al._ (2012) M. Guzzo, J. J. Kas, F. Sottile, M. G. Silly, F. Sirotti, J. J. Rehr, and L. Reining, The European Physical Journal B 85, 1 (2012).
  * Vigil-Fowler  _et al._ (2016) D. Vigil-Fowler, S. G. Louie, and J. Lischner, Physical Review B 93, 235446 (2016).
  * Fujikawa and Arai (2002) T. Fujikawa and H. Arai, Journal of Electron Spectroscopy and Related Phenomena 123, 19 (2002).
  * Fujikawa and Arai (2005) T. Fujikawa and H. Arai, Journal of Electron Spectroscopy and Related Phenomena 149, 61 (2005).
  * Kazama  _et al._ (2014) M. Kazama, H. Shinotsuka, Y. Ohori, K. Niki, T. Fujikawa, and L. Kövér, Physical Review B 89, 045110 (2014).
  * Sébilleau  _et al._ (2011) D. Sébilleau, C. Natoli, G. M. Gavaza, H. Zhao, F. Da Pieve, and K. Hatada, Computer Physics Communications 182, 2567 (2011).
  * Inglesfield (1983) J. Inglesfield, Journal of Physics C: Solid State Physics 16, 403 (1983).
  * Hedin  _et al._ (1998) L. Hedin, J. Michiels, and J. Inglesfield, Physical Review B 58, 15565 (1998).
  * Bechstedt  _et al._ (1983) F. Bechstedt, R. Enderlein, and D. Reichardt, Physica Status Solidi (b) 117, 261 (1983).
  * a (63) “www.questaal.org,” .
  * (17) J. Jackson, L. Petit, A. Mandal, and D. Sébilleau, “Work in progress,” .
  * Atwal and Ashcroft (2002) G. Atwal and N. Ashcroft, Physical Review B 65, 115109 (2002).
  * Pathak and Vashishta (1973) K. Pathak and P. Vashishta, Physical Review B 7, 3649 (1973).
  * De Raedt and De Raedt (1978) H. De Raedt and B. De Raedt, Physical Review B 18, 2039 (1978).
  * 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).
  * Wikborg and Inglesfield (1977) E. Wikborg and J. Inglesfield, Physica Scripta 15, 37 (1977).
  * Götze and Wölfle (1972) W. Götze and P. Wölfle, Physical Review B 6, 1226 (1972).
  * Aulbur  _et al._ (2000) W. G. Aulbur, L. Jönsson, and J. W. Wilkins,  _Quasiparticle Calculations in Solids_ , Solid State Physics, Vol. 54 (Academic Press, 2000) pp. 1–218.
  * Lindhard  _et al._ (1954) J. Lindhard, M. Scharff, and H. Schiøtt, Fys. Medd 28, l954 (1954).
  * Utsumi and Ichimaru (1980) K. Utsumi and S. Ichimaru, Physical Review B 22, 1522 (1980).
  * Kugler (1975) A. Kugler, Journal of Statistical Physics 12, 35 (1975).
  * Mermin (1970) N. D. Mermin, Physical Review B 1, 2362 (1970).
  * Nersisyan and Das (2004) H. B. Nersisyan and A. K. Das, Physical Review E 69, 046404 (2004).
  * Hu and O’Connell (1989) G. Hu and R. O’Connell, Physical Review B 40, 3600 (1989).
  * Arkhipov  _et al._ (2007) Y. V. Arkhipov, A. Askaruly, D. Ballester, A. Davletov, G. Meirkanova, and I. Tkachenko, Physical Review E 76, 026403 (2007).
  * Jindal  _et al._ (1977) V. Jindal, H. Singh, and K. Pathak, Physical Review B 15, 252 (1977).
  * Pines (2018) D. Pines,  _Elementary Excitations in Solids_ (CRC Press, 2018).
  * Arkhipov  _et al._ (2013a) Y. V. Arkhipov, A. Ashikbayeva, A. Askaruly, A. Davletov, and I. Tkachenko, EPL (Europhysics Letters) 104, 35003 (2013a).
  * Tkachenko (2018) I. Tkachenko, Physical Sciences and Technology 5, 16 (2018).
  * Dubovtsev (2019) D. Y. Dubovtsev, PhD Thesis, Universidad Politécnica de Valencia (2019).
  * Arkhipov  _et al._ (2013b) Y. V. Arkhipov, A. Ashikbayeva, A. Askaruly, A. Davletov, and I. Tkachenko, Contributions to Plasma Physics 53, 375 (2013b).
  * Arkhipov  _et al._ (2014) Y. V. Arkhipov, A. B. Ashikbayeva, A. Askaruly, A. E. Davletov, and I. M. Tkachenko, Physical Review E 90, 053102 (2014).
  * Adamjan  _et al._ (1989) S. Adamjan, T. Meyer, and I. Tkachenko, Contributions to Plasma Physics 29, 373 (1989).
  * Yoshida and Takeno (1989) F. Yoshida and S. Takeno, Physics Reports 173, 301 (1989).
  * Umberto Balucani and Tognetti (2003) H. H. L. Umberto Balucani and V. Tognetti, Physics Reports 373, 409 (2003).
  * Kumari and Singh (2020) K. Kumari and N. Singh, European Journal of Physics 41, 053001 (2020).
  * Zwanzig (1961) R. Zwanzig, Physical Review 124, 983 (1961).
  * Mori (1965) H. Mori, Progress in Theoretical Physics 33, 423– (1965).
  * Mountain (1976) R. D. Mountain,  _Advances in Molecular Relaxation Processes_ , 9 (Elsevier, Amsterdam, 1976) Chap. Generalized Hydrodynamics, pp. 225–291.
  * Hansen and McDonald (2005) J.-P. Hansen and I. R. McDonald,  _Theory of Simple Liquids_ (Academic Press, 2005).
  * Reichman and Charbonneau (2005) D. R. Reichman and P. Charbonneau, Journal of Statistical Mechanics: Theory and Experiment 5, P05013 (2005).
  * Janssen (2018) L. M. C. Janssen, Frontiers in Physics 6 (2018).
  * Götze and Sjögren (1995) W. Götze and L. Sjögren, Transport Theory and Statistical Physics 24, 801 (1995).
  * Khamzin  _et al._ (2012) A. A. Khamzin, R. Nigmatullin, and I. I. Popov, Theoretical and Mathematical Physics 173, 1604 (2012).
  * Lemell  _et al._ (2015) C. Lemell, S. Neppl, G. Wachter, K. Tőkési, R. Ernstorfer, P. Feulner, R. Kienberger, and J. Burgdörfer, Physical Review B 91, 241101 (2015).
  * Götze (1978) W. Götze, Solid State Communications 27, 1393 (1978).
  * Snook (2006) I. Snook,  _The Langevin and Generalised Langevin Approach to the Dynamics of Atomic, Polymeric and Colloidal Systems_ (Elsevier, 2006).