First-principles quantum corrections for carrier correlations in double-layer two-dimensional heterostructures

TL;DR

This paper uses first-principles calculations to analyze quantum corrections to charge carrier correlations in two-dimensional heterostructures, revealing the importance of quantum effects in predicting electron-hole puddle behaviors.

Contribution

It provides the first unambiguous first-principles quantum corrections to Thomas-Fermi densities in 2D heterostructures using density-potential functional theory.

Findings

Quantum corrections significantly affect puddle correlations.

Strain potential must exceed impurity potential by tenfold for anti-correlation.

Quantum effects are crucial for matching experimental energy distributions.

Abstract

We present systematic ab initio calculations of the charge carrier correlations between adjacent layers of two-dimensional materials in the presence of both charged impurity and strain disorder potentials using the examples of monolayer and bilayer graphene. For the first time, our analysis yields unambiguous first-principles quantum corrections to the Thomas--Fermi densities for interacting two-dimensional systems described by orbital-free density functional theory. Specifically, using density-potential functional theory, we find that quantum corrections to the quasi-classical Thomas-Fermi approximation have to be taken into account even for heterostructures of mesoscopic size. In order for the disorder-induced puddles of electrons and holes to be anti-correlated at zero average carrier density for both layers, the strength of the strain potential has to exceed that of the impurity potential by at least a factor of ten, with this number increasing for smaller impurity densities. Furthermore, our results show that quantum corrections have a larger impact on puddle correlations than exchange does, and they are necessary for properly predicting the experimentally observed Gaussian energy distribution at charge neutrality.

Tables

Table 1. Table 1: Thomas–Fermi density n TF ​ ( \mathbold ​ r ) subscript 𝑛 TF \mathbold 𝑟 n_{\mathrm{TF}}(\mathbold{r}) and quantum-corrected density n 3 ​ ( \mathbold ​ r ) subscript 𝑛 3 \mathbold 𝑟 n_{3}(\mathbold{r}) for linear and quadratic dispersion in 2D from evaluating Eq. ( 7 ), with degeneracy factor g 𝑔 g , Bessel functions J α ​ ( ) subscript J 𝛼 \mathrm{J}_{\alpha}(\,) , ν ​ ( \mathbold ​ z ) = [ μ − V ​ ( \mathbold ​ z ) ] + 𝜈 \mathbold 𝑧 subscript delimited-[] 𝜇 𝑉 \mathbold 𝑧 {\nu(\mathbold{z})=[\mu-V(\mathbold{z})]_{+}} , K ​ ( \mathbold ​ z ) = ν ​ ( \mathbold ​ z ) / ( ℏ ​ v F ) 𝐾 \mathbold 𝑧 𝜈 \mathbold 𝑧 Planck-constant-over-2-pi subscript 𝑣 F {K(\mathbold{z})=\nu(\mathbold{z})/(\hbar v_{\mathrm{F}})} , σ ​ ( \mathbold ​ z ) = 2 ​ m ​ ν ​ ( \mathbold ​ z ) / ℏ 2 𝜎 \mathbold 𝑧 2 𝑚 𝜈 \mathbold 𝑧 superscript Planck-constant-over-2-pi 2 {\sigma(\mathbold{z})=2m\,\nu(\mathbold{z})/\hbar^{2}} , and [ x ] + subscript delimited-[] 𝑥 [x]_{+} denoting x ​ η ​ ( x ) 𝑥 𝜂 𝑥 x\,\eta(x) . We recover the TF densities from their quantum-corrected successors upon replacing K ​ ( \mathbold ​ r + \mathbold ​ s ) 𝐾 \mathbold 𝑟 \mathbold 𝑠 {K(\mathbold{r}+\mathbold{s})} and σ ​ ( \mathbold ​ r + \mathbold ​ s ) 𝜎 \mathbold 𝑟 \mathbold 𝑠 {\sigma(\mathbold{r}+\mathbold{s})} by their local versions K ​ ( \mathbold ​ r ) 𝐾 \mathbold 𝑟 K(\mathbold{r}) and σ ​ ( \mathbold ​ r ) 𝜎 \mathbold 𝑟 \sigma(\mathbold{r}) , respectively. The case of quadratic dispersion is dealt with in Ref. [ Chau  et al. , 2018 ]; see Appendix  F for the derivations in the case of linear dispersion.

$T\left(\text{mathbold}p\right)$	$n_{\mathrm{TF}}\left(\text{mathbold}r\right)$	$n_3\left(\text{mathbold}r\right)$
$v_F\left|\text{mathbold}p\right|$	$\frac{g}{4\pi }K{\left(\text{mathbold}r\right)}^2$	$\frac{g}{4{\pi }^2}\int \left(d\text{mathbold}s\right)\frac{K{\left(\text{mathbold}r+\text{mathbold}s\right)}^3}{\left|\text{mathbold}s\right|}{J}_1\left(2\left|\text{mathbold}s\right|K\left(\text{mathbold}r+\text{mathbold}s\right)\right)$
$\frac{\text{mathbold}{p}^2}{2m}$	$\frac{g}{4\pi }\sigma \left(\text{mathbold}r\right)$	$\frac{g}{4{\pi }^2}\int \left(d\text{mathbold}s\right)\frac{\sigma \left(\text{mathbold}r+\text{mathbold}s\right)}{{\left|\text{mathbold}s\right|}^2}{J}_2\left(2\left|\text{mathbold}s\right|\sqrt{\sigma \left(\text{mathbold}r+\text{mathbold}s\right)}\right)$

Table 2. Table 2: The correlation measure ξ 𝜉 \xi , given in Eq. ( 18 ), as a function of the ratio R 𝑅 R between strain and disorder strength (cf. Eq. ( 16 )) for one disorder realization of double monolayer graphene in TF approximation, with and without exchange (X). In absolute numbers on the relevant scale of [ − 1 , 1 ] 1 1 [-1,1] for the correlation measures and regarding the spread of ξ 𝜉 \xi due to different disorder realizations, cf. Fig.  5 , we find exchange effects to be insignificant for investigating correlations between the two layers L1 and L2. We made the same observation for different disorder realizations.

$R$	0.1	1	10	20	50
${\xi }^{\mathrm{TF}+X}$	0.492	0.347	-0.013	-0.051	-0.122
${\xi }^{\mathrm{TF}}$	0.480	0.341	-0.015	-0.058	-0.135

Equations

$$E[n,\mu ]$$

$$+\mu \left(N-\int (\mathrm{d}\text{mathbold}r)n(\text{mathbold}r)\right),$$

$$E[V,n,\mu ]$$

$$+E_{\mathrm{int}}[n]+\mu N.$$

$$n(\text{mathbold}r)$$

$$V(\text{mathbold}r)$$

$$N$$

$$E_1[V-\mu ]=\mathrm{tr}\{(H_1-\mu )\eta (\mu -H_1)\},$$

$$\displaystyle n(\mathbold{r})=g\int\limits_{\begin{picture}(16.0,3.0)(-8.0,-3.0)\put(0.0,0.0){\curve(-3,0,-8,0)\curve(3,0,8,0)}\put(8.0,0.0){\curve(0,0,-1.5,1.5)\curve(0,0,-1.5,-1.5)}\put(0.0,0.0){\arc(-3,0){180}}\put(0.0,0.0){\makebox(0.0,0.0)[]{$\cdot$}}\end{picture}}\frac{\mathrm{d}t}{2\pi\mathrm{i}t}\,\mathrm{e}^{\mbox{\footnotesize$\frac{\mathrm{i}t}{\hbar}\mu$}}\langle\mathbold{r}|\mathrm{e}^{\mbox{\footnotesize$-\frac{\mathrm{i}t}{\hbar}H_{1}$}}|\mathbold{r}\rangle,$$

