# Electron-phonon interaction and phonons in 2d doped semiconductors

Francesco Macheda, Thibault Sohier, Paolo Barone, Francesco Mauri

► To cite this version:

Francesco Macheda, Thibault Sohier, Paolo Barone, Francesco Mauri. Electron-phonon interaction and phonons in 2d doped semiconductors. Physical Review B, 2023, pp.094308. 10.1103/Phys-RevB.107.094308. hal-03917628

HAL Id: hal-03917628

https://hal.science/hal-03917628

Submitted on 10 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Electron-phonon interaction and phonons in 2d doped semiconductors

Francesco Macheda,¹ Thibault Sohier,² Paolo Barone,³ and Francesco Mauri⁴,¹

¹Istituto Italiano di Tecnologia, Graphene Labs, Via Morego 30, I-16163 Genova, Italy

²Laboratoire Charles Coulomb (L2C), Université de Montpellier, CNRS, Montpellier, France

³CNR-SPIN, Area della Ricerca di Tor Vergata,
Via del Fosso del Cavaliere 100, I-00133 Rome, Italy

⁴Dipartimento di Fisica, Università di Roma La Sapienza, Roma, Italy

Electron-phonon interaction and phonon frequencies of doped polar semiconductors are sensitive to long-range Coulomb forces and can be strongly affected by screening effects of free carriers, the latter changing significantly when approaching the two-dimensional limit. We tackle this problem within a linear-response dielectric-matrix formalism, where screening effects can be properly taken into account by generalized effective charge functions and the inverse scalar dielectric function, allowing for controlled approximations in relevant limits. We propose complementary computational methods to evaluate from first principles both effective charges – encompassing all multipolar components beyond dynamical dipoles and quadrupoles – and the static dielectric function of doped two-dimensional semiconductors, and provide analytical expressions for the long-range part of the dynamical matrix and the electron-phonon interaction in the long-wavelength limit. As a representative example, we apply our approach to study the impact of doping in disproportionated graphene, showing that optical Fröhlich and acoustic piezoelectric couplings, as well as the slope of optical longitudinal modes, are strongly reduced, with a potential impact on the electronic/intrinsic scattering rates and related transport properties.

## I. INTRODUCTION

The electron-phonon interaction (EPI) is one of the most thoroughly studied topics in solid state physics [1–5] due to the fundamental role it plays in the determination of a variety of physical properties. The prediction and interpretation of e.g., transport experiments [6–16], excited carriers relaxation [17–19] and superconductivity [3, 20], relies on the accurate calculation of the EPI from first-principles, which has become possible in recent years thanks to the development of density functional theory (DFT) [21, 22], density functional perturbation theory (DFPT) [23] and Wannier interpolation technique [24–26], as well as swift progress of computational infrastructures [27–29].

In insulators and undoped semiconductors, long-range Coulomb interactions arise from the charge polarization and lead to non-analytic contributions to phonons and EPIs, including the well-known splitting of longitudinal (LO) and transverse (TO) optical modes [30] as well as the Fröhlich [31] and piezoelectric electron-phonon interactions [32]. In the absence of free carriers and in the long-wavelength limit, the form of such non-analytic contributions is known exactly, enabling a precise evaluation of related effects [23, 32–35]. Due to the strong screening provided by free carriers in partially filled bands, those long-ranged Coulomb interactions are expected to vanish in metals. In doped semiconductors, on the contrary, Coulomb-mediated interactions are only partially screened by the small fraction of added charge carriers. In this intermediate situation between insulator and metal, the long-wavelength behaviour of phonons and EPI may be significantly altered, as recently shown for 3d doped semiconductors [36]. Going beyond phenomenological corrections of the screened quantities

[37], the formalism developed in [36] proposes a clean separation between bare quantities and their screening, thus enabling the evaluation of doping and temperature effects on both of them. In regimes where doping and temperature effects on the bare quantities are negligible, this allows for an efficient yet precise evaluation and interpolation of phonons and EPIs at a given temperature and doping via the modification of screening only.

Given that electrostatically-doped 2D materials are at the heart of the quest for high-efficiency electronic devices [38, 39], the same rigorous treatment of screened Coulomb interactions in quasi-2d systems would be highly beneficial. However, it is now well established that dimensionality alters the dielectric screening behaviour in 2d undoped semiconductors and insulators, and consequently any phonon or EPI involving long-range Coulomb interactions [40]. The implementation of DFPT in 2d boundary conditions [41] has allowed to study those effects; most notably, the LO-TO splitting breaks down at zone center and increases linearly with momentum [42, 43], while the Fröhlich EPI stays finite [44] instead of diverging as the inverse of momentum, as in 3D [32]. Yet, as for 3d materials, the doping effects on the 2d long-range Coulomb interaction are virtually always neglected [45–51], and occasionally included via approximate models [52]. Nonetheless, a direct computation of both phonons and EPIs in the presence of electrostatic doping, as recently implemented in DFPT [41], showed that metallic screening causes a vanishing of the linear-in-momentum LO-TO splitting specific of 2d semiconductors [43], a sizeable effect that might be detected by momentum-resolved electron energy loss spectroscopy in a transmission electron microscope (TEM-EELS) [53]. The intrinsic dimensionality reduction combined with careful Brillouin-zone sampling allowed to
compute transport properties from the Boltzmann equa-
tion formalism in highly-doped 2d semiconductors with-
out resorting to Wannier interpolation methods [54, 55].
However, such a procedure is intrinsically costly and lim-
ited due to the fine sampling of electronic states needed
to account for a small Fermi surfaces, especially at small
dopings. Thus, as typical DFPT calculations scale like
the cube of the number of atoms in the simulation cell,
small dopings, large systems or multiple doping and tem-
perature conditions of the same system remain out of
reach. This latter shortfall in particular prevents any sys-
tematic study of a material transport properties, which
is of the outmost importance for the development of elec-
tronic devices. Recently, a procedure based on the for-
malism of Ref. [40] has been proposed to deal with some
of the above shortfalls [56], but neglecting the effects of
doping and temperature on bare interactions and on the
local-fields components of the response.

Motivated by the above considerations, we extend the
theoretical framework recently outlined for 3d systems
[36] in order to deal with quasi-2d doped semiconduc-
tors. Using a static linear-response dielectric matrix for-
mulation, we introduce screened and unscreened effec-
tive charges. Along with the inverse scalar dielectric
function, these are used to derive general expressions for
the long-range Coulomb contributions to the dynamical
matrix and EPI. At vanishing doping, they reduce to
well-established formulas, both for 3d and 2d systems.
This allows for controlled approximations of screening
effects in appropriate doping and temperature regimes,
whose range of validity can be assessed within the gen-
eral theoretical framework. At variance with the 3d case,
the presence of a non-periodic direction complicates the
electrostatic problem and enforces the dependence of all
the response functions on the out-of-plane component z.
Nonetheless, a great deal of simplification comes when
the wavelength of the interaction is much larger or much
smaller than the typical out-of-plane thickness of the ma-
terial *t*. We denote these regimes as the thin and thick
limits, characterized respectively by a wavevector depen-
dence of the in-plane Coulomb kernel being of the form
$q^{-1}$ and $q^{-2}$. In such regimes, for materials displaying
in-plane mirror symmetry, one can integrate out the out-
of-plane variable via layer-averages to a good approxi-
mation. The resulting long-range components (LRC) of
both the dynamical matrix and EPI can then be con-
nected, as in the 3D case, to a properly modified ver-
sion of the well-known phenomenological theory of Born
and Huang [30], involving in-plane effective charges and
the inverse scalar dielectric function $\epsilon^{-1}$. We further
discuss their range of validity for the two-dimensional
case. The layer-averaging procedure is naturally appro-
priate for those single layers materials where the physical
observables are mainly determined by in-plane electro-
statics. Nonetheless, even when out-of-plane perturba-
tions are important (e.g. when considering remote cou-
plings with surrounding materials), our approach allows
for a simple extraction of phonon perturbations that are

macroscopically unscreened both from the in-plane di-
electric response and from the presence of periodically
repeated images. This makes it suitable to be integrated
in frameworks that account for remote screening in het-
erostructures such as the one proposed in Ref. [57], which
proposes a treatment of the out-of-plane electrostatics a-
posteriori.

Operatively, we propose a fast and precise technique
based on first-principles calculations and Wannier in-
terpolation that is grounded in our general theoretical
formulation and supported by the computation of both
screened and unscreened charge responses. Crucially, in
controlled regimes, accurate quantities can be obtained
from *ab-initio* calculations performed only for the un-
doped setup, while still accounting for the most relevant
doping and temperature effects beyond the state-of-the-
art, and with a great reduction of the computational
workload with respect to brute-force methods. We val-
idate our findings in disproportionated graphene, i.e., a
particular realization of gapped graphene which can be
found in the presence of substrates causing a symmetry-
breaking modulation of potential [58–60] and that has
been recently proposed to host strong polar responses
[61]. Our approach shows that there exists a small doping
regime where the LRCs can be described using the effec-
tive charges value of the undoped setup at zero tempera-
ture and an RPA expression for the screening. The exis-
tence of this regime is expected for any two-dimensional
material. In the strong doping regime, instead, this sim-
plification doesn’t occur and the effective charges are
mostly determined by the appearance of intraband terms
in the electronic polarizability; one shall then resort to
the *ab-initio* calculation of the macroscopically screened
and unscreened effective charge functions in the doped
setup. As a practical example of the implications of our
developments, we show that electronic lifetimes can be
strongly affected by the presence of free carriers, with
likewise critical implications for physical observables. In
particular, the reduction of the electronic lifetimes with
doping opens the way to the engineering of new upper
limits to the carriers mobility in 2d materials via a fine-
tuned choice of the doping and temperature regimes, with
important consequences for the design of efficient field-
effect transistors.

The paper is organized as follows: in Sec. II we present
the general theoretical framework of our approach, start-
ing from the introduction of the response functions and
of the effective charges for two-dimensional materials to
arrive to the expressions of the LRCs of the dynamical
matrix and of the EPI; in Sec. III we discuss the imple-
mentation of our theoretical considerations into an oper-
ative computational approach; in Sec. IV we apply our
developments to the case of disproportionated graphene;
finally, in Sec. V we draw our conclusions; the appen-
dices are devoted instead to the treatment of technical
details which integrate the derivations of the main text.
## II. THEORY

### A. Framework and main formulae

We are interested in the effect of free-carriers on the dynamical matrix C and electron-phonon interaction g as a function of the chemical potential and of the temperature T. To reconcile with the doped semiconductor literature, we prefer to use the free carrier concentration n as a variable, which evaluates to zero in the case of an undoped semiconductor, rather than the chemical potential which is more appropriate for metals.

Our aim is to evaluate the expressions of C and g on fine grids in reciprocal space at any given (n, T), having at disposal with a reasonable effort only their value for the undoped setup at zero Kelvin. Following the same logic as Ref. [36], we exploit the usual separation of short and long-range components of C and g

$$C = C^S + C^L, \quad g = g^S + g^L, \qquad (1)$$

and focus on the description of the long-range components since, as we will show, they are the only ones strongly affected by (n, T). Their description, within well-defined and controlled regimes, requires the knowledge of two main ingredients: i) macroscopically unscreened effective charges and ii) the macroscopic inverse dielectric function. In such regimes, doping and temperature act on the long-range components mostly through the macroscopic inverse dielectric function, which can be approximately expressed as a function of quantities computed in absence of free-carriers. As a result, Eq. 1 may be described using ab-initio techniques only for (n = 0, T = 0). Outside these regimes, we can still compute long-range components at reduced cost on a few selected line, exploiting crystal symmetries.

The theoretical approach that we use to deduce the various components of Eq. 1 is based on the static dielectric matrix formulation of the linear response problem for quasi-2d materials. Despite the generality of the method and the possibility to obtain exact results, we find that employing the static RPA approximation [62] entails a vast simplification of the derivations at a formal level, with at the same time the possibility to (carefully) generalize the conclusions even to the presence of exchange and correlation terms. As done in Ref. [36], we therefore employ the RPA approximation in order to derive the main theoretical results. When computing numerical results, we will reintroduce exchange-correlation effects and defer their theoretical discussion to appendices. In order to simplify the theoretical treatment, we also restrict our arguments to the in-plane electrostatic of quasi-2d systems with in-plane mirror symmetry. Such systems do not mix in-plane and out-of-plane responses at the first order in the expansion of the EPI, while higher orders (such as the piezoelectric coupling) may retain information regarding the out-of-plane components [51]. Nonetheless, the influence of such terms on e.g. the mobility comes mostly

from the coupling with acoustic modes [51], which may be screened statically even at low doping since the plasma and phonon frequencies are comparable. We therefore prefer to keep only the leading order (yet accurate) description of the EPI and discuss its modification is presence of doping.

As anticipated, within these approximations we can obtain relations for the in-plane long-range components of the dynamical matrix and of the EPI involving effective charge functions and the two-dimensional macroscopic inverse dielectric function. In particular, we define the macroscopically unscreened/screened effective charges $\bar{Z}_{s,\alpha}(\mathbf{q})/Z_{s,\alpha}(\mathbf{q})$ starting from the expression of the total charge change arising in a crystal following a collective displacement of the atoms of type s along the Cartesian direction $\alpha$ modulated via an in-plane wavevector $\mathbf{q}$

$$\delta\rho_{s,\alpha}^{\text{tot}}(\mathbf{q}, n, T) := -i \frac{eq}{At} Z_{s,\alpha}(\mathbf{q}, n, T), \qquad (2)$$

$$Z_{s,\alpha}(\mathbf{q}, n, T) := \epsilon^{-1}(\mathbf{q}, n, T) \bar{Z}_{s,\alpha}(\mathbf{q}, n, T), \qquad (3)$$

where $e$ is the electric charge, $A$ the area of the primitive cell of the crystal, $\epsilon^{-1}$ is the two-dimensional macroscopic dielectric function which we will introduce in the next section and $t$ is the typical scale of the electronic response along the out-of-plane direction. The above expressions can be viewed as a generalization of static effective charges tensors to the case of materials with non integer electronic statistical occupations. The long-range components of the dynamical matrix and of the EPI may then be expressed as a function of the effective charges as

$$C_{ss',\alpha\beta}^{\text{L}}(\mathbf{q}, n, T) = \frac{e^2 q^2}{A} v(\mathbf{q}) \bar{Z}_{s,\alpha}^{\text{c.c.}}(\mathbf{q}, n, T) Z_{s',\beta}(\mathbf{q}, n, T) \qquad (4)$$

where $v(\mathbf{q})$ is the two-dimensional Fourier transform of the Coulomb potential, and as

$$g_{\mathbf{q}\nu,mm'}^{\text{L}}(\mathbf{k}) = \frac{ie^2q}{A} v(\mathbf{q}) \langle u_{m\mathbf{k}+\mathbf{q}} | u_{m'\mathbf{k}} \rangle \times \sum_{s\alpha} Z_{s,\alpha}(\mathbf{q}, n, T) e_{s,\alpha}^{\nu}(\mathbf{q}) l_{\mathbf{q}\nu} \left( \frac{M_0}{M_s} \right)^{1/2}, \qquad (5)$$

where $e_{s,\alpha}^{\nu}$ is the eig displacement of the vibrational mode $\nu$ and $l_{\mathbf{q}\nu}$ its zero point motion amplitude, while $u(\mathbf{r}, z)$ are the $\mathbf{r}$-periodic part of the Bloch function and $M_0$ and $M_s$ are, respectively, a reference atomic mass and the atomic mass of the atom s.

Beside the technical interest regarding an improved accuracy in the description of the phononic properties, at a deeper level our framework allows a precise understanding of screening mechanisms at a static level. The consequent extrapolation of the bare unscreened couplings at any given doping and temperature represents a necessary prerequisite to the inclusion of dynamical effects in the description of the screened interactions, which may be of crucial relevance, e.g., in the vicinity of the plasmon resonances.
## B. 2d electrostatic and response functions

Quasi-2d systems are periodic and infinite systems in two dimensions with a finite extension in the third spatial direction—in the case of monolayers this is on the order of the atomic scale, whereas for larger thin films it can reach up to the order of micron. Such systems present a natural distinction between in-plane and out-of-plane properties. Indeed, a periodical quasi-2d system is naturally described by reciprocal space variables in the plane of the material, and a real space variable in the out-of-plane direction. Without loss of generality, the out-of-plane direction is aligned to the $\hat{z}$ direction of a Cartesian coordinate system. 2d quasi-momenta are noted $\mathbf{k}$ for electrons and $\mathbf{q}$ for phonons. To simplify the notation, we will not distinguish between 2d and 3d vectors, whose nature can be inferred from the context. The quasi-2d nature of the problem is reflected in the form of the Bloch theorem for periodic systems

$$ \psi_{m\mathbf{k}}(\mathbf{r}, z) = \frac{e^{i\mathbf{k}\cdot\mathbf{r}}}{\sqrt{N}} u_{m\mathbf{k}}(\mathbf{r}, z), \quad (6) $$

where $m$ is the band index, $u_{m\mathbf{k}}(\mathbf{r}, z)$ is the $\mathbf{r}$-periodic part of the Bloch function and $N$ is the number of cells in the Born-von Karman supercell. The electrostatics of quasi-2d systems is then formulated as a function of $(\mathbf{q}, z)$. In these variables, the Coulomb kernel reads as [40, 41] (see App. A1 for the transform conventions and notations)

$$ v(\mathbf{q}, z - z') = 4\pi \int_{-\infty}^{\infty} \frac{dq_z}{2\pi} \frac{e^{iq_z(z-z')}}{q^2 + q_z^2} = 2\pi \frac{e^{-q|z-z'|}}{q}. \quad (7) $$

The above kernel is involved in the expression of the dielectric response, which in turn determines the LRCs of the dynamical matrix and EPI, as shown in detail in Secs. II D 1 and II D 2. Therefore, all the response functions related to the dielectric one (whose definitions are given in App. A3) need to be expressed as a function of $(\mathbf{q}, z)$. In this spirit, the independent particle polarizability (IPP) of the Khon-Sham system—i.e. the density-density response function of an independent particle system— is written as [63]

$$ \chi^0(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}', z, z') = \frac{2e^2}{A} \sum_{mm'k\mathbf{q}} \frac{f_{m\mathbf{k}} - f_{m'k+\mathbf{q}}}{\epsilon_{m\mathbf{k}} - \epsilon_{m'k+\mathbf{q}}} \int d\mathbf{r} \\ u_{m'k+\mathbf{q}+\mathbf{G}}(\mathbf{r}, z) u_{m\mathbf{k}}^{\text{c.c.}}(\mathbf{r}, z) \int d\mathbf{r}' u_{m\mathbf{k}}(\mathbf{r}', z') u_{m'k+\mathbf{q}+\mathbf{G}'}^{\text{c.c.}}(\mathbf{r}', z'), \quad (8) $$

where $f_{m\mathbf{k}}$ is the Fermi-Dirac occupation of the states with energy $\epsilon_{m\mathbf{k}}$, the factor 2 takes in account spin degeneracy, the spatial integration runs over the unit cell and we have used that

$$ u_{m'k+q+G}(\mathbf{r}', z') = e^{-iG \cdot \mathbf{r}'} u_{m'k+q}(\mathbf{r}', z'), \quad (9) $$

where the $\mathbf{G}$ are reciprocal lattice vectors. The dependence of Eq. 8 on the carrier concentration $n$ and temperature $T$, which is present both in the Fermi-Dirac distributions and in the periodic part of the Bloch wavefunctions (we disregard the latter in this work), has been left implicit not to overburden notation. We will follow the same rule in the rest of this work when possible.

Knowing the expressions for $\chi^0$ and for the Coulomb kernel, we can express the electronic dielectric response matrix, in the RPA for the Khon-Sham ground-state [63], as

$$ \epsilon(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}', z, z') = \delta(z - z') \delta_{\mathbf{G} \mathbf{G}'} \quad (10) \\ -2\pi \int dz'' \frac{e^{-|\mathbf{q}+\mathbf{G}||z-z''|}}{|\mathbf{q}+\mathbf{G}|} \chi^0(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}', z'', z'). $$

The form of Eqs. 8 and 10 allows for effective approximations of the $z, z'$ dependence of the response functions. We first assume that the periodic part of the Bloch's wave functions can be approximated as

$$ u_{m\mathbf{k}}(\mathbf{r}, z) = u_{m\mathbf{k}}(\mathbf{r}) \frac{1}{\sqrt{t}} \theta\left(\frac{t}{2} - |z|\right), \quad (11) $$

where $\theta$ is the Heaviside function and $t$ is defined as the finite layer thickness outside which the electronic cloud vanishes completely. This approximation corresponds to consider the quasi-2d material as an electronically compact homogeneous layer along the out-of-plane direction. One could choose more accurate forms for the $z$ dependence of the wavefunction, but the asymptotic long range expansions of in-plane quantities, that are the focus of this work, do not depend on such choice, as elucidated in App. B.

