Symmetry of superconducting correlations in displaced bilayers of graphene

TL;DR

This paper investigates phonon-mediated superconducting pairing symmetries in displaced bilayer graphene, revealing various possible pairing symmetries and the potential for two-gap superconductivity depending on stacking order and doping.

Contribution

It introduces a comprehensive Green's function approach to analyze superconducting symmetries in displaced bilayer graphene, considering generic coupling potentials and stacking configurations.

Findings

Dominant pairing symmetries include s-wave and p-wave at charge neutrality.

d-wave symmetry can develop at finite doping in AB stacking.

Two-gap superconductivity is possible due to intralayer and interlayer phonon interactions.

Abstract

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.

Tables

Table 1. Table 1: Sublattice positions of BLG for different stacking orders.

Stacking$\to$	AA	AB	AC
${𝝂}_{1,A}=$	(0,0)	(0,0)	(0,0)
${𝝂}_{1,B}=$	$2\left({𝒂}_1+{𝒂}_2\right)/3$	$2\left({𝒂}_1+{𝒂}_2\right)/3$	$2\left({𝒂}_1+{𝒂}_2\right)/3$
${𝝂}_{2,A}=$	(0,0)	$\left({𝒂}_1+{𝒂}_2\right)/3$	$2\left({𝒂}_1+{𝒂}_2\right)/3$
${𝝂}_{2,B}=$	$2\left({𝒂}_1+{𝒂}_2\right)/3$	(0,0)	$\left({𝒂}_1+{𝒂}_2\right)/3$

Table 2. Table 2: Different combinations of momentum functions, appearing in the calculations.

$k_{+}=k_x+i{k}_y$	$k_{+}^2=k_x^2+2i{k}_x{k}_y-k_y^2$	$k_{+}^3=k_x^3+3i{k}_x^2{k}_y-3k_x{k}_y^2-i{k}_y^3$	$k_{+}^4=k_x^4+4i{k}_x^3{k}_y-6k_x^2{k}_y^2-4i{k}_x{k}_y^3+k_y^4$
$k_{-}=k_x-i{k}_y$	$k_{-}^2=k_x^2-2i{k}_x{k}_y-k_y^2$	$k_{-}^3=k_x^3-3i{k}_x^2{k}_y-3k_x{k}_y^2+i{k}_y^3$	$k_{-}^4=k_x^4-4i{k}_x^3{k}_y-6k_x^2{k}_y^2+4i{k}_x{k}_y^3+k_y^4$
$k_{+}{k}_{-}=k_x^2+k_y^2$	$k_{+}^2{k}_{-}^2={\left(k_x^2+k_y^2\right)}^2$	$k_{+}^3{k}_{-}^3={\left(k_x^2+k_y^2\right)}^3$	$k_{+}^4{k}_{-}^4={\left(k_x^2+k_y^2\right)}^4$
$k_{-}^2{k}_{+}=k_x^3-i{k}_x^2{k}_y+k_x{k}_y^2-i{k}_y^3$		$k_{-}^2{k}_{+}^3=k_x^5+i{k}_x^4{k}_y+2k_x^3{k}_y^2+2i{k}_x^2{k}_y^3+k_x{k}_y^4+i{k}_y^5$
$k_{-}^3{k}_{+}=k_x^4-2i{k}_x^3{k}_y-2i{k}_x{k}_y^3-k_y^4$		$k_{+}+k_{-}=2k_x$
$k_{-}^3{k}_{+}^2=k_x^5-i{k}_x^4{k}_y+2k_x^3{k}_y^2-2i{k}_x^2{k}_y^3+k_x{k}_y^4-i{k}_y^5$		$k_{+}-k_{-}=2i{k}_y$

Table 3. Table 3: Phonon mediated electron-electron interaction potentials.

$V_s\propto 1$	$V_{\text{extended}-s}\propto k_x^2+k_y^2$
$V_{p_x}\propto k_x$
$V_{d_{xy}}\propto k_x{k}_y$	$V_{d_{x^2-y^2}}\propto k_x^2-k_y^2$
$V_{f_{y\left(3x^2-y^2\right)}}\propto k_y\left(3k_x^2-k_y^2\right)$
$V_{g_{xy\left(x^2-y^2\right)}}\propto$	$V_{g_{x^4+y^4}}\propto k_x^2\left(k_x^2-3k_y^2\right)$
$k_x{k}_y\left(k_x^2-k_y^2\right)$	$-k_y^2\left(3k_x^2-k_y^2\right)$

Equations

$$H$$

$${\text{K}}_j={\text{K}}_0+{\text{G}}_j,$$

$$\tilde{\mathbf{T}}(\text{k})=\sum _{j=0,1,2}{\text{M}}_j^{\text{XX}}\frac{t_{\perp }({\text{K}}_j+\text{k})}{3}{e}^{i({\text{K}}_j+\text{k})({\text{u}}_2-{\text{u}}_1)},$$

