Dielectric screening by 2D substrates

Two-dimensional (2D) materials, being atomically thin and flexible, are prime candidates as active components in next generation nanoscale devices. Many properties of 2D materials, such as charge density wave instabilities [1], the presence of strongly bound excitons [2], as well as tunable interlayer excitons [3–5], arise from the reduced electronic screening in two dimensions. In fact, electronic screening in 2D systems is non-trivial. Unlike in 3D, where screening can be described by a macroscopic dielectric constant, dielectric screening in 2D is highly nonlocal [6–8]. While electronic screening within 2D materials has been studied extensively, the question still remains of how 2D substrates screen charge perturbations or electronic excitations adjacent to them. This question becomes particularly pertinent with recent interest in incorporating 2D materials as part of larger heterostructures.

One of the major breakthroughs in graphene electronics was the discovery that placing graphene on hexagonal boron nitride (BN) substrates on top of SiO₂ results in superior electron and hole mobilities as compared with graphene placed directly on SiO₂ [9]. Subsequent studies have shown that graphene placed on BN (and BN itself) has significantly lower charge density fluctuations than graphene on SiO₂ [10, 11]. However, it is unclear if the effect results from the migration of charge impurities away from the interface due to the presence of BN, or if it is due simply to a greater spatial separation between graphene and the charge impurities on SiO₂, or if BN itself plays a significant role in screening away the effect of the charge impurities. Electrostatic force microscopy (EFM) experiments on three common 2D materials (graphene, MoS₂ and BN) placed on SiO₂ suggest that electronic screening by atomically thin BN and MoS₂ has relatively weak layer number dependence compared to few layer graphene [12–14]. However, the interpretation of these experimental results is complicated by the presence of impurities and charge transfer from the environment, with some of the thickness-dependent trends being explained by the presence of water molecules [12]. Furthermore, the estimated density of defects in the experiment on BN was lower than that on MoS₂ and graphene (BN: 10⁸ cm⁻²; MoS₂, graphene: 10¹² cm⁻²) [12–14]. Taken together, the importance of screening by 2D materials and the questions raised in the experiments underscore the need for first principles studies on the thickness-dependent intrinsic screening properties of layered 2D materials.

Besides screening of unwanted charge impurities, electronic screening by substrates can also change the quasiparticle (QP) gaps of adjacent materials [15, 16]. This is particularly noticeable for small molecular adsorbates. To understand the effect of screening on the QP gap, consider the energy required to remove an electron from the highest occupied molecular orbital (HOMO) to vacuum. Electronic screening from the substrate will stabilize the resultant hole, moving the HOMO level up. Similarly, screening moves the lowest unoccupied molecular orbital (LUMO) level down. These effects are not present in typical density functional theory (DFT) calculations but are captured by many-body GW calculations [17]. The high cost of GW calculations has motivated simpler models to estimate screening from bulk substrates using a classical image charge model [15, 16]. However, screening by 2D substrates is more subtle, due to the difficulty in defining a macroscopic dielectric constant, and the atomically thin nature of the material. Given the growing interest in organic-2D material heterostructures [18–20], it is timely to explore methods to estimate the energy level alignment at organic-2D material interfaces to facilitate bottom-up design.

In this work, we employ first principles calculations to study the electronic screening of point charge perturbations adjacent to a 2D material. Specifically, we compute ab initio ${{V}_{scr}}(\boldsymbol{r},{\boldsymbol{r}^{\prime}})$, the screening potential at r due to an electron point charge perturbation at r' [21]. The latter represents point charges on an underlying SiO₂ substrate, on which a 2D material is placed, or represents the QP in an adsorbed molecule. We thus simulate ab initio the surface potentials of 2D materials in the presence of underlying point charges, and find that our results compare well with data from EFM experiments [12–14] as well as Kelvin probe [11] experiments. We also provide a reliable back-of-the-envelope way to estimate the HOMO-LUMO gaps for small adsorbed molecules on 2D and 3D substrates, from the band gap of the substrate and the gas phase gap of the molecule.

