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

HAL
open science

# Superconductivity in monolayer and few-layer graphene: III Impurity-induced subgap states and quasi-particle interference patterns

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

► To cite this version:

Emile Pangburn, Louis Haurie, Adeline Crépieux, Oladunjoye A. Awoga, Nicholas Sedlmayr, et al. Superconductivity in monolayer and few-layer graphene: III Impurity-induced subgap states and quasi-particle interference patterns. Physical Review B, 2023, 108, pp.134516. 10.1103/Phys-RevB.108.134516. cea-04258058

HAL Id: cea-04258058

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

Submitted on 25 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.
# Superconductivity in monolayer and few-layer graphene. III. Impurity-induced subgap states and quasiparticle interference patterns

Emile Pangburn,¹ Louis Haurie,¹ Adeline Crépieux,² Oladunjoye A. Awoga,³ Nicholas Sedlmayr,⁴

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 Université, Université de Toulon, CNRS, CPT, Marseille, France

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

⁴Institute of Physics, Maria Curie-Skłodowska University, Plac Marii Skłodowskiej-Curie 1, 20031 Lublin, Poland

⁵Department of Physics and Astronomy, Uppsala University, Box 516, 751 20 Uppsala, Sweden

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

(Received 29 March 2023; revised 29 June 2023; accepted 27 September 2023; published 25 October 2023)

We consider the most energetically favorable symmetry-allowed spin-singlet and spin-triplet superconducting pairing symmetries in monolayer and few-layer graphene and for each calculate the energy spectrum in the presence of a scalar or magnetic impurity. We find that two doubly degenerate subgap states exist for scalar impurities for all types of pairing, except for the spin-singlet s-wave state. For magnetic impurities, two or four subgap states may form depending on the order parameter symmetry. We find that the spin polarization of these states allows one to distinguish between spin-singlet and -triplet pairing; for example, only the spin-triplet states show opposite-energy subgap states with the same spin. We also calculate the quasiparticle interference patterns associated with the subgap states and find that they exhibit features that could distinguish between different types of pairing symmetries, especially a breaking of rotational symmetry for nodal states, which is stronger for the spin-singlet $d_{xy}$ and $d_{x^2-y^2}$ states than for the spin-triplet $p_x$ and $p_y$ states.

DOI: 10.1103/PhysRevB.108.134516

## I. INTRODUCTION

Ever since the discovery of superconductivity in graphene-based systems, such as twisted bilayer graphene [1–12], Bernal (or AB) bilayer graphene [13], and rhombohedral trilayer graphene [14,15], the identification of the pairing symmetries in these unconventional superconductors has been one of the main goals of theoretical and experimental studies. However, the present state of the analysis of microscopic theories for different graphene systems does not allow for a definite answer to this question [16]. Multiple different mechanisms have been proposed based on both phonon-mediated pairing [17,18] and electron-electron interactions [19–21], and they predict different pairing symmetries [8,9,22–37]. Despite the absence of experimental evidence for superconductivity in monolayer graphene, the investigation of pairing symmetry in this system remains an intriguing and relevant topic. Understanding the pairing symmetry in monolayer graphene holds the potential to provide valuable insights into the characteristics of superconductivity in few-layer graphene systems [38].

In two related works [38,39] we examined all expected spin-singlet and spin-triplet superconducting states with the lowest angular momentum ($l \le 3$) in monolayer graphene, as well as in AB-stacked bilayer and ABA- and ABC-stacked trilayer graphene, all without twist. Our goal was to analyze both the basic electronic properties and the topological properties of various superconducting graphene systems by computing their band structure and density of states (DOS), as well as the Chern number and the associated topologically protected edge

