Superconducting Symmetry Phases and Dominant bands in (Ca-) Intercalated AA- Bilayer Graphene

TL;DR

This paper models superconductivity in calcium-intercalated bilayer graphene, revealing dominant d-wave pairing and the possibility of single or dual-gap phases depending on band interactions near the Fermi energy.

Contribution

It introduces a multiband tight-binding model with mirror symmetry to analyze superconducting phases, identifying dominant d-wave pairing and conditions for single or dual-gap superconductivity.

Findings

D-wave superconductivity is dominant in both sectors.

Two flat bands lead to two distinct superconducting gaps.

Transition from single-gap to dual-gap superconductivity is possible between 0-6 K.

Abstract

Built on a realistic multiband tight-binding model, mirror symmetry is used to map a calcium-intercalated bilayer graphene Hamiltonian into two independent single layer graphene-like Hamiltonians with renormalized hopping. The quasiparticles exhibit two types of chirality. Here a quasi-particle consists of two electrons from opposing layers where possess an additional quantum number called "cone index" which can be regarded as the eigenvalue of mirror symmetry operations. To obtain tight-binding parameters, both effective monolayer Schrodinger equations are solved analytically and fitted to first-principles band structure results. Two quasi-particles (four electrons) can team up to build a Cooper pair with even or odd chirality. Treatment of the pairing Hamiltonian leads to two decoupled gap equations. The pairing of quasi-particles with different cone indexes is forbidden. The decoupled gap equations are solved analytically to obtain all the possible superconducting phases. Two nearly "flat bands" crossing the Fermi energy, each related to the graphene-like structures, are responsible for two distinct superconductivity gaps that emerge. Depending on how much these bands are affected by the intercalant and which is closer to the Fermi energy, distorted s-wave or d-wave superconductivity may become dominant. Numerical calculations reveal that d-wave superconductivity is dominant in both sectors. For these two dominant phases, within the range of 0-6 K which superconductivity has been observed, numerically the transition from single-gap to dual-gap superconductivity is possible. Adopting the two-gap viewpoint of superconductivity in C$_6$CaC$_6$, the dominant $d$-wave states should have the same critical temperature. Around $T_c=2K$ these two relations intersect, otherwise, superconductivity has been realized just in one of these two sectors and disappears in the other one.

Tables

Table 1. Table 1: The C-C hopping parameters (in eV) for the six odd symmetry bands of C 6 CaC 6 are denoted by t 0 ​ n ( − ) ​ C ​ C subscript superscript 𝑡 𝐶 𝐶 0 𝑛 t^{(-)CC}_{0n} where the index n 𝑛 n indicates n 𝑛 n -th neighbour.

$n$	$0$	$1$	$2$	$3$	$4$	$5$	$6$	$7$	$8$	$9$
$t_{0n}^{(-)CC}$	${ϵ}_c^{-}=-1.00$	$t_1^{-}=3.04$	$t_1^{{}_{}{}^{\prime }-}=0.92t_1^{-}$	$t_2^{-}=-0.23$	$t_2^{{}_{}{}^{\prime }-}=0.92t_2$	$t_3^{-}=-0.29$	$t_3^{{}_{}{}^{\prime }-}\approx t_3$	$t_4^{-}=-0.02$	$t_4^{{}_{}{}^{\prime }-}\approx t_4^{-}$	$t_5^{-}=-0.05$

Table 2. Table 2: The C-C hopping parameters (in eV) for the seven even-symmetry bands of C 6 CaC 6 are denoted by t 0 ​ n ( + ) ​ C ​ C subscript superscript 𝑡 𝐶 𝐶 0 𝑛 t^{(+)CC}_{0n} where the index n 𝑛 n indicates the n 𝑛 n -th neighbour. In the intercalant plane, Ca-Ca hopping parameters are denoted by t 0 ​ m C ​ a ​ C ​ a superscript subscript 𝑡 0 𝑚 𝐶 𝑎 𝐶 𝑎 t_{0m}^{CaCa} where m 𝑚 m is m 𝑚 m -th Ca neighbor of central Ca. The modified Ca-C hopping parameter is t 0 ​ m C ​ a ​ C superscript subscript 𝑡 0 𝑚 𝐶 𝑎 𝐶 t_{0m}^{CaC} .

$n$	$0$	$1$	$2$	$3$	$4$	$5$	$6$	$7$	$8$	$9$
$t_{0n}^{(+)CC}$	${ϵ}_c^{+}=-0.60$	$t_1^{+}=2.94$	$t_1^{{}_{}{}^{\prime }+}=0.92t_1$	$t_2^{+}=-0.24$	$t_2^{{}_{}{}^{\prime }+}=0.92t_2$	$t_3^{+}=0.27$	$t_3^{{}_{}{}^{\prime }+}\approx t_3^{+}$	$t_4^{+}=-0.02$	$t_4^{+}\approx t_4^{+}$	$t_5^{+}=-0.08$
$m$	$0$	$1$	$2$	$3$	$4$		5
$t_{0m}^{CaCa}$	${ϵ}_{Ca}=1.12$	$-0.35$	$0.06$	$0.06$	$-0.02$	$0.00$
$t_{0m}^{CaC}$	$-$	$0.17$	$-0.14$	$0.08$	$-0.07$	$-0.05$

Equations

$$\hat{H}=-\sum _{i\alpha }\sum _{j\beta ,\sigma }{t}_{i\alpha ,j\beta }^{\sigma ,\sigma }{\hat{c}}_{i\alpha \sigma }^{†}{\hat{c}}_{j\beta \sigma }+\sum _{i\alpha ,\sigma }({ϵ}_{i\alpha }-{\mu }_o){\hat{n}}_{i\alpha \sigma }+\frac{1}{2}\sum _{i\alpha ,\sigma }\sum _{j\beta ,\sigma {}^{\prime }}{U}_{i\alpha ,j\beta }^{\sigma ,\sigma {}^{\prime }}{\hat{n}}_{i\alpha \sigma }{\hat{n}}_{j\beta \sigma {}^{\prime }}={\hat{H}}_N+{\hat{H}}_P.$$

$$\displaystyle\left[{\begin{array}[]{*{20}{c}}H_{11}&\vline&H_{12}&\\ \hline\cr H_{21}&\vline&H_{22}&\end{array}}\right]\left({\begin{array}[]{*{20}{c}}\chi_{1}\\ \chi_{2}\\ \end{array}}\right)=E_{n}\left({\begin{array}[]{*{20}{c}}\chi_{1}\\ \chi_{2}\\ \end{array}}\right),~{}~{}~{}\left({\begin{array}[]{*{20}{c}}\chi_{1}\\ \chi_{2}\\ \end{array}}\right)=\left({\begin{array}[]{*{20}{c}}{\mathscr{C}}_{A_{1}}(\vec{k})\\ {\mathscr{C}}_{B_{1}}(\vec{k})\\ \hline\cr{\mathscr{C}}_{A_{2}}(\vec{k})\\ {\mathscr{C}}_{B_{2}}(\vec{k})\\ \end{array}}\right)$$

$${\mathcal{C}}_{A_1}(\vec{k})=\pm {\mathcal{C}}_{A_2}(\vec{k}),{\mathcal{C}}_{B_1}(\vec{k})=\pm {\mathcal{C}}_{B_2}(\vec{k})$$

$$H^{+}{\chi }_R=E_n^{+}{\chi }_R,H^{-}{\chi }_L=E_n^{-}{\chi }_l,$$

$$\displaystyle\chi_{R}=\frac{1}{\sqrt{2}}(\chi_{1}+\chi_{2})=\left({\begin{array}[]{*{20}{c}}{\mathscr{C}}_{A}^{+}(\vec{k})\\ {\mathscr{C}}_{B}^{+}(\vec{k})\\ \end{array}}\right),~{}~{}~{}\chi_{L}=\frac{1}{\sqrt{2}}(\chi_{1}-\chi_{2})=\left({\begin{array}[]{*{20}{c}}{\mathscr{C}}_{A}^{-}(\vec{k})\\ {\mathscr{C}}_{B}^{-}(\vec{k})\\ \end{array}}\right).$$

$$H=-\sum _{i,j}\sum _{\sigma =\pm }{t}_{i,j}^{\sigma }{a}_{i\sigma }^{†}{b}_{j\sigma }+\sum _{i\sigma }{\epsilon }^{\sigma }(a_{i\sigma }^{†}{a}_{i\sigma }+b_{i\sigma }^{†}{b}_{i\sigma })+h.c$$

$$\displaystyle\left[{\begin{array}[]{*{20}{c}}H_{11}&\vline&H_{12}&\\ \hline\cr H_{12}&\vline&H_{11}&\end{array}}\right]\approx-\left[{\begin{array}[]{*{20}{c}}{\mu}&{t_{1}}f(\vec{k})&\vline&{{\gamma_{0}}}&{\gamma_{1}}f(\vec{k})&\\ {{t_{1}}f^{*}(\vec{k})}&\mu&\vline&{{\gamma_{1}}f^{*}(\vec{k})}&{\gamma_{0}}&\\ \hline\cr{{\gamma_{0}}}&{\gamma_{1}}f(\vec{k})&\vline&\mu&{t_{1}}f(\vec{k})&\\ {{\gamma_{1}}f^{*}(\vec{k})}&{\gamma_{0}}&\vline&{t_{1}}f^{*}(\vec{k})&\mu&\\ \end{array}}\right];~{}~{}H^{\pm}=H_{11}\pm H_{12}=-\left[{\begin{array}[]{*{20}{c}}&{\mu\pm\gamma_{0}}&{t^{\pm}}f(\vec{k})&\\ &{{t^{\pm}}f^{*}(\vec{k})}&{\mu\pm\gamma_{0}}&\\ \end{array}}\right]$$

$$E_n^{\pm }(\vec{k})=\mu +{\eta }_{v-c}(t^{\pm }|f(\vec{k})|\pm {\gamma }_0)$$

$$v_f^{+}-v_f^{-}=\frac{\sqrt{3}a}{\hslash }{\gamma }_1.$$

$${\hat{H}}_N=-\sum _{i\alpha }\sum _{j\beta ,\sigma }{t}_{i\alpha ,j\beta }^{\sigma ,\sigma }{c}_{i\alpha \sigma }^{†}{c}_{j\beta \sigma }+\sum _{i\alpha ,\sigma }({ϵ}_{i\alpha }-{\mu }_o){\hat{n}}_{i\alpha \sigma }.$$

$$\sum _{\beta =0}^{12}{\epsilon }_{\alpha \beta }(\vec{k}){\mathcal{C}}_{\beta }+({ϵ}_{\alpha }-{\mu }_o){\mathcal{C}}_{\alpha }=E(\vec{k}){\mathcal{C}}_{\alpha }\text{mbox}where{\epsilon }_{\alpha \beta }(\vec{k})=-\frac{1}{N}\sum _{ij}{e}^{i\vec{k}.({\vec{r}}_{i\alpha }-{\vec{r}}_{j\beta })}{t}_{i\alpha j\beta }^{\sigma \sigma }.$$

$${\mathcal{C}}_{\alpha }(\vec{k})=\pm {\mathcal{C}}_{\alpha +6}(\vec{k}),\alpha =1,2,...,6$$

$$H^{-}(\vec{k}){\mathcal{C}}_n(\vec{k})=\left(H_{11}-H_{12}\right){\mathcal{C}}_n(\vec{k})=E_n^{-}{\mathcal{C}}_n(\vec{k}),n=1,2,...,6.$$

$$\displaystyle H^{+}_{c}(\vec{k})\left({\begin{array}[]{*{20}{c}}{\mathscr{C}}_{0}(\vec{k})\\ {\sqrt{2}{\bm{\mathcal{C}}_{n}}{\vec{k}})}\end{array}}\right)=\left[{\begin{array}[]{*{20}{c}}{{h_{0}}(\vec{k})}&\vline&{\sqrt{2}{h_{01}}(\vec{k})}\\ \hline\cr{\sqrt{2}{h_{10}}(\vec{k})}&\vline&H^{+}(\vec{k})\end{array}}\right]\left({\begin{array}[]{*{20}{c}}{\mathscr{C}}_{0}(\vec{k})\\ \sqrt{2}{{\bm{\mathcal{C}}_{n}}({\vec{k}})}\end{array}}\right)=E_{n}^{+}\left({\begin{array}[]{*{20}{c}}{\mathscr{C}}_{0}(\vec{k})\\ {\sqrt{2}{\bm{\mathcal{C}}_{n}}({\vec{k}})}\end{array}}\right),$$

$$E_g^{+}=2|t_1^{+}-t_1^{\prime +}|,E_g^{-}=2|t_1^{-}-t_1^{\prime -}|$$

$$E_g^{+}-E_g^{-}=4\Delta {\gamma }_1,E_g^{+}+E_g^{-}=4\Delta t$$

