---
abstract: 'In this work we carry out a theoretical study of the phonon-induced resistivity in twisted double bilayer graphene (TDBG), in which two Bernal-stacked bilayer graphene devices are rotated relative to each other by a small angle $\theta$. We show that at small twist angles ($\theta\sim 1^\circ$) the effective mass of the TDBG system is greatly enhanced, leading to a drastically increased phonon-induced resistivity in the high-temperature limit where phonon scattering leads to a linearly increasing resistivity with increasing temperature. We also discuss possible implications of our theory on superconductivity in such a system, and provide an order of magnitude estimation of the superconducting transition temperature.'
author:
- Xiao Li
- Fengcheng Wu
- 'S. Das Sarma'
bibliography:
- 'Ref-TDBLG.bib'
title: Phonon scattering induced carrier resistivity in twisted double bilayer graphene
---

Introduction
============

Recent experimental discoveries of correlated insulator and superconductivity in twisted bilayer graphene (TBG) [@Cao2018; @Cao2018a; @Yankowitz2019; @Sharpe2019; @Lu2019] have attracted great interest in the community. The fact that the electronic band structure of TBG can become almost flat near the magic angle [@Bistritzer2011] strongly enhances the effect of interactions, making it possible to study novel quantum phases that are otherwise difficult to realize experimentally. Furthermore, the apparent similarities between the phase diagram of TBG and cuprate high-temperature superconductors [@Lee2006] suggest that the study of electron correlations in TBG may provide useful hints for our understanding of the electronic properties in cuprates.

The experimental observation of novel quantum phases in TBG has since stimulated the investigation of other van der Waals heterostructures using the twist angle degree of freedom, including, e.g., the trilayer graphene/h-BN [[[moiré]{}]{}]{} superlattice [@Chen2019; @Chen2019a]. One of the motivations of such studies is to go beyond certain limitations of TBG. For example, although the twist angle offers an unprecedented tuning knob to modify the electron band structure in TBG, it still cannot be changed continuously. To date, properties of TBG are mostly modified by fabricating new devices with different twist angles or by applying hydrostatic pressure [@Carr2018; @Yankowitz2019]. It will thus be advantageous to find a way to modify the band structure of a van der Waals heterostructure continuously near quantum critical points, which will enable a more detailed experimental characterization of the electron correlation effects.

Twisted double bilayer graphene (TDBG) has emerged as a promising platform in this respect, because the band structure of a single Bernal-stacked bilayer graphene [@Neto2009; @McCann2013] can be tuned continuously by an external perpendicular electric field. Consequently, one can expect to adjust the band structure of TDBG continuously by an external electric field. Such a tunability is highly desirable, especially near certain quantum critical points. As a result, TDBG has attracted much attention and rapid experimental [@Shen2019; @Cao2019; @Liu2019] and theoretical [@Lee2019] progress has been made. In particular, the application of an external electric field has indeed given rise to a very rich TDBG phase diagram, including signatures of correlated insulator states as well as superconductivity and possibly ferromagnetism in some cases.

The interesting TDBG physics for small twist angles arises from the same moire flatband physics dominating the extensively studied TBG phenomena near the magic twist angle. Basically, the moire potential for small twist angle strongly flattens the relevant graphene bands, leading to very small band velocities (or very large carrier effective masses), which lead to a great enhancement of all interaction phenomena since typically interaction physics is proportional to the carrier effective mass. In the current work, we use a suitable continuum model TDBG band structure to estimate band flattening effects.

In this work, however, we take a different perspective and study the resistivity in TDBG in the high-temperature limit, when phonon scattering will be the dominant mechanism for resistivity. Thus, instead of focusing on the $T=0$ ground state phase diagram, we investigate the ohmic transport properties of the finite temperature effective metallic TDBG phase above the applicable critical temperatures (or the ground state energy gaps) of the symmetry-broken states where TDBG behaves as a metal. Considering rather clean systems, and focusing specifically on the temperature dependence of carrier resistivity, we neglect effects of disorder, impurities, and defects since the main temperature dependence of metallic resistivity arises from phonon scattering effects. We also ignore all electron-electron interaction effects, and only take into account resistive scattering by acoustic scatttering.

Specifically, the scattering of electrons by acoustic phonons can be generally divided into two regimes: a low-temperature regime ($T<{T_\text{BG}}$) and a high-temperature one ($T>{T_\text{BG}}$). The characteristic temperature ${T_\text{BG}}$ is known as the [[Bloch-Grüneisen]{}]{} (BG) temperature, given by $k_B{T_\text{BG}}= 2\hbar {v_\text{ph}}k_F$ [@Hwang2008; @Min2011], where $k_B$ is the Boltzmann constant, ${v_\text{ph}}$ is the phonon velocity and $k_F$ is the Fermi wave vector. (We note that for regular metals where ${T_\text{BG}}$ is very high, or in any situation where ${T_\text{BG}}>T_D$ with $T_D$ being the Debye temperature, the characteristic temperature defining the low and high temperature phonon scattering regimes is $T_D$ and not ${T_\text{BG}}$.) In the high-temperature regime, the electron resistivity will scale as a linear function of temperature $T$, giving rise to $\rho\sim T$, which has been well understood in the context of graphene devices [@Neto2009; @Peres2010; @Sarma2011]. We are interested in this regime because similar to TBG, a wide range of linear-in-$T$ resistivity has been observed in this temperature range in TDBG [@Shen2019; @Cao2019; @Liu2019]. In the context of TBG, such a behavior is often attributed to the putative ‘strange-metal’ phase [@Cao2019a; @Polshyn2019], although it can be compatible with a phonon-scattering mechanism [@Wu2019], [[albeit]{}]{} with greatly enhanced phonon scattering induced carrier resistivity. In this work, we will theoretically study the phonon-induced resistivity in TDBG in the high-temperature regime, and analyze its compatibility with the experimental observations. In particular, we want to understand whether electron-phonon scattering in TDBG can be a contributing factor for the linear-in-$T$ resistivity seen in recent experiments. Our work can be thought of as the TDBG generalization of Ref. [@Wu2019] or as the small twist angle double-bilayer generalization of Ref. [@Min2011]. The goal is to theoretically obtain the acoustic phonon scattering induced carrier resistivity of TDBG as a function of temperature, twist angle, and carrier density.

