Theory of Coulomb drag for massless Dirac fermions

Coulomb drag between two unhybridized graphene sheets separated by a dielectric spacer has recently attracted considerable theoretical interest. We first review, for the sake of completeness, the main analytical results which have been obtained by other authors. We then illustrate pedagogically the minimal theory of Coulomb drag between two spatially-separated two-dimensional systems of massless Dirac fermions which are both away from the charge-neutrality point. This relies on second-order perturbation theory in the screened interlayer interaction and on Boltzmann transport theory. In this theoretical framework and in the low-temperature limit, we demonstrate that, to leading (i.e. quadratic) order in temperature, the drag transresistivity is completely insensitive to the precise intralayer momentum-relaxation mechanism (i.e. to the functional dependence of the scattering time on energy). We also provide analytical results for the low-temperature drag transresistivity for both cases of"thick"and"thin"spacers and for arbitrary values of the dielectric constants of the media surrounding the two Dirac-fermion layers. Finally, we present numerical results for the low-temperature drag transresistivity in the case in which one of the media surrounding the Dirac-fermion layers has a frequency-dependent dielectric constant. We conclude by suggesting an experiment that can potentially allow for the observation of departures from the canonical Fermi-liquid quadratic-in-temperature behavior of the transresistivity.


I. INTRODUCTION
Electron-electron interactions are a source of coupling between closely spaced nano-electronic circuits. This coupling has commanded a great deal of attention during the past thirty years or so, since it constitutes a potential alternative to the inductive and capacitive couplings of conventional electronics. Early on it was realized 1,2 that "Coulomb mutual scattering" between spatially separated electronic systems provides a mechanism to relax momentum that tends to equalize drift velocities. This intrinsic friction due to electron-electron interactions is modernly referred to as "Coulomb drag" [3][4][5][6][7][8][9] . Early experimental work was carried out by Gramila et al. 10 and by Sivan, Solomon, and Shtrikman 11 in semiconductor double quantum wells.
In these experiments a constant current is imposed on the two-dimensional (2D) electron gas in one of the wells (the "active" or "drive" layer). If no current is allowed to flow in the other well (the "passive" layer), an electric field develops whose associated force cancels the frictional drag force exerted by the electrons in the active layer on the electrons in the passive one. The transresistance ρ D , defined as the ratio of the induced voltage in the passive layer to the applied current in the drive layer, directly measures the rate at which momentum is transferred from the current-carrying 2D electron gas to its neighbor. Coulomb drag is ultimately caused by fluctuations in the density of electrons in each layer since two-dimensional layers with uniformly distributed charge will not exert any frictional forces upon each other 3 .
The study of Coulomb-coupled 2D systems has now been revitalized by advances which have made it possible to prepare robust and ambipolar 2D electron systems (ESs), based on graphene 12 layers or on the surface states of topological insulators (TIs) 13 , that are described by an ultrarelativistic wave equation instead of the nonrelativistic Schrödinger equation.
Single-and few-layer graphene systems can be produced, for example, by mechanical exfoliation of thin graphite 14 or by thermal decomposition of silicon carbide 15 . Isolated graphene layers host massless-Dirac two-dimensional electron systems (MD2DESs) with a four-fold (spin × valley) flavor degeneracy, whereas topologically-protected MD2DESs that have no additional spin or valley flavor labels appear automatically 13 at the top and bottom surfaces of a three-dimensional (3D) TI thin film. The protected surface states of 3D TIs are associated with spin-orbit interaction driven bulk band inversions. 3D TIs in a slab geometry offer two surface states that can be far enough apart to make single-electron tunneling negligible, but close enough for Coulomb interactions between surfaces to be important. Unhybridized MD2DES pairs can be realized in graphene by separating two layers by a dielectric 16 (such as Al 2 O 3 ) or by a few layers of a one-atom-thick insulator such as BN [17][18][19][20] . In both cases interlayer hybridization is negligible and the nearby graphene layers are, from the point of view of single-particle physics, isolated. Isolated graphene layers can be also found on the surface of bulk graphite 21,22 and in "folded graphene" 23 (a natural byproduct of micromechanical exfoliation), or prepared by chemical vapor deposition 22 . We use the term double-layer graphene (DLG) to refer to a system with two graphene layers that are coupled only by Coulomb interactions, avoiding the term bilayer graphene which typically refers to two adjacent graphene layers in the crystalline Bernal-stacking configuration 12 .
DLG and TI thin films are both described at low energies by a Hamiltonian with two MD2DESs 12 coupled only by Coulomb interactions, in the absence of singleparticle tunneling. Coulomb drag between two spatially-separated MD2DESs has recently attracted a great deal of theoretical interest [24][25][26][27][28][29] . The calculations in Refs. [24][25][26][27][28][29] refer to the regime in which both layers are either electron-or hole-doped. The Coulomb drag transresistivity in this case is negative and vanishes like T 2 at low temperatures. Despite the considerable amount of work published on the subject recently [24][25][26][27][28][29] , no clear consensus exists on the dependence of the drag transresistivity in the Fermi-liquid regime on carrier densities in the two layers, on the interlayer distance, and on the dielectric constants of the media surrounding the two layers. The main analytical results obtained earlier by other authors will be summarized below in Sect. II.
The Coulomb drag transresistivity between two MD2DESs in the regime in which one layer is electron doped and the other is hole doped has been recently calculated by Mink et al. 30 . In this intriguing regime, the authors of Ref. 30 have found that ρ D grows logarithmically upon lowering the temperature T towards the critical temperature T c for exciton condensation 31 (condensation of electron-hole pairs in a dipolar condensate).
Coulomb drag between two graphene sheets has been recently measured by Kim et al. 16 . This first experimental study represents an important milestone since the authors of this work have shown that the "strong-coupling" regime, i.e. the regime in which the interlayer distance d is much smaller that the typical separation between two electrons in each layer, is easy to achieve experimentally with two, independently-contacted, graphene sheets 16 . This study has indeed fueled the recent theoretical investigations of Coulomb drag between two MD2DESs mentioned above.
In this Article we present in a pedagogical fashion the minimal theory of Coulomb drag between two spatiallyseparated MD2DESs in the regime in which both layers are either electron-or hole-doped. We will be only concerned with the so-called "Fermi-liquid regime" in which both layers are away from the charge neutrality point. Our theory relies on second-order perturbation theory in the screened interlayer interaction and on Boltzmann transport theory. In this theoretical framework and in the low-temperature limit, we demonstrate that, to leading (i.e. quadratic) order in temperature, the drag transresistivity is completely insensitive to the precise intralayer momentum-relaxation mechanism (i.e. to the functional dependence of the scattering time on energy). This is in contradiction with the findings reported in Refs. 27 and 28. We also provide new analytical results for the low-temperature drag transresistivity in both cases of "thick" and "thin" spacers, correcting in the latter case a mistake contained in Ref. 26. At odds with all the previous literature, our results hold true for arbitrary values of the dielectric constants of the media surrounding the two Dirac-fermion layers. Finally, we present numerical results for the low-temperature drag transresistivity in the case in which one of the media surrounding the two MD2DEs has a strongly-frequency-dependent dielectric constant. We conclude by suggesting an ex- The two layers hosting massless Dirac fermions are located at z = 0 and z = d. In a Coulomb-drag transport setup a constant current flow is imposed in one layer (the bottom one, say). If no current is allowed to flow in the other layer, an electric field develops whose associated force cancels the frictional drag force exerted by the electrons in the bottom layer on the electrons in the top one.
periment with a DLG deposited on SrTiO 3 32 that can pave the way for the observation of departures from the canonical Fermi-liquid quadratic-in-temperature behavior of the transresistivity.
Our manuscript is organized as follows. In Sect. II we report a summary of the main analytical results which have been obtained by other authors. In Sect. III we present the model Hamiltonian and the most important basic definitions. In Sect. IV we present the Kubo formalism approach to the calculation of the drag conductivity, while in Sect. V we present a series of simplifications that lead to the Boltmann-transport expression for the drag conductivity and resistivity. These two Sections do not contain original results but make the paper completely self-contained. Experts can skip Sects. IV-V and go directly to Sects. VI-VII, which contain the most important results of this work and all our original results. In Sect. VIII we summarize our main findings and draw our main conclusions.