We define the screening potential to be ${{V}_{scr}}=W-V$, with V being the bare Coulomb potential and W the static screened Coulomb potential. In real space, the screening potential ${{V}_{scr}}(\boldsymbol{r},{\boldsymbol{r}^{\prime}})$, with r' representing the position of an extra electron point charge and r representing the position of the probe, is given by the formula [22]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{V}_{scr}}(\boldsymbol{r},{\boldsymbol{r}^{\prime}})=&\sum\limits_{\boldsymbol q, \boldsymbol G,{{\boldsymbol G}^{\prime}}}{{{e}^{i(\boldsymbol q+\boldsymbol G)\cdot \boldsymbol r}}\left[\varepsilon _{\boldsymbol G,{{\boldsymbol G}^{\prime}}}^{-1}(\boldsymbol q, \omega =0)-{{\delta}_{\boldsymbol G,{{\boldsymbol G}^{\prime}}}} \right]}\nonumber \\ & \times \frac{4\pi {{e}^{2}}}{\Omega |{\boldsymbol q}+{{\boldsymbol G}^{\prime}}{{|}^{2}}}{{e}^{-i(\boldsymbol q+{{\boldsymbol G}^{\prime}})\cdot {\boldsymbol{r}^{\prime}}}},\nonumber \end{align} \tag{ 1 }$$

where ${{\varepsilon}_{\boldsymbol G,{{\boldsymbol G}^{\prime}}}}(\boldsymbol q, \omega =0)$ is the static dielectric function computed using the random phase approximation (RPA). Details of the calculations are given in the Methods section.

Impurity screening by 2D substrates has recently been studied using EFM experiments [12–14], which effectively measure the difference in surface potentials between the top and bottom of the 2D slabs [14], in the presence of charged impurities on the underlying SiO₂ substrate. Since the difference in surface potentials in non-polar slabs is caused by the induced charge in response to the underlying charged impurity, we can compute this quantity for non-polar slabs of 2D materials using

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Delta {{V}_{scr}}={{V}_{scr}}({\boldsymbol{r}_{top}},{\boldsymbol{r}^{\prime}})-{{V}_{scr}}({\boldsymbol{r}_{bottom}},{\boldsymbol{r}^{\prime}}),\nonumber \end{align} \tag{ 2 }$$

where r' is taken to be ~3 Å beneath the bottom atomic plane (approximating the position of charge impurities on an SiO₂ substrate), and r_top and r_bottom are at the top and bottom atomic planes, respectively. We estimate $\Delta {{V}_{scr}}$ by taking r_top and r_bottom to be directly above r', and averaging over the in-plane coordinates of the perturbing charge, r'. In figure 1, we plot the quantity,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Delta V_{scr}^{dif\!\,\!f}=|\Delta {{V}_{scr}}-\Delta V_{scr}^{bulk}|,\nonumber \end{align} \tag{ 3 }$$

which is the magnitude of the difference between $\Delta {{V}_{scr}}$ computed for the thin film and that for the bulk, as a function of film thickness. This is the quantity that is plotted in experiment (its sign varies depending on the sign of the point charges and the direction of charge transfer), and in figure 1 we have extracted its value for BN [12], ${\rm MoS}_{2}$ [14], and graphene [13] for comparison with our calculated results. For $\Delta V_{scr}^{bulk}$, we use the value of $\Delta {{V}_{scr}}$ computed for the thickest film we have. We find that the qualitative agreement between the computed and experimental values is remarkably good, despite the fact that the values from experiment are influenced by charge transfer effects from the environment [23]. This good agreement implies that the qualitative thickness-dependent screening effects observed in experiment can be explained by intrinsic screening properties of the 2D materials.

Zoom In Zoom Out Reset image size

In [12], the experimental results were explained using a non-linear Thomas-Fermi model that took into account interlayer coupling. According to this model, the measured value of $\Delta V_{scr}^{dif\!f}$ for 1-layer BN would be an outlier, and this was attributed to the presence of much more charge transfer between an underlying dipole water film and the 1-layer BN [12]. However, our ab initio results indicate that the value of $\Delta V_{scr}^{dif\!\,\!f}$ for 1-layer BN fits well with the rest of the data for few-layer BN and the trend observed in experiment agrees well with the trend in our calculations (figure 1(c-d)). We therefore propose that the data for 1-layer BN is not an outlier, and suggest that a non-linear Thomas-Fermi model may not be sufficient for explaining all the non-trivial screening characteristics in 2D substrates.

