Geometric and conventional contribution to superfluid weight in twisted bilayer graphene

TL;DR

This paper calculates the superfluid weight in twisted bilayer graphene, revealing the significant roles of geometric and conventional contributions, and explains how it remains finite despite flat bands, aligning with experimental findings.

Contribution

It introduces a comprehensive calculation of superfluid weight in TBLG considering both geometric and conventional effects, highlighting their relative importance near the magic angle.

Findings

Both contributions exceed expectations for isolated bands.

Geometric contribution dominates at the magic angle.

Superfluid weight can be tuned by chemical potential and twist angle.

Abstract

By tuning the angle between graphene layers to specific "magic angles" the lowest energy bands of twisted bilayer graphene (TBLG) can be made flat. The flat nature of the bands favors the formation of collective ground states and, in particular, TBLG has been shown to support superconductivity. When the energy bands participating in the superconductivity are well-isolated, the superfluid weight scales inversely with the effective mass of such bands. For flat-band systems one would therefore conclude that even if superconducting pairing is present most of the signatures of the superconducting state should be absent. This conclusion is at odds with the experimental observations for TBLG. We calculate the superfluid weight for TBLG taking into account both the conventional contribution and the contribution arising from the quantum geometry of the bands. We find that both contributions are larger than one would expect treating the bands as well-isolated, that at the magic angle the geometric contribution is larger than the conventional one, and that for small deviations away from the magic angle the conventional contribution is larger than the geometric one. Our results show that, despite the flatness of the bands the superfluid weight in TBLG is finite and consistent with experimental observations. We also show how the superfluid weight can be tuned by varying the chemical potential and the twist angle opening the possibility to tune the nature of the superconducting transition between the standard BCS transition and the Berezinskii-Kosterlitz-Thouless transition.

Tables

Table 1. Table S1 : Results obtained from solving the mean-field gap equation.

$\theta {(}^{\circ })$	1.00	1.05	1.10	1.15
$T_c\left(K\right)$	1.2119	1.6305	0.2340	0.0189
${\Delta }_0$	0.4193	0.4346	0.4465	0.4497
${\Delta }_b$	0.2138	0.2118	0.2100	0.2104

Table 2. (a) Solutions at μ = − 0.30 𝜇 0.30 \mu=-0.30 meV

$\theta {(}^{\circ })$	1.00	1.05	1.10	1.15
$T_c\left(K\right)$	1.2119	1.6305	0.2340	0.0189
${\Delta }_0$	0.4193	0.4346	0.4465	0.4497
${\Delta }_b$	0.2138	0.2118	0.2100	0.2104

Table 3. (b) Solutions for θ = 1.05 ∘ 𝜃 superscript 1.05 \theta=1.05^{\circ} .

$\mu \left(meV\right)$	0	-0.10	-0.20	-0.30	-0.40	-0.50
$T_c\left(K\right)$	1.7392	1.7330	1.7043	1.6305	1.4877	1.2455
${\Delta }_0$	0.4346	0.4346	0.4346	0.4346	0.4346	0.4346
${\Delta }_b$	0.2119	0.2119	0.2119	0.2118	0.2118	0.2118

Table 4. (c) Solutions for θ = 1.00 ∘ 𝜃 superscript 1.00 \theta=1.00^{\circ} .

$\mu \left(meV\right)$	0	-0.10	-0.20	-0.30	-0.40	-0.50
$T_c\left(K\right)$	1.2086	1.2097	1.2219	1.2119	1.1561	1.0461
${\Delta }_0$	0.4186	0.4188	0.4191	0.4193	0.4196	0.4197
${\Delta }_b$	0.2138	0.2138	0.2138	0.2138	0.2138	0.2138

Equations

$$H_{\mathbf{K},l}(\mathbf{k})=e^{-il\frac{\theta }{4}{\tau }_z}[\hslash v_F(\mathbf{k}-{\mathbf{\kappa }}_l)\cdot \mathbf{\tau }-\mu {\tau }_0]e^{il\frac{\theta }{4}{\tau }_z},$$

$$H_{TBL,\mathbf{K}}=\left(\begin{array}{cc}{H}_{\mathbf{K},+1}& T(\mathbf{r})\\ T^{†}(\mathbf{r})& H_{\mathbf{K},-1}\end{array}\right),$$

$$\displaystyle H_{BdG}=\left[\begin{array}[]{cc}H_{TBL,{\bf K}}({\bf k})&\hat{\Delta}_{s}\\ \hat{\Delta}_{s}^{\dagger}&-H^{T}_{TBL,{\bf K}^{\prime}}(-{\bf k})\end{array}\right],$$

$$D_{\mu \nu }^s$$

$$D_{\mu \mu }^s$$

$$V_{\pm \mu ij}^d$$

$$V_{\pm \mu ij}^o$$

$$D_{\mu \mu }^{s,\mathrm{conv}}$$

$$D_{\mu \mu }^{s,\mathrm{geom}}$$