II. SUMMARY OF THE MAIN ANALYTICAL RESULTS OBTAINED EARLIER BY OTHER AUTHORS
Despite the large body of theoretical work dedicated to Coulomb drag between two unhybridized MD2DESs [24][25][26][27][28][29] , no consensus appear to exist among different authors. In what follows we summarize the main analytical results that can be found in the existing literature.
1) Tse et al. 24 studied Coulomb drag between two unhybrydized graphene sheets separated by a dielectric by employing Boltzmann transport theory. They neglected the spatial dependence of the dielectric constant in theẑ direction (see Fig. 1) and assumed a momentumindependent scattering time. In the weak-coupling limit, i.e. in the limit in which the interlayer distance d is much larger than the average distance between two electrons in each layer, Tse et al. 24 demonstrated that the low-temperature drag resistivity is given by . (1) In Eq. (1) q TF, is the Thomas-Fermi screening wave vector, which is proportional to k F, , the Fermi wave number in each layer. In the symmetric n 1 = n 2 case the previous equation yields 2) We now move on to summarize the main results by Peres et al. 27 . The authors of this work used Boltzmann transport theory, took into account the momentum dependence of the scattering time and also the spatial dependence of the dielectric constant in theẑ direction. While the main approach followed by these authors is numerical, they do also provide an analytical expression for the Coulomb drag transresistivity in the weak-coupling regime. For an intralayer scattering time that depends linearly on momentum, Peres et al. 27 found 3) Katsnelson 26 used Boltzmann transport theory, took into account, at least partially, the spatial dependence of the dielectric constant in theẑ direction, but did not take into account the momentum dependence of the scattering time. The main results of Ref. 26 are: Note that the weak-coupling result by Katsnelson is in agreement with that of Tse et al. 24 . 4) Hwang et al. 28 used Boltzmann transport theory, took into account the momentum dependence of the scattering time but neglected the spatial dependence of the dielectric constant in theẑ direction. For a momentumindependent intralayer scattering time they find: Note that the weak-coupling result in Eq. (5) differs from the result by Tse et al. 24 .
For an intralayer scattering time that depends linearly on momentum, Hwang et al. 28 By comparing Eq. (5) with Eq. (6), Hwang et al. 28 concluded that the functional dependence of ρ D on n and d is very sensitive to the functional dependence of the intralayer scattering time on momentum. 5) Narozhny et al. 29 have recently presented a systematic study of Coulomb drag between two MD2DESs which is based on perturbation theory in the dimensionless coupling constant α ee = e 2 /(hv), v being the Dirac velocity. The authors of Ref. 29 did not consider the spatial dependence of the dielectric constant along thê z direction and discussed mostly the case of "low doping" (i.e. the regime in which the chemical potential is comparable with or smaller than k B T ). Since the focus of our Article is on the opposite regime, i.e. the Fermiliquid regime, we now provide a short summary of the results of Narozhny et al. 29 in this regime only. The authors of Ref. 29 have demonstrated that the energy dependence of the relaxation time is completely irrelevant in relation with the low-temperature drag transresistivity. In the weak-coupling limit Narozhny et al. 29 found for the drag transresistivity the same doping-and interlayer separation-dependence as in Refs. 24 and 26. As density decreases, Narozhny et al. 29 found that the drag coefficient acquires logarithmic corrections -see Fig. 3 in Ref. 29. In particular, in the limit k F d 1, Eq. (41) in Ref. 29 reads As we will see below, our results agree with those of Narozhny et al. 29 , and represent their generalization to the case of a finite Thomas-Fermi screening length (∝ N f α ee ) and spatially-dependent dielectric constants along theẑ direction.

