Internal screening and dielectric engineering in magic-angle twisted bilayer graphene

TL;DR

This paper investigates how internal screening and dielectric environment influence electronic interactions in magic-angle twisted bilayer graphene, proposing dielectric engineering as a method to explore its quantum many-body states.

Contribution

It reveals the role of interlayer coupling in enhancing internal screening and shows how dielectric environment effects depend on the electronic state, suggesting experimental dielectric tailoring.

Findings

Interlayer coupling enhances internal screening.

Dielectric environment effects depend on electronic state.

Proposes dielectric environment engineering to study quantum states.

Abstract

Magic-angle twisted bilayer graphene (MA-tBLG) has appeared as a tunable testing ground to investigate the conspiracy of electronic interactions, band structure, and lattice degrees of freedom to yield exotic quantum many-body ground states in a two-dimensional Dirac material framework. While the impact of external parameters such as doping or magnetic field can be conveniently modified and analyzed, the all-surface nature of the quasi-2D electron gas combined with its intricate internal properties pose a challenging task to characterize the quintessential nature of the different insulating and superconducting states found in experiments. We analyze the interplay of internal screening and dielectric environment on the intrinsic electronic interaction profile of MA-tBLG. We find that interlayer coupling generically enhances the internal screening. The influence of the dielectric environment on the effective interaction strength depends decisively on the electronic state of MA-tBLG. Thus, we propose the experimental tailoring of the dielectric environment, e.g. by varying the capping layer composition and thickness, as a promising pursuit to provide further evidence for resolving the hidden nature of the quantum many-body states in MA-tBLG.

Figures

Equations

$${\Pi }_0^{\mathbb{G},{\mathbb{G}}^{\prime }}(q)=\frac{g_s{g}_v}{S_{M}N}\sum _{\begin{array}{c}\mathbb{k}\\ \alpha ,\beta \\ {\mathbb{G}}_2,{\mathbb{G}}_2^{\prime }\end{array}}{\mathcal{M}}_{{\mathbb{G}}_2,{\mathbb{G}}_2^{\prime },\mathbb{G},{\mathbb{G}}^{\prime }}^{\alpha \beta }\frac{f_{\mathbb{k}}^{\alpha }-f_{\mathbb{k}+q}^{\beta }}{i\eta +E_{\mathbb{k}}^{\alpha }-E_{\mathbb{k}+q}^{\beta }}.$$

$${\mathcal{M}}_{{\mathbb{G}}_2,{\mathbb{G}}_2^{\prime },\mathbb{G},{\mathbb{G}}^{\prime }}^{\alpha \beta }=(c_{\mathbb{k}}^{\alpha ,{\mathbb{G}}_2}{)}^{†}(c_{\mathbb{k}+q}^{\beta ,{\mathbb{G}}_2^{\prime }-{\mathbb{G}}^{\prime }}{)}^{†}{c}_{\mathbb{k}+q}^{\beta ,{\mathbb{G}}_2^{\prime }}{c}_{\mathbb{k}}^{\alpha ,\mathbb{G}+{\mathbb{G}}_2}.$$

$${\Pi }_0^{\mathrm{Dirac}}(q)=g_s{g}_v{g}_l\frac{q}{16\hslash v_F},$$

$${\Pi }_0^{\mathrm{D}+\mathrm{H}}(q)={\Pi }_0^{\mathrm{Dirac}}(q)+a\frac{q^2}{{\left(1+{\left(\frac{q}{b}\right)}^2\right)}^5}.$$

$$W(q)=\frac{V(q)}{\epsilon (q)},$$

$${\epsilon }_{\text{env}}(q)={\epsilon }_1\frac{1-{\tilde{\epsilon }}_2^2{e}^{-4qd}}{{\left(1+{\tilde{\epsilon }}_2{e}^{-2qd}\right)}^2},$$

$$U=\int \frac{d^{2}q}{(2\pi {)}^2}W(q)={\int }_0^{q_c}\frac{dq}{2\pi }qW(q).$$

$$H\approx \left(\begin{array}{cc}{H}^{+}& 0\\ 0& H^{-}\end{array}\right),$$

$$H^{\xi }=\left(\begin{array}{cc}{h}_1& h_{\theta }^{†}\\ h_{\theta }& h_2\end{array}\right),$$

$$h_l=-\hslash v_F\left[\mathcal{R}(\pm \theta /2)(\mathbb{k}-{\mathbb{K}}_l^{\xi })\right]\cdot (\xi {\sigma }_x,{\sigma }_y),$$

$$h_{\theta }=$$

$$+\left(\begin{array}{cc}u& u^{\prime }{\omega }^{-\xi }\\ u^{\prime }{\omega }^{\xi }& u\end{array}\right)e^{i\xi {\mathbb{G}}_1^M\cdot r}$$

$$+\left(\begin{array}{cc}u& u^{\prime }{\omega }^{\xi }\\ u^{\prime }{\omega }^{-\xi }& u\end{array}\right)e^{i\xi ({\mathbb{G}}_1^M+{\mathbb{G}}_2^M)\cdot r},$$

$${\varphi }_{\mathbb{k}}^{\alpha }(r)=\sum _{\mathbb{G}}{c}_{\mathbb{k}}^{\alpha }{e}^{i(\mathbb{k}+\mathbb{G})\cdot r},$$

$${\Pi }_0^H(q)\sim \frac{q^2}{{\left(1+{\left(\frac{q}{b}\right)}^2\right)}^5},$$

