Electronic cooling via interlayer Coulomb coupling in multilayer epitaxial graphene

In van der Waals bonded or rotationally disordered multilayer stacks of two-dimensional (2D) materials, the electronic states remain tightly confined within individual 2D layers. As a result, electron–phonon interactions occur primarily within layers and interlayer electrical conductivities are low. In addition, strong covalent in-plane intralayer bonding combined with weak van der Waals interlayer bonding results in weak phonon-mediated thermal coupling between the layers. We demonstrate here, however, that Coulomb interactions between electrons in different layers of multilayer epitaxial graphene provide an important mechanism for interlayer thermal transport, even though all electronic states are strongly confined within individual 2D layers. This effect is manifested in the relaxation dynamics of hot carriers in ultrafast time-resolved terahertz spectroscopy. We develop a theory of interlayer Coulomb coupling containing no free parameters that accounts for the experimentally observed trends in hot-carrier dynamics as temperature and the number of layers is varied.

Temperature dynamics when interlayer energy transfer between LD layers (n LD = 10 10 cm −2 ) is ignored. The LD layer temperature T LD (t) resulting from cooling via interlayer Coulombic energy transfer to the HD layers near the substrate (n HD 10 12 cm −2 ) at a constant lattice temperature T L . The distance between the particular LD layer and the HD layers is varied, d HD,LD = 60a, 50a, 40a, 30a, 20a, 10a, a (top to bottom) where a = 3.4 Angstroms. Subfigure (a) shows results for lattice temperature T L = 10 K and Subfigure (b) for T L = 50 K.   Figure 6. Disorder-assisted electron-phonon (supercollision) cooling in LD graphene. Electron temperature dynamics T (t) − T L predicted by the disorder-assisted electron-phonon cooling mechanism for LD graphene with E F = 10 meV at T L = 10 K for variable disorder mean free path l.    Figure 9. Disorder-assisted electron-phonon (supercollision) cooling in HD graphene. Electronic cooling time τ HD in the low electron temperature limit predicted by the disorder-assisted electron-phonon cooling mechanism for HD graphene with E F = 100 meV as a function of T L for variable disorder mean free path l.  Figure 10. Disorder-assisted electron-phonon (supercollision) cooling in LD graphene. Electronic cooling time τ LD in the low electron temperature limit predicted by the disorder-assisted electron-phonon cooling mechanism for LD graphene with E F = 10 meV as a function of T L for variable disorder mean free path l.

SUPPLEMENTARY NOTES
Supplementary Note 1. Normalized differential THz transmission spectra ∆t(ω)/t(ω) Here, we present additional experimental data to further demonstrate that the normalized differential THz transmission spectra are remarkably dispersionless in the detectable frequency range under all experimental conditions. Supplementary Figure 1 and Supplementary Figure 2 show the differential THz transmission spectra normalized to the THz transmission without photoexcitation, ∆t(ω)/t(ω), for variable substrate temperature and variable pump fluence, respectively, for a MEG sample with ∼ 63 layers. The frequencyindependent dynamic THz response of the MEG samples justifies recording the normalized differential THz transmission only at the peak of the THz probe pulse, ∆t/t, as a function of pump-probe delay to map out the relaxation dynamics of the photoexcited carriers. The slight fluctuations in a few of the data scans for frequencies 1.7 THz are due to water vapor absorption arising from very slight fluctuations in the humidity level between the sample and the reference scans. We note that spectra obtained from Fourier-transformed timedomain measurements are much more susceptible to noise than direct spectrally-resolved measurements, because (long-term) fluctuations in the time domain are transferred into uncertainties in the frequency domain during the Fourier transformation process. On the other hand, direct time-resolved THz spectroscopy measurements of the relaxation dynamics are less prone to noise from fluctuations, and a very high signal-to-noise ratio (SNR) can be achieved through sufficient integration.

