# Phonon-limited resistivity of multilayer graphene systems

Seth M. Davis<EMAIL_ADDRESS>Condensed Matter Theory Center and Joint
Quantum Institute, Department of Physics, University of Maryland, College
Park, MD 20742, USA  Yang-Zhi Chou  Condensed Matter Theory Center and Joint
Quantum Institute, Department of Physics, University of Maryland, College
Park, MD 20742, USA  Fengcheng Wu  School of Physics and Technology, Wuhan
University, Wuhan 430072, China  Wuhan Institute of Quantum Technology, Wuhan
430206  Sankar Das Sarma  Condensed Matter Theory Center and Joint Quantum
Institute, Department of Physics, University of Maryland, College Park, MD
20742, USA

###### Abstract

We calculate the theoretical contribution to the doping and temperature ($T$)
dependence of electrical resistivity due to scattering by acoustic phonons in
Bernal bilayer graphene (BBG) and rhombohedral trilayer graphene (RTG). We
focus on the role of nontrivial geometric features of the detailed,
anisotropic $\bm{\mathrm{k}}\cdot\bm{\mathrm{p}}$ band structures of these
systems - e.g. Van Hove singularities, Lifshitz transitions, Fermi surface
anisotropy, and band curvature near the gap - whose effects on transport have
not yet been systematically studied. We find that these geometric features
strongly influence the temperature and doping dependencies of the resistivity.
In particular, the band geometry leads to a nonlinear $T$-dependence in the
high-$T$ equipartition regime, complicating the usual $T^{4}$ to $T$ Bloch-
Grüneisen crossover. Our focus on BBG and RTG is motivated by recent
experiments in these systems that have discovered several exotic low-$T$
superconductivity proximate to complicated hierarchies of isospin-polarized
phases. These interaction-driven phases are intimately related to the
geometric features of the band structures, highlighting the importance of
understanding the influence of band geometry on transport. While resolving the
effects of the anisotropic band geometry on the scattering times requires
nontrivial numerical solution, our approach is rooted in intuitive Boltzmann
theory. We compare our results with recent experiment and discuss how our
predictions can be used to elucidate the relative importance of various
scattering mechanisms in these systems.

##  I Introduction

Rapid progress in the ability to produce clean, stable, 2D layered van der
Walls heterostructures made up of graphene and/or transition metal
dichalcogenides (TMDs) has opened a new subfield of condensed matter physics
Geim and Grigorieva (2013); Novoselov et al. (2006); Bistritzer and MacDonald
(2011); Morell et al. (2010); Li et al. (2019); Kim et al. (2017); Cao et al.
(2018a, b, 2020a, 2021); Yankowitz et al. (2019); Kerelsky et al. (2019); Lu
et al. (2019); Stepanov et al. (2020); Sharpe et al. (2019); Chen et al.
(2020); Rozen et al. (2021); Zhou et al. (2022, 2021a, 2021b); Serlin et al.
(2020); Wu et al. (2018, 2019a); Tschirhart et al. (2022); Polshyn et al.
(2020); Jaoui et al. (2021); Polshyn et al. (2019a); Cao et al. (2020b); Sarma
and Wu (2022); Zhang et al. (2022); Polski et al. (2022); Arora et al. (2020);
Xie and MacDonald (2020); Andrei and MacDonald (2020); Li et al. (2021);
Ghiotto et al. (2021); Pan et al. (2020); Pan and Sarma (2021); Morales-Durán
et al. (2021); Ahn and Sarma (2022); Kerelsky et al. (2021); Khalaf et al.
(2019). The sensitivity of the band structures of these systems to external
control parameters, especially twist angle and displacement field, gives an
unprecedented experimental ability to engineer flat bands and control the
location of geometric band features (e.g. Van Hove singularities and Lifshits
transitions), and thus to tune the relative strength of interaction-driven
physics.

This family of systems has already shown various correlated insulating states
Cao et al. (2018a); Lu et al. (2019), ferromagnetism Sharpe et al. (2019);
Chen et al. (2020), correlation-driven valley and iso-spin polarization Zhou
et al. (2022, 2021a), anomalous quantum hall physics Serlin et al. (2020),
topological insulator physics Wu et al. (2019a); Tschirhart et al. (2022);
Polshyn et al. (2020), metal-insulator transitions Li et al. (2021); Ghiotto
et al. (2021); Pan et al. (2020); Pan and Sarma (2021); Morales-Durán et al.
(2021); Ahn and Sarma (2022), possible “strange metal” resistance scaling at
very low temperature Jaoui et al. (2021); Polshyn et al. (2019a); Cao et al.
(2020b); Sarma and Wu (2022), and most conspicuously, possibly-exotic
superconductivity Cao et al. (2018b); Wu et al. (2018); Zhou et al. (2021b, a,
2022); Zhang et al. (2022); Polski et al. (2022); Arora et al. (2020),
including phases with verified non-spin-singlet pairing Zhou et al. (2022,
2021b). The rich phase diagrams and high experimental control that
characterize these systems has quickly made them into one of the most studied
platforms in condensed matter physics. The above-listed discoveries
demonstrate that geometric band features can have a profound influence on the
effects of interactions on transport properties. In turn, this highlights the
need for a refinement of the basic theories of phonon-limited resistivity as
applied to these materials, accurately taking complex band geometry into
account.

Figure 1: We depict the lattice structure of Bernal bilayer graphene (top) and
rhombohedral trilayer graphene (bottom). The left side of the image shows top-
down views of the xy-plane, labelling atoms with their layer number
($\\{1,2\\}$) and sublattice index ($\\{A,B\\}$). The right side of the figure
shows the stacking from a cross-section view. Figure 2:  Overview of phonon
scattering in Bernal bilayer and rhombohedral trilayer layered graphene
systems. The top row (a,e) depicts qualitatively distinct kinematically-
allowed scattering manifolds for different Bloch states at the same energy. In
these figures, the black curve depicts a Fermi surfaced (near the hole-doped
VHS), the pink dot denotes a reference Bloch state, and the colored points
denote the points in $\bm{\mathrm{k}}$-space that the reference state can
scatter too while conserving conservation of energy and momentum. These are
the “scattering manifolds” [Sec. III.2], which depend on the geometry of the
system. The coloring of the scattering manifolds encodes the transitions rates
from the reference state. These plots demonstrate the nontrivial kinematics
and geometry at play in scattering in these systems. The middle layers
(b,c,f,g) shows the density of states of the two systems, with labels showing
how the Fermi surface geometry changes as the sample is doped. The bottom
layer (d,h) gives the central results of this work, the scaling of the
resistivity with temperature in various regions of $n-T$ space. This is
captured by a heat map of $d\log[\rho(n,T)]/d\log T.$ We emphasize the clear
connection between the scaling behavior and the geometric features in the
density of states. The doted line gives the naive $T^{*}_{BG}$ calculated with
Eq. (1) for an isotropic system. We see that at sufficiently large dopings,
the color contours begin to follow the $\sqrt{|n|}$ profile traced out by the
black dotted $T^{*}_{BG}$ line, reflecting the fact that at large dopings, the
Fermi surface becomes roughly circular. However, our system exhibits a
surprising suppression of $T^{*}_{BG}$ as doping is decreased and the Fermi
surface qualitatively changes. These results are calculated with an inter-
layer potential of $\Delta=0.07\ eV$, and should be compared with Fig. 3,
which treats the zero-field case. The analogous results for a simple Dirac
cone (gapped and ungapped) are given in Fig. 4 for further comparison.  Figure
3:  Resistivity scaling with $T$ due to phonon scattering in ungapped BBG and
RTG, as in the absence of the inter-layer potential, to be compared with Fig.
2. As in Fig. 2, the top layers (a,b,d,e) plot the density of states of the
system for various doping levels, labelled with Fermi surface geometries. The
bottom row (c,f) provides a heat map of $d\log[\rho]/d\log T$ over $n-T$
space, mapping out the various regimes of resistivity scaling. As in the case
of the gapped systems (Fig. 2), we see that the BG transition is strongly
affected by the band geometry. Though this is more subtle without the applied
field, we still see the effects clearly in the RTG case, which still exhibits
an annular Fermi surface over a small doping window. As before, the dotted
line gives the expected BG crossover for an isotropic system; for the ungapped
systems plotted here, this estimate is quite accurate as long as the sample is
sufficiently doped.  Figure 4:  Resistivity scaling with $T$ due to phonon
scattering in simple gapped (right) and ungapped (left) Dirac cones, to be
compared with Figs. 2 and 3. As in those figures, the bottom row provides a
heat map of $d\log[\rho]/d\log T$ over $n-T$ space. The top row plots the
density of states of the system for various doping levels. The dotted line
gives $T^{*}_{BG}$ [Eq. (1)].

In particular, recent experiments in ABC-stacked rhombohedral trilayer
graphene (RTG) and AB-stacked Bernal bilayer graphene (BBG) (Fig. 1) have
discovered superconductivity (SC) proximate to several correlated, iso-spin
polarized phases Zhou et al. (2022, 2021a, 2021b); Zhang et al. (2022) in the
vicinity of Van Hove singularities and Lifshitz transitions in the band
structures. Additionally, there is evidence that some of the superconducting
phases host unconventional, non-spin-singlet pairing. Theories of SC in RTG
and BBG based on Cooper pairing mediated by interaction with acoustic phonons
have been put forth that propose likely explanations for the SC, explaining
the presence of both spin-triplet and spin-singlet phases and providing
roughly accurate transition temperatures Chou et al. (2022a, 2021, b).
Additionally, the proximity of SC phases to various interaction-driven phases
has spurred comparison to strong correlation physics, and several other
explanations centering $e-e$ interactions have been proposed You and
Vishwanath (2021); Ghazaryan et al. (2021); Szabó and Roy (2022); Cea et al.
(2022); Dong and Levitov (2021); Chatterjee et al. (2021); Szabó and Roy
(2022); Qin et al. (2022); Dai et al. (2022).

A time-tested method for ascertaining the relative importance of various
scattering mechanisms in a material is to look for clues in the temperature
dependence of the resistivity. This is because different mechanisms generally
produce various characteristic contributions and finite-$T$ crossovers between
these. For example, this debate is currently unfolding for twisted bilayer
graphene, where it is still unclear whether observed linear-in-$T$ resistivity
dependence is caused by phonons or an interaction-driven strange metal state,
analogous to that famously seen in several highly-correlated systems Sarma and
Wu (2022); Jaoui et al. (2021); Polshyn et al. (2019a); Cao et al. (2020b).

The recent experiments in BBG and RTG show that moiré-induced correlation
effects are not a necessary ingredient for SC in layered graphene systems,
leaving phonon-induced pairing as the de-facto leading candidate for a
universal SC mechanism in these systems. Especially since acoustic phonons
give a consistent theory of SC in both BBG and RTG, it is important to
understand and isolate the contribution to the resistivity that should be
expected due to acoustic phonons in the absence of $e-e$ effects. In
conventional superconductors, electron-phonon couplings extracted from SC tend
to agree well with those extracted from transport measurements. Thus, an
extensive quantitative comparision of the SC data and the transport data is an
important step in elucidating the nature of the SC pairing. Further, since
these systems demonstrate that superconductivity in 2D layered systems can be
intertwined with the nontrivial Fermi surface geometry, they offer an arena to
understand the extent to which these geometric features effect transport
generally.

The general paradigm of acoustic-phonon-limited resistivity in isotropic
(semi)metals is as follows Hwang and Sarma (2008); Min et al. (2011); Wu et
al. (2019b); Li et al. (2020); Hwang and Sarma (2019a); Ziman (1960); Ashcroft
and Mermin (1976). In the low-T regime, where the quantum statistics of the
phonon are important, we expect $\rho\approx T^{d+2}$, where $d$ is the
dimension of the sample. This characterizes the “Bloch-Grüneisen” (BG) regime,
which corresponds to

$\displaystyle k_{B}T\ll k_{B}T^{*}_{BG}=\mathcal{C}_{BG}\cdot(2\hbar
v_{p}k_{F}),$ (1)