III. MODEL HAMILTONIAN AND BASIC DEFINITIONS
We consider two unhybridized layers of massless Dirac fermions, each one described at the noninteracting level by the following single-channel Hamiltonian (h = 1): Here v is the bare electron velocity, k is the k · p momentum, α, β are (sublattice) pseudospin labels, and σ αβ = (σ x αβ , σ y αβ ) is a vector of Pauli matrices which act on the sublattice pseudospin degree-of-freedom. The field operatorψ † k,α, (ψ k,α, ) creates (destroys) an electron with momentum k, pseudospin α, and layer index = 1, 2. In each layer, the Hamiltonian (8) can be easily diagonalized by the matrix where ϕ k is the polar angle of the vector k.
The i-th Cartesian component of the current-density operator (i = x, y) in the -th layer is given bŷ The Fermi wave number in the -th layer is defined by where n > 0 is the excess electron density in the -th layer 33 and N f, is a degeneracy factor (N f, = 4 for a graphene layer, accounting for spin and valley degeneracies, while N f, = 1 for a TI surface state).
In the absence of disorder, the Green's function corresponding to the Hamiltonian in Eq. (8) in the imaginary frequency axis is given by the following 2 × 2 matrix where 1 1 σ is the 2 × 2 identity matrix in the sublatticepseudospin representation and µ is the chemical potential in the -th layer. In the zero-temperature limit µ → ε F, = vk F, , where ε F, is the Fermi energy in the -th layer. The disorder-free Green's function in the eigenstate representation reads instead where ε k,λ = λv|k| are Dirac-band energies.
In what follows we will need the following matrix elements: and The two MD2DESs described by Eq. (8) are coupled electrostatically by long-range Coulomb interactions, which are influenced by the layered dielectric environment (see Fig. 1). The coupling Hamiltonian readŝ is the density-operator for the -th layer and V (q) (with = ) is the 2D Fourier transform of the interlayer Coulomb interaction Here For future purposes, we introduce the dynamically screened interlayer interaction U 12 (q, ω), which, at the random phase approximation (RPA) level, is given by 34,35 where is the RPA dynamical dielectric function. In Eq. (21), is the well-known 37-39 densitydensity (Lindhard) response function of a noninteracting MD2DEs at arbitrary doping n . The Coulomb interaction in the = 1 (top) layer is given by while the Coulomb interaction in the bottom layer, V 22 (q), can be simply obtained from V 11 (q) by interchanging 3 ↔ 1 . Eqs. (18), (19), and (22) have first appeared in Ref. 36 and their explicit derivation has been reported in Ref. 26. Notice that in the "uniform" 1 = 2 = 3 ≡ limit we recover the familiar expressions Most of the previous work on Coulomb drag in DLG has assumed this limit, which rarely applies experimentally.
The aim of this Article is to present a theory of Coulomb drag, which is valid up to second order in the dynamically-screened interaction U 12 (q, ω), for the system described by the Hamiltonian Note that we are not including in Eq. (23) any term describing intralayer electron-electron interactions. These can be treated in an approximate fashion by invoking Landau's theory of normal Fermi liquids 35 , i.e. by renormalizing the microscopic parameters of the intralayer FIG. 2: Second-order Aslamazov-Larkin-type diagrams contributing to the Coulomb drag conductivity (finitetemperature Matsubara formalism). Solid lines denote the single-particle Green's function in the presence of disorder. Wavy lines denote the screened interlayer interaction U12 in Eq. (20). Black dots in the triangular portions of the diagrams denote vertices of current operators in the two layers, which are both proportional to the Pauli matrix σ x , say, if one is interested in the longitudinal drag conductivity. Finally, k± = k ± q, k ± = k ± q, iεn,± = iεn ± iωn, and iε n,± = iε n ± iωn.
HamiltonianĤ . For example, in Eq. (8) one can use the renormalized quasiparticle velocity 40 , v , instead of the bare velocity v. Anyway, a treatment of the impact of intralayer interactions on the Coulomb drag transresistivity is well beyond the scope of this Article.

