Tunneling and Fluctuating Electron-Hole Cooper Pairs in Double Bilayer Graphene

TL;DR

This paper develops a detailed theoretical model explaining the enhancement of tunneling conductance in double bilayer graphene due to fluctuating electron-hole Cooper pairs, predicting temperature and density dependence consistent with experiments.

Contribution

It introduces a comprehensive theory of fluctuating electron-hole Cooper pairs in double bilayer graphene, accounting for multiple pairing channels and their impact on tunneling conductance.

Findings

Fluctuating pairs enhance tunneling conductance via a fluctuational Josephson effect.

The theory predicts a temperature and density dependence of the zero bias conductance peak.

Cleaner samples could achieve higher transition temperatures up to ~50 K.

Abstract

A strong low-temperature enhancement of the tunneling conductance between graphene bilayers has been reported recently, and interpreted as a signature of equilibrium electron-hole pairing, first predicted in bilayers more than forty years ago but previously unobserved. Here we provide a detailed theory of conductance enhanced by fluctuating electron-hole Cooper pairs, which are a precursor to equilibrium pairing, that accounts for specific details of the multi-band double graphene bilayer system which supports several different pairing channels. Above the equilibrium condensation temperature, pairs have finite temporal coherence and do not support dissipationless tunneling. Instead, they strongly boost the tunneling conductivity via a fluctuational internal Josephson effect. Our theory makes predictions for the dependence of the zero bias peak in the differential tunneling conductance on temperature and electron-hole density imbalance that capture important aspects of the experimental observations. In our interpretation of the observations, cleaner samples with longer disorder scattering times would condense at temperatures $T_c$ up to $\sim 50 {\rm K}$, compared to the record $T_c \sim 1.5 $K achieved to date in the experiment.

Figures

Tables

Table 1. Table 1: Coupling constants λ α subscript 𝜆 𝛼 \lambda_{\alpha} for intra-valley Cooper pairs

Channel, $\alpha$	$s$-wave, ${\lambda }^s$	$d$-wave, ${\lambda }^d$
11	$c^4{N}_F{U}_s$	$2c^2{c}^2{N}_F{U}_d$
22	$s^4{N}_F{U}_s$	$2c^2{c}^2{N}_F{U}_d$
12-21	$2s^2{c}^2{N}_F{U}_s$	$\left(c^4+s^4\right)N_F{U}_d$

Table 2. Table 2: Coupling constants λ α subscript 𝜆 𝛼 \lambda_{\alpha} for inter-valley Cooper pairs

Channel, $\alpha$	$s$-wave, ${\lambda }^s$	$d$-wave, ${\lambda }^d$
$12$	$c^2{s}^2{N}_F{U}_s$	$\left(c^4+s^4\right)N_F{U}_d$
$21$	$c^2{s}^2{N}_F{U}_s$	$\left(c^4+s^4\right)N_F{U}_d$
$11$-$22$	$\left(c^4+s^4\right)N_F{U}_s$	$2c^2{s}^2{N}_F{U}_d$

Equations

$$H_0=\sum _{\mathbf{p}}\left[{\hat{\psi }}_{\mathrm{t}\mathbf{p}}^{+}({\hat{h}}_{\mathrm{t}\mathbf{p}}-{\mu }_{\mathrm{t}}){\hat{\psi }}_{\mathrm{t}\mathbf{p}}+{\hat{\psi }}_{\mathrm{b}\mathbf{p}}^{+}({\hat{h}}_{\mathrm{b}\mathbf{p}}+{\mu }_{\mathrm{b}}){\hat{\psi }}_{\mathrm{b}\mathbf{p}}\right].$$

$${\hat{h}}_{\mathrm{t}\mathbf{p}}=\left(\begin{array}{cc}u& \frac{p_{{\bar{v}}_{\mathrm{t}}}^2}{2m}\\ \frac{p_{v_{\mathrm{t}}}^2}{2m}& -u\end{array}\right),{\hat{h}}_{\mathrm{b}\mathbf{p}}=\left(\begin{array}{cc}-u& \frac{p_{v_{\mathrm{b}}}^2}{2m}\\ \frac{p_{{\bar{v}}_{\mathrm{b}}}^2}{2m}& u\end{array}\right).$$

$$|\mathrm{tc}\mathbf{p}⟩=\left(\begin{array}{c}{\mathrm{c}}_{\mathbf{p}}{e}^{-i{v}_{\mathrm{t}}{\varphi }_{\mathbf{p}}}\\ {\mathrm{s}}_{\mathbf{p}}{e}^{i{v}_{\mathrm{t}}{\varphi }_{\mathbf{p}}}\end{array}\right),|\mathrm{bv}\mathbf{p}⟩=\left(\begin{array}{c}{\mathrm{c}}_{\mathbf{p}}{e}^{i{v}_{\mathrm{b}}{\varphi }_{\mathbf{p}}}\\ -{\mathrm{s}}_{\mathbf{p}}{e}^{-i{v}_{\mathrm{b}}{\varphi }_{\mathbf{p}}}\end{array}\right).$$

$$H_{\mathrm{t}}=T^{+}+T=\sum _{\mathbf{p}}\left[{\hat{\psi }}_{\mathrm{b}\mathbf{p}}^{+}{\hat{t}}^{+}{\hat{\psi }}_{\mathrm{t}\mathbf{p}+\mathbf{Q}}+{\hat{\psi }}_{\mathrm{t}\mathbf{p}+\mathbf{Q}}^{+}\hat{t}{\hat{\psi }}_{\mathrm{b}\mathbf{p}}\right].$$

$$H_{\mathrm{int}}=\sum _{\mathbf{p}{\mathbf{p}}^{\prime }\mathbf{q}}\sum _{{\sigma }_{\mathrm{t}}{\sigma }_{\mathrm{b}}}{U}_{\mathbf{q}}{\psi }_{\mathrm{t},\mathbf{p}+\mathbf{q},{\sigma }_{\mathrm{t}}}^{+}{\psi }_{\mathrm{b},{\mathbf{p}}^{\prime }-\mathbf{q},{\sigma }_{\mathrm{b}}}^{+}{\psi }_{\mathrm{b}{\mathbf{p}}^{\prime }{\sigma }_{\mathrm{b}}}{\psi }_{\mathrm{t}\mathbf{p}{\sigma }_{\mathrm{t}}}.$$

$${\hat{\Gamma }}_{l^{\prime }l}=U_l{\delta }_{l^{\prime }l}+\sum _{l^{\prime \prime }}{U}_{l^{\prime }}{\hat{M}}_{l^{\prime }-l^{\prime \prime }}\Pi {\hat{\Gamma }}_{l^{\prime \prime }l}.$$

$$\Pi (\omega ,\mathbf{q})=N_{\mathrm{F}}\left\{\ln \left[\frac{{ϵ}_{\mathrm{c}}}{2\pi T}\right]-\frac{1}{2}\sum _{\zeta =\pm }{⟨\Psi \left(\frac{1}{2}+\frac{i[\omega +\zeta (h+v_{\mathrm{F}}q\cos {\varphi }_{\mathbf{p}})]+\gamma }{4\pi T}\right)⟩}_{{\varphi }_{\mathbf{p}}}\right\}.$$

$${\hat{M}}_{{\sigma }_b^{\prime },{\sigma }_b}^{{\sigma }_t^{\prime },{\sigma }_t}({\varphi }_{\mathbf{p}})=⟨{\sigma }_{\mathrm{t}}^{\prime }|\mathrm{tc}\mathbf{p}+\frac{\mathbf{q}}{2}⟩⟨\mathrm{tc}\mathbf{p}+\frac{\mathbf{q}}{2}|{\sigma }_{\mathrm{t}}⟩⟨{\sigma }_{\mathrm{b}}|\mathrm{bv}\mathbf{p}-\frac{\mathbf{q}}{2}⟩⟨\mathrm{bv}\mathbf{p}-\frac{\mathbf{q}}{2}|{\sigma }_{\mathrm{b}}^{\prime }⟩.$$

$$\hat{M}({\varphi }_{\mathbf{p}})=\left(\begin{array}{cccc}{c}^4& -c^{3}s{e}^{-2i{v}_{\mathrm{b}}{\varphi }_{\mathbf{p}}}& c^{3}s{e}^{-2i{v}_{\mathrm{t}}{\varphi }_{\mathbf{p}}}& -c^2{s}^2{e}^{-2i(v_{\mathrm{t}}+v_{\mathrm{b}}){\varphi }_{\mathbf{p}}}\\ -c^{3}s{e}^{2i{v}_{\mathrm{b}}{\varphi }_{\mathbf{p}}}& c^2{s}^2& -c^2{s}^2{e}^{2i(v_{\mathrm{b}}-v_{\mathrm{t}}){\varphi }_{\mathbf{p}}}& c{s}^3{e}^{-2i{v}_{\mathrm{t}}{\varphi }_{\mathbf{p}}}\\ c^{3}s{e}^{2i{v}_{\mathrm{t}}{\varphi }_{\mathbf{p}}}& -c^2{s}^2{e}^{2i(v_{\mathrm{t}}-v_{\mathrm{b}}){\varphi }_{\mathbf{p}}}& c^2{s}^2& -c{s}^3{e}^{-2i{v}_{\mathrm{b}}{\varphi }_{\mathbf{p}}}\\ -c^2{s}^2{e}^{2i(v_{\mathrm{t}}+v_{\mathrm{b}}){\varphi }_{\mathbf{p}}}& c{s}^3{e}^{2i{v}_{\mathrm{t}}{\varphi }_{\mathbf{p}}}& -c{s}^3{e}^{2i{v}_{\mathrm{b}}{\varphi }_{\mathbf{p}}}& s^4\end{array}\right).$$

$$c^2=\frac{1}{2}\left(1+\frac{u}{{ϵ}_{\mathrm{F}}}\right),s^2=\frac{1}{2}\left(1-\frac{u}{{ϵ}_{\mathrm{F}}}\right).$$

$${\hat{M}}_0=\left(\begin{array}{cccc}{c}^4& 0& 0& 0\\ 0& c^2{s}^2& -c^2{s}^2& 0\\ 0& -c^2{s}^2& c^2{s}^2& 0\\ 0& 0& 0& s^4\end{array}\right).$$

$$\frac{{\hat{\Gamma }}_{00}}{U_{\mathrm{s}}}=\left(\begin{array}{cccc}\frac{1}{L_{11}}& 0& 0& 0\\ 0& \frac{1-c^2{s}^2{U}_{\mathrm{s}}\Pi }{L_{12-21}}& \frac{-c^2{s}^2{U}_{\mathrm{s}}\Pi }{L_{12-21}}& 0\\ 0& \frac{-c^2{s}^2{U}_{\mathrm{s}}\Pi }{L_{12-21}}& \frac{1-c^2{s}^2{U}_{\mathrm{s}}\Pi }{L_{12-21}}& 0\\ 0& 0& 0& \frac{1}{L_{22}}\end{array}\right).$$