$${\chi }_{\mathbf{b}{\mathbf{b}}^{{}^{\prime }}}^{l{l}^{{}^{\prime }}}=\frac{2g_0}{\mathcal{A}}\sum _{\mathbf{q}{n}_1{n}_2}\mathcal{N}(\mathbf{q})[U(\mathbf{q}{)}_{\mathbf{b}l}{]}^{\ast }U(\mathbf{q}{)}_{{\mathbf{b}}^{{}^{\prime }}{l}^{{}^{\prime }}}.$$

$$\chi \mathbf{v}=\mathbf{v},$$

$$H_{BdG}(\mathbf{p})=\left(\begin{array}{cc}{H}_0(\mathbf{p})& \Delta (\mathbf{p})\\ {\Delta }^{†}(\mathbf{p})& -H_0^T(-\mathbf{p})\end{array}\right).$$

$$H_{BdG}(\mathbf{r})=\left(\begin{array}{cc}{H}_0(-i\nabla )& e^{i\frac{\mathbf{k}}{2}\cdot \mathbf{r}}\Delta (-i\nabla )e^{i\frac{\mathbf{k}}{2}\cdot \mathbf{r}}\\ e^{-i\frac{\mathbf{k}}{2}\cdot \mathbf{r}}{\Delta }^{†}(-i\nabla )e^{-i\frac{\mathbf{k}}{2}\cdot \mathbf{r}}& -H_0^T(+i\nabla )\end{array}\right).$$

$$\mathbf{j}=\frac{2e{D}_s}{\hslash }\nabla \varphi ,$$

$$T_{\mathrm{KT}}=\frac{\pi }{2}{D}_s(T_{\mathrm{KT}}).$$

$$\psi (\mathbf{r})=\left(\begin{array}{c}a(\mathbf{p})e^{i\mathbf{p}\cdot \mathbf{r}}{e}^{i\frac{\mathbf{k}}{2}\cdot \mathbf{r}}\\ b(\mathbf{p})e^{i\mathbf{p}\cdot \mathbf{r}}{e}^{-i\frac{\mathbf{k}}{2}\cdot \mathbf{r}}\end{array}\right).$$

$$H_{\mathrm{eff}}(\mathbf{p},\mathbf{k})=\left(\begin{array}{cc}{H}_0(\mathbf{p}+\frac{\mathbf{k}}{2})& \Delta (\mathbf{p})\\ {\Delta }^{†}(\mathbf{p})& -H_0^T(-\mathbf{p}+\frac{\mathbf{k}}{2})\end{array}\right).$$

$$v_{\mu }(\mathbf{p},\mathbf{k})=\frac{\partial H_{\mathrm{eff}}(\mathbf{p},\mathbf{k},\Delta =0)}{\partial p_{\mu }}.$$

$$\mathbf{j}=\frac{1}{L^2}\sum _{\mathbf{p},n}⟨n,\mathbf{p},\mathbf{k}|{\gamma }_z\mathbf{v}(\mathbf{p},\mathbf{k})|n,\mathbf{p},\mathbf{k}⟩n_F(E_n(\mathbf{p},\mathbf{k})),$$

$$\hat{H}=\sum _{\mathbf{p}}{C}_{\mathbf{p},\mathbf{k}}^{†}{H}_{\mathrm{eff}}(\mathbf{p},\mathbf{k})C_{\mathbf{p},\mathbf{k}}={\hat{H}}^0+{\hat{H}}^1,$$

$${\hat{H}}^0=\sum _{\mathbf{p}}{C}_{\mathbf{p},\mathbf{k}}^{†}{H}_{\mathrm{eff}}(\mathbf{p},0)C_{\mathbf{p},\mathbf{k}},$$

$${\hat{H}}^1=\sum _{\mathbf{p}}{C}_{\mathbf{p},\mathbf{k}}^{†}{\gamma }_z\frac{\mathbf{k}}{2}\cdot \mathbf{v}(\mathbf{p},0)C_{\mathbf{p},\mathbf{k}}$$

$$\hat{\mathbf{j}}=\frac{1}{L^2}\sum _{\mathbf{p}}{C}_{\mathbf{p},\mathbf{k}}^{†}{\gamma }_z\mathbf{v}(\mathbf{p},\mathbf{k})C_{\mathbf{p},\mathbf{k}}={\hat{\mathbf{j}}}^0+{\hat{\mathbf{j}}}^1,$$

$${\hat{\mathbf{j}}}^0=\frac{1}{L^2}\sum _{\mathbf{p}}{C}_{\mathbf{p},\mathbf{k}}^{†}{\gamma }_z\mathbf{v}(\mathbf{p},0)C_{\mathbf{p},\mathbf{k}}$$

$$\hat{\mathbf{j}}^{1}=\frac{1}{L^{2}}\sum_{\mathbf{p}}C_{\mathbf{p},\mathbf{k}}^{\dagger}\gamma_{z}\mathbf{k}\cdot\big{[}\nabla_{\mathbf{k}}\mathbf{v}(\mathbf{p},\mathbf{k})\big{]}_{\mathbf{k}=\mathbf{0}}C_{\mathbf{p},\mathbf{k}}.$$