$$\displaystyle\text{M}_{j}^{\text{AA}}=\left(\begin{array}[]{cc}1&e^{-i(2j\pi/3)}\\ e^{+i(2j\pi/3)}&1\end{array}\right),\;\;\;\text{M}_{j}^{\text{AB}}=\left(\begin{array}[]{cc}e^{+i(2j\pi/3)}&1\\ e^{-i(2j\pi/3)}&e^{+i(2j\pi/3)}\end{array}\right),$$

$$\displaystyle\text{M}_{j}^{\text{AC}}=\left(\begin{array}[]{cc}e^{-i(2j\pi/3)}&e^{+i(2j\pi/3)}\\ 1&e^{-i(2j\pi/3)}\end{array}\right).$$

$$i)$$

$$ii)$$

$$g_{\rho {\rho }^{\prime }}^{\sigma {\sigma }^{\prime }}(\tau ,{\tau }^{\prime };\mathbf{r},{\mathbf{r}}^{\prime })=-⟨{\mathcal{T}}_{\tau }{\psi }_{\rho \sigma }(\tau ,{\mathbf{r}}^{\prime }){\psi }_{{\rho }^{\prime }{\sigma }^{\prime }}^{†}({\tau }^{\prime },{\mathbf{r}}^{\prime })⟩,$$

$${\underline{g}}_{\rho {\rho }^{\prime }}^{\sigma {\sigma }^{\prime }}(\tau ,{\tau }^{\prime };\mathbf{r},{\mathbf{r}}^{\prime })=-⟨{\mathcal{T}}_{\tau }{\psi }_{\rho \sigma }^{†}(\tau ,\mathbf{r}){\psi }_{{\rho }^{\prime }{\sigma }^{\prime }}({\tau }^{\prime },{\mathbf{r}}^{\prime })⟩,$$

$$f_{\rho {\rho }^{\prime }}^{\sigma {\sigma }^{\prime }}(\tau ,{\tau }^{\prime };\mathbf{r},{\mathbf{r}}^{\prime })=+⟨{\mathcal{T}}_{\tau }{\psi }_{\rho \sigma }(\tau ,\mathbf{r}){\psi }_{{\rho }^{\prime }{\sigma }^{\prime }}({\tau }^{\prime },{\mathbf{r}}^{\prime })⟩,$$

$$f_{\rho {\rho }^{\prime }}^{\sigma {\sigma }^{\prime }†}(\tau ,{\tau }^{\prime };\mathbf{r},{\mathbf{r}}^{\prime })=+⟨{\mathcal{T}}_{\tau }{\psi }_{\rho \sigma }^{†}(\tau ,\mathbf{r}){\psi }_{{\rho }^{\prime }{\sigma }^{\prime }}^{†}({\tau }^{\prime },{\mathbf{r}}^{\prime })⟩,$$

$$\displaystyle\left(\begin{array}[]{cc}\hat{H}(\mathbf{r})-i\omega_{n}&\hat{\Delta}(\textbf{r})\\ \hat{\Delta}^{*}(\textbf{r})&{\cal T}_{\text{t}}\hat{H}(\mathbf{r}){\cal T}_{\text{t}}^{\dagger}+i\omega_{n}\end{array}\right)\check{g}(i\omega_{n};\mathbf{r},\mathbf{r}^{\prime})=\delta(\mathbf{r}-\mathbf{r}^{\prime}),$$

$$\check{g}(i\omega_{n};\mathbf{r},\mathbf{r}^{\prime})=\left(\begin{array}[]{cc}\hat{g}(i\omega_{n};\mathbf{r},\mathbf{r}^{\prime})&\hat{f}(i\omega_{n};\mathbf{r},\mathbf{r}^{\prime})\\ \hat{f}^{\dagger}(i\omega_{n};\mathbf{r},\mathbf{r}^{\prime})&\hat{\underline{g}}(i\omega_{n};\mathbf{r},\mathbf{r}^{\prime})\end{array}\right).$$

$$\Omega f_{11}^{\uparrow \uparrow }=0,$$

$$\Omega f_{11}^{\uparrow \downarrow }=-{\Delta }_{\text{S}}(k_{-}{k}_{+}+t^2+{\Delta }_{\text{S}}^2+{\omega }_n^2),$$

$$\Omega f_{12}^{\uparrow \uparrow }=-2t{\Delta }_{\text{S}}{k}_{-},$$

$$\Omega f_{12}^{\uparrow \downarrow }=0,$$

$$\Omega f_{11}^{\downarrow \uparrow }=+{\Delta }_{\text{S}}(k_{-}{k}_{+}+t^2+{\Delta }_{\text{S}}^2+{\omega }_n^2),$$

$$\Omega f_{11}^{\downarrow \downarrow }=0,$$

$$\Omega f_{12}^{\downarrow \uparrow }=0,$$

$$\Omega f_{12}^{\downarrow \downarrow }=+2t{\Delta }_{\text{S}}{k}_{+}$$