where $v_{p}$ is the phonon velocity, $k_{F}$ is the Fermi momentum, and
$\mathcal{C}_{BG}\approx\mathcal{O}(1)$ is a material-specific constant.
(Further, $k_{B}T_{BG}\equiv 2\hbar v_{p}k_{F}$ is usually defined.) In the
high-$T$ “equipartition” (EP) regime, $T>>T^{*}_{BG}$, we instead expect
linear-in-T resistivity. We note that single-layer graphene displays these
properties elegantly, with $\mathcal{C}_{BG}\approx 1/6$ Hwang and Sarma
(2008); Efetov and Kim (2010).

The goal of this paper is to give a precise theoretical calculation of the
resistance due to acoustic phonon scattering in BBG and RTG systems in the
presence of an inter-layer potential ($\Delta$). We give concrete predictions
for the doping ($n$) and temperature ($T$) dependence of the resistivity of
these systems in the limit of phonon-dominated transport. The inter-layer
potential (produced by a displacement field) is required to induce SC in BBG,
and tuning this potential can significantly alter the band structure and
control the location of the Van Hove singularities, affecting both the SC and
the interaction-driven phases Zhou et al. (2022, 2021a).

Figure 5:  We plot a heat map of $\log[\rho]$ over doping density and
temperature for the BBG (left), RTG (middle), and Dirac (right) systems for
both the gapless (top) and gapped (bottom) cases. For BBG and RTG, the gapped
systems correspond to a displacement field $\Delta=0.07\ eV$, and the gapped
Dirac system has $M=0.05\ eV$. Several features are prominent. The dark blue
in the lower corners shows the universal features of the BG-transition. All
systems but the ungapped Dirac cone display density dependence throughout the
high-$T$ regime, with a significant increase in resistivity near charge
neutrality. This is especially prominent for the gapped RTG and BBG systems,
where the resistivity spikes to $\rho\approx 3000\Omega$.  Figure 6:  Accurate
approximate calculations of resistivity in the equipartition regime can be
made efficiently with the protocol discussed in Sec. III.3. Here we plot these
approximate results over a large range of $T$, extending the scope of Fig. 5.
We see that the large spike in resistivity near charge neutrality is a
relatively low-$T$ behavior and that the resistivity decreases at higher
temperatures, as $T$ get high enough to excite carriers in the conduction
band. We emphasize that even at high-$T$, the resistivity does not return to
simple linear scaling, but instead asymptotes to a constant value. [See also
Figs. 7 and 16.]

We carry out our calculation in the framework of Boltzmann kinetic theory,
treating the acoustic phonons via the Debye approximation but retaining the
full electronic band structure obtained by the diagonalization of
$\bm{\mathrm{k}}\cdot\bm{\mathrm{p}}$ Hamiltonians Jung and MacDonald (2014);
Zhang et al. (2010). We are able to numerically solve the linearized Boltzmann
equation in the anisotropic band geometry and give quantitative predictions
for the resistance and thus for the BG crossover temperature, $T^{*}_{BG}.$ We
emphasize that accurately treating the non-isotropic band structure is a
significant technical complication, beyond the techniques of prominent earlier
treatments of resistivity in 2D layered graphene structures Hwang and Sarma
(2008); Min et al. (2011); Wu et al. (2019b); Li et al. (2020); Hwang and
Sarma (2019a). Further, these earlier treatments of multi-layer graphene do
not include the effects of the inter-layer potential.

We find that the electronic structure of the layered systems significantly
distorts the BG paradigm explained above. In particular, while the high-$T$
behavior of the scattering rate of an individual Bloch state is linear,
$1/\tau_{\bm{\mathrm{k}}}\propto k_{B}T$, band curvature effects can lead to a
complicated non-linear T-dependence of the resistivity curves. In the cases of
gapped systems, (the displacement field generating $\Delta>0$ opens a gap),
there is a large spike in resistivity near charge neutrality. These band-
curvature effects interfere with the BG crossover, and we find that

Figure 7:  Resistivity data for hole-doped Bernal bilayer (top) and
rhombohedral trilayer (bottom) graphene stacks, evaluated with the
displacement field at $\Delta=0.07\ eV$. Resistivity is given in ohms on a
linear scale. The leftmost two columns give the results of our full numerical
calculation for the resistivity of the two systems at various doping levels up
to $30K$ and $120K$, respectively. The third column gives the same resistivity
curves extended to $800K$, making use of the equipartition assumption
discussed in Sec. III.4. The far right columns indicate the doping levels of
the curves shown in each row. From the full low-$T$ results (left) we may
extract the effective BG crossover. We see that for large dopings, $\rho(T)$
exhibits a BG-EP crossover temperature $T^{*}_{BG}$ as high as $40-60K$;
however, there is a sharp drop in $T^{*}_{BG}$ to around $20K$ upon the
Lifshitz transition to the annular Fermi surface. The $T^{*}_{BG}$ continues
to drop as we approach charge neutrality, dropping as low as $5-10K$. We note
that there is no sharp or discontinuous behavior at the Van Hove
singularities. The high-$T$ equipartition results (center) show how band
curvature effects lead to nonlinearity in $\rho(n,T)$ and how the effects can
complicate the BG crossover. We see also that for dopings close to charge
neutrality, we should expect a large spike in resistivity at moderate
temperatures ($100-300K$).

the approximate power law for resistivity scaling is strongly influenced by
the band structure geometry. This is demonstrated in Figs. 2 and 3 (and most
others in this paper). Further, we note that the anisotropy (i.e., trigonal
warping in graphene systems) in the band structure alters the low-$T$ BG
relaxation rate $T^{4}$ power law to a non-universal,
$\bm{\mathrm{k}}$-dependent $T$-dependence. While this nonlinear-in-$T$
equipartition-regime phonon-limited resistivity is unexpected in the context
of Boltzmann theory, we note that it has been detected experimentally in both
bilayer and trilayer twisted graphene systems Polshyn et al. (2019b);
Siriviboon et al. (2021).

Our paper is organized as follows. In Sec. II, we present an overview of the
main results of the work, emphasizing the most important quantitative aspects
for comparison with experiment and qualitative results that run counter to
common expectations. We then provide a concise review of acoustic phonon
scattering in kinetic theory and present an overview of the calculation of
relaxation times in the BBG and RTG systems in Sec. III. We emphasize the
roles of anisotropy and band curvature, which requires more care than the case
of an isotropic band. The non-linear $T$-dependence we report in the
equipartition regime is unexpected ; Section IV provides more intuition for
these effects. Our concluding discussion is presented in Sec. V.

Some supporting details are relegated to appendices. Appendix A presents the
$\bm{\mathrm{k}}\cdot\bm{\mathrm{p}}$ Hamiltonians used to calculate the band
structure of BBG and RTG. Appendix B discusses the numerical solution for the
relaxation rates in the solution of the linearized Boltzmann equation.
Appendix C discusses the role of the relaxation time approximation in
resistivity calculations. Finally, Appendix D presents additional $\rho(n,T)$
data for the systems of interest, supplementing the results presented in Sec.
II.

##  II Summary of main results

Our central results are the calculations of the doping ($n$) and temperature
($T$) dependence of the resistivity [$\rho(n,T)$] for Bernal bilayer and
rhombohedral trilayer graphene in the presence of a displacement field, under
the assumption that scattering is limited to acoustic phonons (which we treat
in the Debye approximation.) In particular, we give quantitative predictions
for the crossover from the Bloch-Grüneisen regime to the equipartition regime.

We plot $log[\rho(n,T)]$ for low $T$ ($0-120K$) in Fig. 5 and an approximate
extension of these results to higher-$T$ ($0-800K$) in Fig. 6. Individual
curves of $\rho(n,T)$ for fixed $n$ are given in Figs. 7. The most obvious
feature in this data is a strong spike in resistivity ($\rho\approx
3000\Omega$) near to charge neutrality at low $T$. From Figs. 6 and 7, we see
that this is a low-$T$ phenomenon and that resistivity drops and levels out at
higher $T$. However, we note that the high-$T$ resistivity is definitely not
given by a simple $T$-linear power law above the BG regime. In Figs. 2, 3, and
4, we plot $d\log[\rho(n,T)]/d\log(T)$ as an approximate scaling exponent for
the resistivity. These plots act as a sort of “phase diagram” for the various
regimes of $T$-dependence in the resistivity profile. In particular, we find
there is a region where the resistivity curve flattens out to be essentially
constant with $T$, sometimes after a downturn. While this is counter to
high-$T$ phonon expectations, this behavior has been measured in twisted
bilayer Polshyn et al. (2019b) and trilayer Siriviboon et al. (2021) graphene
systems. We stress that this is an effect entirely due to band curvature,
which we discuss further in Sec. IV.

Figures 2, 7, 5, 3, and 4 all demonstrate the BG crossover mentioned in the
introduction. At high dopings, where the Fermi surface is roughly circular, we
find a $T^{*}_{BG}\approx 40-60K$, in line with expectations for a circular
Fermi surface Hwang and Sarma (2008); Min et al. (2011); Wu et al. (2019b); Li
et al. (2020); Hwang and Sarma (2019a). However, we see a sharp drop to around
$20K$ at the Lifshitz transition to an annular Fermi surface, and $T^{*}_{BG}$
continues to drop as we approach charge neutrality. From Figs. 2, 3, and 4, it
is clear that the band geometry created by applying a displacement field
($\Delta$) to the graphene layers causes significant alterations to the
standard BG transition profile. Additionally, the curve-flattening discussed
in the last paragraph can come into effect at $T$ comparable to the crossover
temperature $T^{*}_{BG}$, making the $T^{4}\rightarrow T$ transition difficult
to observe.

Nevertheless, we predict that the phonon contribution to resistivity should
become important at temperatures that vary between $10K$ and $60K$, depending
on the doping, as shown in Fig. 7. This should be compared with what is
currently known from experiment: linear-in-$T$ resistivity dependence has not
been observed under 20$K$ in RTG or under 1.5$K$ in BBG. It is important to
note that the zero-$T$ contribution to resistivity from disorder ranges from
about 30$\Omega$ to 70$\Omega$ in these systems Zhou et al. (2022, 2021a,
2021b); Zhang et al. (2022).

We also report results for phonon scattering in BBG and RTG in the absence of
the displacement field. The data is all given in Fig. 6, and the effective
resistivity power law is extracted in Fig. 3, which should be compared with
Fig. 2. We note that in the absence of the applied field, the high-resistivity
spike near charge neutrality is significantly diminished. However, we still
find high-$T$ nonlinearity in the resistivity curves. Resistivity curves for
the ungapped cases analogous to Fig. 7 can be found in Appendix D.

##  III Resistivity via Boltzmann kinetic theory

The use of Boltzmann kinetic theory to calculate linear response resistivities
due to phonon collisions with Bloch state electrons is well-established
Ashcroft and Mermin (1976); Ziman (1960); Hwang and Sarma (2008); Min et al.
(2011). In this section, we outline the structure of the theory and explain
our calculation, appealing to the Dirac cone of single-layer graphene to
display concepts and highlight departures of our theory from previous work. We
first introduce the model in Sec. III.1, then we state the main results of the
kinetic theory in Sec. III.2 and use these results to give intuition into the
Bloch-Grüneisen crossover in Sec. III.3. Finally, Sec. III.4 discusses the
actual computation of the resistivity.

###  III.1 Model

We use the electronic single-particle Hamiltonian

$\displaystyle H^{e}$
$\displaystyle=\frac{1}{L^{2}}\sum_{\bm{\mathrm{k}}}c^{\dagger}_{\bm{\mathrm{k}}}H^{e}_{\bm{\mathrm{k}}}c_{\bm{\mathrm{k}}},$
(2)