$$j$$

$$T_{\mu \nu }(\mathbf{p})={\gamma }_z\frac{\partial }{\partial p_{\mu }}{v}_{\nu }(\mathbf{p}),$$

$$D_s$$

$$\sum _{\mathbf{p},i}⟨{\psi }_i(\mathbf{p})|{\tau }_z{T}_{\mu \mu }(\mathbf{p})|{\psi }_i(\mathbf{p})⟩n_F(E_i(\mathbf{p}))$$

$$v_{\mu }(\mathbf{p})=\frac{\partial H_{BdG}(\mathbf{p})}{\partial p_{\mu }}$$

$$D_{s}=\frac{\Delta}{4\pi}\bigg{[}\sqrt{1+\frac{\mu^{2}}{\Delta^{2}}}+\frac{\Delta}{|\mu|}\ln\bigg{(}\frac{|\mu|}{\Delta}+\sqrt{1+\frac{\mu^{2}}{\Delta^{2}}}\bigg{)}\bigg{]},$$

Full text

The ability to control accurately the twist angle $\theta$ between two-dimensional crystals forming a van der Waals systems Li et al. (2010); Geim and Grigorieva (2013); Novoselov et al. (2016); Ajayan et al. (2016) has recently emerged as a powerful way to tune the electronic properties of a condensed matter system. The most remarkable example of such tunability has been observed in twisted bilayer graphene (TBLG). For most values of the twist angle between the graphene sheets the systems behave as a normal two-dimensional (2D) semimetal, however for specific “magic angles” Lopes dos Santos et al. (2007); Mele (2010); Morell et al. (2010); Bistritzer and MacDonald (2011); Trambly de Laissardière et al. (2012) the system’s lowest energy bands become almost completely flat and the system may support topological properties Song et al. (2019); Lian et al. (2018); Po et al. (2019); Liu et al. (2019a); Ahn et al. (2019). Quenched kinetic energy in the flat bands increases the importance of interactions and leads to superconductivity and other correlated states Kopnin et al. (2011); Ojajärvi et al. (2018); Wu et al. (2018); Peltonen et al. (2018); Guo et al. (2018); Liu et al. (2018); You and Vishwanath (2019); Zhang and Senthil (2019); Tang et al. (2019); Roy and Juričić (2019); Kozii et al. (2019); Wu (2019); Wu et al. (2019) recently observed in graphene moiré superlattices  Kim et al. (2017); Cao et al. (2018a, b); Yankowitz et al. (2019); Chen et al. (2019); Lu et al. (2019); Shen et al. (2019); Liu et al. (2019b); Cao et al. (2019); Chen et al. (2019). The hallmark signature of the superconducting state is the absence of electrical resistance. For this to happen the superfluid weight $D^{s}_{xx}$ must be nonzero. For an isolated parabolic band at zero temperature $D^{s}_{xx}\propto n/m^{*}$, where $n$ is the electron density, and $m^{*}$ the effective electron’s mass. From this expression one would conclude that the standard signature of superconductivity might be absent for flat bands because one expect $1/m^{*}$ to vanish proportionally to the bandwidth. This is not what happens experimentally in TBLG.

In order to reconcile experimental observations and theory we notice that the above expression for $D^{s}_{xx}$ assumes an isolated band and neglects the interband matrix elements of the current operator. Neither of these assumptions is valid in TBLG. In an isolated band the density of electrons within the band is constant, and therefore when the superconducting transition occurs the chemical potential is renormalized. The superfluid weight depends strongly on the chemical potential and this renormalization is responsible for the appearance $1/m^{*}$ dependence of the intraband (conventional) contribution to the superfluid weight. In a semimetal where both electron- and hole-like bands are present, such as TBLG, the densities in each band are not conserved in the transition separately and the dependence on the chemical potential is weak so that the conventional contribution can be much larger than expected for isolated bands. Moreover, the bandwidth of the low-energy bands, even though very small, is still finite and larger than the superconducting gap. Therefore, the velocity can be large at some points of the Brillouin zone further enhancing the conventional contribution to the weight. On the other hand, the interband matrix elements give rise to the so called geometric part of $D^{s}_{xx}$, which can be large even for completely flat band well isolated from other bands Peotta and Törmä (2015); Liang et al. (2017).

In this work we calculate the superfluid weight of superconducting twisted bilayer graphene (TBLG) taking into account both the conventional and the geometric parts. We assume singlet pairing and use the experimentally measured value of $T_{c}$ to set the value of the coupling constant that enters the mean field gap equation. We obtain the dependence of the superconducting weight on the twist angle and separate the conventional and geometric parts. We find that at one of the “magic angles”, $\theta=1.05^{\circ}$, the geometric contribution is approximately twice as large as the conventional one. However, just off the magic angle the conventional contribution is larger than the geometric one. We also obtain the dependence of the Berezinski-Kosterlitz-Thouless $T_{\mathrm{KT}}$ temperature on $\theta$ and show that its scaling with the chemical potential is different at the magic angle and away from it. Because our calculations take into account the full band structure of TBLG and include both intra- and interband contributions, they can be used for quantitative predictions and they go beyond the models and approximations previously used in deriving bounds for the superfluid weight Hazra et al. (2019); Xie et al. (2019).

