Symmetry of superconducting correlations in displaced bilayers of graphene
Mohammad Alidoust Morten Willatzen Antti-Pekka Jauho





Department of Physics, K.N. Toosi University of Technology, Tehran 15875-4416, Iran

Beijing Institute of Nanoenergy and Nanosystems, Chinese Academy of Sciences,
No. 30 Xueyuan Road, Haidian District, Beijing 100083, China
Department of Photonics Engineering, Technical University of Denmark, Kongens Lyngby, DK-2800, Denmark

Center for Nanostructured Graphene (CNG), DTU Nanotech, Technical University of Denmark, DK-2800 Kongens Lyngby, Denmark



Using a Green's function approach, we study phonon-mediated superconducting pairing symmetries that may arise in bilayer graphene where the monolayers are displaced in-plane with respect to each other. We consider a generic coupling potential between the displaced graphene monolayers, which is applicable to both shifted and commensurate twisted graphene layers; study intralayer and interlayer phonon-mediated BCS pairings; and investigate AA and AB(AC) stacking orders. Our findings demonstrate that at the charge neutrality point, the dominant pairings in both AA and AB stackings with intralayer and interlayer electron-electron couplings can have even-parity s-wave class and odd-parity p-wave class of symmetries with the possibility of invoking equal-pseudospin and odd-frequency pair correlations. At a finite doping, however, the AB (and equivalently AC) stacking can develop pseudospin-singlet and pseudospin-triplet d-wave symmetry, in addition to s-wave, p-wave, f-wave, and their combinations, while the AA stacking order, similar to the undoped case, is unable to host the d-wave symmetry. When we introduce a generic coupling potential, applicable to commensurate twisted and shifted bilayers of graphene, d-wave symmetry can also appear at the charge neutrality point. Inspired by a recent experiment where two phonon modes were observed in a twisted bilayer graphene, we also discuss the possibility of the existence of two-gap superconductivity, where the intralayer and interlayer phonon-mediated BCS picture is responsible for superconductivity. These analyses may provide a useful tool in determining the superconducting pairing symmetries and mechanism in bilayer graphene systems.




Symmetry of superconducting correlations in displaced bilayers of graphene
    Antti-Pekka Jauho
    December 30, 2023
==========================================================================





§ INTRODUCTION



Graphene, a two-dimensional honeycomb arrangement of covalently bound carbon atoms, can be used as a basic building block for novel multilayer systems coupled by interlayer van der Waals forces. These weak interlayer bonding forces allow the possibility to create different ordering stacks that possess dramatically different low-energy properties, offering low-cost tools for rich band structure and material engineering perspectives. The ordering of the stacks can be controlled externally by growth techniques or by the application of tensile stress to laterally shift and/or rotate the layers in plane with respect to each other <cit.>.

Depending on the stacking order, the multilayers can have a number of different configurations. Bilayer graphene (BLG) where the top and bottom layers are mutually rotated is an example of introducing additional features by a simple geometric action. It has been shown that at low angles of the twist, the BLG develops a highly nontrivial band structure <cit.>. An efficient theory tool to approach such a problem and describe its underlying physics is the effective Hamiltonian method. This technique assumes that the two pristine graphene layers are unaffected under the displacement and all changes can be translated into a new coupling potential between the two layers <cit.>.

Physically, one can expect that the charge carriers in BLG experience different interactions due to different phonon modes generated within or between the layers. The two pristine graphenes have identical phonon modes (due to the intralayer covalent bonds) while a different mode appears due to the interlayer coupling (van der Waals bonds). Note that our conclusions depend on having two different coupling strengths (intralayer and interlayer couplings), and the microscopic details of the coupling do not matter.  For example, the interlayer coupling could have a significant covalent component. <cit.>. 
In particular, recent experiments demonstrated the different intralayer and interlayer electron-phonon interactions by using a Raman spectroscopy technique: two resonance peaks in the Raman modes  were observed by tuning the energy of excitation with IR and UV photons <cit.>.


The electron-phonon coupling affects, for instance, the electronic mobility, thermal conductivity, and the (possible) superconducting phases. A recent landmark experiment showed intrinsic superconductivity with a critical temperature of ∼ 1.7 K in a twisted BLG, away from the charge neutrality point μ>0 at a very low carrier density of 10^11 cm^-2. Interestingly, the general features of the phase diagram seemed to share the same phenomenology as high T_c superconductors <cit.>. Soon after, another experiment <cit.> (and recently few others <cit.>) and a very large number of theoretical papers have appeared, many of which focus  on the phonon-mediated single-gap superconductivity, and addressed various aspects of the observed superconducting phase (see, e.g. Refs. [teor_tbg1,teor_tbg2,teor_tbg3,teor_tbg4,teor_tbg5,teor_tbg7,
teor_tbg9,teor_tbg11,teor_tbg12,teor_tbg13,teor_tbg14,teor_tbg15,teor_tbg16,teor_tbg17,teor_tbg18,teor_tbg19,teor_tbg20,teor_tbg21,YWChoi]). However, earlier theories have not discussed the role of deformations in determining the superconducting properties of BLG <cit.>. 

In this paper, we construct the Green's function for various coupling scenarios and present an extensive study of phonon-mediated pairing possibilities in BLG. As mentioned above, a recent experiment has demonstrated two phonon modes in a twisted BLG <cit.>. One mode was attributed to the intralayer vibrations while the other mode was related to interlayer couplings. Therefore, in the same spirit, one can expect that BLG can potentially develop two-gap superconductivity in the phonon-mediated electron-electron coupling picture, where each gap is mainly supported by one of the two phonon modes. Motivated by this picture, we also develop a two-gap superconductivity model for BLG. We employ the effective Hamiltonian method, consider a generic coupling matrix for the interlayer coupling potential of BLG, and study the different superconducting pairings and correlations. The generic expressions we derive can be applied to both twisted and shifted BLG when using the appropriate coupling terms, as will be discussed below.

For concreteness, we study shifted BLG with specific shift directions so that AA stacking turns into AB, AB into AC, and AC into AA in the extreme limits of the shift. We find that for an undoped system, the anomalous Green's function possesses only even-parity s-wave class and odd-parity p-wave class of symmetries for both AA and AB orderings with intralayer/interlayer phonon-mediated electron-electron coupling. For a doped system, however, AB ordering offers a rich variety of pairing symmetries such as s-wave, p-wave, d-wave, f-wave, and symmetry combinations. On the other hand, AA retains its previous even-parity s-wave, odd-parity p-wave (now f-wave class) symmetries. Introducing a generic coupling potential that accommodates in-plane twist and shift between layers, a d-wave symmetry can appear for both AA and AB orderings at μ=0. Therefore, our findings suggest that by shifting the layers of a BLG with respect to each other, an effective switching to d-wave superconducting pairing with specific pseudospin (Pspin) states is accessible. This phenomenon offers a platform to control both superconducting critical temperature and the Pspin degree of freedom of superconducting correlations by experimentally simple actions that can find crucial applications in spintronics in addition to great interest to fundamental researches. The pairing symmetries predicted here can be probed experimentally by point contact tunneling spectroscopy experiments so that the conductance in different directions will be different and reveals these symmetries. It is worth mentioning that similar theory study as the present paper was conducted for black phosphorus monolayer in Ref. BP2019. It was found that a significant transition from effective s-wave (p-wave) to d-wave (f-wave) symmetry class is accessible simply by the exertion of strain into the plane of black phosphorus monolayer.

Using the obtained anomalous Green's function, we also calculate the intralayer and interlayer superconducting gaps as a function of temperature. For a given s-wave electron-electron interaction potential, our results show that the intralayer coupling in AA stacking order has both the largest gap amplitude and critical temperature, while the next largest gap and critical temperature belongs to interlayer coupling in AB stacking order. Also, the shift between the layers has opposite effects on the gap and critical temperature of intralayer and interlayer superconductivity in AA and AB orderings: the displacement applied to AA stacking with intralayer coupling or AB stacking with interlayer coupling reduces both the gap and critical temperature, while it enhances both of them for AA stacking with interlayer coupling and for AB stacking with intralayer coupling. 

We note that our analysis include the influence of chemical potential (away from the charge neutrality point) and are independent of the amplitude of the phonon-mediated electron-electron interaction: Δ. We emphasize that our paper discusses the possible symmetries of superconducting pair correlations away from the magic angle of an incommensurately twisted BLG.

In Sec. <ref>, we first outline the Green's function technique, the effective Hamiltonian approach, and discuss different pairing scenarios that may arise in BLG (some related discussions are presented in Appendix <ref>). In Sec. <ref>, we study the anomalous Green's function both analytically and numerically for different scenarios in AA and AB stacking orders. We first consider an undoped system and derive analytic expressions to the components of anomalous Green's function in Sec. <ref>. Next, for a doped system, we evaluate the Green's functions numerically, and support these results by analytic expressions given in Appendix <ref>. In Sec. <ref>, we consider a generic coupling potential between the two pristine graphene layers and repeat our studies of Sec. <ref>, now in the presence of a small displacement between the two layers. We also study the temperature dependencies of the superconducting gaps. We finally give concluding remarks in Sec. <ref>.




§ LOW ENERGY EFFECTIVE HAMILTONIAN AND GREEN'S FUNCTION APPROACH


We follow the effective Hamiltonian strategy <cit.> where the Hamiltonian governing the low-energy physics of a generic BLG involving an arbitrary in-plane displacement between the layers can be expressed as

    H   =   ∫dk/(2π)^2ψ̂^†(k)H(k)ψ̂(k)
       =   ∫dk/(2π)^2ψ̂^†(k){  H_1( k)ρ_1+H_2( k)ρ_2 + T̃(k)ρ_+ + T̃^†(k)ρ_-   }ψ̂(k),

in which the layer degree of freedom is described in terms of Pauli matrices ρ_0,z,x,y, and we have defined 2ρ_1=ρ_0+ρ_z, 2ρ_2=ρ_0-ρ_z, 2ρ_+=ρ_x+iρ_y, 2ρ_-=ρ_x-iρ_y.  The matrix T̃(k) couples the top (1) and bottom (2) pristine single-layer graphene (SLG) and contains all information about the relative displacement between the layers (see below) <cit.>. The Hamiltonians of the top and bottom layers are denoted by H_1,2( k)=ħ v_F k·σ, respectively, in which σ=(σ_x,σ_y) and σ_x, σ_y are Pauli matrices. The multiplication of Pauli matrices ρ_i and σ_i implies a tensor product so that the Hamiltonian (<ref>) is a 4× 4 matrix and the matrix T̃(k) has an implicit σ (refer to Eqs. (<ref>) and (<ref>) for further details). Here, the two-dimensional momentum k is in the plane of SLG. The field operator associated with the Hamiltonian is given by ψ̂^†(k)=(ψ_1^†,ψ_1^†,ψ_2^†,ψ_2^†) where the first index is 1 (2) for the top (bottom) SLG and the second index is ↑ (↓) for Pspin up -sublattice A- (Pspin down -sublattice B-). Note that, according to the mean-field definition of BCS superconductivity, two electrons with opposite spins and momenta are coupled through an attractive potential. This means that electron and hole excitations should be taken from opposite corners of the Brillouin zone <cit.>. In this case, one can show that the system possesses spin degeneracy <cit.>.  Nonetheless, to denote the Pspin and keep our notation simplified simultaneously, we have dropped the spin and valley indices. A detailed discussion is presented in Appendix <ref>.