where $c^{\dagger}\equiv c^{\dagger}_{s,\xi,\sigma,l,\bm{\mathrm{k}}}$ creates
an electron with crystal momentum $\bm{\mathrm{k}}$ (relative to Dirac point),
spin $s$, valley $\xi$, sublattice $\sigma$, and layer $l$. In our models,
$H^{e}_{\bm{\mathrm{k}}}\equiv\delta_{s,s^{\prime}}\delta_{\xi,\xi^{\prime}}H^{e}_{\sigma,l,\sigma^{\prime},l^{\prime},\bm{\mathrm{k}}}$
is a $\bm{\mathrm{k}}$-dependent matrix coupling together layer and sublattice
degrees of freedom, which are given in Appendix A. This is a
$\bm{\mathrm{k}}\cdot\bm{\mathrm{p}}$ continuum Hamiltonian from Jung and
MacDonald (2014); Zhang et al. (2010), which is very accurate within $1eV$ of
the charge neutrality point. The four degenerate spin-valley flavors remain
decoupled in our calculation and contribute equally to the conductivity
(inverse resistivity).

We are interested in the effects of the electron bands, so we restrict our
model to in-plane longitudinal acoustic phonons and adopt a simple Debye
description. We thus take the phonon Hamiltonian to be

$\displaystyle
H^{p}=\sum_{l,\bm{\mathrm{q}}}\hbar\omega_{\bm{\mathrm{q}}}a_{l,\bm{\mathrm{q}}}^{\dagger}a_{l,\bm{\mathrm{q}}},$
(3)

where $\omega_{\bm{\mathrm{q}}}$ is the phonon dispersion and we use the Debye
approximation $\omega_{\bm{\mathrm{q}}}\approx v_{p}|\bm{\mathrm{q}}|$, where
$v_{p}$ is the phonon velocity. This treatment neglects optical phonons, which
should give a quantitative correction above some temperature. Since optical
phonons have a large excitation gap in graphene, ranging from about $0.15$ to
$0.20eV$ Sohier et al. (2014), they will become relevant at higher
temperatures than we are concerned about here (approximately $1500K$) Sohier
et al. (2014); Xie and Foster (2016); Ghahari et al. (2016). Our neglect of
optical phonons is further justified by the fact that the electron-optical-
phonon couplings are weak in graphene multilayers due to sublattice
polarization Wu et al. (2018).

We couple the electrons to phonons via the well-known deformation potential
coupling Hamiltonian Ziman (1960); Hwang and Sarma (2008); Coleman (2015):

$\displaystyle H^{epc}$
$\displaystyle=\sqrt{\frac{D^{2}\hbar}{2\rho_{M}L^{2}}}\sum_{l,\bm{\mathrm{q}}}\frac{\hat{n}_{\bm{\mathrm{q}},l}}{\sqrt{\omega_{\bm{\mathrm{q}}}}}(-i\bm{\mathrm{q}}\cdot\hat{e}_{\bm{\mathrm{q}}})(a_{\bm{\mathrm{q}},l}+a_{-\bm{\mathrm{q}},l}^{\dagger}).$
(4)

Above, $D$ is the deformation potential, $\rho_{M}$ is the mass density of
monolayer graphene, and $\hat{e}_{\bm{\mathrm{q}}}$ is the desplacement unit
vector of the phonon. Throughout this work, we set $D=25$ eV,
$\rho_{M}=7.6\cdot 10^{-8}g/cm^{2}$, and $v_{p}=2.6\cdot 10^{6}cm/s$ Hwang and
Sarma (2008); Min et al. (2011); Wu et al. (2019b); Efetov and Kim (2010).
Finally, the electron density operator is

$\displaystyle\hat{n}_{\bm{\mathrm{q}},l}\equiv\sum_{\bm{\mathrm{k}}}c^{\dagger}_{(\bm{\mathrm{k}}+\bm{\mathrm{q}}),l}c_{\bm{\mathrm{k}},l}.$
(5)

In Eqs. (2) and (5), sums over unwritten $s,\sigma,\xi,l$ indices are
implicit.

Figure 8:  We plot the kinematically-allowed scattering manifolds for phonon-
scattering on a Dirac cone, given by Eq. (10). The Fermi surface at
$\mu=-0.25eV$ is plotted in black and a reference point on the Fermi surface
is identified with a pink dot. The colored points mark the set of k-space
points that the reference point can scatter too while conserving energy and
momentum. The color coding on the scattering manifold are proportional to the
scattering rate between the two points. The figure is shown for descending
temperatures: $300K$, $100K$, $30K$, $10K$. The Bloch-Grüneisen transition is
demonstrated by the fact that the $300K$ and $100K$ figures (top) only differ
quantitatively by the scale of the color bar, while they are qualitatively
distinct from the lower-temperature versions (bottom).

###  III.2 Kinetic theory

In the so-called “relaxation time approximation” Ashcroft and Mermin (1976)
[also see Appendix C] to Boltzmann kinetic theory, the resistivity tensor
($\rho$) is given by

$\displaystyle[\rho^{ij}(n,T)]^{-1}$
$\displaystyle=\frac{4e^{2}}{TL^{2}}\sum_{\bm{\mathrm{k}}\in
BZ}\tau_{\bm{\mathrm{k}}}v^{i}_{\bm{\mathrm{k}}}v^{j}_{\bm{\mathrm{k}}}f(\varepsilon_{\bm{\mathrm{k}}})[1-f(\varepsilon_{\bm{\mathrm{k}}})],$
(6)

where $T$ is temperature, $L$ is system length, $e$ is the electron charge,
$v_{\bm{\mathrm{k}}}^{j}$ are components of the velocity of the Bloch state
$\bm{\mathrm{k}}$, $f(\varepsilon)$ is the Fermi distribution function, and
the $\tau_{\bm{\mathrm{k}}}$ are the relaxation times of the various Bloch
states. If the band structure and Bloch states are known, the main challenge
in the computation of the resistivity is the computation of the relaxation
times. The leading factor of $4$ follows from the spin and valley degeneracies
of the problem.

In Eq. (6), we have suppressed the band index ($b$) and taken the sum over
$\bm{\mathrm{k}}$ to mean a sum over all Bloch states:
$\bm{\mathrm{k}}\rightarrow(\bm{\mathrm{k}},b)$. We will continue to use this
notation and will explicitly mention when interband excitations or transitions
are important.

Enforcing self-consistency of the relaxation time approximation on the
Boltzmann equation [Appendix (C)], we find that

$\displaystyle\frac{1}{|\bm{\mathrm{v}}_{\bm{\mathrm{k}}}|L^{2}}\sum_{\bm{\mathrm{k^{\prime}}}\in
BZ}\mathcal{W}_{\bm{\mathrm{k}}\rightarrow\bm{\mathrm{k^{\prime}}}}\frac{1-f(\varepsilon_{\bm{\mathrm{k^{\prime}}}})}{1-f(\varepsilon_{\bm{\mathrm{k}}})}\left[l_{\bm{\mathrm{k}}}-l_{\bm{\mathrm{k^{\prime}}}}\cos\theta_{\bm{\mathrm{v}}}\right]$
$\displaystyle=1,$ (7)

where
$l_{\bm{\mathrm{k}}}\equiv|\bm{\mathrm{v}}_{\bm{\mathrm{k}}}|\tau_{\bm{\mathrm{k}}}$
are the “relaxation lengths” (mean free paths), $\theta_{\bm{\mathrm{v}}}$ is
the angle between the Bloch velocities $\bm{\mathrm{v}}_{\bm{\mathrm{k}}}$ and
$\bm{\mathrm{v}}_{\bm{\mathrm{k^{\prime}}}}$, and
$\mathcal{W}_{\bm{\mathrm{k}}\rightarrow\bm{\mathrm{k^{\prime}}}}$ is the
transition rate from state $\bm{\mathrm{k}}$ to $\bm{\mathrm{k^{\prime}}}$. In
the thermodynamic limit, Eq. (7) becomes an integral equation. For a finite-
size system, it is a matrix equation that can be inverted to find the
relaxation lengths. Again, band indices have been suppressed, but
$\bm{\mathrm{k}}$ and $\bm{\mathrm{k^{\prime}}}$ should be taken to stand for
the Bloch states $(\bm{\mathrm{k}},b)$ and
$(\bm{\mathrm{k^{\prime}}},b^{\prime})$.

In the case of the standard deformation potential phonon coupling Hamiltonian,
a standard Fermi’s golden rule calculation gives the transition rates

$\displaystyle\mathcal{W}_{\bm{\mathrm{k}}\rightarrow\bm{\mathrm{k^{\prime}}}}$
$\displaystyle=\frac{\pi
D^{2}}{\rho_{M}v_{p}}|\bm{\mathrm{k^{\prime}}}-\bm{\mathrm{k}}|\Delta(\varepsilon_{\bm{\mathrm{k}}},\varepsilon_{\bm{\mathrm{k^{\prime}}}})\sum_{l}\bigg{|}\langle\psi_{\bm{\mathrm{k^{\prime}}},l}|\psi_{\bm{\mathrm{k}},l}\rangle\bigg{|}^{2}$
$\displaystyle\equiv\hbar
v_{p}|\bm{\mathrm{q}}|\Delta(\varepsilon_{\bm{\mathrm{k}}},\varepsilon_{\bm{\mathrm{k^{\prime}}}})\mathcal{C}_{\bm{\mathrm{k}},\bm{\mathrm{k^{\prime}}}}.$
(8)

with

$\displaystyle\bm{\mathrm{q}}$
$\displaystyle\equiv\bm{\mathrm{k^{\prime}}}-\bm{\mathrm{k}},$ (9)
$\displaystyle\Delta(\varepsilon_{\bm{\mathrm{k}}},\varepsilon_{\bm{\mathrm{k^{\prime}}}})$
$\displaystyle\equiv\begin{aligned}
N_{\bm{\mathrm{q}}}&\delta(\varepsilon^{\prime}-\varepsilon-\hbar
v_{p}|\bm{\mathrm{q}}|)\\\
+(N_{\bm{\mathrm{q}}}+1)&\delta(\varepsilon^{\prime}-\varepsilon+\hbar
v_{p}|\bm{\mathrm{q}}|)\end{aligned},$ (10) $\displaystyle
N_{\bm{\mathrm{q}}}$ $\displaystyle\equiv\frac{1}{\exp(\hbar
v_{p}|\bm{\mathrm{q}}|/k_{B}T)-1}.$ (11)

The Dirac $\delta$-functions in Eq. (10) enforce conservation of energy and
momentum and $N_{\bm{\mathrm{q}}}$ gives the occupation numbers of phonons
available for scattering. The first line in Eq. (10) refers to phonon
absorption processes while the second refers to phonon emission.

For a given energy band geometry, the conservation laws in Eq. (10) determine
a set of scattering manifolds for each Bloch state, corresponding to
absorption and emission of phonons. The summand in Eq. (7) then determines the
rate of transition to each point on the scattering manifold. Written as a sum
over the scattering manifold ($SM$), Eq. (7) takes the form

$\displaystyle\frac{\hbar
v_{p}}{|\bm{\mathrm{v}}_{\bm{\mathrm{k}}}|L^{2}}\sum_{\bm{\mathrm{k^{\prime}}}\in
SM}|\bm{\mathrm{q}}|\mathcal{C}_{\bm{\mathrm{k}},\bm{\mathrm{k^{\prime}}}}\mathcal{F}_{\bm{\mathrm{k}},\bm{\mathrm{k^{\prime}}}}^{\mu,T}\left[l_{\bm{\mathrm{k}}}-l_{\bm{\mathrm{k^{\prime}}}}\cos\theta_{\bm{\mathrm{v}}}\right]$
$\displaystyle=1,$ (12)

with

$\displaystyle\mathcal{F}_{\bm{\mathrm{k}},\bm{\mathrm{k^{\prime}}}}^{\mu,T}\equiv\frac{1-f(\varepsilon_{\bm{\mathrm{k^{\prime}}}})}{1-f(\varepsilon_{\bm{\mathrm{k}}})}\times\left\\{\begin{aligned}
N_{\bm{\mathrm{q}}}\ \ \ \
\varepsilon_{\bm{\mathrm{k^{\prime}}}}>\varepsilon_{\bm{\mathrm{k}}}\\\
N_{\bm{\mathrm{q}}}+1\ \
\varepsilon_{\bm{\mathrm{k^{\prime}}}}<\varepsilon_{\bm{\mathrm{k}}}\end{aligned}\right\\},$
(13)

