<figure><img src="image_1.png" /></figure>

HAL
open science

# Superconductivity in monolayer and few-layer graphene:
I. Review of possible pairing symmetries and basic electronic properties

Emile Pangburn, Louis Haurie, Adeline Crépieux, Oladunjoye A. Awoga, Annica M. Black-Schaffer, Catherine Pépin, Cristina Bena

<figure><img src="image_2.png" /></figure>

► To cite this version:

Emile Pangburn, Louis Haurie, Adeline Crépieux, Oladunjoye A. Awoga, Annica M. Black-Schaffer, et al.. Superconductivity in monolayer and few-layer graphene: I. Review of possible pairing symmetries and basic electronic properties. Physical Review B, 2023, 108, pp.134514. 10.1103/Phys-RevB.108.134514. cea-04259594

HAL Id: cea-04259594

https://cea.hal.science/cea-04259594

Submitted on 26 Oct 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Superconductivity in monolayer and few-layer graphene. I. Review of possible pairing symmetries and basic electronic properties

Emile Pangburn,¹ Louis Haurie,¹ Adeline Crépieux,² Oladunjoye A. Awoga,³ Annica M. Black-Schaffer,⁴ Catherine Pépin,¹ and Cristina Bena¹

¹Institut de Physique Théorique, Université Paris Saclay, CEA CNRS, Orme des Merisiers, 91190 Gif-sur-Yvette Cedex, France

²Aix Marseille Univ, Université de Toulon, CNRS, CPT, Marseille, France

³Solid State Physics and NanoLund, Lund University, Box 118, S-221 00 Lund, Sweden

⁴Department of Physics and Astronomy, Uppsala University, Box 516, S-751 20 Uppsala, Sweden

<figure><img src="image_4.png" /></figure>

(Received 19 January 2023; revised 5 April 2023; accepted 27 September 2023; published 25 October 2023)

We review all symmetry-allowed spin-singlet and spin-triplet superconducting order parameters in graphene (*s*-wave, *d*-wave, *p*-wave, and *f*-wave) generated by generic on-site, nearest-neighbor, and next-nearest-neighbor pairing interactions in a tight-binding model. For each pairing channel, we calculate both the band structure and the dependence of the density of states on energy, chemical potential, and pairing strength. In particular, we distinguish between nodal superconducting states and fully gapped states and study the dependence of gap closing points on the chemical potential and the superconducting pairing strength. We further investigate the difference between mono-, bi-, and trilayer ABC and ABA graphene, including accounting for the effects of trigonal warping.

DOI: 10.1103/PhysRevB.108.134514

## I. INTRODUCTION

The possibility to achieve experimentally isolated graphene layers [1] has kindled widespread interest in making graphene superconducting. Despite giving rise to exotic transport properties [2], the presence of Dirac cones impedes intrinsic superconductivity in undoped monolayer graphene due to the unavailability of electrons for pairing. However, several theoretical studies have shown that when the Fermi surface is far away from the Dirac point, up to and beyond the Van Hove singularity (VHS) where the density of states (DOS) diverges, monolayer graphene can host unconventional superconducting (SC) states, with symmetries ranging from *s* wave and *d* wave to *p* wave and *f* wave [3–15]. Recent studies have shown that doping graphene to such large levels around the VHS is possible [16,17]. Another possible direction to obtain superconductivity in graphene that received a lot of attention recently is to proximize graphene, with reports of observing in spectroscopy experiments both conventional SC features [18,19] as well as more unconventional ones [20–23].

On the other hand, recent experimental discoveries of superconductivity in rhombohedral or ABC-stacked trilayer graphene [24] and in twisted graphene bilayer systems [25] have received a lot of attention [26–57]. A vast range of theoretical proposals have already emerged exploring the mechanisms and symmetries of the SC state in these all-carbon systems [43], including the SC pairing mechanism being both phonon-mediated [28,50,58–60] and electron-interaction mediated [29,34] and various SC spin-singlet and spin-triplet order parameters, ranging from *s*-wave to higher angular momentum symmetries, from exotic chiral *p* + *ip'*-wave to chiral *d* + *id'*-wave symmetries [33,44,47–49,51–53,61–69]. However, at present, there exists no experimental

definite confirmation of a specific mechanism or pairing symmetry, nor an emerging consensus concerning these issues.

While currently achieved critical temperatures in graphene-based systems are only of the order of a few Kelvin [25], they host some tantalizing similarities to the high-temperature cuprate superconductors [70], in particular, a similar phase diagram with multiple regions reminiscent of strongly correlated electron physics, such as Mott insulating [71] and strange metal behavior [72,73]. The importance of strong electron correlations is to be expected due to the normal state hosting low-energy flat-energy bands [74,75], which effectively make even very small electron interactions dominate the kinetic energy. Taken together, understanding the underlying physics of carbon-based superconductors may help unveil the mechanisms at the root of high-temperature SC and thereby also eventually increase the presently accessible SC critical temperatures, which is crucial for technological developments.

In this paper, we focus on revealing all possible symmetry-allowed SC order parameters [3,10,76] in a tight-binding model, without focusing on their possible origins nor on estimating the exact values of the critical temperature, $T_c$, that could be reached in realistic systems. We instead use a *priori* generic values for the pairing amplitudes and review the resulting basic electronic properties, such as the band structure and the density of states. We start for simplicity with monolayer graphene [2], but the goal of our analysis is to also access multilayer systems [77,78], with both ABA and ABC stacking, particularly focusing on the latter, since it has already been demonstrated to exhibit SC properties [24]. We claim that by comparing our results with experimental measurements, one can identify the underlying SC order parameter in a graphitic system.
Overall, our project consists of three parts: in the present paper we focus on calculations of the band structure and density of states, with a particular focus on the gap closing points in the energy spectrum, as well as on evaluating the effects of trigonal warping present in bi- and trilayer graphene. In the upcoming two works we will focus first on the topological properties of the SC state in monolayer and multilayer graphene and on the corresponding edges states [79] and, secondly, on the impurity-induced [21,80–83] subgap states appearing in the presence of impurities [84]. We show that the results of these calculations depend strongly on the underlying SC order parameter, and thus a comparison with experimental measurements would allow us to determine the SC symmetry as well as type of SC pairing.

By using a tight-binding formalism and modeling the SC pairing in real space as on-site (ON), nearest-neighbor (NN), and next-nearest-neighbor (NNN) couplings, in both the spin-singlet and spin-triplet channels, we capture all relevant symmetry possibilities, as classified by group theory [76]. In particular, the possible SC order parameters can be split into spin-singlet pairing with s-wave (both a constant gap and extended s-wave) and d-wave ($d_{xy}$, $d_{x^2-y^2}$, and $d+id'$) spatial symmetry and into spin-triplet pairings with p-wave ($p_x$, $p_y$, and $p+ip'$) and f-wave spatial symmetry. We calculate the full band structure for each pairing possibility and symmetry; in particular, we focus on the dependence of the energy of the lowest (positive energy) band as a function of momentum and distinguish between nodal SC states ($d_{xy}$, $d_{x^2-y^2}$, $p_x$, and $p_y$), which all break the rotational symmetry of the normal state around both the Dirac points and the Brillouin zone center, and fully gapped states (usually ON, $p+ip'$, $d+id'$, and f), which, in general, preserve the spatial symmetries of the normal state. We note that such differences should be feasible to detect using, e.g., angle-resolved photoemission spectroscopy, which would help identify the underlying order parameter.

As expected, the DOS for the nodal states has a linear dependence of energy at low energies (V-shaped), while the fully gapped states produce a U-shaped behavior. However, both exhibit a gap-edge coherence peak in the DOS. We analyze how the DOS evolves for each state as a function of both the chemical potential and the SC pairing strength. In particular, we identify the position of the gap closing points in the parameter space. We find that the spin-triplet $p_x$, $p_y$, and $p+ip'$ states are the most peculiar by exhibiting most of these gap closing points. Among the spin-singlet states, only the $d_{xy}$ state exhibits gap closing points; for this nodal state, we denote the point at which the gap edges merge at zero energy a gap closing point.

Finally, we find that most of the overall physical features are preserved when moving from monolayer to multilayer graphene. The notable difference is a doubling (for bilayer) or tripling (for trilayer) of the number of nodal points, all appearing in close proximity to each other in the Brillouin zone. This also leads to an increase in the number of gap closing points as a function of the chemical potential and the SC pairing strength. We note that, interestingly enough, trigonal warping has a significant effect on the number of gap closing points and can, in fact, also greatly reduce the number of gap closing points.

<figure><img src="image_5.png" /></figure>

$$ \text{ON: } \vec{\Delta}_0 = \Delta_{g_1}^{(0)} $$

$$ \text{NN: } \vec{\Delta}_1 = (\Delta_{g_1}^{(1)}, \Delta_{g_2}^{(1)}, \Delta_{g_3}^{(1)})^T $$

$$ \text{NNN: } \vec{\Delta}_2 = (\Delta_{g_1}^{(2)}, \Delta_{g_2}^{(2)}, \dots, \Delta_{g_5}^{(2)}, \Delta_{g_6}^{(2)})^T $$

<caption>FIG. 1. Honeycomb lattice for monolayer graphene with $g_i^{(1)}$ and $g_i^{(2)}$ corresponding to the NN and NNN vectors. Real-space SC order parameters $\vec{\Delta}_{0,1,2}$ for ON (green), NN (black), and NNN (red) couplings are also indicated in the figure and given in vector form. Honeycomb lattice and its bond vectors are adapted from Ref. [85].</caption>

The remainder of this paper is organized as follows. In Sec. II, we present the tight-binding model for monolayer graphene and all possible SC order parameters up to NNN pairing, as well as the resulting low-energy band structure and DOS. In Sec. III, we present the equivalent model for multilayer graphene and the corresponding modifications to the low-energy band structure and the gap closing points in the DOS. We summarize our results in Sec. IV. Extra information and more calculational details are presented in the Appendices.

## II. SUPERCONDUCTING MONOLAYER GRAPHENE

### A. Tight-binding Hamiltonian

Graphene has a honeycomb hexagonal lattice with two atoms per unit cell, here denoted A and B. We take the three vectors connecting sites A to the NN sites B to be $g_1^{(1)} = (0, -1)$, $g_2^{(1)} = (-\frac{\sqrt{3}}{2}, \frac{1}{2})$ and $g_3^{(1)} = (\frac{\sqrt{3}}{2}, \frac{1}{2})$, see Fig. 1. For simplicity, we assume here that the distance between two carbon atoms, $a_0$, is equal to 1. The standard noninteracting Hamiltonian of graphene [86] involves only NN hopping and can be written as

$$ H_0 = -t \sum_{(i,j)\sigma} [a_{i\sigma}^{\dagger} b_{j\sigma} + b_{j\sigma}^{\dagger} a_{i\sigma}] - \mu \sum_{i\sigma} [a_{i\sigma}^{\dagger} a_{i\sigma} + b_{i\sigma}^{\dagger} b_{i\sigma}], \quad (1) $$

where $t$ is the NN hopping parameter between A and B sites and $\mu$ is the chemical potential. The operators $a_{i\sigma}^{\dagger}(b_{i\sigma}^{\dagger})$ create an electron ON $i$ on sublattice A (B) with spin $\sigma$. Unless otherwise specified, we will consider for simplicity $t=1$. Going to momentum space, we arrive at the well-established Bloch Hamiltonian

$$ H_0(\mathbf{k}) = \begin{pmatrix} \mu & h_0(\mathbf{k}) \\ h_0^*(\mathbf{k}) & \mu \end{pmatrix}, \quad (2) $$

$$ h_0(\mathbf{k}) = -t e^{-ik_y} \left( 1 + 2e^{3ik_y/2} \cos\left(\frac{\sqrt{3}}{2} k_x \right) \right), \quad (3) $$
which for $\mu = 0$ has a vanishing gap at two nonequivalent Dirac points in the Brillouin zone, $\mathcal{K} = (\frac{4\pi}{3\sqrt{3}}, 0)$ and $\mathcal{K}' = (-\frac{4\pi}{3\sqrt{3}}, 0)$. Near these two Dirac points, the electrons obey a linear dispersion relation leading to the famous Dirac cones that control the transport properties of undoped graphene.

## B. Superconducting term

We next consider an interacting term $H_I$, such that the total Hamiltonian is $H = H_0 + H_I$. Here we consider generic four-body interactions extending as far as NNN, which we can write as [85,87]

$$H_I = \frac{1}{2} \sum_{ij, \alpha \beta \delta \gamma} \Gamma_{ij}^{\alpha \beta \delta \gamma} f_{i\alpha}^{\dagger} f_{j\beta}^{\dagger} f_{i\delta} f_{j\gamma}, \quad (4)$$

with

$$\Gamma_{ij}^{\alpha\beta\delta\gamma} = \frac{1}{2} \sum_{\eta} [U_{ij}\sigma_{\alpha\gamma}^{0}\sigma_{\beta\delta}^{0} + J_{ij}\sigma_{\alpha\gamma}^{\eta}\sigma_{\beta\delta}^{\eta}]. \quad (5)$$

Here $\sigma^0$ and $\sigma^\eta$, $\eta \in \{x, y, z\}$ are the Pauli matrices acting in spin space, while $j$ is constrained to be equal to $i$, $i + \bar{\delta}_{NN}$ or $i + \bar{\delta}_{NNN}$, respectively, for ON, NN ($\bar{\delta}_{NN} = g_{1,2,3}^{(1)}$), and NNN interactions ($\bar{\delta}_{NNN} = g_{1,2,3,4,5,6}^{(2)}$). For notational simplicity, we use here the operators $f_{i\sigma}$, which correspond to either $a_{i,\sigma}$ or $b_{i,\sigma}$ depending on whether $i$ corresponds to an electron on an A or B sublattice. By providing the decomposition of the interaction terms in Eq. (5), we separate already at this level the effective Coulomb interactions $U_{ij}$ from the the effective spin-spin interactions $J_{ij}$.

To study superconductivity we need to carry out a mean-field decoupling of $H_I$ into the SC Cooper channels. Following Ref. [85], we introduce the SC order parameters or, equivalently, the mean-field decoupling fields and arrive at

$$H_I^{\text{MF}} = \sum_{\langle ij \rangle, \eta} [\Delta_{ji}^{\eta}(g_{ij}^{\eta})^{\dagger} + (\Delta_{ji}^{\eta})^{*} g_{ij}^{\eta}] \\ + \sum_{\langle\langle ij \rangle\rangle, \eta} [\Delta_{ji}^{\eta}(g_{ij}^{\eta})^{\dagger} + \Delta_{ji}^{\eta}(g_{ij}^{\eta})^{\dagger} + (\Delta_{ji}^{\eta})^{*} g_{ij}^{\eta}] \\ + (\Delta_{ji}^{\eta})^{*} g_{ij}^{\eta}] + 2 \sum_{i, \sigma} [\Delta_{ii}^{\eta}(g_{ii}^{0})^{\dagger} + (\Delta_{ii}^{\eta})^{*} g_{ii}^{0}], \quad (6)$$

where the SC order parameters are given by the self-consistent equations

$$\Delta_{ji}^{\eta} = \Gamma_{ij}^{\eta,+} \langle g_{ij}^{\eta} \rangle, \quad (7)$$

with

$$g_{ij}^{\eta} = \frac{1}{2} \sum_{\alpha\beta} f_{j\alpha} [i\sigma^{\gamma}\sigma^{\eta}]_{\alpha\beta} f'_{i\beta}. \quad (8)$$