states [26,27,40]. The analysis of these properties may help to experimentally distinguish between various order parameters. For example, the DOS, measurable with scanning tunneling microscopy (STM), can, in principle, distinguish between nodal superconductors (SCs) ($d_{xy}$, $d_{x^2-y^2}$, $p_x$, $p_y$ waves), which have a V-shaped DOS, and gapped SCs ($s_{on}$, $s_{ext}$, $p + ip'$, $d + id$, $f$ waves), which have a U-shaped DOS. However, in real experiments it can still be hard to distinguish between these two types of DOSs if the resolution is not sufficient, and disorder may additionally locally perturb the superconducting state [41]. Note that, for simplicity, we choose to assume that all the sources of superconducting pairing for all possible symmetries considered are intralayer, and we do not take into account any interlayer pairing terms, even if to show there is no experimentally justified consensus regarding their effects. The interlayer pairing terms are believed to play a role, for example, in rhombohedral graphite [14], in which superconductivity is observed in the presence of a displacement field. In this final work in the series, we propose another tool to distinguish between different superconducting order parameters in moiréless graphene systems by studying the effects of a single impurity on the local density of states (LDOS) and on the spin-polarized local density of states (SPLDOS). It is already well known from the Anderson theorem [41] that conventional s-wave SCs are not affected by scalar impurities and do not allow the formation of subgap states, whereas scalar impurities usually induce subgap bound (resonance) states in unconventional fully gapped (nodal) SCs [42]. Also, it is well known that a magnetic impurity induces Yu-Shiba-Rusinov subgap states in a SC due to local time reversal symmetry
breaking [43–45]. In what follows we refer to all low-energy impurity states, irrespective of their origin, as simply subgap states.

Although subgap states are often expected, their multiplic- ity and characteristics strongly depend on the underlying sym- metries of the normal state and, especially, the superconduct- ing order parameter [42,46]. Thus, one may hope that studying the features induced by a single impurity would help differen- tiate between the different superconducting states. Motivated by these prospects, we perform an extensive study of the effects of both scalar and magnetic impurities in monolayer, AB bilayer, and trilayer ABA and ABC graphene for the same set of symmetries considered in Refs. [38,39]. Using the well- known T-matrix approach [42], we compute both the spatially averaged LDOS as a function of energy and the Fourier trans- form of the LDOS change induced by the impurity at a given energy (FT-LDOS, also known as the quasiparticle interfer- ence pattern, or QPI). The QPI and the spin-polarized QPI, measurable via STM and spin-polarized STM, respectively, provide a direct connection to angle-resolved photoemission spectroscopy experiments [47–50] and thus contain informa- tion about the band structure of the system. In fact, such measurements have already been used to study impurity scat- tering effects in graphene systems [46,51–57].

We first focus on the energy dependence of the spatially averaged LDOS and SPLDOS. This allows us to study the formation of subgap states, thereby distinguishing between conventional and unconventional SCs. For example, we con- firm that in the presence of scalar impurities the spin-singlet s-wave states do not give rise to subgap states, according to the Anderson theorem [41,46]. On the other hand, magnetic impurities are pair breaking for both spin-singlet and spin- triplet states and thus generate subgap states for all types of pairing symmetries. Here we identify four different subgap states for the fully gapped ( d + id $^{′}$)-, ( p + ip $^{′}$)-, and f -wave symmetries, while the rest exhibit only two subgap states. For the SPLDOS generated by a magnetic impurity we note that pairs of subgap states of the same spin but opposite energies exist for some of the spin-triplet order parameters but none of the spin-singlet order parameters, which thus becomes a clear experimental signature allowing us to identify the existence of a spin-triplet superconducting state. Another difference be- tween the spin-singlet and spin-triplet superconducting states is that the SPLDOS features depend on the impurity spin orientation, while they are automatically independent of the direction of the impurity spin for all spin-singlet supercon- ducting states. We subsequently study the QPI maps and show that the QPI patterns for nodal superconducting states break the sixfold symmetry of the normal state, while the QPI patterns for the gapped superconducting states preserve this symmetry. This establishes that the QPI can distinguish between gapless and nodal order parameters, such as p$_{x}$, p$_{y}$, d$_{xy}$, and d$_{x}$ $^{2}$− y2 waves, from fully gapped order parameters, such as s, s$_{ext}$, p + ip $^{′}$, d + id $^{′}$, and f waves, in supercon- ducting graphene systems. Finally, we note that most of the features of the subgap states are quite generic and unchanged when studying bilayer or trilayer graphene, except for extra subgap states arising in the ABC trilayer case and for a smaller splitting of the features due to the presence of the interlayer coupling.

The rest of this paper is organized as follows. In Sec. II we provide the details of the tight-binding model used and the T-matrix formalism. In Sec. III we focus on monolayer graphene, and we compute the averaged LDOS and SPLDOS for both scalar and magnetic impurities, as well as their mo- mentum dependence, or, equivalently, the QPI. We extend the study to multilayer graphene in Sec. IV, before summarizing our results in Sec. V.

II. MODEL AND METHOD

A. Bulk Hamiltonian

Without trying to justify the pairing mechanism for superconductivity, we consider superconducting graphene de- scribed by a tight-binding Hamiltonian with a pairing term that can take all relevant spin-singlet and spin-triplet symme- tries with the lowest angular momentum. The noninteracting Hamiltonian is given by

$$
\begin{align}
H_0(\mathbf{k}) &= \sum_{\mathbf{k},\alpha} \mu (a_{\mathbf{k}\alpha}^{\dagger} a_{\mathbf{k}\alpha} + b_{\mathbf{k}\alpha}^{\dagger} b_{\mathbf{k}\alpha}) \nonumber \\
&\quad + h_0(\mathbf{k}) a_{\mathbf{k}\alpha}^{\dagger} b_{\mathbf{k}\alpha} + h_0^{*}(\mathbf{k}) b_{\mathbf{k}\alpha}^{\dagger} a_{\mathbf{k}\alpha}, \nonumber \\
h_0(\mathbf{k}) &= -t e^{-ik_y} \left[ 1 + 2e^{3ik_y/2} \cos \left( \frac{\sqrt{3}}{2} k_x \right) \right], \tag{1}
\end{align}
$$

where µ and h$_{0}$ ( k ) are the chemical potential and the kinetic energy, respectively, with t denoting the hopping strength between nearest-neighbor (NN) carbon atoms. Here a † kα ( b † kα ) is the creation operator for an electron with momentum k and spin α in sublattice A ( B ).

We focus primarily on the intralayer NN superconducting pairing, but our results are quite generic and are affected very little if we instead consider intralayer next-nearest-neighbor (NNN) order parameters, like in our earlier works [38,39,58]. This is important to note since self-consistent calculations have shown that the NNN range may be preferred over the NN one for some multilayer graphene configurations [18,58]. The only type of pairing symmetry that cannot be captured by NN pairing is the f -wave state, for which we thus revert to NNN pairing. The SC Hamiltonian can be written as

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

for the spin-singlet channel and

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

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

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

for spin-triplet channels [38,39,58]. Here h η NN ( k ) are the over- all form factors whose expressions depend on both the spin channel chosen and the angular momentum symmetry of the order parameter. Their expressions for the different pairing symmetries are given in Table I.
<table><caption>TABLE I. Expressions for the form factors for different spin-singlet and spin-triplet order parameter symmetries. The overall amplitude for the order parameter is set to $\Delta_0$, the distance between two NN carbon atoms is set to 1, and $\tilde{h}_0(\mathbf{k}) = h_0(\mathbf{k})/t$.</caption><thead><tr><th>η</th><th>Symmetry</th><th>Form factor h<sub>NN</sub><sup>η</sup>(<b>k</b>)</th></tr></thead><tbody><tr><td>0</td><td>s<sub>ext</sub></td><td>h<sub>NN</sub><sup>0,sext</sup>(<b>k</b>) = Δ<sub>0</sub>/√3 h<sub>0</sub>(<b>k</b>)</td></tr><tr><td>0</td><td>d<sub>x<sup>2</sup>-y<sup>2</sup></sub></td><td>h<sub>NN</sub><sup>0,d<sub>x<sup>2</sup>-y<sup>2</sup></sub></sup>(<b>k</b>) = 2Δ<sub>0</sub>/√6 e<sup>-ik<sub>y</sub></sup>[1 - e<sup>-3i/k<sub>y</sub></sup> cos(<sup>√3</sup>&frasl;<sub>2</sub>k<sub>x</sub>)]</td></tr><tr><td>0</td><td>d<sub>xy</sub></td><td>h<sub>NN</sub><sup>0,dxy</sup>(<b>k</b>) = Δ<sub>0</sub>√2i e<sup>i/2 k<sub>y</sub></sup> sin(<sup>√3</sup>&frasl;<sub>2</sub>k<sub>x</sub>)</td></tr><tr><td>x</td><td>p<sub>y</sub></td><td>h<sub>NN</sub><sup>η,p<sub>y</sub></sup>(<b>k</b>) = 2Δ<sub>0</sub>/√6 e<sup>-ik<sub>y</sub></sup>[1 - e<sup>-3i/k<sub>y</sub></sup> cos(<sup>√3</sup>&frasl;<sub>2</sub>k<sub>x</sub>)]</td></tr><tr><td>x</td><td>p<sub>x</sub></td><td>h<sub>NN</sub><sup>η,p<sub>x</sub></sup>(<b>k</b>) = i√2Δ<sub>0</sub>e<sup>i/2 k<sub>y</sub></sup> sin(<sup>√3</sup>&frasl;<sub>2</sub>k<sub>x</sub>)</td></tr><tr><td>x</td><td>f<sub>x</sub>(x<sup>2</sup>-y<sup>2</sup>)</td><td>h<sub>NN</sub><sup>η,f<sub>x</sub></sup>(<b>k</b>) = 2iΔ<sub>0</sub>/√6 [sin(√3k<sub>x</sub>) - 2 sin(<sup>√3</sup>&frasl;<sub>2</sub>k<sub>x</sub>) cos(<sup>3</sup>&frasl;<sub>2</sub>k<sub>y</sub>)]</td></tr></tbody></table>

The explicit matrices corresponding to Eqs. (2)–(5) can be found in Appendix A of Ref. [38].

For the NNN range the above formulas need to be modified such that the pairing terms couple two electrons within the same sublattice. The NNN form factor for the $f_x = f_x(x^2-y^2)$ order parameter, which is the only one considered here, is also given in Table I (we exclude the $f_y(y^2-3x^2)$-wave state because it has multiple nodes and is, as such, highly unfavorable).

For superconducting multilayer graphene, the Hamiltonian is given by

$$H_{\mathbf{k}} = \sum_{\ell=1}^{L} (H_{0}^{(\ell)} + H_{NN}^{(\ell)}) + H_{\text{interlayer}}, \quad (6)$$

where L is the number of layers; $H_0^{(\ell)}$ and $H_{NN}^{(\ell)}$ are the non-interacting and superconducting Hamiltonians, respectively, associated with each layer $\ell$ and given by Eqs. (1)–(5); and $H_{\text{interlayer}}$ is the coupling Hamiltonian between adjacent layers given in Ref. [38]. The interlayer Hamiltonian depends on

three additional parameters, the phase difference $\phi$ between the superconducting state in two adjacent layers and the interlayer couplings $\gamma_1$ and $\gamma_3$ [59], where $\gamma_1$ is the simple interlayer coupling corresponding to hopping between atoms on top of each other and the smaller $\gamma_3$ corresponds to hopping between an atom A in one layer and the neighboring B atoms in the adjacent layer. We have checked that the addition of this trigonal warping $\gamma_3$ with a realistic value $\gamma_3 \le \gamma_1$ does not change our results. We thus set $\gamma_1 = 0.2t$ and $\gamma_3 = 0$ in the rest of the work for simplicity. Note that, as mentioned in the Introduction, we choose to take into account in our model only the intralayer SC pairing terms, and we do not consider the effect of interlayer SC couplings. As it could be argued that, in general, the interlayer couplings are mediated by van der Waals forces and thus are smaller than the intralayer ones [22,60], the interlayer superconducting coupling terms have not been demonstrated experimentally to be negligible, and they may play a role, for example, in rhombohedral graphite [14]. However in this work we assume for simplicity that all SC pairings arise from intralayer couplings.

The band structures of both the normal and superconducting states for various pairing symmetries are illustrated in Ref. [38].

Collating the operators in each layer $\ell$ into a vector, the Hamiltonian can be expressed as

$$H_{\mathbf{k}} = \frac{1}{2} \Psi_{\mathbf{k}}^{\dagger} \hat{H}_{\text{BdG}} \Psi_{\mathbf{k}}, \quad (7)$$

using the basis

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

where $\Psi_{\mathbf{k}}$ thus combines all individual-layer bases $\Psi_{\mathbf{k}\ell}$ and $\hat{H}_{\text{BdG}}$ is the $8L \times 8L$ Bogoliubov–de Gennes (BdG) Hamiltonian matrix. The factor of 8 corresponds to a product of two spins, two sublattices, and the particle-hole doubling of the degrees of freedom in the BdG formalism. Finally, the retarded Green’s function for this system is given by

$$\hat{G}^r(E, \mathbf{k}) = [E + i\delta - \hat{H}_{\text{BdG}}(\mathbf{k})]^{-1}, \quad (9)$$

<figure><img src="image_3.png" /><figcaption>FIG. 1. $\delta\rho(E)$ as a function of energy $E$ and impurity strength $U$ for a scalar impurity. We take $\mu = 0.4t$ and $\Delta_0 = 0.4t$. The dotted lines indicate the SC gap edge, which, as noted in Ref. [38], does not always lie at an energy equal to $\Delta_0$ and may depend on various parameters, including the symmetry of the SC order parameter.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>FIG. 2. $\delta\rho(E)$ as a function of energy $E$ and impurity strength $J_z$ for a z-magnetic impurity. We take $\mu = 0.4t$ and $\Delta_0 = 0.4t$. The dotted lines indicate the gap edge.</figcaption></figure>

with $\delta$ being the quasiparticle lifetime. We set $\delta = 0.03$ in the rest of the work.

## B. Impurity scattering

In this work we are interested in the consequences of introducing a pointlike (scalar or magnetic) impurity. Using the basis in Eq. (8), the Hamiltonian matrix for such pointlike impurity can be written as

$$ \hat{V} = \tau^z \otimes \hat{V}, \quad \hat{V} = \hat{u} \otimes \hat{v}. \qquad (10) $$

Here

$$ \hat{v} = U\sigma^0 + J\sigma^v, \qquad (11) $$

where $\tau^v(\sigma^v)$ are the v-Pauli matrices in the particle-hole (spin) space and $\sigma^0$ is the 2×2 identity matrix, while $\hat{u}$ is a 2L×2L matrix for which all the elements are equal to zero, except for one diagonal element, whose matrix position $i_{imp}$ corresponds to the layer/sublattice of the impurity, which we take to be equal to 1. The parameters $U$ and $J$ are, respectively, the strengths of the scalar and magnetic impurities.

To compute the corresponding variation of the unpolarized (LDOS) and spin-polarized (SPLDOS) local densities of states, we use the T-matrix approach [42]. The T-matrix can be written as

$$ T(E) = \left[ 1_{8L} - \hat{V} \int \frac{d^2\mathbf{k}}{(2\pi)^2} \hat{G}^r(E, \mathbf{k}) \right]^{-1} \hat{V}, \qquad (12) $$

where $1_{8L}$ is an 8L×8L identity matrix. The physical observables (here LDOS and SPLDOS), which can be measured near an impurity, can be expressed directly in terms of this T-matrix if we assume the dilute-limit approximation, such that the impurities are well separated from each other. The Fourier transform of the change in the LDOS induced by the impurity $\delta\rho(\mathbf{q}, E)$ and the same quantity for the SPLDOS $\delta S_v(\mathbf{q}, E)$ can then be written as

$$ \delta\rho(\mathbf{q}, E) = -\frac{1}{2\pi i} \int \frac{d^2\mathbf{k}}{(2\pi)^2} \sum_b \left[ \tilde{g}_{b, \uparrow\uparrow}(E, \mathbf{q}, \mathbf{k}) + \tilde{g}_{b, \downarrow\downarrow}(E, \mathbf{q}, \mathbf{k}) \right], \qquad (13) $$

<figure><img src="image_5.png" /><figcaption>FIG. 3. $\delta S_v(E)$ as a function of energy $E$ and impurity strength $J_z$ for a z-magnetic impurity. We take $\mu = 0.4t$ and $\Delta_0 = 0.4t$. The dotted lines indicate the gap edge.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>FIG. 4. $\delta S_x(E)$ as a function of energy $E$ and impurity strength $J_x$ for an $x$-magnetic impurity. We take $\mu = 0.4t$ and $\Delta_0 = 0.4t$. The dotted lines indicate the gap edge.</figcaption></figure>

$$ \delta S_x(\mathbf{q}, E) = -\frac{1}{2\pi i} \int \frac{d^2\mathbf{k}}{(2\pi)^2} \sum_b [\tilde{g}_{b,\uparrow\downarrow}(E, \mathbf{q}, \mathbf{k}) + \tilde{g}_{b,\downarrow\uparrow}(E, \mathbf{q}, \mathbf{k})], \quad (14) $$

$$ \delta S_y(\mathbf{q}, E) = -\frac{1}{2\pi} \int \frac{d^2\mathbf{k}}{(2\pi)^2} \sum_b [g_{b,\uparrow\downarrow}(E, \mathbf{q}, \mathbf{k}) - g_{b,\downarrow\uparrow}(E, \mathbf{q}, \mathbf{k})], \quad (15) $$

$$ \delta S_z(\mathbf{q}, E) = -\frac{1}{2\pi i} \int \frac{d^2\mathbf{k}}{(2\pi)^2} \sum_b [\tilde{g}_{b,\uparrow\uparrow}(E, \mathbf{q}, \mathbf{k}) - \tilde{g}_{b,\downarrow\downarrow}(E, \mathbf{q}, \mathbf{k})], \quad (16) $$

where the index $b$ runs over all electron bands (the hole bands are not taken into account since, experimentally, only the available electron density of states is measured) and

$$ g/\tilde{g}(E, \mathbf{q}, \mathbf{k}) = G^r(E, \mathbf{q})T(E)G^r(E, \mathbf{q}+\mathbf{k}) \pm [G^r(E, \mathbf{k}+\mathbf{q})]^*T^*(E)[G^r(E, \mathbf{q})]^*. \quad (17) $$

At $\mathbf{q} = 0$, the quantities $\delta\rho(\mathbf{q} = 0, E) \rightarrow \delta\rho(E)$ and $\delta S_v(\mathbf{q} = 0, E) \rightarrow \delta S_v(E)$ correspond to the spatially averaged disorder-induced LDOS and SPLDOS, respectively. In the next two sections, we plot and analyze $\delta\rho(E)$ and $\delta S_v(E)$ as a function of energy and impurity strength to establish the formation of subgap states. Furthermore, at constant energy, the QPI patterns described by Eqs. (13)–(16) provide a map in reciprocal space of the possible scattering processes. Experimentally, the QPI patterns are obtained by performing a fast Fourier transform of the STM measurements of the LDOS in real space [46,54].

### III. MONOLAYER GRAPHENE

#### A. Unpolarized and spin-polarized local densities of states

We first consider the spatially averaged LDOS, i.e., $\delta\rho(E)$ and $\delta S_v(E)$, in the presence of both scalar and magnetic

impurities. If and when subgap states form, these quantities will display clear peaks inside the superconducting gap. The position of the peaks may depend on various parameters, such as the impurity strength, the amplitude of the SC order parameter and its symmetry, and the chemical potential. The gap edge is identified by the energy of the superconducting coherence peak. We find it numerically in the impurity-free DOS $\rho_0(E)$ in the same manner for both fully gapped and nodal superconducting states.

In Fig. 1 we plot $\delta\rho(E)$ for all the order parameter symmetries as a function of energy and impurity strength in the presence of a scalar impurity. We first note that for the on-site (ON) and $s_{ext}$ symmetries, there is no impurity subgap peak. This is consistent with previous observations in the literature: conventional s-wave SCs are unaffected by the presence of nonmagnetic or scalar impurities [41,42,46] and thus do not show any subgap states since these impurities do not break time reversal symmetry [61]. This reasoning can also be applied to extended s-wave SCs: as long as the chemical potential is chosen such that the order parameter is almost constant along the Fermi surface, the phenomenology is approximately the same [62].

In contrast, the spin-singlet nodal, $d_{xy}$- and $d_{x^2-y^2}$-wave, and the fully gapped, chiral ($d_{xy} + id_{x^2-y^2}$)-wave ($d + id'$) states show spin degenerate subgap states in the presence of a scalar impurity. For these states the Anderson theorem does not forbid the presence of subgap states, even for a nonmagnetic impurity. The physical interpretation is that scattering by an impurity disturbs the phase distribution for some particular directions of the momenta in all these nontrivial superconducting states [42]. This has also been noted in d-wave SCs on the square lattice, modeling the cuprate SCs [42], and is also in agreement with former studies on the chiral ($d + id'$)-wave superconducting state in graphene [46,63]. Similarly, we find subgap states for all spin-triplet states, both the nodal $p_x$- and $p_y$-wave states, and the fully gapped, chiral ($p_x + ip_y$)-wave ($p + ip'$) superconducting state, as well as the fully gapped $f_x$-wave state. For all these subgap states we find that their energies evolve with the impurity strength such that the states cross zero energy at given, but different, impurity strengths. A similar observation was also noted in Ref. [64].

We next study the effect of a magnetic impurity on the formation of subgap states. First, we plot in Fig. 2 $\delta\rho(E)$ as a function of energy and impurity strength for all the order parameter symmetries considered. We first note that for a magnetic impurity we find subgap states for all types of pairing symmetries, including the s-wave states. We also note that for the fully gapped ($d + id'$)-, ($p + ip'$)-, and f-wave states, the spin degeneracy has been lifted, and we now have four distinct subgap states rather than two pairs of degenerate ones. On the other hand, the nodal superconducting $d_{xy}$-, $d_{x^2-y^2}$-, $p_x$-, and $p_y$-wave states show the same number of subgap states, i.e., two distinct states, for both magnetic and scalar impurities; however, a pair of extra impurity states often arises outside the gap. As a consequence, the number of subgap states can be used as a simple tool to experimentally distinguish between various order parameters; that is, if we can identify four distinct subgap states, then we can be sure to have a ($d + id'$)-, ($p + ip'$)-, or f-wave pairing. We note, however, that the reverse may not always work
<figure><img src="image_7.png" /><figcaption>FIG. 5. $|\delta\rho(\mathbf{q})|$ at the values of impurity strength $U$ and energy provided in Table II for a scalar impurity. We take $\mu = 0.4t$ and $\Delta_0 = 0.4t$. The Brillouin zone is indicated by dashed lines.</figcaption></figure>

since the four states may be too close together to distinguish experimentally.

In order to get a better understanding of what happens for magnetic impurities, we also look at the spin polarization of the induced subgap states. Here we find that each impurity gives rise to a nonzero spin polarization only in the spin channel parallel to its spin direction; thus, for an α-magnetic impurity we plot only the α-magnetic component of the SPL-DOS, with $\alpha = x, y, z$. We first plot in Fig. 3 the SPLDOS $\delta S_z(E)$ as a function of magnetic impurity strength and energy for all types of pairings for a z-magnetic impurity. Here we

find that the SPLDOS shows even more clearly the difference between the two and four subgap scenarios, as it clearly differentiates between the different states. We can next ask if the interplay between the direction of the impurity spin and the choice of the triplet channel, which we fixed to $x$ above, influences the results. We have checked that for this particular choice a y-magnetic impurity yields exactly the same behavior as the z-magnetic impurity. However, for an x-magnetic impurity, the LDOS is unchanged, but the x-SPLDOS differs in the spin-triplet channel. In Fig. 4 we illustrate this by plotting the x-SPLDOS for all spin-triplet states.

<figure><img src="image_8.png" /><figcaption>FIG. 6. $|\delta\rho(\mathbf{q})|$ at zero energy and the corresponding impurity strength values $U$ in Table III for a scalar impurity. We take $\mu = 0.4t$ and $\Delta_0 = 0.4t$. The Brillouin zone is indicated by dashed lines.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>FIG. 7. $|\delta\rho(\mathbf{q})|$ at zero energy and the corresponding impurity strength $J_z = J$ values in Table III. We take $\mu = 0.4t$ and $\Delta_0 = 0.4t$. The Brillouin zone is indicated by dashed lines.</figcaption></figure>

We further note in Fig. 3 that for the spin-singlet superconducting states, opposite-energy states have opposite spin, while for the spin-triplet states the opposite-energy states have the same spin. This, however, seems to be a feature dependent on the direction of the impurity spin: for the different impurity directions considered in Fig. 4, we find that the states with opposite energy also have opposite spin for spin-triplet pairing. We thus conclude that the peculiar occurrence of having the same spin for subgap states with opposite energies is a distinguishing characteristic of a spin-triplet pairing state ($p_x, p_y, p + ip'$, and $f$ waves) and could be used as an experimental signature to identify a spin-triplet order parameter. Moreover, the dependence of the SPLDOS on the direction of the impurity spin is also a characteristic unique to spin-triplet pairing, which additionally could be used to distinguish between spin-singlet and spin-triplet triplet order parameters.

## B. Quasiparticle interference

Figures 1 and 2 describe the dependence of the average LDOS change induced by an impurity as a function of energy and impurity strength and thus tell us at which energy the subgap states form. In what follows we are interested in the spatial dependence of the subgap states. In particular, we study the Fourier transform of the LDOS and of the SPLDOS at a given subgap energy peak as a function of momentum. We primarily focus on two different peak energies, $E = 0$ and $E \neq 0$. In Tables II and III we provide the values of the scalar

<table><caption>TABLE II. Values of impurity strength and energy used to generate the QPI plot for scalar (U) and magnetic (J) impurities for all considered order parameter symmetries in Fig. 5.</caption><thead><tr><th>Energy</th><th>Scalar</th><th>Magnetic</th><th>Symmetry</th></tr></thead><tbody><tr><td>0.02</td><td></td><td>J = 3.5</td><td>s<sub>ext</sub></td></tr><tr><td>0.1</td><td></td><td>J = 2.5</td><td>d<sub>x<sup>2</sup>-y<sup>2</sup></sub></td></tr><tr><td>0.1</td><td>U = 2.5</td><td>J = 2.5</td><td>d<sub>xy</sub></td></tr><tr><td>0.1</td><td>U = 3</td><td>J = 3</td><td>d<sub>x<sup>2</sup>-y<sup>2</sup></sub> + id<sub>xy</sub></td></tr><tr><td>0.1</td><td>U = 1.5</td><td>J = 1.5</td><td>p<sub>x</sub></td></tr><tr><td>0.1</td><td>U = 1.5</td><td>J = 1.5</td><td>p<sub>y</sub></td></tr><tr><td>0.1</td><td>U = 2</td><td>J = 2</td><td>p<sub>x</sub> + ip<sub>y</sub></td></tr><tr><td>0.2</td><td>U = 1.5</td><td>J = 3</td><td>f<sub>x</sub>(x<sup>2</sup>-y<sup>2</sup>)</td></tr><tr><td>0.2</td><td>U = 3</td><td>J = 1.5</td><td>s<sub>ON</sub></td></tr></tbody></table>

<table><caption>TABLE III. Values of impurity strength corresponding to a zero-energy subgap state used to generate the QPI for scalar (U) and magnetic (J) impurities for all considered order parameter symmetry in Figs. 6–9.</caption><thead><tr><th>Energy</th><th>Scalar</th><th>Magnetic</th><th>Symmetry</th></tr></thead><tbody><tr><td>0</td><td></td><td>J = 2</td><td>s<sub>ext</sub></td></tr><tr><td>0</td><td></td><td>J = 5</td><td>d<sub>x<sup>2</sup>-y<sup>2</sup></sub></td></tr><tr><td>0</td><td>U = 5</td><td>J = 5</td><td>d<sub>xy</sub></td></tr><tr><td>0</td><td>U = 6</td><td>J = 6</td><td>d<sub>x<sup>2</sup>-y<sup>2</sup></sub> + id<sub>xy</sub></td></tr><tr><td>0</td><td>U = 5</td><td>J = 5</td><td>p<sub>x</sub></td></tr><tr><td>0</td><td>U = 5</td><td>J = 5</td><td>p<sub>y</sub></td></tr><tr><td>0</td><td>U = 6</td><td>J = 6</td><td>p<sub>x</sub> + ip<sub>y</sub></td></tr><tr><td>0</td><td>U = 6</td><td>J = 6</td><td>f<sub>x</sub>(x<sup>2</sup>-y<sup>2</sup>)</td></tr><tr><td>0</td><td>U = 2.5</td><td>J = 2.5</td><td>s<sub>ON</sub></td></tr></tbody></table>

impurity strength U and of the magnetic impurity strength J and the corresponding peak energies for each order parameter symmetry. All energies are given in units of $t$; that is, we set $t = 1$.

In the following we also plot only the absolute values of $\delta\rho(\mathbf{q})$ and $\delta S_\alpha(\mathbf{q})$, $\alpha = x, y, z$. This is because, for the hexagonal structure of the lattice, these are generally complex quantities that have both nonzero real and imaginary parts. However, at present it is very hard to distinguish experimentally between the real and imaginary parts: experiments calculating the QPI patterns based on fast Fourier transform cannot precisely keep track of either the phase or the sign. Moreover, spin-polarized STM experiments with corresponding QPIs are still in their infancy. Thus, in order to avoid the overload of information, we focus on only the absolute Fourier transform values. If more accurate experimental data become available, this study could easily be extended and refined to take into account separately the real and imaginary parts of the Fourier transforms of both the LDOS and SPLDOS.

Similar to what we did above, we start by considering a scalar impurity ($U \neq 0$ and $J = 0$) and calculate the QPI patterns for both nonzero energy in Fig. 5 and zero energy in Fig. 6. Since the $s_{ON}$-wave and $s_{ext}$-wave order parameters do not exhibit any subgap states for a scalar impurity, we do not include them in the scalar impurity QPI analysis. Overall, the QPI patterns are dominated by a central feature at the center of the Brillouin zone (Γ point), corresponding to intranodal scattering (in the normal-state band structure) of the electrons by the impurity, and by six features localized at the corners of the Brillouin zone (K points), corresponding to internodal scattering.

We further note that the gapless nodal states clearly produce a QPI pattern that breaks the sixfold symmetry, while the gapped states, i.e., the $f$-wave and the chiral $(d + id')$- and $(p + ip')$-wave states, all show QPI patterns that preserve the full rotation symmetry of the lattice. This is fully consistent with the symmetries of the superconducting order parameter, modulo the order parameter phase that might change sign but which seemingly does not affect the QPI, in contrast to the case of $d$-wave cuprates [65]. It is also fully consistent with the fact that these states have a symmetry-preserving superconducting band structure [38] for all gapped states. QPI would thus be a good experimental tool to distinguish between
<figure><img src="image_10.png" /><figcaption>FIG. 8. $|\delta S_z(\mathbf{q})|$ at zero energy and corresponding impurity strength values $J_z = J$ in Table III. We take $\mu = 0.4t$ and $\Delta_0 = 0.4t$. The Brillouin zone is indicated by dashed lines.</figcaption></figure>

nodal states that break rotation symmetry and the gapped states, which do not. A similar observation was already made in Refs. [46,63] in the comparison of the nodal d-wave states and the chiral d-wave states. Here we establish that this also holds for both spin-singlet d-wave and spin-triplet p-wave symmetries in graphene.

We next analyze the QPI patterns generated by a magnetic impurity ($J \neq 0$ and $U = 0$). We here plot only the QPI corresponding to the zero-energy peaks, as we find that the effect for the nonzero energy subgap states is very similar. We further find that in the presence of a z-magnetic impurity we recover the same features for $|\delta\rho(\mathbf{q})|$ as those for the scalar

impurity depicted in Fig. 6. The main difference is that, in the presence of a magnetic impurity, the states with $s_{ON}$-wave and $s_{ext}$-wave order parameters also exhibit subgap states. Thus, to avoid repetition, in Fig. 7 we plot $|\delta\rho(\mathbf{q})|$ for only these two s-wave states, for which the corresponding values of the z-magnetic impurity strength, $J = J_z$, are given in Table III.

While $|\delta\rho(\mathbf{q})|$ does not show any significant differences, $|\delta S_z(\mathbf{q})|$, i.e., the spin-polarized LDOS, shows more interesting features, which we plot in Fig. 8 for the same z-magnetic impurity. The main differences from the QPI of the scalar impurity are a ringlike feature arising in the center of the Brillouin zone in the d-wave superconducting states and a reduction in the asymmetry for the K-point features.

<figure><img src="image_11.png" /><figcaption>FIG. 9. $|\delta S_x(\mathbf{q})|$ at zero energy and corresponding impurity strength values $J_x = J$ in Table III for an x-magnetic impurity. We take $\mu = 0.4t$ and $\Delta_0 = 0.4t$. The Brillouin zone is indicated by dashed lines.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>FIG. 10. $\delta\rho(E)$ as a function of energy and impurity strength $U$ for a scalar impurity at $\Delta_0 = 0.4$ and $\mu = 0.4$ for the nodal SC states in ABC-stacked trilayer graphene. The dotted lines indicate the gap edge.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>FIG. 11. $\delta\rho(E)$ as a function of energy and impurity strength $J_z$ for a z-magnetic impurity at $\Delta_0 = 0.4$ and $\mu = 0.4$ for the nodal SC states in ABC-stacked trilayer graphene. The dotted lines indicate the gap edge.</figcaption></figure>

We also note that for the spin-triplet pairing states, the x-SPLDOS QPI for an x-magnetic impurity is different from the y and z ones and becomes more reminiscent of the spin-singlet ones, which have no dependence on x-, y-, and z-magnetic impurity directions. For example, as depicted in Fig. 9, the spin-polarized QPIs for $p_{x-}$, $p_{y-}$, and $(p + ip')$-wave states at zero energy and the corresponding impurity strength values $J_x = J$ provided in Table III acquire more similarities to the $d_{xy-}$, $d_{x^2-y^2-}$, and $(d + id')$-wave states, in that they show a ring of high intensity for the feature in the center of the

<figure><img src="image_14.png" /><figcaption>FIG. 12. $\delta S_z(E)$ as a function of energy and impurity strength $J_z$ for a z-magnetic impurity at $\Delta_0 = 0.4$ and $\mu = 0.4$ for the nodal SC states in ABC-stacked trilayer graphene. The dotted lines indicate the gap edge.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>FIG. 13. $\delta S_x(E)$ as a function of energy and impurity strength $J_x$ for an x-magnetic impurity at $\Delta_0 = 0.4$ and $\mu = 0.4$ for the spin-triplet nodal SC states in ABC-stacked trilayer graphene. The dotted lines indicate the gap edge.</figcaption></figure>

Brillouin zone and increased asymmetry for the features at the corners of the Brillouin zone.

## IV. MULTILAYER GRAPHENE

We next consider both AB-stacked bilayer graphene and ABC- or ABA-stacked trilayer graphene. We first note that for most of the order parameter symmetries, the number of subgap states, their impurity-strength dependence, and their spin dependences are quite universal, generic features and do not depend on the number of layers or the stacking. In what follows, to avoid redundancy, we present the LDOS and SPLDOS results only when there is a difference from the generic case. In particular, we find differences for ABC-stacked trilayer graphene in the presence of gapless $d_{xy-}$, $d_{x^2-y^2-}$, $p_{x-}$, or $p_{y}$-wave order parameters.

Figures 10 and 11 show $\delta\rho(E)$ in ABC-stacked trilayer graphene for all nodal order parameters in the presence of a scalar impurity and a z-magnetic impurity, respectively. Since, in our calculations, the LDOS and SPLDOS are averaged over all atoms in all layers, the results do not depend on the position of the impurity chosen, so we here arbitrarily consider an impurity in the top layer, located on an atom that does not sit on top of any atoms in the neighboring layer. We note that for ABC-stacked trilayer graphene extra subgap states appear besides the two subgap states observed in all the other graphene systems. We have checked using tight-binding calculations that even when these states are close to zero energy for extended parameter ranges, they do not appear to correspond to Majorana zero modes. Furthermore, Fig. 12 plots $\delta S_z(E)$ for the z-magnetic impurity. For completeness we show in Fig. 13 the effect of changing the spin orientation for the spin-triplet nodal states by plotting $\delta S_x(E)$ for an x-magnetic impurity. Like before, we find that the x-magnetic impurity shows a different behavior than y- and z-magnetic impurities for spin-triplet order parameters.

We next turn to the QPI patterns. For simplicity we focus first on only the zero-energy plots for AB-bilayer graphene in the presence of a scalar impurity. The corresponding impurity strength values are quasi-identical to those for the monolayer, and thus, we use the same values as those presented in Table III. These results can be generalized to the other configurations because we find the differences from the monolayer analysis to be quite generic. We here choose to calculate
<figure><img src="image_16.png" /><figcaption>FIG. 14. $|\delta\rho(\mathbf{q})|$ at zero energy and the corresponding impurity strength values $U$ in Table III, evaluated in the top layer for a scalar impurity placed in the bottom layer on an A-sublattice atom. We take $\Delta_0 = 0.4$ and $\mu = 0.4$. The Brillouin zone is indicated by dashed lines.</figcaption></figure>

only the contribution to the LDOS from the top-layer atoms since that is what is measured experimentally [56,57]. Like in Refs. [56,57], the QPI patterns are different depending on whether the impurity is placed in the top layer or in the bottom layer. However, the measured QPI for a given sample becomes an average between all possible contributions given a random distribution of impurities between the atoms in the two layers. In Fig. 14 we plot the QPI resulting from a bottom-layer A-sublattice impurity, while in Fig. 15 we plot the QPI from a top-layer A-sublattice impurity. Note that here in the top layer the A atom is the atom that does not sit on

top of another atom, while the A atom in the bottom layer is the one sitting directly underneath another atom. We find that the main difference for a bottom-layer impurity, compared to a top-layer impurity, consists of having a more equal intensity between the central feature at the $\Gamma$ point (corresponding to intranodal scattering) and the features at the corners of the Brillouin zone (corresponding to internodal scattering). Thus, for the A bottom-layer impurity the corner features appear to be sharper. Note also that the features exhibit an extra split due to the interband effect introduced by the interlayer hopping compared to monolayer graphene.

<figure><img src="image_17.png" /><figcaption>FIG. 15. $|\delta\rho(\mathbf{q})|$ at zero energy and the corresponding impurity strength values $U$ in Table III, evaluated in the top layer for a scalar impurity placed in the top layer on an A-sublattice atom. We take $\Delta_0 = 0.4$ and $\mu = 0.4$. The Brillouin zone is indicated by dashed lines.</figcaption></figure>
For ABA- and ABC-stacked trilayer graphene we obtain overall similar QPI patterns (not shown), underlining the generic features that we observe for all QPI patterns: (i) breaking of the sixfold rotational symmetry for the nodal *d*-wave and *p*-wave states, which is more pronounced for the *d*-wave symmetries; (ii) the dependence of the details of the QPI patterns on the nature of the impurity (magnetic or scalar), as well as its spin direction for the spin-triplet states; and (iii) spin-polarized measurements that help us distinguish better between various order parameter symmetries. Overall, these features could help us to experimentally identify the symmetry of the superconducting states in different mono- and few-layer graphene systems.

V. CONCLUSION

We calculated the impurity-induced LDOS and SPLDOS, as well as their Fourier transforms (through the QPI patterns), for SC monolayer, AB-stacked bilayer, and ABA- and ABC-stacked trilayer graphene for all expected SC order parameters resulting from NN pairing (NNN pairing in the case of *f*-wave symmetry). We assumed only SC intralayer pairing terms in our model and neglected the role of interlayer SC couplings. We analyzed the formation of subgap states as a function of energy and impurity strength and found that the number of subgap bound states depends on the type of order parameter. For a scalar impurity we found no subgap states for *s* waves, both on-site and extended *s* waves, while two spin-degenerate subgap states appear for all other order parameter symmetries. For a magnetic impurity we found two subgap states for order parameters with *s*-wave symmetries and for nodal states with *d*<sub>xy</sub>-, *d*<sub>x<sup>2</sup>-y<sup>2</sup></sub>-, *p*<sub>x</sub>-, and *p*<sub>y</sub>-wave symmetries, while four subgap states exist for the fully gapped (*d* + *id*')-, (*p* + *ip*')-, and

*f*-wave states. The spin polarization of the impurity states is also different depending on whether one has a spin-singlet or spin-triplet order parameter and could thus be used to distinguish between the two. In particular, the spin-triplet superconducting states are the only ones for which the opposite-energy subgap states may have an identical spin polarization and for which the spin structure of the subgap states may depend on the direction of the impurity spin. These observations could provide an experimental test to distinguish unambiguously, via spin-polarized STM, between spin-singlet and spin-triplet SC order parameters, as well as between gapped and nodal pairings. The analysis of the QPI patterns additionally showed a breaking of the sixfold symmetry for nodal states, while the gapped states preserve this crystalline symmetry, in agreement with the observation that these states also have a symmetry-preserving superconducting band structure [38]. Except in a few peculiar situations, our results do not change significantly for bilayer or trilayer graphene, such that we can easily extend our conclusions to multilayer graphene, and thus, the features described here are quite generic and independent of the number of layers or the graphene layer stacking.

ACKNOWLEDGMENTS

A.M.B.-S. acknowledges 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. N.S. would like to thank the National Science Centre (NCN, Poland) for funding under Grant No. 2018/29/B/ST3/01892. O.A.A. acknowledges funding from NanoLund.

[1] 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).

[2] 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).

[3] X. Lu, P. Stepanov, W. Yang, M. Xie, M. A. Aamir, I. Das, C. Urgell, K. Watanabe, T. Taniguchi, G. Zhang, A. Bachtold, A. H. MacDonald, D. K. Efetov, Superconductors, orbital magnets and correlated states in magic-angle bilayer graphene, Nature (London) **574**, 653 (2019).

[4] 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).

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