$$E_{sh;m,l}^{\pm }(t_i,\overrightarrow{{\xi }_i},\vec{k})={\mu }_m^{\pm }(\vec{k})-{\mu }_o+\frac{1}{2}\left[{\epsilon }_A+{\epsilon }_B+{(-1)}^l\sqrt{{\left({\epsilon }_A-{\epsilon }_B\right)}^2+4w_m^{\pm }(t_i,\overrightarrow{{\xi }_i},\vec{k})}\right],m=1,2,3;l=1,2$$

$${\hat{H}}_N={\hat{H}}_N^{+}\oplus {\hat{H}}_N^{-}=\sum _{ij\sigma }\sum _{\alpha ,\beta =0}^6{t}_{i\alpha \sigma ,j\beta \sigma }^{+}{\hat{c}}_{i\alpha \sigma }^{+†}{\hat{c}}_{j\beta \sigma }^{+}\oplus \sum _{ij\sigma }\sum _{\alpha ,\beta =1}^6{t}_{i\alpha \sigma ,j\beta \sigma }^{-}{\hat{c}}_{i\alpha \sigma }^{-†}{\hat{c}}_{j\beta \sigma }^{-}$$

$${\hat{c}}_{i\alpha \sigma }^{(\pm )†}=\frac{1}{\sqrt{2}}({\hat{c}}_{i\alpha \sigma }^{†}\pm {\hat{c}}_{i,\alpha +6,\sigma }^{†}),t_{i\alpha \sigma ,j\beta \sigma }^{\pm }=t_{i\alpha \sigma ,j\beta \sigma }^{inter}\pm t_{i\alpha \sigma ,j\beta +6\sigma }^{intra},\alpha ,\beta =1,...,6.$$

$$\displaystyle{\hat{H}}_{su}=\sum_{\vec{k}}\Lambda^{{\dagger}}(\vec{k})\left({\begin{array}[]{*{20}{c}}{{H_{su}^{+}}(\vec{k})}&0\\ 0&{H_{su}^{-}(\vec{k})}\end{array}}\right)\Lambda(\vec{k})$$

$$E_{m,s}^{Q+}(\vec{k})=s\left(E_m^{+}(\vec{k})+\sum _{n=1}^7\frac{{\left|{\Delta }_{mn}^{+}(\vec{k})\right|}^2}{E_m^{+}(\vec{k})+E_n^{+}(\vec{k})}\right),{\Delta }_{mn}^{+}(\vec{k})=\sum _{\alpha =1}^9{\Omega }_{mn}^{+\alpha }(\vec{k}){\Delta }_{+}^{\alpha }s=\pm 1,m=1,2,...7$$

$$E_{m,s}^{Q-}(\vec{k})=s\left(E_m^{-}(\vec{k})+\sum _{n=8}^{13}\frac{{\left|{\Delta }_{mn}^{-}(\vec{k})\right|}^2}{E_m^{-}(\vec{k})+E_n^{-}(\vec{k})}\right),{\Delta }_{mn}^{-}(\vec{k})=\sum _{\alpha =1}^9{\Omega }_{mn}^{-\alpha }(\vec{k}){\Delta }_{-}^{\alpha }s=\pm 1,m=8,9,...13$$

$${\Omega }_{mn}^{1\pm }(\vec{k})$$

$${\Omega }_{mn}^{2\pm }(\vec{k})$$

$${\Omega }_{mn}^{3\pm }(\vec{k})$$

$${\Omega }_{mn}^{4\pm }(\vec{k})$$

$${\Omega }_{mn}^{5\pm }(\vec{k})$$

$${\Omega }_{mn}^{6\pm }(\vec{k})$$

$${\Omega }_{mn}^{7\pm }(\vec{k})$$

$${\Omega }_{mn}^{8\pm }(\vec{k})$$

$${\Omega }_{mn}^{9\pm }(\vec{k})$$

$${\Delta }_{mn}^{\pm }(\vec{k})=g_0^{\pm }⟨{\hat{d}}_m^{\pm \uparrow }(\vec{k}){\hat{d}}_n^{\pm \downarrow }(\vec{k})⟩$$

$$\displaystyle\left[\begin{array}[]{*{20}{c}}A^{\pm}&B^{\pm}&B^{\pm}\\ B^{\pm}&C^{\pm}&D^{\pm}\\ B^{\pm}&D^{\pm}&C^{\pm}\end{array}\right]\left(\begin{array}[]{*{20}{c}}\Sigma_{i}\pm\Sigma_{i}^{{}^{\prime}}\\ \Pi_{i}\pm\Pi_{i}^{{}^{\prime}}\\ \Delta_{i}\pm\Delta_{i}^{{}^{\prime}}\end{array}\right)=-\frac{1}{g_{0}^{\pm}}\left(\begin{array}[]{*{20}{c}}\Sigma_{i}\pm\Sigma_{i}^{{}^{\prime}}\\ \Pi_{i}\pm\Pi_{i}^{{}^{\prime}}\\ \Delta_{i}\pm\Delta_{i}^{{}^{\prime}}\end{array}\right)$$

$$\displaystyle A^{\pm}=\left[\begin{array}[]{*{20}{c}}\Gamma^{\pm}_{11}&\Gamma^{\pm}_{12}&\Gamma^{\pm}_{12}\\ \Gamma^{\pm}_{12}&\Gamma^{\pm}_{11}&\Gamma^{\pm}_{12}\\ \Gamma^{\pm}_{12}&\Gamma^{\pm}_{12}&\Gamma^{\pm}_{11}\end{array}\right],\;C^{\pm}=\left[\begin{array}[]{*{20}{c}}\Gamma^{\pm}_{44}&\Gamma^{\pm}_{45}&\Gamma^{\pm}_{45}\\ \Gamma^{\pm}_{45}&\Gamma^{\pm}_{44}&\Gamma^{\pm}_{45}\\ \Gamma^{\pm}_{45}&\Gamma^{\pm}_{45}&\Gamma^{\pm}_{44}\end{array}\right],\;\;B^{\pm}=\left[\begin{array}[]{*{20}{c}}\Gamma_{14}&\Gamma^{\pm}_{15}&\Gamma^{\pm}_{15}\\ \Gamma^{\pm}_{15}&\Gamma^{\pm}_{14}&\Gamma^{\pm}_{15}\\ \Gamma^{\pm}_{15}&\Gamma^{\pm}_{15}&\Gamma^{\pm}_{14}\end{array}\right],\;\;D^{\pm}=\left[\begin{array}[]{*{20}{c}}\Gamma^{\pm}_{47}&\Gamma^{\pm}_{48}&\Gamma^{\pm}_{48}\\ \Gamma^{\pm}_{48}&\Gamma^{\pm}_{47}&\Gamma^{\pm}_{48}\\ \Gamma^{\pm}_{48}&\Gamma^{\pm}_{48}&\Gamma^{\pm}_{47}\end{array}\right].$$

Full text

I Introduction

Discovery of new superconducting phases, often at low temperature, has been one of the active achievements in recent decades, sometimes overshadowed by the high profile effort to push toward higher temperature superconductors. Remarkable progress in the synthesis or fabrication of nanostructures in the recent years is opening new horizons in engineering the physical properties of new classes of materials. Two dimensional (2D) phases, including superconducting ones,brun2017 ; Saito2016 received renewed attention from the discovery that monolayer graphene can be peeled off, and subsequently found that it can be synthesized from the gas phase. Dirac (charge neutrality) points, has become of rising interest. While graphene, with its massless Dirac points, hosts many unusual quantum properties with potential applications, pristine graphene does not support superconductivity due to its vanishing density of states (DOS) at the Dirac points.

Despite that, discussion of superconductivity in doped graphene has been profuse with theoretical modelsUchoa2007 ; Nandkishore2012 ; Nandkishore2014 ; Black-Schaffer2007 ; Kiesel2012 ; Ma2014 . With the many pictures raising various possibilities (such as chiral d+id pairingNandkishore2012 ; Nandkishore2014 ; Black-Schaffer2007 , coexistence of both chiral d+id and f-waveKiesel2012 and p+ipMa2014 ,) but little of a certain nature. However, superconductivity in lithium decorated graphene around 6K was predicted theoreticallyProfeta2012 and subsequently reported experimentally.Ludbrook2015 Graphite intercalation compound (GIC) superconductors have been known and studied for decades.Weller-naturephysics2005 Recent advances in few layer graphene (FLG) fabrication methods have caused theoretical attention to their electronic physics. Due to its strong two dimensionality, bilayer graphene (BlG) has provided an attractive platform for studying 2D electron correlation effects.Nori2016 Because of the weak interlayer van der Waals forces, the layers of the graphene bilayer can rotate relative to each other and form different ordering stakes i.e. AA, AB (Bernal phase), or even “twisted bilayer” form.

Following the indication of a Mott insulating phase in twisted bi-layered grapheneCao-Mott2018 and the observation of superconductivity upon doping,Cao2018 these types of systems have revived interest and may realize a new class of superconductor. It seems that bilayer graphene can host superconductivity, magnetism, and other unusual phases. In addition to relativistic character of quasi particles in single layer graphene, interlayer coupling causes the bilayer graphene to exhibit behaviors that are not observed in the single layer graphene. Most fascinating behaviors of TBlG are played by interlayer hybridization of nearby Dirac K points of opposite layers via interlayer tunneling in a spatially periodic wayMacDonald2019 . Interlayer coupling causes quasi particles in the AB-stacked bilayer graphene behave as massive chiral quasi-particle with parabolic dispersion near the Dirac Fermi points. In this manuscript we will see that quasi-particles in AA-stacked BlG can exhibit extra aspects of such behaviors.

Beside TBlG, when the electron-electron correlations effects are taken into consideration in AA and AB-stacking structures, theory suggests a variety of instabilities with potential technological applications, including unusual quantum Hall effects, antiferromagnetic phases, and tunable band gap opening at the charge neutral points (for a review see [Nori2016, ]), some of which have been confirmed experimentally.Velasco2012 . Although it is difficult to distinguish AA-stacked BlG from single layer graphene but Some authors claimed to fabricate AA-BlG experimentally. However, it has been received much less attention than more stable AB-stacked phase. Unlike Fermi points in monolayer graphene and AB-stacked, the well nested Fermi surface of pristine AA-BlG consists of small electron and hole pockets of equal area. This feature has drastic consequences, tending the system toward electron instabilities such as antiferromagnetism at zero doping and bilayer exciton condensation when doped.Nori2016 ; Rakhmanov2012 ; Akzyanov2013 ; Sboychakov2013 ; SboychakovMetal2013

Superconducting instabilities in doped or gated AB-stacking phases have been predicted. An effective two-band Hamiltonian with attractive interactions was used in RefVucicevic2012 to investigate the possibility of a time reversal breaking $d+id$ phase in moderately doped AB-stacked BlG. Using a weak-coupling renormalization group formalism, the possibility of unconventional superconducting orders from repulsive interaction on doped AB-stacked honeycomb bilayer in $d$-wave, $f$-wave, and pair density wave channels were discussed in by James et al.James2014 Spin triplet $s$-wave pairing could also arise.Hosseini2012 Superconductivity reported in Ca-intercalated bilayer graphene represents the thinest limit of graphite intercalation compounds (GICs) Kanetani2012 , at 4KSatoru2016 and around 6.4K in Ca-doped graphene laminatesChapman2016 . Even so, the superconducting phase of AA-graphene has been addressed very little in the literatures very rarely.

Based on an effective Hamiltonian with an attractive potential between inter- and intra-layer near neighbor sublattices Alidoust et al. Alidoust2019 to study phonon-mediated superconducting pairing symmetries that may arise in AA, AB (and AC)-stacking bilayer graphene at the charge neutral point and beyond (by varying chemical potential). They claimed that at a finite doping, AB stacking can develop singlet and triplet $d$-wave symmetry beside $s$-wave, $p$-wave and $f$-wave that can be achieved at the charge neutral point, while the AA-stacked phase, similar to the undoped case, is unable to accept $d$-wave pairing.

Motivated by experimental observation of superconductivity in Ca doped bilayer graphene, a more realistic model has been introduced here to obtain all possible superconducting symmetries which can be arise analytically. In this manuscript, we follow the notion that Calcium doped bilayer graphene as mentioned in the ref.Kanetani2012 as the thinest limit of graphite intercalation compounds for which the structure Consists of Ca intercalated bilayer AA- stacked graphene (Fig.1(a)). However, recently, another possibility has been raised experimentally Endo2019 .

Based on mean-field treatment of an extended Hubbard model, a realistic tight-binding model has been used where its parameters are determined by a fitting to DFT band structure. By adding an effective attractive interactions between interlayer and intralayer electrons, all possible superconductivity pairing symmetry characters of C6CaC6 has been studied in details ( which can be applied to any related graphene-like structures such as B3N3CaB3N3). We will take advantage of the observation that, mathematically one can use mirror symmetry properties of Bloch coefficients of intercalated AA -stacked bilayer graphene and interpret their Hamiltonian as two independent single layer pseudo-graphene structures (even and odd sectors) where one of them (even symmetry sector) is decorated with calcium layer. This notion leads to the emergence of a topological number called cone index ($c=\pm 1$). Conservation of cone index during Klein tunneling across an n-p junction is one of the interesting unique behaviors in AA-stacked bilayer where it raises the possibility of ‘cone-tronic’ devices based on AA-BlG.Sanderson2013