To model the TBLG we use the approach described in Ref. Bistritzer and MacDonald, 2011; Wu et al., 2018. The low-energy states of the isolated single layers of graphene are located at the ${\bf K}$ and ${\bf K}^{\prime}=-{\bf K}$ valleys of the Brillouin zone (BZ). Close to ${\bf K}$ the Hamiltonian for each layer $l=\pm 1$ is

[TABLE]

where $v_{F}=10^{6}$ m/s is graphene’s Fermi velocity, $\mu$ is the chemical potential, and ${\tau_{i}}$ ($i=0,1,2,3$) are the $2\times 2$ Pauli matrices in sublattice space. Because of the rotation of each layer by angle $\theta/2$ the Dirac cone position in layer $l$ is shifted to ${\boldsymbol{\kappa}}_{l}$. We choose moiré BZ in which ${\boldsymbol{\kappa}}_{l}$ are located at the corners and refer to the center of this BZ as the $\gamma$ point. This leads to a Hamiltonian for TBLG around the ${\bf K}$ point

[TABLE]

with periodically varying interlayer tunneling terms $T({\bf r})=w[T_{0}+e^{-i{\bf b}_{2}\cdot{\bf r}}T_{+1}+e^{-i({\bf b}_{2}-{\bf b}_{1})\cdot{\bf r}}T_{-1}]$, where $T_{j}=\tau_{0}+\cos(2\pi j/3)\tau_{x}+\sin(2\pi j/3)\tau_{y}$, ${\bf b}_{1}=(\sqrt{3}Q,0)$ and ${\bf b}_{2}=(\sqrt{3}Q/2,3Q/2)$ are reciprocal basis vectors, $Q=\frac{8\pi}{3a_{0}}\sin(\theta/2)$, $a_{0}$ is the lattice constant of graphene and $w=118$ meV Jung et al. (2014); Wu et al. (2018). $H_{{\bf K}^{\prime}}$ is obtained from $H_{{\bf K}}$ via time-reversal.

We leave $d$-wave pairingWu et al. (2018) for the future studies and focus on $s$-wave pairing. In the presence of superconductivity the mean field theory in Nambu space is described by the Bogoliubov-de-Gennes Hamiltonian

[TABLE]

and $\hat{\Delta}_{s}=\Delta\tau_{0}\sum_{{\bf b}}\Delta_{\bf b}e^{i{\bf b}\cdot{\bf r}}$, where $\Delta$ is the overall amplitude of the superconducting gap, and $\Delta_{\bf b}$ is the normalized coefficient of the ${\bf b}=m_{1}{\bf b}_{1}+m_{2}{\bf b}_{2}$ ($m_{i}\in Z$) Fourier component. In the remainder we assume $\Delta=1.764k_{B}T_{c}$, and determine $T_{c}$ and the coefficients $\Delta_{\bf b}$ by solving the linearized gap equationWu et al. (2018); SM .

Using standard linear response theory we can obtain the expression for the superconducting weight Peotta and Törmä (2015); Liang et al. (2017); SM

[TABLE]

where $L\times L$ is the size of the two dimensional system, $n(E)$ is the Fermi distribution function, $E_{i}$, $|\psi_{i}({\bf k})\rangle$ are the eigenvalues and eigenvectors of $H_{BdG}$, and $\mu,\nu=x,y$ represent the directions. In the remainder, we focus on the case $\nu=\mu$. We have the velocity operators $\hat{v}_{\mu}({\bf k})=\partial H_{BdG}/\partial k_{\mu}$, $\hat{v}_{cf,\mu}({\bf k})=(1/2)\gamma_{z}\partial H_{BdG}/\partial k_{\mu}$, ($\gamma_{z}$ is the Pauli matrix acting in Nambu space).

Let $H_{+}$ and $H_{-}$ be the particle and hole Hamiltonians, respectively, of $H_{BdG}$, $|\psi_{\pm m}\rangle$ the eigenstates of $H_{\pm}$, $w_{\pm im}\equiv\langle\psi_{\pm m}|\psi_{i}\rangle$, and $v^{+}_{\mu}\equiv\partial_{\mu}{H_{+}}$, $v^{-}_{\mu}\equiv-\partial_{\mu}{H_{-}}$. In terms of these quantities we have Liang et al. (2017):

[TABLE]

where $m,n$ and $p,q$ index the particle and hole bands. The matrix elements with $m\neq n$ and $p\neq q$ in Eq. (7) represent pure interband contribution. By defining

[TABLE]

we can separate Eq. (7) into a conventional and a geometric part:

[TABLE]

Below we show that both the conventional Scalapino et al. (1992, 1993) and the geometric contribution Peotta and Törmä (2015); Liang et al. (2017) are important for the superfluid weight in TBLG.

