# Ab-initio Van der Waals electrodynamics: polaritons and electron scattering from plasmons and phonons in BN-capped graphene

Francesco Macheda, Francesco Mauri, Thibault Sohier

► To cite this version:

Francesco Macheda, Francesco Mauri, Thibault Sohier. Ab-initio Van der Waals electrodynamics: polaritons and electron scattering from plasmons and phonons in BN-capped graphene. Physical Review B, 2024, 110 (11), pp.115407. 10.1103/PhysRevB.110.115407. hal-04723775

HAL Id: hal-04723775

https://hal.science/hal-04723775v1

Submitted on 7 Oct 2024

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
# Ab-initio Van der Waals electrodynamics: polaritons and electron scattering from plasmons and phonons in BN-capped graphene

Francesco Macheda,¹ Francesco Mauri,¹,² and Thibault Sohier³

¹Dipartimento di Fisica, Università di Roma La Sapienza, Piazzale Aldo Moro 5, I-00185 Roma, Italy

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

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

Plasmons and polar phonons are elementary electrodynamic excitations of matter. In 2d and at long wavelengths, they couple to light and act as the system polaritons. They also dictate the scattering of charged carriers. Van der Waals heterostructures offer the opportunity to couple excitations from different layers via long-range Coulomb interactions, modifying both their dispersion and their scattering of electrons. Even when the excitations do not couple, they are still influenced by the screening from all layers, leading to complex dynamical interactions between electrons, plasmons and polar phonons. We develop an efficient ab initio model to solve the dynamical electric response of Van der Waals heterostructures, accompanied by a formalism to extract relevant spectroscopic and transport quantities. Notably, we obtain scattering rates for electrons of the heterostructure coupling remotely with electrodynamic excitations. We apply those developments to BN-capped graphene, in which polar phonons from BN couple to plasmons in graphene. We study the nature of the coupled excitations, their dispersion and their coupling to graphene's electrons. Regimes driven by either phonons or plasmons are identified, as well as a truly hybrid regime corresponding to the plasmon-phonon-polariton at long wavelengths. Those are studied as a function of the graphene's Fermi level and the number of BN layers. In contrast with standard descriptions in terms of surface-optical phonons, we find that the electron-phonon interaction stems from several different modes. Moreover, the dynamical screening of the coupling between BN's LO phonons and graphene's electrons crosses over from inefficient to metal-like depending on the relative value of the phonons' frequency and the energetic onset of interband transitions. While the coupling is significant in general, the associated scattering of graphene's carriers is found to be negligible in the context of electronic transport.

## I. INTRODUCTION

Van der Waals heterostructures¹,² (vdWH) have emerged as some of the most promising ways to explore and exploit the properties of materials at the nanoscale³. Stacking different layers of 2d materials offers the opportunity to tailor novel properties. Unsurprisingly, this comes with some challenges on the modelling side. The electrodynamic response of a material, i.e. the response to a momentum and frequency-dependent electric potential $V(q, \omega)$, dictates a wide range of physical properties. This is evident from the fundamental importance of the inverse screening function⁴ in describing optical and electronic properties. The electrodynamic elementary excitations of the system are intrinsic, self-sustained collective modes driven by long-range electromagnetic interactions⁵. This work will focus on plasmons and polar phonons, that is collective oscillations of electronic or atomic plasma. Other excitations such as excitons or magnons are out of scope.

In 2d and at long wavelengths, longitudinal plasmons and phonons couple to light⁶⁻⁸. The resulting light-matter quasiparticles, polaritons, are usually discussed at momenta close to the light cone ($q_0 \gtrsim \omega/c$, where c is the speed of light). However, highly-confined plasmon- and phonon-polaritons⁹⁻¹⁶ in 2d materials have been observed at momenta 10 to 100 times larger. Those can be studied without accounting for relativistic retardation effects, as will be done here. Electrodynamic collective modes in

general are also a major source of scattering, via remote interactions between layers¹⁷,¹⁸. A good understanding and predictive models for the emerging interlayer couplings between collective modes thus lead to the exciting prospect of engineering light-matter interactions⁶,⁷ and electron scattering¹⁹ in future photonic and electronic devices based on vdWH.

To obtain polariton dispersions, one can resort to empirical and analytical approaches such as continuum electrodynamics²⁰⁻²², ultimately relying on parametrized models for the dielectric functions. As a consequence, electron scattering involving coupled plasmons and phonons systems²³⁻²⁷ is not treated from a fully microscopic point of view.

In order to be predictive in complex systems, a microscopic ab-initio approach that directly simulates the electrodynamic response of vdWH is highly desirable. Unfortunately, direct ab initio calculations of vdWH including plasmonic and phononic excitations are costly. Indeed, lattice mismatches imply large supercells which in addition to the number of layers increases considerably the number of atoms in the simulation cell. Another strategy is to coarse-grain the layer properties and focus on interlayer interactions. The single layer quantities are easily obtained from first principles, and used in a model to build the properties of the multilayer. In the current framework, this amounts to extracting single layer dynamical responses from ab-initio, and model the multilayer response.
We build on a previous work based on this strategy in the static limit¹⁹ and add dynamical electronic and atomic contributions. With respect to similar efforts in the literature²⁸⁻³¹, the Van der Waals electrodynamics (VED) model presented here offers both technical and conceptual improvements. Technical improvements include the full momentum-dependent phononic response, rather than resorting to the long-wavelength limit of Born effective charge models²⁹, which notably misses the correct relative weight of electronic and phononic polarizabilities. This is key to match our results to density functional perturbation theory (DFPT) in the appropriate limits, in terms of both dispersion and intensity of the collective modes. Conceptual progress stems from a simple yet powerful setup and exploitation of the formalism. We first build a matrix describing the response of given layer when a potential perturbation is applied to another layer. The imaginary part of this matrix yields the elementary excitations of the vdWH. From those, the response to an arbitrary external potential or charge density perturbation can be reconstructed and then projected on a given probe. When perturbation and probe are uniform over the layers, the model yields the basic spectroscopic properties of the slab, such as the Electronic Energy Loss Spectroscopy (EELS) cross Section in transmission. More complex quantities describing surface-probing experiments can also be extracted. Importantly, by perturbing and probing with a charge density localized on a given layer, we extract scattering rates for electrons coupling with plasmons, phonons, and hybrid excitations²³,²⁶,²⁷,³². Such ab-initio calculation of interlayer electron scattering mechanism, including dynamical screening from electrons and phonons, opens exciting perspectives for the study of vdWH. Indeed, they can be used for transport calculations²⁴,²⁷,³³⁻³⁹, Raman scattering⁴⁰⁻⁴², excited carriers relaxation¹⁷,⁴³⁻⁴⁷, and superconductivity⁴⁸⁻⁵⁰.

We focus on graphene encapsulated with multilayer hexagonal boron nitride, a prototypical system both for polaritonics and electronic transport. It hosts both charge plasmons from graphene and polar phonons from BN. Their coexistence in the same energy range leads to their mixing and anti-crossing¹²,²⁶,⁵¹. We present an extensive exploitation of the VED framework, discussing spectroscopic spectral functions and electron scattering including contributions from plasmons, phonons, and hybrid excitations. We analyse the contributions of different layers to gain insight on the nature of the excitations (plasmon, phonon or hybrid), and propose a method to separate the corresponding scattering strengths. We study the impact of two main parameters for BN-capped graphene, i.e. the number of BN layers and the Fermi level (electron doping) of graphene, that can be tuned experimentally⁵²,⁵³ to change the properties of the plasmons, phonons, and their scattering of electrons. We focus notably on the nature of the polaritonic states, and the dynamical screening of the remote coupling between graphene's electrons and BN's LO phonons. Gaining

access to the microscopic origin of the interactions, we progress over state-of-the-art modelling and show that in realistic systems the remote coupling is mediated by several different LO phonon branches, and assess its impact on the graphene's resistivity value at room temperature. The paper is structured as follows. Section II describes the theoretical and computational frameworks. Section III describes the collective modes of the BN/graphene systems. Section IIIA considers graphene's plasmons and BN's phonons before coupling, validating the method against known results and DFPT calculations. Section IIIB explores the nature of the excitations in the coupled system, as well as their dependence on the number of BN layers and graphene's Fermi level. Finally, Section IV studies the coupling of graphene's electrons with those collective modes, and compares the associated scattering rates to those associated with graphene's intrinsic optical modes.

## II. ELECTRODYNAMIC RESPONSE OF VAN DER WAALS HETEROSTRUCTURES

The density-density response of a system at momentum $q$ and frequency $\omega$, $\chi(q, \omega)$, enters the determination of many physical quantities, from spectroscopic responses to intrinsic scattering mechanisms. This Section first shows how the Van der Waals ElectroDynamics (VED) model builds the response of a layered heterostructure from single layer responses. Then, it shows how it can be used to obtain the collective modes dispersions, their EELS and surface-probes responses, and their coupling to electrons to deduce scattering rates.

### A. Van der Waals ElectroDynamics (VED) Model

The VED is the dynamical generalisation of the model developed in Ref.¹⁹. The general idea, sketched in Fig. 1, is to build the electrodynamic response of a layered Van der Waals heterostructure (vdWH) from the response of single layers. The following methodology is general, but we will specialize to BN-capped graphene as a prototypical system. The response of each different layer is simulated and parametrized from first principles. The static, clamped-ions electronic response is directly extracted from density functional perturbation theory (DFPT) calculations, therefore including localfields and exchange-correlation effects. The dynamical contribution from free carriers is modelled within the Random Phase Approximation (RPA), equivalent to a time-dependent Hartree theory without local-fields. The dynamical atomic response is evaluated from DFPT ingredients. The ab-initio simulation of the whole vdWH, very costly even in the static Born-Oppenheimer approximation of DFPT, is avoided. Direct ab initio simulations of the interplay of dynamical contributions from plasmons and phonons would further require time-dependent
density-functional theory, extremely costly in large supercells.

Here, the dynamical response of the full vdWH is obtained by solving a self-consistent set of equations that determine the potentials and charges on each layer. The computational burden of our method is limited to the single layer static DFPT response calculations needed for each different kind of layer, plus a small overhead to solve the coupled equations of the VED model⁵⁴. The single layer responses are easily parametrized and stored, such that for all vdWH combining layers in the database, only the VED coupled equations are left to solve. The main approximation of this methodology is to neglect the wavefunction overlap between layers. This essentially comes down to an interlayer coarse-graining approximation where different layers see each other only through average macroscopic potentials. This approximation holds remarkably well in vdWH, where the interlayer bonds are weak by definition.

<figure><img src="image_2.png" /><figcaption>FIG. 1: Sketch of the Van der Waals electrodynamics (VED) model. Single layer monopole ($f^0$) and dipole ($f^1$) density responses to even ($\phi^0$) and odd ($\phi^1$) potentials are extracted from DFPT. This is then used to build the multilayer response to arbitrary potential or charge perturbation.</figcaption></figure>

## 1. General definitions

The general aim is to simulate the macroscopic charge density response of the vdWH to a macroscopic external potential. When working with screening/dielectric properties in 2d, it is preferable to Fourier transform only the in-plane space variables ($x, y \to q$), while keeping the full out-of-plane $z$ dependence of the response. The external potential perturbation $V_{\text{ext}}(q, \omega, z)$ is thus periodic in-plane and monochromatic at momentum $q$ and frequency $\omega$. The total charge density response $\rho(q, \omega, z)$ is determined by both the rearrangement of the electronic cloud and the displacement of the point-like ions. For both electrons and atoms, the responses are assumed to be isotropic in the plane, depending only on $q = |q|$. This is usually a good approximation for many 2d material (e.g. hexagonal materials like BN and graphene). Within linear response theory, the charge density is obtained via

the density-density response function $\chi$ as