In the above equations, we have an ON pairing determined by $\Delta_{ii}^{\eta}$, as well as pairings between NNs $\langle ij \rangle$ and NNNs $\langle\langle ij \rangle\rangle$, both denoted by $\Delta_{ji}^{\eta}$. As spin-orbit coupling is negligible in graphene [88] (carbon atoms are light with an atomic number of $Z = 6$), we can separate the treatment of spin-singlet SC, where $\eta = 0$, so we have the interactions $\Gamma_{ij}^{0,+} = U_{ij} - 3J_{ij}$, from spin-triplet SC, where $\eta = x, y, z$ and $\Gamma_{ij}^{\eta,+} = U_{ij} + J_{ij}$. We use here the notation $\Gamma_{ij}^{\eta,+}$ for interactions acting in the

<table><caption>TABLE I. Spin-singlet SC symmetries ($\eta = 0$) for ON, NN, and NNN pairing in the form of the $\bar{\Delta}$ order parameter (up to an overall amplitude $\Delta_0$) and their spatial symmetries in reciprocal space.</caption><thead><tr><th>Range</th><th>Δ̄</th><th>Symmetry</th></tr></thead><tbody><tr><td>ON</td><td>1</td><td>s<sub>ON</sub></td></tr><tr><td>NN</td><td>(1,1,1)<sup>T</sup> / √3</td><td>s<sub>ext</sub></td></tr><tr><td>NN</td><td>(2,-1,-1)<sup>T</sup> / √6</td><td>d<sub>x<sup>2</sup>-y<sup>2</sup></sub></td></tr><tr><td>NN</td><td>(0,-1,1)<sup>T</sup> / √2</td><td>d<sub>xy</sub></td></tr><tr><td>NNN</td><td>(1,1,1,1,1)<sup>T</sup> / √6</td><td>s<sub>ext</sub></td></tr><tr><td>NNN</td><td>(-1,2,-1,-1,2,-1)<sup>T</sup> / √12</td><td>d<sub>x<sup>2</sup>-y<sup>2</sup></sub></td></tr><tr><td>NNN</td><td>(-1,0,1,-1,0,1)<sup>T</sup> / 2</td><td>d<sub>xy</sub></td></tr></tbody></table>

SC channels, but a similar treatment can be done for $\Gamma_{ij}^{\eta,+}$ to capture putative mean-field magnetic or charge-ordered states. As we will shortly see, this model captures all reasonable pairing symmetries, including phonon-mediated superconductivity as well as states very reminiscent of the cuprate high-temperature superconductors [87]. More specifically, an ON coupling necessarily implies spin-singlet superconductivity and $i = j$, such that $\Delta_{ji}^{\eta} = \bar{\Delta}_{ON}$ is a one-component vector. For NN pairing, we can have both spin-singlet $\eta = 0$ and spin-triplet $\eta = x, y, z$ superconductivity. Here $\Delta_{ji}^{\eta}$ has three nonzero components, one on each NN bond, and can thus be written compactly as a three-component vector $\bar{\Delta}_1$, see Fig. 1. Similarly, for NNN we create a 6-component vector $\bar{\Delta}_2$, see Fig. 1, but where only three components are technically independent. Because the spin-orbit coupling is negligible in graphene, $H^{\text{MF}}$ does not break spin rotation symmetry. Consequently, the three spin-triplet order parameters are degenerate. Therefore, when it comes to the basic electronic properties, such as the band structure and DOS, it is enough to focus on only one triplet channel. In what follows, we will generally use the $\eta = x$ channel.

If we know the interaction parameters $U$ and $J$ we could use Eq. (7) to calculate self-consistently all SC order parameters $\bar{\Delta}$. But we very rarely do, and currently we have very little notion of these interaction parameters in SC graphene systems. Instead, we concentrate in this paper on the symmetries and overall properties of all relevant SC states. As such, we assume a given overall amplitude $\Delta_0$, for either $\bar{\Delta}_{ON}$, $\bar{\Delta}_1$, or $\bar{\Delta}_2$. Moreover, we know that at the critical temperature, the SC order parameter is forced to belong to one of the irreducible representations of the crystal symmetry group [76], in this case the hexagonal lattice group $D_{6h}$. This means we can readily quantify all allowed symmetries of the SC state by using the irreducible representations of the crystal symmetry group. Starting with ON pairing, we further can have only spin-singlet s-wave symmetry as the pairing is fully localized in real space, and hence uniform in momentum space. This state belongs to the trivial irreducible representation. For NN and NNN pairing, the possibilities include both spin-singlet and spin-triplet states, with the lowest spatial harmonics going up to f-wave symmetry (higher harmonics induce more nodes and will thus not be energetically favored in general). In Table I, we summarize the allowed $\bar{\Delta}$ and their symmetries
<table><caption>TABLE II. Spin-triplet SC symmetries ($\eta \in \{x, y, z\}$) for ON, NN, and NNN pairing in the form of the $\tilde{\Delta}$ order parameter (up to an overall amplitude $\Delta_0$) and their spatial symmetries in reciprocal space.</caption><thead><tr><th>Range</th><th>Δ̃</th><th>Symmetry</th></tr></thead><tbody><tr><td>NN</td><td>(2, -1, -1)<sup>T</sup></td><td>p<sub>y</sub></td></tr><tr><td>NN</td><td>(0, -1, 1)<sup>T</sup></td><td>p<sub>x</sub></td></tr><tr><td>NNN</td><td>(1, -1, -1, -1, -1)<sup>T</sup></td><td>f<sub>x</sub>(x<sup>2</sup> - 3y<sup>2</sup>)</td></tr><tr><td>NNN</td><td>(-1, -2, -1, 1, 2, 1)<sup>T</sup></td><td>p<sub>x</sub></td></tr><tr><td>NNN</td><td>(-1, 0, 1, 1, 0, -1)<sup>T</sup></td><td>p<sub>y</sub></td></tr></tbody></table>

for ON, NN, and NNN pairing for spin-singlet pairing ($\eta = 0$) and in Table II for spin-triplet pairing ($\eta \in \{x, y, z\}$) [10,85]. We note that the extracted symmetries are those obtained for the intraband pairing in reciprocal space and around the Brillouin zone center, see further Sec. II D.

### C. Bogoliubov–de-Gennes Hamiltonian

To proceed, we rewrite the total Hamiltonian $H = H_0 + H^MF$ in momentum space in the compact Bogoliubov–de-Gennes (BdG) matrix form. Writing the momentum-space Hamiltonian in a matrix form will later allow us to diagonalize $H$ and obtain the band structure. For this purpose, we introduce the spinor

$$ \psi_{\mathbf{k}} = \{a_{\mathbf{k}\uparrow}, b_{\mathbf{k}\uparrow}, a_{\mathbf{k}\downarrow}, b_{\mathbf{k}\downarrow}, a_{-\mathbf{k}\uparrow}^{\dagger}, b_{-\mathbf{k}\uparrow}^{\dagger}, a_{-\mathbf{k}\downarrow}^{\dagger}, b_{-\mathbf{k}\downarrow}^{\dagger}\}^{\mathrm{T}}, \quad (9) $$

with $a_{\mathbf{k}\sigma}$ and $b_{\mathbf{k}\sigma}$ the usual annihilation operators, but now $\mathbf{k}$ belonging to the first Brillouin zone. With this notation, we express the total Hamiltonian:

$$ H = \sum_{\mathbf{k}} \psi_{\mathbf{k}}^{\dagger} H_{\text{BdG}}(\mathbf{k}) \psi_{\mathbf{k}}. \quad (10) $$

Note that the dimension of the matrix associated to the BdG Hamiltonian is doubled compared to the standard BdG Hamiltonian, since we take into account separately the up and down spins, as well as the electrons and holes.

The normal-state terms within the BdG form can be written as (see Appendix A for details)

$$ H_0 = \sum_{\mathbf{k},\sigma} h_0(\mathbf{k}) (a_{\mathbf{k}\sigma}^{\dagger} b_{\mathbf{k}\sigma} - b_{-\mathbf{k},\sigma} a_{-\mathbf{k},\sigma}^{\dagger}) \\ + \sum_{\mathbf{k},\sigma} h_0^*(\mathbf{k}) (b_{\mathbf{k}\sigma}^{\dagger} a_{\mathbf{k}\sigma} - a_{-\mathbf{k},\sigma} b_{-\mathbf{k},\sigma}^{\dagger}) \\ - \mu \sum_{\mathbf{k}\sigma} (a_{\mathbf{k},\sigma}^{\dagger} a_{\mathbf{k},\sigma} - a_{-\mathbf{k},\sigma} a_{-\mathbf{k},\sigma}^{\dagger} + b_{\mathbf{k},\sigma}^{\dagger} b_{\mathbf{k},\sigma} - b_{-\mathbf{k},\sigma}^{\dagger} b_{-\mathbf{k},\sigma}^{\dagger}). \quad (11) $$

We next write all the SC terms within the BdG form. The details of the calculations are given in Appendix A. Writing out everything explicitly for NN pairing, we have in real space for spin-singlet pairing ($\eta = 0$):

$$ H_{\mathrm{NN}}^{0} = \sum_{\langle ij \rangle} \Delta_{ij}^{0} (a_{i\uparrow}^{\dagger} b_{j\downarrow}^{\dagger} - a_{i\downarrow}^{\dagger} b_{j\uparrow}^{\dagger}) + \mathrm{H.c.}, \quad (12) $$

<table><caption>TABLE III. Form factors $h_{\mathrm{NN}}^0(\mathbf{k})$ and $h_{\mathrm{NNN}}^0(\mathbf{k})$ for each spin-singlet SC symmetry ($\eta = 0$) for NN and NNN pairing. Here $\tilde{h}_0(\mathbf{k}) = \frac{\tilde{h}_0(\mathbf{k})}{i}$.</caption><thead><tr><th>Symmetry</th><th>Range</th><th>Form factor</th></tr></thead><tbody><tr><td>s<sub>ext</sub></td><td>NN</td><td>h<sub>NN</sub><sup>0,ext</sup>(<b>k</b>) = <sup>&Delta;<sub>0</sub></sup>&frasl;<sub>&radic;3</sub> h&#773;<sub>0</sub>(<b>k</b>)</td></tr><tr><td>d<sub>x<sup>2</sup>-y<sup>2</sup></sub></td><td>NN</td><td>h<sub>NN</sub><sup>0,d<sub>x<sup>2</sup>-y<sup>2</sup></sub></sup>(<b>k</b>) = <sup>2&Delta;<sub>0</sub></sup>&frasl;<sub>&radic;6</sub> e<sup>-ik<sub>y</sub></sup> [1 - e<sup>-3i/k<sub>y</sub></sup> cos(<sup>&radic;3</sup>&frasl;<sub>2</sub>k<sub>x</sub>)]</td></tr><tr><td>d<sub>xy</sub></td><td>NN</td><td>h<sub>NN</sub><sup>0,dxy</sup>(<b>k</b>) = &Delta;<sub>0</sub>&radic;2i e<sup>i/k<sub>y</sub></sup> sin(<sup>&radic;3</sup>&frasl;<sub>2</sub>k<sub>x</sub>)</td></tr><tr><td>s<sub>ext</sub></td><td>NNN</td><td>h<sub>NNN</sub><sup>0,ext</sup>(<b>k</b>) = <sup>2&Delta;<sub>0</sub></sup>&frasl;<sub>&radic;6</sub> [cos(&radic;3k<sub>x</sub>) + 2 cos(<sup>&radic;3</sup>&frasl;<sub>2</sub>k<sub>x</sub>) cos(<sup>3</sup>&frasl;<sub>2</sub>k<sub>y</sub>)]</td></tr><tr><td>d<sub>x<sup>2</sup>-y<sup>2</sup></sub></td><td>NNN</td><td>h<sub>NNN</sub><sup>0,d<sub>x<sup>2</sup>-y<sup>2</sup></sub></sup>(<b>k</b>) = <sup>4&Delta;<sub>0</sub></sup>&frasl;<sub>&radic;12</sub> [cos(&radic;3k<sub>x</sub>) - cos(<sup>&radic;3</sup>&frasl;<sub>2</sub>k<sub>x</sub>) cos(<sup>3</sup>&frasl;<sub>2</sub>k<sub>y</sub>)]</td></tr><tr><td>d<sub>xy</sub></td><td>NNN</td><td>h<sub>NNN</sub><sup>0,dxy</sup>(<b>k</b>) = -2&Delta;<sub>0</sub> sin(<sup>3</sup>&frasl;<sub>2</sub>k<sub>y</sub>) sin(<sup>&radic;3</sup>&frasl;<sub>2</sub>k<sub>x</sub>)</td></tr></tbody></table>

and for spin-triplet pairing ($\eta = x, y, z$) we have

$$ H_{\mathrm{NN}}^{x} = \sum_{\langle i,j \rangle} \Delta_{ij}^{x} (a_{i\uparrow}^{\dagger} b_{j\uparrow}^{\dagger} - a_{i\downarrow}^{\dagger} b_{j\downarrow}^{\dagger}) + \mathrm{H.c.}, \quad (13) $$

$$ H_{\mathrm{NN}}^{y} = \sum_{\langle i,j \rangle} \Delta_{ij}^{y} (a_{i\uparrow}^{\dagger} b_{j\uparrow}^{\dagger} + a_{i\downarrow}^{\dagger} b_{j\downarrow}^{\dagger}) + \mathrm{H.c.}, \quad (14) $$

$$ H_{\mathrm{NN}}^{z} = \sum_{\langle i,j \rangle} \Delta_{ij}^{z} (a_{i\uparrow}^{\dagger} b_{j\downarrow}^{\dagger} + a_{i\downarrow}^{\dagger} b_{j\uparrow}^{\dagger}) + \mathrm{H.c.} \quad (15) $$

By Fourier transforming these into reciprocal space, we obtain

$$ H_{\mathrm{NN}}^{0} = \sum_{\mathbf{k}} h_{\mathrm{NN}}^{0}(\mathbf{k}) (a_{\mathbf{k}\uparrow}^{\dagger} b_{-\mathbf{k}\downarrow}^{\dagger} - a_{\mathbf{k}\downarrow}^{\dagger} b_{-\mathbf{k}\uparrow}^{\dagger}) + \mathrm{H.c.} \quad (16) $$

and

$$ H_{\mathrm{NN}}^{x} = \sum_{\mathbf{k}} h_{\mathrm{NN}}^{x}(\mathbf{k}) (a_{\mathbf{k}\uparrow}^{\dagger} b_{-\mathbf{k}\uparrow}^{\dagger} - a_{\mathbf{k}\downarrow}^{\dagger} b_{-\mathbf{k}\downarrow}^{\dagger}) + \mathrm{H.c.}, \quad (17) $$

$$ H_{\mathrm{NN}}^{y} = i \sum_{\mathbf{k}} h_{\mathrm{NN}}^{y}(\mathbf{k}) (a_{\mathbf{k}\uparrow}^{\dagger} b_{-\mathbf{k}\uparrow}^{\dagger} + a_{\mathbf{k},\downarrow}^{\dagger} b_{-\mathbf{k},\downarrow}^{\dagger}) + \mathrm{H.c.}, \quad (18) $$

$$ H_{\mathrm{NN}}^{z} = \sum_{\mathbf{k}} h_{\mathrm{NN}}^{z}(\mathbf{k}) (a_{\mathbf{k}\uparrow}^{\dagger} b_{-\mathbf{k}\uparrow}^{\dagger} + a_{\mathbf{k},\downarrow}^{\dagger} b_{-\mathbf{k},\downarrow}^{\dagger}) + \mathrm{H.c.} \quad (19) $$

Here $h_{\mathrm{NN}}^{\eta}(\mathbf{k})$ are the overall form factors whose expressions depend on both the spin channel and the spatial symmetry of the order parameter. Their general expression is

$$ h_{\mathrm{NN}}^{\eta}(\mathbf{k}) = \Delta_{\mathrm{NN}}^{\eta,d=1} e^{-ik_y} + \Delta_{\mathrm{NN}}^{\eta,d=2} e^{\frac{i}{2}k_y - \frac{\sqrt{3}i}{2}k_x} \\ + \Delta_{\mathrm{NN}}^{\eta,d=3} e^{\frac{i}{2}k_y + \frac{\sqrt{3}i}{2}k_x}, \quad (20) $$