$$\Omega =k_{-}^2{k}_{+}^2+2k_{-}{k}_{+}({\Delta }_{\text{S}}^2-t^2+{\omega }_n^2)+(t^2+{\Delta }_{\text{S}}^2+{\omega }_n^2{)}^2.$$

$$\Omega f_{11}^{\uparrow \uparrow }=0,$$

$$\Omega f_{11}^{\uparrow \downarrow }=-{\Delta }_{\text{S}}(k_{-}{k}_{+}+{\Delta }_{\text{S}}^2+{\omega }_n^2),$$

$$\Omega f_{12}^{\uparrow \uparrow }=-it{\Delta }_{\text{S}}{\omega }_n,$$

$$\Omega f_{12}^{\uparrow \downarrow }=0,$$

$$\Omega f_{11}^{\downarrow \uparrow }=+{\Delta }_{\text{S}}(k_{-}{k}_{+}+{\Delta }_{\text{S}}^2+{\omega }_n^2),$$

$$\Omega f_{11}^{\downarrow \downarrow }=0,$$

$$\Omega f_{12}^{\downarrow \uparrow }=t{\Delta }_{\text{S}}(k_{-}-k_{+}),$$

$$\Omega f_{12}^{\downarrow \downarrow }=+it{\Delta }_{\text{S}}{\omega }_n,$$

$$\Omega =(k_{-}{k}_{+}+{\Delta }_{\text{S}}^2{)}^2+(2k_{-}{k}_{+}+t^2+2{\Delta }_{\text{S}}^2){\omega }_n^2+{\omega }_n^4.$$

$$\Omega f_{11}^{\uparrow \uparrow }=0,$$

$$\Omega f_{11}^{\uparrow \downarrow }=-2it{\Delta }_{\text{B}}{\omega }_n,$$

$$\Omega f_{12}^{\uparrow \uparrow }=0,$$

$$\Omega f_{12}^{\uparrow \downarrow }=-{\Delta }_{\text{B}}(k_{-}{k}_{+}-t^2+{\Delta }_{\text{B}}^2+{\omega }_n^2),$$

$$\Omega f_{11}^{\downarrow \uparrow }=2it{\Delta }_{\text{B}}{\omega }_n,$$

$$\Omega f_{11}^{\downarrow \downarrow }=0,$$

$$\Omega f_{12}^{\downarrow \uparrow }={\Delta }_{\text{B}}(k_{-}{k}_{+}-t^2+{\Delta }_{\text{B}}^2+{\omega }_n^2),$$

$$\Omega f_{12}^{\downarrow \downarrow }=0,$$

$$\Omega =(k_{-}{k}_{+}-t^2+{\Delta }_{\text{B}}^2{)}^2+2(k_{-}{k}_{+}+t^2+{\Delta }_{\text{B}}^2){\omega }_n^2+{\omega }_n^4.$$

$$\Omega f_{11}^{\uparrow \uparrow }=0,$$

$$\Omega f_{11}^{\uparrow \downarrow }=-t{\Delta }_{\text{B}}{k}_{+},$$

$$\Omega f_{12}^{\uparrow \uparrow }=0,$$

$$\Omega f_{12}^{\uparrow \downarrow }=-{\Delta }_{\text{B}}(k_{-}{k}_{+}+{\Delta }_{\text{B}}^2+{\omega }_n^2),$$

$$\Omega f_{11}^{\downarrow \uparrow }=t{\Delta }_{\text{B}}{k}_{+},$$

$$\Omega f_{11}^{\downarrow \downarrow }=0,$$

$$\Omega f_{12}^{\downarrow \uparrow }={\Delta }_{\text{B}}(k_{-}{k}_{+}+t^2+{\Delta }_{\text{B}}^2+{\omega }_n^2),$$

$$\Omega f_{12}^{\downarrow \downarrow }=0,$$

Full text

I 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 Butz ; Kisslinger ; Hass ; Borysiuk ; Liu .

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 tbg-2 ; tbg-1 ; tbg0 ; tbg1 ; tbg2 ; tbg3 ; tbg4 ; tbg5 ; tbg6 ; tbg8 ; tbg9 ; tbg10 ; tbg11 ; tbg12 ; tbg13 ; tbg14 ; 2ph-el_exp . 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 ray1 ; ray2 ; ray_slip .

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. tbg2 ; tbg3 ; tbg4 ; tbg5 ; tbg6 ; tbg8 ; tbg9 ; tbg10 ; tbg11 ; tbg12 ; tbg13 ; tbg14 ; 2ph-el_exp . 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 2ph-el_exp .

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 $\sim 1.7$ K in a twisted BLG, away from the charge neutrality point $\mu>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 herrero . Soon after, another experiment C.Dean_exp (and recently few others grexp1 ; grexp2 ; grexp3 ) 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 Uchoa ; Hosseini1 ; Hosseini2 ; Hesselmann .

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 2ph-el_exp . 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 $\mu=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: $\Delta$. 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. II, 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 A). In Sec. III, 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. III.1. Next, for a doped system, we evaluate the Green’s functions numerically, and support these results by analytic expressions given in Appendix B. In Sec. IV, we consider a generic coupling potential between the two pristine graphene layers and repeat our studies of Sec. II, 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. V.