where $\bm{\mathrm{k}}\in SM$ indicates a summation over the scattering
manifold of states picked out by the delta functions in Eq. (10). We emphasize
that all implicit dependence of the relaxation lengths on the temperature or
chemical potential are due to
$\mathcal{F}_{\bm{\mathrm{k}},\bm{\mathrm{k^{\prime}}}}^{\mu,T}$.

From Eqs. (6-8) we see that resistance scales linearly with $D^{2}/\rho_{M}$,
so our results are easy to adjust for different values of these parameters.
The dependence on $v_{p}$ is more involved, since it also affects the geometry
of the scattering manifolds.

Fig. 8 shows how scattering rates can vary across the scattering manifold,
using a Dirac cone as a simple example. It also visually demonstrates the
transition between the BG and EP regimes, which we discuss next.

Figure 9:  We plot the log of the relaxation lengths
$\\{l_{\bm{\mathrm{k}}}\\}$ for a Dirac cone band structure for various
temperatures. We fix the chemical potential at $\mu=-0.25\ eV$ and perform the
calculation at $T=10\ K$ (top), $T=30\ K$ (middle), and $T=100\ K$ (bottom).
The far left column shows the relaxation lengths plotted over energy, while
the central and rightmost columns give a heat map of the relaxation lengths in
momentum space for the bottom (center) and top (right) bands, respectively. We
emphasize that at low T, states near the Fermi level become long-lived.
Further, states near the Dirac point are always long-lived due to a vanishing
scattering manifold.

###  III.3 Bloch-Grüneisen and Equipartition regimes

The low-$T$ BG regime is best understood in the case of an isotropic
($l_{\bm{\mathrm{k}}}\rightarrow l_{\varepsilon_{\bm{\mathrm{k}}}}$ and
$\bm{\mathrm{v}}_{\bm{\mathrm{k}}}\parallel\bm{\mathrm{k}}$) and quasi-elastic
($\varepsilon^{\prime}\approx\varepsilon$) system, such as graphene Hwang and
Sarma (2008). In this case, we can replace the velocity angle with the
momentum angle ($\theta_{\bm{\mathrm{v}}}=\theta_{\bm{\mathrm{k}}}$) and Eq.
(7) simplifies to a direct formula for the relaxation time:

$\displaystyle\frac{1}{\tau_{\bm{\mathrm{k}}}}$ $\displaystyle=\frac{\hbar
v_{p}}{L^{2}}\sum_{\bm{\mathrm{k^{\prime}}}\in
FS}|\bm{\mathrm{q}}|\mathcal{C}_{\bm{\mathrm{k}},\bm{\mathrm{k^{\prime}}}}\mathcal{F}_{\bm{\mathrm{k}},\bm{\mathrm{k^{\prime}}}}^{\mu,T}\left[1-\cos\theta_{\bm{\mathrm{k}}}\right],$
(14)

where $\bm{\mathrm{k}}\in FS$ indicates a summation over the Fermi surface,
which is taken to be indistinguishable from the scattering manifold in the
quasi-elastic approximation.

For small $\bm{\mathrm{q}}$,
$1-\cos\theta_{\bm{\mathrm{k}}}\approx|\bm{\mathrm{q}}|^{2}$ and
$\mathcal{C}_{\bm{\mathrm{k}},\bm{\mathrm{k^{\prime}}}}\approx 1$, and the
summand of Eq. (14) scales with $\bm{\mathrm{q}}$ roughly as
$|\bm{\mathrm{q}}|^{3}$. For low $T$, the Fermi functions
$1-f(\varepsilon_{\bm{\mathrm{k^{\prime}}}})$ and the phonon occupation
function $N_{\bm{\mathrm{q}}}$ effectively restrict the sum in Eq. (14) to
$\bm{\mathrm{k^{\prime}}}$ with $\hbar v_{p}|\bm{\mathrm{q}}|\leq k_{B}T$.
Summing $|\bm{\mathrm{q}}|^{3}$ over the portion of the [$(d-1)$-dimensional]
scattering manifold within a radius proportional to $T$ gives the famous
power-law defining the BG regime:

$\displaystyle\frac{1}{\tau_{\bm{\mathrm{k}}}}\propto T^{d+2}.$ (15)

However, if we do not assume isotropy, then we must restore

$\displaystyle 1-\cos\theta_{\bm{\mathrm{k}}}\rightarrow
1-\frac{l_{\bm{\mathrm{k^{\prime}}}}}{l_{\bm{\mathrm{k}}}}\cos\theta_{\bm{\mathrm{v}}}$
(16)

in Eq. (14). The small-$|\bm{\mathrm{q}}|$ limit of the right hand side of Eq.
(16) is not necessarily proportional to $|\bm{\mathrm{q}}|^{2}$, since it
depends on the way $l_{\bm{\mathrm{k^{\prime}}}}\rightarrow
l_{\bm{\mathrm{k}}}$ and
$\bm{\mathrm{v}}_{\bm{\mathrm{k^{\prime}}}}\rightarrow\bm{\mathrm{v}}_{\bm{\mathrm{k}}}$
as $\bm{\mathrm{k^{\prime}}}\rightarrow\bm{\mathrm{k}}.$ We therefore expect
anisotropy to introduce non-universal, $\bm{\mathrm{k}}$-dependent
modifications of the BG power law in the $T$-dependence of each relaxation
time $\tau_{\bm{\mathrm{k}}}.$

In the high-$T$ limit, expanding in small $\Delta\varepsilon/T$, we find

$\displaystyle\mathcal{F}_{\bm{\mathrm{k}},\bm{\mathrm{k^{\prime}}}}^{\mu,T}=\frac{k_{B}T}{\hbar
v_{p}|\bm{\mathrm{q}}|}+\mathcal{O}(\Delta\varepsilon/T),$ (17)

and inserting into Eq. (7) gives

$\displaystyle\frac{k_{B}T}{|\bm{\mathrm{v}}_{\bm{\mathrm{k}}}|L^{2}}\sum_{\bm{\mathrm{k^{\prime}}}\in
SM}\mathcal{C}_{\bm{\mathrm{k}},\bm{\mathrm{k^{\prime}}}}\left[l_{\bm{\mathrm{k}}}-l_{\bm{\mathrm{k^{\prime}}}}\cos\theta_{\bm{\mathrm{v}}}\right]$
$\displaystyle=1+\mathcal{O}(\Delta\varepsilon/T)^{2}.$ (18)

Solving Eq. (18) order-by-order in $1/T$, we see that the high-$T$ form of the
relaxation length is

$\displaystyle l_{\bm{\mathrm{k}}}$
$\displaystyle=\frac{c_{\bm{\mathrm{k}}}}{k_{B}T}+\mathcal{O}(\Delta\varepsilon/T)^{3}.$
(19)

We note that the $\mathcal{O}(1)$ term in the $\Delta\varepsilon/T$ expansion
of $\mathcal{F}_{\bm{\mathrm{k}},\bm{\mathrm{k^{\prime}}}}^{\mu,T}$ in Eq.
(17) rather remarkably vanishes, preventing a
$\mathcal{O}(\Delta\varepsilon/T)^{2}$ term in Eq. (19). This implies that the
high-$T$ scattering rate (due to phonons) of a given Bloch state should be
purely linear, going to zero in the $T\rightarrow 0$ extrapolation.

The equipartition regime is the range of temperature for which Eq. (19) holds
for all Block states $\bm{\mathrm{k}}$. Unlike the case in the BG regime, the
linear-in-$T$ power law for the relaxation rate of the EP regime is not
affected by anisotropy - all band structure information is encoded in the
“length constants” $c_{\bm{\mathrm{k}}}$.

###  III.4 Resistivity computation

Equations (6-11) combined with knowledge of the Bloch states give all the
tools necessary to make a resistivity prediction. We solve Eqs. (7) for
scattering lengths for each Bloch state [see Appendix B for discussion.] We
emphasize that in general, the relaxation lengths $\\{l_{\bm{\mathrm{k}}}\\}$
implicitly depend on temperature and chemical potential through the Fermi
functions and phonon occupation number ($N_{\bm{\mathrm{q}}}$) in Eq. (7).
Once the $\\{l_{\bm{\mathrm{k}}}\\}$ are known for a given pair $(n,T)$, the
resistivity can be computed through Eq. (6). We plot the relaxation lengths
for a Dirac cone band structure in Fig. 9, keeping $\mu$ fixed as we vary $T$.
These results illustrate that states near the Fermi surface become long-lived
at low $T$.

In Secs. III.2-III.3, we have suppressed the band index in summations over
Bloch states. The Bernal bilayer and rhombohedral trilayer
$\bm{\mathrm{k}}\cdot\bm{\mathrm{p}}$ Hamiltonians have four and six bands,
respectively, while the Dirac cone model has two. In each case, we have two
“low energy” bands near charge neutrality: a “valence” (hole) band and a
“conduction” (particle) band. The higher energy bands, when present, are over
$3.5eV$ from charge neutrality. We note that the Fermi distributions in Eq.
(6) suppress excitations in these higher energy bands for the temperatures and
dopings we are interested in. However, it is important to keep both the
conduction and valence bands as charge carriers may be excited in both bands,
especially in the gapless systems. In all the models we study, interband
transitions between the conduction and valence bands are forbidden by
kinematics (i.e. the phonon velocity is too low). Interband transitions into
higher energy bands are kinematically allowed, but thermally irrelevant.

It is important to note that as we scan $T$ for fixed $n$, $\mu(n,T)$ can
change, and this can be quite drastic near a gap. We must therefore calculate
$\mu(n,T)$ self-consistently via

$\displaystyle n=\frac{4}{L^{2}}\sum_{\bm{\mathrm{k}}\in
BZ}f(\varepsilon_{\bm{\mathrm{k}}}).$ (20)

The prefactor 4 above follows from the spin and valley degeneracies. We stress
that accurately computing the $T$-dependence of $\mu(n,T)$ near the band edge
requires keeping both the valence and conduction bands, even if $T$ is far too
low to excite carriers across the gap.

Figure 10:  We plot heat maps of the relaxation length in k-space for the
Bernal bilayer band structure. We set $\mu=-0.066\ eV$ (left) and $\mu=-0.058\
eV$ (right). For each $\mu$, we do the calculation for $T=10\ K$ (top), $T=30\
K$ (middle), and $T=100\ K$ (bottom). and We see that at low $T$, the specific
geometry of the Fermi surface is very important to relaxation, but that this
information tends to get washed out at higher $T$.  Figure 11:  We plot the
equipartition “length constants” over a wide range of energy for the various
systems under study. The top row gives Bernal bilayer graphene, the middle row
gives rhombohedral trilayer graphene, and the bottom row gives standard Dirac
cone graphene. Gapless systems are on the left and gapped systems are on the
right. With these values stored we can efficiently compute the resistivity of
these systems up to very high temperatures using Eq. 19, though our results
will miss the low-$T$ BG physics, as discussed in Sec. III.3. This is how we
generate the high-T results in Figs. 7, 6. We emphasize that the ungapped
Dirac cone graphene, which has diverging $c_{\bm{\mathrm{k}}}$ near charge
neutrality due to a vanishing scattering manifold, is the outlier here. All
other systems we consider have band curvature effects near charge neutrality
that more than compensate for the vanishing scattering manifolds and suppress
the divergence of $c_{\bm{\mathrm{k}}}.$ Figure 12:  We plot the
“$g$-function”, [$g(\varepsilon)$], defined in Eqs. (21-22) for Bernal bilayer
(top), rhombohedral trilayer (center), and Dirac cone graphene (bottom). We
plot $g(\varepsilon)$ for both ungapped (left) and gapped (right) cases. These
graphs demonstrate clearly why single ungapped Dirac cone graphene has such a
robust high-$T$ linear resistivity and why all the other systems display
nonlinear resistivity effects at high-$T$.  Figure 13:  We plot heat maps of
$d\log[\rho(n,T)]/d\log T$ for the various systems using the resistivities
calculated in the high-$T$ equipartition calculation (the data from Fig. 6).
This should be compared with Figs. 2,3, and 4. Since this is based on
resistivity data from the equipartition regime, it does not contain any BG
physics, and any deviations from linear scaling are due to band curvature
effects. We emphasize that these figures show linear-in-$T$ scaling for
dopings away from charge neutrality, but then demonstrate a flattening of the
$\rho(T)$. We further emphasize that the flattening seen near the van-Hove
singularities is present here as well, indicating it is an effect of band
geometry in the thermal averaging and not a transition in the nature of the
scattering.