$$V_{\mathrm{ext}}(\text{mathbold}r)=V_{\mathrm{C}}(\text{mathbold}r)+V_{\mathrm{S}}(\text{mathbold}r)$$

$$V_{+}(\text{mathbold}r)=2(\mu +\Delta )-V_{-}(\text{mathbold}r)$$

$$E_{\mathrm{H}}[n]=\frac{W}{2}\int (\mathrm{d}\text{mathbold}r)(\mathrm{d}\text{mathbold}{r}^{\prime })\frac{n(\text{mathbold}r)n(\text{mathbold}{r}^{\prime })}{\max (|\text{mathbold}r-\text{mathbold}{r}^{\prime }|,b)},$$

$$V_{\mathrm{H}}(\text{mathbold}r)=W\int (\mathrm{d}\text{mathbold}r)\frac{n(\text{mathbold}{r}^{\prime })}{\max (|\text{mathbold}r-\text{mathbold}{r}^{\prime }|,b)}$$

$$W=\frac{{\mathrm{e}}_{\mathrm{cgs}}^2}{\kappa {ϵ}_0}=\hslash v_{\mathrm{F}}{r}_{\mathrm{s}}\approx \frac{2.187}{\kappa }ul$$

$$W\approx \frac{0.1596}{\kappa }{u}^{\prime }l$$

$$V_{\mathrm{X}}(\text{mathbold}r)$$

$$\times \left[\frac{1}{4}\log \left(\frac{4\sqrt{3}}{\sqrt{\pi |n(\text{mathbold}r)|}}\right)-0.5757\right]$$

$$V_{\mathrm{C}}(\text{mathbold}r)=W\int (\mathrm{d}\text{mathbold}{r}^{\prime })\frac{C(\text{mathbold}{r}^{\prime })}{\sqrt{|\text{mathbold}r-\text{mathbold}{r}^{\prime }{|}^2+d^2}},$$

$$V_{\mathrm{S}}(\text{mathbold}r)=RW\int (\mathrm{d}\text{mathbold}{r}^{\prime })\frac{\tilde{C}(\text{mathbold}{r}^{\prime })}{\sqrt{|\text{mathbold}r-\text{mathbold}{r}^{\prime }{|}^2+d^2}},$$

$$V_{\mathrm{ext}}(L2)(\text{mathbold}r)={V_{\mathrm{C}}(\text{mathbold}r)|}_{d=\tilde{d}}+\int (\mathrm{d}\text{mathbold}{r}^{\prime })\frac{W{n}_1(\text{mathbold}r+\text{mathbold}{r}^{\prime })}{\sqrt{|\text{mathbold}{r}^{\prime }{|}^2+(\tilde{d}-d{)}^2}}.$$

$$\xi [n_1,n_2]=\frac{⟨(n_1-⟨n_1⟩)(n_2-⟨n_2⟩)⟩}{n_{1,\mathrm{rms}}{n}_{2,\mathrm{rms}}}$$

$${\xi }_c[n_1,n_2]=\frac{⟨(n_1-⟨n_1{⟩}_c)(n_2-⟨n_2{⟩}_c){⟩}_c}{n_{1,\mathrm{rms},c}{n}_{2,\mathrm{rms},c}},$$

$$⟨A{⟩}_c=\frac{{\sum }_{i=1}^{\Lambda }{\chi }_c(i)A(i)}{{\sum }_{i=1}^{\Lambda }{\chi }_c(i)}$$

$$V_{-}^{(i+1)}=V_{\mathrm{ext}}+{\frac{\delta E_{\mathrm{int}}[n]}{\delta n}|}_{n_{-}=n_{-}^{(i)}}$$

$$||n^{(i+1)}-n^{(i)}||$$

$${\chi }_{ij}$$

$$t^{(i+1)}(\text{mathbold}{r}_j)$$

$$s^{(i+1)}$$

$${\mathrm{MAX}}_i$$

$${\mathrm{MIN}}_i$$

$$n(\text{mathbold}r)$$

$$⟨\text{mathbold}r|U_{\nu }|\text{mathbold}r⟩$$

$$=:\prod _{i=1}^{\nu /2}\int (\mathrm{d}\text{mathbold}{r}_i)(\mathrm{d}\text{mathbold}{p}_i)c_i{\mathrm{e}}^{-\frac{\mathrm{i}t}{\hslash }A}$$

$$⟨\text{mathbold}r|\eta (\mu -H_1)|\text{mathbold}r⟩$$

$$n_{\mathrm{TF}}(\text{mathbold}r)$$

Full text

PhySH: Two-dimensional electron system, Density functional approximations, Heterostructures

I Introduction

The simulation of two-dimensional (2D) materials and prediction of their properties has become a mainstay of materials science over the past decade, with the promise and realization of valuable applications in both industrial technology and fundamental research Geim and Grigorieva (2013). Theoretical and computational methods for 2D materials have been advanced into a sophisticated machinery that enables researchers to deal with ever more realistic settings Rodriguez-Vega et al. (2014). The widely used Kohn-Sham density functional theory (KS-DFT) Kohn and Sham (1965); Dreizler and Gross (1990) presents one particularly popular ab initio approach with the capability of accurately handling hundreds of interacting particles (and up to thousands of atoms in cases where linear-scaling methods apply Goedecker (1999); Soler et al. (2002); Bowler and Miyazaki (2012); Fox et al. (2014); Aarons et al. (2016); Cole and Hine (2016)). The development of functionals for 2D systems has been lagging behind that of their 3D counterparts for various reasons: Some of the most heavily used 3D KS-DFT functionals are not bounded from below Chiodo et al. (2012); Kaplan et al. (2018) when the 2D limit is approached. This stems from the improper scaling behaviour of the density Chiodo et al. (2012). Furthermore, even consistent first-order gradient corrections of the kinetic energy density (used in developing meta-generalized gradient approximations for KS-DFT) were unknown for 2D fermion systems until recently Holas et al. (1991); Brack and Bhaduri (2003); Trappe et al. (2016). Nonetheless, considerable progress has been made and alternative derivations of some of those 2D functionals have been obtained since the turn of the millenium Pollack and Perdew (2000); Reimann and Manninen (2002); Pittalis et al. (2007, 2008); Constantin et al. (2008); Pittalis et al. (2009a, b); Pittalis and Räsänen (2010); Räsänen et al. (2010); Vilhena et al. (2014).