II Low energy effective Hamiltonian and green’s function approach

We follow the effective Hamiltonian strategy ray1 ; ray2 ; ray_slip 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

[TABLE]

in which the layer degree of freedom is described in terms of Pauli matrices $\rho_{0,z,x,y}$, and we have defined $2\rho_{1}=\rho_{0}+\rho_{z}$, $2\rho_{2}=\rho_{0}-\rho_{z}$, $2\rho_{+}=\rho_{x}+i\rho_{y}$, $2\rho_{-}=\rho_{x}-i\rho_{y}$. The matrix $\tilde{\bm{T}}(\textbf{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) ray1 ; ray2 ; ray_slip . The Hamiltonians of the top and bottom layers are denoted by $H_{1,2}({\bf k})=\hbar v_{\textbf{F}}{\bf k}\cdot{\bm{\sigma}}$, respectively, in which ${\bm{\sigma}}=(\sigma_{x},\sigma_{y})$ and $\sigma_{x},\sigma_{y}$ are Pauli matrices. The multiplication of Pauli matrices $\rho_{i}$ and $\sigma_{i}$ implies a tensor product so that the Hamiltonian (II) is a $4\times 4$ matrix and the matrix $\tilde{\bm{T}}(\textbf{k})$ has an implicit $\bm{\sigma}$ (refer to Eqs. (3) and (8) 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 $\hat{\psi}^{\dagger}(\textbf{k})=(\psi_{1\uparrow}^{\dagger},\psi_{1\downarrow}^{\dagger},\psi_{2\uparrow}^{\dagger},\psi_{2\downarrow}^{\dagger})$ where the first index is 1 (2) for the top (bottom) SLG and the second index is $\uparrow$ ($\downarrow$) 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 RMP-2008-Beenakker ; RMP-2009-Neto ; PRL-2013-Halter . In this case, one can show that the system possesses spin degeneracy RMP-2008-Beenakker ; RMP-2009-Neto ; PRL-2013-Halter . 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 A.

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 ${\bm{\nu}}_{A}=(0,0),{\bm{\nu}}_{B}=2({\bm{a}}_{1}+{\bm{a}}_{2})/3$ and two lattice vectors ${\bm{a}}_{1}=a(1,0),{\bm{a}}_{2}=a(1,\sqrt{3})/2$. Using this notation, the different stackings can be described by the sublattice positions given in Table 1. Figure 1(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. 1(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 1).

The reciprocal lattice of a BLG can be described by the lattice vectors ${\bm{b}}_{1}=2a^{-1}\pi(1,-1/\sqrt{3}),{\bm{b}}_{2}=2a^{-1}\pi(0,2/\sqrt{3})$ in reciprocal space as shown in Fig. 1(c). Since we are interested in the low-energy physics, the Hamiltonian can be expanded around reciprocal sublattice points:

[TABLE]

in which $\textbf{K}_{0}=(2{\bm{b}}_{1}+{\bm{b}}_{2})/3,\textbf{G}_{0}={\bm{0}},\textbf{G}_{1}=-{\bm{b}}_{1},\textbf{G}_{2}=-{\bm{b}}_{1}-{\bm{b}}_{2}$ (see Fig. 1(b)). The periodic arrangement of atoms in SLG is spanned by $\mathbf{r}=\textbf{R}_{j}+{\bm{\nu}}_{\alpha}$, where $\textbf{R}_{j}$ are integer combinations of the primitive lattice vectors ${\bm{a}}_{j}$, and the orbitals centered at positions ${\bm{\nu}}_{\alpha}$. In the following, $\mathbf{r}$ 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. 1(b). The displacement enters the coupling matrix $\tilde{\bm{T}}$ via a phase factor ray1 ; ray2 ; ray_slip ,

[TABLE]

in which $t_{\perp}(\textbf{K}_{j}+\textbf{k})$ is the interlayer hopping amplitude, given by $t_{\perp}({\bf q})=V_{u}^{-1}\int d{\bf r}t_{\perp}({\bf r})e^{i{\bf q}{\bf r}}$, t_{\perp}({\bf r})=\sum_{i,j}t_{2,j}^{1,i}\Big{\langle}{\bf r},2\Big{|}c^{\dagger}_{2,j}c_{1,i}\Big{|}{\bf r}^{\prime},1\Big{\rangle} with $c^{\dagger},c$ being the creation and annihilation operators and ${\bf r}^{\prime}={\bf r}+{\bm{\delta}}$ in which ${\bm{\delta}}$ denotes the distance between hopping sites. $V_{u}$ is the volume of unit cell, $i,j$ run over the lattice sites, $\textbf{u}_{1,2}$ are the displacements of the top and bottom layers, and $\text{M}_{j}^{\text{XX}}$ defines stacking order:

[TABLE]

In what follows, we consider two types of couplings that may lead to phonon-mediated superconductivity in a BLG 2ph-el_exp . 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 2ph-el_exp . 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:

[TABLE]

where $\Delta_{\text{S}}$ and $\Delta_{\text{B}}$ are the BCS spin-singlet phonon-mediated electron-electron coupling within each layer and between the layers, respectively. For more details see Appendix A.

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; BP2019

[TABLE]

where ${\cal T}_{\tau}$ is the time ordering operator, and $\tau,\tau^{\prime}$ are the imaginary times. Here, $\sigma,\sigma^{\prime}$ and $\rho,\rho^{\prime}$ denote the Pspin and layer indices and $\langle...\rangle$ is the thermodynamics averaging. In particle-hole space, the Green’s function satisfies:

[TABLE]

in which $\omega_{n}=\pi(2n+1)k_{B}T$ is the Matsubara frequency, $n\in{Z}$, $k_{B}$ is the Boltzman constant, $T$ is temperature, ${\cal T}_{\text{t}}$ is the time-reversal operator, and $\hat{\Delta}(\mathbf{r})$ is the superconducting gap. Note that $\hat{H}(\mathbf{r})$ is obtained by replacing ${\bf k}\rightarrow-i\bm{\nabla}$ in $H(\textbf{k})$, in Eq. (II). The matrix form of the Green’s function is given by;

[TABLE]

We denote $4\times 4$ matrices by ‘hat’ symbol, $\hat{\square}$, and $8\times 8$ matrices by ‘check’ symbol, $\check{\square}$.

We now proceed to derive the anomalous Green’s function $f$ for the different coupling scenarios (Eq. (26) below gives the self-consistency condition these functions must satisfy).BP2019

III 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., $\Delta_{\text{S}}^{1}\equiv\Delta_{\text{S}}^{2}$. In what follows, we use a compact notation for the different combinations of momenta (see Table 2) 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}_{\rho\rho^{\prime}}^{\sigma\sigma^{\prime}}$ correspond to pairing correlations that may arise in the system whereas the normal Green’s function ${g}_{\rho\rho^{\prime}}^{\sigma\sigma^{\prime}}$ contains information about local density of states. Hence, in what follows, we focus on the anomalous Green’s function ${f}_{\rho\rho^{\prime}}^{\sigma\sigma^{\prime}}$ and do not present explicitly the normal Green’s function ${g}_{\rho\rho^{\prime}}^{\sigma\sigma^{\prime}}$ 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_{\pm}$(=$k_{x}\pm ik_{y}$) $\rightarrow$ $p$-wave, $k_{\pm}^{2}$ $\rightarrow$ $d$-wave, and $k_{\pm}^{2}k_{\mp}$ $\rightarrow$ $f$-wave, etc. (see also Tables 2 and 3).

In Secs. III.1 and III.2, we consider undoped $\mu=0$ and doped $\mu\neq 0$ BLG, respectively. Each of these sections studies both AA and AB orderings with both intralayer $\Delta_{\text{S}}$ and interlayer $\Delta_{\text{B}}$ phonon-mediated electron-electron interaction.

III.1 Undoped bilayer graphene: $\mu=0$

We consider two cases separately: spin singlet intralayer phonon-mediated electron-electron coupling $\Delta_{\text{S}}$ and spin-singlet interlayer phonon-mediated electron-electron coupling $\Delta_{\text{B}}$. In both cases, we assume a constant coupling potential between the two monolayers given by: $t=t^{\dagger}$. The components of coupling potential are defined by Eq. (3).

III.1.1 Intralayer coupling: $\Delta_{\text{S}}\neq 0,\Delta_{\text{B}}=0$

$\bullet$AA stacking: In this case, the components of the anomalous Green’s function Eqs. (13) are:

[TABLE]

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.

$\bullet$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:

[TABLE]

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$.

III.1.2 Interlayer BCS coupling: $\Delta_{\text{S}}=0,\Delta_{\text{B}}\neq 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. (12b). For AA and AB stack orderings of the two monolayer graphene, we find the following anomalous components of the Green’s function, Eqs. (13);

$\bullet$AA stacking:

[TABLE]

$\bullet$AB stacking:

[TABLE]

Comparing with the Green’s function components of the same stacking orders but with intralayer phonon-mediated electron-electron coupling discussed in Subsec. III.1.1, 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. III.1.1. The same symmetry correspondence can be seen between AB stacking of this current subsection and Sec. III.1.1. In that case, BLG develops odd-parity even-frequency correlations. Note that, unlike the intralayer phonon-mediated electron-electron coupling scenario, Sec. III.1.1, 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_{\pm}$. See Tables 2 and 3). 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.

III.2 Bilayer graphene with finite doping $\mu\neq 0$

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

III.2.1 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 B 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\sim 0.01T_{c}$). We also set a ‘representative’ ratio to the chemical potential over coupling strength $\mu/t=1.5$ throughout the numerical study without loss of generality. Figure 2 illustrates the momentum space map of the real and imaginary parts of the Green’s function components in AA stacking order. In Fig. 2(a), the phonon-mediated electron-electron coupling is of the intralayer type $\Delta_{\text{S}}\neq 0,\Delta_{\text{B}}=0$ while in Fig. 2(b) we consider only interlayer coupling $\Delta_{\text{S}}=0,\Delta_{\text{B}}\neq 0$. For the intralayer coupling, Pspin-singlet $f_{11}^{\uparrow\downarrow},f_{12}^{\uparrow\downarrow},f_{11}^{\downarrow\uparrow},f_{12}^{\downarrow\uparrow}$ components are real and have $s$-wave ($\propto k_{x}^{2}+k_{y}^{2}$) symmetry while the other components with triplet Pspin state show $p$-wave ($\propto k_{\pm}$) and $f$-wave ($\propto k_{\pm}^{2}k_{\mp}$) 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 ($\propto k_{x}^{2}+k_{y}^{2}$) symmetry class. The interlayer Pspin-singlet correlations, $f_{12}^{\uparrow\downarrow}$ and $f_{12}^{\downarrow\uparrow}$, 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. (31) and (33) in Appendix B. According to the momentum combinations presented in Table 2, Eqs. (31) and (33) contain no term of type single $k_{\pm}^{2}$ that results in $d$-wave symmetry class. The allowed momentum combinations are either odd-parity or possess $s$-wave symmetry.

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