In the real space a quasi-particle consist of two electrons from opposing layers with the same symmetrical position where they possess an additional quantum cone index beside their chirality nature near the Dirac points. The cone index concept is a unique feature in describing quasi-particles behaviors in AA-BlG with respect to single layer graphene. Two quasi particles (i.e. “four electrons”) can team up to build a Cooper pair as one can see in Fig. 1(b,c). We will see that only quasiparticles with the same cone index can be paired.

It will be shown that the question of superconducting phases in metal-intercalated bilayer graphene such as C6CaC6 can be decoupled to two independent gap equations corresponding to each of the even and odd sectors which they can be solved analytically (or nearly so) to obtain all possible pairing symmetry phases which can be probed experimentally.

The two gap nature of superconductivity that is one of the unique feature of MgB2 Putti2006 can be inspected here similarly. We have numerically predicted that in the temperature range of 0-6 K the phase transition from d-wave single-gap to dual-gap d-wave superconductivity could be observed. Using ab initio anisotropic Migdal-Eliashberg theory including Coulomb interaction, Margine et al.Margine2016 concluded that C6CaC6 should support phonon mediated superconductivity with a critical temperature $T_{c}=6-7$K, within the range of observations, and it exhibits two distinct superconducting gaps on the electron and hole Fermi surface pockets which is in agreement with the result has been obtained in the present manuscript.

The rest organization of the paper is as follows. Section II introduces the model Hamiltonian that we study, with Sec. III setting the stage by obtaining mostly analytic diagonalization of the non-interacting system. In Sec. IV, the treatment of pairing and presentation of superconducting phases is presented. In Sec. V, analytic insight are complemented by numerical solutions, followed by a discussion and summary in Sec. VI. Many of the analytic expressions are delegated to Appendices.

II Model Hamiltonian

The system we consider, illustrated in Fig.1(a), consists of $AA$ stacked bilayer graphene intercalated by Ca metal layer in which intercalant atoms are located on the central symmetry plane of bilayer graphene at the center of neighboring carbon hexagons. The distance between the graphene layers is calculated to be $h=4.63\text{\AA}$ in the case of Ca intercalation. The nearest in-plane Ca-Ca distance is $\xi=4.26\text{\AA}$. Charge transfer from Ca to the graphene layers leads to breaking symmetries of hopping amplitudes, and of C-C bond lengths similarly to those of Li decorated monolayer graphene[Rouhollah2018 ]. The attractive interaction between metal cations and C atoms after charge transfer contracts the Ca-C distance and reduces the C-C bond lengths in the Ca-centered hexagon to $a_{1}=1.419\text{\AA}$. As a result the bond length of neighboring C atoms in different hexagons is somewhat larger, at $a_{2}=1.423\text{\AA}$. Also in this “shrunken bilayer graphene”Rouhollah2018 the hopping integrals between short-bond inter- and intra-layer carbons are respectively $t^{11}_{1}$ and $t^{12}_{1}$, while those between stretched carbon sites will be denoted $t^{{}^{\prime}11}_{1}$ and $t^{{}^{\prime}12}_{1}$. The lattice then becomes a two dimensional hexagonal Bravais lattice with thirteen atomic sites. The sites of $i$-th cell will be labeled as $A_{i1}^{m},A_{i2}^{m},A_{i3}^{m},B_{i1}^{m},B_{i2}^{m},B_{i3}^{m}$ and $Ca$, where $m$ is layer index and takes $m=1,2$.

The Hamiltonian of this system is

[TABLE]

where $\alpha$ and $\beta$ run over the sublattice orbitals $A_{i}^{m}p_{z}$, $B_{i}^{m}p_{z}$ and Ca $s$. Here ${\hat{c}}_{i\alpha\sigma}^{\dagger}$, ${\hat{c}}_{i\alpha\sigma}$ are creation and annihilation operators of an electron with spin $\sigma$ on subsite $\alpha$ of $i$th lattice site, and $\hat{n}_{i\sigma}={\hat{c}}^{\dagger}_{i\sigma}{\hat{c}}_{i\sigma}$ is the electron number operator. The chemical potential is $\mu_{0}$ and $t_{i\alpha,j\beta}$ is the hopping integral from $\alpha$ subsite of $i$th site to the $\beta$ subsite of $j$th site. Here $U$ is an effective negative interaction between electrons in the extended (negative $U$) Hubbard model that allows the possibility of superconductivity.

III The non-interacting System

In this section a thirteen band tight binding model, consisting of twelve $p_{z}$ C orbitals and the Ca $s$ orbital, for Ca-intercalated bilayer graphene is constructed, to be applied to study superconducting states of this system within BdG theory. The non-interacting system Hamiltonian is invariant under mirror symmetry, which leads to division of the intercalated AA-BlG band structure into two sectors characterized by eigenvalues of mirror operation. Here we take advantage of the mirror symmetry. We first apply this reduction to the simple case of pristine AA-BlG.

III.1 Reducible Tight Binding Model for Pristine AA-staked Bilayer Graphene

The unit cell of AA-BlG, illustrated in Fig. 2(a), consists of four atoms $A_{1}$, $B_{1}$ in the top layer and $A_{2}$, $B_{2}$ in the bottom layer. The Schrödinger equation for this system in terms of Bloch coefficients is given by

[TABLE]

where $2\times 1$ column matrices $\chi_{1}$ and $\chi_{2}$ are components of AA-BlG iso-spinors, $H_{11}=H_{22}$ and $H_{12}=H_{21}$ are intralayer and interlayer hopping matrix respectively. The 4-component ”Dirac spinor” (Dirac representation) is reducible. Using the unitary transformation U=\frac{1}{\sqrt{2}}\left(\begin{array}[]{*{20}{c}}I_{2}&~{}~{}I_{2}&\\ I_{2}&-I_{2}&\end{array}\right), it decomposes into two irreducible representations, acting only on two 2-component right and left hand “Weyl spinors.” Similar transformations decouple the 4-component iso-spinor of AA-BlG into two 2-component chiral iso-spinors. This connection is an pedagogically useful mathematical similarity between the Schrödinger equation of AA-BlG and the Dirac equation, because it can give a insight into the prediction of quasiparticle behavior in AA-BlG in analogy with relativistic particles such as neutrinos. Mirror symmetry leads to

[TABLE]

Inserting Eq. 15 into Eq. 14 and defining new single layer graphene-like Hamiltonians $H^{\pm}=H_{11}\pm H_{12}$, the four band Schrödinger equation converts into two decoupled single layer graphene two band Schrödinger equations of the form

[TABLE]

wherein 2-component iso-spinors (has been shown in Fig. 2c) are given by

[TABLE]

The $\pm$ sign appearing in the even and odd sector Hamiltonians are the $h=\pm 1$ eigenvalues of the mirror symmetry operator \hat{S}_{h}=\left(\begin{array}[]{*{20}{c}}0&~{}~{}I_{2}&\\ I_{2}&0&\end{array}\right) which is the same as fifth Dirac gamma matrix $\hat{\gamma}^{5}$ that reflects the right ($\chi_{R}$) or left hand ($\chi_{L}$) chirality of quasiparticles in the relativistic quantum field theory. This additional $\pm$ topological index called by Sanderson and Ang (AS) as “cone index.”Sanderson2013 . SA have shown that quasiparticles in AA-BlG are not only chiral but are also characterized by “cone index,”. So we will refer to the cone index as h-chirality index and to usual helicity appears again as v-chirality. According to this notion, one can describe Dirac cones in AA-BlG shown in the Fig. 2(b) with two kind of chirality with respect to asymmetric in such a way that the structure and its vertical (“v-chirality”) and horizontal (“h-chirality”) mirror image are not superimposable. This chirality is a general aspect of AA-BlG quasi-particles that holds for general hoppings and over the entire Brillouin zone. SA show that electron transport across a barrier must conserve the cone index, a consequence of the Klein tunneling behavior in AA-stacked BlG. In the following sections we will extend the consequence of the cone index notion to superconductivity pairing.

A quasi-particle in the even sector (odd sector) consists of two electron (2$e$ charge) from opposing layers with the same spatial symmetry which possess $h=+1$ ($-1$) cone index and on-site energy of $\varepsilon^{+}=+\gamma_{0}$ ($-\gamma_{0}$) respectively. Hopping of a quasi-particle from A subsite to nearest neighbor B subsite in the even-sector (odd-sector) changes the energy by $t^{+}=t_{1}+\gamma_{1}$ ($t^{-}=t_{1}-\gamma_{1}$). This decoupling is more than just a simple mathematical diagonalization and symmetry characterization. One can interpret the Hamiltonian of AA-BlG as a single layer honeycomb lattice Hamiltonian with two types of charge carriers described by

[TABLE]

where $a_{i\sigma}^{\pm{\dagger}}$ ($b_{i\sigma}^{\pm{\dagger}}$) is the creation operator of a quasiparticle in the A(B) subsite of the $i$th site with $\sigma=\pm 1$ wherein $\sigma$ represents h-chiral Pseudo-spin (h-Pspin). where $a_{i\sigma}^{\pm{\dagger}}$ ($b_{i\sigma}^{\pm{\dagger}}$) is the creation operator of a quasiparticle in the A(B) subsite of the $i$th site with $\sigma=\pm 1$ wherein $\sigma$ represents $h$-chiral Pseudo-spin ($h$-Pspin).

In the intralayer and interlayer nearest neighbor hopping approximation, viz. $t_{1}$, $\gamma_{0}$ and $\gamma_{1}$ the system Hamiltonian is given by

[TABLE]

Here $t^{\pm}=t_{1}\pm\gamma_{1}$ and $t_{1}$ ($\gamma_{1}$) is the nearest neighbor intralayer(interlayer) hopping between A and B sublattices and the direct interlayer hopping (i.e. $A_{1}$ to $A_{2}$ or $B_{1}$ to $B_{2}$) is given by $\gamma_{0}$. Here $\mu$ is chemical potential, also $f(\vec{k})=\sum_{i=1}^{3}e^{i\vec{k}.\vec{\delta_{i}}}$ wherein $\vec{\delta_{i}}$ vectors connect A subsite to three in-plane nearest neighbor B subsites.

The decoupled Schrödinger equation Eq. 16 has eigenvalues

[TABLE]

where $\eta_{v-c}=\pm 1$. The four bands of AA-BlG separates into two independent up-down shifted single layer graphene bands where they are referred to as the even and odd sector, respectively.

Near the Dirac points, the dispersion energies $E_{n}^{\pm}(\vec{k})=\hbar v_{f}^{\pm}|\vec{k}|\pm\gamma_{0}$ of odd and even sectors are shown by two up-down shifted Dirac cones in Fig. 2(b) where $v_{f}^{\pm}=\frac{\sqrt{3}t^{\pm}a}{2\hbar}$ is Fermi velocity of $h=\pm 1$ quasi particles. Generalizing the tight binding model to include further neighbor hopping terms can highlight some hidden aspects of the AA-stacked Dirac cone quasi-particles. For second neighbor interlayer hopping, $\gamma_{1}$ taken into regard one can distinguish quasi-particles with the same chirality (v-chirality) and different cone index (h-chirality) from their velocities which could be inspected experimentally. The Fermi velocity of Dirac cone with $h=-1$ chirality decreases as interlayer hopping increases, while the velocity of $h=+1$ quasi-particle increases:

[TABLE]

In the strong inter layer coupling $t_{1}\to(-)\gamma_{1}$ Fermi velocity $v_{f}^{-}(v_{f}^{+})\to 0$ and odd(even) sector bandwidth tends to zero. As shown in the next subsection, interlayer coupling $\gamma_{1}$ may be considerable, so that inequality could be considerable.

III.2 Analytic Tight Binding Model for Intercalated Bilayer Graphene

In this subsection we will generalize the previous procedures to include the case of experimentally observed structures such as Li or Ca intercalated bilayer graphene. We will follow the notion it believe that these structure are intercalated AA-stacked bilayer graphene as has been shown in the Fig. 1. From the beginning, the Hamiltonian is generalized to incorporate several broken symmetries, including the on-site energies, hopping integrals, and bonds lengths (geometry). Due to this generalization, it can be used to obtain analytic dispersion energies of not only C6CaC6, but also related graphene-like structures such as B3N3CaB3N3. The Hamiltonian of such a non-interacting system is

[TABLE]

where $\alpha$ and $\beta$ run over sublattice orbitals $A_{i}^{m}$, $B_{i}^{m}$ and the intercalated atom (e.g. $Ca$) orbital. The Schrödinger equation for this system in terms of Bloch coefficients in $\vec{k}$ space becomes

[TABLE]