Figure 1 shows the dispersion of the lower flat band without superconductivity. It is plotted in the moiré BZ for two different values of $\theta$: $\theta=1.05^{\circ}$, the magic angle, and $\theta=1.00^{\circ}$. For each angle we see a sharp feature in the dispersion at the $\gamma$ point, which is away from zero energy. For $\theta=1.05^{\circ}$ the bandwidth of the nearly flat moiré band is about 2 meV, whereas for $\theta=1.00^{\circ}$ it is around 5 meV. The bandwidth of the lowest energy bands and the value of the magic angle can differ between experiments Kerelsky et al. (2019), a fact that can be taken into account by tuning $w$, see SM . For $\theta=1.00^{\circ}$ we see that the bands exhibit deep and narrow valleys, green regions emanating from the $\gamma$ point. Around these valleys the quasiparticle energy $\epsilon({\bf k})$ varies rapidly with ${\bf k}$ producing high local velocity despite the fact that the bandwidth is only few meVs.

Figures 2 (a), (b) show the profile of $\epsilon({\bf k})$ for the lowest excitation in the presence of s-wave pairing, for $\theta=1.05^{\circ}$, and $\theta=1.00^{\circ}$. The amplitude and Fourier components of the superconducting gap are obtained by solving the mean field gap equation SM . We see that also in the presence of a superconducting gap the bands exhibit the same qualitative features as the bands with no pairing Fig 1.

Figures 2 (c), (d) show the momentum space profile of the integrand, $d^{s,\rm conv}_{xx}({\bf k})$, that enters the expression (8) for $D^{s,\rm conv}_{xx}$ for $\theta=1.05^{\circ},1.00^{\circ}$, respectively. We see that for $\theta=1.05^{\circ}$ $d^{s,\rm conv}_{xx}$ is peaked at $\gamma$ point, and is otherwise quite uniform and small. At bands crossings $d_{xx}^{s,{\rm geom}}$ is expected to be large as long as the Berry curvature is not zero, regardless of the nature of the crossing Rhim and Yang (2019). For $\theta=1.00^{\circ}$ $d^{s,\rm conv}_{xx}$ is strongly peaked at the position of the valleys that we identified in Fig. 1 (b). This clearly shows that the conventional contribution to $D^{s}$ can depend very strongly on the twist angle and in general cannot be assumed to be negligible despite the smallness of the bandwidth. The reason is that even for narrow bands, the expectation value of the velocity operators can be non-negligible. Figure. 2 (e), (f) show the profile of the integrand, $d^{s,\rm geom}_{xx}({\bf k})$, that enters the expression of $D^{s,\rm geom}_{xx}$ for the same conditions used to obtain panels (c) and (d). For $\theta=1.05^{\circ}$ $d^{s,\rm geom}_{xx}({\bf k})$ is strongly peak at the $\gamma$ point and on average is larger than the conventional term. This shows that at the magic angle the geometric contribution to $D^{s}_{xx}$ is significant and larger than the conventional contribution. For $\theta=1.00^{\circ}$, however, $d^{s,\rm conv}_{xx}({\bf k})$ is large in most of the moiré BZ so that the conventional contribution to $D^{s}_{xx}$ is larger than the geometric one. As the bands become flatter the conventional contribution, for fixed electron’s density, decreases and so we can expect its importance to decrease relative to the geometric contribution. Fig. 2 (g), (h) show the sum $d^{s,\rm conv}_{xx}({\bf k})+d^{s,\rm geom}_{xx}({\bf k})$. It is worth pointing out that the spin Chern number, $C$, of the lowest energy bands is zero, but in general $D^{s,\rm geom}_{xx}$ is nonzero even when $C=0$ Liang et al. (2017).

We continue by obtaining the dependence of $D^{s,\rm conv}_{xx}$, and $D^{s,\rm geom}_{xx}$ on the chemical potential. From the initial discussion we expect $D^{s,\rm conv}_{xx}$ to increase with the electron density and therefore with $|\mu|$. The scaling of $D^{s,\rm geom}_{xx}$ with respect to $\mu$ depends on the details of the quantum metric of the bands Liang et al. (2017). Fig. 3 shows the evolution of $D^{s,\rm conv}_{xx}$, $D^{s,\rm geom}_{xx}$, and $D^{s}_{xx}$ with $\mu$ for the cases of $\theta=1.05^{\circ}$ and $\theta=1.00^{\circ}$. To obtain these results the superconducting gap is obtained for each value of $\mu$. The results of Fig. 3 confirm the expectation that $D^{s,\rm conv}_{xx}$ increases with $|\mu|$, for both the magic angle and $\theta=1.00^{\circ}$. They also show that for both angles the geometric contribution decreases with $|\mu|$. Considering that $D^{s}_{xx}$ controls the critical temperature, $T_{\mathrm{KT}}$, for the Berezinskii-Kosterlitz-Thouless (BKT) phase transition Berezinskiǐ (1971); Kosterlitz and Thouless (1973), the results of Fig. 3 show that in TBLG it could be possible in principle to tune the nature of the transition, BCS, or BKT by simply tuning the chemical potential.