Figure 3 shows the pairing correlations in a BLG with AB ordering of the top and bottom pristine layers. In Fig. 3(a), the BCS electron-electron coupling is of the intralayer type, i.e., $\Delta_{\text{S}}\neq 0,\Delta_{\text{B}}=0$, while Fig. 3(b) contains results with interlayer BCS electron-electron coupling, i.e., $\Delta_{\text{S}}=0,\Delta_{\text{B}}\neq 0$. Clearly, in Fig. 3(a), components $f_{11}^{\uparrow\downarrow},f_{12}^{\uparrow\uparrow},f_{11}^{\downarrow\uparrow},f_{12}^{\downarrow\downarrow}$ 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_{\pm}^{2}$) symmetries can be confirmed by the analytic expressions in Eqs. (32). In Fig. 3(b), for the interlayer BCS coupling, the only component that possesses effective $d$-wave symmetry is $f_{11}^{\downarrow\downarrow}$. Checking through Eqs. (34), we see that $f_{11}^{\downarrow\downarrow}$ expression contains single $k_{+}^{2}$ which produces the $d$-wave symmetry. The components $f_{11}^{\uparrow\uparrow},f_{12}^{\uparrow\downarrow},f_{12}^{\downarrow\uparrow}$ are of effective $s$-wave type while $f_{11}^{\uparrow\downarrow},f_{12}^{\uparrow\uparrow},f_{11}^{\downarrow\uparrow},f_{12}^{\downarrow\downarrow}$ show $p$-wave symmetry beside terms proportional to $k_{+}k_{-}^{2}k_{+}^{2}$ and $k_{+}^{2}k_{-}$. Mathematically, the reason that $f_{11}^{\downarrow\downarrow}$ gains $d$-wave symmetry is the specific combination of the interlayer coupling $\tilde{\bm{T}}$ which has off-diagonal terms, and the phonon-mediated electron-electron coupling $\Delta_{\text{B}}$ that is on the off-diagonal entries of Eq. (16).

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 arxiv_nonlinearMSE .