[6] 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).

[7] Y. Cao, D. Rodan-Legrain, J. M. Park, N. F. 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).

[8] 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).

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

[10] 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).

[11] 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).

[12] A. Fischer, Z. A. 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).

[13] H. Zhou, L. Holleis, Y. Saito, L. Cohen, W. Huynh, C. L. Patterson, F. Yang, T. Taniguchi, K. Watanabe, and
A. F. Young, Isospin magnetism and spin-polarized superconductivity in Bernal bilayer graphene, *Science* **375**, 774 (2022).

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

[15] 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).

[16] E. Lake, A. S. Patri, and T. Senthil, Pairing symmetry of twisted bilayer graphene: A phenomenological synthesis, *Phys. Rev. B* **106**, 104506 (2022).

[17] B. Lian, Z. Wang, and B. A. Bernevig, Twisted bilayer graphene: A phonon-driven superconductor, *Phys. Rev. Lett.* **122**, 257002 (2019).

[18] 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).

[19] 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).

[20] A. Jimeno-Pozo, H. Sainz-Cruz, T. Cea, P. A. Pantaleón, and F. Guinea, Superconductivity from electronic interactions and spin-orbit enhancement in bilayer and trilayer graphene, *Phys. Rev. B* **107**, L161106 (2023).

[21] P. A. Pantaleón, A. Jimeno-Pozo, H. Sainz-Cruz *et al.*, Superconductivity and correlated phases in non-twisted bilayer and trilayer graphene, *Nat. Rev. Phys.* **5**, 304 (2023).

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

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