The different dielectric screening from the environment and differing defect densities in the different experiments [12–14] make it difficult to deduce from experiment the relative intrinsic thickness-dependent screening properties of the various materials. Using our ab initio data, we see that $\Delta V_{scr}^{dif\!\,\!f}$ changes with thickness most quickly for few-layer graphene. However, the differences we predict are much less drastic than those observed in the experiments [12–14].

We address the quantitative difference between our ab initio results and the experimental data by including the effects of environmental screening using a dielectric constant $\varepsilon$ (figure 1). The value of $\varepsilon$ is deduced by fitting our ab initio results in (3) to the experimental data.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Delta V_{scr}^{model}(1L)&=\frac{1}{\varepsilon}\left(V_{scr}^{ab{{\rm \text -}}initio}(1L)-V_{scr}^{ab{{\rm \text -}}initio}(bulk) \right)\nonumber \\ & ={{V}^{\exp}}(1L)-{{V}^{\exp}}(bulk)\nonumber \end{align} \tag{ 4 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Delta V_{scr}^{model}({{n}_{\max}}L)&=\frac{1}{\varepsilon}\left(V_{scr}^{ab{{\rm \text -}}initio}({{n}_{\max}}L)-V_{scr}^{ab{{\rm \text -}}initio}(bulk) \right)\nonumber \\ & ={{V}^{\exp}}({{n}_{\max}}L)-{{V}^{\exp}}(bulk).\nonumber \end{align} \tag{ 5 }$$

The model results compare very well with the experimental data (figure 1). We note that our definition of $\Delta {{V}_{scr}}$ is likely to overestimate the difference in surface potentials between the top and bottom of the slabs, for typical defect densities, and therefore, our fitted values of $\varepsilon$ should be taken to be an upper bound to the actual value of $\varepsilon$ in the experiment. Nonetheless, we comment that the value of $\varepsilon$ found for the BN data (56.16) is consistent with the presence of water [24], which is believed to exist in the BN experiment [12]. On the other hand, the value of $\varepsilon$ found for the graphene data (2.32) agrees well with the dielectric constant estimated in [13].

It is well known that graphene devices placed atop BN substrates on SiO₂ have much higher mobility than graphene devices placed directly atop SiO₂ [9]. This has been partially attributed to the atomically flat nature of BN. At the same time, experiments also provide evidence of reduced charge potential fluctuations on top of BN, compared to SiO₂ [10, 11]. For example, in [11], the autocorrelation function, C(r), computed for the measured surface potential using Kelvin probe microscopy, was smaller in magnitude on top of BN compared to SiO₂, from which it was deduced that BN reduced charge potential fluctuations. It is unclear if this effect results from the migration of charge impurities away from the interface due to the presence of BN, or if it is due simply to a greater spatial separation between graphene and the charge impurities on SiO₂, or if BN itself plays a significant role in screening away the effect of the charge impurities. We directly address the role of BN in reducing charge potential fluctuations by computing the 2D autocorrelation function, C(r), of the screened and bare Coulomb potentials in different planes above a point charge perturbation, which is taken to be ~3 Å below BN, similar to what is expected for charge impurities on an underlying SiO₂ substrate. In order to simulate C(r) for large enough r, we have to obtain the value of $W(\boldsymbol{r},{\boldsymbol{r}^{\prime}})$ for large $\left|\boldsymbol r-{\boldsymbol{r}^{\prime}} \right|$, which implies that we need a very dense q-point mesh for computing ${{\varepsilon}_{\boldsymbol G,{{\boldsymbol G}^{\prime}}}}(\boldsymbol q,\omega =0)$. In this work, we have used a q-point mesh of 73  ×  73  ×  1 for BN when computing ${{\varepsilon}_{\boldsymbol G,{{\boldsymbol G}^{\prime}}}}(\boldsymbol q,\omega =0)$ (corresponding equivalently to a real-space supercell of 73  ×  73  ×  1 in the expansion of $W(\boldsymbol{r},{\boldsymbol{r}^{\prime}})$), and limited our study to monolayers. From figure 2(a), we can see that C(r) for the screened Coulomb potential W in a plane 1 Å above the BN surface is smaller, and drops off more slowly with r, than C(r) for the bare Coulomb potential V in a plane 1 Å above the impurity (below the BN layer). These results for C(r) are in qualitative agreement with the experimental results in [11], where the BN slab is thicker. Furthermore, both the rate at which C(r) drops off with r and the magnitude of C(r) are smaller for W than for V in the same plane above BN. These features imply that intrinsic screening from BN is important for reducing the potential fluctuations created by underlying charge impurities. In figure 2(b), we show C(r) for W and V above 1L BN, 1L MoS₂, and graphene, with the supercell sizes of the latter two chosen to closely match the area of the BN supercell. We see that, while C(r) for V above BN and graphene are the same (the small inconsistency in the plot is due to the slightly different areas of the cells), C(r) for W of graphene is both flatter and greatly reduced in magnitude compared to BN. This observation is consistent with the larger screening expected from semi-metallic graphene, compared to the large band gap insulator, BN. The smaller magnitude of C(r) for V above MoS₂ compared to BN, meanwhile, is indicative of the greater thickness of the MoS₂ layer. Nevertheless, C(r) for W above MoS₂ shows potential fluctuation reductions beyond those of the corresponding V curve, demonstrating that MoS₂ also effectively screens the charge impurity. Thus, we can conclude that the non-trivial dielectric screening is a general feature of 2D materials, and that the thickness of the screening substrate is important. However, the effective screening of impurities by BN in particular, together with the fact that BN itself is relatively free of dangling bonds and charge impurities, makes it an ideal substrate to improve charge mobility in graphene [9–11].

Zoom In Zoom Out Reset image size

Interestingly, the quantity ${{V}_{scr}}(\boldsymbol{r},{\boldsymbol{r}^{\prime}})$ is also useful for estimating the HOMO-LUMO gaps of small molecules adsorbed on substrates. In particular, the renormalization of molecular levels for molecules physisorbed on a substrate is given by the static polarization integral [15],

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{P}_{j}}=-\frac{1}{2}\iint{d {\boldsymbol{r}} d{\boldsymbol{r}^{\prime}}{{\phi}_{j}}(\boldsymbol{r})\phi _{j}^{*}({\boldsymbol{r}^{\prime}})\Delta W(\boldsymbol{r},{\boldsymbol{r}^{\prime}})}{{\phi}_{j}}({\boldsymbol{r}^{\prime}})\phi _{j}^{*}(\boldsymbol{r}),\nonumber \end{align} \tag{ 6 }$$