IV 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 $\theta/2$, one clockwise and the other counter clockwise. Therefore, the uncoupled single layers can be described by $H_{1,2}^{\text{S}}(\pm\theta/2,\theta_{\text{k}})=\Theta_{z}(\pm\theta/2)H^{\text{S}}(\theta_{\text{k}})\Theta_{z}^{\dagger}(\pm\theta/2)$, where $\theta_{\text{k}}$ is the direction of the particles’ momentum in each layer and $\Theta_{z}(\theta)$ 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 $\tilde{\bm{T}}({\theta,\mathbf{r}},\textbf{k})$ from Eq. (3), now with the phase factor $\exp(i\textbf{g}\cdot{\mathbf{r}})$ in which $\textbf{g}_{j}=\Theta_{z}(+\theta/2)\textbf{G}_{j}-\Theta_{z}(-\theta/2)\textbf{G}_{j}$. For commensurate rotation angles, e.g., $\theta=21.787^{\circ},38.213^{\circ}$ one finds an effective $4\times 4$ Hamiltonian in the momentum space for the twisted BLG tbg1 ; tbg-2 . 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 tbg15 ; tbg16 ; tbg17 ; tbg18 ; tbg19 ; teor_tbg8 ; teor_tbg10 ; teor_tbg11 ; teor_tbg12 ; teor_tbg13 ; teor_tbg14 . Also, an AA-stacked bilayer can be transformed into an AB stacking order by keeping layer $1$ fixed, i.e. $\textbf{u}_{1}=0$, and shifting layer $2$ by $\textbf{u}_{2}=\epsilon({\bm{a}}_{1}+{\bm{a}}_{2})/3$ in Eq. (3). One can easily show that this transformation changes AA order into AB, AB into AC, AC into AA when $\epsilon$ is equal to unity.

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

