Effect of quasiparticle excitations and exchange-correlation in Coulomb drag in graphene

Coulomb drag in double layer graphene systems separated by an h-BN interlayer allows probing of the electron-electron interactions in the effective limit of zero layer separation. Although these interactions can be influenced by plasmons, phonons and exchange and correlation effects, these excitations have never been studied altogether, missing the effects of their coupling on the drag physics. Here we study theoretically the effects of these quasiparticles and their coupling, including also the effects of the electronic exchange and correlation, and demonstrate that the drag resistivity can attain a maximum value at room temperature and beyond, where hybridized plasmon-phonon modes contribute significantly. In particular, the hybridization of the plasmons with the hyperbolic phonons of h-BN, confined within the reststrahlen bands, enhance the drag resistivity. This study paves the way for the exploration of novel many-body physics phenomena in systems coupled through emerging 2D hyperbolic materials. Coulomb drag describes a phenomenon where long-range Coulomb interactions occur between charge carriers in two electrically isolated systems, and a current applied to one system can induce a current (or voltage) in the other. Here, the authors theoretically investigate the contribution of plasmons, phonons, and exchange-correlation to Coulomb drag in a graphene/h-BN/graphene heterostructures.

C oulomb drag is a phenomenon which occurs in two 2dimensional (2D) electrical conductors placed close to each other but still electrically isolated, such that an electric current in one conductor, called active layer, drags carriers inducing a closed-circuit current (or an open-circuit voltage) in the other conductor, called passive layer. This phenomenon is caused by the transfer of momentum between carriers in the different layers due to the interlayer electron-electron interactions [1][2][3][4] . Coulomb drag is of special interest because it is identically zero unless these interactions are present, thus becoming an ideal experimental probe. This tool provides a detailed insight into the many-particle interactions in low-dimensional systems, which play a central role in a wide range of condensed matter phenomena like the fractional quantum Hall effect, the hightemperature superconductivity, the Wigner crystallization, the exciton condensates, and the Mott transition. The first experimental demonstration of Coulomb drag was carried out using 2D electron gases (2DEGs) in GaAs double quantum well structures 5,6 and, during two decades, they were the most used material systems for experiments. This situation changed though with the recent experiments in double-layer graphene systems 7 , specially those using few-layer hexagonal boron nitride (h-BN) as the insulating interlayer 8,9 , and the intense theoretical effort that has accompanied them [10][11][12][13][14][15][16][17][18][19][20][21] .
Graphene is a 2D semimetal that naturally hosts a 2DEG, whose carriers are Dirac-like massless fermions with a linear dispersion near the Dirac point. The atomically flat surface without dangling bonds of the 2D h-BN crystal provides graphene with the largest carrier mobility as compared to any other insulating substrate. Moreover, an h-BN interlayer thickness d of only 1 nm (a few monolayers) is sufficient to ensure the electrical isolation between two graphene layers, thus providing a strong interlayer interaction. In addition, the graphene/h-BN/graphene (G/h-BN/G) system allows to move from the Fermi-liquid regime to the Dirac point by tuning in-situ the Fermi energy (E F ), i.e. the carrier density (of either electrons (n) or holes (p)), through electrostatic gating. Thus, with n(p) in the 10 11 -10 13 cm À2 range, and with the carrier scattering length l given by l % nðpÞ À1=2 , the d=l ( 1 condition can be achieved easily 8 . Coulomb drag can be enhanced by plasmonic excitations supported by the 2DEGs in the active and passive layers. Moreover, it can also be assisted by phonon exchange through the dielectric interlayer. The contribution of both plasmons and phonons was first addressed in doped polar semiconductor heterostructures. The enhancement of the drag by plasmons was shown to become significant for T ! 0:2T F 22-24 , where T F is the Fermi temperature. In the case of phonons, at very low temperatures, these materials only showed a contribution from the acoustic phonons [25][26][27][28] . However, at higher temperatures, the electron-energy scales (Fermi and plasma energies) become comparable to the energy of the longitudinal optical (LO) phonon 29 , so that coupled plasmon-LO phonon modes were shown to contribute further to the Coulomb drag by virtual phonon exchange 30 .
Graphene supports plasmons in the THz to mid-infrared regime and their effect in Coulomb drag has been evaluated in double-layer systems 31 as well as in lateral waveguides 32 . In both cases, it has been shown that the drag resistivity can be enhanced by the plasmon contribution for T ! 0:15T F . Moreover, the weak coupling of the electrons to the graphene acoustic phonons 33 allows to study the interlayer phonon contribution to the Coulomb drag at high electron temperatures. However, none of these studies consider the full frequency-dependent permittivity of the dielectric in between the graphene layers or waveguides but just a static permittivity, thus lacking the contribution to the Coulomb drag from the phonon exchange mediated by the dielectric material. In G/h-BN/G systems, Amorim et al. 34 studied the contribution to the Coulomb drag of the virtual phonon exchange through the h-BN interlayer. Although the role of the h-BN anisotropy on the phonon-mediated mechanism was already pointed out, the hyperbolic nature of the phonons in this material was not described till a couple of years later [35][36][37] . Moreover, this study only took into account the contribution from the intraband region, neglecting that from the interband region where the plasmons lie. Therefore, there is a lack of a comprehensive study of the effect of plasmons and phonons altogether on the Coulomb drag in graphene.
The drag resistivity is typically calculated in the framework of the random phase approximation (RPA), which assumes that the electrons respond to the external field plus the time-dependent Hartree potential determined by the external and the polarization charges. The RPA neglects the short-range effects due to exchange and correlation (XC), which can modify the Coulomb interactions 38 . XC effects can be taken into account by employing local field factors (LFFs), which describe the repulsion hole around an electron due to the exchange repulsion and correlation. LFFs become more important for large wave vectors and small carrier densities, and they have been shown to have a significant impact on the plasmon dispersion and the drag resistivity of 2D electron bilayers 39 .
In this article, we study the Coulomb drag mechanism in G/h-BN/G structures along with the effects of the contribution of plasmons and phonons, and their coupling, by going beyond RPA and including the XC effects by means of LFFs. We show that the inclusion of the LFFs strongly modifies the drag resistivity becoming an essential ingredient for comparing theory and experiment. We also show the effect of the hyperbolic phonons of the h-BN interlayer on the drag resistivity, which has never been discussed so far according to our knowledge. In particular, we demonstrate that the hybridization of plasmons and phonons leads to a strong enhancement of the drag resistivity by the hyperbolic plasmon-phonon (HPP) modes in the intraband region, whereas the contribution from the surface plasmon-phonon (SPP) modes in the interband region reduces as compared to that from the unhybridized plasmon modes. We also demonstrate that the maximum value of the drag resistivity can be obtained within a physically realizable experimental regime.

Results
Plasmon-phonon hybridization and exchange-correlation effects. Here we consider a stacked structure of h-BN/G/h-BN/G/ h-BN. The top and bottom h-BN layers are thick enough to be considered semi-infinite, whereas the thin h-BN layer between the two graphene layers has a thickness d. Both graphene layers are doped, and carriers in one layer interact with those in the other layer by both Coulomb interaction and phonon exchange through the h-BN interlayer that leads to dressed interactions. The two graphene layers can be modelled using their polarizability, which can be calculated within the RPA framework. Here we consider the full finite temperature form of the graphene polarizability 40 . A transport-dominated regime, with τ >> τ ee , is assumed, where τ and τ ee are, respectively, the carrier transport and the electron-electron interaction times. h-BN is an anisotropic wide bandgap insulator whose in-and out-of-plane frequencydependent relative permittivities, ϵ x ðωÞ and ϵ z ðωÞ respectively, have opposite sign in the two reststrahlen bands, thus supporting a series of hyperbolic phonon modes that are confined within these bands (see Supplementary Note 1 for mathematical expressions and graphs of ϵ xðzÞ ðωÞ) 36,41 . To model the h-BN layer, we first define the expressions ϵ s ðωÞ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ϵ z ðωÞϵ x ðωÞ p and ϵ a ðωÞ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ϵ z ðωÞ=ϵ x ðωÞ p , that are plotted in Fig. 1a, b respectively.
The bare intra-(i ¼ j) and interlayer (i ≠ j) Coulomb interactions in the graphene layers (i, j = 1, 2) are then given by 34,42,43 where ϵ s ðωÞ acts as an effective permittivity and ϵ a ðωÞ is a factor responsible for the origin of the HPP in the two restrahlen bands of h-BN with opposite phase. These Coulomb interactions get modified by the XC effects, thus leading to the effective Coulomb interactions given by where G ij ðq; ωÞð¼ G ji ðq; ωÞÞ are the intra-(i ¼ j) and interlayer (i ≠ j) LFFs (see Supplementary Note 2 for more details and graphs) and V eff ij ðq; ωÞ ¼ V eff ji ðq; ωÞ. The total dielectric screening function is given by where χ ð0Þ l ðq; ω; TÞ is the temperature-dependent polarizability or irreducible density-density correlation function of layer l (l = 1, 2) 40 . Clearly, the XC effects are cancelled out and the RPA is recovered by setting G ij ðq; ωÞ ¼ 0.
The bare interactions will be screened by charge carriers in the two graphene layers. These screened interactions can be described by the coupled Dyson equations and written in matrix form as Uðq; ω; TÞ ¼ Vðq; ωÞ þ Vðq; ωÞ:χ ð0Þ ðq; ω; TÞ:Uðq; ω; TÞ; ð4Þ where Uðq; ω; TÞ is the 2 × 2 matrix of the screened interaction. From Eq. (4), one can find that the dressed interlayer interaction is given by Similarly, the reducible density-density correlation function Xðq; ω; TÞ can also be defined in terms of χ ð0Þ ðq; ω; TÞ and written in matrix form as Xðq; ω; TÞ ¼ Vðq; ωÞ þ Vðq; ωÞ:χ ð0Þ ðq; ω; TÞ:Xðq; ω; TÞ: ð6Þ In single graphene monolayer systems, one can define the electron loss function as the À=½1=ϵ T ðq; ω; TÞ to obtain the hybridized plasmon-phonon dispersion. However, in the doublelayer graphene system, this function changes sign for the acoustic branches. In oder to avoid this inconvenience, here we plot alternatively the imaginary part of the trace of the matrix Xðq; ω; TÞ, defined as the loss function which characterizes the spectral density of the hybridized quasiparticle excitations: Coulomb drag. For the double-layer graphene system, the drag resistivity is a 2 2 matrix given by In this equation, σ 11 ðσ 22 Þ is the temperature-dependent intralayer conductivity of the top (bottom) layer, given by the expression e 2 μ D τ=ðπ_ 2 Þ, where τ is the graphene carrier transport time, which is assumed to be energy independent, and μ D is the effective Drude weight, which depends on the carrier density and the temperature (see Supplementary Note 3 for more details on the graphene conductivity, chemical potential, and Drude weight), and σ D σ 12 ¼ σ 21 is the interlayer drag conductivity, which, up to second-order perturbation in the interlayer interaction, is given by where ω and q are the transferred energy and momentum from active to passive layer at temperature T with Γ l ðq; ω; TÞ being the rectification function or non-linear susceptibility of layer l 15 . The final expression for the drag resistivity can be then written as where Iðq; ω; TÞ is the drag intensity given by TOðLOÞ , are the in-(out-of-) plane transverse (longitudinal) optical phonon frequencies of h-BN. Both ϵ s ðωÞ and ϵ a ðωÞ vary rapidly around the two restrahlen bands. The y axis is the absolute magnitude of ϵ s ðωÞ and ϵ a ðωÞ, both of which are complex functions of ω, and the filling colour is the argument. The phase remains the same for ϵ s ðωÞ in the two restrahlen bands but it is of opposite sign for ϵ a ðωÞ, hence leading to a variation in the total phase of the exponential term of Eq. (1) within these two restrahlen bands, which gives rise to the hyperbolic modes confined within these bands.
Numerical results. Figure 2a, b shows the loss function, Lðq; ω; TÞ, calculated at T = 0 within the RPA only and with the inclusion of the XC effects, respectively. For simplicity, two identical (same carrier density) graphene layers separated by 6 nm have been considered to be placed in a homogeneous medium of dielectric constant 4.74 (i.e. the ϵ s ðω ¼ 0Þ of h-BN). The Lðq; ω; TÞ, which is related to the dynamical structure factor Sðq; ω; TÞ that accounts for the spectral strength of the various elementary excitations, characterizes the collective modes (plasmons) of the coupled 2DEGs of the two graphene layers. There are two such modes, one where the electron densities in the two layers oscillate in phase, called as the optical mode (O, ω $ ffiffi ffi q p ), and one where the oscillations are out of phase, called as the acoustic mode (A, ω $ q). Since both modes lie within the Pauliblocked interband region (v F q < ω < v F ð2k F À qÞ), i.e. outside of the particle-hole continuum (PHC), they cannot decay via intraband excitations. However, both modes get damped as soon as they enter into the region ω < v F q and ω > v F ð2k F À qÞ, as they decay via interband PHC transitions (Landau damping) that are now allowed by the energy-momentum conservation. It has to be also noted that the acoustic plasmon becomes degenerate with the PHC as d ! 0. The comparison of the colour scale bars of Fig. 2a,  b indicates that the maximum value of Lðq; ω; 0Þ in the XC case is almost double that in the RPA case. Moreover, the two plasmon branches shift to higher q values for a fixed ω and to lower ω values for a fixed q, as shown, respectively, in Fig. 2c, d. In addition, the inclusion of the XC reduces the interlayer Coulomb interactions compared to those provided by the RPA alone, as shown in Fig. 2e, f for a fixed value of both ω and q, respectively. Thus, the XC increases the ϵ T and decreases the V eff 12 , which in turn decreases the dressed interlayer Coulomb interaction, U 12 ðq; ω; TÞ, given by Eq. (5). As a result, the drag resistivity, which is directly proportional to j U 12 ðq; ω; TÞ ð j 2 , decreases with the inclusion of XC (see Supplementary Note 2 for more details and plots). Figure 3a shows the loss function, Lðq; ω; TÞ, whereas Fig. 3b, c shows the Coulomb drag intensity, Iðq; ω; TÞ, without and with phonons, respectively, all calculated for the h-BN/G/h-BN/G/h-BN system for a temperature of 50 K. Figure 3d-f plots the same magnitudes, respectively, as Fig. 3a-c but for a temperature of 300 K. For simplicity, both graphene layers are considered to have identical carrier density (with E F1 ¼ E F2 ¼ 0:064 eV) and the h-BN interlayer thickness is d ¼ 6 nm. The effect of h-BN phonons and their hybridization with the plasmons in graphene can be eliminated by setting ϵ s ðω ¼ 0Þ and ϵ a ðω ¼ 0Þ in the expressions of V eff 11ð22Þ and V eff 12ð21Þ only (see Supplementary Note 1 for more details and graphs). In this way, the direct comparison of the results with and without phonons provides the net contribution of the phonon exchange to the drag intensity. In the Lðq; ω; TÞ plot, both plasmon modes, optical and acoustic, hybridize with the phonons of h-BN. In the case of the higher temperature value (300 K, Fig. 3d), both optical and acoustic modes start to get damped. The reason for this can be found in the imaginary part of the graphene temperature-dependent polarizability =½χ ð0Þ l ðq; ω; TÞ, shown in Fig. 4, that represents the PHC. In this figure, the hybridized SPP modes are indicated by the superimposed black lines. At low temperature, the plasmon  7), calculated for the case of a RPA and b XC for the double-layer graphene system with background dielectric constant ϵ ¼ 4:74. The solid lines separate the regions with Landau damping due to intraband (ω < v F q) and interband (ω > v F q, v F ð2k F À qÞ) particle-hole continuum and the Pauli-blocked interband region (v F q < ω < v F ð2k F À qÞ). Cross section of Lðq; ω; 0Þ for c ω=E F ¼ 0:9 and d q=k modes, which lie in the region v F q < ω < v F ð2k F À qÞ, are undamped as they do not overlap with the PHC, as shown in Fig. 4a for T ¼ 50 K. Nevertheless, at high temperature, these plasmon modes overlap with the PHC, as shown in Fig. 4b for T ¼ 300 K, thus reducing their lifetime as they now can decay via particle-hole excitations. The drag intensity, which is directly proportional to the drag resistivity, ρ D , allows us to discriminate the different quasiparticle excitations contributing to the Coulomb drag. Figure 3b, c indicates that, at low temperature (50 K), the major contribution to the drag resistivity comes from the intraband region, ω < v F q, and it is identical for both cases of with and without phonons.
This has a twofold explanation. On the one hand, the contribution to the drag resistivity is restricted to low frequency values as the temperature T in the sinh term in the denominator of Eq. (11) sets an effective upper cut-off in frequency. Thus, at low temperature, the HPP modes, appearing at higher frequencies (the two reststrahlen bands lie at around 100 and 180 meV), do not contribute to the Coulomb drag. On the other hand, the drag intensity is proportional to the non-linear susceptibility, which vanishes at low temperature in the region where the SPP modes lie, v F q < ω < v F ð2k F À qÞ, as shown in Fig. 5a. In this figure, the hybridized SPP modes are indicated by the superimposed black lines. Thus, at low temperature, these modes do not contribute either to the Coulomb drag.
The drag intensity changes dramatically at high temperature, as shown, for example, in Fig. 3e, f for T ¼ 300 K. In the case when phonons of h-BN are not considered, Fig. 3e, the drag intensity is markedly enhanced due to the contribution from the optical and acoustic plasmon modes in the interband region. When the effect of the h-BN phonons is taken into account, as shown in Fig. 3f, the strong enhancement of the drag intensity is caused by the contribution of both hybridized surface and HPP modes. These contributions clearly mimic the dispersion shown by the loss function in Fig. 3d, where the SPP modes are contained in the interband region while the HPP modes propagate confined within the two reststrahlen bands of h-BN even into the intraband region.
Apart from the general higher effective upper cut-off in frequency, the specific reason for the enhancement of the Coulomb drag intensity at higher temperatures lies in the fact that the scattering matrix element given by the screened Coulomb interaction U 12 ðq; ω; TÞ is inversely proportional to ϵ T ðq; ω; TÞ (see Eq. (5)). The latter becomes very small (vanishing at T ¼ 0) for certain ωðqÞ corresponding to the plasmons of the system, here optical and acoustic modes. This makes U 12 ðq; ω; TÞ very with Landau damping due to intraband (ω < v F q) and interband (ω > v F q, v F ð2k F À qÞ) particle-hole continuum and the Pauli-blocked interband region (v F q < ω < v F ð2k F À qÞ), whereas the dashed lines correspond to the bulk phonons of h-BN that define the two reststrahlen bands. The coupling of the graphene plasmons (optical and acoustic) and the h-BN phonons leads to hybridized surface and hyperbolic plasmon-phonon modes (SPP and HPP, respectively); the superscript in HPP denotes the upper (U) and lower (L) reststrahlen bands. Density plot of the drag intensity, Iðq; ω; TÞ, calculated at T ¼ 50 K without phonons (b), T ¼ 50 K with phonons (c), T ¼ 300 K without phonons (e), and T ¼ 300 K with phonons (f). The temperature and the screening condition strongly influence the Coulomb drag. Parameters are E F1 ¼ E F2 ¼ 0:064 eV and d ¼ 6 nm. The scale bars in a and b also apply, respectively, to d and c, e and f. Fig. 4 Polarizability of graphene. Density plot of the imaginary part of the graphene polarizability, =½χ ð0Þ l ðq; ω; TÞ, for a temperature T ¼ 50 K (a) and T ¼ 300 K (b) with E F ¼ 0:064 eV. The dashed lines separate the regions with Landau damping due to intraband (ω < v F q) and interband (ω > v F q, v F ð2k F À qÞ) particle-hole continuum and the Pauli-blocked interband region (v F q < ω < v F ð2k F À qÞ). The superimposed black curves are the hybridized surface plasmon-phonon (SPP) modes of the h-BN/graphene/ h-BN/graphene/h-BN system. large at the ωðqÞ of these plasmons, which greatly increase the total electron-electron scattering in the system. At low temperature, it has already been explained above that the contribution of the plasmons to the Coulomb drag is suppressed because the nonlinear susceptibility, Γ l ðq; ω; TÞ, vanishes in the region v F q < ω < v F ð2k F À qÞ where the plasmons lie, Fig. 5a. However, at high temperature, Γ l ðq; ω; TÞ acquires a non-vanishing finite value in this region that overlaps with the plasmon-phonon modes, as shown in Fig. 5b. In this figure, the hybridized SPP modes are indicated by the superimposed black lines. Therefore, the plasmon-phonon modes start to contribute to the Coulomb drag intensity, as pointed out in Fig. 3e, f. It should also be noted that the plasmon enhancement of the Coulomb drag is primarily due to acoustic mode, since it is lower in energy than the optical mode and runs closer to the PHC, thus facilitating its excitation at lower temperatures than the optical mode. Nonetheless, the optical mode also contributes to the drag, being its contribution larger the higher the temperature.
The drag resistivity, ρ D , has been numerically calculated as a function of the three parameters that govern it, namely the temperature T, the carrier density (defined here in terms of E F ), and the h-BN interlayer thickness d. In order to isolate the effects of the different quasiparticle excitations, ρ D has been calculated with and without phonons, for the whole ω-q space as well as for the inter-and intraband regions separately. Figure 6a-c presents ρ D as a function of T calculated for E F1 ¼ E F2 ¼ 0:064 eV and d ¼ 4, 6, and 9 nm, respectively, showing a good agreement with the experimental data of drag resistivity 8 . In all cases, ρ D is initially enhanced with increasing temperature, peaking at around 300 K. Beyond that threshold, ρ D starts to decrease. In the G/h-BN/G system considered here, this peak gets considerably more intense when the plasmon-phonon hybridization is considered (solid lines), as compared to the case where only plasmons are taken into account (dashed lines). Thus, the full consideration of the hybridized SPP and HPP modes leads to a 10%, 16%, and 25% percent higher values of ρ D for d ¼ 4, 6, and 9 nm, respectively. Although the overall magnitude of ρ D decreases with increasing d, the larger difference between the cases with and without phonons for the thicker h-BN interlayer points out the strong contribution arising from the HPP modes confined in the h-BN reststrahlen bands, whose number increases with increasing d. The final reduction of ρ D beyond the peak temperature is a consequence of the merging of the PHC with plasmon modes as shown in Fig. 4b, as explained earlier. It is important to remark that these results show that ρ D can be enhanced by various collective excitations even at room temperature and beyond.
In Fig. 7a, b, the ρ D of Fig. 6b is calculated separately for the inter-and intraband regions, respectively. In the interband region, Fig. 7a, ρ D slightly decreases when considering the h-BN phonons because the two plasmon branches split into multiple hybridized SPP modes which have finite lifetime. This splitting occurs around the two reststrahlen bands of h-BN. Since only HPP modes are allowed to propagate in the restrahlen bands and they primarily lie in the intraband region, the spectral weight of hybridized SPP modes in the interband region is smaller than that of their counterparts (plasmon modes) when h-BN phonons are not considered. Therefore, there is a net slight reduction in the interband contribution to ρ D with inclusion of h-BN phonons. In the intraband region, Fig. 7b, the difference in the ρ D when h-BN phonons are considered is the result of the enhancement of ρ D by the HPP modes, which become more prominent even at room temperature and beyond. Figure 7c shows the qualitative evolution of ρ D with respect to the temperature (normalized by T F ) due to the contributions of the SPP and HPP modes, which are shown separately as shaded colour regions. For T ≲ 0:01T F , ρ D $ T 3 as shown in the inset. Above that range and for T ≲ 0:07T F ðT PHC Þ, ρ D $ T 2:2 is dominated by the PHC. When the temperature is increased further, the SPP modes start to contribute with ρ D $ T 1:75 till T SPP ($ 0:3T F ). In the absence of any other excitations, ρ D is expected to decrease with increasing temperature. However, Fig. 7c shows that ρ D keeps increasing Fig. 5 Non-linear susceptibility of graphene. Density plot of the non-linear susceptibility of graphene, Γ l ðq; ω; TÞ, for a temperature T ¼ 50 K (a) and T ¼ 300 K (b) with E F ¼ 0:064 eV. The dashed lines separate the regions with Landau damping due to intraband (ω < v F q) and interband (ω > v F q, v F ð2k F À qÞ) particle-hole continuum and the Pauli-blocked interband region (v F q < ω < v F ð2k F À qÞ). The superimposed black curves are the hybridized surface plasmon-phonon (SPP) modes of the h-BN/graphene/ h-BN/graphene/h-BN system. shifting the maximum value of ρ D towards a higher temperature, T HPP ($ 0:45T F ) and provinding a ρ D $ T 0:35 dependence due to the contribution from the HPP modes. Beyond T HPP , the system enters finally into the Boltzmann regime where ρ D decreases as $1=T 2 .
The calculated dependence of ρ D on E F and d is shown in Fig. 8a-c along with the experimental data from ref. 8 . The comparison shows an excellent agreement with the theory. The ρ D values decrease with increasing E F , as shown in Fig. 8a, b for T ¼ 240 K and 300 K, respectively. The functional form of ρ D varies as ρ D $ 1=E β F (or $1=n β=2 ), with a β value of 2.3, 2.7, and 2.7 for d ¼ 4, 6, and 9 nm respectively, which agrees well with the reported experimental data 8 . Similarly, ρ D decreases with increasing d as the two graphene layers start to decouple, as shown in Fig. 8c for both T ¼ 240 K and 300 K. In this case, the functional form of ρ D has been assessed separately for the k F d < 1 and k F d > 1 regimes, with ρ D varying as 1=d γ with γ ¼ 1:5 and 1.75, respectively. It should be noted that, close to the Dirac point, the effects of the charge density fluctuations, known as electronhole puddles, in the two graphene layers need to be taken into account 44 . In the calculations presented here, far from the Dirac point, these effects have been neglected.