The main result of this work is the application of the above analysis to
Bernal bilayer and rhombohedral trilayer graphene stacks. These results are
presented and discussed in Sec. II. We use
$\bm{\mathrm{k}}\cdot\bm{\mathrm{p}}$ Hamiltonians for these systems Jung and
MacDonald (2014); Zhang et al. (2010), which we provide in Appendix A. The
band structure further gives the density of states and Fermi surface
geometries depicted in Figs. 2,3.

The nontrivial band geometry of these systems gives scattering manifolds that
depend qualitatively on not only the Fermi level, but also the specific Bloch
state in question, as depicted in Fig. 2. Since the bands are not isotropic
and the phonon scattering cannot be considered “quasi-elastic” Hwang and Sarma
(2008), we need to find the full solution of Eq. (7). Solving Eq. (7) for the
$\\{l_{k}\\}$ repeatedly for many values of $n$ and $T$, we calculate the
resistivity data given in Figs. 5,7. Data showing how scattering lengths vary
throughout the band structure are given in Fig. 10.

The equipartition regime scaling coefficients, $c_{\bm{\mathrm{k}}}$, are
given for all the models of interest in Fig. 11. In the case of Dirac cone
graphene, we see that there is a divergence of $c_{\bm{\mathrm{k}}}$ at the
Dirac point, arising from a vanishing set of scattering states. However, in
all the other models under consideration, band curvature effects near charge
neutrality more than compensate for the vanishing scattering manifolds and
suppress $c_{\bm{\mathrm{k}}}$.

##  IV Nonlinear $T$-dependence of resistivity

The “common knowledge” of high-$T$ phonon scattering is that the resistivity
scales linearly with $T$ above the BG crossover regime Ashcroft and Mermin
(1976); Ziman (1960); Hwang and Sarma (2008); Hwang and Sarma (2019a). While
it is true that each individual relaxation length has the high-$T$ scaling of
Eq. (19), the $T$-dependence of the resistivity itself can be quite nonlinear.
Indeed, our calculations for BBG and RTG predict a nonlinear $T$-dependence of
the resistivity, especially in the vicinity of the gap. [See Figs. 7,6.]

Our calculations predict that the phonon scattering will crossover from the BG
regime to the EP regime at an effective BG crossover temperature that varies
from as high as $40\ K$ at high doping to as low as $<10\ K$ near charge
neutrality. However, in the EP regime, we start to see sharp reductions in
slope of the resistivity at temperatures as low as $40\ K$ [See Fig. 7]. For
dopings closer to charge neutrality, we see the resistivity peak and drop
precipitously at $T\approx 300\ K$. This behavior has already been observed in
twisted bilayer graphene, Polshyn et al. (2019b), at temperatures and
resistivity values qualitatively consistent with our results here.

We note that the non-linear $T$-dependence resembles the same sort of
resistivity profiles that have been characterized as “resistivity saturation”
Poniatowski et al. (2021); Sarkar et al. (2018); Hwang and Sarma (2019b);
Emery and Kivelson (1995) and are sometimes associated with a breakdown of
kinetic theory at the Mott-Ioffe-Regel limit Gurvitch (1981); Millis et al.
(1999); Emery and Kivelson (1995); Hussey et al. (2004). However, we stress
that our results are fully in the Boltzmann framework. The possibility that
the apparent resistivity saturation type effect could arise purely from the
electron-phonon coupling effects was pointed out in the literature before
Millis et al. (1999), but the physics of this apparent saturation in the
current work is qualitatively different, arising not from non-Boltzmann strong
coupling physics, but from subtle band structure effects as discussed in our
paper.

In the rest of this section, we provide some intuition for the non-linear
$T$-dependence of the resistivity. As discussed above, the high-$T$ relaxation
lengths are given in terms of the $n,T$-independent constants
$\\{c_{\bm{\mathrm{k}}}\\}.$ We can gain an understanding of the non-linearity
of $\rho(T)$ by considering the function

$\displaystyle\delta^{ij}g(\varepsilon)$
$\displaystyle\equiv\frac{1}{L^{2}}\sum_{\bm{\mathrm{k}}}\frac{v^{i}_{\bm{\mathrm{k}}}v^{j}_{\bm{\mathrm{k}}}}{|v_{\bm{\mathrm{k}}}|}c_{\bm{\mathrm{k}}}\delta[\varepsilon-\varepsilon_{\bm{\mathrm{k}}}].$
(21)

In the EP regime, we can write the high-$T$ conductivity in terms of
$g(\varepsilon)$:

$\displaystyle\sigma=\frac{1}{(k_{B}T)^{2}}\int d\varepsilon
g(\varepsilon)f(\varepsilon)[1-f(\varepsilon)]$ (22)

We see that even in the equipartition regime, we only have linear scaling of
the resistivity if the integral over $\varepsilon$ in Eq. (22) scales linearly
with $T$, which will be true as long as $g(\varepsilon)$ has a good linear
approximation in a window of width $k_{B}T$ about $\mu(n,T).$ The functions
$g(\varepsilon)$ are plotted in Fig. 12. We see that gapless Dirac cone
graphene has a perfectly flat $g(\varepsilon),$ (though our figure shows
finite-size effects near the Dirac point), giving the familiar, perfectly
linear resistivity in the equipartition regime. On the other hand, we see that
gapless Dirac cone graphene is the exception - all of the other systems
studied exhibit band curvature that manifests nonlinearity in
$g(\varepsilon).$ The gapless bilayer and trilayer systems exhibit
$g(\varepsilon)$ that can be roughly approximated as linear over small
$\varepsilon$-windows when sufficiently doped. However, we expect a
qualitative change when $T\approx|\mu(n,T)|,$ and the integral in Eq. (22)
crosses the zero-energy point, where we expect the scaling of the integral in
Eq. (22) to crossover from linear-in-$T$ to quadratic-in-$T$. This would
result in a crossover to a roughly $T$-independent resistivity when
$T\approx\mu(n,T),$ which is indeed what we see in Fig. (3). All three gapped
systems exhibit more curvature in $g(\varepsilon),$ even when far from charge
neutrality, but may still be linearly approximated in a small $T$-window.
However, sharp qualitative changes in $g(\varepsilon)$ occur at a band edge,
so we expect sharp qualitative changes in the resistivity scaling when
$T\approx|\mu(n,T)|-E_{gap}/2$ and again when $T\approx|\mu(n,T)|+E_{gap}/2$.
For a system without a band gap, we would only expect a single kink. Figures 7
and 16 \- 15 demonstrate this intuition. Two distinct kinks are visible in
many resistivity curves in Figs. 7 and 16, which plot the data for gapped
systems, while curves in Figs. 14 and 15 tend to have a single kink.

We emphasize that the nonlinear $T$-dependence of the resistivity is in
general due to the curvature of the bands and not necessarily related to
interband excitations Polshyn et al. (2019b). For instance, in the hole-doped
systems in Fig. 7, with the potential difference at $\Delta=0.07\ eV$, the gap
is approximately $0.1\ eV$ wide. However, nonlinear $T$-dependence is seen at
temperatures as low as $40\ K$, which is far to cold to excite appreciable
states in the conduction band.

##  V Discussion and conclusions

We have calculated the electrical DC resistivity of Bernal bilayer and
rhombohedral trilayer graphene systems, due to scattering off of acoustic
phonons. We extend previous study by using a detailed
$\bm{\mathrm{k}}\cdot\bm{\mathrm{p}}$ band structure and focusing our
attention on the roles of geometric features of the band structure of these
systems, including those affected by a displacement field.

We develop a thoroughly nontrivial transport theory for carrier resistivity
due to electron-acoustic phonon interaction in experimentally relevant RTG and
BBG multilayer graphene systems. The theory, while using the standard graphene
acoustic phonons and the conventional electron-phonon deformation potential
coupling, includes the full effects of RTG and BBG band structures (even
including an applied electric displacement field) non-perturbatively by
employing a full $\bm{\mathrm{k}}\cdot\bm{\mathrm{p}}$ description. The
qualitative importance of the van Hove singularities and the anisotropies in
the graphene band structures are exactly incorporated in the theory by
iteratively solving the integral Boltzmann transport equation. This leads to
several qualitatively new features in the resistivity (e.g. inapplicability of
the simple Bloch-Grüneisen criteria for linear versus non-linear resistivity
in temperature, apparent resistivity saturation behavior at higher
temperatures, and other features as discussed in this paper), which have not
been discussed in the transport literature of electronic materials before in
any context. We provide concrete predictions for the doping and temperature
dependence of resistivity in RTG and BBG multilayers, finding that simple
considerations for a Bloch-Grüneisen temperature separating the linear-in-$T$
high-temperature resistivity from the non-linear low-temperature resistivity
does not apply because of the band geometry introducing strong modifications
of the resistivity behavior.

Our results are important in two contexts. First, the experiments in BBG and
RTG have shown that the exotic superconductivity and the various interaction-
driven correlated states are closely related to the nontrivial geometric
features of the band structures, including Fermi-surface reshapes and Van Hove
singularities. This spotlights the enhanced effects of band geometry on
scattering processes in complex 2D systems. As 2D layered heterostructures are
currently ascendant in condensed matter physics, it is important to study the
relationship between band geometry and transport directly and to modify
intuitions gained in three dimensions. Second, it is crucial in the
investigation of the origin of the superconductivity in moiréless layered
graphene systems to understand the relative importance of various scattering
mechanisms. Our work provides a clear and concrete picture of how the
resistivity should behave in a phonon-dominated system. If strong deviations
from these results are seen in experiment, that could serve as evidence that
the scattering mechanisms other than phonons play dominant roles in transport.
This would point to directions for non-phonon pairing in the observed
superconductivity.

The doping and temperature dependence of the resistivity of these systems
behave similarly and with many interesting features. We find that the BG
crossover in the qualitative $T$-dependence of the scattering rates varies as
a function of doping from as low as $5K$ to as high as $60K$. However, we note
that this crossover temperature depends strongly on the geometric features of
the band structure, and is sharply reduced by the emergence of the annular
Fermi surface, which is related to the observed SC. Further, we find that band
curvature effects also give rise to a non-linear $T$-dependence of the
resistivity at temperatures in the intermediate range of $60K-300K$. While our
results show an interesting sensitivity to changes in Fermi surface geometry,
they are remarkably smooth at the Van Hove singularity.

Our results are qualitatively compatible with what is currently known in
experiment Zhou et al. (2022, 2021a, 2021b). We have not yet seen evidence of
the high-$T$ nonlinear equipartition resistivity in BBG or RTG, but very
similar effects have been observed in twisted bilayer Polshyn et al. (2019b)
and trilayer Siriviboon et al. (2021) graphene systems. While the BG crossover
to linear scattering has not yet been observed in these systems at low
temperatures, our results show that current experiment cannot rule out the
possibility that these systems are dominated by phonon scattering. In
particular, no linear-in-$T$ region has been observed below $20K$ in RTG and
the zero temperature resistivity varies from $20-70\Omega$ (c.f. Fig.S6 in
Zhou et al. (2021a)). Our predictions are compatible with these experimental
results. However, our results do make it clear that extensive experimental
resistivity data over wider ranges of doping and temperature (from $0$ to
$300K$) should be sufficient to tell if there are strong deviations from the
phonon-dominated picture. Comparison of our results with future additional
experimental resistivity data could be a crucial step in discovering the
origin of SC in these systems. Further, the low Fermi velocities and high
density of states at the Van Hove singularity should enhance the effects of
electron-electron interactions. Since our calculations do not predict sharp
features to emerge at the Van Hove singularities in a purely phonon picture,
observations of such features in the resistivity could serve as evidence for
strong-coupling physics that could underlie the systems’ superconductivity.