[24] 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).

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

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

[27] 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).

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

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

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

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

[32] B. Roy and V. Juričić, Unconventional superconductivity in nearly flat bands in twisted bilayer graphene, *Phys. Rev. B* **99**, 121407(R) (2019).

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

[34] A. L. Szabó and B. Roy, Extended Hubbard model in undoped and doped monolayer and bilayer graphene: Selection rules and organizing principle among competing orders, *Phys. Rev. B* **103**, 205135 (2021).

[35] A. L. Szabó and B. Roy, Competing orders and cascade of degeneracy lifting in doped Bernal bilayer graphene, *Phys. Rev. B* **105**, L201107 (2022).

[36] S. Wolf, T. Gardener, K. Le Hur, and S. Rachel, Topological superconductivity on the honeycomb lattice: Effect of normal state topology, *Phys. Rev. B* **105**, L100505 (2022).

[37] S. Wolf, T. L. Schmidt, and S. Rachel, Unconventional superconductivity in the extended Hubbard model: Weak-coupling renormalization group, *Phys. Rev. B* **98**, 174515 (2018).

[38] E. Pangburn, L. Haurie, A. Crépieux, O. A. Awoga, A. M. Black-Schaffer, C. Pépin, and C. Bena, this issue, Superconductivity in monolayer and few-layer graphene. I. Review of possible pairing symmetries and basic electronic properties, *Phys. Rev. B* **108**, 134514 (2023).

