Theory of Coulomb Drag in Spatially Inhomogeneous Materials

Coulomb drag between parallel two-dimensional electronic layers is an excellent tool for the study of electron-electron interactions. In actual experiments, the layers display spatial charge density fluctuations due to imperfections such as external charged impurities. However, at present a systematic way of taking these inhomogeneities into account in drag calculations has been lacking, making the interpretation of experimental data problematic. On the other hand, there exists a highly successful and widely accepted formalism describing transport within single inhomogeneous layers known as effective medium theory. In this work, we generalize the standard effective medium theory to the case of Coulomb drag between two inhomogeneous sheets and demonstrate that inhomogeneity in the layers has a strong impact on drag transport. In the case of exciton condensation between the layers, we show that drag resistivity takes on a value determined by the amplitude of density fluctuation. Next we consider drag between graphene sheets, in which the existence of spatial charge density fluctuations is well-known. We show that these inhomogeneities play a crucial role in explaining existing experimental data. In particular, the temperature dependence of the experimentally observed peaks in drag resistivity can only be explained by taking the layer density fluctuations into account. We also propose a method of extracting information on the correlations between the inhomogeneities of the layers. The effective medium theory of Coulomb drag derived here is general and applies to all two-dimensional materials.


I. INTRODUCTION
The effects of electron-electron interactions in transport measurements are usually a small correction to predictions from models of non-interacting electrons. Coulomb drag is special because it is identically zero unless interactions are present [1], making it an ideal experimental probe of electron-electron interactions [2]. A typical experiment measuring drag involves driving a current in one (active) layer and measuring the induced potential drop in a physically separated (passive) layer caused by the Coulomb force as shown in Fig. 1. The corresponding induced electric field is then divided by the current density in the active layer to yield the drag resistivity. Studies of this effect now have a history of almost thirty years. The first experiments [3][4][5] were performed using double layer two-dimensional electronic gases (i.e. GaAlAs heterostructures), followed by a series of associated theoretical works [6][7][8][9]. A subject of special interest in drag studies is the formation and detection of interlayer exciton condensates [10] which hold the potential for application in low power electronics (see Ref. [11] and references therein). Drag measurements are the standard method for detecting exciton condensation since it shows strong signatures in the drag resistivity [12,13].
Recent years have seen a sustained experimental effort to understand Coulomb drag in two dimensional materials such as graphene [14][15][16][17] and its bilayer [18][19][20][21] due to their high level of tunability and the ability to reach unprecedentedly small separations between the lay-ers while still keeping them electrically isolated, both of which are favourable to the formation of exciton condensates. These two-dimensional layers however also come with a drawback. Due to their two-dimensional nature, they tend to possess charge density fluctations [22] that arise either due to external charged impurities [23] or corrugations in the topography of the sheets [24], as shown in Fig. 1. As such inhomogeneity is known to play a role in the single layer carrier transport of two-dimensional materials [23], it is natural to expect that they will also play a role in double layer drag transport. Up till now however, there has not been a method for systematically including the inhomogeneity of the layers in calculations of drag resistivity. The situation is very different when it comes to single layer resistivity, for which there exists a well-known formalism that successfully describes charge and heat transport in an inhomogeneous layer known as effective medium theory (EMT) [25][26][27]. In this work, we close the gap by generalizing EMT for the first time to the case of Coulomb drag between two inhomogeneous sheets and demonstrate the importance of inhomogeneity in drag transport by applying the resulting formalism to two examples. First, we show that in the case of interlayer exciton condensation [10,28], where a divergent drag resistivity is expected at zero temperature [29], charge density fluctuations yield a finite value determined by the amplitude of spatial density fluctuations. Next, we apply the drag EMT formalism to drag between two inhomogeneous graphene sheets. The standard homogeneous theory predicts drag resistivity peaks that decrease as temperature increases, in contradiction with experiment by Gorbachev et al. [16] where the opposite is seen. We show that upon inclusion of density inhomogeneity, the drag resistivity peaks increase with temperature within the range of experiment, thus resolving the contradiction. Gorbachev et al. also report measuring anomalously straight drag resistivity isolevels whereas standard homogeneous theory predicts curved ones. We demonstrate that these straight isolevels are in fact caused by the presence of charge density fluctuations. Lastly, there is an ongoing controversy surrounding the nature of correlations between the layers' fluctuations. As we discuss in detail later, there exist arguments that they are correlated [30], anti-correlated [16] or simply uncorrelated. We demonstrate using drag EMT that it is possible to deduce the nature of the correlations by measuring drag resistivity along different lines in the two-layer density parameter space. The plan of this paper is as follows. Sec. II. presents the derivation of Coulomb drag EMT, of which several example applications will be given in the next two sections. Sec. III investigates excitonic drag in the presence of density fluctuations and Sec. IV studies the impact of these fluctuations on drag between graphene sheets. Sec. V concludes with a discussion of this work and the problems that may be pursued in future based on it.