$${\Pi }_0^{\mathbb{G},{\mathbb{G}}^{\prime }}(q)=\frac{g_s{g}_v}{S_{M}N}\sum _{\begin{array}{c}\mathbb{k}\\ \alpha ,\beta \\ {\mathbb{G}}_2,{\mathbb{G}}_2^{\prime }\end{array}}{\mathcal{M}}_{{\mathbb{G}}_2,{\mathbb{G}}_2^{\prime },\mathbb{G},{\mathbb{G}}^{\prime }}^{\alpha \beta }\frac{f_{\mathbb{k}}^{\alpha }-f_{\mathbb{k}+q}^{\beta }}{i\eta +E_{\mathbb{k}}^{\alpha }-E_{\mathbb{k}+q}^{\beta }}[1-\Theta (\Delta -|E_{\mathbb{k}+q}^{\beta }-\mu |)\Theta (\Delta -|E_{\mathbb{k}}^{\alpha }-\mu |)].$$

Full text

Introduction. Stacking two graphene layers at a twist angle $\theta$ on top of each other leads to twisted bilayer graphene (tBLG) featuring a moiré pattern with an intricate emergent low-energy electronic structure. For small twist angles $\theta<2\,^{\circ}$, the resulting superlattices host several thousand atoms per unit cell. In this situation, the electronic bands around the charge neutrality point (CNP) become very flat Suárez Morell et al. (2010); Bistritzer and MacDonald (2011), which facilitates strong correlation effects. Recent experiments Cao et al. (2018a, b); Yankowitz et al. (2019); Choi et al. (2019) reported the emergence of possibly unconventional superconducting and insulating states in magic-angle tBLG (MA-tBLG) at different levels of doping. The insulating states occur for commensurate fillings at both electron and hole dopings Cao et al. (2018a); Yankowitz et al. (2019); Lu et al. (2019), signaling a possible Mott-Hubbard origin Koshino et al. (2018); Rademaker and Mellado (2018); Pizarro et al. (2019). Around these insulating states, superconductivity emerges Cao et al. (2018b); Yankowitz et al. (2019); Lu et al. (2019), resembling the phase diagram of high-$T_{c}$ cuprate superconductors Keimer et al. (2015); Kordyuk (2015) and other unconventional superconductors Kurosaki et al. (2005); Sipos et al. (2008). Different models Po et al. (2018); Roy and Juričić (2019); Dodaro et al. (2018); Liu et al. (2018); Rademaker and Mellado (2018); Kennes et al. (2018); Isobe et al. (2018); Koshino et al. (2018); Pizarro et al. (2019); Wu et al. (2018); Thomson et al. (2018); Ochi et al. (2018); Guinea and Walet (2018); Po et al. (2018); Laksono et al. (2018); Venderbos and Fernandes (2018); Choi and Choi (2018); Kozii et al. (2018); Kang and Vafek (2018); Seo et al. (2018); Xie and MacDonald (2018) have been proposed to understand the physics behind these insulating and superconducting states.

Electronic correlations in ultrathin systems such as MA-tBLG depend decisively on the effective electron-electron interaction profile, which is determined by a delicate interplay of screening processes taking place in the dielectric environment and the material itself Wehling et al. (2011); Rösner et al. (2015); Schönhoff et al. (2016). The profile of the effective electron-electron interactions presents a central uncertainty in the current understanding of MA-tBLG.

In this Letter, we provide a quantitative study of the internal polarizability and effects of the dielectric environment on the effective electron-electron interaction in MA-tBLG. We consider MA-tBLG in dielectric surroundings, see Fig. 1(a-b), directly resembling different experimental setups: (a) MA-tBLG is separated by possible hexagonal boron nitride (hBN) encapsulation from a dielectric environment ($\varepsilon_{2}$) at distance $d$ below and above. The case of a metallic gate at a certain distance $d$ from MA-tBLG is included in this analysis for $\varepsilon_{2}\rightarrow\infty$. (b) The special case of the dielectric $\varepsilon_{2}$ in direct contact with MA-tBLG as realized for $d=3.35\>\text{\AA}$. Our main findings are: $(1)$ intrinsic screening in MA-tBLG is enhanced by interlayer coupling and is larger than previously assumed in the model of uncoupled graphene double layers Stauber and Kohler (2016); Cao et al. (2018a); Koshino et al. (2018); Zhang (2018); Liu et al. (2018); Guinea and Walet (2018); González and Stauber (2019); Sboychakov et al. (2018); Fu et al. (2018); González et al. (1999), $(2)$ in a possible Mott insulating phase, the effective interaction can be modulated on the order of $50\%$ by changing the dielectric environment, and $(3)$ metallic states of MA-tBLG are mostly insensitive to the dielectric environment. Finally, a generically applicable, realistic, and yet simple model of the effective electron-electron interaction in MA-tBLG setups as depicted in Fig. 1 is provided.