where $d = 1, 2, 3$ correspond to the three NN bonds, following the convention of Fig. 1 and Tables I and II. We summarize $h_{\mathrm{NN}}^{\eta}(\mathbf{k})$ for each symmetry in Table III for spin-singlet pairing and in Table IV for spin-triplet pairing. In these tables, we also present the corresponding results for $h_{\mathrm{NNN}}^{\eta}(\mathbf{k})$ (see Appendix A for details).

For completeness, we write out the full BdG Hamiltonian in each spin channel in their matrix form. For the spin-singlet
<table><caption>TABLE IV. Form factors $h_{NN}^{\eta}(\mathbf{k})$ and $h_{NNN}^{\eta}(\mathbf{k})$ for each spin-triplet SC symmetry ($\eta = x, y, z$) for NN and NNN pairing.</caption><thead><tr><th>Symmetry</th><th>Range</th><th>Form factor</th></tr></thead><tbody><tr><td>$p_y$</td><td>NN</td><td>$h_{NN}^{\eta, p_y}(\mathbf{k}) = \frac{2\Delta_0}{\sqrt{6}}e^{-ik_y}[1 - e^{\frac{3i}{2}k_y}\cos(\frac{\sqrt{3}}{2}k_x)]$</td></tr><tr><td>$p_x$</td><td>NN</td><td>$h_{NN}^{\eta, p_x}(\mathbf{k}) = i\sqrt{2}\Delta_0 e^{\frac{i}{2}k_y}\sin(\frac{\sqrt{3}}{2}k_x)$</td></tr><tr><td>$f_{x}(x^2-3y^2)$</td><td>NNN</td><td>$h_{NNN}^{\eta, p_x}(\mathbf{k}) = \frac{2i\Delta_0}{\sqrt{6}}[\sin(\sqrt{3}k_x) - 2\sin(\frac{\sqrt{3}}{2}k_x)\cos(\frac{\sqrt{3}}{2}k_y)]$</td></tr><tr><td>$p_x$</td><td>NNN</td><td>$h_{NNN}^{\eta, p_x}(\mathbf{k}) = \frac{4i\Delta_0}{\sqrt{12}}[\cos(\sqrt{3}k_x) + \cos(\frac{3}{2}k_y)]\sin(\frac{\sqrt{3}}{2}k_x)$</td></tr><tr><td>$p_y$</td><td>NNN</td><td>$h_{NNN}^{\eta, p_y}(\mathbf{k}) = -2i\Delta_0\sin(\frac{3}{2}k_y)\cos(\frac{\sqrt{3}}{2}k_x)$</td></tr></tbody></table>

pairing ($\eta = 0$), we obtain

$$ \begin{pmatrix} H_0(\mathbf{k}) & \mathbf{0}_{2 \times 2} & \mathbf{0}_{2 \times 2} & -h_{\Delta}^0(\mathbf{k}) \\ \mathbf{0}_{2 \times 2} & H_0(\mathbf{k}) & h_{\Delta}^0(\mathbf{k}) & \mathbf{0}_{2 \times 2} \\ \mathbf{0}_{2 \times 2} & (h_{\Delta}^0(\mathbf{k}))^{\dagger} & -H_0(\mathbf{k}) & \mathbf{0}_{2 \times 2} \\ -(h_{\Delta}^0(\mathbf{k}))^{\dagger} & \mathbf{0}_{2 \times 2} & \mathbf{0}_{2 \times 2} & -H_0(\mathbf{k}) \end{pmatrix}, \quad (21) $$

and for $\eta = x$ spin-triplet pairing, we have

$$ \begin{pmatrix} H_0(\mathbf{k}) & \mathbf{0}_{2 \times 2} & -h_{\Delta}^x(\mathbf{k}) & \mathbf{0}_{2 \times 2} \\ \mathbf{0}_{2 \times 2} & H_0(\mathbf{k}) & \mathbf{0}_{2 \times 2} & h_{\Delta}^x(\mathbf{k}) \\ -(h_{\Delta}^x(\mathbf{k}))^{\dagger} & \mathbf{0}_{2 \times 2} & -H_0(\mathbf{k}) & \mathbf{0}_{2 \times 2} \\ \mathbf{0}_{2 \times 2} & (h_{\Delta}^x(\mathbf{k}))^{\dagger} & \mathbf{0}_{2 \times 2} & -H_0(\mathbf{k}) \end{pmatrix}, \quad (22) $$

where $\mathbf{0}_{2 \times 2}$ is a $2 \times 2$ null matrix, $H_0(\mathbf{k})$ is the normal state Hamiltonian matrix given in Eq. (2), and

$$ h_{\Delta}^{0}(\mathbf{k}) = \frac{1}{2} \begin{pmatrix} \Delta_{\text{ON}} + h_{\text{NNN}}^{0}(-\mathbf{k}) & h_{\text{NN}}^{0}(\mathbf{k}) \\ h_{\text{NN}}^{0}(-\mathbf{k}) & \Delta_{\text{ON}} + h_{\text{NNN}}^{0}(-\mathbf{k}) \end{pmatrix}, $$

$$ h_{\Delta}^{x}(\mathbf{k}) = \frac{1}{2} \begin{pmatrix} h_{\text{NN}}^{x}(-\mathbf{k}) & -h_{\text{NN}}^{x}(\mathbf{k}) \\ h_{\text{NN}}^{x}(-\mathbf{k}) & h_{\text{NN}}^{x}(-\mathbf{k}) \end{pmatrix} $$

are the corresponding SC order parameter matrices. The SC order parameter matrices for the $\eta = y, z$ spin-triplet pairing are given in Appendix B. Diagonalizing the matrices in Eqs. (21) and (22) yields both the full energy spectrum and the eigenstates of the various SC states.

### D. Intraband pairing symmetries

While the matrix forms in Eqs. (21) and (22) allow for straightforward numerical diagonalization to easily find, e.g., the energy spectrum, it is still useful to first analyze the order parameter in more detail. For this purpose, it is beneficial to not work in the sublattice basis with the operators $a_{k\sigma}$, $b_{k\sigma}$, but to instead work in the band basis with the operators $c_{k\sigma}$, $d_{k\sigma}$ for the two bands, where the normal-state Hamiltonian $H_0$ is diagonal, with the band energies $\epsilon_c(\mathbf{k})$, $\epsilon_d(\mathbf{k})$ as diagonal entries. In particular, in the band basis we know that $\epsilon_c(\mathbf{k}) = 0$ (or $\epsilon_d(\mathbf{k}) = 0$, depending on the value of $\mu$) on the Fermi surface. In graphene, the Fermi surface forms circles around the Dirac points at $\mathcal{K}$ and $\mathcal{K}'$ with increasing radius as the chemical potential increases from 0. At $|\mu| = 1$, these two

Dirac Fermi surfaces meet and instead form a separatrix line joining the M points in the Brillouin zone. Finally, for $|\mu| > 1$, the Fermi surface becomes centered around the Brillouin zone center $\Gamma$.

Setting $\alpha_{\mathbf{k}} = (a_{k\sigma}, b_{k\sigma})^T$ and $\chi_{\mathbf{k}} = (c_{k\sigma}, d_{k\sigma})^T$, we find that the basis change can be expressed as

$$ \alpha_{\mathbf{k}} = \hat{U}(\mathbf{k})\chi_{\mathbf{k}}, \quad (23) $$

with the unitary matrix

$$ \hat{U}(\mathbf{k}) = \frac{1}{\sqrt{2}} \begin{pmatrix} -e^{-i\varphi_k} & e^{-i\varphi_k} \\ 1 & 1 \end{pmatrix}, \quad (24) $$

where $\varphi_k = \arg(h_0(\mathbf{k}))$. Using $\hat{U}(\mathbf{k})$ to transform also the SC terms from the sublattice basis of the previous subsections into the band basis, we arrive for spin-singlet ON pairing at

$$ H_{\text{ON}}^{0} = \Delta_{0} \sum_{\mathbf{k}} (c_{\mathbf{k}\uparrow}^{\dagger} c_{-\mathbf{k}\downarrow}^{\dagger} + d_{\mathbf{k}\downarrow}^{\dagger} d_{-\mathbf{k}\uparrow}^{\dagger}) + \text{H.c.}, \quad (25) $$

while for spin-singlet NN pairing we get

$$ H_{\text{NN}}^{0} = \Delta_{0} \sum_{\mathbf{k},d} \Delta_{d} [\cos(\mathbf{k} \cdot \mathbf{g}_{d}^{(1)} - \varphi_{k}) (d_{\mathbf{k}\downarrow}^{\dagger} d_{-\mathbf{k}\uparrow}^{\dagger} - c_{\mathbf{k}\uparrow}^{\dagger} c_{-\mathbf{k}\downarrow}^{\dagger}) \\ + i \sin(\mathbf{k} \cdot \mathbf{g}_{d}^{(1)} - \varphi_{k}) (c_{\mathbf{k}\uparrow}^{\dagger} d_{-\mathbf{k}\downarrow}^{\dagger} - d_{\mathbf{k}\uparrow}^{\dagger} c_{-\mathbf{k}\downarrow}^{\dagger})] + \text{H.c.}, \quad (26) $$

and for spin-triplet NN pairing with $\eta = z$:

$$ H_{\text{NN}}^{z} = \Delta_{0} \sum_{\mathbf{k},d} \Delta_{d} [i \sin(\mathbf{k} \cdot \mathbf{g}_{d}^{(1)} - \varphi_{k}) (d_{\mathbf{k}\downarrow}^{\dagger} d_{-\mathbf{k}\uparrow}^{\dagger} - c_{\mathbf{k}\uparrow}^{\dagger} c_{-\mathbf{k}\downarrow}^{\dagger}) \\ + \cos(\mathbf{k} \cdot \mathbf{g}_{d}^{(1)} - \varphi_{k}) (c_{\mathbf{k}\uparrow}^{\dagger} d_{-\mathbf{k}\downarrow}^{\dagger} - d_{\mathbf{k}\uparrow}^{\dagger} c_{-\mathbf{k}\downarrow}^{\dagger})] + \text{H.c.} \quad (27) $$

Here $d = 1, 2, 3$ mark the three NN bonds, following the convention of Fig. 1, and $\Delta_d$ is the $d$th component of the NN $\bar{\Delta}$ bond order parameter, given in Table I.

From Eqs. (26) and (27), we see that in the band basis superconductivity consists of both intraband pairing through the terms $c_{\mathbf{k}\uparrow}^{\dagger}c_{-\mathbf{k}\downarrow}^{\dagger}$ and $d_{\mathbf{k}\uparrow}^{\dagger}d_{-\mathbf{k}\downarrow}^{\dagger}$ for bands c and d, respectively, and interband pairing through terms such as $c_{\mathbf{k}\uparrow}^{\dagger}d_{-\mathbf{k}\downarrow}^{\dagger}$. Generally, the existence of both intra- and interband pairing means that the energy spectrum of the BdG Hamiltonian cannot simply be expressed as $E = \pm\sqrt{\epsilon^2 + |\Delta|^2}$, but instead both the intra- and interband pairing terms contribute to the energy E. Here we are not concerned with deriving the exact expressions, but instead note that for a qualitative understanding it has in the past often been enough to only consider the intraband pairing and ignore the interband term [3,10,11]. This can be understood by noticing that only one of the bands crosses zero energy (except at the Dirac point), and therefore contributes to the formation of the Fermi surface. Thus, pairing within that band is the most important contribution, while pairing in the other band and in-between bands are often of less importance. It is also the symmetry of the intraband order parameters in the Brillouin zone that we use for naming the different possible symmetries as s, d, p, and f waves in Tables I–IV.

Based on the above arguments, we analyze the intraband SC order parameters for both NN and NNN pairing. In Fig. 2,
<figure><img src="image_6.png" /><figcaption>FIG. 2. Intraband SC order parameter form factors for NN ($s_{ext}$-, $d$-, $p$-wave) pairing. Top (bottom) row corresponds to spin-singlet (spin-triplet) pairing. Cyan hexagon denotes the first Brillouin zone, while yellow (black) circles centered around the $K, K'$ points ($\Gamma$ point) are the normal-state Fermi surfaces for $\mu = 0.4t, (\mu = 1.2t)$. White regions correspond to nodes, or zeros, of the SC order parameter.</figcaption></figure>

we plot the intraband SC order parameter symmetry for NN pairing or, equivalently, the form factors in front of the intraband terms in Eqs. (26) and (27). We also mark the Fermi surface for both $\mu = 0.4$ (yellow lines) and $\mu = 1.2$ (black lines). We first note that the extended $s$ wave follows the symmetry of the normal state, with nodal points only at the Dirac points. Assuming a simplified energy dispersion, $E = \pm\sqrt{\epsilon^2 + |\Delta|^2}$, with $\Delta$ only given by the intraband component, we would then expect a fully gapped energy spectrum for any finite $\mu$. Moving on to the $d_{x^2-y^2}$- and $d_{xy}$-wave order parameters, we note that, for large chemical potential $|\mu| > 1$, when the Fermi surface is centered around $\Gamma$ (black lines), the SC state has $d$-wave symmetry, and we thus expect four nodal points in the band structure, one at each of the intersections between the Fermi surface and the zero-gap lines (white in Fig. 2). In contrast, around $K, K'$, the $d_{x^2-y^2}$-wave order has an effective $p_y$-wave symmetry, while the $d_{xy}$-wave order has an effective $p_x$-wave symmetry. As such, for $|\mu| < 1$, when the Fermi surface is centered around $K, K'$ (yellow lines), one expects two nodal points per Dirac cone. Here we note that the two Dirac cones have opposite signs on their effective $p$-wave order, such that the spin-singlet state still has an overall even spatial symmetry, as required by Fermi-Dirac statistics. We can also combine the two $d$-wave order parameters into a chiral combination $d_{x^2-y^2+ixy} = d + id'$, resulting in a fully gapped order parameter with restored rotational symmetry as seen in Fig. 2. It is only at the $\Gamma, K, K'$ points that the order parameter is equal to zero. This full gap is the reason why the chiral $d + id'$-wave combination is often found to be most stable among all the $d$-wave states [3, 4, 7, 8, 12, 89, 90]. Moving on to the NN spin-triplet order parameters, we find that the $p_y$- and $p_x$-wave order parameters have a $p$-wave symmetry both around the $\Gamma$ point and the $K, K'$ Dirac points. This results in two nodal points per Fermi surface when it is centered around the Dirac points and also two nodal points

for a Fermi surface centered around $\Gamma$. We also note that the $p$-wave states host additional nodal points, in particular, at the $M$ points. This is a consequence of adopting a $p$-wave symmetry to the sixfold symmetric honeycomb lattice with its hexagonal first Brillouin zone. Finally, the chiral combination $p_x + ip_y = p + ip'$ [13] is again fully gapped, except at the $\Gamma, K, K'$, as well as at the $M$ points where both individual $p$-wave components have nodes. Thus, we expect the $p + ip'$ SC to be fully gapped except when the chemical potential $|\mu| = 0, 1, 3$.