The structure of the paper is the following. In Section \[Section:Model\] we set up a continuum model for TDBG and demonstrate that its low-energy bands become almost flat (i.e. very large effective mass or equivalently very small effective velocity) at small twist angles ($\theta\sim 1^\circ$). In Section \[Section:Resistivity\] we explain the theoretical framework we use to evaluate phonon-induced resistivity in TDBG and present our numerical results. In Section \[Section:Discussions\] we provide some additional discussions. In particular, we will comment on the possible implications of our theory on superconductivity in TDBG, and provide a rough estimate of the superconducting transition temperature $T_c$ arising from the enhanced electron-phonon coupling. Finally, in Section \[Section:Conclusions\] we provide a brief summary of our results.

Continuum description of a twisted double bilayer graphene \[Section:Model\]
============================================================================

We start by introducing the continuum model of a TDBG. We consider two Bernal-stacked bilayer graphene (BLG) rotated relative to each other by a small angle $\theta$, as shown in Fig. \[Fig:Bandstructure\](a). In particular, we adopt the convention that the top BLG will be rotated by an angle of $\theta/2$, while the bottom one will be rotated by $-\theta/2$. As a result, the continuum description of TDBG near valley $+K$ can be written as $$\begin{aligned}
    \mathcal{H}_{+} = {\begin{pmatrix} h_t({\bm{k}}) & T(\bm{r}) \\ T^\dagger(\bm{r}) & h_b({\bm{k}}) \end{pmatrix}}, \label{Eq:TheModel}\end{aligned}$$ which is given in the basis of {$A_1$, $B_1$, $A_2$, $B_2$, $A_3$, $B_3$, $A_4$, $B_4$}. Here $A$ and $B$ denote the two sublattices of BLG and indices $1$-$4$ denote the four atomic layers, with $1$-$2$ belonging to the top BLG and $3$-$4$ to the bottom one.

In the continuum model Eq. , $h_{t(b)}$ denotes the Hamiltonian for the isolated top (bottom) BLG, given by $$\begin{aligned}
    h_{\lambda}({\bm{k}}) = 
    \begin{pmatrix}
        0 & \hbar vk_{\lambda}^\ast e^{il_{\lambda}\theta/2} & 0 & 0  \\
        \hbar vk_{\lambda}e^{-il_{\lambda}\theta/2} & 0 & \gamma_1 & 0\\
        0 & \gamma_1 & 0 & \hbar vk_{\lambda}^\ast e^{il_{\lambda}\theta/2}\\
        0 & 0 & \hbar vk_{\lambda}e^{-il_{\lambda}\theta/2} & 0 
    \end{pmatrix}. \label{Eq:BareBLG}\end{aligned}$$ In the above equation $\lambda = t,b$ denotes the top and bottom BLG, and $l_{t(b)}=+1 (-1)$. In addition, $k_{t(b)}\equiv k_x^{(t(b))} + ik_y^{(t(b))}$ denotes the (complex) in-plane momentum measured from the Brillouin zone corner of the top (bottom) BLG. In addition, $v = \SI{1e6}{m/s}$ is the bare Dirac velocity of monolayer graphene, while $\gamma_1$ is the interlayer coupling energy of an isolated BLG. Note that the value of $\gamma_1$ in the literature varies widely from $\SI{300}{meV}$ to $\SI{400}{meV}$ [@McCann2013]. In this work we will take $\gamma_1 = \SI{380}{meV}$, but note that results do depend on the specific choice of the $\gamma_1$ band parameter.

In TDBG the [[[moiré]{}]{}]{} potential arising from the twist angle between the two BLGs only induces direct coupling between atomic layers $2$ and $3$. As a result, the [[[moiré]{}]{}]{} potential term in the continuum model Eq.  can be written as $$\begin{aligned}
    T(\bm{r}) = 
    \begin{pmatrix}
        0 & 0 \\ \mathbbm{t}(\bm{r}) & 0 
    \end{pmatrix}, \end{aligned}$$ where $\mathbbm{t}(\bm{r}) = w \displaystyle\sum_{j=1}^3 T_j e^{i\bm{Q}_j\cdot \bm{r}}$. Here $w\simeq\SI{118}{meV}$ [@Wu2018] and $$\begin{aligned}
    T_j = \sigma_0 + \cos(2\pi j/3)\sigma_x + \sin(2\pi j/3)\sigma_y, \; (j = 1, 2, 3), \end{aligned}$$ where $\sigma_i$ are the Pauli matrices. The three vectors $\bm{Q}_j$ read as $$\begin{aligned}
    \bm{Q}_1 = K_\theta \left(\dfrac{\sqrt{3}}{2}, \dfrac{1}{2}\right), \; 
    \bm{Q}_2 = K_\theta \left(-\dfrac{\sqrt{3}}{2}, \dfrac{1}{2}\right), \; 
    \bm{Q}_3 = K_\theta (0, -1), \notag \end{aligned}$$ with $K_\theta = 4\pi/(3a_M)$. Here $a_M = a_0/[2\sin(\theta/2)]$ is the lattice constant of TDBG, and $a_0$ is the lattice constant of monolayer graphene.