Internal screening. As a first step, we calculate the intrinsic polarization function of MA-tBLG in the random phase approximation (RPA) using the low-energy continuum model for the electronic bands of tBLG from Refs. Bistritzer and MacDonald (2011); Koshino et al. (2018). In this model SM , two sets of Dirac electron bands originating from the lower and the upper graphene layers hybridize due to the interlayer coupling, which is modulated with the periodicity of the moiré superlattice (c.f. Fig. 1(c)). In reciprocal space, there is correspondingly a coupling between Dirac electron states at each $\mathbb{k}$ vector in one layer with states at $\mathbb{k}+\mathbb{G}$ in the other layer, where $\mathbb{G}$ is a reciprocal lattice vector associated with the moiré superlattice: $\mathbb{G}=m\mathbb{G}_{1}^{M}+n\mathbb{G}_{2}^{M}$ with $m$, $n$ integers and $\mathbb{G}_{1}^{M}$ and $\mathbb{G}_{2}^{M}$ spanning the reciprocal lattice of tBLG, see Fig. 1(d). Then, the electronic states of tBLG are expanded in terms of coupled two-dimensional Dirac spinors up to a certain plane wave cutoff $G_{c}=8G_{M}$, where $G_{M}=|\mathbb{G}_{1}^{M}|=|\mathbb{G}_{2}^{M}|$ SM . MA-tBLG is then realized when the twisting angle is set to $1.05\>^{\circ}$.

In reciprocal space, the polarization operator $\Pi_{0}$ is a function of the scattering vector $\mathbb{q}$ and a matrix indexed by reciprocal lattice vectors $\mathbb{G}$ and $\mathbb{G}^{\prime}$. At zero transferred frequency, $\Pi_{0}$ takes the form Vonsovsky and Katsnelson (1989)

[TABLE]

Here, $g_{s}=g_{v}=2$ are the spin and valley degeneracy factors, $S_{M}=\sqrt{3}\lambda^{2}/2$ is the moiré unit cell area, where $\lambda\approx 134\>\text{\AA}$ is the moiré lattice constant, $N=100$ is the number of $\mathbb{k}$ points, $\alpha$ and $\beta$ encode the sublattice ($A_{1}$, $B_{1}$, $A_{2}$ and $B_{2}$) and band indices SM , $\mathbb{G}_{2}$, $\mathbb{G}^{\prime}_{2}$, $\mathbb{G}$ and $\mathbb{G}^{\prime}$ are reciprocal lattice vectors, $f_{\mathbb{k}}^{\alpha}$ and $E_{\mathbb{k}}^{\alpha}$ are, respectively, the Fermi function and the band energy of $\alpha$ at $\mathbb{k}$, $\eta=0.5\cdot 10^{-6}-10^{-8}\>\text{eV}$ is the broadening parameter, and we consider the temperature $T\approx 50$ K (inverse temperature $\beta=200\>\text{eV}^{-1}$) to calculate the Fermi functions $f_{k}^{\alpha}$ in Eq. (1) SM . $\Pi_{0}$ is thus given in units of $\text{eV}^{-1}\,\text{\AA}^{-2}$. The overlap matrix $\mathcal{M}$ results from the Dirac spinor plane wave expansion coefficients $c_{\mathbb{k}}^{\alpha,\mathbb{G}}$ via

[TABLE]

We then study the internal screening in MA-tBLG for two possible scenarios: first, screening according to the RPA, which assumes MA-tBLG to be in a conventional metallic state; second, a scenario resembling insulating states which are modeled within the so-called constrained RPA (cRPA) Aryasetiawan et al. (2004). In the cRPA, polarization processes taking place inside the low-energy flat bands are excluded from Eq. (1) by setting $f_{\mathbb{k}}^{\alpha}=\frac{1}{2}$ for all states $\alpha$ from this low-energy sector. The resulting cRPA partially screened interactions are those which should be considered in effective Hamiltonians solely dealing with the low-energy flat bands, and correspond also to the screened interactions that would be expected if electronic correlations suppress low-energy polarization processes like in a Mott insulator van Loon et al. (2014); Mot .

The results presented in the following are obtained at charge neutrality with the continuum model from Ref. Koshino et al. (2018) at a twist angle of $\theta=1.05\>^{\circ}$. We show in the supplemental material SM that our results are robust against changes in twist angle, doping, and interlayer coupling parameters.

Fig. 2 compares the RPA and cRPA polarization functions SM ; GGp for MA-tBLG to the case of uncoupled tBLG, where the interlayer coupling is neglected. Uncoupled tBLG hosts 8 flavors of ideal Dirac fermions due to spin $g_{s}=2$, valley $g_{v}=2$, and layer degeneracy $g_{l}=2$ resulting in a linearly-$q$ dependent polarization function Katsnelson (2012):

[TABLE]

where $v_{F}$ is the Fermi velocity. For $q<G_{M}$, the intrinsic polarization functions of MA-tBLG, both in RPA and cRPA, systematically exceed the uncoupled tBLG model. Intrinsic screening is, thus, significantly larger in MA-tBLG than previously assumed Stauber and Kohler (2016); Cao et al. (2018a); Koshino et al. (2018); Zhang (2018); Liu et al. (2018); Guinea and Walet (2018); González and Stauber (2019); Sboychakov et al. (2018); Fu et al. (2018); González et al. (1999). Only at larger momentum transfer $q\gtrsim 1.5G_{M}\approx 0.08\>\text{\AA}^{-1}$, the RPA and cRPA polarization functions of MA-tBLG approach values corresponding to the uncoupled model. This behavior can be understood from the interplay of intra- and interlayer coupling: at sufficiently large momentum transfer $q>G_{M}$, the intralayer coupling $\hbar v_{F}|q|\gtrsim 0.3\>\text{eV}$ dominates over the interlayer one. At $q<G_{M}$, however, the interlayer coupling cannot be neglected.