With the approximation of Eq. 11, Eq. 8 becomes

$$ \chi^0(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}', z, z') = \frac{1}{t^2} \theta\left(\frac{t}{2} - |z|\right) \theta\left(\frac{t}{2} - |z'|\right) \times \\ \chi^0(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}'), \quad (12) $$

where $\chi^0(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}')$ corresponds exactly to the IPP of a two dimensional system [64, 65] (see App. A5). In other words, the approximation of Eq. 11 implies that the IPP of a quasi-2d material can be modelled for the 2d case and then extended uniformly along the $z$ direction inside the layer thickness, while the presence of the $1/t^2$ pre-factor assures the correct dimensionality of the response.

Next, the out-of-plane variables are integrated out of the response functions, via an average along the out-of-plane direction. We define the layer averaged dielectric matrix as (see App. A3)

$$ \tilde{\epsilon}(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}') = \int_{-\frac{t}{2}}^{\frac{t}{2}} dz \frac{dz'}{t} \epsilon(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}', z, z'). \quad (13) $$
This quantity relates to the 2D IPP as follows (see App. A 5):

$$ \tilde{\epsilon}(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}') = \delta_{\mathbf{G}, \mathbf{G}'} - \tilde{v}(\mathbf{q} + \mathbf{G})\chi^0(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}'), \quad (14) $$

where

$$ \tilde{v}(\mathbf{q} + \mathbf{G}) = \frac{4\pi}{|\mathbf{q} + \mathbf{G}|^2 t} \times \left( 1 - \frac{2}{|\mathbf{q} + \mathbf{G}|t} e^{-|\mathbf{q} + \mathbf{G}|\frac{t}{2}} \sinh\left(|\mathbf{q} + \mathbf{G}|\frac{t}{2}\right) \right). \quad (15) $$

Eq. 14 and 15 are particularly pleasant because we can deduce the asymptotic behaviour of the dielectric response function in relevant limits where the Coulomb kernel assumes the simple expression

$$ \tilde{v}(\mathbf{q} + \mathbf{G}) = \begin{cases} \frac{2\pi}{|\mathbf{q} + \mathbf{G}|} & |\mathbf{q} + \mathbf{G}|t \ll 1 \\ \frac{4\pi}{|\mathbf{q} + \mathbf{G}|^2 t} & |\mathbf{q} + \mathbf{G}|t \gg 1 \end{cases} . \quad (16) $$

The above limits for a quasi-2d material are the aforementioned thin limit—$|\mathbf{q} + \mathbf{G}|t \ll 1$—and the thick limit—$|\mathbf{q} + \mathbf{G}|t \gg 1$, where $\mathbf{q}$ is still intended to be small in order to allow Taylor expansions. In these limits Eq. 16 shows that the dependence of the Coulomb kernel upon in-plane components of the wavevector assumes the formal expression typical of, respectively, a 2d and a 3d system.

This observation is in line with what was already noted in Ref. [41], i.e. there exists a scale that discriminates between the two and three dimensional character of the response functions. The same can be easily shown also for the interacting polarizability $\chi$ by considering its Dyson equation. Of course, for a realistic material there will be a crossover between the two regimes as a function of $|\mathbf{q}|$ or $|\mathbf{G}|$. For example, if we consider systems where $t$ is of the order of the lattice parameter, as done in this work, then in the long wavelength limit $\epsilon(\mathbf{q}, \mathbf{q}) := \epsilon(\mathbf{q})$— called the ‘head’ of the dielectric matrix—will behave as for the case of a 2d material while $\epsilon(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}')$— the ‘body’—will mostly follow the 3d behaviour (see also App. A 6 for terminology). If instead the lattice parameter is very large, we can expect the full matrix to show 2d behaviours. Conversely, if we consider a multilayer structure with a small in-plane unit cell, then all the elements of the response are expected to become substantially 3d with the increase of the number of layers. When all the relevant elements of the response are in the thin and/or thick limit, then our layer averaging procedure is well justified, as explained in App. B. Intermediate crossover regimes, existing in ranges that depend on material-dependant internal parameters, are far more complicated to treat since general asymptotic formulae cannot be deduced. Brute-force first-principles methods then have to be used to evaluate the response functions.

From now on, unless otherwise stated, we drop the tilde notation for the layer-averaged quantities, which can be easily recognized from the context.

### C. Effective charges

The electrostatic response of materials to external perturbing potentials can be effectively described in terms of the effective charges, which are macroscopic quantities in the sense that they do not depend on $\mathbf{G}$ vectors explicitly. They do contain, however, all the information regarding the microscopic response of the material, which instead explicitly depends on the $\mathbf{G}$ vectors components (the so-called local-fields). To obtain an expression for the effective charges defined in Sec. II A, we first introduce the reciprocal space expression of the screened Coulomb potential [66]

$$ w(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}') = \epsilon^{-1}(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}')v(\mathbf{q} + \mathbf{G}'). \quad (17) $$

Explicit formulae relating the matrix elements of $w$ and its inverse $w^{-1}$ are given in App. A 6; the long wavelength expansion of $w$ is given in A 6 a, and the limits for its macroscopic components are

$$ w(\mathbf{q}) = \epsilon^{-1}(\mathbf{q})v(\mathbf{q}) = \begin{cases} 2\pi \frac{1}{\mathbf{q}+\mathbf{q}'\cdot\mathbf{B}'\cdot\mathbf{q}} & \text{thin} \\ \frac{4\pi}{t} \frac{1}{\mathbf{q}\cdot\mathbf{B}''\cdot\mathbf{q}} & \text{thick} \end{cases}, \quad (18) $$

were $B'$ and $B''$ are defined from the asymptotic form of Eq. A45. They are, respectively, the tensorial generalization of the effective dielectric screening length $r_{eff}$ as defined in Ref. [41] for a material with $\epsilon_{eff}^0 \approx 1$, and of the electronic dielectric constant $\epsilon^\infty$.

Then, we consider a collective displacement of the atoms of type $s$ along the Cartesian direction $\alpha$ modulated via a wavevector $\mathbf{q}$. From the electronic point of view, this may be regarded as an external charge density perturbation, which can be expressed for in-plane displacements as (see App. A 2 and Ref. [67])

$$ t\delta\rho_{s,\alpha}^{\text{ext}}(\mathbf{q} + \mathbf{G}) = -i\frac{Z_s e}{A} [(q_\alpha + G_\alpha)e^{-i\mathbf{G}\cdot\boldsymbol{\tau}_s}], \quad (19) $$

where $Z_s e$ is the atomic charge and $\boldsymbol{\tau}_s$ indicates the position of the atom in the unit cell. Its relation to the total electrostatic potential $\delta V_{s,\beta}^{\text{tot}}(\mathbf{q})$, obtained as the sum of the induced and the external potentials, in the thin and thick limits can simply be expressed in terms of the $w$ tensor (see App. C) as

$$ \delta V_{s,\alpha}^{\text{tot}}(\mathbf{q}) = \sum_{\mathbf{G}'} w(\mathbf{q}, \mathbf{q} + \mathbf{G}') t \delta\rho_{s,\alpha}^{\text{ext}}(\mathbf{q} + \mathbf{G}'). \quad (20) $$

As shown in App. A 7 in the RPA, we can define the macroscopically unscreened density response $\delta\bar{\rho}^{\text{tot}}$:

$$ \delta\rho_{s,\alpha}^{\text{tot}}(\mathbf{q}) = \epsilon^{-1}(\mathbf{q})\delta\bar{\rho}_{s,\alpha}^{\text{tot}}(\mathbf{q}), \quad \delta\bar{\rho}_{s,\alpha}^{\text{tot}}(\mathbf{q}) = \frac{\delta V_{s,\alpha}^{\text{tot}}(\mathbf{q})}{w(\mathbf{q})} \quad (21) $$

and conveniently write

$$ \frac{\delta V_{s,\alpha}^{\text{tot}}(\mathbf{q})}{w(\mathbf{q})} = -i \frac{e q}{A} \bar{Z}_{s,\alpha}(\mathbf{q}). \quad (22) $$
We obtain at last

$$ \bar{Z}_{s,\alpha}(\mathbf{q}) = i \frac{Aq}{e} \sum_{\mathbf{G}'} \frac{\epsilon^{-1}(\mathbf{q}, \mathbf{q} + \mathbf{G}')}{\epsilon^{-1}(\mathbf{q})} \delta V_{s,\alpha}^{\text{ext}}(\mathbf{q} + \mathbf{G}'). \quad (23) $$

From the above equation it is evident that both $\bar{Z}_{s,\alpha}(\mathbf{q})$ and $Z_{s,\alpha}(\mathbf{q})$ include the local-fields components of the response that contribute to the variation of the macroscopic potential. They differ in the inclusion of the long-range, macroscopic components of the screening response through Eq. 3. As anticipated, the macroscopic character of Eq. 22 is evident from the absence of any explicit reference to the lattice vectors $\mathbf{G}$ (see also App. C). As shown in App. A 7, $\delta\bar{\rho}_{s,\alpha}^{\text{tot}}(\mathbf{q})$ is analytical and allows for a Taylor expansion. Therefore, $\bar{Z}_{s,\alpha}(\mathbf{q})$ can be written as (expliciting the dependence on the doping level $n$ and $T$ dependence)

$$ \bar{Z}_{s,\alpha}(\mathbf{q}, n, T) = \frac{i}{q} M_{s,\alpha}(n, T) + \frac{q\beta}{q} Z_{s,\alpha\beta}^{*}(n, T) + \frac{i}{2} \frac{q\beta}{q} q_{\gamma} Q_{s,\alpha\beta\gamma}(n, T) + \dots. \quad (24) $$

Since the charge density change (and the effective charge functions) is a real quantity in direct space, Eq. 24 comprises alternating imaginary and real terms in the reciprocal space expansion. The coefficients of the expansion are site-dependent tensorial quantities (with rank proportional to the order in $q$ of the expansion), and as such they comply with the site symmetries, transforming as the totally symmetric irreducible representation. The rank-1 tensor $M_{s,\alpha}(n, T)$ - a polar vector, transforming as a force - is strictly zero in insulators/semiconductors where $n = 0$, as it arises from intraband terms. In the presence of free carriers, it is allowed only for those atoms whose Wyckoff positions are not fixed by symmetry, i.e., that can be subject to forces that do not lower the crystallographic symmetries, and it is indeed related to the presence of a Fermi energy shift (see Eq. (79) of Ref. [23] and the discussion in Ref. [36]). The second term corresponds to Born effective charge tensors, the third to dynamical effective quadrupole tensors and so on (sums over repeated indexes are intended). In Eq. 24 we highlighted the doping and temperature dependence that mostly comes from intraband contribution in the IPP for $n \neq 0$, as shown in App. A 5.

## D. Long range components

### 1. Dynamical matrix

We are now ready to provide the asymptotic form for the LRC of the dynamical matrix of a generic 2d material, valid for insulators, semiconductors (doped or not) and metals, and which becomes exact in the thin and thick limits. We start from the component of the force constants matrix that gives rise to LRCs, expressed as

a function of the inverse dielectric screening in an all-electron formalism [68]

$$ C_{ss',\alpha\beta}(\mathbf{R}_p, \mathbf{R}_{p'}) = \frac{\partial^2}{\partial(\mathbf{R}_p + \boldsymbol{\tau}_s)_\alpha \partial(\mathbf{R}_{p'} + \boldsymbol{\tau}_{s'})_\beta} \int d\mathbf{r} \int_{-\infty} d\epsilon \epsilon^{-1} (\mathbf{R}_p + \boldsymbol{\tau}_s, \mathbf{r}, \tau_{sz}, z) \frac{\mathcal{Z}_s \mathcal{Z}_{s'} e^2}{|\mathbf{r} - \mathbf{R}_{p'} - \boldsymbol{\tau}_{s'}|}, \quad (25) $$

where $s, s'$ are atomic indexes, $\alpha, \beta$ are Cartesian indexes and the spatial integration runs over the whole crystal, while $\mathbf{R}_p$ and $\boldsymbol{\tau}_s$ indicate the position of the atom $s$ in the cell $p$, as detailed in App. A 1. The above expression is amenable for in-plane derivatives (i.e. for $\alpha, \beta = x, y$) while the out-of-plane direction is more cumbersome and model-dependent. We restrict to the study of in-plane derivatives, i.e. to in-plane modes. The LRCs of the out-of-plane modes are generally less important on the final dispersion [40]; if needed, a more refined approach based on Ref. [40] should be developed. Further restricting to single layer materials with mirror symmetry [69], we can then fix the out-of-plane atomic coordinate $\tau_{sz} = 0$ and can rewrite the LRC, as shown in App. A 8, as

$$ C_{ss',\alpha\beta}^{\text{L}}(\mathbf{q}) = \frac{\mathcal{Z}_s \mathcal{Z}_{s'} e^2}{A} \sum_{\mathbf{G}, \mathbf{G}'} (q_\alpha + G_\alpha) (q_\beta + G'_\beta) \times w(\mathbf{q}, \mathbf{q} + \mathbf{G}') \frac{w^{\text{c.c.}}(\mathbf{q}, \mathbf{q} + \mathbf{G})}{w(\mathbf{q})} e^{i\mathbf{G} \cdot \boldsymbol{\tau}_s - i\mathbf{G}' \cdot \boldsymbol{\tau}_{s'}}. \quad (26) $$

If we now recast Eq. 26 as a function of macroscopic physical quantities alone, i.e. that do not depend on $\mathbf{G}$, as firstly done by Born and Huang starting from a phenomenological theory valid only for the case of undoped semiconductors [70], we end up with Eq. 4. One may also find useful to rewrite Eq. 4 in a more symmetric form, i.e. as a function of the screened Coulomb potential as

$$ C_{ss',\alpha\beta}^{\text{L}}(\mathbf{q}, n, T) = \frac{e^2 q^2}{A} w(\mathbf{q}) \bar{Z}_{s,\alpha}^{\text{c.c.}}(\mathbf{q}, n, T) \bar{Z}_{s',\beta}(\mathbf{q}, n, T). \quad (27) $$

Eq. 4 is non-analytical for semiconductors and insulators due to the asymptotic form of the dielectric screening [36]. The non-analyticity is cured by the presence of free-carriers for metals and doped semiconductors, for which $\lim_{\mathbf{q} \to 0} \epsilon^{-1}(\mathbf{q}) = 0$. For the latter, however, this happens in a vanishingly small region around $\Gamma$ in the limit of zero doping.

To connect to well known formulae, we rewrite Eq. 4 at the leading order in the $\mathbf{q}$ expansion for the case of an insulator or an undoped semiconductor

$$ C_{ss',\alpha\beta}^{\text{L}}(\mathbf{q}) = \frac{e^2}{A} w(\mathbf{q}) \sum_{\gamma} q_{\gamma} Z_{s,\gamma\alpha}^{*} \sum_{\gamma'} q_{\gamma'} Z_{s',\gamma'\beta}^{*}; \quad (28) $$

one can see that the thin limit of Eq. 28 is equivalent to the in-plane component of Eqs. 45-46 of Ref. [40] or to Eq. 4 of Ref. [43], while the thick limit is the standard textbook version of the LRC of the dynamical
<figure><img src="image_2.png" /><figcaption>FIG. 1: Upper) ab-initio band structure and its Wannier interpolation in the vicinity of the Fermi level, taken as energy reference at 0 eV; inset: schematics of disproportionated graphene, with the cartesian reference systems used in this work. Lower) Phonon band structure as obtained from Fourier interpolation.</figcaption></figure>

matrix given by Born and Huang [30] (Eq. 18 of [71]). The crossover between the thin and the thick limit of Eq. 28 can be observed not just as a function of $\mathbf{q}$, where the crossover scale is simply given by the effective dielectric screening length $r_{eff}$, but also as a function of $t$ for those layered materials that can be brought continuously from the one layer setup to the bulk form (as for example increasing the number of layers of h-BN in the AA stacking [42, 43, 72]). As argued in App. B, in this case the crossover scale given by $t$ is to be more appropriately intended as the dielectric thickness of the material.

## 2. EPI

After the derivation of the LRC of the dynamical matrix in Sec. IID 1, we can proceed along the same line and derive the LRC of the EPI. In this case, we start from the all-electron expression [5]

$$g_{\mathbf{q}\nu}(\mathbf{r}, z) = \int d\mathbf{r}'dz' \epsilon^{-1}(\mathbf{r}, \mathbf{r}', z, z') g_{\mathbf{q}\nu}^{\mathbf{b}}(\mathbf{r}', z'), \quad (29)$$

where $g^{\mathbf{b}}$ is the bare electron-phonon coupling induced by a unit displacement of wavevector $\mathbf{q}$ of the atoms along a mode $\nu$

$$g_{\mathbf{q}\nu}^{\mathbf{b}}(\mathbf{r}, z) = e^{-i\mathbf{q} \cdot \mathbf{r}} \Delta_{\mathbf{q}\nu} V^{\text{e-n}}(\mathbf{r}, \{\mathbf{R}_p + \boldsymbol{\tau}_s\}), \quad (30)$$

where

$$V^{\text{e-n}}(\mathbf{r}, \{\mathbf{R}_p + \boldsymbol{\tau}_s\}) = - \sum_{ps} \frac{\mathcal{Z}_s e^2}{|\mathbf{r} - \mathbf{R}_p - \boldsymbol{\tau}_s|}. \quad (31)$$

$V^{\text{e-n}}$ is therefore the cell-periodic electron-nuclei interaction and $\Delta_{\mathbf{q}\nu}$ indicates the variation following a phonon displacement, as more precisely defined in App. A 4. We bracket Eq. 29 with the Bloch functions of Eq. 11, leading to

$$\langle \psi_{m\mathbf{k}+\mathbf{q}} | g_{\mathbf{q}\nu}(\mathbf{r}, z) | \psi_{m'\mathbf{k}} \rangle = \int_{-\frac{t}{2}}^{\frac{t}{2}} \frac{dz'}{t} dz' \sum_{\mathbf{G}\mathbf{G}'\mathbf{G}''} u_{m\mathbf{k}+\mathbf{q}}^{\text{c.c.}}(\mathbf{G}) \times \\ u_{m'\mathbf{k}}(\mathbf{G}') \epsilon^{-1}(\mathbf{q} + \mathbf{G}' - \mathbf{G}, \mathbf{q} + \mathbf{G}'', z, z') g_{\mathbf{q}\nu}^{\mathbf{b}}(\mathbf{G}'', z'), \quad (32)$$

where

$$g_{\mathbf{q}\nu}^{\mathbf{b}}(\mathbf{G}'', z') = i \frac{2\pi e^2}{A} \sum_{s\alpha} \frac{\mathcal{Z}_s e^{-|\mathbf{q}+\mathbf{G}''||z'|}}{|\mathbf{q}+\mathbf{G}''|} \times \\ (q_\alpha + G''_\alpha) e^{-i\mathbf{G}'' \cdot \boldsymbol{\tau}_s} e_{s,\alpha}^{\nu}(\mathbf{q}) l_{\mathbf{q}\nu} \left( \frac{M_0}{M_s} \right)^{1/2}, \quad (33)$$

and $l_{\mathbf{q}\nu}$ is the zero point motion amplitude (App. A 4). Following exactly the same procedures as in Sec. IID 1 we obtain the LRC as

$$g_{\mathbf{q}\nu,mm'}^{\mathbf{L}}(\mathbf{k}) = \frac{i2\pi e^2 |\mathbf{q}|}{A} \sum_{\mathbf{G}\mathbf{G}'} u_{m\mathbf{k}+\mathbf{q}}^{\text{c.c.}}(\mathbf{G}) u_{m'\mathbf{k}}(\mathbf{G}') \times \\ w^{\text{c.c.}}(\mathbf{q}, \mathbf{q} + \mathbf{G}' - \mathbf{G}) \sum_{s,\alpha} \bar{Z}_{s,\alpha}(\mathbf{q}) e_{s,\alpha}^{\nu}(\mathbf{q}) l_{\mathbf{q}\nu} \left( \frac{M_0}{M_s} \right)^{1/2}. \quad (34)$$

Differently from the dynamical matrix case, we haven't yet isolated the head of the $w$ tensor in the above expression. In general, we cannot set $\mathbf{G}' = \mathbf{G}$ blindly because we would lose corrections to the non-leading order expansion of the EPI due to the dependence of the wings of $w$ on $\mathbf{q}$. We will however set $\mathbf{G}' = \mathbf{G}$ and verify a posteriori that this approximation is correct. We do expect the terms coming from $\mathbf{G} \neq \mathbf{G}'$ to be small because $w(\mathbf{q}, \mathbf{q} + \mathbf{G}' - \mathbf{G})$ contributes to a further power of $q$ with respect to the terms coming from $w(\mathbf{q})$ (see also the discussion in Ref. [35]). Eq. 34 then becomes Eq. 5, which may be rewritten as a function of $w(\mathbf{q})$ as

$$g_{\mathbf{q}\nu,mm'}^{\mathbf{L}}(\mathbf{k}) = \frac{ie^2q}{A}w(\mathbf{q}) \langle u_{m\mathbf{k}+\mathbf{q}}|u_{m'\mathbf{k}} \rangle \times \\ \sum_{s\alpha} \bar{Z}_{s,\alpha}(\mathbf{q}, n, T) e_{s,\alpha}^{\nu}(\mathbf{q}) l_{\mathbf{q}\nu} \left( \frac{M_0}{M_s} \right)^{1/2}, \quad (35)$$

where we have supposed that the unperturbed Bloch functions and the phonon polarizations do not change appreciably with $n$ and $T$. To connect with well known formulae, at the leading order for undoped semiconductors we find

$$g_{\mathbf{q}\nu,mm'}^{\mathbf{L}}(\mathbf{k}) = i w(\mathbf{q}) \langle u_{m\mathbf{k}+\mathbf{q}}|u_{m'\mathbf{k}} \rangle \times \\ \sum_{s,\alpha\beta} q_{\beta} Z_{s,\alpha\beta}^{*} e_{s,\alpha}^{\nu}(\mathbf{q}) l_{\mathbf{q}\nu} \left( \frac{M_0}{M_s} \right)^{1/2}. \quad (36)$$
One can see that Eq. 36 in the thin limit is equivalent to Eq. 8 of Ref. [44], while in the thick limit it is equivalent to Eq. 4 of [33]—considering the phase difference as explained App. A1— and to the long wavelength expansion of Eq. 9 of [34] once the matrix element between the periodic part of the Bloch functions is approximated to $\delta_{mn}$.

### III. COMPUTATIONAL APPROACH

#### A. Connection with theory and general strategy

The theoretical framework developed in the previous sections requires the knowledge of the effective charge functions, whose expression is given in Eqs. 3 and 24, in order to determine the LRCs of the dynamical matrix and the EPI. The central quantity is the layer-averaged total charge density change per unit surface, $t\delta\rho_{s,\alpha}^{\text{tot}}(\mathbf{q}, n, T)$, that in the RPA approximation is connected to the macroscopically unscreened and screened effective charges $\bar{Z}_{s,\alpha}(\mathbf{q}, n, T)$ and $Z_{s,\alpha}(\mathbf{q}, n, T)$. More precisely, in App. A7 we show that in RPA it is possible to deduce the macroscopically unscreened density change $\delta\bar{\rho}_{s,\alpha}^{\text{tot}}(\mathbf{q})$, connected to $\bar{Z}_{s,\alpha}(\mathbf{q}, n, T)$ by Eq. A49, by simply imposing that the macroscopic component of the layer-averaged electrostatic potential is zero. We also remark and stress that in RPA the connection between screened and unscreened quantities is simply attained through the macroscopic inverse dielectric function.

To connect to the theoretical derivations with a deeper insight, we can sum up the above observations saying that the layer-averaging procedure performed in the theoretical section and App. A7 was engineered to separate short and long range components of the dynamical matrix and of the EPI, and to mathematically demonstrate the analiticity of the expansion Eq. 24, which terms can be deduced from the computation of $\delta\bar{\rho}_{s,\alpha}^{\text{tot}}(\mathbf{q})$. However, in ab-initio calculations we can in principle retain the out-of-plane dependence in the solution of the response problem for those quantities that we demonstrated to be analytic functions of in-plane momentum. In fact, the reintroduction of the z-dependence here does not impact the analiticity of the expressions.

From an intuitive point of view, this corresponds to treat our slab of material as a compact layer when it responds to the long-range macroscopic electrostatic, but not disregarding its local-field dependence (even the out-of-plane) when looking at the analytical unscreened charges and potentials. We detail in Sec. IIIC how this can be done in practice. We anticipate that within this treatment, for both in-plane and out-of-plane perturbations, our computational approach can be used to deduce macroscopic charge densities and potentials with the correct out-of-plane dependence, which are macroscopically unscreened both from the in-plane dielectric response and from the presence of periodically repeated images. This is of paramount importance if one wants to extract ‘bare’

couplings to be inserted within formalisms to compute remote EPI and their screening in van der Waals heterostructures [57].

We conclude this section by noting that for realistic DFT calculations we will need to include exchange and correlation terms; in this case the connection between screened and unscreened quantities is not as straightforward as in the RPA case. The procedure of unscreening itself is not unique, as detailed in App. D. In the framework of this work, it suffices to say that a good approximation for the ab-initio evaluation of the long range components is attained by computing the inverse dielectric function in the RPA+xc approximation (i.e. RPA plus the response of the approximated DFT exchange-correlation functional [5]) and computing the unscreened effective charge functions by setting to zero the total macroscopic electrostatic+xc potential—more details are given in Sec. IIIC and App. D.

#### B. Evaluation of the macroscopic inverse dielectric function

A practical way to compute the layer-averaged inverse dielectric matrix for a given doping and temperature is to resort to ab-initio calculations. This approach is based on the methodology developed in Ref. [64], and consists in evaluating the inverse dielectric screening as

$$ \tilde{\epsilon}^{-1}(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}') = \frac{\delta\tilde{V}^{\text{tot}}(\mathbf{q} + \mathbf{G})}{\delta\tilde{V}^{\text{ext}}(\mathbf{q} + \mathbf{G}')}, \quad (37) $$

i.e. as the variation of the layer-averaged Khon-Sham potential $\delta\tilde{V}^{\text{tot}}$ as a consequence of the variation of the external potential $\delta\tilde{V}^{\text{ext}}$. In 2d-materials the above equation is valid provided that the Coulomb potential used in the DFT framework has been modified with the Coulomb cutoff technique [41]. Indeed, Eq. 37 correctly takes in account the local-fields (LFs) corrections to $\epsilon^{-1}(\mathbf{q}, n, T)$ if $\delta V^{\text{tot}}(\mathbf{q}, n, T)$ is the total electronic potential change at the end of the DFPT cycle. If instead we stop the DFPT cycle after the first iteration, i.e. when the response of the system is still non-interacting, and use Eq. 37 with the corresponding potential change, we obtain $\epsilon^{-1}(\mathbf{q}, n, T)$ without LFs corrections. We remind that, typically, LFs corrections are on the order of 10% on $\epsilon^{\infty}$ for 3d materials [73–75] or on the inverse dielectric function for 2d materials [64]; we will show in Sec. IVC that such corrections are small also in the system under study in this work.

At the RPA level, the macroscopic dielectric function can also be obtained (through Eq. 14) from the evaluation of the IPP, which is amenable to simplifying approximations in the doped case at finite temperature. As done in Ref. [36], one can split the IPP into two contributions
as:

$$ \chi^{0, \text{dop.}}(\mathbf{q}, n, T) = \chi^{0, \text{undop.}}(\mathbf{q}) + \delta\chi^0(\mathbf{q}, n, T) \quad (38) $$

$$ \chi^{0, \text{undop.}}(\mathbf{q}) = \frac{2}{A} \sum_{mm'k} \frac{\theta(\epsilon_{m'k}) - \theta(\epsilon_{mk})}{\epsilon_{mk} - \epsilon_{m'k+q}} |\langle u_{mk} | u_{m'k+q} \rangle|^2 $$

$$ \delta\chi^0(\mathbf{q}, n, T) = \frac{2}{A} \sum_{mm'k} \frac{\delta f_{mk} - \delta f_{m'k+q}}{\epsilon_{mk} - \epsilon_{m'k+q}} |\langle u_{mk} | u_{m'k+q} \rangle|^2. $$

where $\delta f = f^{\text{dop.}} - f^{\text{undop.}}$ and $\epsilon_{mk}, \epsilon_{m'k}$ are taken with respect to the chemical potential. We can then write

$$ \epsilon(\mathbf{q}, n, T) \approx \frac{1}{\epsilon^{-1, \text{undop.}}(\mathbf{q})} - e^2 v(\mathbf{q}) \delta\chi^0(\mathbf{q}, n, T), \quad (39) $$

and approximate $\epsilon^{-1}(\mathbf{q}, n, T) \approx 1/\epsilon(\mathbf{q}, n, T)$, which corresponds to neglecting LFs. We have also used that, within a rigid-band approximation, band energies and overlaps can be obtained from DFT calculations in the undoped setup. A practical simplification comes from neglecting the full $\mathbf{q}$ dependence of the dielectric response and using the asymptotic expressions

$$ 1/\epsilon^{-1, \text{undop.}}(\mathbf{q}) \approx \begin{cases} 1 + r_{\text{eff}}q & \text{thin} \\ \epsilon^{\infty} & \text{thick} \end{cases}. \quad (40) $$

Further neglecting the sum over valence/conduction states in $\delta\chi^0$ for electron/hole doping, and taking $\langle u_{mk} | u_{m'k+q} \rangle = \delta_{mm'}$ [64] in $\delta\chi^0$, one finally gets for the macroscopic inverse dielectric function:

$$ \epsilon^{-1}(\mathbf{q}, n, T) \approx \begin{cases} [1 + r_{\text{eff}}q - 2\pi e^2/q \delta\chi^0(\mathbf{q}, n, T)]^{-1} & \text{thin} \\ [\epsilon^{-\infty} - \frac{4\pi}{t}e^2/q^2 \delta\chi^0(\mathbf{q}, n, T)]^{-1} & \text{thick} \end{cases} \quad (41) $$

$$ \delta\chi^0(\mathbf{q}, n, T) \approx \frac{2}{A} \sum_{mk} \frac{\delta f_{mk} - \delta f_{mk+q}}{\epsilon_{mk} - \epsilon_{mk+q}}. $$

The above equation depends on quantities that can be computed directly in the undoped case ($r_{\text{eff}}$ and $\epsilon^{\infty}$) or that can be evaluated via Wannier interpolation ($\epsilon_{n\mathbf{k}}$). For this second case, restricting the sum over $m$ only to valence/conduction bands is particularly convenient and allows to consider only a subset of bands which are situated near the chemical potential. Given the above, the $n$ and $T$ dependence of Eq. 41 stems only from the occupation functions $f_{n\mathbf{k}}$, which makes it easy to implement.

### C. Computation of effective charges

Within the Coulomb cutoff technique of Ref. [41], the 2D system is still replicated periodically along the $\hat{z}$ direction, implying a discrete Fourier transform with reciprocal space vectors $G_z$. Standard DFPT calculations readily provide all the $\{G, G_z\}$ components of screened $t\delta\rho^{\text{tot}}(\mathbf{q} + \mathbf{G}, G_z, n, T)$, from which we select in particular the macroscopic and layer-averaged $\mathbf{G} = 0, G_z = 0$

component. To obtain the macroscopically unscreened $\delta\rho^{\text{tot}}(\mathbf{q}, n, T)$ and the related $\tilde{Z}(\mathbf{q}, n, T)$, we perform the DFPT calculation while setting to zero the $\{\mathbf{G} = 0, G_z\}$ components of i) the change of the local part of the pseudopotential and ii) change of the Hartree and exchange-correlation potentials. This corresponds to the condition $\delta V^{\text{tot}}(\mathbf{q}) = 0$. In particular, setting all the $G_z$ components to zero in the DFPT problem (as opposed to disregarding only the $G_z = 0$ component) is coherent with the procedure of layer-averaging the Maxwell's equation, as presented in App. C, and with the formal derivation of the effective charge expansion (App. A 7). We observe that, even though the theory of Sec. II has been developed from an all-electron perspective, the implementation of pseudopotentials in the calculation does not spoil the conclusion regarding the long wavelength expansions presented in Sec. II. This is because the approximations introduced by pseudopotentials may in general fail to reproduce the all-electron results on the response function in the opposite limit, i.e. for $\mathbf{q} \to \infty$ (when core contributions may become visible).