$$ \rho(q, \omega, z) = \int dz' \chi(q, \omega, z, z') V_{\text{ext}}(q, \omega, z'). \quad (1) $$

The integration in the out-of-plane variable is meant to cover all the space, since $z$ is free of periodic boundary conditions. The connection between the total and the external potential is then given by the inverse dielectric matrix

$$ \epsilon^{-1}(q, \omega, z, z') = \delta(z - z') + \frac{2\pi e^2}{q} \int dz'' e^{-q|z-z''|} \chi(q, \omega, z'', z'). \quad (2) $$

The possibility of having any kind of external potential has to be reduced in order to have a tractable problem, while at the same time maintaining the main physical features of the perturbations. Therefore, we find it useful to restrict the $z$ dependence of the external potential to a simplified functional form. As formalized in the following Section, we follow the approach of Ref. ¹⁹, which is an approximation of the exact method of Ref. ⁵⁵.

## 2. Dual basis set

To write the response problem in an easily solvable matrix form, the first step is to approximate the continuous out-of-plane variables of the response function of Eq. 1 with discrete indices over the layers. For each layer, $V_{\text{ext}}(z)$ is expanded over $N_b$ elements composing a subset of a complete basis set. The induced density can itself be expressed as a linear combination of (different) $N_b$ basis vectors. Eq. 1 then transforms into a matrix equation with densities and potentials expanded over a dual basis set (in the same spirit as Ref. ²⁸). We start by defining the mapping between the out-of-plane coordinate $z$ and the layer index. For $N_l$ layers, $k \in \{1, \dots, N_l\}$ and

$$ z \in [z_k - d/2, z_k + d/2] \to k. \quad (3) $$

We assume that each layer has a finite thickness $d$ in the out-of-plane direction around its central coordinate $z_k$, within which the density response is fully contained. $d$ depends on the kind of layer but for simplicity, here, we take it to be layer independent. Within the ‘boxes’ defined in Eq. 3, we define the restricted ($N_b = 2$) basis set for the external potential in the form

$$ \phi_k^i(z) = \begin{cases} (z - z_k)^i & z \in k \\ 0 & z \notin k \end{cases}, \quad (4) $$

with $i = 0, 1$. Analogously, we define the restricted basis set for the densities as

$$ f_k^i(z) = \begin{cases} \mathcal{F}_k^i(z - z_k) & z \in k \\ 0 & z \notin k \end{cases}, \quad (5) $$
which we refer to as profile functions. $\mathcal{F}_k^i(z)$ are the normalized out-of-plane profiles of the density response to the basis potentials of Eq. 4 when the $k$ layer is placed at $z_k = 0$, as depicted in Fig. 1. They separate the density response into monopolar (z-symmetric, $i = 0$) and dipolar (z-antisymmetric, $i = 1$) components. In principle, they might depend on $q$ and $\omega$. However, we have observed negligible dependence on $q$ and assume the same for $\omega$, as discussed in Section II A 6. The profile functions (Eq. (A2) of Ref. 19) are normalized as

$$ \int dz f_k^i(z) \phi_l^j(z) = \delta_{ij} \delta_{kl}. \qquad (6) $$

The admitted potentials in our problem can now be written as

$$ V(q, \omega, z) = \sum_{ik} V_k^i(q, \omega) \phi_k^i(z), \qquad (7) $$

$$ V_k^i(q, \omega) = \int dz f_k^i(z) V(q, \omega, z). \qquad (8) $$

Similarly, for the density

$$ \rho_k(q, \omega, z) = \sum_i \rho_k^i(q, \omega) f_k^i(z), \qquad (9) $$

$$ \rho_k^i(q, \omega) = \int dz \phi_k^i(z) \rho(q, \omega, z). \qquad (10) $$

Notice that in Eq. 1 and in the Random Phase Approximation (RPA), the induced density is the sum of two contributions due to the purely electronic polarization ('el') and the atomic (phonon) mediated one ('ph'). Since the observables that we study in this work are integrated quantities along $z$, it is convenient to approximate atomic polarization to have the same $z$-profile of the electronic one, with a small error as discussed later. In formulae

$$ \rho_k^i(q, \omega) = \rho_{k,el}^i(q, \omega) + \rho_{k,ph}^i(q, \omega). \qquad (11) $$

### 3. Single layer response function

We here consider the basic building block of our methodology, i.e. the response of each single layer to potential perturbations. As shown in Ref. 55, the out-of-plane dependence of the single layer density-density response function can be written in a separable form, i.e. with functions of only $z$ or $z'$. Accordingly, in this work we express the density-density response function for a single layer in the dual basis set as

$$ [\chi_{1L}]_k(q, \omega, z, z') = \sum_i Q_k^i(q, \omega) f_k^i(z) f_k^i(z'), \qquad (12) $$

$$ Q_k^i(q, \omega) = \int dz dz' \phi_k^i(z) [\chi_{1L}]_k(q, \omega, z, z') \phi_k^i(z'), \qquad (13) $$

where $Q_k^0$ and $Q_k^1$ are the amplitudes of the 'monopole' and 'dipole' response functions$^{56}$. Note that the off-diagonal matrix elements of $\chi_{1L}$ (i.e. $f_k^i(z) f_k^j(z')$ with

$i \neq j$ terms) are null since for the single layers considered in this work we have in-plane mirror symmetry. To simplify notations, we will now drop the $(q, \omega)$ dependence when not essential. The density response for a single layer embedded in an heterostructure reads

$$ \rho_k^i = Q_k^i V_{k,\text{eff}}^i, \qquad (14) $$

where $V_{k,\text{eff}}^i$ is the effective potential felt by the embedded single layer. It is the sum of the actual external potential and of the tails of the macroscopic Hartree potential induced by all the other layers. $V_{k,\text{eff}}^i$ reduces to just the external potential if the single layer is suspended in vacuum.

### 4. Interaction between layers

The process of obtaining the multilayer response from the single layer response functions above is described in Ref. 19. Only the main steps are summarized here. The objective is to reduce the problem to a finite and limited number of elements, i.e. the layers, interacting with each other via simple interlayer couplings. Those interlayer couplings contain all the information about the out-of-plane structure of the system. The starting point is the Poisson equation for the potential induced by the charge variations

$$ V_{\text{ind}}(q, \omega, z) = \frac{2\pi e^2}{q} \int dz' e^{-q|z-z'|} \rho(q, \omega, z'). \qquad (15) $$

In the dual basis set introduced in this work, the matrix elements of the Coulomb interaction read

$$ v_{kl}^{ij}(q) = \frac{2\pi e^2}{q} F_{kl}^{ij}(q), \qquad (16) $$

$$ F_{kl}^{ij}(q) = \int f_k^i(z) \int e^{-q|z-z'|} f_l^j(z') dz dz', $$

which essentially describe the potential generated by a unit monopole ($j = 0$) or dipole ($j = 1$) charge-density on layer $l$ and projected on layer $k$. Eq. 15 then becomes

$$ V_{k,\text{ind}}^i = \sum_{jl} v_{kl}^{ij} \rho_l^j. \qquad (17) $$

Notice that in the above equations the Coulomb potential is assumed to be frequency-independent. This assumption is valid if we restrict to wavevectors much larger than the light cone $q \gg \frac{\omega}{c}$, i.e. if we can disregard relativistic retardation effects. The collective modes discussed in this work, including highly-confined plasmon-phonon-polaritons, are far from the light cone. In general though, relativistic retardation effects are important when discussing phonon-polaritons in 3D materials$^{57}$, meaning that there will be a critical threshold for the thickness of any heterostructure after which the approximation of Eq. 16 fails to describe the physics of phonon-polaritons, i.e. their wavevectors get near to the light cone.
## 5. Multilayer response function

Given the coupling of Eq. 16, one can build a set of linear equations describing the collective response of the layers¹⁹. As anticipated, each individual layer $k$ is assumed to respond to an effective external potential $V_{\text{eff}}$ which is the sum of the applied external potential plus the sum of the potentials induced by all other layers through Eq. 17, i.e. to

$$V_{k,\text{eff}}^i = V_{k,\text{ext}}^i + \sum_{l \neq k, j=0,1} v_{kl}^{ij} \rho_l^j. \quad (18)$$

Using Eq. 14, we can write a self-consistent system of equations for $\rho_k^i$:

$$\rho_k^i = Q_k^i \left[ V_{k,\text{ext}}^i + \sum_{l \neq k, j=0,1} v_{kl}^{ij} \rho_l^j \right]. \quad (19)$$

Eq. 19 is composed of $2 \times N_l$ equations to solve for $\rho_k^i$ at each $q$ and $\omega$, where $i = 0, 1$ and $k = 0, \dots, N_l$.

Once the system is solved, one can determine the density-density response function of the full heterostructure. In fact, the matrix-form of the multi-layer density-density response function of Eq. 1 is

$$\chi_{kl}^{ij}(q, \omega) = \int dzdz' \phi_k^i(z) \chi(q, \omega, z, z') \phi_l^j(z'), \quad (20)$$

and its value can be deduced from

$$\rho_k^i = \sum_{jl} \chi_{kl}^{ij} V_{l,\text{ext}}^j. \quad (21)$$

Mind that the above equation is the matrix form of Eq. 1 written in the dual basis set introduced in Sec. II A 2, and it is different from the single layer density-density response function of Eq. 12, because it relates the induced density on each layer to the external potential applied to the full heterostructure, rather than to the effective potential felt by each layer.

## 6. Computation of single layer response

While the basis set in which potentials are decomposed are simple analytical functions, the profile functions used as a basis set for densities are material dependent and computed in DFPT, as explained in App. B. As anticipated, their dependence in $q$ is very mild. Including it for graphene and BN brings less than 1% impact on the observables treated in this work, with respect to taking the $q = 0$ value. We therefore neglect it.

The fundamental building block of the VED method is the single layer density response function to a potential perturbation. As anticipated, in RPA it is split in purely electronic (clamped-ion) and atomic-mediated contributions (total minus clamped-ion) as

$$Q^i(q, \omega) = Q_{\text{el}}^i(q, \omega) + Q_{\text{ph}}^i(q, \omega), \quad (22)$$

where we have dropped the layer index for simplicity, as we will do for the rest of this Section. In this work, we will only consider dynamical contributions from plasmons in metals (doped graphene) and longitudinal polar-optical LO phonons in semiconductors or insulators (BN). They contribute to the monopolar component of the electronic and atomic contributions to the response, respectively. The dipolar atomic contribution is neglected, while the electronic dipolar contribution is evaluated in its static limit. In other words, we assume:

$$Q^0(q, \omega) \approx \begin{cases} Q_{\text{el}}^0(q, \omega = 0) + Q_{\text{ph}}^0(q, \omega) & (\text{semiconductors}) \\ Q_{\text{el}}^0(q, \omega) & (\text{metals}) \end{cases}, \quad (23)$$

$$Q^1(q, \omega) \approx Q_{\text{el}}^1(q, \omega = 0). \quad (24)$$

In the general case, the dipolar atomic contribution stemming from the out-of-plane response of ZO phonons would appear in Eq. 24. Neglecting it is a safe approximation in the main energy range of interest for this work, i.e. around the BN's LO phonons energy. Indeed, $Q^1(q, \omega_{\text{LO}}) = Q_{\text{el}}^1(q, 0) \gg Q_{\text{ph}}^1(q, \omega_{\text{LO}})$, as numerically shown in App. B 1. For energies close to the BN's ZO phonon ($\omega_{\text{ZO}} \sim \omega_{\text{LO}}/2$), the approximation does not hold, since the response is dominated by ZO phonons. We choose to neglect ZO phonons in this work. However, one might wonder if the ZO phonon's response could lead to significant alterations of the plasmon response. This is not the case for the doping levels studied in this work, since the ZO coupling to graphene's plasmons is much weaker than for LO, as also noted in Ref. 29. This weaker coupling can be understood as follows. The electric field generated by graphene's plasmon extend to other layers as $e^{-q|z-z_{\text{Gr}}|}$. Since the plasmon dispersion is at relatively small momenta, this generated potential is fairly flat and its projections on the dipole profile functions of BN layers are small (compared to the monopole ones). The ZO mode being associated with the dipole part of the atomic response $Q_{\text{ph}}^1$, it will only be weakly coupled to the plasmon. Neglecting the atomic dipolar term in Eq. 24 is therefore justified also for the determination of the plasmon dispersion.

The static ($\omega = 0$) electronic contributions of Eqs. 23 and 24 are parametrized from DFPT calculations, as detailed in App. B and Ref. 19, while the dynamical ones are modelled as follows. For graphene, we need $Q_{\text{el}}^0(q, \omega)$. One can include the dynamical contribution from free carriers by using the RPA through the non-interacting (irreducible) polarizability, as

$$Q_{\text{el}}^{0}(q, \omega) = \frac{\chi_{\text{irr}}(q, \omega)}{1 - v_{\text{Gr},\text{Gr}}^{00} \chi_{\text{irr}}(q, \omega)}, \quad (25)$$

$$\chi_{\text{irr}}(q, \omega) = \frac{2}{(2\pi)^2} \sum_{nm} \int d^2\mathbf{k} \frac{n_{\varepsilon_{n\mathbf{k}}}^{\text{FD}} - n_{\varepsilon_{m\mathbf{k}+\mathbf{q}}}^{\text{FD}}}{\varepsilon_{n\mathbf{k}} - \varepsilon_{m\mathbf{k}+\mathbf{q}} + \hbar\omega + i\hbar\eta_{\text{pl}}} \quad (26)$$

$$\times |\langle u_{n\mathbf{k}} | u_{m\mathbf{k}+\mathbf{q}} \rangle|^2,$$
where $n_{\epsilon n_k}^{\text{FD}} = \frac{1}{e^{(\epsilon_{n_k} - \epsilon_F)/k_B T} + 1}$ is the Fermi-Dirac occupation function for state of momentum $\mathbf{k}$ in band $n$ with energy $\epsilon_{n_k}$ and Bloch periodic function $u_{n_k}$, $\epsilon_F$) being the Fermi energy, $k_B$ the Boltzmann constant and $T$ the temperature. There is a factor 2 for spin degeneracy. The plasmon linewidth $\eta_{pl}$ is an external parameter in the model. In the case of graphene considered in this work, we simply compute the above $\chi_{\text{irr}}$ semi-numerically in the Dirac cone model, with a Fermi velocity that includes GW corrections$^{58}$. The wavefunction overlap is $|\langle u_{n_k} | u_{m_k+q} \rangle|^2 = \frac{1+nm \cos\theta}{2}$ with $n, m = \pm 1$ for $\pi, \pi^*$ bands and $\theta$ the angle between $\mathbf{k}$ and $\mathbf{k}+\mathbf{q}$. Deviations from the Dirac cone model are assumed negligible for the frequencies and doping levels studied in the following. The effect of other band levels beside $\pi$ and $\pi^*$ is also neglected because their contribution to metallic screening is very weak$^{59}$.

For semiconducting materials with polar phonons, long-range Coulomb interactions induce macroscopic electric fields in the crystal. The microscopic responses to such fields add up in a contribution $Q_{\text{ph}}^0(q, \omega)$ that can be expressed in terms of macroscopic quantities only$^{55}$. We deduce $Q_{\text{ph}}^0(q, \omega)$ from the atomic contribution to the inverse dielectric function

$$[\epsilon_{1L}^{-1}]^i(q, \omega) = \int dzdz' \phi^i(z) \epsilon_{1L}^{-1}(q, \omega, z, z') \phi^i(z'). \quad (27)$$

For BN's single layer, the monopolar contribution may be written as (see App. A 2 b)

$$[\epsilon_{1L}^{-1}]^0(q, \omega) = [\epsilon_{1L,el}^{-1}]^0(q) \left( 1 + \frac{\mathbf{e}_{LOq} \mathcal{D}^L(q) \mathbf{e}_{LOq}^*}{(\hbar\omega + i\hbar\eta_{LO})^2 - \hbar^2\omega_{LOq}^2} \right), \quad (28)$$

where $\mathbf{e}_{LOq} \mathcal{D}^L \mathbf{e}_{LOq}^*$ is the long-range, in-plane polar contribution to the dynamical matrix projected on the eigenvector of mode $LO^{60}$, and

$$[\epsilon_{1L,el}^{-1}]^0(q) = 1 + v_{\text{BN,BN}}^{00}(q) Q_{\text{el}}^0(q) \quad (29)$$

is the static electronic dielectric function. The expression for $\mathbf{e}_{LOq} \mathcal{D}^L(q) \mathbf{e}_{LOq}^*$ where the profile of the atomic polarizations are approximated with the electronic ones is given in Eq. A26. We found that a simpler yet more accurate approximation (discussed in App. A 2 b) is to deduce $\mathbf{e}_{LOq} \mathcal{D}^L(q) \mathbf{e}_{LOq}^*$ from the difference of the squared frequencies of the LO phonon with respect to the TO phonon

$$\mathbf{e}_{LOq} \mathcal{D}^L(q) \mathbf{e}_{LOq}^* = \hbar^2 \omega_{LOq}^2 - \hbar^2 \omega_{TOq}^2. \quad (30)$$

where $\omega_{LOq}$ and $\omega_{TOq}$ are computed via direct DFPT calculations. In particular, for BN the single layer LO mode has a DFPT frequency of $\omega_{LOq} \sim 1500 \text{ cm}^{-1}$. The phonon linewidth $\eta_{LO}$ is an external parameter for our model, that we take as a constant ($\sim 1$ meV). We finally get

$$[\epsilon_{1L}^{-1}]^0(q, \omega) = 1 + v_{\text{BN,BN}}^{00}(q) [Q_{\text{el}}^0(q) + Q_{\text{ph}}^0(q, \omega)], \quad (31)$$

with

$$Q_{\text{ph}}^{0}(q, \omega) = \frac{\hbar^2 \omega_{\text{LOq}}^2 - \hbar^2 \omega_{\text{TOq}}^2}{(\hbar\omega + i\hbar\eta_{\text{LO}})^2 - \hbar^2 \omega_{\text{LOq}}^2} \times \left[ \epsilon_{1\text{L,el}}^{-1} \right]^{0}(q) \frac{1}{v_{\text{BN,BN}}^{00}(q)}. \quad (32)$$

The above expression contains all in-plane multipolar orders of the charge density expansions, going beyond the leading order from Born effective charges at $q \to 0$, an improvement with respect to Ref.$^{29}$. Notice that Eq. 32, as well as Eq. 25, contain the Coulomb kernel in the profile basis function, instead of its $q \to 0$ form $2\pi e^2/q$ used in Ref.$^{29}$. While asymptotically equivalent at $q \to 0$, this has significant quantitative impacts on the atomic response, phonon dispersions and electron-phonon couplings as $q$ increases, as shown in Fig. 16.

As a final comment, we remind that the disregard of dynamical effects in Eq. 24 means that we cannot study the out-of-plane plasmons of the heterostructures, which we leave for future studies.

## B. Measurable quantities

### 1. Density-density response

The density-density response function of the heterostructure $\chi$ defined in the previous Section can be used to evaluate several fundamental response functions. For example, in the transmission setup of Electronic Energy Loss Spectroscopy (EELS), one is interested in the scattering of an electron-beam with all the possible excitations of the system. An electron in the beam of energy $E_i = \hbar^2 q_i^2 / 2m_e$ (with a typical order of magnitude of $\sim 10$ keV) is scattered by the material, that absorbs a quantity of energy $\omega$ through the excitation of internal degrees of freedom. The Stokes EELS cross-section is then proportional to$^{61,62}$

$$\frac{d^2\sigma}{d\Omega d\omega}(q, \omega) \propto \frac{-1}{(q^2 + q_z^2)^2} \times [1 + n_{\omega}^{\text{BE}}] \text{Im} \left[ \int dz dz' e^{iq_z(z-z')} \chi(q, z, z', \omega) \right], \quad (33)$$

where $q_z = q_z(q, \omega) = q_i - \sqrt{q_i^2 - q^2 - 2m_e\omega/\hbar}$, $m_e$ being the electron mass, and $n_{\omega}^{\text{BE}} = \frac{1}{e^{\omega/k_B T} - 1}$ is the Bose-Einstein statistical distribution. The $q$-dependent prefactors of Eq. 33 are not important to the aims of this work, and we will drop them in the following. We will also drop $n_{\omega}^{\text{BE}}$, since it can be easily re-inserted a posteriori in the evaluation of the cross-section. Then, in the approximation where $q_z$ is negligible—as in typical EELS experiments on heterostructures that are not too
<figure><img src="image_3.png" /><figcaption>FIG. 2: Sketch of the computation of the macroscopic $\chi_M$ and of the coupling between electrons and collective modes in the VED framework. On the top, only the monopole part of the density-density response $\chi_M$ is used, since the projection of the dipole part on the probe vanishes. On the bottom, the induced potential response is plotted, which drives the coupling of electrons. The plotted curves are computed with the VED method in 2BN/Gr/2BN for relatively large momentum at $\omega \sim \omega_{LO}$.</figcaption></figure>

thick—the cross section is proportional to

$$ \frac{d^2\sigma}{d\Omega d\omega}(q, \omega) \propto -\operatorname{Im}[\chi_M(q, \omega)], \quad (34) $$

$$ \chi_M(q, \omega) = \int dz dz' \chi(q, \omega, z, z'). \quad (35) $$

'M' stands for 'macroscopic', since $\chi_M$ is the in-plane average of the full heterostructure's density-density response integrated along the out-of-plane coordinates. Eq. 35 is rewritten, in our formalism, as

$$ \chi_M(q, \omega) = \sum_{kl} \chi_{kl}^{00}(q, \omega), \quad (36) $$

i.e. as the response to an external potential perturbation of the form $V_{\text{ext}}(q, \omega, z) = \sum_k \phi_k^0(z)$. Consistently with the neglect of q-dependent prefactors, in this work we will always normalize $\chi$-related quantities at each q to

their maximum, unless otherwise stated. In doing so, the plots of $-\operatorname{Im}\chi_M$ or EELS scattering cross Sections become equivalent$^{63}$.

The peaks of $-\operatorname{Im}\chi_M$ in the $(q, \omega)$ plane determine the collective modes of the system that are symmetrical with respect to the out-of-plane centre of the heterostructure. The antisymmetric modes instead average to zero when integrated in the out-of-plane direction, and are therefore dark in an EELS experiment. The antisymmetric modes are though visible in the total spectral function of the density-density response, defined as

$$ \chi_{\text{Tr}}(q, \omega) = \sum_{ik} \chi_{kk}^{ii}(q, \omega), \quad (37) $$

Minus the imaginary part of Eq. 37 in fact contains all the collective excitations of the system, and it is well defined even where modes cross.

More complex response functions can be defined from the $\chi$ matrix, in particular for perturbations that are non uniform over the layers or containing finite projections on the dipole response. As an example, one could consider evanescent waves as a way to reveal plasmons of graphene on a substrate$^{9,10}$. We here propose a simple modellization for such kind of measurement. We assume that we have a tip that can both induce exponentially suppressed potentials in the heterostructure and measure averages of the induced potential on a given surface layer S. In our formalism, the form of the evanescent ('ev') external potential has components

$$ V_{k,\text{ev}}^{i}(q, \omega) = \int dz f_{k}^{i}(z) e^{-q|z-z_{\text{S}}|}, \quad (38) $$

where $z_S$ is the central coordinate of the surface layer. Using Eq. 17, the average of the potential on the surface layer, i.e. its monopolar component, is proportional to

$$ V_{S,\text{ind}}^{0} = \sum_{ik} v_{Sk}^{0i}(q) \rho_{k}^{i}(q, \omega) = \sum_{ijkl} v_{Sk}^{0i} \chi_{kl}^{ij}(q, \omega) V_{l,\text{ev}}^{j}(q, \omega). \quad (39) $$

From this induced potential we can define a 'local' density-density response via the Coulomb kernel

$$ \chi_{\text{loc}}(q, \omega) = V_{S,\text{ind}}^{0}(q, \omega) / v_{\text{SS}}^{00}(q). \quad (40) $$

More realistic and complicated descriptions are possible, but they go beyond the scope of this work. We end this Section reminding that from $\chi$ one can also derive the sheet optical conductivity, to define figures of merit such as the propagation quality factor, i.e. the distance (in number of wavelength) traveled before decay, and reflectivity$^{16,20}$.

## 2. Scattering rates

The electrons of the heterostructure experience scattering from collective modes. Understanding and quantifying those processes is essential to characterise the relaxation of the electron energy and momentum. In the single
layer case, the squared modulus of the electron-phonon and/or electron-plasmon couplings describing this scattering are obtained$^{23,64}$ via the product of the imaginary part of the inverse total dielectric function, containing electronic and atomic contributions, and the Coulomb kernel, as shown in App. A3. For the multilayer, if we only assume that the matrix elements of $(z-z_k)$ between Bloch states is small for every $k$, then the scattering rate for an electron in the Bloch state $nk$ reads (App. A3)

$$ \tau_{nk}^{-1} = \frac{2\pi}{\hbar} \frac{1}{N_q} \sum_{q} \int_{-\infty}^{\infty} d\omega A(\omega) \delta(\hbar\omega + \epsilon_{mk+q} - \epsilon_{nk}) \times \\ \sum_{kl} \langle u_{nk} | \phi_k^0 | u_{mk+q} \rangle \langle u_{mk+q} | \phi_l^0 | u_{nk} \rangle g_{kl}^2(q, \omega), \\ A(\omega) = \left[ n_{\omega}^{\text{BE}} + \frac{1}{2} + \text{sign}(\omega) \frac{1}{2} - \text{sign}(\omega) n_{\epsilon_{mk+q}}^{\text{FD}} \right], \quad (41) $$

where $N_q$ is the number of points in the $q$-grid, while $g_{kl}^2$ is

$$ g_{kl}^2(q, \omega) = \text{sign}(\omega) \sum_{ijk'l'} \frac{-1}{A\pi} \text{Im} \left[ v_{kk'}^{0i}(q) \chi_{k'l'}^{ij}(q, \omega) v_{l'l}^{j0}(q) \right]. \quad (42) $$

$A$ is the area of the unit cell of the single layer, taken equal for all the layers in this work for simplicity. Notice that the sign($\omega$) term is included in the definition of $g^2$ so that it doesn't become negative with $\omega \to -\omega$, due to the property Eq. A2. Scattering times are therefore always positive definite. Classically, the coupling $g_{kl}^2$ is the imaginary part of the potential projected on layer $k$ when layer $l$ is perturbed with a normalized charge density perturbation, as sketched in the bottom panel of Fig. 2.

The interaction of electrons of the heterostructure with phonons, plasmons, and any hybrid collective modes is entirely contained in the response of the system. In practice, it is mediated via every kind of collective modes contained inside $\chi$. Only in certain limits, i.e. away from plasmon-phonon hybridization, can the above coupling be formally identified with electron-phonon or electron-plasmon couplings, see discussion in App. A3.

Note also that $g_{kl}^2$ is the interaction stripped of the Bloch functions overlaps, while the scattering rate expression includes the overlaps. For BN, Bloch overlaps can be approximated to unity or zero at small $q$. For graphene, this would be notably wrong. Indeed, the wavefunctions overlap strongly depends on the angle between the momenta of the two electrons. We will nonetheless still study the coupling $g_{kl}^2$ in the following, and use Eq. 41 with the overlap given under Eq. 26 while studying the impact of remote coupling on graphene's electron scattering rates.

### III. COLLECTIVE MODES

In this section, the VED formalism is used to obtain the electrodynamic excitations of graphene and BN, and

their interplay within the heterostructure. The more general term "collective mode" is also used, with the understanding that in the current context, they are the one driven by long-range electromagnetic interactions. The treatment of the coupling between electrons and collective modes is treated in Sec. IV.

#### A. Graphene's plasmons and BN's phonons before coupling

In this Section, the collective modes of multilayer BN and graphene are considered separately, before any coupling between the two, as shown in Fig. 3. We show that the present model reproduces the features of multilayer BN and doped graphene, at the cost of a simple DFPT calculation for a single layer of each material. Before coupling, those systems are fairly well-known, such that the VED method can be compared to literature and direct DFPT calculations. At the same time, it brings new insight on certain aspects like the relative intensity of the different modes in the EELS intensity.

##### 1. graphene's plasmons

Graphene is a very promising medium for plasmonics$^{65,66}$, due to the strength and tunability of the plasmons. Analytical models$^{67-69}$ of their dispersion are obtained by considering the polarizability of doped graphene—Eq. (17) of Ref.$^{67}$— and the RPA relation to obtain the plasmon dispersion. The dielectric environment may be taken into account by a phenomenological background dielectric constant $\kappa$. In the small momentum limit, the plasmon dispersion is found to be

$$ \hbar\omega_{\text{pl}}(q) \sim_{q \to 0} \sqrt{\frac{2e^2\epsilon_F}{\kappa}} \sqrt{q}. \quad (43) $$

Two important features of the plasmon dispersion are shown here: the $\sqrt{q}$ asymptotic behavior, and a dispersion that scales with the Fermi surface, since $\hbar\omega_{\text{pl}}/\epsilon_F \propto \sqrt{q/k_F}$, with a proportionality constant that depends only on the dielectric environment.

The left panel of Fig. 3 shows $-\text{Im}[\chi_M(q, \omega)]$ for isolated graphene with a Fermi level of $\epsilon_F = 0.2$ eV and at room temperature. The only chosen parameter is the plasmon linewidth, fixed at $\eta_{\text{pl}} = 5$ meV, consistent in order of magnitude with theoretical$^{70,71}$ and experimental$^{11,13}$ investigations. The maxima of $-\text{Im}[\chi_M(q, \omega)]$ agree well with the asymptotic expression Eq. 43, as represented by the blue continuous line for $\kappa = 1$, at least up to $q/k_F \sim 0.2$. The same results (in units of $\hbar\omega/\epsilon_F$ and $q/k_F$) are obtained for different Fermi levels. The main features of the dispersion are thus recovered. Other spectral properties are also consistent with the literature. Notably, the spectral weight scales with $\epsilon_F^{67}$.
<figure><img src="image_4.png" /></figure>

FIG. 3: (Left) graphene's plasmonic excitation compared to the analytical model of Eq. 43 with $\kappa = 1$ (blue curve). The colourmap represents $-\text{Im}\chi_M$, normalized by the global maximum of the intensity. The dashed red curve is the extrapolated maximum of $-\text{Im}\chi_M$. The black dashed lines separate the following zones. I: no electron-hole transition and undamped/dispersed plasmon, II: interband electron-hole continuum, III: intraband electron-hole continuum, IV: no electron-hole transition. (Middle) Polar LO phonon modes of 5-BN from the resonances of the full spectral function $-\text{Im}\chi_{\text{Tr}}$, normalized at each momentum. The relative signs of the layer components of the phonon polarizations (corresponding to the sign of $\text{Im}\rho_k^M$, defined in the text) are as follows, from highest to lowest energy mode: $(+++++)$, $(- - 0 + +)$, $(- + + + -)$, $(+ - 0 + -)$, $(+ - + - +)$. Red crosses are phonon frequencies computed in DFPT. (Right) $\text{Im}\chi_M$ of 5-BN, normalized at each momentum. The highest LO mode is by far the most intense. Increasing the signal by a factor 100 in the inset reveals that two other modes are active, although much weaker.

In the left panel of Fig. 3, different zones are delineated with dashed lines, corresponding to the boundaries of the intra- and interband particle-hole continua. In practice, energy and momentum conservation allow intra- and interband electronic transitions only in zones III and II, respectively. Note that those boundaries are only strict at zero temperature. Otherwise they are smeared by the Fermi-Dirac occupation function.

## 2. BN's polar-optical phonons

Multilayer BN is ubiquitous in 2d devices$^{72}$. It is often used as an encapsulator, to protect the active layer from the environment and thus get closer to ideal intrinsic properties. By itself, it is a very promising platform for various photonics applications$^{73}$. In particular, BN's phonon-polaritons$^{15,20,74,75}$ have attracted a lot of interest based on their ability to shape and control light in matter. In 2d, phonon-polaritons are simply the polar-optical phonons of the system$^8$, and near-field microscopy techniques developed to probe polaritons in 2d materials are one of the best ways to probe BN's phonon dispersions, e.g. versus number of layers$^{20,74}$. In addition to

a detailed microscopic understanding, the VED model provides valuable insights to interpret and predict the results of such experiments. In particular, it clarifies which modes are active and their relative intensities.

The middle panel of Fig. 3 shows the total spectral function $-\text{ImTr}\chi$ for 5-layers BN (5-BN). This is computed in the VED model, using inputs from single-layer BN DFPT calculations. The only input parameter is the BN's phonon linewidth $\eta_{\text{LO}} = 0.6$ meV, chosen fairly small here in order to obtain separate peaks. Only phonon excitations are present, and the red crosses indicate phonon frequencies computed directly in DFPT. We observe a fairly good agreement. The error likely comes from modelling approximations, as discussed in App. B; also, in DFPT calculation each BN single layer does not possess anymore mirror symmetry as instead it is enforced in our vdWH.

Many interesting characteristics of polar phonons in 2d materials are recovered. The highest branch displays the dimensionality signature of the LO-TO splitting$^{76-78}$, i.e. the non-analytic and vanishing splitting at $q \to 0$, marked by the linear increase of the dispersion at small but finite momentum. The total number of LO modes is equal to the number of layers. The highest mode cor-
responds to in-phase contributions from all the layers, while the lower energy modes corresponds to various out-of-phase combinations. The slopes of the latter vanish at $q \to 0$.

The right panel of Fig. 3 shows $-\text{Im}\chi_M$ for 5-layers BN (5-BN). As already mentioned, among all the LO modes, only the symmetric ones appear in $-\text{Im}\chi_M$. Indeed, the perturbation and probe are symmetric in $z$, implying that the anti-symmetric modes are inactive. The symmetry of the modes is obtained via the relative sign of the imaginary part of the density response $\rho_k^M = \sum_l \chi_{kl}^{00}$ to a uniform external perturbation at the different resonances, and it is indicated in the caption of Fig. 3. The relative amplitude of the layers' density response decreases from the center to the outside of the stack, consistent with the relative amplitude of phonon displacement obtained in DFPT. The VED framework provides solid grounds to explain or predict the relative intensity of the modes in, e.g., EELS experiments. In the case of uniform perturbation and probe, note that the intensity of the peaks decreases with energy, and the highest one dominates largely.

## B. Interplay of plasmons and phonons

Beyond the interest of graphene and BN taken individually, it is essential to investigate their interactions. Since BN is often used as an encapsulator to protect graphene and other 2d materials, it is important to understand exactly how it can interact and modify the intrinsic properties of the encapsulated layer. The coupling between the two can also be exploited to tailor the polaritons$^{12,51}$, thus opening new pathways for electromagnetic waves manipulation; or to better control the energy and momentum relaxation of graphene's electrons$^{17}$. For BN-capped graphene, the main tuning parameters are the Fermi level of graphene and the number of BN layers. We study the interplay of graphene's plasmons and BN's LO phonons by focusing on two aspects: their dispersion, discussed here, and their scattering of graphene's electrons, discussed in Sec. IV. On top of previous analytical works on plasmon-phonon interactions$^{23,26,27,71}$, we bring the quantitative and predictive insight of microscopic ab initio simulations.

Fig. 4 shows $-\text{Im}\chi_M$ for the prototypical system made of graphene with a Fermi level of $\epsilon_F = 0.2$ eV encapsulated by 10-BN on each side. The plasmon and phonon linewidths are set to $\eta_{\text{pl}} = 5$ meV, $\eta_{\text{LO}} = 1$ meV (for numerical efficiency and plotting purposes, within experimental order of magnitude). We clearly observe an anti-crossing that indicates the plasmon-LO phonon coupling. Let us already name $(q, \omega)$ regions of interest outside of the anti-crossing, as annotated in Fig. 4: (1) the plasmon region covers the plasmon dispersion away from the typical LO phonon frequency, $\omega \neq \omega_{\text{LO}}$; (2) the phonon region follows the phonon dispersion at large momenta, to the right of the plasmon dispersion $q/k_F > \hbar\omega_{\text{LO}}/\epsilon_F$; (3)

<figure><img src="image_5.png" /><figcaption>FIG. 4: (Top) Macroscopic $-\text{Im}\chi_M$ for 10-BN/Gr@0.2eV/10-BN, normalized at each momentum. Blue and red dashed lines are the uncoupled excitations: LO dispersion of 20 BN layers and plasmon of free-standing graphene, respectively.</figcaption></figure>

the plasmon-phonon hybrid region covers the small momenta $q \to 0$ around the LO phonon frequency $\omega \sim \omega_{\text{LO}}$.

In the plasmon region (1), the plasmon is screened by the presence of BN, pushing its dispersion towards the $\hbar\omega/\epsilon_F = q/k_F$ dashed line, as expected from external screening in Eq. 43. Along their dispersion, the LO phonons experience different screening regimes from graphene. In the phonon region (2), they experience metallic screening, suppressing the LO-TO splitting and pushing the dispersion downward. Indeed, in the presence of static metallic screening, the slope of the highest mode dispersion at small momenta would vanish. Here, however, in the polariton region (3) the free-carrier screening from graphene is inefficient and we recover the finite slope one would obtain in BN alone for $q \to 0$ and $\omega \sim \omega_{\text{LO}}$, and a stiffening of the slope at small but finite momenta.

Those results are developed in the following. First, the nature of the excitations is carefully studied, in order to disentangle, when possible, phonons and plasmons. Second, we consider the case of a surface probing setup to access other collective modes. Third, the effect of the number of BN layers is analysed looking at 1, 10, and 30 BN layers on each side. Finally, the effect of graphene's Fermi level is considered, with two other representative doping: $\epsilon_F = 0.1$ and 0.3 eV.

### 1. Nature of the excitations and layer contributions

Here, we look at the layers from which the response originates to unravel the nature of the excitations. This procedure provides valuable insights into the interplay of phonons and plasmons, showing which mode drives the
<figure><img src="image_6.png" /><figcaption>FIG. 5: 10BN/Gr@0.2eV/10BN. Contributions to $-Im\chi_M$ from graphene and BN, $-Im\chi_{Gr}$ and $-Im\chi_{BN}$ as described by Eq. 45, normalized at each $q$ point. The inset focuses on the hybrid plasmon-phonon excitation at small $q$ and around $\hbar\omega_{LO}$, with the intensity magnified by a factor of 3 for better visualization. Blue and red dashed lines are the uncoupled excitations: LO dispersion of 20 BN layers and isolated graphene's plasmon, respectively.</figcaption></figure>

excitation, and how electrons screen it.

In Fig. 5, a unit potential perturbation is applied to all layers and the response coming from either graphene or BN are separated as follows

$$ \chi_M = \chi_{BN} + \chi_{Gr}, \quad (44) $$

$$ \chi_{BN} = \sum_{k=BN,l} \chi_{kl}^{00}, \quad \chi_{Gr} = \sum_{k=Gr,l} \chi_{kl}^{00}. \quad (45) $$

Notice that the sum over BN runs over all the BN layers. Fig. 5 shows $-Im\chi_{BN}$ and $-Im\chi_{Gr}$. It clearly reveals 3 types of excitations, clarifying the nomenclature of the 3 regions of $(q, \omega)$ space defined in the previous section. First, in the phonon region (2), we observe a positive contribution from BN layers and a negative one from graphene. Note that the sum of the two is always positive, so that here $|-Im\chi_{BN}| > |-Im\chi_{Gr}|$. Therefore, the excitation is driven by BN's polar phonons, and graphene is simply responding with (static) free-carrier screening. Second, in the plasmon region (1), we observe instead a positive contribution from graphene, and a negative one from BN layers. Thus, the plasmon is driving the excitation, and BN's electrons respond by screening it. Third, in the plasmon-phonon hybrid region (3) we observe positive contributions from both graphene and BN. This is a true hybrid polariton, driven by both plasmonic and phononic responses. The insets zoom in on this plasmon-phonon hybrid. The dispersions of the modes before the coupling are also plotted, i.e. graphene's plasmon without BN and the highest LO phonon mode corresponding to the total number of BN layers in absence of graphene. The slope of the hybrid excitation dispersion (i.e. the

group velocity of the polariton) approaches that of the isolated 20-BN system when $q \to 0$. In that limit, the polariton gradually becomes pure LO phonon, and contextually $-Im\chi_{Gr}$ vanishes. At small but finite $q$, however, the slope remains fairly high, rather than flattening like that of the intrinsic BN's phonons. This is thanks to a gradual shift of the plasmon-phonon hybrid from phonon to plasmon, as $-Im\chi_{BN}$ decreases with increasing $q$ (eventually changing sign) while $-Im\chi_{Gr}$ increases. Eventually, following this branch, we fall back on the plasmon-driven excitation. We thus have a clear interpretation of all the excitations and the role played by the electrons, phonons, and plasmons.