In the RPA model, i.e. in the metallic case, the enhancement of screening for $q<G_{M}$ is due to the high density of states (DOS) around the Fermi level originating from the flat bands. In the cRPA screening processes inside flat bands are excluded, which explains why the polarization function is smaller than in the RPA case. Still, the cRPA polarizability is clearly enhanced as compared to the uncoupled case for $q<G_{M}$. This can be understood as follows: While polarization processes taking place entirely inside the flat band manifold are excluded in cRPA, gapped transitions between states associated with peaks in the DOS, e.g. between the flat bands and higher energy states with corresponding peaks in the MA-tBLG DOS around $E\approx\pm 20$ meV and $\pm 60\,$meV SM are possible. These gapped transitions are reminiscent of atomic systems. The hydrogen atom, as simplest example, would yield a q-dependent polarizability of the form $\Pi_{0}^{\rm{H}}(q)\sim\frac{q^{2}}{\left(1+\left(\frac{q}{b}\right)^{2}\right)^{5}}$ Yang et al. (1991), where the parameter $b$ is an effective inverse orbital radius. Superimposing this quasi-atomic model $\Pi_{0}^{\rm{H}}$ with a Dirac electron background capturing all higher energy processes leads to the ansatz:

[TABLE]

Fig. 2 shows that Eq. 4 fits the cRPA numerics very well, upon choosing $a\approx 4.1\>\text{eV}^{-1}$ and $b\approx 0.088\>\text{\AA}^{-1}$. Thus, we suggest to use Eq. (4) for calculations of cRPA screened Coulomb interaction matrix elements in MA-tBLG.

While screening in MA-tBLG has been often approximated in terms of uncoupled tBLG Cao et al. (2018a); Koshino et al. (2018); Zhang (2018); Liu et al. (2018); Guinea and Walet (2018); González and Stauber (2019); Sboychakov et al. (2018); Fu et al. (2018); González et al. (1999), an alternative point of view is that MA-tBLG is a patchwork of AA and AB stacked bilayer graphene (BLG) regions. For $q\gtrsim G_{M}/2$ the MA-tBLG polarizability fits, indeed, into the range marked by the low-energy density of states (DOS) of AA BLG, $N^{AA}(E_{F})=0.0158\>\text{eV}^{-1}\,\text{\AA}^{-2}$, and AB BLG, $N^{AB}(E_{F})=0.0042\>\text{eV}^{-1}\,\text{\AA}^{-2}$ Hwang and Das Sarma (2008). Further below we show that the AB BLG model gives good estimates for the effective local interactions when compared with the cRPA. On physical grounds, our results for the polarization function show that MA-tBLG can be seen as a flat band system embedded in an almost metallic background generated by the higher energy bands, which are separated by approximately $15\>\text{meV}$ from the low-energy flat bands SM .

Dielectric engineering. We now assess the possibilities for dielectric engineering of MA-tBLG in experimental setups as shown in Fig. 1(a-b). The screened interaction is

[TABLE]

where $V(q)=2\pi e^{2}/q$ is the bare interaction and $\varepsilon(q)=\varepsilon_{\text{env}}(q)+V(q)\Pi_{0}(q)$ is the dielectric function Katsnelson (2012). The latter accounts for the screening resulting from the electrons in the bands of MA-tBLG via $\Pi_{0}(q)$ plus the screening produced by the dielectric surrounding and higher energy bands of MA-tBLG summarized in $\varepsilon_{\text{env}}(q)$. In the setups from Fig. 1, the background dielectric function reads Rösner et al. (2015)

[TABLE]

where we approximated $\varepsilon_{1}\approx 4$ Katsnelson (2012) and $\tilde{\varepsilon}_{2}=(\varepsilon_{1}-\varepsilon_{2})/(\varepsilon_{1}+\varepsilon_{2})$. Together with the polarizability $\Pi_{0}(q)$ obtained either numerically in cRPA from Eq. (1) or conveniently from the fit of Eq. (4), Eqs. (5) and (6) specify the appropriately screened interaction of MA-tBLG in different dielectric surroundings, which should be incorporated in interacting electron models of the low-energy flat bands.

To illustrate the influence of the dielectric surrounding on the effective interaction in MA-tBLG in terms of a single descriptor, we consider an effective local interaction $U$ obtained by Fourier transformation of $W(q)$:

[TABLE]

Here, the cutoff $q_{c}$ relates to the 2D radial spread $r_{WF}$ of the Wannier functions (WF) constructed from the MA-tBLG low-energy bands: $q_{c}\approx\pi/r_{WF}$. We choose $r_{WF}$ to be half of the moiré lattice constant $\lambda$, hence $q_{c}\approx 0.047\>\text{\AA}^{-1}$.

Fig. 3(a) shows the influence of the dielectric environment $\varepsilon_{2}$ on $U$, where we assume that the surrounding dielectric is in direct contact with MA-tBLG (c.f. Fig. 1(b)). Assuming a metallic state of MA-tBLG as in the RPA model, the dielectric environment has barely any effects, and $U$ remains almost constant around $10\>\text{meV}$. This insensitivity of metallic states in MA-tBLG to the dielectric surrounding is due to a high intrinsic polarizability which masks any environmental polarizabilities (c.f. RPA curves in Fig. 2 and Fig. S3 of the supplemental material SM ).

On the other hand, the cRPA estimated $U$ varies from $40\>\text{meV}$ to $\approx 20\>\text{meV}$ when $\varepsilon_{2}$ increases. Interestingly, $40\>\text{meV}$ is approximately the splitting of the upper and lower Hubbard bands when a gap is opened at the CNP, as recently measured in scanning tunneling spectroscopy (STS) Choi et al. (2019). Besides its relevance for low-energy effective Hamiltonians, the cRPA solution resembles the fully screened interaction in certain insulating states, where screening from the low-energy states is suppressed. In this situation the dielectric environment can modify the effective interaction $U$ up to $\approx 40-50\,\%$.

