The dielectric genome of van der Waals heterostructures
Kirsten Andersen
Simone Latini
Kristian S. Thygesen
Vertical stacking of two-dimensional (2D) crystals, such as graphene and hexagonal boron nitride, has recently lead to a new class of materials known as van der Waals heterostructures (vdWHs) with unique and highly tunable electronic properties. Ab-initio calculations should in principle provide a powerful tool for modeling and guiding the design of vdWHs, but in their traditional, form such calculations are only feasible for commensurable structures with a few layers. Here we show that the dielectric properties of realistic, incommensurable vdWHs comprising hundreds of layers can be calculated with ab-initio accuracy using a multi-scale approach where the dielectric functions of the individual layers (the dielectric building blocks) are coupled simply via their long-range Coulomb interaction. We use the method to illustrate the 2D-3D dielectric transition in multi-layer MoS$_2$ crystals, the hybridization of quantum plasmons in large graphene/hBN heterostructures, and to demonstrate the intricate effect of substrate screening on the non-Rydberg exciton series in supported WS$_2$.
paper.tex
\documentclass[twocolumn]{revtex4-1} %journal = nalefd, manuscript = letter,

\usepackage{graphicx}
\usepackage[export]{adjustbox}
\usepackage{amsmath}
\usepackage{subfig}
\usepackage[justification=justified]{caption}
\usepackage{verbatim}
\usepackage{color}
\usepackage{pdfpages}

\newcommand{\G}{\mathbf{G}}
\newcommand{\rr}{\mathbf{r}}
\newcommand{\chit}{\tilde{\chi}}
\newcommand{\R}{\mathbf{R}}
\newcommand{\q}{\mathbf{q}}
\newcommand{\qparf}{\mathbf{q}_\parallel}
\newcommand{\qpar}{\mathbf{q}_\parallel}
\newcommand{\Gpar}{\mathbf{G}_\parallel}
\newcommand{\rpar}{\mathbf{r}_\parallel}
\newcommand{\mr}{\mathrm}
\newcommand{\m}{\mathrm}
\newcommand{\mb}{\mathbf}
\newcommand{\ve}{\mathbf{v}_e}
\newcommand{\chim}{\boldsymbol{\chi}}
\newcommand{\chitm}{\boldsymbol{\tilde{\chi}}}
\newcommand{\Vtm}{\mathbf{\tilde{V}}}
\newcommand{\Vim}{\mathbf{V^I}}
\newcommand{\Vm}{\mathbf{V}}
\newcommand{\chiom}{\boldsymbol{\chi^0}}
\usepackage[abs]{overpic}

\begin{document}
\title{The dielectric genome of van der Waals heterostructures}

\author{Kirsten Andersen}
\email{kiran@fysik.dtu.dk}
\affiliation{Center for Atomic-scale Materials Design, Department of
Physics, Technical University of Denmark, DK - 2800 Kgs. Lyngby, Denmark}

\author{Simone Latini}
\affiliation{Center for Atomic-scale Materials Design, Department of Physics, Technical University of Denmark, DK - 2800 Kgs. Lyngby, Denmark}
\affiliation{Center for Nanostructured Graphene, Technical University of Denmark, DK - 2800 Kgs. Lyngby, Denmark}

\author{Kristian S. Thygesen}
\email{thygesen@fysik.dtu.dk}
%\affiliation{Center for Atomic-scale Materials Design, Department of
%Physics,  and Center for Nanostructured Graphene, Technical University of Denmark, DK - 2800 Kgs. Lyngby, Denmark}

\affiliation{Center for Atomic-scale Materials Design, Department of Physics, Technical University of Denmark, DK - 2800 Kgs. Lyngby, Denmark}
\affiliation{Center for Nanostructured Graphene, Technical University of Denmark, DK - 2800 Kgs. Lyngby, Denmark}

\keywords{van der Waals heterostructures, 2D materials, density-functional theory, dielectric function, excitons, plasmons}



\begin{abstract}
Vertical stacking of two-dimensional (2D) crystals, such as graphene and hexagonal boron nitride, has recently lead to a new class of materials known as van der Waals heterostructures (vdWHs) with unique and highly tunable electronic properties. Ab-initio calculations should in principle provide a powerful tool for modeling and guiding the design of vdWHs, but in their traditional, form such calculations are only feasible for commensurable structures with a few layers. Here we show that the dielectric properties of realistic, incommensurable vdWHs comprising hundreds of layers can be calculated with ab-initio accuracy using a multi-scale approach where the dielectric functions of the individual layers (the dielectric building blocks) are coupled simply via their long-range Coulomb interaction. We use the method to illustrate the 2D-3D dielectric transition in multi-layer MoS$_2$ crystals, the hybridization of quantum plasmons in large graphene/hBN heterostructures, and to demonstrate the intricate effect of substrate screening on the non-Rydberg exciton series in supported WS$_2$.
\end{abstract}
\maketitle