Here the $\beta$=0, 1, 2, …, 12 subscripts refer respectively to intercalant Ca, $A_{1}^{1}$, $A_{2}^{1}$, $A_{3}^{1}$, $B_{1}^{1}$, $B_{2}^{1}$, $B_{3}^{1}$, $A_{1}^{2}$, $A_{2}^{2}$, $A_{3}^{2}$, $B_{1}^{2}$, $B_{2}^{2}$ & $B_{3}^{2}$ and N is the number of unit cells. Here $\epsilon_{A_{i}}=\epsilon_{A}$ and $\epsilon_{B_{i}}=\epsilon_{B}$.

Mirror symmetry of this system result in the relations

[TABLE]

which reflects the mirror symmetry through the Ca plane that separates even and odd states (i.e. $h=\pm 1$ h-chirality). By inserting Eq. 34 into Eq. 33, with more detail given in Appendix A, one obtains two independent Schrödinger equations corresponding to eigenvectors ${\Psi_{i}^{-}}(\vec{k})^{T}=(0,~{}~{}~{}\bm{\mathcal{C}}_{i}(\vec{k}),~{}-\bm{\mathcal{C}}_{i}(\vec{k}))$, ${\Psi_{j}^{+}}(\vec{k})^{T}=({\mathscr{C}}_{0}(\vec{k}),~{}~{}\bm{\mathcal{C}}_{j}(\vec{k}),~{}+\bm{\mathcal{C}}_{j}(\vec{k}))$ respectively. For the odd eigensystem, the Schrödinger Eq.99 reduces to following $6\times 6$ matrix eigenvalue problem

[TABLE]

The Schrödinger Eq. 35 can be solved analytically, with the six eigenvalues presented in Appendix I, Eq. 150. These expressions are unaffected by the intercalant layer due to the separation of even and odd mirror symmetries, but the presence of Ca will renormalize parameters. For the even mirror sector, the Schrödinger equation Eq.99 reduces to the following $7\times 7$ matrix eigenvalue problem

[TABLE]

where $n=7,...,13$. The other seven bands of the Schrödinger equation Eq. 99 can be obtained from solving the new Schrödinger Eq. 44. The k-dependent part of corresponding matrix components of $H_{11}$ and $H_{12}$ are identical in form, differing only in hopping parameters, hence $H^{\pm}(\vec{k})$ can be considered as a shrunken graphene monolayer Hamiltonian with renormalized hopping parameters. It follows that, Similar to the pristine bilayer graphene (see Nori2016 ), intercalated bilayer graphene can be interpreted as two independent pseudo-graphene monolayers where one of them is dressed by a modified hopping Ca layer. The thirteen bands of bilayer graphene divide into two groups, six bands group (odd symmetries) corresponding to $H^{-}$ Hamiltonian and seven bands group (even symmetries) which are eigenvalues of $H_{c}^{+}$ matrix.

Mathematically many of the results obtained in ref. [Rouhollah2018, ] can be generalized to these graphene like structures but with renormalized hopping parameters. For general $\vec{k}$, except at $\Gamma$, it is challenging to obtain an exact analytical solution of the Schrödinger equations of Eq. 44. At the graphene Dirac points which have become folded back to the this supercell $\Gamma$ point, symmetry breaking of nearest neighbor intra- and interlayer hopping parameters i.e. $t_{11},~{}{t^{\prime}}_{11}$ and $t_{12},~{}{t^{\prime}}_{12}$ (Fig. 1b) results in two unequal small gaps with different centers, corresponding to six (odd-sector) and seven (even-sector) band pseudo-graphene Hamiltonians, given by

[TABLE]

where $t_{1}^{\pm}=t_{11}\pm t_{12}$ and $t_{1}^{\prime\pm}=t^{\prime}_{11}\pm t^{\prime}_{12}$. For the use of graphene as field effect transistors, it is necessary to create an tunable gap. Tunable and sizable band gap can be constructed in single layerFarjam2009 by decoration and in the bilayer graphe by intercalation as can be seen from Eq. 45.

The effect of symmetry braking of the the inter-layer coupling parameter i.e. $\Delta\gamma_{1}=t_{12}-{t^{\prime}}_{12}$ leads to inequality of even and odd sector gaps. Knowing the size of these energy gaps, one can find the difference in the first nearest neighbor intra and inter-layer hopping parameters symmetry breaking (suppose $\Delta t>0$, $\Delta\gamma_{1}>0$ and $\Delta t>\Delta\gamma_{1}$),

[TABLE]

where $\Delta t=t_{11}-{t^{\prime}}_{11}$. These two gaps are characteristic of AA-IBlG. In the case of Li- intercalated BlG, experimental ARPES spectra (Fig. 4 Ref.Caffrey2016 ) shows two distinct gaps of wide $E_{g}^{-}=0.20eV$ and $E_{g}^{+}=0.46eV$. In this reference the authors has equated the ratio of these two gaps as the ratio of the interlayer skew coupling parameters ($\gamma_{2}$ and $\gamma_{2}^{\prime}$ in their notation). Equation 46 slightly correct the discussion that has been stated in the [Sec. III, sub-sec.C of Ref.Caffrey2016 ] about this ratio.

Eq. 44 results in the event that intercalant layer hopping parameters to graphene sheets are negligible, as in the case of Li-decorated graphene where Li atoms fully ionize and the Li-associate band lies above the Fermi energy, so Li-C hopping effects are negligible. In this particular case the odd (-) and the even (+) nontrivial eigenvalues of $H^{\pm}$ matrix, are given by

[TABLE]

wherein $\mu_{m}^{\pm}(\vec{k})$, $w_{m}^{\pm}(t_{i},\vec{\xi_{i}},\vec{k})$ are defined in Eqs. 151 and 160 respectively. Details for obtaining these results are presented in Appendix A and ref. [Rouhollah2018, ].

Now further notation is established. Similar to the case of pristine AA-BlG investigated in the previous subsection, this transformation recast the noninteracting Hamiltonian Eq. 32 as the direct sum of two single layer pseudo-graphene structures with renormalized hopping integrals in the right and left hand chiral representation of the form

[TABLE]

wherein, as illustrated in the Fig.2(c), we introduced quasiparticle creation operators and hopping integrals in real space,

[TABLE]

The creation operator $\hat{c}^{(\pm)\dagger}_{i\alpha\sigma}$ creates an quasi-particle with $h=\pm 1$ h-chirality and spin ($\sigma$) at $i\alpha$ atomic subsite of each of these two graphene-like structures, participate directly in the formation of superconducting phases.

The separation of thirteen bands of intercalated bilayer graphene into groups of six bands with odd-symmetry and seven bands with even symmetry has strong advantages. The tight-binding model band structures can be fit to DFT band structure results with greater accuracy and simplicity, for example. Especially when the pairing interactions are introduced, this transformation reduces the speed of numerical calculations significantly and provides additional insight into physical properties of bilayer graphene.

III.3 Fit to DFT band structures

This formalism has been applied and divided into two separated effective single layer (shrunken)- pseudo-graphene models. DFT calculation is used to obtain the electronic structure data. Gnu-plot of DFT data has been shown with blue thin dashed lines in the background of Fig. 3(a),(b). The six odd bands and seven even bands were fit to DFT bands with results shown with color lines in the Fig. 3. The main problem that emerged here is that odd and even sector bands are not separated in the DFT data. But by inspecting the DFT band structure and to be careful in analytical calculations, knowing that odd sector does not affected by Ca-C coupling the odd and even sector flat bands can easily be distinguished. Emerging of two distinct gaps in the Dirac point of both sectors is the other guide to perform fitting. The reduced fitting parameters are given in Tables 1 and 2.

We follow the model presented in ref. [Rouhollah2018, ]. As illustrated in Fig. 1(a) on each of bilayer graphene sheets, the $A_{1}$ sublattice site of the central unit cell is chosen as the origin labeled by [math], and the $B_{1}$ site in the adjacent hexagon is considered as the second C atom neighbor. While just slightly longer than the nearest neighbor atoms $B_{2}$ and $B_{3}$ in the same hexagon, this neighbor is labeled by $n=2$, and so on the further neighbors are labeled. In Fig. 1 (a), the big hexagon included up to nine intra-plane neighbors but for the pristine graphene it is surrounded by five neighbors. C-C hopping from [math]-subsite to intra-plane $n$th neighbor($t_{i0jn}^{intra}$) plus (minus) hopping from [math]-subsite to inter-plane $n$th neighbor ($t_{i0jn}^{inter}$) has been shown by $t^{{\pm}CC}_{0n}$. In-plane Ca-Ca hoppings $t^{Ca-Ca}_{0m}$ are included up to $m=4$ neighbors. Modified Ca to C hopping integrals in Eq. 44, which are defined as $\sqrt{2}$ times the hopping from central Ca to $m$th neighbor C atoms, are denoted by $t^{CaC}_{0m}$ and obtained up to $m=5$ neighbors.

The six odd bands and seven even bands specified by Eqs. 35 and 44 respectively, have reduced hopping integrals given by $t^{\pm}_{im\sigma,jn\sigma}=t^{inter}_{im,jn}\pm t^{intra}_{im,jn}$, DFT calculated bands were fitted to tight binding odd bands of Eq. 47, with results presented in Fig. 3(a) and Table I. The even bands which are solutions of Eq. 44 are obtained numerically and fitted to the DFT bands. The results are illustrated in Fig. 3(b) and Table II.

There are two flat bands with $d$-wave Bloch character: one in each of the odd and even sectors. The opposite signs of nearest neighbor inter- and intra-layer hopping amplitudes, $t_{1}^{11}$ and $t_{1}^{12}$, leads to reduced bandwidths of the even states (+ sign) in Eq. 44. Larger interlayer hopping $t_{1}^{12}$ leads to smaller bandwidth ($H^{\pm}=H_{11}\pm H_{12}$), while for the other six odd-symmetry bands, the bandwidth can increase as can be seen in Fig. 3. On the other hand, the bandwidth of the even sector flat band is reduced, again due to the calcium to carbon hopping while the odd sector is not affected by Ca-C hoppings. For this reason, the flat band belonging to the even bands group plays a major role in superconductivity.

IV Superconducting Pairing and States

IV.1 Bogoliubov-de Gennes Transformation

We treat the thirteen band Hubbard model in mean field approximation to investigate superconductivity in intercalated bilayer graphene. Singlet pairing is considered and, as illustrated in Fig.  1, pairing interactions are pictured in real space as interactions between nearest neighbors on inter- and intra-layer carbon atoms. This superconducting Hamiltonian can be transformed, as for the non-interacting case, to the direct sum of two independent superconducting Hamiltonian corresponding to odd and even symmetries pseudo-graphene structures

[TABLE]

where $\Lambda^{\dagger}(\vec{k})=\left([\hat{c}^{\dagger}_{0\uparrow}(\vec{k})\hat{c}^{+\dagger}_{1\uparrow}(\vec{k})...\hat{c}^{+\dagger}_{6\uparrow}(\vec{k})\;\hat{c}_{0\downarrow}(-\vec{k})\hat{c}_{1\downarrow}^{+}(-\vec{k})...\hat{c}_{6\downarrow}^{+}(-\vec{k})]~{}~{}~{}[\hat{c}^{-\dagger}_{1\uparrow}(\vec{k})\;\hat{c}^{-\dagger}_{2\uparrow}(\vec{k})...\hat{c}^{-\dagger}_{6\uparrow}(\vec{k})\;\hat{c}_{1\downarrow}^{-}(-\vec{k})\;\hat{c}_{2\downarrow}^{-}(-\vec{k})...\hat{c}_{6\downarrow}^{-}(-\vec{k})]\right)$, in which $\hat{c}^{(\pm)\dagger}_{m\sigma}(\vec{k})=\frac{1}{\sqrt{2}}(\hat{c}^{\dagger}_{m\sigma}(\vec{k})\pm\hat{c}^{\dagger}_{m+6,\sigma}(\vec{k}))$ and ${H_{su}^{+}}$ and ${H_{su}^{-}}$ are Hamiltonians of even and odd symmetry pseudo-graphene structures respectively; for more information see Appendix C. Decoupling of these Hamiltonians means there is no effective pairing between an electron in the even sector with one in the odd sector. Using the fact that the gap is small on the electronic scale, applying perturbation up to second order gives quasiparticle energies from Eq.52 (see ref. [Rouhollah2018, ]) as

[TABLE]

[TABLE]