We stress that we only modify the electrostatic problem for the $\{\mathbf{G} = 0, G_z\}$ component, and not for all the $\mathbf{G}$. This means that for $\mathbf{G} \neq 0$ we are still keeping the out-of-plane dependence of the response problem, which translates in retaining the correct out-of-plane behaviour of (unscreened) quantities that we demonstrated to be analytic functions of in-plane momenta. From a formal perspective, this may be thought as restoring the out-of-plane variable when evaluating expressions such as Eq. A48, an operation that does not spoil the analyticity of the expressions. Secondly, for in-plane perturbations our condition is asymptotically equivalent to the procedure used to deduce in-plane effective charge functions via DFPT at $\mathbf{q} = 0$, as shown in the results sections.

### D. Evaluation of the LRCs

We now discuss the general procedure that shall be used within our framework to practically implement an efficient interpolation of the long range components of the dynamical matrix and of the EPI, for general $n$ and $T$.

For insulators and undoped semiconductors Eqs. 4 and 5 (usually approximated using $\langle u_{mk+q} | u_{m'k} \rangle \sim \delta_{mm'}$) are non-analytical asymptotic formulae in reciprocal space deduced in the $\mathbf{q} \to 0$ limit and can therefore be used only in the neighborhoods of $\Gamma$. Such non-analyticity implies that their expression is not prone to be transformed in real space (in order to be then back interpolated on the full BZ) using a finite set of Wannier functions or plane waves. To overcome this problem, within the state-of-the-art methodology one performs an ansatz for the real space transform of Eq. 4 and 5 that is able to recover the $\mathbf{q} \to 0$ asymptotic behaviours of Eqs. 4 and 5 (up to a certain order) once they are transformed back in reciprocal space on the full BZ [34, 76]. Notice that
the ansatz is not uniquely determined. Alternatively, one can perform the ansatz directly in Fourier space by extending each order expansion of Eqs. 4 and 5 at $q+G$ vectors in such a way that periodicity and continuity at zone boundary is fulfilled [34]. This extended expression for the LRCs on the full BZ (which we will refer to as xLRCs or xLs) is then used to perform the following interpolation scheme for a pair of points $(k_0, q_0)$: one first interpolates the well-behaved differences $C^S = C - C^{xL}$ and $g^S = g - g^{xL}$, where C and g are obtained ab-initio on a coarse grid defined on the full BZ, via Fourier or Wannier interpolation [25]; then, the xLRC are re-added evaluating their closed-form expression at $(k_0, q_0)$.

For the case of doped semiconductors, as already discussed, Eq. 4 and 5 are in principle analytical, even though the extension of the region where the non analyticity is cured depends on the magnitude of the carrier concentration. Practically, for small doping levels it is not guaranteed that Wannier interpolation can be performed without isolating the LRCs, but instead we should repeat the above described procedure for each n and T. The xLRCs would then be formally obtained from their expression in the undoped setup substituting $Z_{s,\alpha}^{\text{undop.}}(q) \rightarrow Z_{s,\alpha}(q, n, T)$ and $\bar{Z}_{s,\alpha}^{\text{undop.}}(q) \rightarrow \bar{Z}_{s,\alpha}(q, n, T)$. In practice, for the doping levels and temperatures studied in this work, we will assume (and verify a posteriori) that the analytical, short-range components of the matrix elements are the same as the undoped setup at zero temperature [77], i.e. that $C = C^{\text{S,undop.}} + C^{\text{xL,dop.}}(n, T)$ and $g = g^{\text{S,undop.}} + g^{\text{xL,dop.}}(n, T)$. In this case we can perform only one ab-initio calculation of C and g on the coarse grid for the undoped setup at zero temperature, interpolate $C^{\text{S,undop.}} = C - C^{\text{xL,undop.}}$ and $g^{\text{S,undop.}} = g - g^{\text{xL,undop.}}$ on fine meshes, and finally add $C^{\text{xL,dop.}}(n, T)$ and $g^{\text{xL,dop.}}(n, T)$ with the correct n and T dependencies on the fine grid. Notice that if we are interested in obtaining the correct expressions for C and g only in a small region in the neighborhood of $\Gamma$, then the last step of the above procedure can be simplified since the xLRCs in reciprocal space can be evaluated using directly the asymptotic expansions Eqs. 4 and 5.

# IV. RESULTS

## A. Disproportionated graphene

The sublattice symmetry of graphene can be broken via the interaction with substrates as, e.g., SiC [58–60]. Here, every second carbon atom has a neighbor in the bottom layer so that a different potential is felt by atoms belonging to different sublattices, thus breaking the equivalence between carbon atoms, with a consequent reported band splitting up to 0.5 eV between the valence and conduction bands at the K point. In this setup, so-called disproportionated graphene is equivalent to monolayer h-BN from a symmetry point of view, and non-zero Born effective charge tensors and piezoelectric

coefficients arise in the system [61]. This system may be efficiently simulated creating an imbalance of valence charge between two otherwise equivalent neighbouring carbon atoms, more precisely by adding/subtracting a fractional value of valence electrons (and a compensating ion charge) to their pseudopotentials. The Born effective charge tensor, defined as

$$Z_{s,\alpha\beta}^* = \frac{A}{e} \frac{\partial P_{\alpha}}{\partial v_{s,\beta}(\mathbf{0})} \Big|_{\mathbf{E}=\mathbf{0}} \quad (42)$$

where $P_\alpha$ is a Cartesian component of the polarization vector per unit surface and $v_{s,\beta}(\mathbf{0})$ is the displacement of the atom s along the cartesian direction $\beta$ (as already defined in Sec. IIC and App. A1), is diagonal with two independent components $Z_{\parallel}^*$ and $Z_{\perp}^*$ describing in-plane and out-of-plane Born effective charges [78]. Since we focus on the in-plane components of the Born effective charges only, which are two order of magnitude bigger than the out-of-plane ones in disproportionated graphene, we will adopt the following notation $Z_{s,xx}^* = Z_{s,yy}^* := Z_s^* := (-1)^s Z^*$. This is consistent with the charge neutrality condition for the undoped system (for the doped system, as we will see, the sum over the atom of the Born effective charges may assume non-zero values). The piezoelectric tensor

$$\bar{e}_{\alpha\beta\gamma} = \frac{\partial P_{\alpha}}{\partial \epsilon_{\beta\gamma}} \Big|_{\mathbf{E}=\mathbf{0}} \quad (43)$$

with $\epsilon_{\beta\gamma}$ being the strain tensor, has instead a single independent coefficient $\bar{e}_{yyy} = -\bar{e}_{yxx} = -\bar{e}_{xyx} = -\bar{e}_{xxy}$, where $\hat{x}$ and $\hat{y}$ are the Cartesian directions of the reference frame of the inset of Fig. 1. In the frozen-ion approximation, the piezoelectric tensor is expressed as a function of the dynamical quadrupole charges as [78]

$$\bar{e}_{\alpha\beta\gamma}^{\text{FI}} = -\frac{1}{2A} \sum_{s} (Q_{s,\beta\alpha\gamma} + Q_{s,\alpha\beta\gamma} - Q_{s,\gamma\alpha\beta}) = \frac{i}{A} \left[ \frac{\partial(q\bar{Z}_{s,\beta})}{\partial q_{\alpha}} + \frac{\partial(q\bar{Z}_{s,\alpha})}{\partial q_{\beta}} - \frac{\partial(q\bar{Z}_{s,\gamma})}{\partial q_{\alpha}} \right]_{\mathbf{q}=0} \quad (44)$$

which accordingly admit the non-zero components $Q_{s,yyy} = -Q_{s,yxx} = -Q_{s,xyx} = -Q_{s,xxy} := Q_s$; notice that the definition of $\bar{e}_{\alpha\beta\gamma}^{\text{FI}}$ in terms of second derivatives of $\bar{Z}$ allows the generalization to the case of finite doping and temperature. We also have $\bar{e}_{yyy}^{\text{FI}} = -\bar{e}_{yxx}^{\text{FI}} = -\bar{e}_{xyx}^{\text{FI}} = -\bar{e}_{xxy}^{\text{FI}} := \bar{e}^{\text{FI}}$.

In disproportionated graphene the values of the Born effective charges and of the quadrupoles are topological properties that, adopting a simple tight-binding model consisting in a graphene Dirac Hamiltonian plus an on-site energy diagonal term $\pm\Delta$, do not depend on the magnitude of the opened gap [61] and may be expressed in terms of the Berry curvature of the band manifold. Nonetheless, it is found that along with the reduction of the on-site energy and the band gap, the region of k-point where the Berry curvature is relevantly different from 0 is
closer and closer to **K** [61]. It follows that for this system, even for small dopings at finite temperature, important modifications to the values of $Z_s^*$ and $Q_s$ may appear, differently from what found for 3d 3C-SiC [36]. For this reason, we will study the different qualitative behaviours of the dynamical matrix and EPI in a weak and strong doping regimes (WDR and SDR respectively), as well as intermediate regimes. Disproportionated graphene is a model candidate to benchmark the whole of our theoretical developments and compare exact results (within the DFT framework) with controlled approximations.