[39] 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).

[40] M. Z. Hasan and C. L. Kane, Colloquium: Topological insulators, *Rev. Mod. Phys.* **82**, 3045 (2010).

[41] P. W. Anderson, Theory of dirty superconductors, *J. Phys. Chem. Solids* **11**, 26 (1959).

[42] A. V. Balatsky, I. Vekhter, and J.-X. Zhu, Impurity-induced states in conventional and unconventional superconductors, *Rev. Mod. Phys.* **78**, 373 (2006).

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

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

[45] A. I. Rusinov, On the theory of gapless superconductivity in alloys containing paramagnetic impurities, *Sov. Phys. JETP* **29**, 1101 (1969).

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

[47] L. Simon, C. Bena, F. Vonau, M. Cranney, and D. Aubel, Fourier-transform scanning tunnelling spectroscopy: The possibility to obtain constant-energy maps and band dispersion using a local measurement, *J. Phys. D* **44**, 464010 (2011).

[48] J. Hoffman, K. McElroy, D.-H. Lee, K. Lang, H. Eisaki, S. Uchida, and J. Davis, Imaging quasiparticle interference in $Bi_2Sr_2CaCu_2O_{8+\delta}$, *Science* **297**, 1148 (2002).

[49] K. McElroy, R. Simmonds, J. Hoffman, D.-H. Lee, J. Orenstein, H. Eisaki, S. Uchida, and J. Davis, Relating atomic-scale electronic phenomena to wave-like quasiparticle states in superconducting $Bi_2Sr_2CaCu_2O_{8+\delta}$, *Nature* (London) **422**, 592 (2003).

