Effects of dynamical screening on the BCS-BEC crossover in double bilayer graphene: Density functional theory for exciton bilayers

We derive a gap equation for bilayer excitonic systems based on density functional theory and benchmark our results against quantum Monte-Carlo simulations and recent experiments on double bilayer graphene. The gap equation has a mean-field form but includes a consistent treatment of dynamical screening. We show that the gap survives at much higher densities than previously thought from mean-field estimates which gives strong indications that the double-bilayer graphene systems at zero magnetic field can be used as model systems to investigate the BCS-BEC crossover. Furthermore, we show that Josephson-like transfer of pairs can be substantial for small band gaps and densities.

In a recent experiment [1] on two bilayer graphene sheets separated by a WSe 2 insulating barrier an enhanced tunneling rate was measured, which indicates the formation of an exciton condensate. While exciton formation has previously been achieved in similar systems in the quantum Hall regime at high magnetic fields [2,3] these were the first experimental indications of exciton formation in bilayer systems at zero magnetic field. In bilayer graphene systems the electron and hole concentrations are tunable by metallic gates, which opens up for the possibility to tune the device from the Bose-Einstein condensate (BEC) to the Bardeen-Cooper-Schrieffer (BCS) regime with possible applications in novel electronic devices [4].
Theoretically, bilayer excitonic systems are typically described within mean-field theory [5][6][7][8][9][10] or using Quantum Monte-Carlo (QMC) simulations [11][12][13]. While the latter is very accurate the calculations are restricted to simplified single-band systems which partly limit the predictive power. The mean-field calculations, on the other hand, can be applied to realistic systems but neglect dynamical screening as well as higher order self-energy diagrams and are therefore, as we will argue below, only accurate in the low-density limit. The role of dynamical screening for single bilayer graphene systems and topological insulators was previously investigated in Ref. [14] by solving the imaginary axis gap equation. Similar to the present work it was shown that a static mean-field description severely underestimates the gap in a large parameter regime. In this letter we derive a new gap equation, based on density functional theory, which yields an effective account of dynamical screening at a sufficiently low computational cost to be applied to realistic materials. The method is applied to double bilayer graphene, which is the most promising bilayer candidate for achieving exciton condensation without an applied magnetic field. To our knowledge this is the first time that the effect of dynamical screening is studied for realistic multiband models of this material. We begin by giving a brief account of the derivation (the complete derivation is given in the supplemental material (SM)), then we bench-mark our method against QMC simulations for simplified models [13] and finally consider realistic systems where we compare to the experimental results in Ref. [1] and mean-field calculations in Ref. [10]. We show that, contrary to mean-field theories, the results from our effective gap equation qualitatively agree both with QMC simulations and experimental measurements. Due to its low computational cost the theory can provide an efficient means to predict exciton condensation in novel materials and an important complement to existing methods.
We start with the Kohn-Sham Hamiltonian for the electron-hole system H KS = drψ † (r)H e (r)ψ(r) + drφ † (r)H h (r)φ(r) − drdr ∆(r, r )ψ † (r)φ † (r ) + h.c. (1) Here ψ(r) (φ(r)) are the electron (hole) field operators and ∆(r, r ) is the Kohn-Sham electron-hole potential. Assuming band diagonal phase-coherent singlet coupling (i.e. electrons of spin σ, ,momentum k and band γ couple to holes with spin −σ momentum −k and band γ) the Hamiltonian can be rewritten as Here ε e/h kγ are the eigenenergies of H e/h ,b † andâ † are the corresponding creation operators and ∆ kγ is the matrix element of ∆(r, r ) in the electron and hole eigenstates (see SM). This approximation is the analogue of the decoupling approximation in superconductor density functional theory (scDFT). For graphene bilayers in ABstacking a tight-binding fit to the low-energy band structure is given by [10,15,16] 2 where Γ k = t 2 1 + (2 vk) 2 + (2E g vk) 2 /t 2 1 , Unless otherwise specified we use the bilayer bandstructure above with the intercell distance a = 0.246nm, intralayer hopping t 0 = 3.16eV and interlayer hopping t 1 = 0.38eV [10,17] for both the electron and hole layers. For simplicity we will assume ε e kγ = ε h −kγ = ε kγ in the derivation below, although it is straight-forward to generalize the derivation to different electron and hole dispersions. All calculations were performed at a temperature of 1.5 K. The Hamiltonian in Eq. 2 is diagonalized by a Bogoliubov transformation The amplitudes u kγ and v kγ are subject to the equations Utilizing the the decoupling approximation and Fourier transforming to frequency space yields Using the Sham-Schlüter connection [18][19][20] with the same approximations as in scDFT [21][22][23] one can then derive the gap equation for the Kohn-Sham system where Σ a is the exchange-correlation part of the anomalous self-energy and we have assumed that the effect of the normal self-energies are included in the bare dispersion. A generalized gap equation which include also the electron-electron and hole-hole correlations is given in the SM. It is straight-forward to show that a static exchange approximation to the self-energy yields the usual meanfield gap equation The geometrical form-factor F is related to the overlap of the single particle states |Ψ kγ ≡ |f kγ e ikr as F γγ kk = f k,γ |f k ,γ f k ,γ |f k,γ [24] (See the SM for details).
In order to include the effect of dynamical screening we compute Σ a in Eq. 15 within the GW -approximation The screened electron-hole interaction for the bilayer system is given by [5,10,14] where the k and Matsubara indexes are implicit for all quantities and Unless otherwise specified we use the dielectric constant = 2 for bilayer graphene sheets in a few layers of hexagonal boron nitride [10,25]. Π n = Π (n) q (iν n ) and Π a = Π (a) q (iν n ) are the normal and anomalous polarizations which are computed within the random-phase approximation (See SM).
In order to have an efficient implementation at low temperatures all Matsubara sums have to be done analytically. This is straightforward for the polarizations, but for the self-energy in Eq. 17 it is not as clear. In this work we adopt a similar strategy as Ref. 26 and fit W eh k (ν l ) to a multiple plasmon pole model We use up to 100 poles in order to get the correct shape of the peaks in ImW eh . The position of the main pole and the relative weight of the poles are determined by the analytical continuation of W eh (iν n ) to the real axis, while the overall weight of the poles is determined by a least square fit of the imaginary axis data. We evaluated all Matsubara sums in Eqs. 15 and 17 analytically using the MatsubaraSum package for Mathematica [27]. The details of the fitting procedure are given in the SM.
To benchmark our method we first consider a model with single parabolic electron and hole bands and compare with the Quantum Monte-Carlo simulations by Rìos et. al [13]. In Fig. 1 we compare the results from solving the full DFT gap equation (Eq. 15), the standard mean-field gap equation (Eq. 16) with the static value of the screened interaction (screened MF-equation) and the same equation but with the unscreened value of the interaction (bare MF-equation) with the corresponding results from the QMC simulations in Ref. 13. The screened MF equation was benchmarked against QMC previously, with encouraging results [8]. However, in these calculations a different form of the screened interaction was used and the formfactor was ignored in the gap equation. As we will see below, a correct treatment of the interaction and formfactor yields a larger screening and substantially worsens the performance of the screened MF equation. We begin by considering the maximum value of the gap. To facilitate the comparison we adopt the scaled units of Ref. 13 with the effective energy unit Ha * = m * e /m e κ −2 Ha and distances in units of a * 0 = κ(m e /m * e )a 0 . The density n is defined through r s /a * 0 where r s = 1/(πn). We use the same interlayer distance d/a * 0 = 0.4. In Ref. 13 the gap was estimated by fitting the excitation energies E(k) obtained from the QMC simulations to an effective mean-field gap-equation with k-independent gap E(k) = (k 2 /2m * − µ) 2 + ∆ 2 with m * , µ and ∆ as fitting variables. In our simulations we use the effective mass which gave the optimal fit for each density (extracted from Fig. 3 in Ref 13) and consider the maximum value of the k-dependent gap.
At low densities (high r s /a * 0 ) all mean-field calculations yield the same gap (left panel in Fig. 1). As the density is increased the screening reduces the gap of the screened MF result compared to the bare MF results substantially and at r s /a * 0 ≈ 4.5 the screened meanfield gap rapidly goes to zero. The additional correction in the DFT gap equation reduces the screening and increases the gap compared to the static MF results, while it is always smaller than the corresponding value in the bare MF calculation. Hence the DFT gap equation incorporates the frequency dependence of the interaction in an effective static theory, similar to the Z-factor approach for Hubbard models [28]. Due to the decrease of the effective screening the full gap goes to zero around r s /a * 0 ≈ 1.2 − 2, which corresponds to much higher densities than the static MF gap and is in qualitative agree-ment with the QMC results. At low and intermediate densities (r s /a * 0 > 4) all mean-field calculations underestimate the gap compared to the QMC data. Since the bare MF calculations sets the upper limit of ∆ obtainable by any mean-field approximation relying on the first-order diagram of Σ, this indicates that the differences at low density may originate from many-body effects beyond a first-order mean-field description. Both the screened mean-field and DFT gap equation are restricted to the first-order term of the self-energy expansion in W eh . This is a good approximation as long as W eh is sufficiently small or equivalently the screening is large, which explains why the DFT gap-equation agrees better with the QMC data for high densities than for intermediate and low densities. Thus, the poor performance of the screened MF equation for high densities is related to the neglect of the frequency dependence of the interaction, which is to a large extent cured by the effective treatment in the DFT gap-equation. The mean-field treatments give reasonable results also in the low-density limit, which is expected since the mean-field treatment becomes exact in the limit n → 0. It should be noted that the DFT gap-equation (Eq. 15) is not restricted to the GW -approximation for the self-energy, and hence in principle higher order self-energy terms can be accounted for within our framework.
In the right panel of Fig 1 we compare the condensate fractions c. In our mean-field results the condensate fractions can be estimated from the Bogoliubov amplitudes as while in the QMC calculations it was estimated by interpolating the two-body rotationally averaged density matrix. The transition to the BCS phase is characterized by a value of c < 0.2. Contrary to the QMC calculations the static MF simulations predicts c to vanish at r s /a * 0 ≈ 4.5, before the system has entered the BCS phase. The bare MF calculation on the other hand predicts a too high value of c for small r s /a * 0 . The full calculations agree qualitatively well with the QMC results, especially for small values of r s /a * 0 where both calculations predicts a transition to the BCS phase at around r s /a * 0 ≈ 2. For low and intermediate densities (intermediate-high r s /a * 0 ) the full gap equation overestimates c compared to the QMC data and predicts the crossover to the BEC phase with c ≈ 0.8 at r s /a * 0 ≈ 4, which can be compared to a transition at r s /a * 0 ≈ 7 in the QMC simulations. The overestimation of c at intermediate and low densities is inherent to all mean-field approximations and is again related to the neglect of higher order self-energy terms. However, the inclusion of the dynamical screening in the full gap equation improves the static MF results substantially and yields results which at least qualitatively agree with the QMC simulations in Ref. 13. It should be noted that in more realistic multiband models the screening is larger which, according to the discussion above, indicates that the DFT gap equation is accurate for a larger density range for these systems than for the single-band model discussed above. So far we have benchmarked our method against QMC simulations for simple single-band models. In the remainder of this letter we will focus on more realistic models for double bilayer graphene with the single-particle dispersion relation in Eq. 3. These systems were extensively studied in Ref. 10 within a static mean-field approximation for band gaps ≤ 210 meV, which is close to the maximum gap obtainable for these systems (≈ 300 meV) [16,[29][30][31][32]. In Fig. 2 we show the superfluid gaps (∆ ± ) for the conduction band (γ = +1) and valence band (γ = −1) channels. The strength of the interband screening depends intrinsically on the bandgap and the larger the band gap the smaller the screening. Therefore ∆ is larger and survives at higher densities for large band gaps. However, for large band gaps the condensate fraction ( Fig. 3) also decreases more slowly as function of the density. Due to these counteracting tendencies the static mean-field approximation predicts that ∆ vanishes in the crossover region with 0.2 < c + < 0.8 and that no system, independent of the band gap, reaches the BCS state at c + ≈ 0.2. The full calculations, based on the DFT gap equation, yield a fundamentally different conclusion. In these calculations the effect of screening is reduced substantially which yields larger superfluid gaps that survive at higher densities. The picture of a relatively constant ∆ that rapidly goes to zero at a critical density where the screening sets in is replaced by a slowly decaying gap for high densities. Looking at c + in Fig. 3 one can see that the BCS state is reached before the gap vanishes, for all band gaps considered. Furthermore, the slowly decaying tail of the gap function and the condensate fractions in Figs. 2-3 suggest that, for the ideal system, it may be possible to stabilize a condensate with a small gap (≈1 meV) and a small condensate fraction (≈ 1%) for high densities also in systems with small band gaps.
Our static mean-field calculations agree with Ref. [10], and similar to these calculations ∆ + >> ∆ − for all bandgaps and densities, which suggest that the two gap equations are decoupled and that a Josephson-like trans- fer of pairs is negligible at all densities. In the full calculations the reduced screening yields much larger gaps. For small bandgaps ∆ ≈ E g , which yields a strong coupling of the gap equations and ∆ + ≈ ∆ − . Hence, contrary to the MF calculations we predict that Josephson-like transfer of pairs is substantial for small bandgaps.
In Ref. 1 a double bilayer graphene system with WSe 2 dielectric barrier of 1.4nm thickness was investigated experimentally. The carrier populations and effective bandgap were controlled by external gates and the tunneling current (I int ) was measured. Around n e = n h = 7·10 11 cm −2 the differential tunneling conductance g int = dI int /dV int was strongly enhanced, which cannot be explained with a single-band tunneling model and hence signals the formation of an interbilayer exciton condensate. The enhanced tunneling was measured at gate voltages corresponding to a bandgap of 6 meV in the single-band model, although the actual bandgap can vary slightly from this value due to many-body effects [33]. The system we consider in Fig. 2 has a thinner barrier and an hBN dielectric, which has a smaller dielectric constant than the WSe 2 barrier in the experimental setup. Hence the calculations in this paper can be considered to provide an upper limit for the superfluid gap in the experimental setup. Within the static mean-field description the gap vanishes already around the density 0.3·10 11 cm −2 for a bandgap of 6 meV, and therefore these calculations qualitatively fails to explain the experimental observations. The full calculation on the other hand has a nonvanishing gap at n = 7 · 10 11 cm −2 , although the maximum gap is at lower densities.
To summarize we have derived a new gap equation for bilayer excitonic systems which includes an effective treatment of the frequency dependence of the electronhole interaction. The frequency dependence has a dramatic effect on the superfluid gap, which is enhanced substantially and survives at much higher densities. Contrary to the usual static mean-field treatment our results qualitatively agree with both to QMC simulations for simplified systems [13] and experimental measurements for real systems [1]. Our method therefore provides an important complement to existing methods to describe bilayer excitonic systems. For double bilayer graphene at zero magnetic field our calculations predict non-vanishing superfluid gaps across the entire BEC-BCS crossover. Double bilayer graphene may therefore provide a convenient solid-state systems that can be used to investigate the new emerging physics in the crossover regime as well as to design new devices. Contrary to static mean-field calculations we predict that the superfluid gaps are substantial even for small band gaps (< 10 meV) and that Josephson-like transfer of pairs can be substantial.
F.N. and F.A. acknowledge financial support from the Knut and Alice Wallenberg Foundation and the Swedish Research Council (Vetenskapsrådet, VR). The computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at LUNARC. F. N. would like to thank N. Prasad and W. Burg for helpful discussions regarding the experimental results in Ref. [1] Supplemental information for "Effects of dynamical screening on the BCS-BEC crossover in double bilayer graphene: Density functional theory for exciton bilayers" To understand why the frequency dependence of the interaction has such a dramatic effect on the gap we consider W (ω), obtained by Padé analytic continuation, in Fig 1 at k = k F and n = 5 · 10 11 cm −2 for different bandgaps. A static approximation to W is sufficient if ReW (ω) has a weak frequency dependence for low frequencies. From Fig.  1 it is clear that this is not the case for these systems. On the contrary, the pole in ImW occurs at low frequencies around 500meV which yields a corresponding Kramer's Kronig feature in the real part. Therefore these systems are not well described by a purely static approximation to W .

A. Density Functional Theory
DFT for excitonic bilayers can be derived in close anology to DFT for superconductors (scDFT) 1-3 The electron-hole Hamiltonian of the bilayer system can in second-quantized form be expressed aŝ whereψ (r) (φ (r)) are the electron (hole) field operators and all terms that can be incorporated as a shift of the chemical potential have been ignored. The electrons and holes are assumed to be confined on their respective layers and electron-hole recombination is ignored. The generalized coordinate r includes both the spatial r and spin σ andT e andT h are the electron and hole kinetic energy operators, respectively. In analogy with Ref. 1 this formulation is a multicomponent formulation of DFT that relies on three densities, 1. The electron density n e (r) = σ ψ † (r, σ)ψ(r, σ) .
According to the Hohenberg-Kohn theorem there is, in thermal equilibrium, a one-to-one mapping between the densities n e (r),n h (r), χ(r, r ) and their conjugate potentials V e (r), V h (r), and ∆(r, r ). The grand canonical potential Ω is minimized by the equilibrium densities where the universal function F is defined as S is the entropy 4 andρ 0 the grand canonical density matrix Now we want to map this to a non-interacting Kohn-Sham system with the same ground-state density as the interacting system. The Kohn-Sham system is described by the following thermodynamic potential: The Hartree energy contains three contributions

B. Bogoliubov equations
The Kohn-Sham Hamiltonian is given bŷ where H e/h = ∓∇ 2 e/h /2 + v e/h KS . Following Refs. 1 and 5 we then definê where the operators γ follow Fermionic anti-commutation relations: Furthermore the transformation is chosen to diagonalize H KŜ where E g is the ground state energy. This implies that the commutators between the Hamiltonian and the operators are given by: From the Fermionic anti-commutation relations of the electron and hole field operators the commutator between the hole and electron annihilation operators and the Hamiltonian are given by: [ψ(r),Ĥ KS ] = H e (r)ψ(r) + dr ∆(r, r )φ † (r ) (26) [φ(r),Ĥ KS ] = H h (r)φ(r) − dr ∆(r, r )ψ † (r ) (27) Now inserting the expressions for the operators in terms of theγ:s, making use of the commutation relation for theγ:s on the left hand side and equating equal coefficients of theγ:s yields the Kohn-Sham Bogoliubov-de Gennes (KS-BdG) equations:

C. Decoupling approximation
Since the electron and hole states are solutions to separate Hamiltonians we need to specify which states that should couple in the decoupling approximation. In practice we will consider band-diagonal singlet phase-coherent coupling, i.e. electrons from band γ with momentum k and spin σ couple to holes in band γ with momentum −k and spin −σ. This is similar to Ref 6 which is the appropriate choice for s-wave band diagonal pairing. This yields The solutions of Eqs. 32-33 are where we have defined E kγ = ξ 2 The corresponding amplitudes u kγ and v kγ are subject to the equations

D. Green's functions
The following Kohn-Sham propagators enter the Feynman diagrams: 1. The usual electron Green's function: 2. The hole Green's function:

The anomalous Green's functions
By inserting the Bogoliubov transformation in the expression for the propagators and also utilizing Tγγ = Tγ †γ † = 0 (47) with the corresponding Fourier transforms Within the decoupling approximation the propagators are given by

E. Nambu-Gorkov Green's functions
Before we proceed we have to define the Nambu-Gorkov propagators. First we define the Nambu-Gorkov field operators asΨ The Green's function is then the tensor product The Kohn-Sham Nambu Gorkov Green's function obeys the equation of motion We will divide the self-energy into the static part which comes from the exchange-correlation potentials and the restΣ (rτ, r τ ) =Σ xc (rτ, r τ ) +Ū (r, r), . (69) The first order self-energy diagram in W (corresponding to the GW -approximation) is given bȳ

G. Sham-Schlüter connection
The The Sham-Schlüter equations 7-9 and the corresponding equations for superconducting systems developed by Marques et al. 3 make use of the equality between the Kohn-Sham and interacting densities: Using the division of the self-energy Eq. 68 we get (henceforth dropping the index eh in F with the implicit notation of F = F eh and F † = F he † ) 1 β ωn σσ1σ2 dr 1 dr 2 [Ḡ s σσ1 (rr 1 ω n )Σ xc σ1σ2 (r 1 , r 2 ω n )Ḡ σ2σ (r 2 rω n )] 11 = 2 β ωn dr 1 v e xc (r 1 )G e s (rr 1 ; ω n )G e (r 1 r; ω n ) + v h xc (r 1 )F s (rr 1 ; ω n )F † (r 1 r; ω n ) + 2 β ωn dr 1 dr 2 ∆ * xc (r 1 r 2 )F s (rr 1 ; ω n )G e (r 2 r; ω n ) + ∆ xc (r 1 , r 2 )G e s (rr 1 ; ω n )F † (r 2 r; ω n ) (78) where we also have used the fact that ωn f (−ω n ) = ωn f (ω n ). Now we write out the first line in each equation explicitly assuming that the self-energy does not couple different spin-channels in Nambu-Gorkov space 1 β ωn dr 1 dr 2 Σ xc 11 (r 1 , r 2 , ω n )G e s (rr 1 ; ω n )G e (r 2 r; ω n ) + Σ xc 22 (r 1 , r 2 )F s (rr 1 ; ω n )F † (r 2 r; ω n ) + 1 β ωn dr 1 dr 2 Σ xc 21 (r 1 r 2 ω n )F s (rr 1 ; ω n )G e (r 2 r; ω n ) + Σ xc 12 (r 1 , r 2 , ω n )G e s (rr 1 ; ω n )F † (r 2 r; ω n ) = 1 β ωn dr 1 v e xc (r 1 )G e s (rr 1 ; ω n )G e (r 1 r; ω n ) + v h xc (r 1 )F s (rr 1 ; ω n )F † (r 1 r; ω n ) + 1 β ωn dr 1 dr 2 ∆ * xc (r 1 r 2 )F s (rr 1 ; ω n )G e (r 2 r; ω n ) + ∆ xc (r 1 , r 2 )G e s (rr 1 ; ω n )F † (r 2 r; ω n ) (82) 1 β ωn dr 1 dr 2 Σ xc 22 (r 1 , r 2 , ω n )G h s (rr 1 ; ω n )F † (r 2 r ; ω n ) + Σ xc 11 (r 1 , r 2 , ω n )F † s (rr 1 ; ω n )G e (r 2 r ; ω n ) + 1 β ωn dr 1 dr 2 Σ xc 12 (r 1 , r 2 , ω n )F † s (rr 1 ; ω n )F † (r 2 r ; ω n ) + Σ xc 21 (r 1 , r 2 , ω n )G h s (rr 1 ; ω n )G e (r 2 r ; ω n ) = 1 β ωn dr 1 v h xc (r 1 )G h s (rr 1 ; ω n )F † (r 1 r ; ω n ) + v e xc (r 1 )F † s (rr 1 ; ω n )G e (r 1 r ; ω n ) + 1 β ωn dr 1 dr 2 ∆ xc (r 1 r 2 )F † s (rr 1 ; ω n )F † (r 2 r ; ω n ) + ∆ * xc (r 1 , r 2 )G h s (rr 1 ; ω n )G e (r 2 r ; ω n ) To simplify the equations we make use of the decoupling approximation and take the matrix-elements of the equations. We will assume time-reversal symmetry, which implies that the LDA eigenstates are real. Furthermore, to simplify the notation we will ignore the band-index in the derivation below. However, it should be noted that the band-index enters in the same way as the k-index, so that the final gap equation is valid also for the multiorbital models. We can then write the Green's function as Furthermore we define and so on for the different electron-electron, electron-hole, hole-hole and hole-electron components. Correspondingly Any matrix element with v e xc including hole wavefunctions is set to zero and vice versa. Now we insert these expressions into the equations above. In the diagonal components, which are derived from the equations of the normal density we integrate over r and use the orthonormality of the wave-functions to get the equations for the single components. Furthermore in the off-diagonal components, which are derived from the anomalous density we multiply by ψ k (r) and φ * −k (r ) and integrate over r and r in the third equation and the other way around in the last equation to get the equations for the single components.
Eliminating v e xc and v h xc yields a linear equation for ∆ This reduces to the mean-field result if we set Σ 22 = Σ 11 = 0 and assume that Σ 21 is a static exchange term, where F γγ kk is the geometrical formfactor (defined below).

H. Further approximations
For practical calculations we will assume that the bare dispersion already includes the effect of the normal selfenergies and therefore set Σ 11 = Σ 22 = 0 in Eq. 108. Furthermore we will approximate the interacting propagators in the gap equation by their Kohn-Sham counterparts and evaluate the anomalous self-energy within the GW -approximation.
We compute the eigenvectors corresponding to the eigenvalues kγ numerically.

J. The screened interaction and self-energy
In the superfluid state the effective electron-hole interaction is given by 6,13,14 where the k and Matsubara indexes are implicit for all quantities and is the effective dielectric constant, d the separation between the bilayers and e the electronic charge. We will now specialize to the case with equal electron hole dispersions so that 1m = 2m and the electron and hole polarizations are equal. For this case the normal and anomoulous polarizations are given by Here the factor 4 comes from the spin and valley degeneracy. In order to evaluate the frequency sums in the gap equation analytically we need an analytic expression for the self-energy. Therefore we have to fit W to a rational function. We choose the method of Ref. 15 and fit W to a plasmon pole functionW where the ω ik are the position of the poles which are obtained by a least-square fitting procedure as described below. In practice we compute W along the imaginary axis using eq. 112 and continue it to the real axis using Padé approximant. In order to be able to probe low temperatures we use an effective frequency mesh for the fitting procedure of W (iν m ) which is sparser than the actual frequency mesh corresponding to the temperature of the calculations. However, all quantities (including the Fermi-Dirac distributions in the polarizations and the final frequency sums) are computed for the correct temperature. Therefore this approximation corresponds to an effective temperature smearing in the analytic continuation of W and for the least-square fit. The resulting parameters a i;nkn k and ω i;nkn k are not sensitive to the choice of mesh.
For the bilayer bandstructure considered in this work ImW k (ω) only contains a single peak (Fig. 2). However, since this peak is asymmetric with respect to the peak maximum it is not well represented by a single plasmon pole approximation. In order to remedy this we approximate W by a multiple plasmon pole approximation with up to 100 poles. We use two trains with equidistant poles, one starting at some lowest frequency ω min and ending just before the plasma frequency and the other starting at or after the plasma frequency and ending at the frequency ω max . ω min and ω max are determined as the frequencies ω where max(ImW (ω))/ImW (ω ) = 1 30 The relative weight between the peaks within each train is determined by the strength of ImW , while the overall weight of each train is determined by the least square fit. We also have an additional peak at the plasma frequency which is given a separate weight in the least square fit. Thus, we determine three parameters in the least square fit ofW (iν m ). The procedure is sketched in Fig. 2. Then the self-energy is given by: Σ 12 γ (q, iω n ) = k,γ ,νl F γ,γ q,q−k F eh γ (q − k, iω n − iν l )W eh (k, ν l ) = k,γ ,νl F γ,γ q,q−k (α kγ l − β kγ l ) (119) k,γ ,νl where n F is the Fermi-dirac distribution and n B the Bose-Einstein distribution. The self-energies contain four different kinds of terms: These expressions can be inserted in the gap equation and the Matsubara sums evaluated analytically using the MatsubaraSum Mathematica package 16 .

K. k-point sampling
We use circular coordinates for the k-point sampling. In the radial k-direction we use an exponential mesh, which is densest around k = 0 and |k| = k F where k F = √ πn. For |k| < k F /2 it is given by The mesh is then mirrored to k F and mirrored again to 3k F /2. After 3k F /2 the mesh is given by x i = (i − N a k ) * α 1 /(N k − N a k ) C = k cut − 1.5 · k F /(exp(α 2 ) − 1).
where N a k are the number of k-points below 3kF 2 and N k the total number of k-points. We use 150-200 k-points with N a k ≈ 3 4 N k , α 1 ≈ 3.5 and α 2 ≈ 5.5. For the angular part θ we use an equidistant mesh with 50-100 points.