Theory of Coulomb drag in spatially inhomogeneous 2D materials

Coulomb drag is a favored experimental probe of Coulomb interactions between layers of 2D materials. In reality, these layers display spatial charge density fluctuations known as puddles due to various imperfections. A theoretical formalism for incorporating density inhomogeneity into calculations has however not been developed, making the understanding of experiments difficult. Here, we remedy this by formulating an effective medium theory of drag that applies in all 2D materials. We show that a number of striking features at zero magnetic field in graphene drag experiment which have not been explained by existing literature emerge naturally within this theory. Applying the theory to a phenomenological model of exciton condensation, we show that the expected divergence in drag resistivity is replaced by a peak that diminishes with increasing puddle strength. Given that puddles are ubiquitous in 2D materials, this work will be useful for a wide range of future studies. An understanding of the interlayer electronic interactions in two dimensional heterostructures is required to advance their potential applications in low-power electronics. The authors develop a transport theory incorporating charge inhomogeneities in order to explain the behavior of Coulomb drag observed experimentally in double layer heterostructures.

I ndirect excitons are quasiparticles that form when an electron from one layer and a hole from another (see Fig. 1) bind together due to Coulomb interaction. Under the right conditions [1][2][3][4] , these may condense into a superfluid phase 5,6 below a critical condensation temperature T c . This phase is believed to hold great promise for low-power electronics 7 and there has been a sustained experimental pursuit of its realization and detection in various two-dimensional material heterostructures [8][9][10][11][12][13][14][15][16] . The tool of choice for detecting a condensate in these experiments is the Coulomb drag measurement 17 which involves driving a current in one (active) electronic layer and measuring the induced potential drop in the other (passive) layer as shown in Fig. 1. The corresponding induced electric field is then divided by the current density in the active layer to yield drag resistivity. Tunneling between the two layers is prevented during this process by a dielectric spacer between them so that the induced electric field is solely due to the interlayer Coulomb interaction. The exciton condensate is then expected to reveal itself as a divergence 18 or enhancement 19 of the drag resistivity.
There exists a large body of theoretical literature on Coulomb drag in 2D systems including a comprehensive review paper 20 . Most of these theoretical works describe drag between two ideal electronic layers each possessing a uniform charge density. In reality, electronic layers are known to possess charge density fluctations known as puddles [21][22][23] as shown in Fig. 1 arising from either charged impurities 24 or local strain in the layer 25 and there is no known method of eliminating them completely. In spite of this, a controlled way of calculating their effects in actual drag measurements has never been developed, creating difficulties in interpreting drag measurements searching for exciton condensates and other many-body phases.
In this work, we formulate a rigorous framework accounting for puddles in Coulomb drag calculations by generalizing Bruggeman's effective medium theory (EMT) approximation [26][27][28][29] for the first time to the case of Coulomb drag between two inhomogeneous sheets (see Fig. 2). The EMT is a machinery that averages over inhomogeneities in a macroscopically inhomogeneous medium in a self-consistent manner in order to obtain its effective behavior. It has been applied in a wide array of contexts, including magnetoresistance in polycrystalline metals, superconductivity in granular materials and hydrodynamics in suspensions to name but a few 29 . This work establishes that this powerful machinery may be brought to bear on the Coulomb drag problem and sheds considerable light on existing experimental data. The resulting drag EMT applies in all 2D materials. We illustrate the utility of this formalism through two examples. First, we consider drag between graphene sheets as studied in the experiment by Gorbachev et al. 8 , focusing only on drag in zero magnetic field. This experiment shows a number of observed features in qualitative contradiction with the predictions of standard theory ignoring puddles, with an unexpected positive peak in drag resistivity at the double neutrality point (DNP) being the one that has received the most attention. The most widely accepted theoretical explanation of the zero-field experimental data invokes a puddle-induced thermoelectric 'energy drag' 30 effect. The work of Song et al. 30 however does not consider the non-zero 'momentum drag' [31][32][33][34][35][36] resistivity that will inevitably be induced at the DNP by these same puddles, leaving the question of what exactly causes the peak still open. A second experimental feature that has not been reproduced from theoretical microscopic calculations is the increase of the finite density drag resistivity peaks as a function of temperature. Microscopic calculations of drag in homogeneous samples in fact yield finite density drag peaks that decrease with increasing temperature (see for instance Fig. 8b of Narozhny et al. 36 ). In fact, no theoretical work to date has successfully reproduced both the above experimental features. We show here that once puddles with negative interlayer correlations are taken into account using drag EMT, momentum drag alone is able to reproduce all these features with nearly quantitative agreement. As a second application, we demonstrate using a phenomenological model that puddles cause the divergence in drag resistivity expected upon exciton condensation to be re-normalized to a peak that goes down as the reciprocal of the puddle strength. These two example applications show that puddles must be taken into account in order to understand existing drag experiments on graphene and suggest that they will also play a role in the behavior of drag resistivity in an exciton condensate. Since height corrugations are ubiquitous in 2D materials 37 and these in turn lead to puddles 25 , we expect that the present formalism will be useful for explaining many future 2D Coulomb drag experiments.