We also perform a similar analysis for NNN pairing, for which we in fact find that the interband pairing terms are identically zero. Thus, the energy $E$ in the SC state is exactly given by $E = \pm\sqrt{\epsilon^2 + |\Delta|^2}$, with $\epsilon$ the energy of the normal state and $\Delta$ the intraband order parameter. This can be understood by noting that the NNN pairing couples electrons on the same sublattice and thus there is no sublattice mixing or any interband pairing terms. In Fig. 3, we plot the resulting NNN intraband order parameter symmetries. Overall, we find very similar results to the NN pairing with only a few exceptions: The extended $s$-wave state has a peculiar nodal structure centered around $\Gamma$. However, unless $\mu$ is fine-tuned such that the Fermi surface exactly hits this node, the SC state will be fully gapped. Also, for NNN pairing, an $f$-wave symmetry is allowed. This order parameter displays the required $f_x = f_{x(x^2-3y^2)}$-wave symmetry, with a total of six nodal points for a Fermi surface centered around $\Gamma$ (black line), but is notably fully gapped on Fermi surfaces centered around the $K, K'$ Dirac points. Technically, there also exists another, $f_y = f_{y(y^2-3x^2)}$-wave, state that has nodal lines going through the $K, K'$ Dirac points, such that the energy for $|\mu| < 1$ has six nodes per Dirac cone. However, due to its high number of nodes per Dirac cone, this state will be much less energetically favorable at all doping levels $|\mu| < 1$ and we do not consider it in this paper.
<figure><img src="image_7.png" /><figcaption>FIG. 3. Intraband SC order parameter form factors for NNN ($s_{ext}$, d-, p-, $f_x$-wave) pairing. Top (bottom) row corresponds to spin-singlet (spin-triplet) pairing. Cyan hexagon denotes the first Brillouin zone, while yellow (black) circles centered around the $K, K'$ points ($\Gamma$ point) are the normal-state Fermi surfaces for $\mu = 0.4$, ($\mu = 1.2$). White region corresponds to the nodes, or zeros, of the SC order parameter.</figcaption></figure>

### E. Lowest energy bands

Having analyzed the symmetry of the intraband order parameter in detail in the previous subsection, we now turn to the complete solution, attained by diagonalizing the BdG Hamiltonian in Eqs. (21) and (22) for spin-singlet and spin-triplet pairing, respectively. In this subsection, we are interested in the low-energy band structure, which we analyze by plotting the lowest energy band ($E > 0$) as a function of momentum. These results directly tell us about both the existence of nodal points in the energy spectrum and the overall symmetry of the SC state. We focus primarily on ON s-wave, NN s-, d-, and p-wave, and NNN f-wave SC states. Due to the similarities between the remaining NNN order parameters and the corresponding NN ones, we expect to cover all relevant behavior with this selection.

In Fig. 4, we plot the lowest energy band as a function of $k_x$ and $k_y$ for $\mu = 0.4$. This value corresponds to a normal-state Fermi surface that is a small circle centered around the Dirac points $K, K'$, see Fig. 2. As expected, and consistent with the intuitive analysis in the previous subsection, we find that the ON s-wave and NN $s_{ext}$-wave SC states preserve the rotational symmetry in reciprocal space around the $K, K'$ Dirac points. Also, these SC states give rise to a fully gapped band structure, even though the band gap is small for $s_{ext}$ for $\Delta_0 = 0.4$ and $\mu = 0.4$ as shown in Fig. 4(b). This small gap at low $\mu$ comes from the fact that $s_{ext}$ has the same symmetry as the normal state On the other hand, we note that the $d_{xy}$, $d_{x^2-y^2}$, $p_x$, and $p_y$-wave SC states all break the rotation symmetry around the $K, K'$ points. In particular, we observe a nodal energy spectrum with two nodal points per Fermi surface arranged in agreement with the nodal structure of the intraband order parameters in Fig. 2. These nodal points can be viewed as the split of the original normal-state Dirac point into two points in the SC state. Thus the $d_{xy}$, $d_{x^2-y^2}$, $p_x$, and $p_y$-wave SC states all correspond to nodal superconductors.

By continuing with analyzing the chiral $d + id'$- and $p + ip'$-wave SC states, we find fully gapped energy spectra and a restored rotational symmetry around the $K, K'$ points (up to the same threefold symmetry as for the s-wave state). We note that the convention to generate these states is to consider both a $d_{xy}$ and $d_{x^2-y^2}$ amplitude of $\Delta_0$ for the $d + id'$ state, and similarly for the $p + ip'$ state. This convention corresponds to an effective coupling $\sqrt{2}$ stronger for the chiral states than for the nodal ones. Finally, for the NNN $f_x(x^2-3y^2)$-wave state, we find a fully gapped state with preserved rotation symmetry around the $K, K'$ points. Again, all these results are consistent with the intraband order parameter picture discussed in Fig. 2.

In Fig. 4, we focus on Fermi surfaces centered around the $K, K'$ points. For completeness, we report in Appendix C the equivalent results for $\mu = 1$, where the Fermi level sits at the VHS with the Fermi surface forming a separatrix line connecting the M points. Similarly, we report the results for NNN pairing in Appendix D, with the conclusion that NN and NNN pairing host very similar low-energy band structures.

### F. Density of states and gap closing points

Having studied the lowest energy bands, we next turn to evaluating the DOS at low energies. For this purpose, we define the Green's function of the Bloch Hamiltonian $H_{BdG}(q)$ as

$$G(q, E) = (E - H_{BdG}(q) + i0^{+})^{-1}. \qquad (28)$$

We then integrate $G(q, E)$ over the first Brillouin zone to get the DOS $\rho_0(E)$ [91],

$$\rho_0(E) = -\frac{1}{\pi} \text{Tr}_{\text{el}} \left[ \text{Im} \int_{\text{BZ}} dq G(q, E) \right], \qquad (29)$$