Here $(\Delta^{1}_{\pm}~{}~{}\Delta^{2}_{\pm}~{}~{}\Delta^{3}_{\pm})=[(g_{1}\Sigma_{1}\pm g_{1}^{{}^{\prime}}\Sigma_{1}^{{}^{\prime}})~{}~{}(g_{1}\Sigma_{2}\pm g_{1}^{{}^{\prime}}\Sigma_{2}^{{}^{\prime}})~{}~{}(g_{1}\Sigma_{3}\pm g_{1}^{{}^{\prime}}\Sigma_{3}^{{}^{\prime}})]$, $(\Delta^{4}_{\pm}~{}~{}\Delta^{5}_{\pm}~{}~{}\Delta^{6}_{\pm})=[(g_{0}\Delta_{1}\pm g_{0}^{{}^{\prime}}\Delta_{1}^{{}^{\prime}})~{}~{}(g_{0}\Delta_{2}\pm g_{0}^{{}^{\prime}}\Delta_{2}^{{}^{\prime}})~{}~{}(g_{0}\Delta_{3}\pm g_{0}^{{}^{\prime}}\Delta_{3}^{{}^{\prime}})]$ and $(\Delta^{7}_{\pm}~{}~{}\Delta^{8}_{\pm}~{}~{}\Delta^{9}_{\pm})=[(g_{0}\Pi_{1}\pm g_{0}^{{}^{\prime}}\Pi_{1}^{{}^{\prime}})~{}~{}(g_{0}\Pi_{2}\pm g_{0}^{{}^{\prime}}\Pi_{2}^{{}^{\prime}})~{}~{}(g_{0}\Pi_{3}\pm g_{0}^{{}^{\prime}}\Pi_{3}^{{}^{\prime}})]$, where $<ij>$ subscript has been dropped for brevity. Also

[TABLE]

Here ${\mathscr{C}_{m}}({E_{n}^{-}})$ is the $m$th component of $n$th column eigenvector of $H^{-}$ and ${\mathscr{C}_{m}}({E_{n}^{+}})$ is $(m+1)$th component of $n$th eigenvector of $H_{c}^{+}(\vec{k})$. Band order parameters $\Delta_{mn}^{+}(\vec{k})$ are defined such that first electron is in the $m$th band and second electron is in the $n$th band of $H_{c}^{+}(\vec{k})$, also $\Delta_{mn}^{-}(\vec{k})$ is defined such that first electron is in the $m$th band and second electron is in the $n$th band of $H^{-}(\vec{k})$. Note that an electron in the $m$th band of $H_{c}^{+}(\vec{k})$ and an electron in $n$th band of $H^{-}(\vec{k})$ cannot be paired; i.e for this case $\Delta_{mn}^{\pm}(\vec{k})=0$.

The Bogoliubov-de Gennes transformation used in Eq. 52 shows that pairing amplitudes should be $\Delta_{\pm}^{\alpha}=<\hat{c}^{\pm}_{\alpha,i}\hat{c}^{\pm}_{\alpha,j}>$ which implies that all inter- and intra-layer pairing amplitudes in real space are equal, $g_{0}=g_{0}^{{}^{\prime}}$ and $g_{1}=g_{1}^{{}^{\prime}}$. This restriction makes the matrix gap equations hermitian and implies that band order parameters, $\Delta_{mn}^{\pm}(\vec{k})$ can be interpreted physically as pairing of electrons in different bands with pairing interaction $g_{0}^{\pm}$. In this limit $\Delta_{mn}^{\pm}(\vec{k})$ is equal to the product of band Green function and $g_{0}$,

[TABLE]

where $\hat{d}^{\pm\sigma}_{i}(\vec{k})=\sum^{7}_{m=1}{\mathscr{C}}^{\pm*}_{m}(E_{i}(\vec{k}))\hat{c}^{\sigma}_{m}(\vec{k})$ annihilates an electron with spin $\sigma$ in the $i$th even or odd bands with energy $E_{i}^{\pm}(\vec{k})$.

IV.2 Two Gap Superconducting Pairings and States

The linearized gap equation can be decoupled by minimizing free energy with respect to nearest neighbor pairing, or equivalently with respect to $\Delta_{\pm}^{\alpha}$, for more detail see Appendix C. Minimization of free energy with respect to $\Delta_{\pm}^{\alpha}$ gives

[TABLE]

in which $A$, $B$, $C$ and $D$ matrices have been introduced as

[TABLE]

wherein, $\Gamma$ matrix elements are given by

[TABLE]

Equation 66 can be interpreted as two independent gap equations for odd (minus sign) and even (plus sign) pseudo-graphene systems. The impact is that superconductivity can be established independently in two distinct sectors this system. In the next section we numerically inspect which of these pseudo-graphene sectors, odd or even, play major rules of superconductivity.

Gholami et al.Rouhollah2018 solved such gap equations for Li-decorated single layer graphene. The gap equations in Eq. 234 are the same form so they can be solved similarly. The $A$, $B$, $C$, and $D$ matrices have identical structures, hence they share eigenvectors: $V_{s}^{T}=(1\;\;1\;\;1)$, $V_{d_{xy}}^{T}=(1\;\;-1\;\;0)$, and $V_{d_{x^{2}-y^{2}}}^{T}=(1\;\;1\;\;-2)$, where the latter two are degenerate. Their eigenvalues, in obvious notation, are

[TABLE]

Similar to decorated single layer graphene Rouhollah2018 , for each of the gap equations given by Eq. 234 there are nine independent solutions. The first three superconducting states with island (localized) character can be expressed in compact form as

[TABLE]

where $V_{sy}$ refers to one of the $V_{s}$, $V_{d_{xy}}$ or $V_{d_{x^{2}-y^{2}}}$-wave symmetries. Pairing in these phases cannot propagate. The other six superconducting states of Eq. 234 have the explicit form

[TABLE]

where

[TABLE]

Here $c_{sy}=c_{s}^{\pm},c_{d}^{\pm}$, $b_{sy}=b_{s}^{\pm},b_{d}^{\pm}$, $d_{sy}=d_{s},d_{d}^{\pm}$ and $J_{sy}^{l,\pm}=-\frac{1}{g_{0}^{\pm}}$ for each symmetry, and $+$ superscript refer to the even sector and $-$ superscripts to the odd sector. In each of above categories, ${d_{{x^{2}}-{y^{2}}}}$ and ${d_{xy}}$ phases are degenerate. Similar to decorated single layer graphene, only three of solutions for which $l=2$ are physically reachable in the framework of mean field theory. In the limit of pristine bilayer graphene, these three states convert to usual s-wave and d-wave symmetries. In the sec. V, for odd and even sectors we illustrate from numerical solutions which of these three phases are dominant.

IV.3 Flat band(s) Superconductivity: Strong interlayer coupling

To make a rough estimate and provide mathematical insight into the physics, one can diagonalize normal state Hamiltonian of pristine bilayer graphene in the mini-Brillouin zone of C6CaC6. As shown in Fig. 3, two conduction bands corresponding to odd and even sector are weakly dispersive near the Fermi energy along $\Gamma\to M$, which seems to play a major role in the formation of superconducting Cooper’s pairs.

In the case of pristine bilayer graphene, odd (-) and even (+) so called flat bands are the minimum of ($E_{\alpha,2}^{\pm}(\vec{k}),E_{\beta,2}^{\pm}(\vec{k})$) along different high symmetry paths that are given by Eqs. 146 and 147, where their Bloch wave function, viz. Eq. 148, are similar to those of ref. Rouhollah2018 . They have linear combination of $d_{x^{2}-y^{2}}$ and $d_{xy}$ character and are responsible for superconducting pairing $d_{x^{2}-y^{2}}$ and $d_{xy}$.

One can ask: what is so special about these flat bands? To address this question, we return to the matrix gap equation of Eq. 80. The right hand side contains the product of a form factor given by $\Omega_{ni}^{\alpha}(\vec{k})\Omega_{ni}^{*\beta}(\vec{k})+\Omega_{ni}^{\beta}(\vec{k})\Omega_{ni}^{*\alpha}(\vec{k})$ and the thermal occupation factor over the energy denominator i.e. $\frac{\tanh(\frac{{E_{n}}}{2k_{B}T})}{E_{n}(\vec{k})+E_{i}(\vec{k})}$. The form factor is a function of the Bloch wave coefficients of normal state Hamiltonian. Using Eqs. 148, 149 one can investigate that in the limited case of pristine bilayer graphene at the nearest neighbor approximation, these Bloch wave coefficients are the same for both sectors and this is almost for the next neighbor approximation. As such, it is independent of chemical potential $\mu$. Thus the form factor is the same for the both odd and even sector of band structures. Since $\frac{tanh(x)}{x}\to 1$ as $x\to 0$, when one of the conduction odd or even flat bands and so their corresponding valance bands becomes completely flat at the Fermi level then $\frac{\tanh(\frac{{\beta E_{i}}}{2})}{E_{i}(\vec{k})+E_{j}(\vec{k})}\to\frac{\beta}{4}$, where $E_{i}(\vec{k})$ or $E_{j}(\vec{k})$ are one of these flat bands. In this case the dominant contribution comes from these mutual conduction and valence flat bands, and one can show that all gap equation block matrix elements in Eq. 79 are equal to $A^{\pm}.$ In this event, depending on whether the flat bands belong to the odd or even sector, one can use Eqs. 55, 80, and 148 to show that

[TABLE]

Cooper pair interaction potentials $g_{0}$ of d-wave symmetry, i.e. $g_{0}^{d}$ and s-wave symmetry $g_{0}^{s}$, are given by

[TABLE]

In this case $\Gamma_{12}^{\pm}<0$ and $g_{0}^{d}$ is less than $g_{0}^{s}$, so $d$-wave symmetry is dominant, with an extraordinary decrease in pairing potential interaction proportional to the critical temperature. This “ultra” decrease of pairing interaction can explain the importance of the flat bands in the formation of Cooper pairs in twisted bilayer graphene. Here, another point that can be deduced from mathematical calculations is that in the limit of strong interlayer hopping when inter-layer hoppings tends to minus (plus) of intra-layer hoppings, as one can see from Appendix D, all of the six bands of even (odd) sector become flat while the other sector bandwidths increases. Then one can show that the gap matrix elements are

[TABLE]

and so all possible superconducting symmetries are degenerate with pairing potential $g_{0}=k_{B}T_{c}$.

V Numerical Results

V.1 General features

To know in the variety of doping regimes which of the pairing symmetries (distorted $s$-wave or $d$-wave) are dominant, and also to inspect in the which sectors of the band structure Cooper pairs with the lowest pairing potential can constructed, superconducting gap equations of odd and even sectors i.e. Eq. 66 are solved numerically. The result is shown in Fig. 4.

Similar to Li intercalated single layer graphene (ref. [Rouhollah2018, ]), at moderate doping, $d$-wave superconductivity always dominates in both sectors of the C6CaC6 band system. Distorted $s$-wave symmetry only survives at high doping levels. For odd and even sector flat bands, the density of states peaks at the $M$ critical point. This point, for the odd flat band, is located about 0.5 eV above the Fermi level, and about 0.44 eV below the Fermi level for the case of the even flat band. Exploration of each of these flat bands can be engineered by applying a gate voltage on the bilayer via a change in the chemical potential $\mu$. As illustrated in Fig. 4, around $\mu=0$ at $T_{c}=1K$ i.e. C6CaC6 is not affected by gate voltage. Dominant $d$-wave symmetry pairing phases are degenerate, so pairing can arise in both odd and even flat bands. This occurs when the averaged density of states of both bands are the same near the Fermi level. When the odd flat band is dominant near Fermi level due to electron doping, $d$-wave pairing dominates in this band, with the minimum of the pairing interaction energy ($g_{0}^{-}=0.35$) corresponding to the $M$ point near $\mu=0.5eV$.

V.2 Hole doping

Hole doping by gating, by Ca$\rightarrow$Na substitution, or by Ca deficiency, leads to the situation where the even flat band dominates. Dominant $d$-wave symmetry occurs in this band with the lowest pairing energy ($g_{0}^{+}=0.12$), again at the critical $M$ point at $\mu=-0.44eV$. From Fig. 4 it can be seen that when the even flat band is dominant (hole doping), the pairing potential $g_{0}$ of $d$-wave symmetry emerging from this band is less than the case that dominant $d$-wave symmetry occurs in the odd flat band (by electron doping). For example, with a factor of one-third at their critical $M$ point i.e. $g_{0}^{+}=\frac{1}{3}g_{0}^{-}$, that means that when the even flat band reaches near the Fermi level, reduction of bandwidth due to both interlayer C-C interaction($H_{12}$) and C-Ca layer interactions lead to a sharp increase in the density of states. While both C-C and Ca-C interlayer interaction decrease the pairing potential $g_{0}$, one can numerically inspect that reduction of the bandwidth due to graphene interlayer interaction, more affected the energy of pairing in the even flat band than Ca-C layers interaction.

In the case of Ca intercalated bilayer graphene, for a given $T_{c}$ the pairing interaction potential $g_{0}$ (proportional to superconducting gap energy, $|\Delta|^{2}$) for dominant $d$-wave phases of the odd and even superconducting gaps are illustrated in Fig. 5. It is evident that for a given critical temperature $T_{c}$, superconductivity can be single gap or two gap and dominant superconducting pairings can occur between electrons in the odd $H^{-}$, or even $H^{+}_{c}$, sector.

VI Discussion and Summary