Discussion
In this report, we have studied the Coulomb drag mechanism in an h-BN/G/h-BN/G/h-BN system and have evaluated the contribution of plasmons and phonons, as well as their hybrization, thereto. We have also included the effects of XC in the standard RPA, and have found that these effects further modify the plasmon modes of the system and thus drag resistivity. XC allows us to replicate the experimental data to a very good accuracy. In order to discriminate the different quasiparticles and their effects, the contributions from the inter-and intraband regions with and without h-BN phonons have been evaluated. We have showed that apart from the drag enhancement led by the graphene plasmons, their hybridization with the hyperbolic phonons of h-BN, which extend into the intraband region, can further enhance the drag resistivity even at room temperature and beyond. This contrasts with the surface-like hybridized plasmon-phonon modes in the interband region, whose contribution to the drag resistivity reduces, as compared to the unhybridized plasmon modes, due to their splitting imposed by the h-BN restrahlen bands. Moreover, the calculated maximum of the Coulomb drag resistivity takes place at room temperature for physically realizable experimental conditions. This result is relevant, as Coulomb drag experiments in G/h-BN/G systems have not been yet explored for temperatures larger than 240 K 8 .
It is also worth mentioning that the XC effects introduced in this study via LFFs obtained from the quantum Monte Carlo (QMC) method are applicable for the Coulomb drag in singlelayer graphene, where the interaction parameter (r s ) has a small value. However, under different experimental conditions, such as Fig. 8 Fermi energy-and interlayer thickness-dependence of drag resistivity. Drag resistivity, ρ D , calculated as a function of the Fermi energy E F (E F ¼ E F1 ¼ E F2 ) for T ¼ 240 K (a) and T ¼ 300 K (b) calculated for an h-BN thickness d ¼ 4, 6, and 9 nm. c ρ D as a function of the h-BN interlayer thickness d (with E F ¼ E F1 ¼ E F2 ¼ 0:064 eV and for T ¼ 240 K and T ¼ 300 K). The functional forms of ρ D on E F and d are also indicated. The comparison with the data from the experiment reported in ref. 8 (circles in a and c) shows a very good agreement with the theory. Fig. 7 Contributing mechanisms to temperature-dependence of drag resistivity. Drag resistivity, ρ D , same as in Fig. 6b but calculated separately for the a interband and b intraband regions as a function of temperature with (solid line) and without phonons (dashed line). The inclusion of the phonons of h-BN enhances ρ D in the intraband region while reduces it in the interband region. c Dependence of ρ D in Fig. 6b on T=T F with the contributions from the particle-hole continuum (PHC) and the different quasiparticles, namely hybridized surface and hyperbolic plasmon-phonon modes (SPP and HPP, respectively). The vertical dashed lines indicate their characteristic temperatures, whereas the temperature functional form of ρ D for the different regions is indicated along its line shape. The inset shows a magnified view of the low T=T F values.
for example single-layer graphene in a large magnetic field, or in other systems, such as bilayer graphene (BLG), r s can take large values and thus it might be necessary to use other forms of the LFFs (beyond the QMC method) and/or express r s in terms of ϵ T ðq; ω; TÞ instead of ϵ s ðωÞ to fit the data (see Supplementary Note 2 for details).
In addition, the described mechanisms of the enhancement of the Coulomb drag in the double-layer graphene system G/h-BN/G, led by the graphene plasmons hybridized with the hyperbolic phonons of h-BN, are also foreseen to have an important role in the more complex double BLG systems such as BLG/WSe 2 /BLG 45 and BLG/h-BN/BLG 46 , where Coulomb drag has already been measured at room temperature 47,48 . In this case, it is necessary though to take into account the different plasmon dispersion of the BLG 49 . Furthermore, the impact of these results may also extend to the analysis of the Coulomb drag in other 2D material systems with a h-BN interlayer, such as phosphorene/h-BN/phosphorene 50 , or more generally to any future double 2D layer systems including a naturally hyperbolic 2D interlayer material, such as, for example, black phosphorous 51 or transition metal dichalcogenides 52 . This can open new ways for the exploration of novel many-body physics phenomena. For example, superposing graphene on black phosphorous has been shown to lead to strained superlattices, which in turn create large pseudomagnetic fields 53 that are shown to have substantial effects on the Coulomb drag 8 .

Data availability
Data that support the findings of this study are available from the authors on reasonable request.