## 2. Probing other modes

The macroscopic response function $\chi_M$ only probes symmetric modes because perturbation and probe are uniform over the layers. This is not the case of surface spectroscopic techniques, like scattering-type Scanning Near-field Optical Microscopy$^{20}$ (sSNOM) or 2d High-resolution electron energy loss spectroscopy$^{79}$ (2d - HREELS). A quantitatively accurate simulation of those experimental setups is out of scope here. However, it is informative to consider a setup where perturbation and probe are localised on the surface, and qualitatively assess the consequences on the activity and relative intensity of the collective modes. Using the model explained in Sec. II B 1, we expect the physical observables to be proportional to $\chi_{loc}(q, \omega)$, plotted in Fig. 6.
<figure><img src="image_7.png" /><figcaption>FIG. 6: $-\text{Im}\chi_{\text{loc}}$ normalized at each $q$ in 10-BN/Gr@0.2eV/10-BN, representing the results of a probe measurement at the surface, with the model explained in Sec. II B 1. The inset zooms in on the anti-crossing region, and shows the spectral function $-\text{Im}\chi_{\text{Tr}}$.</figcaption></figure>

The relative intensities of the modes is very different with respect to the plots of $-\text{Im}\chi_M$, with the out-of-phase LO modes more visible. Interestingly, the most intense phonon-like mode is not the highest LO mode, but the mode corresponding to the BN layers moving in-phase with respect to layers on the same side of the graphene layer, but out-of-phase with respect to the BN layers on the other side of graphene. This mode is unaffected by screening from graphene, since any electric field from phonons on one side of graphene would be equal and opposite to the field from the other side. Thus, the net potential from phonons felt by graphene vanishes, and it does not respond. This mode actually crosses the dispersion of the in-phase LO mode, becoming higher in energy when the latter is strongly screened by graphene. This raises the question of which mode is actually most seen in surface-probing experiments. Near field spectroscopy would probe the vicinity of the light cone, dominated by the plasmon-phonon hybrid at vanishing momentum. Techniques probing larger momenta like 2d-HREELS$^{79}$, however, might very well be dominated by this anti-symmetric mode.

### 3. Effects of number of BN layers

The number of BN layers may vary significantly in different realisation of similar devices. Whether it is a parameter one can control precisely to tailor the properties of the system, or a fluctuating number depending on the experimental realization of the sample, it is important to understand its impact on the collective modes of BN-

capped systems. Beyond the fairly obvious increase of the plasmon-phonon coupling, the VED framework enables the quantification and a detailed interpretation of the changes in the collective modes dispersions.