The influence of the metallic gate ($\varepsilon_{2}\rightarrow\infty$) at a distance $d$ from MA-tBLG (c.f. Fig. 1(b)) on the effective local interaction $U$ is shown in Fig. 3(b). In an insulating state (cRPA case), the influence of the metallic gate is strong if the distance from MA-tBLG is on the order of $d\lesssim 20\>\text{\AA}$. In this case, $U$ is reduced by a factor of two when approaching the gate from $d\sim 30\>\text{\AA}$ to $d\sim 10\>\text{\AA}$. The fact that the metallic gate has to be as close as $d\lesssim 20\>\text{\AA}$, which is much smaller than the moiré lattice constant ($134\>\text{\AA}$), in order to see any substancial effect on $U$ even in the cRPA / insulating case is due to MA-tBLG being a flat-band system embedded in an already almost metallic background. In the RPA case, the metallic gate has again a much smaller effect on $U$. Only for metallic gates very close to MA-tBLG ($d\approx 3.35\>\text{\AA}$), $U$ is reduced from $10\>\text{meV}$ to $5\>\text{meV}$. Thus, in metallic states of MA-tBLG (RPA case), metallic gates can only influence the electron-electron interactions due to external screening if they are very close to the MA-tBLG, and a possible hBN spacer separating the gate from the MA-tBLG should be as thin as a monolayer. The dependence of the effective local cRPA screened interaction $U$ on the dielectric surrounding $\varepsilon_{2}$ and its distance $d$ from MA-tBLG is summarized in Fig. 3(c). $U$ can be reduced from $40\>\text{meV}$ (free standing case) to $10\>\text{meV}$ (strong screening environment) if the surrounding dielectric gets sufficiently close $d\lesssim 20\>\text{\AA}$.

Additionally to the RPA and cRPA numerics according to Eq. (1), we also propose simple approximate forms to calculate screened Coulomb interaction matrix elements. First, by choosing the expression of Eq. (4) for the polarizability, we reproduce also the effective cRPA screened local interaction quite accurately. Even more simplistically, Fig. 3 also shows that screening based on the DOS of $AB$ BLG provides us with a good estimate for $U$ in the cRPA case.

Summary. We have shown that the sensitivity of the interaction profile in MA-tBLG to the dielectric environment depends decisively on the quantum state of the electrons in the low-energy flat-band manifold. Currently, the nature of the insulating states at half-filling, of superconductivity near half-filling, and of the gaps measured in STS at charge neutrality are major open problems in the understanding of MA-tBLG. What would an experimentally detected significant influence of the dielectric surrounding on spectral, thermodynamic, or transport properties in any of these phases imply? Based on our results any strong influence of the external dielectric means that charge screening in the low-energy bands is substantially suppressed as compared to a conventional metallic state as assumed in RPA. That is hardly possible for insulating or any kind of ordered states resulting from a weak coupling instability, which affect only a small sector of states in the vicinity of the non-interacting Fermi surface. Thus, significant impact of external dielectrics on correlated insulating states implies that they are at strong coupling in the sense that states from the entire mini-BZ contribute significantly. Similarly, any impact of the dielectric surrounding on superconductivity implies strong coupling physics at work, and possibly hints at an unconventional origin.

Currently, there are arguments in favor and against strong coupling physics governing different parts of the MA-tBLG phase diagram. On the pro side, is the ratio of effective interactions and bandwidth: in case of partial screening as described in cRPA, the effective interaction largely exceeds the bandwidth, which is suggestive of strong coupling. In line with a strong coupling scenario, close to charge neutrality, STS further measured ”gaps” on the scale of several $10\,$meV Choi et al. (2019) and transport experiments reported suppression of conductance around charge neutrality being robust against temperature and magnetic fields Yankowitz et al. (2019); Cao et al. (2018a); Lu et al. (2019). At and around half-filling, the situation is more complex: the generic shape of the phase diagram in this region appears to exhibit similarities to the cuprate high-$T_{c}$ superconductors, which might be seen as an indication pro strong coupling physics. At half filling, however, the insulating states are fragile, i.e. they are easily destroyed by magnetic fields ($B\sim 5\>\text{T}$) and temperature ($T\sim 4\>\text{K}$), in line with gaps that are small ($\Delta\sim 0.3\>\text{meV}$) Yankowitz et al. (2019); Cao et al. (2018a); Lu et al. (2019). All these energies are at least an order of magnitude smaller than the expected width of the flat bands, and would be consistent with a weak coupling scenario. Experimental studies of the response of MA-tBLG to changes in the surrounding dielectric will thus be useful to distinguish between strong- and weak-coupling scenarios, and to pinpoint the nature of the different quantum many-body states in MA-tBLG.

We thank D. A. Abanin, E. Bascones, F. Guinea for useful conversations. JMP and TW acknowledge funding from DFG-RTG 2247 (QM3) and the European Graphene Flagship. RT is supported by DFG-SFB 1170 Tocotronics (project B04), DFG-SPP 1666, and further acknowledges financial support from the DFG through the Würzburg-Dresden Cluster of Excellence on Complexity and Topology in Quantum Matter – ct.qmat (EXC 2147, project-id 39085490). RV acknowledges funding from DFG-TRR 49.

I SUPPLEMENTAL MATERIAL