###### Acknowledgements.

We thank Christopher David White, Jiabin Yu, Jay D. Sau, and Matthew S. Foster
for helpful discussion. This work is supported by the Laboratory for Physical
Sciences (S.M.D, Y.-Z.C., and S.D.S). It is also partially funded by JQI-NSF-
PFC (Y.-Z.C.). F.W. is supported by startup funding of Wuhan University.

##  Appendix A Hamiltonians for stacked graphene systems

To calculate the band structure for the Bernal bilayer graphene stack, we use
in Eq. (2) the Hamiltonian introduced in Jung and MacDonald (2014), and used
also in Zhou et al. (2022); Chou et al. (2022a, b):

$\displaystyle
H^{e}=\begin{bmatrix}-\Delta&v_{0}\bar{k}&-v_{4}\bar{k}&-v_{3}k\\\
v_{0}k&\Delta_{2}-\Delta&t_{1}&-v_{4}\bar{k}&\\\
-v_{4}k&t_{1}&\Delta_{2}+\Delta&v_{0}\bar{k}\\\
-v_{3}\bar{k}&-v_{4}k&v_{0}k&\Delta\end{bmatrix},$ (23)

where we use the dimensionless, valley-dependent, (anti)holomorphic momenta
$k\equiv a_{0}(\xi k_{x}+ik_{y})$ and $\bar{k}=a_{0}(\xi k_{x}-ik_{y})$ for
valley $\xi\in\\{\pm\\}$, where $a_{0}$ is the lattice constant for graphene
($a_{0}=0.246\ nm$). The parameters take the following values (all quantities
in $eV$): $\Delta_{2}=0.015,t_{1}=0.361,v_{0}=2.261,v_{3}=0.245,v_{4}=0.12$.
The interlayer potential is $\Delta$, and in our calculations this is either
set to $0.07\ eV$ or $0\ eV$. The basis for this matrix is
$\\{1A,1B,2A,2B\\},$ where $A,B$ correspond to sublattice and $1,2$ correspond
to layer.

For the rhombohedral trilayer stack, we use the Hamiltonian introduced in
Zhang et al. (2010), and used also in Zhou et al. (2021a, b); Chou et al.
(2021, 2022b):

$\displaystyle
H^{e}=\begin{bmatrix}\Delta_{2}+\Delta+\delta&\gamma_{2}/2&v_{0}\bar{k}&v_{4}\bar{k}&v_{3}k&0\\\
\gamma_{2}/2&\Delta_{2}-\Delta+\delta&0&v_{3}\bar{k}&v_{4}k&v_{0}k\\\
v_{0}k&0&\Delta_{2}+\Delta&\gamma_{1}&v_{4}\bar{k}&0\\\
v_{4}k&v_{3}k&\gamma_{1}&-2\Delta_{2}&v_{0}\bar{k}&v_{4}\bar{k}\\\
v_{3}\bar{k}&v_{4}\bar{k}&v_{4}k&v_{0}k&-2\Delta_{2}&\gamma_{1}\\\
0&v_{0}\bar{k}&0&v_{4}k&\gamma_{1}&\Delta_{2}-\Delta\end{bmatrix},$ (24)

where we use the same notation ($k,\bar{k}$) as in Eq. (23) and the following
parameters (all quantities in $eV$):
$\Delta_{2}=-0.0023,\delta=-0.0105,v_{j}=\gamma_{j}\sqrt{3}/2,\gamma_{0}=3.1,\gamma_{1}=0.38,\gamma_{3}=-0.29,\gamma_{4}=-0.141$.
Again, the interlayer potential is $\Delta$, and in our calculations this is
either set to $0.07\ eV$ or $0\ eV$. The basis for the RTG Hamiltonian is
$\\{1A,3B,1B,2A,2B,3A\\}.$

##  Appendix B Numerical implementation of resistivity calculation

In our numerical calculations for $\\{l_{\bm{\mathrm{k}}}\\}$, we usually
retain approximately $10^{6}$ Bloch states, and must solve a rather large
linear system [Eq. (7)] for each pair of values $(n,T)$. In our main results
[Figs.7,5], we do this on a 50-by-60 grid in $n-T$-space. This is necessary to
understand the low-$T$ physics, but the EP regime can be studied much more
efficiently since the $\\{c_{\bm{\mathrm{k}}}\\}$ defined in Eq. (19) are
independent of both $n,T$. Once we solve directly for the
$c_{\bm{\mathrm{k}}}$, calculating the EP approximation to the resistivity is
as simple as computing $\mu(n,T)$ via Eq. (20) and then using Eq. (19) in Eq.
(6). This is how we compute the EP resistivity in Figs.7 and 6.

We discuss the numerical solution of Eq. (7) in the main text. In order to
discuss the existence and uniqueness of solutions to Eq. (7), as well as the
convergence of iterative methods, we will re-cast this in the traditional
notation of a linear operator problem. Letting $\bm{\mathrm{k}}$ in the
Brillouin zone act as a vector index, we define the vector $\hat{b}$ and the
matrices $\hat{A},\hat{D}$, indexed by $\bm{\mathrm{k}}\in BZ$.

$\displaystyle\hat{A}_{\bm{\mathrm{k}},\bm{\mathrm{k^{\prime}}}}$
$\displaystyle=\mathcal{W}_{\bm{\mathrm{k}}\rightarrow\bm{\mathrm{k^{\prime}}}}[1-f^{0}(\varepsilon_{\bm{\mathrm{k^{\prime}}}})]\cos\theta_{\bm{\mathrm{v}}}$
(25) $\displaystyle\hat{D}_{\bm{\mathrm{k}},\bm{\mathrm{k^{\prime\prime}}}}$
$\displaystyle=\delta_{\bm{\mathrm{k}},\bm{\mathrm{k^{\prime\prime}}}}\sum_{\bm{\mathrm{k^{\prime}}}}\mathcal{W}_{\bm{\mathrm{k}}\rightarrow\bm{\mathrm{k^{\prime}}}}[1-f^{0}(\varepsilon_{\bm{\mathrm{k^{\prime}}}})]$
(26) $\displaystyle\hat{b}_{\bm{\mathrm{k}}}$
$\displaystyle=|\bm{\mathrm{v}}_{\bm{\mathrm{k}}}|L^{2}[1-f^{0}(\varepsilon_{\bm{\mathrm{k}}})]$
(27)

With this notation, Eq. (7) takes the form

$\displaystyle(\hat{D}-\hat{A})\hat{l}=\hat{b}.$ (28)

The solution for the relaxation lengths is then a matrix inversion problem. A
unique solution exists if $\det[\hat{D}-\hat{A}]\neq 0$, which is always true
in this case due to the diagonal dominance of $\hat{D}-\hat{A}$. Since our
problem is large and we compute the matrix elements only as needed in the
computation, Eq. (28) is most effectively solved via an iterative method. We
set

$\displaystyle\hat{l}^{i+1}\leftarrow\hat{D}^{-1}(\hat{A}\hat{l}^{i}+\hat{b})$
(29)

repeatedly until convergence. This is simply a case of Gauss-Seidel iteration,
which is guaranteed to converge to the unique solution. (This guarantee is
again provided by diagonal dominance.)

Explicitly, in the $(i+1)^{\text{th}}$ iteration ($i\geq 0$), we define
$\\{l^{(i+1)}_{\bm{\mathrm{k}}}\\}$ in terms of
$\\{l^{(i)}_{\bm{\mathrm{k}}}\\}$ via

$\displaystyle
l^{(i+1)}_{\bm{\mathrm{k}}}=\frac{|\bm{\mathrm{v}}_{\bm{\mathrm{k}}}|L^{2}-\hbar
v_{p}\sum_{\bm{\mathrm{k^{\prime}}}}|\bm{\mathrm{q}}|\mathcal{C}_{\bm{\mathrm{k}},\bm{\mathrm{k^{\prime}}}}\mathcal{F}_{\bm{\mathrm{k}},\bm{\mathrm{k^{\prime}}}}^{\mu,T}\cos\theta_{\bm{\mathrm{v}}}l^{(i)}_{\bm{\mathrm{k^{\prime}}}}}{\hbar
v_{p}\sum_{\bm{\mathrm{k^{\prime}}}}|\bm{\mathrm{q}}|\mathcal{C}_{\bm{\mathrm{k}},\bm{\mathrm{k^{\prime}}}}\mathcal{F}_{\bm{\mathrm{k}},\bm{\mathrm{k^{\prime}}}}^{\mu,T}}.$
(30)

In order to optimize for quick convergence, we initialize the procedure using
the explicit formula for an isotropic system with quasi-elastic scattering:

$\displaystyle l^{(0)}_{\bm{\mathrm{k}}}=\left[\frac{\hbar
v_{p}}{|\bm{\mathrm{v}}_{\bm{\mathrm{k}}}|L^{2}}\sum_{\bm{\mathrm{k^{\prime}}}}|\bm{\mathrm{q}}|\mathcal{C}_{\bm{\mathrm{k}},\bm{\mathrm{k^{\prime}}}}\mathcal{F}_{\bm{\mathrm{k}},\bm{\mathrm{k^{\prime}}}}^{\mu,T}\left(1-\cos\theta_{\bm{\mathrm{v}}}\right)\right]^{-1}.$
(31)

In practice, we find very quick convergence and only use two Gauss-Seidel
iterations. We emphasize that our iterative algorithm is a numerical approach
to solving the full BTE, as given in Eqs. (7,12), which is different from yet
equivalent to another commonly-employed technique of “iterating the collision
integral”.

Additionally, to numerically solve Eq. (12) on a discrete momentum grid, we
must broaden the delta functions defining the scattering manifold [see Eq.
(10)]. In practice, we do this by broadening the delta function to a finite-
width step function of a certain small “tolerance”. We then check that our
results are independent of the tolerance variable. We note that our results
are very insensitive to reasonable variation of the tolerance. We also
emphasize that this procedure reproduces the known analytical results for a
single Dirac cone with great accuracy.

##  Appendix C Relaxation time approximation in non-isotropic systems

In the case of elastic scattering and an isotropic band structure, it is well-
known that the solution to the relaxation time approximation to the Boltzmann
equation is also a solution to the full (linearized) Boltzmann equation
Ashcroft and Mermin (1976). In our case, we assume neither isotropy nor
(quasi-)elasticity, which are both present in earlier treatments Hwang and
Sarma (2008); Min et al. (2011); Wu et al. (2019b); Li et al. (2020); Hwang
and Sarma (2019a). In this appendix, we discuss the extent to which the
relaxation time approach holds for our systems.

The canonical “relaxation time approximation” to the Boltzmann equation is the
replacement of the collision integral for the scattering out of state
$\bm{\mathrm{k}}$ with the expression

$\displaystyle\mathcal{I}_{c}[F_{\bm{\mathrm{k}}}]\rightarrow\mathcal{I}^{\text{RT}}_{c}[F_{\bm{\mathrm{k}}}]\equiv\frac{-1}{\tau_{\bm{\mathrm{k}}}}\left[f(\varepsilon_{\bm{\mathrm{k}}})-F_{\bm{\mathrm{k}}}\right],$
(32)

where $F_{\bm{\mathrm{k}}}$ is the full non-equilibrium distribution function
on the set of Bloch states and $f(\varepsilon_{\bm{\mathrm{k}}})$ is the Fermi
distribution function. This introduces the relaxation times as timescales for
the occupation of state $\bm{\mathrm{k}}$ to reach equilibrium.

In the absence of temperature gradients or external magnetic fields, the non-
equilibrium distribution function may be written to linear order in $E$ in
terms of the relaxation times as