$${\hat{M}}_0=\left(\begin{array}{cccc}{c}^4& 0& 0& -c^2{s}^2\\ 0& c^2{s}^2& 0& 0\\ 0& 0& c^2{s}^2& 0\\ -c^2{s}^2& 0& 0& s^4\end{array}\right).$$

$$\frac{{\hat{\Gamma }}_{00}}{U_{\mathrm{s}}}=\left(\begin{array}{cccc}\frac{1-s^4{U}_{\mathrm{s}}\Pi }{L_{11-22}}& 0& 0& \frac{-c^2{s}^2{U}_{\mathrm{s}}\Pi }{L_{11-22}}\\ 0& \frac{1}{L_{12}}& 0& 0\\ 0& 0& \frac{1}{L_{21}}& 0\\ \frac{-c^2{s}^2{U}_{\mathrm{s}}\Pi }{L_{11-22}}& 0& 0& \frac{1-c^4{U}_{\mathrm{s}}\Pi }{L_{11-22}}\end{array}\right).$$

$$\ln \left[\frac{T_0}{{\bar{T}}_0}\right]+\frac{1}{2}\sum _{\zeta =\pm }{⟨\Psi \left(\frac{1}{2}+\frac{i\zeta (h+v_{\mathrm{F}}{Q}_0\cos {\varphi }_{\mathbf{p}})+\gamma }{4\pi T}\right)-\Psi \left(\frac{1}{2}\right)⟩}_{{\varphi }_{\mathbf{p}}}=0.$$

$$L_0^{-1}(\omega ,\mathbf{q})=\frac{1}{{\lambda }_0}\frac{1}{i\omega \tau -{ϵ}_{\mathbf{q}}},{ϵ}_{\mathbf{q}}=ϵ+\frac{{\xi }^2{\mathbf{q}}^2}{{\hslash }^2}.$$

$$\tau =\frac{\hslash {\Psi }^{\prime }\left(\frac{1}{2}+\frac{\gamma }{4\pi T}\right)}{4\pi T},\xi =\frac{\hslash v_F|{\Psi }^{\prime \prime }\left(\frac{1}{2}+\frac{\gamma }{4\pi T}\right){|}^{\frac{1}{2}}}{8\pi T}.$$

$$G_{\mathrm{T}}(V)=\frac{8A{e}^2}{\hslash }\frac{\mathrm{Im}[\chi (eV,\mathbf{Q})]}{eV}.$$

$$\chi (\tau ,\mathbf{q})=-⟨T_{\mathrm{M}}T(\tau ,\mathbf{q})T^{+}(0,\mathbf{q})⟩.$$

$${\chi }_0(\omega ,\mathbf{q})={\hat{t}}^{+}{\hat{M}}_0\Pi \hat{t},$$

$${\hat{t}}_l^{\prime }=t{\delta }_{l0}+\sum _{l^{\prime }}{U}_l{\hat{M}}_{l-l^{\prime }}\Pi {\hat{t}}_{l^{\prime }}^{\prime }.$$

$${\chi }^0=(c^4|t_{11}{|}^2+s^4|t_{22}{|}^2+|t_{12}-t_{21}{|}^2{c}^2{s}^2)\Pi$$

$$\chi =\left(\frac{c^4|t_{11}{|}^2}{L_{11}}+\frac{s^4|t_{22}{|}^2}{L_{22}}+\frac{|t_{12}-t_{21}{|}^2{c}^2{s}^2}{L_{12-21}}\right)\Pi .$$

$${\chi }_0=(c^2{s}^2|t_{12}{|}^2+c^2{s}^2|t_{21}{|}^2+|c^2{t}_{11}-s^2{t}_{22}{|}^2)\Pi .$$

$$\chi =\left(\frac{c^2{s}^2|t_{12}{|}^2}{L_{12}}+\frac{c^2{s}^2|t_{21}{|}^2}{L_{21}}+\frac{|c^2{t}_{11}-s^2{t}_{22}{|}^2}{L_{11-22}}\right)\Pi .$$

$$G_{\mathrm{T}}^0=\frac{8A{e}^2}{\hslash }\frac{s^4|t_{22}{|}^2\mathrm{Im}[\Pi (eV,\mathbf{Q})]}{eV}.$$

$$G_{\mathrm{T}}=\frac{8A{e}^2}{\hslash }\left(\frac{s^4|t_{22}{|}^2}{|L_{22}{|}^2}+\frac{|t_{12}-t_{21}{|}^2{c}^2{s}^2}{|L_{12-21}{|}^2}\right)\frac{\mathrm{Im}[\Pi (eV,\mathbf{Q})]}{eV}.$$

$$G_{\mathrm{T}}=\frac{8A{e}^2}{\hslash }\frac{s^4|t_{22}{|}^2}{|L_{11-22}{|}^2}\frac{\mathrm{Im}[\Pi (eV,\mathbf{Q})]}{eV}.$$

$$F(eV,\mathbf{Q})=\frac{\mathrm{Im}[\Pi ]}{eV|L_0{|}^2}=\frac{\tau }{(eV\tau {)}^2+{ϵ}_{\mathbf{Q}}^2}.$$

$${\Gamma }_{{\sigma }_{\mathrm{b}}^{\prime }{\sigma }_{\mathrm{b}}}^{{\sigma }_{\mathrm{t}}^{\prime }{\sigma }_{\mathrm{t}}}(i{p}_n,{\mathbf{p}}^{\prime },\mathbf{p})=U_{{\mathbf{p}}^{\prime }-\mathbf{p}}{\delta }_{{\sigma }_{\mathrm{b}}^{\prime }{\sigma }_{\mathrm{b}}}^{{\sigma }_{\mathrm{t}}^{\prime }{\sigma }_{\mathrm{t}}}+T\sum _{i{\omega }_n{\mathbf{p}}^{\prime \prime }}{U}_{{\mathbf{p}}^{\prime }-{\mathbf{p}}^{\prime \prime }}{G}_{{\sigma }_t^{\prime }{\sigma }_t^{\prime \prime }}^{\mathrm{t}}(i{p}_n+i{\omega }_n,{\mathbf{p}}_{+}^{\prime \prime })G_{{\sigma }_{\mathrm{b}}^{\prime \prime }{\sigma }_{\mathrm{b}}^{\prime }}^{\mathrm{b}}(i{\omega }_n,{\mathbf{p}}_{-}^{\prime \prime }){\Gamma }_{{\sigma }_{\mathrm{b}}^{\prime \prime }{\sigma }_{\mathrm{b}}}^{{\sigma }_{\mathrm{t}}^{\prime \prime }{\sigma }_{\mathrm{t}}}(i{p}_n,{\mathbf{p}}^{\prime \prime },\mathbf{p}).$$

Full text

††preprint: APS/123-QED

I I. Introduction

The possibility of Cooper pairing in a system with spatially separated electrons and holes in semiconductor quantum wells was first anticipated more than forty years ago Lozovik and Yudson (1975); Shevchenko (1976). According to theory, strong Coulomb interactions allow pairing at elevated temperatures, which would provide a physical realization of dipolar superfluidity that is potentially relevant for applications. The paired state is fragile however, and can be suppressed by disorder Efimkin et al. (2011); Bistritzer and MacDonald (2008), or by Fermi-line mismatches due to the differences between electron and hole anisotropies Efimkin and Lozovik (2011); Seradjeh (2012) that are always present in conventional semiconductors Conti et al. (1998). In fact equilibrium pairing has until recently been observed only in the presence of strong magnetic fields that quench the kinetic energies of electrons and holes and drive the system to the regime of strong correlations Eisenstein and MacDonald (2004); Eisenstein (2014).

Recent progress in fabricating single-atomic-layer two-dimensional materials has renewed interest in electron-hole pairing Lozovik and Sokolik (2008); Min et al. (2008); Zhang and Joglekar (2008); Efimkin et al. (2012); Sodemann et al. (2012); Mink et al. (2011); Lozovik et al. (2012); Lozovik and Sokolik (2010); Perali et al. (2013); Conti et al. (2017); Neilson et al. (2014); Shumway and Gilbert (2012); Su and MacDonald (2017); Debnath et al. (2017); Abergel et al. (2013). Graphene-based two-dimensional electron systems not only have high mobility and almost perfect electron-hole symmetry but they make it possible to fabricate closely-spaced, and therefore strongly interacting, independently gated and contacted double layer structures. Very recently low-temperature enhancement of the tunneling conductance between graphene bilayers has been observed at matched concentrations of electrons and holes Burg et al. (2018). A typical conductance trace is presented in Fig. 1, where we see a tunneling conductance that appears to diverge at $T_{0}\approx 1.5\;\hbox{K}$, signaling equilibrium pair condensation. This observation provides the first clear experimental signature of equilibrium electron-hole pair condensation in the absence of magnetic field Com (a, b).

Enhanced tunneling conductance has been observed previously in semiconductors bilayers in the strong field quantum Hall regime, Spielman et al. (2000) and has been interpreted as an internal Josephson effect Eisenstein and MacDonald (2004). The differential conductance, does not diverge however, and instead has a sharp peak at zero bias. The property that the conductance peak width is smaller than temperature, and smaller than the single-electron scattering rate (i.e. the Landau level width) nevertheless points to a collective origin of the peak. Bilayers in the quantum Hall regime are predicted to support dissipationless Josephson-like tunneling currents in the presence of long-range electron-hole coherence Wen and Zee (1993); Ezawa and Iwazaki (1993). The development of a quantitative theory of enhanced tunneling in quantum Hall systems Stern et al. (2001); Hyart and Rosenow (2011); Fogler and Wilczek (2001); Jack et al. (2004) that fully explains the peak width has been challenged by the importance of inhomogeneity and disorder, and by strong interactions in the presence of dispersionless Landau levels. Phase fluctuations that are inevitably present due to the two-dimensional Berezinskii-Kosterlitz-Thouless nature of the phase transition Berezinskii (1971, 1972); Kosterlitz and Thouless (1973); Kosterlitz (1974) also play a role. The theory of enhanced tunneling is simpler at zero magnetic field, at least in the weak coupling regime where the electron-hole pairing energy is small compared to the Fermi energy, allowing experiments to be explained more fully as we demonstrate below.