[TABLE]

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, $\mu=0$, and only include the intralayer coupling: $\Delta_{\text{S}}\neq 0$, $\Delta_{\text{B}}=0$. The results for the anomalous Green’s function are:

[TABLE]

[TABLE]

An examination of these equations leads to the following conclusions. Similar to Sec. III.1, the generic coupling matrix leads to vanishing $f_{11}^{\uparrow\uparrow},f_{11}^{\downarrow\downarrow}$. However, the symmetries are no longer similar to those explored in Sec. III.1. As seen, here, $\Omega$ includes terms of single $k_{\pm}^{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 ($\Delta_{S}=0,\Delta_{B}\neq 0$ at $\mu=0$). The pairing correlations are given in Eqs. (25):

[TABLE]

We see that here, also $f_{11}^{\uparrow\uparrow},f_{11}^{\downarrow\downarrow}$ are zero. Similarly to the previous case, $\Omega$ includes single $k_{\pm}^{2}$ terms, carrying $d$-wave symmetry. Note that these single $k_{\pm}^{2}$ terms are not vanishing here because all the entries of coupling matrix $\tilde{\bm{T}}$ are assumed nonzero with no specific symmetry. In all the previous cases, with specific stacking orders, $\Omega$ includes powers of $k_{\pm}k_{\mp}$ 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 4 and 5 are shifted ($\epsilon=0.2$) counterparts of Figs. 2 and 3. When comparing these figures, we see that, for AA stacking, the shift has induced imaginary parts into $f_{11}^{\uparrow\downarrow},f_{11}^{\downarrow\uparrow},f_{12}^{\uparrow\downarrow},f_{12}^{\downarrow\uparrow}$ 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}^{\uparrow\downarrow},f_{12}^{\downarrow\uparrow}$ with $\Delta_{\text{S}}\neq 0$ and nonzero imaginary parts into $f_{11}^{\uparrow\uparrow},f_{12}^{\uparrow\downarrow}$ with $\Delta_{\text{B}}\neq 0$. In all cases, the displacement tends to first destroy the symmetries and, subsequently, at large values of $\epsilon$ 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 stateherrero ; C.Dean_exp , 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 $\Delta_{\text{S}}$ and $\Delta_{\text{B}}$ are nonzero. The results of this two-gap superconductivity scenario are given by Eqs. (35) of the Appendix. In deriving Eqs. (35) we consider bulk properties , i.e., $\Delta=\Delta^{\dagger}$ and assume that the two layers develop identical superconducting gap $\Delta_{\text{S}}$, which is unequal to the interlayer superconducting gap $\Delta_{\text{B}}$. To be specific, we consider twisted bilayers by an amount $\theta$ with zero doping, i.e., $\mu=0$. The intralayer Pspin-triplet components remains zero. Similarly to the single gap superconductivity case, $\Omega$ contains $k_{\pm}^{2}$ terms due to nonzero coupling potential terms considered.

IV.1 Temperature dependence of superconducting gaps

Recalling Eqs. (31)-(34), 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({\textbf{k},\textbf{k}^{\prime}})$ 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. (16):

[TABLE]

in which $\alpha,\alpha^{\prime},\beta,\beta^{\prime}$ 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 3.