[50] T. Zhang, P. Cheng, X. Chen, J.-F. Jia, X. Ma, K. He, L. Wang, H. Zhang, X. Dai, Z. Fang, X. Xie, and Q.-K.
Xue, Experimental demonstration of topological surface states protected by time-reversal symmetry, Phys. Rev. Lett. **103**, 266803 (2009).

[51] C. Bena and S. A. Kivelson, Quasiparticle scattering and local density of states in graphite, Phys. Rev. B **72**, 125432 (2005).

[52] C. Bena, Effect of a single localized impurity on the local density of states in monolayer and bilayer graphene, Phys. Rev. Lett. **100**, 076601 (2008).

[53] G. M. Rutter, J. Crain, N. Guisinger, T. Li, P. First, and J. Stroscio, Scattering and interference in epitaxial graphene, Science **317**, 219 (2007).

[54] I. Brihuega, P. Mallet, C. Bena, S. Bose, C. Michaelis, L. Vitali, F. Varchon, L. Magaud, K. Kern, and J.-Y. Veuillen, Quasiparticle chirality in epitaxial graphene probed at the nanometer scale, Phys. Rev. Lett. **101**, 206802 (2008).

[55] P. Mallet, I. Brihuega, S. Bose, M. M. Ugeda, J. M. Gomez-Rodriguez, K. Kern, and J.-Y. Veuillen, Role of pseudospin in quasiparticle interferences in epitaxial graphene probed by high-resolution scanning tunneling microscopy, Phys. Rev. B **86**, 045444 (2012).