Supplementary Note 2. Interlayer energy transfer linearized in ∆T
Here, we derive an approximation for the interlayer Coulombic energy transfer rate between a lightly doped (LD) layer and a highly doped (HD) layer in MEG systems. We neglect electron tunneling between layers. Energy transfer between layers can nevertheless occur when an electron in one layer scatters off an electron in a remote layer. To first order in the temperature difference δT = T LD − T HD , we find that the energy transfer rate between thermal electron distributions in two Coulomb-coupled layers is given by: (1) Equation 1 can be derived by summing over all interlayer electron collision processes and combining a Fermi golden-rule expression for the transition rates with a random phase approximation expression for the electron-electron scattering amplitudes. We focus on the low temperature limit (i.e. T ≪ T F,LD , T F,HD ), where several physical approximations can be made. In this case, it is clear that the ω → 0 limit dominates the integrand of Equation 1. We note that in this low frequency limit, is the density of states at the Fermi energy. This linear dependence on frequency is also found in parabolic band 2DEG's, and describes how the particle-hole excitation spectrum vanishes with decreasing excitation energy ω [1]. Additionally, in the degenerate regime interlayer particle scattering allows a maximum change in electron wavevector of ∆q max = 2k F , a fact reflected in: The interlayer Coulomb interaction, proportional to e −qd , naturally places the limit q 1/d, where ν j is the density of states in the j'th HD layer and κ is the background dielectric function of a thin film on SiC (κ ≈ (10 + 1)/2 = 5.5). Because of the large difference in carrier density between layers near and far from the substrate, allowed q values are much smaller than q MEG TF and the screened interlayer interaction reduces to (2πe 2 )/(κq MEG TF ). The remaining integrals in Equation 1 contain a logarithmic divergence at zero wavevector. This is removed using a cutoff of ω/v F which reflects the zero value of Im[χ(q, ω)] above the intraband particle-hole continuum.
Finally, we find that the linearized interlayer energy transfer rate per area between a pair of graphene layers in a MEG system (where one is HD and the other LD) is: where The density of states in the HD layer (LD layer) is denoted by ν HD (ν LD ). The sum over index j runs over all HD layers in the MEG system, accounting for the approximation that the dominant screening effect originates from the majority fraction of carriers in the HD layers (see the main text). T F,LD = E F,LD /k B is the Fermi temperature of the LD layer. Equation 4 becomes exact in the degenerate limit T ≪ T F,LD , T F,HD .

Supplementary Note 3. Interlayer energy transfer to leading order in T
Here, we derive in detail the leading order in temperature formula for the interlayer Coulombic energy transfer rate between one lightly doped (LD) layer of graphene and one highly doped (HD) layer of graphene. This follows similar steps as in the previous section, but without the assumption of infinitesimal temperature separation between the LD and HD layers. As described in the main text, the interlayer energy transfer rate per area between an LD and an HD layer is: where n B (T ) is a Bose distribution function and v sc LD,HD = v q /ǫ RPA (q, ω, T LD , T HD ) is the screened interaction between electrons in opposing layers within the random phase approximation (RPA). Given that the separation of the two layers is a distance d, the RPA dielectric function is given by: where v q = 2πe 2 /κq and χ i (q, ω, T ) is the temperature dependent non-interacting densityresponse function of the i'th layer of graphene. We next assume the temperature of the electrons in the HD layer is pinned to the lattice temperature, and approximate this as zero relative to the high temperature in the LD layer, i.e. T HD = T L → 0. Relabeling T LD as T we have: In the degenerate limit, k B T << E F,LD , the Bose distribution function limits the important frequencies to approximately ω ≤ E F,LD / , and reveals that it is the leading order in ω which will yield the leading order in T . With this motivation we make use of the low frequency limit, Im[χ i (q, ω → 0)] = −ν i ω/(v F q), where ν i = 2E F,i /(π 2 v 2 F ) is the density of states at the Fermi energy. Similarly, the dielectric function can be reduced to: where the Thomas-Fermi wavevector of the i'th layer is defined as q TF i = qv q ν i . Making use of the fact that in our MEG samples k F,LD d << 1 and q TF HD /k F,LD >> 1, we simplify the interlayer energy transfer rate per area to: where we have introduced the dimensionless wavevector, Q = q/k F,LD , the dimensionless frequency, Ω = ω/E F,LD , and the dimensionless temperature, t = k B T /E F,LD . The wavevector integration here diverges logarithmically. To remedy this, we identify that Im[χ LD (q, ω, 0)] vanishes for ω > v F q, which cuts off the wavevector integration for Q < Ω. Carrying out the remaining integrals, we obtain the leading order in temperature energy loss rate per area of the LD electrons:

Supplementary Note 4. Interlayer energy transfer with no LD-LD layer coupling
Here, we illustrate the distance dependence of the interlayer Coulombic energy transfer rate by calculating the temperature dynamics of an individual LD layer coupled to the HD layers in MEG systems, when we neglect energy transfer between pairs of LD layers.
Although Coulomb coupling between LD layers is very strong [2], by temporarily forcing Q el LD,LD ′ → 0 we can gain insight into the distance dependence of the interlayer energy transfer rate. The temperature dynamics of each individual LD layer coupled to the HD layers are then independent and obey: If we also approximate the heat capacity in the LD layers by the neutral graphene formula Here, we present some simple calculations in support of our assumption that acoustic phonon cooling is capable of pinning the electronic temperature in the HD layers at the lattice temperature while these layers act as a heat sink for energy dissipation from the remaining hot LD layers. Heuristically, we also note that previous ultrafast optical spectroscopy experiments [3] have observed thermal relaxation times in the HD layers of MEG of ∼ 10 ps, much faster than the relaxation times in the LD layers of ∼ 100 − 500 ps that we report here.
Acoustic phonon cooling has previously been investigated in the context of hot-carrier cooling in disorder-free monolayer graphene [4,5]. In these single layer systems, as a result of the relatively large optical phonon energy in graphene ( ω op ≈ 200 meV [6]), acoustic phonons serve as the primary intrinsic cooling mechanism over a large temperature window extending up towards T ∼ 250 K [4]. To estimate the ability of acoustic phonon cooling to take away the electronic energy that is Coulomb-transfered from LD to HD layers, we use Equation 14 in Bistritzer and MacDonald [4] to calculate the total energy transfer rate to the the lattice. Using the carrier density profile of the four most highly doped layers measured in experiment (E F,1−4 = 360, 218, 140, and 93 meV), we find that Q ph HD = 9.09(T HD − T L ) W cm −2 can be transferred to the lattice via carrier-phonon scattering in the HD layers.
We can calculate the ratio of Q ph to the interlayer energy transfer rate Q el from the N LD layers to these four HD layers using Equation 13 in the main text. Supplementary Figure 4 suggests that for both N = 30 and N = 60 layer MEG systems, the acoustic phonon cooling power is sufficient to keep the HD layers pinned to the lattice temperature while absorbing energy from the distant LD layers.