Discovery of new superconducting phases, often at low temperature, has been one of the active achievements in recent decades. Superconductivity in the lithium-coated single layer graphene with respect to to the case of calcium-decorated single layer is more capable Ludbrook2015 and also reported experimentally, whereas this situation is opposing in the bilayer graphene. Superconductivity has been reported in the Ca-intercalated bilayer graphene around $T_{c}=4K$ while Li-intercalated bilayer graphene is not superconductor Satoru2016 ; Chapman2016 . Experimental fabrication of Li and Ca-intercalated bilayer graphene has been reported in the ref.Caffrey2016 and ref.Kanetani2012 respectively. Li and Ca atoms are suggested to intercalate between graphene layers with an ordered structures similar to that of bulk GICs like LiC6 i.e. two graphene layers are AA-stacking. That it demonstrates the vitally important role of intercalant inter-band in the formation of superconductivity Cooper pairs.

Most theoretical microscopic models of pristine honeycomb bilayer superconductivity concentrate on the more stable AB stacking of bilayer graphene.Vucicevic2012 ; James2014 . To our knowledge, there are few studies that focus on pristine AA stacked bilayer graphene, Alidoust2019 and no analogous studies that concentrate on intercalated bilayer graphene. Based on ab initio calculations of electron-phonon coupling, anisotropic Migdal-Eliashberg theory has been applied by some authors Margine2016 ; Durajski2019 to give strong evidence that this system is a phonon mediated two-gap superconductor with predicted $T_{c}$ around 7K. Recently unconventional superconductivity up to $T_{c}=1.7$K has been reported in gated twisted bilayer graphene where the layers are rotated relative to each other by a magic angle of $1.1^{\circ}$. Superconductivity in this low doping regime of band filling cannot be addressed within the framework of conventional electron-phonon coupling based on Migdal’s adiabatic approximation. This discovery has opened speculation that this superconducting behavior may shed light on other systems in which superconductivity arises from an insulating phase.Cao2018 This development also highlights studies such as ours, which does not rely on the mechanism, but instead on more general pairing concepts and the specific electronic structure.

Many theorists have suggested that exotic superconductivity gaps arise in some materials are related to a peculiarity of the normal-state band structure. In this kind of issue angle-resolved photoemission spectroscopy (ARPES) extensively has been applied to analysis of the normal state band structure. To determine structural and electronic properties of material, tight binding model in addition to DFT calculation has been used to interpret experimental results achieved from ARPES. Following this point of view, extended Hubbard model has been used here to address superconductivity of Ca- intercalated bilayer graphene.

The main results are achieved in two steps: first, for the normal-state (non-interacting part) a more realistic effective tight-binding model with two decoupled symmetry sectors is derived, with the parameters determined by a fitting to the DFT band structure; second, the dominant superconducting pairing channels are discussed based on a mean-field treatment of a Hubbard model obtained by adding (effective) attractive interactions between the electrons. The summary, results and comparisons are presented below.

VI.1 Non interacting Part: Normal state

In the first part of this manuscript we have taken the advantages of mirror symmetry operation on the AA-stacked BlG through the central plane and generalized it to include intercalated bilayer graphene.

VI.1.1 AA-stacked pristine bilayer graphene: Two kinds of quasiparticles

The honeycomb lattice structure makes quasi-particles in the single layer graphene behave as a massless Dirac particles at low energies that provide a proper platform to examine the characteristic effects of QED, such as the Klein paradox and Zitterbewegung, which were never observed in particle physics. In addition to relativistic nature of quasiparticles in single layer graphene, they exhibit extra aspects of such behaviors in AA-stacked BlG. Interlayer coupling causes the bilayer graphene to exhibit properties that are not observed in the single layer graphene. In QED, the 4-component ”Dirac spinor” (Dirac representation) decomposes into two irreducible representations, acting only on two 2-component right and left hand “Weyl spinors.” There is an pedagogically useful mathematical similarity between the Schrödinger equation of AA-BlG in these two representations and the Dirac equation. The non-interacting AA-stacked Hamiltonian is invariant under mirror symmetry, leads to division of the AA-BlG band structure into two even and odd sectors characterized by eigenvalues of mirror operation $h=\pm 1$ (analogous to two decouple Weyl equations for massless relativistic chiral particles). Each of these sectors describes a graphene-like structure i.e. $H^{+}$ and $H^{-}$. In this notion as has been shown in the Fig. 2c, up and down pseudo-spin (h-Pspin) of irreducible blocks of AA-BlG Hamiltonian viz. $H^{+}$ and $H^{-}$, consists of two electrons with the same spin which each one located at similar sub-sites in the opposing layers. These quasi-particles describe by an additional index has been called “cone index”. Here we refer to this index as h-chirality index. According to this notion, one can describe Dirac cones in AA-BlG shown in the Fig. 2(b) with two kind of chirality with respect to asymmetric in such a way that the structure and its vertical (“v-chirality”) and horizontal (“h-chirality”) mirror image are not superimposable. This chirality (h-chirality) is a general aspect of AA-BlG quasi-particles that holds for general hoppings and over the entire Brillouin zone and it is unrelated to the helicity operation. This is in contrast to the famous graphene chirality (helicity) that occurs just at low energies near the Dirac cones.

Physically, AA-stacked BlG can be interpreted as a “single layer honeycomb lattice” that instead of $1e^{-}$ charge carriers, there are two types of fermionic quasi-particles with $2e^{-}$ charge, moving through it that differ in a quantum number called cone index (h-chirality). Quasi-particles with different $h$-chirality don’t interact but move independently. Also $(\pm 1)$ $h$-chiral quasi-particles have $\pm\gamma_{0}$ on-site energies (similar to the positive and negative energy of the particles and the anti-particles in the QED). Hopping of quasi-particles with (+1) $h$-chirality constructs the even sector of band structure while the odd sector made by (-1) $h$-chiral quasiparticles. Near the Dirac cone points, quasi-particles with ($\pm 1$) $h$-chirality are moving with Fermi velocities $v_{f}^{\pm}$. One can distinguish quasi-particles with the same chirality (v-chirality) and different cone index (h-chirality) from their velocities which in the case of strong interlayer coupling could be observed experimentally.

VI.1.2 Intercalated AA-staked bilayer graphene: Mirror symmetry operation advantage

Based on mean-field treatment of an extended Hubbard model, a realistic thirteen band tight binding model, has been constructed to include the case of experimentally observed structures such as Ca intercalated bilayer graphene, where its parameters are determined by a fitting to DFT band structure. We followed the notion that Calcium doped bilayer graphene as the thinest limit of graphite intercalation compounds (Fig.1(a)).

In our previous work, the effects of Li-decoration on structural and electronic band structure of single layer graphene has been demonstrated in details and also symmetry character of the band-branches illustrated as well as the possible superconducting phases of lithium decorated single layer graphene LiC6 were obtained analytically and analyzed.Rouhollah2018 . The Brillouin zone (BZ) of these structure is one third of that of graphene, with the Dirac points folded back to the $\Gamma$ point. In this mini-BZ, the two $\pi$ bands of (pristine) graphene folds to six branches and their different symmetries (d+id, s,…) are also separated as illustrated in Fig. 2 of ref. [Rouhollah2018, ]. Generalization of these results to include intercalated bilayer graphene (IBlG) has strong advantages and provides additional insight into its physical properties. This is possible through decoupling of normal and superconducting Hamiltonians of IBlG into two independent corresponding single layer pseudo-graphene Hamiltonians, coupled only by a common chemical potential.

Similar to the pristine AA-staked BlG, accounting for symmetries of Bloch wave coefficients, the $13\times 13$ Hamiltonian of IBlG converts, by mirror symmetry, to two decouple sectors: an $7\times 7$ even symmetry sector $H_{c}^{+}$ with involvement of the intercalant (coated single layer pseudo-graphene) and the $6\times 6$ odd sector $H^{-}$, for which the intercalant provides only renormalized hopping amplitudes and break down symmetry of hopping integrals (six bands shrunken single layer pseudo-graphene). Therefore, all previous discussions about $2e^{-}$ charge, h-Pspin and chirality of quasiparticles in pristine AA-stacked BlG are extended to include IBlG.

Periodic perturbation of graphene layers potential due to ordered Intercalant atoms causes hopping integrals symmetry to break and so two distinct gaps of size $E_{g}^{\pm}=2|t_{11}^{\pm}-{t^{\prime}}_{11}^{\pm}|$ open at the Dirac point (folded to $\Gamma$ point) of each of the even and odd sector pseudo- graphene structures. These two gaps are characteristic of AA-IBlG. Knowing the size of these energy gaps, one can find the difference in the first nearest neighbor intra and inter-layer hopping parameters symmetry breaking i.e. ($\Delta t=t_{11}-{t^{\prime}}_{11}$) and ($\Delta\gamma=t_{12}-{t^{\prime}}_{12}$). In the case of Li- intercalated BlG, experimental ARPES spectra (Fig. 4 Ref.Caffrey2016 ) shows two distinct gaps of wide $E_{g}^{-}=0.20eV$ and $E_{g}^{+}=0.46eV$. We slightly correct the discussion that has been stated in the [Sec. III, sub-sec.C of Ref.Caffrey2016 ] about the relation between these two gaps and symmetry breaking interlayer coupling parameters.

In the case of Li- intercalated BlG, Li-s orbital is fully ionized and Li-C hybridization is negligible, so the odd and even sector band structure are similar to the band structure of pristine shrunken- graphene $C_{6}$. The difference of course is due to $\pm$ sign that appears in the even and odd sectors between intra and inter layer hopping terms which leads to different bandwidth. The even and odd sectors Schrödinger equation solved analytically (or nearly so). From the beginning, the Hamiltonian is generalized to incorporate several broken symmetries, including the on-site energies, hopping integrals, and bonds lengths (geometry). Due to this generalization, it can be used to obtain analytic dispersion energies of not only C6CaC6, but also related graphene-like structures such as B3N3CaB3N3.

VI.1.3 Tight Binding Parametrization of Ca- intercalated bilayer graphene from DFT

Dividing the thirteen bands into seven even-symmetry bands and six odd-symmetry bands, considerably facilitates tight Binding Parametrization from DFT. We used up to nine neighbor approximation tight binding model considering the symmetry breaking of bond length and hopping integral parameters across different direction of hexagons. The main problem with DFT data is that the odd and even sector data bands are not separated. But by inspecting the DFT band structure and to be careful in analytical calculations (knowing that odd sector does not affected by Ca-C coupling directly and so remain graphen like), the odd and even sector DFT flat bands can easily be distinguished. Emerging of two distinct gaps in the Dirac point of both sectors is the other guide to perform fitting. The reduced fitting parameters are given in Tables 1 and 2 with results are shown in the Fig. 3.

VI.2 Interacting Part: Superconductivity

SA show that electron transport across a barrier must conserve the cone index, a consequence of the Klein tunneling behavior in AA-stacked BlG. Here it is discussed that Cone index footprint can be traced also in the formation of superconducting Cooper pairs. Due to this index, the salient differences are emerged between the Cooper pairs in single layer graphene and AA-stacked bilayer graphene. The two types of even and odd symmetry superconductivity are predictable in AA-stacked bilayer graphene.

VI.2.1 Odd/ Even Superconducting Gap equations and Symmetry Phases

Similar to the normal state Hamiltonian, the superconducting Hamiltonian is also block diagonalized into two sectors. Each of these two blocks represents the superconducting Hamiltonian of the even and odd single layer pseudo-graphene structures. The impact is that superconductivity in AA-stacked BlG can be established independently in two distinct band structure sectors. Two quasi particles (i.e. “four electrons”) just with the same h-chirality (i.e. cone index) can team up to build a Cooper pair (Fig. 2c). In the other words, pairing in bilayer graphene arises between quasiparticles inside of the coated single layer pseudo-graphene $H^{+}_{c}$ band structure (even sector-superconductivity) or uncoated shrunken single layer pseudo-graphene $H^{-}$ (odd sector-superconductivity) separately; even-odd pairing is impossible without further symmetry breaking.

Two distinct superconductivity gap equations corresponding to $H_{c}^{+}$ and $H^{-}$ single layer pseudo-graphene structures emerged from minimization of the free energy. These behaviors show that general aspects of superconductivity in the (Li-)decorated single layer grapheneRouhollah2018 and (Ca-) intercalated bilayer graphene are similar, and their behaviors are different primarily in the probability of two gap superconductivity in the bilayer structures. A difference of course is that interlayer interaction becomes a key factor; the Cooper pairs in AA-stacked BlG instead of $2e^{-}$ charge have $4e^{-}$ charge and additionally have an right or left hand h-chirality index. This theoretical prediction requires empirical inspect.