The enhancement of the tunneling conductance in the double bilayer graphene system has been observed over a wide density range 4\cdot 10^{10}$$\,\sim\,$$10^{12}\,\hbox{cm}^{-2}, where electronic correlations vary from moderate to weak López Ríos et al. (2018). The Bardeen-Cooper-Schrieffer (BCS) theory of electron-hole pairing has greater validity at weaker pairing. True internal Josephson behavior in this case occurs only below a critical temperature $T_{0}$ and is preceded by enhancement of the tunneling conductance that diverges as $T_{0}$ is approached as illustrated in Fig. 1. This critical behavior has been predicted by one of the authors Efimkin and Lozovik (2013a) and has been interpreted as a fluctuational internal Josephson effect. It originates from partly coherent fluctuating electron-hole Cooper pairs Efimkin and Lozovik (2013a, b); Rist et al. (2013) that are a precursor of equilibrium pairing and reminiscent of Aslamazov-Larkin and related effects in superconductors Varlamov et al. (2018); Larkin and Varlamov (2005); Skocpol and Tinkham (1975). Above $T_{0}$ fluctuating Cooper pairs have a finite coherence time Perali et al. (2013, 2002); Mueller (2017); Franz (2007) and cannot support a dissipationless tunneling current. While the recent observations do qualitatively agree with earlier theory, the double bilayer graphene system has some important differences compared to the double parabolic electron gas models considered previously. These are related to the system’s well known $2\pi$ momentum space Berry phases. We show here that accounting properly for these differences provides a better account of the low-temperature tunneling anomalies.

In the present work we have developed a theory of the fluctuational internal Josephson effect in a system with closely spaced graphene bilayers. As we show below, the presence of valley and sublattice degrees of freedom provides three competing electron-hole channels for both intra-valley and inter-valley Cooper pairs. We show that three channels are nearly independent and have have different condensation temperatures and different sublattice structure. The experimental enhancement of the tunneling conductance by fluctuating Cooper pairs can be explained only by the presence of competing channels that dominate in different temperature ranges. In the vicinity of $T_{0}$, the tunneling conductance at zero bias is predicted to have a critical divergence $G_{\mathrm{T}}\sim(T-T_{0})^{-2}$, that matches the experimental data well as we see in Fig. 1. The calculated dependence of the tunneling conductivity on inter-layer voltage bias and carrier-density imbalance also match the experimental data Burg et al. (2018) reasonably. We conclude that the observed enhancement of the tunneling conductance in double bilayer graphene is well explained by our fluctuational internal Josephson effect theory.

The paper is organized as follows. In Sec. II we introduce a model that describes the low-energy physics of two closely spaced graphene bilayers. Sec. III is devoted to a description of fluctuating Cooper pairs above the critical temperature $T_{0}$. In Sec. IV we use these results as a starting point for a theory of the tunneling conductance. In Sec. V we compare our calculations with the recent experimental data. Finally in Sec. VI we discuss limitations of our theory and aspects of the experimental data that are still not well understood, and present our conclusions.

II II. Model

II.1 A. Weak and strong coupling regimes

The system of interest contains two graphene bilayers separated by an insulator as sketched in Fig. 2-a. An external electric field perpendicular to the bilayers induces an excess of electrons in the top bilayer ($\mathrm{t}$) and their deficit in the bottom bilayer ($\mathrm{b}$). Inevitably, it also results in gaps $2|u|$ in isolated bilayer graphene spectra. The latter are in Fig. 2-b in the case of matched concentrations for electrons and holes, the most favorable regime for the Cooper pairing. Double bilayer physics is very rich and has been considered in a number of recent papers Perali et al. (2013); Conti et al. (2017); Neilson et al. (2014); Su and MacDonald (2017); López Ríos et al. (2018) that aim to provide realistic predictions of the critical temperature $T_{0}$ based on the microscopic model. Here we follow a different route and consider a minimal phenomenological model that accounts for disorder and describes the instability of the system towards electron-hole pairing in the weak coupling regime and can explain the observed enhancement of the tunneling conductance between graphene bilayers Burg et al. (2018).

The physics of electron-hole Copper pairing depends on three dimensionless parameter: Wigner-Seitz interaction strength parameter $r_{\mathrm{s}}=me^{2}/\kappa\hbar^{2}k_{\mathrm{F}}$ that scales the ratio of interactions and kinetic energy in an individual bilayer Com (c), a parameter $k_{\mathrm{F}}d$ that scales the inter-layer Coulomb interactions respect to the intra-layer one, and a parameter $|u|/\epsilon_{\mathrm{F}}$ that determines the low-energy spectrum of individual bilayers. Here $k_{\mathrm{F}}$ and $\epsilon_{\mathrm{F}}$ are Fermi wave vector and energy for electrons and holes, $m$ is their effective constant, $\kappa$ and $d$ are a dielectric constant for the spacer between bilayers and its thickness. If $r_{\mathrm{s}}\ll 1,$ $k_{\mathrm{F}}d\gg 1$ and $|u|/\epsilon_{\mathrm{F}}\ll 1$ the system is in the weak coupling regime with pairing correlations only in the vicinity of Fermi level for electrons and holes. This regime can be described by the BCS theory for electron-hole pairing. The nature of the strong coupling regime ($r_{\mathrm{s}}\gg 1$ and $k_{\mathrm{F}}d\sim 1$) depends on the ratio $|u|/\epsilon_{\mathrm{F}}$. If $|u|\gg\epsilon_{\mathrm{F}}$ the state is a Bose-Einstein condensate (BEC) of indirect excitons that represent a bound state of electron and hole. It has been argued that the most favorable conditions for observation of electron-hole condensation are reached near the mid-point of the BEC-BCS crossover Perali et al. (2013); Conti et al. (2017); Neilson et al. (2014); Gortel and Swierkowski (1996); López Ríos et al. (2018). In the opposite case the paired state is a multi-band BCS-like paired state Sodemann et al. (2012); Mink et al. (2011); Lozovik et al. (2012); Lozovik and Sokolik (2010); Conti et al. (2019) where pairing correlations also span to remote bands (valence band in the layer with excess of electrons and conduction band in the layer with excess of holes).

Enhanced conductance has been seen over a range of electron and hole densities that covers $4\cdot 10^{10}\sim 10^{12}\;\hbox{cm}^{-2}$ and corresponds to $r_{\mathrm{s}}=3.1\sim 0.72$ and $k_{\mathrm{F}}d=0.07\sim 0.35$ Est . This range suggests to the transition from moderate to weak coupling regime occurs with increasing of electron and hole concentrations. In the experimental setup Burg et al. (2018) the Fermi energy $\epsilon_{\mathrm{F}}$ of charge carriers and the gap $2|u|$ in the spectrum of individual bilayers are not controlled independently, but with the same gate (that produces the electric field perpendicular to bilayers). While the density can be tuned within a wide range, the ratio $|u|/\epsilon_{\mathrm{F}}\approx d_{\mathrm{BG}}/d\approx 0.15$ is approximately fixed. Here $d_{\mathrm{BG}}$ is the thickness of individual bilayers. This favors the multi-band BCS-like state at strong coupling and extends the region of the applicability for the phenomenological weak coupling BCS theory that will be employed below.

II.2 B. Phenomenological model in the weak coupling regime

The spectrum sketched in Fig. 2-b has spin and valley degeneracy. The spin degrees of freedom simply added a factor of $2$ in the tunneling conductance when the condensed state preserves spin-invariance and do not need to be treated explicitly. The low energy states in bilayer graphene are concentrated around two inequivalent valleys $K$ ($v=1$) and $K^{\prime}$ ($v=-1$) situated at the corners of the first Brillouin zone, and are described by the two-band Hamiltonian McCann and Koshino (2013)

[TABLE]

where $\hat{\psi}_{\mathbf{p}}=\{\psi_{\mathrm{l}1\mathbf{p}},\psi_{\mathrm{l}2\mathbf{p}}\}$ is a spinor of annihilation operators for electrons in both layers with sublattice indexes $\sigma=1,2$, that are numerated according to the sketch in Fig. 2-a. We assume that there is a deficit of electrons in the bottom layer, but it is instructive not to perform the transformation to field operators of holes. $\mu_{\mathrm{t}}=\epsilon_{\mathrm{F}}+h$ and $\mu_{\mathrm{b}}=\epsilon_{\mathrm{F}}-h$ characterize the electric potentials in the top and bottom bilayers. Here $\epsilon_{\mathrm{F}}$ is the average of the electron and hole Fermi energies at neutrality, while $2h\ll\epsilon_{\mathrm{F}}$ is their difference. The matrices $\hat{h}_{\mathrm{t}\mathbf{p}}$ and $\hat{h}_{\mathrm{b}\mathbf{p}}$ are

[TABLE]

Here $v=\pm 1$ and $\bar{v}=\mp 1$ are valley indexes, $m$ is the electron mass, and $p_{v}=p_{x}+ivp_{y}$. The electric field perpendicular to bilayers opens a gap $2|u|$ separating conduction $\epsilon_{\mathrm{c}\mathbf{p}}=\epsilon_{\mathbf{p}}$ and valence $\epsilon_{\mathrm{v}\mathbf{p}}=-\epsilon_{\mathbf{p}}$ bands in each bilayer with $\epsilon_{\mathbf{p}}=\sqrt{u^{2}+(p^{2}/2m)^{2}}$. In the weak coupling regime the pairing correlations appear in the vicinity of Fermi lines for electrons and holes, and the presence of remote bands (valence band in the layer with excess of electrons and conduction band in the layer with excess of holes) can be neglected. In this regime only the conduction band of the top bilayer and the valence band of the bottom bilayer are relevant, and the corresponding spinor wave functions are

[TABLE]

Here $\phi_{\mathbf{p}}$ is the polar angle; $c_{\mathbf{p}}=\cos(\vartheta_{\mathbf{p}}/2)$ and $s_{\mathbf{p}}=\sin(\vartheta_{\mathbf{p}}/2)$ with $\cos(\vartheta_{\mathbf{p}})=u/\epsilon_{\mathbf{p}}$. The spinors have valley dependent chirality $\pm v_{\mathrm{t}(\mathrm{b})}\phi_{\mathbf{p}}$ that defines a sublattice structure of fluctuating electron-hole Cooper pairs as will be shown below. We introduce disorder with the help of phenomenological scattering rates $\gamma_{\mathrm{t}(\mathrm{b})}$. It is important in this theory to observe that because the electron and hole components of the Cooper pair are spatially separated and have opposite charges, both short-range and long-range Coulomb disorder lead to pair-breaking Efimkin et al. (2011); Bistritzer and MacDonald (2008).

In experiment Burg et al. (2018) the relative angle $\theta$ between graphene bilayers can be adjusted. Since valleys $K$ and $K^{\prime}$ reside at the corners of the first Brillouin zone, the valley momenta in the two layers do not match in the presence of a twist. For momentum conserving tunneling, the current is maximized when the layers are aligned ($\theta=0$) or twisted by $\theta=n\pi/3$. For even $n$, valley K (K’) in one layer is aligned with valley K (K’) in the other layer whereas for $n$ odd valley K (K’) in one layer is aligned with valley K’ (K) in the other layer. When states are labeled by their momenta relative to the Brillouin-zone corners, the tunneling Hamiltonian for $\theta$ close to $n\pi/3$