We next define the parameters that characterize a BLG system with a relative displacement between the top and bottom SLG. The monolayer pristine graphene can be divided into two sublattices A and B that can be described by two sublattice positions in real space ν_A = (0,0), ν_B =2( a_1+ a_2)/3 and two lattice vectors a_1 = a(1,0),  a_2 = a(1,√(3))/2. Using this notation, the different stackings can be described by the sublattice positions given in Table <ref>. Figure <ref>(a) displays graphene lattice in real space with sublattices A and B marked by dark red and dark blue circles, respectively. Also, we have illustrated AB ordering in Fig. <ref>(b) where only the sublattice B of top monolayer graphene is aligned with the sublattice A of bottom layer. In AA case, A and B sublattices of the top layer are aligned with A and B sublattices of the bottom layer (see Table <ref>).

The reciprocal lattice of a BLG can be described by the lattice vectors b_1 = 2a^-1π(1,-1/√(3)),  b_2 = 2a^-1π(0,2/√(3)) in reciprocal space as shown in Fig. <ref>(c). Since we are interested in the low-energy physics, the Hamiltonian can be expanded around reciprocal sublattice points:

    K_j=K_0+G_j,

in which K_0=(2 b_1+ b_2)/3, G_0= 0, G_1=- b_1, G_2=- b_1- b_2 (see Fig. <ref>(b)).
The periodic arrangement of atoms in SLG is spanned by 𝐫=R_j+ν_α, where R_j are integer combinations of the primitive lattice vectors a_j, and the orbitals centered at positions ν_α. In the following, 𝐫 is the location at which the coupling between the top and bottom layers takes place. For example, the coupling of A and B sublattices shown in Fig. <ref>(b). The displacement enters the coupling matrix T̃ via a phase factor <cit.>,


    T̃(k)=∑_j=0,1,2M_j^XXt_⊥(K_j+k)/3e^i(K_j+k)(u_2-u_1),

in which t_⊥(K_j+k) is the interlayer hopping amplitude, given by t_⊥( q)=V_u^-1∫ d r t_⊥( r)e^i q r, t_⊥( r)=∑_i,jt_2,j^1,i⟨ r,2|c^†_2,jc_1,i| r',1⟩ with c^†,c being the creation and annihilation operators and r'= r+δ in which δ denotes the distance between hopping sites. V_u is the volume of unit cell, i,j run over the lattice sites, u_1,2 are the displacements of the top and bottom layers, and M_j^XX defines stacking order:

    M_j^AA=( [           1 e^-i(2jπ/3); e^+i(2jπ/3)           1 ]),   M_j^AB=( [ e^+i(2jπ/3)           1; e^-i(2jπ/3) e^+i(2jπ/3) ]),
       M_j^AC=( [ e^-i(2jπ/3) e^+i(2jπ/3);           1 e^-i(2jπ/3) ]).





In what follows, we consider two types of couplings that may lead to phonon-mediated superconductivity in a BLG <cit.>. The top and bottom SLG are assumed to possess identical phonon modes while the coupling between the two SLG induces a new and different phonon mode <cit.>. Thus, the BCS electron coupling within each layer is called `intralayer phonon-mediated electron-electron coupling' and the electron coupling between the two layers is called `interlayer phonon-mediated electron-electron coupling' throughout this paper. In this basis, the two types of pairing mechanisms are characterized by the following two-electron amplitudes: 


    i)     Δ_S⟨ψ^†_1ψ^†_1⟩+h.c.  and Δ_S⟨ψ^†_2ψ^†_2⟩+H.c., 
        ii)     Δ_B⟨ψ^†_1ψ^†_2⟩+H.c.  ,


where Δ_S and Δ_B are the BCS spin-singlet phonon-mediated electron-electron coupling within each layer and between the layers, respectively. For more details see Appendix <ref>.