As a conclusion to this section, we show in Fig. 1 a schematics of disproportionated graphene and the Cartesian reference frame that we adopt in this work, alongside with its *ab-initio* electronic band structure (first-principles and Wannier interpolation obtained with the computational parameter explained in Sec. IV B) and with the Fourier interpolated phonon dispersion.

## B. Computational parameters

To perform the first-principles calculations we use a private version of the QUANTUM ESPRESSO (QE) code [27], firstly developed in Ref. [53] to compute EELS cross sections via the extraction of the screened density response to external Cartesian perturbations. We use a lattice parameter of $a = 4.6487$ Bohr, PBE-GGA functionals [79] and norm-conserving pseudopotentials where the valence charge of the two carbon atoms has been altered, respectively by $\pm 0.1e$, in order to mimic the sublattice symmetry breaking of disproportionated graphene [61]. The result is a gap opening at the **K**-point of around $E_G \sim 0.4$ eV, as shown in Fig. 1. We sample the BZ with telescopic grids—following Ref. [80]—that in the densest region are equivalent to Monkhorst-Pack grids of dimension $256^2$ in order to well describe the region around the chemical potential, and use a energy cutoff for the plane wave basis set of 90 Ry.

To perform the *ab-initio* calculations in the doped setup we simulate hole densities from to $n \sim 9.5 \times 10^9 \text{cm}^{-2}$ to $n \sim 9.5 \times 10^{12} \text{cm}^{-2}$ at finite temperature. We define an effective Fermi momentum $k_F$ as the wavevector where the hole occupation halves from its value at the top of the valence band. The doping is simulated in a double-gate setup where the carrier concentration introduced in the graphene layer is compensated by two equally distanced gates that are positioned $\pm 8.3704$ Bohrs away from the graphene plane along the $\hat{z}$ direction. The Coulomb kernel is treated within the 2d Coulomb-cutoff technique as developed in Ref. [41], using an interlayer distance between graphene periodic images of $c = 37.7941$ Bohr.

As already anticipated and shown later, the methodology developed in this work requires to perform Wannier interpolation of electronic and vibrational quantities only for the undoped setup. In this case, the Wannier interpolation of the electronic properties is obtained using

the SCDM method for the determination of the starting projections [81, 82] as implemented in Wannier90 [83] and EPW [29], using 5 Wannier functions and a Gaussian entanglement with $\mu$ equal to the top of the valence band and $\sigma = 10$ eV. Wannierization proves to be of fundamental importance when describing the electronic properties of disproportionated graphene. In fact, simple tight-binding models cannot fully account for the complete orbital composition of the Bloch functions implying poor accuracy in the determination of, e.g., the angular dependency of the EPI. The Wannier interpolation of the electron-phonon matrix elements and of the dynamical matrix is performed using a private version of EPW, adapted for the study of the current work, using coarse meshes of $48 \times 48 \times 1$ for the electrons and $12 \times 12 \times 1$ for the phonons. The convergence of the electronic inverse scattering times (see Sec. IV F) is obtained using fine grids of **q** points of dimensions $700^2$ and a Gaussian smearing for the Dirac delta functions of 15 meV.

## C. Inverse dielectric function

The accurate description of the inverse dielectric function is a focal point of our approach. We first note that LFs have a quantitatively small effect on the asymptotic small **q** value of $\epsilon^{-1}(\mathbf{q})$, as shown in Fig. 2. It is evident that in both doped and undoped cases, LFs have an impact on $\epsilon^{-1}$ of the order of some percent for wavevectors smaller than $\mathbf{q} \sim 0.02 \text{Bohr}^{-1}$. Thus, they may be neglected in this region, which is the most interesting in practice. Incidentally, this shows that the inclusion of *xc* effects in the inverse dielectric response is of little importance to our aims. Importantly, we notice that the approximation of the inverse dielectric function via Eq. 41 works very well.

## D. Effective charges

<table><caption>TABLE I: Reciprocal space expansion of $Z_{s,\alpha}(\mathbf{q})$, for the case of undoped disproportionated graphene along different directions of the BZ; $\epsilon^{-1}(\mathbf{q})$ is approximated as in Eq. 46.</caption><thead><tr><th>Z<sub>s,α</sub></th><th colspan="2">Reciprocal space direction</th></tr></thead><tbody><tr><td></td><td>(q, 0)</td><td>(0, q)</td></tr><tr><td>Z<sub>s,x</sub></td><td>(-1)<sup>s</sup>Z<sup>*</sup> + O(q<sup>2</sup>)<br>1+r<sub>eff</sub>q + O(q<sup>2</sup>)</td><td>O(q<sup>2</sup>)</td></tr><tr><td>Z<sub>s,y</sub></td><td>(1/2Q<sub>s</sub>q + O(q<sup>2</sup>))<br>1+r<sub>eff</sub>q + O(q<sup>2</sup>)</td><td>(-1)<sup>s</sup>Z<sup>*</sup>q - 1/2Q<sub>s</sub>q + O(q<sup>2</sup>)<br>1+r<sub>eff</sub>q + O(q<sup>2</sup>)</td></tr></tbody></table>

Now, we can study the macroscopically screened and unscreened layer-averaged effective charges, $Z_{s,\alpha}(\mathbf{q})$ and $\bar{Z}_{s,\alpha}(\mathbf{q})$, for undoped disproportionated graphene along a specific direction in the reciprocal space. We focus on the
<figure><img src="image_3.png" /><figcaption>FIG. 2: Upper) Dielectric response, for disproportionate graphene in the undoped case, from the model of Eq. 46, compared with the DFT calculation (RPA+xc) obtained using Eq. 37 with or without the inclusion of local-fields (LFs) corrections, as explained in Sec. III, along the reciprocal space Cartesian line form (q, 0). The Cartesian reference frame used in this work is depicted in the inset of Fig. 1. Lower) Same quantities for the case of $n \sim 9.5 \times 10^9 \text{cm}^{-2}$ and a temperature of 4 meV, compared with the RPA approximation of Eq. 41. The effective Fermi wavevector $k_F$, defined in the text, is represented by a black vertical dashed line.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIG. 3: Upper) Real part of the macroscopically unscreened ('unscr') and screened ('scr') effective charge for undoped disproportionated graphene at zero temperature, for the carbon atom with increased/decreased valence charge $C^\uparrow/C^\downarrow$. Lower) Same notation for the imaginary part divided by $q$, summed over the atoms in the lowest panel to obtain $2A\tilde{e}^{\text{FI}}$ ('piezo') as intercept with the vertical axis (see Eq. 44). In the upper panel, the continuous horizontal lines correspond to $Z_s^{*,\text{QE}}$, while in the lower one to the fitted values of the quadrupoles $Q_s$ and of $2A\tilde{e}^{\text{FI}}$. The lines passing through the screened charge data are, respectively, $\frac{Z_s^{*,\text{QE}}}{1+r_{\text{effq}}}$, $\frac{Q_s}{1+r_{\text{effq}}}$ and $\frac{2A\tilde{e}^{\text{FI}}}{1+r_{\text{effq}}}$ where $r_{\text{eff}} \sim 97$ Bohr is computed using Eq. 47.</figcaption></figure>

in-plane components of the effective charges. Within this framework, $M_{s,\alpha}(n, T) = 0$ by symmetry considerations.

The choice of the line is performed in order to isolate the different terms of the expansion of Eq. 24. In fact, the symmetry properties of the effective charges tensors enforce specific behaviours, as shown in Tab. I for two possible directions. We choose the (q, 0) line since the leading order coefficients of $\bar{Z}_{s,x}(\mathbf{q})$ and $\bar{Z}_{s,y}(\mathbf{q})$ are respectively the Born effective charges and the dynamical quadrupole tensors. Since there are no lines where the leading order of the expansion of $\bar{Z}_{s,\alpha}(\mathbf{q})$ is represented by the octupole term, we stop our expansion at the quadrupole level. We then plot in the upper panel of Fig. 3 the real part of $\bar{Z}_{s,x}(\mathbf{q})$ and $\bar{Z}_{s,y}(\mathbf{q})$, alongside with the value of Born effective charges calculated

from ab-initio DFPT calculation at $\mathbf{q} = 0$, which we denote by $Z_s^{*,\text{QE}}$. For small wavevectors, both $\text{Re}\bar{Z}_{s,x}(\mathbf{q})$ and $\text{Re}Z_{s,x}(\mathbf{q})$ tend to $Z_s^{*,\text{QE}}$. For larger wavevectors ($q \sim 0.02 \text{ Bohr}^{-1}$), the unscreened charge departs from its asymptotic value. Indeed, the identification of $\bar{Z}_{s,x}(\mathbf{q})$ with the effective charges starts to degrade because we are exiting from the thin limit—in fact, if we evaluate $t \sim 6$ Bohr, which is the typical interlayer distance in graphite, at $0.02 \text{ Bohr}^{-1}$ we have $qt \sim 0.1$.

As it regards $\text{Re}Z_{s,x}(\mathbf{q})$ instead, we notice that it follows
lows the simple model from Ref. [44]

$$ \mathrm{Re}Z_{s,x}(\mathbf{q}) = \epsilon^{-1}(\mathbf{q})Z_s^{*,\mathrm{QE}} \quad (45) $$

$$ \epsilon^{-1}(\mathbf{q}) = \frac{1}{1 + r_{\mathrm{eff}}q}, \quad (46) $$

$$ r_{\mathrm{eff}} = (\epsilon^{\infty, \mathrm{QE}} - 1) \frac{c}{2}, \quad (47) $$

way beyond the expected regime of validity, where $c \gg t$ is the separation of the periodic images used in the simulation. Dividing the Born effective charges with the simplest thin limit asymptotic expression of $\epsilon^{-1}(\mathbf{q})$ thus works surprisingly well.

This seems to be a feature common to different materials [43, 44, 57], due to cancellation of errors between the wrong estimations of the effective charge by its asymptotic value and of the screening with respect to the local field inclusion at finite wavevectors. All the above conclusions can be drawn also analyzing the quadrupole and piezoelectric tensors, as done in the lower panel of Fig. 3. At difference with the Born effective charge case, the ab-initio calculation of quadrupole from DFPT in QE is not implemented at the time of writing. We therefore estimate the quadrupole values directly from the asymptotic values of $-\mathrm{Im}2\bar{Z}_{s,y}(\mathbf{q})/q$. In the same way, we estimate also the value of the frozen ion piezoelectric tensor, as defined in Eq. 44. We find that $\bar{e}^{\mathrm{FI}} = 6.1 \times 10^{-10} \mathrm{C/m}$, in perfect agreement with the values found in Ref. [61].

We now consider the same quantities but in the doped case, starting with the WDR and SDR at, respectively, $n \sim 9.5 \times 10^9 \mathrm{cm}^{-2}$ and $n \sim 1.9 \times 10^{12} \mathrm{cm}^{-2}$. In the first case, to mimic the zero temperature regime we use a Fermi-Dirac (FD) occupation with a temperature of 4 meV, while in the second we use a Methfessel-Paxton (MP) distribution [84] with a smearing of 40 meV—we will come back on elucidating this choice later. We plot in Fig. 4 the same quantities as Fig. 3, respectively for the WDR on the left column and for the SDR on the right column. In the WDR regime we observe a reduction of the asymptotic value of $\mathrm{Re}\bar{Z}_{s,x}(\mathbf{q}, n, T)$ with respect to the undoped case that only slightly impacts on the change of $\mathrm{Re}Z_{s,x}(\mathbf{q}, n, T)$ at small wavevector, which is instead mostly determined by the influence of metallic screening. We also notice that in the regime where metallicity fades out, i.e. at large $\mathbf{q}$, the charge response returns to the typical value for the semiconductor response.

The same conclusions cannot seemingly be drawn for the quadrupole term, since peculiar behaviours arise at small $\mathbf{q}$. This is not a serious issue since the relevant macroscopic parameter which enters the EPI in the long wavelength limit is the piezoelectric tensor, which is well-behaved. Indeed, in general at the quadrupole order we see that the added carriers substantially alter the unscreened part of the response through the $C(\mathbf{q}, n, T)$ term of Eq. A46, which originates mainly from the intraband contributions to the IPP coming from partially filled bands. In the charge expression, this means that $-\mathrm{Im}2\bar{Z}_{s,x}(\mathbf{q}, n, T)/q$ strongly departs from the behaviour shown in the undoped case. Nonetheless, if we look at

the value of the piezoelectric tensor in the WDR as defined by the second derivative of the $\bar{Z}$ quantities (Eq. 44), we notice that this is very similar to $\bar{e}^{\mathrm{FI}}$, because the sum over atoms of the term coming from $C(\mathbf{q}, n, T)$ results to cancel out almost completely. It follows that $-\mathrm{Im}2\sum_s Z_{s,x}(\mathbf{q}, n, T)/q \approx 2A\epsilon^{-1}(\mathbf{q}, n, T)\bar{e}^{\mathrm{FI}}$. The independence of the frozen ion piezoelectric tensor from doping in the WDR is interesting since the acoustic EPI is directly proportional to $\bar{e}^{\mathrm{FI}}$ in the long wavelength limit. This justifies neglecting the complicated doping and temperature dependence coming from the $C(\mathbf{q}, n, T)$ term in the evaluation of the leading orders of the EPI.

In the SDR the situation is very different. The value of the Born effective charges is completely suppressed and so is the value of the dynamical quadruple tensors, due to the strong modifications of the wings of the dielectric response at high doping concentrations. Notice that now the sum rule for the Born effective charges $\sum_s Z_s^* = 0$ is strongly broken. This is allowed because the total dynamical matrix is expressed as a function of Eq. A54 as $C_{ss',\alpha\beta}^{\mathrm{tot}}(\mathbf{q}) = C_{ss',\alpha\beta}(\mathbf{q}) - \delta_{ss'} \sum_{s''} C_{ss'',\alpha\beta}(\mathbf{0})$, and the relevant translational invariance condition $\lim_{\mathbf{q} \to 0} \sum_{ss'} C_{ss',\alpha\beta}^{\mathrm{tot}}(\mathbf{q}) = C_{ss',\alpha\beta}^{\mathrm{tot}}(\mathbf{0})$ [68] is trivially satisfied in the metallic regime of the dielectric response even if $\sum_s Z_s^* \neq 0$.

To conclude, the qualitative difference between the WDR and the SDR is that in the first regime the $n, T$ dependencies of the Born effective charges and of the piezoelectric tensor can be considered, to a good approximation, to enter only in the head of the inverse dielectric screening. This is highly important because it implies a simple operative procedure to determine the LRCs of the dynamical matrix and EPIs, disregarding the effects of doping on unscreened effective charges (the C term in Eq. A46). To quantify this observation, Fig. 5 shows the ab-initio screened charges of Fig. 4 together with the expressions $\epsilon^{-1}(\mathbf{q}, n, T)Z_s^{*\mathrm{QE}}$, $\epsilon^{-1}(\mathbf{q}, n, T)Q_s$ and $2A\epsilon^{-1}(\mathbf{q}, n, T)\bar{e}^{\mathrm{FI}}$, where $\epsilon^{-1}(\mathbf{q}, n, T)$ is evaluated via Wannier interpolation of Eq. 41. We also perform the same comparison in Fig. 6 for $\mathrm{Re}\bar{Z}(\mathbf{q}, n, T)$ for several different densities, using a MP smearing of 40 meV to evaluate the statistical occupations. MP smearing has the drawback of losing the correspondence with a physical temperature, but it requires coarser samplings of the BZ in order to get converged results, which is of extreme importance in order to perform numerous calculations. In general, we find a good agreement between ab-initio calculation and the RPA modelling but, as expected, this increasingly deteriorates while approaching the SDR.

## E. Frequencies and EPI

Now that the behaviour of $\bar{Z}_{s,\alpha}(\mathbf{q}, n, T)$ and $Z_{s,\alpha}(\mathbf{q}, n, T)$ is understood, we turn to the interpolation of the dynamical matrix and EPI, using the procedure exposed in Sec. III D. In principle, Eqs. 4 and 5 are rigorously valid only in the thin and thick limits. Nonetheless,
<figure><img src="image_5.png" /><figcaption>FIG. 4: Same quantities and same symbol-colour scheme as in Fig. 3, for the weak and strong doping regimes (left and right columns). The horizontal lines are the values of the asymptotic quantities deduced for the undoped setup.</figcaption></figure>

we can expect the validity of Eq. 5 to be extended like that of $Z_{s,\alpha}(\mathbf{q}, n, T)$ in Eq. 46. For Eq. 4 the considerations are a bit different since the LRC of the dynamical matrix depends both on the unscreened and the screened charge density changes. Here, the validity of Eq. 4 is extended using for $\tilde{Z}_{s,\alpha}(\mathbf{q}, n, T)$ the leading order asymptotic value obtained for the undoped setup, i.e. $Z_s^{*,QE}$. This is coherent with the fact that at large $\mathbf{q}$ the LRC perform a crossover to the value of the undoped setup. At sufficiently small wavevectors instead, the correct numerical value of $\tilde{Z}_{s,\alpha}(\mathbf{q}, n, T)$ is not thoroughly needed since the whole LRC is suppressed via screening.

As discussed in the previous section, the evaluation of $Z_{s,\alpha}(\mathbf{q}, n, T)$ may be performed exploiting the RPA approximation of the head of the inverse dielectric screening. If this is not precise enough, as in the SDR, $Z_{s,\alpha}(\mathbf{q}, n, T)$ may be evaluated *ab-initio*. In that case, we shall perform calculations on the minimum number of lines that are required by a symmetry analysis to describe all the relevant independent components of $Z_{s,\alpha}(\mathbf{q}, n, T)$,

and then use symmetry relations to complete the description on the full BZ. Also, we can sample $Z_{s,\alpha}(\mathbf{q}, n, T)$ on each needed line on a relatively coarse set of points and then perform an inexpensive linear interpolation on finer sets of points. All the various steps that compose our strategy carefully avoid to resort to brute-force calculations of all the matrix elements of the dynamical matrix and EPI on the full BZ. Of course, for special regimes which fall very far from the thin or thick limits, our machinery is no more applicable.

The previous observations are quantified in Fig. 7, interpolating on the line $(q, 0)$ for the dynamical matrix and EPI for the WDR and SDR (left and right columns). In both cases the LRCs are evaluated using directly the $Z_{s,\alpha}(\mathbf{q}, n, T)$ coming from *ab-initio* calculations. The long-range features of the coupling and of the phonon dispersion are strongly suppressed in the region near $k_F$, and an excellent agreement between *ab-initio* calculations and interpolated quantities is obtained. Finally, as for effective charges, frequencies and EPI values
<figure><img src="image_6.png" /><figcaption>FIG. 5: Same quantities and same symbol-colour scheme as in Fig. 3, for the weak doping regime, compared to the results obtained screening $Z_s^{*,QE}$, $Q_s$ and $\bar{e}^{FI}$ with the $\epsilon^{-1}(\mathbf{q}, n, T)$ evaluated via Wannier interpolation of Eq. 41 (dashed lines).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIG. 6: Comparison of Re$Z_{s,x}(\mathbf{q})$ obtained with ab-initio calculations and the RPA approximation as explained in the text (upper and lower), for different values of hole concentration per unit cell, reported in units of $10^{-3}$e units in the legend.</figcaption></figure>

for different dopings are plotted in Fig. 8. Ab-initio results are compared with those obtained computing the LRCs with the RPA approximation to $\epsilon^{-1}(\mathbf{q}, n, T)$. As for the charges, we find that the general trends agree in the WDR while they are increasingly different while approaching the SDR.

## F. Lifetimes

We can finally quantify the impact that the correct descriptions of the EPI and phonon frequencies have on physical observables; in particular, we compute the electronic inverse lifetimes evaluated in the Self Energy Relaxation Time Approximation (SERTA) [12] as

$$ \tau_{m\mathbf{k}}^{-1} = \frac{2\pi}{\hbar} \sum_{m'\nu} \int \frac{d\mathbf{q}}{A_{\text{BZ}}} |g_{\mathbf{q}\nu, m'm}(\mathbf{k})|^2 \times \left[ (1 - f_{m'\mathbf{k}+\mathbf{q}} + n_{\nu\mathbf{q}})\delta(\epsilon_{m\mathbf{k}} - \epsilon_{m'\mathbf{k}+\mathbf{q}} - \hbar\omega_{\nu\mathbf{q}}) + (f_{m'\mathbf{k}+\mathbf{q}} + n_{\nu\mathbf{q}})\delta(\epsilon_{m\mathbf{k}} - \epsilon_{m'\mathbf{k}+\mathbf{q}} + \hbar\omega_{\nu\mathbf{q}}) \right], \quad (48) $$