[TABLE]

Here $\mathbf{Q}=\mathbf{Q}_{\theta}+\mathbf{Q}_{\mathbf{B}}$ is the momentum splitting between valleys in the different layers. The twist contribution at small relative angle $\theta\ll 1$ can be approximated as $\mathbf{Q}_{\theta}=-[\mathbf{q}_{\mathrm{K}}\times\mathbf{e}_{\mathrm{z}}]\theta$ where $\mathbf{q}_{\mathrm{K}}$ is the momentum for the Brillouin-zone corner in bilayer graphene. It has opposite directions for valleys $K$ and $K^{\prime}$, and its magnitude is $q_{\mathrm{K}}=4\pi\hbar/3a_{0}$ with $a_{0}$ the corresponding Bravais lattice period. The contribution induced by an in-plane magnetic field $\mathbf{B}_{||}$ is the same for two valleys and is equal to $Q_{\mathbf{B}}=ed[\mathbf{B}_{||}\times\mathbf{e}_{\mathrm{z}}]/\hbar c$. Because each bilayer is represented by a two-band model, the matrix $\hat{t}$ has four matrix elements which we treat in a phenomenological way below, with the expectation that since $t_{22}$ corresponds to the tunneling between adjacent sublayers while $t_{11}$ to the tunneling between remote sublayers, $|t_{11}|\ll|t_{12}|\approx|t_{21}|\ll|t_{22}|$.

III III.Fluctuating Cooper pairs

III.1 A. The Cooper instability

Due to Coulomb interactions between electrons and holes, the double-bilayer system is unstable towards Cooper pairing. Here we omit repulsive interactions within each graphene bilayer since its main effect in the considered weak coupling regime is a simple renormalization of the quasiparticle spectra. The inter-bilayer attraction

[TABLE]

Here $U_{\mathbf{q}}$ is the screened Coulomb potential estimated in our previous work Efimkin et al. (2012); Sodemann et al. (2012); Mink et al. (2011). Below we employ a multipole decomposition of the interaction and set the momenta magnitudes to the Fermi momentum so that $U_{\mathbf{q}}$ reduces to a constant $U_{l}$ for each orbital angular momentum channel $l$. The set of $U_{l}$ parameters are also treated as phenomenological parameters with the expectation that the $s$-wave moment $U_{\mathrm{s}}\equiv U_{0}$ is largest. Since valleys are well separated in momentum space we neglect inter-valley scattering for electrons and holes.

The instability of double layer system towards electron-hole Cooper pairing is signaled by divergence of the electron-hole channel scattering vertex $\Gamma^{\sigma^{\prime}_{\mathrm{t}}\sigma_{\mathrm{t}}}_{\sigma^{\prime}_{\mathrm{b}}\sigma_{\mathrm{b}}}(\omega,\mathbf{p}^{\prime},\mathbf{p},\mathbf{q})$ at zero frequency $\omega$. Note that in our approximation scattering conserves valley indices for electrons ($v_{\mathrm{t}}$) and holes ($v_{\mathrm{b}}$). The $\Gamma$-vertex satisfies the Bethe-Saltpeter equation presented in Fig. 3-a, and is algebraic within the multipole approximation. It is instructive to combine $\sigma_{\mathrm{t}}$ and $\sigma_{\mathrm{b}}$ into a single index $|\sigma_{\mathrm{t}}\sigma_{\mathrm{b}}\rangle$ that varies between $1$ and $4$ as $|1\rangle=|11\rangle$, $|2\rangle=|12\rangle$, $|3\rangle=|21\rangle$ and $|4\rangle=|22\rangle$. With this definition the Bethe-Saltpeter equation can be written (See Appendix A for details) in a compact matrix form:

[TABLE]

Here momentum and frequency dependence are suppressed, $l$ ($l^{\prime}$) is the orbital momentum for the relative motion of two particles before (after) scattering, and $\hat{\Gamma}_{l^{\prime}l}$ is the corresponding scattering matrix. We have separated a factor of $\Pi(\omega,\mathbf{q})$ which also appears as the single-step pair propagator in the Cooper ladder sum of a bilayer system without sub lattice degrees of freedom:

[TABLE]

Here $\epsilon_{\mathrm{c}}$ is an energy cutoff that is required for momentum-independent interactions. The average $\langle\cdot\cdot\cdot\rangle_{\phi_{\mathbf{p}}}$ is calculated respect to a polar angle $\phi_{\mathbf{p}}$. $\Psi(x)$ is the logarithmic derivative of $\Gamma$-function (or the digamma function), while $\zeta=\pm 1$ is the summation index. $\gamma=\gamma_{\mathrm{t}}+\gamma_{\mathrm{b}}$ is the pair breaking rate, which is the sum of the scattering rates for electrons $\gamma_{\mathrm{t}}$ and holes $\gamma_{\mathrm{b}}$. The expression for $\Pi(\omega,\mathbf{q})$ in Eq. (7) is well known Varlamov et al. (2018); Larkin and Varlamov (2005) from previous work on systems without layer degrees of freedom. The chiral nature of the bilayer graphene charge carriers is captured by the nontrivial matrix form-factor $\hat{M}_{l}$ in the Bethe-Salpheter Eq. (6). The matrix form-factor $\hat{M}_{l}$ is defined as the multipole moment of the two-particle matrix element

[TABLE]

Corrections to the form-factor $\hat{M}_{l}$ due to finite Cooper pair momentum $q$, $|\Delta M|=q^{2}|u|/4p_{\mathrm{F}}^{2}\epsilon_{\mathrm{F}}$, are negligible in the weak coupling regime since $q\ll p_{\mathrm{F}}$. As a result, the matrix (8) can be approximated as follows

[TABLE]

Here the coefficient $c$ and $s$ correspond to the coherence factors $c_{\mathbf{p}}$ and $s_{\mathbf{p}}$ in (3) evaluated at the average Fermi energy $\epsilon_{\mathrm{F}}$ for electrons and holes and are given by

[TABLE]

The matrix form-factor (9) shapes the sublattice structure of fluctuating electron-hole Cooper pairs in double-bilayer graphene. Note that it couples scattering channels $\hat{\Gamma}_{l^{\prime}l}$ with different orbital momenta. Importantly, $\hat{M}_{l}$ has only even harmonics $l=0,\pm 2,\pm 4$ that forbid scattering between states with even and odd orbital momenta. For isotropic Coulomb interactions, the $s$-wave moment $U_{\mathrm{s}}$ is expected to be largest. It is instructive to start by neglecting all other moments. In that case only $\hat{\Gamma}_{00}$ is nonzero. The latter depends on the $s$-wave moment of the form-factor $\hat{M}_{0}$ that has a different form for intra-valley and inter-valley Cooper pairs. We discuss these two case separately below.

III.2 B. Intra-valley Cooper pairs

For intra-valley ($v_{\mathrm{t}}=1$ and $v_{\mathrm{b}}=1$) electron-hole Cooper pairs the s-wave moment of the form-factor $\hat{M}_{0}=\langle\hat{M}(\phi_{\mathbf{p}})\rangle_{\phi_{\mathbf{p}}}$ is given by

[TABLE]

The scattering problem decouples into the three channels identified in Ref. [Neilson et al., 2014] and the corresponding scattering vertex is given by

[TABLE]

Here $L_{\alpha}=1-\lambda^{\mathrm{s}}_{\alpha}\Pi/N_{\mathrm{F}}$ is a dimensionless inverse Cooper propagator for channel $\alpha$ and $\lambda^{\mathrm{s}}_{\alpha}$ is the corresponding coupling constant specified in Tab. (1). $L_{\alpha}$ vanishes at the critical temperature $T_{\alpha}$ for the electron-hole pairing instability in channel $\alpha$.

In the absence of disorder and electron-hole density imbalances, the critical temperatures are given by $\bar{T}_{\alpha}=2e^{\mathrm{C}}\Lambda\exp[-1/\lambda^{\mathrm{s}}_{\alpha}]/\pi$, where $C=0.577$ is the Euler constant. Although the coupling constant values $\lambda_{\alpha}$ can be fine-tuned by the displacement field, as it is presented in Fig. 4-a, their hierarchy is universal for the case $|u|\ll\epsilon_{\mathrm{F}}$. The coupling constant $\lambda_{12\hbox{-}21}^{\mathrm{s}}\approx 1/2$ is almost twice as large as the constants $\lambda_{11(22)}^{\mathrm{s}}\approx 1/4$, ensuring domination of the mixed channel $12\hbox{-}21$. Physically the presence of two sub lattice combinations ($|12\rangle$ and $|21\rangle$) doubles the number of states that take part in the Cooper pairing.

III.3 C. Inter-valley Cooper pairs

For inter-valley ($v_{\mathrm{t}}=1$ and $v_{\mathrm{b}}=-1$) electron-hole Cooper pairs the s-wave moment of the form-factor $\hat{M}_{0}=\langle\hat{M}(\phi_{\mathbf{p}})\rangle_{\phi_{\mathbf{p}}}$ is given by

[TABLE]

This case also decouples into three independent channels, with scattering vertex

[TABLE]

Interestingly, the sublattice structure of the Cooper pairs and the corresponding coupling constants are different in intra-valley and inter-valley cases. The latter are presented in the Tab. (2) and their dependence on displacement field parameter $u$ is shown in Fig. 4-b. In that case the hierarchy between coupling constants is universal and does not depend on the ratio between $|u|$ and $\epsilon_{\mathrm{F}}$. The mixed channel $11\hbox{-}22$ has the highest critical temperature.

The separation of scattering problem into three channels is not an artifact of the $s$-wave truncation, but is maintained when higher multipole momenta of interactions $U_{\mathrm{l}}$ are taken into account. When the $d$-wave interaction $U_{\mathrm{d}}\equiv U_{\pm 2}$ is nonzero ($s$- and $p$-wave momenta are decoupled and the latter is irrelevant) the scattering matrix $\Gamma_{l^{\prime}l}$ is nonzero for $l=-2,0,2$. We show below that the only effect of the $d$-wave momentum $U_{\mathrm{d}}$ on the tunneling conductance between graphene bilayers is a renormalization of coupling constants $\lambda_{\alpha}=\lambda_{\alpha}^{\mathrm{s}}+\lambda_{\alpha}^{\mathrm{d}}$. The $d$-wave coupling constants $\lambda_{\alpha}^{\mathrm{d}}$ are summarized in Tabs. 1 and  2.

III.4 D. Disorder and density imbalances