However, a systematic ab inito methodology that is universally applicable and scales favorably with particle number, thereby enabling high-throughput computations of mesoscopic systems, is not yet available — with orbital-free density functional theory (OF-DFT) being the suspected saviour for almost a century Thomas (1927); Fermi (1927); Cangi et al. (2011); Burke (2012); Karasiev et al. (2013); Pribram-Jones et al. (2015); Witt et al. (2018). While KS-DFT scales cubically with particle number in generic settings, OF-DFT scales linearly, and sub-linear scaling can be achieved in special cases. Functional development, in particular concerning the kinetic energy functional Wang and Teter (1992); Della Sala et al. (2015); Xia and Carter (2015), and implementations of OF-DFT have gained momentum in recent years Karasiev et al. (2012, 2013); Chen et al. (2015); Das et al. (2015); Chen et al. (2016); Mi et al. (2016); Constantin et al. (2018); Witt et al. (2018), also in conjunction with other techniques like ab initio molecular dymanics González and González (2008). Chemical accuracy is approached in selected cases Xia et al. (2012); Borgoo et al. (2014); Espinosa Leal et al. (2015). If quantum effects play a minor role or if the considered system is largely homogeneous, OF-DFT can also be used in its most basic form, the Thomas–Fermi (TF) approximation. For instance, the effect of exchange on large disordered systems with long-range interactions was studied in Ref. [Rossi and Das Sarma, 2008] using OF-DFT in TF approximation. The TF model is not only of historical significance, but presents, as an exact constraint for homogeneous systems and in the limit of infinite nuclear charges Lieb and Siman (1977), an important base line for benchmarking proposed systematic density functional improvements. To what extent then do corrections to the TF approach play a crucial role or dominate over exchange effects in 2D materials? (They do, indeed, for a number of relevant fermionic systems, ranging from atomic Fermi gases to molecules and single atoms.)

The most severe obstacle for OF-DFT in taking over as the workhorse of theoretical chemistry and materials science is the lack of accurate, reliable, systematic, and preferably universal quantum corrections to the quasi-classical TF approximation, in particular for the kinetic energy of low-dimensional systems Holas et al. (1991); Brack and Bhaduri (2003); Salasnich (2007); Trappe et al. (2016, 2017). While ad-hoc corrections to the quasi-classical limit and heuristic approximations are available for kinetic energy and particle density of low-dimensional systems van Zyl et al. (2013, 2014), successful derivations of systematic and consistent corrections are scarce Ribeiro et al. (2015); Trappe et al. (2015, 2016, 2017); Chau et al. (2018). One promising route towards systematic orbital-free quantum corrections is provided by density-potential functional theory (DPFT) Englert (1988, 1992); Cinal and Englert (1993); Trappe et al. (2016, 2017); Chau et al. (2018), a more flexible reformulation of the original Hohenberg-Kohn DFT Hohenberg and Kohn (1964); Dreizler and Gross (1990), which circumvents the need for an explicit kinetic-energy density-functional and provides natural ways for systematic semiclassical expansions.

In this article we explore the applicability of DPFT for 2D materials by assessing quantum-corrections to the TF approximation for double-layer heterostructures of mono- and bilayer graphene. Of particular interest to us are situations that are not easily tackled with orbital-based techniques, for example 2D material sheets of mesoscopic size that are subjected to aperiodic disorder potentials. Such situations are for example of current interest in studies on Coulomb drag Narozhny and Levchenko (2016) where there exists an unsettled controversy as to whether the behavior of drag measured in experiment Gorbachev et al. (2012) is due to correlation Song and Levitov (2012) or anti-correlation Ho et al. (2018) between the density fluctuations of the layers.

Our work contributes in several ways to answering some of the questions raised above. Sections II and III provide the computational framework for obtaining quantum-corrected carrier-densities of 2D materials using DPFT. The expressions for the semiclassical particle densities developed here enable us to decide whether or not the quasi-classical TF approximation is sufficient for describing at least conglomerate properties like average inter-layer correlations of heterostructures. Section IV introduces the generic double-layer system, with both layers subjected to one layer of charged impurities, while only one of the layers is strained. Charge and strain disorder potentials are expected to compete in creating correlated (from charged impurities) and anticorrelated (from strain) carrier densities in the two layers. Our model setup is designed to extract the strain strengths required for switching between correlation and anticorrelation. In Sec. V we apply our new approach to double-monolayer graphene and double-bilayer graphene. We discuss whether or not the electron-hole puddles of both layers, interacting electrostatically, require a self-consistent inter-layer treatment. Finally, we analyse the effects of quantum-corrections and exchange energy on the correlations with the aid of phase diagrams that chart the correlation measures as functions of impurity density, carrier density, and ratios of strain and charge disorder. The appendix gathers background information on the units, system parameters, correlation measures, and numerical procedures employed here.

II Density-potential functional theory

Instead of resorting to the computationally demanding orbital-based Kohn-Sham DFT, we make use of orbital-free density-potential functional theory (DPFT) Englert (1988, 1992). It is formally equivalent to the Hohenberg-Kohn formulation, but makes systematic improvements upon the TF approximation technically feasible — in particular for low-dimensional systems.

Specifically, by Legendre-transforming $E_{\mathrm{kin}}[n]$, the kinetic energy functional of the particle density $n(\mathbold{r})$, w.r.t. the new variable ${V(\mathbold{r})=\mu-\frac{\delta E_{\mathrm{kin}}[n]}{\delta n(\mathbold{r})}}$, we recast the total energy of an interacting quantum system with interaction energy $E_{\mathrm{int}}[n]$,

[TABLE]

as the density-potential functional

[TABLE]

Here, the external potential $V_{\mathrm{ext}}(\mathbold{r})$ yields the external energy $E_{\mathrm{ext}}[n]$, and the particle number $N$ is enforced via the Lagrange multiplier $\mu$, viz. the chemical potential. From Eq. (II) we obtain the ground-state solutions of the three variables $V$, $n$, and $\mu$ by self-consistently solving

[TABLE]

Equation (5) is obtained by combining $\partial E[V,n,\mu]/\partial\mu$ with Eq. (3) and reveals the particle number constraint in Eq. (II).

Equations (3)–(5) are exact and reminiscent of the KS scheme, but without the need of orbitals. However, the noninteracting case aside, we have to approximate the unknown potential functional ${E_{1}[V-\mu]}$; here the subscript indicates that $E_{1}$ can be written in terms of a single-particle trace over a function of the single-particle Hamilton operator Englert (1992). We also have to provide the equally important interaction energy $E_{\mathrm{int}}[n]$ as an explicit functional of the particle density $n$. Approximate particle densities follow directly from approximations of ${E_{1}[V-\mu]}$ (or, rather, its functional derivative) for any given potential $V$. As is evident from Eq. (4), $V$ constitutes an effective single-particle potential with interaction effects effectively included for any given density $n$

Following Refs. [Englert and Schwinger, 1984; Englert, 1992; Cinal and Englert, 1993; Trappe et al., 2016, 2017; Chau et al., 2018] we approximate $E_{1}$ by its noninteracting version as the single-particle trace

[TABLE]

where ${H_{1}=H_{1}(\mathbold{R},\mathbold{P})=T(\mathbold{P})+V(\mathbold{R})}$ is a single-particle Hamiltonian with dispersion relation $T$ and potential energy $V$, while the trace includes the degeneracy factor $g$. For example, ${g=4}$ accounts for the spin and valley multiplicity of unpolarized charge carriers in the cases of mono- and bilayer graphene. $\mathbold{R}$ and $\mathbold{P}$ are the position and momentum operators, respectively, and $\eta(\,)$ denotes the step function.

The explicit expression for $E_{1}$ in Eq. (6) results in an explicit expression for the particle density in terms of arbitrary functions $V(\mathbold{r})$ via Eq. (3). The approximate nature of Eq. (6) aside, the exact particle density including all quantum corrections is thereby obtained for any specified interaction energy $E_{\mathrm{int}}[n]$ and without reference to orbitals. Specifically, Eqs. (3) and (6), together with the Fourier transform of the step function, yield the particle density111The contour integration circumvents the singularity at ${t=0}$ in the lower half-plane.

[TABLE]