where ${{\phi}_{j}}$ represents the j th molecular orbital, and $\Delta W$ is the change in screened Coulomb potential in the molecule upon adsorption. Thus,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle E_{gap}^{{{V}_{scr}}}=E_{gap}^{gas-phase}-({{P}_{{\rm HOMO}}}+{{P}_{{\rm LUMO}}}).\nonumber \end{align} \tag{ 7 }$$

For small molecules, we may approximate $\Delta W(\boldsymbol{r},{\boldsymbol{r}^{\prime}})$ by $\Delta W({\boldsymbol r},{\boldsymbol r})\approx \Delta {{V}_{scr}}({\boldsymbol r},{\boldsymbol r})$ [15]. It has been shown that for benzene physisorbed on bulk metallic or semiconducting substrates, $\Delta {{V}_{scr}}({\boldsymbol r},{\boldsymbol r})$ can be approximated by a classical image potential model. Such an approximation assumes that $\Delta {{V}_{scr}}({\boldsymbol r},{\boldsymbol r})$ is given by the screening potential of the substrate alone, and gives good agreement with benchmark GW results for the HOMO-LUMO gap of benzene on various bulk substrates [15, 16]. For 2D materials, the classical image potential model cannot be directly applied because of the ambiguity in defining a macroscopic dielectric constant, as discussed above. By approximating $\Delta {{V}_{scr}}({\boldsymbol r},{\boldsymbol r})$ as the screening potential ${{V}_{scr}}({\boldsymbol r},{\boldsymbol r})$ of the 2D material substrate, with r being the center of the molecule (point charge model), we can estimate the HOMO-LUMO gap of benzene physisorbed on 2D materials without having to do a full GW calculation for the interface, using the following equation:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle E_{gap}^{{{V}_{scr}}}=E_{gap}^{gas-phase}+{{V}_{scr}}(\boldsymbol r,\boldsymbol r),\nonumber \end{align} \tag{ 8 }$$