II. COULOMB DRAG EFFECTIVE MEDIUM THEORY
We consider the standard drag setup -two parallel 2D sheets of identical size separated by some finite distance, with a current flowing through the active layer while the passive layer remains an open circuit. To model the presence of inhomogeneity, we assume that the active and passive layers are each made up of N patches (commonly referred to as 'puddles') each with its own conductivity, σ A i and σ P i respectively, where i = 1, · · · , N . We assume that the puddles in both layers are circles of radius a, with the ith puddle of the active layer lying exactly atop the ith puddle of the passive, as in Fig. 2. Both puddles are assumed to be of equal area. This assumption is of general applicability because we are allowed to define as many circular patches and make them as small as we wish. We do not make any further assumptions about the nature of these puddles and our derivation is applicable regardless of whether the puddles are correlated, anti-correlated or uncorrelated.
Taking the continuum limit, we obtain the final results of Eqs. (22) and (23). We now begin the detailed derivation. Our first task is to solve for the field inside the ith (where i is arbitrary) puddle that has been embedded in the effective medium as described above. We assume that the puddles are regions of uniform 2D polarization in a direction parallel to the effective electric field of the layer. These polarizations are denoted by M A = M A e x and M P = M P e x for We have thus solved for the electric fields within the ith pair of puddles. In order for these fields to be combined with the self-consistency equations in a useful manner however, we must first be able to write the electric field inside each puddle as a function of only its own layer's effective medium field (i.e., instead of being a function of the effective medium fields of both layers). We achieve this by requiring that just given the two effective medium layers without puddles (i.e., Fig. 2 with the puddles taken out), there will be no current flow in the passive layer. Physically, this means that E A 0 and E P 0 represent uniform fields that effectively model the spatially fluc-tuating fields in the two layers. The expression for this condition is given by where we have made use of the fact that i f i = 1. Equation (18) is the well-known discrete single layer EMT equation applied to the active layer. Equation (19) is more complicated and comprises two terms summing to zero. Since the two unknowns σ D E and σ P E cannot be determined by this single equation, we require another condition. Since the drag conductivity is very small compared to the in-plane conductivity within either layer, we may obtain this condition by approximating that the interlayer interaction has a negligible effect on the passive layer conductivity. This then implies that the standard single layer EMT equation (cf. (18) ) applies to the passive layer, and the two terms of Eq. (19) must both be individually zero. Equating the first term to zero and solving for σ D E yields while setting the second term to zero yields The former is a newly derived discrete EMT drag equation while the latter is just the already-known discrete single layer EMT equation applied to the passive layer. Note that in the case of there being only two puddles in each layer (relevant only at double charge neutrality) so that f i = 1 2 , i = 1, 2, we recover the results of Ref. [31] which considered drag between inhomogeneous layers each consisting of only two areal components. Generalizing Eq. (20) to the continuum limit, we obtain the EMT drag conductivity equation where n A and n P denote the charge densities in the active and passive layers and σ E D is the effective medium theory averaged drag conductivity obtained by solving the equation. P bi (n A , n P ) is the joint probability distribution of finding two points on the layers with one point lying directly above the other having charge densities (n A , n P ). Doing the same for Eqs. (18) and (21) where i = A, P denotes the layer index and σ i (n i ) is the homogeneous conductivity of layer i at uniform density n i . P mono (n i ) is the single layer probability density of finding a point on layer i with charge density n i .
Equation (22) is the main result of this work and represents the first generalization of EMT to the drag problem. We emphasize that it is general and applies to drag between any two sheets of two-dimensional material. The drag resistivity is given by solving Eqs. (22) and (23) for the three conductivities σ E D , σ A E , σ P E and inserting them into The standard homogeneous theory [8,9,32] is recovered in the limit of n (A,P) rms → 0. To perform actual calculations, one must choose specific probability distributions for P mono and P bi . We choose for the former the usual monovariate 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 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. A value of η = 1 (−1) corresponds to perfectly correlated (anti-correlated) charge density fluctuations within the two layers, while a value of η = 0 corresponds to uncorrelated fluctuations. Mathematically, η is defined by where the angular brackets refer to averaging over the areas of the two layers. Static charged impurities in the surroundings of two sheets held close together can lead to correlated fluctuations, whereas random strain in the sheets together with strong Coulomb coupling between them can lead to anti-correlated fluctuations. If the separation between layers is fairly large, the fluctuations will tend to be uncorrelated due to each sheet seeing a potential of different origin. All these scenarios may be modeled in Eq. (27) by choosing the value of η accordingly. In an earlier version of this work, we claimed that Onsager reciprocity relation is violated for η = 0 but this has since been found to be due to a numerical error. We show in Appendix B a proof that Onsager reciprocity is obeyed in the drag EMT formalism.
Before moving on to applications, we discuss the physical conditions under which inhomogeneities strongly affect transport and the drag EMT must be applied. There are three energy scales that need to be considered. First, the temperature of the system k B T is important. Second, we introduce the typical layer Fermi energy E F (n) wherē n ≡ √ n A n P is the typical average layer density. Lastly, we define an inhomogeneity energy scale E F (n * ) where n * ≡ √ n A rms n P rms is the typical root mean square fluctuation of density in the layers. Generally speaking, the impact of fluctuations (i.e. the amount by which drag changes after including density fluctuations) is strong when is satisfied. The reason is as follows. The drag EMT essentially yields an average of the homogeneous ρ D over a region of density space centered at average densities (n A , n P ) with an area on the order of (n * ) 2 , with the exact shape and orientation of the averaging region determined by η and the relative magnitudes of n A rms and n P rms . Such a coarse-graining procedure makes a big difference when performed over regions in which the function being averaged possesses turning points. In the case of drag resistivity, these occur at the double neutrality point, and in the regions |Ē F | ∼ k B T where the finite density drag peaks occur, and they will be inside the region of averaging when Eq. (28) is true. Inhomogeneity is negligible if E * F max(Ē F , k B T ) because the averaging encompasses a region over which ρ D does not change much. Lastly, if E * F max(Ē F , k B T ), then ρ E D will be approximately zero everwhere since the averaging window includes the high density regions where drag has already gone to zero for all intents and purposes. Under the right conditions, electrons and holes in the two layers are expected to bind together forming stable bosonic excitons, which can condense into a superfluid exciton condensate (see for instance Refs. [33] and [34]). In particular, the drag resistivity of an exciton condensate diverges as temperature approaches zero [29] because the magnitude of drag conductivity approaches that of the single layer conductivity. We demonstrate that this divergence is suppressed in the presence of charge density inhomogeneity. We model the exciton condensate at zero temperature using the following phenomenological expressions for the various conductivities. The monolayer conductivity is given by and the drag conductivity by In the above, i = A, P is a layer index, and σ 0 = e 2 h , n 0 = 10 10 cm −2 . A and α are phenomenological coefficients that can be given arbitrary values depending on the specific material under consideration. We note that the above expressions produce the correct behavior in various limits. When the passive layer is in open-circuit configuration and no current flows in it, they lead to equal electric fields in the two layers, as expected from Ref. [29]. When the passive layer is short-circuited so that there is no electric field across it, they yield equal charge currents in the two layers.
In the absence of inhomogeneity (ie. n (A,P) rms = 0), it is clear from Eq. (24) that a divergence in ρ D occurs at perfectly matched opposite densities n A = −n P . To investigate the effect of density fluctuations on exciton drag, we substitute the above conductivity expressions into the EMT equations (22) and (23) and use the probability distributions in Eqs. (26) and (27) with various values of n (A,P) rms . As shown in Fig. 3, density fluctuations suppress the divergence in drag resistivity. Furthermore, the magnitude of drag at perfectly matched densities goes inversely as n (A,P) rms as shown in the inset. Here, we have used A = 5, α = 1 and η = 0 but our numerics suggest that these statements still apply for arbitrary values of A and all positive powers α. They also apply regardless of the value of correlation coefficient η. This finding demonstrates that sample inhomogeneity is important when using drag resistivity as a probe of exciton condensation and should be of great interest in the ongoing search for exciton condensation [18,21,35].