Disorder and electron-hole density imbalances both reduce the critical temperature $\bar{T}_{0}$ (index [math] corresponds to the channel $\alpha$ that has the largest critical temperature of electron-hole pairing). The latter can stabilize the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) pairing state Fulde and Ferrell (1964); Larkin and Ovchinnikov (1965); Seradjeh (2012); Efimkin and Lozovik (2011); Pieri et al. (2007); Subas ı et al. (2010) with finite Cooper pair momentum $\mathbf{Q}_{0}$ (Ref. Pieri et al. (2007); Subas ı et al. (2010) also outlines a stabilization of the Sarma-phases Sarma (1963) within the BEC-BCS crossover. These phases are unstable in the weak-coupling regime considered here). The critical temperature $T_{0}$ and the instability momentum $\mathbf{Q}_{0}$ for a channel with highest critical temperature satisfy the equation $L_{0}(0,\mathbf{Q}_{0})=0$, which can be recast as follows:

[TABLE]

The dependence of the critical temperature $T_{0}$ on the pair-breaking rate $\gamma$ and the electron-hole imbalance expressed as a difference in Fermi energies, $h$, are illustrated in Fig. 5-a. The electron-hole pair instability is suppressed when the pair-breaking rate exceeds a critical value $\gamma\approx 1.78\,\bar{T}_{0}$. If the rate does not exceed $\gamma\approx T_{0}$ the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) state with finite Cooper pair momentum $\mathbf{Q}_{0}$ is stabilized by finite imbalance. The dependence of the instability momentum $\mathbf{Q}_{0}$ is presented in Fig.5-b. Its magnitude can be approximated by $v_{\mathrm{F}}Q_{0}\approx h$, corresponding to the difference between the Fermi momenta for electrons and holes.

III.5 E. Coherence time and length of fluctuating Cooper pairs

In the absence of electron-hole imbalance, the Cooper propagator $L_{0}^{-1}$ of the dominating channel at small frequencies and momenta simplifies to

[TABLE]

Here $\epsilon=\ln\left[T/T_{0}\right]$ is the energy scale that is required to create a uniform fluctuating Cooper pair. It vanishes at the critical temperature $T_{0}$ and is linear $\epsilon\approx\Delta T/T_{0}$ in its vicinity $\Delta T\approx T_{0}$. Here $\Delta T=T-T_{0}$. $\tau$ and $\xi$ are characteristic time and spatial scales for Cooper pairs that are given by

[TABLE]

They are connected to the coherence time and length of Cooper pairs by $\tau^{*}=\tau/2\epsilon$ and $\xi^{*}=\xi/\sqrt{\epsilon}$ that diverge at the critical temperature for electron-hole condensation $T_{0}$. The Cooper propagator (16) has its only pole on the imaginary frequency axis at $\omega_{\mathbf{q}}=-i/2\tau^{*}_{\mathbf{q}}$ with $\tau_{\mathbf{q}}^{*}=\tau/2\epsilon_{\mathbf{q}}$, reflecting the dissipative nature of Cooper pairs dynamics. Due to the presence of a finite temporal coherence time $\tau^{*}$, fluctuating Cooper pairs do not provide a dissipationless Josephson current, but do strongly enhance the tunneling conductance at zero voltage bias.

IV IV. Tunneling conductivity

IV.1 A. Linear response theory for the tunneling conductance

When inter-layer tunneling amplitudes are treated in leading order of perturbation theory, the inter bilayer tunneling conductance at finite voltage bias $V$ is Mahan (1990); Zheng and MacDonald (1993); Jungwirth and MacDonald (1996)

[TABLE]

Eq. (18) accounts for the fourfold degeneracy due to the presence of valley and spin degrees of freedom and $A$ is the sample area. Here $\chi(\Omega)$ is the retarded correlation function corresponding to the imaginary-time ordered correlation function constructed from tunneling operators

[TABLE]

$\chi(\omega,\mathbf{q})$ can be constructed from $\chi(\tau,\mathbf{q})$ by the usual Fourier transform and analytical continuation steps. Without Coulomb interactions between electrons and holes the response function $\chi(\omega,\mathbf{q})$ corresponds to the single electron-hole loop diagram depicted in Fig. 3-b and is given by

[TABLE]

Here $\hat{t}=\{t_{11},t_{12},t_{21},t_{22}\}$ is the vector of tunneling matrix elements in the compact representation, and $\Pi\equiv\Pi(\omega,\mathbf{q})$ is a single step pair propagator in the Cooper ladder defined in Eq. (7). In the presence of Coulomb interactions, the single-particle Green functions and the tunneling vertex need to be renormalized. The renormalization of Green functions results in a dip of the density of states at the Fermi level Rist et al. (2013), that does not produce any singularities in the tunneling conductivity and thus is unimportant and can be neglected. The renormalized tunneling vertex $\hat{t}\rightarrow\hat{t}^{\prime}\equiv\hat{t}^{\prime}(\omega,\mathbf{q})$ diverges at the critical temperature $T_{0}$ and is responsible for the drastic enhancement of the tunneling conductivity in its vicinity. The renormalization of tunneling vertex is presented in Fig. 3-d and the corresponding equation for $\hat{t}^{\prime}$ can be written as

[TABLE]

The matrix form-factor $\hat{M}$ couples even orbital channels and we neglect all multipole moments except for $s$- and $d$. (The $p$-wave multipole moment is decoupled from the $s$-wave one and is irrelevant.) Since the form-factors $\hat{M}$ are different for intra-valley and inter-valley Cooper pairs, we again consider these two cases separately.

IV.2 B. Intra-valley tunneling

Without Coulomb interactions between electrons and holes the response function $\chi(\omega,\mathbf{q})$ is a sum of three non-interfering terms that correspond to three channels $\alpha$ introduced in Sec. III and identified in Ref. [Neilson et al., 2014]:

[TABLE]

The three channels are not coupled by Coulomb interactions and are renormalized independently as follows:

[TABLE]

Here $L_{\alpha}$ is the inverse Cooper propagator for each channel and $\lambda_{\alpha}=\lambda_{\alpha}^{\mathrm{s}}+\lambda_{\alpha}^{\mathrm{d}}$ is the corresponding coupling constant. Importantly the only role of the $d$-wave interaction moment is the renormalization the coupling constant $\lambda_{\alpha}$.

IV.3 C. Inter-valley tunneling

When opposite valleys are aligned by a twist angle between bilayers close to $\theta=\pi/3$, the response function $\chi(\omega,\mathbf{q})$ for noninteracting electrons and holes is

[TABLE]

It again is a sum of three non-interfering terms that correspond to three channels $\alpha$. Coulomb interactions do not couple the channels but renormalize as follows

[TABLE]

Each channels acquires its own Cooper propagator $L_{\alpha}$ with the coupling constant $\lambda_{\alpha}=\lambda_{\alpha}^{\mathrm{s}}+\lambda_{\alpha}^{\mathrm{d}}$.

IV.4 D. Tunneling conductance

Due to the remarkable conservation of momentum for electron tunneling, the conductance $G_{\mathrm{T}}$ is observable only if graphene bilayers are aligned with twist angle $\theta=0$ or $\theta=\pi/3$. In the former case the tunneling is intra-valley, while in the latter case it is inter-valley. Therefore we refer to the cases $\theta=0$ and $\theta=\pi/3$ as to intra-valley and inter-valley alignments. As we explain below, there important differences in the enhancement of tunneling conductance by fluctuating Cooper pair in these two cases.

Within the model of non-interacting electrons and holes the tunneling conductance is dominated by tunneling between adjacent sublayers $t_{22}$ and can be approximated (both for intra- and inter-valley alignments) as follows

[TABLE]

In the presence of interactions, the enhancement of tunneling conductance by fluctuating Cooper pairs works very differently for intra- and inter-valley alignments. The reason is a drastic difference in the sublattice structure for the pairing correlations.

For intra-valley fluctuating Cooper pairs, the dominating channel $12\hbox{-}21$ does not involve pairing correlations at adjacent sublayers. As result, the corresponding contribution of the channel $12\hbox{-}21$ is weaker than that of $22$, since the latter involves tunneling between adjacent sublayers, except in the vicinity of the critical temperature $T_{0}$. The tunneling conductance is governed by the competition of two channels and can be approximated as follows

[TABLE]

As we explain in the next section, the competition between channels is essential for the explanation of the experimental data Burg et al. (2018), that in our interpretation is strongly depend on mutual orientation of graphene bilayers.

For the inter-valley alignment, the dominant channel for fluctuating Cooper pairs involves pairing correlations at adjacent sublayers. As a result, the tunneling conductance can be well approximated by a single term corresponding to the channel $11\hbox{-}22$ and is given by

[TABLE]

The inter-valley alignment case has not been studied experimentally yet, but according to our theory is most favorable for observations of the fluctuational internal Josephson effect.

IV.5 E. Critical behavior of tunneling conductance

At matched concentrations of electrons and holes the voltage dependence of tunneling conductance $G_{\mathrm{T}}$ in the vicinity of the critical temperature $T_{0}$ acquires a Lorentzian shape which is governed by the factor

[TABLE]

Here $\epsilon_{\mathbf{Q}}=\epsilon+\xi^{2}\mathbf{Q}^{2}/\hbar^{2}$ and $\epsilon=\ln(T/T_{0})$ can be interpreted as an energy of fluctuating Cooper pairs. In the absence of a valley splitting $\mathbf{Q}=0$, the amplitude of the peak has a long high temperature tail $F(0,0)=\tau/\ln^{2}[T/T_{0}]$. It diverges near the critical temperature $T_{0}$ as a function of $\Delta T=T-T_{0}$ in the critical manner as $F(0,0)\approx\tau T_{0}^{2}/\Delta T^{2}$ with the index $2$. Its width at half maximum $eV_{\mathrm{HM}}=1/\tau^{*}=2\Delta T/T_{0}\tau$ is equal to the inverse coherence time $\tau^{*}$ of fluctuating Cooper pairs and vanishes linearly at $T_{0}$. In the presence of a valley splitting that can be induced by in-plane magnetic field or relative twist, fluctuating Cooper pairs with finite momentum $\mathbf{Q}$ are probed in tunneling experiments. Temperature dependence of the peak width is modified as $eV_{HM}=\left(1+(\xi^{*}\mathbf{Q})^{2}\right)/\tau^{*}$. It has a universal form as a function of coherence time $\tau_{*}=2\tau/\epsilon$ and length $\xi^{*}=\xi/\sqrt{\epsilon}$ for fluctuating Cooper pairs that allow to extract them form the experimental data in the presence of in-plane magnetic field.