with r in the center of the molecule. We use the above method to predict the HOMO-LUMO (π  −  π^*) gaps for single benzene molecules adsorbed on 1- and 2-layer BN, 1-layer MoS₂, and graphene. We find that the benzene molecules adsorb at a height of ~3.25 Å above the top atomic plane and use this adsorption height in the calculations. We note that in this regime, small amounts of hybridization between the benzene LUMO and the electronic states of the substrate are present for BN and MoS₂ (table S1). This is in contrast to previous literature utilizing equation (7), where there is no hybridization between benzene and the substrate, either because of the substrate material or the chosen adsorption heights [15, 16]. Nevertheless, the benzene molecules do not interact strongly with the substrates and it is instructive to see how equations (7) and (8) fare in these more realistic adsorption geometries. We compare the predictions of equations (7) and (8) with benchmark GW calculations of the same system in 3  ×  3 supercells (as well as a 4  ×  4 supercell for 1L BN) (figure 3). Notably, the 2D materials lead to significant renormalization of the benzene HOMO-LUMO gap from its gas phase value (10.6 eV). The agreement between $E_{gap}^{{{V}_{scr}}}$ and the GW values is reasonable and shows that ${{V}_{scr}}(\boldsymbol r,\boldsymbol r)$ captures most of the physics of screening by the 2D substrates, despite the small amount of hybridization present. We note that the effects of hybridization can be included in our estimate of the HOMO-LUMO gap using more involved methodologies [25, 26]. We can further improve our estimation in equation (8) by using Bader charges [27], centered at each atom of the molecule, instead of using a point charge model. These results are given in table 1.

Zoom In Zoom Out Reset image size

Table 1. Predicted π  −  π^* gaps for benzene on 2D substrates. The numbers are computed using the point charge model. Provided in square brackets are the numbers computed using the Bader charge model.

	GW π  −  π* gap of full interface system (eV)	$E_{gap}^{{{V}_{scr}}}$ (eV)
Benzene/1L BN (3  ×  3)	8.35	9.04 [9.21]
Benzene/1L BN (4  ×  4)	8.73	9.04 [9.18]
Benzene/2L BN (3  ×  3)	8.17	8.79 [9.08]
Benzene/1L MoS2 (3  ×  3)	7.66	8.14 [8.36]
Benzene/graphene (3  ×  3)	7.51	7.81 [7.95]