We next consider the Green's functions for the generic BLG system in the presence of superconductivity. The normal Green's functions g and the anomalous Green's function f are defined as follows; <cit.>


    g _ρρ'^σσ'(τ,τ';  𝐫,  𝐫') = -⟨ T_τψ_ρσ (τ, 𝐫') ψ_ρ'σ'^†(τ', 𝐫')  ⟩,       
          g_ρρ'^σσ'(τ,τ';  𝐫,  𝐫') = - ⟨ T_τψ^†_ρσ (τ, 𝐫) ψ_ρ'σ'(τ', 𝐫')  ⟩,       
           f_ρρ'^σσ'(τ,τ';  𝐫,  𝐫') = + ⟨ T_τψ_ρσ (τ, 𝐫) ψ_ρ'σ'(τ', 𝐫')  ⟩,       
          f_ρρ'^σσ'†(τ,τ';  𝐫,  𝐫') = + ⟨ T_τψ^†_ρσ (τ, 𝐫) ψ^†_ρ'σ'(τ', 𝐫')  ⟩,


where T_τ is the time ordering operator, and τ, τ' are the imaginary times. Here, σ,σ' and ρ,ρ' denote the Pspin and layer indices and ⟨ ... ⟩ is the thermodynamics averaging.
In particle-hole space, the Green's function satisfies:

    ( [          Ĥ( 𝐫)-iω_n               Δ̂(r);             Δ̂^*(r) T_tĤ( 𝐫) T_t^†+iω_n ])ǧ(iω_n; 𝐫, 𝐫')
    =δ( 𝐫- 𝐫'),

in which ω_n=π (2n+1) k_B T is the Matsubara frequency, n∈ Z, k_B is the Boltzman constant, T is temperature, T_t is the time-reversal operator, and Δ̂( 𝐫) is the superconducting gap. Note that Ĥ( 𝐫) is obtained by replacing k→ -i∇ in H(k), in Eq. (<ref>). 
The matrix form of the Green's function is given by;

    ǧ(iω_n; 𝐫, 𝐫')=(  [    ĝ(iω_n; 𝐫, 𝐫')   f̂(iω_n; 𝐫, 𝐫'); f̂^†(iω_n; 𝐫, 𝐫')    ĝ(iω_n; 𝐫, 𝐫') ]).

We denote 4× 4 matrices by `hat' symbol, □̂, and 8× 8 matrices by `check' symbol, □̌.

We now proceed to derive the anomalous Green's function f for the different coupling scenarios (Eq. (<ref>) below gives the self-consistency condition these functions must satisfy).<cit.>




§ BILAYER GRAPHENE WITHOUT ANY DEFORMATIONS


To begin, we consider spin-singlet intralayer and interlayer phonon-mediated electron-electron couplings and find possible superconducting pair correlations in a BLG with zero doping. We thus assume that the top and bottom SLG host identical phonon-mediated electron-electron interactions, i.e., Δ_S^1≡Δ_S^2. In what follows, we use a compact notation for the different combinations of momenta (see Table <ref>) that arise during the calculations to simplify the resulting expressions. We first consider the simplest cases that result in analytic solutions and gradually add parameters to reach more complicated low-energy effective Hamiltonians.



This section is split into two subsections, where we study the undoped and doped cases, respectively, considering both AA and AB stacking orders. For symmetry reasons, the AC stacking order leads to similar results as AB stacking and, thus, we omit presenting results for the AC stacking order. We also assume that the system is translationally invariant in the xy plane (the plane of BLG) and thus the corresponding Hamiltonian is diagonal in momentum space as a function of k_x and k_y. Consequently, we deal with a set of algebraic equations in momentum space for the Green's function. The anomalous components of Green's function f_ρρ'^σσ' correspond to pairing correlations that may arise in the system whereas the normal Green's function g_ρρ'^σσ' contains information about local density of states. Hence, in what follows, we focus on the anomalous Green's function f_ρρ'^σσ' and do not present explicitly the normal Green's function g_ρρ'^σσ' although we have derived all the components in our actual calculations. 

In the following, inspired by the atomic orbital symbols, we use the following nomenclature: terms proportional to k_±(=k_x± ik_y) → p-wave, k_±^2 → d-wave, and k_±^2k_∓ → f-wave, etc. (see also Tables <ref> and <ref>).

In Secs. <ref> and <ref>, we consider undoped μ=0 and doped μ≠ 0 BLG, respectively. Each of these sections studies both AA and AB orderings with both intralayer Δ_S and interlayer Δ_B phonon-mediated electron-electron interaction. 



 §.§ Undoped bilayer graphene: μ=0


We consider two cases separately: spin singlet intralayer phonon-mediated electron-electron coupling Δ_S and spin-singlet interlayer phonon-mediated electron-electron coupling Δ_B. In both cases, we assume a constant coupling potential between the two monolayers given by: t=t^†. The components of coupling potential are defined by Eq. (<ref>). 



  §.§.§ Intralayer coupling: Δ_S 0, Δ_B=0


∙AA stacking:
In this case, the components of the anomalous Green's function Eqs. (<ref>) are: 






    Ω f_11^= 0,
          Ω f_11^=-Δ_S(k_-k_++t^2+Δ_S^2+ω_n^2), 
          Ω f_12^=-2tΔ_S k_-, 
          Ω f_12^=0, 
          Ω f_11^=+Δ_S(k_-k_++t^2+Δ_S^2+ω_n^2), 
          Ω f_11^=0, 
          Ω f_12^=0, 
          Ω f_12^=+2tΔ_S k_+
          Ω=k_-^2k_+^2+2k_-k_+(Δ_S^2-t^2+ω_n^2)+(t^2+Δ_S^2+ω_n^2)^2.


As seen, the intralayer Pspin-triplet and interlayer Pspin-singlet correlations are not induced while, as expected, the intralayer Pspin-singlet pairings are nonzero and have even-parity even-frequency symmetry. Interestingly, the Pspin-triplet interlayer correlations are nonzero, proportional to the coupling strength between the layer and possess odd-parity even-frequency symmetry.


∙AB stacking: Next, we consider AB stacking order and intralayer opposite-spin phonon-mediated electron-electron interactions. Constructing the anomalous Green's function, we find: 


    Ω f_11^= 0, 
          Ω f_11^=-Δ_S(k_-k_++Δ_S^2+ω_n^2), 
          Ω f_12^=- itΔ_Sω_n, 
          Ω f_12^=0, 
          Ω f_11^=+Δ_S(k_-k_++Δ_S^2+ω_n^2), 
          Ω f_11^=0, 
          Ω f_12^=tΔ_S (k_- - k_+), 
          Ω f_12^=+itΔ_Sω_n, 
          Ω=(k_-k_++Δ_S^2)^2+(2k_-k_++t^2+2Δ_S^2)ω_n^2+ω_n^4.


Similar to the AA case, the intralayer Pspin-triplet correlations are zero while the Pspin-singlet pairings are nonzero with even-parity even-frequency. Due to the specific coupling matrix in the AB stacking order, one of the Pspin-singlet correlation components is odd-parity and even in frequency while the Pspin-triplet ones are even in parity and odd in frequency. Note that these unconventional correlations are proportional to the interlayer coupling strength t.



  §.§.§ Interlayer BCS coupling: Δ_S=0, Δ_B 0


Here, we consider a situation where one electron from the top layer and one from the bottom layer are coupled through the interlayer phonons which is given by Eq. (<ref>). For AA and AB stack orderings of the two monolayer graphene, we find the following anomalous components of the Green's function, Eqs. (<ref>);


∙AA stacking: 


    Ω f_11^= 0, 
          Ω f_11^=-2itΔ_Bω_n, 
          Ω f_12^=0, 
          Ω f_12^=-Δ_B (k_-k_+-t^2+Δ_B^2+ω_n^2), 
          Ω f_11^=2itΔ_Bω_n, 
          Ω f_11^=0, 
          Ω f_12^=Δ_B (k_-k_+-t^2+Δ_B^2+ω_n^2), 
          Ω f_12^=0, 
          Ω=(k_-k_+-t^2+Δ_B^2)^2+2(k_-k_++t^2+Δ_B^2)ω_n^2+ω_n^4.




∙AB stacking: 


    Ω f_11^= 0, 
          Ω f_11^=-tΔ_B k_+, 
          Ω f_12^=0, 
          Ω f_12^=-Δ_B (k_-k_++Δ_B^2+ω_n^2), 
          Ω f_11^=tΔ_B k_+, 
          Ω f_11^=0, 
          Ω f_12^=Δ_B (k_-k_++t^2+Δ_B^2+ω_n^2), 
          Ω f_12^=0, 
          Ω=k_-^2k_+^2+2k_-k_+(Δ_B^2+ω_n^2)+(Δ_B^2+ω_n^2)(Δ_B^2+t^2+ω_n^2).




Comparing with the Green's function components of the same stacking orders but with intralayer phonon-mediated electron-electron coupling discussed in Subsec. <ref>, we readily find that the AA stacking in the current scenario develops even-parity and odd-frequency superconducting correlations similar to the AB stacking order of Sec. <ref>. The same symmetry correspondence can be seen between AB stacking of this current subsection and Sec. <ref>. In that case, BLG develops odd-parity even-frequency correlations. Note that, unlike the intralayer phonon-mediated electron-electron coupling scenario, Sec. <ref>, here all components of the anomalous Green's function are spin-singlet. Also, the odd-parity components are of p-wave type (i.e., proportional to k_±. See Tables <ref> and <ref>). Thus, for the undoped case, the only possible pairing correlations with the intralayer and interlayer phonon-mediated electron-electron couplings are of s-wave and p-wave class regardless of the stacking order.












 §.§ Bilayer graphene with finite doping μ≠ 0


Next, we consider finite doping, and reinvestigate the superconducting correlations given by the anomalous Green's function as discussed in the previous section. 



  §.§.§ AA stacking  order: phonon-mediated intralayer and interlayer electron-electron couplings


In the presence of a finite doping, analytic expressions for the anomalous Green's function become cumbersome. However, they allow one to deduce and analyze the pairing symmetries that may arise. We have presented them in Appendix <ref> and exploit them as a benchmark for our numerical calculations.

To evaluate the effective momentum dependencies, we set the Matsubara frequency fixed at the lowest mode with n=0 for the individual components of the anomalous Green's function evaluated at low temperatures (typically, T∼ 0.01T_c). We also set a `representative' ratio to the chemical potential over coupling strength μ/t=1.5 throughout the numerical study without loss of generality. Figure <ref> illustrates the momentum space map of the real and imaginary parts of the Green's function components in AA stacking order. In Fig. <ref>(a), the phonon-mediated electron-electron coupling is of the intralayer type Δ_S≠ 0, Δ_B=0 while in Fig. <ref>(b) we consider only interlayer coupling Δ_S= 0, Δ_B≠ 0. For the intralayer coupling, Pspin-singlet f_11^,f_12^,f_11^,f_12^ components are real and have s-wave (∝ k_x^2+k_y^2) symmetry while the other components with triplet Pspin state show p-wave (∝ k_±) and f-wave (∝ k_±^2k_∓) symmetries with both nonzero imaginary and real parts. The interlayer coupling changes this picture slightly. Yet, the correlations with triplet Pspin state are of effective f-wave class while the intralayer Pspin-singlet correlations are purely imaginary with s-wave (∝ k_x^2+k_y^2) symmetry class. The interlayer Pspin-singlet correlations, f_12^ and f_12^, are of s-wave class and even in frequency. 

In short, the AA stacking order with and without doping and in the presence of intralayer or interlayer phonon-mediated electron-electron BCS coupling does not support d-wave symmetry class. Hence, AA stacking order can only yields s-wave, p-wave, and f-wave type superconducting correlations with the phonon-mediated scenarios considered in this paper. This can be fully confirmed by examining the analytic expressions given by Eqs. (<ref>) and (<ref>) in Appendix <ref>. According to the momentum combinations presented in Table <ref>, Eqs. (<ref>) and (<ref>) contain no term of type single k_±^2 that results in d-wave symmetry class. The allowed momentum combinations are either odd-parity or possess s-wave symmetry.






  §.§.§ AB stacking order: phonon-mediated intralayer and interlayer electron-electron couplings


Figure <ref> shows the pairing correlations in a BLG with AB ordering of the top and bottom pristine layers. In Fig. <ref>(a), the BCS electron-electron coupling is of the intralayer type, i.e., Δ_S≠ 0, Δ_B=0, while Fig. <ref>(b) contains results with interlayer BCS electron-electron coupling, i.e., Δ_S= 0, Δ_B≠ 0. Clearly, in Fig. <ref>(a), components f_11^,f_12^,f_11^,f_12^ show a combination of d-wave and s-wave symmetries both in imaginary and real parts. The rest of the pairing components are, however, of the p-wave and f-wave type similar to the AA ordering case we considered before. These observed odd-parity and d-wave (proportional to single k_±^2) symmetries can be confirmed by the analytic expressions in Eqs. (<ref>). In Fig. <ref>(b), for the interlayer BCS coupling, the only component that possesses effective d-wave symmetry is f_11^. Checking through Eqs. (<ref>), we see that f_11^ expression contains single k_+^2 which produces the d-wave symmetry. The components f_11^,f_12^,f_12^ are of effective s-wave type while f_11^,f_12^,f_11^,f_12^ show p-wave symmetry beside terms proportional to k_+k_-^2k_+^2 and k_+^2k_-. Mathematically, the reason that f_11^ gains d-wave symmetry is the specific combination of the interlayer coupling T̃ which has off-diagonal terms, and the phonon-mediated electron-electron coupling Δ_B that is on the off-diagonal entries of Eq. (<ref>).

Our prediction of the symmetry changes and generation of Pspin superconducting correlations can be confirmed experimentally. A relevant experiment that can reveal distinctive and direct evidence for the prediction of symmetry change is a high-resolution angular point-contact tunneling spectroscopy experiment. The angular tunneling spectroscopy experiment can determine the angular dependence of superconducting correlations before and after introducing the displacement and twist into the layers. Therefore, one should be able to distinguish between angle-independent and angle-dependent symmetry classes. Another method that might be utilized in investigating the change in the symmetry profile of superconducting correlations discussed above is nonlinear Meissner effect to image the symmetry profile of pairing correlations <cit.>. 




§ SHIFTED AND TWISTED BILAYER GRAPHENE


To study twisted BLG, one can assume that the top and bottom layers are rotated mutually by an amount of θ/2, one clockwise and the other counter clockwise. Therefore, the uncoupled single layers can be described by H_1,2^S(±θ/2,θ_k)=Θ_z(±θ/2)H^S(θ_k)Θ_z^†(±θ/2), where θ_k is the direction of the particles' momentum in each layer and Θ_z(θ) is the rotation operator around the z axis perpendicular to the graphene plane. In this case, the twist affects the coupling matrix as well so that one finds T̃(θ, 𝐫,k) from Eq. (<ref>), now with the phase factor exp (ig· 𝐫) in which g_j=Θ_z(+θ/2)G_j - Θ_z(-θ/2)G_j. 
For commensurate rotation angles, e.g., θ=21.787^∘, 38.213^∘ one finds an effective 4× 4 Hamiltonian in the momentum space for the twisted BLG <cit.>. As the twisting angle becomes smaller, the superlattice cell gets larger. Note that for the noncommensurate angles, the bilayer system is inhomogeneous in real space and a very large basis set is needed to describe it <cit.>. Also, an AA-stacked bilayer can be transformed into an AB stacking order by keeping layer 1 fixed, i.e. u_1=0, and shifting layer 2 by u_2= ϵ( a_1 +  a_2)/3 in Eq. (<ref>). One can easily show that this transformation changes AA order into AB, AB into AC, AC into AA when ϵ is equal to unity.

Nevertheless, in what follows, we consider a generic matrix with four different components for the coupling term, namely:

    T̃=([ t_11 t_12; t_21 t_22 ]),

and derive the anomalous Green's function, and finally, investigate a specific case of shifted BLG. Note that our analytic results below are applicable to twisted BLG where the low-energy physics can be described by the effective Hamiltonian approach. One simply needs to replace appropriate coupling terms (i.e., t_11, t_12, t_21, t_22), describing a twist.


As an illustrative example, we consider the pristine case,  μ=0, and only include the intralayer coupling: Δ_S≠ 0, Δ_B= 0. The results for the anomalous Green's function are:


    Ω f_11^= 0,
    
          Ω f_11^=-Δ_S(Δ_S ^2+k_- k_++t_11
       t_22-t_12 t_21+ω_n^2),
    
          Ω f_12^=-Δ_S  k_- (t_22+t_22^*)-i Δ_Sω_n
       (t_12+t_21^*),
    
          Ω f_12^=-Δ_S  (-k_- t_21+k_+ t_21^*-i ω_n
       (t_11-t_22^*)),
    
          Ω f_11^=+Δ_S(Δ_S ^2+k_- k_++t_11
       t_22-t_12 t_21+ω_n^2),
    
          Ω f_11^=0,
    
          Ω f_12^=-Δ_S  (-k_- t_12^*+k_+ t_12-i ω_n
       (t_11^*-t_22)),
    
          Ω f_12^=+Δ_S  k_+ (t_11+t_11^*)+i Δ_Sω_n
       (t_12^*+t_21),



    Ω=k_-^2 (k_+^2-t_12^* t_21)+k_-
       (-k_+ (-2 (Δ_S ^2+ω_n^2)+ .. 
          .. t_11 t_22^*+t_11^* t_22)-i
       ω_n (t_12^* (t_11+t_22)+  t_21
       (t_11^*+t_22^*)))  
            -k_+^2 t_12 t_21^*-i
       k_+ ω_n (t_21^*
       (t_11+t_22)+t_12
       (t_11^*+t_22^*))+ 
          ω_n^2 (2 Δ_S
       ^2+t_11 t_11^*+t_12 t_12^*+t_21
       t_21^*+    t_22 t_22^*)-   
          (-Δ_S ^2-t_11
       t_22+t_12 t_21)   (Δ_S ^2+t_11^*
       t_22^*-t_12^* t_21^*)+ω_n^4.

An examination of these equations leads to the following conclusions. Similar to Sec. <ref>, the generic coupling matrix leads to vanishing f_11^, f_11^. However, the symmetries are no longer similar to those explored in Sec. <ref>. As seen, here, Ω includes terms of single k_±^2 that makes substantial changes. These terms have d-wave symmetry and thus the resultant Green's functions and corresponding superconducting correlations can now carry this symmetry as well. 

Next, we consider the interlayer coupling (Δ_S=0, Δ_B 0 at μ=0). The pairing correlations are given in Eqs. (<ref>):


    Ω f_11^= 0,
    
          Ω f_11^=Δ_B  (k_- t_21+k_+ t_12+i ω_n
       (t_11+t_22)),
    
          Ω f_12^=i Δ_B  (t_21^* t_22-t_12 t_22^*),
    
          Ω f_12^=-Δ_B(Δ_B ^2+k_- k_+-t_11
       t_22^*+t_21 t_21^*+ω_n^2),
    
          Ω f_11^=Δ_B  (k_- t_21+k_+ t_12+i ω_n
       (t_11+t_22)),
    
          Ω f_11^=0,
    
          Ω f_12^=Δ_B(Δ_B ^2+k_- k_+-t_11^*
       t_22+t_12 t_12^*+ω_n^2),
    
          Ω f_12^=-i Δ_B  (t_11 t_12^*-t_11^* t_21),
    
          Ω=-Δ_B ^4-k_-^2 k_+^2+k_-^2 t_12^* t_21-2
       Δ_B ^2 k_- k_+-  
          ω_n^2 (2 Δ_B ^2+ 2
       k_- k_++t_11 t_11^*+t_12
       t_12^*+t_21 t_21^*+t_22 t_22^*)+   
            i
       ω_n (k_- t_12^*
       (t_11+t_22)+  k_- t_21
       (t_11^*+t_22^*)+   k_+ t_21^*
       (t_11+t_22)+ 
            k_+ t_12
       (t_11^*+t_22^*))+k_- k_+ t_11
       t_22^*+  k_- k_+   t_11^* t_22+k_+^2
       t_12 t_21^*-  
            t_11 t_11^* t_22
       t_22^*+ t_11 t_12^* t_21^* t_22+Δ_B ^2
       t_11 t_22^*+ t_11^* t_12 t_21 t_22^*+Δ_B
       ^2 t_11^* t_22
           -Δ_B ^2 t_12 t_12^*- t_12
       t_12^* t_21 t_21^*-Δ_B ^2 t_21
       t_21^*-ω_n^4.


We see that here, also f_11^, f_11^ are zero. Similarly to the previous case, Ω includes single k_±^2 terms, carrying d-wave symmetry. Note that these single k_±^2 terms are not vanishing here because all the entries of coupling matrix T̃ are assumed nonzero with no specific symmetry. In all the previous cases, with specific stacking orders, Ω includes powers of k_± k_∓ that carry s-wave symmetry.




At finite doping, the expressions become cumbersome, and we resort to numerical investigations. To examine the effect of displacement on the effective symmetry profile of pairing correlations, we study the shifted BLG scenario numerically. Figures <ref> and <ref> are shifted (ϵ=0.2) counterparts of Figs. <ref> and <ref>. When comparing these figures, we see that, for AA stacking, the shift has induced imaginary parts into f_11^, f_11^, f_12^, f_12^ pairings with intralayer coupling and real parts for the interlayer coupling. Also, in the case of AB stacking, this displacement gives rise to nonzero real parts into f_12^, f_12^ with Δ_S≠ 0 and nonzero imaginary parts into f_11^, f_12^ with Δ_B≠ 0. In all cases, the displacement tends to first destroy the symmetries and, subsequently, at large values of ϵ recovers either AA, AB, or AC symmetries depending on the initial ordering we start with. In the next section, we see how this deformation affects the critical temperature and superconducting gap and link them to the superconducting correlations discussed so far. 

As mentioned in passing, inspired by recent experiments in normal state<cit.>, BLG may develop a “two-gap” superconductivity if the electron-electron coupling strength within the layers is unequal to the electron-electron coupling between the two layers. Following this idea, we have calculated the anomalous Green's function when both Δ_S and Δ_B are nonzero. The results of this two-gap superconductivity scenario are given by Eqs. (<ref>) of the Appendix. In deriving  Eqs. (<ref>) we consider bulk properties , i.e., Δ=Δ^† and assume that the two layers develop identical superconducting gap Δ_S, which is unequal to the interlayer superconducting gap Δ_B. To be specific, we consider twisted bilayers by an amount θ with zero doping, i.e., μ=0. The intralayer Pspin-triplet components remains zero. Similarly to the single gap superconductivity case, Ω contains k_±^2 terms due to nonzero coupling potential terms considered.




 §.§ Temperature dependence of superconducting gaps 


Recalling Eqs. (<ref>)-(<ref>), the AB-stacked graphene bilayer with finite doping can potentially support a number of different pairings: s-wave, p-wave, d-wave, f-wave, and their combinations. In addition to the anomalous Green's function, an interaction potential V(k,k') enters the equations, determining the superconducting gap function and critical temperature. In general, the superconducting gap function can be momentum dependent and it satisfies the following equation that should be solved self-consistently together with Eq. (<ref>):

    Δ^αα'_ββ'(k)=-∑_k' V^αα',σσ'_ββ',ρρ'(k,k') f^σσ'_ρρ'(k'),

in which α, α', β, β' are Pspin and layer indexes, respectively. The interaction potential can be both spin and layer dependent. Therefore, one can expand the interaction potential in spherical harmonics and in terms of the atomic orbital symbols introduced in Table <ref>.


The type and details of phonon-mediated electron-electron interactions can be obtained through optical absorption and Raman spectroscopy experiments <cit.>. With the relevant anomalous Green's functions at hand, Eqs. (<ref>)-(<ref>), it is sufficient to insert the interaction potential into Eq. (<ref>) to obtain a self-consistent equation for the gap function, and solve it at different temperatures. As an example, we consider s-wave pairing (i.e., V_s in Table <ref>) and plot the superconducting gap as a function of temperature for different scenarios presented in Fig. <ref>. Note that we have used an identical interaction potential for intralayer and interlayer coupling potential in AA and AB orderings and set μ/t=1.5 throughout our calculations (increasing the density of charged particles enhances the superconducting gap and critical temperature). With these assumptions, we find that the intralayer phonon-mediated electron-electron coupling in AA ordering hosts the largest superconducting gap that, from energy point of view, results in the lowest ground-state energy. We thus normalize all other superconducting gaps and temperature by the intralayer zero temperature superconducting gap of AA ordering, i.e., Δ_S^AA(T=0) and its critical temperature T_c,S^AA, respectively. Figure <ref> illustrates that the interlayer phonon-mediated electron-electron coupling in AA ordering results in the smallest critical temperature (Fig. <ref>(a)) while its AB ordering counterpart Fig. <ref>(d) has the closest gap amplitude and critical temperature to the intralayer AA ordering (the largest one). We have also applied a small displacement, i.e., ϵ=0.2, according to the discussion in the previous subsection. Compairing Figs. <ref>(a) and <ref>(d), we find that the intralayer and interlayer couplings in AA and AB stacking orders share similar behavior in all features in the presence and absence of the displacement. While the displacement in intralayer coupling of AA stacking, Fig. <ref>(a), decreases the gap amplitude and critical temperature, it enhances these quantities in AB ordering, Fig. <ref>(c). This finding is reversed for the interlayer BCS coupling scenario in AA configuration. The real parts of f_11^, f_11^ and f_12^, f_12^ play a crucial role resulting in dramatic changes to the gap function and critical temperature. The underlying reason for these changes can be understood by compairing Figs. <ref> and <ref> with Figs. <ref> and <ref>. The displacement ϵ induces real parts to f_11^ and f_12^ with interlayer and intralyer phonon-mediated electron-electron coupling of AA and AB orderings, respectively, Figs. <ref>(b) and <ref>(c). However, a nonzero displacement ϵ applied to AA and AB orderings with intralayer and interlayer BCS coupling, respectively, induces nonzero imaginary parts to f_11^ and f_12^ that changes slightly the gap function and critical temperature, Figs. <ref>(a) and <ref>(d).

One can repeat the study presented in Fig. <ref> for other interaction types given in Table <ref> and the given Green's functions in Eqs. (<ref>)-(<ref>). Nevertheless, we postpone such a study to future works when more details of interactions or relevant experimental data are available. It is, however, apparent that considering d-wave or g-wave interactions in Table <ref> for the electron-electron coupling potential V(k,k') with its spin state being singlet, the ratio of superconducting gap and critical temperature increases (according to the traditional calculations of such interaction potentials <cit.>) for cases where the Green's function itself possesses d-wave symmetry, as discussed in the previous section. Therefore, such a discussion should be relevant for the recent experiment <cit.> that observed a puzzling behavior combining the features of high T_c and conventional superconductors. Accordingly, by deforming the BLG, d-wave spin-singlet pairing may become dominant and increase the ratio of the superconducting gap and T_c.


In our calculations above, we have simply followed the BCS picture of superconductivity where the coupling of two particles with opposite spins and momenta are building blocks of superconductivity. Additionally, we have assumed that the influence of twist and displacement in superconducting BLG are encoded in the coupling terms of Eq. (<ref>) and the electron-electron coupling interaction potentials remain unchanged. Nevertheless, if twist and displacement change the amplitude of electron-electron coupling interaction potentials only, the above findings remain intact. However, if twist and displacement introduce momentum-dependent changes to the electron-electron coupling interaction potentials, similar to those presented in Table <ref>, one should repeat the above calculations with these momentum-dependent potentials. 




§ SUMMARY


In summary, using an effective Hamiltonian approach to discuss in-plane displacements between graphene monolayers in a bilayer graphene system, we investigate the superconducting pairing correlations by deriving the anomalous Green's function. Motivated by a recent experiment, we consider both intralayer and interlayer phonon-mediated electron-electron couplings in AA and AB stacking orders. Our results reveal that both AA and AB configurations at the charge neutrality point, μ=0, can only develop even-parity s-wave and odd-parity p-wave superconducting correlations. At a finite doping, μ≠ 0, this finding remains intact in a AA system with the addition of odd-parity f-wave symmetry while AB ordering, exclusively, can host even-parity d-wave symmetries. Introducing a generic coupling potential between graphene monolayers in bilayer graphene, we show that displacement of graphene monolayers can induce d-wave symmetry at μ=0 as well. Our results suggest that a switching to d-wave symmetry can be achieved by a simple displacement of the two coupled pristine graphene layers at a finite density of charged carriers. We also discuss the possible appearance of pseudospin-triplet and odd-frequency pairings. Finally, we consider an s-wave interaction potential and study the superconducting gap function and critical temperature in a specific scenario of displacement (in-plane shifting of the monolayers of graphene) in AA and AB stacking orders. We find that AA stacking with an intralayer electron-electron coupling most favorably hosts s-wave superconductivity while AA ordering with interlayer electron-electron coupling is the next desirable platform. Also, the exertion of a slight in-plane shift in bilayer graphene can increase both the amplitude of the superconducting gap and critical temperature in AA and AB orderings with phonon-mediated interlayer and intralayer electron-electron couplings.



M.A. is supported by Iran's National Elites Foundation (INEF). Center for Nanostructured Graphene is supported by the Danish National Research Foundation (Project No. DNRF103).







§ HAMILTONIAN OF SLG WITH SUPERCONDUCTIVITY


A monolayer graphene can be described by two valleys in k-space: K =2π/3a(1,+1/√(3)) and K'=2π/3a(1,-1/√(3)). At low energies, the dispersion relation near the corners of the Brillouin zone, i.e., K and K' are linear and can be described by a Dirac Hamiltonian: H_ K(k) = ħ v_Fk·σ.
The corresponding wavefunction is a two-component spinor, namely,

    ħ v_F([        0 k_x-ik_y; k_x+ik_y        0 ] )([ ψ_A K(k); ψ_B K(k) ] )=ε([ ψ_A K(k); ψ_B K(k) ] ).
 
Here, A and B denote the sublattices of graphene honeycomb lattice. Note that the Hamiltonian above couples these two sublattices. This appears because the nearest neighbor site of each sublattice site belongs to the other sublattice. This can be clearly seen in Fig. <ref>(a), where sublattices are marked by blue and red colors. Hence, the total spin-less Hamiltonian can be expressed by:

    H(k) = ([  H_ K(k)        0;        0 H_ K'(k) ] ),
 
where H_ K'(k) = ħ v_Fk·σ^*. Next, we can add spin to this Hamiltonian by introducing Pauli matrices acting in spin-space: τ=(τ_x,τ_y,τ_z). In this case, the spin-full Hamiltonian can be expressed by H_tot=H(k)τ_0 where τ_0 reflects the absence of spin coupling. In SLG the two valleys are decoupled in the absence of an external potential<cit.>. Note that the form Eq. (<ref>) applies even in the presence of elastic impurities and disorder. Therefore, the Hamiltonian becomes degenerated and reduces to a 4× 4 Hamiltonian. Now we incorporate superconductivity by introducing a two-electron amplitude,

    Δ_S⟨ψ^†_A↑ Kψ^†_B↓ K'⟩+H.c.,

in which Δ_S is the gap representing BCS spin-singlet phonon-mediated electron-electron coupling between sublattices A, B and valleys K, K'. Note that there are several options to introduce a two-electron amplitude in the presence of superconductivity, including Pspin-triplet coupling that couples two electrons in the same sublattice <cit.>. In the main text, we have used the above two-electron amplitude Eq. (<ref>) in our calculations. It is worth mentioning that the hole-excitation block of the Hamiltonian in the particle-hole space, when introducing superconductivity, is equivalent to the Hamiltonian of K' valley. Therefore, superconductivity can be viewed as a mean to introduce valley coupling in SLG. When two normal SLGs are coupled (making a BLG), more options for incorporating the two-electron amplitudes in its superconducting phase are generated. To clarify the presence of the second SLG, we have introduced indices 1,2 that label top and bottom SLGs as described in the main text. In this case, to be consistent with the previous case, we consider the following two-electron amplitude:   

    Δ_B⟨ψ^†_1A↑ Kψ^†_2B↓ K'⟩+H.c.,

where Δ_B is the gap representing BCS spin-singlet phonon-mediated electron-electron coupling between sublattices A, B, valleys K, K', and layers 1,2. To simplify our notation throughout the presentation in the main text, due to the spin-degeneracy<cit.>, we have dropped spin and valley K, K' indices and only keep the sublattice indices. Since we discuss AA and AB stacking orders, from now on, we change our notation of A,B sublattices to ↑,↓ and call them pseudospin (Pspin). 
 
 
 

§ GREEN'S FUNCTION OF DISPLACED BLG WITH SUPERCONDUCTIVITY

In this part, we present expressions derived for the anomalous Green's function in the presence of a finite doping, considering intralayer and interlayer phonon-mediated electron-electron couplings.




∙Intralayer BCS coupling, Δ_S≠ 0, Δ_B= 0, in AA stacking order


    Ω f_11^= +2 Δ_S  k_- μ((Δ_S ^2+μ ^2+ω_n^2)^2+k_-^2 k_+^2+2 k_- k_+ (Δ_S ^2-μ
       ^2+t^2+ω_n^2)-  3 t^4-2 t^2 (Δ_S ^2-μ ^2+ω_n^2)),
    
       Ω f_11^=-Δ_S(k_-^3 k_+^3-k_-^2 k_+^2 (-3 Δ_S ^2+μ ^2+t^2-3
       ω_n^2)+k_- k_+ ((Δ_S ^2+μ ^2+ω_n^2) (3 Δ_S ^2-μ ^2+3 ω_n^2)-t^4+2 t^2 (Δ_S
       ^2+5 μ ^2+ω_n^2))+ . 
       . (Δ_S ^2+(t-μ )^2+ω_n^2) (Δ_S ^2+μ ^2+t^2+ω_n^2) (Δ_S ^2+(μ
       +t)^2+ω_n^2)),
    
       Ω f_12^=-2 Δ_S  k_- t ((Δ_S ^2-3 μ ^2+ω_n^2) (Δ_S
       ^2+μ ^2+ω_n^2)+k_-^2 k_+^2+2 k_- k_+
       (Δ_S ^2+μ ^2-t^2+ω_n^2)+t^4+2 t^2 (Δ_S ^2+μ
       ^2+ω_n^2)),
    
       Ω f_12^=-2 Δ_Sμ  t (-3 k_-^2 k_+^2+2 k_- k_+ (-Δ_S ^2+μ
       ^2+t^2-ω_n^2)+(Δ_S ^2+(t-μ )^2+ω_n^2)
       (Δ_S ^2+(μ +t)^2+ω_n^2)),
    
       Ω f_11^=-Δ_S(-k_-^3 k_+^3+k_-^2 k_+^2 (-3 Δ_S ^2+μ ^2+t^2-3
       ω_n^2)+k_- k_+ (-(Δ_S ^2+μ ^2+ω_n^2) (3 Δ_S ^2-μ ^2+3 ω_n^2)+t^4-2 t^2 (Δ_S
       ^2+5 μ ^2+ω_n^2))-  . 
       .  (Δ_S ^2+(t-μ )^2+ω_n^2) (Δ_S ^2+μ ^2+t^2+ω_n^2) (Δ_S ^2+(μ
       +t)^2+ω_n^2)),
    
       Ω f_11^=-2 Δ_S  k_+ μ((Δ_S ^2+μ ^2+ω_n^2)^2+k_-^2
       k_+^2+2 k_- k_+ (Δ_S ^2-μ ^2+t^2+ω_n^2)-3
       t^4-2 t^2 (Δ_S ^2-μ ^2+ω_n^2)),
    
       Ω f_12^=+2 Δ_Sμ  t (-3 k_-^2 k_+^2+2 k_- k_+ (-Δ_S ^2+μ
       ^2+t^2-ω_n^2)+(Δ_S ^2+(t-μ )^2+ω_n^2)
       (Δ_S ^2+(μ +t)^2+ω_n^2)),
    
       Ω f_12^=+2 Δ_S  k_+ t ((Δ_S ^2-3 μ ^2+ω_n^2) (Δ_S
       ^2+μ ^2+ω_n^2)+k_-^2 k_+^2+2 k_- k_+
       (Δ_S ^2+μ ^2-t^2+ω_n^2)+t^4+2 t^2 (Δ_S ^2+μ
       ^2+ω_n^2)),
       Ω=(k_-^2 k_+^2+2 k_- k_+ (Δ_S ^2-(t-μ )^2+ω_n^2)+(Δ_S ^2+(t-μ )^2+ω_n^2)^2)
       (k_-^2 k_+^2-2 k_- k_+ (-Δ_S ^2+(μ
       +t)^2-ω_n^2)+(Δ_S ^2+(μ +t)^2+ω_n^2)^2),




∙Intralayer BCS coupling, Δ_S≠ 0, Δ_B= 0, in AB stacking order


    Ω f_11^= +2 Δ_Sμ(k_-^3 k_+^2+2 k_-^2 k_+ (Δ_S ^2-μ
       ^2+ω_n^2)+k_- (Δ_S ^2+μ ^2+ω_n^2)^2-k_+ t^2 (μ ^2+ω_n^2)),
    
       Ω f_11^=-Δ_S(k_-^3 k_+^3+k_-^2 k_+^2 (3 Δ_S ^2-μ ^2+3
       ω_n^2)+k_- k_+ ((Δ_S ^2+μ ^2+ω_n^2) (3 Δ_S ^2-μ ^2+3 ω_n^2)+t^2 (μ
       ^2+ω_n^2))-   2 i μ  t^2 ω_n (Δ_S
       ^2+k_+^2+μ ^2)-.
       . 2 k_+^2 μ ^2 t^2+ω_n^4 (3
       (Δ_S ^2+μ ^2)+t^2)+ω_n^2 (3 (Δ_S ^2+μ
       ^2)^2+Δ_S ^2 t^2)-   2 i μ  t^2 ω_n^3+(Δ_S ^2+μ
       ^2) (Δ_S ^2+μ  (μ -t)) (Δ_S ^2+μ  (μ
       +t))+ω_n^6),
    
       Ω f_12^=-Δ_S  t (-2 k_-^3 k_+ μ +k_-^2 (2 μ(-Δ_S ^2+(μ
       +i ω_n)^2)+k_+^2 (μ +i ω_n))+2 i k_-
       k_+   ω_n (Δ_S ^2+μ ^2+ω_n^2)- .
       . (μ -i
       ω_n) (Δ_S ^2+μ ^2-t (μ +i ω_n)+ω_n^2) (Δ_S ^2+μ  (μ +t)+i t ω_n+ω_n^2)),
    
       Ω f_12^=+2 Δ_Sμ  t (-k_-^2 k_+ (μ -i ω_n)+k_- (μ +i
       ω_n) (Δ_S ^2+k_+^2+μ ^2+ω_n^2)-k_+
       (μ -i ω_n) (Δ_S ^2+μ ^2+ω_n^2)),
    
       Ω f_11^ =-Δ_S(-k_-^3 k_+^3+k_-^2 k_+^2 (-3 Δ_S ^2+μ ^2-3
       ω_n^2)-k_- k_+ ((Δ_S ^2+μ ^2+ω_n^2) (3 Δ_S ^2-μ ^2+3 ω_n^2)+t^2 (μ
       ^2+ω_n^2))-   2 i μ  t^2 ω_n (Δ_S
       ^2+k_+^2+μ ^2).
       . +2 k_+^2 μ ^2 t^2-ω_n^4 (3
       (Δ_S ^2+μ ^2)+t^2)-ω_n^2 (3 (Δ_S ^2+μ
       ^2)^2+Δ_S ^2 t^2)-   2 i μ  t^2 ω_n^3-(Δ_S ^2+μ
       ^2) (Δ_S ^2+μ  (μ -t)) (Δ_S ^2+μ  (μ
       +t))-ω_n^6),
    
       Ω f_11^=-2 Δ_Sμ(k_-^2 k_+^3+k_- k_+^2 (2 (Δ_S ^2-μ
       ^2+ω_n^2)+t^2)+Δ_S ^2 k_- t^2-k_+^3 t^2-k_+
       (Δ_S ^2+μ ^2+ω_n^2) (-Δ_S ^2-μ ^2+t^2-ω_n^2)),
    
       Ω f_12^=-Δ_S  t (2 i μω_n (k_-+k_+) (Δ_S ^2+k_-
       k_++μ ^2)+2 i μω_n^3
       (k_-+k_+)+(k_+-k_-) ((Δ_S ^2+k_-
       k_+)^2  -μ ^4+μ ^2 t^2)- .
       .  ω_n^2 (k_–k_+)
       (2 Δ_S ^2+2 k_- k_++t^2)+ω_n^4
       (k_+-k_-)),
    
       Ω f_12^=-Δ_S  t (k_-^2 k_+^2 (μ -i ω_n)-2 k_- k_+
       (k_+^2 μ +i ω_n (Δ_S ^2+μ ^2+ω_n^2))-2 k_+^2 μ(Δ_S ^2-(μ -i ω_n)^2)- .
       .  (μ +i ω_n) (Δ_S ^2+μ  (μ +t)-i t ω_n+ω_n^2) (Δ_S ^2+μ ^2-μ  t+i t
       ω_n+ω_n^2)),
       Ω=Δ_S ^8+4 Δ_S ^6 (μ ^2+ω_n^2)+k_-^4 k_+^4+4
       k_-^3 k_+^3 (Δ_S ^2-μ ^2+ω_n^2)+2 k_-^2
       (k_+^2 (3 Δ_S ^4-2 Δ_S ^2 μ ^2+3 μ ^4+ω_n^2 (6
       Δ_S ^2-2 μ ^2+t^2)-μ ^2 t^2+3 ω_n^4) . 
       . -2 Δ_S ^2 μ ^2
       t^2)+4 k_- k_+ (Δ_S ^2+μ ^2+ω_n^2)
       ((Δ_S ^2+ω_n^2)^2-μ ^4+t^2 (μ ^2+ω_n^2))+4 Δ_S ^2 (-k_+^2 μ ^2 t^2+(μ ^2+ω_n^2)^3+t^2 (ω_n^4-μ ^4))  
        +2 Δ_S ^4 (3
       (μ ^2+ω_n^2)^2+t^2 (ω_n^2-μ
       ^2))+(μ ^2+ω_n^2)^2 ((t-μ )^2+ω_n^2) ((μ +t)^2+ω_n^2),




∙Interlayer BCS coupling, Δ_S= 0, Δ_B≠ 0, in AA stacking order


    Ω f_11^= 8 i Δ_B  k_- μ  t ω_n (Δ_B ^2+k_- k_++μ
       ^2-t^2+ω_n^2),
    
       Ω f_11^=2 i Δ_B  t ω_n (k_-^2 k_+^2+2 k_- k_+
       (Δ_B ^2+3 μ ^2-t^2+ω_n^2)+2 ω_n^2 (Δ_B
       ^2+μ ^2+t^2)+(Δ_B ^2+μ ^2-t^2)^2+ω_n^4),
    
       Ω f_12^=2 Δ_B  k_- μ((Δ_B ^2+μ ^2+ω_n^2)^2+k_-^2
       k_+^2+2 k_- k_+ (Δ_B ^2-μ ^2-t^2+ω_n^2)+t^4-2 t^2 (Δ_B ^2+μ ^2+3 ω_n^2)),
    
       Ω f_12^=-Δ_B(Δ_B ^2+k_- k_++μ ^2-t^2+ω_n^2)
       (k_-^2 k_+^2+2 k_- k_+ (Δ_B ^2-μ
       ^2-t^2+ω_n^2)+2 ω_n^2 (Δ_B ^2+μ
       ^2+t^2)+(Δ_B ^2+μ ^2-t^2)^2+ω_n^4),
    
       Ω f_11^=2 i Δ_B  t ω_n (k_-^2 k_+^2+2 k_- k_+
       (Δ_B ^2+3 μ ^2-t^2+ω_n^2)+2 ω_n^2 (Δ_B
       ^2+μ ^2+t^2)+(Δ_B ^2+μ ^2-t^2)^2+ω_n^4),
    
       Ω f_11^=-8 i Δ_B  k_+ μ  t ω_n (Δ_B ^2+k_- k_++μ
       ^2-t^2+ω_n^2),
    
       Ω f_12^=Δ_B(Δ_B ^2+k_- k_++μ ^2-t^2+ω_n^2)
       (k_-^2 k_+^2+2 k_- k_+ (Δ_B ^2-μ
       ^2-t^2+ω_n^2)+2 ω_n^2 (Δ_B ^2+μ
       ^2+t^2)+(Δ_B ^2+μ ^2-t^2)^2+ω_n^4),
    
       Ω f_12^=-2 Δ_B  k_+ μ((Δ_B ^2+μ ^2+ω_n^2)^2+k_-^2 k_+^2+2 k_- k_+ (Δ_B ^2-μ
       ^2-t^2+ω_n^2)+t^4-2 t^2 (Δ_B ^2+μ ^2+3 ω_n^2)),
       
       Ω=k_-^4 k_+^4+4 k_-^3 k_+^3 (Δ_B ^2-μ ^2-t^2+ω_n^2)+2 k_-^2 k_+^2 (3 Δ_B ^4-2 Δ_B ^2 μ ^2+3 μ ^4+3
       t^4-2 ω_n^2 (-3 Δ_B ^2+μ ^2+t^2)-6 Δ_B ^2 t^2+2 μ ^2
       t^2+3 ω_n^4)
       +4 k_- k_+ ((Δ_B ^2-μ
       ^2+ω_n^2) (Δ_B ^2+μ ^2+ω_n^2)^2-t^6+t^4
       (3 Δ_B ^2+μ ^2-ω_n^2)+t^2 (-3 Δ_B ^4-2 Δ_B ^2
       (μ ^2+ω_n^2)+..
       ..μ ^4+10 μ ^2 ω_n^2+ω_n^4))+(2 ω_n^2 (Δ_B ^2+μ
       ^2+t^2)+(Δ_B ^2+μ ^2-t^2)^2+ω_n^4)^2,




∙Interlayer BCS coupling, Δ_S= 0, Δ_B≠ 0, in AB stacking order


    Ω f_11^= -2 Δ_Bμ  t ((Δ_B ^2+μ ^2+ω_n^2)^2-k_-^2
       k_+^2),
    
       Ω f_11^=-Δ_B  k_+ t (k_-^2 k_+^2+2 k_- k_+ (Δ_B ^2+μ
       ^2+ω_n^2)+(Δ_B ^2+μ ^2+ω_n^2)
       (Δ_B ^2-3 μ ^2+t^2+ω_n^2)),
    
       Ω f_12^=2 Δ_B  k_- μ((Δ_B ^2+μ ^2+ω_n^2)^2+k_-^2
       k_+^2+k_- k_+ (2 (Δ_B ^2-μ ^2+ω_n^2)+t^2)),
    
       Ω f_12^=-Δ_B(Δ_B ^2+k_- k_++μ ^2+ω_n^2)
       (k_-^2 k_+^2+2 k_- k_+ (Δ_B ^2-μ ^2+ω_n^2)+(Δ_B ^2+μ ^2+ω_n^2) (Δ_B ^2+μ
       ^2+t^2+ω_n^2)),
    
       Ω f_11^=Δ_B  k_+ t (k_-^2 k_+^2+2 k_- k_+ (Δ_B ^2+μ
       ^2+ω_n^2)+(Δ_B ^2+μ ^2+ω_n^2)
       (Δ_B ^2-3 μ ^2+t^2+ω_n^2)),
    
       Ω f_11^=-2 Δ_B  k_+^2 μ  t (2 (Δ_B ^2-μ ^2+ω_n^2)+2
       k_- k_++t^2),
    
       Ω f_12^=Δ_B(k_-^3 k_+^3+k_-^2 k_+^2 (3 Δ_B ^2-μ ^2+t^2+3
       ω_n^2)+k_- k_+ (Δ_B ^2+μ ^2+ω_n^2) (3 Δ_B ^2-μ ^2+3 t^2+ 3 ω_n^2)+.
       . (Δ_B
       ^2+μ ^2+ω_n^2) (Δ_B ^2+(t-μ )^2+ω_n^2)
       (Δ_B ^2+(μ +t)^2+ω_n^2)),
    
       Ω f_12^=-2 Δ_B  k_+ μ((Δ_B ^2+μ ^2+ω_n^2)^2+k_-^2 k_+^2+k_- k_+ (2 (Δ_B ^2-μ
       ^2+ω_n^2)+t^2)),
       Ω=k_-^4 k_+^4+4 k_-^3 k_+^3 (Δ_B ^2-μ ^2+ω_n^2)+2 k_-^2 k_+^2 (3 Δ_B ^4-2 Δ_B ^2 (μ ^2-3
       ω_n^2)+3 μ ^4-2 μ ^2 ω_n^2+t^2 (Δ_B ^2-μ
       ^2+ω_n^2)+3 ω_n^4)+
        4 k_- k_+
       (Δ_B ^2+μ ^2+ω_n^2)^2 (Δ_B ^2-μ
       ^2+t^2+ω_n^2)+(Δ_B ^2+μ ^2+ω_n^2)^2
       (Δ_B ^2+(t-μ )^2+ω_n^2) (Δ_B ^2+(μ
       +t)^2+ω_n^2),




∙Intralayer and interlayer BCS coupling, Δ_S≠ 0, Δ_B≠ 0, with a generic coupling matrix


    Ω f_11^= 0,
    
       Ω f_11^=-Δ_B^2 Δ_S-Δ_B (k_-^- t_21+k_+^- t_12+i
       ω_n (t_11+t_22))-  Δ_S(Δ_S^2+k_-^-
       k_+^-+t_11 t_22-t_12 t_21+ω_n^2),
    
       Ω f_12^=-(Δ_S k_-^+ t_22 - Δ_B t_21^* t_22 + Δ_S k_-^- t_22^* + Δ_B t_12 t_22^* +
     i (Δ_B (-k_-^- + k_-^+) + Δ_S (t_12 + t_21^*)) ω_n),
    
       Ω f_12^=-(Δ_B^3 + Δ_S (-k_-^+ t_21 + k_+^- t_21^* - i (t_11 - t_22^*) ω_n) +
     Δ_B (Δ_S^2 + k_-^+ k_+^- + t_21 t_21^* - t_11 t_22^* + ω_n^2)),
    
       Ω f_11^=Δ_B^2 Δ_S + Δ_B (k_+^- t_12 + k_-^- t_21 + i (t_11 + t_22) ω_n) +
     Δ_S (Δ_S^2 + k_-^- k_+^- - t_12 t_21 + T1 t_22 + ω_n^2),
    
       Ω f_11^=0,
    
       Ω f_12^=Δ_B^3 + Δ_S (-k_+^+ t_12 + k_-^- t_12^* + i (t_11^* - t_22) ω_n) +
     Δ_B (Δ_S^2 + k_-^- k_+^+ + t_12 t_12^* - t_11^* t_22 + ω_n^2),
    
       Ω f_12^=Δ_S k_+^+ t_11 + Δ_S k_+^- t_11^* - Δ_B t_11 t_12^* + Δ_B t_11^* t_21 +
     i (Δ_B (-k_+^- + k_+^+) + Δ_S (t_12^* + t_21)) ω_n,
       Ω=Δ_B^4 + Δ_S^4 + k_-^- k_-^+ k_+^- k_+^+ - k_-^- k_-^+ t_12^* t_21 - k_+^- k_+^+ t_12 t_21^* +
     t_12 t_12^* t_21 t_21^* +
     Δ_BΔ_S ((k_-^- - k_-^+) (t_12^* + t_21) + (k_+^- - k_+^+) (t_12 + t_21^*)) - 
       
     k_-^+ k_+^- t_11^* t_22 - t_11 t_12^* t_21^* t_22 - k_-^- k_+^+ t_11 t_22^* - t_11^* t_12 t_21 t_22^* +
     t_11 t_11^* t_22 t_22^* -
     i (k_+^- t_11^* t_12 + k_-^- t_11 t_12^* +   k_-^+ t_11^* t_21 + k_+^+ t_11 t_21^* +  k_-^+ t_12^* t_22 + 
       
        k_+^- t_21^* t_22 + k_+^+ t_12 t_22^* + k_-^- t_21 t_22^*) ω_n + (k_-^- k_+^- +
        k_-^+ k_+^+ + t_11 t_11^* + t_12 t_12^* + t_21 t_21^* +
        t_22 t_22^*) ω_n^2 + ω_n^4 +
     Δ_B^2 (2 Δ_S^2 + k_-^+ k_+^- + k_-^- k_+^+ + 
        t_12 t_12^* + t_21 t_21^* - t_11^* t_22 -
        t_11 t_22^* + 2 ω_n^2) +
     Δ_S^2 (k_-^- k_+^- + k_-^+ k_+^+ - t_12 t_21 - t_12^* t_21^* + t_11 t_22 + t_11^* t_22^* +
        2 ω_n^2).


Here we have defined k_±^+=kcos (θ_k+θ/2)± iksin (θ_k+θ/2) and k_±^-=kcos (θ_k-θ/2)± iksin (θ_k-θ/2) according to the main text for a twisted BLG by an amount θ.



000


Butz B. Butz, C. Dolle, F. Niekiel, K. Weber, D. Waldmann, H. B. Weber, B. Meyer, and E. Spiecker, Dislocations in bilayer graphene,  http://dx.doi.org/10.1038/nature12780Nature 505, 533 (2014).

Hass J. Hass, F. Varchon, J. E. Millan-Otoya, M. Sprinkle, N. Sharma, W. A. de Heer, C. Berger, P. N. First, L. Magaud, and E. H. Conrad, Why Multilayer Graphene on 4H-SiC(0001) Behaves Like a Single Sheet of Graphene, http://link.aps.org/doi/10.1103/PhysRevLett.100.125504Phys. Rev. Lett. 100, 125504 (2008).

Kisslinger F. Kisslinger, C. Ott, C. Heide, E. Kampert, B. Butz, E. Spiecker, S. Shallcross, and H. B. Weber, Linear magnetoresistance in mosaic-like bilayer graphene, http://dx.doi.org/10.1038/nphys3368Nat. Phys. 11, 650 (2015).


Borysiuk J. Borysiuk, J. Soltys, and J. Piechota, Stacking sequence dependence of graphene layers on SiC Experimental and theoretical investigation, http://aip.scitation.org/doi/abs/10.1063/1.3585829Journal of Applied Physics 109(9), 093523 (2011).


Liu Z. Liu, K. Suenaga, P. J. F. Harris, and S. Iijima, Open and Closed Edges of Graphene Layers, http://link.aps.org/doi/10.1103/PhysRevLett.102.015501Phys. Rev. Lett. 102, 015501 (2009).

tbg1 J. M. B. Lopes dos Santos, N. M. R. Peres, and A. H. Castro Neto, Graphene Bilayer with a Twist: Electronic Structure
, https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.99.256802Phys. Rev. Lett. 99, 256802 (2007).

tbg-2 E. J. Mele, Commensuration and interlayer coherence in twisted bilayer graphene, https://link.aps. org/doi/10.1103/PhysRevB.81.161405Phys. Rev. B 81, 161405 (2010).


tbg-1 E. Morell, J. D. Correa, P. Vargas, M. Pacheco, and Z. Barticevic, Flat bands in slightly twisted bilayer graphene: Tight-binding calculations, https://link.aps.org/doi/10.1103/ PhysRevB.82.121407Phys. Rev. B 82, 121407 (2010).

tbg0 G. T. Laissardiere, D. Mayou, and L. Magaud, Numerical studies of confined states in rotated bilayers of graphene, https://link.aps.org/doi/10.1103/PhysRevB.86. 125413Phys. Rev. B 86, 125413 (2012).



tbg2  K. Sato, R. Saito, C. Cong, T. Yu, and M. S. Dresselhaus, Zone folding effect in Raman G-band intensity of twisted bilayer graphene, https://journals.aps.org/prb/abstract/10.1103/PhysRevB.86.125414Phys. Rev. B 86, 125414 (2012).

tbg3 P. Moon, M. Koshino, Optical absorption in twisted bilayer graphene, https://journals.aps.org/prb/abstract/10.1103/PhysRevB.87.205404Phys. Rev. B 87, 205404 (2013).

tbg4  A. Righi, P. Venezuela, H. Chacham, S.D. Costa, C. Fantini, R.S. Ruoff, L. Colombo, W.S. Bacsa, M.A. Pimenta, Resonance Raman spectroscopy in twisted bilayer graphene, https://www.sciencedirect.com/science/article/pii/S0038109813002548Solid State Commun. 175, 13 (2013).

tbg5 J. Campos-Delgado, L. G. Canclado, A.C. Achete, A. Jorio, and J.-P. Raskin, Raman scattering study of the phonon dispersion in twisted bilayer graphene, https://link.springer.com/article/10.1007/s12274-013-0304-zNano Res. 6, 269 (2013).

tbg6  Y. Wang, Z. Su, W. Wu, S. Nie, N. Xie, H. Gong, Y. Guo, J. Hwan Lee, S. Xing, X. Lu, H. Wang, X. Lu, K. McCarty, S.-s. Pei, F. Robles-Hernandez, V. G. Hadjiev, and J. Bao, Resonance Raman spectroscopy of G-line and folded phonons in twisted bilayer graphene with large rotation angles, http://aip.scitation.org/doi/abs/10.1063/1.4821434Appl. Phys. Lett. 103, 123101 (2013).

tbg8  C. C. Lu, Y. C. Lin, Z. Liu, C. H. Yeh, K. Suenaga, P. W. Chiu, Twisting bilayer graphene superlattices, http://pubs.acs.org/doi/abs/10.1021/nn3059828ACS Nano 7, 2587 (2013).

tbg9  Z. Ni, L. Liu, Y. Wang, Z. Zheng, L.J. Li, T. Yu, Z. Shen, G-band Raman double resonance in twisted bilayer graphene: evidence of band splitting and folding, https://journals.aps.org/prb/abstract/10.1103/PhysRevB.80.125404Phys. Rev. B 80, 125404 (2009).

tbg10  K. Kim, S. Coh, L. Z. Tan, W. Regan, J. M. Yuk, E. Chatterjee, M. F. Crommie, M. L. Cohen, S. G. Louie, and A. Zettl, Raman spectroscopy study of rotated double-layer graphene: misorientation-angle dependence of electronic structure, https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.108.246103Phys. Rev. Lett. 108, 246103 (2012).

tbg11  R. W. Havener, H. Zhuang, L. Brown, R. G. Hennig, and J. Park, Angle- resolved Raman imaging of interlayer rotations and interactions in twisted bilayer graphene, http://pubs.acs.org/doi/abs/10.1021/nl301137kNano. Lett. 12, 3162 (2012).

tbg12
P. Venezuela, M. Lazzeri, F. Mauri, Theory of double-resonant Raman spectra in graphene: intensity and line shape of defect-induced and two-phonon bands, https://journals.aps.org/prb/abstract/10.1103/PhysRevB.84.035433Phys. Rev. B 84, 035433 (2011).

tbg13  J. Maultzsch, S. Reich, C. Thomsen, H. Requardt, and P. Ordejan, Phonon dispersion in graphite, https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.92.075501Phys. Rev. Lett. 92, 075501 (2004).

tbg14
A.I. Cocemasov, D.L. Nika, A.A. Balandin, Phonons in twisted bilayer graphene, https://journals.aps.org/prb/abstract/10.1103/PhysRevB.88.035428Phys. Rev. B 88, 035428 (2013).

2ph-el_exp G. S. N. Eliel, M. V. O. Moutinho, A. C. Gadelha, A. Righi, L. C. Campos, H. B. Ribeiro, Po-Wen Chiu, K. Watanabe, T. Taniguchi, P. Puech, M. Paillet, T. Michel, P. Venezuela and M. A. Pimenta, Intralayer and interlayer electron–phonon interactions in twisted graphene heterostructures, https://www.nature.com/articles/s41467-018-03479-3Nat. Comm. 9, 1221 (2018).


ray1 D. Weckbecker, S. Shallcross, M. Fleischmann, N. Ray, S. Sharma, and O. Pankratov, Low-energy theory for the graphene twist bilayer, http://link.aps.org/doi/10.1103/PhysRevB.93. 035452Phys. Rev. B 93, 035452 (2016).

ray2 N. Ray, M. Fleischmann, D. Weckbecker, S. Sharma, O. Pankratov, and S. Shallcross, Electron-phonon scattering and in-plane electric conductivity in twisted bilayer graphene, http://link. aps.org/doi/10.1103/PhysRevB.94.245403Phys. Rev. B 94, 245403 (2016).

ray_slip N. Ray, F. Rost, D. Weckbecker, M. Vogl, S. Sharma, R. Gupta, O. Pankratov, and S. Shallcross, Going beyond k.p theory: a general method for obtaining effective Hamil-tonians in both high and low symmetry situations, https://arxiv.org/abs/1607.00920arXiv:1607.00920 (2016).


herrero  Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras and P. Jarillo-Herrero, Unconventional superconductivity in magic-angle graphene superlattices, https://www.nature.com/articles/nature26160Nature 556, 43 (2018).

C.Dean_exp M. Yankowitz, S. Chen, H. Polshyn, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, C. R. Dean, Tuning superconductivity in twisted bilayer graphene, https://arxiv.org/abs/1808.07865arXiv:1808.07865(2018).

grexp1 E. Codecido, Q. Wang, R. Koester, S. Che, H. Tian, R. Lv, S. Tran, K. Watanabe, T. Taniguchi, F. Zhang, M. Bockrath, C. N. Lau, Correlated Insulating and Superconducting States in Twisted Bilayer Graphene Below the Magic Angle, https://arxiv.org/abs/1902.05151arXiv:1902.05151.


grexp2 S. Moriyama, Y. Morita, K. Komatsu, K. Endo, T. Iwasaki, S. Nakaharai, Y. Noguchi, Y. Wakayama, E. Watanabe, D. Tsuya, K. Watanabe, T. Taniguchi, Observation of superconductivity in bilayer graphene/hexagonal boron nitride superlattices, https://arxiv.org/abs/1901.09356arXiv:1901.09356.

grexp3 X. Lu, P. Stepanov, W. Yang, M. Xie, M. A. Aamir, I. Das, C. Urgell, K. Watanabe, T. Taniguchi, G. Zhang, A. Bachtold, A. H. MacDonald, D. K. Efetov, Superconductors, Orbital Magnets, and Correlated States in Magic Angle Bilayer Graphene, https://arxiv.org/abs/1903.06513arXiv:1903.06513.

teor_tbg1  C. Xu and L. Balents, Topological Superconductivity in Twisted Multilayer Graphene, https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.121.087001Phys. Rev. Lett. 121, 087001 (2018).

teor_tbg2  N. F. Q. Yuan and L. Fu, A Model for Metal-Insulator Transition in Graphene Superlattices and Beyond, https://journals.aps.org/prb/abstract/10.1103/PhysRevB.98.045103Phys. Rev. B 98, 045103 (2018).

teor_tbg3  H. C. Po, L. Zou, A. Vishwanath, and T. Senthil, Origin of Mott insulating behavior and superconductivity in twisted bilayer graphene, https://journals.aps.org/prx/abstract/10.1103/PhysRevX.8.031089
Phys. Rev. X 8, 031089 (2018).

teor_tbg4  G. Baskaran, Theory of Emergent Josephson Lattice in Neutral Twisted Bilayer Graphene (Moire is Different), https://arxiv.org/abs/1804.00627arXiv:1804.00627 (2018).

teor_tbg5  J. F. Dodaro, S. A. Kivelson, Y. Schattner, X. Q. Sun,
and C. Wang, Phases of a phenomenological model of twisted bilayer graphene, https://journals.aps.org/prb/abstract/10.1103/PhysRevB.98.075154Phys. Rev. B 98, 075154 (2018).


teor_tbg7  J. Kang and O. Vafek, Symmetry, maximally localized Wannier states, and low energy model for the twisted bilayer graphene narrow bands, https://journals.aps.org/prx/abstract/10.1103/PhysRevX.8.031088Phys. Rev. X 8, 031088 (2018).


teor_tbg9  H. Isobe, N. F. Q. Yuan, and L. Fu, Superconductivity and Charge Density Wave in Twisted Bilayer Graphene, https://journals.aps.org/prx/abstract/10.1103/PhysRevX.8.041041Phys. Rev. X 8, 041041
(2018).


teor_tbg11  Y.-Z. You and A. Vishwanath, Superconductivity from Valley Fluctuations and Approximate SO(4) Symmetry in a Weak Coupling Theory of Twisted Bilayer Graphene, http://arxiv.org/abs/1805.06867arXiv:1805.06867 (2018).

teor_tbg12  F. Wu, A. H. MacDonald, and I. Martin, Theory of phonon-mediated superconductivity in twisted bilayer graphene, https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.121.257001Phys. Rev. Lett. 121, 257001 (2018).

teor_tbg13  D. M. Kennes, J. Lischner, and C. Karrasch, Strong Correlations and d+id Superconductivity in Twisted Bilayer Graphene, https://journals.aps.org/prb/abstract/10.1103/PhysRevB.98.241407 Phys. Rev. B 98, 241407 (2018).

teor_tbg14  C.-C. Liu, L.-D. Zhang, W.-Q. Chen, and F. Yang, Chiral SDW and d + id superconductivity in the magic-angle twisted bilayer-graphene, https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.121.217001Phys. Rev. Lett. 121, 217001 (2018).

teor_tbg15  J. M. Murray and O. Vafek, Excitonic and superconducting orders from repulsive interaction on the doped honeycomb bilayer,  https://link.aps.org/doi/10.1103/PhysRevB.89.205119 Phys. Rev. B 89, 205119 (2014).

teor_tbg16  F. Guinea and B. Uchoa, Odd-momentum pairing and superconductivity in vertical graphene heterostructures, https://link.aps.org/doi/10.1103/PhysRevB.86.134521  Phys. Rev. B 86, 134521 (2012).

teor_tbg17  Y.-H. Zhang, D. Mao, Y. Cao, P. Jarillo-Herrero, and T. Senthil, Moire Superlattice with Nearly Flat Chern Bands: Platform for (Fractional) Quantum Anomalous Hall Effects and Unconventional Superconductivity, http://arxiv.org/abs/1805.08232 arXiv: 1805.08232 (2018) .

teor_tbg18  M. Fidrysiak, M. Zegrodnik, and J. Spaek, Unconventional topological superconductivity and phase diagram for a two-orbital model of twisted bilayer graphene, https://journals.aps.org/prb/abstract/10.1103/PhysRevB.98.085436 Phys. Rev. B 98, 085436 (2018).

teor_tbg19 B. Lian, Z. Wang, B. Andrei Bernevig, Twisted Bilayer Graphene: A Phonon Driven Superconductor, https://arxiv.org/abs/1807.04382arXiv:1807.04382 (2018).

teor_tbg20 E. Laksono, J. N. Leaw, A. Reaves, M. Singh, X. Wang, S. Adam, X. Gu, Singlet superconductivity enhanced by charge order in nested twisted bilayer graphene Fermi surfaces, https://www.sciencedirect.com/science/article/pii/S0038109818304903Solid State Commun. 282, 38 (2018).

teor_tbg21 Y. Su and S.-Z. Lin, Pairing symmetry and spontaneous vortex-antivortex lattice in superconducting twisted-bilayer graphene: Bogoliubov-de Gennes approach, https://journals.aps.org/prb/abstract/10.1103/PhysRevB.98.195101Phys. Rev. B 98, 195101 (2018).

YWChoi Y. W. Choi, H. J. Choi, Electron-phonon interaction in magic-angle twisted bilayer graphene, https://journals.aps.org/prb/abstract/10.1103/PhysRevB.98.241412Phys. Rev. B 98, 241412 (2018)

Uchoa  B. Uchoa and A. H. Castro Neto, Superconducting States of Pure and Doped Graphene, https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.98.146801Phys. Rev. Lett. 98, 146801 (2007).

BP2019 M. Alidoust, M. Willatzen, and A.-P. Jauho, Control of superconducting pairing symmetries in monolayer black phosphorus, https://journals.aps.org/prb/abstract/10.1103/PhysRevB.99.125417Phys. Rev. B 99, 125417 (2019).

Hosseini1  M.V. Hosseini, M. Zareyan, Unconventional superconducting states of interlayer pairing in bilayer and trilayer graphene, https://journals.aps.org/prb/abstract/10.1103/PhysRevB.86.214503Phys. Rev. B 86, 214503 (2012).

Hosseini2 M.V. Hosseini, M. Zareyan, Model of an exotic chiral superconducting phase in a graphene bilayer, https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.108.147001Phys. Rev. Lett. 108, 147001 (2012).

Hesselmann S. Hesselmann, D. D. Scherer, M. M. Scherer, and S. Wesse, Bond-ordered states and f-wave pairing of spinless fermions on the honeycomb lattice, https://journals.aps.org/prb/abstract/10.1103/PhysRevB.98.045142Phys. Rev. B 98, 045142
(2018).

RMP-2008-Beenakker C. W. J. Beenakker, Andreev reflection and Klein tunneling in graphene, https://journals.aps.org/rmp/abstract/10.1103/RevModPhys.80.1337Rev. Mod. Phys. 80, 1337 (2008).

RMP-2009-Neto A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, The electronic properties of graphene, https://journals.aps.org/rmp/abstract/10.1103/RevModPhys.81.109Rev. Mod. Phys. 81, 109 (2009).

PRL-2013-Halter K. Halterman, O. T. Valls, and M. Alidoust, Spin-Controlled Superconductivity and Tunable Triplet Correlations in Graphene Nanostructures, https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.111.046602Phys. Rev. Lett. 111, 046602 (2013).

arxiv_nonlinearMSE S. Bae, Y. Tan, A. P. Zhuravel, L. Zhang, S. Zeng, Y. Liu, T. A. Lograsso, R. Prozorov, A. T. Venkatesan, S. M. Anlage, Dielectric Resonator Method For Determining Gap Symmetry Of Superconductors Through Anisotropic Nonlinear Meissner Effect, https://arxiv.org/abs/1901.08762arXiv:1901.08762.

tbg15  H. C. Po, L. Zou, T. Senthil, A. Vishwanath, Faithful Tight-binding Models and Fragile Topology of Magic-angle Bilayer Graphene,  https://arxiv.org/abs/1808.02482arXiv:1808.02482 (2018).

tbg16  L. Huder, A. Artaud, T. L. Quang, G. T. Laissardiere, A.M. Jansen, G. Lapertot, C. Chapelier, and V. T. Renard, Electronic Spectrum of Twisted Graphene Layers under Heterostrain, https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.120.156405 Phys. Rev. Lett. 120, 156405 (2018).

tbg17  H. Patel, L. Huang, C.-J. Kim, J. Park, M. W. Graham, θ-Tunable Photoluminescence from Interlayer Excitons in Twisted Bilayer Graphene,  https://arxiv.org/abs/1806.06312 	arXiv:1806.06312 (2018).


tbg18  M. P. Zaletel, S. Geraedts, Z. Papic, and E. H. Rezayi, Evidence for a topological exciton Fermi sea in bilayer graphene,  https://journals.aps.org/prb/abstract/10.1103/PhysRevB.98.045113 Phys. Rev. B 98, 045113 (2018).

tbg19  Z. Song, Z. Wang, W. Shi, G. Li, C. Fang, B. A. Bernevig, All "Magic Angles" Are "Stable" Topological, https://arxiv.org/abs/1807.10676 arXiv:1807.10676 (2018).

teor_tbg8  L. Rademaker and P. Mellado, Charge-transfer insulation in twisted bilayer graphene, https://arxiv.org/abs/1805.05294arXiv:1805.05294 (2018).

teor_tbg10  M. Koshino, N. F. Q. Yuan, M. Ochi, K. Kuroki, and
L. Fu, Maximally-localized Wannier orbitals and the extended Hubbard model for the twisted bilayer graphene, https://journals.aps.org/prx/abstract/10.1103/PhysRevX.8.031087Phys. Rev. X
8, 031087 (2018).

teor_tbg11 B. Padhi, P. Phillips, Pressure-Induced Metal-Insulator Transition in Twisted Bi-layer Graphene, https://arxiv.org/abs/1810.00884arXiv:1810.00884.

teor_tbg12 B. Padhi, C. Setty, and P. W. Phillips, Doped twisted bilayer graphene near magic angles: Proximity to Wigner crystallization, not Mott insulation, https://doi.org/10.1021/acs.nanolett.8b02033Nano Lett. 18, 6175 (2018).

teor_tbg13 A. V. Rozhkov, A. O. Sboychakov, A. L. Rakhmanov, and F. Nori, Electronic properties of graphene-based bilayer systems, https://doi.org/10.1016/j.physrep.2016.07.003Phys. Rep. 648, 1 (2016).

teor_tbg14 A. V. Rozhkov, A. O. Sboychakov, A. L. Rakhmanov, and F. Nori, Single-electron gap in the spectrum of twisted bilayer graphene, https://doi.org/10.1103/PhysRevB.95.045119Phys. Rev. B 95, 045119 (2017).

mineev V.P. Mineev, K Samokhin, Introduction to Unconventional Superconductivity, CRC Press, (1999).

BP1-prb2018 M. Alidoust, M. Willatzen, A.-P. Jauho, Strain-engineered Majorana zero energy modes and φ_0 Josephson state in black phosphorus, https://link.aps.org/doi/10.1103/PhysRevB.98.085414Phys. Rev. B 98, 085414 (2018).

BP2-prb2018 M. Alidoust, M. Willatzen, A.-P. Jauho, Fraunhofer response and supercurrent spin switching in black phosphorus with strain and disorder, https://journals.aps.org/prb/abstract/10.1103/PhysRevB.98.184505Phys. Rev. B 98, 184505 (2018).

mm P. Marconcini and M. Macucci, The k.p method and its application to graphene, carbon nanotubes and graphene nanoribbons: the Dirac equation
, https://www.sif.it/riviste/sif/ncr/econtents/2011/034/08-09/article/0Rivista del Nuovo Cimento 34, 489 (2011).















