Accessing topological superconductivity via a combined STM and renormalization group analysis

The search for topological superconductors has recently become a key issue in condensed matter physics, because of their possible relevance to provide a platform for Majorana bound states, non-Abelian statistics, and fault-tolerant quantum computing. We propose a new scheme which links as directly as possible the experimental search to a material-based microscopic theory for topological superconductivity. For this, the analysis of scanning tunneling microscopy, which typically uses a phenomenological ansatz for the superconductor gap functions, is elevated to a theory, where a multi-orbital functional renormalization group analysis allows for an unbiased microscopic determination of the material-dependent pairing potentials. The combined approach is highlighted for paradigmatic hexagonal systems, such as doped graphene and water-intercalated sodium cobaltates, where lattice symmetry and electronic correlations yield a propensity for a chiral singlet topological superconductor state. We demonstrate that our microscopic material-oriented procedure is necessary to uniquely resolve a topological superconductor state.


I. INTRODUCTION
The search for topological states of matter has recently generated a flurry of broad activity in the field of superconductivity, both experimentally and theoretically (for several paradigmatic directions, see e.g. Refs. [1][2][3][4][5][6][7][8]. A topological superconductor (SC) is an unprecedented state of quantum matter, which possesses a full pairing gap in the bulk but gapless exotic surface states if a surface termination exists, as well as non-trivial vortex bound states 2,5,9 . A chiral superconductor 1,2,10 with broken time-reversal symmetry (TRS) may be considered the superconducting analogue of the quantum Hall phase as characterized by a non-trival Chern number of its Bogoliubov bands 11 , whereas a topological superconductor conserving TRS 6 is closely related to the quantum spin Hall phase along with its non-trivial Z 2 invariant 12 . Recently, chiral SCs have enjoyed significant attention, exhibiting a variety of exotic phenomena based on their non-trivial topology 1 , such as hosting Majorana vortex bound states 2,9 and gapless chiral edge modes, that carry quantized thermal or spin currents (see e.g. Ref. 13 ). The Majorana bound states can be interpreted as elusive fermionic "particles" equivalent to their own "antiparticles", and have potential applications in fault-tolerant topological quantum computation 14 .
In view of these striking properties, it would be desirable to develop guidelines to identify materials with the potential to host chiral topological SC states. As it turns out, the interplay of lattice symmetry, the shape and the size of the Fermi surface (fermiology), multi-orbital character, and electron-electron (e-e) interaction are decisive for an unconventional chiral pairing mechanism.
Material-specific research into this direction first concentrated on the perovskite Sr 2 RuO 4 , where experimental evidence points to a chiral odd-parity p-wave SC state 15 . Thus, Sr 2 RuO 4 has become a hallmark candidate for a quasi-two-dimensional electronic analogue of superfluid 3 He 16 . Due to the square-lattice symmetry in this quasi-two-dimensional material, a chiral pwave order parameter can be induced already at the instability level, lifting the degeneracy of the p x -and p y -wave SC channels: the chiral ("p x + ip y ")-wave order parameter combination, which spontaneously breaks TRS, then maximizes the condensation energy of the SC state below T c (for recent studies see e.g. 17,18 ). However, the topologically protected Majorana edge modes have -so far -not uniquely been identified, despite strong experimental efforts 19 . It suggests that the odd-parity paring, along with its non-trivial spin dependence, induces challenges which in terms of complexity even overshadow the original task, i.e. to identify a material with topological chiral SC.
It is hence one central strategic suggestion of this work to focus the attention on evenparity topological chiral SC. This happens for a "d+id" topological state, where the Chern number for the Bogoliubov bands is given by C = ±2. While it is a bit more subtle to obtain unambiguous Majorana bound states in the "d+id" vortex cores as compared to the "p+ip" case, this should still be possible by the application of a weak Zeeman field that does not affect the superconducting gap 20 . In the absence of spin-orbit coupling and in the presence of SU(2) spin invariance, this state corresponds to class C in the periodic table of topological insulators and superconductors [21][22][23] . The difficulty in realizing such a (d + id) singlet state is that generic fermiology and interactions, which might yield a d-wave state on a square lattice where most unconventional superconductors are found, favor d x 2 −y 2 over d xy pairing. This changes for hexagonal lattices, where the SC instability for d x 2 −y 2 and d xy is degenerate, as both form one two-dimensional irreducible lattice representation. Similar to the "p+ip" case for the square lattice, it is hence conceivable that there is a natural propensity for chiral singlet superconductivity on a hexagonal lattice in order to maximize condensation energy 24 .
In contrast to "p+ip", however, the singlet character of the unconventional pairing should make its emergence more generic, as it stems from electron-mediated pairing where large-q particle-hole fluctuations tend to drive singlet superconductivity. This has been recently addressed in several theoretical scenarios of a d+id state such as for doped graphene 25-31 , waterintercalated sodium cobaltates 32 , and the pnictide SC SrPtAs 33 . ARPES data on chemicallydoped graphene shows that the highly-doped regime close to the van Hove singularity, albeit with a currently high degree of disorder, is accessible 34 . Still, superconductivity in graphene has not yet been experimentally confirmed. From this perspective, water-intercalated sodium cobaltates and the pnictide SrPtAs may be more promising because, in both of these compounds, superconductivity has been already discovered 35,36 . Furthermore, some indications of unconventional superconductivity were observed in Knight shift data on cobaltates 37 , as well as in muon-spin rotation/relaxation measurements and nuclear quadrupole experiments on SrPtAs 38,39 . While no unambiguous experimental confirmation of chiral d-wave SC for these materials exists so far, we believe that further experimental attention is certainly warranted and could lead to the first unambiguous identification of chiral topological SC.
Along this path, a characteristic challenge in the search for unconventional chiral SC is the pronounced competition between different orders, such as spin density-wave (SDW) and different SC orders. On the hexagonal lattice, there exists a pronounced competition between the TRS-broken ("d + id") SC state and an f-wave state with TRS 24,40 , as exemplified in Fig. 1.
This clearly calls for methods which are capable of distinguishing the competing channels at the instability level, i.e. at the low energies of order k B T C of the SC gap features. This is the main strength of renormalization-group methods, in particular, of the functional renormalization group (fRG) method (for reviews see e.g. Refs. 24,41 ). The fRG allows for a systematic transition (i.e. renormalization) from a " high-energy" starting Hamiltonian (obtained e.g. from a Density-Functional Theory (DFT) and corresponding calculation of multi-orbital e-e interaction parameters) down to the low-energy k B T C effective theory, where the SC channels can be resolved 24 .
The resulting pairing functions stemming from the fRG calculations are one crucial ingredient of our analysis, which is subsequently employed as a refined microscopic input for our spectroscopic analysis. For doped graphene, the fRG results are schematically summarized in the phase diagram of Fig. 1a. This phase diagram 30 plots the critical instability scale Λ C of the fRG results, which gives an upper bound for the transition temperature T C (see below) as a function of doping at, and away from the van-Hove singularity (vHs). At this vHs (orange area), chiral (d+id) pairing competes with, but wins over the spin-density wave channel (details for the underlying band structure and the Hamiltonian are summarized in the Methods Section).
Away from the vHS (blue area), Λ C ∼ T C drops and whether the (d+id) or the competing f-wave instability is preferred depends on the range of electron-electron interactions (see again Meth. Sect.). Our corresponding fRG results for the possible phases in doped cobaltates are displayed in Fig. 1b  Knowing the precise functional form of the pairing function is of fundamental importance to make contact with experimental signatures. As seen e.g. for the sodium cobaltate example, the "d+id" state might exhibit a strongly anisotropic gap function 32 , which of course significantly affects the spectroscopic measures. It is thus the main task of our work to elevate spectrosopic signatures to provide evidence in favor of a possible chiral topological SC state. To this, we combine here the microscopic theory, i.e. fRG method, with the theory of scanning tunneling spectroscopy (STM). As evidenced, for example, by The conventional STM theory has already been successful in a variety of situations, where it has been viewed as a phenomenology, assuming a certain ansatz for an anisotropic order parameter [42][43][44] . Employing a normal metal-insulator-anisotropic SC (N-I-S) junction set-up, the assumed, effective pairing potential, felt by the quasiparticles 42,43 , was introduced into the Blonder-Tinkham-Klapwijk (BTK) formula 45 , which is a much tested and reliable formula for the conductance and dI/dV-characteristic of N-I-S junctions. This phenomenological procedure gave valuable insight, if one had to distinguish between e. g. a "structureless" s-wave and an anisotropic, for example, d-wave SC state 42 . However, precisely in the chiral superconductors as mentioned above, competing anisotropic SC channels such as (d+id), (p+ip) and f-wave may appear. Thus, to set up a microscopic criterion and, simultaneously, a reliable test, an as much as possible unbiased microscopic theory is required as an "input" into the BTK formula, which is what we elaborate on in the following discourse. where SC has been observed experimentally.
In the second step the differential conductance of these quasi-2D SC is calculated using an N-I-S setup with the pairing potentials from the fRG calculation. Such a N-I-S junction, formulated for a δ-function barrier model 42 , is known to imitate very well an experimental STM setup. This is documented by a variety of applications, where the pairing potentials have been used as a phenomenological "ansatz" for distinguishing between different SC symmetry channels 19,[42][43][44] . The effective barrier height is represented by Z 0 = mH 2 k F N , where m, k F N and H denote the electron mass, the Fermi momentum on the normal side, and the barrier potential, respectively. Both in-plane and out-of-plane setups are considered, where in the first (second) case the STM lies in (is oriented perpendicular to) the SC plane. Solving, as in the previous phenomenological studies, the BdG equations with the appropriate boundary conditions, but now with the momentum dependent microscopic pairing potentials, the coefficients (probabilities) for Andreev reflection r A and normal reflection r N are obtained. They determine, via the BTK-formula 45 the conductance, i.e.
for the quasiparticle injection with energy E = −eV , where V is the bias voltage, and θ is the incident angle with respect to the interface (Fig. 2). More details for the calculation of the normalized differential conductance can be found in the Methods Sect. Here, it suffices to note that the conductance contains two distinct pairing potentials ∆ + and ∆ − . They correspond to the effective pairing potentials for transmitted electron-like quasiparticles (ELQ) and hole-like quasiparticles (HLQ), respectively, as shown in Fig. 2. The total conductance is given by integrating the angle-resolved conductance of Eq. (1) over all transverse momenta, which are the independent modes. Except in the low-barrier (small Z 0 -) limit, which is not of interest here (our choice of Z 0 = 5 corresponds to a high barrier), the conductance peak corresponds to the Andreev bound-state (ABS) energy level, which is formed at the edge of the superconductor (see Fig. 3 of Ref. 42 ).

Differential conductance spectra for the in-plane setup
In the previous sections a theory for the tunneling spectroscopy of a N-I-S junction, combined with the microscopic fRG derivation of the underlying pairing potentials, was proposed as a new and efficient tool for identifying, in particular, a chiral SC state with broken TRS. Revealing the corresponding edge states of topological superconductors is crucially important because various novel features, such as non-Abelian anyons and Majorana fermions can be directly seen in STM spectroscopy 3,19,46 .
The differential conductance for the in-plane setup is shown in polar plots (left parts of

Figs. 3 and 4)
, where the radial axis is the quasiparticle excitation energy normalised by the superconducting energy (gap) scale ∆, which is a common energy scale obtained from fRG.
α denotes the angle between the interface normal (n) and the k x -direction (see Fig. 2). The absolute value of the pairing potential is also shown in a polar plot in order to compare it to the differential conductance. For each phase, cross sections at three different angles α are given in the upper right panel. The chosen angles are indicated by black lines in the dI/dV characteristics.
More specifically, the differential conductance spectra, obtained in the in-plane setup simulating an STM experiment, are presented for the (d+id)-and f-pairing phases of graphene ( Fig. 3) and for water-intercalated sodium cobaltates (Fig. 4). Comparing the conductance spectra (dI/dV) for the (d+id)-pairing phase ( Fig. 3a) with the f-pairing phase (Fig. 3b), it becomes immediately obvious why the fRG + STM calculation is an ideal tool to distinguish a (d+id)-wave TRS-broken phase from a f-wave TRS conserving phase. In graphene, for the chiral (d+id) case, the polar dI/dV intensity plot displays, in addition to the outer structure (which maps the local density of states with the peaks appearing at the local pairing potentials -compare with the pairing potential in Fig. (3a)) an inner structure for energies smaller than the superconducting band gap, which is called "sunflower structure" in what follows. It is this inner colored "sunflower structure" which immediately signals the presence of a TRS-broken SC state. A similar, but even more complicated inner structure for the cobaltates (below the band gap E=6.7 ∆) is again the sign of broken TRS (see Fig. 4a).
In order to understand the physical content of these differential conductance spectra for the (d+id) and f-pairing phases in more detail, let us summarize what kind of physical quantity the tunneling spectroscopy is detecting. The energy level giving the conductance peak is determined by a quantum condition of the bound quasiparticles (QPs) in a "pseudo quantum well". This quantum well is formed by the N-I-S junction (Fig. 2), where an electron injected from the N-side is transmitted into the superconductor (S) ejecting an Andreev hole-like quasiparticle.
This hole-like quasi-particle with its wave vector k − F S scatters into an electron-like quasiparticle with k + F S after reflection from the I-S interface, with a corresponding change in the effective pair potentials from ∆ − = ∆(k − F S ) to ∆ + = ∆(k + F S ) at the insulator (I) (see 42 and, in particular, Fig. 3 in this reference). The bound states form at the interface and tunneling electrons in the N-I-S junction flow via these bound states. These bound states have been shown in earlier work by Tanaka and coworkers 47 to converge to the edge states of the superconductor (trivial or non-trivial) in the large barrier-height limit (which is considered here).
The analogy with the quasiparticles and their Andreev reflection in a "pseudo quantum well" is useful for a transparent physical understanding of the inner "sunflower structure" appearing in the STM spectra of Figs 3a and 4a, as discussed below. The quantum-well analogy has been suggested by Kashiwaya et. al. 42 and shows the equivalence of the bound QP condition of the N-I-S junction with that of QPs in a S-N-S structure (with thickness of N → 0 and no difference of the superconducting phase across the junction), in which the pair potentials of the two superconductors are ∆ + and ∆ − , respectively. The QPs in the pseudo quantum well (normal region N) are confined if their energies are less than the amplitudes of both pair potentials, i.e. E < min(|∆ + | , |∆ − |). The bound QPs travel along a closed path by repeating Andreev reflections at both N-S and S-N interfaces. This is, then, equivalent to the ∆ + and ∆ − scattering processes at the I-S interface, as schematically shown in Fig. 2.
The corresponding peak condition for the bound QPs, which was first reported by Kashiwaya et.al. 42 , is given in Eq. (7) of the Methods Section. Here, we consider its general properties, which help to understand the occurrence of the outer peaks and that of the inner "sunflower structure" in Figs. 3 and 4.
For ∆ − = ∆ + , the peak condition, i.e. Eq. (7), is fulfilled if the energy E of the injected particle is E = |∆ ± |. Consequently, a peak at the bandgap occurs, i.e. at the outer structure in Figs. 3, 4. We call this peak the local bandgap peak, since as shown in these figures, the bandgap value depends on the (local) polar angle φ. The polar angle is defined as φ = arctan ky(φ) kx(φ) , with k x (φ) = k F (φ) cos φ and k y (φ) = k F (φ) sin φ and is given for ELQ by φ + = θ S − α and for HLQ by φ − = π − θ S − α (see Fig. 2). α denotes the angle between the interface normal (n) and the k x -direction and θ S is the angle between the interface normal and the momentum of the ELQ in the superconductor.
The celebrated zero-energy Andreev bound state (ZEBS) with E = 0 occurs if the phase difference of the pair potentials ∆ + and ∆ − , i. e. Φ + − Φ − (where ∆ ± = |∆ ± | exp(iΦ ± )), in Eq. (7) is ±π. Φ ± simplifies to Φ ± = φ ± , if k F does not depend on φ and a d + id pair potential with equal mixing of d x 2 −y 2 and d xy on the square lattice is considered. The novel aspect of more harmonics and an angle dependent k F leads to a more complicated peak structure in the differential conductance curves (see Figs. 3, 4 and Supplementary Information). One possibility to fulfill Φ + − Φ − = ±π is ∆ + = −∆ − for purely real pair potentials. The resulting ZEBS are seen in our fRG + STM calculations for the f-pairing phase in Fig. 3b, as well as in Fig. 4b.
However, for a complex (i.e. "d+id") pairing potential, the condition ∆ + = −∆ − is not sufficient. This is exactly the situation encountered for the inner "sunflower structure" of the (d+id)-order parameter shown in Fig 3a: Since ∆ + = ∆ x 2 −y 2 (θ S − α) + i∆ xy (θ S − α) and ∆ − = ∆ x 2 −y 2 (π−θ S −α)+i∆ xy (π−θ S −α), the phase difference between the two pair potentials is not restricted to multiples of π. Thus, the peak position moves between 0 and min(|∆ + | , |∆ − |), depending on the relation of ∆ + , ∆ − and the angle α. This peak is also called double-split peak, since the zero energy conduction peak is split into two peaks positioned symmetrically at positive and negative finite energies. This confirms the usefulness of the analogy with the pseudo quantum well and implies that indeed the quasiparticles in the quantum well are only confined if their energies are less than the amplitudes of both pair potentials.
In contrast to the conductance curves of the d+id pairing phases, that reveal clear signatures of broken TRS, the conductance curves of the f pairing phases contain zero energy peaks, which are present due to conserved TRS. These zero energy peaks (ZEBS) are seen in the cuts.
Additionally, the inset in Fig. 3b shows a zoom for small quasiparticle energies, displaying the ZEBS in white. As discussed already above, these peaks originate form an antisymmetric pairing ∆ + = −∆ − , and are a signature of conserved TRS. Furthermore, the differential conductance seems to be "rotated" with respect to the pair potential (see Figs. 3b and 4b). The physical reason for this effect is a lack of inversion symmetry for the f order parameter with respect to the origin of the k x -k y plane (let us mention that this inversion is preserved for both real and imaginary parts of the d+id order parameter). This lack of inversion symmetry and the corresponding "rotation" of the f order parameter can be easily understood, considering for example graphene (Fig. 3b) at an angle α = π 6 . Then, ELQ exhibit a pair potential of ∆ + = ∆(θ S − π 6 ) and HLQ of ∆ − = ∆( 5π 6 − θ S ). Using the pair potential (the signs of the pair potential for a phase difference of π are opposite in Fig. 3b), we find ∆ + = −∆ − , giving rise to a zero energy peak, which is indeed found in the conductance spectrum at α = π 6 . Let us add a few more details, concerning the results in Figs. 3 and 4. Figure 3 shows In general, the height of the peaks depends on how many incident angles θ contribute to the resonance for a given α-direction. In the case of the cobaltates, the gap is very anisotropic and more harmonics contribute than for graphene. Consequently, the number of incident angles contributing to a resonance is smaller, giving a smaller peak height.
Summarizing, the in-plane setup is sensitive to the magnitude and the phase of the pair potential, and allows to distinguish the different pairing phases (d+id and f). While the f-wave paring phase gives zero energy peaks (ZEBS) typical for order parameters with conserved TRS, the signature of broken TRS can be clearly seen in the beautiful " sunflower structure" in the differential conductance spectrum of the d+id paring phase of our proto-type examples.
Differential conductance spectra for the out-of-plane setup Figure 5 shows the differential conductance obtained in the out-of-plane setup for the d + The "high-energy" starting point, i.e. the Hamiltonian has to be accurately determined, for the "band-structure" part H 0 (typically taken from a fit to an a-priori DFT calculation) and for the interaction H int (taken e.g. from a "cRPA" determination of the various terms). Details of the parameter choices can be found in Ref. 24 .
For the cobaltates, the Hamiltonian includes three hybridized orbitals per site (d xy , d yz , d zx ) and reads where c † iασ denotes the electron creation operator with spin σ =↑, ↓ and orbital α at site i. The occupation number is defined as n iασ = c † iασ c iασ . In addition, t represents the hopping mediated by O pz orbitals and t' corresponds to a direct Co-Co-hopping, D is the crystal-field splitting, and µ the chemical potential. These parameters are set to t = 0.1eV , and interorbital Coulomb interactions, respectively. The remaining interaction parameters are J H = J p = 0.07eV for Hund's rule coupling J H and pair hopping J p .
In graphene, the tight-binding Hamiltonian H 0 is where n = i,σ n i,σ = i,σ c † i,σ c i,σ . c † i,σ is the creation operator for an electron with spin σ at site i, µ denotes the chemical potential, and t 1···3 is the hopping strength for nearest neighbour (1), second nearest neighbour (2) and third nearest neighbour (3) hopping. Coulomb interaction is included by a long-range Hubbard-type Hamiltonian H int with where U 0···2 gives the Coulomb repulsion scale from on-site (0) to second nearest neighbour (2) interactions, respectively.
The near-degeneracy between SC and density-wave orders is strongly influenced by a subtle interplay between deviations from perfect nesting (taken into account in H 0 of Eqs. The 4PF V Λ (k, −k, q, −q) in the Cooper channel is, then, decomposed into different eigenmode contributions 24 .
where i is a symmetry decomposition index. The leading instability of that channel corresponds to an eigenvalue w SC i (Λ) first diverging under the flow of Λ. f SC i (k) is the SC form factor of pairing mode i, which tells us about the SC pairing symmetry and hence gap structure associated with it. In the fRG, from the final Cooper channel 4PFs, this quantity is computed along the discretized FSs. f SC i (k) is the gap function, which then enters the BTK formula (1) for the conductance. Details of the functional RG-procedure can be found in Ref. 24 .
For the STM-part of the calculation, we use the NIS junction setup and solve the Bogoliubovde Gennes equations for the normal and the superconducting parts. The insulator is modelled by a delta-Dirac potential barrier with strength H 42 . We apply the approximation, that the quasiparticle excitation energy E and the maximum absolute value of the pairing potential max{|∆|} are much smaller than the Fermi energy E F . Consequently, the pairing is only relevant close to the Fermi surface. Within these approximations, the wave function obtained by solving the BdG equations (see Supplementary Information) is independent of the dispersion relation of the material (we assume to have a quadratic term in the dispersion). We use the continuity of the wave function at the interface and the matching of the derivative of the wave functions of the normal and the superconducting side to obtain the coefficients for Andreev (r A ) and normal (r N ) reflection. The total conductance of the NIS junction is obtained by integrating the BTK conductance over all independent contributions, which is integrating over all k y in our in-plane setup. We normalise the conductance by the conductance of a NIN junction in the same geometrical setup.
A peak in the conductance spectrum is obtained, if the following condition, first reported by Kashiwaya et. al. 42 is fulfilled where , Φ ± = arg(∆ ± ) and E denotes the quasiparticle excitation energy.
In general, Φ ± = φ ± because k F depends on φ. However, for the generic d + id pair potential, Supplementary Information: Accessing topological superconductivity via a combined STM and renormalization group analysis Normalised differential conductance The wave functions for the normal part (N) and the superconducting part (S) are obtained solving the Bogoliubov-de Gennes Hamiltonian.
and Φ ± = arg (∆ ± ). ∆ + (∆ − ) denotes the pairing potential exhibited by the electron-like (hole-like) quasiparticles (ELQ/HLQ), where we define ∆ ± = ∆(φ ± ). In the in-plane setup φ + = θ S − α and φ − = π − θ S − α. θ S denotes the angle of an ELQ with respect to the interface normal and α gives the mismatch between the interface normal and the k x -direction of the pairing potential (see also Fig. 2, main part). If the Fermi momenta of the normal region and the superconductor are identical, i.e. k F S = k F N , then θ S = θ N =θ. In the out-of plane setup φ + = φ − = φ, θ + = θ S and θ − = π − θ S , where the ± index refers to ELQ and HLQ, respectively. The insulator is modelled with a delta-Dirac barrier with effective, incident angle dependent, barrier height Z = Z 0 cosθ , where Z 0 = mH 2 k F N . We use the continuity of the wave function at the interface and the matching condition of the derivative of the wave functions of the normal and the superconducting side to obtain the coefficients for Andreev (r A ) and normal (r N ) reflection. The total conductance of the NIS junction is obtained by integrating the BTK conductance over all independent contributions, which is integrating over all k y in our in-plane setup. We normalise the conductance by the conductance of a NIN junction in the same geometrical setup. The normalised differential conductance in the in-plane setup σ in is In the out-of-plane setup, we have to integrate over all momenta lying in the plane of the superconductor. The normalised differential conductance in the out-of-plane setup σ out is Harmonics for the pairing potential of graphene with short range Coulomb interactions The following four harmonics were used to fit the gap for the d+id pairing phase of graphene with short range Coulomb interactions.
Fitting these basis functions to the numerical data gives an analytical expression for the pairing potential.
The fitting factors are c 1 = c 2 = 1.68 and c 3 = c 4 = 0.32. ∆ is an arbitrarily chosen superconducting energy scale, which is the same for all phases and allows to compare the energy scale of the phases. It is used for the normalisation of the energy in the plots. The same harmonics were used to obtain the expression for the long-range Coulomb interaction pairing phase. For the cobaltates, 12 harmonics were used to acchieve a fit. In the appendix of Ref. 24 a general scheme describing how to obtain the lattice harmonics is given.
Differential conductance for graphene with long-range Coulomb interactions Fig. 3 in the main part of the manuscript displays the in-plane setup differential conductance spectra for the short-range Coulomb interactions in graphene at the van-Hove singularity. and 4a of the main manuscript, the local bandgap peak occurs if the energy E of the injected particle is E = |∆ ± |. As before, this finding is in analogy to the direct mapping of the maximum value of the pairing potential in certain α-directions, as discussed for the short-range d + id case in the main text.
For the f-wave pairing at the same doping x = 0.15, we chose U 0 = 10eV , U 1 = 5eV and U 2 = 3eV (see Eq. (5) in Sect. IV of the main text). Again our general statement holds, i.e.
that the physics of the TRS-preserving f-pairing case, for the longer-range Coulomb situation, is crucially different from the chiral d + id case: In the f-wave scenario, we only find the outer local bandgap peaks (rotated by an angle π 6 in respect to the pairing potential) plus the zero energy bound state peaks.

3.
FIG. S1: Differential conductance spectra (dI/dV) for the in-plane setup and corresponding pairing potentials for the long-range Coulomb interaction pairing phase of graphene. The left panels show the differential conductance spectra for the d + id pairing phase (a) and the f pairing phase (b). The corresponding pairing potentials are shown on the right as a polar plot of the absolute value of the pairing potential. The differential conductance spectra, obtained in the in-plane setup simulating an STM experiment, are plotted as a function of the quasiparticle energy E (radial axis, normalised by the reference band gap ∆) and the angle α between the interface normal and the kx-direction of the pairing potential (polar axis). Differences in brightness of the colors are due to the plot style. Note, that the scale changes with the height due to the perspective.