As ${{V}_{scr}}(\boldsymbol r,\boldsymbol r)$ evaluated 3.25 Å away from the substrate is a good proxy for the π  −  π^* gap renormalization in adsorbed benzene, we evaluate this quantity for ~20 2D and 3D substrates (the 3D substrates are represented by slabs that are at least 10 Å thick). Both 2D and 3D substrates are chosen to cover a wide range of DFT Kohn–Sham gaps, $E_{gap}^{substrate}$. Interestingly, both 2D and 3D substrates obey the same approximately linear relation between $E_{gap}^{substrate}$ and the screening response to a point charge 3.25 Å above the substrate (figure 4). We note that the value of ${{V}_{scr}}$ is rather sensitive to the choice of the surface plane, which, for non-planar surfaces such as CF, LiCl, TiO₂, and SiO₂, is not obvious. Here we simply choose the surface plane to be the average of the coordinates of the surface atoms, though a more rigorous approach merits investigation. The observed linear behavior is in stark contrast to the well-known fact that screening within the 2D material is much weaker than that within 3D materials [6, 28–30]. It also questions the validity of the common practice of computing an effective dielectric constant for a composite system by averaging dielectric constants of individual components [31, 32] —for an atomically thin material, the space occupied by vacuum would be much larger than that by the material, giving an effective dielectric constant close to 1 instead. In fact, using a classical electrostatics model, we can explain why 2D substrates can screen adjacent electronic excitations as well as 3D substrates with similar QP gaps/dielectric constants (see supplementary information (stacks.iop.org/TDM/6/035036/mmedia)). Our fully quantum mechanical calculations also show that the induced charge [21], ${{\rho}_{scr}}(\boldsymbol{r},{\boldsymbol{r}^{\prime}})$, which generates the screening response ${{V}_{scr}}(\boldsymbol{r},{\boldsymbol{r}^{\prime}})$, is primarily located at the surface (see supplementary information; figure 5). The localization of induced charge at the surface is in agreement with classical electrostatics, which tells us that the polarization-induced charges must lie only at the surface where there is a discontinuity in the dielectric environment. However, the quantum mechanical calculations point to a small distribution of the induced charge above and below the surface plane, with out-of-plane orbitals enabling the atomically thin 2D materials like BN to better screen the perturbing charge. We remark that the linear relation in figure 4 enables one to estimate the HOMO-LUMO gaps of small molecules physisorbed on substrates (2D and 3D), knowing only the gas phase HOMO-LUMO gaps and the DFT Kohn–Sham band gaps of the substrates. Specifically, this can be done by using the equation for the best-fit line in figure 4, as follows:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle E_{gap}^{{{V}_{scr}}}=E_{gap}^{gas-phase}+0.29E_{substrate}^{gap}-3.04.\nonumber \end{align} \tag{ 9 }$$

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The predicted HOMO-LUMO gaps, shown in figure 3 (hollow symbols), give remarkably good agreement with the benchmark GW calculations, given the simplicity of equation (9). To enable predictions of HOMO-LUMO gaps for small molecules with different adsorption heights, we provide the corresponding equations of best fit in table S3 for typical adsorption heights of 3.2–3.6 Å. These results give a more general version of equation (9) that depends on the adsorption height:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle E_{gap}^{{{V}_{scr}}}=&~E_{gap}^{gas-phase}+(-0.13d+0.71)E_{substrate}^{gap} \nonumber \\ & +(-0.82{{d}^{2}}+7.03d-17.20).\nonumber \end{align} \tag{ 10 }$$

Using equation (10), we obtain values of $E_{gap}^{{{V}_{scr}}}$ within 0.2% of those from equation (9) for d  =  3.25 Å. Equivalent equations using the GW QP substrate gap instead of the DFT Kohn–Sham band gap of the substrate are provided in the supplementary information.

In summary, we have used first principles calculations to study dielectric screening by 2D substrates, explicitly computing the electronic screening response at r, ${{V}_{scr}}(\boldsymbol{r},{\boldsymbol{r}^{\prime}})$, due to a point charge perturbation at r'. These calculations enable one to ascertain unambiguously the intrinsic role of 2D materials in screening charged impurities—a task that is not possible experimentally. Our results shed light on the interpretation of recent experiments [10–14] aimed at understanding the effects of screening by 2D substrates, and show that BN layers play an important role in screening out the effect of charged impurities, such as those present on underlying SiO₂ substrates. Atomically thin BN monolayers also screen remarkably well, enabled by the out-of-plane orbitals, which allow for the distribution of induced charge to generate a screening response. We also provide simple equations that allow one to estimate rather reliably the HOMO-LUMO gaps of small molecules physisorbed on 2D and 3D substrates, using only the gas phase HOMO-LUMO gap, the substrate band gap, and the adsorption height of the molecule. Taken together, our results provide important insights into the nature of 2D screening that will help both experimentalists and theorists to better understand and interpret the non-trivial interactions between 2D materials and adjacent charges.

Calculations were performed using BerkeleyGW [33] as well as our own scripts to compute the non-local screened Coulomb potential and the induced charge from the BerkeleyGW output. These scripts can be found at the following website: https://github.com/quek-group/wrrp. The mean-field wavefunctions and eigenvalues used as starting points for the BerkeleyGW calculations were computed using DFT as implemented in Quantum ESPRESSO [34, 35].