In section A of this Supplemental Material, we give a short description of the low-energy continuum model developed in Refs. Bistritzer and MacDonald (2011); Koshino et al. (2018) for tBLG and employed here. We show the band structure and the DOS around the charge neutrality point (CNP) for the first magic angle $1.05\>^{\circ}$. In section B, we discuss the matrix structure of the polarization operator given in Eq. (1) of the main text as obtained in RPA and cRPA. Section C analyzes the dielectric function $\varepsilon(q)$ of MA-tBLG in different dielectric environments. Section D discusses the temperature dependence of the polarization functions in RPA and cRPA. In section E, we study the dependencies of internal screening in tBLG on twist angle, doping and variations in vertical coupling parameters used in the continuum model. In section F, we investigate how gaps on the order of $0.1$ meV to $1$ meV opened at the CNP and at half-filling of the flat bands affect the polarization function of MA-tBLG.

I.1 LOW-ENERGY CONTINUUM MODEL FOR tBLG

We outline here the low-energy continuum model which we used for the description of the electronic band structure of MA-tBLG. The model has been introduced in Bistritzer and MacDonald (2011); Koshino et al. (2018). By neglecting the intervalley coupling, the total low-energy Hamiltonian can be written as

[TABLE]

where $H^{\xi}$, $\xi=\pm$, describe the two valleys separately. The intravalley blocks are written in the space of the sublattices $(A,B)$ and layer $(1,2)$ degrees of freedom, $(A_{1},B_{1},A_{2},B_{2})^{\xi}$

[TABLE]

where $h_{l}$ with $l=1,2$ are the intralayer terms, and $h_{\theta}$ the interlayer ones. Then, following derivations in Koshino et al. (2018), the intralayer term can be written as

[TABLE]

where $\pm$ stands for $l=1$ and $2$, respectively, $\hbar v_{F}/a=2.1354\>\text{eV}$, $\mathcal{R}(\theta)$ is the 2D rotation matrix, $\sigma_{x}$ and $\sigma_{y}$ are the 2D Pauli matrices and $\mathbb{K}_{l}^{\xi}$ are the rotated graphene $K$ points of each layer $l$ in each valley $\xi$. For the interlayer term we use

[TABLE]

with $\omega=e^{i2\pi/3}$, the AA interlayer coupling $u=0.0797\>\text{eV}$ and the AB interlayer coupling $u^{\prime}=0.0975\>\text{eV}$. The inequality $u\neq u^{\prime}$ accounts for possible corrugation effects in MA-tBLG Koshino et al. (2018). Uncoupled tBLG calculations are performed by switching off the interlayer terms, i.e. $h_{\theta}\equiv 0$. Due to the valley and spin degeneracy, we simply diagonalize $H^{+}$ to calculate the intrinsic static polarization $\Pi_{0}^{\mathbb{G},\mathbb{G}^{\prime}}$ and including degeneracy factors $g_{v}=2$ and $g_{s}=2$ in Eq. (1).

The diagonalization of $H^{+}$ is performed in the $k$-space. Then, there is a coupling between Dirac electron states at each $\mathbb{k}$ vector in one layer with states at $\mathbb{k}+\mathbb{G}$ in the other layer, where $\mathbb{G}=m\mathbb{G}_{1}^{M}+n\mathbb{G}_{2}^{M}$ is the reciprocal vector associated with the moiré superlattice, $m$ and $n$ are integers, and $\mathbb{G}_{1}^{M}$ and $\mathbb{G}_{2}^{M}$ span the reciprocal lattice of tBLG. The electronic states are

[TABLE]

where $\alpha$ encodes the sublattice / layer $X=A_{1},B_{1},A_{2},B_{2}$ and band $n$ indices Koshino et al. (2018). These eigenstates are expanded up to a certain plane wave cutoff $G_{c}$. In this work, we use $G_{c}=8$ for the calculation of $\Pi_{0}(q)$ in Fig. 2, and $G_{c}=2$ in the calculations of the matrices $\Pi_{0}^{\mathbb{G},\mathbb{G}^{\prime}}(\mathbb{q})$ in Fig. S2.

In Fig. S1, we show the band structure and the DOS of MA-tBLG for the first magic angle $1.05\>^{\circ}$, as well as its comparison with uncoupled tBLG. The flat bands with bandwidth $W\approx 7\>\text{meV}$, as well as the higher energy bands separated by gaps of the order of $10\>\text{meV}$ from the flat bands are well reproduced.

In cRPA, polarization processes taking place inside the flat band manifold are excluded. However, gapped transitions between states associated with peaks in the DOS, e.g. between the flat bands and higher energy DOS peaks around, are still possible. These gapped transitions are reminiscent of the atomic systems such as the hydrogen atom. The q-dependent polarizability from transitions between hydrogen $s$ and $p$-orbitals reads Yang et al. (1991)

[TABLE]

where $b$ is an effective inverse orbital radius. Superimposing this quasi-atomic model with a Dirac electron model capturing all higher energy processes leads to the ansatz of Eq. (4) from the main text. As shown in the main text, this ansatz yields a very good fit of the polarization function in the cRPA case.

I.2 POLARIZATION MATRIX