For the case of the inter-valley alignment ($\theta=\pi/3$), the temperature dependence of the tunneling conductance (28) is well approximated by Eq. (29) in a wide temperature range. For the case of intra-valley alignment ($\theta=0$) channels $12\hbox{-}21$ and $22$ compete with each other. As a result, tunneling conductance Eq. (27) can not be approximated by a simple analytical expression over as wide temperature range. It is well approximated by the critical behavior given by Eq. (29) only in the narrow temperature range where the contribution of the channel $21\hbox{-}21$ dominates.

V V. Comparison with experiment

In Ref. [Burg et al., 2018] the twist angle between graphene layers can in principle be tuned to access the intra-valley ($\theta=0$) and inter-valley ($\theta=\pi/3$) tunneling cases, although experimental results for tunneling conductance $G_{\mathrm{T}}$ have so far been reported only for one alignment. A preliminary analysis of the data suggests that a divergent zero bias peak appears on the top of a background with a weaker temperature dependence. This behavior can be explained by the competition between channels $12\hbox{-}21$ and $22$. This is not surprising since top and bottom graphene bilayers originate from the same flake, and are aligned with the twist angle $\theta=0$. Calculations for the intra-valley alignment are presented below and compared with the experimental data, while ones for the inter-valley one are presented in Appendix B.

To fit the experimental data, the tunneling conductance $G_{\mathrm{T}}$ has been calculated with the help of Eq. (27). The model has a large number of fitting parameters. To adjust their values and compare our results with the experimental data we use the following strategy that involves six steps.

1. The experimental data has already been carefully analyzed Burg et al. (2018) within a model of noninteracting charge carriers. In the wide range of concentrations the non-interacting model explains the tunneling data very well except in the case of opposite polarity charge carriers with nearly equal electron and hole densities. Based on the fits to experimental data away from matched electron and hole densities we can confidently assign values for the adjacent layer tunneling amplitude, $|t_{22}|=30\;\mu\hbox{eV}$, and the disorder-broadening energies $\gamma_{\mathrm{t}(\mathrm{b})}=4\;\hbox{meV}$. Note that the measured disorder broadening parameter $\gamma_{\mathrm{t}(\mathrm{b})}$ is much larger than $T_{0}$, where $T_{0}$ is the temperature at which the tunneling conductance appears to diverge experimentally. It is immediately clear therefore, even before performing a detailed analysis, that the condensation temperature must be strongly suppressed by disorder.

2. The colossal enhancement of the tunneling conductance has been observed in the wide density range $4\cdot 10^{10}\sim 10^{12}\hbox{cm}^{-2}$. The strength of the interactions in the system decreases with doping and we have chosen the elevated doping level $n=7.4\cdot 10^{11}\;\hbox{cm}^{-2}$ because of the detailed experimental data is available in this case. This doping level corresponds to the average concentration of electrons and holes $n=(n_{\mathrm{e}}+n_{\mathrm{h}})/2$, while the electron-hole imbalance $\Delta n\ll n$ can be present.

3. The transport experiments with similar double bilayer graphene samples and with similar gating geometry have already been performed Lee et al. (2017). The doping level of charge carriers $\epsilon_{\mathrm{F}}$ and the gap $2|u|$ in the electronic spectrum of bilayer graphene are independent and are controlled by the same gate. At the considered average density $n=7.4\cdot 10^{11}\;\hbox{cm}^{-2}$, the effective mass parameter of bilayer graphene can be approximated as $m\approx 0.04\;m_{0}$, while the gap is equal to $2|u|\approx 6.6\;\hbox{meV}$ and is much smaller than the corresponding Fermi energy $\epsilon_{\mathrm{F}}\approx 20\;\hbox{meV}$. It follows that $|u|/\epsilon_{\mathrm{F}}\approx 0.16$, which implies sublayer polarization within the bilayers $c^{2}\approx 0.58$ and $s^{2}\approx 0.42$ to be modest.

4. The bare critical temperature without disorder $\bar{T}_{0}\approx 50\;\hbox{K}$ can be recalculated from the actual critical temperature at the considered doping level $T_{0}=1.5\;\hbox{K}$ and the Cooper pair-breaking rate $\gamma=\gamma_{\mathrm{t}}+\gamma_{\mathrm{b}}=8\;\hbox{meV}$ that has been chosen above with help of Eq. (15). The system is in the regime of strong pair breaking $T_{0}\ll\bar{T}_{0}\sim\gamma$ and the value $\gamma/T_{0}\approx 1.74$ is very close to the critical value $1.78$. Experimentally, singular behavior of the tunneling conductance is observed only in the cleanest samples.

5. The bare critical temperature $T_{0}$ of Cooper pairing in the weak-coupling regime is given by $\bar{T}_{0}=2e^{\mathrm{C}}\Lambda\exp[-1/\lambda_{0}]/\pi$, where $C=0.577$ is the Euler constant. Here $\lambda_{0}=\lambda_{0}^{\mathrm{s}}+\lambda_{0}^{\mathrm{d}}$ corresponds to the mixed-channel $12\hbox{-}21$ that dominates in the considered regime $|u|\ll\epsilon_{\mathrm{F}}$. Approximating the high energy cutoff as $\epsilon_{\mathrm{c}}\approx 2\epsilon_{\mathrm{F}}$ and employing the values of $\bar{T}_{0}$ and $\epsilon_{\mathrm{F}}$ chosen above we get $\lambda_{0}\approx 0.44$. The applicability condition for the weak-coupling approach $\lambda_{0}\ll 1$ is not well satisfied, but this value for $\lambda_{0}$ corresponds to moderate coupling regime $\lambda_{0}\lesssim 1$ that justifies the approximations used in our theory.

6. The relative contribution to the tunneling conductance $G_{\mathrm{T}}$ of channels $12\hbox{-}21$ and $22$ can be characterized by a dimensionless parameter $r=|t_{12}-t_{21}|^{2}/|t_{22}|^{2}$ that we treat in the phenomenological way. Along with the sample area $A$, $r$ and $A$ are only free parameters of the model that have not yet been assigned. We adjust them by fitting the measured temperature dependence for tunneling conductance at zero voltage bias $V=0$ (and also at matched concentrations of electrons and holes and zero in-plane magnetic field), that is presented in Fig. 1, with the theory that incorporates the effect of fluctuating electron-hole Cooper pairs (Eq. 27). The area $A$ can be be obtained by matching the high temperature behavior of $G_{\mathrm{T}}$ since the contribution of $12\hbox{-}21$ in this case is negligible small. It results in $A\approx 397\;\mu\hbox{m}^{2}$. The value of $r\approx 7.6\cdot 10^{-5}$ is obtained by fitting the singular behavior of the tunneling conductance in the vicinity of the critical temperature. The corresponding theoretical curve is also presented in Fig. 1 and matches with the experimental data reasonably well over a wide temperature range.

In Fig. 1 we also present calculations within the model of noninteracting electrons and holes (Eq. (26)). This model severely underestimates the tunneling conductance $G_{\mathrm{T}}$ in the case of matched electron and hole concentrations, where interactions are crucial to explain enhanced tunneling conductance at low temperature and singular behavior in the vicinity of the critical temperature $T_{0}$ for Cooper pair condensation. We will keep all parameters chosen above in further calculations and investigate an impact of finite voltage bias between bilayers $V$, electron-hole imbalance $\Delta n$, and in-plane magnetic field $B$ at tunneling conductance. We will also present only results of the theory that incorporates the effect of fluctuating Cooper pairs.

A comparison between theory and experiment for the voltage dependence of the tunneling conductance between graphene bilayers $G_{\mathrm{T}}$ is presented in Fig. 6. The zero bias peaks emerge with a decreasing temperature on a top of a smooth background that corresponds to the channel $22$. The width of the background $eV_{\mathrm{HM}}$ is governed by the single-particle energy scales $2\pi T$ and $\gamma$ and is approximately equal to the largest of them. The width of the zero bias peak in the vicinity of $T_{0}$ is much smaller than the single-particle disorder scale, that demonstrates its collective origin. While the temperature dependence of peak height is well fit by the theory, the width dependence is captured only qualitatively and is overestimates by the factor of 2. The experimental data also exhibit a voltage asymmetry that becomes more prominent at low temperatures. Within our phenomenological model an asymmetrical voltage dependence of the tunneling conductance can be obtained if the scattering rates of electrons and holes $\gamma_{\mathrm{t}(\mathrm{b})}$ that define Cooper pair scattering rate as $\gamma=\gamma_{\mathrm{t}}+\gamma_{\mathrm{b}}$ are energy-dependent. The simple linear dependence $\gamma=\gamma+\gamma^{\prime}\omega$ with a phenomenological parameter $\gamma^{\prime}$ does not capture the observed asymmetry however. A quantitative understanding of the asymmetry requires a microscopic understanding of disorder mechanisms that is outside the scope of the present work.

The comparison between theory and experiment for the dependence of the tunneling conductance $G_{\mathrm{T}}$ on the electron-hole imbalance $\Delta n$ and the voltage bias $V$ are summarized in Fig. 7 and a zero-bias-voltage cut of the comparison made in this color plot at zero voltage bias is presented in Fig. 8. The theoretical curves again agree reasonably well with the data. A density imbalance splits the Fermi lines of the electrons and holes and disfavors their Cooper pairing. As seen in Fig. 5, the FFLO state with a finite Cooper pair momentum can not be stabilized in the strong pair breaking regime realized in the experiment. The critical temperature $T_{0}$ is maximal for Cooper pairs with zero Cooper pair momentum and decreases monotonically in the presence of imbalance and vanishes if the latter exceeds the critical temperature. As a result, the dependence of tunneling conductance on imbalance $\Delta n$ and in-plane magnetic filed $B$ that is presented in Fig. 9 is smooth and featureless. Fig. 9 presents theoretical curves since the experimental data for magnetic field and density-balance dependence at this temperature is not yet available. The theory qualitatively explains the decrease of peak height with magnetic field studied experimentally at lower temperatures $T_{0}\approx 1.5\;\hbox{K}$, but considerably overestimates its effect. In this case the system is in the paired state whose behavior lies outside the range of validity of the present theory. Other less fundamental limitations might also explain this discrepancy as we discuss in more details in the next section. We conclude, that our theory of the fluctuational internal Josephson effect, combined with specific features of the multiple-channel structure of pairing in double bilayer graphene provides a reasonable overall description of experiment.

VI VI. Discussions