DFT wavefunctions and eigenvalues were computed within the local density approximation (LDA) for all systems with the exception of those involving MoS₂, where the PBE exchange-correlation functional [36] was used in conjunction with optimized norm-conserving Vanderbilt (ONCV) pseudopotentials [37, 38]. The ONCV pseudopotentials allow for the inclusion of Mo semicore electrons without a prohibitively high kinetic energy cutoff for wavefunctions. In MoS₂ systems, the inclusion of semicore states in the Mo pseudopotential has a negligible effect on the value of V_scr; however, the effect on the QP gap is significant (table S4). The kinetic energy cutoff for wavefunctions was set to 80 Ry (60 Ry) for LDA (PBE) pseudopotentials.

The static dielectric and self-energy matrices were computed using a 16 Ry cutoff for matrix elements. A 10 Ry cutoff for the sum over bands was imposed for all systems except for the benzene/BN (4  ×  4  ×  1 supercell) interface, which used a 6 Ry cutoff. All GW gaps were computed using the modified static remainder method [39] in order to ensure convergence with respect to the sum over unoccupied bands. The π  −  π^* gaps presented in table 1 were calculated at the Γ-point.

${{V}_{scr}}(\boldsymbol{r},{\boldsymbol{r}^{\prime}})$ was calculated following equation (1), where ${{\varepsilon}_{\boldsymbol G,{{\boldsymbol G}^{\prime}}}}(\boldsymbol q,\omega =0)$ is the static RPA dielectric matrix computed using BerkeleyGW. (We note that our code can handle dielectric matrices produced with or without the non-uniform neck sampling (NNS) method [40] for q-point sampling.) In our implementation, the inverse of the dielectric matrix is multiplied in reciprocal space by the bare Coulomb interaction to obtain ${{W}_{\boldsymbol G,{{\boldsymbol G}^{\prime}}}}(\boldsymbol q,\omega =0)$, from which the bare Coulomb interaction is subtracted to give ${{V}_{scr,\boldsymbol G,{{\boldsymbol G}^{\prime}}}}(\boldsymbol q,\omega =0)$. We then perform a 6D inverse Fourier transform to get ${{V}_{scr}}(\boldsymbol{r},{\boldsymbol{r}^{\prime}};\boldsymbol q)$. Finally, we calculate the screening potential as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{V}_{scr}}(\boldsymbol{r},{\boldsymbol{r}^{\prime}})=\sum\limits_{\boldsymbol q}{{{e}^{i{\boldsymbol q}\cdot (\boldsymbol r-{{\boldsymbol r}^{\prime}})}}}{{V}_{scr}}(\boldsymbol{r},{\boldsymbol{r}^{\prime}};\boldsymbol q).\nonumber \end{align} \tag{ 11 }$$

The bare Coulomb interaction used in both the GW and ${{V}_{scr}}(\boldsymbol{r},{\boldsymbol{r}^{\prime}})$ calculations was truncated [41] in the direction normal to the interface to avoid spurious interactions between periodic images.

The induced charge [21], ${{\rho}_{scr}}(\boldsymbol{r},{\boldsymbol{r}^{\prime}})$, which generates the screening response ${{V}_{scr}}(\boldsymbol{r},{\boldsymbol{r}^{\prime}})$, is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{V}_{scr}}(\boldsymbol{r},{\boldsymbol{r}^{\prime}})=-\int{d{\boldsymbol{r}^{\prime\prime}}\frac{{{e}^{2}}}{|\boldsymbol{r}-{\boldsymbol{r}^{\prime\prime}}|}{{\rho}_{scr}}({\boldsymbol{r}^{\prime\prime}},{\boldsymbol{r}^{\prime}})}.\nonumber \end{align} \tag{ 12 }$$

${{\rho}_{scr}}(\boldsymbol{r},{\boldsymbol{r}^{\prime}})$ is computed by first constructing the full RPA polarizability ${{\chi}_{\boldsymbol G,{{\boldsymbol G}^{\prime}}}}(\boldsymbol q,\omega =0)$ [42],

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\chi}_{\boldsymbol G,{{\boldsymbol G}^{\prime}}}}(\boldsymbol q,\omega =0)=v_{\boldsymbol G}^{-1}(\boldsymbol q)\left(\varepsilon _{_{\boldsymbol G,{{\boldsymbol G}^{\prime}}}}^{-1}(\boldsymbol q,\omega =0)-{{\delta}_{\boldsymbol G,{{\boldsymbol G}^{\prime}}}} \right),\nonumber \end{align} \tag{ 13 }$$