Within this continuum description, the band structure of TDBG can be obtained by diagonalizing a large matrix in the momentum space which connects each ${\bm{k}}$ point in the first superlattice [[[moiré]{}]{}]{} Brillouin zone (MBZ) to three other points ${\bm{k}}+ \bm{Q}_j$ ($j = 1,2,3$) in an adjacent MBZ. The resulting band structures for TDBG at three different twist angles are shown in Fig. \[Fig:Bandstructure\](b). For a large twist angle ($\theta=5.0^\circ$), the band structure near the $\kappa_{-}$ point in the MBZ is close to that of pristine BLG, although a small band gap opens up due to the broken inversion symmetry. In contrast, for much smaller twist angles ($\theta = 1.5^\circ$ and $\theta = 1.31^\circ$), the bands near the $\kappa_{-}$ point become very flat, giving rise to a large density of states (DOS) near the band bottom. In Fig. \[Fig:Bandstructure\](c)-(d) we show the DOS per spin per valley $\nu(\varepsilon_F)$ in TDBG, which indeed becomes quite large for small twist angles ($\theta\sim1^\circ$). In addition, the flattened bands also lead to a much reduced Fermi velocity, as shown in Fig. \[Fig:FermiVelocity\]. This physics strongly enhances phonon scattering as we discuss later in the paper.

Approximate zero-energy eigenstates
-----------------------------------

In order to obtain the full band structure of TDBG one must numerically diagonalize a large matrix. However, near the $\kappa_{\pm}$ points in the MBZ it is possible to obtain an approximate analytical expression for the two lowest-energy eigenstates. Such a method was first developed in Ref. [@Bistritzer2011] to obtain an approximate two-band model for TBG, and was used to obtain approximate lowest-energy eigenstates in TDBG in Ref. [@Lee2019].

Specifically, one can truncate the Hamiltonian $\mathcal{H}_+$ in Eq.  by retaining only four momentum points ${\bm{k}}_t$, and ${\bm{k}}_b^{(j)} \equiv {\bm{k}}_t + \bm{Q}_j$ ($j = 1, 2, 3$), and obtain a $8$-by-$8$ Hamiltonian $H(k)$ $$\begin{aligned}
    H(k) = 
    \begin{pmatrix}
        {\tilde{h}}_0(k) & T_1(k) & T_2(k) & T_3(k) \\
        T_1^{\dagger}(k) & {\tilde{h}}_1(k) & 0 & 0 \\
        T_2^{\dagger}(k) & 0 & {\tilde{h}}_2(k) & 0 \\
        T_3^{\dagger}(k) & 0 & 0 & {\tilde{h}}_3(k) \\
    \end{pmatrix}, \label{Eq:8-bandModel}\end{aligned}$$ where [@Lee2019] $$\begin{aligned}
    {\tilde{h}}_0(k) &=
    \begin{pmatrix}
        0 & -v^2(k_\theta^\ast)^2/\gamma_1 \\
        -v^2(k_\theta)^2/\gamma_1 & 0
    \end{pmatrix}, \notag\\
    {\tilde{h}}_j(k) &= 
    \begin{pmatrix}
        0 & -v^2\left[(k+Q_j)_{-\theta}^\ast\right]^2/\gamma_1 \\
        -v^2\left[(k+Q_j)_{-\theta}\right]^2/\gamma_1 & 0
    \end{pmatrix}, \notag\\
    T_j(k) &= 
    \begin{pmatrix}
        -t_Mvk_\theta^\ast/\gamma_1 & t_Mv^2\lambda_j^\ast k_\theta^\ast(k+Q_j)_{-\theta}^\ast/\gamma_1^2\\
        t_m\lambda_j & -t_Mv(k+Q_j)_{-\theta}^\ast/\gamma_1
    \end{pmatrix}. \end{aligned}$$ In the above results, we have defined $\lambda_j = e^{2\pi ij/3}$ ($j = 1, 2, 3$). In addition, we have introduced a short-hand notation that $k_\theta \equiv (k_x + ik_y)e^{-i\theta/2}$.

One can verify that in the $k\to0$ limit the two zero-energy eigenstates of $H(k)$ in Eq.  can be approximately written as $$\begin{aligned}
    \ket{\Psi^{(\alpha)}} = S_\alpha 
    \begin{pmatrix}
        \psi_0^{(\alpha)} \\ \psi_{1}^{(\alpha)}\\ \psi_{2}^{(\alpha)} \\ \psi_3^{(\alpha)}
    \end{pmatrix}, 
    \quad \alpha = A, B, \label{Eq:ZeroEnergyState}\end{aligned}$$ where $S_\alpha$ is the normalization factor, $\psi_0^{(\alpha)}$ is a two-component spinor, and $$\begin{aligned}
    \psi_j^{(\alpha)} = \dfrac{t_M}{vK_\theta^4}{\begin{pmatrix} 0 & K_\theta^2(Q_j)_\theta \\ 0 & \dfrac{\gamma_1}{v}e^{ij\phi}\left[(Q_j)_\theta^2\right]^\ast \end{pmatrix}} \psi_0^{(\alpha)} \equiv M_j \psi_0^{(\alpha)}, \quad j = 1, 2, 3. \notag \end{aligned}$$ The wave function normalization can be determined from the following condition, $$\begin{aligned}
    1 &= \inner{\Psi^{(\alpha)}}{\Psi^{(\alpha)}} 
    = \lvert S_\alpha\rvert^2 \left(\psi_0^{(\alpha)}\right)^\dagger \left[1 + \sum_{j=1}^{3} M^\dagger_j M_j\right]\psi_0^{(\alpha)} \notag\\
    &= \lvert S_\alpha\rvert^2 \left(\psi_0^{(\alpha)}\right)^\dagger {\begin{pmatrix} 1 & 0 \\ 0 & 1+3\Delta \end{pmatrix}} \psi_0^{(\alpha)}, \end{aligned}$$ where $\Delta = \dfrac{\gamma_1^2t_M^2}{v^4K_\theta^4}\left(1+\dfrac{v^2K_\theta^2}{\gamma_1^2}\right)$. Now if we adopt the natural choice of $\psi_0^{(A)} = \begin{pmatrix} 1\\0\end{pmatrix}$, and $\psi_0^{(B)} = \begin{pmatrix} 0\\ 1 \end{pmatrix}$, we find that $$\begin{aligned}
    S_A = 1, \quad S_B = \dfrac{1}{\sqrt{1+3\Delta}}. \label{Eq:Normalization}\end{aligned}$$ This approximate analytical expression for the lowest-energy states in TDBG will be helpful for our analysis of phonon-induced resistivity in the next section.