An increase of $\Delta$, keeping $\mu$ fixed, is expected to cause an increase of $D^{s}_{xx}$. This is confirmed by the results of Fig. 4. Again, we can see at $\theta=1.05^{\circ}$ the geometric contribution is significant while at $\theta=1.00^{\circ}$ the conventional contribution dominates.

In Figs. 5 (b), (d) we show the BKT transition temperature as a function of $\mu$ obtained from the equation $k_{B}T_{\mathrm{KT}}=\pi D^{s}(\Delta(T_{\mathrm{KT}}),T_{\mathrm{KT}})$, assuming $\Delta(T)=1.764k_{B}T_{c}\sqrt{1-T/T_{c}}$. The pre-factor on the r.h.s. of the equation for $T_{\mathrm{KT}}$ is twice $\pi/2$ due to the valley degeneracy. In Figs. 5 (a), (c) the curves with solid circles show the evolution of $D^{s}(\Delta(T),T)$ with $T$ for different values of $\mu$ for $\theta=1.05^{\circ}$ and $\theta=1.00^{\circ}$, respectively. The intersection of these curves with the solid line $k_{B}T$ returns the values of $T_{\mathrm{KT}}(\mu)$. We note that $T_{\mathrm{KT}}$ is fairly close to $T_{c}$, well above the lower bound set by previous studies Xie et al. (2019).

We obtained the value of $D^{s}_{xx}$ ($D^{s,\rm conv}_{xx}$, $D^{s,\rm geom}_{xx}$) for different twist angles using the corresponding values of $T_{c}$. The results are shown in Fig. 6 (a). We see that despite the fact that $T_{c}$ is lower for $\theta=1.10^{\circ}$ than for $\theta=1.05^{\circ}$ the superconducting weight is larger for $\theta=1.10^{\circ}$. This is because for $\theta=1.10^{\circ}$ the conventional contribution to $D^{s}_{xx}$ is much larger than at the magic angle. The results of Figs. 6(a) clearly show that $D^{s}_{xx}$ varies strongly with the twist angle, and that, as a function of $\theta$, the dominant contribution to $D^{s}_{xx}$ can either be the conventional one or the geometric. It is somewhat surprising, that even for twist angles as small as $1.00^{\circ}$, corresponding to a bandwidth of the lowest energy bands of just 5 meV, the conventional contribution is larger than the geometric one.

Figure 6 (b) shows the dependence of $T_{c}$ and $T_{\mathrm{KT}}$ on the twist angle. We see that both $T_{c}$ and $T_{\mathrm{KT}}$ are maximum at the magic angle and decrease rapidly for $\theta$ larger than the magic angle. The results of Fig. 6 (b) suggests that it may be possible to tune $T_{\mathrm{KT}}$ by tuning the twist angle. Taking into account finite size effects, this can change the nature of the normal-superconductor phase transition.

In summary, we have shown that in twisted bilayer graphene, despite the flatness of the low energy bands, the superconducting weight $D^{s}_{xx}$ is finite and large enough to explain the experimental observation of superconducting behavior in these systems. We find that the share of the geometric and conventional contributions to $D^{s}_{xx}$ depends on the twist angle: at the magic angle the geometric contribution dominates, for angles slightly away from the magic angle the conventional contribution dominates. This qualitative difference is also reflected in the scaling of $D^{s}_{xx}$ with $\mu$: at the magic angle $D^{s}_{xx}$ decreases with $|\mu|$, a somewhat surprising result due to the conventional contribution to $D^{s}_{xx}$ being almost independent of $\mu$ whereas the geometric, large, contribution decreases with $|\mu|$ at the magic angle. Away from the magic angle we find the more conventional behavior of $D^{s}_{xx}$ growing with $|\mu|$ as the conventional contribution, that grows with $|\mu|$, dominates. This has the important implication that at the magic angle, by simply increasing $|\mu|$ twisted bilayer graphene can be tuned into a regime for which the Berezinzki-Kosterlits-Thouless transition is significantly smaller than the BCS critical temperature. This result shows that twisted bilayer graphene is an exceptional system in which the nature, BKT or BCS, of the superconducting transition can be tuned and experimentally studied in unprecedented ways. We stress that the superfluid weight is one of the few observable signatures of the Berry phase in the Brillouin zone. Measurements of the superfluid weight with the known experimental techniques Hebard and Fiory (1980); Turneaure et al. (2000); Bert et al. (2011, 2012); Kapon et al. (2019) can directly test our predictions of its parametric dependencies.

Note: In the process of completing the manuscript, we became aware of a related recent preprint by Aleksi Julku et al. Julku et al. (2019).

I Supplemental material

I.1 I. Solutions of the mean field gap equation

To obtain an estimate of pairing parameters we solve the linearized gap equation Wu et al. (2018) $\Delta_{{\bf b}l}=\sum_{{\bf b}^{{}^{\prime}}l^{{}^{\prime}}}\chi_{{\bf b}{\bf b}^{{}^{\prime}}}^{ll^{{}^{\prime}}}\Delta_{{\bf b}^{{}^{\prime}}l^{{}^{\prime}}}$, where the pairing susceptibility is