where v_G(q) is the Coulomb potential in reciprocal space. Next, assuming the external potential to be a point charge at r', we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\rho}_{scr,\boldsymbol G,{{\boldsymbol G}^{\prime}}}}(\boldsymbol q,\omega =0)={{\chi}_{\boldsymbol G,{{\boldsymbol G}^{\prime}}}}(\boldsymbol q,\omega =0){{v}_{{{\boldsymbol G}^{\prime}}}}(\boldsymbol q).\nonumber \end{align} \tag{ 14 }$$

We then perform a 6D inverse Fourier transform to get ${{\rho}_{scr}}(\boldsymbol{r},{\boldsymbol{r}^{\prime}};{\boldsymbol{q}})$ and obtain

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\rho}_{scr}}(\boldsymbol{r},{\boldsymbol{r}^{\prime}})=\sum\limits_{\boldsymbol q}{{{e}^{i \boldsymbol q\cdot (\boldsymbol r-{{\boldsymbol r}^{\prime}})}}}{{\rho}_{scr}}(\boldsymbol{r},{\boldsymbol{r}^{\prime}};\boldsymbol q).\nonumber \end{align} \tag{ 15 }$$

Interfaces were created by placing benzene 3.25 Å atop 3  ×  3  ×  1 substrate supercells (BN, 2L BN, MoS₂, and graphene), as well as a 4  ×  4  ×  1 supercell (BN). The Brillouin zones (BZ) of the cells were sampled via a uniform 6  ×  6  ×  1 grid of q-points, with the exception of the 3  ×  3  ×  1 benzene/BN cell, which used an 8  ×  8  ×  1 sampling to ensure commensurability with the corresponding 4  ×  4  ×  1 supercell. In all cases except benzene/2L BN we used the recently-implemented nonuniform neck sampling (NNS) method [40] to ensure convergence with respect to q-point sampling. The omission of NNS in the benzene/2L BN system is justified by tests on BN, which show that the use of NNS has a minimal effect on its band gap and V_scr for a 6  ×  6  ×  1 sampling of the BZ. The BZ of the bare substrates were sampled via 18  ×  18  ×  1 q-points and also make use of NNS.

All systems feature 15 Å of vacuum in order to separate periodic images normal to the surface. In separate tests, we verified that increasing the vacuum by a further 10 Å does not alter the value of V_scr at 3.25 Å above the surface.

Bader charges were calculated with the Bader Charge Analysis tool [43] using charge densities calculated in VASP [44].

The 2D autocorrelation, C(r), for a function $f(\bf\boldsymbol{r})$, was computed as in [45]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle C(r)=\frac{1}{2\pi}\int_{0}^{2\pi}{d\phi \left\langle \left\langle\,f({\boldsymbol r})\,f({\bf 0}) \right\rangle \right\rangle},\nonumber \end{align} \tag{ 16 }$$

where the angular brackets denote an average over the area and the integration averages over angular orientations. For the autocorrelation of the screened Coulomb potential, we note that the dielectric matrices used to construct W were computed using the unit cell with a Brillouin Zone sampling equivalent to the supercell size in order to prevent artificial periodicity within the cell area. For BN we used a 73  ×  73  ×  1 supercell, corresponding to a charge impurity density of 3.47  ×  10¹¹ cm⁻². The supercell sizes of MoS₂ and graphene, 58  ×  58  ×  1 and 73  ×  73  ×  1, respectively, were chosen to closely match the area, and therefore the impurity density, of the BN supercell.

We acknowledge support from Grant NRF-NRFF2013-07 from the National Research Foundation, Singapore and Grant MOE2016-T2-2-132 from the Ministry of Education, Singapore. We thank Miguel Dias Costa for his help in parallelizing codes used in this work and Shaffique Adam for discussions on [11]. Computations were performed on the NUS CA2DM cluster and National Supercomputing Centre Singapore (NSCC) under project IDs: 1100446 and 1100949. We acknowledge support from the Singapore National Research Foundation, Prime Minister's Office, under its medium-sized centre program.

The authors declare no competing financial interest. Keian Noori and Nicholas Lin Quan Cheng contributed equally to this work.