These two even and odd superconducting gap equations were solved analytically to obtain the relations between the superconducting pairing potential and resulting ordered phases. The two sets of gap equations have solutions similar to those obtained in our previous work, for decorated single layer graphene.Rouhollah2018 Seven hybridized orbitals in pseudo- coated single layer graphene support nine possible bond pairing amplitudes. There are nine superconducting phases with $p_{x}$, $p_{y}$, $f$, $s^{\pm}$, $d_{xy}^{\pm}$ and $d_{x^{2}-y^{2}}^{\pm}$ atomic orbital-like symmetries corresponding to each of these even(+)/odd(-) gap equations. Only three of them are physically reachable, denoted by $\Psi^{{\pm}}_{2,s}$, $\Psi^{{\pm}}_{2,{d{xy}}}$, and $\Psi^{{\pm}}_{2,{d_{x^{2}-y^{2}}}}$. These symmetries almost preserve properties from a two band model of pristine graphene (Fig. 4 Ref.Rouhollah2018 ). The $d$-wave solutions are degenerate and so it can support chiral $d_{x^{2}-y^{2}}+id_{xy}$ superconductivity in each of these sectors. These three phases are distorted by intercalation. In fact, the significant difference which appears between two bands pristine $C_{2}$ pairing symmetries and shrunken graphene $C_{6}$(decorated graphene) is a skewness factor i.e. $\alpha_{sy}^{l,\pm}\neq 1$ in front of the self consistent gaps solutions. In each of these even or odd sectors, one band is weakly dispersive near the Fermi energy along $\Gamma\to M$ where its Bloch wave function has linear combination of $d_{x^{2}-y^{2}}$ and $d_{xy}$ character, and is responsible for $d_{x^{2}-y^{2}}$ and $d_{xy}$ pairing with lowest pairing energy in our model(see Ref.Rouhollah2018 ). Because of the high density of states of carriers in this band, d-wave superconductivity is more robust against disorder than s-wave.

VI.2.2 Dominant Bands: Possibility of Two Gap superconductivity

Superconductivity could be established in the odd or even sector of intercalated AA-stacked or simultaneously in both sectors. Even/odd sectors are coupled just via chemical potential. Two nearly “flat bands” with d-wave symmetry Bloch character, crossing the Fermi energy, each related to the graphene-like structures, are responsible for two distinct d-wave superconductivity gaps that could be emerged. The distorted s-wave superconductivity is constructed between quasi-particles in upper bands of both sector that have s-wave symmetry character. At moderate doping, distorted d-wave superconductivity is dominant in both sectors while in high doping, distorted s-wave becomes preferable. Superconductivity with different phases in each of these sectors e.g. distorted s-wave in one sector and d-wave phase symmetry in the other, is not so possible.

To know whether superconductivity in IBlG is single gap or multi gap, depending on what type of intercalant is used, one can inspect numerically which sector will prevail. Hybridization of Carbon ($CP_{z}$) and intercalant ($Is$) orbitals, electron or hole doping factor (chemical potential), nearest neighbor hopping symmetry breaking ($\Gamma$ gap opening) and interlayer coupling ($H_{12}$) are important factors in specifying superconductivity pairing symmetry and dominant bands. While just the even sector is under the influence of C-I orbital hybridization and odd sector bands does not affected directly (except an small gap opening), but interlayer coupling has a dual effects. It reduces the bandwidth of one sector (e.g. the even sector) but simultaneously it increases the bandwidth of the other sector(e.g. odd sector). It therefore plays a crucial role in electronic correlation effects.

Mathematical analysis shows that in the limit of strong interlayer hopping, so that inter- and intra- layer hoppings are the same (up to a minus sign), then bandwidths of the odd (even) sector become completely flat while the bandwidths of the other sector are doubled. In this case superconductivity established in the flat band sector, even with a small Cooper pairing potential $g_{0}$, leads to a high critical temperature i.e. $k_{B}T_{c}=g_{0}$. In this limit all of the possible superconducting phases, i.e. $d,p$ and $s$-wave symmetries, are degenerate. This observation suggests there may be some aspect of unconventional superconductivity in bilayer AA- graphene related to inter- versus intralayer hopping effects be available under high pressure.

The best conditions to induced superconductivity in IBlG are those that interlayer coupling be strong and band structure of even or odd sector or both slightly be deviated from pristine graphen-like structures. Under these circumstances if electron or hole doping causes the nearly flat bands of odd or even sector meet the Fermi surface then chiral $d+id$ odd or even superconductivity may induce. An inadequate inspecting of the band structure of different metal intercalated $C_{6}MC_{6}$ (M=Li, Na, K, Rb, Cs, Gr, Be, Mg, Ca,Sr and Ba) such as that can be seen in the Fig. 3(a) of Ref. Tomoaki2017 , shows that odd sector flat band (not affected by I-C orbital coupling) always lay on top of even sector flat band, that means the interlayer coupling $t_{12}$ in all of them is positive. The structures in which interlayer (IL) band is empty have negligible I-C coupling and so both even and odd sector have the structure similar to shrunken graphene $C_{6}$. The even band of Li intercalated BlG meet the Fermi Surface and interlayer coupling is stronger than the others. These is an interesting structure that could host even-superconducting beside normal odd-Dirac quasi-particles. But since IL band is empty, symmetry not allow out of plane phonon vibration to trigger superconductivity. However it may exhibit richer correlation effect than single layer graphene under pressure, gating or proximity effects. K and Rb intercalated BlG have potential to exhibit superconductivity, while in Gr intercalated BlG it seems that simultaneous odd and even Dirac cones coexist.

Experimental evidence for superconductivity in Ca- intercalated BlG has been reported. Motivated by this observation we performed numerical calculation for $C_{6}CaC_{6}$. Numerical calculations show that for both even and odd sector gap equations, $d$-wave phases, i.e. $\Psi_{{\pm},{d_{xy}}}^{2}$ and $\Psi_{{\pm}{d_{x^{2}-y^{2}}}}^{2}$, are dominant (smaller $g_{0}$ means less interaction energy for pairing) and slightly distorted by intercalation i.e. $\alpha_{d}^{2,\pm}\approx 1$ while $s$-wave symmetry i.e. $\Psi_{{\pm},s}^{2}$ require greater energy and are significantly distorted, $\alpha_{s}^{2,\pm}\neq 1$. Phase Transition from d-wave single-gap to d-wave dual-gap Superconductivity in calcium intercalated bilayer graphene is possible. Although $T_{c}$ experimentally around 4K and theoretically calculated near 6K are reported, also using Raman spectroscopy, possibility of distinguishing intralayer and interlayer electron-phonon interactions in samples of twisted bilayer graphene has been reported by ref.[Eliel2018 ] relying on these results, from Fig. 5 it can be seen that both even and odd $d$-wave phases are nearly degenerate at 2K, consistent with this system being a two gap superconductor around $T_{c}\approx 2K$. Our results support two d-wave gap superconductivity that has been proposed in ref.[Margine2016, ] (Fig. 2), although different sectors were not separated in their studies.

Relying on pre-mentioned properties, AA-stacked bilayer graphene may exhibit feature-rich electronic properties than singlelayer graphen. It seems that study of superconductivity in pristine and intercalated AA-stacking BlG could tend to interesting experimental achievements- such as (even and odd) chiral superconducting $d+id$ pairing which has been predicted primarily in pristine single layer graphene at van Hove singularity point, and also simultaneous coexistence of different phases e.g. superconductivity and normal Dirac quasiparticles or two gap superconductivity (with different chirality) in different branches of their band-structures.

VII Acknowledgments

R. Gholami acknowledges support that allowed an extended visit to the University of California Davis during part of this work. W.E.P. was supported by NSF grant DMR-1207622.

VIII References

IX Author contributions

R.G. and R.M. conceived of the presented idea. R.G. developed the theory and performed the computations. R.M. and W. E. P verified the analytical methods. R.M. supervised the project. R.G, R. M. and W. E. P. discussed the results and commented on the manuscript. S. M performed DFT calculation and designed the figures.

X Additional information

Competing interests: The authors declare no competing interests.

XI Appendix A: Accurate Tight Binding Model

In our previous work we used realistic multiband tight binding model for decorated monolayer graphene and obtained its band structure analytically.Rouhollah2018 Here we follow and generalize that method and find analytic solutions for the intercalated bilayer graphene spectrum in general form. We consider Bloch ket state of Eq.32 as

[TABLE]

in which ${\vec{r}}_{n\alpha}={\vec{r}}_{n}+{\vec{d}}_{\alpha}$ and ${\vec{r}}_{n}$ is $n$th Bravais lattice site vector position and ${\vec{d}}_{\alpha}$ is vector position of the $\alpha$-th subsite with respect to unit cell $n$. The Ca sublattice is labeled by $\alpha=0$ also $A_{1}^{1}$, $A_{2}^{1}$, $A_{3}^{1}$, $B_{1}^{1}$, $B_{2}^{1}$, $B_{3}^{1}$, $A_{1}^{2}$, $A_{2}^{2}$, $A_{3}^{2}$, $B_{1}^{2}$, $B_{2}^{2}$, $B_{3}^{2}$ subsites are labeled by $\alpha=1,...,12$ respectively. $\left|\phi_{n\alpha}\right\rangle=\left|\phi_{n\alpha}({\vec{r}}-{\vec{r}}_{n}-{\vec{d}}_{\alpha})\right\rangle$ is the atomic $\pi$ electron ket state of subsite $\alpha$ of site $n$. The Schrödinger equation for this system is

[TABLE]

Symmetries of this system imply that $\mathscr{C}_{\alpha}(\vec{k})=\pm\mathscr{C}_{\alpha+6}(\vec{k}).$ The Schrödinger equation Eq.89 can be written in the following $13\times 13$ matrix form eigenvalue problem

[TABLE]

where the column matrix $\bm{\mathcal{C}}({\vec{k}})$ is $\bm{\mathcal{C}}({\vec{k}})=({\mathscr{C}}_{1}({\vec{k}})\;{\mathscr{C}}_{2}({\vec{k}})\;...\;{\mathscr{C}}_{6}({\vec{k}}))^{T}$ and the dispersion matrices satisfy $H_{11}=H_{22}$, $H_{12}=H_{21}$ and $h_{01}=h_{02}=h_{10}^{\dagger}=h_{20}^{\dagger}.$ The Ca-C dispersion row matrices $h_{01}(\vec{k})=h_{02}(\vec{k})=(h_{CaA}(\vec{k})~{}~{}h_{CaB}(\vec{k}))$ are given by

[TABLE]

Here $t^{CaC}_{i}$ is the hopping amplitude from Ca to $i$th neighbor C atom. The Ca-Ca dispersion is $h_{0}(\vec{k})={{\varepsilon_{CaCa}}(\vec{k})+{\epsilon_{Ca}}-{\mu_{o}}}$ where

[TABLE]

The interlayer dispersion matrices $H_{11}$, $H_{22}$ and interlayer dispersion matrices $H_{12}$ and $H_{21}$ are given by

[TABLE]

where the off-diagonal carbon-carbon dispersion matrices are

[TABLE]

[TABLE]

in which $m$ and $n$ are layer index. Shorthand notation has been introduced as follows:

[TABLE]

and it has been supposed that $w_{t}=\frac{t^{{}^{\prime}mn}_{1}}{t_{1}^{mn}}=\frac{t^{{}^{\prime}mn}_{2}}{t_{2}^{mn}}=...$ and $\vec{\xi_{i}}=\vec{\tau_{i}}+2\vec{\delta_{i}}$. The $\tau$ and $d$-functions are given by

[TABLE]

Using the following unitary transformation one can separate the bilayer graphene Hamiltonian Eq. 99 into two decoupled single layer pseudo-graphene Hamiltonians, where one of them is decorated with the intercalant layer

[TABLE]

Here $H^{\pm}=H_{11}(\vec{k})\pm H_{12}$ in matrix notation is

[TABLE]

in which $\vec{k}$-dependent on-site energies have been defined as $\varepsilon_{1}^{\pm}({\vec{k}})=\epsilon_{A}-\mu_{0}+\alpha^{\pm}(\vec{k})$ and $\varepsilon_{2}^{\pm}({\vec{k}})=\epsilon_{B}-\mu_{0}+\alpha^{\pm}(\vec{k}).$ Also the following shorthand notation has been introduced

[TABLE]

Unitary transformation of Eq. 136 divides thirteen bands of intercalated bilayer graphene into, six and seven bands groups. Following the approachRouhollah2018 that has been applied to monolayer decorated graphene, an exact analytical solution of the six-band group can be found in general case. These bands are eigenvalues of $H^{-}$ matrix and are not affected directly by the intercalant band.

In the special case of pristine bilayer graphene in which $\gamma^{\pm*}(\vec{k})=\theta^{\pm}(\vec{k})=\beta^{\pm}(\vec{k})$, $\varepsilon_{1}^{\pm}(\vec{k})=\varepsilon_{2}^{\pm}(\vec{k})$ and $\tau_{i}^{\pm}(\vec{k})=d_{i}^{\pm}(\vec{k})$, Eq. 143 easily can be diagonalized to find eigenvalues and also eigenvectors. The eigenvalues are given by

[TABLE]

[TABLE]

[TABLE]

with eigenvectors are given by replacing $m$ in the following equation with $m=0,1,2$ respectively

[TABLE]

wherein

[TABLE]