see Refs. [Golden, 1957, 1960; Light and Yuan, 1973; Lee and Light, 1975; Chau et al., 2018].

We seek to approximate the time evolution operator ${U=\mathrm{e}^{-\frac{\mathrm{i}t}{\hbar}H_{1}}}$ systematically via split-operator methods, for example of the Suzuki-Trotter type Hatano and Suzuki (2005); Chau et al. (2018). The quasi-classical approximation ${U\approx U_{2}=\mathrm{e}^{-\frac{\mathrm{i}t}{\hbar}T(\mathbold{P})}\mathrm{e}^{-\frac{\mathrm{i}t}{\hbar}V(\mathbold{R})}}$ recovers the quasi-classical TF density $n_{\mathrm{TF}}$, while ${U_{3}=\mathrm{e}^{-\frac{\mathrm{i}t}{2\hbar}T(\mathbold{P})}\mathrm{e}^{-\frac{\mathrm{i}t}{\hbar}V(\mathbold{R})}\mathrm{e}^{-\frac{\mathrm{i}t}{2\hbar}T(\mathbold{P})}}$ produces the first quantum-corrected density $n_{3}$ in a series of expressions that utilize higher-order factorizations Chau et al. (2018) of $U$. We give the corresponding 2D densities for linear and quadratic dispersion in Table 1 and outline the derivation of $n_{3}$ for the case of linear dispersion in Appendix F. In contrast to $n_{\mathrm{TF}}(\mathbold{r})$, which is restricted to classically allowed regions of the potential and only depends on the local value $V(\mathbold{r})$, $n_{3}(\mathbold{r})$ samples $V$ in an extended region and exhibits evanescent tails beyond the quantum-classical border.

III Self-consistent simulation of disordered 2D materials

We target 2D systems with chemical potential in the vicinity of the Dirac point (the point where valence and conduction bands touch) for graphene (effective bilayer graphene). The usual tight binding approach absorbs the lattice structure in an effective Hamiltonian and yields noninteracting quasiparticles in a homogeneous (${V_{\mathrm{ext}}(\mathbold{r})=0}$) environment. The electronic structures, viz. atoms, of the materials are thus not modeled explicitly here. The single-particle energies, viz. band structure, associated with such quasiparticles is the dispersion relation $T(\mathbold{p})$ whose operator version appears in Eq. (6). $T(\mathbold{p})$ can be an arbitrary function, but for the purpose of this work we shall restrict ourselves to the analytically more tractable cases of linear and quadratic dispersion222Hence, the scope of our investigation also extends to energies further away from the Dirac point (band gap) as long as $T(\mathbold{p})$ can be considered linear (quadratic)..

In the following we outline the procedures involved for arriving at the ground state solutions of Eqs. (3)–(5); further details are provided in Appendix E. Upon adding external potentials $V_{\mathrm{C}}(\mathbold{r})$ and $V_{\mathrm{S}}(\mathbold{r})$ that model charged impurites and strain, respectively, we initiate the self-consistent loop of Eqs. (3)–(5) by evaluating the density (denoted $n_{-}(\mathbold{r})$ for the quasiparticles that follow the dispersion of the conduction band) with the external potential

[TABLE]

for these conduction quasiparticles, see Appendix C for details. Since no interactions are included at this stage, the effective potential is ${V_{-}(\mathbold{r})=V_{\mathrm{ext}}(\mathbold{r})}$. The density $n_{+}(\mathbold{r})$ of valence quasiparticles, which follow the inverted dispersion, e.g. $T(\mathbold{p})=-v_{\mathrm{F}}|\mathbold{p}|$ in the case of graphene, is built from the same density expression as $n_{-}(\mathbold{r})$ but takes as an input the inverted potential $2\,(\mu+\Delta)-V_{\mathrm{ext}}(\mathbold{r})$, with an optional bandgap $\Delta$ (in case of mono- and bilayer graphene, we have ${\Delta=0}$). The effective potential for the valence quasiparticles reads

[TABLE]

and equals the external potential for the valence quasiparticles if interactions are omitted. The such obtained carrier density ${n(\mathbold{r})=n_{-}(\mathbold{r})-n_{+}(\mathbold{r})}$ updates the effective potential via the interaction contribution in Eq. (4).

As an approximate interaction energy $E_{\mathrm{int}}[n]$ for the quasiparticles, we employ the regularized Hartree term for the Coulomb energy,

[TABLE]

where $b$ is half the lattice constant of the numerical implementation and $W=\hbar v_{\mathrm{F}}r_{\mathrm{s}}$, with the ratio $r_{\mathrm{s}}$ of Coulomb potential energy and kinetic energy for graphene. Upon functional differentiation, Eq. (10) leads to the Hartree potential333Equation (11) can be regarded as a simplified version of the smooth models for the Coulomb potential discussed in Ref. [González-Espinoza et al., 2016].

[TABLE]

as an approximate interaction contribution in Eq. (4). Eyeing means of comparison and higher accuracy, we may supplement $E_{\mathrm{H}}$ with an exchange energy Dirac (1930), leading to the exchange potential $V_{\mathrm{X}}$; see Appendix B for details.

The updated effective potential $V$ then determines a new quasiparticle density $n_{-}(\mathbold{r})$ via Eq. (3), thereby closing the self-consistent loop. This process is repeated until a predefined relative precision is reached (we find $10^{-6}$ to be sufficient) when comparing the local densities of subsequent loop iterations. The chemical potential is adjusted in each iteration to enforce a given particle number, viz. average carrier density. Figure 1 highlights the differences in the converged quasiparticle densities $n_{\mathrm{TF}}$ and $n_{3}$ of a single graphene layer with and without exchange. Both exchange and quantum corrections tend to decrease the peaks of the density landscape. This effect is well-known in the case of exchange Rossi and Das Sarma (2008), and a smoothening of the carrier density in an external potential $V_{\mathrm{C}}$ is to be expected when tunneling starts to play a role with the inclusion of quantum corrections. In fact, when comparing the corresponding densities in Fig. 1 we find that quantum corrections considered in this work can dominate over exchange effects.

The analysis of disorder averages reveals a striking instance of this observation. As observed in Fig. 1 for a single disorder realization, $n_{3}$ and $n_{\mathrm{TF}}$ differ in their density distribution function. For $n_{3}$ the integration over a finite region in the disorder landscape, see Table 1, tends to result in smoother densities compared with $n_{\mathrm{TF}}$. This effect can be quantified by density histograms compiled from many disorder realizations. The density histograms in Fig. 2, calculated for a graphene monolayer on SiO2 with ${\bar{n}_{\mathrm{imp}}=10^{12}/\mathrm{cm}^{2}}$, corroborate the snapshot of one disorder realization in Fig. 1. Figures 1 and 2 highlight the (a priori) importance of not only addressing exchange but also quantum corrections beyond the TF approximation for quantitatively viable investigations of 2D materials via orbital-free DFT: Both the inclusion of exchange and of quantum corrections indicate less pronounced peaks of the carrier-density landscape. Compared to the TF approximation, the quantum corrections exhibit an averaging effect with broader density distributions that have relatively more weight on intermediate values of the density rather than a strong distribution maximum at zero density; cf. Fig. 2 (left). Although exchange also shows visibly less pronounced densities in Fig. 1, this effect stems from a global reduction in density variance rather than a redistribution of densities from very small towards intermediate values. The distributions of the quantum-corrected densities (including exchange) are captured by their Gaussian (‘G’) or Lorentzian (‘L’) fits (with offset) more accurately than the distributions resulting from the TF approximation. This observation is in line with experimental results for graphene that point towards Gaussian distributions of density and energy in the presence of charged disorder Martin et al. (2008); Xue et al. (2011). In Fig. 2 we compare our calculations ‘$n_{\mathrm{TF}}$(+X)’ and ‘$n_{3}$(+X)’, with their Gaussian fits ‘G(TF)’ and ‘G’, directly with the experimental data ‘exp.’ from scanning tunnelling spectroscopy; cf. Ref. [Xue et al., 2011]. We find our quantum-corrected approach to predict the experimental data much better than what can be obtained from the TF approximation — in both qualitative and quantitative terms. In view of this stark improvement over the TF approximation, we want to stress again that our quantum-corrected density expressions are based on first principles without adjustable parameters or fits, and rely solely on controlled approximations to quantum mechanics.