IV. KUBO-FORMULA APPROACH TO THE CALCULATION OF THE DRAG CONDUCTIVITY
Coulomb drag starts at second order in a perturbative expansion for the "drag conductivity" σ D in powers of the interlayer Coulomb interaction. To avoid infrared pathologies (stemming from the long-range nature of the Coulomb interaction) of certain integrals that appear in the theory, Coulomb interactions must be screened: from now on we will work with the dynamically-screened interaction U 12 introduced in Eq. (20) rather than with the bare potential V 12 in Eq. (18).
The two Aslamazov-Larkin-type diagrams 41 that contribute to σ D up to second order in U 12 are depicted in Fig. 2. In this figure, solid lines represent disordered propagators in the sublattice-pseudospin representation, while wavy lines represent the screened Coulomb interaction U 12 . Finally, the vertices (black dots in the triangular diagrams) are current operators, one in each layer, evaluated at q = 0 from the very beginning since we are interested in the drag conductivity in the uniform limit. Both current operators can be taken along thex direction since we are interested in the current response of the "passive" layer in thex direction to an electric field applied in the same direction in the "active" layer, i.e. we are interested in the longitudinal drag conductivity.
Straightforward algebraic manipulations lead to the following compact expression for the sum of the two diagrams in Fig. 2: where β = (k B T ) −1 is the usual thermal factor, Ω m and ω n are bosonic Matsubara frequencies and is the so-called "non-linear susceptibility". Here ε n is a fermionic Matsubara frequency, the index = 1, 2 refers to the layer degree of freedom, and the symbol Tr[. . .] denotes a trace over sublattice-pseudospin indices. We now need to perform the analytic continuation iΩ m → Ω + i0 + in Eq. (24) and then take the d.c. Ω → 0 limit. It is very well known that the analytic continuation must be carried out after performing the sum over the Matsubara frequencies ω n . Let us assume that U 12 (q, z), seen as a function of the complex frequency z, is analytic, i.e. that it has no branch cuts. Under this assumption, the function Γ (q, z, z ) has two branch cuts when z or z are on the real axis and is analytic elsewhere, including the points z = 0 and z = 0. We then introduce and perform the Matsubara sum in Eq. (24) using a standard procedure: where n B (z) = (e βz − 1) −1 is the Bose-Einstein occupation factor on the complex plane and C is an appropriate contour which excludes the branch cuts of the integrand (see Fig. 3), located at m(z) = 0, −iΩ m . Notice that the contour C as defined in Fig. 3 includes the points z = 0 and z = −iΩ m , where the integrand is analytic. We obtain since only the integrals around the branch cuts contribute to J (iΩ m ), while the integral over the circle vanishes when its radius is sent to infinity. Here ω ± = ω ± i0 + and P stands for the Cauchy principal value, i.e. the point ω = 0 is excluded from the integration.
We are now ready to perform the analytical continuation to real frequencies. After some algebra we find We finally take the d.c. limit Ω → 0 in Eq. (29) and we obtain since the terms in the last line of Eq. (29) vanish at least like Ω 2 . Thus the drag conductivity in the d.c. limit reads 5 This is the central result of this Section and was first obtained in the context of Coulomb drag between ordinary 2D parabolic-band electron gases -see e.g. Ref.

5.
A. The non-linear susceptibility Let us go back to the definition of the non-linear susceptibility Γ given in Eq. (25). We introduce the following definition so that Eq. (25) reads We here remind the reader that ω 1 and ω 2 (ε n ) are bosonic (fermionic) Matsubara frequencies. Let us now concentrate on the first term on the r.h.s. of Eq. (33), i.e. on J g (q, iω 1 , iω 2 ). For the sake of simplicity, in what follows we omit to indicate explicitly the q-dependence of the functions g (q, iε, iε , iε ) and J g (q, iω 1 , iω 2 ) every time we write equations that involve them. The complete notation will be restored only at the end of the mathematical manipulations that we perform on Eq. (33). We first perform the fermionic Matsubara sum in the first term on the r.h.s. of Eq. (33) by the standard pro-cedure: where n F (z) = (e βz + 1) −1 is the Fermi-Dirac occupation factor on the complex plane. The function g(z, z+iω 2 , z+ iω 1 ) has three branch cuts in the complex plane located at m(z) = 0, −iω 1 , −iω 2 . Choosing a suitable contour of integration C, which encircles only the poles of n F (z) and leaves outside the branch cuts (in complete analogy with the contour drawn in Fig. 3), we find Here ε ± = ε ± i0 + and we used that n F (ε + iω m ) = n F (ε) if ω m is a bosonic Matsubara frequency. Notice that we are allowed to remove the index ± from ε ± when it is summed to iω 1 or iω 2 (since they will be analytically continued to ω ± ) but not when it is summed to iω 1 −iω 2 . Analytically continuing iω 1 → ω + 1 and iω 2 → ω − 2 [to obtain e.g. Γ 1 (q, ω + , ω − ) in the integrand on the r.h.s. of Eq. (31)] we finally obtain We observe that the second term on the second line and the first term on the third line of Eq. (36) are exactly zero since they involve products of three Green's functions with poles on the same half of the complex plane. After some algebraic manipulations we get We now take the limit ω 1 , ω 2 → ω and recast Eq. (34) in the following form: where we have introduced the retarded (advanced) Green's function G R(A) (k, ε). The second term on the right hand of Eq. (33) can be treated in an analogous manner. The final expression for the non-linear susceptibility reads Eq. (39) together with Eq. (31) represent the most important results of this Section and are the starting point for the calculation of the drag transresistivity ρ D .