IV. DRAG IN GRAPHENE SETUPS
In this section, we demonstrate that just as in single layer graphene transport, inhomogeneity plays an important role in graphene drag transport and it is necessary to take inhomogeneity into account in order to explain the experimental data in the literature. We begin with a brief review of drag calculation.

Review of Coulomb drag theory
Coulomb drag has been studied theoretically in graphene monolayers in several works [32,[36][37][38][39]. The drag conductivity σ D (n A , n P ) between two sheets at uniform densities n A and n P respectively can be derived diagramatically [8,9] as where T is temperature, and d the interlayer spacer width. V (d) is the dynamically screened interlayer Coulomb interaction [40], and Γ x refers to the xcomponent of the nonlinear susceptibility in monolayer graphene as given in Ref. [32]. In our calculations here, we choose parameters based on the experimental setup in Gorbachev et al. [16] where both graphene sheets are encapsulated in hexagonal boron nitride (hBN) and separated by a hBN spacer. Further details on these quantities may be found in Appendix C. In all our calculations, we use the following parameters unless otherwise specified-the interlayer spacing is d = 9nm, distance of the active (passive) layer from the charged impurity plane is 20nm (10nm) and the impurity concentration on the charged impurity plane is n imp = 15 × 10 10 cm −2 .
The in-plane conductivity of layer i at uniform density n i is given by where i = A, P , D(E) = 2|E|/(π 2 v 2 F ) is the density of states and f (E) the Fermi function, given by f (E, and τ i (E) is the intralayer transport scattering time. The chemical potential µ i is determined by n i and T (see Appendix C). We assume that electroncharged impurity scattering is the dominant scattering mechanism and neglect all others. The previous two equations together with Eq. (24) yield drag resistivity at any point in density space (n A , n P ), assuming completely homogeneous graphene sheets.

Temperature Dependence of drag resistivity peaks
Using the above expressions, we calculate drag resistivity along the line of oppositely matched densities n A = −n P and show the results in Fig. 4(a). At each temperature, the drag resistivity follows the standard density dependence. At charge neutrality, it is zero due to electron-hole symmetry. As we move away from charge neutrality, drag increases to a peak (henceforth referred to simply as the 'drag peak') and subsequently goes to zero at high density as screening between the layers effectively decouples them. However, a surprise occurs as we tune temperature. The drag peaks decrease as temperature increases, in direct contradiction with the experimental observations of Gorbachev et al. This puzzling contradiction appears to have gone unnoticed in the literature.
The resolution lies in including inhomogeneity in the calculation of drag using EMT. This is done by substituting Eqs. (31) and (32) into the EMT equations (22) and (23) to obtain the effective conductivities σ D E , σ A E and σ P E , which are then substituted into Eq. (24) to obtain the effective drag resistivity. We assume uncorrelated density fluctuations η = 0 in Eq. (27) because interlayer correlations at finite densities away from the double neutrality point n A , n P = 0 are expected to be weak due to screening. The drag peaks thus obtained by EMT are smaller than their values assuming perfect homogeneity because EMT essentially averages ρ D over some region in density space and this can only result in maxima decreasing in height. The drag peaks also follow a non-monotonic temperature dependence, increasing with temperature to a peak at some temperature set by the inhomogeneity strength, and decreasing thereafter. We choose inhomogeneity strength n (A,P) rms = (7, 14) × 10 10 cm −2 so that the highest peak occurs at T = 240K (i.e. the highest temperature studied in Gorbachev et al.) and show the results of our calculation in Fig. 4(b). At high densities and temperatures such that v F √ πn A,P and k B T are much greater than v F πn (A,P) rms , the homogeneous and EMT drag calculations yield essentially the same values, as expected since the disorder energy scale is now the smallest energy scale of the problem. As a corollary of this, the decrease in drag peak at temperatures above 240K occurs because inhomogeneity becomes unimportant and the homogeneous behavior re-emerges. The decrease in drag peak above 240K is a concrete prediction of our theory that should be testable in existing experimental setups.
The drag resistivities calculated in Fig. 4 are smaller than those of the experiment, possibly due to enhancement mechanisms such as dielectric inhomogeneity [41] and virtual phonons [42] which we have not included in our calculations. We defer a detailed study of these effects to a future work. For now, we take them into account phenomenologically by comparing our theoretically calculated drag resistivity with experiment at high density and temperature (where inhomogeneity plays a negligible role) and using the extracted discrepancy factor to scale up our calculated ρ D values by hand. Comparing the experimental ρ D at n A = −n P = 60 × 10 10 cm −2 and T = 240K (i.e. the largest density and temperature in the data of Gorbachev et al.) with the corresponding ρ D calculated using the homogeneous theory, we find a dis-  crepancy of 3.96. A similar discrepancy was also found by Gorbachev et al. We thus define a 'dressed' drag resistivityρ D as the calculated ρ D multiplied by a factor of 3.96. Comparing the dressed EMT drag resistivity peaks at various inhomogeneity strengths with experiment in Fig.  5, we see that the curve for n (A,P ) rms = (7, 14) × 10 10 cm −2 agrees well with experiment. We choose a larger root mean square density fluctuation for one layer because in general one layer is expected to be nearer the impurity plane than the other. However, we note that our numerics do not show much difference when the individual layer root mean square fluctuations are changed so long as the total fluctuation n (A) rms + n (P ) rms is unchanged.

Drag isolevels in the presence of inhomogeneity
We show in Figs. 6(a) and (b) our calculations of ρ D as a function of n A and n P using homogeneous theory and EMT respectively. The drag isolevels are strongly concave in the homogeneous case, whereas they are essentially straight in the presence of inhomogeneity. This is due to the fact that EMT basically does a form of averaging of the conductivities in density space. At high densities much greater than n (A,P) rms , the contour lines start becoming more concave since the influence of inhomogeneity becomes increasingly unimportant at high densities. These straight contour lines have also been observed in the experimental works of Ref. [16] and [15].

Drag resistivity with correlated inhomogeneities
Various proposals have been made concerning the existence of correlations in the density fluctuations of the active and passive layers. Gorbachev et al. argue for the existence of anti-correlated fluctuations in which each hole (electron) puddle lies predominantly above an electron (hole) puddle. This occurs if the density fluctuations arise from strain-induced corrugations in the graphene sheets [24] and the sheets deform themselves so as to minimize the electrostatic potential energy. On the other hand, correlated fluctuations in which each hole (electron) puddle lies predominantly above another hole (electron) puddle have also been suggested [30]. This situation occurs if the fluctuations arise from charged impurities in the surrounding environment [23] since the puddles in the two layers experience potentials arising from one and the same set of charges. Lastly, the density fluctuations tend toward being completely uncorrelated as the interlayer spacing becomes large.
We study the effect of all three possible types of correlation on ρ D and compare with the homogeneous theory. samples and three different values of correlation coefficient η respectively. Strictly speaking, η changes as a function of density (and temperature) but we assume it to be constant since we are interested mainly in the qualitative effect of correlations that have more to do with the sign of η rather than its exact value. With the exception of the double neutrality point (DNP), inhomogeneity always causes a decrease in magnitude of ρ D , regardless of the nature of correlation. This is because drag EMT performs an averaging of the homogeneous ρ D in density space, as mentioned at the end of Sec. II. This averaging can only lead to an increase in magnitude of ρ D at minima of |ρ D |, and the only one such minima that occurs is located at the double neutrality point [32].
The presence of non-zero correlations also leads to a finite ρ D at the double neutrality point, with correlation (anti-correlation) leading to a negative (positive) drag. This is to be expected since drag between sheets of the same (opposite) sign of charge density is negative (positive). Gorbachev et al. report measuring positive ρ D at the double neutrality point in their experiment and that this positive drag always appeared in the regime of n A , n P ∼ n rms . There are two proposed explanations for this. Song and Levitov [30] have demonstrated that energy exchange between the two layers in the presence of correlated fluctuatons yields positive ρ D due to thermoelectric effects and suggest that this might explain experiment. They refer to this effect as 'energy drag' and to the standard Coulomb force-mediated drag considered in this work (i.e. Eq. (31)) and many others [32,[36][37][38][39] as 'momentum drag'. Song and Levitov ignore the negative contribution of correlated momentum drag at the DNP (cf. Fig. 7(c)). On the other hand, Gorbachev et al. ignore energy drag and explain their experiment by considering only momentum drag at the DNP and arguing (without the use of drag EMT) that anti-correlated fluctuations arising from random strain in the graphene sheets give rise to positive drag (cf. Fig. 7(d)). In the presence of nonzero correlations, momentum and energy drag compete at the DNP and a full analysis involving both must be performed in order to predict ρ D . Such an analysis has not been performed and the exact origin of positive drag at the DNP remains an open problem. The drag EMT in this work does not include energy drag and is thus unable to make quantitative predictions of drag at the DNP. However, it might still be possible to shed some light on the nature of correlations between the two sheets away from the DNP, where n A , n P n rms and energy drag is negligible. This is done by measuring ρ D along the two lines n A = ±n P and comparing the magnitude of the drag peaks (i.e. the peaks at finite density away from DNP) along the two lines. As shown in Figs. 7(c) and (d), drag EMT predicts that correlated (anti-correlated) fluctuations lead to drag peaks of larger magnitude along the n A = n P (n A = −n P ) line than the n A = −n P (n A = n P ) line. Fig. 7(b) shows that uncorrelated fluctuations lead to peaks of equal magnitude along the two lines. This suggests a means of deducing the nature of correlations experimentally. We note that this should be done at low temperatures since inhomogeneity effects become weak at high temperatures.

V. DISCUSSION
Coulomb drag is a favored probe for studying electronelectron interactions but a standard framework for studying drag in the presence of sample inhomogeneity has been absent till this work. We have generalized effective medium theory to include Coulomb drag, thus providing for the first time a systematic means of incorporating spatial density fluctuations into drag calculations. The importance of the formalism was demonstrated through several examples. Standard theory neglecting inhomogeneity suggests that an indirect exciton condensate possesses infinite drag resistivity at zero temperature [29]. We showed that the presence of density fluctuations yields a finite drag resistivity with a value determined by the amplitude of the fluctuations. In the case of drag between graphene sheets, we demonstrated that inhomogeneity is crucial for explaining the experimental observations made in the Manchester experiment of Ref. [16]. We showed that the rise of the drag resistivity peaks with temperature is a direct result of inhomogeneity and cannot be explained by the standard (homogeneous) theory. Our calculations also yield isodrag lines that bear a striking resemblance to those from the experiment. Lastly, there is an ongoing controversy concerning the sign of correlation between the density fluctuations of the two layers. We showed that this may be deduced from drag measurements along different lines in the (n A , n P ) density space.
Several future avenues of research arise from this work. First, we have considered only the Coulomb-mediated momentum exchange between the layers and ignored interesting proposals of possible thermoelectric effects [30,43]. In principle, it is possible to formulate a generalized effective medium theory of drag incorporating thermoelectricity along the lines of Ref. [44]. Second, it would be of great interest to explore the effects of inhomogeneity in double layer graphene in the hydrodynamic regime and also in different drag setups such as bilayer graphene [19,20], topological insulators [45], two-dimensional electron gases [35] and hybrid graphenetwo-dimensional gas setups [46]. Third, the present work may be generalized to obtain a magnetodrag EMT (i.e. an EMT for Coulomb drag in the presence of a magnetic field). Such a theory is at present still lacking in the literature and would be especially important because to date, successful observations of exciton condensation have typically involved non-zero magnetic fields [10,18,21]. Last, it would be interesting to implement the dielectric inhomogeneity [41] and virtual phonon [42] enhancement mechanisms in theoretical calculations to obtain precise agreement with experiment. We have also assumed for the sake of simplicity that η and n (A,P ) rms are effectively constant as functions of density and temperature. This assumption does not affect the qualitative accuracy of the results but nonetheless has a quantitative effect. One possible way to remedy this is to use the rigorous theory presented in Ref. [47] for calculating η and n (A,P ) rms as a function of the system parameters. We leave these as problems for future work. For completeness, we also review the derivation of the monolayer EMT result, Eq. (23). Useful reviews of monolayer EMT may be found in Consider a sheet of 2D material which is made up of a patchwork of N areas (i.e. puddles) each with differing conductivities σ i , where i = 1, · · · , N . The areal fraction of the ith puddle is denoted by f i with i f i = 1. We wish to calculate the effective conductivity of this sheet. We first imagine that each puddle is embedded in a homogeneous effective medium of conductivity σ E and through which permeates a uniform electric field E 0 . Next, we work out the electric field E i inside each puddle in the form Figure 8. A circular region of inhomogeneity embedded inside a homogeneous effective medium of conductivity σE. E i = (· · · ) E 0 . Consider the ith puddle embedded in the effective medium as shown in Fig. 8. For simplicity, we assume that all puddles are circles of radius a. We find the electric field inside this puddle.
The conductivity of the medium is σ E and it is permeated by a uniform field E 0 = E 0 e x corresponding to an external field that we will refer to as the primary field. We denote the field inside the puddle as E i and the field outside as E E . We assume that the puddle possesses a uniform polarization M pointing in the same direction as the external field, so that M = M e x . From electrostatics, such a circular puddle is associated with an electric field in the region outside it given by with a corresponding potential Here, the subscript 's' is short for 'secondary' and we placed the origin of our radial coordinate system at the center of the circular puddle. Hence, the total field outside (r > a) is, taking the e r components, = E 0 cos(θ) + a 2 2 0 r 2 M cos(θ) and the potential is Next, we consider the field inside the puddle. We guess that the field inside is proportional to the externally applied field E 0 . Hence, We now use boundary conditions to solve for the unknown C, and obtain the desired field inside the puddle. The boundary conditions that must be obeyed at the boundary of the puddle are the continuity of potential, and continuity of radial current density. In equations, they are and We can make use of these two boundary conditions to solve for the two unknowns M and C. This yields The fact that solutions for C and M exist validates our earlier guess in Eq. (A6), which we know is the unique physical solution due to the uniqueness theorem. We are not interested in M here. Our objective was to find E i , which we have successfully done by determining C. Explicitly, we have found that where E i is the magnitude of E i . We substitute the above into the EMT self-consistency condition and simplify to obtain Generalizing to a continuum of puddles and denoting the continuous puddle label by n, one obtains Eq. (23).

APPENDIX B: PROOF OF ONSAGER RECIPROCITY
In the context of drag, Onsager reciprocity [48] predicts that there should be no difference in drag resistivity when the active and passive layers are switched.
We prove that Onsager reciprocity is obeyed by the drag EMT derived in this work. Mathematically, proving Onsager reciprocity amounts to interchanging all 'A' and 'P' superscripts (denoted henceforth by 'A ↔ P') and showing that drag resistivity remains the same. Since it is clear that interchanging σ E A and σ E P in the drag resistivity expression Eq. (24) makes no difference, our remaining task is to prove smilar invariance for σ E D . We start by considering the discretized version of σ E D in Eq. (20), which may be rewritten as where and .

(B2)
The numerator Eq. (B1) is invariant under A ↔ P because of the invariance of the standard (homogeneous) drag conductivity, as seen from Eq. (31). To prove that the denominator is also invariant, we take the difference of the monolayer EMT equations (18) and (21), Simplifying and rearranging, we obtain .
(B3) This proves that Eq. (B2) is also invariant under A ↔ P. We have thus established that the discrete version of σ E D obeys Onsager reciprocity. Taking the continuum limit, we conclude that the continuous version Eq. (22) does too and σ E D remains unchanged when the roles of the two layers are swapped. This completes the proof of Onsager reciprocity.
APPENDIX C: HOMOGENEOUS DRAG THEORY

Drag Conductivity Expressions
The dynamically screened interlayer Coulomb interaction is 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. We assume d = 9nm in all our calculations here, corresponding to the sample measured in Fig 3(a) of Gorbachev et al. Π i is the dynamical polarizability of layer i as given in Ref. [40]. The physical meaning of the nonlinear susceptibility is that it is a response function relating the voltage felt across each layer with the current it induces within it [32] via i(ω) = dr 1 dr 2 Γ(ω, r 1 , r 2 )V (r 1 )V (r 2 ).
The x-component of the nonlinear susceptibility of layer i is given by = Γ i (ω, q, µ i /k B T ) cos(θ q ), where we convert charge density to chemical potential using the methods detailed in the next subsection. It is convenient at this stage to switch to dimensionless notation, (C5) With this notation, we follow the approach detailed in Ref. [32] to obtain the expressions Γ i (ω, q, µ i k B T ) = − 4e v FΓ i (ω,q,μ i ), Γ i (ω,q,μ i ) = 1 4π G(ω,q,μ i )q, G(ω,q,μ i ) = K i (z,q,ω), |ω| <q, where τ i is the transport scattering time of layer i. Note however that the dimensionless frequency and momenta defined in this work differ from that in Ref. [32] by a factor of 1/2. In this paper we assume that electron-charged impurity scattering dominates over all other scattering mechanisms. In this case, τ i is given by where q = |k − k | and θ k,k is the angle between the initial and final wave vectors k and k in a scattering event. n (i) imp is the areal concentration of charged impurities in an impurity plane located at a distance d (i) imp away from the graphene sheet. V imp (q) = 2πe 2 /(κq) exp(−qd imp is the distance between the graphene sheet and the charged impurities which are assumed to lie in a single plane. In our calculations, we assume that both layers see one and the same impurity plane so that n A imp = n P imp ≡ n imp . In the experiment of Gorbachev et al., the impurity plane is expected to be near the SiO 2 wafer on which the drag heterostructure rests. Hence, we assume d P imp = 10nm and d A imp = 20nm corresponding to the fact that the passive layer lies nearer the SiO 2 wafer in Gorbachev et al.
The single layer dielectric function s is given by where Π(q, T ) is the static polarizability of graphene. As pointed out in Ref. [32], the nonlinear susceptibility contains a logarithmic divergence along the lineq =ω so long as τ has an E-dependence. In this work, we prevent σ D from diverging in calculations by using the dynamical polarizability of Ref. [40] in the dielectric function. This introduces a divergence in the denominator of σ D which cancels that in the numerator, leaving behind a finite and well-defined quantity. We shall assume the value of the graphene 'fine structure constant' is r s = 0.568, corresponding to an estimated value of κ = 3.5 for hexagonal boron nitride. The above expressions constitute all the ingredients one needs to calculate the drag homogeneous conductivity.

Relation between Charge Density and Chemical Potential
Here we review the one-to-one correspondence between charge density and chemical potential given a fixed temperature. The charge density is defined as n = n e − n h where n e and n h refer to the electron and hole densities respectively. These densities are obtained from the chemical potential µ and temperature T via where Li 2 refers to the dilogarithm function. This equation allows us to find the charge density given the chemical potential and temperature. We can also find the chemical potential given the charge density and temperature. This is done by noting that E F = v F π|n| sign(n) and using the relation where