IV Double layer setup

The treatment of monolayers in the previous section forms our basis for the description of more complicated heterostructures. Figure 3 illustrates a two-layer system, where both layers L1 and L2 are sandwiched between h-BN and subjected to a charged impurity layer from the SiO2 substrate. We expose L2 to the same charge disorder that affects L1, though at a larger separation, but refrain from adding disorder on L2 in order to avoid confusing inter-layer correlation effects with effects from independent disorder on L2. For the same reason we model the strain of L1 and the charge disorder by the identical type of disorder, albeit in different realizations. The layer separation of $5\,$nm suffices to justify a merely classical electrostatic interaction between L1 and L2, i.e., inter-layer tunneling of charge carriers can be neglected — in contrast to intra-layer tunneling through the disorder potential landscape. The latter is missed by densities in TF approximation but captured (in part) via the higher-order Suzuki-Trotter factorizations.

In what follows we address the puddle correlations between L1 and L2 as a function of the ratio $R$ between the disorder strength of the strain and that of the charged impurities. Further details are provided in Appendix C. The sum of $V_{\mathrm{C}}$ and $V_{\mathrm{S}}$ results in electron-hole puddles within the first layer L1, whose electrostatic potential adds to the external potential $V_{\mathrm{C}}$ for the charge carriers in the second layer444Strictly speaking, the electrostatic effects of L1 and L2 on each other have to be treated self-consistently, but the layer separation of 5nm ensures that the backreaction of L2 on L1 is of the order of $10^{-6}$ of the external potential of L1 and therefore negligible at the level of precision we are aiming for in this work.. We expect maximal puddle correlation if ${R=0}$, that is, when no strain can obscure the then dominating effect of the charged impurities on both layers: Owing to the intra-layer Coulomb interaction, the puddles in L1 are much reduced in weight compared with the case of noninteracting carriers. That is, the tendency of a puddle in L1 to electrostatically induce a puddle of opposite charge in L2 is overcompensated by the charge disorder, which exhibits the tendency to induce a puddle of the same charge555In the (unphysical) limit of vanishing layer separation these tendencies become certainties, but a finite layer separation in concert with the disordered potential landscape allows for local variations of correlation, thereby affecting global correlation measures quantitatively.. Following the same line of reasoning, we expect maximal anticorrelation if ${R=\infty}$, with the transition from correlation to anticorrelation occuring at some value ${R>1}$. In the following section we substantiate these claims with quantitative predictions for graphene and bilayer graphene.

V Density correlations in double-layers of mono- and bilayer graphene

We quantify the inter-layer correlations of electron-hole-puddles of the double-layer system described in Sec. IV by solving Eqs. (3) and (4) self-consistently666During the self-consistent loop the chemical potential $\mu$ is adjusted such that the mean carrier densities ${\left<\right.{n_{j}}\left.\right>}$ are kept at a fixed value (zero, unless stated otherwise). and by comparing the converged carrier densities ${n_{1}=n(\mathrm{L1})}$ and ${n_{2}=n(\mathrm{L2})}$ of layers L1 and L2 locally. To that end we calculate the two correlation measures $\xi[n_{1},n_{2}]$ and $\xi_{c}[n_{1},n_{2}]$, which yield a value of one for perfectly correlated electron-hole puddles (i.e., if the density distribution $n_{2}$ is proportional to $n_{1}$ and their values have the same sign at each position $\mathbold{r}$), minus one for perfect anticorrelation (i.e., if $n_{2}$ is proportional to $-n_{1}$), and are designed for tracking the transition between these two extremes; see Appendix D for details.

Figure 4 depicts potentials and densities for graphene, viz. linear dispersion, calculated for mean carrier densities ${\left<\right.{n_{j}}\left.\right>=0}$ and equal strengths of strain and charge disorder (ratio ${R=1}$). Due to the screening effects of the Coulomb interaction within L1, the effective potential for L1 exhibits less variability then the total external potential, i.e., the sum of the disorder potentials $V_{\mathrm{C}}$ and $V_{\mathrm{S}}$. The quasiparticle density $n_{1}$ of L1 , which can be viewed as resulting from this effective potential, induces an external electrostatic potential for the carrier density $n_{2}$ of L2. However, for ${R=1}$ the charge disorder potential dominates the total external potential for $n_{2}$ with magnitudes by a factor of more than 50 larger than those of the electrostatic potential caused by $n_{1}$. For the setting that leads to Fig. 4, the magnitudes of the total external potential for L2 are smaller than those of L1 by a factor of 3–5.

As a result, the density fluctuations of $n_{2}$ are diminished compared with those of $n_{1}$ by a factor of 5–10. It is therefore well justified to refrain from a self-consistent treatment of the electrostatically induced potentials $V_{\mathrm{ind}}$ of both layers, and to consider only $V_{\mathrm{ind}}$ stemming from $n_{1}$. As is evident from the bottom row of Fig. 4, the spatial distributions of $n_{1}$ and $n_{2}$ are correlated rather than anticorrelated for ${R=1}$.

Repeating the calculation which yields the results illustrated in Fig. 4 for different values of $R$, we find the critical value $R_{0}$ at the transition from inter-layer correlation to anticorrelation. The correlation measure $\xi$ signifies inter-layer puddle-correlation (anticorrelation) by taking on positive (negative) values. The corresponding diagram in Fig. 5, calculated with charged impurities density $\bar{n}_{\mathrm{imp}}=10^{12}\,\mathrm{cm}^{-2}$, exhibits ${R_{0}\approx 10\pm 2}$ when $n_{3}$ is used, and a similar value in the case of $n_{\mathrm{TF}}$. Here, we take into account exchange effects and report a rough error estimate simply based on our numerical findings from five disorder realizations777At this point we do not calculate disorder averages, but rather showcase the predictions of a few disorder realizations (also bearing in mind the limited number of samples usually investigated in actual experiments).. Evidently, the disparity in electron-hole puddle landscapes between $n_{3}$ and $n_{\mathrm{TF}}$ as seen in Fig. 1 does not translate into an appreciable difference between $\xi^{(3)+\mathrm{X}}$ (for $n_{3}$) and $\xi^{\mathrm{TF+X}}$ (for $n_{\mathrm{TF}}$).