where $A_{\text{BZ}}$ is the area of the BZ, and $n_{\nu\mathbf{q}}$ is the Bose-Einstein occupation factor of the phonon of frequency $\omega_{\nu\mathbf{q}}$.

We compare the cases where the asymptotic expressions of the LRCs are evaluated using the effective charges of the undoped setup screened by the undoped dielectric screening or with correct screening dependence on $n$ and $T$, within the RPA approximation to $\epsilon^{-1}(\mathbf{q}, n, T)$. We choose a setup tailored to highlight the differences between the two approaches. We fix the chemical potential at the top of valence band and set $T = 300$ K, so that the numerical difference in the approaches is only due to screening. We do not take into account the change in the effective charges values due to $n$ and $T$ since this is highly material dependent and the current aim is to compare the present treatment of screening with state-of-the-art methods, more than a refined calculation of the electronic lifetimes. As shown in Fig. 9, discarding the doping dependence of the screening, as currently done in most state-of-the-art first-principles calculations within the rigid-band approximation, implies a
<figure><img src="image_8.png" /><figcaption>FIG. 7: Values for the TO and LO phonon frequencies (upper row) and of the EPI (central and lower rows) for the undoped (circles) and doped (squares) setups, in the WDR (left column) and in the SDR (right column). The central row shows the Frölich coupling with the LO mode (LO), while the lower one is relative to the coupling between the TA mode and the dynamical effective quadrupoles. The continuous and dashed lines represent the Wannier interpolation of, respectively, the undoped and the doped setups. The TO modes and their EPI are not showed for the doped case as they are substantially left invariant by doping.</figcaption></figure>

very strong overestimation of the inverse scattering times in the region near the top of the valence band. This is mostly attributed to a wrong estimation of the piezoelectric coupling between electrons and TA phonons. For electronic states that are below 200 meV from the top of the valence band, we also observe an overestimation of the scattering times which are mostly determined by the Γ optical phonons, due to the overestimation of the Frölich coupling. As mentioned before, in the case of disproportionated graphene even stronger reductions of the electronic inverse lifetimes may come from the change of the effective charge tensors with doping and temperature, which are routinely not taken in account in state-of-the-art calculations. Notice that the value of the scattering times does not change further relevantly with increasing the number of carriers. This means that our lifetimes are

one typical of the high doping regime where dynamical effects are negligible. Also, such dynamical effects may not prevent the screening of piezoelectric coupling since plasma and acoustic frequencies may be of the same order of magnitude.

V. CONCLUSIONS

We analyse the dependence on doping and temperature of effective charges, EPI and phonon frequencies in quasi-two-dimensional doped semiconductors, within a linear-response dielectric-matrix formulation that allows for controlled approximations of the effect of electronic screening. We further propose a fast and accurate interpolation method based on Wannier functions
<figure><img src="image_9.png" /><figcaption>FIG. 8: Same quantities as in Fig. 7, compared between ab-initio calculations (left column) and the RPA approach as described in the text (right column).</figcaption></figure>

that enables a quantitative analysis at a feasible com-
putational cost. We show that neglecting free-carrier
screening on the piezoelectric and Fröhlich interactions,
as done in state-of-the-art computational approaches,
leads to a substantial overestimation of scattering rates
in specific doping-temperature regimes, which may have
a strong impact on the determination of transport prop-
erties. However, our general formulation is not limited
to those couplings, but applies to any EPI that are ac-
cessible within DFPT. The proposed approach for deal-
ing with electronic screening lays the foundation for fur-
ther extensions tackling other less conventional types of
EPI, such as the vector coupling or the electron two-
phonon scattering, that may play an important role in
polar metals and doped quantum paraelectrics [85–88].
As a concluding remark, we mention that the extension
to finite-frequencies dependence within a time-dependent
approach can also give access to nonadiabatic effects on
effective charges and, hence, on lattice dielectric proper-
ties [80, 89] which contributes to the shapes of the dy-

namical structure factor probed by EELS [53] and other
inelastic scattering experiments [90].

Acknowledgements—We acknowledge the European Union’s Horizon 2020 research and innovation program under grant agreements no. 881603-Graphene Core3 and the MORE-TEM ERC-SYN project, grant agreement No 951215. We acknowledge that the results of this research have been achieved using the DECI resource Mahti CSC based in Finland at https://research.csc.fi/-/mahti with support from the PRACE aisbl; we also acknowledge PRACE for awarding us access to Joliot-Curie Rome at TGCC, France. This paper is dedicated to the living memory of Dr. Nicola Bonini.
<figure><img src="image_10.png" /><figcaption>FIG. 9: Comparison of inverse lifetime scattering computed along high symmetry lines using Eq. (48) where the LRCs are computed using Eqs. (4)-(5) with and without n and T dependence: the inverse lifetimes are given in absolute magnitude in the right panel, and as linewidths in the left and central ones. The light blue dots and the central panel are computed with the method proposed in this work, while grey dots and the left panel are obtained using state-of-the-art methods; the black dotted line indicates the position of the chemical potential.</figcaption></figure>

# Appendix A: Conventions, definitions and derivations

## 1. Fourier transforms

We assume Born-von Karman cyclic boundary conditions and consider a supercell made up of N primitive cells, whose lattice vectors are $T_1 = N_1t_1$, $T_2 = N_2t_2$ with $N_1N_2 = N$, where $t_1, t_2$ are the direct lattice vectors (and $g_1, g_2$ the reciprocal ones). We define the Fourier transform in the Born-von Karman supercell as

$$f(\mathbf{q} + \mathbf{G}) = \int f(\mathbf{r})e^{-i(\mathbf{q}+\mathbf{G})\cdot\mathbf{r}}d\mathbf{r}, \quad (A1)$$

$$f(\mathbf{r}) = \frac{1}{NA} \sum_{\mathbf{q}\mathbf{G}} f(\mathbf{q} + \mathbf{G})e^{i(\mathbf{q}+\mathbf{G})\cdot\mathbf{r}}, \quad (A2)$$

where A is the area of the unit cell, the integration is intended to run over the Born-von Karman supercell and $\mathbf{q} = m_1/N_1\mathbf{g}_1 + m_2/N_2\mathbf{g}_2$ where $m_1, m_2$ are integer numbers, while $\mathbf{G}$ are defined as reciprocal lattice vectors. In the case of quantities which are cell-periodic, the above expression is evaluated at $\mathbf{q} = 0$, and the domain of integration is reduced to the primitive cell while taking $N = 1$. Analogously, the Fourier transform of quantities

dependent on two real space indexes is written as

$$f(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}') = \frac{1}{NA} \int d\mathbf{r} d\mathbf{r}' f(\mathbf{r}, \mathbf{r}') \times e^{-i[(\mathbf{q}+\mathbf{G})\cdot\mathbf{r}-(\mathbf{q}+\mathbf{G}')\cdot\mathbf{r}]}, \quad (A3)$$

$$f(\mathbf{r}, \mathbf{r}') = \frac{1}{NA} \sum_{\mathbf{q}\mathbf{G}\mathbf{G}'} f(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}') \times e^{i[(\mathbf{q}+\mathbf{G})\cdot\mathbf{r}-(\mathbf{q}+\mathbf{G}')\cdot\mathbf{r}]}, \quad (A4)$$

where we have used that $\mathbf{q}$ is the same in both reciprocal space arguments as a consequence of $f(\mathbf{r}, \mathbf{r}') = f(\mathbf{r} + \mathbf{R}, \mathbf{r}' + \mathbf{R})$ for any $\mathbf{R}$ belonging to the direct lattice— when $\mathbf{G} = \mathbf{G}' = 0$ we will often shorten $f(\mathbf{q}, \mathbf{q}) := f(\mathbf{q})$. Since we need to transform also along the out-of-plane direction to obtain the expression for the Coulomb kernel of Eq. 6, we define

$$f(\mathbf{q} + \mathbf{G}, q_z) = \int f(\mathbf{r}, z)e^{-i(\mathbf{q}+\mathbf{G})\cdot\mathbf{r}}e^{-iq_z z}d\mathbf{r}dz, \quad (A5)$$

$$f(\mathbf{r}, z) = \frac{1}{2\pi NA} \int \sum_{\mathbf{q}\mathbf{G}} f(\mathbf{q} + \mathbf{G}, q_z) e^{i(\mathbf{q}+\mathbf{G})\cdot\mathbf{r}} e^{iq_z z} dq_z. \quad (A6)$$

With these definitions, the 3d transform of the Coulomb kernel is

$$v(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}', q_z) = \frac{4\pi e^2}{|\mathbf{q} + \mathbf{G}|^2 + q_z^2} \delta_{\mathbf{G}\mathbf{G}'}, \quad (A7)$$
from which Eq. 6 follows using the residue theorem on the antitransform to the z variable.

For the transform of the force constants matrix, we define

$$C_{ss',\alpha\beta}(\mathbf{q}) = \sum_{p} C_{ss',\alpha\beta}(\mathbf{R}_p) e^{i\mathbf{q} \cdot (\mathbf{R}_p + \boldsymbol{\tau}_s - \boldsymbol{\tau}_{s'})}, \\ C_{ss',\alpha\beta}(\mathbf{R}_p) = \frac{1}{N} \sum_{\mathbf{q}} C_{ss',\alpha\beta}(\mathbf{q}) e^{-i\mathbf{q} \cdot (\mathbf{R}_p + \boldsymbol{\tau}_s - \boldsymbol{\tau}_{s'})}, \quad (\text{A8})$$

where $\mathbf{R}_p$ is the vector that identifies the cell $p$ and $\boldsymbol{\tau}_s, \boldsymbol{\tau}_{s'}$ are the basis vectors of the atoms $s, s'$ in the unit cell. Accordingly, a phonon of polarization $e_{s,\alpha}^{\nu}(\mathbf{q})$ induces a displacement $v$ of the atom $s$ along the direction $\alpha$ in the cell $p$ that can be written as

$$v_{sp,\alpha}^{\nu}(\mathbf{q}) = e_{s,\alpha}^{\nu}(\mathbf{q}) e^{i\mathbf{q} \cdot (\mathbf{R}_p + \boldsymbol{\tau}_s)}. \quad (\text{A9})$$

The above definition differentiates from the one of Refs. [25, 71, 91, 92] for the explicit presence of $\boldsymbol{\tau}_s$ at the exponential, while it is the same of Refs. [67, 68]—the advantage of this definition is that the response to such monochromatic perturbation behaves as a scalar under a reference frame change. We also define

$$v_{s,\alpha}^{\nu}(\mathbf{q}) = e_{s,\alpha}^{\nu}(\mathbf{q}) \sum_{p} e^{i\mathbf{q} \cdot (\mathbf{R}_p + \boldsymbol{\tau}_s)}. \quad (\text{A10})$$

## 2. Definition of charge density change

In general, for the monochromatic linear response problem we may write the external perturbation on the atoms as [67] $v_{s,\alpha}(\mathbf{q}) = \sum_{p} \lambda_{\alpha} e^{i\mathbf{q} \cdot (\mathbf{R}_p + \boldsymbol{\tau}_s)}$ where $\lambda$ is the dimensional amplitude of the atomic displacement. For a generic charge density or potential $f$ that depends on the set of atomic positions $\{\mathbf{R}_p + \boldsymbol{\tau}_s\}$, we define its variation in the linear response regime with respect to the external perturbation $v_{s,\alpha}(\mathbf{q})$ as

$$\left. \frac{\partial f(\mathbf{r}, \{\mathbf{R}_p + \boldsymbol{\tau}_s\})}{\partial \lambda_{\alpha}} \right|_{\lambda_{\alpha}=0} = \sum_{\mathbf{G}} \delta f_{s,\alpha}(\mathbf{q} + \mathbf{G}, z, \{\mathbf{R}_p + \boldsymbol{\tau}_s\}) e^{i(\mathbf{q} + \mathbf{G}) \cdot \mathbf{r}}; \quad (\text{A11})$$

$\sum_{\mathbf{G}} \delta f_{s,\alpha}(\mathbf{q} + \mathbf{G}, z, \{\mathbf{R}_p + \boldsymbol{\tau}_s\}) e^{i\mathbf{G} \cdot \mathbf{r}}$ is clearly a cell-periodic function and, using a notation similar to the one of Ref. [23], may be as well indicated with $f^{\mathbf{q}}(\mathbf{r})$; throughout the text we shorten the notation by defining $\delta f_{s,\alpha}(\mathbf{q} + \mathbf{G}, z) := \delta f_{s,\alpha}(\mathbf{q} + \mathbf{G}, z, \{\mathbf{R}_p + \boldsymbol{\tau}_s\})$. As an example, the ionic charge density change induced by $v_{s,\alpha}(\mathbf{q})$ is written, given the unperturbed expression for the ionic charge density $\rho^{\text{ion}}(\mathbf{r}) = \sum_{ps} \mathcal{Z}_s e \delta(\mathbf{r} - \mathbf{R}_p - \boldsymbol{\tau}_s)$, as

$$\left. \frac{\partial \rho^{\text{ion}}(\mathbf{r}, \{\mathbf{R}_p + \boldsymbol{\tau}_s\})}{\partial \lambda_{\alpha}} \right|_{\lambda_{\alpha}=0} = \sum_{\mathbf{G}} \delta \rho_{s,\alpha}^{\text{ion}}(\mathbf{q} + \mathbf{G}, z) e^{i(\mathbf{q} + \mathbf{G}) \cdot \mathbf{r}}. \quad (\text{A12})$$

Since $\delta \rho_{s,\alpha}^{\text{ion}}$ is itself the perturbing external charge density of the electrostatic problem, we rename it $\delta \rho_{s,\alpha}^{\text{ext}}$ for clarity. For a 2d material where all the atoms are disposed on one single layer and are subjected only to in-plane perturbations of the atoms, which is the matter of study of this work, we obtain

$$\delta \rho_{s,\alpha}^{\text{ext}}(\mathbf{q} + \mathbf{G}, z) = -i \frac{\mathcal{Z}_s e}{A} [(q_{\alpha} + G_{\alpha}) e^{-i\mathbf{G} \cdot \boldsymbol{\tau}_s}] \delta(z), \quad (\text{A13})$$

which, when layer-averaged as defined in Eq. A20, gives Eq. 19. The same definitions may be extended to the variation of potentials. Notice that in case of the use of pseudopotentials, following a collective atomic displacement the local part of the pseudopotential produces an external charge in a form similar Eq. A13, but which includes a form factor coming from the shape of the pseudopotential. The form factor is though relevantly different from zero for $q \gg \frac{1}{R_c}$, where $R_c$ is the core radius of the pseudopotential; similarly, the non-local part of the pseudopotential is irrelevant in the $q \ll \frac{1}{R_c}$ limit.

## 3. Response functions

We define the following linear response functions in real space:

$$\Delta \rho^{\text{ind}}(\mathbf{r}) = \int d\mathbf{r}' \chi(\mathbf{r}, \mathbf{r}') \Delta V^{\text{ext}}(\mathbf{r}'), \quad (\text{A14})$$

$$\Delta \rho^{\text{ind,H}}(\mathbf{r}) = \int d\mathbf{r}' \Pi(\mathbf{r}, \mathbf{r}') \Delta V^{\text{tot,H}}(\mathbf{r}'), \quad (\text{A15})$$

$$\Delta \rho^{\text{ind}}(\mathbf{r}) = \int d\mathbf{r}' \chi^0(\mathbf{r}, \mathbf{r}') \Delta V^{\text{tot}}(\mathbf{r}'), \quad (\text{A16})$$

$$\Delta V^{\text{tot}}(\mathbf{r}) = \int d\mathbf{r}' \epsilon^{-1}(\mathbf{r}, \mathbf{r}') \Delta V^{\text{ext}}(\mathbf{r}'), \quad (\text{A17})$$

where $\Delta V^{\text{ext}}$ is an external perturbing potential, $\Delta V^{\text{ind}}$ is the change of potential in the crystal induced by the external perturbation, $\Delta V^{\text{tot,H}} := \Delta V^{\text{ind,H}} + \Delta V^{\text{ext}}$ and $\Delta V^{\text{tot}} := \Delta V^{\text{ind}} + \Delta V^{\text{ext}}$ where the difference between $\Delta V^{\text{ind,H}}$ and $\Delta V^{\text{ext}}$ is respectively the neglect or the inclusion of the exchange correlation terms in the response (we use the same superscript meaning for the charge density $\rho$) [93]. Notice that $\chi^0$ is the IPP of the Khon-Sham system and is expressed in Eq. 8 with wavefunctions that have been computed with the inclusion of the $xc$ potential [63], while $\Pi$ would not include the $xc$ contributions in the self consistent field calculation; in our work we use the RPA approximation in which the total response $\chi$ is approximated as

$$\chi^{-1} = \chi^{0,-1} - v - K_{xc} \rightarrow \chi^{-1,\text{RPA}} = \chi^{0,-1} - v; \quad (\text{A18})$$

from this approximation and $\epsilon^{-1,\text{RPA}} = \mathcal{I} + v\chi^{\text{RPA}}$, it follows that $\epsilon^{-1,\text{RPA}}\chi^{-1,\text{RPA}}\chi^0 = \mathcal{I}$ which means that $\epsilon^{\text{RPA}} = \mathcal{I} - v\chi^0$, i.e. Eq. 10—see App. D for considerations beyond RPA. The definitions from Eq. A14 to
Eq. A17 are given for a general real space perturbation, and therefore in reciprocal space they can be used with the substitution $\Delta \to \delta$ where $\delta$ assumes the meaning given in App. A2. In this case, to connect the external perturbation to the external charge change, we use the electrostatic relation

$$ \delta V_{s,\alpha}^{\text{ext}}(\mathbf{q} + \mathbf{G}, z) = 2\pi \int dz' \frac{e^{-|\mathbf{q}+\mathbf{G}||z-z'|}}{|\mathbf{q}+\mathbf{G}|} \delta\rho_{s,\alpha}^{\text{ext}}(\mathbf{q} + \mathbf{G}, z'). \quad (A19) $$

Also, to our aims it is useful to define the averaging of a generic function $f(\star, z, z')$ over the layer width $t$ (defined in Sec. II B) as

$$ \tilde{f}(\star) = \int_{-\frac{t}{2}}^{\frac{t}{2}} dz \frac{dz'}{t} f(\star, z, z'), \quad (A20) $$

which we use to get rid of the $z, z'$ dependence of, e.g., the response functions. With the above definition, the layer-averaged response functions satisfy

$$ \int dz'' \sum_{\mathbf{G}''} \tilde{f}(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}'') \times \quad (A21) \\ \tilde{f}^{-1}(\mathbf{q} + \mathbf{G}'', \mathbf{q} + \mathbf{G}') = \delta_{\mathbf{G}\mathbf{G}''}. $$

Notice that, from a more formal perspective, the layer-averaging procedure may be interpreted as keeping the leading order of the hyperbolic cosine in the appropriate expressions of Ref. [40].

### 4. Definition of $\Delta_{q\nu}$

$\Delta_{q\nu}$ is defined as in Ref. [5], but with phase convention coherent with Eq. A10, i.e. as

$$ \Delta_{q\nu} V = e^{i\mathbf{q}\cdot\mathbf{r}} \Delta_{q\nu} v, \quad (A22) $$

$$ \Delta_{q\nu} v = l_{q\nu} \sum_{s\alpha} \left(\frac{M_0}{M_s}\right)^{\frac{1}{2}} e^{\nu s\alpha} (\mathbf{q}) \partial_{s\alpha, \mathbf{q}} v, \quad (A23) $$

$$ \partial_{s\alpha, \mathbf{q}} v = \sum_{p} e^{-i\mathbf{q} \cdot (\mathbf{r} - \mathbf{R}_p - \tau_s)} \frac{\partial V}{\partial \tau_{s\alpha}} \Big|_{\mathbf{r} - \mathbf{R}_p - \tau_s}, \quad (A24) $$

where $\Delta_{q\nu}v$ is a cell-periodic function and again we have restricted the derivation to be performed wit respect to in-plane displacements of the atoms. We also define

$$ l_{q\nu} = \left[ \frac{\hbar}{2M_0 \omega_{q,\nu}} \right]^{1/2}. \quad (A25) $$

Notice that

### 5. Independent particle polarizability and dielectric matrix

Starting from the approximation of Eq. 11, Eq. 8 becomes Eq. 12:

$$ \begin{aligned} \chi^0(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}', z, z') &= \frac{2e^2}{At^2} \sum_{mm'k} \frac{f_{mk} - f_{m'k+\mathbf{q}}}{\epsilon_{mk} - \epsilon_{m'k+\mathbf{q}}} \times \\ & \langle u_{m'k+\mathbf{q}+\mathbf{G}} | u_{mk} \rangle \langle u_{mk} | u_{m'k+\mathbf{q}+\mathbf{G}'} \rangle \theta\left(\frac{t}{2} - |z|\right) \theta\left(\frac{t}{2} - |z'|\right) \\ &= \frac{1}{t^2} \theta\left(\frac{t}{2} - |z|\right) \theta\left(\frac{t}{2} - |z'|\right) \chi^0(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}'). \end{aligned} \quad (A26) $$

The brakets stand for integration on the in-plane real space variables and the expression $\chi^0(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}')$, as already anticipated in Sec. II, is the IPP of a two dimensional Khon-Sham system [64, 65]. To obtain an expression for the dielectric tensor, we plug Eq. A26 inside Eq. 10 obtaining

$$ \begin{aligned} \epsilon(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}', z, z') = \delta(z - z') \delta_{\mathbf{G}\mathbf{G}'} + & \quad (A27) \\ - \frac{2\pi}{t^2} \frac{P(\mathbf{q} + \mathbf{G}, z)}{|\mathbf{q} + \mathbf{G}|} \chi^0(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}') \theta\left(\frac{t}{2} - |z'|\right), & \end{aligned} $$