Phonon-induced resistivity \[Section:Resistivity\]
==================================================

In the previous section we have shown that for small twist angles ($\theta\sim1^\circ$), the Fermi velocity of TDBG near the MBZ corners can become quite small. As we will show in this section, such a substantial reduction in Fermi velocity can give rise to a much enhanced phonon-induced resistivity in the high-temperature ($T\gg{T_\text{BG}}$) limit as happens also for TBG at small twist angles [@Wu2019]. In particular, in such a limit the resistivity scales linearly with temperature, $\rho\approx CT$, and our goal is to estimate the coefficient $C$ and explain how it increases substantially at small twist angles. To verify the validity of our theory, we also numerically evaluate the resistivity in the full temperature range (i.e. $T\gg{T_\text{BG}}$ as well as $T\ll{T_\text{BG}}$), and estimate the crossover temperature above which this enhanced linear-in-$T$ resistivity regime applies.

Resistivity from Boltzmann transport theory
-------------------------------------------

To begin with, we recall that in the Boltzmann transport theory the energy-averaged scattering time $\tau$ in monolayer graphene in the limit of $k_BT\ll \varepsilon_F$ is given by [@Min2011] $$\begin{aligned}
    {\langle\tau\rangle}^{-1} = \dfrac{2\pi}{\hbar}\nu_0\lvert W(k_F)\rvert^2 I, \label{Eq:ScatteringTime}\end{aligned}$$ where $\nu_0$ is the DOS per spin and valley at the Fermi energy, and $\lvert W(k_F)\rvert^2 = D^2\hbar k_F/(2\rho_m {v_\text{ph}})$ is the squared matrix element for acoustic phonon scattering. Here $\rho_m=\SI{7.6e-8}{g/cm^2}$ is the mass density of a single graphene sheet, ${v_\text{ph}}=\SI{2.6e6}{cm/s}$ is the phonon velocity in monolayer graphene, $D=\SI{25}{eV}$ is the acoustic phonon deformation potential [@Min2011], and $v_F$ is the Fermi velocity. The integral $I$ has the following form, $$\begin{aligned}
    I = \int \dfrac{\phi}{2\pi} \dfrac{F(q)(1-\cos\phi)}{\epsilon^2(q)} \dfrac{2q}{k_F}\beta \hbar\omega_q N_q(N_q+1), \label{Eq:IntegralI}\end{aligned}$$ where $q = 2k_F \sin(\phi/2)$ is the magnitude of the acoustic phonon wave vector, $\beta = 1/(k_B T)$, and $F(q)$ is the chiral factor defined as the square of the wave function overlap between incoming and scattered electrons. In addition, $N_q = (e^{\beta\hbar\omega_q}-1)^{-1}$ is the phonon occupation number, with $\omega_q = {v_\text{ph}}q$ being the frequency of the acoustic phonon. Finally, $\epsilon(q)$ is the dielectric function, which takes into account the screening effect at wave vector $q$. In this work we will only consider the unscreened limit, so we will take $\epsilon(q) = 1$. The reason for neglecting screening, which is easy to include, is that there is no experimental evidence that the electron-acoustic phonon resistive scattering gets screened in graphene as a direct comparison between theory [@Hwang2008] and experiment [@Efetov2010] supports the unscreened approximation. We thus do not believe that screening plays any role in TDBG (or TBG) phonon scattering.

We can apply the above formalism to the case of TDBG, and obtain the electron resistivity as $\rho = \sigma^{-1}$, where $\sigma$ is the electron conductivity, given by $$\begin{aligned}
  \sigma = g_s g_v e^2\nu(\varepsilon_F)\dfrac{v_F^2}{2}{\langle\tau\rangle}. \label{Eq:Conductivity}\end{aligned}$$ In the above equation $g_s = 2$ and $g_v = 2$ are the degeneracies due to electron spin and valley degrees of freedom, respectively, while $\nu(\varepsilon_F)$ is the DOS per spin per valley in TDBG shown in Fig. \[Fig:Bandstructure\]. It is worth noting that when evaluating the scattering time ${\langle\tau\rangle}$ in TDBG using Eq. , we should replace the DOS $\nu_0$ there by $\nu(\varepsilon_F)/2$, for the following reasons. In this work we only consider electron densities below the van Hove singularity in TDBG, in which case the topology of the Fermi surface consists of two disconnected mini-valleys ($\kappa_{\pm}$) in the MBZ near the $+K$ valley in the original Brillouin zone of BLG. As a result, the scattering matrix element $W(k_F)$ in Eq.  is only appreciable for electrons within the same mini-valley. The above observation leads us to conclude that in the low-density regime we are interested in, only half of electrons at the Fermi surface contribute to the scattering time. Consequently, we need to substitute $\nu_0$ by $\nu(\varepsilon_F)/2$ in Eq.  [^1]. Putting everything together, we finally obtain the following expression for the resistivity in TDBG: $$\begin{aligned}
 \rho =\dfrac{1}{2g_sg_v}\left(\dfrac{h}{e^2}\right)\left(\dfrac{D^2k_F I}{\hbar\rho_m{v_\text{ph}}v_F^2}\right). \label{Eq:Resistivity}\end{aligned}$$