Although it is not surprising per se that integrated quantities like $\xi$ are less sensitive to local differences between $n_{3}$ and $n_{\mathrm{TF}}$, it cannot be assumed a priori. Our quantitative analysis shows that a quasi-classical approach to inter-layer correlations in monolayer-graphene heterostructures is justified in case of rather large impurity densities like of $10^{12}\,\mathrm{cm}^{-2}$ used for Fig. 5. Commonly, however, the TF approximation is less reliable for smaller particle numbers, and quantum corrections can be expected to play a more dominant role as the carrier density is reduced. The disorder potential for smaller $\bar{n}_{\mathrm{imp}}$ is less pronounced and gives rise to puddles that exhibit smaller carrier densities on average. Indeed, with ${\bar{n}_{\mathrm{imp}}=10^{11}/\mathrm{cm}^{2}}$ employed for Fig. 6, the quantum-corrected $\xi^{(3)+\mathrm{X}}$ can be clearly distinguished from the quasi-classical $\xi^{\mathrm{TF+X}}$. Our data shown in Fig. 6 point to the critical value ${R_{0}^{\mathrm{TF+X}}\approx 15\pm 2}$ with an error estimate similar to that in Fig. 5. An increased $R_{0}$ at lower $\bar{n}_{\mathrm{imp}}$ can be understood from the following simplified picture. With typical values $\left.v_{\mathrm{C}}\right|_{d}$ and $v_{\mathrm{S}}$ of the charge and strain potentials for L1 and neglecting interactions, we have typical values {n_{\mathrm{TF}}(\mathrm{L1})\propto\big{(}\mu-v_{\mathrm{ext}}(\mathrm{L1})\big{)}^{2}\approx(\left.v_{\mathrm{C}}\right|_{d}+v_{\mathrm{S}})^{2}} for ${\mu\approx 0}$, i.e., approximately ${n_{\mathrm{TF}}(\mathrm{L1})\propto v_{\mathrm{S}}^{2}}$ at ${R\approx 10}$. Then, the typical TF density in L2 is determined by ${v_{\mathrm{ext}}(\mathrm{L2})=\left.v_{\mathrm{C}}\right|_{\tilde{d}}+v_{\mathrm{ind}}}$, where the typical values ${\left.v_{\mathrm{C}}\right|_{\tilde{d}}\approx\frac{1}{2}\left.v_{\mathrm{C}}\right|_{d}}$ at ${\tilde{d}-d\approx 5}\,$nm and ${v_{\mathrm{ind}}\propto v_{\mathrm{S}}^{2}}$ are roughly equal if ${\bar{n}_{\mathrm{imp}}=10^{12}/\mathrm{cm}^{2}}$ and ${R=10}$. For a scaled ${\bar{n}_{\mathrm{imp}}=\lambda\times 10^{12}/\mathrm{cm}^{2}}$ and the same ${R=10}$, $v_{\mathrm{ext}}(\mathrm{L1})$ scales with $\lambda$ as well, but the typical values of ${n_{\mathrm{TF}}(\mathrm{L1})}$ then scale like $\lambda^{2}\,v_{\mathrm{S}}^{2}$, such that ${v_{\mathrm{ext}}(\mathrm{L2})=\lambda\left.v_{\mathrm{C}}\right|_{\tilde{d}}+\lambda^{2}\,v_{\mathrm{ind}}}$. For the case of ${\lambda=1/10}$, as represented by Fig. 6, the strength of $v_{\mathrm{ind}}$, relative to $\left.v_{\mathrm{C}}\right|_{\tilde{d}}$, is diminished by a factor of ten. As strain feeds into $v_{\mathrm{ind}}$, not into $\left.v_{\mathrm{C}}\right|_{\tilde{d}}$, (relatively) more strain is required to counteract the effect of $\left.v_{\mathrm{C}}\right|_{\tilde{d}}$, implying a larger critical $R_{0}$ at lower impurity densities. Disregarding the uncertainties for the critical $R$, one could estimate ${R_{0}^{(3)+\mathrm{X}}\approx 30}$ for the quantum-corrected correlations. However, $\xi^{(3)+\mathrm{X}}$ rather exhibits an extended cross-over regime, where the magnitude of strain can be varied substantially with little effect on the average correlation $\xi^{(3)+\mathrm{X}}$. Our main result is thus that corrections to the TF approximation can become important even for integrated or averaged observables of 2D materials — given the proper conditions, for example, small carrier densities.

We briefly turn to the impact of exchange effects in inter-layer correlations. Since quantum corrections modify the TF carrier distributions of the individual layers more profoundly than exchange does, cf. Fig. 1, we expect that exchange plays a minor role in determining carrier correlations of double-layer systems. This is confirmed with Table 2 in Appendix B, where $\xi$ is determined in TF approximation for a single disorder realization, with and without exchange. In view of the magnitude of uncertainties displayed in Figs. 5 and 6 we find negligible quantitative differences when omitting or including exchange.

Our results on double bilayer graphene in Fig. 7 show a slight quantitative difference between the transitions of $\xi^{(3)}$ and $\xi^{\mathrm{TF}}$ for the here employed density of charged impurities ${\bar{n}_{\mathrm{imp}}=10^{12}\,\mathrm{cm}^{-2}}$. However, the overall trend for both approximations is still very similar and the spread of $\xi$, originating in multiple disorder realizations like in Figs. 5 and 6, still implies ${R_{0}^{(3)}\approx R_{0}^{\mathrm{TF}}}$ within the depicted uncertainties. Extrapolating from the double-monolayer case the double-bilayer graphene can be expected to exhibit more pronounced differences between $\xi^{(3)}$ and $\xi^{\mathrm{TF}}$ for smaller density of charged impurities.

Figure 8 provides another angle on the transition from correlation to anticorrelation for double bilayer graphene: For fixed ${\left<\right.{n_{1}}\left.\right>=0}$ we chart $\xi^{\mathrm{TF}}$ as a function of $\left<\right.{n_{2}}\left.\right>$ and find a qualitative change in the shape of the curves $\xi(\left<\right.{n_{2}}\left.\right>)$ as $R$ transits across ${R_{0}^{\mathrm{TF}}\approx 8}$, consistent with our findings in Fig. 7.

VI Conclusions and perspectives

The main result of this work is thus that both exchange and quantum corrections beyond the Thomas–Fermi approximation are essential for quantitatively reliable density distributions in the context of 2D materials. While exchange plays a minor role for some integrated quantities like inter-layer correlations of double-layer systems, quantum corrections become important for smaller average carrier-densities, viz. cleaner systems. In other words, systems like graphene on boron nitride samples with ${\bar{n}_{\mathrm{imp}}\approx 10^{10}\mathrm{cm}^{-2}}$ call for a more sophisticated description than the Thomas–Fermi approximation can provide. We showed that the regime of high carrier densities may still be described quasiclassically in the computationally highly efficient Thomas–Fermi approximation. For specific systems, however, quantum corrections have to be taken into account at the quantitative level. In this work we provided one such example when analyzing the crossover from correlation to anticorrelation in bilayer heterostructures of 2D materials. We focussed on describing double-layers of mono- and bilayer graphene via recently developed techniques of calculating quantum corrections, but our approach can be easily adapted for other heterostructures.