V. BOLTZMANN-TRANSPORT LIMIT
Eqs. (31) and (39) need to be simplified for any practical purpose and for a quantitative estimate of the Coulomb drag transresistivity.
We first switch from the off-diagonal sublatticepseudospin representation of the Green's functions to the diagonal representation. In the latter representation Eq. (39) reads where we used the definitions in Eqs. (14)- (15) for the density and current vertices.
We now show that the terms with λ = λ do not contribute to Γ (q, ω + , ω − ).
Let us indeed consider all the terms in the triple sum in Eq. (40) in which λ = −λ ≡λ . We first recall that σ x k+qλ ,k+qλ = iλ sin(ϕ k+q ) and We then consider the term on the third line of Eq. (40) and perform the following changes of variables: k → −k and ε → −ε. Using the fact that sin(ϕ −k−q ) = − sin(ϕ k+q ) and sin(ϕ −k − ϕ −k−q ) = sin(ϕ k − ϕ k+q ), we can write the off-diagonal (OD) contribution (λ = −λ ) to Γ OD (q, ω + , ω − ) as Since n F (z) + n F (−z) = 1, the term in square brackets is identically zero. The OD contribution to Γ (q, ω + , ω − ) thus vanishes. From now on we can set λ = λ in Eq. (40) and sum only over λ and λ . We can further simplify Eq. (40) by assuming that, in the presence of weak disorder, the Green's function in the -th layer can be well approximated by the expression, where τ (k) represents a momentum-dependent scattering time in the -th layer and ξ ( ) k,λ ≡ λv|k|−µ are Diracband energies measured from the chemical potential µ of the -th layer.
Note that here τ (k) should not be interpreted as a one-particle scattering time, but rather as the transport scattering time 6 . This statement can be justified within the Kubo formalism by including disorder-related vertex corrections, or, in a much more transparent way, by using the Boltzmann transport equation, see Appendix A.
Below, we will assume that the transport scattering time τ (k) is isotropic, i.e. that it depends only on k = |k|, but allow different scattering times in the two layers. The specific functional dependence of τ on k depends on a particular model of intralayer impurity scat-tering. In the simple case of a momentum-independent scattering time we recover the usual "relaxation time approximation". A scattering time which depends linearly on momentum, τ (k) ∝ k, is a very popular model in the graphene literature. Early on it was understood 42 that such functional dependence of τ on k is needed to explain the linear-in-carrier-density d.c. conductivities that are experimentally measured in samples on dielectric substrates such as SiO 2 (and is typically attributed to charged impurities located close to the graphene sheet). Below, we will not assume any specific functional dependence of τ on k. Our low-temperature analytical results for the Coulomb drag transresistivity do not depend on the particular scattering model that one chooses.
At this point, it is useful to introduce the one-particle spectral function corresponding to Eq. (43), and the identity Using these definitions in Eq. (40) we obtain the following approximate expression for the non-linear susceptibility: From now on we will introduce the simplified notation J x k,λ ≡ vσ x kλ,kλ = vλ cos(ϕ k ), where in the last equality we have used Eq. (15). The definition of J x k,λ should not be confused with the definition of the current-density operatorĴ i q, in second quantization given in Eq. (10). In this Article we are interested in the limit in which intralayer scattering is weak, which is, for example, the most relevant regime for high-quality DLG samples. In this limit we can approximate the spectral functions in Eq. (46) with δ functions, Straightforward algebraic manipulations of Eq. (46) with the use of Eq. (47) yield the following Boltzmanntransport (BT) expression for the non-linear susceptibility: We have dubbed Eq. (48) "BT expression" for the nonlinear susceptibility since after inserting it in Eq. (31) one obtains precisely an expression for the Coulomb drag conductivity that can be derived within a Boltzmannequation approach. This will be shown in Appendix A. Note that, in the usual "relaxation-time approximation" (τ independent of k), the BT expression for the non-linear susceptibility is simply proportional to the intralayer scattering time, Γ BT (q, ω) ∝ τ .
A. Coulomb drag transresistivity in the Boltzmann-transport limit We now focus our attention on the quantity which is actually measured in experiments, i.e. the drag transresistivity ρ D , which is precisely the ratio between the voltage drop in the passive layer and the current in the active layer is indeed ρ D . This quantity can be easily found by inverting the 2 × 2 conductivity matrix, where σ is the intralayer conductivity and the last approximation in Eq. (49) holds true only if σ D σ . Within BT theory the intralayer conductivity at finite temperature is given by 43 while the drag conductivity σ D is given by Eq. (31) with the expression (48) for the non-linear susceptibility. The final expression for the BT drag resistivity reads Note that in the relaxation-time approximation the BT drag transresistivity does not depend on the transport scattering times in the two layers (even if these are different), since, as noted earlier in Sect. V, Γ BT is proportional to τ in this approximation and so is σ .
The situation turns out to be much more complicated in the case in which τ is an arbitrary function of k = |k|. A full numerical treatment of the three-dimensional integral in Eq. (51) at any finite temperature is beyond the scope of the present work and will be the subject of a forthcoming publication.
In the next Section, however, we demonstrate that a dramatic simplification occurs in the low-temperature limit. We will indeed prove that when k B T is much smaller than the Fermi energies in the two layers, k B T min (ε F, ), ρ BT D is insensitive to the precise functional dependence τ = τ (k) of the intralayer scattering time. Since for typical values of doping the Fermi energy in a graphene sheet is very large (≈ 1400 K for a carrier density on the order of 10 12 cm −2 ), analytical results in the low-temperature limit are very useful when both layers are sufficiently away from the charge neutrality point.
In the intriguing regime in which one layer (both layers) lies (lie) at the charge neutrality point BT theory is inapplicable. Coulomb drag in this regime is dominated by unavoidable electron-hole puddles 44 and is certainly very interesting but well beyond the scope of the present Article. For reasons of symmetry, lowest-order Boltzmann transport theory applied to the regime in which one of the two layers is at the neutrality point gives ρ BT D = 0 24,25 .