In order to calculate the resistivity in TDBG using the above equation, we need to evaluate the integral $I$, whose explicit form is given in Eq. . It can be simplified by setting $x = q/(2k_F) = \sin(\phi/2)$, which yields $$\begin{aligned}
  I = \dfrac{16}{\pi}\int_{0}^{1} dx \dfrac{F(2k_Fx)}{\sqrt{1-x^2}} \dfrac{z_\text{BG} x^4 e^{z_\text{BG} x}}{(e^{z_\text{BG} x}-1)^2}, \end{aligned}$$ where $z_\text{BG} = {T_\text{BG}}/T$. In the high-temperature limit ($T\gg{T_\text{BG}}$) we are interested in, we find that $I\approx z_{\infty}/z_\text{BG}$, where $$\begin{aligned}
    {z_{\infty}}= \dfrac{16}{\pi}\int_0^{1} dx \dfrac{x^2F(2k_Fx)}{\sqrt{1-x^2}}, \end{aligned}$$ and therefore the resistivity in Eq.  becomes $\rho \approx CT$, where the coefficient $C$ is given by $$\begin{aligned}
    C = \dfrac{\pi D^2k_Bz_{\infty}}{2g_s g_v e^2\hbar\rho_m{v_\text{ph}}^2v_F^2}. \label{Eq:Coefficient_C}\end{aligned}$$ Therefore, phonon-induced electron resistivity becomes linear in $T$ in the high-temperature ($T\gg{T_\text{BG}}$) limit, a regime we focus on in this work. In addition, from the above result one can see that the quantity $z_\infty$ is a key quantity in this calculation, which depends solely on the chiral form factor $F(q)$ \[or equivalently, $F(\phi)$\]. Thus, we will discuss this quantity first.

The chiral form factor
----------------------

We will use three different approximations to evaluate the chiral form factor $F(\phi)$ and hence ${z_{\infty}}$ for low-energy conduction-band states in TDBG. Specifically, we will use the two-band and four-band description for a pristine BLG, as well as a low-energy two-band description for TDBG. We will see that they capture different aspects of the band structure. A more accurate estimate of $F(\phi)$ necessitates a full numerical evaluation, which we leave for future studies. We comment that, given the simplified nature of our TDBG band structure model, it is unclear that a full numerical calculation of the form factor is warranted.

### Two-band model for bilayer graphene

We start with the simplest case, where a pristine BLG is described by a two-band model, given by $$\begin{aligned}
    H_\text{BLG-2band} = -\dfrac{\hbar^2v^2}{\gamma_1}{\begin{pmatrix} 0 & (k_x+ik_y)^2 \\ (k_x-ik_y)^2 & 0 \end{pmatrix}}. \end{aligned}$$ As a result, the conduction band eigenstates are given by $$\begin{aligned}
    \ket{\psi_{+}({\phi_{{\bm{k}}}})} = \dfrac{1}{\sqrt{2}}\begin{pmatrix}e^{-2i{\phi_{{\bm{k}}}}}\\1\end{pmatrix}, \end{aligned}$$ and the chiral form factor is given by $$\begin{aligned}
    F(\phi) = \lvert\inner{\psi_{+}({\phi_{{\bm{k}}}}+\phi)}{\psi_{+}({\phi_{{\bm{k}}}})}\rvert^2. \label{Eq:2BandFormFactor}\end{aligned}$$ In order to evaluate ${z_{\infty}}$, we note that $x = \sin(\phi/2)$, and thus $F(\phi) = (1-2x^2)^2$. It follows that $$\begin{aligned}
    {z_{\infty}}^\text{(BLG-2band)} = \dfrac{16}{\pi}\int_0^1 dx \dfrac{x^2}{\sqrt{1-x^2}}(1-2x^2)^2 = 2. \label{Eq:zinf-BLG-2band}\end{aligned}$$ As a result, within the two-band model of BLG, ${z_{\infty}}$ is a constant, independent of either electron density or twist angle.

### Four-band model for bilayer graphene