where we have defined

$$ P(\mathbf{q} + \mathbf{G}, z) = \frac{2}{|\mathbf{q} + \mathbf{G}|} \times \quad (A28) \\ \times \begin{cases} e^{-|\mathbf{q} + \mathbf{G}|z} \sinh\left(|\mathbf{q} + \mathbf{G}|\frac{t}{2}\right) & z > \frac{t}{2} \\ 1 - e^{-|\mathbf{q} + \mathbf{G}|\frac{t}{2}} \cosh\left(|\mathbf{q} + \mathbf{G}|z\right) & |z| \le \frac{t}{2} \\ e^{|\mathbf{q} + \mathbf{G}|z} \sinh\left(|\mathbf{q} + \mathbf{G}|\frac{t}{2}\right) & z < -\frac{t}{2} \end{cases}. $$

We concentrate on the domain $z, z' \in [-\frac{t}{2}, \frac{t}{2}]$ because it is the one where we are interested to evaluate $\epsilon^{-1}(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}', z, z')$. In this region Eqs. A27 and A28 can be further simplified using the definition of layer-averaged dielectric function of Eq. 13. Plugging Eqs. A27 and A28 in Eq. 13, we finally obtain Eq. 14 of the main text.

The appropriate long wavelength expansions of Eq. 14 are obtained analyzing the asymptotic behaviours of $\chi^0$. For insulators and undoped semiconductors at zero temperature we can write for the leading orders