VI. THE LOW-TEMPERATURE LIMIT
In this Section we use Eq. (51) to calculate ρ BT D analytically in the low-temperature limit, for arbitrary values of the dielectric constants i in Fig. 1, and for a generic scattering time τ (k).
As already anticipated above, the low-temperature regime is readily identified by the inequality k B T min (ε F, ). In this limit we can: (i) use the low-temperature expression for the intralayer conductivity 42,43 ii) evaluate the non-linear BT susceptibilities Γ BT 1 (q, ω) and Γ BT 2 (q, ω) at zero temperature and only to lowest order in ω in the low-frequency ω/ min (ε F, ) → 0 limit; and iii) replace the dynamically-screened interlayer interaction with the much simpler statically-screened interaction U 12 (q, 0). Indeed, the thermal factor sinh −2 (βω/2) in Eq. (51) represents an effective cut off on the values of ω in the integral in Eq. (51): it must be ω 2/β.
To find the asymptotic behavior of Γ BT (q, ω) in the limit ω → 0 we first re-write Eq. (48) in the form Clearly Γ BT (q, 0) = 0. We thus need to calculate the derivative of Γ BT (q, ω) with respect to ω at ω = 0 and in the zero-temperature limit. We find where the last equality is valid only in the zero-temperature limit. More explicitly, we find Eq. (55) is one of the most important results of this Section. An explicit expression for lim T →0 Γ BT (q, ω → 0) will be given below in Eq. (58). Notice that the two δ-functions in the above expression impose that only intraband excitations (λ = λ = +1) contribute to Γ BT (q, ω → 0) in the low-temperature limit. Moreover, the two δ-functions together pin the absolute values of k and k + q to be equal to k F, . This implies that in the low-temperature limit the transport scattering times τ (k + q) and τ (k) inside the square bracket in Eq. (55) must be both evaluated on the Fermi surface of the -th layer, i.e. |k + q| = |k| = k F, . Since τ (k) depends only on the absolute value of its argument (and not on the polar angle of k), both τ (k + q) and τ (k) can be factorized out of the integral in Eq. (55). This is precisely the reason why at the end of Sect. V A we claimed that in the low-temperature limit ρ BT D is completely insensitive to the precise intralayer scattering mechanism.
More mathematically, we have: We now employ the following identity, which can be easily derived by using the aforementioned condition |k + q| = |k| = k F, . The two integrals in Eq. (56) can be easily carried out analytically: we find where we have used that q x = q cos(ϕ q ). Eq. (58) is one of the most important results of this Section. We stress again that the applicability of this equation is not limited to the case of a momentum-independent scattering time. Rather, Eq. (58) applies to a generic intralayer scattering time τ = τ (k). The physical meaning of this result is that in the low-temperature limit the non-linear susceptibility is determined by what happens close to the Fermi surface. What matters thus is only the magnitude of τ (k) evaluated at the Fermi momentum k F, of theth layer. This conclusion is in agreement with Narozhny et al. 29 .
Using Eq. (58) and the following integral we finally arrive at the desired result for the lowtemperature BT drag transresistivity: where q max ≡ min(2k F,1 , 2k F,2 ). Once again, we stress that the result in Eq. (60) does not depend on the precise functional form of τ (k).