$\displaystyle F_{\bm{\mathrm{k}}}\approx
f(\varepsilon_{\bm{\mathrm{k}}})+\frac{1}{T}f(\varepsilon_{\bm{\mathrm{k}}})[1-f(\varepsilon_{\bm{\mathrm{k}}})](e\bm{\mathrm{E}}\cdot\bm{\mathrm{v}}_{\bm{\mathrm{k}}})\tau_{\bm{\mathrm{k}}}\equiv
F^{1}_{\bm{\mathrm{k}}}$ (33)

The distribution function in Eq. (33) is a solution of the Boltzmann equation
under the approximation Eq. (32) and calculating the current from the
distribution function in Eq. (33) gives Eq. (6) in the main text.

The relaxation time approximation is generally uncontrolled, and the true
collision integral in the Boltzmann equation is

$\displaystyle\mathcal{I}_{c}[F_{\bm{\mathrm{k}}}]$
$\displaystyle=-\sum_{\bm{\mathrm{k^{\prime}}}}\mathcal{W}_{\bm{\mathrm{k}}\rightarrow\bm{\mathrm{k^{\prime}}}}F_{\bm{\mathrm{k}}}[1-F_{\bm{\mathrm{k^{\prime}}}}]-\mathcal{W}_{\bm{\mathrm{k^{\prime}}}\rightarrow\bm{\mathrm{k}}}F_{\bm{\mathrm{k^{\prime}}}}[1-F_{\bm{\mathrm{k}}}].$
(34)

However, if there exist $\\{\tau_{\bm{\mathrm{k}}}\\}$ such that for
$F^{1}_{\bm{\mathrm{k}}}$ given by Eq. (33), we have
$\mathcal{I}_{c}[F^{1}_{\bm{\mathrm{k}}}]=\mathcal{I}^{RT}_{c}[F^{1}_{\bm{\mathrm{k}}}]$
to first order in $\bm{\mathrm{E}}$, then Eq. (33) is in fact a solution to
the full linearized Boltzmann equation.

Evaluating Eq. (34) on the distribution function $F^{1}_{\bm{\mathrm{k}}}$ and
using the principle of detailed balance, one may see that

$\displaystyle\mathcal{I}_{c}[F^{1}_{\bm{\mathrm{k}}}]$
$\displaystyle=\frac{-e\bm{\mathrm{E}}}{T}\cdot\sum_{\bm{\mathrm{k^{\prime}}}}\mathcal{W}_{\bm{\mathrm{k}}\rightarrow\bm{\mathrm{k^{\prime}}}}\left[\begin{aligned}
f(\varepsilon_{\bm{\mathrm{k}}})[1-f(\varepsilon_{\bm{\mathrm{k^{\prime}}}})]\\\
\times\big{(}\tau_{\bm{\mathrm{k}}}\bm{\mathrm{v}}_{\bm{\mathrm{k}}}-\tau_{\bm{\mathrm{k^{\prime}}}}\bm{\mathrm{v}}_{\bm{\mathrm{k^{\prime}}}}\big{)}\end{aligned}\right].$
(35)

Comparing with

$\displaystyle\mathcal{I}^{RT}[F^{1}_{\bm{\mathrm{k}}}]$
$\displaystyle=-e\bm{\mathrm{E}}\cdot\bm{\mathrm{v}}_{\bm{\mathrm{k}}}\frac{1}{T}f(\varepsilon_{\bm{\mathrm{k}}})[1-f(\varepsilon_{\bm{\mathrm{k}}})],$
(36)

we find that Eq. (7) is necessary and sufficient for Eq. (33) to be a solution
to the linearized Boltzmann equation.

##  Appendix D Additional data

In this final appendix, we compile additional data for the temperature and
doping dependencies of the resistivity for BBG and RTG. We provide the zero
displacement field ($\Delta=0$) counterparts to Fig. 7, as well as particle-
doped data complementing Fig. 7.

Figure 14:  We plot resistivity data over temperature for hole-doped and
electron-doped Bernal bilayer systems in the absence of a displacement field
($\Delta=0$). This should be compared with the Bernal bilayer data in Figs. 7
and 16. As with Fig. 7, the leftmost two columns give the results of our full
numerical calculation for the resistivity of the two systems at various doping
levels up to $30K$ and $120K$, respectively. The third column gives the
high-$T$ results in the EP regime. The far-right column denotes the doping
values corresponding to the resistivity curves.

Figure 15:  We plot resistivity data over temperature for hole-doped and
electron-doped rhombohedral trilayer systems in the absence of a displacement
field ($\Delta=0$). This should be compared with the rhombohedral trilayer
data in Figs. 7 and 16. As with Fig. 7, the leftmost two columns give the
results of our full numerical calculation for the resistivity of the two
systems at various doping levels up to $30K$ and $120K$, respectively. The
third column gives the high-$T$ results in the EP regime. The far-right column
denotes the doping values corresponding to the resistivity curves.

Figure 16:  We plot resistivity data over temperature for electron-doped
samples, complementing the hole-doped data presented in Fig. 7. We emphasize
that the results are analogous to the hole-doped side. As with Fig. 7, the
leftmost two columns give the results of our full numerical calculation for
the resistivity of the two systems at various doping levels up to $30K$ and
$120K$, respectively. The third column gives the high-T results in the EP
regime. The far-right column denotes the doping values corresponding to the
resistivity curves.