Turning to the issue of Coulomb drag Gorbachev et al. (2012), our analysis reveals that the strain potential has to be at least ten times stronger than the impurity potential in order for anti-correlation between the layers to occur in both monolayer and bilayer graphene. Any strain potential weaker than this will result in correlation. The DPFT framework presented here thus provides the following path towards an experimental verification of the nature of correlations in graphene heterostructures, a non-trivial task since it is impossible to directly measure the local densities of encapsulated 2D layers via local scanning tunneling microscopy. The strength of the strain potential Gibertini et al. (2012) may first be extracted from local height fluctuation data obtained via atomic force microscopy Dean et al. (2010), followed by the impurity strength from conductivity measurements Adam et al. (2007). The resulting ratio $R$ may then be compared against Figs. 5 to 7 to determine the nature of correlation between the layers. This would be especially useful for clarifying the source of the unexpected sign changes Li et al. (2016); Lee et al. (2016); Simonet et al. (2017) in Coulomb drag experiments that have been explained by earlier works simply asserting correlation Song and Levitov (2012); Schütt et al. (2013) or anti-correlation Gorbachev et al. (2012); Ho et al. (2018) without proof. In the case of bilayer graphene, Ref. [Zarenia et al., 2018] theoretically demonstrated that a multiband mechanism based on the thermal smearing of Fermi energy in both layers explains the unexpected sign changes at high densities where puddles may be ignored. A complete analysis at the double neutrality where both puddles and multiband effects are important is however still lacking. In future, this may be addressed by generalizing the DPFT framework presented here to incorporate multiband effects Zarenia et al. (2018) in the determination of puddle-induced quantities such as the interlayer correlation and density fluctuations.

Finally, the findings of this work are useful in the study of van der Waals heterostructures Geim and Grigorieva (2013) consisting of a pair of two-dimensional electronic layers separated by thin dielectric spacers. These structures serve as ideal platforms for the study of a range of interesting physical effects such as exciton condensation Li et al. (2017); Liu et al. (2017); Burg et al. (2018); López Ríos et al. (2018) and strong light-matter interaction Britnell et al. (2013); Yu et al. (2013); Withers et al. (2015) and are thus an area of intense research activity. It is likely that the inter-layer puddle correlations calculated within this work plays a role in the above physics.

Acknowledgements.

We acknowledge the support of the National Research Foundation of Singapore under its Fellowship program (NRF-NRFF2012-01) and the National University of Singapore Young Investigator Award (R-607-000-094-133). We also acknowledge the use of the dedicated computational facilities at the Centre for Advanced 2D Materials.

Appendix A Units and system parameters

Throughout this work we use cgs units, measure lengths in ${l=\sqrt{3}a\approx 2.46\,\mbox{\r{A}}}$, where $a$ is the lattice constant of graphene, energies in ${u=\hbar v_{\mathrm{F}}/l\approx 2.678\,\mathrm{eV}}$ (for monolayer graphene, with the Fermi velocity ${v_{\mathrm{F}}=10^{6}}\,$m/s) and ${u^{\prime}=\hbar^{2}/(m\,l^{2})\approx 36.71}\,$eV (for bilayer graphene with quasiparticle mass Raza ${m\approx 0.39\,\mathrm{eV}/(2\,v_{\mathrm{F}}^{2})}$). With the elementary charges in cgs- and SI-units connected via ${\mathrm{e}_{\mathrm{SI}}=\sqrt{4\pi}\mathrm{e}_{\mathrm{cgs}}}$ and the dielectric constants related via ${\epsilon=\kappa\epsilon_{0}}$, the Coulomb coupling strength in Eqs. (10) and (B) reads

[TABLE]

for graphene and

[TABLE]

for bilayer graphene. The static limit Geick et al. (1966) for h-BN amounts to setting ${\kappa=5.09}$.

We restrict the numerical evaluations to square-shaped sheets of 2D materials with edge length of ${L=200\,\mathrm{nm}=406.5\,l}$, unless stated otherwise. A graphene sheet of $(200\,\mathrm{nm})^{2}$ consists of approximately 1.5 million carbon atoms, while the number of charge carriers, derived from the $n_{3}$-type realizations of $n_{\pm}(\mathbold{r})$, typically ranges between 100 and 10000 for layer L1 — depending on the employed disorder strengths. We confirmed the convergence of the correlation measure $\xi$ by increasing the resolution — here a modest grid size of $105^{2}$ grid points proved enough for converging $\xi$ sufficiently.

The main contribution to the quantum-corrected density $n_{3}(\mathbold{r})$ stems from the effective potential $V$ in the vicinity of $\mathbold{r}$, see Table 1, but the spatial integral for $n_{3}$ is a priori indefinite. We restrict this spatial integral to the square-shaped sheet of area $L^{2}$. Hence, the bulk of the sheet is modeled adequately, but the region close to the edges requires further consideration. Two scenarios suggest themselves. First, the Coulomb tails of the charged impurities could be taken into account beyond the here employed sheet area $L^{2}$ (which would increase the computational cost substantially, with little effect on $\xi$). For this case of a truly finite sheet also the linear (quadratic) dispersion for graphene (bilayer graphene) would require a modification due to the non-periodic boundary conditions. Second, we could view the sheet as a representative part of a larger sample whose disorder potential is unknown beyond the employed sheet area888An extension of the position space support by employing appropriate reflections of the disordered sheet at its edges, would counteract this inaccuracy, but likely introduce spurious effects due to the artificial periodicity, along with a substantial increase in computational cost.. That is, both scenarios render the $n_{3}$ densities less precise near the edges (and particularly near the corners) of the sheet when omitting the spatial regions beyond the employed area $L^{2}$ (which amounts to setting ${V(\mathbold{z})=\mu}$ for ${\mathbold{z}\notin L^{2}}$). We deem this inaccuracy acceptable since the majority of the sheet is adequately taken into account, and the effect on $\xi$ is likely irrelevant given the observed variation due to different disorder realizations.

Appendix B Exchange

Exchange potentials for graphene are available, for example in Ref. [Rossi and Das Sarma, 2008],

[TABLE]

in the units used here. Figure 9 illustrates the exchange effects on the converged effective potential $V$, and Table 2 demonstrates that exchange is negligible as far as inter-layer correlations are concerned.

An approximate exchange potential for other materials, viz.  other dispersion relations, can for instance be obtained by calculating the according Dirac exchange energy for 2D systems and adding correlation in the spirit of Ref. [Vosko et al., 1980]. It would also be interesting to see if exchange-correlation functionals developed for the 2D electron gas like in Ref. [Pittalis and Räsänen, 2010] can be used for systems with quadratic dispersion differing in the particle mass, e.g. for effective bilayer graphene. This topic is, however, beyond the scope of this work.

Appendix C Disorder potentials

We separate the sheet (with optional strain) of the first layer L1 from the plane that holds the charged impurities by ${d=1}\,$nm, cf. Fig. 3. The impurity-induced Coulomb potential in a sheet at distance $d$ reads

[TABLE]

where $C$ takes on values ${c_{i}\in\{-1,0,1\}}$ randomly on the grid points ${i\in\{1,\dots,\Lambda\}}$. We employ an average density of charged impurities $n_{\mathrm{imp}}$ (e.g., $10^{12}\,\mathrm{cm}^{-2}$, corresponding to 400 impurities on $L^{2}$), with half of the impurities positively charged, such that the net charge in the impurity plane is zero. We model the strain as

[TABLE]

with $R$ controlling the relative impact of strain and charged impurities, while $\tilde{C}$ and $C$ represent different disorder realizations of the same type. The sum of Eqs. (15) and (16) is the external potential $V_{\mathrm{ext}}$, cf. Eq. (8), for a monolayer setup as well as for the first sheet of the double-layer structure. The second sheet of the double-layer structure is not affected by $V_{\mathrm{S}}$, but only by the charged impurities at a distance of ${\tilde{d}=6}\,$nm, cf. Fig. 3, and by the electrostatic potential from the converged carrier distribution $n_{1}$ on the first sheet:

[TABLE]

Appendix D Correlation measures

In Section V we use a normalized version of the interlayer correlation measure Rodriguez-Vega et al. (2014)

[TABLE]

to quantify the degree of (anti-)correlation between the density distributions $n_{1}$ and $n_{2}$. Here, $\left<\right.{n_{j}}\left.\right>=\frac{1}{\Lambda}\sum_{i=1}^{\Lambda}n_{j}(i)$ is the mean of the discretized density of layer $j$ on the grid points ${i\in\{1,\dots,\Lambda\}}$, and the root mean square of the density fluctuations in layer $j$ is ${n_{j,\mathrm{rms}}=\sqrt{\left<\right.{(n_{j}-\left<\right.{n_{j}}\left.\right>)^{2}}\left.\right>}}$.

To check for spurious artefacts of the correlation measure $\xi$ in Eq. (18) and provide an error estimate for the (anti-)correlations, we also consider the alternative measure

[TABLE]

which yields stronger contrast between situations of correlation and anticorrelation — at the expense of larger variance; cf. Fig. 10. Here we define the average

[TABLE]

of a quantity $A$ and the corresponding root mean square $n_{j,\mathrm{rms},c}=\sqrt{\left<\right.{(n_{j}-\left<\right.{n_{j}}\left.\right>_{c})^{2}}\left.\right>_{c}}$. The characteristic function $\chi_{c}(i)$ at grid point $i$ equals one if ${n_{j}(i)\notin I_{j}}$ (for both ${j=1,2}$) and zero otherwise999All grid points are taken into account, i.e., ${\chi_{c}(i)=1}$ for all $i$, if ${c=1/2}$. Then, Eqs. (18) and (19) coincide. For ${c=0}$ all data is dismissed. — with the cutoff interval $I_{j}=[\mathrm{min}_{i}\{n_{j}(i)\}+c\,\Delta_{j},\mathrm{max}_{i}\{n_{j}(i)\}-c\,\Delta_{j}]$ and the density spread ${\Delta_{j}=\mathrm{max}_{i}\{n_{j}(i)\}-\mathrm{min}_{i}\{n_{j}(i)\}}$. We choose ${c=0.4}$ to dismiss the fluctuations of correlation at small local carrier density. Owing to its increased variance, we refrain from depicting $\xi_{c}$ in Sec. V.

Appendix E Self-consistent ground-state variables

In this appendix we outline the numerical procedures for obtaining the ground-state variables of DPFT, viz. stationary solutions, $n$, $V$, and $\mu$ of Eqs. (3)–(5), adapted to the quasi-particle system described in Sec. III. Figure 11 illustrates the main steps for implementing the self-consistent loop, with a detailed description below.

Initialization ${(\mbox{iteration }i=0)}$:

• (I)

Calculate ${n_{\pm}^{(0)}=n_{\pm}^{(i=0)}=n_{\nu}[V_{\pm}^{(0)}-\mu^{(0)}]}$ using ${V_{-}^{(0)}=V_{\mathrm{ext}}}$ from Eq. (8) and $V_{+}^{(0)}$ from Eq. (9). The initial choice of the chemical potential is arbitrary, e.g., ${\mu^{(0)}=0}$. We omit $\mathbold{r}$-dependences in this appendix for the sake of notational simplicity. ${n_{\nu}[V-\mu](\mathbold{r})}$ are functionals of ${V-\mu}$ with parametric dependence on $\mathbold{r}$, for example, the density expressions given in Table 1.

• (II)

Adapt $\mu^{(0)}$ until ${N=\int(\mathrm{d}\mathbold{r})\,n^{(0)}=\int(\mathrm{d}\mathbold{r})\,(n_{-}^{(0)}-n_{+}^{(0)})}$. For instance, we demand ${N=0}$ to be reached within an absolute accuracy of $10^{-4}$.

Self-consistent loop:

• (III)

Update effective potential

[TABLE]

and determine effective potential ${V_{+}^{(i+1)}}$ via Eq. (9).

• (IV)

Update densities ${\tilde{n}_{\pm}^{(i+1)}:=n_{\nu}[V_{\pm}^{(i+1)}-\mu^{(i+1)}]}$, starting with ${\mu^{(i+1)}=\mu^{(i)}}$, and adjusting $\mu^{(i+1)}$ until ${N=\int(\mathrm{d}\mathbold{r})\,\tilde{n}^{(i+1)}=\int(\mathrm{d}\mathbold{r})\,(\tilde{n}_{-}^{(i+1)}-\tilde{n}_{+}^{(i+1)})}$.

• (V)

Mix old and new densities with mixing parameter ${\theta\in(0,1)}$: ${n^{(i+1)}=(1-\theta)\,n^{(i)}+\theta\,\tilde{n}^{(i+1)}}$. More sophisticated density mixing like Pulay or Broyden mixing can be expected to improve the convergence behaviour.

• (VI)

Check convergence with an appropriate norm criterion ${||n^{(i+1)}-n^{(i)}||<\epsilon}$. We use ${\epsilon=10^{-6}}$ and

[TABLE]

with grid positions $\mathbold{r}_{j}$,

[TABLE]

If densities have converged to within $\epsilon$, the stationary solutions $n$, $V$, and $\mu$ of Eqs. (3)–(5) are found; otherwise return to step (III) for the next iteration.

The number of self-consistent iterations $i$ required for convergence is largely determined by the density-mixing parameter $\theta$, which is highly system-dependent. Here, we can choose ${\theta=0.5}$ for weak disorder, but need to approach the equilibrium very gently in case of strong disorder ($\theta\lesssim 0.01$). In our implementation $i$ usually reaches values on the order of $\sim 10/\theta$.

Appendix F The quantum-corrected density $n_{3}$

for linear dispersion

In this appendix, we outline the derivation of the density expressions for linear dispersion given in Table 1. Equations (3) and (6) yield Eq. (7) with the aid of the Fourier transform of the step function $\eta(\cdot)$:

[TABLE]

where the degeneracy factor $g$ is part of the trace in Eq. (6). Generally, a Suzuki-Trotter factorization $U_{\nu}$ of the time-evolution operator with

[TABLE]

and real-valued functions $c_{i}$ and $A$ fixed by the here inserted completeness relations in position and momentum space, implies

[TABLE]

For example, the case $U_{3}$ in Sec. II is formally obtained by setting ${\alpha_{1}=\alpha_{2}=-\frac{\mathrm{i}t}{2\hbar}}$, ${\beta_{1}=-\frac{\mathrm{i}t}{\hbar}}$, and ${\beta_{2}=0}$ in $U_{4}$. Note that $U_{3}$ is denoted as $U_{3^{\prime}}$ in Ref. [Chau et al., 2018].

We obtain the TF density by choosing ${\nu=2}$ and ${\alpha_{1}=\beta_{1}=-\frac{\mathrm{i}t}{\hbar}}$. Specifying linear dispersion, ${T(\mathbold{p})=v_{\mathrm{F}}|\mathbold{p}|}$, we find

[TABLE]

whose analytic result is easily obtained and given in Table 1. The computation of $n_{3}$ requires the insertion of three completeness relations, and the so introduced two-dimensional momentum integrations are reduced to one-dimensional integrals via appropriate coordinate transformations. The subsequent analytic evaluations of the remaining one-dimensional integrals lead to the expression for $n_{3}$ in Table 1, leaving one spatial integral to be evaluated numerically.