<figure><img src="image_8.png" /><figcaption>FIG. 7: $-\text{Im}\chi_M$ for 1-BN/Gr@0.2eV/1-BN (top) and 30BN/Gr@0.2/30BN (bottom) normalized at each $q$. The blue curve and dashed red line are the uncoupled modes (LO dispersion of total number of BN layers in blue, isolated graphene's plasmon in red)</figcaption></figure>

Fig. 7 shows $-\text{Im}\chi_M$ for graphene at $\epsilon_F = 0.2$ eV encapsulated by 1 or 30 layers of BN on each side, on the top and bottom panels respectively. This complements the dispersion of Fig. 4, with 10-BN on each side. As the number of BN layers increases, the features already noted in the previous Sections are enhanced. Due to overall stronger phonon responses, the plasmon-phonon coupling and the size of the anti-crossing increase. Also, the stronger electronic screening from BN pushes the plasmonic dispersion closer to the $\hbar\omega/\epsilon_F = q/k_F$ line. The slope of the polariton at $q \to 0$ and $\omega \sim \omega_{LO}$ increases,
following that of the intrinsic LO phonon, as the polar
contributions from all the layers add up$^{78}$.

We remind that, as discussed after Eq. 16, our frame-
work can be applied until the slope of the longitudinal
phonon is not so steep to enter in the light-cone region,
an event that happens for a certain limiting width of the
layer (way beyond what is considered here). Increasing
the width further causes relativistic retardation effects
to become important and we therefore enter the realm of
three-dimensional phonon-polaritonics. Notice that this
gives a rigorous criteria to distinguish between a slab and
a proper three-dimensional materials from the polariton-
ics point of view.

4. *Effects of graphene's Fermi level*

Graphene's Fermi level is often tuned electrostatically$^{80-82}$. Within a field-effect transistor geometry, free carriers are accumulated in graphene, typically up to the order of $10^{13}$ cm$^{-2}$. This corresponds to Fermi levels of a few hundreds of meV. Thus, in the following, three Fermi levels $\epsilon_F = 0.1, 0.2, 0.3$ eV are considered and the consequences on the collective modes of the BN-encapsulated graphene are studied. Since the effect of doping on the plasmon are relatively well-known, we focus on the rich and yet unresolved discussion around the dynamical screening of phonons and their coupling to electrons.

Fig. 8 shows $-\text{Im}\chi_M$ for graphene encapsulated by 10 BN, at Fermi levels $\epsilon_F = 0.1$ eV and $\epsilon_F = 0.3$ eV, in addition to the $\epsilon_F = 0.2$ eV case from Fig. 4. In all cases, some trends persist. The plasmon-phonon hybrid follows the intrinsic LO phonon dispersion in the $q \to 0$ limit, while the LO-TO splitting is reduced due to free-carrier screening from graphene in the intraband electron-hole continuum, that is for $q > \frac{\hbar\omega_{LO}}{\epsilon_F}k_F = \frac{\omega_{LO}}{v_F}$.
The behavior between those two regimes ($q \to 0$ and $q > \frac{\omega}{v_F}$), however, is quite sensitive to the Fermi level. In the case of $\epsilon_F = 0.2$ eV, where phonon energy and Fermi level are similar, the plasmon-phonon anti-crossing happens around the middle of the triangular region I. For $\epsilon_F = 0.3$ eV, this is shifted towards smaller momenta, and so is the region where we can observe the intrinsic LO phonon dispersion. Indeed, assuming the plasmon dispersion to be of the form Eq. 43, the momentum $q$ where the plasmon and phonon cross is proportional to $k_F$. For both $\epsilon_F = 0.2$ eV and $\epsilon_F = 0.3$ eV, we observe the intrinsic and screened LO dispersion on the left and right of the anti-crossing, respectively. The qualitative behavior starts to change significantly when $\hbar\omega_{LO} > \epsilon_F$, and the crossing of plasmon and phonon dispersions fall into the interband electron-hole continuum. There is no clear anti-crossing then. Instead, a collective mode emerges that roughly follows the intrinsic LO phonon dispersion. A clear free carrier screening only starts in the intraband continuum. Qualitatively, the difference between interband and intraband comes from the possibility in

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

<figure><img src="image_10.png" /><figcaption>FIG. 8: $-\text{Im}\chi_M$ for 10-BN/Gr@0.1eV/10-BN (top) and 10-BN/Gr@0.3eV/10-BN (bottom), normalized at each $q$, in the plotted $\omega$ window. See also Fig. 4 for the case of $\epsilon_F = 0.2$ eV. Dashed black lines indicate intra- and inter-band electron-holes continuum. The dashed red line corresponds to the plasmon dispersion in isolated graphene. The blue line indicates the phonon dispersion of the highest LO mode in 20-BN without graphene.</figcaption></figure>

the latter case to have Lindhard/Thomas-Fermi type of
screening$^{57,83}$, ultimately coming down to the non van-
ishing Bloch overlaps in the $q \to 0$ limit$^{60}$.

IV. ELECTRON SCATTERING

The VED model provides the coupling between elec-
trons and the various collective modes, to be included in
the calculations of scattering times. Above, the study of
quantities like $\chi_M$ or $\chi_{Tr}$ informed us about the disper-
sion and nature of the modes. In this Section we inves-
tigate how those modes may scatter the electrons of the heterostructure.

A. Validation in multilayer BN

<figure><img src="image_11.png" /><figcaption>FIG. 9: Total coupling, for a 5-BN heterostructure, compared to DFPT calculations performed as described in App. C. The coupling is entirely due to the interaction between electrons and phonons, statically screened by the BN's dielectric function.</figcaption></figure>

We first validate the method by comparing the VED results on multilayer BN to the DFPT ones, obtained as described in App. C. Indeed, since multilayer BN is a large gap semiconductor, all electronic contributions to the response are treated in the static limit within the VED model in that case. This corresponds to the adiabatic Born-Oppenheimer approximation of DFPT, and the results should thus differ only due to the approximations of the VED model (see App. B).

Fig. 9 shows the electron-phonon coupling with the LO modes in 5-BN. Practically, using Eq. 42 we compute

$$g^2(q, \omega) = \frac{1}{N_l^2} \sum_{kl} g_{kl}^2(q, \omega), \quad (46)$$

which physically is the square of the electron-phonon coupling since other collective modes are absent. We also define

$$\bar{g}^2(q) = \int_0^\infty g^2(q, \omega) d\omega. \quad (47)$$

In the following, the bar sign will be used consistently on all other $g^2$ quantity to signify integration over $\omega$ from zero to infinity. Notice that the sums run on all the BN layers. In fact, in the vdW approximation the Bloch electronic state of the full heterostructure layers are degenerate so that electrons are delocalized over the full structure. The agreement with DFPT is excellent, with the remaining error likely due to aforementioned modelling approximations.

Note that we do not force any phonon perturbation singularly in the VED method. The external perturbation itself excites all the possible phonons modes. Therefore, we get the sum of the couplings with all phonon modes, so we compare with the sum of the (squared) couplings over DFPT modes. Nevertheless, as seen in DFPT and expected in VED, the highest mode largely dominates the coupling. This result also validates the spectral weights of the peaks in the VED method as a function of momentum, which is a significant progress with respect to other ab-initio works$^{29}$.

B. Scattering from coupled plasmon and phonons

We now continue with the prototypical system consisting of graphene sandwiched by $n$ layers of BN on each side. Here we consider the coupling with graphene's electrons, i.e. we compute

$$g_{\text{Gr}}^2(q, \omega) = \sum_{ijk'l'} \frac{-1}{A\pi} \text{Im} \left[ v_{\text{Gr}k'}^{0i}(q) \chi_{k'l'}^{ij}(q, \omega) v_{l'\text{Gr}}^{j0}(q) \right]. \quad (48)$$

The above equation can be interpreted as using normalized charge densities as perturbation and probes of the system, as sketched in Fig. 2. The coupling of Eq. 48 contains contributions from phonons and plasmons. We plot it in the central column of Fig. 10 in the full $(q, \omega)$ plane, normalized at each $q$, for different $n$ (1,10 and 30 BN per side) at a doping level of $\epsilon_F = 0.2$eV. The left and right panels of Fig. 10 show $g_{\text{Gr}}^2(\omega)$ and its running integral $\int_0^\omega g_{\text{Gr}}^2(\omega')d\omega'$, $q = 1.5k_F$ and $q = 3.8k_F$. In both cases, the number peaks around the phonon energy increases with the number of BN layers. Clearly, graphene's electrons couple to several LO modes with different polarizations. This implies that considering a single surface-optical phonon for the remote coupling between electrons and phonons, as routinely done for bulk substrates like $SiO_2$$^{84}$, is not valid for BN encapsulation and vdWH in general. At $q = 3.8k_F$, up to $\hbar\omega/\epsilon_F \sim 1.75$, no electron-hole excitation is possible, and the coupling is totally determined by the electron-phonon interaction. For $q = 1.5k_F$, however, $g_{\text{Gr}}^2(\omega)$ contains significant contributions from the low-energy spectrum of the intraband electron-hole continuum for $\hbar\omega/\epsilon_F < 1$, from the different phononic branches of BN around $\hbar\omega/\epsilon_F \sim 1$, and finally from the dispersed plasmon dispersion in the interband region for $\hbar\omega/\epsilon_F > 1.5$.

While it contributes to electronic lifetimes detected e.g. via ARPES, the scattering of graphene's electrons coming from graphene's own electronic excitations (plasmon and electron-hole excitations) does not necessarily contribute to the momentum relaxation responsible for resistivity (depending on the magnitude and origin of plasmon dampening). Indeed, momentum is exchanged with essentially the same entity, i.e. graphene's electrons$^{26}$. Electron-phonon scattering, on the other hand, does. It
<figure><img src="image_12.png" /><figcaption>FIG. 10: (Central) $g_{Gr}^2(q, \omega)$, normalized at each q, for an heterostructure of n-BN/Gr@0.2eV/n-BN, with n=1,10,30. The blue and orange vertical lines represent the q values for which $g_{Gr}^2(\omega)$, and its running integral, are represented in the (left) and (right) panels. In the left panel, the contributions to $g_{Gr}^2(\omega)$ come from the intraband electron-hole continuum ($\hbar\omega/\epsilon_F < 1$), from the electron-phonon coupling $\hbar\omega/\epsilon_F \sim 1$, and from the coupling with the dispersed plasmon in the interband region $\hbar\omega/\epsilon_F > 1.5$. In the right panel, up to $\hbar\omega/\epsilon_F > 1.75$ only electron-phonon is present. Notice that the electron-phonon coupling gets contribution from many different modes, contrary to the picture of one single surface phonon responsible for all the interaction.</figcaption></figure>

is then desirable to separate the different couplings arising from the interaction with different collective modes.

1. Separation of the coupling with different modes

In Section III B 1 we separated the contributions from graphene and BN to gain insight on the nature of the excitations. Fig. 11 shows the results of a similar procedure applied to interaction between electrons and collective mode interactions. We separate the contributions as follows

$$g_{Gr}^2 = g_{Gr}^{2,\text{BN}} + g_{Gr}^{2,\text{Gr}}, \quad (49)$$

$$g_{Gr}^{2,\text{BN}}(q, \omega) = \frac{-1}{A\pi} \sum_{k=\text{BN}} \sum_{ijl} v_{\text{Gr},k}^{0,i} \text{Im} \left[ \chi_{kl}^{ij} \right] v_{l,\text{Gr}}^{j0},$$

$$g_{Gr}^{2,\text{Gr}}(q, \omega) = \frac{-1}{A\pi} \sum_{k=\text{Gr}} \sum_{ijl} v_{\text{Gr},k}^{0i} \text{Im} \left[ \chi_{kl}^{ij} \right] v_{l,\text{Gr}}^{j0}. \quad (50)$$

Note that even though each of the above terms may be negative, $g_{Gr}^2$ is always positive.

We observe similar trends as for the macroscopic $-\text{Im}\chi_M$ of Fig. 5. Three types of coupling emerge: one driven by BN's phonons and screened by graphene when $g_{Gr}^{2,\text{BN}} > 0$ and $g_{Gr}^{2,\text{Gr}} < 0$ in region (2); one driven by graphene's plasmon and screened by BN when $g_{Gr}^{2,\text{Gr}} > 0$ and $g_{Gr}^{2,\text{BN}} < 0$.
<figure><img src="image_13.png" /><figcaption>FIG. 11: Contributions to the coupling of graphene's electrons in 10BN/Gr@0.2eV/10BN, i.e. $g_{Gr}^{2,BN}(q, \omega)$ and $g_{BN}^{2,Gr}(q, \omega)$ as described by Eq. 50, normalized at each $q$ point. We also represent the phonon contribution to the coupling, i.e. $g_{Gr}^{2,BN}$ of Eq. 51.</figcaption></figure>

0 in region (1); and one corresponding to a true plasmon-phonon hybrid, when both $g_{Gr}^{2,BN} > 0$ and $g_{Gr}^{2,Gr} > 0$ in region (3).

Using the information contained within this separation of the contributions, one may define a dynamically screened remote electron-phonon interaction corresponding to the scattering of graphene's electrons by BN's phonons. The phonon-driven contribution to scattering can be systematically extracted, with exact asymptotic limits, as

$$g_{Gr}^{2,ph} = \begin{cases} 0 & \text{if } g_{Gr}^{2,BN} < 0, \\ g_{Gr}^{2,BN} & \text{if } g_{Gr}^{2,BN} > 0 \text{ and } g_{Gr}^{2,Gr} > 0, \\ g_{Gr}^{2} & \text{if } g_{Gr}^{2,BN} > 0 \text{ and } g_{Gr}^{2,Gr} < 0. \end{cases} \quad (51)$$

The above procedure of extracting an electron-phonon coupling is arbitrary, but it assures the coupling to be a positive definite quantity. Disentangling plasmon and phonon scattering as such is useful when considering momentum scattering for electronic transport. Despite providing the main contribution, the current procedure does not necessarily provide a full description of the momentum relaxation process. To do so, one could for example consider a system of coupled Boltzmann equations for electrons and collective modes²⁷ for dynamical couplings, and carefully extract a momentum relaxation time. We leave this for future studies.

## 2. Plasmon and phonon contributions

Fig. 12 shows the coupling of graphene's electrons with the collective modes for different setups. In the upper

panel we plot the total dynamical coupling $\tilde{g}_{Gr}^2$ (black dashed line), compared to electron-phonon part of the coupling $\tilde{g}_{Gr}^{2,ph}$ (blue line), and to $\tilde{g}_{Gr}^2$ with $Q_{ph,BN}^0 = 0$ (red line), which represents the plasmon coupling screened by BN's electrons.

The lower panel of Fig. 12 shows the expected asymptotic behaviors for the electron-phonon contribution. At large momenta, $\tilde{g}_{Gr}^{2,ph}$ coincides with $\tilde{g}_{Gr}^2$ computed with graphene's static electronic response $Q_{el,Gr}^0(\omega = 0)$. Here, the dynamical nature of the plasmon is not felt, and graphene's free carriers screen the LO phonon statistically. We thus recover the static, Born-Oppenheimer approximation of DFPT (see App. C for a direct comparison of VED and DFPT electron-phonon couplings in a BN/Graphene system in this regime). In the small momentum limit, the total coupling between electrons and collective modes can be qualitatively divided into two contributions. The dominant part to the total coupling is due to the plasmon ($\tilde{g}_{Gr}^2$ with $Q_{ph,BN}^0 = 0$). The phonon contribution, $\tilde{g}_{Gr}^{2,ph}$, is seen to approach $\tilde{g}_{Gr}^2$ with $Q_{el,Gr}^i = 0$, which represents the electron-phonon coupling obtained by completely suppressing the response from graphene but keeping the static electronic screening from BN. In other words, the phonon is too fast for graphene's electrons to respond. We will refer to the coupling in this region as 'intrinsic BN'.

We remind that, following the discussion of App. A 3, only in these asymptotic regimes the electron-phonon interaction deduced in this work may be rigorously separated from the electron-plasmon scattering and used inside a Fermi golden rule approach to compute electronic scattering rates due to phonons. In between the asymptotic regimes, following our description there is a
<figure><img src="image_14.png" /></figure>

FIG. 12: Coupling between graphene's electrons and collective modes, integrated over positive frequencies, for different setups. (Upper) $\bar{g}_{Gr}^2$: in presence of both dynamical phononic and dynamical plasmonic contributions. $\bar{g}_{Gr}^2$ with $Q_{ph, BN}^0 = 0$: in the absence of the dynamical phonon contribution from BN, but including the BN's static electronic response. This corresponds to electron-plasmon interaction screened by BN environment. $\bar{g}_{Gr}^{2,ph}$: the contribution to the total coupling due to electron-phonon when the responses are all dynamical. (Lower) $\bar{g}_{Gr}^2$ with $Q_{el, Gr}^0(\omega = 0)$: in the absence of the dynamical contribution from Gr, i.e. by using the static screening from graphene's electrons. This corresponds to remote electron-phonon coupling screened from graphene in the Born-Oppenheimer approximation, as in DFPT calculations. $\bar{g}_{Gr}^2$ with $Q_{el, Gr}^i = 0$: the response of graphene layer is set to zero, i.e. it does not respond to any perturbation. This corresponds to remote electron-phonon coupling screened only by BN's dielectric function. The vertical dashed line indicates the passage from the region where no electron-hole excitation is possible, to the intraband region.

region where there is no clear phonon-driven coupling ($\bar{g}_{Gr}^{2, BN} < 0$), which explains the dip in the electron-phonon coupling curve. Nevertheless, the full coupling with the hybrid plasmon-phonon mode could induce some

momentum relaxation, which we leave for future studies. Now that we discussed the division of the couplings, we can study the effect of the number of layers and of the Fermi level on the electron-phonon coupling.

### 3. Effects of number of BN layers

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

FIG. 13: Evolution of the electron-phonon coupling as a function of the number of BN layers.

Fig. 13 shows the evolution of the electron-phonon coupling as a function of the number of layers. The ‘intrinsic BN’ $q \to 0$ limit of the electron-phonon coupling increases as $\frac{n}{\sqrt{n}} = \sqrt{n}$, where n is the number of BN layers. This numerator comes from the sum of the dipole fields over the layers, while the denominator accounts for the the highest LO phonon eigenvector normalization$^{78}$. The coupling decreases with the number of layers at large $q$. This can be traced back to the intrinsic screening from BN, of the form $\epsilon \sim 1 + N\alpha q + o(q^2)$. There is a saturation in this decrease of the coupling at large $q$. Indeed, the reach of the interlayer Coulomb interactions goes as $e^{-qd}$ where $d$ is the distance between two given layers. At a given $q$, graphene only feels BN layers at distance such that $qd \sim 1$. At large $q$, the threshold $d$ decreases and BN layers added beyond that won’t interact with graphene’s electrons.

For the system with only 1 BN on each side and at small $q$, the coupling can actually be stronger than the ‘intrinsic BN’ coupling, just before reaching it in the $q \to 0$ limit. The phonon part of the plasmon-phonon hybrid is thus stronger than the intrinsic phonon itself. One may consider that the plasmon is driving the phonon part beyond the normal self-sustained phonon excitation. This is likely due to the relative strength of the plasmon with respect to the phonons.
<figure><img src="image_16.png" /><figcaption>FIG. 14: (Upper) Evolution of the electron-phonon coupling as a function of doping. The vertical dash-dot lines represents the limit of the electron-hole continuum. For transport, only the momenta larger than this limit mater (not a hard limit, modulated by temperature). (Lower) Asymptotic limits of the electron-phonon coupling, as in Fig. 12, but for $\varepsilon_F = 0.1$eV. The vertical dashed lines indicate the passage from the region where no electron-hole excitation is possible, to the interband and then intraband regions.</figcaption></figure>

### 4. Effects of graphene's Fermi level

We verified that the plasmon part of the coupling, which is dominant in the $q \to 0$ limit, reaches a maximum value independent of $\varepsilon_F$. Also, removing the atomic contribution to the screening, the electron-plasmon coupling was checked to be mostly independent of the Fermi level in units of $q/k_F$. Thus, differences in the total coupling at different doping levels are mostly due to the phonon contribution, which we study here.

The upper panel Fig. 14 shows the electron-phonon coupling as a function of the Fermi level in graphene. The VED framework offers an interesting opportunity to elucidate how the electron-phonon interaction is screened by free-carriers as a function of doping. As discussed further in Section IV D, this is a non-trivial question in the case of doped 2d semiconductors.

Neglecting screening from graphene is only valid in the $q \to 0$ limit, where the electron-phonon coupling tends to the intrinsic LO phonon coupling (see lower panels of Figs. 12 and 14). The vertical dashed lines in Fig. 14 show the momenta at which the electron-hole continuum is entered. For $\varepsilon_F = 0.2$ and $0.3$ eV, this coincides with the onset of strong metallic screening as obtained in the static limit of graphene's response. Only in the case of $\varepsilon_F = 0.1$ eV does the coupling maintain a large value for a considerable interval of wavevectors when entering the electron-hole continuum.

The reason can be traced back to the crossing the interband region II of the electron-hole continuum. In fact, for $\varepsilon_F = 0.1$ eV (lower panel of Fig. 14) region I is encountered when $0 < q/k_F < 0.4$, region II when $0.4 < q/k_F < 2$ and region III when $q/k_F > 2$. Instead, for $\varepsilon_F = 0.2$ (lower panel of Fig. 12) we have that region II is not encountered, and the phonons directly cross from region I to III at $q/k_F \sim 0.9$. At large momenta, for all the Fermi energies the coupling is statically screened by graphene's free carriers. The difference between the curves of Fig. 14 at $q > 2k_F$ is due to $q/k_F$ scale, since graphene's screening is constant and the coupling is independent of $k_F$.

We conclude that within the electron-hole continuum, the screening of the electron-phonon interaction by graphene's free carriers is always well approximated by its static limit in the intraband continuum. Only when the phonon dispersion crosses the interband continuum does it significantly weakens. In that situation, as shown in Fig. 8, the anti-crossing disappears to be replaced by a collective mode which mostly follows the intrinsic LO phonon dispersion. Note that the case of semi-metallic graphene is peculiar, since interband transitions are easily accessible, as discussed in Sec. IV D in details.

### C. Remote versus intrinsic scattering of graphene's electrons

As a final demonstration of the insight brought by the VED framework, we compare the remote scattering from BN's phonons to the intrinsic scattering from graphene's optical phonons. We focus here on momentum scattering, which contributes to the resistivity of graphene. Following Ref. 58, in graphene, one may define an isotropic and energy-dependent transport scattering time $\tau^{\text{Tr}}$ from which the resistivity may be expressed as:

$$ \frac{1}{\rho} = \frac{e^2}{\pi \hbar^2} \int_{-\infty}^{\infty} d\varepsilon |\varepsilon| \tau^{\text{Tr}}(\varepsilon) \left( -\frac{\partial n_{\varepsilon}^{\text{FD}}}{\partial \varepsilon} \right). \quad (52) $$

In general, $\tau^{\text{Tr}}(\varepsilon)$ is obtained by solving the full Boltzmann transport equations, i.e. a set of coupled equations relating $\tau^{\text{Tr}}$ at different energies$^{58}$. All phonons contribute to each of the coupled equations, and it is not straightforward to separate their respective contributions to $\tau^{\text{Tr}}(\varepsilon)$. However, within the relaxation time
approximation ($\tau^{\text{Tr}}(\varepsilon)$ weakly dependent on energy), the equations decouple and the scattering rate $\tau^{\text{Tr},-1}(\varepsilon)$ (i.e. the inverse of the scattering time) is trivially expressed as a sum over phonon modes. The quantitative impact on the intrinsic resistivity is negligible at room temperature and in the range of Fermi levels studied here (see Fig.6 of Ref.⁵⁸). Within this approximation, it is thus reasonable to compute and compare scattering rates from each type of phonons. $\tau^{\text{Tr},-1}$ then differs from the self-energy/ARPES scattering rates $\tau^{-1}$ (Eq. 41) by an additional $(1 - \cos\theta)$ term inside the angular integration, where $\theta$ is the scattering angle between $\mathbf{k}$ and $\mathbf{k} + \mathbf{q}$. That term gives more importance to the back-scattering transitions⁵⁸ ($\theta = \pi$). Finally, notice that those approximations preserve the electron-hole symmetry of scattering rates.

For intrinsic scattering mechanisms in graphene, we use results from App. B of Ref.⁵⁸. We focus on the scattering rates coming from LO and TO phonons at zone center ($q \sim \Gamma$) with frequency $\omega_{\text{LO}\Gamma} = \omega_{\text{TO}\Gamma} = \omega_O \sim 0.2$eV, and $A_1$ phonons at zone border ($q \sim K$) with frequency $\omega_K \sim 0.15$eV. Keeping only phonon-dependent quantities, the transport scattering rate from the sum of LO and TO phonons is proportional to

$$ \tau_{\text{O}}^{\text{tr},-1}(\varepsilon) = \frac{\hbar\beta_{\text{O}}^2}{2M\omega_{\text{O}}} W_{\pm}(\hbar\omega_{\text{O}}, \varepsilon), \quad (53) $$

$$ W_{\pm}(\hbar\omega_{\text{O}}, \varepsilon) = \left( n_{\omega_{\text{O}}}^{\text{BE}} + \frac{1}{2} \mp \frac{1}{2} \right) \times \frac{|\varepsilon \pm \hbar\omega_{\text{O}}|}{(\hbar v_{\text{F}})^2 (1 - n_{\varepsilon}^{\text{FD}})} (1 - n_{\varepsilon \pm \hbar\omega_{\text{O}}}) \quad (54) $$

where the summation over both the $\pm$ terms is intended ('+' stands for absorption and '-' for emission). $\beta_O$ is the strength of the optical electron-phonon coupling, with $\frac{\hbar\beta_O^2}{2M\omega_O} = 0.11$ eV² as obtained from GW calculations. $M$ is the mass of carbon atoms, and $n^{\text{BE}}$ and $n^{\text{FD}}$ are taken at room temperature here. The scattering rate for the zone-border phonon at K reads

$$ \tau_{\text{K}}^{\text{tr},-1}(\varepsilon) = \frac{\hbar\beta_{\text{K}}^2}{2M\omega_{\text{K}}} W_{\pm}(\hbar\omega_{\text{K}}, \varepsilon) \int \frac{d\theta}{2\pi} (1 - s_{\text{K}} \cos\theta) \times (1 - \cos\theta) = \frac{\hbar\beta_{\text{K}}^2}{2M\omega_{\text{K}}} W_{\pm}(\hbar\omega_{\text{K}}, \varepsilon)(1 + s_{\text{K}} \frac{1}{2}), \quad (55) $$

where $s_K = \text{sign}(\varepsilon)\text{sign}(\varepsilon \pm \hbar\omega_K)$ is positive (negative) for intraband (interband) transitions, and $\frac{\hbar\beta_K^2}{2M\omega_K} = 0.18$ eV² as obtained from GW calculations at large doping⁸⁵. The calculation of the correct value for $\beta_K$, in particular at small doping, is still matter of active theoretical and experimental research⁴¹,⁸⁶,⁸⁷. The value used here is a lower limit, such that the scattering from those phonons is likely much stronger in practice. The first factor in the integral comes from the angular dependence of the wavefunctions bracketing the electron-phonon interaction, with a different sign of the cosine for intraband

<table><caption>TABLE I: Ratio of remote transport scattering rates $\tau_{\text{BN}}^{\text{tr},-1}(\varepsilon_{\text{F}})$ and the intrinsic zone-center $\tau_{\text{O}}^{\text{tr},-1}(\varepsilon_{\text{F}})$ or zone-boundary $\tau_{\text{K}}^{\text{tr},-1}(\varepsilon_{\text{F}})$ ones, for a temperature of $T = 300$K.</caption><thead><tr><th>ε<sub>F</sub> (eV)</th><th>n</th><th>τ<sub>BN</sub><sup>tr,-1</sup>(ε<sub>F</sub>)/τ<sub>O</sub><sup>tr,-1</sup>(ε<sub>F</sub>)</th><th>τ<sub>BN</sub><sup>tr,-1</sup>(ε<sub>F</sub>)/τ<sub>K</sub><sup>tr,-1</sup>(ε<sub>F</sub>)</th></tr></thead><tbody><tr><td rowspan="3">0.1</td><td>1</td><td>2.24</td><td>0.20</td></tr><tr><td>10</td><td>2.80</td><td>0.25</td></tr><tr><td>30</td><td>2.23</td><td>0.20</td></tr><tr><td rowspan="3">0.2</td><td>1</td><td>0.87</td><td>0.052</td></tr><tr><td>10</td><td>0.92</td><td>0.055</td></tr><tr><td>30</td><td>0.92</td><td>0.055</td></tr><tr><td rowspan="3">0.3</td><td>1</td><td>0.67</td><td>0.04</td></tr><tr><td>10</td><td>0.67</td><td>0.04</td></tr><tr><td>30</td><td>0.78</td><td>0.05</td></tr></tbody></table>

and interband transitions. The second is the usual angular term in transport giving more importance to the back-scattering transitions.

For remote scattering from BN's phonons, we compute the scattering time as explained in App. A4. In particular, we start from Eq. 41 and insert the $(1 - \cos\theta)$ term. Then, as for the intrinsic scattering, we assume that the relevant BN's phonons are dispersionless phonons of frequency $\sim \omega_{\text{BN}} = 0.18$eV. We then consider only contributions coming from electron-phonon scattering, i.e. we consider only the coupling $\tilde{g}_{\text{Gr}}^{2,\text{ph}}$, which is isotropic in momentum space to a very good approximation. Energy conservation is imposed, and the sum over $\mathbf{q}$ points transform in an angular integral. Since energy conservation changes from emission to absorption, we define the admitted scattering wavevectors, at a given $\theta$, as $q_{\theta}^{\pm}$ (see Fig. 15). As detailed in App. A4 we obtain that the remote phonon scattering rate is proportional to

$$ \tau_{\text{BN}}^{\text{tr},-1}(\varepsilon) = W_{\pm}(\hbar\omega_{\text{BN}}, \varepsilon) \int \frac{d\theta}{2\pi} \tilde{g}_{\text{Gr}}^{2,\text{ph}}(q_{\theta}^{\pm}) \frac{1 + s_{\text{BN}} \cos\theta}{2} \times (1 - \cos\theta), \quad (56) $$

where $s_{\text{BN}} = \text{sign}(\varepsilon)\text{sign}(\varepsilon \pm \hbar\omega_{\text{BN}})$. The value of $\tilde{g}_{\text{Gr}}^{2,\text{ph}}(q_{\theta}^{\pm})$ is taken from Fig. 14. Note that the sign of the cosine in the first term of the angular integral of Eq. 56, still depending on whether the transition is intra- or interband, is reversed with respect to Eq. 55. This is easily understood noting that the zone border phonon at K leads to an intervalley transition, while BN's phonons scatter within the same K valley.

We evaluate the transport scattering rates for $\varepsilon = \varepsilon_{\text{F}} = 0.1, 0.2, 0.3$ eV for heterostructures with n-BN layers on each side, n=1,10,30, at room temperature (300K). The results are presented is Tab. I. For $\varepsilon_{\text{F}} = 0.2, 0.3$ eV, we find that $\tau_{\text{BN}}^{\text{tr},-1}(\varepsilon_{\text{F}})/\tau_{\text{O}}^{\text{tr},-1}(\varepsilon_{\text{F}}) < 1$ while $\tau_{\text{BN}}^{\text{tr},-1}(\varepsilon_{\text{F}})/\tau_{\text{K}}^{\text{tr},-1}(\varepsilon_{\text{F}}) \ll 1$. For these doping levels $\varepsilon_{\text{F}} > \hbar\omega_{\text{BN}}$, therefore the admitted scattering is intraband. The remote coupling with BN's LO phonon is always screened by graphene. To understand why in this case
remote scattering is smaller than the intrinsic zone center
optical scattering, we remark that the electron-phonon
coupling for graphene's intrinsic LO and TO phonons is
$\frac{\hbar\beta_O^2}{2M\omega_O} = 0.11 \text{ eV}^2$, independent of the scattering angle$^{58}$.
Fig. 14 shows coupling strengths around the same value
for BN's polar phonons. However, the scattering rates
also include graphene's electronic wavefunction overlap.
Since for intraband transitions this factor vanishes in
back-scattering geometry ($s_{BN} = 1$), i.e. the most im-
portant for transport, one obtains a contribution smaller
than graphene's zone center optical phonons in trans-
port, themselves scattering much less than zone border
$A_1$ phonons. This holds also for large number of BN lay-
ers, since the coupling is essentially independent of the
number of layers for large enough momenta, see Fig. 13.
Increasing doping further only weakens the relative con-
tribution from BN's phonons.

For $\epsilon_F = 0.1$ eV, the situation is different. We find that $\tau_{\text{BN}}^{\text{tr},-1}(\epsilon_F)/\tau_{\text{O}}^{\text{tr},-1}(\epsilon_F) > 1$, but still $\tau_{\text{BN}}^{\text{tr},-1}(\epsilon_F)/\tau_{\text{K}}^{\text{tr},-1}(\epsilon_F) < 1$. To understand why now remote phonons are more important than intrinsic zone-center phonons, we notice that at small doping the possibility of unscreened interband transitions increases the scattering rate due to a strong coupling, as shown in Fig. 14. Further, for interband process the overlap of graphene's wavefunctions now reaches a maximum in the back-scattering geometry ($s_{\text{BN}} = -1$). The contribution of BN's LO phonons to graphene's resistivity is thus larger and comparable to the intrinsic contribution from graphene's LO and TO phonons. However, we find that the scattering from BN's phonons remains smaller than the contribution of intrinsic optical $A_1$ phonons at $\mathbf{K}$ at room temperature, despite the overall larger coupling. This is due to the strong reduction of the remote phonon scattering at small $\theta$, and therefore small $q$, due to the $(1 - \cos\theta)^2$ angular factor, and to the role played by phonon occupation factors at room temperature. For the intraband contribution (related to phonon absorption), the final electronic states are at $\epsilon_F = 0.3$ eV, and the norm of the corresponding momenta go from $2k_F$ to $4k_F$, where the remote coupling is screened. Let us underscore the use of a conservative estimate for coupling to the optical $A_1$ phonons at $\mathbf{K}$, such that they would scatter *at least* 4 times more than remote BN phonons.

In conclusion, remote electron-phonon coupling with BN's LO phonons does not seem to explain the discrepancies seen between ab initio simulations of intrinsic phonon-limited transport in graphene and experiments in BN-capped graphene at room temperature and higher$^{58,88,89}$. However, while we don't expect the dominance of the $A_1$ phonons at $\mathbf{K}$ to be challenged, a more accurate assessment of BN's remote contributions to graphene's resistivity is left for future investigations. Indeed, BN's out-of-plane ZO phonons would also participate to electron scattering and, as mentioned before, the definition of electron-phonon coupling proposed here does not describe the full momentum relaxation process from the plasmon-phonon excitation.

D. Discussion: screening of electron-phonon
interactions in doped 2d semiconductors

In this work we specifically quantify the remote scat-
tering of graphene’s electrons from BN’s phonons. More
generally, qualitative conclusions can be drawn on the
scattering from long-range polar phonons in doped 2d
semiconductors. In particular, the VED framework sheds
some light on the screening of those phonons from the free
carriers added by doping.

Without computing the full dynamical response as done
here, there are two easily accessible limits used in the
literature: static screening or no screening. It is gen-
eral knowledge that static screening should be used at
high doping when the plasma frequency (which increases
with free carrier density) is expected to be well above
the phonon frequency, i.e. phonons are much slower than
free carriers. In the opposite small doping limit, one may
neglect free-carrier screening if the plasma frequency is
much smaller than the phonon, i.e. phonons are much
faster than free carriers. In the case of 2d materials, the
universal square root behavior of the plasmon dispersion
$\omega_{\text{pl}} \propto \sqrt{q}$ (as opposed to a constant in 3D) somewhat
complicates the picture, as there is no single plasma fre-
quency to compare to. The systems studied in this work
are then an instructive example for the 2d case.

We only consider scattering rates that involve real
states, i.e. within the electron-hole continua. Since the
scattering mechanisms sensitive to screening mostly come
from polar-optical phonons, we can assume the existence
of a constant typical phonon frequency $\omega_{\text{ph}}$, despite the
presence of a linear dispersion at small momenta. If $\omega_{\text{ph}}$
is lower than the onset of the interband electron-hole con-
tinuum $\omega_{\text{inter}}$, i.e. if $\omega_{\text{ph}} < \omega_{\text{inter}}$, then the screening of
scattering from polar phonons by free carriers is always
effective, as in the case of the lower panel of Fig. 12. In-
deed, in this case the phonon dispersion only crosses the
electron-hole continuum in the intraband Landau damp-
ing region. When $\omega_{\text{ph}} > \omega_{\text{inter}}$, the situation is simi-
lar to the lower panel of Fig. 14. While screening re-
mains effective in the intraband continuum, the phonon
also crosses the interband continuum, where free-carrier
screening vanishes or at least is weaker.

Graphene being gapless, $\omega_{\text{inter}}$ can be made very small with respect to $\omega_{\text{ph}}$ by lowering the doping. Typical high-mobility doped semiconductors, however, have a well-isolated transport band$^{90}$, and $\omega_{\text{inter}}$ remains finite even at vanishing doping. Further, other bands are far enough that $\omega_{\text{ph}} < \omega_{\text{inter}}$, in which case the strong screening from free carriers is effective. Of course, the assessment of a realistic situations is more complex, with multiple phonon frequencies and multiple onsets of interband transitions to be considered.

Nevertheless, given the above arguments, we conclude
that simply neglecting free-carrier screening in 2d doped
semiconductors is not valid, even in the low doping limit.
Yet, this approximation is often done in state-of-the-art
first-principles electronic transport calculations$^{91-99}$. In-
stead, without access to the full dynamical screening,
using the static free-carrier screening might be a more
reasonable and robust approximation<sup>19,60,90,100-102</sup>.

V. CONCLUSIONS

In summary, we have developed an ab initio computa-
tional technique capable of describing plasmons, phonons
and their mixing in Van der Waals heterostructures. We
avoid the computational effort of simulating the full het-
erostructure. Instead, single layer response functions are
computed and inserted in a system of equations that
takes into account the electrodynamic interactions be-
tween different layers. Despite the reduced computa-
tional cost, our method is able to reproduce the results
obtained by full ab initio calculations in the static limit
of electronic responses, and extends to the dynamical
case. Those developments are applied to BN-capped
graphene, where we study the features of the collective
excitations and their evolution when changing the num-
ber of BN layers and the doping level. In particular, we
show that the phonon and plasmon mode mixing gener-
ates interesting observable results, such as the stiffening
of the plasmon-phonon-polariton group velocity due to
both vanishing remote screening of BN's phonons from
graphene and the interplay with the plasmon. We also
obtain the coupling of the collective modes with elec-
trons of the heterostructure. In particular, we investi-
gate the role played by the plasmon in screening the re-
mote electron-phonon coupling between graphene's elec-
trons and BN's LO phonons. We quantitatively describe
the crossover between the regions where the coupling is
screened by graphene's free carriers (intraband region),
to the one where it is not (interband region). This implies
a contribution to graphene's low doping resistivity from
remote BN's LO phonons that is similar to graphene's
zone-center optical phonons, but it is still negligible with
respect to the scattering from intrinsic zone-border opti-
cal phonons. Thus, it is not sufficient to explain the large
discrepancy between experiments and state-of-the-art ab-
initio simulations<sup>58,88,89,103</sup>. For a general heterostruc-
ture, one expects the interplay between plasmons and
phonons to increase in complexity. This may be due to a
complex band structure of the metallic layer presenting
e.g. subbands<sup>104</sup>, intersubband or multi-valley<sup>105</sup> plas-
mons, or to an increased importance of cross-talks be-
tween in-plane and out-of-plane electrostatics<sup>55</sup>. Both of
these cases are exciting perspectives for this work.

ACKNOWLEDGMENTS

We thank A. Guandalini, P. Barone, M. Stengel, M.
Gibertini and G. Cassabois for useful discussions. F.
Mauri and F. Macheda acknowledge the MORE-TEM
ERC-SYN project, grant agreement no. 95121. F. Mauri
and F. Macheda acknowledge the EuroHPC Joint Under-

taking for awarding this project access to the EuroHPC
supercomputer LUMI, hosted by CSC (Finland) and the
LUMI consortium through a EuroHPC Regular Access
call. First-principles calculations were performed with
the support of MESO@LR-Platform at the University of
Montpellier.
# Appendix A: Theoretical aspects of dispersions, couplings and scattering rates

## 1. Basic definitions

In these Sections we restore the possible anisotropy of the system using the vectori symbols (as e.g. **q** and **r**) for momentum and space related quantities. We define a generic retarded response as

$$f^R(\mathbf{r}, \mathbf{r}', t) = -\frac{i}{\hbar} \theta(t) \langle [\hat{A}(\mathbf{r}, t), \hat{A}(\mathbf{r})] \rangle \quad (A1)$$

where $\langle \rangle$ indicates the average over the ground state (or the thermal average), $\theta$ is the Heaviside function and $\hat{A}$ is a generic operator. Its Fourier transform satisfies the property

$$f^R(\mathbf{q}, \omega) = [f^R(\mathbf{q}, -\omega)]^* . \quad (A2)$$

The retarded density-density response function in the random phase approximation (RPA) is expressed as a function of the irreducible independent-particle density-density response $\chi_{\text{irr}}^R$ as in Eq. 26. The non-interacting retarded Green function for a phonon of branch $\mu$ at zero temperature may be written as a function of the harmonic oscillator ladder operators $a, a^\dagger$ in real time as

$$D_{\mu}^{R}(\mathbf{q}, t) = -i\theta(t) \langle a_{\mu\mathbf{q}}^{\dagger}(t) a_{\mu\mathbf{q}} + a_{\mu-\mathbf{q}}(t) a_{\mu-\mathbf{q}}^{\dagger} \rangle, \quad (A3)$$

and in the frequency space as

$$\begin{aligned} D_{\mu}^{R}(\mathbf{q}, \omega) = & \frac{2\hbar\omega_{\mu\mathbf{q}}}{(\hbar\omega + i\hbar\eta)^2 - \hbar^2\omega_{\mu\mathbf{q}}^2} = \\ & -\frac{1}{\hbar\omega + \hbar\omega_{\mu\mathbf{q}} + i\hbar\eta} + \frac{1}{\hbar\omega - \hbar\omega_{\mu\mathbf{q}} + i\hbar\eta}. \end{aligned} \quad (A4)$$

As shown later, the expression of the total inverse dielectric function contains the density-density response and the phonon propagator, both in the retarded flavour. We drop the retarded apex everywhere in the text, so that any Green/response function is to be intended as retarded. We also remind that to deduce the finite temperature expression of the retarded Green/response functions, one computes their time ordered expression in the Matsubara formalism and then performs the substitution $i\hbar\omega_n \rightarrow \hbar\omega + i\hbar\eta$ $^{64}$. As a consequence of using retarded Green functions, all the computed quantities will be analytical and will therefore respect the Kramers-Kronig relations$^{106}$.

## 2. Long-range atomic-mediated inverse dielectric function/screened Coulomb

The goal of this Section is to justify the form of Eq. 28 for a wide gap semiconductor such as BN. Therefore, we here perform the derivation of the long-range atomic-mediated inverse dielectric function and the related screened Coulomb interaction. We start from the

derivation for the 3d case. In this case, following Ref. 107 we use the RPA approximation for the electronic response and obtain that the atomic-mediated inverse dielectric function can be expressed in terms of effective charges. We then use the conclusions and procedures of Ref. 60 to extend the results to the case of a 2d single layer of small width.

### a. 3d expression for the inverse dielectric function

The screened electronic interaction is defined as$^4$

$$v(\mathbf{r}) = \frac{4\pi e^2}{|\mathbf{r}|}, \quad (A6)$$

$$w(41) = \int d\mathbf{r}'' \epsilon^{-1}(4, \mathbf{r}'', t_1) v(\mathbf{r}'' - \mathbf{r}_1), \quad (A7)$$

where the compact notation $1 = (\mathbf{r}_1, t_1)$ has been introduced. The total inverse dielectric function—and consequently the Coulomb potential—is usually split into ‘electronic’ (clamped-ion) and ‘phononic’ (total minus clamped-ion) contributions as

$$\epsilon^{-1} = \epsilon_{\text{el}}^{-1} + \epsilon_{\text{ph}}^{-1}, \quad w = w_{\text{el}} + w_{\text{ph}}. \quad (A8)$$

Such separation bring great simplification if vertex corrections that couple the systems are small$^{64}$, as assumed here. The phononic part of the screened Coulomb interaction may be written, following Refs. 108–110, as

$$\begin{aligned} \epsilon_{\text{ph}}^{-1}(12) = & \sum_{s\alpha p \atop s'\alpha' p'} \int d(34) \epsilon_{\text{el}}^{-1}(13) (\nabla_{\mathbf{r}_3})_{\alpha} V_s(\mathbf{r}_3 - \mathbf{u}_s^p) \times \\ & D_{s\alpha p, s'\alpha' p'}(t_3 t_4) \frac{\nabla_{\mathbf{r}_4}^2}{4\pi e^2} \epsilon_{\text{el}}^{-1}(42) (\nabla_{\mathbf{r}_4})_{\alpha'} V_{s'}(\mathbf{r}_4 - \mathbf{u}_{s'}^{p'}). \end{aligned} \quad (A9)$$

where $D$ is the dressed phonon propagator. $\mathbf{u}_s^p = \mathbf{R}_p + \tau_s$ are the atomic equilibrium positions, which in reciprocal space representation read as

$$\mathbf{u}_s^q = \sum_p e^{i\mathbf{q} \cdot (\mathbf{R} + \tau_s)} \mathbf{u}_s^p. \quad (A10)$$

$V_s$ indicates, in the all-electron formalism, the ion-electron interaction due to the ion at site $s$, and its Fourier transform is

$$V_{s\alpha}(\mathbf{q} + \mathbf{G}) = \frac{1}{e} v(\mathbf{q} + \mathbf{G}) \rho_{s\alpha}^{\text{ext}}(\mathbf{q} + \mathbf{G}), \quad (A11)$$

$$\rho_{s\alpha}^{\text{ext}}(\mathbf{q} + \mathbf{G}) = -i \frac{\mathcal{Z}_s e}{V} [(q_\alpha + G_\alpha)e^{-i\mathbf{G} \cdot \tau_s}], \quad (A12)$$

where $\mathcal{Z}_s$ are the ionic charges. We assume to be looking at times large enough that all the response functions depend only on relative time, the retardation condition restoring the causality at the end. We express the phonon propagator as
$$D_{s\alpha p, s'\alpha' p'}(\omega) = \frac{1}{\sqrt{M_s M_{s'}}} \sum_{\mu} \int \frac{d\mathbf{q}}{V_{\text{BZ}}} \frac{e_{\mu\mathbf{q}}^{s\alpha} e^{i\mathbf{q}\cdot(\mathbf{R}_p + \tau_s)} e_{\mu\mathbf{q}}^{s'\alpha'} e^{-i\mathbf{q}\cdot(\mathbf{R}_{p'} + \tau_{s'})}}{(\hbar\omega + i\hbar\eta)^2 - \hbar^2\omega_{\mu\mathbf{q}}^2}, \quad (\text{A13})$$

where $e_{\mu\mathbf{q}}^{s\alpha}$ are eigenvectors of the dynamical matrix and $\omega_{\mu\mathbf{q}}$ its eigenvalues. Notice that we have written the phonon propagator in the 'clean limit', i.e. we assume a delta-like spectral function. In BN (bulk or single layer) the clean limit is respected. Notice also that the dynamical frequencies coincide with the adiabatic ones, since it is a large gap semiconductor. We now introduce the in-plane screened and unscreened charge densities, and related effective charges, in the RPA approximation for $\epsilon_{\text{el}}^{-160,107,111-114}$. These read as

An explicit expression for the effective charges is found as

$$\rho_{s\alpha}(\mathbf{q}, \omega) = \epsilon_{\text{el}}^{-1}(\mathbf{q}, \omega) \bar{\rho}_{s\alpha}(\mathbf{q}, \omega), \quad (\text{A14})$$

$$\bar{Z}_{s\alpha}(\mathbf{q}, \omega) = i \frac{V}{eq} \bar{\rho}_{s\alpha}(\mathbf{q}, \omega). \quad (\text{A15})$$

$$\bar{Z}_{s\alpha}(\mathbf{q}, \omega) = \frac{\bar{Z}_s}{qv(\mathbf{q})} \times \sum_{\mathbf{G}} \frac{\epsilon_{\text{el}}^{-1}(\mathbf{q}, \mathbf{q} + \mathbf{G}, \omega)}{\epsilon_{\text{el}}^{-1}(\mathbf{q}, \omega)} (q_{\alpha} + G_{\alpha}) v(\mathbf{q} + \mathbf{G}) e^{-i\mathbf{G} \cdot \tau_s}, \quad (\text{A16})$$

where $\mathbf{G}$ are reciprocal lattice vectors. Using all the above expressions, and Fourier transforming Eq. A9 we find, for the atomic-mediated macroscopic component of the Coulomb interaction

$$w_{\text{ph}}(\mathbf{q}, \omega) = -\frac{1}{V} q^2 [w_{\text{el}}(\mathbf{q}, \omega)]^2 \frac{1}{\sqrt{M_s M_{s'}}} \sum_{\mu} \frac{(\sum_{s\alpha} \bar{Z}_{s\alpha}(\mathbf{q}, \omega) e_{\mu\mathbf{q}}^{s\alpha}) (\sum_{s'\alpha'} \bar{Z}_{s'\alpha'}(-\mathbf{q}, \omega) e_{\mu-\mathbf{q}}^{s'\alpha'})}{(\hbar\omega + i\hbar\eta)^2 - \hbar^2\omega_{\mu\mathbf{q}}^2} = \\ \frac{V}{e^2} [w_{\text{el}}(\mathbf{q}, \omega)]^2 \frac{1}{\sqrt{M_s M_{s'}}} \sum_{\mu} \frac{(\sum_{s\alpha} \bar{\rho}_{s\alpha}(\mathbf{q}, \omega) e_{\mu\mathbf{q}}^{s\alpha}) (\sum_{s'\alpha'} \bar{\rho}_{s'\alpha'}(-\mathbf{q}, \omega) e_{\mu-\mathbf{q}}^{s'\alpha'})}{(\hbar\omega + i\hbar\eta)^2 - \hbar^2\omega_{\mu\mathbf{q}}^2}, \quad (\text{A17})$$

where

$q$, one can rewrite Eq. A17 as

$$e_{\mu-\mathbf{q}}^{s'\alpha'} = e_{\mu\mathbf{q}}^{s'\alpha'*}. \quad (\text{A18})$$

For an isolant such as BN, for frequencies much smaller than the band gap, we have

$$w_{\text{ph}}(\mathbf{q}, \omega) = \frac{1}{V} \sum_{s\alpha\beta, s'\alpha'\beta'} [w_{\text{el}}(\mathbf{q})]^2 \frac{1}{\sqrt{M_s M_{s'}}} \times \\ \sum_{\mu} \frac{\left(\sum_{s\alpha\beta} \bar{Z}_{s\alpha\beta} q\beta e_{\mu\mathbf{q}}^{s\alpha}\right) \left(\sum_{s'\alpha'\beta'} \bar{Z}_{s'\alpha'\beta'} q\beta' e_{\mu-\mathbf{q}}^{s'\alpha'}\right)}{(\hbar\omega + i\hbar\eta)^2 - \hbar^2\omega_{\mu\mathbf{q}}^2}, \quad (\text{A20})$$

$$\epsilon_{\text{el}}^{-1}(\mathbf{q}, \omega) \sim \epsilon_{\text{el}}^{-1}(\mathbf{q}), \quad \bar{Z}(\mathbf{q}, \omega) \sim \bar{Z}(\mathbf{q}). \quad (\text{A19})$$

The unscreened effective charges have an expansion of the form$^{60,107}$

$$\bar{Z}_{s\alpha}(\mathbf{q}) = \frac{q\beta}{q} \bar{Z}_{s\alpha\beta} - \frac{i q\beta}{2q} q\gamma Q_{s\alpha\beta\gamma} + \frac{1}{3!} \frac{q\beta}{q} q\gamma q\delta O_{s\alpha\beta\gamma\delta} + \mathcal{O}(q^3),$$

where the various orders of the expansion are the Born effective charges, dynamical quadrupoles and octupoles, as so on. With the above expansion, at leading order in

which is consistent with the form obtained in Ref.$^4$. Eq. A20 includes only the leading order Frölich coupling$^{115-119}$, whereas Eq. A17 include quadrupoles and higher multipolar order expansions$^{92,120,121}$. Using the form of Eq. (28) of Ref.$^{60}$, one can rewrite the equivalent of Eq. A17 for the inverse dielectric matrix in the static case as a function of the long-range component of the dynamical matrix:

$$\epsilon_{\text{ph}}^{-1}(\mathbf{q}, \omega) = \epsilon_{\text{el}}^{-1}(\mathbf{q}) \sum_{\mu} \frac{\mathbf{e}_{\mu\mathbf{q}} \mathcal{D}^{\text{L}}(\mathbf{q}) \mathbf{e}_{\mu\mathbf{q}}^{*}}{(\hbar\omega + i\hbar\eta)^2 - \hbar^2\omega_{\mu\mathbf{q}}^2}. \quad (\text{A21})$$
b. 2d case for a single BN single layer

As in the rest of this work, for a BN single layer we only consider in-plane displacements of atoms, and suppose that there are no cross-talks between in-plane and out-of-plane components of the response functions. The error of this assumption, since our single layers always present in-plane mirror symmetry, is of order $O(q^2z^2)$. As explained in Ref.⁶⁰, if one also assume that the material thickness is much smaller than the wavelength of the perturbation (thin limit), one can then forget about the out-of-plane direction of the electrostatic problem and, for the formulae of App. A 2 a, replace the three-dimensional Coulomb kernel with its two dimensional counterpart, and volume-intensive quantities with the area-intensive ones. In this work, we extend the validity of 2d formulae beyond the thin limit by using Eq. 16 to describe the 2d Coulomb interaction. We finally obtain for the LO phonon, assuming a small finite lifetime

$$[\epsilon_{1L,ph}^{-1}]^0(\mathbf{q}, \omega) = [\epsilon_{1L,el}^{-1}]^0(\mathbf{q}) \frac{\mathbf{e}_{LOq} \mathcal{D}^L(\mathbf{q}) \mathbf{e}_{LOq}^*}{(\hbar\omega + i\hbar\eta_{LO})^2 - \hbar^2\omega_{LOq}^2}, \quad (A22)$$

$$[\epsilon_{1L,el}^{-1}]^0(\mathbf{q}) = 1 + v_{BN, BN}(\mathbf{q})Q_{el}^0(\mathbf{q}), \quad (A23)$$

$$[\epsilon_{1L,ph}^{-1}]^1(\mathbf{q}, \omega) = 0. \quad (A24)$$

The above equations are the one used, for the isotropic case, in Sec. II A 6. The explicit expression of the numerator of Eq. A22 for the thin limit is⁶⁰

$$\mathbf{e}_{LOq} \mathcal{D}^L(\mathbf{q}) \mathbf{e}_{LOq}^* = \frac{2\pi A}{q} \times \left( \sum_{s\alpha} \rho_{s\alpha}(\mathbf{q}) e_{LOq}^{s\alpha} \right) \left( \sum_{s'\alpha'} \bar{\rho}_{s'\alpha'}(\mathbf{q}) e_{LOq}^{s'\alpha'} \right)^*. \quad (A25)$$

In the approximation of this work that the atomic profile density is the same of the electronic one, we extend the above expression to

$$\mathbf{e}_{LOq} \mathcal{D}^L(\mathbf{q}) \mathbf{e}_{LOq}^* = A v_{BN,BN}^{00} \times \left( \sum_{s\alpha} \rho_{s\alpha}(\mathbf{q}) e_{LOq}^{s\alpha} \right) \left( \sum_{s'\alpha'} \bar{\rho}_{s'\alpha'}(\mathbf{q}) e_{LOq}^{s'\alpha'} \right)^*. \quad (A26)$$

Using the Poisson equation we obtain an equivalent formula

$$\mathbf{e}_{LOq} \mathcal{D}^L(\mathbf{q}) \mathbf{e}_{LOq}^* = \frac{A}{v_{BN,BN}^{00}(\mathbf{q})} \int f_{BN}^0(z) \bar{V}_{LO}(\mathbf{q}, z) dz \quad (A27) \\ \times \int f_{BN}^0(z) V_{LO}^*(\mathbf{q}, z) dz.$$

In the above equation we have used the abbreviation $\bar{V}_{LO}(\mathbf{q}, z) = \sum_{s,\alpha} \frac{e_{LOq}^{s,\alpha}}{\sqrt{M_s}} \frac{\partial \bar{V}_{KS}(\mathbf{q}, z)}{\partial u_{s,\alpha}}$ is the projection on the LO mode of $\frac{\partial \bar{V}_{KS}(\mathbf{q}, z)}{\partial \mathbf{u}^3}$ (and the same for non-barred potential; 'KS' indicates the Khon-Sham potential). $\bar{V}_{LO}$ is computed in this work as explained in Ref.⁶⁰. The approximation of A27 can be bettered assuming that the difference between longitudinal and tranverse optical phonon frequencies is all attributable to the long-range component of the dynamical matrix—i.e. short range interactions differences between LO and TO are negligible. This amounts to writing

$$\mathbf{e}_{LOq} \mathcal{D}^L(\mathbf{q}) \mathbf{e}_{LOq}^* = \hbar^2 \omega_{LOq}^2 - \hbar^2 \omega_{TOq}^2, \quad (A28)$$

as it is customarily done e.g. for 3d cases. This is the same approximation included in Eq. 30 in the isotropic case. The difference between using Eq. A28 and Eq. A27 is commented in App. B. Inserting Eq. A28 in Eq. A22, we end up with

$$[\epsilon_{1L,ph}^{-1}]^0(\mathbf{q}, \omega) = [\epsilon_{1L,el}^{-1}]^0(\mathbf{q}) \frac{\hbar^2 \omega_{LOq}^2 - \hbar^2 \omega_{TOq}^2}{(\hbar\omega + i\hbar\eta_{LO})^2 - \hbar^2 \omega_{LOq}^2}, \quad (A29)$$

which in the isotropic case directly leads to Eq. 32. If both the positive and negative phonon spectral functions are very sharp, for positive frequencies we have

$$-\mathrm{Im} \left[ \epsilon_{1L,ph}^{-1} \right]^0 (\mathbf{q}, \omega) \sim_{\omega \sim \omega_{LOq}} \pi \left[ \epsilon_{1L, el}^{-1} \right]^0 (\mathbf{q}) \times \frac{\hbar^2 \omega_{LOq}^2 - \hbar^2 \omega_{TOq}^2}{2\hbar\omega_{LOq}} \delta(\hbar\omega - \hbar\omega_{LOq}). \quad (A30)$$

In the above derivations we have always treated the phonon system as non-interacting; this hypothesis may be lifted by dressing the phonon propagator with the appropriate self-energies. We simplified the treatment by using for $\eta$ a finite constant value, i.e. $\eta_{LO}$.

c. Rewriting atomic-mediated screened Coulomb as a function of electron-phonon in the static case

In the same approximations of App. A 2 a and A 2 b, Eq. A17 can be rewritten in an interesting fashion for the static case via the Frölich electron-phonon coupling stripped of the Bloch functions overlaps. In fact, using Eq. (36) of Ref.⁶⁰ we define

$$g_{LOq}^{\mathrm{Fr}} = \frac{i}{A} w_{\mathrm{el}}(\mathbf{q}) \sum_{s\alpha\beta} q_{\beta} \bar{Z}_{s\alpha\beta} e_{LOq}^{s\alpha} \times \left( \frac{\hbar}{2M_s \omega_{LOq}} \right)^{1/2}. \quad (A31)$$

Eq. A17 then becomes

$$w_{\mathrm{ph}}(\mathbf{q}, \omega) = A |g_{\mathrm{LOq}}^{\mathrm{Fr}}|^2 \frac{2\hbar\omega_{\mathrm{LOq}}}{(\hbar\omega + i\hbar\eta_{\mathrm{LO}})^2 - \hbar^2\omega_{\mathrm{LOq}}^2}, \quad (A32)$$
with

$$-\mathrm{Im}[w_{\mathrm{ph}}(\mathbf{q}, \omega)] \underset{\omega \sim \omega_{\mathrm{LOq}}}{\sim} \pi A |g_{\mathrm{LOq}}^{\mathrm{Fr}}|^2 \delta(\hbar\omega - \hbar\omega_{\mathrm{LOq}}). \qquad (\mathrm{A33})$$

This means that the weight of the poles of $w_{ph}(\mathbf{q}, \omega_{LO})$ is proportional to the squared modulus of the Frölich coupling.

In the general case where the static approximation for the electronic screening cannot be performed, the imaginary part of the screened Coulomb interaction will contain both electron-electron and electron-phonon contributes. The interplay of the electronic system and the phononic one determines the dispersion of plasmon and phonon poles, and their coupling.

### 3. Scattering rates and $g^2$

The scattering rate $\tau^{-1}$ of electrons relates to the retarded electronic self-energy as

$$-\mathrm{Im}\Sigma_{\mathrm{el}} = \frac{\hbar}{2}\tau^{-1}. \qquad (\mathrm{A34})$$

Following Refs. 23 and 27, one can deduce the depopulation rate of electrons, i.e. scattering times, using the structure factor and the fluctuation-dissipation theorem, obtaining an expression proportional to $-\mathrm{Im}[\epsilon^{-1}(\mathbf{q})]v(\mathbf{q})$. We generalize such 3d expressions to the 2d case. We take into account the matrixial ordering of the inverse dielectric matrix and Coulomb kernel (not relevant in the purely 3d case) by noting that the electronic self-energy reads 122

$$\Sigma_{\mathrm{el}}(12) = i\hbar \int d(34)G(13)\Gamma(324)w(41^+). \qquad (\mathrm{A35})$$

The scattering rate for a Bloch electron $nk$ is then obtained as

$$k_B T \frac{\partial f_{nk}}{\partial \epsilon_{nk}} \tau_{nk}^{-1} = \frac{2\pi}{\hbar} \frac{1}{AN_q} \sum_q \int_{-\infty}^{\infty} d\omega \int d\mathbf{r} dz d\mathbf{r}' dz' d\bar{z} \times \\ n_{\epsilon_{nk}}^{\mathrm{FD}} (1 - n_{\epsilon_{m\mathbf{k}+\mathbf{q}}}^{\mathrm{FD}}) (1 + n_{\omega}^{\mathrm{BE}}) \times \\ u_{nk}^{*}(\mathbf{r}, z) u_{nk}(\mathbf{r}', z') u_{m\mathbf{k}+\mathbf{q}}(\mathbf{r}, z) u_{m\mathbf{k}+\mathbf{q}}^{*}(\mathbf{r}', z') \times \\ \mathrm{Im} \left[ -\epsilon^{-1}(\mathbf{q}, z, \bar{z}, \omega) \frac{v(\mathbf{q}, \bar{z}, z')}{\pi} \right] \delta(\hbar\omega + \epsilon_{m\mathbf{k}+\mathbf{q}} - \epsilon_{nk}). \qquad (\mathrm{A36})$$

where we have assumed that final states are represented by non-interacting electrons. Depending on the derivation, the overlap factor comes from bringing the self-energy in the basis of the Bloch basis functions, and taking its diagonal components (see e.g. Eq. (157) of Ref. 122), or from the connection between the scattering rate and the structure factor (see e.g. Eq. (5) of Ref. 23). The integral over positive frequencies will produce the

scattering due to phonon/plasmon emission, while the one over negative frequencies is related to absorption because $1 + n_{\omega}^{\mathrm{BE}} = -n_{\omega}^{\mathrm{BE}}$, where the minus sign is compensated by the property A2 applied on the inverse dielectric matrix. To connect with usual formulae, by using the algebraic relations discussed in the Supplementary Material of Ref. 38, one finds

$$\tau_{nk}^{-1} = \frac{2\pi}{\hbar} \frac{1}{AN_q} \sum_q \int_{-\infty}^{\infty} d\omega \int d\mathbf{r} dz d\mathbf{r}' dz' d\bar{z} \mathrm{sign}(\omega) A(\omega) \times \\ u_{nk}^{*}(\mathbf{r}, z) u_{nk}(\mathbf{r}', z') u_{m\mathbf{k}+\mathbf{q}}(\mathbf{r}, z) u_{m\mathbf{k}+\mathbf{q}}^{*}(\mathbf{r}', z') \times \\ \mathrm{Im} \left[ -\epsilon^{-1}(\mathbf{q}, z, \bar{z}, \omega) \frac{v(\mathbf{q}, \bar{z}, z')}{\pi} \right] \delta(\hbar\omega + \epsilon_{m\mathbf{k}+\mathbf{q}} - \epsilon_{nk}), \\ A(\omega) = \left[ n_{\omega}^{\mathrm{BE}} + \frac{1}{2} + \mathrm{sign}(\omega) \frac{1}{2} - \mathrm{sign}(\omega) n_{\epsilon_{m\mathbf{k}+\mathbf{q}}}^{\mathrm{FD}} \right]. \qquad (\mathrm{A37})$$

In the dual basis set used in this work we have

$$\int d\bar{z} \epsilon^{-1}(\mathbf{q}, z, \bar{z}, \omega) v(\mathbf{q}, \bar{z}, z') = \sum_{\substack{ii'j'j \\ kk'l'l}} \phi_k^i(z) v_{kk'}^{ii'}(\mathbf{q}) \chi_{k'l'l}^{i'j'}(\mathbf{q}, \omega) v_{l'l}^{j'j}(\mathbf{q}) \phi_l^j(z'). \qquad (\mathrm{A38})$$

In the scattering rate expression we have pairs of cell-periodic Bloch functions which are building blocks for the density, which gets coherently bracketed with the basis set of Eq. 4. In fact, Eq. A37 becomes

$$\tau_{nk}^{-1} = \sum_{\substack{ii'j'j \\ kk'l'l}} \frac{2\pi}{\hbar} \frac{1}{AN_q} \sum_q \int_{-\infty}^{\infty} d\omega \langle u_{nk} | \phi_k^i | u_{m\mathbf{k}+\mathbf{q}} \rangle \times \\ \langle u_{m\mathbf{k}+\mathbf{q}} | \phi_l^j | u_{nk} \rangle \mathrm{sign}(\omega) A(\omega) \times \\ \mathrm{Im} \left[ -\frac{1}{\pi} v_{kk'}^{ii'}(\mathbf{q}) \chi_{k'l'l}^{i'j'}(\mathbf{q}, \omega) v_{l'l}^{j'j}(\mathbf{q}) \right] \delta(\hbar\omega + \epsilon_{m\mathbf{k}+\mathbf{q}} - \epsilon_{nk}). \qquad (\mathrm{A39})$$

#### a. BN single layer in the static case: electron-phonon

The scattering rate due to electron-phonon coupling for a BN single layer is obtained by using that the density-density response function is of the form of Eq. 12, where the atomic-mediated contribution to the response is obtained from Eqs. A22 and A24. Due to in-plane mirror symmetry, the unperturbed Bloch wavefunctions are even with respect to the layer position, so that in Eq. A39 one has the selection rule $i = j = 0$. The scattering rate may then be written as

$$\tau_{nk}^{-1,\mathrm{LO}} = \frac{1}{A} \frac{2\pi}{\hbar} \int_{-\infty}^{\infty} \hbar d\omega |\langle u_{nk} | u_{m\mathbf{k}+\mathbf{q}} \rangle|^2 \mathrm{sign}(\omega) A(\omega) \times \\ \mathrm{Im} \left[ -[\epsilon_{1\mathrm{L,ph}}^{-1}]^0 (\mathbf{q}, \omega) \frac{v_{\mathrm{BN,BN}}^{00}(\mathbf{q})}{\pi} \right] \delta(\hbar\omega + \epsilon_{m\mathbf{k}+\mathbf{q}} - \epsilon_{nk}), \qquad (\mathrm{A40})$$
which using Eqs. A32 and A5 becomes the standard expression

$$ \tau_{n\mathbf{k}}^{-1,\text{LO}} = \frac{2\pi}{\hbar} \frac{1}{N_{\mathbf{q}}} \sum_{\mathbf{q}} |g_{\text{LOq}}^{\text{Fr}}|^2 |\langle u_{n\mathbf{k}} | u_{m\mathbf{k}+\mathbf{q}} \rangle|^2 \times \\ \left[ \left(1 + n_{\omega_{\text{LOq}}}^{\text{BE}} - n_{\varepsilon m\mathbf{k}+\mathbf{q}}^{\text{FD}}\right) \delta(\hbar\omega_{\text{LOq}} + \varepsilon_{m\mathbf{k}+\mathbf{q}} - \varepsilon_{n\mathbf{k}}) + \right. \\ \left. \left(n_{\omega_{\text{LOq}}}^{\text{BE}} + n_{\varepsilon m\mathbf{k}+\mathbf{q}}^{\text{FD}}\right) \delta(-\hbar\omega_{\text{LOq}} + \varepsilon_{m\mathbf{k}+\mathbf{q}} - \varepsilon_{n\mathbf{k}}) \right]. \quad (\text{A41}) $$

b. Heterostructures

In general, one can express the imaginary part of the screened Coulomb as a squared modulus times a delta (with a minus sign), i.e. as in Eq. A33, only when one the electronic dielectric function can be treated as static. In the heterostructures studied in this work, this is not always the case, due to the coupling of plasmons and phonons. Nevertheless, there exist limits where one can use $w_{el}(\mathbf{q}, \omega) \sim w_{el}(\mathbf{q}, 0)$. They correspond to the cases where the coupling is screened statically by graphene and BN (large $|\mathbf{q}|$ limit), or just by BN (small $|\mathbf{q}|$ limit). Only in these asymptotic limits, decoupling the electron-phonon and electron-plasmon contribution to the total coupling is not ambiguous. E.g., one could determine determine a screened electron-phonon coupling directly using $\epsilon_{el}^{-1123}$. The same conclusions can be obtained looking at the electronic contribution to the screened Coulomb kernel, for a system without atomic polarization. In that case, one finds that the pole residue is related to the electron-plasmon coupling$^{124}$.

The general expression Eq. A39 can nevertheless be simplified in some situations. Bloch overlaps are important to understand selection rules and the localization of the indexes $k$ and $l$. If e.g. a Bloch electron is delocalized in the out-of-plane direction over the whole structure, such as in the case of a purely BN heterostructure, brackets with $k \neq l$ are all relevant. If instead we are interested in the electron-phonon coupling between graphene's Bloch electrons, which are localized on the graphene's plane, then we must have $k = l = \text{Gr}$, since the overlap vanish for the other cases. Moreover, as done for the single layer, we still assume the unperturbed Bloch wavefunctions are even with respect to the layer position, so that in Eq. A39 one has $i = j = 0$. The result of this operation is

$$ \tau_{n\mathbf{k}}^{-1} = \frac{2\pi}{\hbar} \frac{1}{AN_{\mathbf{q}}} \sum_{\mathbf{q}} \sum_{kk'l'l} \langle u_{n\mathbf{k}} | \phi_k^0 | u_{m\mathbf{k}+\mathbf{q}} \rangle \langle u_{m\mathbf{k}+\mathbf{q}} | \phi_l^0 | u_{n\mathbf{k}} \rangle \times \\ \int_{-\infty}^{\infty} d\omega \operatorname{sign}(\omega) \mathcal{A}(\omega) \operatorname{Im} \left[ -\frac{1}{\pi} v_{kk'}^{0i'}(\mathbf{q}) \chi_{k'l'}^{i'j'}(\mathbf{q}, \omega) v_{l'l}^{j'0}(\mathbf{q}) \right] \times \\ \delta(\hbar\omega + \varepsilon_{m\mathbf{k}+\mathbf{q}} - \varepsilon_{n\mathbf{k}}). \quad (\text{A42}) $$

In the above expression, in the asymptotic limits where

we can decouple the electron-phonon from the electron-electron interaction, $-\operatorname{Im}\nu\chi\nu$ reduces to the squared modulus of the electron-phonon coupling. Motivated by this observation, we recast the above expression with the following notation

$$ \tau_{n\mathbf{k}}^{-1} = \frac{2\pi}{\hbar} \frac{1}{N_{\mathbf{q}}} \sum_{\mathbf{q}} \int_{-\infty}^{\infty} d\omega \mathcal{A}(\omega) \delta(\hbar\omega + \varepsilon_{m\mathbf{k}+\mathbf{q}} - \varepsilon_{n\mathbf{k}}) \times \\ \sum_{kl} \langle u_{n\mathbf{k}} | \phi_k^0 | u_{m\mathbf{k}+\mathbf{q}} \rangle \langle u_{m\mathbf{k}+\mathbf{q}} | \phi_l^0 | u_{n\mathbf{k}} \rangle g_{kl}^2(\mathbf{q}, \omega), \quad (\text{A43}) $$

where $g_{kl}^2(\mathbf{q}, \omega)$ is given by Eq. 42.

4. Practical expression for the remote phonon transport scattering rates

Eq. A39 is the most general expression for the scattering rate of a heterostructure. Neglecting the scalar product of Bloch functions with $\phi^1$, one obtains Eq. A43.

<figure><img src="image_17.png" /><figcaption>FIG. 15: Scattering process for an electron with an initial energy $\epsilon$. Graphical representation of the isoenergetic surface at $\epsilon$ and $\epsilon + \hbar\omega_{ph}$, on which the electron ends up after an absorption of a phonon of frequency $\omega_{ph}$. In this case, we are representing an electron-doped Dirac cone. Several different setups of the Fermi surface and the isoenergetic surfaces are possible in dependence of the relative magnitude of $\epsilon$ and $\omega_{ph}$, and considering also the emission case.</figcaption></figure>

Here, we are interested in the scattering rates where the initial and final states are electrons of graphene, and the scattering is mediated by remote phonons. In the Dirac cone approximation, we take $\mathbf{k}$ to be a vector around $\mathbf{K}$. Without loss of generality due to isotropy, we consider a $\mathbf{k}$ on the horizontal axis, as depicted in Fig. 15. In Fig. 15 we consider initial electrons belonging to isosurface $\epsilon$. We then define the $\theta$ angle as the scattering angle between $\mathbf{k}$ and $\mathbf{k} + \mathbf{q}$, i.e. $\theta = \theta(\mathbf{q})$. The transport scattering time due to remote phonon scattering is then
obtained as

$$ \tau_{n\mathbf{k},-1} = \frac{2\pi}{\hbar} \frac{1}{N_q} \sum_{\mathbf{q}} \int_{-\infty}^{\infty} d\omega A(\omega) \delta(\hbar\omega + \epsilon_{m\mathbf{k}+\mathbf{q}} - \epsilon_{n\mathbf{k}}) \times \\ \sum_{kl} \langle u_{n\mathbf{k}} | \phi_k^0 | u_{m\mathbf{k}+\mathbf{q}} \rangle \langle u_{m\mathbf{k}+\mathbf{q}} | \phi_l^0 | u_{n\mathbf{k}} \rangle g_{Gr}^{2,ph}(\mathbf{q}, \omega) [1 - \cos \theta(\mathbf{q})], \qquad (A44) $$

where $g_{Gr}^{2,ph}$ is defined by Eq. 51. To have more manageable expressions, we now approximate the phonon frequency to be dispersionless, i.e. we assume that $g_{Gr}^{2,ph}$ is peaked at $\omega = \pm \omega_{ph}$ and approximate

$$ g_{Gr}^{2,ph}(\mathbf{q}, \omega) \sim g_{Gr}^{2,ph}(q) \delta(\omega \pm \omega_{ph}). \qquad (A45) $$

We remind that the above coupling is of long-range nature, and it is in very good approximation isotropic, as we used in the right hand side of the equation. Given this coupling, if we consider the special case $\epsilon = \epsilon_F$ the two energy conserving Dirac delta functions of Eq. A44 (for positive and negative frequencies, i.e. emission of absorption) imply that $\mathbf{k} + \mathbf{q}$ lies on the isoenergetic lines $|\epsilon_F \pm \hbar\omega_{ph}|$. As presented in Fig. 15 for an absorption process, $q$ and $\theta$ are now related geometrically, i.e. we can write $q_\theta$. Since for a given $\theta$ we have two different values of $q$ for emission and absorption, we will call them $q_\theta^\pm$ respectively. They respect the following relations

$$ q_{\theta=0}^{\pm} \le q_{\theta}^{\pm} \le q_{\theta=\pi}^{\pm}, \qquad (A46) $$

$$ \frac{q_{\theta=0}^{\pm}}{k_F} = \frac{| |\epsilon_F \pm \hbar\omega_{ph}| - |\epsilon_F| |}{|\epsilon_F|}, \qquad (A47) $$

$$ \frac{q_{\theta=\pi}^{\pm}}{k_F} = \frac{| |\epsilon_F| + |\epsilon_F \pm \hbar\omega_{ph}| |}{|\epsilon_F|}, \qquad (A48) $$

Now, the sum over $\mathbf{q}$ can then be transformed in just an angular integration, keeping in mind that the volume element $qd\mathbf{q}$ will bring a factor $|\epsilon_F \pm \hbar\omega_{ph}|$. The Bloch overlaps of graphene can be written in the form given below Eq. 26. With algebraic rewriting of the statistical occupation factors, at last we obtain Eq. 54, if we choose $\omega_{ph} = \omega_{BN}$.

# Appendix B: Computational details

We here detail the numerical computation of the ingredients of the VED method, and compare our results to previous approaches in Fig. 16.

## a. DFT calculations

DFT and DFPT calculations of ground states, dielectric responses and phonons are carried out with the QUANTUM ESPRESSO package$^{125,126}$, with 2d boundary conditions$^{127}$, and optimized norm-conserving Vanderbilt pseudopotentials$^{128}$ from the pseudo-DOJO library$^{129}$.

## b. Single-layer static electronic response from DFPT

The parametrization of single layer static electronic responses is detailed in$^{19}$. In the framework of DFPT, the layer is perturbed with a monochromatic in-plane periodic function of the form $V^i(q, z) = V_0\phi^i$. The density response $\rho^i$ is extracted from DFPT, renormalized by $V_0$ and integrated to obtain $Q^i(q)$. Examples of plots of $Q^i(q)$ can be found in Fig. 9 of Ref.$^{19}$. The profile function is then defined as $\rho^i(q, z)/Q^i(q)$. A small $q$-dependence of the profile function exists in practice, but has a small impact on the observables.

## c. In-plane contribution to atomic response

We here compare the approximation of Eq. A27 to the one of Eq. 30. In Eq. A27 the potential is computed using a procedure similar to Ref.$^{60}$—an exact procedure may be developed following Ref.$^{55}$, but we leave it for further developments. As anticipated in the main text, this assumes that the LO potentials are generated by a charge density with the same profile as the electronic response $f_{BN}^0$. In practice that is not exactly true. We find that the difference of implementing Eq. A27 instead of Eq. 30 (which is expected to be a better approximation) is order of 5%.

## 1. Neglect of out-of-plane atomic contribution

An upper limit estimation of the atomic contribution to Eq. 24 in the long wavelength limit is given by

$$ Q_{ph}^{1}(q, \omega) \sim \frac{1}{AM^{*}(\omega + i\eta_{ZO})^{2} - \omega_{ZO}^{2}}, \qquad (B1) $$

where $M^*$ is the reduced mass, we have neglected the wavevector dependence of the phonon frequency $\omega_{ZO} \sim 800\text{cm}^{-1}$, and $|Z^\perp| = 0.245$ is the absolute value of the Born effective charges of BN as computed from DFPT calculations performed at $q = 0$. With respect to the expression for the in-plane atomic response Eq. 32, in Eq. B1 we have approximated the out-of-plane dielectric function to unity, thus expecting an overestimation of $Q_{ph}^1(q, \omega)$, in line with what discussed in the previous Section in Fig. 16. The smearing is taken as at an unphysical value of $\eta_{ZO} \sim 80\text{cm}^{-1}$, for representation aims. In Fig. 17 we plot the real part of the ratio $Q_{ph}^1(q, \omega)/|Q_{el}^1(q)|$.

As seen from the plot, Eq. B1 has a 5% influence on the total value of $Q^1$ at $\omega = \omega_{LO}$ (vertical red line), meaning Eq. 24 is a good approximation for the BN single layer. Nevertheless, contributions may add up when adding several layers in a vdWH, i.e. going towards the hexagonal bulk limit. Even in this case though, the atomic contri-
<figure><img src="image_18.png" /><figcaption>FIG. 16: (Top) Comparison of $\tilde{g}^2$ in single layer BN, obtained via the VED model and the QEH model of Ref.²⁹, reducing to the LO phonon coupling. The discrepancy is attributed to the different evaluation of $-\text{Im}[\int_0^\infty Q^0(q, \omega)d\omega]$ (Bottom). Notice that the QEH error is up to 38% at the largest value of $q$ with respect to DFPT values, that are on top of the results of this work. There main source of error is the use, in the QEH method, of the strictly 2d Coulomb kernel $\frac{2\pi}{q}$ rather than the more appropriate $v_{\text{BN,BN}}^{00}$ to deduce the interacting response function $Q^0(q, \omega)$ from the irreducible response functions of electrons and phonons. Correcting this error, the QEH methodology is very near to our results.</figcaption></figure>

bution to $Q_{\text{ph}}^1$ is negligible. In fact, we can write

$$ \epsilon_0^{zz} = \epsilon_\infty^{zz} + \frac{4\pi e^2 |Z^{zz}|^2}{V M^* \omega_{\text{ZO}}^2}, \quad (\text{B2}) $$

where $|Z^{zz}| = 0.76$ and $\epsilon_\infty^{zz} = 2.729$ as from DFPT calculations at zone centers, concluding that $(\epsilon_0^{zz} - \epsilon_\infty^{zz})/\epsilon_\infty^{zz} \sim 7\%$. Even in this case, the effect of the ZO degrees of freedom is small on the total value of the dielectric function. Since the ZO phonon frequency is roughly equidistant from $\omega = 0$ and $\omega_{\text{LO}}$, the approximation of Eq. 24 remains good for the determination of LO dispersions even when stacking multiple BN layers.

<figure><img src="image_19.png" /><figcaption>FIG. 17: $Re(Q_{\text{ph}}^1(q, \omega))/|Q_{\text{el}}^1(q)|$ as a function of $\omega$, for the wavevectors relevant for this work. For the ZO mode at $\sim 800\text{cm}^{-1}$ (black line), $\eta_{\text{ZO}} \sim 80\text{cm}^{-1}$ is employed for representation purposes, the values at $\omega = \omega_{\text{LO}}$ (red line) not being influenced by this choice. We notice that the estimated upper limit of the ratio between the atomic and electronic contribution reaches a small value of $\sim 5\%$ at $\omega = \omega_{\text{LO}}$.</figcaption></figure>

## Appendix C: Comparison of screened electron-phonon coupling in VED and DFPT

When considering only the static limit of electronic responses, the VED method reproduces electron-phonon interactions obtained in DFPT calculations in the Born-Oppenheimer approximation. The $g$ coupling is extracted from DFPT for a given mode $\mu$ on layer k as

$$ g_{\mu,k}^{\text{DFPT}}(\mathbf{q}) = \left| \int dz f_k(z) V_{\mu}(\mathbf{q}, z) \right|, \quad (\text{C1}) $$

where $\bar{V}_\mu(\mathbf{q}, z) = \sum_s \frac{e_{\mu q}^s}{\sqrt{M_s}} \cdot \frac{\partial \bar{V}_{\text{KS}}(\mathbf{q}, z)}{\partial \mathbf{u}^q}$. In the above expression we can take the modulus because the coupling is always positive. From a practical point of view, it is implemented to remove the phase uncertainty of the phonon polarization vectors coming out from first-principles calculations. Notice that in DFPT the perturbation is explicitly inserted in the system via a displacement of the atoms, therefore we cannot choose the layer on which the external perturbation is applied. This is the reason why Eq. C1 presents only one layer index, that is treated as the first index (probe index) of $g_{kl}^2$.

In this Section, we consider a simple system made of BN and graphene. Ground state DFT calculations are performed with a non-uniform grid of momenta to properly sample the small Fermi surface of graphene and be able to use an electronic smearing corresponding to the Fermi-Dirac distribution at room temperature. The non-uniform grid varies from 12 × 12 to the equivalent of 196 × 196 around the Fermi surface. Graphene is doped with an extra density of electrons $n = 10^{13} \text{cm}^{-2}$, corresponding to a Fermi level $\sim 0.3 \text{ eV}$. Phonons are simulated at 5 momenta in the $\Gamma-M$ direction. We then select
the electron-phonon coupling coming from the (only) LO
mode of BN.

To compute the comparable quantity in the VED
method, we use the static limit of the electronic response
of graphene's electrons, computed directly in the sin-
gle layer graphene's system with the same parameters
as the BN/Graphene system. We then compute $\bar{g}_{\text{Gr}}^2(q)$
and $\bar{g}_{\text{BN}}^2(q)$ from Eq. 48.

Fig. 18 shows a good agreement between direct DFPT
and the VED method for the main quantity of interest in
this work $\bar{g}_{\text{Gr}}^2(q)$. Along with Fig. 9, this confirms that
the VED method reproduces DFPT in the static limit
of electronic responses even when metallic screening is
present. A larger error is observed for $\bar{g}_{\text{BN}}^2(q)$ in this case,
with respect to Fig. 9. This may be explained by the
fact that the total effective potential, i.e. the potential
generated by LO phonons and screened by electrons, dis-
plays much sharper variations in the out-of-plane direc-
tion due to the strong metallic screening from graphene.
This would enhance errors related to the approximation
of atomic profile functions with the electronic ones in the

VED model.

<figure><img src="image_20.png" /><figcaption>FIG. 18: Comparison of electron-phonon coupling matrix elements in VED ($\bar{g}_{\text{Gr}}^2$ and $\bar{g}_{\text{BN}}^2$) and DFPT results in BN/Gr. In the VED method, the static limit of graphene's electronic response is used, so that the coupling is only of electron-phonon type.</figcaption></figure>

1 A. Castellanos-Gomez, X. Duan, Z. Fei, H. R. Gutierrez, Y. Huang, X. Huang, J. Quereda, Q. Qian, E. Sutter, and P. Sutter, Nature Reviews Methods Primers **2**, 58 (2022).

2 K. S. Novoselov, A. Mishchenko, A. Carvalho, and A. H. Castro Neto, Science **353**, aaac9439 (2016), arxiv:1411.1235v1.

3 Y. Liu, N. O. Weiss, X. Duan, H.-C. Cheng, Y. Huang, and X. Duan, Nature Reviews Materials **1**, 16042 (2016).

4 L. V. Keldysh, A. A. Maradudin, and D. A. Kirzhnitz, *The Dielectric Function of Condensed Systems* (Elsevier Science, Amsterdam, 2014).

5 P. Nozières and D. Pines, *Theory Of Quantum Liquids: Superfluid Bose Liquids* (Westview Press, New York, 1999).

6 D. N. Basov, M. M. Fogler, and F. J. García De Abajo, Science **354**, aag1992 (2016).

7 T. Low, A. Chaves, J. D. Caldwell, A. Kumar, N. X. Fang, P. Avouris, T. F. Heinz, F. Guinea, L. Martin-Moreno, and F. Koppens, Nature Materials **16**, 182 (2017).

8 N. Rivera, T. Christensen, and P. Narang, Nano Letters **19**, 2653 (2019).

9 J. Chen, M. Badioli, P. Alonso-González, S. Thongratanasiri, F. Huth, J. Osmond, M. Spasenović, A. Centeno, A. Pesquera, P. Godignon, A. Zurutuza Elorza, N. Camara, F. J. G. De Abajo, R. Hillenbrand, and F. H. L. Koppens, Nature **487**, 77 (2012).

10 Z. Fei, A. S. Rodin, G. O. Andreev, W. Bao, A. S. McLeod, M. Wagner, L. M. Zhang, Z. Zhao, M. Thiemens, G. Dominguez, M. M. Fogler, A. H. C. Neto, C. N. Lau, F. Keilmann, and D. N. Basov, Nature **487**, 82 (2012).

11 H. Yan, T. Low, W. Zhu, Y. Wu, M. Freitag, X. Li, F. Guinea, P. Avouris, and F. Xia, Nature Photonics **7**, 394 (2013).

12 V. W. Brar, M. S. Jang, M. Sherrott, S. Kim, J. J. Lopez, L. B. Kim, M. Choi, and H. Atwater, Nano Letters **14**, 3876 (2014).

13 A. Woessner, M. B. Lundeberg, Y. Gao, A. Principi, P. Alonso-González, M. Carrega, K. Watanabe, T. Taniguchi, G. Vignale, M. Polini, J. Hone, R. Hillenbrand, and F. H. L. Koppens, Nature Materials **14**, 421 (2015).

14 E. Yoxall, M. Schnell, A. Y. Nikitin, O. Txoperena, A. Woessner, M. B. Lundeberg, F. Casanova, L. E. Hueso, F. H. L. Koppens, and R. Hillenbrand, Nature Photonics **9**, 674 (2015).

15 N. Li, X. Guo, X. Yang, R. Qi, T. Qiao, Y. Li, R. Shi, Y. Li, K. Liu, Z. Xu, L. Liu, F. J. García De Abajo, Q. Dai, E.-G. Wang, and P. Gao, Nature Materials **20**, 43 (2021).

16 M. Chen, Y. Zhong, E. Harris, J. Li, Z. Zheng, H. Chen, J.-S. Wu, P. Jarillo-Herrero, Q. Ma, J. H. Edgar, X. Lin, and S. Dai, Nature Communications **14**, 4782 (2023).

17 K. J. Tielrooij, N. C. Hesp, A. Principi, M. B. Lundeberg, E. A. Pogna, L. Banszerus, Z. Mics, M. Mascotte, P. Schmidt, D. Davydovskaya, D. G. Purdie, I. Goykhman, G. Soavi, A. Lombardo, K. Watanabe, T. Taniguchi, M. Bonn, D. Turchinovich, C. Stampfer, A. C. Ferrari, G. Cerullo, M. Polini, and F. H. Koppens, Nature Nanotechnology **13**, 41 (2018), arxiv:1702.03766.

18 W. Yang, S. Berthou, X. Lu, Q. Wilmart, A. Denis, M. Rosticher, T. Taniguchi, K. Watanabe, G. Fève, J.-M. Berroir, G. Zhang, C. Voisin, E. Baudin, and B. Plaçais, Nature Nanotechnology **13**, 47 (2018).

19 T. Sohier, M. Gibertini, and M. J. Verstraete, Physical Review Materials **5**, 024004 (2021), arxiv:2011.04961.

20 S. Dai, Z. Fei, Q. Ma, A. S. Rodin, M. Wagner, A. S. McLeod, M. K. Liu, W. Gannett, W. Regan, K. Watanabe, T. Taniguchi, M. Thiemens, G. Dominguez, A. H. C. Neto, A. Zettl, F. Keilmann, P. Jarillo-Herrero, M. M. Fogler, and D. N. Basov, Science **343**, 1125 (2014).

21 R. Fuchs and K. L. Kliewer, Phys. Rev. **140**, A2076 (1965).
22 E. N. Economou, Phys. Rev. **182**, 539 (1969).

23 M. E. Kim, A. Das, and S. D. Senturia, Physical Review B **18**, 6890 (1978).

24 B. A. Sanborn, Physical Review B **51**, 14256 (1995).

25 S. Fratini and F. Guinea, Physical Review B **77**, 195415 (2008), arxiv:0711.1303v3.

26 Z.-Y. Ong and M. V. Fischetti, Physical Review B **86**, 165422 (2012).

27 A. Hauber and S. Fahy, Physical Review B **95**, 045210 (2017).

28 K. Andersen, S. Latini, and K. S. Thygesen, Nano Letters **15**, 4616 (2015), arxiv:1506.02463.

29 M. N. Gjerding, L. S. R. Cavalcante, A. Chaves, and K. S. Thygesen, The Journal of Physical Chemistry C **124**, 11609 (2020).

30 L. Sponza and F. Ducastelle, arXiv:2011.07811 [cond-mat] (2020), arxiv:2011.07811 [cond-mat].

31 Y. in't Veld, M. I. Katsnelson, A. J. Millis, and M. Rösner, 2D Materials **10**, 045031 (2023).

32 G.-Q. Hai, N. Studart, and G. E. Marques, Physical Review B **55**, 1554 (1997).

33 S. Poncé, F. Macheda, E. R. Margine, N. Marzari, N. Bonini, and F. Giustino, Physical Review Research **3**, 043022 (2021).

34 M. V. Fischetti, Physical Review B **44**, 5527 (1991).

35 V. A. Jhalani, J.-J. Zhou, J. Park, C. E. Dreyer, and M. Bernardi, Physical Review Letters **125**, 136602 (2020), arxiv:2002.08351.

36 F. Macheda and N. Bonini, Physical Review B **98**, 201201 (2018).

37 G. Brunin, H. P. C. Miranda, M. Giantomassi, M. Royo, M. Stengel, M. J. Verstraete, X. Gonze, G.-M. Rignanese, and G. Hautier, Physical Review Letters **125**, 136601 (2020), arxiv:2002.00628.

38 F. Macheda, S. Poncé, F. Giustino, and N. Bonini, Nano Letters **20**, 8861 (2020).

39 S. Gopalan, M. L. Van De Put, G. Gaddemane, and M. V. Fischetti, Physical Review Applied **18**, 054062 (2022).

40 S. Piscanec, M. Lazzeri, F. Mauri, A. C. Ferrari, and J. Robertson, Physical Review Letters **93**, 185503 (2004).

41 T. Venanzi, L. Graziotto, F. Macheda, S. Sotgiu, T. Ouaj, E. Stellino, C. Fasolato, P. Postorino, V. Mišekis, M. Metzelaars, P. Kögerler, B. Beschoten, C. Coletti, S. Roddaro, M. Calandra, M. Ortolani, C. Stampfer, F. Mauri, and L. Baldassarre, Physical Review Letters **130**, 256901 (2023).

42 A. Grüneis, J. Serrano, A. Bosak, M. Lazzeri, S. L. Molodtsov, L. Wirtz, C. Attaccalite, M. Krisch, A. Rubio, F. Mauri, and T. Pichler, Physical Review B - Condensed Matter and Materials Physics (2009), 10.1103/Phys-RevB.80.085423, arxiv:0904.3205.

43 R. Kim, V. Perebeinos, and P. Avouris, Physical Review B **84**, 075449 (2011).

44 G. Marini and M. Calandra, Physical Review B **104**, 144103 (2021).

45 M. Harb, R. Ernstorfer, T. Dartigalongue, C. T. Hebeisen, R. E. Jordan, and R. J. D. Miller, The Journal of Physical Chemistry B **110**, 25308 (2006).

46 M. Bernardi, D. Vigil-Fowler, J. Lischner, J. B. Neaton, and S. G. Louie, Physical Review Letters **112**, 257402 (2014).

47 A. C. Betz, S. H. Jhang, E. Pallecchi, R. Ferreira, G. Fève, J.-M. Berroir, and B. Plaçais, Nature Physics **9**, 109 (2013).

48 J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Physical Review **108**, 1175 (1957).

49 J. R. Schrieffer, *Theory of Superconductivity*, rev. printing, 4. print ed., Advanced Book Classics (Westview Press, Bolder, Colo, 1999).

50 F. Marsiglio, Annals of Physics **417**, 168102 (2020).

51 S. Dai, Q. Ma, M. K. Liu, T. Andersen, Z. Fei, M. D. Goldflam, M. Wagner, K. Watanabe, T. Taniguchi, M. Thiemens, F. Keilmann, G. C. Janssen, S. E. Zhu, P. Jarillo-Herrero, M. M. Fogler, and D. N. Basov, Nature Nanotechnology **10**, 682 (2015), arxiv:1501.06956.

52 J. Sonntag, K. Watanabe, T. Taniguchi, B. Beschoten, and C. Stampfer, Physical Review B **107**, 075420 (2023).

53 M. Yankowitz, Q. Ma, P. Jarillo-Herrero, and B. J. LeRoy, Nature Reviews Physics **1**, 112 (2019).

54 The code will be made public on gitlab after acceptance.

55 M. Royo and M. Stengel, Physical Review X **11**, 041027 (2021).

56 We can now clarify the meaning of the dual basis set in mathematical terms: at each fixed $z$, $\chi_{1L}$ is an element of the dual space of the potentials. The separability of $\chi_{1L}(q, \omega, z, z') = (\chi_{1L}(q, \omega, z', z))^*$ imply that the basis set of the densities is the same as of the dual space, justifying our nomenclature choice.

57 G. Grosso and G. Pastori Parravicini, Solid State Physics (Academic Press, San Diego, 2000).

58 T. Sohier, M. Calandra, C.-H. C.-H. Park, N. Bonini, N. Marzari, and F. Mauri, Physical Review B **90**, 125414 (2014), arxiv:1407.0830.

59 T. Sohier, M. Calandra, and F. Mauri, Phys. Rev. B **91**, 165428 (2015).

60 F. Macheda, T. Sohier, P. Barone, and F. Mauri, Physical Review B **107**, 094308 (2023).

61 V. U. Nazarov, New Journal of Physics **17**, 073018 (2015).

62 K. Sturm, Zeitschrift für Naturforschung A **48**, 233 (1993).

63 R. Senga, K. Suenaga, P. Barone, S. Morishita, F. Mauri, and T. Pichler, Nature **573**, 247 (2019).

64 G. D. Mahan, *Many-Particle Physics*, 2nd ed., Physics of Solids and Liquids (Plenum Press, New York, 1990).

65 A. N. Grigorenko, M. Polini, and K. S. Novoselov, Nature Photonics **6**, 749 (2012).

66 A. Politano and G. Chiarello, Nanoscale **6**, 10927 (2014).

67 E. H. Hwang and S. Das Sarma, Physical Review B **75**, 205418 (2007).

68 M. Polini, R. Asgari, G. Borghi, Y. Barlas, T. Pereg-Barnea, and A. H. MacDonald, Physical Review B **77**, 081411 (2008).

69 a. Hill, S. a. Mikhailov, and K. Ziegler, EPL (Europhysics Letters) **87**, 27005 (2009).

70 A. Principi, G. Vignale, M. Carrega, and M. Polini, Physical Review B **88**, 195405 (2013).

71 A. Principi, M. Carrega, M. B. Lundeberg, A. Woessner, F. H. L. Koppens, G. Vignale, and M. Polini, Physical Review B **90**, 165408 (2014).

72 K. Zhang, Y. Feng, F. Wang, Z. Yang, and J. Wang, Journal of Materials Chemistry C **5**, 11992 (2017).

73 J. D. Caldwell, I. Aharonovich, G. Cassabois, J. H. Edgar, B. Gil, and D. N. Basov, Nature Reviews Materials **4**, 552 (2019).

74 Z. Shi, H. A. Bechtel, S. Berweger, Y. Sun, B. Zeng, C. Jin, H. Chang, M. C. Martin, M. B. Raschke, and F. Wang, ACS Photonics **2**, 790 (2015).

75 S. Dai, W. Fang, N. Rivera, Y. Stehle, B.-Y. Jiang,
J. Shen, R. Y. Tay, C. J. Ciccarino, Q. Ma, D. Rodan-Legrain, P. Jarillo-Herrero, E. H. T. Teo, M. M. Fogler, P. Narang, J. Kong, and D. N. Basov, Advanced Materials **31**, 1806603 (2019).

76 D. Sánchez-Portal and E. Hernández, Physical Review B **66**, 235415 (2002).

77 K. H. Michel and B. Verberck, Physical Review B **80**, 224301 (2009).

78 T. Sohier, M. Gibertini, M. Calandra, F. Mauri, and N. Marzari, Nano Letters **17**, 3758 (2017), arxiv:1612.07191.

79 J. Li, L. Wang, Z. Tao, W. Zhong, S. Xue, G. Miao, W. Wang, J. Guo, and X. Zhu, "Observation of the Breakdown of Optical Phonon Splitting in a Two-dimensional Polar Monolayer," (2023), arxiv:2306.17458 [cond-mat].

80 K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Science **306**, 666 (2004).

81 F. Schwierz, Nature Nanotechnology **5**, 487 (2010).

82 A. C. Ferrari, F. Bonaccorso, V. Fal'ko, K. S. Novoselov, S. Roche, P. Bøggild, S. Borini, F. H. L. Koppens, V. Palermo, N. Pugno, J. A. Garrido, R. Sordan, A. Bianco, L. Ballerini, M. Prato, E. Lidorikis, J. Kivioja, C. Marinelli, T. Ryhänen, A. Morpurgo, J. N. Coleman, V. Nicolosi, L. Colombo, A. Fert, M. Garcia-Hernandez, A. Bachtold, G. F. Schneider, F. Guinea, C. Dekker, M. Barbone, Z. Sun, C. Galiotis, A. N. Grigorenko, G. Konstantatos, A. Kis, M. Katsnelson, L. Vandersypen, A. Loiseau, V. Morandi, D. Neumaier, E. Treossi, V. Pellegrini, M. Polini, A. Tredicucci, G. M. Williams, B. Hee Hong, J.-H. Ahn, J. Min Kim, H. Zirath, B. J. Van Wees, H. Van Der Zant, L. Occhipinti, A. Di Matteo, I. A. Kinloch, T. Seyller, E. Quesnel, X. Feng, K. Teo, N. Rupesinghe, P. Hakonen, S. R. T. Neil, Q. Tannock, T. Löfwander, and J. Kinaret, Nanoscale **7**, 4598 (2015).

83 G. Giuliani and G. Vignale, *Quantum Theory of the Electron Liquid* (Cambridge University Press, 2005).

84 S. Q. Wang and G. D. Mahan, Phys. Rev. B **6**, 4517 (1972).

85 C. Attaccalite, L. Wirtz, M. Lazzeri, F. Mauri, and A. Rubio, Nano letters **10**, 1172 (2010), 20222744.

86 L. Graziotto, F. Macheda, T. Sohier, M. Calandra, and F. Mauri, Phys. Rev. B **109**, 075420 (2024).

87 L. Graziotto, F. Macheda, T. Venanzi, G. Marchese, S. Sotgiu, T. Ouaj, E. Stellino, C. Fasolato, P. Postorino, M. Metzelaars, P. Kögerler, B. Beschoten, M. Calandra, M. Ortolani, C. Stampfer, F. Mauri, and L. Baldassarre, Nano Letters **24**, 1867 (2024).

88 L. Banszerus, T. Sohier, A. Epping, F. Winkler, F. Lيبisch, F. Haupt, K. Watanabe, T. Taniguchi, K. Müller-Caspary, N. Marzari, F. Mauri, B. Beschoten, and C. Stampfer, ArXiv:1909.09523 (2019), arxiv:1909.09523.

89 L. Wang, I. Meric, P. Y. Huang, Q. Gao, Y. Gao, H. Tran, T. Taniguchi, K. Watanabe, L. M. Campos, D. a Muller, J. Guo, P. Kim, J. Hone, K. L. Shepard, and C. R. Dean, Science **342**, 614 (2013), 24179223.

90 T. Sohier, M. Gibertini, and N. Marzari, 2D Materials **8**, 015025 (2020), arxiv:2007.16110.

91 K. Kaasbjerg, K. S. Thygesen, and K. W. Jacobsen, Physical Review B Condensed Matter and Materials Physics **85**, 115317 (2012), arxiv:1201.5284.

92 S. Poncé, M. Royo, M. Gibertini, N. Marzari, and M. Stengel, Phys. Rev. Lett. **130**, 166301 (2023).

93 C. Zhang, R. Wang, H. Mishra, and Y. Liu, Phys. Rev. Lett. **130**, 087001 (2023).

94 L. Cheng, C. Zhang, and Y. Liu, Computational Materials Science **194**, 110468 (2021).

95 X. Duan, T. Wang, Z. Fu, L. Liu, and J.-Y. Yang, Phys. Chem. Chem. Phys. **25**, 10175 (2023).

96 J. Ma, D. Xu, R. Hu, and X. Luo, Journal of Applied Physics **128**, 035107 (2020), https://pubs.aip.org/aip/jap/article/pdf/doi/10.1063/5.0015102/15251722/035107_1_online.pdf.

97 L. Cheng and Y. Liu, Journal of the American Chemical Society **140**, 17895 (2018).

98 S. Song, Y. Sun, S. Liu, J.-H. Yang, and X.-G. Gong, Phys. Rev. B **107**, 155101 (2023).

99 W. Li, S. Poncé, and F. Giustino, Nano Letters **19**, 1774 (2019).

100 N. Ma and D. Jena, Physical Review X **4**, 011043 (2014), arxiv:1310.7157.

101 T. Sohier, M. Gibertini, D. Campi, G. Pizzi, and N. Marzari, Nano Letters **19**, 3723 (2019).

102 T. Sohier, D. Campi, N. Marzari, and M. Gibertini, Phys. Rev. Mater. **2**, 114010 (2018).

103 C.-H. Park, N. Bonini, T. Sohier, G. Samsonidze, B. Kozinsky, M. Calandra, F. Mauri, and N. Marzari, Nano Letters **14**, 1113 (2014).

104 S.-C. Lee and I. Galbraith, Physical Review B **59**, 15796 (1999).

105 T. Sohier, E. Ponomarev, M. Gibertini, H. Berger, N. Marzari, N. Ubrig, and A. F. Morpurgo, Physical Review X **9**, 031019 (2019), arxiv:1901.08012.

106 R. De L. Kronig, Journal of the Optical Society of America **12**, 547 (1926).

107 F. Macheda, P. Barone, and F. Mauri, Physical Review Letters **129**, 185902 (2022).

108 L. Hedin and S. Lundqvist, in *Solid State Physics*, Vol. 23 (Elsevier, 1970) pp. 1-181.

109 O. V. Dolgov, D. A. Kirzhnits, and E. G. Maksimov, Reviews of Modern Physics **53**, 81 (1981).

110 A. E. Karakozov and E. G. Maksimov, Sov. Phys. - JETP (Engl. Transl.); (United States) **47**:2 (1978).

111 L. Binci, P. Barone, and F. Mauri, Physical Review B **103**, 134304 (2021).

112 G. Marchese, F. Macheda, L. Binci, M. Calandra, P. Barone, and F. Mauri, Nature Physics (2023), 10.1038/s41567-023-02203-3.

113 C. E. Dreyer, S. Coh, and M. Stengel, Phys. Rev. Lett. **128**, 095901 (2022).

114 C.-Y. Wang, S. Sharma, E. K. U. Gross, and J. K. Dewhurst, Physical Review B **106**, L180303 (2022).

115 H. Fröhlich, Advances in Physics **3**, 325 (1954), https://doi.org/10.1080/00018735400101213.

116 J. Sjakste, N. Vast, M. Calandra, and F. Mauri, Phys. Rev. B **92**, 054307 (2015).

117 C. Verdi and F. Giustino, Phys. Rev. Lett. **115**, 176401 (2015).

118 P. Vogl, Phys. Rev. B **13**, 694 (1976).

119 R. M. Pick, M. H. Cohen, and R. M. Martin, Physical Review B **1**, 910 (1970).

120 M. Royo and M. Stengel, Phys. Rev. X **11**, 041027 (2021).

121 S. Poncé, M. Royo, M. Stengel, N. Marzari, and M. Gibertini, Phys. Rev. B **107**, 155424 (2023).

122 F. Giustino, Reviews of Modern Physics **89**, 015003 (2017), arxiv:1603.06965.

123 C. Verdi, F. Caruso, and F. Giustino, Nature Commun
cations **8**, 15769 (2017).

<footnote n="124">¹²⁴ F. Caruso and F. Giustino, Phys. Rev. B **94**, 115208 (2016).</footnote>

<footnote n="125">¹²⁵ P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. Dal Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, R. M. Wentzcovitch, A. D. Corso, G. Fratesi, S. D. Gironcoli, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, Journal of Physics: Condensed Matter **21**, 395502 (2009), arxiv:0906.2569.</footnote>

<footnote n="126">¹²⁶ P. Giannozzi, O. Andreussi, T. Brumme, O. Bunau, M. Buongiorno Nardelli, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, M. Cococcioni, N. Colonna, I. Carnimeo, A. Dal Corso, S. De Gironcoli, P. Delugas, R. A. Distasio, A. Ferretti, A. Floris, G. Fratesi, G. Fugallo, R. Gebauer, U. Gerstmann, F. Giustino, T. Gorni, J. Jia, M. Kawamura, H. Y. Ko, A. Kokalj, E. Küçükbenli, M. Lazzeri, M. Marsili, N. Marzari, F. Mauri, N. L. Nguyen, H. V. Nguyen, A. Otero-DeLa-Roza, L. Paulatto, S. Poncé, D. Rocca, R. Sabatini, B. Santra, M. Schlipf, A. P. Seitsonen, A. Smogunov, I. Timrov, T. Thonhauser, P. Umari, N. Vast, X. Wu, and S. Baroni, Journal of Physics Condensed Matter **29**, 465901 (2017), arxiv:1709.10010.</footnote>

<footnote n="127">¹²⁷ T. Sohier, M. Calandra, and F. Mauri, Physical Review B **96**, 075448 (2017), arxiv:1705.04973.</footnote>

<footnote n="128">¹²⁸ D. R. Hamann, Physical Review B **88**, 085117 (2013).</footnote>

<footnote n="129">¹²⁹ M. J. van Setten, M. Giantomassi, E. Bousquet, M. J. Verstraete, D. R. Hamann, X. Gonze, and G. M. Rignanese, Computer Physics Communications **226**, 39 (2018), arxiv:1710.10138.</footnote>