The class of 2D materials which started with graphene is rapidly expanding and now includes metallic and semiconducting transition metal dichalcogenides\cite{Wang2012} in addition to group III-V semi-metals, semiconductors and insulators\cite{Sahin2009}. These atomically thin materials exhibit unique opto-electronic properties with high technological potential\cite{Britnell2013,Choi2013, Shih2014,Woessner2014,Withers2015}. However, the 2D materials only form the basis of a new and much larger class of materials consisting of vertically stacked 2D crystals held together by weak van der Waals forces. In contrast to conventional heterostructures which require complex and expensive crystal-growth techniques to epitaxially grow the single-crystalline semiconductor layers, vdWHs can be stacked in ambient conditions with no requirements of lattice matching. The latter implies a weaker constraint, if any, on the choice of materials that can be combined into vdWHs.

The weak inter-layer binding suggests that the individual layers of a vdWH largely preserve their original 2D properties modified only by the long range Coulomb interaction with the surrounding layers. Turning this argument around, it should be possible to predict the overall properties of a vdWH from the properties of the individual layers. In this Letter we show that this can indeed be achieved for the dielectric properties. Conceptually, this extends the Lego brick picture used by Geim and Grigorieva\cite{Geim2013} for the atomic structure of a vdWH, to its dielectric properties. Specifically, we develop a semi-classical model which takes as input the dielectric functions of the individual isolated layers computed fully quantum mechanically and condensed into the simplest possible representation, and couple them together via the Coulomb interaction, see Figure \ref{scheme}. Despite the complete neglect of interlayer hybridization, the model provides an excellent account of both the spatial and dynamical dielectric properties of vdWHs. The condensed representation of the dielectric function of all isolated 2D crystals can thus be regarded as the dielectric genome of vdWHs. 

In addition to its conceptual value, our approach overcomes a practical limitation of conventional first-principles methods. Such methods are not only computationally demanding, but also rely on periodic boundary conditions which are incompatible with the incommensurable interfaces found in vdWHs. In fact, for many purposes, an in-plane lattice mismatch between neighbouring 2D crystals is preferred because it reduces the interlayer coupling and thus minimises the risk of commensurate-incommensurate transitions\cite{Woods2014}, and formation of Moire patterns\cite{Kang2013} and associated band structure reconstructions\cite{Lu2014} which are typical for systems with similar lattice constants. This emphasises the need for alternative approaches for modelling vdWHs.


%%%%%%%%%
%METHOD
%%%%%%%%%

\begin{figure*}
  \includegraphics[width = 0.9\linewidth]{Figure1.pdf}
   \caption{Schematic of the QEH model. (a) The density response function and dielectric function of the heterostructure are calculated from the dielectric building blocks of the individual layers assuming a purely electrostatic interaction between the layers. The dielectric building blocks are calculated ab-initio for the isolated layers. They comprise the monopole and dipole components of the density response function, $\tilde \chi_{M/D}$, together with the spatial shape of the electron density, $\rho_{M/D}(z)$, induced by a constant and linear applied potential, respectively. (b) Monopole and dipole induced densities (blue) together with the associated potentials (red) for monolayer MoS$_2$.}
 \label{scheme}
\end{figure*}

The dielectric function is one of the most important material response functions. It determines the effective interaction between charged particles in the material, contains information about the collective oscillations of the electron gas (plasmons)\cite{Pitarke2007}, and enters as a fundamental ingredient in many-body calculations of e.g. excitons and quasiparticle band structures\cite{Onida2002,Hybertsen1986}. 