A. Dimensionless variables
It is now useful to introduce dimensionless variables. We scale the wave number q in Eq. (60) with k F,1 k F,2 by introducing x = q/ k F,1 k F,2 and the effective interaction U 12 with e 2 / k F,1 k F,2 by introducingŪ 12 = U 12 k F,1 k F,2 /e 2 .
In these reduced units Eq. (60) reads (for physical rea-sons we restore Planck's constant from now on) where x max = min(2 k F,1 /k F,2 , 2 k F,2 /k F,1 ). In Eq. (61) we have introduced the dimensionless coupling constant α ee = e 2 /(hv), which has a value ≈ 2.2 in DLG and To proceed further, we write in an explicit manner the statically-screened dimensionless interlayer interac-tionŪ 12 (x, 0): where we have introduced the crucially-important dimensionless parameter and the dimensionless form factors with The only ingredient we have used to write Eq. (62) is the well-known behavior of the static Lindhard function of a doped MD2DES for wave numbers q smaller than twice the Fermi momentum k F, : Two asymptotic limits can now be easily inferred from the general low-temperature theory in Eqs. (61)-(67): 1) the weak-coupling ("large interlayer distance" and/or "high density") limit, i.e. the limit in which ξ 1, and 2) the strong-coupling ("small interlayer distance" and/or "low density") limit, i.e. the limit in which ξ 1.
As discussed briefly above, the applicability of BT theory in the low-density regime is questionable: for the BT theory to remain valid, one should think of entering the strong-coupling regime by reducing the interlayer separation d while keeping fixed the factor k F,1 k F,2 . Note also that, in the strong-coupling regime, diagrams of higher order with respect to those presented in Fig. 2 might become relevant. The applicability of RPA -Eq. (21) -is also questionable in this regime. Only a comparison with experimental results can tell us about the importance of these subtle issues of many-body theory.
B. The weak-coupling limit In the limit ξ 1, straightforward algebraic manipulations of Eqs. (61)-(67) yield the following result for the low-temperature BT drag transresistivity: where we have introduced a new reduced variable, y = xξ. The quadrature with respect to the variable y can be carried out analytically yielding the following result: where ζ(x) is the Riemann zeta function and ζ(3) ≈ 1.2. Note that this quadrature does not depend neither on 1 nor on 3 .
In summary, we find The dependence of ρ BT D on temperature, densities and interlayer distance displayed by Eq. (71) is in agreement with the results by Tse et al. 24 and by Katsnelson 26 , even though the numerical prefactor we find is four times smaller (we remind the reader that for DLG N f,1 = N f,2 = 4). It is important to observe that the weak-coupling and low-temperature BT Coulomb drag transresistivity is sensitive only to the dielectric constant 2 between the two MD2DESs (see Fig. 1).
C. The strong-coupling limit In the limit ξ → 0 we can expand all the functions f ij (xξ) that appear in the effective interaction (62) is powers of their argument since x is bounded from above, 0 ≤ x ≤ x max , and ξ → 0. Following this procedure we find the following asymptotic expression forŪ 12 (x, 0) in the limit ξ → 0: where we have introduced the Thomas-Fermi screening wave number q TF, = N f, α ee k F, . Using this result in Eq. (61) we find the final expression for the low-temperature BT drag transresistivity in the strongcoupling limit: Eq. (73) simplifies considerably if we assume equal densities and degeneracies in the two layers: n 1 = n 2 ≡ n, k F,1 = k F,2 ≡ k F , ε F,1 = ε F,2 ≡ ε F , and N f,1 = N f,2 ≡ N f . In this case Eq. (73) reduces to where (75) Note that in the strong-coupling limit lim T →0 ρ BT D does not depend on the interlayer distance d and scales like 1/n. As far as the dielectric constants i are concerned, the strong-coupling and low-temperature BT Coulomb drag transresistivity depends only on 1 + 3 but not on 2 .
The quadrature in Eq. (75) can be easily carried out analytically. The result is Note that F( 1 + 3 , N f α ee ) diverges logarithmically in the weak-screening N f α ee → 0 limit, in agreement with Eq. (41) of Ref. 29.

VII. DEVIATIONS FROM THE LOW-TEMPERATURE QUADRATIC BEHAVIOR DUE TO A FREQUENCY-DEPENDENT SUBSTRATE DIELECTRIC CONSTANT
In this Section we illustrate how deviations from the Fermi-liquid quadratic-in-temperature dependence of the Coulomb drag transresistivity can occur in the situation in which the dielectric constant i of one of the media surrounding the two Dirac-fermion layers has a strong frequency dependence.
Consider, for example, a DLG system on a substrate like SrTiO 3 . Following Ref. 45, we model the dielectric constant of this substrate by a frequency-dependent function of the form At room temperature ∞ = 5.2 , 0 = 310, ω 0 /(2π) = 2.7 THz, and γ/(2π) = 1.3 THz. We would like to understand what is the qualitative role played by the change 3 → 3 (ω) in the Coulomb drag transresistivity. For simplicity, we assume that the two graphene sheets comprising the double layer have the same density n 1 = n 2 = n (k F,1 = k F,2 ≡ k F ). We also assume that the Fermi energy ε F corresponding to n is the largest energy scale in the problem: under this assumption, we are still allowed to expand the BT nonlinear susceptibility Γ BT (q, ω) to lowest order in the small parameterhω/ε F .

Despitehω
ε F , we can have two distinct regimes since 3 (ω) brings in a new frequency scale, i.e. ω 0 : i) ω ω 0 ε F /h and ii) ω 0 ω ε F /h. In the first regime what matters is obviously 3 (ω = 0) = 0 . In the second regime, instead, we have to retain the full frequency dependence of 3 (ω). The existence of this second regime ensures the possibility to observe deviations from ρ BT D ∝ T 2 above a certain temperature scale. The low-temperature Coulomb drag transresistivity for the situation described above is given by Eq. (51) with Γ BT (q, ω) given by the expression in Eq. (58) and σ given by the expression in Eq. (52): Here the notation U 12 (q, 0)| 3→ 3 (ω) means that from the point of view of the electronic system we still have to use the statically-screened interlayer interaction, while from the point of view of the substrate we have to take into account the frequency dependence of 3 through the use of Eq. (77).
The double integral in Eq. (78) can be easily performed numerically and some illustrative results are summarized in Fig. 4. For temperatures T < ∼h ω 0 /k B we find the usual Fermi-liquid quadratic-in-temperature behavior, −ρ BT D /ρ 0 = aT 2 where ρ 0 = h/e 2 and a is a numerical coefficient whose actual value depends on the carrier density n and on the interlayer distance d. However, we clearly see from Fig. 4 that, for temperatures T > ∼h ω 0 /k B , the low-temperature Coulomb drag transresistivity deviates from the canonical Fermi-liquid behavior.
Further deviations from the Fermi-liquid temperature dependence are induced by the explicit (and strong) temperature dependence of the dielectric constant of SrTiO 3 -see, for example, Ref. 32 -which here has been neglected for simplicity.