[TABLE]

Here $\mathcal{N}({\bf q})=\frac{1-n_{F}[\epsilon_{n_{1}}({\bf q})]-n_{F}[\epsilon_{n_{2}}(\bf q)]}{\epsilon_{n_{1}}({\bf q})+\epsilon_{n_{2}}({\bf q})-2\mu}$, $U({\bf q})_{{\bf b}l}=\langle u_{n_{1}}({\bf q})|e^{i{\bf b}\cdot{\bf r}}|u_{n_{2}}({\bf q})\rangle_{l}$, $g_{0}$ is the pairing coupling constant due to electron-phonon interaction, $\mathcal{A}$ is the total area of the sample, $\epsilon_{n_{1,2}}({\bf q})$ and $u_{n_{1,2}}$ are eigenvalue and wavefunctions of the non-superconducting Hamiltonian $H_{TBL,{\bf K}}({\bf k})$ and $l(l^{{}^{\prime}})$ represents the layer index. We keep the reciprocal basis vectors whose lengths are no larger than twice of $b\equiv|{\bf b}_{1}|$.

When the temperature approaches $T_{c}$, the maximum eigenvalues of $\chi$ approach one. Correspondingly, the linearized gap equation becomes

[TABLE]

here ${\bf v}$ contains $\Delta_{0}$, $\Delta_{b}$ and higher order Fourier components of $\Delta({\bf r})$. Because the higher order terms are small, we can neglect them and keep only $\Delta_{0}$ and $\Delta_{b}$. The results are listed in table (1(c)).

Everywhere in the text we have asssumed $\Delta=1.764k_{B}T_{c}$ at zero temperature, where 1.764 is the prefactor for weak-coupling theory in metals. This assumption is made for simplicity because the actual prefactor in flat band systems depends on the details of the model.

I.2 II. Calculation of the superfluid weight

Our starting point is a BdG Hamiltonian for a singlet superconductor

[TABLE]

Here $H_{0}(\mathbf{p})$ is $n\times n$ normal state Hamiltonian and $\Delta(\mathbf{p})$ is the superconducting order parameter.

We assume that the order parameter supports a gradient of phase $\phi=\mathbf{k}\cdot\mathbf{r}$

[TABLE]

leading to a current

[TABLE]

where $D_{s}$ is the superfluid weight. We choose the units so that $e=1$ and $\hbar=1$ everywhere. The BKT transition temperature is given by

[TABLE]

We point out that in the literature several different conventions are used. In particular, one often uses $e^{i\phi(\mathbf{r})/2}\Delta(-i\nabla)e^{i\phi(\mathbf{r})/2}\to e^{i\phi(\mathbf{r})}\Delta(-i\nabla)e^{i\phi(\mathbf{r})}$ in Eq. (S4) and $2e\to e$ in Eq. (S5). With this convention $D_{s}$ becomes 4 times larger and $\pi/2$ is replaced by $\pi/8$, but this of course does not affect the predictions for $T_{\rm KT}$.

To diagonalize the Hamiltonian (S4) we use an ansatz

[TABLE]

This way we find that $a(\mathbf{p})$, $b(\mathbf{p})$ and corresponding eigenenergies $E(\mathbf{p})$ can be solved using the effective Hamiltonian

[TABLE]

The operator for quasiparticle velocity in the basis described by the coefficients $a(\mathbf{p})$, $b(\mathbf{p})$ is given by

[TABLE]

The current density is given by

[TABLE]

where $\gamma_{z}$ is the Pauli matrix in the Nambu space, $|n,\mathbf{p},\mathbf{k}\rangle$ are the eigenstates and $E_{n}(\mathbf{p},\mathbf{k})$ the corresponding eigenenergies of Hamiltonian $H_{\rm eff}(\mathbf{p},\mathbf{k})$.

The Hamiltonian (in Nambu space $C_{\mathbf{p},\mathbf{k}}^{\dagger}$) can be written as

[TABLE]

where

[TABLE]

and

[TABLE]

The current density operator can be written as

[TABLE]

where

[TABLE]

and

[TABLE]

We assume that the current flows parallel to the gradient of phase $\mathbf{k}=k\mathbf{e}_{\mu}$, $\langle\hat{\mathbf{j}}\rangle=j\mathbf{e}_{\mu}$ and by using linear response theory to obtain

[TABLE]

where

[TABLE]

$E_{i}(\mathbf{p})$ and $\psi_{i}(\mathbf{p})$ are the eigenenergies and eigenstates of $H_{BdG}(\mathbf{p})$. We have also used a shorthand notation $v_{\mu}(\mathbf{p})=v_{\mu}(\mathbf{p},\mathbf{k}=0)$. Therefore, we can identify

[TABLE]

Additionally we can simplify the expression using

[TABLE]