The (inverse) dielectric function is related to the electron density response function, $\chi$, via 
\begin{equation}\label{eq:epsilon}
\epsilon^{-1}(\rr,\rr', \omega) = \delta(\rr-\rr')+\int \frac{1}{\vert \rr-\rr''\vert } \chi(\rr'', \rr', \omega)d\rr''.
\end{equation}
In our quantum-electrostatic heterostructure (QEH) model the calculation of the dielectric function is divided into two parts. In the first part the in-plane averaged density response function of each of the freestanding layers, $\chi_i(z, z', \qpar, \omega)$, are obtained from ab-initio calculations. In practice we treat the in-plane momentum transfer, $\qpar$, as a scalar since  most 2D materials are isotropic within the plane. From $\chi_i$ we calculate the magnitude of the monopole/dipole component of the density induced by a potential with a constant/linear variation across the layer and in-plane variation $\exp(i \rpar \cdot \qpar)$:
\begin{align}\label{eq:chi}
&\tilde \chi_{i\alpha}(\qpar, \omega) = \int z^\alpha \chi_i(z, z', \qpar, \omega) z'^\alpha dz dz'.
\end{align}
Here $\alpha=0,1$ for the monopole and dipole components, respectively. In addition we calculate the spatial form of the induced density, $\rho_{i\alpha}(z,\qpar)$. With a proper normalization of $\rho_{i\alpha}$ we can then write 
\begin{align}\label{eq:rho}
\int \chi_i(z, z',\qpar,\omega) z'^\alpha dz' = \tilde{\chi}_{i\alpha}(\qpar,\omega) \rho_{i\alpha}(z,\qpar)
\end{align}
We have found that while $\tilde \chi_{i\alpha}$ depends strongly on frequency, $\rho_{i\alpha}$ does not. The data set $(\tilde \chi_{i\alpha},\rho_{i\alpha})$ with $\alpha=0,1$ or equivalently $\alpha=M,D$ constitutes the dielectric building block of layer $i$, as illustrated in Figure \ref{scheme}. According to Eq. (\ref{eq:rho}) the dielectric building block allows us to obtain the density induced in the (isolated) layer $i$ by a constant/linear potential. It is straightforward to extend the dielectric building blocks to account for higher-order moments in the induced density described by $\alpha>1$, but we have found the dipole approximation to be sufficient in all cases considered.

In the second part of the QEH model, the density response function of the vdWH in the discrete monopole/dipole representation is obtained by solving a Dyson equation that couples the dielectric building blocks together via the Coulomb interaction. The Dyson equation for the full density response function giving the magnitude of the monopole/dipole density on layer $i$ induced by a constant/linear potential applied to layer $j$, reads (omitting the $\qpar$ and $\omega$ variables for simplicity)
\begin{align}\label{eq:dysonhet}
\chi_{i\alpha, j\beta}& =
\tilde \chi_{i\alpha} \delta_{i\alpha,j\beta} \: + \tilde \chi_{i\alpha} \sum_{k\neq i ,\gamma} V_{i\alpha,k\gamma}  \, \chi_{k \gamma, j\beta}.
\end{align}
The Coulomb matrices are defined as
\begin{align}\label{eq:coulomb}
V_{i\alpha,k\gamma} (\qpar)=\int \rho_{i\alpha}(z,\qpar) \Phi_{k\gamma}(z,\qpar) dz
\end{align}
where $\Phi_{k\gamma}$ is the potential associated with the induced density, $\rho_{k\gamma}$, which we calculate on a uniform grid by solving a 1D Poisson equation. Note that we leave out the self-interaction terms in Eq. (\ref{eq:dysonhet}) since the intralayer Coulomb interaction is already accounted for by the uncoupled $\tilde \chi_{i\alpha}$. The (inverse) dielectric function of Eq. (\ref{eq:epsilon}) in the monopole/dipole basis becomes
\begin{align}\label{eq:eps}
\epsilon^{-1}_{i\alpha, j\beta}&(\qpar,\omega) =
\delta_{i\alpha,j\beta}  +  \sum_{k \gamma} V_{i\alpha,k\gamma} (\qpar) \, \chi_{k \gamma, j\beta}(\qpar,\omega).
\end{align}
More details on the method and computations are provided in the supporting information.

A database containing the dielectric building blocks of a large collection of 2D materials has been constructed, and is available from our website \cite{database}. It presently contains more than 50 transition metal dichalcogenides and oxides, graphene at different doping levels, and hBN, and more materials are being added. From here the data files can be downloaded together with a Python module for calculating the dielectric function and associated properties of any combination of these materials. QEH model calculations for vdWHs containing a few hundred layers can be performed on a standard PC. In Figure \ref{database} we show the $\qpar$-dependent static dielectric functions of the monolayer transition metal dichalcogenides and -oxides presently contained in our database (for a complete overview of the materials see Ref.~\cite{rasmussen_computational_2015}).   All the dielectric functions show the same qualitative form, in particular they become 1 for $\qpar\to0$ and $\qpar \to \infty$, however there is quite some variation in their magnitude. As expected the size of the dielectric function correlates well with the size of the band gap of the material indicated by the colour. 

\begin{figure}
  \includegraphics[width = \linewidth]{Figure2.pdf}
   \caption{The static dielectric function $\epsilon(\qpar,\omega=0)$ of the 51 transition metal dichalcogenides and oxides included in the database, that span a large range of magnitudes. The relation to the quasiparticle G$_0$W$_0$ band gap of the materials, calculated in Ref.~\cite{rasmussen_computational_2015} is shown in color.}
 \label{database}
\end{figure}


First-principles calculations were performed with the GPAW code\cite{GPAW2,Yan2011a}. Single-particle wave functions and energies were calculated within the local density approximation (LDA) using 400 eV plane wave cut-off and at least $45\times45$ sampling of the 2D Brillouin zone. Density response functions and dielectric functions were calculated within the random phase approximation (RPA). The RPA does not include electron-hole interactions, but generally yields good results for the static dielectric properties of semi-conductors and dynamical response of metals. Except for MoS$_2$ bulk, we included at least 15 \AA~of vacuum in the super cells perpendicular to the layers and applied a truncated Coulomb kernel to avoid long range screening between periodically repeated structures. All response functions were calculated in a plane wave basis including reciprocal lattice vectors up to at least 50 eV. A similar cut off was used for the sum over empty states and convergence was carefully checked. The frequency dependence of the response functions was represented on a non-linear frequency grid ranging from 0 to 35 eV, with an initial grid spacing of 0.02 eV. All details of the calculations and atomic structure geometries are provided in the supporting information.

 \begin{figure*}
\includegraphics[width = \linewidth]{Figure3.pdf}
  	  \caption{2D-3D transition of the dielectric function. (a) Atomic structure of a 20 layer MoS$_2$ slab. (b) The macroscopic static dielectric function $\epsilon_M(\qpar, \omega = 0)$ as a function of the in-plane momentum transfer for different number of layers, $N$. The macroscopic dielectric function relates the total potential averaged over the width of the slab to an external potential of the form $V_{ext}(r_\Vert,z)=\exp(ir_\Vert \qpar)$. The dielectric functions increase monotonically with $N$ converging slowly towards the dielectric function of bulk MoS$_2$ obtained from an ab-initio calculation. Excellent agreement between the QEH model and the ab initio results are seen for $N=1,2$. The slow convergence towards the bulk result is due to the strong spatial variation of the induced potential in the surface region of the slabs. This can be seen in panel (c) which shows $V_{ext}/V_{tot}(z)$, i.e.~the local dielectric function, for an external potential constant across the slab and with in-plane wave vector $\qpar=0.036$\AA$^{-1}$ for $N=50$.}
  \label{mos2macro}
\end{figure*}

As a first application of the QEH model, we study how the (static) dielectric function of a 2D material evolves as the layer thickness increases towards the bulk. One of the most characteristic differences between 2D and 3D materials is the behaviour of the dielectric function in the long wave length limit: For a bulk semiconductor, the dielectric function $\epsilon(\qpar)$ tends smoothly to a constant value larger than unity as $\qpar \to 0$. In contrast $\epsilon(\qpar)=1+O(\qpar)$ for a 2D semiconductor implying a complete absence of screening in the long wave length limit \cite{Cudazzo2011,Huser2013}. 

Ab initio calculations were performed for the dielectric function of MoS$_2$ monolayer, bilayer, and bulk, and the QEH model was used for multilayer structures up to 100 layers. Figure \ref{mos2macro} (b) shows the dielectric functions averaged over the slabs, i.e.~the macroscopic dielectric function, as function of the in-plane momentum transfer. For large $\qpar$ the dielectric functions show similar behavior. However, whereas $\epsilon(0)=14$ for the bulk, the dielectric functions of the slabs decrease sharply to 1 for small $\qpar$. This demonstrates that the dielectric properties of a vdWH of thickness $L$ are 2D like for $\qpar \ll 1/L$ and 3D like for $\qpar \gg 1/L$. Interestingly, also the result for bulk MoS$_2$ shows reminiscence of the 2D nature of the constituent layers, where the magnitude of the dielectric function has a slight drop when $\qpar \to 0$.

The QEH model describes the change in the dielectric function from mono- to bilayer very accurately in spite of the well known differences between the mono- and bilayer band structures\cite{Cheiwchanchamnangij2012}. This shows that hybridisation driven band structure effects, i.e.~quantum confinement, have negligible influence on the dielectric properties of a vdWH and is the main reason for the success of the QEH model. The model result seems to converge towards the ab initio bulk result, however, convergence is not fully reached even for $N=100$. The slow convergence towards the bulk result is mainly due to the spatial variation of the induced potential across the slab. In Figure \ref{mos2macro} (c) we show the $z$-dependent dielectric function defined as $\epsilon(z)=V_{ext}/V_{tot}(z)$, for a constant (along $z$) external potential with a long wavelength in-plane variation for $N=50$. Although $\epsilon(z)$ is close to the ab-initio bulk value (dashed line) in the middle of the slab, screening is strongly suppressed in the surface region. Increasing the slab thickness beyond 50 layers brings the QEH result even closer to the bulk result in the middle of the slab, but a small underestimation remains originating from the difference in the band structures of the monolayer and bulk systems. The suppressed screening in the surface region is a direct consequence of the anisotropic nature of the layered MoS$_2$ crystals which limits the screening of perpendicular fields relative to in-plane fields, and is expected to be a general property of vdWHs.

The model can also be used to calculate the response to fields polarized along the $z$-direction, i.e. perpendicular to the layers. In this case the perpendicular component, $\epsilon_{zz}(q_z=0)$, can be calculated by applying an external potential with a linear variation along $z$. In the discrete basis of the QEH model, such a field is represented by a vector with 0 for all monopole components and 1 for all dipole components. Comparing the averaged slope of the total potential to the slope of the applied linear potential for a slab of N=100 layers of MoS2 yields $\epsilon_{zz}= 7.8$.  This value is somewhat larger than the bulk value of 6.03, however, due to long range surface effects it is not necessarily to be expected that the two numbers should coincide. In fact, we find excellent agreement between the QEH model and full ab-initio calculation of $\epsilon_{zz}$ for a four layer MoS$_2$ slab (see supporting information).

\begin{figure*}
 \includegraphics[width = \linewidth]{Figure4.pdf}
  \caption{Plasmons in graphene/hBN heterostructures. (a) Two graphene sheets separated by three layers of hBN.  (b) Eigenvalues of the heterostructure dielectric function $\epsilon(\omega)$. Only the two eigenvalue curves that fullfill the plasmon condition $\mr{Re}\epsilon_n(\omega_P)=0$ are shown. (c) The eigen-potential, $\psi(\omega_P)$, and associated density, $\rho(\omega_P)$, of the plasmon modes. The plasmons  correspond to the antisymmetric ($+-$) and symmetric ($++$) combinations of the isolated graphene plasmons. (d) Plasmon dispersion for heterostructures containing 1 and 3 layers of hBN. Full lines denote the QEH model while ab-initio results are denoted by symbols. (e+f) Energy and weight of the plasmon modes for up to 20 layers hBN between the graphene sheets. Results for equivalent structures with vacuum filling the gap are also shown. Dashed black lines indicate the plasmon energy and weight in an isolated graphene sheet. Overall, the QEH model is in excellent agreement with the full ab initio calculations performed for up to 3 layers hBN.}
  \label{graphene-BN}
\end{figure*}

Next, we consider the hybridisation of plasmons in graphene sheets separated by hBN, see Figure \ref{graphene-BN}(a). Plasmons in graphene on hBN were recently found to propagate with low loss \cite{Woessner2014}, and the close to perfect lattice match between the two layers enables full ab initio calculations for the thinnest heterostructures. Here we use doped graphene that has a finite density of states at the Fermi level, giving rise to two-dimensional plasmons with energies in the regime 0-2 eV. The plasmon energies goes to zero in the optical limit, $\qpar \rightarrow 0$ as is characteristic for plasmons in 2D metals\cite{Hwang2007, Shin2011}. We calculate the effect of hBN on the plasmons using the QEH model for up to 20 layers of hBN and compare to full ab-initio calculations for 1-3 layers of hBN.

To identify the plasmons of the heterostructure we follow Ref. \cite{Andersen2012}. In brief, we compute the eigenvalues, $\epsilon_n(\omega)$, of the heterostructure dielectric function for each frequency point and identify a plasmon energy, $\hbar \omega_P$, from the condition $\text{Re}\epsilon_n(\omega_P)=0$, see Figure \ref{graphene-BN}(b). The corresponding eigenvector, $\phi_n(\omega_P)$, represents the potential associated with the plasmon oscillation, see panel (c). This analysis identifies two plasmons corresponding to the symmetric ($++$) and antisymmetric ($+-$) combinations of the graphene plasmons as previously found for two freestanding graphene sheets \cite{Hwang2009}. For 1-3 hBN layers, the QEH model perfectly reproduces the ab-initio results for the dielectric eigenvalues, plasmon energy, and weight. The latter was defined as the area under the peaks in the loss function $-\mr{Im}\epsilon^{-1}(\qpar,\omega)$, see panel (b). The densities and potentials of the plasmon eigenmodes shown in panel (c) are also reproduced fairly accurately by the model, where the qualitative differences for the induced densities, $\rho(z)$, are due to the use of a limited basis of the monopole and dipole response for each layer. In panels (e-f) the result of full ab-initio calculations are shown by symbols while the QEH results are shown by continuous lines. The effect of the hBN buffer (dashed lines) is to red shift and damp the plasmons compared to the result for two graphene sheets separated by the same amount of vacuum (full lines). This is also reflected by the relatively large amount of electron density located on the hBN during the plasma oscillation, see panel (c).

\begin{figure*}
  \includegraphics[width = \linewidth]{Figure5.pdf}
  \caption{Excitons in supported WS$_2$. (a) Monolayer WS$_2$ adsorbed on a h-BN substrate. (b) The screened interaction between an electron and a hole localised within a WS$_2$ monolayer adsorbed on hBN. For comparison the unscreened $1/r$ potential is shown. The radial probability distribution of the first five excitons, $r\vert F(r)\vert^2$, are also shown (arbitrary normalization). (c) The calculated binding energies of the lowest five excitons in freestanding WS$_2$ (green dashed) and WS$_2$ on hBN (blue) and MoS$_2$ (cyan). Experimental values from Ref. \cite{Chernikov2014} for WS$_2$ on SiO$_2$ are shown in red. The 2D hydrogen model with a $1/\epsilon r$ potential is shown for $\epsilon=1.7$. (d) The dielectric function of the WS$_2$ layer defined as $\epsilon(q)=V(q)/W(q)$, where $V(q)$ and $W(q)$ are the bare and screened interaction in the WS$_2$ layer, respectively. (e) The screened interaction in the WS$_2$ layer as function of $\log(r)$. (f) The relative difference between the screened interaction in the supported and freestanding WS$_2$. Inset shows the relative difference between $E_b$ for the supported and freestanding WS$_2$.}
  \label{fig:excitons}
\end{figure*}

Finally, we explore some characteristic features of excitons in freestanding and supported 2D semiconductors. A straight forward generalisation of the well known Mott-Wannier model\cite{Wannier1937} leads to the following eigenvalue equation for the excitons of a 2D semiconductor\cite{Cudazzo2011, Berkelbach2013}:
%
%
%
\begin{equation}
\left[-\frac{\nabla_{2D}^2}{2\mu_{ex}}+W(\rr)\right]F(\rr)=E_bF(\rr),
\label{eq:MWHamiltonian}
\end{equation}
%
%
%
where $E_b$ is the exciton binding energy, $F(\rr)$ is the wave function, $\mu_{ex}$ is the effective mass, and $W(\rr)$ is the screened electron-hole interaction. Assuming that the electron and hole are localised in layer 1, the Fourier transformed screened electron-hole interaction is obtained from the static ($\omega=0$) response function Eq. (\ref{eq:dysonhet}) and Coulomb interaction matrix Eq. (\ref{eq:coulomb}) of the QEH model,
%
%
%
\begin{equation}
W(\qpar) = V_{1M,1M}+\sum_{i \alpha,j\beta} V_{1 M,j \beta}(\qpar) \chi_{j \beta, i \alpha}(\qpar)V_{i \alpha, 1 M} (\qpar).
\label{eq:W_q}
\end{equation}
%
%
%
The first term is the bare, i.e. unscreened, electron-hole interaction in layer 1 under the assumption that the electron and hole densities can be represented by the induced monopole density, $\rho_{1M}(z)$. The second term describes the screening from the surrounding layers and layer 1 itself. Note that the above equation can be easily generalised to describe the screened interaction between charges localised in different layers (relevant for indirect excitons). 

In Ref. \cite{Chernikov2014} Chernikov et al. observed a peculiar non-hydrogenic Rydberg series for the excitons in a single layer of WS$_2$ adsorbed on a SiO$_2$ substrate. Here we use the QEH model to calculate the screened electron-hole interaction within the WS$_2$ layer from the dielectric function of the full heterostructure. Since the QEH is applicable only to layered materials we place WS$_2$ on a 100 layer thick slab of hBN which has dielectric constant very similar to that of SiO$_2$ (both around 4). For comparison we performed similar calculations using MoS$_2$ as substrate (dielectric constant larger than SiO$_2$). Figure \ref{fig:excitons} (c) shows the five lowest $s$-excitons calculated from Eq. (\ref{eq:MWHamiltonian}) for both freestanding and supported WS$_2$. For freestanding WS$_2$, we obtain $E_b=0.59$ eV for the lowest exciton in good agreement with previous ab-initio calculations\cite{Shi2013}. The enhanced screening from the substrate lowers the exciton binding energies bringing the entire series closer to the experimental values (red), in particular for the hBN substrate.  

The dielectric function of the WS$_2$ layer defined as $\epsilon(q)=V(q)/W(q)$, where $V(q)$ and $W(q)$ are the bare and screened interaction in the WS$_2$ layer, respectively. Figure \ref{fig:excitons} (d) shows that the dielectric function of the supported WS$_2$ layer exceeds unity in the $\qpar\to 0$ limit. For structures of finite width, $L$, the dielectric function will in practice tend to unity for very small $\qpar \ll 1/L$. Here the result have been extrapolated to infinite substrate thickness, where $\epsilon(\qpar)$ tends to a value larger than unity. This means that the nature of the screening within the layer is not strictly 2D because the bulk substrate is able to screen the long wave length fields. In real space, the screened potentials diverge as $\log(r)$ for small $r$ and decay as $1/r$ for large $r$, see panel (e). In panel (f) we show how the substrate affects $W(r)$: The relative deviation from $W(r)$ of the freestanding layer vanishes for small and large $r$ but becomes significant at intermediate distances. As a consequence, the substrate-induced change in the exciton binding energy is relatively larger for intermediate exciton sizes. These results clearly demonstrate the profound, nonlocal influence of substrates on the dielectric screening and excitations in 2D materials.  

In conclusion, we have demonstrated that the spatial and dynamical dielectric properties of a vdWH can be accurately and efficiently obtained from the dielectric properties of its constituent 2D crystals. The presented quantum-electrostatic heterostructure model (QEH) exploits this feature and enables the calculation of the dielectric properties and collective electronic excitations of realistic incommensurable heterostructures with ab-initio precision. The dielectric building blocks for more than fifty different 2D materials are available in an open database allowing 2D materials researchers to efficiently predict and design the dielectric properties of realistic vdWHs.

\section*{acknowledgement}
The authors thank Karsten Jacobsen for inspiring discussions. The authors acknowledge support from the Danish Council for Independent Research's Sapere Aude Program through grant no. 11-1051390. The Center for Nanostructured Graphene (CNG) is sponsored by the Danish National Research Foundation, Project DNRF58.


\section*{Supporting information available}
Detailed description of our quantum-electrostatic heterostructure (QEH) model and the computational details for all the ab-initio calculations are given in the supporting information.  


\bibliography{bibtex}

\cleardoublepage
\onecolumngrid
\section*{}
\includepdf[pages={1-6}]{supinfo.pdf}
\end{document}




paper.bbl
\providecommand*\mcitethebibliography{\thebibliography}
\csname @ifundefined\endcsname{endmcitethebibliography}
  {\let\endmcitethebibliography\endthebibliography}{}
\begin{mcitethebibliography}{30}
\providecommand*\natexlab[1]{#1}
\providecommand*\mciteSetBstSublistMode[1]{}
\providecommand*\mciteSetBstMaxWidthForm[2]{}
\providecommand*\mciteBstWouldAddEndPuncttrue
  {\def\EndOfBibitem{\unskip.}}
\providecommand*\mciteBstWouldAddEndPunctfalse
  {\let\EndOfBibitem\relax}
\providecommand*\mciteSetBstMidEndSepPunct[3]{}
\providecommand*\mciteSetBstSublistLabelBeginEnd[3]{}
\providecommand*\EndOfBibitem{}
\mciteSetBstSublistMode{f}
\mciteSetBstMaxWidthForm{subitem}{(\alph{mcitesubitemcount})}
\mciteSetBstSublistLabelBeginEnd
  {\mcitemaxwidthsubitemform\space}
  {\relax}
  {\relax}

\bibitem[Wang et~al.(2012)Wang, Kalantar-Zadeh, Kis, Coleman, and
  Strano]{Wang2012}
Wang,~Q.~H.; Kalantar-Zadeh,~K.; Kis,~A.; Coleman,~J.~N.; Strano,~M.~S.
  \emph{Nat. Nanotechnol.} \textbf{2012}, \emph{7}, 699--712\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem[{\c{S}}ahin et~al.(2009){\c{S}}ahin, Cahangirov, Topsakal, Bekaroglu,
  Akturk, Senger, and Ciraci]{Sahin2009}
{\c{S}}ahin,~H.; Cahangirov,~S.; Topsakal,~M.; Bekaroglu,~E.; Akturk,~E.;
  Senger,~R.~T.; Ciraci,~S. \emph{Phys. Rev. B} \textbf{2009}, \emph{80},
  155453\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem[Britnell et~al.(2013)Britnell, Ribeiro, Eckmann, Jalil, Belle,
  Mishchenko, Kim, Gorbachev, Georgiou, Morozov, Grigorenko, Geim, Casiraghi,
  {Castro Neto}, and Novoselov]{Britnell2013}
Britnell,~L.; Ribeiro,~R.~M.; Eckmann,~a.; Jalil,~R.; Belle,~B.~D.;
  Mishchenko,~a.; Kim,~Y.-J.; Gorbachev,~R.~V.; Georgiou,~T.; Morozov,~S.~V.;
  Grigorenko,~a.~N.; Geim,~a.~K.; Casiraghi,~C.; {Castro Neto},~a.~H.;
  Novoselov,~K.~S. \emph{Science} \textbf{2013}, \emph{340}, 1311--4\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem[Sup~Choi et~al.(2013)Sup~Choi, Lee, Yu, Lee, Hwan~Lee, Kim, Hone, and
  Jong~Yoo]{Choi2013}
Sup~Choi,~M.; Lee,~G.-H.; Yu,~Y.-J.; Lee,~D.-Y.; Hwan~Lee,~S.; Kim,~P.;
  Hone,~J.; Jong~Yoo,~W. \emph{Nat. Commun.} \textbf{2013}, \emph{4},
  1624\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem[Shih et~al.(2014)Shih, Wang, Son, Jin, Blankschtein, and
  Strano]{Shih2014}
Shih,~C.~J.; Wang,~Q.~H.; Son,~Y.; Jin,~Z.; Blankschtein,~D.; Strano,~M.~S.
  \emph{ACS Nano} \textbf{2014}, \emph{8}, 5790--5798\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem[Woessner et~al.(2015)Woessner, Lundeberg, Gao, Principi,
  Alonso-Gonz{\'a}lez, Carrega, Watanabe, Taniguchi, Vignale, Polini, Hone,
  Hillenbrand, and Koppens]{Woessner2014}
Woessner,~A.; Lundeberg,~M.~B.; Gao,~Y.; Principi,~A.; Alonso-Gonz{\'a}lez,~P.;
  Carrega,~M.; Watanabe,~K.; Taniguchi,~T.; Vignale,~G.; Polini,~M.; Hone,~J.;
  Hillenbrand,~R.; Koppens,~F. H.~L. \emph{Nat Mater} \textbf{2015}, \emph{14},
  421--425\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem[Withers et~al.(2015)Withers, Del Pozo-Zamudio, Mishchenko, Rooney,
  Gholinia, Watanabe, Taniguchi, Haigh, Geim, Tartakovskii, and
  Novoselov]{Withers2015}
Withers,~F.; Del Pozo-Zamudio,~O.; Mishchenko,~A.; Rooney,~A.~P.; Gholinia,~A.;
  Watanabe,~K.; Taniguchi,~T.; Haigh,~S.~J.; Geim,~A.~K.; Tartakovskii,~A.~I.;
  Novoselov,~K.~S. \emph{Nat. Mater.} \textbf{2015}, \emph{14}, 301--306\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem[Geim and Grigorieva(2013)Geim, and Grigorieva]{Geim2013}
Geim,~a.~K.; Grigorieva,~I.~V. \emph{Nature} \textbf{2013}, \emph{499},
  419--25\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem[Woods et~al.(2014)Woods, Britnell, Eckmann, Ma, Lu, Guo, Lin, Yu, Cao,
  Gorbachev, Kretinin, Park, Ponomarenko, Katsnelson, Gornostyrev, Watanabe,
  Taniguchi, Casiraghi, Gao, Geim, and Novoselov]{Woods2014}
Woods,~C.~R. et~al.  \emph{Nature Phys.} \textbf{2014}, \emph{10},
  451--456\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem[Kang et~al.(2013)Kang, Li, Li, Xia, and Wang]{Kang2013}
Kang,~J.; Li,~J.; Li,~S.~S.; Xia,~J.~B.; Wang,~L.~W. \emph{Nano Lett.}
  \textbf{2013}, \emph{13}, 5485--5490\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem[Lu et~al.(2014)Lu, Li, Watanabe, Taniguchi, and Andrei]{Lu2014}
Lu,~C.-P.; Li,~G.; Watanabe,~K.; Taniguchi,~T.; Andrei,~E. \emph{Phys. Rev.
  Lett.} \textbf{2014}, \emph{113}, 156804\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem[Pitarke et~al.(2007)Pitarke, Silkin, Chulkov, and
  Echenique]{Pitarke2007}
Pitarke,~J.~M.; Silkin,~V.~M.; Chulkov,~E.~V.; Echenique,~P.~M. \emph{Rep.
  Prog. in Phys.} \textbf{2007}, \emph{70}, 1--87\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem[Onida et~al.(2002)Onida, Reining, and Rubio]{Onida2002}
Onida,~G.; Reining,~L.; Rubio,~A. \emph{Rev. Mod. Phys.} \textbf{2002},
  \emph{74}, 601--659\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem[Hybertsen and Louie(1986)Hybertsen, and Louie]{Hybertsen1986}
Hybertsen,~M.~S.; Louie,~S.~G. \emph{Phys. Rev. B} \textbf{1986}, \emph{34},
  5390\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem[dat()]{database}
The dielectric building blocks and QEH software can be downloaded at
  https://cmr.fysik.dtu.dk/vdwh/vdwh.html\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem[Rasmussen and Thygesen(2015)Rasmussen, and
  Thygesen]{rasmussen_computational_2015}
Rasmussen,~F.~A.; Thygesen,~K.~S. \emph{J. Phys. Chem. C} \textbf{2015},
  Accepted\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem[Enkovaara and et~al.(2010)Enkovaara, and et~al.]{GPAW2}
Enkovaara,~J.; et~al., \emph{J. Phys. Condens. Matte} \textbf{2010}, \emph{22},
  253202\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem[Yan et~al.(2011)Yan, Mortensen, Jacobsen, and Thygesen]{Yan2011a}
Yan,~J.; Mortensen,~J.~J.; Jacobsen,~K.~W.; Thygesen,~K.~S. \emph{Phys. Rev. B}
  \textbf{2011}, \emph{83}, 245122\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem[Cudazzo et~al.(2011)Cudazzo, Tokatly, and Rubio]{Cudazzo2011}
Cudazzo,~P.; Tokatly,~I.~V.; Rubio,~A. \emph{Phys. Rev. B} \textbf{2011},
  \emph{84}, 085406\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem[H\"{u}ser et~al.(2013)H\"{u}ser, Olsen, and Thygesen]{Huser2013}
H\"{u}ser,~F.; Olsen,~T.; Thygesen,~K.~S. \emph{Phys. Rev. B} \textbf{2013},
  \emph{87}, 1--14\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem[Cheiwchanchamnangij and Lambrecht(2012)Cheiwchanchamnangij, and
  Lambrecht]{Cheiwchanchamnangij2012}
Cheiwchanchamnangij,~T.; Lambrecht,~W. R.~L. \emph{Phys. Rev. B} \textbf{2012},
  \emph{85}, 1--4\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem[Hwang and Sarma(2007)Hwang, and Sarma]{Hwang2007}
Hwang,~E.; Sarma,~S. \emph{Phys. Rev. B} \textbf{2007}, \emph{75}, 205418\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem[Shin et~al.(2011)Shin, Kim, Kim, Kim, Noh, Kim, and Chung]{Shin2011}
Shin,~S.~Y.; Kim,~N.~D.; Kim,~J.~G.; Kim,~K.~S.; Noh,~D.~Y.; Kim,~K.~S.;
  Chung,~J.~W. \emph{Appl. Phys. Lett.} \textbf{2011}, \emph{99}, 082110\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem[Andersen et~al.(2012)Andersen, Jacobsen, and Thygesen]{Andersen2012}
Andersen,~K.; Jacobsen,~K.~W.; Thygesen,~K.~S. \emph{Phys. Rev. B}
  \textbf{2012}, \emph{86}, 245129\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem[Hwang and {Das Sarma}(2009)Hwang, and {Das Sarma}]{Hwang2009}
Hwang,~E.; {Das Sarma},~S. \emph{Phys. Rev. B} \textbf{2009}, \emph{80},
  205405\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem[Chernikov et~al.(2014)Chernikov, Berkelbach, Hill, Rigosi, Li, Aslan,
  Reichman, Hybertsen, and Heinz]{Chernikov2014}
Chernikov,~A.; Berkelbach,~T.~C.; Hill,~H.~M.; Rigosi,~A.; Li,~Y.;
  Aslan,~O.~B.; Reichman,~D.~R.; Hybertsen,~M.~S.; Heinz,~T.~F. \emph{Phys.
  Rev. Lett.} \textbf{2014}, \emph{113}, 076802\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem[Wannier(1937)]{Wannier1937}
Wannier,~G.~H. \emph{Phys. Rev.} \textbf{1937}, \emph{52}, 191\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem[Berkelbach et~al.(2013)Berkelbach, Hybertsen, and
  Reichman]{Berkelbach2013}
Berkelbach,~T.~C.; Hybertsen,~M.~S.; Reichman,~D.~R. \emph{Phys. Rev. B}
  \textbf{2013}, \emph{88}, 045318\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\bibitem[Shi et~al.(2013)Shi, Pan, Zhang, and Yakobson]{Shi2013}
Shi,~H.; Pan,~H.; Zhang,~Y.-W.; Yakobson,~B.~I. \emph{Phys. Rev. B}
  \textbf{2013}, \emph{87}, 155304\relax
\mciteBstWouldAddEndPuncttrue
\mciteSetBstMidEndSepPunct{\mcitedefaultmidpunct}
{\mcitedefaultendpunct}{\mcitedefaultseppunct}\relax
\EndOfBibitem
\end{mcitethebibliography}