VIII. SUMMARY OF OUR MAIN RESULTS AND CONCLUSIONS
In summary, we have presented in a complete and as much pedagogical as possible manner the minimal theory of Coulomb drag between two spatially-separated twodimensional (2D) systems of massless Dirac fermions, which are both away from the charge-neutrality point. Our theory relies on second-order perturbation theory in the screened interlayer interaction and on Boltzmann transport theory.
In this well-established theoretical framework, which has also been adopted by other authors in earlier works [24][25][26][27][28][29] , we have clearly demonstrated that the precise functional dependence of the intralayer scattering time on momentum plays absolutely no role in determining the low-temperature Coulomb drag transresistivity in the Fermi-liquid regime. This is in contradiction with the findings reported in Refs. 27 and 28 while it agrees with the conclusions reached by Narozhny et al. 29 .
For two layers with identical degeneracies (N f,1 = N f,2 ≡ N f ) and densities (n 1 = n 2 ≡ n), we find the following results for the low-temperature Coulomb drag transresistivity: 1) In the weak-coupling limit -Eq. (71) -we find The dependence of this result on layer-separation d and doping n agrees with that found in Refs. 24,26 and 29.
2) In the strong-coupling limit -Eq. (74) -we find where the explicit expression of the function F( 1 + 3 , N f α ee ) is reported in Eq. (76). The independence of this result on layer-separation d and the dependence on doping n (∝ 1/n) agree with the findings of Ref. 29, even though the authors of this work captured only the weak-screening N f α ee → 0 asymptotic behavior of the function F( 1 + 3 , N f α ee ) and did not take into account the spatial dependence of the dielectric constant along theẑ direction.
General results for n 1 = n 2 and N f,1 = N f,2 can be found in Eqs. (71) and (74).
Finally, we have shown how deviations from the Fermiliquid temperature dependence can occur when one of dielectric constants of the media surrounding the Diracfermion layers depends strongly on frequency and temperature. Double-layer graphene systems deposited, for example, on SrTiO 3 , a very well-known insulator close to a ferroelectric instability 32,45 , can be used as a testbed for this idea.
In the future, we plan to present a systematic numerical study of Eq. (51) at arbitrary temperatures and to investigate more deeply the strong-coupling limit by studying i) the impact of diagrams of order higher than two in the screened interlayer interaction and ii) beyond-RPA corrections. equation by writing × 1 + λλ cos(ϕ k − ϕ k+q ) 2 1 + λ λ cos(ϕ k − ϕ k −q ) 2 . (A11) To obtain Eq. (A11) we used the identities ∂f (0) (k, λ) ∂ξ ( ) k,λ variables (k ↔ k +q, k ↔ k −q, λ ↔ λ , and λ ↔ λ ) and thus the integral is multiplied by an additional factor 1/2. From Eq. (A19) we can derive the longitudinal drag conductivity σ D , which reads as following: where the function Γ BT (q, ω) coincides with the one defined in Eq. (48) of the main text. In writing Eq. (A20) we have introduced an auxiliary variable ω to disentangle k from k :