[56] V. Kaladzhyan, F. Joucken, Z. Ge, E. A. Quezada-Lopez, T. Taniguchi, K. Watanabe, J. Velasco, Jr., and C. Bena, Quasiparticle interference patterns in bilayer graphene with trigonal warping, Phys. Rev. B **104**, 235425 (2021).

[57] F. Joucken, Jr., C. Bena, Z. Ge, E. Quezada-Lopez, S. Pinon, V. Kaladzhyan, T. Taniguchi, K. Watanabe, A. Ferreira, and J. Velasco, Jr., Direct visualization of native defects in graphite

and their effect on the electronic properties of Bernal-stacked bilayer graphene, Nano Lett. **21**, 7100 (2021).

[58] 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).

[59] L. M. Malard, J. Nilsson, D. C. Elias, J. C. Brant, F. Plentz, E. S. Alves, A. H. Castro Neto, and M. A. Pimenta, Probing the electronic structure of bilayer graphene by Raman scattering, Phys. Rev. B **76**, 201401(R) (2007).

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

[61] M. Skvortsov and M. Feigelman, Subgap states in disordered superconductors, J. Exp. Theor. Phys. **117**, 487 (2013).

[62] M. I. Salkola, A. V. Balatsky, and D. J. Scalapino, Theory of scanning tunneling microscopy probe of impurity states in a *D*-wave superconductor, Phys. Rev. Lett. **77**, 1841 (1996).

[63] 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).

[64] 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).

[65] T. Pereg-Barnea and M. Franz, Magnetic-field dependence of quasiparticle interference peaks in a *d*-wave superconductor with weak disorder, Phys. Rev. B **78**, 020509(R) (2008).