Next, we consider pristine BLG in the four-band description. The corresponding Hamiltonian is given by $$\begin{aligned}
    H_\text{BLG-4band} = 
    \begin{pmatrix}
        0 & \hbar vke^{-i{\phi_{{\bm{k}}}}} & 0 & 0\\
        \hbar vke^{i{\phi_{{\bm{k}}}}} & 0 & \gamma_1 & 0\\
        0 & \gamma_1 & 0 & \hbar vke^{-i{\phi_{{\bm{k}}}}}\\
        0 & 0 & \hbar vke^{i{\phi_{{\bm{k}}}}} & 0
    \end{pmatrix}, \end{aligned}$$ which is written in the {$A_1$, $B_1$, $A_2$, $B_2$} basis. We consider the lower conduction band of BLG, whose energy is $$\begin{aligned}
    E_c({\bm{k}}) = \dfrac{1}{2}\left[\sqrt{4\hbar^2v^2k^2+\gamma_1^2} - \gamma_1\right], \end{aligned}$$ and the corresponding wave function is given by $$\begin{aligned}
    \psi&({\phi_{{\bm{k}}}})^\dagger = \\
    &\begin{bmatrix}
        e^{2i{\phi_{{\bm{k}}}}}\dfrac{-\sqrt{1+\eta}}{2}, & e^{i{\phi_{{\bm{k}}}}}\dfrac{-\sqrt{1-\eta}}{2}, & e^{i{\phi_{{\bm{k}}}}}\dfrac{\sqrt{1-\eta}}{2}, & \dfrac{\sqrt{1+\eta}}{2}
    \end{bmatrix}. \notag\end{aligned}$$ In the above expression, $n$ is the electron density and $\eta = \left(1 + \frac{n}{n_0}\right)^{-1/2}$, with $n_0 = k_0^2/\pi$ and $\hbar vk_0 = \gamma_1/2$. From this wave function, we can obtain the chiral form factor $F(\phi)$ as follows, $$\begin{aligned}
    F(\phi) = \dfrac{1}{4}\left[(1-\eta) + (1+\eta)\cos\phi\right]^2. \end{aligned}$$ Note that in the low-density limit ($\eta\to1$) the above form factor reduces to $\cos^2\phi$, the result derived from the two-band model given in Eq. , as expected.

The expression for ${z_{\infty}}$ derived from the four-band model for pristine BLG has an appreciable electron density dependence. In particular, we find that $$\begin{aligned}
    {z_{\infty}}^\text{(BLG-4band)} = \dfrac{1}{2}(5\eta^2-2\eta+1). \label{Eq:zinf-BLG-4band}\end{aligned}$$

### Two-band model for TDBG

Finally, we consider a low-energy two-band description of TDBG. Because we are only interested in the chiral form factor, we do not need the exact two-band model for TDBG. Instead, we know from symmetry considerations that to leading order the two-band model for TDBG must be of the form $$\begin{aligned}
    H_\text{TDBG} = \mathcal{A} {\begin{pmatrix} 0 & (k_x+ik_y)^2 \\ (k_x-ik_y)^2 & 0 \end{pmatrix}}, \label{Eq:t-DBLG}\end{aligned}$$ where the coefficient $\mathcal{A}$ depends on band structure details, which we do not need. However, we do need explicit expressions for the basis states of this two-band Hamiltonian at $k=0$, which were already given in Eq.  as $\ket{\Psi^{(A)}}$ and $\ket{\Psi^{(B)}}$. With this knowledge, we can write down general expressions for the eigenstates of this two-band model at small ${\bm{k}}$ as follows, $$\begin{aligned}
    \ket{\zeta, {\bm{k}}} &= \dfrac{1}{\sqrt{2}}\left(\ket{\Psi^{(A)}} + \zeta e^{-2i{\phi_{{\bm{k}}}}}\ket{\Psi^{(B)}} \right) \notag\\
    &\equiv 
    \begin{pmatrix} 
        \ket{\zeta, {\bm{k}}}_0, & \ket{\zeta, {\bm{k}}}_1, & \ket{\zeta, {\bm{k}}}_2, &\ket{\zeta, {\bm{k}}}_3
    \end{pmatrix}, \label{Eq:TwobandEigenstate}\end{aligned}$$ where $\zeta = \pm1$ is the band index, and $$\begin{aligned}
    \ket{\zeta, {\bm{k}}}_n = \dfrac{1}{\sqrt{2}}\left(S_A\ket{\psi_n^{(A)}} + S_B\zeta e^{-2i{\phi_{{\bm{k}}}}}\ket{\psi_n^{(B)}} \right), \end{aligned}$$ where $n = 0, 1, 2, 3$. One can verify that $\ket{\zeta = \pm1, {\bm{k}}}$ are indeed the two eigenstates of the two-band Hamiltonian Eq. .