$$ \lim_{\mathbf{q} \to 0} \chi^0(\mathbf{q}) = \frac{2q_\alpha q_\beta}{A} \sum_{mm'k} \frac{\theta(\epsilon_{m'k}) - \theta(\epsilon_{mk})}{\epsilon_{mk} - \epsilon_{m'k}} \times \quad (A29) $$

$$ \langle u_{mk} | \partial_{k\alpha} u_{m'k} \rangle^* \langle u_{mk} | \partial_{k\beta} u_{m'k} \rangle := \mathbf{q} \cdot \mathbf{A}(\mathbf{G}'), \quad (A30) $$

$$ \lim_{\mathbf{q} \to 0} \chi^0(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}') = \frac{2}{A} \sum_{mm'k} \frac{\theta(\epsilon_{m'k}) - \theta(\epsilon_{mk})}{\epsilon_{mk} - \epsilon_{m'k}} \times \quad (A31) $$

where it is intended that Eq. A30 is valid for $\mathbf{G}' \neq 0$ and Eq. A31 for $\mathbf{G} \neq 0$ and $\mathbf{G}' \neq 0$. To obtain the
above expressions we have used that $\langle u_{mk} | u_{m'k'} \rangle = 0$ if $m$ is a valence state and $m'$ a conduction state (or viceversa) and that there are no intraband transitions that can vanish the energy denominator, so that the $q$ dependence of the limits is coming only from the long-wavelength expansion of the periodic part of the Bloch functions. For the dielectric matrix we find

$$ \epsilon(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}') = \frac{1}{t} \begin{cases} 1 - v(q)\mathbf{q} \cdot \mathbf{B} \cdot \mathbf{q} & \mathbf{G}, \mathbf{G}' = 0 \\ -v(q)\mathbf{q} \cdot \mathbf{A}^{\text{c.c.}}(\mathbf{G}') & \mathbf{G} = 0 \\ -v(G)\mathbf{q} \cdot \mathbf{A}(\mathbf{G}) & \mathbf{G}' = 0 \\ \delta_{\mathbf{G}\mathbf{G}'} - v(G)C(\mathbf{G}, \mathbf{G}') & \end{cases}, \quad (\text{A32}) $$

where $v$ is defined in Eq. 15. For metals or doped semiconductors the asymptotic limits for the IPP instead read

$$ \lim_{\mathbf{q} \to 0} \chi^0(\mathbf{q}) = \frac{2}{A} \sum_{m\mathbf{k}} \frac{\delta f_{m\mathbf{k}} - \delta f_{m\mathbf{k}+\mathbf{q}}}{\epsilon_{m\mathbf{k}} - \epsilon_{m\mathbf{k}+\mathbf{q}}} + \mathbf{q}B\mathbf{q}, \quad (\text{A33}) $$

$$ \lim_{\mathbf{q} \to 0} \chi^0(\mathbf{q}, \mathbf{q} + \mathbf{G}') = \frac{2}{A} \sum_{m\mathbf{k}} \frac{\delta f_{m\mathbf{k}} - \delta f_{m\mathbf{k}+\mathbf{q}}}{\epsilon_{m\mathbf{k}} - \epsilon_{m\mathbf{k}+\mathbf{q}}} \times \langle u_{m\mathbf{k}} | u_{m\mathbf{k}+\mathbf{G}'} \rangle + \mathbf{q} \cdot \mathbf{A}(\mathbf{G}'), \quad (\text{A34}) $$

$$ \lim_{\mathbf{q} \to 0} \chi^0(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}') = \frac{2}{A} \sum_{mm'\mathbf{k}} \frac{\delta f_{m\mathbf{k}} - \delta f_{m'\mathbf{k}+\mathbf{q}}}{\epsilon_{m\mathbf{k}} - \epsilon_{m'\mathbf{k}+\mathbf{q}}} \times \langle u_{m\mathbf{k}} | u_{m'\mathbf{k}+\mathbf{G}'} \rangle^{\text{c.c.}} \langle u_{m\mathbf{k}} | u_{m'\mathbf{k}+\mathbf{G}'} \rangle + C(\mathbf{G}, \mathbf{G}'), \quad (\text{A35}) $$

where we have used the approximation that $\langle u_{m\mathbf{k}} | u_{m'\mathbf{k}+\mathbf{q}} \rangle = \delta_{mm'}$ and therefore the sum over $m$ is restricted to the states whose occupation is significantly different from 1 or 0. In this case the response functions depend on the carrier concentration $n$ and temperature $T$ through the occupation factors, i.e. $\epsilon(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}') \to \epsilon(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}', n, T)$. Using a purely qualitative argument, the expression of the dielectric matrix as a function of $n$ can be obtained in the degenerate limit performing the substitution $\mathbf{q} \cdot \mathbf{B} \cdot \mathbf{q} \to k_{TF}^2(n) + \mathbf{q} \cdot \mathbf{B} \cdot \mathbf{q}$ in Eq. A32, where $k_{TF}^2(n)$ is the Thomas-Fermi screening wavevector of the material, or using $\mathbf{q} \cdot \mathbf{B} \cdot \mathbf{q} \to k_D^2(n, T) + \mathbf{q} \cdot \mathbf{B} \cdot \mathbf{q}$ in the non-degenerate limit, where $k_D(n, T)$ is the Debye screening wavevector [94]; we are here supposing that the wings of the dielectric response matrix are not affected by doping—as shown in Sec. IV this is not rigorously true, but as a matter of fact wings are less sensitive to doping than the head. The above substitution cures the non-analiticity of the LRCs of the dynamical matrix and of the EPI typical of the insulator case (see Secs. IID1 and IID2), smoothing their singular behaviours in a small region around $\Gamma$. Quantitatively, the above description of the doping effect on the screening is too rough since this is strongly dependent on the microscopic details of the material and temperature, and therefore more refined strategies, as the one proposed in Eq. 41, have to be implemented.

## 6. Inversion of $w^{-1}$

We define the inverse screened Coulomb potential $w^{-1}$ as the inverse of Eq. 17:

$$ w^{-1}(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}') = \frac{\epsilon(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}')}{v(\mathbf{q} + \mathbf{G})}. \quad (\text{A36}) $$

The $w^{-1}$ matrix is Hermitian because of the Hermiticity of $\chi^0$ and Eq. 14. For a generic Hermitian matrix $w^{-1}$, following [67, 68], we write:

$$ w^{-1} = \begin{pmatrix} P & Q \\ Q^{\dagger} & S \end{pmatrix}, \quad (\text{A37}) $$

$$ w = \begin{pmatrix} W & X \\ X^{\dagger} & Z \end{pmatrix}, \quad (\text{A38}) $$

where $P = w^{-1}(\mathbf{q})$ is the head of the matrix, $Q = w^{-1}(\mathbf{q}, \mathbf{q} + \mathbf{G}')$ is the wing and $S = S^{\dagger} = w^{-1}(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}')$ is the body, and the same goes for $W, X, Z$. We have

$$ W = (P - QS^{-1}Q^{\dagger})^{-1}, \quad (\text{A39}) $$

$$ X = -WQS^{-1}, \quad (\text{A40}) $$

$$ Z = S^{-1} + X^{\dagger}W^{-1}X. \quad (\text{A41}) $$

$w$ may then be rewritten as:

$$ w = \begin{pmatrix} W & X \\ X^{\dagger} & S^{-1} + X^{\dagger}W^{-1}X \end{pmatrix}. \quad (\text{A42}) $$

We also define

$$ \hat{w} = \begin{pmatrix} W & X \\ X^{\dagger} & X^{\dagger}W^{-1}X \end{pmatrix}, \quad (\text{A43}) $$

which can be rewritten as

$$ \hat{w}(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}') = \frac{\hat{w}(\mathbf{q}, \mathbf{q} + \mathbf{G}')\hat{w}(\mathbf{q} + \mathbf{G}, \mathbf{q})}{\hat{w}(\mathbf{q})} = \frac{\hat{w}(\mathbf{q}, \mathbf{q} + \mathbf{G}')\hat{w}^{\text{c.c.}}(\mathbf{q}, \mathbf{q} + \mathbf{G})}{\hat{w}(\mathbf{q})}. \quad (\text{A44}) $$

For the case of undoped semiconductors, the $\hat{w}$ tensor contains all the non-analytical terms that give rise to the LRCs of the EPI and of the dynamical matrix; notice that the head and the wings of the $w$ and the $\hat{w}$ tensors coincide.

### a. Asymptotic expansion of the $w$ tensor

To evaluate Eq. 22, we need to know the asymptotic expressions of the $w$ tensor, given the inversion formulae of App. A6. For insulators and undoped semiconductors
these are

$$w(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}') = \begin{cases} w(\mathbf{q}) & \mathbf{G}, \mathbf{G}' = 0 \\ w(\mathbf{q}) [\mathbf{q} \cdot \mathbf{A}^{\prime \text{c.c.}}(\mathbf{G}')] & \mathbf{G} = 0 \\ w(\mathbf{q}) [\mathbf{q} \cdot \mathbf{A}'(\mathbf{G})] & \mathbf{G}' = 0 \\ L(\mathbf{G}, \mathbf{G}') + S^{-1}(\mathbf{G}, \mathbf{G}') & \end{cases} \quad (A45)$$

$$w(\mathbf{q}) = \begin{cases} \frac{2\pi}{|\mathbf{q}+\mathbf{q}' \cdot \mathbf{q}|} & \text{thin} \\ \frac{4\pi}{t \mathbf{q} \cdot \mathbf{B}'' \cdot \mathbf{q}} & \text{thick} \end{cases}$$

where c.c. stands for complex conjugate and the expressions for $\mathbf{A}'$, $L$, $H$, which are different for the thin and thick limits, can be easily worked out. For metals or doped semiconductors, the above expressions have to be generalized to include the terms coming from the expansions Eqs. A33, A34 and A35. The generalization can be easily worked out and is not reported here. We though notice that if we are just interested in the ratio between the wing and the head of the $w$ entering Eq. 26 and Eq. 22, then we can recast it as

$$\frac{w(\mathbf{q}, \mathbf{q} + \mathbf{G}', n, T)}{w(\mathbf{q}, \mathbf{q}, n, T)} = \begin{cases} 1 & \mathbf{G}' = 0 \\ \mathbf{q} \cdot \mathbf{A}^{\prime \text{c.c.}}(\mathbf{G}') + \mathcal{C}(\mathbf{q}, n, T) & \mathbf{G}' \neq 0 \end{cases}, \quad (A46)$$

where the term $\mathcal{C}(\mathbf{q}, n, T)$ stems formally from the change of $\chi^0$ due to (marginally) the modification of the interband terms and to (substantially) the appearance of an intraband contribution to $\chi^0$ in presence of doping at finite temperature (see Eqs. A33 and A34 for the limiting values); practically, we associate the presence of $\mathcal{C}(\mathbf{q}, n, T)$ exclusively to the existence of intraband contributions. Such term can be detected in the density response in the region where the dielectric response acquires strong metallic features (see Sec. IV).

## 7. Proof of Eq. 24

To prove the expansion of Eq. 24, we consider the electrostatic problem imposing that the $\mathbf{G} = 0$ component of the change of the total potential, $\delta V_{s,\alpha}^{\text{tot}}(\mathbf{q})$, is null for any possible charge density perturbation—notice that we are considering the electrostatic problem where each term of the Maxwell’s equation has been layer-averaged independently (see also App. C). This request can be satisfied only if in Eq. A38 we have $W = \eta$ and $X = \eta$ where $\eta$ is an infinitesimally small number. We notice that $\lim_{\eta \to 0} \hat{w} = 0$, so that we are left only with the short range components of the $w$ matrix. The total potential may now be written, for $\mathbf{G} \neq 0$ as

$$\delta \bar{V}_{s,\alpha}^{\text{tot}}(\mathbf{q} + \mathbf{G}) = \sum_{\mathbf{G}' \neq 0} S^{-1}(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}') t \delta \rho_{s,\alpha}^{\text{ext}}(\mathbf{q} + \mathbf{G}'); \quad (A47)$$

the bar over a quantity is here used to indicate that it is computed imposing $\delta V_{s,\alpha}^{\text{tot}}(\mathbf{q}) = 0$—the bar notation for $\bar{Z}_{s,\alpha}$ introduced in Sec. II is not casual, as we will see in a moment. By the definition of $X^0$ we obtain

$$\begin{aligned} t \delta \bar{\rho}_{s,\alpha}^{\text{ind}}(\mathbf{q}) &= \sum_{\substack{\mathbf{G} \neq 0 \\ \mathbf{G}' \neq 0}} \chi^0(\mathbf{q}, \mathbf{q} + \mathbf{G}) S^{-1}(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}') \times \\ & \qquad t \delta \bar{\rho}_{s,\alpha}^{\text{ext}}(\mathbf{q} + \mathbf{G}'). \end{aligned} \quad (A48)$$

To obtain the total charge density change in this particular setup, $\delta \bar{\rho}_{s,\alpha}^{\text{tot}}(\mathbf{q})$, we sum $\delta \rho_{s,\alpha}^{\text{ext}}(\mathbf{q})$ to $\delta \bar{\rho}_{s,\alpha}^{\text{ind}}(\mathbf{q})$. Using the asymptotic expressions for $w$, $\chi^0$ and $S^{-1}$ (for the first see App. A6a, for the second see App. A5, while for the third we just get $S^{-1}(\mathbf{G}, \mathbf{G}')$) both in the thin and thick limits one finds that $t \delta \bar{\rho}_{s,\alpha}^{\text{tot}}(\mathbf{q})$ is exactly equivalent to the asymptotic expansion of the left hand side of Eq. 22, i.e.

$$\delta \bar{\rho}_{s,\alpha}^{\text{tot}}(\mathbf{q}) = -i \frac{eq}{At} \bar{Z}_{s,\alpha}(\mathbf{q}). \quad (A49)$$

This means that $-i \frac{eq}{A} \bar{Z}_{s,\alpha}(\mathbf{q})$ corresponds to the Fourier transform of the total charge density change generated by $t \delta \rho_{s,\alpha}^{\text{ext}}(\mathbf{q})$ once that we have imposed the absence of macroscopic electric fields, i.e. $t \delta \bar{\rho}_{s,\beta}^{\text{tot}}(\mathbf{q})$—exactly as in the three dimensional case studied in Ref. [36]; the reason for the use of the bar notation for the unscreened effective charges is now evident. Eq. A48 also shows us a fundamental property: $t \delta \bar{\rho}_{s,\beta}^{\text{tot}}(\mathbf{q})$ is manifestly analytic and as such may be expanded in a Taylor series, therefore justifying Eq. 24.

From Eqs. A49,22, 2, 3 and the RPA electrostatic relation between charges and potentials we deduce Eqs. 21, from which the definition of $Z_{s,\alpha}(\mathbf{q})$ starting from $\bar{Z}_{s,\alpha}(\mathbf{q})$ (Eq. 3) is now naturally clear, and is coherent with Eq. 2. Notice that the relation between screened and unscreened charges, explicitly deduced here for the first order of the expansion, is valid at all orders, as shown by the alternative derivation given in App. D.

The above argument is fundamental to prove that Eq. 24 contains the in-plane dynamical effective charges. For example, the Born effective charge is defined as in Eq. 42, which in our notation (see App. A2) is equivalent to write

$$Z_{s,\alpha\beta}^* = \left. \frac{A}{e} \frac{\partial P_{\alpha}}{\partial \lambda_{\beta}} \right|_{\lambda=0, \mathbf{q}=0, \mathbf{E}=0}; \quad (A50)$$

notice that the star notation here does not mean complex conjugate, but has always been historically used to identify the Born effective charges, that are properly real quantities. We now define the cell-periodic layer-averaged charge density change on a single atom (SA) exploiting the superposition of the charge density change in the linear regime for the response

$$\delta \bar{\rho}_{s,\beta}^{\text{SA}}(\mathbf{r} - \mathbf{R}_p - \tau_s) = \frac{A}{(2\pi)^2} \int d\mathbf{q} \delta \bar{\rho}_{s,\beta}^{\text{q,tot}}(\mathbf{r}) e^{-i\mathbf{q} \cdot (\mathbf{R}_p - \tau_s - \mathbf{r})}, \quad (A51)$$
with the inverse transform given by

$$ \delta\bar{\rho}_{s,\beta}^{\text{q,tot}}(\mathbf{r}) = \sum_{p} \delta\bar{\rho}_{s,\beta}^{\text{SA}}(\mathbf{r} - \mathbf{R}_{p} - \boldsymbol{\tau}_{s})e^{i\mathbf{q}\cdot(\mathbf{R}_{p}+\boldsymbol{\tau}_{s}-\mathbf{r})}. \quad (\text{A52}) $$

Then, we write

$$ \begin{gathered} \left. \frac{\partial P_{\alpha}}{\partial \lambda_{\beta}} \right|_{\lambda=0, \mathbf{q}=0, \mathbf{E}=0} = \frac{1}{NA} \sum_{p} \int d\mathbf{r} (r_{\alpha} - R_{p,\alpha} - \tau_{s,\alpha}) \times \\ \delta\bar{\rho}_{s,\beta}^{\text{SA}}(\mathbf{r} - \mathbf{R}_{p} - \boldsymbol{\tau}_{s}) = i \left. \frac{\partial}{\partial q_{\alpha}} \delta\bar{\rho}_{s,\beta}^{\text{tot}}(\mathbf{q}) \right|_{\mathbf{q}=0}, \end{gathered} \quad (\text{A53}) $$

where the integration is intended to run over the whole crystal. The last equality of Eq. A53 proves the identification between effective charges and the expansion of the unscreened charge density change; the same can be done also for higher orders. We need to stress here that these identifications are exact only in the thin and thick limits, i.e. in the regimes in which the LRCs can be matched to phenomenological expressions that involve only effective charges and the head of the inverse screening matrix. In these cases, these identifications are crucial in order to have an operative simple method to extract from ab-initio calculations the value of the dynamical effective charges without passing from the explicit calculation of all the structure of the w tensor, which may depend on the adopted model.

## 8. Derivation of Eq. 26

For a single layer we can write for the dynamical matrix in reciprocal space

$$ C_{ss', \alpha \beta}(\mathbf{q}) = \frac{2\pi Z_s Z_{s'} e^2}{A} \sum_{\mathbf{G} \mathbf{G}'} e^{i\mathbf{G} \cdot \boldsymbol{\tau}_s - i\mathbf{G}' \cdot \boldsymbol{\tau}_{s'}} \int_{-\infty}^{\infty} dz (q_{\alpha} + G_{\alpha}) (q_{\beta} + G'_{\beta}) \epsilon^{-1}(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}', 0, z) \frac{e^{-|\mathbf{q} + \mathbf{G}'||z|}}{|\mathbf{q} + \mathbf{G}'|} \quad (\text{A54}) $$

To obtain asymptotic formulae, we now wish to get rid of the z dependence of the above equation. To do so, we replace $\epsilon^{-1}(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}', 0, z)$ with its layer average contextually restricting $z \in [-\frac{t}{2}, \frac{t}{2}]$. In fact, we can approximate $\epsilon^{-1}$ assuming that the induced density is uniform inside the slab of width $t$, and vanishes outside, obtaining an expression similar to Eq. A27—see App. B for a discussion on this approximation. We also define

$$ Q(\mathbf{q} + \mathbf{G}) = \frac{4\pi}{|\mathbf{q} + \mathbf{G}|^2 t} \left(1 - e^{-|\mathbf{q} + \mathbf{G}| \frac{t}{2}}\right), \quad (\text{A55}) $$

which let us rewrite Eq. A54 as

$$ C_{ss', \alpha \beta}(\mathbf{q}) = \frac{Z_s Z_{s'} e^2}{A} \sum_{\mathbf{G} \mathbf{G}'} (q_{\alpha} + G_{\alpha}) (q_{\beta} + G'_{\beta}) \times \\ w(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}') \frac{Q(\mathbf{q} + \mathbf{G}')}{v(\mathbf{q} + \mathbf{G}')^2} e^{i\mathbf{G} \cdot \boldsymbol{\tau}_s - i\mathbf{G}' \cdot \boldsymbol{\tau}_{s'}}. \quad (\text{A56}) $$

In the thin and thick limits we have $\frac{Q(|\mathbf{q} + \mathbf{G}'|)}{v(\mathbf{q} + \mathbf{G}')^2} \to 1$; since we are interested in these limits we can set the ratio of the two functions to 1, obtaining eventually Eq. 26 when isolating the non analytical terms coming from the w tensor. For intermediate regimes the behaviour of the dynamical matrix (and EPI) will be, in general, strongly model dependent and therefore the division into short and long range components is not well-defined. The evaluation of Eq. A56 passes through the knowledge of the w matrix; the inversion of the $w^{-1}$ matrix can be performed symbolically, in line with Refs. [67, 68], as already shown in App. A6.

## Appendix B: Layer thickness and layer-averaging

The physical working hypothesis assumed in Sec. II is that, in appropriate limits, we can approximate a quasi two-dimensional system as a electronically compact homogeneous layer of width $t$ along the $\hat{z}$ direction, following Eq. 11. More specifically, we assume that whenever a $z$-dependent quantity enters an exact relation, we can substitute its value with its layer-average given by Eq. A20. The validity of this approximation is debatable, and indeed we show in Figs. 10 and 11 that the unperturbed electronic density along the $\hat{z}$ direction is not a window function, but it decays exponentially on a certain characteristic length, in contradiction with what assumed in Eq. 11. Nonetheless, we now argue that even in the exact treatment of the $z, z'$ dependence of our problem, we would obtain the same closed forms of the long-range components involving in-plane effective charges as if starting from Eq. 11, with the only assumption that there exists a scale $\hat{t}$ which can be used to define the thin and thick limits.

We start from the expression of the dielectric tensor

$$ \epsilon(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}', z, z', n, T) = \delta(z - z')\delta_{\mathbf{G} \mathbf{G}'} + \quad (\text{B1}) \\ -2\pi \int dz'' \frac{e^{-|\mathbf{q} + \mathbf{G}'||z-z''|}}{|\mathbf{q} + \mathbf{G}'|} \chi^0(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}', z'', z', n, T). $$

In the thin limit, we can define a scale $\hat{t}$ corresponding to the typical size of the extension of the induced charge along the out-of-plane direction; with this definition, we can approximate $e^{-|\mathbf{q} + \mathbf{G}'||z-z'|} \sim 1$ for $|\mathbf{q} + \mathbf{G}'|\hat{t} \ll 1$, to obtain

$$ \epsilon(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}', z, z', n, T) \sim \delta(z - z')\delta_{\mathbf{G} \mathbf{G}'} + \quad (\text{B2}) \\ - \frac{2\pi}{|\mathbf{q} + \mathbf{G}'|} \int dz'' \chi^0(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}', z'', z', n, T). $$

A similar relation holds between $\epsilon^{-1}$ and $\chi^{-1}$. Now, since to get to the formulae of the LRC of the dynamical matrix all the relevant approximations involve products of $\epsilon^{-1}$ with exponentials in the $z, z'$ coordinates, that can be approximated to unity, the final form of the LRC depends only on the layer-averaged value of the inverse dielectric matrix that can be defined from Eq. B2 using Eq. A21
over the length $\hat{t}$. The asymptotic formal expressions of the LRC of the dynamical matrix are therefore equal to the one of our model, since the layer-averaging does not introduce in general changes to the asymptotic limits of the inverse dielectric matrix; of course, this is true only if stick to the in-plane components of the expansions, otherwise further terms coming from the out-of-plane direction may be present as shown in Ref. [40]. For the LRC of the EPI modifications of the asymptotic formulae may come from the out-of-plane integration of the product of the inverse dielectric matrix with the periodic parts of the Bloch functions at $\mathbf{k}$ and $\mathbf{k}+\mathbf{q}$. Nonetheless, in the general case we expect such integration to tend to a constant non zero value in the limit $\mathbf{q} \to 0$ and therefore to recover the form of the LRCs of our model even in this case. Even if the form of the asymptotic models are the same, one would need to obtain the correct prefactors coming from the different kinds of approximation: we overcome this obstacle in our work by taking the values of the effective charges and of the dielectric constants directly from ab-initio calculations.

Analogously, for $|\mathbf{q} + \mathbf{G}|\hat{t} \gg 1$ we can use that the integral of Eq. B1 is relevantly different from zero only when in the IIP we have $z'' \approx z$:

$$ \epsilon(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}', z, z') \sim \delta(z - z') \delta_{\mathbf{G} \mathbf{G}'} + \quad (B3) $$

$$ -2\pi \frac{P(\mathbf{q} + \mathbf{G}, z)}{|\mathbf{q} + \mathbf{G}|} \chi^0(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}', z, z') \sim $$

$$ \delta(z - z') \delta_{\mathbf{G} \mathbf{G}'} - 4\pi \frac{1}{|\mathbf{q} + \mathbf{G}|^2} \chi^0(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}', z, z'). \quad (B4) $$

In the derivation of the LRCs for every occurrence of the inverse dielectric function under the integration sign together with $e^{-|\mathbf{q}+\mathbf{G}||z-z'|}$ we can use the same approximation of Eq. B3. Again, the asymptotic in-plane limits are not modified and the LRCs have the same formal expression as for our model.

From the above considerations it follows that the typical scale $\hat{t}$ shall be regarded as the scale where the charge response of the system is substantially different from zero.

## Appendix C: Layer-averaged Maxwell's equations

The Maxwell's equations are classical equations which relate macroscopic physical densities and currents to electromagnetic fields. In this appendix, we discuss the form of the first Maxwell's equation in the thin limit within the RPA approximation. We start from the RPA electrostatic relation between charge and potential

$$ \delta V^{\text{tot}}(\mathbf{q} + \mathbf{G}, z) = 2\pi \int dz' \frac{e^{-|\mathbf{q} + \mathbf{G}||z - z'|}}{|\mathbf{q} + \mathbf{G}|} \delta\rho^{\text{tot}}(\mathbf{q} + \mathbf{G}, z'). \quad (C1) $$

We now substitute each term of the above equation with its layer-average, as defined in App. A3 and discussed in

App. B; we have

$$ \delta V^{\text{tot}}(\mathbf{q} + \mathbf{G}) = v(\mathbf{q} + \mathbf{G})t\delta\rho^{\text{tot}}(\mathbf{q} + \mathbf{G}); \quad (C2) $$

if we now specify $\delta\rho^{\text{tot}} = \delta\rho^{\text{ind}} + \delta\rho^{\text{ext}}$ and use the relation of Eq. A15 between $\delta\rho^{\text{ind}}$ and $\delta V^{\text{tot}}$, we obtain

$$ \sum_{\mathbf{G}'} \epsilon(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}') \frac{\delta V^{\text{tot}}(\mathbf{q} + \mathbf{G}')}{v(\mathbf{q} + \mathbf{G})} = t \delta\rho^{\text{ext}}(\mathbf{q} + \mathbf{G}), \quad (C3) $$

which using Eq. A36 becomes

$$ \sum_{\mathbf{G}'} w^{-1}(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}') \delta V^{\text{tot}}(\mathbf{q} + \mathbf{G}') = t \delta \rho^{\text{ext}}(\mathbf{q} + \mathbf{G}). \quad (C4) $$

For the thin limit in particular we can write

$$ \sum_{\mathbf{G}'} \epsilon(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}') \delta V^{\text{tot}}(\mathbf{q} + \mathbf{G}') = 2\pi \frac{t \delta \rho^{\text{ext}}(\mathbf{q} + \mathbf{G})}{|\mathbf{q} + \mathbf{G}|}, \quad (C5) $$

where $t\delta\rho^{\text{ext}}(\mathbf{q}+\mathbf{G})$ is dimensionally a 2d density and the left hand side is not explicitly dependent on $t$. The above relation, which follows from Eqs. C1 and A16, relates the total potential change to the external charge perturbation, and includes all the in-plane non-local-fields component of the response.

The form of the classical Maxwell's equation instead follows from the relation between the total potential change and the macroscopically unscreened total charge change through the macroscopic inverse dielectric tensor, as shown in Sec. II D 1 and in particular in Eq. 22. Such equation may be rewritten for clarity as a function of the density as

$$ \delta V_{s,\alpha}^{\text{tot}}(\mathbf{q}) = w(\mathbf{q}) t \delta \rho_{s,\alpha}^{\text{tot}}(\mathbf{q}). \quad (C6) $$

Notice that if we compare the above expression with Eq. C2 evaluated at $\mathbf{G} = 0$ we obtain the relation of Eq. 21 between macroscopically screened and unscreened charges. Physically, this means that $\delta\bar{\rho}_{s,\alpha}^{\text{tot}}(\mathbf{q})$ acts electrostatically as an external macroscopic perturbing potential that has been renormalized by the local fields on the system under consideration.

## Appendix D: xc and response functions

The many body theoretical treatment of the response functions has been performed up to now directly in the RPA approximation, therefore neglecting any effect stemming from the presence of an exchange-correlation kernel $K_{xc}$ in the response. We will now try to extend the treatment beyond such an approximation. For the sake of clarity, we will restrict to the thick limit only and will take $t = 1$.
<figure><img src="image_11.png" /><figcaption>FIG. 10: Heat-maps of the in-plane electronic density of disproportionated graphene, expressed in atomic units, for several values of the z coordinate. The graphene layer is situated at z = 0. The heat-maps are taken at (upper left) z = 0 Bohr, (upper right) z = 0.4648 Bohr, (lower left) z = 0.9297 Bohr and (lower right) z = 1.3946 Bohr.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIG. 11: Slice of electronic density for a line along the z direction whose projection on the graphene layer is situated at the point $\mathbf{r} = (2.3244, 0)$ Bohr.</figcaption></figure>

In absence of $K_{xc}$, the screening function and its inverse are written as

$$ \epsilon^{-1, RPA} = \mathcal{I} + v\chi^{RPA}, \qquad (\text{D1}) $$

$$ \epsilon^{RPA} = \mathcal{I} - v\chi^0; \qquad (\text{D2}) $$

We define the screened Coulomb interaction as

$$ w = \epsilon^{-1, RPA} v; \qquad (\text{D3}) $$

if we now take Eq. A47 and A48, using Eq. A40 we obtain

$$ \delta \bar{\rho}_{s,\alpha}^{\text{ind}}(\mathbf{q}) + \delta \rho_{s,\alpha}^{\text{ext}}(\mathbf{q}) = \sum_{\mathbf{G}'} \frac{w(\mathbf{q}, \mathbf{q} + \mathbf{G}')}{w(\mathbf{q})} \delta \rho_{s,\alpha}^{\text{ext}}(\mathbf{q} + \mathbf{G}'). \qquad (\text{D4}) $$

Now, one uses Eq. 22 (where it is only used that $V^{\text{ext}} = v\rho^{\text{ext}}$) and eventually it follows

$$ \delta V_{s,\alpha}^{\text{tot}}(\mathbf{q}) = w(\mathbf{q}) \delta \bar{\rho}_{s,\alpha}^{\text{tot}}(\mathbf{q}), \qquad (\text{D5}) $$

$$ \delta \rho_{s,\alpha}^{\text{tot}}(\mathbf{q}) = \epsilon^{-1}(\mathbf{q}) \delta \bar{\rho}_{s,\alpha}^{\text{tot}}(\mathbf{q}), \qquad (\text{D6}) $$

$$ Z_{s,\alpha}(\mathbf{q}) = \epsilon^{-1}(\mathbf{q}) \bar{Z}_{s,\alpha}(\mathbf{q}). \qquad (\text{D7}) $$

A first crucial observation is that, numerically, in the undoped case the value of the unscreened Born effective charges is the same when computed in the RPA or in the RPA+xc with an LDA kernel, while this is not true for the value of $\epsilon^{-1}(\mathbf{q})$ [95]; with Eq. A48 in mind, we use this finding to state that at small $\mathbf{q}$ the values of all $S^{-1}$ elements are fairly independent of $xc$ effects in the undoped case [96].

In presence of $K_{xc}$, the screening function and its inverse (which cannot be called ‘dielectric’ anymore) are written as

$$ \epsilon^{-1} = \mathcal{I} + (v + K_{xc})\chi, \qquad (\text{D8}) $$

$$ \epsilon = \mathcal{I} - (v + K_{xc})\chi^0; \qquad (\text{D9}) $$

properly speaking $K_{xc}$ should present a long-wavelength behaviour proportional to $1/q^2$, in the most common approximations $K_{xc}$ is independent of $\mathbf{q}$ (in LDA) [97] or
dependent on positive powers of **q** (in GGA) [98]. We still define the screened Coulomb interaction as

$$w = \epsilon^{-1}v, \qquad (\text{D10})$$

even though the matrix is now not Hermitian and the rules for inversion change. For a generic non Hermitian matrix $w^{-1}$, we write [68]:

$$w^{-1} = \begin{pmatrix} P & Q \\ R & S \end{pmatrix}, \qquad (\text{D11})$$

$$w = \begin{pmatrix} W & X \\ Y & Z \end{pmatrix}, \qquad (\text{D12})$$

where the inversion rules are now

$$W = (P - QS^{-1}R)^{-1}, \qquad (\text{D13})$$

$$X = -WQS^{-1}, \qquad (\text{D14})$$

$$Y = -S^{-1}RW, \qquad (\text{D15})$$

$$Z = S^{-1} + YW^{-1}X, \qquad (\text{D16})$$

where we notice that Eq. D14 is equal to Eq. A40. The wing of $w^{-1}$ is now (assuming local $K_{xc}$ with no local field dependence for the sake of simplicity of the argument)

$$Q(\mathbf{q}, \mathbf{q} + \mathbf{G}) = - \frac{v(\mathbf{q}) + K_{xc}(\mathbf{q})}{v(\mathbf{q})} \chi^0(\mathbf{q}, \mathbf{q} + \mathbf{G}). \qquad (\text{D17})$$

As already discussed before, in the undoped case the presence of $K_{xc}$ in the response affects mainly the terms of $w$ which depend on $W$. Also, $xc$ terms are not expected to cure the non analytical behaviours typical of the RPA case, so that we are still induced to split the $w$ matrix in a short range component involving only the body $S^{-1}$, which is mainly RPA-like, and a remainder which is non analytical. To isolate the short-range component, we are therefore induced to adopt the same prescription of the RPA case to nullify the total macroscopic potential, comprised now also of the $xc$ component; this explain the procedure described in Sec. III C. Also, within our interpolating procedure explained in Sec. III D, this splitting means that in the doped case we are approximating $S^{-1}$ as the one of the undoped case; the goodness of the approximation is therefore connected to the magnitude of $xc$ contributions to the local fields in presence of doping. With the above prescription, Eq. A48 can now be written as

$$\delta \bar{\rho}_{s,\alpha}^{\text{ind}}(\mathbf{q}) = - \sum_{\substack{\mathbf{G} \neq 0 \\ \mathbf{G}' \neq 0}} \frac{v(\mathbf{q})}{v(\mathbf{q}) + K_{xc}(\mathbf{q})} Q(\mathbf{q}, \mathbf{q} + \mathbf{G}) \times \\ S^{-1}(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}') \delta \rho_{s,\alpha}^{\text{ext}}(\mathbf{q} + \mathbf{G}'). \qquad (\text{D18})$$

We use Eq. A40 obtaining

$$\delta \bar{\rho}_{s,\alpha}^{\text{ind}}(\mathbf{q}) = \sum_{\mathbf{G}' \neq 0} \frac{w(\mathbf{q}, \mathbf{q} + \mathbf{G}')}{\epsilon^{-1}(\mathbf{q})[v(\mathbf{q}) + K_{xc}(\mathbf{q})]} \delta \rho_{s,\alpha}^{\text{ext}}(\mathbf{q} + \mathbf{G}') \qquad (\text{D19})$$

$$\delta \bar{\rho}_{s,\alpha}^{\text{tot}}(\mathbf{q}) = \delta \bar{\rho}_{s,\alpha}^{\text{ind}}(\mathbf{q}) + \delta \rho_{s,\alpha}^{\text{ext}}(\mathbf{q}) \qquad (\text{D20})$$

As evident, in principle we lose the correspondence between the above equation and Eq. 22, and therefore barred quantities are no more related to the unbarred ones simply via the screening function. In practice, for the undoped case or small dopings and for the ordinary approximations to $K_{xc}$, we find that the correspondence is respected to a satisfactory precision. Notice that our considerations imply the well-known consequence that in the undoped case the Born-Huang formulae works even when their ingredients are computed in a RPA+xc approach; in fact, $\epsilon^{-1}$ in the above equations now contains the effects due to the presence of $K_{xc}$ in the response, while as already mentioned the effective charge functions are fairly insensitive to $xc$ effects. Notice also that, formally, with the current prescription the definition of, e.g., the Born effective charge tensors of Eq. 42 has to be intended as evaluated at total macroscopic field equal to zero, and not only the electrostatic one.

In a more general treatment, we find it interesting to rewrite the responses as

$$\epsilon^{-1} = \mathcal{I} + v\tilde{\chi} \qquad \epsilon = \mathcal{I} - v\tilde{\chi}^0, \qquad (\text{D21})$$

$$\tilde{\chi} = v^{-1}[v + K_{xc}]\chi \qquad \tilde{\chi}^0 = v^{-1}[v + K_{xc}]\chi^0. \qquad (\text{D22})$$

In this case, we can interpret $\tilde{\chi}^0$ as the polarizability response that enters in the Maxwell's equations. We now define a new induced density, namely

$$\delta \bar{\rho}_{s,\alpha}^{\text{ind}}(\mathbf{q}) = \sum_{\substack{\mathbf{G} \neq 0 \\ \mathbf{G}' \neq 0}} \tilde{\chi}^0(\mathbf{q}, \mathbf{q} + \mathbf{G}) S^{-1}(\mathbf{q} + \mathbf{G}, \mathbf{q} + \mathbf{G}') \times \\ \delta \rho_{s,\alpha}^{\text{ext}}(\mathbf{q} + \mathbf{G}'); \qquad (\text{D23})$$

summing the above expression to the external density, we obtain

$$\delta \bar{\rho}_{s,\alpha}^{\text{ind}}(\mathbf{q}) + \delta \rho_{s,\alpha}^{\text{ext}}(\mathbf{q}) = \sum_{\mathbf{G}'} \frac{w(\mathbf{q}, \mathbf{q} + \mathbf{G}')}{w(\mathbf{q})} \delta \rho_{s,\alpha}^{\text{ext}}(\mathbf{q} + \mathbf{G}') \qquad (\text{D24})$$

and therefore we can identify it with Eq. 22; we therefore have

$$\delta V_{s,\alpha}^{\text{tot}}(\mathbf{q}) = w(\mathbf{q}) \delta \bar{\rho}_{s,\alpha}^{\text{tot}}(\mathbf{q}) := -i \frac{eq}{A} w(\mathbf{q}) \tilde{Z}_{s,\alpha}(\mathbf{q}). \qquad (\text{D25})$$

Interestingly, the expression of $g^L$ is the same as Eq. 5 with the replacement $Z \to \tilde{Z}$, while in the expression for $C^L$ we can approximate $\tilde{Z}$ at its undoped value and send $Z \to \tilde{Z}$.

We end this section noting that the definition of effective charges given in presence of $xc$ in the context of this work, i.e. in the context of interpolation of short and long range components of the dynamical matrix and of the EPI, is not physical for the case of doped semiconductors or metals. The physical value would be in fact obtained defining the effective charges from the total density change obtained at null macroscopic electrostatic potential, leaving the $xc$ component of the total potential
unaltered; such density would be automatically analyti-

cal and would respect translational invariance. We leave
the investigation of such effective charges to future works.

[1] J. Ziman, *Electrons and Phonons: The Theory of Transport Phenomena in Solids*, International series of monographs on physics (OUP Oxford, 2001).

[2] G. Grimvall, *The electron-phonon interaction in metals* (North-Holland, Amsterdam, 1981).

[3] J. Schrieffer, *Theory Of Superconductivity*, Advanced Books Classics (Avalon Publishing, 1999).

[4] G. Mahan, *Many-Particle Physics*, Physics of Solids and Liquids (Springer US, 1990).

[5] F. Giustino, Electron-phonon interactions from first principles, Rev. Mod. Phys. **89**, 015003 (2017).

[6] F. Macheda and N. Bonini, Magnetotransport phenomena in *p*-doped diamond from first principles, Phys. Rev. B **98**, 201201 (2018).

[7] G. Brunin, H. P. C. Miranda, M. Giantomassi, M. Royo, M. Stengel, M. J. Verstraete, X. Gonze, G.-M. Rignanese, and G. Hautier, Phonon-limited electron mobility in si, gaas, and gap with exact treatment of dynamical quadrupoles, Phys. Rev. B **102**, 094308 (2020).

[8] F. Macheda, S. Poncé, F. Giustino, and N. Bonini, Theory and computation of hall scattering factor in graphene, Nano Letters **20**, 8861 (2020).

[9] S. Poncé, D. Jena, and F. Giustino, Hole mobility of strained gan from first principles, Phys. Rev. B **100**, 085204 (2019).

[10] S. Poncé, F. Macheda, E. R. Margine, N. Marzari, N. Bonini, and F. Giustino, First-principles predictions of hall and drift mobilities in semiconductors, Phys. Rev. Research **3**, 043022 (2021).

[11] N.-E. Lee, J.-J. Zhou, H.-Y. Chen, and M. Bernardi, Ab initio electron-two-phonon scattering in gaas from next-to-leading order perturbation theory, Nature Communications **11**, 1607 (2020).

[12] S. Poncé, W. Li, S. Reichardt, and F. Giustino, First-principles calculations of charge carrier mobility and conductivity in bulk semiconductors and two-dimensional materials, Reports on Progress in Physics **83**, 036501 (2020).

[13] J.-J. Zhou and M. Bernardi, Ab initio electron mobility and polar phonon scattering in gaas, Phys. Rev. B **94**, 201201 (2016).

[14] F. Macheda and N. Bonini, *Ab initio prediction of thermomagnetic and thermoelectric transport phenomena in 3d and 2d materials*, Ph.D. thesis (2021).

[15] J. Ma, A. S. Nissimagoudar, and W. Li, First-principles study of electron and hole mobilities of si and gaas, Phys. Rev. B **97**, 045201 (2018).

[16] C. Di Paola, F. Macheda, S. Laricchia, C. Weber, and N. Bonini, First-principles study of electronic transport and structural properties of cu₁₂sb₄s₁₃ in its high-temperature phase, Phys. Rev. Research **2**, 033055 (2020).

[17] R. Kim, V. Perebeinos, and P. Avouris, Relaxation of optically excited carriers in graphene, Phys. Rev. B **84**, 075449 (2011).

[18] M. Harb, R. Ernstorfer, T. Dartigalongue, C. T. Hebeisen, R. E. Jordan, and R. J. D. Miller, Carrier re-

laxation and lattice heating dynamics in silicon revealed by femtosecond electron diffraction, The Journal of Physical Chemistry B **110**, 25308 (2006).

[19] A. C. Betz, S. H. Jhang, E. Pallecchi, R. Ferreira, G. Fève, J.-M. Berroir, and B. Plaçais, Supercollision cooling in undoped graphene, Nature Physics **9**, 109 (2013).

[20] J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Theory of superconductivity, Phys. Rev. **108**, 1175 (1957).

[21] P. Hohenberg and W. Kohn, Inhomogeneous electron gas, Phys. Rev. **136**, B864 (1964).

[22] W. Kohn and L. J. Sham, Self-consistent equations including exchange and correlation effects, Phys. Rev. **140**, A1133 (1965).

[23] S. Baroni, S. de Gironcoli, A. Dal Corso, and P. Giannozzi, Phonons and related crystal properties from density-functional perturbation theory, Rev. Mod. Phys. **73**, 515 (2001).

[24] N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vanderbilt, Maximally localized wannier functions: Theory and applications, Rev. Mod. Phys. **84**, 1419 (2012).

[25] F. Giustino, M. L. Cohen, and S. G. Louie, Electron-phonon interaction using wannier functions, Phys. Rev. B **76**, 165108 (2007).

[26] M. Calandra, G. Profeta, and F. Mauri, Adiabatic and nonadiabatic phonon dispersion in a wannier function approach, Phys. Rev. B **82**, 165111 (2010).

[27] P. Giannozzi, O. Baseggio, P. Bonfa, D. Brunato, R. Car, I. Carnimeo, C. Cavazzoni, S. de Gironcoli, P. Delugas, F. Ferrari Ruffino, A. Ferretti, N. Marzari, I. Timrov, A. Urru, and S. Baroni, Quantum espresso toward the exascale, The Journal of Chemical Physics **152**, 154105 (2020), https://doi.org/10.1063/5.0005082.

[28] G. Pizzi, V. Vitale, R. Arita, S. Blügel, F. Freimuth, G. Géranton, M. Gibertini, D. Gresch, C. Johnson, T. Koretsune, J. Ibañez-Azpiroz, H. Lee, J.-M. Lihm, D. Marchand, A. Marrazzo, Y. Mokrousov, J. I. Mustafa, Y. Nohara, Y. Nomura, L. Paulatto, S. Poncé, T. Ponderer, J. Qiao, F. Thöle, S. S. Tsirkin, M. Wierzbowska, N. Marzari, D. Vanderbilt, I. Souza, A. A. Mostofi, and J. R. Yates, Wannier90 as a community code: new features and applications, Journal of Physics: Condensed Matter **32**, 165902 (2020).

[29] S. Poncé, E. Margine, C. Verdi, and F. Giustino, Epw: Electron-phonon coupling, transport and superconducting properties using maximally localized wannier functions, Computer Physics Communications **209**, 116 (2016).

[30] M. Born and K. Huang, Dynamical theory of crystal lattices, American Journal of Physics **23**, 474 (1955), https://doi.org/10.1119/1.1934059.

[31] H. Fröhlich, Electrons in lattice fields, Advances in Physics **3**, 325 (1954), https://doi.org/10.1080/00018735400101213.

[32] P. Vogl, Microscopic theory of electron-phonon interaction in insulators or semiconductors, Phys. Rev. B **13**,
694 (1976).

[33] C. Verdi and F. Giustino, Fröhlich electron-phonon ver- tex from first principles, Phys. Rev. Lett. 115, 176401 (2015).

[34] J. Sjakste, N. Vast, M. Calandra, and F. Mauri, Wan- nier interpolation of the electron-phonon matrix elements in polar semiconductors: Polar-optical coupling in gaas, Phys. Rev. B 92, 054307 (2015).

[35] G. Brunin, H. P. C. Miranda, M. Giantomassi, M. Royo, M. Stengel, M. J. Verstraete, X. Gonze, G.-M. Rignanese, and G. Hautier, Electron-phonon beyond fröhlich: Dy- namical quadrupoles in polar and covalent solids, Phys. Rev. Lett. 125, 136601 (2020).

[36] F. Macheda, P. Barone, and F. Mauri, Electron-phonon interaction and longitudinal-transverse phonon splitting in doped semiconductors, Phys. Rev. Lett. 129, 185902 (2022).

[37] Q. Ren, C. Fu, Q. Qiu, S. Dai, Z. Liu, T. Masuda, S. Asai, M. Hagihala, S. Lee, S. Torri, T. Kamiyama, L. He, X. Tong, C. Felser, D. J. Singh, T. Zhu, J. Yang, and J. Ma, Establishing the carrier scattering phase di- agram for zrnisn-based half-heusler thermoelectric mate- rials, Nature Communications 11, 3142 (2020).

[38] Q. H. Wang, K. Kalantar-Zadeh, A. Kis, J. N. Cole- man, and M. S. Strano, Electronics and optoelectronics of two-dimensional transition metal dichalcogenides, Na- ture Nanotechnology 7, 699 (2012).

[39] M. Chhowalla, D. Jena, and H. Zhang, Two-dimensional semiconductors for transistors, Nature Reviews Materials 1, 16052 (2016).

[40] M. Royo and M. Stengel, Exact long-range dielectric screening and interatomic force constants in quasi-two- dimensional crystals, Phys. Rev. X 11, 041027 (2021).

[41] T. Sohier, M. Calandra, and F. Mauri, Density func- tional perturbation theory for gated two-dimensional het- erostructures: Theoretical developments and application to flexural phonons in graphene, Phys. Rev. B 96, 075448 (2017).

[42] K. H. Michel and B. Verberck, Phonon dispersions and piezoelectricity in bulk and multilayers of hexagonal boron nitride, Phys. Rev. B 83, 115328 (2011).

[43] T. Sohier, M. Gibertini, M. Calandra, F. Mauri, and N. Marzari, Breakdown of optical phonons' splitting in two-dimensional materials, Nano Letters 17, 3758 (2017).

[44] T. Sohier, M. Calandra, and F. Mauri, Two-dimensional fröhlich interaction in transition-metal dichalcogenide monolayers: Theoretical modeling and first-principles calculations, Phys. Rev. B 94, 085415 (2016).

[45] J. Ma, D. Xu, R. Hu, and X. Luo, Examining two- dimensional fröhlich model and enhancing the elec- tron mobility of monolayer inse by dielectric engineer- ing, Journal of Applied Physics 128, 035107 (2020), https://doi.org/10.1063/5.0015102.

[46] K. Kaasbjerg, K. S. Thygesen, and K. W. Jacobsen, Phonon-limited mobility in *n*-type single-layer mos<sub>2</sub> from first principles, Phys. Rev. B 85, 115317 (2012).

[47] W. Li, S. Poncé, and F. Giustino, Dimensional crossover in the carrier mobility of two-dimensional semiconduc- tors: The case of inse, Nano Letters 19, 1774 (2019).

[48] T. Deng, G. Wu, W. Shi, Z. M. Wong, J.-S. Wang, and S.- W. Yang, Ab initio dipolar electron-phonon interactions in two-dimensional materials, Phys. Rev. B 103, 075410 (2021).

[49] M. M. Khatami, M. L. Van de Put, and W. G. Vanden-

berghe, First-principles study of electronic transport in germanane and hexagonal boron nitride, Phys. Rev. B 104, 235424 (2021).

[50] S. Poncé, M. Royo, M. Gibertini, N. Marzari, and M. Stengel, Drift and hall mobility of two-dimensional materials from first principles (2022).

[51] S. Poncé, M. Royo, M. Stengel, N. Marzari, and M. Gibertini, Long-range electrostatic contribution to the electron-phonon couplings and mobilities of two- dimensional materials (2022).

[52] K. Kaasbjerg, K. S. Thygesen, and A.-P. Jauho, Acoustic phonon limited mobility in two-dimensional semiconduc- tors: Deformation potential and piezoelectric scattering in monolayer mos<sub>2</sub> from first principles, Phys. Rev. B 87, 235312 (2013).

[53] R. Senga, K. Suenaga, P. Barone, S. Morishita, F. Mauri, and T. Pichler, Position and momentum mapping of vi- brations in graphene nanostructures, Nature 573, 247 (2019).

[54] T. Sohier, D. Campi, N. Marzari, and M. Gibertini, Mo- bility of two-dimensional materials from first principles in an accurate and automated framework, Phys. Rev. Ma- terials 2, 114010 (2018).

[55] T. Sohier, M. Gibertini, and N. Marzari, Profiling novel high-conductivity 2D semiconductors, 2D Materials, 1 (2020), arXiv:2007.16110.

[56] C. Zhang and Y. Liu, Phonon-limited transport of two- dimensional semiconductors: Quadrupole scattering and free carrier screening, Phys. Rev. B 106, 115423 (2022).

[57] T. Sohier, M. Gibertini, and M. J. Verstraete, Remote free-carrier screening to boost the mobility of fröhlich- limited two-dimensional semiconductors, Phys. Rev. Ma- terials 5, 024004 (2021).

[58] S. Y. Zhou, G.-H. Gweon, A. V. Fedorov, P. N. First, W. A. de Heer, D.-H. Lee, F. Guinea, A. H. Castro Neto, and A. Lanzara, Substrate-induced bandgap opening in epitaxial graphene, Nature Materials 6, 770 (2007).

[59] K. Novoselov, Mind the gap, Nature Materials 6, 720 (2007).

[60] M. S. Nevius, M. Conrad, F. Wang, A. Celis, M. N. Nair, A. Taleb-Ibrahimi, A. Tejeda, and E. H. Conrad, Semi- conducting graphene from highly ordered substrate in- teractions, Phys. Rev. Lett. 115, 136802 (2015).

[61] O. Bistoni, P. Barone, E. Cappelluti, L. Benfatto, and F. Mauri, Giant effective charges and piezoelectricity in gapped graphene, 2D Materials 6, 045015 (2019).

[62] D. Bohm and D. Pines, A collective description of elec- tron interactions. i. magnetic interactions, Phys. Rev. 82, 625 (1951).

[63] G. Giuliani, G. Vignale, A. Eringen, and C. U. Press, *Quantum Theory of the Electron Liquid*, Masters Series in Physics and Astronomy (Cambridge University Press, 2005).

[64] T. Sohier, M. Calandra, and F. Mauri, Density-functional calculation of static screening in two-dimensional materi- als: The long-wavelength dielectric function of graphene, Phys. Rev. B 91, 165428 (2015).

[65] E. H. Hwang and S. Das Sarma, Dielectric function, screening, and plasmons in two-dimensional graphene, Phys. Rev. B 75, 205418 (2007).

[66] Notice that the screened Coulomb potential depends on two reciprocal space index, as typical for the case of non- homogeneous electron gas; the homogeneous part of the interaction is indeed represented by the diagonal compo-
nents.

[67] M. Stengel, Flexoelectricity from density-functional perturbation theory, Phys. Rev. B **88**, 174106 (2013).

[68] R. M. Pick, M. H. Cohen, and R. M. Martin, Microscopic theory of force constants in the adiabatic approximation, Phys. Rev. B **1**, 910 (1970).

[69] In the case of more layers, one can model the response of the material as constant along the out-of-plane direction within the slab of the material, or one can introduce more refined models where each layer is treated separately, and remote couplings are considered.

[70] M. Born and K. Huang, Dynamical theory of crystal lattices (1954).

[71] P. Giannozzi, S. de Gironcoli, P. Pavone, and S. Baroni, Ab initio calculation of phonon dispersions in semiconductors, Phys. Rev. B **43**, 7231 (1991).

[72] P. Barone and F. Mauri, To be published, (2022).

[73] R. Resta and A. Baldereschi, Dielectric matrices and local fields in polar semiconductors, Phys. Rev. B **23**, 6615 (1981).

[74] S. Baroni and R. Resta, Ab initio calculation of the macroscopic dielectric constant in silicon, Phys. Rev. B **33**, 7017 (1986).

[75] Lannoo, M., Local field effects on the static dielectric constant of crystalline semiconductors or insulators, J. Phys. France **38**, 473 (1977).

[76] X. Gonze, J.-C. Charlier, D. Allan, and M. Teter, Interatomic force constants from first principles: The case of α-quartz, Phys. Rev. B **50**, 13035 (1994).

[77] At a given (small) doping, this is certainly valid in the usual case that the initial grid is coarse enough that the xLRCs of the undoped and doped setups coincide on that grid. However, if the initial coarse grid is fine enough to sample points within the metallic region of the response, then the different forms of the xLRCs invalidates the identity of the SRCs between the undoped and the doped setups on the coarse grid. This is when we make an approximation.

[78] J. Hong and D. Vanderbilt, First-principles theory of frozen-ion flexoelectricity, Phys. Rev. B **84**, 180101 (2011).

[79] J. P. Perdew, K. Burke, and M. Ernzerhof, Generalized gradient approximation made simple, Phys. Rev. Lett. **77**, 3865 (1996).

[80] L. Binci, P. Barone, and F. Mauri, First-principles theory of infrared vibrational spectroscopy of metals and semimetals: Application to graphite, Phys. Rev. B **103**, 134304 (2021).

[81] A. Damle, L. Lin, and L. Ying, Compressed representation of kohn-sham orbitals via selected columns of the density matrix, Journal of Chemical Theory and Computation **11**, 1463 (2015).

[82] V. Vitale, G. Pizzi, A. Marrazzo, J. R. Yates, N. Marzari, and A. A. Mostofi, Automated high-throughput wannierisation, npj Computational Materials **6**, 66 (2020).

[83] G. Pizzi, V. Vitale, R. Arita, S. Blügel, F. Freimuth, G. Géranton, M. Gibertini, D. Gresch, C. Johnson, T. Koretsune, J. Ibañez-Azpiroz, H. Lee, J.-M. Lihm, D. Marchand, A. Marrazzo, Y. Mokrousov, J. I. Mustafa, Y. Nohara, Y. Nomura, L. Paulatto, S. Poncé, T. Ponweiser, J. Qiao, F. Thöle, S. S. Tsirkin, M. Wierzbowska, N. Marzari, D. Vanderbilt, I. Souza, A. A. Mostofi, and J. R. Yates, Wannier90 as a community code: new features and applications, Journal of Physics: Condensed Matter **32**, 165902 (2020).

[84] M. Methfessel and A. T. Paxton, High-precision sampling for brillouin-zone integration in metals, Phys. Rev. B **40**, 3616 (1989).

[85] M. N. Gastiasoro, J. Ruhman, and R. M. Fernandes, Superconductivity in dilute sr dio3: A review, Annals of Physics **417**, 168107 (2020), eliashberg theory at 60: Strong-coupling superconductivity and beyond.

[86] A. Kumar, V. I. Yudson, and D. L. Maslov, Quasiparticle and nonquasiparticle transport in doped quantum paraelectrics, Phys. Rev. Lett. **126**, 076601 (2021).

[87] A. Kumar, P. Chandra, and P. A. Volkov, Spin-phonon resonances in nearly polar metals with spin-orbit coupling, Phys. Rev. B **105**, 125142 (2022).

[88] J. Ruhman and P. A. Lee, Superconductivity at very low density: The case of strontium titanate, Phys. Rev. B **94**, 224515 (2016).

[89] C. E. Dreyer, S. Coh, and M. Stengel, Nonadiabatic born effective charges in metals and the drude weight, Phys. Rev. Lett. **128**, 095901 (2022).

[90] S. K. Sinha, Theory of inelastic x-ray scattering from condensed matter, Journal of Physics: Condensed Matter **13**, 7511 (2001).

[91] X. Gonze, First-principles responses of solids to atomic displacements and homogeneous electric fields: Implementation of a conjugate-gradient algorithm, Phys. Rev. B **55**, 10337 (1997).

[92] X. Gonze and C. Lee, Dynamical matrices, born effective charges, dielectric permittivity tensors, and interatomic force constants from density-functional perturbation theory, Phys. Rev. B **55**, 10355 (1997).

[93] In other notations in the context of DFT or TDDFT $\chi^0$ is referred to as $\chi^{KS}$ or $\chi^{ir}$.

[94] H. Ehrenreich, Screening effects in polar semiconductors, Journal of Physics and Chemistry of Solids **8**, 130 (1959).

[95] In two dimensions the difference between the dielectric screenin in RPA and RPA+$xc$ for undoped semiconductors is anyway small. The same conclusions are reached though also for the three dimensional case of Ref. [36], where the difference is more marked.

[96] The effects of the exchange-correlation terms in $\chi^0$ are present through the periodic part of the Bloch function and the electronic energies; nonetheless, such effects are considered also in the RPA approximation.

[97] M. S. Hybertsen and S. G. Louie, Ab initio static dielectric matrices from the density-functional approach. i. formulation and application to semiconductors and insulators, Phys. Rev. B **35**, 5585 (1987).

[98] P. Ghosez, X. Gonze, and R. W. Godby, Long-wavelength behavior of the exchange-correlation kernel in the kohn-sham theory of periodic systems, Phys. Rev. B **56**, 12811 (1997).