## References

  * Geim and Grigorieva (2013) A. K. Geim and I. V. Grigorieva, Nature 499, 419 (2013), URL https://doi.org/10.1038%2Fnature12385.
  * Novoselov et al. (2006) K. S. Novoselov, E. McCann, S. V. Morozov, V. I. Fal’ko, M. I. Katsnelson, U. Zeitler, D. Jiang, F. Schedin, and A. K. Geim, Nature Physics 2, 177 (2006), URL https://doi.org/10.1038%2Fnphys245.
  * Bistritzer and MacDonald (2011) R. Bistritzer and A. H. MacDonald, Proceedings of the National Academy of Sciences 108, 12233 (2011), URL https://doi.org/10.1073%2Fpnas.1108174108.
  * Morell et al. (2010) E. S. Morell, J. D. Correa, P. Vargas, M. Pacheco, and Z. Barticevic, Physical Review B 82 (2010), URL https://doi.org/10.1103%2Fphysrevb.82.121407.
  * Li et al. (2019) X. Li, F. Wu, and A. H. MacDonald, _Electronic structure of single-twist trilayer graphene_ (2019), URL https://arxiv.org/abs/1907.12338.
  * Kim et al. (2017) K. Kim, A. DaSilva, S. Huang, B. Fallahazad, S. Larentis, T. Taniguchi, K. Watanabe, B. J. LeRoy, A. H. MacDonald, and E. Tutuc, Proceedings of the National Academy of Sciences 114, 3364 (2017), URL https://doi.org/10.1073%2Fpnas.1620140114.
  * Cao et al. (2018a) Y. Cao, V. Fatemi, A. Demir, S. Fang, S. L. Tomarken, J. Y. Luo, J. D. Sanchez-Yamagishi, K. Watanabe, T. Taniguchi, E. Kaxiras, et al., Nature 556, 80 (2018a), URL https://doi.org/10.1038%2Fnature26154.
  * Cao et al. (2018b) Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Nature 556, 43 (2018b), URL https://doi.org/10.1038%2Fnature26160.
  * Cao et al. (2020a) Y. Cao, D. Chowdhury, D. Rodan-Legrain, O. Rubies-Bigorda, K. Watanabe, T. Taniguchi, T. Senthil, and P. Jarillo-Herrero, Phys. Rev. Lett. 124, 076801 (2020a), URL https://link.aps.org/doi/10.1103/PhysRevLett.124.076801.
  * Cao et al. (2021) Y. Cao, J. M. Park, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Nature 595, 526 (2021).
  * Yankowitz et al. (2019) M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean, Science 363, 1059 (2019), URL https://doi.org/10.1126%2Fscience.aav1910.
  * Kerelsky et al. (2019) A. Kerelsky, L. J. McGilly, D. M. Kennes, L. Xian, M. Yankowitz, S. Chen, K. Watanabe, T. Taniguchi, J. Hone, C. Dean, et al., Nature 572, 95 (2019), URL https://doi.org/10.1038%2Fs41586-019-1431-9.
  * Lu et al. (2019) X. Lu, P. Stepanov, W. Yang, M. Xie, M. A. Aamir, I. Das, C. Urgell, K. Watanabe, T. Taniguchi, G. Zhang, et al., Nature 574, 653 (2019), URL https://doi.org/10.1038%2Fs41586-019-1695-0.
  * Stepanov et al. (2020) P. Stepanov, I. Das, X. Lu, A. Fahimniya, K. Watanabe, T. Taniguchi, F. H. Koppens, J. Lischner, L. Levitov, and D. K. Efetov, Nature 583, 375 (2020).
  * Sharpe et al. (2019) A. L. Sharpe, E. J. Fox, A. W. Barnard, J. Finney, K. Watanabe, T. Taniguchi, M. A. Kastner, and D. Goldhaber-Gordon, Science 365, 605 (2019), URL https://doi.org/10.1126%2Fscience.aaw3780.
  * Chen et al. (2020) G. Chen, A. L. Sharpe, E. J. Fox, Y.-H. Zhang, S. Wang, L. Jiang, B. Lyu, H. Li, K. Watanabe, T. Taniguchi, et al., Nature 579, 56 (2020), URL https://doi.org/10.1038%2Fs41586-020-2049-7.
  * Rozen et al. (2021) A. Rozen, J. M. Park, U. Zondiner, Y. Cao, D. Rodan-Legrain, T. Taniguchi, K. Watanabe, Y. Oreg, A. Stern, E. Berg, et al., Nature 592, 214 (2021).
  * Zhou et al. (2022) H. Zhou, L. Holleis, Y. Saito, L. Cohen, W. Huynh, C. L. Patterson, F. Yang, T. Taniguchi, K. Watanabe, and A. F. Young, Science 375, 774 (2022), URL https://doi.org/10.1126%2Fscience.abm8386.
  * Zhou et al. (2021a) H. Zhou, T. Xie, T. Taniguchi, K. Watanabe, and A. F. Young, Nature 598, 434 (2021a), URL https://doi.org/10.1038%2Fs41586-021-03926-0.
  * Zhou et al. (2021b) H. Zhou, T. Xie, A. Ghazaryan, T. Holder, J. R. Ehrets, E. M. Spanton, T. Taniguchi, K. Watanabe, E. Berg, M. Serbyn, et al., Nature 598, 429 (2021b), URL https://doi.org/10.1038%2Fs41586-021-03938-w.
  * Serlin et al. (2020) M. Serlin, C. L. Tschirhart, H. Polshyn, Y. Zhang, J. Zhu, K. Watanabe, T. Taniguchi, L. Balents, and A. F. Young, Science 367, 900 (2020), URL https://doi.org/10.1126%2Fscience.aay5533.
  * Wu et al. (2018) F. Wu, A. MacDonald, and I. Martin, Physical Review Letters 121 (2018), URL https://doi.org/10.1103%2Fphysrevlett.121.257001.
  * Wu et al. (2019a) F. Wu, T. Lovorn, E. Tutuc, I. Martin, and A. H. MacDonald, Phys. Rev. Lett. 122, 086402 (2019a), URL https://link.aps.org/doi/10.1103/PhysRevLett.122.086402.
  * Tschirhart et al. (2022) C. L. Tschirhart, E. Redekop, L. Li, T. Li, S. Jiang, T. Arp, O. Sheekey, T. Taniguchi, K. Watanabe, K. F. Mak, et al., _Intrinsic spin hall torque in a moire chern magnet_ (2022), URL https://arxiv.org/abs/2205.02823.
  * Polshyn et al. (2020) H. Polshyn, J. Zhu, M. A. Kumar, Y. Zhang, F. Yang, C. L. Tschirhart, M. Serlin, K. Watanabe, T. Taniguchi, A. H. MacDonald, et al., Nature 588, 66 (2020), URL https://doi.org/10.1038%2Fs41586-020-2963-8.
  * Jaoui et al. (2021) A. Jaoui, I. Das, G. Di Battista, J. Díez-Mérida, X. Lu, K. Watanabe, T. Taniguchi, H. Ishizuka, L. Levitov, and D. K. Efetov, _Quantum critical behavior in magic-angle twisted bilayer graphene_ (2021), URL https://arxiv.org/abs/2108.07753.
  * Polshyn et al. (2019a) H. Polshyn, M. Yankowitz, S. Chen, Y. Zhang, K. Watanabe, T. Taniguchi, C. R. Dean, and A. F. Young, Nature 15 (2019a), URL https://doi.org/10.1038/s41567-019-0596-3.
  * Cao et al. (2020b) Y. Cao, D. Chowdhury, D. Rodan-Legrain, O. Rubies-Bigorda, K. Watanabe, T. Taniguchi, T. Senthil, and P. Jarillo-Herrero, Phys. Rev. Lett. 124, 076801 (2020b), URL https://link.aps.org/doi/10.1103/PhysRevLett.124.076801.
  * Sarma and Wu (2022) S. D. Sarma and F. Wu, _Strange metallicity of moiré twisted bilayer graphene_ (2022), URL https://arxiv.org/abs/2201.10270.
  * Zhang et al. (2022) Y. Zhang, R. Polski, A. Thomson, . Lantagne-Hurtubise, C. Lewandowski, H. Zhou, K. Watanabe, T. Taniguchi, J. Alicea, and S. Nadj-Perge, _Spin-orbit enhanced superconductivity in bernal bilayer graphene_ (2022), URL https://arxiv.org/abs/2205.05087.
  * Polski et al. (2022) R. Polski, Y. Zhang, Y. Peng, H. S. Arora, Y. Choi, H. Kim, K. Watanabe, T. Taniguchi, G. Refael, F. von Oppen, et al., _Hierarchy of symmetry breaking correlated phases in twisted bilayer graphene_ (2022), URL https://arxiv.org/abs/2205.05225.
  * Arora et al. (2020) H. S. Arora, R. Polski, Y. Zhang, A. Thomson, Y. Choi, H. Kim, Z. Lin, I. Z. Wilson, X. Xu, J.-H. Chu, et al., Nature 583, 379 (2020), URL https://doi.org/10.1038%2Fs41586-020-2473-8.
  * Xie and MacDonald (2020) M. Xie and A. MacDonald, Physical Review Letters 124 (2020), URL https://doi.org/10.1103%2Fphysrevlett.124.097601.
  * Andrei and MacDonald (2020) E. Y. Andrei and A. H. MacDonald, Nature Materials 19, 1265 (2020), URL https://doi.org/10.1038%2Fs41563-020-00840-0.
  * Li et al. (2021) T. Li, S. Jiang, L. Li, Y. Zhang, K. Kang, J. Zhu, K. Watanabe, T. Taniguchi, D. Chowdhury, L. Fu, et al., Nature 597, 350 (2021), URL https://doi.org/10.1038%2Fs41586-021-03853-0.
  * Ghiotto et al. (2021) A. Ghiotto, E.-M. Shih, G. S. S. G. Pereira, D. A. Rhodes, B. Kim, J. Zang, A. J. Millis, K. Watanabe, T. Taniguchi, J. C. Hone, et al., Nature 597, 345 (2021), URL https://doi.org/10.1038%2Fs41586-021-03815-6.
  * Pan et al. (2020) H. Pan, F. Wu, and S. D. Sarma, Physical Review B 102 (2020), URL https://doi.org/10.1103%2Fphysrevb.102.201104.
  * Pan and Sarma (2021) H. Pan and S. D. Sarma, Physical Review Letters 127 (2021), URL https://doi.org/10.1103%2Fphysrevlett.127.096802.
  * Morales-Durán et al. (2021) N. Morales-Durán, A. H. MacDonald, and P. Potasz, Physical Review B 103 (2021), URL https://doi.org/10.1103%2Fphysrevb.103.l241110.
  * Ahn and Sarma (2022) S. Ahn and S. D. Sarma, Physical Review B 105 (2022), URL https://doi.org/10.1103%2Fphysrevb.105.115114.
  * Kerelsky et al. (2021) A. Kerelsky, C. Rubio-Verdú, L. Xian, D. M. Kennes, D. Halbertal, N. Finney, L. Song, S. Turkel, L. Wang, K. Watanabe, et al., Proceedings of the National Academy of Sciences 118 (2021).
  * Khalaf et al. (2019) E. Khalaf, A. J. Kruchkov, G. Tarnopolsky, and A. Vishwanath, Phys. Rev. B 100, 085109 (2019), URL https://link.aps.org/doi/10.1103/PhysRevB.100.085109.
  * Chou et al. (2022a) Y.-Z. Chou, F. Wu, J. D. Sau, and S. D. Sarma, Physical Review B 105 (2022a), URL https://doi.org/10.1103%2Fphysrevb.105.l100503.
  * Chou et al. (2021) Y.-Z. Chou, F. Wu, J. D. Sau, and S. D. Sarma, Physical Review Letters 127 (2021), URL https://doi.org/10.1103%2Fphysrevlett.127.187001.
  * Chou et al. (2022b) Y.-Z. Chou, F. Wu, J. D. Sau, and S. D. Sarma, _Acoustic-phonon-mediated superconductivity in moiréless graphene multilayers_ (2022b), URL https://arxiv.org/abs/2204.09811.
  * You and Vishwanath (2021) Y.-Z. You and A. Vishwanath, _Kohn-luttinger superconductivity and inter-valley coherence in rhombohedral trilayer graphene_ (2021), URL https://arxiv.org/abs/2109.04669.
  * Ghazaryan et al. (2021) A. Ghazaryan, T. Holder, M. Serbyn, and E. Berg, Phys. Rev. Lett. 127, 247001 (2021), URL https://link.aps.org/doi/10.1103/PhysRevLett.127.247001.
  * Szabó and Roy (2022) A. L. Szabó and B. Roy, Physical Review B 105 (2022), URL https://doi.org/10.1103%2Fphysrevb.105.l081407.
  * Cea et al. (2022) T. Cea, P. A. Pantaleón, V. o. T. Phong, and F. Guinea, Phys. Rev. B 105, 075432 (2022), URL https://link.aps.org/doi/10.1103/PhysRevB.105.075432.
  * Dong and Levitov (2021) Z. Dong and L. Levitov, _Superconductivity in the vicinity of an isospin-polarized state in a cubic dirac band_ (2021), URL https://arxiv.org/abs/2109.01133.
  * Chatterjee et al. (2021) S. Chatterjee, T. Wang, E. Berg, and M. P. Zaletel, _Inter-valley coherent order and isospin fluctuation mediated superconductivity in rhombohedral trilayer graphene_ (2021), URL https://arxiv.org/abs/2109.00002.
  * Szabó and Roy (2022) A. L. Szabó and B. Roy, Phys. Rev. B 105, L201107 (2022), URL https://link.aps.org/doi/10.1103/PhysRevB.105.L201107.
  * Qin et al. (2022) W. Qin, C. Huang, T. Wolf, N. Wei, I. Blinov, and A. H. MacDonald, arXiv preprint arXiv:2203.09083 (2022).
  * Dai et al. (2022) H. Dai, R. Ma, X. Zhang, and T. Ma, arXiv preprint arXiv:2204.06222 (2022).
  * Hwang and Sarma (2008) E. H. Hwang and S. D. Sarma, Physical Review B 77 (2008), URL https://doi.org/10.1103%2Fphysrevb.77.115449.
  * Min et al. (2011) H. Min, E. H. Hwang, and S. D. Sarma, Physical Review B 83 (2011), URL https://doi.org/10.1103%2Fphysrevb.83.161404.
  * Wu et al. (2019b) F. Wu, E. Hwang, and S. D. Sarma, Physical Review B 99 (2019b), URL https://doi.org/10.1103%2Fphysrevb.99.165112.
  * Li et al. (2020) X. Li, F. Wu, and S. D. Sarma, Physical Review B 101 (2020), URL https://doi.org/10.1103%2Fphysrevb.101.245436.
  * Hwang and Sarma (2019a) E. H. Hwang and S. D. Sarma, Physical Review B 99 (2019a), URL https://doi.org/10.1103%2Fphysrevb.99.085105.
  * Ziman (1960) J. M. Ziman, _Electrons and Phonons_ (Oxford University Press, 1960), ISBN 9780198507796.
  * Ashcroft and Mermin (1976) N. W. Ashcroft and N. D. Mermin, _Solid State Physics_ (Harcourt College Publishers, 1976), ISBN 9780030839931\.
  * Efetov and Kim (2010) D. K. Efetov and P. Kim, Physical Review Letters 105 (2010), URL https://doi.org/10.1103%2Fphysrevlett.105.256805.
  * Jung and MacDonald (2014) J. Jung and A. H. MacDonald, Phys. Rev. B 89, 035405 (2014), URL https://link.aps.org/doi/10.1103/PhysRevB.89.035405.
  * Zhang et al. (2010) F. Zhang, B. Sahu, H. Min, and A. H. MacDonald, Physical Review B 82 (2010), URL https://doi.org/10.1103%2Fphysrevb.82.035409.
  * Polshyn et al. (2019b) H. Polshyn, M. Yankowitz, S. Chen, Y. Zhang, K. Watanabe, T. Taniguchi, C. R. Dean, and A. F. Young, Nature Physics 15, 1011 (2019b).
  * Siriviboon et al. (2021) P. Siriviboon, J.-X. Lin, X. Liu, H. D. Scammell, S. Liu, D. Rhodes, K. Watanabe, T. Taniguchi, J. Hone, M. S. Scheurer, et al., _A new flavor of correlation and superconductivity in small twist-angle trilayer graphene_ (2021), URL https://arxiv.org/abs/2112.07127.
  * Sohier et al. (2014) T. Sohier, M. Calandra, C.-H. Park, N. Bonini, N. Marzari, and F. Mauri, Phys. Rev. B 90, 125414 (2014), URL https://link.aps.org/doi/10.1103/PhysRevB.90.125414.
  * Xie and Foster (2016) H.-Y. Xie and M. S. Foster, Phys. Rev. B 93, 195103 (2016), URL https://link.aps.org/doi/10.1103/PhysRevB.93.195103.
  * Ghahari et al. (2016) F. Ghahari, H.-Y. Xie, T. Taniguchi, K. Watanabe, M. S. Foster, and P. Kim, Phys. Rev. Lett. 116, 136802 (2016), URL https://link.aps.org/doi/10.1103/PhysRevLett.116.136802.
  * Coleman (2015) P. Coleman, _Introduction to Many-Body Physics_ (Cambridge University Press, 2015), ISBN 9780521864886\.
  * Poniatowski et al. (2021) N. R. Poniatowski, T. Sarkar, S. D. Sarma, and R. L. Greene, Physical Review B 103 (2021), URL https://doi.org/10.1103%2Fphysrevb.103.l020501.
  * Sarkar et al. (2018) T. Sarkar, R. L. Greene, and S. D. Sarma, Physical Review B 98 (2018), URL https://doi.org/10.1103%2Fphysrevb.98.224503.
  * Hwang and Sarma (2019b) E. H. Hwang and S. D. Sarma, Physical Review B 99 (2019b), URL https://doi.org/10.1103%2Fphysrevb.99.085105.
  * Emery and Kivelson (1995) V. J. Emery and S. A. Kivelson, Phys. Rev. Lett. 74, 3253 (1995), URL https://link.aps.org/doi/10.1103/PhysRevLett.74.3253.
  * Gurvitch (1981) M. Gurvitch, Phys. Rev. B 24, 7404 (1981), URL https://link.aps.org/doi/10.1103/PhysRevB.24.7404.
  * Millis et al. (1999) A. J. Millis, J. Hu, and S. D. Sarma, Physical Review Letters 82, 2354 (1999), URL https://doi.org/10.1103%2Fphysrevlett.82.2354.
  * Hussey et al. (2004) N. E. Hussey, K. Takenaka, and H. Takagi, Philosophical Magazine 84, 2847 (2004), URL https://doi.org/10.1080%2F14786430410001716944.