Electron-phonon interaction and pairing mechanism in superconducting Ca-intercalated bilayer graphene

Using the ab initio anisotropic Eliashberg theory including Coulomb interactions, we investigate the electron-phonon interaction and the pairing mechanism in the recently-reported superconducting Ca-intercalated bilayer graphene. We find that C6CaC6 can support phonon-mediated superconductivity with a critical temperature Tc = 6.8–8.1 K, in good agreement with experimental data. Our calculations indicate that the low-energy Caxy vibrations are critical to the pairing, and that it should be possible to resolve two distinct superconducting gaps on the electron and hole Fermi surface pockets.

, however predictive ab initio calculations of the critical temperature have not yet been reported. In this work we elucidate the role of the electron-phonon interaction in the normal and in the superconducting state of C 6 CaC 6 by performing state-of-the-art ab initio calculations powered by electron-phonon Wannier interpolation 25,26 . For the normal state we study the electron self-energy and spectral function in the Migdal approximation; for the superconducting state we solve the anisotropic Migdal-Eliashberg equations including Coulomb interactions from first principles. Our main findings are: (i) superconductivity in C 6 CaC 6 can be explained by a phonon-mediated pairing mechanism; (ii) in contrast to bulk CaC 6 , low-energy Ca vibrations are responsible for the majority of the EPC in the superconducting state; (iii) unlike bulk CaC 6 , C 6 CaC 6 should exhibit two superconducting gaps. For clarity all technical details are described in the Methods. Figure 1(a) shows a ball-and-stick model of Ca-intercalated bilayer graphene, while Fig. 1(b-d) show the corresponding band structure, Brillouin zone, and Fermi surface, respectively. Two sets of bands cross the Fermi level around the Γ point. The bands labeled as α * and β * in Fig. 1(b) represent π * states, and are obtained by folding the Dirac cone of graphene from K to Γ , following the superstructure induced by Ca. The band labeled as IL is the Ca-derived nearly-free electron band, which disperses upwards from about 0.5 eV below the Fermi energy at Γ 1,23,24 . In Fig. 1(d) we see two sets of Fermi surface sheets: triangular hole pockets around the K′ points, corresponding to the α ⁎ 1 states; and a bundle of small electron pockets around Γ , which arise from the other π * states (α ⁎ 2 , β ⁎ 1 , and β ⁎ 2 ) and from the IL band. In order to quantify the strength of EPCs for each of these bands we calculate the spectral function in the normal state using the Migdal approximation 27 . Figure 1(e-g) show that the EPC induces sudden changes of slope in the bands, which are referred to as 'kinks' in high-resolution ARPES experiments. A pronounced kink is seen in all bands at a binding energy of 180 meV. This feature can be assigned to the coupling with the in-plane C xy stretching modes 28 (see Supplementary Fig. S1 for the phonon dispersion relations). The occurrence of this high-energy kink in the π * bands is consistent with the observed broadening of the quasiparticle peaks in the ARPES spectra of bulk CaC 6 in the same energy range 13 . This feature results from a sharp peak at 180 meV in the real part of the electron self-energy of the π * electron pockets, as shown in Supplementary Fig. S2. A second kink is clearly seen at a binding energy of 70 meV, and is mostly visible for the π * bands defining the hole pockets [magenta line in Fig. 1(e)] and for the interlayer band [cyan line in Fig. 1(g)]. This low-energy kink corresponds to a second, smaller peak at the same energy in the real part of the electron self-energy ( Supplementary Fig. S2), and arises from the coupling with the out-of-plane C z modes of the graphene sheets. Closer inspection reveals also a third kink around 12 meV, however this feature is hardly discernible and is unlikely to be observed in ARPES experiments. This faint structure arises from the coupling to the in-plane Ca xy vibrations, and can be seen more clearly in the real and the imaginary parts of the electron self-energy in Supplementary Fig. S2.
From the calculated electron self-energy we can extract the electron-phonon mass enhancement parameter λ F for each band using the ratio between the bare and the renormalized Fermi velocities. Along the Γ K′ direction we obtain λ F = 0.53 (α ⁎ 1 ), 0.48 (α ⁎ 2 ), 0.30 (β ⁎ 1 and β ⁎ 2 ) for the π * bands, and λ F = 0.68 for the IL band. Along the other two directions Γ M′ and M′ K′ the mass enhancement parameters are up to a factor of three smaller than the corresponding values along Γ K′ , suggesting a rather anisotropic EPC. Our findings are in agreement with ARPES studies on bulk KC 6 and CaC 6 superconductors 10,22 and also with recent work on Rb-intercalated bilayer graphene 3 .
We now move from the normal state to the superconducting state, and consider an electron-phonon pairing mechanism in analogy with bulk CaC 6 . Figure 2(a) shows a comparison between the isotropic Eliashberg spectral functions, α 2 F(ω), and the cumulative total EPC, λ(ω), calculated for C 6 CaC 6 and for bulk CaC 6 . In both cases we can distinguish three main contributions to the EPC, to be associated with the Ca xy vibrations (~10 meV), the out-of-plane C z vibrations (~70 meV), and the in-plane C xy modes (~180 meV). The Eliashberg functions of C 6 CaC 6 and CaC 6 look similar in shape, and the total EPC is λ = 0.71 in both cases. However the relative contributions of each set of modes differ considerably. The low-energy Ca modes are slightly softer in C 6 CaC 6 , leading to a larger contribution to the EPC than in CaC 6 . These modes account for 60% of the total coupling in the bilayer, while they only account for less than 30% of the coupling in bulk CaC 6 . In both cases the EPC strength associated with the in-plane C-C stretching modes is too weak (15% of the total) to make a sizable contribution to the superconducting pairing. This is somewhat counterintuitive, given that the C xy modes lead to the most pronounced kinks in the spectral function in Fig. 1. The smaller contribution of the out-of plane C z modes to the EPC and the softening of the in-plane Ca xy vibrations, obtained when going from the bulk to the bilayer, are similar to the results found for Li-and Ca-decorated monolayer graphene 29 . In the monolayer case, the removal of quantum confinement causes a shift of the IL wave function farther away from the graphene layer as compared to bulk, and, therefore, the EPC coupling between π * and IL states mediated by C z vibrations is reduced. Although in the bilayer case the IL state is strongly localized around the Ca atom, the fact that the interlayer charge density is only present between the graphene layers gives rise to a weaker coupling of the out-of-plane C z modes with the interlayer electrons, and, therefore, a lower contribution to the global EPC.
To check whether the softening of the low energy Ca xy modes is due to the difference in the structural parameters between the bilayer and the bulk, we recalculate the vibrational spectrum of bilayer C 6 CaC 6 after setting the in-plane lattice constant and the interlayer distance to the values of the bulk. As shown in Supplementary  Fig. S3, the in-plane C xy modes (ω > 105 meV) soften due to the increase in the in-plane lattice constant, and the out-of-plane C z and Ca z modes (20 < ω < 45 meV) harden due to the decrease in the interlayer distance. On the other hand, the two lowest-lying modes involving mainly in-plane Ca xy vibrations are considerably less affected, although a closer inspection reveals an overall slight hardening of approximately 1 meV, particularly along the Γ M′ and M′ K′ directions. This suggests that the observed softening of the low-energy Ca xy modes from bulk to bilayer is most likely caused by changes in the electronic structure which is consistent with an increased density of states at the Fermi level.
In order to determine the superconducting critical temperature of C 6 CaC 6 we solve the anisotropic Eliashberg equations 26,30,31 , with the Coulomb pseudopotential μ * calculated from first principles (the calculation of μ * is discussed below). In Fig. 2(b,c) we plot the energy-dependent distribution of the superconducting gap, separated into contributions corresponding to the two sets of Fermi surfaces centered around the Γ and K′ points. We see that two distinct gaps open on the Γ -centered electron pockets, with average values Δ 1 = 0.55 meV and Δ 2 = 1.29 meV at zero temperature. The Δ 1 gap is characterized by a very narrow energy profile and the EPC on these pockets is essentially isotropic, resulting primarily from the coupling with the Ca xy phonons [ Fig. 2(d)]. The Δ 2 gap exhibits a much broader energy profile (0.82 < Δ 2 < 1.75 meV) and originates mainly from the coupling to Ca modes and out-of-plane C z phonons [ Fig. 2(e)]. In between the two Γ -centered gaps, a third gap with an average value Δ 3 = 0.99 meV opens on the triangular hole pockets around K′ [α ⁎ 1 states] [ Fig. 2(f)]. These states couple primarily to C xy phonons and the gap has an anisotropic character with a large spread in energy (0.73 < Δ 3 < 1.25 meV). In Supplementary Fig. S4 we show the anisotropic EPC parameters λ k leading to this superconducting gap structure. For completeness, we compare our results with the gap structure of bulk CaC 6 . In the latter case, only one superconducting gap is predicted (see Supplementary Fig. S5), in agreement with previous theoretical studies 20,32 and experiments 33,34 . Although multiple sheets of the Fermi surface contribute to the superconducting gap, there is no separation into distinct gaps, giving rise to a smeared multigap structure 20,32 . This situation is similar to the Δ 2 and Δ 3 gaps in the bilayer case. Based on our results we suggest that in C 6 CaC 6 high-resolution ARPES experiments might be able to resolve two distinct gaps on the electron pockets, corresponding to Δ 1 and Δ 2 , but only one gap on the triangular hole pockets, corresponding to Δ 3 .
The calculations of the gap function and the superconducting critical temperature in Fig. 2(b,c) require the knowledge of the Coulomb pseudopotential μ * . In order to determine this parameter we first calculate the dimensionless electron-electron interaction strength within the random-phase approximation, obtaining μ = 0.254 (see Methods). Then we renormalize this interaction using μ * = μ/[1 + μlog(ω pl /ω ph )] 35 , where ω pl and ω ph are characteristic electron and phonon energy scales, respectively. We set ω pl = 2.5 eV, corresponding to the lowest plasmon resonance in GICs 36,37 , and ω ph = 200 meV, corresponding to the highest phonon energy in C 6 CaC 6 . Figure 2(g) shows the variation of the Coulomb pseudopotential μ k across the Fermi surface. The repulsive interaction is strongest on the electron pockets, and weakest on the hole pockets. Since the EPC exhibits a very similar anisotropy, the net coupling strength λ µ − ⁎ k k is only moderately anisotropic and can be replaced by a single, average μ * in the Eliashberg equations. From Fig. 2(g) we obtain μ * = 0.155, and this is the value employed in Fig. 2(b,c). As we can see in Fig. 2(b,c), our Eliashberg calculations yield a superconducting critical temperature T c = 8.1 K, which is only slightly higher than the experimental value = T 4 c exp K reported in ref. 4 . In order to check the role of anisotropic Coulomb interactions we repeat the calculations by considering a Coulomb pseudopotential resolved over the electron and hole pockets. In this case we find the decrease in T c to be very small, Δ T c = 0.3 K. Furthermore we explore the sensitivity of the calculated T c to the choice of the characteristic phonon energy ω ph . To this aim we solve the Eliashberg equations again, this time by setting ω ph equal to the Matsubara frequency cutoff (5 × 200 meV, see Methods). This alternative choice leads to μ * = 0.207 and T c = 6.8 K (Supplementary Fig. S6), which is in even better agreement with experiment. For completeness in Supplementary Fig. S7 we also show the dependence of T c on the characteristic energy ω c as obtained using the standard McMillan formula 38 . Consistent with our Eliashberg calculations, we find that large variations of ω c only change T c by a few K's. These additional tests show that our results are solid, therefore we can safely claim that the ab initio Eliashberg theory yields T c = 6.8-8.1 K for C 6 CaC 6 . The close agreement between these values and experiment supports the notion that Ca-intercalated bilayer graphene is a conventional phonon-mediated superconductor.
In conclusion, we studied entirely from first principles the electron-phonon interaction and the possibility of phonon-mediated pairing in the newly-discovered superconducting C 6 CaC 6 . We showed that the Ca vibrations play an important role in the pairing but do not carry a sharp signature in the normal-state band structure; conversely the high-frequency in-plane C xy vibrations lead to pronounced photoemission kinks but have a small contribution to the pairing. The good agreement between the critical temperature calculated here and the recent experiments of ref. 4 indicate that Ca-intercalated bilayer graphene is an electron-phonon superconductor. The present work calls for high-resolution spectroscopic investigations, as well as for calculations based on alternative ab initio methods 39,40 , in order to test our prediction of two distinct superconducting gaps in C 6 CaC 6 .