Supplementary Note 6. Disorder-assisted electron-phonon (supercollision) cooling
Here, we present an application of the recently proposed disorder-assisted electron-phonon (supercollision) cooling mechanism [7] to MEG. We investigate the qualitative and quantitative differences between disorder-assisted electron-phonon cooling in HD and LD graphene and we clearly demonstrate that this cooling mechanism cannot alone explain electronic cooling in high quality MEG.
The electron temperature dynamics of a single graphene layer in the framework of the disorder-assisted electron-phonon cooling mechanism is given by: where Q sc = ∂ t E is the electronic cooling rate due to supercollisions, C = ∂ T E is the electronic heat capacity and E is the electronic energy density. The electronic cooling rate due to supercollisions in the degenerate limit when k B T ≪ E F is given by [7]: with a rate coefficient A = 9.62g 2 ν 2 (E F )k 3 B /( k F l), where g = D/ 2ρv 2 s is the electronphonon coupling constant and ν(E F ) = E F /(2π 2 v 2 F ) is the density of states at the Fermi level per one spin and one valley flavor. The rest of the parameters are the Fermi velocity v F , the sound velocity v s , the deformation potential D, the mass density ρ and the disorder mean free path l. The electronic cooling rate in the non-degenerate limit when k B T ≫ E F is given by [8]: with a rate coefficient B = ((4k 2 B )/(E 2 F ))A. Because both the supercollision cooling rate and the electronic heat capacity have different functional dependence on the electron temperature at high and at low doping densities, we need to consider the two cases separately.
For HD graphene, the heat capacity is given by C = αT = ((2πE F k 2 B )/(3 2 v 2 F ))T . By substituting this expression and Equation 13 in Equation 12, we obtain [9]: The electronic cooling timescale is governed by the rate coefficient A/α ∝ E F /k F l ∝ 1/l. In the high and low electron temperature limits, Equation 15 reduces to: with solutions given by: For LD graphene, the heat capacity is approximated by C = βT 2 = ((18ζ(3)k 3 B )/(π 2 v 2 F ))T 2 . By substituting this expression and Equation 14 in Equation 12, we obtain: The electronic cooling timescale is governed by the rate coefficient In the high and low electron temperature limits, Equation 20 reduces to: with solutions given by: We observe that the different functional form of the cooling rate and the heat capacity at high and at low doping densities results in qualitatively different electron temperature relaxation dynamics. In the low electron temperature limit, in particular, the temperature dynamics follows an exponential form with a lattice temperature dependent characteristic time τ HD = ((3A/α)T L ) −1 ∝ l/T L for HD graphene, and with a lattice temperature dependent characteristic time τ LD = ((5B/β)T 2 L ) −1 ∝ E F l/T 2 L for LD graphene. To obtain quantitative estimates of the full electron temperature dynamics predicted by the disorder-assisted electron-phonon cooling mechanism, we solve Equation 15 and Equation 20 numerically. In all calculations, we use v F = 1 × 10 6 m s −1 , v s = 2.1 × 10 4 m s −1 , ρ = 7.6 × 10 −7 kg m −2 and D = 20 meV. Supplementary Figure 5 and Supplementary Figure 6 show the calculated electron temperature dynamics for HD graphene with E F = 100 meV and for LD graphene with E F = 10 meV, respectively, at T L = 10 K for variable disorder mean free path l. The large uncertainty in the value of the disorder length scale leads to a very large spread in the calculated electron temperature dynamics.
The experimental measurement of the precise value of the disorder mean free path between supercollisions is a challenging task and reliable estimates are lacking in the literature.
One reliable way to characterize the degree of disorder which can potentially be related to supercollision cooling is directly from the width of the Dirac cone near the Dirac point in the graphene band structure as measured in high-resolution angle-resolved photoemission spectroscopy (ARPES). We use such recent ARPES measurements to estimate the disorder mean free path value appropriate for our MEG samples [10]. In Ref. 10 ( Figure S1), correlation lengths of ∼ 1 − 3 nm in exfoliated graphene and correlation lengths exceeding ∼ 50 nm (limited only by the instrument resolution, but are expected to be even longer) in C-face MEG have been reported. Supplementary Figure 7 shows the calculated electron temperature dynamics for HD graphene with E F = 100 meV and disorder mean free path l = 5 nm for variable T L . We note that these calculations are consistent with recent experimental results based on photocurrent measurements on HD chemical-vapor-deposited (CVD) graphene [9]. Supplementary Figure 8 shows the calculated electron temperature dynamics for LD graphene with E F = 10 meV and disorder mean free path l = 50 nm for variable T L .
We also present the electronic cooling times in the low electron temperature limit predicted by the disorder-assisted electron-phonon cooling mechanism that are more straightforward to compare to the experiment. Supplementary Figure 9 and Supplementary Figure 10 show the calculated electronic cooling times τ HD for HD graphene with E F = 100 meV and τ LD for LD graphene with E F = 10 meV, respectively, as a function of T L for variable disorder mean free path l. We observe that the predictions of the supercollision cooling model applied to the HD layers of MEG are inconsistent both in magnitude and lattice temperature dependence with the ultrafast time-resolved THz spectroscopy experiments. We also observe that the model predictions for the LD layers of MEG can roughly capture the order of magnitude and the lattice temperature dependence for some disorder mean free path, but not the layer number dependence. Because the quality of our MEG samples is expected to be even higher than we have conservatively assumed in these calculations, we conclude that the disorder-assisted electron-phonon (supercollision) cooling can provide a parallel cooling channel, but it is not dominant in MEG.