where $G(\omega_{n},\mathbf{p})=(i\omega_{n}-H_{BdG}(\mathbf{p}))^{-1}$ and $\omega_{n}$ are Matsubara frequencies. By assuming also that $\Delta(\mathbf{p})$ is independent of momentum $\mathbf{p}$ we obtain

[TABLE]

and arrive at the expression used in the main text.

For single isolated parabolic band at zero temperature the expressions above give $D_{s}\propto n/m^{*}$, where $n$ is the density and $m^{*}$ is the effective mass of the electrons. In this calculation one needs to take into account that the density of electrons within the band is constant, and therefore when the superconducting transition occurs the chemical potential is renormalized. In particular, in the limit of small density and large $m^{*}$ the chemical potential is renormalized well below the band. This can be understood as a BCS to BEC crossover driven by decreasing $n/(m^{*}\Delta)$. In a semimetal where both electron- and hole-like bands are present, such as TBLG, similar renormalization of the chemical potential below the conduction band (above the valence band) cannot occur because it would lead to appearance of large density of holes (electrons) in the valence band (conduction band). The densities in each band are not conserved in the transition separately (only the total density is conserved) and therefore the conventional contribution to the superfluid weight is necessarily much larger at small densities than expected for isolated bands.

In this manuscript we always consider a system consisting of two valleys ${\bf K}$ and ${\bf K}^{\prime}$ and calculate $D_{s}$ only for a single valley. The valley degeneracy gives an additional factor of 2, so that $\pi/2$ in Eq. (S6) has to be replaced with $\pi$.

To check our numerical codes we have calculated the superfluid weight for a single valley of TBLG numerically as a function of the effective interlayer tunneling amplitude $w$ in the case of uniform $s$-wave pairing [see Fig.S1 (a)]. Analytically, for $w=0$ one obtains

[TABLE]

where the first term is the conventional superfluid weight and the second term is the geometricLiang et al. (2017). The numerics demonstrate that the superfluid weight quickly approaches the analytical result when the coupling between layers is decreased. By increasing $w$ we find that $D_{s}/\Delta$ first increases and then decreases. Qualitatively similar results are obtained also in the case of spatially dependent $\Delta(\mathbf{r})$ with parameters from Table 1(c). However, in this calculation we have not self-consistently taken into account the influence of $w$ on the superconducting order parameter. Experimentally the effect of $w$ can be studied by applying pressure.

I.3 III. The tight binding model of twisted bilayer graphene

Following the approach described in Ref. Bistritzer and MacDonald, 2011; Wu et al., 2018, we can construct the tight binding model of TBLG. In TBLG, the conservation of the crystal momentum requires ${\bf k}_{b}={\bf k}_{t}+({\bf G}_{t}-{\bf G}_{b})$ where ${\bf G}_{i}$ is the reciprocal lattice wave-vector in layer $i$, with $i=t,b$ representing the top or bottom layer. For small twist angles a fairly accurate description is obtained by just keeping the tunneling processes for which $|{\bf k}_{b}-{\bf k}_{t}|=|{\bf G}_{t}-{\bf G}_{b}|=2K\sin(\theta/2)$, where $K=4\pi/3a_{0}$. There are three vectors ${\bf Q}_{i}={\bf G}_{t}-{\bf G}_{b}$ ($i=1,2,3)$ for which $Q=2K\sin(\theta/2)$ and so all the tunneling processes for which $|{\bf k}_{b}-{\bf k}_{t}|=Q$ are taken into account by keeping all the recursive tunneling processes on a honeycomb structure constructed in momentum space with nearest neighbor points connected by the vectors ${\bf Q}_{i}$. The honeycomb arrangement can be thought of as a triangular lattice. By keeping ${\bf b}=0,{\bf b}_{1},{\bf b}_{2}$ we can write out the tight binding Hamiltonian around the ${\bf K}$ point

[TABLE]

where $\boldsymbol{\kappa}_{+1}=(0,0)$ and $\boldsymbol{\kappa}_{-1}=(0,Q)$, and the basis is

[TABLE]

with $\phi_{{\bf k}}=(c_{{\bf k}A},c_{{\bf k}B})$. Here $c_{{\bf k}A,B}$ is the electron annihilation operator with momentum ${\bf k}$ at sublattice A,B. Similarly, We can include more ${\bf b}$ in this Hamiltonian.

I.4 IV. Dependence of the superfluid weight on the number of bands

It is interesting to study how the superfluid weight depends on the number of bands that are included in the calculation. Figure S2 shows the dependence of both the conventional and geometric part of $D^{s}$ on the number of bands. We see that the conventional part depends only weakly on the number of bands, $n_{\rm bands}$, but that the geometric part depends very strongly on $n_{\rm bands}$. We see that keeping only the two nearly flat bands is not enough to get accurate estimates of the geometric contribution. However, we find that when $n_{\rm bands}=10$, the $D^{s}_{\rm geom}$ is already very close (less than 2% away) to the value obtained keeping as many as 196 bands. This seems consistent with recent results that suggest that a minimal model for TBLG might require a minimum of 10 bands Po et al. (2019)