Results
Model. The derivation of drag EMT parallels that of the usual EMT for inhomogeneous single layer transport (see Supplementary Note 1). Consider the Coulomb drag setup of Fig. 1-two parallel 2D electronic layers with charge density fluctuations separated by some finite distance with an insulating spacer material, with a current flowing through the active layer while the passive layer remains open. We model the inhomogeneities by In reality, the twodimensional layers possess density fluctuations (i.e., puddles) as shown but their effects are largely ignored theoretically, leading to contradiction with experiment. We solve this by generalizing the well-established effective medium theory [27][28][29] to the case of Coulomb drag. The scale bar shows local density in units of 10 10 cm −2 Fig. 2 Embedding procedure of effective medium theory. An arbitrary ith pair of inclusions in the inhomogeneous layers are each embedded inside an effective medium. The effective in-plane conductivities of the inhomogeneous samples are given by σ E A and σ E P and the effective drag conductivity by σ E D assuming the active and passive layers are each made up of N inclusions each with its own constant local conductivity, σ A i and σ P i , respectively, where i = 1, …, N. In other words, we model the puddles using a collection of inclusions. These inclusions are assumed to be concentric circles of radius a (i.e., the ith inclusion of the active layer is assumed to lie exactly atop the ith inclusion of the passive). The local drag conductivity σ i D corresponding to each ith pair of inclusions is also assumed to be constant. In order for the assumption of constant local conductivities within the inclusions to be justified, the typical conductivity fluctuations within each inclusion must be small. More precisely, σ r ð Þ À σ i ð Þ =σ i j j ( 1 must be obeyed for all points r within each inclusion, and σ i is the average conductivity of the area within the inclusion. Here σ i refers generically to the intralayer conductivities, as well as the drag conductivity. We do not make any further assumptions about the nature of the inclusions. Choosing an arbitrary ith pair of inclusions, one from each layer, we regard them each as being embedded within an effective medium as shown in Fig. 2. The two effective media have within them effective uniform fieldsẼ A;0 andẼ P;0 representing the total field from all other inclusions, as well as the fieldsẼ A;S andẼ P;S arising from the embedded inclusions. Drag conductivity between the two effective media shall be denoted by σ E D , while the effective media possess effective conductivities σ E A and σ E P , respectively. Within the EMT approximation, σ E A , σ E P and σ E D are the intralayer and drag conductivities that will be obtained in experimental measurements on inhomogeneous substrates. We solve for the electric fields inside the inclusionsẼ i A andẼ i P (see Supplementary Note 2 for details) and substitute them into the EMT selfconsistency equations P i f iẼi A ¼Ẽ A;0 and P i f iẼi P ¼Ẽ P;0 where f i refers to the areal fraction of the ith inclusion relative to the whole layer. These equations enforce the condition within each layer that upon averaging the intra-inclusion field over all inclusions, one should recover the effective uniform fieldẼ 0 . Note that the effective fieldsẼ A;0 andẼ P;0 drop out from the equation once the above substitutions are made as is expected since the above condition makes no demands on the specific values of the effective fields. Modeling the densities of the inclusions with continuous probablity distribution functions, we obtain expressions for the effective drag conductivity σ E D and in-plane conductivities σ E A and σ E P . The effective drag conductivity is then obtained as where n ′ A denotes the local charge density in a small region on the active layer and n ′ P denotes the same in a small region directly beneath it on the passive layer, and σ D n ′ A ; n ′ P À Á denotes the drag conductivity between two sheets of uniform densities n ′ A and n ′ P . P bi n ′ A ; n ′ P À Á is the joint probability distribution of having two such local regions with given charge densities n ′ A ; n ′ P À Á . The effective in-plane conductivities σ E i on the other hand are given by the solutions of where i = A, P henceforth denotes the layer index (instead of inclusion index) and σ i n ′ i À Á is the conductivity of layer i at uniform density n ′ i . P mono n ′ i À Á is the single layer probability density of a point on layer i having charge density n ′ i . This is the well-known effective conductivity expression from single layer EMT 27-29. To calculate drag resistivity in the presence of puddles, we take the clean (i.e., puddle-free) conductivities calculated as a function of n ′ A ; n ′ P À Á and substitute them into the EMT Eqs. (1) and (2). This yields the effective conductivities σ E D , σ E A , σ E P in the presence of puddles. The latter step may be understood as an effective averaging of the homogeneous drag resistivity over densities in a manner prescribed by electrostatics. Finally, the effective drag resistivity ρ E D is obtained by inserting the effective conductivities into The above Eqs. (1)-(3) constitute the generalization of Bruggeman's historic EMT expression 26 to the drag problem. We emphasize that they apply to drag between two layers of arbitrary 2D electronic material. Having completed the derivation, let us make some further caveats. First, so long as the inclusion size a is chosen much larger than the electronic mean free path l in both layers, one may assume the intralayer and drag conductivities σ i n ′ i À Á and σ D n ′ A ; n ′ P À Á are the same as those derived assuming infinite samples (i.e., no finite size effects). Second, we have ignored the drag effect between neighboring inclusions (i.e., the ith inclusion in each layer has a Coulomb drag effect on only the ith inclusion of the other layer and has no effect on other inclusions). This is again justified provided that a ) l in both layers. Last, the above derivation assumes that transport is dominated by scattering within the inclusions instead of across the boundaries, thus ignoring scattering processes at the boundaries of the inclusions. We note that an earlier work 38 has shown that the above conditions are satisfied in graphene.
The probability distributions P mono and P bi are characterized by the average layer densities n A,P , the root mean square density fluctuations n A;P rms , and the interlayer correlation coefficient η which ranges between −1 for perfect negative correlation of puddles and 1 for perfect positive correlation. The last three quantities n A rms , n P rms and η are parameters of the theory. The n A;P rms values are typically extracted from experimentally obtained intralayer conductivities using Eqs. (2). In the case where a microscopic theory has been formulated for the puddle correlations between the layers 39 , the latter information is already sufficient to determine the correlation function η as a function of the experimental phase space comprised by average densities (n A , n P ) and temperature. In the absence of a microscopic theory for puddle correlations, η may be extracted from experiment using a single slice of drag resistivity in density space (n A , n P ) at fixed temperature and a second slice of drag measured as a function of temperature at fixed density (preferably low density since correlations weaken with increasing temperature). This leaves the rest of the large parameter space available for independent fit parameter-free verification of the theory. The drag EMT formalism may be used as a probe of the depth and the correlations of puddle fluctuations in future applications of 2D van der Waals heterostructures 40 .
Drag in graphene. We consider 'normal' drag (i.e., without exciton condensation) in graphene layers and show that a number of distinct features observed in the experiment 8 can only be explained using drag EMT. For simplicity, we restrict our considerations to temperatures below 200K so that virtual phononinduced enhancements to the drag resistivity that become noticeable starting at about 150 K 41 may be neglected while still maintaining a reasonable level of accuracy. The microscopic theory in the disorder-scattering dominated regime for drag between graphene layers ignoring charge puddles and considering screened interlayer interactions up to second order has been worked out in several papers [31][32][33][34][35][36] and will henceforth be referred to as the standard momentum drag theory or 'standard' theory for brevity. We show here that this theory shows striking differences with the experimental data of Gorbachev et al. 8 performed on graphene encapsulated in hexagonal boron nitride (hBN). These discrepancies do not imply that standard momentum drag theory 31-36 is wrong. Instead, they arise because the samples of Gorbachev et al. 8 possess significant charge puddle fluctuations whereas the standard theory applies to spatially homogeneous systems.
We note that theories for Coulomb drag extending into the carrier-carrier collision dominated (i.e., hydrodynamic) limit have also been developed [42][43][44][45] . We focus here on the disorderdominated regime and begin by using the standard theory [31][32][33][34][35][36] (see Methods for detailed expressions) to calculate the drag resistivity ρ D of graphene encapsulated in hBN. We take into account enhancements to ρ D due to random dielectric inhomogeneity 46 known to occur in the encapsulating hBN substrate 47 . The particular device studied in Fig. 3a of Gorbachev et al. 8 displays dielectric inhomogeneity enhancement by a factor of about 3.6 based on a fit to the experiment (see details in Methods) and we denote the enhanced drag resistivity asρ D 3:6ρ D . We note that this work may be modified to include various effects such as dielectric spatial anisotropy 41 and frequency dependence 34 .
In Fig. 3a, we plot drag resistivityρ D calculated using standard theory as a function of oppositely matched layer densities at several temperatures. These curves compared to the experiment of Gorbachev et al 8 . shown in Fig. 3b show qualitative differences in addition to a general quantitative disagreement. First, the experiment sees the peaks at finite density (i.e., the 'outer' peaks) increase as temperature increases, whereas standard theory predicts that they decrease. Next, a strong positive central peak at the DNP is seen in the experiment, while standard theory predicts zero drag due to electron-hole symmetry. Last, standard theory predicts drag resistivity peaks bunched up relatively near the DNP which rapidly decay to zero as density increases. The experiment shows a more smeared out structure with the outer peaks occurring at noticeably higher densities and decaying more slowly to zero. To give a sense of the behavior of drag resistivity throughout density space, we plot in Fig. 3c a color map ofρ D calculated as a function of layer densities. We note that a similar plot results for drag in the hydrodynamic regime 43 . To demonstrate another difference between standard theory and experiment, Fig. 3d shows the contour lines ofρ D in the density space. Each contour line comes in the form of a triangle with two sides lying almost atop the axes and the third 'outward' line being rather concave. Gorbachev et al. 8 point out that this disagrees with the contour lines obtained experimentally which take the form of triangles with two sides slightly shifted away from the axes, and with the 'outward' line being noticeably straighter (see Fig. 1c and Fig. 5 in the supplementary material of Gorbachev et al. 8 ).
Existing works. Before presenting our drag EMT calculations, it is useful to briefly review existing theoretical literature to place our findings in context. The problem of the central drag peak has drawn three explanations. First, it has been theoretically shown by Schütt et al. 43 that a third-order interlayer interaction term gives rise to a nonzero drag resistivity at the DNP. This however cannot be the explanation of the experimental data because third order drag remains finite (see Fig. 3 of Schütt et al.) as temperature goes to zero whereas the peak in experiment 8 vanishes in this limit. The more 'popular' explanation is the energy drag mechanism 30 that gives rise to positive drag at the DNP due to interlayer energy exchange 48 in the presence of positively correlated puddles (η > 0)  8 . c Drag resistivity calculated as a function of layer densities at T = 70K. d Contour lines for the plot in c. Contour lines are drawn every 0.5 Ω from ±1 Ω to ±3 Ω. Interlayer spacing d = 9 nm throughout. Scale bars representρ D in units of Ω caused by charged impurities (a similar effect is mentioned in Schütt et al. 43 ). Since Song et al. 30 neglects the sizable negative contribution of momentum drag in the presence of positively correlated puddles at the DNP while taking them into account only in energy drag, one cannot definitely conclude that the latter is the cause of the central peak.
The last proposed explanation is by Gorbachev et al. 8 who argue that local sheet corrugations (known as ripples) give rise to spatially varying potentials which induce puddles 25 . At the DNP, the oppositely charged electron and hole puddles in the layers attract one another, creating negative interlayer correlations (η < 0) which lead to a positive peak at the DNP since momentum drag between oppositely charged carriers is positive (see upper left and lower right quadrants of Fig. 3c).
Numerical results. Having shown that the zero-field data of Gorbachev et al. 8 is still not well understood, we present our drag EMT calculations in Fig. 4. The axes n A and n P refer to the average layer densities deduced from gate voltages or Hall resistivity measurements in experiment. Figure 4a shows drag along the line of oppositely matched densities in the presence of negatively correlated puddles and agrees remarkably well with experiment (see Fig. 3b) while Fig. 4b shows the same but in the presence of positively correlated puddles. The agreement of the experiment with Fig. 4a supports the presence of negatively correlated puddles as argued in Gorbachev et al. 8 . The experimental increase of the outer peaks with temperature and the drag peaks at the DNP (see Supplementary Fig. 2 for temperature dependence) are both successfully reproduced in Fig. 4a with an almost quantitative level of agreement with experiment. The nonzero drag at the DNP is possible because electron-hole symmetry is broken by the puddles in the presence of nonzero correlations. Puddles also cause the outer peaks to move out to densities where This general smearing of the drag resistivity throughout the density space is demonstrated in Fig. 4c. Figure 4d shows the contour lines of drag EMT resistivity calculated as a function of n A and n P . The previously concave contour lines have straightened and the other lines have moved away from the axes, in agreement with experiments 8,9 (see Supplementary Fig. 3). The above calculations were performed by substituting the standard theory expressions [31][32][33][34][35][36] (see Methods section) into the EMT equations (1) and (2) and numerically evaluating the result with the layer densities modeled by the bivariate normal distribution in Eq. (7) with average densities n A , n P and root mean square fluctuations n drag in our calculations and discuss this later. As mentioned earlier, n A;P rms should be determined from measurements of intralayer conductivity. Since such measurements are not available in the experiment of Gorbachev et al. 8 , we have chosen n A;P rms ¼ 6 10 10 cm À2 by fitting to the experimental drag resistivities at the outer peaks, where correlations are assumed to be weak so η may be set to zero during this fit. These n rms values are comparable with typical measured density fluctuations in graphene on hBN 23,22 . In all panels of Fig. 4 with the exception of panel (b), we made use of the negative correlation function η given by Eq. (16). Since there is currently no microscopic theory available for arbitrary interlayer correlations, we used a form of η which we believe to be reasonable. First, η was required to decrease in magnitude with increasing density and temperature since screening (which suppressess correlations) in graphene increases with these two quantities. The functional form of η in Eq. (16) was chosen to satisfy this requirement and its numerical prefactors were obtained by fitting to two slices of drag resistivity in density and temperature parameter space (see Methods for more details on this fitting procedure). Since we have only used three lines in the parameter space of n A , n P and T for our fits, this leaves the rest of the large three-dimensional parameter space for fit parameter-free verification of our predictions. q . This suggests that cleaner samples in future will show outer peaks that decrease with increasing temperature. We stress that the temperature trend of the outer peaks depends only on the magnitude of n rms and not on η (see Supplementary Fig. 5 for the uncorrelated η = 0 case).
The above changes from homogeneous behavior are explained by noting that the EMT calculation basically averages (in a manner prescribed by the laws of electrostatics) the homogeneous drag resistivity over density space n ′ A ; n ′ P À Á with the center of the averaging region given by the average densities (n A , n P ) and the size of the averaging region given by the root mean square fluctuations n A;P rms . In the case of uncorrelated puddles η = 0, the averaging region is circular in shape, while in the case of negatively (positively) correlated puddles η < 0 (η > 0), the averaging region is an ellipse with its major axis parallel to the line n ′ P ¼ Àn ′ A (n ′ P ¼ n ′ A ). It is easily verified that at the DNP, negative (positive) correlations involve an average involving mainly regions of positive (negative) drag resistivity (see for instance Fig. 3c). This explains the sign of the drag at DNP in Fig. 4a, b. The change in temperature dependence of the outer peaks occurs because the homogeneous peaks at low temperature are narrower in width (see Fig. 3a) and thus more strongly diminished by the averaging procedure. A more physical way to understand Fig. 4 is to note that drag increases (decreases) with increasing temperature at fixed density in the low (high) temperature regime 36 . Since the presence of puddles causes the local densities in the layers to saturate in magnitude at a lower bound of about n A;P rms (with the relative signs of n A and n P in these local regions given by η), if n A;P rms > n A j j; n P j j, the system can be in the low temperature regime locally (which determines the drag characteristics) even if n A and n P are in the high temperature regime (i.e., k B T> E A;P F ). This causes the low temperature behavior of increasing with T to occur even at average densities that are low relative to temperature. As temperature goes to high values where k B T ) E ðA;PÞ F;rms , drag in this n A;P rms > n A j j; n P j j regime goes to zero everywhere with increasing temperature because all local regions are in the high temperature regime. We now discuss the implications of our results for the role of energy drag 30 . The energy drag resistivity values reported in Song et al. 30 in fact give contributions to drag resistivity at the DNP that are roughly equal in magnitude but opposite in sign to the central peaks shown in Fig. 4a, b. This implies that momentum and energy drag will cancel each other to give close to zero drag at the DNP, in contradiction with the experiment. There are however reasons suggesting that energy drag (see details in Supplementary Note 4) is far smaller than momentum drag at the DNP and the latter alone causes the large central peak. The fact that our calculations using momentum drag alone agrees well with experiment supports this conclusion (see also Supplementary Fig. 6).
So far, we have used drag EMT to explain features that have already been reported experimentally. We now predict a new unreported feature. Our calculations predict that drag resistivity does not go to zero along the axis lines n A = 0 and n P = 0 in the presence of nonzero interlayer puddle correlations. Instead, there is a transition in the sign of drag as one tunes the average density of one layer while keeping the other fixed at zero. To be specific, we fix average density of the passive layer n P = 0 and tune the average density of the active layer n A from −∞ to +∞ as shown in Fig. 5 while using the negative-valued correlation function η given in Eq. (16). The double sign change occurs because the EMT average in the case of negative correlations involves an average of the homogeneous drag resistivity over a slanted ellipse. Consistent with this, Gorbachev et al. 8 report that the lines of zero drag in density space follow "closely but not exactly" the neutrality points of the individual layers (see Fig. 4c and Fig. 5). We note that if interlayer puddle correlations were positive, measuring drag along this same line would result in Fig. 5 multiplied by an overall minus sign.
Puddles and drag upon exciton condensation. Drag resistivity is expected to jump discontinuously or diverge as temperature approaches T c 18 . This occurs because the drag conductivity in Eq.
(3) approaches the values of the single layer in-plane conductivities. The exact value of T c in graphene is a controversial topic 49,50 , and it is not known whether a finite T c is even possible. Our objective here is not to discuss this subtle question but to investigate the impact of puddles on drag resistivity in a system where T c is assumed to be finite (without the need for a magnetic field) and has already been reached. We describe this hypothetical Here we calculate effective drag resistivity as a function of n A along the line n P = 0 using effective medium theory with negatively correlated puddles. We use the same parameters as in Fig. 4a. We predict this set of crossovers of the sign of drag may be seen in future experiments system using the following physically motivated phenomenological expressions. The monolayer conductivity is given by and the drag conductivity by Here, i = A, P is the layer index, and σ 0 ¼ e 2 h , n 0 = 10 10 cm −2 . B, β and C ( B are phenomenological coefficients, the values of which can be chosen quite arbitrarily to fit the material being considered. The statements that we make here do not depend on their particular values. The expressions above are chosen to satisfy the following physical requirements. When the layer densities are perfectly matched (i.e., equal in magnitude but opposite in sign) and the passive layer is open, the electric fields in the layers should be equal based on a minimum dissipation premise 51 , whereas in the case of a short-circuited passive layer with zero electric field, the layers must contain charge currents of equal magnitude in opposite directions. One easily verifies that Eqs. (4) and (5) satisfy these requirements. Whenever the layers have densities of the same sign, the usual non-excitonic drag that is much weaker than excitonic drag is to be expected. This is represented by the second term with C ( B.
We study the effect of puddles in this model by substituting Eqs. (4) and (5) into the EMT equations and calculating ρ E D as a function of charge density n A = −n P for various n A rms ¼ n P rms in Fig. 6. We have used the particular values B = 5, β = 1 and η = 0 but stress that we find numerically the following statements apply for arbitrary η and all positive β, B and C so long as C ( B and n A;P rms ( n A;P . Figure 6 shows that the infinite drag resistivity at perfectly matched densities gets suppressed by puddles to a finite value that decreases with increasing puddle strength. The inset shows that the magnitude of drag at perfectly matched densities goes inversely as n A;P rms . In actual experiment where several samples of differing quality will typically be studied, a drag resistivity that goes as 1/n rms may be interpreted as a sign of exciton condensation.
Discussion. We have generalized the powerful formalism of effective medium theory to the problem of Coulomb drag between 2D materials possessing puddle fluctuations with arbitrary interlayer correlation. Applying it to graphene drag at zero magnetic field reproduces the striking features observed experimentally that have yet to be fully understood 8 . Our work shows that these features arise due to momentum drag in the presence of negative interlayer puddle correlations.
There is currently strong experimental [52][53][54][55] and theoretical 56,57 interest in indirect exciton formation in transition metal dichalcogenide heterostructures. Some experimental drag studies along this line have already begun in heterostructures of MoS 2 58 , a material in which the existence of puddles has been experimentally verified 59,60 , Since the drag EMT formalism developed here applies to arbitrary 2D materials, we believe that it will prove important in the understanding of these and future drag experiments.
Several venues of investigation open up following the present work. The drag EMT formalism may be generalized to include finite magnetic fields 10 and interlayer thermoelectric effects 61 along the lines of earlier theoretical work on single layer graphene thermoelectricity 62 . The present work and its extensions will also be useful for assessing and designing potential technological devices such as those in Refs 15,63 .

Methods
Charge density distributions. We model the single layer charge density using a Gaussian distribution, where n i without the prime superscript denotes the average charge density of layer i set by the external gate voltage. n i rms is the root mean square density fluctuation about the average caused by charged impurities or corrugations and quantifies the strength of inhomogeneity in the sample.
We model the double layer distribution using the bivariate normal probability distribution P bi n ′ A ; n ′ P À Á P bi n ′ A ; n ′ P ; n A ; n P ; n A rms ; n P rms ; η where the interlayer correlation coefficient η quantifies the charge density fluctuations between the two layers. Mathematically, η is defined by where the angular brackets refer to averaging over the areas of the two layers. A value of η = 1 (−1) corresponds to perfectly positively (negatively) correlated charge density fluctuations within the two layers, while a value of η = 0 corresponds to uncorrelated fluctuations. All forms of interlayer correlation may be modeled in Eq. (7) by choosing the value of η accordingly.
Graphene conductivity expressions. The drag conductivity σ D (n A , n P ) between two sheets at uniform densities n A and n P , respectively can be derived diagramatically 64,65 , or from a Boltzmann equation approach 66 as where T is temperature, and d the interlayer spacer width. V(q, ω,d) is the dynamically screened interlayer Coulomb interaction given by where the double-layer dielectric function ε D is given by  with bare interlayer and intralayer Coulomb potentials V 12 ðq; dÞ ¼ V 21 ðq; dÞ ¼ 2πe 2 expðÀqdÞ=κq and V 11 (q) = V 22 (q) = 2πe 2 /κq, respectively, where κ is the dielectric constant of the material encapsulating the graphene sheets and d is the interlayer spacing. Π i is the dynamical polarizability of layer i as shown in existing theoretical literature 67 . Γ x refers to the x-component of the nonlinear susceptibility in monolayer graphene, where we convert charge density to chemical potential using the methods detailed in Supplementary Note 5. We make use of the dimensionless notation, With this notation, we follow the approach detailed in Narozhny et al. 36 to obtain the expressions Gω;q;μ i À Á ¼ À where τ i is the transport scattering time of layer i, assumed here to be energyindependent. This is done only for computational efficiency. We have performed preliminary numerical checks with energy-dependent scattering time corresponding to charged impurities (as seen for example in Eq. (10) of Hwang et al. 68 ) and found that it makes no qualitative difference to the behavior of drag resistivity. Note that the dimensionless frequency and momenta defined in this work differ from that in Narozhny et al. 36 by a factor of 1/2.
The in-plane conductivity of layer i at uniform density n i is given by where i = A, P, DðEÞ ¼ 2jEj=ðπ h 2 v 2 F Þ is the density of states and f(E) the Fermi function, given by f E; μ i À Á ¼ exp . In the case of completely uniform sheets, we have σ E D ¼ σ D n A ; n P ð Þand σ E A;P ¼ σ A;P n A;P . In the absence of exciton condensation, σ D ( σ A;P , and hence Eq. (3) becomes ρ E D % À σ E D σ E A σ E P , in which case τ A,P drop out of the equation and drag resistivity becomes independent of the scattering mechanism. We assume d = 9nm in all our calculations here, corresponding to the sample in Fig. 3(a) of the Manchester experiment 8 . We assume the dielectric constant of hBN to be κ = 3.5 69 .
Fitting procedure for dielectric enhancement. The dielectric inhomogeneity enhancement factor of 3.6 is obtained by dividing the 130 K drag resistivity seen in the experiment of Gorbachev et al. 8 by that predicted by homogeneous momentum drag theory in the high density regime n A = −n P = 6 × 10 11 cm −2 , where inhomogeneities have a negligible effect on transport.
Fitting procedure for correlation function. The correlation function η used in Fig. 4a, c, d is chosen in the following way. Based on the physical requirement that correlations decrease with increasing density and temperature due to enhanced screening, we consider the following functional form η ¼ F n A ; n P ; T ð Þ Θ ÀF n A ; n P ; T ð Þ ð Þ ; where F n A ; n P ; T ð Þ¼B þ C T À T 0 ð Þ Θ T À T 0 ð Þ ½ e ÀD ffiffiffiffiffiffiffiffiffi ffi Here Θ denotes the Heaviside stepfunction. We determine the numerical coefficients by considering two slices of drag resistivity obtained in the experiment. First, we fit our EMT drag resistivity to the behavior of the experimental drag at the DNP as a function of temperature to obtain B = −1 (i.e., negative correlations), C = 0.005 K −1 and T 0 = 70 K. Next, we fit to the experimental drag resistivity measured as a function of n A along the line n A = −n P at T = 70 K (see Fig. 3(b)) and obtain D = 0.5 × 10 10 cm 2 . This correlation function gives perfect negative correlation between the layer puddles at low temperature and decreases in magnitude as density and temperature increase, reflecting that screening gets stronger with increasing density and temperature. We note that this fit was performed by trying different functions by hand. It is likely that even better agreement with experiment will be obtained using more rigorous optimization procedures.
Code availability. Codes used to generate the data shown here are available at https://github.com/derekhoyh/Coulomb-drag-EMT.