Our theory of the internal fluctuational Josephson effect does not account for interactions between fluctuating Cooper pairs. The Gaussian nature of the theory we employ is more clearly seen within an alternate derivation of the tunneling conductance. The latter employs the auxiliary field approach and is presented in Appendix C. Interactions between fluctuating pairs can be safely omitted in the wide range of temperatures $\Delta T\sim T_{0}$, and are important only in the critical region $\Delta T_{Gi}\lesssim Gi\,T_{0}$ where fluctuations are large and strongly interfere with each other. Here $Gi=T_{0}/E_{F}$ is the Ginzburg number calculated in Appendix D. The double bilayer graphene system studied experimentally Burg et al. (2018) is close to the weak coupling regime, and the critical region $\Delta T_{\mathrm{Gi}}=T_{0}^{2}/E_{\mathrm{F}}\approx 10\;\hbox{mK}$ is much smaller than the temperature range $\Delta T\approx 4\;\hbox{K}$ where the tunneling conductance is strongly enhanced. The picture of noninteracting fluctuating Cooper pairs is therefore well justified to address the basic phenomenon identified in experiment.

Fluctuating Cooper pairs in conventional superconductors alter the thermodynamics of the normal state only in the critical region $\Delta T_{Gi}.$ The conductivity and the magnetic susceptibility Varlamov et al. (2018) are however not only singular at the critical temperature, but have long high temperature tails Skocpol and Tinkham (1975). The high temperature tail for the diamagnetic susceptibility $\chi\sim\chi_{\mathrm{L}}/\ln^{2}[T/T_{0}]$, where $\chi_{\mathrm{L}}$ is the Landau diamagnetic susceptibility in the normal state, was predicted Aslamazov and A.I. (1975); Bulaevskii (1974); Maki (1973) theoretically and observed in experiments Skocpol and Tinkham (1975). The reason is the paired state is superconducting and mediates the perfect diamagnetism that makes even a small number of fluctuating Cooper pairs important. Similarly the equilibrium paired state of spatially separated electrons and holes provides in fluctuational internal Josephson effect that colossally enhances the inter-layer tunneling. That is why even a small density of fluctuating Cooper pairs can make a strong impact on the tunneling conductance above the critical temperature $T_{0}$ which also has the high-temperature tail proportional to $1/\ln^{2}[T/T_{0}]$ that is clearly seen in Fig. 1.

In the vicinity of the critical temperature $T_{0}$ zero bias peak shape is universal and governed by the factor $F(eV,\mathbf{Q})$ given by Eq. (29). It can be rewritten as $F(eV,\mathbf{Q})=\mathrm{Im}[L_{0}^{-1}(eV,\mathbf{Q})]/eV$. Thus an imaginary part of the Cooper propagator $L_{0}^{-1}(eV,\mathbf{Q})$, that can be interpreted as Cooper pair susceptibility Com (d), is directly probed in tunneling experiments Com (e). It should be noted that Cooper pair-susceptibility of a superconductor can also be probed in tunneling Josephson junction in which one side is near its critical temperature while the other is well below its critical temperature Ferrell (1969); Scalapino (1970). The junction does not support dissipationless Josephson tunneling current, but the tunneling current at finite voltage is strongly enhanced by fluctuating Cooper pairs that grow in the vicinity of $T_{0}$. The latter has been observed experimentally Anderson et al. (1972); Aspen and Goldman (1981).

The enhancement of intra-valley tunneling that corresponds to the alignment ($\theta=0$) has only been reported so far. In this case the instability happens in one channel, while the main contribution to the tunneling conductance comes form the different one. As a result the divergent contribution to the conductance due to fluctuating Cooper pairs appears on the top of background with weak temperature dependence that dominates at higher temperatures. We predict an enhancement of inter-valley tunneling ($\theta=\pi/3$) to be much more profound because Cooper pairs in the dominating pairing channel involve electrons and holes settled at adjacent sublayers. We present calculations for this case in Appendix B.

The Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) state with finite Cooper pair momentum has been predicted in conventional superconductors more than sixty years ago Fulde and Ferrell (1964); Larkin and Ovchinnikov (1965). It requires the splitting of Fermi surfaces/lines for pairing electrons with opposite spins. The splitting can be induced by a magnetic field provided that its paramagnetic effect is larger than its diamagnetic one (Chandrasekhar-Clogston limit  Chandrasekhar (1962); Clogston (1962)). This condition is rarely satisfied even in layered conventional superconductors subjected to in-plane magnetic field. There are few observations in heavy-fermion and organic superconductors where FFLO state signatures have been claimed but are still debated (See Ref. Matsuda and Shimahara (2007) and Refs. Kinnunen et al. (2018); Gubbels and Stoof (2013) for reviews of progress in solid state and cold atoms systems). So far the FFLO state has not been unambiguously identified. In double bilayer graphene the densities of electrons and holes can be controlled separately in a way that opens the FFLO state up for experimental study in a condensed matter system Efimkin and Lozovik (2011); Seradjeh (2012). The FFLO can be unambiguously identified if it appears from the dependence of the zero-bias peak on imbalance and in-plane magnetic field, since the latter makes it possible to probe Cooper pairs with finite momentum. In the vicinity of the instability to the uniform paired state, the tunneling conductance monotonically decreases with in-plane magnetic field (as presented in Fig. 9). In the vicinity of an instability to the FFLO state the tunneling conductance achieves a maximum at finite field-induced momentum shift $\mathbf{Q}_{\mathrm{B}}=\mathbf{Q}_{0}$ where $\mathbf{Q}_{0}$ is the corresponding momentum of Cooper pairs. We discuss how distinguish these states in more details in Appendix B, where calculations for the inter-valley tunneling are presented.

The sensitivity of Cooper pairing to a disorder opens a possibility of a granular electron-hole state in the presence of its strong long-range variations. In this state the pairing happens in disconnected or weakly coupled regions with minimal amount of disorder and does not support the spatial coherence. It makes the transport properties of the system including Coulomb drag effect to be different from ones in the uniform paired state. A tunneling conductance in the granular state is still colossally enhanced since the latter requires temporal coherence of Cooper pairs but not the spatial one.

The interpretation of experiment provided by our theory suggested that the pairing critical temperature would be substantial if samples with weaker disorder could be fabricated. This finding is perhaps a bit surprising since the experiments are for the most part conducted in the weak to moderate coupling regime ($r_{\mathrm{s}}=0.72\sim 3.3$) where some researchers have argued that critical temperatures should be strongly suppressed by screening, especially accounting for spin and valley degeneracy Kharitonov and Efetov (2008); Lozovik et al. (2012) (See also arguments that this approach considerably underestimates the critical temperature Sodemann et al. (2012); Mink et al. (2011); Perali et al. (2013); Conti et al. (2017)). Our theory also suggests that high pairing temperatures should be achievable in double single-layer graphene systems, since there is nothing in its structure that puts single layers at a disadvantage relative to bilayer. Moreover, for the linear spectrum in monolayer graphene the Wigner-Seitz radius is defined in a different way $r_{\mathrm{s}}=2.19/\kappa$ and could achieve even larger value $r_{s}=1.1$ if hBN is used as a spacer between graphene sheets. Here $\kappa$ is the corresponding dielectric constant. Future experimental work which seeks to weaken pair-breaking by disorder and which explores double single-layer graphene systems as well, is therefore important.

In summary, the theory of the fluctuational internal Josephson effect developed here explains the anomalies in the tunneling conductance between graphene bilayers observed experimentally at equal electron and hole densities, including their dependence on temperature, bias voltage bias and electron-hole imbalance. Some aspects of the observations are nevertheless not understood. First of all the observed dependence of the tunneling conductance on bias voltage has an asymmetry between positive and negative bias that becomes more prominent with decreasing temperature. At first glance the asymmetry is unexpected and surprising since the electronic spectrum of two graphene bilayers with matched concentrations of electrons and holes is symmetric, as it is clearly seen in Fig. 2. The symmetry can be broken by Coulomb impurities if most of them are of the same charge. For example positive charges (ionized donors) provide repulsive scattering for holes and attractive scattering for electrons. Our model takes the scattering rates for electrons and holes $\gamma_{\mathrm{t}(\mathrm{b})}$ to be momentum and energy independent and ignores these common complications. An asymmetry can be introduced can be introduced in phenomenological way by making the assumption that Cooper pair-breaking time is energy dependent $\gamma(\omega)$. Approximating it by a linear function does predict an asymmetry of the tunneling conductance that grows with decreasing temperature, but the shape of the experimental curves is not captured by this simple ansatz. To clarify whether or not the observed asymmetry can be explained by the presence of charge impurities, a more microscopic description of their scattering characteristics is needed and this is outside of the scope of the present work. Secondly it is not clear whether or not our theory can capture the dependence of tunneling conductance on magnetic field since more experimental data is needed. Comparison with data obtained at $T\approx 1.5\;\hbox{K}$ suggests that the theory considerably overestimates the effect of magnetic field. Nevertheless, at such temperatures the system is in the paired state or in the critical regime that is outside of the applicability range of the theory of Gaussian fluctuations. This discrepancy can be due to other reasons. Bilayer graphene as other two-dimensional systems have long-range density variations and if the corresponding length is smaller than $\hbar/Q_{\mathrm{B}}$ the effect of the magnetic can not be reduced just to the relative shift of dispersion for electrons from different layers. To better understand capabilities of the theory more experimental data is needed.

To conclude, we have developed a theory of the fluctuational internal Josephson effect in the system of two closely spaced graphene bilayers. The presence of valley and sublattice degrees of freedom provides three competing electron-hole channels for both intra-valley and inter-valley Cooper pairs. We show that three channels are nearly independent and have different critical temperatures of the condensation and sublattice structures. The observed enhancement of the tunneling conductance can be explained only by the presence of competing channels that dominate in different temperature ranges. The theory reasonably captures the dependence of the conductance on temperature, voltage bias between bilayers and electron-hole imbalance. We also argue that the enhancement is much stronger for inter-valley tunneling than for the intra-valley one that has been reported recently. We also discuss how to distinguish the uniform state and the FFLO state with finite Cooper pair momentum that can be stabilized in the system by an electron-hole imbalance.

Acknowledgements.

This work was supported by the Army Research Office under Award W911NF-17-1-0312, by ARO MURI 3004628717, and by the Welch Foundation under grant F-1473. DKE acknowledges support from the Australian Research Council Centre of Excellence in Future Low-Energy Electronics Technologies (FLEET).

VII Appendixes

Appendix A A. Bethe-Salpether equation

Here we present a detailed derivation of Eq. (6) from the paper. The Bethe-Saltpeter equation that is illustrated in Fig. 3 can be written as follows

[TABLE]

Here $\mathbf{p}_{\pm}=\mathbf{p}\pm\mathbf{q}/2$; $\omega_{n}=(2n+1)\pi/T$ and $p_{n}=2n\pi/T$ are fermionic and bosonic Matsubara frequencies. The electron Green function $\hat{G}^{\mathrm{t}(\mathrm{b})}(\i\omega_{n},\mathbf{p})$ in the sublattice space can be presented as follows

[TABLE]