When calculating the chiral form factor, we will consider phonon scattering in the two layers independently. In particular, note that the first component of the four-component eigenstate $\ket{\zeta = \pm1, {\bm{k}}}$ in Eq.  resides in the top BLG, while the other three reside in the bottom one. As a result, the chiral form factor should be evaluated as follows, $$\begin{aligned}
    F(\phi) &= \lvert _0\inner{\zeta', {\bm{k}}'}{\zeta, {\bm{k}}}_0\rvert^2  + \bigg\lvert \sum_{j=1}^{3} {}_j\inner{\zeta', {\bm{k}}'}{\zeta, {\bm{k}}}_j \bigg\rvert^2 \notag\\
    &\equiv |\mathcal{F}_1(\phi)|^2 + |\mathcal{F}_2(\phi)|^2. \end{aligned}$$ The above derivations lead to the following results, $$\begin{aligned}
    \mathcal{F}_1(\phi) = \dfrac{1}{2}\left(1+\dfrac{\zeta\zeta'}{1+3\Delta}e^{2i\phi}\right), \; 
    \mathcal{F}_2(\phi) = \dfrac{\zeta\zeta'}{2}\dfrac{3\Delta}{1+3\Delta}e^{2i\phi}, \end{aligned}$$ which then gives rise to the following form factor for the two-band model of TDBG, $$\begin{aligned}
    F(\phi) = \dfrac{1}{4}\left[1+\dfrac{9\Delta^2+1}{(3\Delta+1)^2} + \dfrac{2}{3\Delta+1}\cos(2\phi) \right]. \label{Eq:TDBG-FormFactor}\end{aligned}$$ Such a chiral form factor yields the following result for ${z_{\infty}}$, $$\begin{aligned}
    {z_{\infty}}^\text{(TDBG)} = 1 + \dfrac{9\Delta^2+1}{(3\Delta + 1)^2}. \label{Eq:zinf-t-dBLG}\end{aligned}$$

Some numerical results for ${z_{\infty}}$ under different approximations are given in Fig. \[Fig:FormFactor\]. One can see that the three different approximations of ${z_{\infty}}$ are of the same order, although they capture different aspects of the band structure. In the rest of the paper, we will use both ${z_{\infty}}^\text{(TDBG)}$ and ${z_{\infty}}^\text{(BLG-4band)}$ to calculate the phonon-induced resistivity. Note that, in order to evaluate ${z_{\infty}}$ for TDBG accurately, one has to resort to full numerical evaluations from the band structure. Although such a calculation is beyond the scope of this work, we expect that the exact value of ${z_{\infty}}$ is still within the same order of magnitude as the ones we used in this work.

Phonon-induced resistivity: High temperature limit
--------------------------------------------------

After explaining the calculations of ${z_{\infty}}$, we are now ready to evaluate the phonon-induced resistivity explicitly. In this subsection we will consider the high-temperature limit first, when the resistivity is a linear function of temperature, and then present results for the full temperature range in the next subsection. Before showing our results, however, we make a few comments on our numerical evaluation of the coefficient $C$ using Eq. . First, we will use both ${z_{\infty}}^\text{(TDBG)}$ and ${z_{\infty}}^\text{(BLG-4band)}$ to approximate ${z_{\infty}}$, and demonstrate how different approximations affect the final value of $C$. Second, the Fermi velocity $v_F$ will be extracted directly from the full numerical band structure of TDBG, instead of just from the two-band effective model in Eq. . Finally, all of our calculations are limited to carrier densities below the van Hove singularities in the band structure, because our theory will break down for higher carrier densities, when the topology of the Fermi surface is different from our assumptions because of the complications arising from van Hove singularities. Thus, our theory is explicitly limited to low carrier densities where the Fermi level stays below the van Hove singularities.

Some numerical results for the coefficient $C$ are given in Fig. \[Fig:DensityDependence\]. In (a)-(b) the results for two different twist angles are shown. One can see that as the twist angle decreases from $5.0^\circ$ to $1.31^\circ$, the coefficient $C$ increases substantially. Such a trend is also apparent in panel (c), which shows how the coefficient $C$ depends on the twist angle at a fixed electron density. We note that such an angular dependence with resistivity increasing strongly with decreasing twist angle is consistent with recent resistivity measurements in TDBG [@Shen2019; @Cao2019; @Liu2019].

In addition, we find from Fig. \[Fig:DensityDependence\] that within our theory the coefficient $C$ has a strong density dependence. This feature in our theory arises from the fact that the coefficient $C$ is inversely proportional to the Fermi velocity, which has a strong density dependence for parabolic bands. When we compare our theory of phonon-induced resistivity in TDBG with the experimental results, we find that we cannot capture the very weak density dependence of $C$ observed in the experiment. In particular, at low carrier densities and small twist angles ($\theta\sim 1^\circ$) experimental results show that the coefficient $C$ has almost no dependence on the carrier density [@Shen2019; @Cao2019; @Liu2019]. Our current understanding of this disprepancy is that the electron correlation may play an important role in this limit, which strongly suppresses the density dependence of $C$, an effect not accounted for in our theory. We expect that for larger twist angles ($\theta\gtrsim 2.0^\circ$), when the band width becomes large, electron-phonon scattering will overcome the electron correlation effects, and become the dominant mechanism for resistivity in the high-temperature limit. Our theory should be more applicable in that regime, and the coefficient $C$ should exhibit a substantial carrier density dependence in the experimental measurements. It will thus be interesting to carry out an experiment to resolve the crossover between these two regimes, which will help us better understand the role of electron correlation in TDBG.

It is also interesting to draw a comparison between TBG and TDBG in this context. In particular, note that such a density dependence in $C$ is absent in the case of TBG, even within the framework of phonon-induced resistivity and at small twist angles [@Wu2019]. The underlying reason is straightforward: the low-energy electronic states in TBG has a linear dispersion. As a result, the Fermi velocity and hence the coefficient $C$ in TBG does not depend on the carrier density. It is worth noting that such an absence of density dependence in $C$ is consistent with the experimental observations in TBG [@Cao2018a]. By contrast, TDBG bands are parabolic, and hence one expects a density dependence in the temperature coefficient of the resistivity.

[We do mention, however, that our calculated TDBG resistivity approximately agrees with recent measurements [@KimPrivate] where we obtain $d\rho/dT$ $\sim$ $\SI{95}{\Omega/K}$ at a twist angle of $\theta\sim 1.24^\circ$ and a carrier density of $\SI{3.0e12}{cm^{-2}}$, compared with the experimental values of $d\rho/dT$ $\sim$ $\SI{75}{\Omega/K}$.]{} But any definitive agreement between our current theory and the mesured TDBG temperature dependent resistivity awaits a careful experimental study of temperature, twist angle, and density dependence of TDBG resistivity, which is unavailable right now.

Phonon-induced resistivity: Full temperature range
--------------------------------------------------

Finally, we evaluate the phonon-induced resistivity in the full temperature range. Such a calculation will not only allow us to present a complete result for the phonon-induced resistivity in TDBG, but also help us determine the temperature range in which the resistivity is linear in $T$.

To begin with, we consider the integral $I$ in the low-temperature limit ($T\ll{T_\text{BG}}$), which can be evaluated by introducing $y = z_\text{BG} x$, yielding $$\begin{aligned}
I \approx \dfrac{16}{\pi z_\text{BG}^4} \int_0^{+\infty} \dfrac{y^4 e^{y}}{(e^y-1)^2} dy = \dfrac{16}{\pi} \times \dfrac{4! \zeta(4)}{z_\text{BG}^4} \propto T^4, \end{aligned}$$ where we have used $F(0)\approx 1$, and $\zeta(s)$ is the Riemann-$\zeta$ function. As a result, in the low-temperature limit the resistivity depends on the temperature as $\rho \propto T^4$. This is the so-called [[Bloch-Grüneisen]{}]{} regime of phonon scattering, which in 3D metals produces a $T^5$ power law for the temperature-dependent resistivity. In addition, the resistivity becomes independent of the chiral form factor $F(\theta)$ in this low-temperature limit. In Fig. \[Fig:FullResistivity\] we show the resistivity in TDBG across the full temperature range for a carrier density of $\SI{3.0e12}{cm^{-2}}$ and two different twist angles. These results are obtained by using the chiral form factor for the two-band model of TDBG in Eq. . One can clearly observe that the resistivity has a $\rho \propto T^4$ behavior in the low-temperature range, while a $\rho\propto T$ behavior in the high-temperature limit.

After obtaining the resistivity in the full temperature range, it is instructive to examine the crossover temperature above which the resistivity becomes linear in temperature. It has been established previously that the linear-in-$T$ behavior already kicks in at a characteristic temperature $T\gtrsim T_L\approx{T_\text{BG}}/4$ [@Hwang2008; @Min2011]. In Fig. \[Fig:TL\] we plot the crossover temperature $T_L$ in TDBG as a function of carrier density for three different twist angles. We find that $T_L$ is below $\SI{11}{K}$ for almost all carrier densities and twist angles we considered. As a result, our analysis of the phonon-induced resistivity in the high temperature limit should be applicable above $\sim \SI{11}{K}$. In fact, recent resistivity measurements in TDBG indeed show a linear-in-$T$ behavior for temperatures between $\SI{10}{K}$ and $\SI{30}{K}$ [@KimPrivate], a temperature range where our theory is applicable. Therefore, our theory will be relevant for the understanding of the linear-in-$T$ resistivity observed in recent experiments in TDBG [@Shen2019; @Cao2019; @Liu2019]. We note that in the low-temperature [[Bloch-Grüneisen]{}]{} regime, the $T^4$ power law in the resistivity may not be easy to discern because of other resistive scattering contributions such as electron-impurity and electron-electron interactions which are neglected in our theory.

Phonon-Mediated Superconductivity\[Section:Discussions\]
========================================================

We now discuss possible implication of our theory on superconductivity in TDBG. The electron-acoustic phonon coupling mediates an effective attractive electron-electron interaction with a strength given by $g_0$ $=$ $D^2/(4 \rho_m {v_\text{ph}}^2) \approx \SI{50}{meV\cdot nm^2}$ [@Wu2019]. The dimensionless electron-phonon coupling constant is determined by $\lambda^* = g_0 \nu(\varepsilon_F)$, where $\nu(\varepsilon_F)$ is the DOS per spin and valley. Because of the narrow bandwidth for small twist angle ($\sim 1^{\circ}$), $\lambda^*$ in TDBG can reach order of $0.25$ given the DOS shown in Fig. \[Fig:Bandstructure\](d). The superconducting transition temperature $T_c$ can be roughly estimated as $k_B T_c = \Lambda \exp(-1/\lambda^*)$ within a BCS-type theory, where $\Lambda$ is a cutoff energy approximately given by the flatband bandwidth ($\sim 5$ meV). Therefore, $T_c$ can be of order $\SI{1}{K}$ from electron-phonon interactions. Moreover, the effective attractive interactions mediated by acoustic phonons have an enlarged SU(2) $\times$ SU(2) symmetry, namely, each valley has its own spin rotational symmetry. Therefore, acoustic phonons mediate both spin singlet and spin triplet pairings [@Wu2019], and can account for the spin triplet superconductivity experimentally identified in TDBG [@Liu2019]. We note that the possibility of the electron-phonon mediated superconductivity at low temperatures ($<\SI{1}{K}$) and the large phonon-induced linear-in-$T$ resistivity at high temperatures ($>\SI{10}{K}$) are closely connected, both arising from the strongly enhanced electron-phonon coupling induced by flatband moire physics, as has already been emphasized in the context of TBG in Ref. [@Wu2019].

Conclusion \[Section:Conclusions\]
==================================

To summarize, in this work we developed a theory to calculate the phonon-induced resistivity in twisted double bilayer graphene in the high-temperature ($T>{T_\text{BG}}$) limit, where the resistivity $\rho$ scales linearly with temperature $T$, $\rho \approx CT$. We present a quantitative analysis of the coefficient $C$ and showed that it increases substantially as the twist angle $\theta$ is reduced. However, since we did not account for electron correlation effects on the resistivity, we expect that our predictions are likely applicable for devices in which the twist angle is relatively large ($\theta\gtrsim 2^\circ$). The main qualitative conclusion of our theory is that for $T>\SI{10}{K}$ or so, TDBG should manifest a very large linear-in-$T$ resistivity arising from phonon scattering at small twist angles. The linear coefficient should manifest a strong density dependence, which is not seen in current experiments for reasons not obvious right now.

Acknowledgment {#acknowledgment .unnumbered}
==============

This work is supported by Microsoft and Laboratory of Physics. X.L. also acknowledges support from City University of Hong Kong (Project No. 9610428).

[^1]: Note that we should still use the full DOS $\nu(\varepsilon_F)$ in Eq. , because all electrons at the Fermi surface contribute to the conductivity.