where the trace Tr$_{el}$ is performed over the electron modes only (first half of the diagonal of the Green's function matrix).
<figure><img src="image_8.png" /><figcaption>FIG. 4. Lowest energy band for $\mu = 0.4$, $\Delta_0 = 0.4$ and for SC states with ON s-wave, NN $s_{ext-}$, $d_{xy-}$, $d_{x^2-y^2-}$, $p_x-$, $p_y-$, $p + ip'$, $d + id'$-wave, as well as NNN $f_{x(x^2-3y^2)-}$-wave symmetries. Dark blue corresponds to zero energy. The inset of (b) presents the energy dispersion along a line-cut across the K point, as illustrated by the black dashed line. It shows that $s_{ext}$-wave pairing is fully gapped.</figcaption></figure>

We expect the DOS for the nodal states to have a linear dependence of energy at low energies (V-shaped), while the fully gapped states produce U-shaped behavior. Indeed, this difference is visible in Fig. 5, where we plot the DOS as a function of energy for the ON s-wave and the NN $d_{xy}$-wave SC states. The s-wave state displays a full gap giving rise to an overall U-type DOS profile. Note that the position of the coherence peaks corresponds exactly to the value of the SC pairing term $\Delta_0 = 0.4$. The nodal d-wave state, on the other hand, exhibits no hard gap, but a V-type DOS centered around $E = 0$, which is due to the existence of nodal quasi-particles. Note that the DOS also exhibits coherence peaks, but at values slightly smaller than the coupling term $\Delta_0 = 0.4$, indicating that the values of the SC coupling and the energy of the coherence peaks, while proportional, are not always exactly the same, as it can also be seen in Fig. 6. Moreover, we note an additional near V-type feature at higher energies, centered around $E = \mu = 0.4$. This corresponds to the normal-state Dirac point appearing at finite energy due to the finite $\mu$ [21].

In Fig. 6, we plot the DOS as a function of the overall strength, or amplitude, of the SC order parameter $\Delta_0$ for the same chemical potential $\mu = 0.4$ as in Fig. 4. The overall particle-hole asymmetry, especially strong at higher energies, is due to a pronounced particle-hole asymmetry in the normal state for finite $\mu$. Focusing primarily on low energies, as relevant for superconductivity, we note that for the ON s-wave state, the energy gap is roughly linear in $\Delta_0$. We observe a similar linear dependence for the gap edges (or coherence peaks) as a function of $\Delta_0$ for the NN $s_{ext-}$, $d_{xy-}$, $d_{x^2-y^2-}$, $p_x-$, and $p_y$-wave SC states. Note, however, that for the $d_{xy-}$, $d_{x^2-y^2-}$, $p_x-$, and $p_y$-wave states, the intensity is not zero inside the gap, which is consistent with the observations in the previous sections that all these states exhibit nodes in the band structure.

In Fig. 6, we also note that the overall size of the gap or, equivalently, half the distance between the coherence peaks for the nodal states, is no longer equal to the overall amplitude of the SC order parameter $\Delta_0$, but generally lower, except for the ON s-wave state and the NNN f-wave state. The
<figure><img src="image_9.png" /><figcaption>FIG. 5. DOS as a function of energy for ON s-wave (a) and NN $d_{xy}$-wave (b) SC states at $\Delta_0 = 0.4$ and $\mu = 0.4$ illustrating the difference between the U-shaped and V-shaped behavior, corresponding to a full gap and a nodal spectrum, respectively.</figcaption></figure>

two exceptions have the normal-state Fermi surface being gapped out by an SC order parameter constant along the Fermi surface, which generates a gap equal to $\Delta_0$. For all the other states, the low-energy spectrum is strongly modified by having the SC order parameter vary along the normal-state Fermi surface. Still, even in these cases we can identify a gap that is initially linearly increasing with $\Delta_0$, albeit with a smaller than 1 coefficient. At larger $\Delta_0$, this linearity breaks down for some of the symmetries. This is not surprising, as $\Delta_0 \sim 1$ corresponds to the energy scale of the VHS point, $\mu = 1$, which marks the energy where the normal-state band structure dramatically changes. For example, in Fig. 6 we noticed a gap closing point for the $p_x$, $p_y$ and $p + ip'$ at $\Delta_0 \sim 1$.

We also plot in Fig. 7 the DOS as a function of chemical potential $\mu$, while instead keeping the SC amplitude $\Delta_0$ fixed. Normally, changing the chemical potential should not affect the DOS much as it only alters the underlying normal-state band structure. Indeed, this is what we see for the ON s-wave which exhibits a constant gap for a large range of $\mu$. For some of the other states, however, the intraband SC order parameters exhibit nodes at the $K, K'$ Dirac points, such that when $\mu = 0$, both the SC order parameter and the normal state energy is zero at the Dirac points. This leads to a gapless system at $\mu = 0$ and to a linear dependence of the gap for small $\mu$.

Around $\mu = 1$, we observe further substantial changes in the gap, or, equivalently, in the distance between the gap coherence peaks, for several SC states ($f_x, p_x, p_y, d_{xy}, p + ip'$). These changes consist primarily in the existence of points in parameter space where the energy gap is closing. To track these points, we plot in Fig. 8 the energy gap in the band

structure for the chiral $p + ip'$-wave state as a function of both the chemical potential $\mu$ and the SC amplitude $\Delta_0$. We note that a continuous gap closing line indeed occurs connecting $\mu = 1$ with $\Delta_0 \sim 1.2$. These gap closing points are important for experimentally distinguishing between different SC states. They are also important from the point of view of topology, as we will discuss in detail in the next part of this paper. We can trace the gap closing points in all the p-wave states to the existence of nodes in the intraband SC order parameter at the M points in the Brillouin zone, as depicted in Fig. 2. This leads to both a vanishing band energy and intraband order parameter at $\mu \approx 1$. To verify this behavior, we plot in Fig. 9 the lowest energy band for the chiral $p + ip'$-wave state at three different values of $\mu$. Keeping $\Delta_0 = 0.4$, we find that the gap closing occurs at the M points of the Brillouin zone and at $\mu \approx 0.9$, corresponding to the middle panel of Fig. 9. The difference from the intuitively expected value of $\mu \approx 1$ comes from interband pairing due to the sizable $\Delta_0$. For $\mu$ below this gap closing at M, the topology of the low-energy contours corresponds to having an origin at the $K, K'$ Dirac points, while for larger $\mu$'s the contours are centered around the $\Gamma$ point. We also note that the $d_{xy}$-wave intraband order parameter has nodes at two M points, which generates the gap closing point at $\mu = 1$ in Fig. 7. However, the $d_{x^2-y^2}$-wave state, and consequently the chiral $d + id'$-wave state, host no nodes as a function of $\mu$. This different behavior offers a clear way to distinguish between all p- and d-wave states. Overall these results show that, even if the total energy is also dependent on the interband pairing, many features can be understood simply by reasoning, using the intraband order parameter and its nodes.

To provide a comprehensive treatment, we also report the DOS for the NNN pairing in Appendix E. Here we find that s-ext-wave and d-wave states have exactly the same type of behavior as their NN counterparts, while the NNN p-wave states only show a gap closing as a function of $\mu$, at $\mu = 1$, but not as a function of the SC pairing $\Delta_0$. The difference can be explained by the influence of the interband pairing term for NN pairing but not for NNN pairing. As such, the gap closing stays at $\mu = 1$ for all values of $\Delta_0$.

### III. SUPERCONDUCTING MULTILAYER GRAPHENE

We next turn to the treatment of SC multilayer graphene. Here the 2D graphene sheets are stacked together with interlayer coupling primarily through van der Waals interaction. The naive configuration of placing two layers directly on top of each other is not energetically stable. The most stable form of stacking is instead known as AB, or Bernal stacking. Figure 10(b) illustrates this stacking: the two layers are stacked with a relative relative translation of $a_0\hat{x}$. For bilayer graphene, there is only this one choice. For trilayer graphene, we have two choices, we can either perform a $a_0\hat{x}$ translation between the first and second layers, then also a $a_0\hat{x}$ translation between the second and third layers, giving rise to a rhombohedral trilayer system, or ABC stacking, see Fig. 10(d). We can also perform a translation of $-a_0\hat{x}$ for the third layer, which yields a Bernal trilayer system, or ABA stacking, see Fig. 10(c). ABA is the energetically most stable stacking, while ABC is technically only metastable. Still, ABC-stacked
<figure><img src="image_10.png" /><figcaption>FIG. 6. DOS as a function of SC amplitude $\Delta_0$ for $\mu = 0.4$ for SC states with ON s-wave, NN $s_{ext}$, $d_{xy}$, $d_{x^2-y^2}$, $p_x$, $p_y$, $p + ip'$, $d + id'$-wave, as well as NNN $f_{x(x^2-3y^2)}$-wave symmetries. Dark blue represents zero DOS.</figcaption></figure>

multilayer graphene has been constructed using layer-by-layer deposition in multiple experiments [92,93]. ABC stacking has also been produced by imposing mechanical constraints on Bernal graphite [94].

To be able to model the normal state in multilayer graphene, we need to introduce interlayer hopping terms. We denote the two closest hoppings between sites A and B in neighboring layers by $\gamma_1$ and $\gamma_3$. The parameter $\gamma_1$ is the direct hopping term between planes, coupling sites A and B that are spatially on top of each other. Since this is a purely vertical hopping, it appears as a constant in the tight-binding Hamiltonian. In contrast, $\gamma_3$ parametrizes the nondirect interlayer hopping by coupling sites A and B that are not spatially on top of each other. Due to its finite intralayer reach, it has a k-space modulation and acts as a trigonal deformation term. The values of these parameters have previously been determined, for example, in Ref. [96] to be $\gamma_1 = 0.38$ eV and $\gamma_3 = 0.38$ eV.

We assume that the interlayer SC coupling is negligible compared to the intralayer pairing terms. Several arguments can be brought to support this assumption, for example, the interaction between the graphene layers are van der Waals forces, and therefore much smaller than the intralayer

couplings. Moreover, the in-plane interactions in each graphene layer, given by Eqs. (4) and (5), have been shown to be strong enough to make graphene almost strongly correlated [3,97]. Still, the normal state band structure is strongly influenced at the K points by $\gamma_1$ and $\gamma_3$, the interaction between the graphene layers, and an interlayer pairing has also been shown to be possible in multilayer graphene [64]. However, with the SC state only larger at higher doping levels, away from the K points, we in the following choose to limit ourselves to intralayer SC pairing terms.

To write the BdG Hamiltonian for SC bilayer and trilayer graphene, we use $\alpha \in \{1, 2, 3\}$ to describe the bottom (1), middle (2), or top (3) layers, respectively, and we write

$$c_k^\alpha = \{a_{k\uparrow}^\alpha, b_{k\uparrow}^\alpha, a_{k\downarrow}^\alpha, b_{k\downarrow}^\alpha\}. \qquad (30)$$

Thus, the bases to use for bilayer and trilayer graphene are, respectively,

$$\Psi_k^{\text{bilayer}} = \{c_k^1, c_k^2, c_{-k}^{1\dagger}, c_{-k}^{2\dagger}\}^T \text{ and} \qquad (31)$$

$$\Psi_k^{\text{trilayer}} = \{c_k^1, c_k^2, c_k^3, c_{-k}^{1\dagger}, c_{-k}^{2\dagger}, c_{-k}^{3\dagger}\}^T. \qquad (32)$$

Below we treat each case separately.
<figure><img src="image_11.png" /><figcaption>FIG. 7. DOS as a function of chemical potential $\mu$ for SC amplitude $\Delta_0 = 0.4$ for SC states with ON s-wave, NN $s_{ext}^-$, $d_{xy}^-$, $d_{x^2-y^2}^-$, $p_{x^-}$, $p_{y^-}$, $p + ip'$, $d + id'$-wave, as well as NNN $f_{x(x^2-3y^2)}$-wave symmetries. Dark blue represents zero DOS.</figcaption></figure>

A. Bilayer Hamiltonian

The BdG Hamiltonian for AB bilayer graphene can be expressed as

$$H_{\text{AB}} = \begin{pmatrix} H_0^{\text{AB}} & H_{\text{SC}} \\ H_{\text{SC}}^{\dagger} & -H_0^{\text{AB}} \end{pmatrix}, \qquad (33)$$

<figure><img src="image_12.png" /><figcaption>FIG. 8. DOS as a function of chemical potential $\mu$ and SC amplitude $\Delta_0$ for the NN $p + ip'$-wave state. Dark blue represents zero DOS and illustrates the gap closing points.</figcaption></figure>

with

$$H_0^{\text{AB}} = \begin{pmatrix} H_1 & H_{12} \\ H_{12}^{\dagger} & H_2 \end{pmatrix}, \quad H_{\text{SC}} = \begin{pmatrix} H_{\text{SC}}^1 & 0 \\ 0 & H_{\text{SC}}^2 \end{pmatrix}. \qquad (34)$$

The components of the normal-state tight-binding Hamiltonian $H_0^{\text{AB}}$ are given by $H_1 = H_2 = \mathbf{1}_{2 \times 2} \otimes H_0(\mathbf{k})$ which describe the intralayer hoppings in each layer 1 and 2, while the interlayer terms are given by [77,78],

$$H_{12} = \mathbf{1}_{2 \times 2} \otimes H_{\perp}(\mathbf{k}), \quad H_{\perp} = \begin{pmatrix} 0 & \gamma_3 \tilde{h}_0^*(\mathbf{k}) \\ \gamma_1 & 0 \end{pmatrix}, \qquad (35)$$

where $H_0(\mathbf{k})$ and $h_0(\mathbf{k})$ are, respectively, given by Eqs. (2) and (3) in Sec. IIC. For details of the physics of each element of $H_\perp$ in relation to the normal state of bilayer graphene, see Refs. [77,78].

The SC state is represented through $H_{\text{SC}}^{1,2}$, which are $4 \times 4$ matrices describing the intralayer SC pairing in layers 1 and 2, respectively. Due to the absence of interlayer SC pairing, $H_{\text{SC}}$ is diagonal in the layer index. Furthermore, making the natural assumption of having the same SC mechanism in each layer, we consider that both layers have the same SC symmetry, leading to $h_{\text{SC}} = \text{SC}^1 = H_{\text{SC}}^2$, and $h_{\text{SC}}$ is identical to the SC monolayer graphene matrix derived in Sec. IIC.
<figure><img src="image_13.png" /><figcaption>FIG. 9. Lowest energy band structure for SC amplitude $\Delta_0 = 0.4$ at three different values of the chemical potential, $\mu = 0.4$ (a), 0.9 (b), and 1.4 (c) for the NN $p + ip'$-wave state. Only the middle panel (b) hosts zero-energy states.</figcaption></figure>

## B. Trilayer Hamiltonian

For ABA/ABC trilayer graphene, the BdG Hamiltonian is given by

$$H_{\text{ABA/ABC}} = \begin{pmatrix} H_0^{\text{ABA/ABC}} & H_{\text{SC}} \\ H_{\text{SC}}^{\dagger} & -H_0^{\text{ABA/ABC}} \end{pmatrix}, \quad (36)$$

where $H_{\text{SC}}$ is a block diagonal matrix with the three blocks given by ($H_{\text{SC}}^1, H_{\text{SC}}^2, H_{\text{SC}}^3$). As in the bilayer case, $H_{\text{SC}}$ does not couple the different layers, such that the matrix is diagonal in the layer space and is also characterized for each layer by the same SC Hamiltonian, $H_{\text{SC}}^1 = H_{\text{SC}}^2 = H_{\text{SC}}^3$. For the normal-state Hamiltonian, we have [77]

$$H_0^{\text{ABC/ABA}} = \begin{pmatrix} H_1 & H_{12} & 0 \\ H_{12}^{\dagger} & H_2 & H_{23}/H_{23}^{T} \\ 0 & H_{23}^{\dagger}/(H_{23}^{T})^{\dagger} & H_3 \end{pmatrix}. \quad (37)$$

<figure><img src="image_14.png" /><figcaption>FIG. 10. Illustration of bilayer and trilayer graphene adapted from Ref. [95]. (a) monolayer graphene. (b) AB or Bernal-stacked bilayer graphene. (c) ABA or Bernal-stacked trilayer graphene. (d) ABC, or rhombohedral-stacked, trilayer graphene. Red arrows correspond to intralayer hopping $t$, magenta to interlayer hopping $\gamma_1$, and cyan to trigonal warping terms $\gamma_3$.</figcaption></figure>

Here the intralayer terms are identical, $H_1 = H_2 = H_3$ and are given by Eq. (3) above. Furthermore, for the ABC trilayer $H_{23} = H_{12}$, with $H_{12}$ given by the bilayer expression in Eq. (35). This is the consequence of having the same interlayer coupling between the bottom and middle layers as between the middle and top layers. In contrast, for the ABA trilayer, the type of pairs of sites connecting the top two layers is opposite to that connecting the bottom two layers: $\gamma_1$ connects a site A in the middle layer to a site B in the top layer, while $\gamma_3$ connects a site B in the middle layer to a site A in the top layer. This swap in the roles of $\gamma_1$ and $\gamma_3$ in the interlayer Hamiltonian is captured by the transpose operation $T$ on $H_{23}$ in the ABA Hamiltonian.

## C. Lowest energy bands

Having derived the Hamiltonians for bilayer and trilayer graphene, we continue with exact diagonalization to find the energy spectrum. In the following, we focus mainly on bilayer graphene, aiming to extract the generic differences between the monolayer and any multilayer band structures. First, we describe briefly the effect of trigonal warping in the normal state. To this end, we plot in Fig. 11 the lowest energy band in bilayer graphene for zero chemical potential and vanishing SC pairing. When the trigonal warping $\gamma_3$ is set to zero, we obtain at very low energy (dark blue region) a quadratic and isotropic dispersion relation with circular isoenergy contours around the $K, K'$ Dirac points [78]. When including also the

<figure><img src="image_15.png" /><figcaption>FIG. 11. Lowest energy band for undoped bilayer graphene ($\mu = 0$) in the normal state ($\Delta_0 = 0$) when $\gamma_3 = 0$ (a) and $\gamma_3 = 0.2$ (b). Trigonal warping gives rise to deformations in the lowest-energy isoenergy contours.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>FIG. 12. Lowest energy band structure for $\mu = 0.4$, $\Delta_0 = 0.4$, $\gamma_1 = 0.2$, $\gamma_3 = 0$, and for SC states with ON s-wave, NN $s_{ext}^-$, $d_{xy}^-$, $d_{x^2-y^2}^-$, $p_x^-$, $p_y^-$, $p + ip'$-, and $d + id'$-wave, as well as $f_{x(x^2-3y^2)}$-NNN symmetries. Only the region in the Brillouin zone close to the Dirac point $\mathcal{K}$ is displayed. Dark blue corresponds to zero energy. The inset of (b) presents the energy dispersion along a line cut across the K point, as illustrated by the black dashed line. It shows that $s_{ext}$-wave pairing is fully gapped.</figcaption></figure>

interlayer $\gamma_3$ coupling, we note that the behavior of the isoen-
ergy contours around the $\mathcal{K}$, $\mathcal{K}'$ points is significantly altered
and acquires a trigonal distortion also at the lowest energies,
hence the name trigonal warping term for $\gamma_3$.

Moving on to the SC states, we plot in Fig. 12 the lowest
energy bands for the same SC symmetries as in monolayer
graphene (cf. Fig. 4), but now zoomed in around the Dirac
point $\mathcal{K}$ to capture most clearly the effects of the interlayer
coupling. We here first use $\gamma_3 = 0$. We find that for the fully
gapped states, the lowest energy bands are similar to the
monolayer case. However, for all the nodal states, we find a
significant difference from the monolayer, in that all the nodal
points double by exhibiting a small splitting in momentum
space; see dark blue color gradient. Thus, the nodal $d$- and
$p$-wave states now have four nodal points, two on each side
of the Dirac point. Next, we turn on the trigonal warping
and find that it has a significant effect on the formation of the
nodes. In Fig. 13, we plot the lowest band as a function of $k_y$
for $k_x = 0$ for the illustrative case of NN $p_y$-wave symmetry.

We note that the two nodal points in monolayer graphene split
and give rise to four nodal points for bilayer graphene with no
trigonal warping. When turning on $\gamma_3$, we find that this term
slightly gaps two of the four nodes.

D. Density of states and gap closing points

It is also interesting to study the DOS in bilayer and tri-
layer graphene, in particular, elucidating the existence of gap
closing points. Using a similar procedure as in the previous
section, we calculate the DOS as a function of energy, chem-
ical potential, and SC order parameter amplitude. We find
that it has a very similar dependence to that of monolayer
graphene. For completeness, we provide plots of the DOS
for the most interesting case, the $p + ip'$-wave state, in Ap-
pendix F for both bilayer and trilayer graphene. The main
features observed for the monolayer are preserved, except
that the number of gap closing points is different for multi-
layer graphene. Focusing on these gap closing points for the
<figure><img src="image_17.png" /><figcaption>FIG. 13. Lowest energy band structure for the NN $p_y$-wave SC state as a function of $k_y$ for $k_x = 0$. We compare monolayer graphene (a), bilayer graphene with no trigonal warping, $\gamma_1 = 0.2$ and $\gamma_3 = 0$ (b), and bilayer graphene with finite trigonal warping, $\gamma_1 = 0.2$ and $\gamma_3 = 0.2$ (c). Here $\Delta_0 = 0.4$ and $\mu = 0.4$.</figcaption></figure>

chiral $p + ip'$-wave SC state, we plot in Fig. 14 the gap in the energy spectrum as a function of the chemical potential and the SC amplitude for bilayer graphene without trigonal warping ($\gamma_3 = 0$). We find that there are now two gap closing lines for each $\Delta_0$ as a function of $\mu$, compared to the single gap closing line in monolayer graphene in Fig. 8. Adding a finite $\gamma_3$ changes this, and we find again only one gap closure line. Thus, the trigonal warping has a strong influence on the number of gap closing points. Indeed, this is confirmed in the right panel in Fig. 14, where we plot the DOS as a function of both $\mu$ and $\gamma_3$. We find a similar effect in trilayer graphene, with the number of gap closing points oscillating between one and three. A more detailed analysis is presented in Appendix F.

Finally, we note the formation of an interesting additional gap closing point at small $\mu$ and $\Delta_0$, not present in the monolayer case; see the dark blue line close to $\mu \approx 0.2$ in Fig. 14. Since any realistic system will have a rather small $\Delta_0$, this is particularly alluring from an experimental point of view. It appears that this gap closing point is not overly sensitive to the value of $\gamma_3$. In Fig. 15, we plot the lowest energy band for the parameters corresponding to this gap closing point and establish that this interlayer-induced gap closing occurs at the $K, K'$ Dirac points, in contrast to the previously described gap closing points that occur due to a nodal intraband pairing at the $M$ points. We attribute this gap closing to the combination of nodal points at the $K, K'$ points in the intraband SC order parameter and the bottom of the second graphene band touching the $K, K'$ points when the chemical potential $\mu$ is equal to the interlayer coupling $\gamma_1 = 0.2$. As such, it is intimately tied to the multilayer aspect.

## IV. CONCLUSIONS

In this paper, we reviewed all symmetry-allowed spin-singlet and spin-triplet SC states in monolayer, bilayer, and trilayer graphene with different stacking. By allowing ON pairing, as well as pairing between NNs and NNNs, we captured all possible spin-singlet and -triplet order parameters, from s- and d-wave to p- and f-wave states, including the chiral $d + id'-$ and $p + ip'$-wave states. To analyze the properties of these states, we calculated the low-energy band structure, as well as the DOS as a function of the chemical potential and the SC pairing strength. The different SC states can be classified in two large classes, the fully gapped states, such as the $s-$, $f_{x(x^2-3y^2)-}$, $d + id'-$, and $p + ip'$-wave states, that exhibit a U-shaped DOS and the nodal ones, such as the $d_{xy}$-, $d_{x^2-y^2}$-, $p_x-$, and $p_y$-wave states that exhibit nodal points in the band dispersion and a V-shaped DOS. Moreover, we focused on the existence of gap closing points in the DOS when changing the physical parameters and found that many of these points can be understood by a careful examination of the symmetry of the intraband SC order parameter in monolayer graphene. For bilayer and trilayer graphene, we found that the interlayer coupling splits the nodal points in the band structure, as well as the gap closing points. We also analyzed the effect of the trigonal warping present in bilayer and trilayer graphene on the formation of gap closing points. By distinguishing between nodal and fully gapped SC states, our paper provides an experimentally viable tool to differentiate between different SC symmetries in all carbon-based superconductors. In future works, we will discuss the topology and formation of edge states in these systems, as

<figure><img src="image_18.png" /><figcaption>FIG. 14. DOS as a function of chemical potential $\mu$ and SC amplitude $\Delta_0$ for the NN $p + ip'$-wave state in bilayer graphene for $\gamma_3 = 0$ (a), $\gamma_3 = 0.2$ (b), with $\gamma_1 = 0.2$. (c) DOS as a function of $\gamma_3$ and $\mu$ for $\Delta_0 = 0.4$ and $\gamma_1 = 0.2$. Dark blue corresponds to a vanishing DOS and to the gap closing points.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>FIG. 15. Lowest energy band for bilayer graphene for the NN $p + ip'$-wave SC state for $\mu = 0.1$, $\Delta_0 = 0.1$, $\gamma_1 = 0.2$, and $\gamma_3 = 0$. The dark blue regions correspond to zero energy.</figcaption></figure>

well as the formation of Shiba states, in the hope to provide additional tools to further distinguish experimentally between various order parameters arising in SC graphene materials.

## ACKNOWLEDGMENTS

We thank M. Alvarado and T. Löthman for useful discussions. We acknowledge financial support from the Swedish Research Council (Vetenskapsrådet Grant No. 2018-03488) and the Knut and Alice Wallenberg Foundation through the Wallenberg Academy Fellows program.

# APPENDIX A: DERIVATION OF THE MOMENTUM SPACE HAMILTONIAN

In this Appendix, we provide additional details for the derivation of the BdG Hamiltonian in momentum space.

## 1. Normal state Hamiltonian

For the tight-binding term in momentum space, we have

$$H_{\text{TB}}(k) = \sum_{k,\sigma} (h_0(k)a_{k,\sigma}^{\dagger} b_{k,\sigma} + h_0^{*}(k)b_{k,\sigma}^{\dagger} a_{k,\sigma}). \quad (A1)$$

In the BdG form, and also taking into account the doubling of the number of modes to include separately the electron and the hole spectrum for both up and down spins, this becomes

$$H_{\text{TB}}(k) = \sum_{k,\sigma} h_0(k)(a_{k\sigma}^{\dagger}b_{k\sigma} - b_{-k,\sigma}^{\dagger}a_{-k,\sigma}) + \sum_{k,\sigma} h_0^{*}(k)(b_{k\sigma}^{\dagger}a_{k\sigma} - a_{-k,\sigma}b_{-k,\sigma}^{\dagger}), \quad (A2)$$

with a form factor

$$h_0(k) = -t \left[ e^{-ik_y} + 2e^{i\frac{k_y}{2}} \cos \left( \frac{\sqrt{3}}{2} k_x \right) \right]. \quad (A3)$$

For the chemical potential term, we have

$$H_{\mu} = -\mu \sum_{k\sigma} (a_{k,\sigma}^{\dagger} a_{k,\sigma} + b_{k,\sigma}^{\dagger} b_{k,\sigma}). \quad (A4)$$

In the BdG form, this becomes

$$H_{\mu} = -\mu \sum_{k\sigma} (a_{k,\sigma}^{\dagger} a_{k,\sigma} - a_{-k,\sigma} a_{-k,\sigma}^{\dagger} + b_{k,\sigma}^{\dagger} b_{k,\sigma} - b_{-k,\sigma}^{\dagger} b_{-k,\sigma}^{\dagger}). \quad (A5)$$

## 2. On-site pairing

In real space, the ON pairing term takes the following form:

$$H_{\text{ON}} = 2 \left[ \sum_{i} \Delta_{\text{ON}} a_{i\uparrow}^{\dagger} a_{i\downarrow}^{\dagger} + \sum_{i} \Delta_{\text{ON}} b_{i\uparrow}^{\dagger} b_{i\downarrow}^{\dagger} + \text{H.c.} \right], \quad (A6)$$

where the factor of 2 come from the doubling of the number of modes to take into account separately the electrons/holes with spin up and spin down. The Fourier transform of this term is

$$H_{\text{ON}} = 2 \left[ \sum_{k\sigma} \Delta_{\text{ON}} a_{k\uparrow}^{\dagger} a_{-k\downarrow}^{\dagger} + \sum_{k\sigma} \Delta_{\text{ON}} b_{k\uparrow}^{\dagger} b_{-k\downarrow}^{\dagger} + \text{H.c.} \right], \quad (A7)$$

which can be written in the BdG form as

$$H_{\text{ON}} = \sum_{k\sigma} \Delta_{\text{ON}} (a_{k\uparrow}^{\dagger} a_{-k\downarrow}^{\dagger} - a_{-k\downarrow}^{\dagger} a_{k\uparrow}^{\dagger}) + \sum_{k\sigma} \Delta_{\text{ON}} (b_{k\uparrow}^{\dagger} b_{-k\downarrow}^{\dagger} - b_{-k\downarrow}^{\dagger} b_{k\uparrow}^{\dagger}) + \text{H.c.} \quad (A8)$$

## 3. NN pairing

For NN pairing, in the spin-singlet channel ($\eta = 0$) we have

$$H_{\text{NN}}^{0} = \sum_{\langle i,j \rangle} \Delta_{ij}^{\eta=0} (a_{i\uparrow}^{\dagger} b_{j\downarrow}^{\dagger} - a_{i\downarrow}^{\dagger} b_{j\uparrow}^{\dagger}) + \text{H.c.} \quad (A9)$$

while in the spin-triplet channel we find

$$H_{\text{NN}}^{x} = \sum_{\langle i,j \rangle} \Delta_{ij}^{\eta=x} (a_{i\uparrow}^{\dagger} b_{j\uparrow}^{\dagger} - a_{i\downarrow}^{\dagger} b_{j\downarrow}^{\dagger}) + \text{H.c.}, \quad (A10)$$

$$H_{\text{NN}}^{y} = \sum_{\langle i,j \rangle} \Delta_{ij}^{\eta=y} (a_{i\uparrow}^{\dagger} b_{j\uparrow}^{\dagger} + a_{i\downarrow}^{\dagger} b_{j\downarrow}^{\dagger}) + \text{H.c.}, \quad (A11)$$

$$H_{\text{NN}}^{z} = \sum_{\langle i,j \rangle} \Delta_{ij}^{\eta=z} (a_{i\uparrow}^{\dagger} b_{j\downarrow}^{\dagger} + a_{i\downarrow}^{\dagger} b_{j\uparrow}^{\dagger}) + \text{H.c.} \quad (A12)$$

We can now perform the Fourier transform as before. In the spin-singlet channel, we obtain

$$H_{\text{NN}} = \sum_{k} h_{\text{NN}}^{0} (a_{k\uparrow}^{\dagger} b_{-k\downarrow}^{\dagger} - a_{k\downarrow}^{\dagger} b_{-k\uparrow}^{\dagger}) + \text{H.c.}, \quad (A13)$$

and in the spin-triplet channel

$$H_{\text{NN}}^{x} = \sum_{k} h_{\text{NN}}^{y} (a_{k\uparrow}^{\dagger} b_{-k\uparrow}^{\dagger} - a_{k\downarrow}^{\dagger} b_{-k\downarrow}^{\dagger}) + \text{H.c.}, \quad (A14)$$

$$H_{\text{NN}}^{y} = i \sum_{k} h_{\text{NN}}^{z} (a_{k\uparrow}^{\dagger} b_{-k\uparrow}^{\dagger} + a_{k\downarrow}^{\dagger} b_{-k\downarrow}^{\dagger}) + \text{H.c.}, \quad (A15)$$

$$H_{\text{NN}}^{z} = \sum_{k} h_{\text{NN}}^{y} (a_{k\uparrow}^{\dagger} b_{-k\downarrow}^{\dagger} + a_{k\downarrow}^{\dagger} b_{-k\uparrow}^{\dagger}) + \text{H.c.} \quad (A16)$$
Here $h_{NN}^{\eta}$ are the form factors, whose expressions depend on both the spin channel and the symmetry of the order parameter. Their general expression is

$$h_{NN}^{\eta} = \Delta_{NN}^{\eta,d=1} e^{-ik_y} + \Delta_{NN}^{\eta,d=2} e^{\frac{i}{2}k_y - \frac{\sqrt{3}i}{2}k_x} + \Delta_{NN}^{\eta,d=3} e^{\frac{i}{2}k_y + \frac{\sqrt{3}i}{2}k_x}. \quad (A17)$$

Here, $d=1, 2, 3$ correspond to the three NNs, following the convention of Fig. 1. The values of $\Delta^{\eta,d}$ are also detailed in Tables I and II for each symmetry. By replacing these values, we obtain the form factors for NN pairing ($h_{NN}^{0,sext}$, $h_{NN}^{0,dxy}$, $h_{NN}^{0,d_{xy-2}^2}$, $h_{NN}^{\eta,p_x}$, $h_{NN}^{\eta,p_y}$, $h_{NN}^{\eta,f_x}$) in Tables III and IV.

We subsequently express the above NN pairing term in the BdG form

$$ \begin{aligned} H_{NN}^{0} &= \frac{1}{2} \sum_{k} h_{NN}^{0}(k) (a_{k\uparrow}^{\dagger} b_{-k\downarrow}^{\dagger} - b_{-k\downarrow}^{\dagger} a_{k\uparrow}^{\dagger} - a_{k\downarrow}^{\dagger} b_{-k\uparrow}^{\dagger} \\ &\quad + b_{-k\uparrow}^{\dagger} a_{k\downarrow}^{\dagger}) + \text{H.c.}, \\ &= \frac{1}{2} \sum_{k} (h_{NN}^{0}(k) a_{k\uparrow}^{\dagger} b_{-k\downarrow}^{\dagger} - h_{NN}^{0}(-k) b_{k\downarrow}^{\dagger} a_{-k\uparrow}^{\dagger} \\ &\quad - h_{NN}^{0}(k) a_{k\downarrow}^{\dagger} b_{-k\uparrow}^{\dagger} + h_{NN}^{0}(-k) b_{k\uparrow}^{\dagger} a_{-k\downarrow}^{\dagger}) + \text{H.c.} \end{aligned} \quad (A18) $$

for the singlet channel and

$$ H_{NN}^{x} = \frac{1}{2} \sum_{k} (h_{NN}^{x}(k) a_{k\uparrow}^{\dagger} b_{-k\downarrow}^{\dagger} - h_{NN}^{x}(-k) b_{k\uparrow}^{\dagger} a_{-k\downarrow}^{\dagger} \\ - h_{NN}^{x}(k) a_{k\downarrow}^{\dagger} b_{-k\downarrow}^{\dagger} + h_{NN}^{x}(-k) b_{k\downarrow}^{\dagger} a_{-k\downarrow}^{\dagger}) + \text{H.c.}, $$

$$ H_{NN}^{y} = \frac{i}{2} \sum_{k} (h_{NN}^{y}(k) a_{k\uparrow}^{\dagger} b_{-k\downarrow}^{\dagger} - h_{NN}^{y}(-k) b_{k\uparrow}^{\dagger} a_{-k\downarrow}^{\dagger} \\ + h_{NN}^{y}(k) a_{k\downarrow}^{\dagger} b_{-k\downarrow}^{\dagger} - h_{NN}^{y}(-k) b_{k\downarrow}^{\dagger} a_{-k\downarrow}^{\dagger}) + \text{H.c.}, $$

$$ H_{NN}^{z} = \frac{1}{2} \sum_{k} (h_{NN}^{z}(k) a_{k\uparrow}^{\dagger} b_{-k\downarrow}^{\dagger} - h_{NN}^{z}(-k) b_{k\downarrow}^{\dagger} a_{-k\uparrow}^{\dagger} \\ + h_{NN}^{z}(k) a_{k\downarrow}^{\dagger} b_{-k\uparrow}^{\dagger} - h_{NN}^{z}(-k) b_{k\uparrow}^{\dagger} a_{-k\downarrow}^{\dagger}) + \text{H.c.} \quad (A19) $$

for the triplet channel.

### 4. NNN pairing

For NNN pairing in the spin-singlet channel ($\eta=0$), we have

$$ H_{NNN}^{0} = \sum_{\langle\langle i, j \rangle\rangle} \Delta_{ij}^{\eta=0} (a_{i\uparrow}^{\dagger} a_{j\downarrow}^{\dagger} - a_{i\downarrow}^{\dagger} a_{j\uparrow}^{\dagger} \\ + b_{i\uparrow}^{\dagger} b_{j\downarrow}^{\dagger} - b_{i\downarrow}^{\dagger} b_{j\uparrow}^{\dagger}) + \text{H.c.}, \quad (A20) $$

and in the spin-triplet channel

$$ H_{NNN}^{x} = \sum_{\langle\langle i, j \rangle\rangle} \Delta_{ij}^{\eta=0} (a_{i\uparrow}^{\dagger} a_{j\uparrow}^{\dagger} - a_{i\downarrow}^{\dagger} a_{j\downarrow}^{\dagger} \\ + b_{i\uparrow}^{\dagger} b_{j\uparrow}^{\dagger} - b_{i\downarrow}^{\dagger} b_{j\downarrow}^{\dagger}) + \text{H.c.}, \quad (A21) $$

$$ H_{NNN}^{y} = \sum_{\langle\langle i, j \rangle\rangle} \Delta_{ij}^{\eta=0} (a_{i\uparrow}^{\dagger} a_{j\uparrow}^{\dagger} + a_{i\downarrow}^{\dagger} a_{j\downarrow}^{\dagger} + b_{i\uparrow}^{\dagger} b_{j\uparrow}^{\dagger} \\ + b_{i\downarrow}^{\dagger} b_{j\downarrow}^{\dagger}) + \text{H.c.}, \quad (A22) $$

$$ H_{NNN}^{0} = \sum_{\langle\langle i, j \rangle\rangle} \Delta_{ij}^{\eta=0} (a_{i\uparrow}^{\dagger} a_{j\downarrow}^{\dagger} + a_{i\downarrow}^{\dagger} a_{j\uparrow}^{\dagger} \\ + b_{i\uparrow}^{\dagger} b_{j\downarrow}^{\dagger} + b_{i\downarrow}^{\dagger} b_{j\uparrow}^{\dagger}) + \text{H.c.}, \quad (A23) $$

where the sums are now performed over the NNN pairs. Exactly as before, performing the Fourier transform gives

$$ H_{NNN}^{0} = \frac{1}{2} \sum_{k} h_{NNN}^{\eta=0}(k) (a_{k\uparrow}^{\dagger} a_{-k\downarrow}^{\dagger} - a_{k\downarrow}^{\dagger} a_{-k\uparrow}^{\dagger} \\ + b_{k\uparrow}^{\dagger} b_{-k\downarrow}^{\dagger} - b_{k\downarrow}^{\dagger} b_{-k\uparrow}^{\dagger}) + \text{H.c.} \quad (A24) $$

in the spin-singlet channel and

$$ H_{NNN}^{x} = \frac{1}{2} \sum_{k} h_{NNN}^{x}(k) (a_{k\uparrow}^{\dagger} a_{-k\uparrow}^{\dagger} - a_{k\downarrow}^{\dagger} a_{-k\downarrow}^{\dagger} \\ + b_{k\uparrow}^{\dagger} b_{-k\uparrow}^{\dagger} - b_{k\downarrow}^{\dagger} b_{-k\downarrow}^{\dagger}) + \text{H.c.}, \quad (A25a) $$

$$ H_{NNN}^{y} = \frac{i}{2} \sum_{k} h_{NNN}^{y}(k) (a_{k\uparrow}^{\dagger} a_{-k\uparrow}^{\dagger} + a_{k\downarrow}^{\dagger} a_{-k\downarrow}^{\dagger} \\ + b_{k\uparrow}^{\dagger} b_{-k\uparrow}^{\dagger} + b_{k\downarrow}^{\dagger} b_{-k\downarrow}^{\dagger}) + \text{H.c.}, \quad (A25b) $$

$$ H_{NNN}^{z} = \frac{1}{2} \sum_{k} h_{NNN}^{z}(k) (a_{k\uparrow}^{\dagger} a_{-k\downarrow}^{\dagger} + a_{k\downarrow}^{\dagger} a_{-k\uparrow}^{\dagger} \\ + b_{k\uparrow}^{\dagger} b_{-k\downarrow}^{\dagger} + b_{k\downarrow}^{\dagger} b_{-k\uparrow}^{\dagger}) + \text{H.c.}, \quad (A25c) $$

in the spin-triplet one. Note that we have a $\frac{1}{2}$ factor appearing since one takes twice into account the links between the NNNs when summing over all the atoms. Following the NNN conventions of Fig. 1, we write

$$ \begin{aligned} h_{NNN}^{\eta} = \Delta_{NNN}^{\eta,d=1} e^{-\frac{\sqrt{3}i}{2}k_x + \frac{3i}{2}k_y} + \Delta_{NNN}^{\eta,d=2} e^{-i\sqrt{3}k_x} \\ + \Delta_{NNN}^{\eta,d=3} e^{-\frac{\sqrt{3}i}{2}k_x - \frac{3i}{2}k_y} + \Delta_{NNN}^{\eta,d=4} e^{\frac{\sqrt{3}i}{2}k_x - \frac{3i}{2}k_y} \\ + \Delta_{NNN}^{\eta,d=5} e^{i\sqrt{3}k_x} + \Delta_{NNN}^{\eta,d=6} e^{\frac{\sqrt{3}i}{2}k_x + \frac{3i}{2}k_y}. \end{aligned} \quad (A26) $$

The explicit expressions of the resulting NNN form factors ($h_{NNN}^{0,sext}$, $h_{NNN}^{0,dxy}$, $h_{NNN}^{0,d_{xy-2}^2}$, $h_{NNN}^{\eta,p_x}$, $h_{NNN}^{\eta,p_y}$, $h_{NNN}^{\eta,f_x}$) are given in Tables III and IV.

## APPENDIX B: SPIN-TRILET PAIRING HAMILTONIAN

In this Appendix, we derive the Hamiltonian for spin-triplet pairing with $\eta = y, z$. The result for $\eta = x$ is given in the main text. For the $\eta = y$ spin-triplet channel, the Hamiltonian reads

$$ \begin{pmatrix} H_0(\mathbf{k}) & \mathbf{0}_{2 \times 2} & h_{\Delta}^{y}(\mathbf{k}) & \mathbf{0}_{2 \times 2} \\ \mathbf{0}_{2 \times 2} & H_0(\mathbf{k}) & \mathbf{0}_{2 \times 2} & h_{\Delta}^{y}(\mathbf{k}) \\ (h_{\Delta}^{y}(\mathbf{k}))^{\dagger} & \mathbf{0}_{2 \times 2} & -H_0(\mathbf{k}) & \mathbf{0}_{2 \times 2} \\ \mathbf{0}_{2 \times 2} & (h_{\Delta}^{y}(\mathbf{k}))^{\dagger} & \mathbf{0}_{2 \times 2} & -H_0(\mathbf{k}) \end{pmatrix}, \\ h_{\Delta}^{y}(\mathbf{k}) = \frac{i}{2} \begin{pmatrix} -h_{NNN}^{y}(-\mathbf{k}) & h_{NN}^{y}(\mathbf{k}) \\ -h_{NN}^{y}(-\mathbf{k}) & -h_{NNN}^{y}(-\mathbf{k}) \end{pmatrix}, \quad (B1) $$
<figure><img src="image_20.png" /><figcaption>FIG. 16. Lowest energy band structure for $\mu = 1$, $\Delta_0 = 0.4$, and for SC states with ON s-wave, NN $s_{ext}$, $d_{xy}$, $d_{x^2-y^2}$, $d + id'$, $p_x$, $p_y$, $p + ip'$-wave, as well as NNN $f_{x(x^2-3y^2)}$ symmetries. Dark blue corresponds to zero energy.</figcaption></figure>

while for the $\eta = z$ spin-triplet channel we find

$$ \begin{pmatrix} H_0(\mathbf{k}) & \mathbf{0}_{2 \times 2} & \mathbf{0}_{2 \times 2} & -h_{\Delta}^{z}(\mathbf{k}) \\ \mathbf{0}_{2 \times 2} & H_0(\mathbf{k}) & h_{\Delta}^{z}(\mathbf{k}) & \mathbf{0}_{2 \times 2} \\ \mathbf{0}_{2 \times 2} & (h_{\Delta}^{z}(\mathbf{k}))^{\dagger} & -H_0(\mathbf{k}) & \mathbf{0}_{2 \times 2} \\ -(h_{\Delta}^{z}(\mathbf{k}))^{\dagger} & \mathbf{0}_{2 \times 2} & \mathbf{0}_{2 \times 2} & -H_0(\mathbf{k}) \end{pmatrix}, \qquad h_{\Delta}^{z}(\mathbf{k}) = \frac{1}{2} \begin{pmatrix} h_{\text{NNN}}^{z}(-\mathbf{k}) & -h_{\text{NN}}^{z}(\mathbf{k}) \\ h_{\text{NN}}^{z}(-\mathbf{k}) & h_{\text{NNN}}^{z}(-\mathbf{k}) \end{pmatrix}, \quad (B2) $$

where $\mathbf{0}_{2 \times 2}$ is a 2 × 2 null matrix, $H_0(\mathbf{k})$ is the normal state Hamiltonian matrix given in Eq. (2), and $h_{\Delta}^y, h_{\Delta}^z$ are the form factor for spin-triplet SC symmetry.

## APPENDIX C: LOWEST ENERGY BAND FOR $\mu = 1$

In this Appendix, we provide plots of the lowest energy bands at $\mu = 1$ in Fig. 16. This is to be compared to the results for $\mu = 0.4$ given in Fig. 4 in the main text. For a doping close to the VHS, $\mu = 1$, we thus

plot the lowest energy band as a function of $k_x$ and $k_y$. We note that in this case the bands for some of the SC symmetries reach zero energy in the vicinity of the M points (corresponding to the dark blue regions). This gives rise to a gap closing point in the DOS around the value of $\mu = 1$.

## APPENDIX D: BAND STRUCTURE FOR NNN PAIRING

In this Appendix, we provide plots of the lowest energy bands for NNN pairing in Fig. 17. This is to be compared to results for NN pairing given in Fig. 4 in the main text. We thus consider both singlet and triplet NNN pairing with various symmetries. We note similarities with the bands resulting from NN pairing Fig. 4, such as the nodal points in the structure for the $d_{xy}$, $d_{x^2-y^2}$, $p_x$, and $p_y$-wave (doped in dark blue) and a fully gapped structure for the rest.
<figure><img src="image_21.png" /><figcaption>FIG. 17. Lowest energy band structure for $\mu = 0.4$, $\Delta_0 = 0.4$, and for SC states with NNN pairing with $s_{ext-}$, $d_{xy-}$, $d_{x^2-y^2-}$, $d + id'$, $p_x-$, $p_y-$, $p + ip'$, and $f_{x(x^2-3y^2)-}$wave symmetries. Dark blue corresponds to zero energy.</figcaption></figure>

## APPENDIX E: DOS FOR NNN PAIRING

In this Appendix, we provide information on the DOS for NNN pairing. First we plot the DOS for NNN pairing as a function of energy and SC pairing strength in Fig. 18 for the same parameter values as those considered for the NN pairing described in the main text and in Fig. 6. We note that the main difference from the NN coupling case is that the dependence with $\Delta_0$ is roughly linear, and in contrast to the NN case, there is no critical value of $\Delta_0$ for which we observe additional gap closings, for any of the symmetries. The dependence of the DOS with energy and the chemical potential is very similar to that for the NN couplings as seen when comparing Fig. 19 with Fig. 7 in the main text. Furthermore, as mentioned in the main text, the gap closing for the NNN $p + ip'$-wave state occurs at $\mu = 1$, independent of the value of the SC coupling, as seen in Fig. 20.

## APPENDIX F: DOS FOR NN $p + ip'$-WAVE PAIRING IN MULTILAYER GRAPHENE

In this Appendix, we provide additional information of the DOS for NN $p + ip'$-wave pairing in multilayer graphene. For a vanishing trigonal warping, $\gamma_3 = 0$, we obtain the DOS as a function of energy and chemical potential in Fig. 21 for different multilayer configurations. We note that these results are very similar to monolayer graphene, except for the number of gap closing points. As discussed in the main text, this number is also strongly affected by the trigonal warping, which we illustrate in Fig. 22, where we plot the DOS in the presence of non-zero trigonal warping $\gamma_3 = 0.2$.

As opposed to the monolayer case for which we have a single gap closing point as a function of the chemical potential, for trilayer graphene we find that this number oscillates between one and three. This is illustrated for ABA trilayer graphene in Fig. 23 and for ABC trilayer graphene in Fig. 24. We find that the effect of trigonal warping is quite pronounced by modifying the gap closing points.
<figure><img src="image_22.png" /><figcaption>FIG. 18. DOS as a function of SC amplitude $\Delta_0$ for $\mu = 0.4$ for SC states with NNN pairing with $s_{ext-}$, $d_{xy-}$, $d_{x^2-y^2-}$, $d+id'-$, $p_x-$, $p_y-$, $f_{x(x^2-3y^2)-}$ and $p+ip'$-wave symmetries.</figcaption></figure>
<figure><img src="image_23.png" /><figcaption>FIG. 19. DOS as a function of $\mu$ for $\Delta_0 = 0.4$ for SC states with NNN pairing with $s_{ext}^-$, $d_{xy}^-$, $d_{x^2-y^2}^-$, $d+id'$, $p_x^-$, $p_y^-$, $f_{x(x^2-3y^2)}^-$ and $p+ip'$-wave symmetries.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>FIG. 20. Gap value as a function of the chemical potential $\mu$ and the SC amplitude $\Delta_0$ for the NNN $p+ip'$-wave state.</figcaption></figure>
<figure><img src="image_25.png" /><figcaption>FIG. 21. DOS as a function of energy and $\mu$ for $\Delta_0 = 0.4$, $\gamma_1 = 0.2$, $\gamma_3 = 0$ in bilayer (a), trilayer ABA (b), and trilayer ABC (c) graphene for the NN $p + ip'$-wave SC state.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>FIG. 22. DOS as a function of energy and $\mu$ for $\Delta_0 = 0.4$, $\gamma_1 = 0.2$, $\gamma_3 = 0.2$ in bilayer (a), trilayer ABA (b), and trilayer ABC (c) graphene for the NN $p + ip'$-wave SC state.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>FIG. 23. DOS as a function of the chemical potential $\mu$ and the SC amplitude $\Delta_0$ for the NN $p + ip'$-wave state in ABA trilayer graphene for $\gamma_3 = 0$ (a) and $\gamma_3 = 0.2$ (b) with $\gamma_1 = 0.2$. (c) DOS as a function of $\gamma_3$ and $\mu$ for $\Delta_0 = 0.4$ and $\gamma_1 = 0.2$. Dark blue corresponds to a vanishing DOS and to the gap closing points.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>FIG. 24. DOS as a function of the chemical potential $\mu$ and the SC amplitude $\Delta_0$ for the NN $p + ip'$-wave state in ABC trilayer graphene for $\gamma_3 = 0$ (a) and $\gamma_3 = 0.2$ (b) with $\gamma_1 = 0.2$. (c) DOS as a function of $\gamma_3$ and $\mu$ for $\Delta_0 = 0.4$ and $\gamma_1 = 0.2$. Dark blue corresponds to a vanishing DOS and to the gap closing points.</figcaption></figure>
[1] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V. Grigorieva, and A. A. Firsov, Electric field effect in atomically thin carbon films, *Science* **306**, 666 (2004).

[2] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, The electronic properties of graphene, *Rev. Mod. Phys.* **81**, 109 (2009).

[3] A. M. Black-Schaffer and S. Doniach, Resonating valence bonds and mean-field *d*-wave superconductivity in graphite, *Phys. Rev. B* **75**, 134512 (2007).

[4] C. Honerkamp, Density waves and cooper pairing on the honeycomb lattice, *Phys. Rev. Lett.* **100**, 146404 (2008).

[5] J. Linder, A. M. Black-Schaffer, T. Yokoyama, S. Doniach, and A. Sudbø, Josephson current in graphene: Role of unconventional pairing symmetries, *Phys. Rev. B* **80**, 094522 (2009).

[6] F. M. Pellegrino, G. G. Angilella, and R. Pucci, Pairing symmetry of superconducting graphene, *Eur. Phys. J. B* **76**, 469 (2010).

[7] R. Nandkishore, L. Levitov, and A. V. Chubukov, Chiral superconductivity from repulsive interactions in doped graphene, *Nat. Phys.* **8**, 158 (2012).

[8] M. L. Kiesel, C. Platt, W. Hanke, D. A. Abanin, and R. Thomale, Competing many-body instabilities and unconventional superconductivity in graphene, *Phys. Rev. B* **86**, 020507(R) (2012).

[9] W.-S. Wang, Y.-Y. Xiang, Q.-H. Wang, F. Wang, F. Yang, and D.-H. Lee, Functional renormalization group and variational Monte Carlo studies of the electronic instabilities in graphene near $1/4$ doping, *Phys. Rev. B* **85**, 035414 (2012).

[10] A. M. Black-Schaffer and C. Honerkamp, Chiral *d*-wave superconductivity in doped graphene, *J. Phys.: Condens. Matter* **26**, 423201 (2014).

[11] A. M. Black-Schaffer, W. Wu, and K. Le Hur, Chiral *d*-wave superconductivity on the honeycomb lattice close to the Mott state, *Phys. Rev. B* **90**, 054521 (2014).

[12] T. Löthman and A. M. Black-Schaffer, Defects in the *d* + *id*-wave superconducting state in heavily doped graphene, *Phys. Rev. B* **90**, 224504 (2014).

[13] T. Ma, F. Yang, H. Yao, and H.-Q. Lin, Possible triplet *p* + *ip* superconductivity in graphene at low filling, *Phys. Rev. B* **90**, 245114 (2014).

[14] A. M. Alsharari and S. E. Ulloa, Inducing chiral superconductivity on honeycomb lattice systems, *J. Phys.: Condens. Matter* **34**, 205403 (2022).

[15] J. B. Curtis, N. R. Poniatowski, Y. Xie, A. Yacoby, E. Demler, and P. Narang, Stabilizing fluctuating spin-triplet superconductivity in graphene via induced spin-orbit coupling, *Phys. Rev. Lett.* **130**, 196001 (2023).

[16] P. Rosenzweig, H. Karachachian, S. Link, K. Muller, and U. Starke, Tuning the doping level of graphene in the vicinity of the Van Hove singularity via ytterbium intercalation, *Phys. Rev. B* **100**, 035445 (2019).

[17] P. Rosenzweig, H. Karachachian, D. Marchenko, K. Küster, and U. Starke, Overdoping graphene beyond the Van Hove singularity, *Phys. Rev. Lett.* **125**, 176403 (2020).

[18] H. B. Heersche, P. Jarillo-Herrero, J. B. Oostinga, L. M. K. Vandersypen, and A. F. Morpurgo, Bipolar supercurrent in graphene, *Nature (London)* **446**, 56 (2007).

[19] L. Bretheau, J. I.-J. Wang, R. Pisoni, K. Watanabe, T. Taniguchi, and P. Jarillo-Herrero, Tunnelling spectroscopy of Andreev states in graphene, *Nat. Phys.* **13**, 756 (2017).

[20] A. D. Bernardo, O. Millo, M. Barbone, H. Alpern, Y. Kalheim, U. Sassi, A. K. Ott, D. D. Fazio, D. Yoon, M. Amado, A. C. Ferrari, J. Linder, and J. W. A. Robinson, *p*-wave triggered superconductivity in single-layer graphene on an electron-doped oxide superconductor, *Nat. Commun.* **8**, 14024 (2017).

[21] O. A. Awoga and A. M. Black-Schaffer, Probing unconventional superconductivity in proximitized graphene by impurity scattering, *Phys. Rev. B* **97**, 214515 (2018).

[22] D. Perconte, K. Seurre, V. Humbert, C. Ulysse, A. Sander, J. Trastoy, V. Zatko, F. Godel, P. R. Kidambi, S. Hofmann, X. P. Zhang, D. Bercioux, F. S. Bergeret, B. Dlubak, P. Seneor, and J. E. Villegas, Long-range propagation and interference of *d*-wave superconducting pairs in graphene, *Phys. Rev. Lett.* **125**, 087002 (2020).

[23] D. Perconte, D. Bercioux, B. Dlubak, P. Seneor, F. S. Bergeret, and J. E. Villegas, Superconducting proximity effect in *d*-wave cuprate/graphene heterostructures, *Ann. Phys.* **534**, 2100559 (2022).

[24] H. Zhou, T. Xie, T. Taniguchi, K. Watanabe, and A. F. Young, Superconductivity in rhombohedral trilayer graphene., *Nature (London)* **598**, 434 (2021).

[25] Y. Cao, V. Fatemi, S. Fang, K. Watanabe, T. Taniguchi, E. Kaxiras, and P. Jarillo-Herrero, Unconventional superconductivity in magic-angle graphene superlattices, *Nature (London)* **556**, 43 (2018).

[26] H. C. Po, L. Zou, A. Vishwanath, and T. Senthil, Origin of Mott insulating behavior and superconductivity in twisted bilayer graphene, *Phys. Rev. X* **8**, 031089 (2018).

[27] F. Guinea and N. R. Walet, Electrostatic effects, band distortions, and superconductivity in twisted graphene bilayers, *Proc. Natl. Acad. Sci. USA* **115**, 13174 (2018).

[28] F. Wu, A. H. MacDonald, and I. Martin, Theory of phonon-mediated superconductivity in twisted bilayer graphene, *Phys. Rev. Lett.* **121**, 257001 (2018).

[29] D. M. Kennes, J. Lischner, and C. Karrasch, Strong correlations and *d* + *id* superconductivity in twisted bilayer graphene, *Phys. Rev. B* **98**, 241407(R) (2018).

[30] 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, and D. K. Efetov, Superconductors, orbital magnets and correlated states in magic-angle bilayer graphene, *Nature (London)* **574**, 653 (2019).

[31] M. Yankowitz, S. Chen, H. Polshyn, Y. Zhang, K. Watanabe, T. Taniguchi, D. Graf, A. F. Young, and C. R. Dean, Tuning superconductivity in twisted bilayer graphene, *Science* **363**, 1059 (2019).

[32] G. Chen, A. L. Sharpe, P. Gallagher, I. T. Rosen, E. J. Fox, L. Jiang, B. Lyu, H. Li, K. Watanabe, T. Taniguchi, J. Jung, Z. Shi, D. Goldhaber-Gordon, Y. Zhang, and F. Wang, Signatures of tunable superconductivity in a trilayer graphene moiré superlattice, *Nature (London)* **572**, 215 (2019).

[33] M. Alidoust, M. Willatzen, and A.-P. Jauho, Symmetry of superconducting correlations in displaced bilayers of graphene, *Phys. Rev. B* **99**, 155413 (2019).

[34] L. Balents, C. R. Dean, D. K. Efetov, and A. F. Young, Superconductivity and strong correlations in moiré flat bands, *Nat. Phys.* **16**, 725 (2020).

[35] M. Oh, K. P. Nuckolls, D. Wong, R. L. Lee, X. Liu, K. Watanabe, T. Taniguchi, and A. Yazdani, Evidence for
unconventional superconductivity in twisted bilayer graphene, *Nature* (London) **600**, 240 (2021).

[36] E. Y. Andrei and A. H. MacDonald, Graphene bilayers with a twist, *Nat. Mater.* **19**, 1265 (2020).

[37] M. Christos, S. Sachdev, and M. S. Scheurer, Superconductivity, correlated insulators, and Wess-Zumino-Witten terms in twisted bilayer graphene, *Proc. Natl. Acad. Sci. USA* **117**, 29543 (2020).

[38] D. V. Chichinadze, L. Classen, and A. V. Chubukov, Nematic superconductivity in twisted bilayer graphene, *Phys. Rev. B* **101**, 224513 (2020).

[39] Y. Cao, D. Rodan-Legrain, J. M. Park, N. F. Q. Yuan, K. Watanabe, T. Taniguchi, R. M. Fernandes, L. Fu, and P. Jarillo-Herrero, Nematicity and competing orders in superconducting magic-angle graphene, *Science* **372**, 264 (2021).

[40] X. Wu, W. Hanke, M. Fink, M. Klett, and R. Thomale, Harmonic fingerprint of unconventional superconductivity in twisted bilayer graphene, *Phys. Rev. B* **101**, 134517 (2020).

[41] T. Yu, D. M. Kennes, A. Rubio, and M. A. Sentef, Nematicity arising from a chiral superconducting ground state in magic-angle twisted bilayer graphene under in-plane magnetic fields, *Phys. Rev. Lett.* **127**, 127001 (2021).

[42] E. Khalaf, S. Chatterjee, N. Bultinck, M. P. Zaletel, and A. Vishwanath, Charged skyrmions and topological origin of superconductivity in magic-angle graphene, *Sci. Adv.* **7**, eabf5299 (2021).

[43] T. Cea and F. Guinea, Coulomb interaction, phonons, and superconductivity in twisted bilayer graphene, *Proc. Natl. Acad. Sci. USA* **118**, e2107874118 (2021).

[44] Y.-Z. Chou, F. Wu, J. D. Sau, and S. Das Sarma, Correlation-induced triplet pairing superconductivity in graphene-based moiré systems, *Phys. Rev. Lett.* **127**, 217001 (2021).

[45] A. Fischer, Z. A. H. Goodwin, A. A. Mostofi, J. Lischner, D. M. Kennes, and L. Klebl, Unconventional superconductivity in magic-angle twisted trilayer graphene, *npj Quantum Mater.* **7**, 5 (2022).

[46] H. Kim, Y. Choi, C. Lewandowski, A. Thomson, Y. Zhang, R. Polski, K. Watanabe, T. Taniguchi, J. Alicea, and S. Nadj-Perge, Evidence for unconventional superconductivity in twisted trilayer graphene, *Nature (London)* **606**, 494 (2022).

[47] T. Löthman, J. Schmidt, F. Parhizgar, and A. M. Black-Schaffer, Nematic superconductivity in magic-angle twisted bilayer graphene from atomistic modeling, *Commun. Phys.* **5**, 1 (2022).

[48] T. Cea, P. A. Pantaleón, V. o. T. Phong, and F. Guinea, Superconductivity from repulsive interactions in rhombohedral trilayer graphene: A Kohn-Luttinger-like mechanism, *Phys. Rev. B* **105**, 075432 (2022).

[49] S. Chatterjee, T. Wang, E. Berg, and M. P. Zaletel, Inter-valley coherent order and isospin fluctuation mediated superconductivity in rhombohedral trilayer graphene, *Nat. Commun.* **13**, 6013 (2022).

[50] Y.-Z. Chou, F. Wu, J. D. Sau, and S. Das Sarma, Acoustic-phonon-mediated superconductivity in rhombohedral trilayer graphene, *Phys. Rev. Lett.* **127**, 187001 (2021).

[51] A. L. Szabó and B. Roy, Metals, fractional metals, and superconductivity in rhombohedral trilayer graphene, *Phys. Rev. B* **105**, L081407 (2022).

[52] Y.-Z. You and A. Vishwanath, Kohn-luttinger superconductivity and intervalley coherence in rhombohedral trilayer graphene, *Phys. Rev. B* **105**, 134524 (2022).

[53] A. Ghazaryan, T. Holder, M. Serbyn, and E. Berg, Unconventional superconductivity in systems with annular Fermi surfaces: Application to rhombohedral trilayer graphene, *Phys. Rev. Lett.* **127**, 247001 (2021).

[54] H. Dai, R. Ma, X. Zhang, and T. Ma, Quantum Monte Carlo study of superconductivity in rhombohedral trilayer graphene under an electric field, *Phys. Rev. B* **107**, 245106 (2023).

[55] D.-C. Lu, T. Wang, S. Chatterjee, and Y.-Z. You, Correlated metals and unconventional superconductivity in rhombohedral trilayer graphene: A renormalization group analysis, *Phys. Rev. B* **106**, 155115 (2022).

[56] A. Ghazaryan, T. Holder, E. Berg, and M. Serbyn, Multilayer graphenes as a platform for interaction-driven physics and topological superconductivity, *Phys. Rev. B* **107**, 104502 (2023).

[57] T. Cea, Superconductivity induced by the intervalley coulomb scattering in a few layers of graphene, *Phys. Rev. B* **107**, L041111 (2023).

[58] T. J. Peltonen, R. Ojajärvi, and T. T. Heikkilä, Mean-field theory for superconductivity in twisted bilayer graphene, *Phys. Rev. B* **98**, 220504(R) (2018).

[59] Y.-Z. Chou, F. Wu, J. D. Sau, and S. Das Sarma, Acoustic-phonon-mediated superconductivity in Bernal bilayer graphene, *Phys. Rev. B* **105**, L100503 (2022).

[60] Y.-Z. Chou, F. Wu, J. D. Sau, and S. Das Sarma, Acoustic-phonon-mediated superconductivity in moiréless graphene multilayers, *Phys. Rev. B* **106**, 024507 (2022).

[61] J. Vučičević, M. O. Goerbig, and M. V. Milovanović, d-wave superconductivity on the honeycomb bilayer, *Phys. Rev. B* **86**, 214505 (2012).

[62] M. V. Hosseini and M. Zareyan, Model of an exotic chiral superconducting phase in a graphene bilayer, *Phys. Rev. Lett.* **108**, 147001 (2012).

[63] M. Vahid Hosseini and M. Zareyan, Graphene-based ferromagnetic superconductors, *Appl. Phys. Lett.* **101**, 252602 (2012).

[64] M. V. Hosseini and M. Zareyan, Unconventional superconducting states of interlayer pairing in bilayer and trilayer graphene, *Phys. Rev. B* **86**, 214503 (2012).

[65] M. V. Hosseini, Inhomogeneous superconductivity in the presence of time-reversal symmetry, *Europhys. Lett.* **110**, 47010 (2015).

[66] E. Thingstad, A. Kamra, J. W. Wells, and A. Sudbø, Phonon-mediated superconductivity in doped monolayer materials, *Phys. Rev. B* **101**, 214513 (2020).

[67] M. Alidoust, A.-P. Jauho, and J. Akola, Josephson effect in graphene bilayers with adjustable relative displacement, *Phys. Rev. Res.* **2**, 032074(R) (2020).

[68] H. Dai, J. Hou, X. Zhang, Y. Liang, and T. Ma, Mott insulating state and d + *id* superconductivity in an abc graphene trilayer, *Phys. Rev. B* **104**, 035104 (2021).

[69] B. Roy and I. F. Herbut, Unconventional superconductivity on honeycomb lattice: Theory of kekule order parameter, *Phys. Rev. B* **82**, 035429 (2010).

[70] A. Kerelsky, L. J. McGilly, D. M. Kennes, L. Xian, M. Yankowitz, S. Chen, K. Watanabe, T. Taniguchi, J. Hone, C. Dean et al., Maximized electron interactions at the magic angle in twisted bilayer graphene, *Nature (London)* **572**, 95 (2019).
[71] Y. Saito, J. Ge, K. Watanabe, T. Taniguchi, and A. F. Young, In-dependent superconductors and correlated insulators in twisted bilayer graphene, Nat. Phys. **16**, 926 (2020).

[72] Y. Cao, D. Chowdhury, D. Rodan-Legrain, O. Rubies-Bigorda, K. Watanabe, T. Taniguchi, T. Senthil, and P. Jarillo-Herrero, Strange metal in magic-angle graphene with near Planckian dissipation, Phys. Rev. Lett. **124**, 076801 (2020).

[73] R. Lyu, Z. Tuchfeld, N. Verma, H. Tian, K. Watanabe, T. Taniguchi, C. N. Lau, M. Randeria, and M. Bockrath, Strange metal behavior of the hall angle in twisted bilayer graphene, Phys. Rev. B **103**, 245424 (2021).

[74] H. Wang, J.-H. Gao, and F.-C. Zhang, Flat band electrons and interactions in rhombohedral trilayer graphene, Phys. Rev. B **87**, 155116 (2013).

[75] T. T. Heikkilä and G. E. Volovik, Flat bands as a route to high-temperature superconductivity in graphite, in *Basic Physics of Functionalized Graphite* (Springer, New York, 2016), pp. 123–143.

[76] M. Sigrist and K. Ueda, Phenomenological theory of unconventional superconductivity, Rev. Mod. Phys. **63**, 239 (1991).

[77] A. Grüneis, C. Attaccalite, L. Wirtz, H. Shiozawa, R. Saito, T. Pichler, and A. Rubio, Tight-binding description of the quasi-particle dispersion of graphite and few-layer graphene, Phys. Rev. B **78**, 205425 (2008).

[78] E. McCann and M. Koshino, The electronic properties of bilayer graphene, Rep. Prog. Phys. **76**, 056503 (2013).

[79] A. Crépieux, E. Pangburn, L. Haurie, O. A. Awoga, A. M. Black-Schaffer, N. Sedlmayr, C. Pépin, and C. Bena, this issue, Superconductivity in monolayer and few-layer graphene. II. Topological edge states and Chern numbers, Phys. Rev. B **108**, 134515 (2023).

[80] L. Yu, Bound state in superconductors with paramagnetic impurities, Chin. J. Phys. (Peking), Engl. Transl. **21**, 75 (1965).

[81] H. Shiba, Classical spins in superconductors, Prog. Theor. Phys. **40**, 435 (1968).

[82] A. I. Rusinov, On the Theory of Gapless Superconductivity in Alloys Containing Paramagnetic Impurities, Sov. Phys. JETP **29**, 1101 (1969).

[83] Y.-Q. Li and T. Zhou, Impurity effect as a probe for the pairing symmetry of graphene-based superconductors, Front. Phys. **16**, 1 (2021).

[84] E. Pangburn, L. Haurie, A. Crépieux, O. A. Awoga, N. Sedlmayr, A. M. Black-Schaffer, C. Pépin, and C. Bena, this issue, Superconductivity in monolayer and few-layer

graphene. III. Impurity-induced subgap states and quasiparticle interference patterns, Phys. Rev. B **108**, 134516 (2023).

[85] O. A. Awoga, T. Löthman, and A. M. Black-Schaffer, Superconductivity and magnetism in the surface states of ABC-stacked multilayer graphene, Phys. Rev. B **108**, 144504 (2023).

[86] K. Saito, *Chemical Physics of Molecular Condensed Matter* (Springer, New York, 2020).

[87] T. Löthman and A. M. Black-Schaffer, Universal phase diagrams with superconducting domes for electronic flat bands, Phys. Rev. B **96**, 064505 (2017).

[88] A. H. Castro Neto and F. Guinea, Impurity-induced spin-orbit coupling in graphene, Phys. Rev. Lett. **103**, 026804 (2009).

[89] A. M. Black-Schaffer, Edge properties and Majorana fermions in the proposed chiral *d*-wave superconducting state of doped graphene, Phys. Rev. Lett. **109**, 197001 (2012).

[90] O. A. Awoga, A. Bouhon, and A. M. Black-Schaffer, Domain walls in a chiral *d*-wave superconductor on the honeycomb lattice, Phys. Rev. B **96**, 014521 (2017).

[91] V. Kaladzhyan, C. Bena, and P. Simon, Characterizing unconventional superconductors from the spin structure of impurity-induced bound states, Phys. Rev. B **93**, 214514 (2016).

[92] K. Yan, H. Peng, Y. Zhou, H. Li, and Z. Liu, Formation of bilayer Bernal graphene: layer-by-layer epitaxy via chemical vapor deposition, Nano Lett. **11**, 1106 (2011).

[93] R. H. Vervuurt, W. M. Kessels, and A. A. Bol, Atomic layer deposition for graphene device integration, Adv. Mater. Interfaces **4**, 1700232 (2017).

[94] M. Yankowitz, J. I.-J. Wang, A. G. Birdwell, Y.-A. Chen, K. Watanabe, T. Taniguchi, P. Jacquod, P. San-Jose, P. Jarillo-Herrero, and B. J. LeRoy, Electric field control of soliton motion and stacking in trilayer graphene, Nat. Mater. **13**, 786 (2014).

[95] Y. Shan, Y. Li, D. Huang, Q. Tong, W. Yao, W.-T. Liu, and S. Wu, Stacking symmetry governed second harmonic generation in graphene trilayers, Sci. Adv. **4**, (2018).

[96] A. B. Kuzmenko, I. Crassee, D. van der Marel, P. Blake, and K. S. Novoselov, Determination of the gate-tunable band gap and tight-binding parameters in bilayer graphene using infrared spectroscopy, Phys. Rev. B **80**, 165406 (2009).

[97] T. O. Wehling, E. Şaşıoğlu, C. Friedrich, A. I. Lichtenstein, M. I. Katsnelson, and S. Blügel, Strength of effective coulomb interactions in graphene and graphite, Phys. Rev. Lett. **106**, 236805 (2011).