The type and details of phonon-mediated electron-electron interactions can be obtained through optical absorption and Raman spectroscopy experiments tbg2 ; tbg3 ; tbg4 ; tbg5 ; tbg6 ; tbg8 ; tbg9 ; tbg10 ; tbg11 ; tbg12 ; tbg13 ; tbg14 ; 2ph-el_exp . With the relevant anomalous Green’s functions at hand, Eqs. (31)-(35), it is sufficient to insert the interaction potential into Eq. (26) 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 3) and plot the superconducting gap as a function of temperature for different scenarios presented in Fig. 6. Note that we have used an identical interaction potential for intralayer and interlayer coupling potential in AA and AB orderings and set $\mu/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., $\Delta_{\text{S}}^{\text{AA}}(T=0)$ and its critical temperature $T_{c,\text{S}}^{\text{AA}}$, respectively. Figure 6 illustrates that the interlayer phonon-mediated electron-electron coupling in AA ordering results in the smallest critical temperature (Fig. 6(a)) while its AB ordering counterpart Fig. 6(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., $\epsilon=0.2$, according to the discussion in the previous subsection. Compairing Figs. 6(a) and 6(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. 6(a), decreases the gap amplitude and critical temperature, it enhances these quantities in AB ordering, Fig. 6(c). This finding is reversed for the interlayer BCS coupling scenario in AA configuration. The real parts of $f_{11}^{\uparrow\downarrow},f_{11}^{\downarrow\uparrow}$ and $f_{12}^{\uparrow\downarrow},f_{12}^{\downarrow\uparrow}$ 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. 2 and 3 with Figs. 4 and 5. The displacement $\epsilon$ induces real parts to $f_{11}^{\uparrow\downarrow}$ and $f_{12}^{\uparrow\downarrow}$ with interlayer and intralyer phonon-mediated electron-electron coupling of AA and AB orderings, respectively, Figs. 6(b) and 6(c). However, a nonzero displacement $\epsilon$ applied to AA and AB orderings with intralayer and interlayer BCS coupling, respectively, induces nonzero imaginary parts to $f_{11}^{\uparrow\downarrow}$ and $f_{12}^{\uparrow\downarrow}$ that changes slightly the gap function and critical temperature, Figs. 6(a) and 6(d).

One can repeat the study presented in Fig. 6 for other interaction types given in Table 3 and the given Green’s functions in Eqs. (31)-(35). 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 3 for the electron-electron coupling potential $V({\textbf{k},\textbf{k}^{\prime}})$ with its spin state being singlet, the ratio of superconducting gap and critical temperature increases (according to the traditional calculations of such interaction potentials mineev ) 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 herrero 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. (22) 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 3, one should repeat the above calculations with these momentum-dependent potentials.

V 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, $\mu=0$, can only develop even-parity $s$-wave and odd-parity $p$-wave superconducting correlations. At a finite doping, $\mu\neq 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 $\mu=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.

Acknowledgements.

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).

Appendix A Hamiltonian of SLG with superconductivity

A monolayer graphene can be described by two valleys in $k$-space: K =$\frac{2\pi}{3a}(1,+1/\sqrt{3})$ and K′=$\frac{2\pi}{3a}(1,-1/\sqrt{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_{\bf K}(\textbf{k})=\hbar v_{\textbf{F}}\textbf{k}\cdot{\bm{\sigma}}$. The corresponding wavefunction is a two-component spinor, namely,

[TABLE]

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. 1(a), where sublattices are marked by blue and red colors. Hence, the total spin-less Hamiltonian can be expressed by:

[TABLE]

where $H_{{\bf K}^{\prime}}(\textbf{k})=\hbar v_{\textbf{F}}\textbf{k}\cdot{\bm{\sigma}}^{*}$. Next, we can add spin to this Hamiltonian by introducing Pauli matrices acting in spin-space: ${\bm{\tau}}=(\tau_{x},\tau_{y},\tau_{z})$. In this case, the spin-full Hamiltonian can be expressed by $H_{tot}=H(\textbf{k})\tau_{0}$ where $\tau_{0}$ reflects the absence of spin coupling. In SLG the two valleys are decoupled in the absence of an external potentialmm . Note that the form Eq. (28) applies even in the presence of elastic impurities and disorder. Therefore, the Hamiltonian becomes degenerated and reduces to a $4\times 4$ Hamiltonian. Now we incorporate superconductivity by introducing a two-electron amplitude,

[TABLE]

in which $\Delta_{\text{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 BP1-prb2018 ; BP2-prb2018 ; BP2019 . In the main text, we have used the above two-electron amplitude Eq. (29) 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:

[TABLE]

where $\Delta_{\text{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-degeneracyRMP-2008-Beenakker ; RMP-2009-Neto ; PRL-2013-Halter , 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 $\uparrow$,$\downarrow$ and call them pseudospin (Pspin).

Appendix B 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.

$\bullet$Intralayer BCS coupling, $\Delta_{\text{S}}\neq 0,\Delta_{\text{B}}=0$, in AA stacking order

[TABLE]

$\bullet$Intralayer BCS coupling, $\Delta_{\text{S}}\neq 0,\Delta_{\text{B}}=0$, in AB stacking order

[TABLE]

$\bullet$Interlayer BCS coupling, $\Delta_{\text{S}}=0,\Delta_{\text{B}}\neq 0$, in AA stacking order

[TABLE]

$\bullet$Interlayer BCS coupling, $\Delta_{\text{S}}=0,\Delta_{\text{B}}\neq 0$, in AB stacking order

[TABLE]

$\bullet$Intralayer and interlayer BCS coupling, $\Delta_{\text{S}}\neq 0,\Delta_{\text{B}}\neq 0$, with a generic coupling matrix

[TABLE]

Here we have defined $k_{\pm}^{+}=k\cos(\theta_{\text{k}}+\theta/2)\pm ik\sin(\theta_{\text{k}}+\theta/2)$ and $k_{\pm}^{-}=k\cos(\theta_{\text{k}}-\theta/2)\pm ik\sin(\theta_{\text{k}}-\theta/2)$ according to the main text for a twisted BLG by an amount $\theta$.