However, except at $\Gamma$ point it is challenging (and unhelpful) to obtain an exact forms of the seven-bands group analytically. These bands are eigenvalues of the $H^{+}_{c}$ matrix, and analytic expressions for them can be obtained just in the particular case of no hopping between intercalant layer and graphene sheet, similar to the case for lithium intercalated bilayer graphene where intercalant band is empty (no Li-C hopping). In these cases nontrivial solutions are eigenvalues of $H^{+}$ matrix. Eigenvalues of $H^{-}$ and $H^{+}$ matrices are given by Rouhollah2018

[TABLE]

in which $\vec{k}$ dependent chemical potentials are defined as,

[TABLE]

The $\Pi_{0}^{\pm}(t_{2},\vec{\xi_{i}},\vec{k})$ function is introduced as

[TABLE]

and, ${c_{0}^{\pm}}(t_{2},\vec{\xi_{i}},\vec{k})={\gamma^{\pm*}}(\vec{k})[{\beta^{\pm}}(\vec{k}){\theta^{\pm}}(\vec{k})]+{\gamma^{\pm}}(\vec{k})[{\beta^{\pm}}(\vec{k}){\theta^{\pm}}(\vec{k})]^{*}$ and ${c_{1}^{\pm}}(t_{2},\vec{\xi_{i}},\vec{k})={\left|\beta^{\pm}(\vec{k})\right|^{2}}+{\left|\theta^{\pm}(\vec{k})\right|^{2}}+{\left|\gamma^{\pm}(\vec{k})\right|^{2}}$. Also, $w_{m}^{\pm}(t_{i},\vec{\xi_{i}},\vec{k})$ are eigenvalues of the following matrices

[TABLE]

The eigenvalues can be obtained as

[TABLE]

wherein

[TABLE]

[TABLE]

One can write the $H_{D}$ matrix in $H^{\pm}$ bases i. e. $H^{{}^{\prime}}_{D}=U^{\dagger}H_{D}U$,

[TABLE]

The upper left portion of Eq.176 can be obtained sufficiently well by perturbation theory.

XII Appendix B: Bogoliubov-de Gennes Transformation

The interacting Hamiltonian $H_{su}$ in matrix representation is $14\times 14$ matrix,

[TABLE]

where $\hat{\Psi}^{\dagger}(\vec{k})=(\hat{c}^{\dagger}_{0\uparrow}(\vec{k})\hat{c}^{\dagger}_{1\uparrow}(\vec{k})\;\hat{c}^{\dagger}_{2\uparrow}(\vec{k})\;...;\hat{c}^{\dagger}_{12\uparrow}(\vec{k})\;\hat{c}_{0\downarrow}(-\vec{k})\hat{c}_{1\downarrow}(-\vec{k})\;\hat{c}_{2\downarrow}(-\vec{k})\;...\;\hat{c}_{12\downarrow}(-\vec{k}))$. $H_{N}$ is Hamiltonian of the normal state and $H_{p}$ is the pair interaction matrix. The full matrix must be diagonalized to obtain the quasiparticle spectrum.

The mean field superconducting Hamiltonian of Eq. 1 in Nambu space is ${{\hat{H}}_{su}}=\sum_{\vec{k}}{{{\hat{\Psi}}^{\dagger}}(\vec{k})}{H_{su}}(\vec{k})\hat{\Psi}(\vec{k})$ where $H_{su}$ in matrix representation is,

[TABLE]

and $\hat{\Psi}^{\dagger}(\vec{k})=(\hat{c}^{\dagger}_{0\uparrow}(\vec{k})\hat{c}^{\dagger}_{1\uparrow}(\vec{k})\;\hat{c}^{\dagger}_{2\uparrow}(\vec{k})\;...;\hat{c}^{\dagger}_{12\uparrow}(\vec{k})\;\hat{c}_{0\downarrow}(-\vec{k})\hat{c}_{1\downarrow}(-\vec{k})\;\hat{c}_{2\downarrow}(-\vec{k})\;...\;\hat{c}_{12\downarrow}(-\vec{k}))$ where $H_{11}(\vec{k})=H_{22}(\vec{k})$. The interlayer pairing matrices are ${H_{11}^{P}(\vec{k})}={H_{22}^{P}(\vec{k})}$ and interlayer pairing matrices are ${H_{12}^{P}(\vec{k})}={H_{21}^{P}(\vec{k})}$. The pairing matrices are given by

[TABLE]

where $m$ and $n$ are layer index which can take 1 or 2. The order parameters accordingly in Fourier space are

[TABLE]

where $<ij>$ subscript indicate nearest neighbor pairing amplitude in real space as illustrated in Fig.1. Introducing the following unitary transformation matrix,

[TABLE]

one can transform Eq.210. Eq.196 can be transformed to the block diagonalize form

[TABLE]

in which

[TABLE]

and, $\Lambda^{\dagger}(\vec{k})=\left([\hat{c}^{\dagger}_{0\uparrow}(\vec{k})\hat{c}^{+\dagger}_{1\uparrow}(\vec{k})...\hat{c}^{+\dagger}_{6\uparrow}(\vec{k})\;\hat{c}_{0\downarrow}(-\vec{k})\hat{c}_{1\downarrow}^{+}(-\vec{k})...\hat{c}_{6\downarrow}^{+}(-\vec{k})]~{}~{}~{}[\hat{c}^{-\dagger}_{1\uparrow}(\vec{k})\;\hat{c}^{-\dagger}_{2\uparrow}(\vec{k})...\hat{c}^{-\dagger}_{6\uparrow}(\vec{k})\;\hat{c}_{1\downarrow}^{-}(-\vec{k})\;\hat{c}_{2\downarrow}^{-}(-\vec{k})...\hat{c}_{6\downarrow}^{-}(-\vec{k})]\right),$ where $\hat{c}^{(\pm)\dagger}_{m\sigma}(\vec{k})=\frac{1}{\sqrt{2}}(\hat{c}^{\dagger}_{m\sigma}(\vec{k})\pm\hat{c}^{\dagger}_{m+6,\sigma}(\vec{k})).$ New pairing matrices $H_{p}^{+}(\vec{k})$ and $H_{p}^{-}(\vec{k})$ are defined as,

[TABLE]

From Eq. 213 it can be seen that the superconducting Hamiltonian $H_{su}$ can be diagonalized into two decoupled new superconducting Hamiltonian $H^{+}_{su}$ and $H^{-}_{su}$. Thus electrons just can be paired within the seven bands sector $H_{c}^{+}$ or within the six bands sector $H^{-}$, without coupling between the sectors. Thus superconductivity in bilayer graphene cane be interpreted as two decoupled monolayer graphene-like systems with independent behaviors.

XIII Appendix C: Two Superconducting Gap Equations

The linearized superconducting gap equation are obtained by minimizing the quasiparticle free energy with respect to the nearest neighbor order parameter, or equivalently with respect to $\Delta_{\pm}^{\alpha}$. The free energy of system is

[TABLE]

For $F^{+}$ the summation runs over $n=1,...,7$ giving $E_{n}^{Q}=E_{n,s}^{Q+}$; for $F^{-}$ the summation takes $n=8,...,13$ values giving $E_{n}^{Q}=E_{n,s}^{Q-}$, with $E_{n,s}^{Q\pm}$ introduced in Eqs. 53 and 54.

Minimization of the free energy with respect to $\Delta_{+}^{\alpha}$ gives

[TABLE]

giving independent gap equations for the seven bands odd-symmetry graphene-like Hamiltonian $H_{c}^{+}$. Minimizing the free energy with respect to $\Delta_{-}^{\alpha}$ gives

[TABLE]

where $\Delta^{\alpha}$ as illustrated in Fig.1(b) covers all possible nearest neighbor inter- and intra-layer C-C pairing amplitudes. Eqs. 258 and Eq. 224 in matrix form written as

[TABLE]

Equivalently, Eqs. 258 and 224 can be combined in the following non-Hermitian eigenvalue problem,

[TABLE]

in which $\kappa_{0}=\frac{g_{0}}{g_{0}^{{}^{\prime}}}$ and

[TABLE]

$\kappa=\frac{g_{1}}{g_{0}}=\frac{g_{1}^{{}^{\prime}}}{g_{0}^{{}^{\prime}}}$, $V_{1}^{T}=(\Sigma_{1}~{}\Sigma_{2}~{}\Sigma_{3})$, $V_{2}^{T}=(\Pi_{1}~{}\Pi_{2}~{}\Pi_{3})$, and $V_{3}^{T}=(\Delta_{1}~{}\Delta_{2}~{}\Delta_{3})$. Also $V_{1}^{{}^{\prime}T}=(\Sigma_{1}^{{}^{\prime}}~{}\Sigma_{2}^{{}^{\prime}}~{}\Sigma_{3}^{{}^{\prime}})$, $V_{2}^{{}^{\prime}T}=(\Pi_{1}^{{}^{\prime}}~{}\Pi_{2}^{{}^{\prime}}~{}\Pi_{3}^{{}^{\prime}})$, and $V_{3}^{{}^{\prime}T}=(\Delta_{1}^{{}^{\prime}}~{}\Delta_{2}^{{}^{\prime}}~{}\Delta_{3}^{{}^{\prime}})$. Equation 241 is in fact the matrix representation of gap equation resulting from minimization of the free energy with respect to nearest neighbor order parameters instead of $\Delta_{\pm}^{\alpha}$, which can be solved to obtain the differing superconductivity phases and pairing interaction potentials.

In the limiting case $\kappa_{0}\rightarrow 1$, the inter- and intra-layer pairing amplitudes in real space are equal, i.e. $g_{0}=g_{0}^{{}^{\prime}}$ and $g_{1}=g_{1}^{{}^{\prime}}$. This restriction makes the matrix gap equation hermitian and implies that band order parameters $\Delta_{mn}^{\pm}(\vec{k})$ can be interpreted physically as pairing of electrons in different bands with pairing interaction $g_{0}^{\pm}$. Just in this limit $\Delta_{mn}^{\pm}(\vec{k})$ is equal to the product of band Green function and $g_{0}$,

[TABLE]

Here $\hat{d}^{\pm\sigma}_{i}(\vec{k})=\sum^{7}_{m=1}{\mathscr{C}}^{\pm*}_{m}(E_{i}(\vec{k}))\hat{c}^{\sigma}_{m}(\vec{k})$ annihilates an electron with spin $\sigma$ in the $i$th six (odd) or seven (even) sector bands with energy $E_{i}^{\pm}(\vec{k})$. In this limit, Eq. 241 has two solutions $\Psi_{\Delta}=\Psi_{\Delta}^{{}^{\prime}}$ or $\Psi_{\Delta}=-\Psi_{\Delta}^{{}^{\prime}}$. The corresponding gap equations Eqs. 234 and 241 become decoupled gap equations corresponding to the even or odd sector of the graphene-like systems.

[TABLE]

XIV Appendix D: flat band(s) Superconductivity

Mirror symmetry transformation rearranges the noninteracting Hamiltonian Eq. 32 as the direct sum of two single layer pseudo-graphene structures $\hat{H}_{N}=\hat{H}^{+}_{N}\oplus\hat{H}^{-}_{N}$ (even sector (+ sign) and odd sector (- sign))

[TABLE]

with renormalized hopping integrals of the form $t^{\pm}_{i\alpha\sigma,j\beta\sigma}=t^{inter}_{i\alpha\sigma,j\beta\sigma}\pm t^{intra}_{i\alpha\sigma,j\beta+6\sigma}$. In the limit case of strong interlayer hopping wherein $t^{inter}_{i\alpha\sigma,j\beta\sigma}\to\pm t^{intra}_{i\alpha\sigma,j\beta+6\sigma}$ one can see that odd (or even) sector bandwidths completely become flat while the other sector bandwidth doubles. In this limit, thermal weight factor of Eq. 80 $\frac{\tanh(\frac{{\beta E_{i}^{\pm}}}{2})}{E_{j}^{\pm}(\vec{k})+E_{i}^{\pm}(\vec{k})}\to\frac{\beta}{4}$ and so $\Gamma$ matrix elements are given by

[TABLE]

These elements are linked to the normal state Bloch coefficients via the $\Omega_{ij}^{\pm\beta}(\vec{k})$ factors that given by Eq. 55. For the case of pristine bilayer graphene, Eq. 258 can be determined analytically. In this limit, normal state Bloch coefficients are given by

[TABLE]

wherein $e^{i\phi_{m}^{\pm}(\vec{k})}=\frac{\eta_{m}^{*\pm}}{|\eta_{m}^{\pm}|}$ and $\eta_{m}^{\pm}(\vec{k})=d_{2}^{\pm}(\vec{k})+u_{m}d_{1}^{\pm}(\vec{k})+u_{m}^{*}d_{3}^{\pm}(\vec{k})$. $\Omega_{ij}^{\pm\beta}(\vec{k})$ factors can be calculated by substituting Bloch coefficients Eq. 259 in the Eq. 55. For instant one can show

[TABLE]

By calculating a large number of these factors and replacing them in the Eq. 258 one can obtain

[TABLE]