Dielectric screening by 2D substrates

Two-dimensional (2D) materials are increasingly being used as active components in nanoscale devices. Many interesting properties of 2D materials stem from the reduced and highly non-local electronic screening in two dimensions. 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. Thickness-dependent dielectric screening properties have recently been studied using electrostatic force microscopy (EFM) experiments. However, it was suggested that some of the thickness-dependent trends were due to extrinsic effects. Similarly, Kelvin probe measurements (KPM) indicate that charge fluctuations are reduced when BN slabs are placed on SiO2, but it is unclear if this effect is due to intrinsic screening from BN. In this work, we use first principles calculations to study the fully non-local dielectric screening properties of 2D material substrates. Our simulations give results in good qualitative agreement with those from EFM experiments, for hexagonal boron nitride (BN), graphene and MoS2, indicating that the experimentally observed thickness-dependent screening effects are intrinsic to the 2D materials. We further investigate explicitly the role of BN in lowering charge potential fluctuations arising from charge impurities on an underlying SiO2 substrate, as observed in the KPM experiments. 2D material substrates can also dramatically change the HOMO-LUMO gaps of adsorbates, especially for small molecules, such as benzene. We propose a reliable and very quick method to predict the HOMO-LUMO gaps of small physisorbed molecules on 2D and 3D substrates, using only the band gap of the substrate and the gas phase gap of the molecule.