Fig. S2 shows the polarization matrix $\Pi_{0}^{\mathbb{G},\mathbb{G}^{\prime}}(\mathbb{q})$ for $\mathbb{q}\rightarrow\mathbb{\Gamma}$, $\mathbb{K}$ and $\mathbb{M}$ as obtained in RPA and cRPA using $G_{c}=2$. In the cRPA case (top row of Fig. S2), the diagonal elements are the leading elements for all $\mathbb{q}$. In the RPA case (bottom row of Fig. S2), local field effects, i.e. off-diagonal matrix elements of $\Pi_{0}^{\mathbb{G},\mathbb{G}^{\prime}}(\mathbb{q})$ with $\mathbb{G}\neq\mathbb{G}^{\prime}$, become progresively more important when $\mathbb{q}$ is close to the mini BZ boundaries. However, the diagonal elements are still largest for small $q$, which is the decisive region of $q$-space to determine the impact of substrate screening effects.

I.3 DIELECTRIC FUNCTION FOR tBLG

The screening of Coulomb interactions is encoded in the dielectric function $\varepsilon(q)=\varepsilon_{\text{env}}(q)+V(q)\Pi(q)$, where $\varepsilon_{\text{env}}(q)$ as given in Eq. (6) of the main text accounts for the different dielectric environments (c.f. Fig. 1 of the main text). We neglected local field effects in the study of the dielectric function which is well justified at all $q$ in the cRPA case and at small $q$, which are most relevant here, also in the RPA case (c.f. section 2 of the supplement). Fig. S3 shows the dielectric functions for the cRPA (top row) and RPA (bottom row) for MA-tBLG in different dielectric environments $\varepsilon_{2}$ at different distances $d=3.35$, $10$ and $20\>\text{\AA}$.

In the RPA case, the intrinsic screening is strong. Hence, the dielectric environment has a small effect, as seen in Fig. S3. In the cRPA case, the dielectric environment affects the screening in a considerable $q$-range covering e.g. the whole first mini BZ if the dielectric / metallic gate is as close as $d\lesssim 10\>\text{\AA}$. These behaviors connect with Fig. 3 of the main text, where we have shown that the effective local interaction remains almost unaffected by the dielectric surrounding in the metallic state (RPA case) but can vary up to $40-50\%$ in the Mott insulating state (cRPA case).

I.4 TEMPERATURE DEPENDENCE OF THE POLARIZATION FUNCTION

Fig. S4 shows the temperature dependence of the polarization function $\Pi_{0}(q)$ as calculated in the RPA and cRPA according to Eq. (1) of the main text. We compare $T\approx 50$ K (inverse temperature $\beta=200$ eV*-1*) and $T\approx 10$ K ($\beta=1000$ eV*-1*). In the cRPA, $\Pi_{0}(q)$ is essentially temperature independent in this range due to the fact that cRPA models an effectively gapped system. In contrast, the RPA polarizability is strongly temperature dependent: in RPA, $\Pi_{0}(q)$ increases as the temperature is lowered, which is understandable from the band width of the low energy flat bands being on the order of only a few meV. These results confirm the conclusions drawn in the main text: in cRPA, the system is sensitive to the dielectric environment essentially independently of the temperature. In RPA, on the other hand, the system is not affected by the surrounding dielectrics due to strong internal screening. This statement holds already at $T=50$ K as discussed in the main text and becomes even more pronounced at smaller temperatures as can be seen in Fig. S4(a).

I.5 DEPENDENCE OF INTERNAL SCREENING ON TWIST ANGLES, DOPING, AND INTERLAYER COUPLING PARAMETERS

We assess in the following the robustness of the conclusions drawn in the main text with respect to variations in the experimental conditions (twist angle and doping) and details of the model employed (interlayer coupling parameters). We study the polarization function of tBLG for twist angles away from the first magic angle $1.05\>^{\circ}$, finite doping, and different ratios of AA to AB interlayer coupling $u/u^{\prime}$, which corresponds to different assumptions on vertical relaxations Bistritzer and MacDonald (2011); Koshino et al. (2018); Tarnopolsky et al. (2019).

I.5.1 Twist angle dependence

The twist angle dependence of the internal screening in tBLG can be understood from Fig. S5. The polarization functions obtained in RPA are shown in Fig. S5(a) for twist angles $1.00\>^{\circ}\leq\theta\leq 2.00\>^{\circ}$ and in Fig. S5(b) for $\theta\geq 2.00\>^{\circ}$. Within the range $1.00\>^{\circ}\leq\theta\leq 1.20\>^{\circ}$ the RPA polarization functions behave qualitatively similar to the magic angle case of $\theta=1.05\>^{\circ}$ in the sense that the RPA polarization functions exceed by far the model of uncoupled Dirac fermions. For $\theta>2\>^{\circ}$ the situation changes and the RPA screening recovers the linear-$q$ dependence as expected from Eq. (3), but with a renormalized Fermi velocity $v_{F}^{*}$, where $v_{F}^{*}(\theta=3\>^{\circ})/v_{F}\approx 0.882$, $v_{F}^{*}(\theta=4\>^{\circ})/v_{F}\approx 0.895$, $v_{F}^{*}(\theta=5\>^{\circ})/v_{F}\approx 0.927$ and $v_{F}^{*}(\theta=6\>^{\circ})/v_{F}\approx 0.958$. In the cRPA case (Fig. S5(c)), the low $q$ behavior is similar for all twist angles $\theta<2\>^{\circ}$. Thus, the conclusions on the interplay of external versus internal screening drawn in the main text will be valid not only at the magic angle but in an extented range $1.00\>^{\circ}\leq\theta\leq 1.20\>^{\circ}$ around the magic angle.