Methods
The calculations are performed within the local density approximation to density-functional theory 41,42 using planewaves and norm-conserving pseudopotentials 43,44 , as implemented in the Quantum-ESPRESSO suite 45 . The planewaves kinetic energy cutoff is 60 Ry and the structural optimization is performed using a threshold of 10 meV/Å for the forces. C 6 CaC 6 is described using the × 3 3 R30° supercell of graphene with one Ca atom per unit cell, and periodic images are separated by 15 Å. The optimized in-plane lattice constant and interlayer separation are a = 4.24 Å and d = 4.50 Å, respectively. Bulk CaC 6 is described using the rhombohedral lattice, and the optimized lattice constant and rhombohedral angle are a = 5.04 Å and α = 50.23°, respectively. The electronic charge density is calculated using an unshifted Brillouin zone mesh with 24 2 and 8 3 k-points for C 6 CaC 6 and CaC 6 , respectively, and a Methfessel-Paxton smearing of 0.02 Ry. The dynamical matrices and the linear variation of the self-consistent potential are calculated within density-functional perturbation theory 46 on the irreducible set of regular 6 2 (C 6 CaC 6 ) and 4 3 (CaC 6 ) q-point grids.
The electron self-energy, spectral function, and superconducting gap are evaluated using the EPW code 25,26,47 . The electronic wavefunctions required for the Wannier-Fourier interpolation 48,49 in EPW are calculated on uniform unshifted Brillouin-zone grids of size 12 2 (C 6 CaC 6 ) and 8 3 (CaC 6 ). The normal-state self-energy is calculated on fine meshes consisting of 100,000 inequivalent q-points, using a broadening parameter of 10 meV and a temperature T = 20 K. For the anisotropic Eliashberg equations we use 120 × 120 (40 × 40 × 40) k-point grids and 60 × 60 (20 × 20 × 20) q-point grids for C 6 CaC 6 (CaC 6 ). The Matsubara frequency cutoff is set to five times the largest phonon frequency (5 × 200 meV), and the Dirac delta functions are replaced by Lorentzians of widths 100 meV and 0.5 meV for electrons and phonons, respectively. The technical details of the Eliashberg calculations are described extensively in ref. 26 .
The electron-electron interaction strength is obtained as μ = N F 〈 〈 V k,k′ 〉 〉 FS , where N F is the density of states at the Fermi energy, , and W is the screened Coulomb interaction within the random phase approximation 50 . Here 〈 〈 ⋅ 〉 〉 FS indicates a double Fermi-surface average, and k stands for both momentum and band index. The screened Coulomb interaction is calculated using the Sternheimer approach 51,52 . Linear-response equations are solved using 26 × 26 Brillouin-zone grids, corresponding to 70 inequivalent points. The energy cutoff of the dielectric matrix is 10 Ry (815 planewaves).