3 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][4][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][7][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 SiO2 results in superior electron and hole mobilities as compared with graphene placed directly on SiO2. 9 Subsequent studies have shown that graphene placed on BN (and BN itself) has significantly lower charge density fluctuations than graphene on SiO2. 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 SiO2, 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, MoS2 and BN) placed on SiO2 suggest that electronic screening by atomically thin BN and MoS2 has relatively weak layer number dependence compared to few layer graphene. [12][13][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 4 water molecules. 12 Furthermore, the estimated density of defects in the experiment on BN was lower than that on MoS2 and graphene (BN: 10 8 cm -2 ; MoS2, graphene: 10 12 cm -2 ). [12][13][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][19][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 ( , ') scr V r r , the screening potential at r due to an electron point charge perturbation at r'. 21 The latter represents point charges on an underlying SiO2 substrate, on which a 2D material is placed, or represents the 5 quasiparticle 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][13][14] as well as Kelvin probe 11 experiments. We also provide a reliable back-of-the-envelope way to estimate the HOMO-LUMO gap 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 scr V W V =−, with V being the bare Coulomb potential and W the static screened Coulomb potential. In real space, the screening potential ( , ') scr V r r , 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 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 electrostatic force microscopy (EFM) experiments, [12][13][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 SiO2 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 where r' is taken to be ~3 Å beneath the bottom atomic plane (approximating the position of charge impurities on an SiO2 substrate), and rtop and rbottom are at the top and bottom atomic planes, 6 respectively. We estimate scr V  by taking rtop and rbottom to be directly above r', and averaging over the in-plane coordinates of the perturbing charge, r'. In Figure 1, we plot the quantity, which is the magnitude of the difference between scr V  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 Figure 1b). 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.

7
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 diff scr V  changes with thickness most quickly for few-layer graphene, and least quickly for BN.
However, the differences we predict are much less drastic than those observed in the experiments. [12][13][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  ( Figure 1). The value of  is deduced by fitting our ab initio results in (3) to the experimental data.
( ) The model results compare very well with the experimental data ( Figure 1). We note that our definition of scr V  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  should be taken to be an upper bound to the actual value of  in the experiment. Nonetheless, we comment that the value of  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  found for the graphene data (2.32) agrees well with the dielectric constant estimated in Ref. 13.  (4) and (5)).
It is well known that graphene devices placed atop BN substrates on SiO2 have much higher mobility than graphene devices placed directly atop SiO2. 9  W r r ), and limited our study to monolayers. From Figure 2a, 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 Ref. 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 2b, we show C(r) for W and V above 1L BN, 1L MoS2, 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),  Interestingly, the quantity ( , ') scr V r r is also useful for estimating the HOMO-LUMO gap 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 where j  represents the jth molecular orbital, and W  is the change in screened Coulomb potential in the molecule upon adsorption. Thus, For small molecules, we may approximate 15 It has been shown that for benzene physisorbed on bulk metallic or semiconducting substrates, 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 MoS2, 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 MoS2 (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 substrate and it is instructive to see how equations (7)(8) fare in these more realistic adsorption geometries. We compare the predictions of equations (7-8) with benchmark GW calculations of the same system in 3 x 3 supercells (as well as a 4x4 supercell for 1L BN) ( Figure 3). Notably, the 2D materials lead to significant renormalization of the benzene HOMO-  Table 1. Figure 3. GW π-π* gap of benzene on BN, MoS2 and graphene versus the π-π* gap predicted using Vscr. (equation (8)). The filled symbols denote the values predicted using the ab initio Vscr

13
(with a point charge model), while the hollow symbols denote values predicted using Vscr obtained from the best fit line in Figure 4.
GW π-π* gap of full  Figure S4. The same absolute isovalue is used for all plots. Blue denotes depletion of charge, while red denotes accumulation of charge.
As ( , ) scr V r 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 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) Using equation (10), we obtain values of scr V gap E within 0.2% of those from equation (9)  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 gap 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.  (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 (4x4x1 supercell) interface, which used a 6 Ry cutoff. All GW gaps were computed using the modified static remainder method 40 in order to ensure convergence with respect to the sum over unoccupied bands. The π-π* gaps presented in Table 1 were calculated at the Γ-point. , ' Then we have , , ' , ' ' We then perform a six-dimensional inverse Fourier transform to get ( , '; ) scr r r q supercell. In all cases except benzene/2L BN we used the recently-implemented nonuniform neck sampling (NNS) method 41 to ensure convergence with respect to k-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 Vscr for a 6x6x1 sampling of the BZ. The BZ of the bare substrates were sampled via 18x18x1 k-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 Vscr at 3.25 Å above the surface.
Bader charges were calculated with the Bader Charge Analysis tool 44 using charge densities calculated in VASP. 45 The 2D autocorrelation, C(r), for a function () fr, was computed as in Ref. 46: 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 73x73x1 supercell, corresponding to a charge impurity density of 3.47 x 10 11 cm -2 . The supercell sizes of MoS2 and graphene, 58x58x1 and 73x73x1, respectively, were chosen to closely match the area, and therefore the impurity density, of the BN supercell.
Author Contributions § Keian Noori, and Nicholas Lin Quan Cheng contributed equally to this work. Figure S1. Electric field as a function of the distance z from the surface of the material.      Tables   Table S1. Projection of isolated benzene monolayer HOMO and LUMO onto benzene/substrate interface states.

Classical electrostatic screening of a point charge
For a point charge of charge q at a distance d above the z-axis, the electric field is simply given by ( ) = 4πϵ 0 − | − | 3 (S1) For a point charge of charge q at a distance d above a semi-infinite conductor (z < 0), the potential and the electric field (for z > 0) is solved by adding an image charge of charge −q on the other side of the plane. The electric field is For a point charge of charge q at a distance d above a semi-infinite dielectric (z < 0) of dielectric constant , the image charge is instead The corresponding electric field is For a point charge of charge q at a distance d above the top surface (z = 0) of a dielectric slab with dielectric constant  and thickness t, the electric field is given to be 2 We plot the above 4 quantities in Figure S1, taking  = 5, d = 3 Å and t = 6 Å. This models the electric field due to a point charge perturbation 3 Å above a 2D material. Here, the electric field 28 for a thin dielectric slab is almost the same as that of a semi-infinite dielectric with the same dielectric constant, which implies that a thin dielectric slab can screen the electric field due to a point charge perturbation as well as a semi-infinite dielectric with the same dielectric constant.
This implies that the similarity in the behavior of 2D and 3D substrates in screening point charge perturbations adjacent to it can be explained using a simple classical model. Figure S1. Electric field as a function of the distance z from the surface of the material.

Relation of Vscr to the dielectric constant
Previous work on 3D semiconducting substrates demonstrates that there is a clear linear trend between the HOMO-LUMO π-π* gap of benzene and the band gap of the substrate (computed using LDA). 3 Interestingly, from our results, we find that both 2D and 3D substrates obey the same approximately linear relation between the substrate DFT gap and the screening response to a point charge 3.25 Å above the substrate (Figure 4). The dielectric screening capability of a material is most commonly described by its macroscopic dielectric constants, //  (out-of-plane) and  ⊥ (inplane). In Figures S2-3   The induced charge ( , ') scr rr  represents the charge rearrangement at r in the system due to a point charge perturbation at r'. We fix a negative point charge perturbation ~3.25 Å away from the surface and analyze the spatial extent of the induced charge by integrating the induced charge over a cylindrical region of radius r and height d as illustrated in Fig. S5. We maximize the radius of the cylindrical region to r = 6 Å and integrate the induced charge over the cylindrical region as a function of d. At d = 4.5 Å, most of the induced charge is captured in the cylindrical region, which indicated that a large amount of the induced charge is located close to the surface (Fig. S6). Interestingly, the amount of induced charge in the cylinder scales linearly as the planar-averaged Vscr(r,r) at 3.25 Å above the surface (Fig. S7). This observation is consistent with plots of the charge response for 2D and 3D surfaces (Fig. 5). Thus, the fact that a 2D material can screen nearby electronic excitations just as well as a 3D material with a similar GW gap may be explained by the observation that screening is predominantly a surface effect. From a classical viewpoint, it is known that the polarization-induced charges must lie only at the surface where there is a discontinuity in the dielectric environment. These quantum mechanical calculations, however, 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.   given in eV.

Tables
The best-fit equations for Vscr using the GW -corresponding to equations (9) and (10) in the main text -are:  Table S4. Effect of inclusion of semicore electrons and GW self-energy q→0 behavior on GW gap and Vscr in MoS2 layers. The treatment of q→0 in the GW self-energy (∑) is done via either a Monte Carlo (MC) averaging scheme or the nonuniform neck sampling (NNS) scheme. 20 All pseudopotentials use the LDA exchange-correlation functional.