In Fig. S5(d), we show impact of dielectric engineering on the effective local interaction $U$ as obtained in RPA and cRPA for twisting angles $1.00\>^{\circ}\leq\theta\leq 1.20\>^{\circ}$. We considered a dielectric $\varepsilon_{2}$ in direct contact with the tBLG corresponding to the setup depicted in Fig. 1(b) of the main text, and calculated $U$ according to Eq. (7) of the main text with $q_{c}$ being rescaled in all cases according to the twist-angle dependent moiré lattice constants. In the entire range of twist angles considered, here, the RPA effective interactions $U$ are insensitive to the dielectric environment, while the effective interactions in cRPA change with $\varepsilon_{2}$ in a similar manner for all twisting angles. Therefore, we can conclude that our main results are robust against variations of the twist angle around the magic angle in the range $1.00\>^{\circ}\leq\theta\leq 1.20\>^{\circ}$. This range of twist angles is far larger than uncertainties of twist angles in experiments of Refs. Cao et al., 2018a, b; Lu et al., 2019 and covers the range of angles within which possible unconventional superconducting and insulating states have been observed therein.

I.5.2 Doping dependence

In Fig. S6, we show RPA/cRPA polarization functions for different levels of doping. The band filling $\nu=n/n_{0}$ is given relative to the CNP $\nu=0$ in multiples of the density $n_{0}$ needed to fill one moiré superlattice band. Due to spin and valley degeneracy, $\nu=\pm 4$ then describes fully empty / filled flat bands and $\nu=\pm 2$ the half-filled flat bands. Around half-filling of the positive and negative flat bands (i.e. $\nu=\pm 2$), we find that RPA screening is of the same order as in the undoped case. Significant changes occur in the RPA case only in the limit of completely empty / filled flat bands ($\nu\approx\pm 4$). The cRPA results are virtually independent of the doping level. Thus, our results hold for the range of dopings at which unexpected superconducting and insulating states have been observed in Refs. Cao et al., 2018a, b; Lu et al., 2019.

I.5.3 Interlayer coupling dependence

The continuum model employed, here, depends sensitively on the ratios of AA to AB interlayer coupling $u/u^{\prime}$, which encodes vertical relaxation effects. Different $u/u^{\prime}$ ratios have been employed in the literature Bistritzer and MacDonald (2011); Koshino et al. (2018); Tarnopolsky et al. (2019). The work of Ref. Bistritzer and MacDonald, 2011 assumed equal AA and AB interlayer coupling $u/u^{\prime}=1$, which neglects the larger interlayer distance in the AA regions as compared to the AB regions. To account for these relaxation effects, $u/u^{\prime}\approx 0.82$ has been considered in Ref. Koshino et al., 2018 and also used for all calculations shown in the main text of this work. In the following, we study the impact of variations in the $u/u^{\prime}$ ratio on our results: Fig S7 shows RPA/cRPA polarization functions for MA-tBLG obtained with $u/u^{\prime}=0.70$, $0.82$ Koshino et al. (2018), $0.90$ and $0.95$ . We obtain almost the same polarization functions as obtained from the model considered in the main text ($u/u^{\prime}=0.82$ Koshino et al. (2018)). Thus, our conclusions are also robust against variations of the $u/u^{\prime}$ ratio in the low-energy continuum model.

Taken together, we have shown that the conclusions drawn on the interplay of internal and external screening in MA-tBLG are robust against variations of system parameters like twist angle, doping, and details of the model parametrization.

I.6 INFLUENCE OF ENERGY GAPS ON INTERNAL SCREENING OF MA-tBLG

In the experiments of Refs. Cao et al. (2018a, b); Yankowitz et al. (2019); Lu et al. (2019) the opening of energy gaps at integer filling factors $\nu$ has been reported. In particular, gaps on the order of $0.1$ meV - $0.3$ meV at half-filling and other integer fillings $\nu\neq 0$ were found experimentally in Refs. Cao et al. (2018a, b); Yankowitz et al. (2019); Lu et al. (2019). Gaps found at the CNP ($\nu=0$) are likely larger Choi et al. (2019); Lu et al. (2019). Currently, the nature of all these gaps is unclear. If (some of) these gaps arise from a weak coupling scenario, it can be expected that screening processes from a low energy window on the order of the gap are suppressed, while all higher energy contributions to the polarization function remain essentially unaffected. To mimic such a situation, we study how exclusion of polarization processes taking place inside an energy window of $\pm\Delta$ around the chemical potential affects the polarization function. To this end, we modify Eq. (1) from the main text to calculate the polarization function in the following manner:

[TABLE]

This formula reduces to Eq. (1) from the main text for $\Delta\to 0$ and corresponds to a flavor of cRPA, where all states from a certain energy window $\pm\Delta$ are assumed to form the correlated subspace. In this flavor of RPA, $f_{\mathbb{k}}^{\alpha}$ denote the non-modified Fermi functions belonging to band $\alpha$ at crystal momentum $\mathbb{k}$.

In Fig. S8, we show the polarization function of MA-tBLG calculated according to Eq. (14) with gaps in the range of $\Delta=0.05\>\text{meV}$ to $1\>\text{meV}$ being included in the simulations. For $\Delta\leq 0.15\>\text{meV}$, the polarization functions are very similar to the ones obtained within full RPA. Thus, weak coupling instabilities opening gaps in the $\lesssim 0.3\>\text{meV}$ range can be expected to barely affect the internal polarization function of tBLG. Correspondingly, we expect that the RPA results shown in the main text are representative for conventional metallic states and gapped states resulting from weak coupling mechanism, as assumed in the main text. For gaps opened by strong coupling mechanisms, the situation can be very different as explained in the main text and illustrated in the $\Delta=1\>\text{meV}$ case in Fig. S8.