Here $\gamma_{\mathrm{t}(\mathrm{b})}$ are scattering rates for electrons (holes). We have also neglected the presence of the valence band in the layer with excess of electrons and the presence of the conduction band in the layer with the excess of holes. The product of Green functions that appears in (30) can be written with as follows

[TABLE]

Here the matrix form-factor

[TABLE]

reflects the chiral nature of charge carriers in bilayer graphene, while $C(ip_{n},\mathbf{q},i\omega_{n},\mathbf{p})$ contains information only about energy spectrum for electrons and holes and is given by

[TABLE]

In the weak coupling regime that we consider in the paper pairing correlations do appear in the vicinity of Fermi lines for electrons and holes. As a result the vertex $\hat{\Gamma}(ip_{n},\phi_{\mathbf{p}^{\prime}},\phi_{\mathbf{p}})$, the Fourier transform of interactions $U(\phi_{\mathbf{p}^{\prime}}-\phi_{\mathbf{p}})$, and the form-factor $\hat{M}(\phi_{\mathbf{p}})$ can be safely approximated by their values at the Fermi level ($|\mathbf{p}|=p_{\mathrm{F}}$ and $|\mathbf{p}^{\prime}|=p_{\mathrm{F}}$) and depend only on the corresponding polar angles ($\phi_{\mathbf{p}}$ and $\phi_{\mathbf{p}^{\prime}}$). After decomposition over multipole momenta $\hat{\Gamma}_{l^{\prime},l}$, $U_{l}$ and $\hat{M}_{l}$ the equation (30) becomes algebraic and can be presented in a compact form

[TABLE]

Here $\Pi_{l}(ip_{n},\mathbf{q})$ is given by

[TABLE]

Its $s$-wave component $\Pi_{0}(ip_{n},\mathbf{q})\equiv\Pi(ip_{n},\mathbf{q})$ has the ultraviolet logarithmic divergence that is usually present in the weak coupling pairing theories and coincides with the single-step pair propagator in the Cooper ladder sum of a bilayer system without sublattice degrees of freedom. Its detailed derivation can be found in textbooks Varlamov et al. (2018); Larkin and Varlamov (2005) and its explicit expression is presented in the main part as Eq. (7). Therefor all information about the chiral nature of the bilayer graphene charge carriers is hidden in the nontrivial matrix form-factor $\hat{M}_{l}$. At finite $l$ the value $\Pi_{l}(ip_{n},\mathbf{q})$ is nonzero only at finite Cooper pair momentum and are much smaller than $\Pi_{0}(ip_{n},\mathbf{q})$ and can be neglected. As a result, the Eq. (33) reduces to Eq. (6) from the main text of the paper.

Appendix B B. Inter-valley alignment and identification of the FFLO state

In the main text of the paper calculations for intra-valley tunneling ($\theta=0$) are discussed, while ones for the inter-valley alignment ($\theta=\pi/3$) are presented here. The tunneling conductance is dominated by the channel $11\hbox{-}22$ because of fluctuating Cooper pairs in the dominating channel involve correlations at adjacent sublayers. The tunneling conductance can be approximated as Eq. (28). We use the same set of parameters that has been used above to fit the experimental data except for pair-breaking rate $\gamma$. For the latter we use $\gamma=2,\,4,\,\hbox{and}\;8\;\hbox{meV}$. The first two values correspond to cleaner samples compared to ones that have been studied experimentally Burg et al. (2018). The temperature dependence of tunneling conductance is presented in Fig. 10. Its enhancement of the conductance is considerably stronger than that for the inter-valley tunneling. The pair-breaking rate $\gamma$ determines the critical temperature of pair condensation, but weakly influences the temperature dependence of tunneling conductance above it.

For the pair-breaking rate $\gamma=2\;\hbox{meV}$ the critical temperature is $T_{0}\approx 42\;\hbox{K}$. According to the phase diagram Fig. 5 the ratio $\gamma/T_{0}\approx 0.52$ is small enough to stabilize the Fulde-Ferrell-Larkin-Ovchinnikov (FFLO) state by the electron-hole imbalance. In-plane magnetic filed results in a relative shift of electronic dispersions between layers and makes it possible to probe fluctuating Cooper pairs with finite momentum $Q_{\mathrm{B}}=edB_{||}/\hbar c$. The dependence of tunneling conductance at zero voltage bias on the electron-hole imbalance and magnetic field at temperature $T=55\;\hbox{K}$ is presented in Fig. 11-a. The tunneling conductance has a monotonic dependence demonstrating that the system is far away from the instability to the FFLO state and Cooper pair fluctuations with zero momentum $\mathbf{Q}=0$ are the most intensive. The dependence at $T=25\;\hbox{K}$ is presented in Fig. 11-b with a dashed line that denotes a phase boundary of a paired state. Since the critical imbalance required to suppress the pairing instability decreases with magnetic field the paired state is the uniform BCS one. Nevertheless, the tunneling conductance achieves maximum at finite value of magnetic field that shows that the systems is close to the FFLO state and fluctuations with finite Cooper momentum are the most intensive. The dependence at $T=5\;\hbox{K}$ is presented in Fig. 11-c. The dashed line that denotes the phase boundary is non-monotonous and has a kink at finite value of magnetic field. It clearly demonstrates that the equilibrium FFLO state is stabilized by the electron-hole imbalance.

Appendix C C. Bosonic picture of the fluctuational internal Josephson effect

Here we present a bosonic picture of the fluctuational internal Josephson effect that can be employed with a help of the field integral formalism. We consider the contact interactions $U$ between electrons and holes that correspond to the truncation of all multipole momenta of interactions except the $s$-wave one. In the main text of the paper we have demonstrated that they are unimportant. The correlation function of tunneling operators $\chi(\omega,\mathbf{q})$ that defines the tunneling conductance $G_{\mathrm{T}}$ according to Eq. (18) can be extracted from the corresponding imaginary time correlation function $\chi(\tau,\mathbf{r})$ given by

[TABLE]

Here $Z[\Lambda^{+}_{\tau,\mathbf{r}},\Lambda_{\tau,\mathbf{r}}]$ is the statistical sum with an auxiliary bosonic field $\Lambda_{\tau,\mathbf{r}}$ introduced to the action $S$ of the system as follows

[TABLE]

Here $\hat{\psi}_{\mathrm{t}}\equiv\hat{\psi}_{\mathrm{t}\tau\mathbf{r}}$ and $\hat{\psi}_{\mathrm{b}}\equiv\hat{\psi}_{\mathrm{b}\tau\mathbf{r}}$ are spinor fermionic fields for electrons from top (t) and bottom (b) layers with labeling described in Sec. IIII. If they are integrated out and the corresponding action is expanded in the lowest order in tunneling matrix elements $\hat{t}$, the Fourier transform $\chi_{q}$ do appears in the action as follows $S=-\sum_{q}\chi_{q}|\Lambda_{q}|^{2}$. Here $q=\left\{q_{n},\mathbf{q}\right\}$ with bosonic Matsubara frequency $q_{n}=2\pi Tn$. For noninteracting electrons and holes calculations are straightforward and result in

[TABLE]

After analytical continuation we get the Eq. (20) from the main text. In the case of interacting electrons and holes it is instructive start with the Hubbard-Stratonovich transformation. It eliminates interactions but introduces the bosonic field $\hat{\Delta}\equiv\hat{\Delta}_{\tau\mathbf{r}}$ corresponding to electron-hole Cooper pairs. The action $S$ is modified as follows

[TABLE]

where $\hat{\Delta}^{\prime}_{\tau\mathbf{r}}=\hat{\Delta}_{\tau\mathbf{r}}+\hat{t}\,\Lambda_{\tau\mathbf{r}}$. Above the critical temperature $T_{0}$ of the electron-hole pairing the saddle point of the action is trivial $\langle\hat{\Delta}\rangle=0$ and the field $\hat{\Delta}$ corresponds to Cooper pair fluctuations Com (f). In the wide temperature range $\Delta T_{Gi}\lesssim\Delta T\sim T_{0}$ outside the critical regime $\Delta T\lesssim\Delta T_{Gi}$ fluctuations can be approximated by the noniteracting Gaussian theory. Here $\Delta T_{Gi}=Gi\;T_{0}$ with Ginzburg number $Gi=T_{0}/E_{\mathrm{F}}$ calculated in Appendix D. Integrating out fermions and expanding the action up to the second order in the bosonic field $\hat{\Delta}^{\prime}$ results in

[TABLE]

Here all matrices are in the compact representation, and $\hat{\Gamma}_{q}$ is the scattering vertex calculated in the Sec. III of the paper. Some of its components vanish at the critical temperature $T_{0}$ of the electron-hole Cooper pairing. The last term corresponds to the response $\chi^{0}_{q}$ function for noninteracting electrons and holes that is given by Eq.( 35). The action represents the bosonic picture of the Josephson effect and is valid outside the weak coupling regime. The action (37) is quadratic in the bosonic field $\Delta_{q}$ and after its integration we get the tunneling response function $\chi_{q}$ to be given by

[TABLE]

With a help of Eqs. (12) and (11) we get the tunneling response function for the intra-valley tunneling (23). In the same way with a help of Eqs. (14) and (13) we get the response function for the inter-valley one (25).

Appendix D D. Ginzburg criterion

The developed theory of the fluctuational internal Josephson effect implies that Cooper pair fluctuations are Gaussian and do not interact with each other. The interactions can be safely omitted in the wide range of temperatures $\Delta T\sim T_{0}$ except the critical region $\Delta T\lesssim\Delta T_{\mathrm{Gi}}$ where fluctuations are overgrown are strongly interfere with each other. The range $\Delta T_{\mathrm{Gi}}$ can be estimated from the Ginzburg criterion Ginzburg (1960) that compares the contribution of Gaussian fluctuations to the heat capacity $C_{\mathrm{FL}}=T_{0}/4\pi\xi^{2}\Delta T$ with predictions of the mean-field theory below the critical temperature $C_{\mathrm{MF}}=N_{\mathrm{F}}/gT_{0}=\epsilon_{\mathrm{F}}/4\pi\xi^{2}T_{0}$. Here $g=4\xi^{2}/(\hbar v_{\mathrm{F}})^{2}$ is the strength of contact interactions between fluctuating Cooper pairs neglected so far. The contribution of fluctuations $C_{\mathrm{MF}}$ grows with decreasing of temperatures and dominates in the temperature range $\Delta T_{Gi}=GiT_{0}$ with Ginzburg number $Gi=T_{0}/E_{\mathrm{F}}$. It does not depend explicitly on the pair-breaking rate $\gamma$ but only on the critical temperature $T_{0}$. It should be noted that the Ginzburg criterion can be derived microscopically in a more strict way by explicit analysis of the role of interactions between fluctuating Cooper pairs  Dogson (2001).