High-resolution real-space evaluation of the self-energy operator of disordered lattices: Gade singularity, spin-orbit effects and p-wave superconductivity

Disorder is a key factor influencing the behavior of condensed states of matter, however the true extent of its impact is generally difficult to determine due to the prominent roles played by quantum interference, entanglement between spin and orbital degrees of freedom and proximity to quantum critical points. Here we show that the one-particle disorder self-energy -- a direct probe of the renormalization of low-energy excitations due to defects and impurities distributed randomly in a crystal -- can be obtained by means of unbiased spectral expansions of lattice Green's functions in a computationally expedient manner. Our scheme provides a powerful framework to map out the frequency and wavevector dependence of electronic excitations in unprecedented large tight-binding systems, up to $10^9$ orbitals, with energy resolution only limited by the mean level spacing. We demonstrate the versatility of our approach in 3 distinct problems: (i) the Gade singularity in honeycomb layers with dilute topological defects; (ii) the rich landscape of impurity resonances in a spin-orbit-coupled ferromagnet; and (iii) the tailoring of emergent s-wave and p-wave superconducting phases in graphene via atomic defects. These examples reveal rich features in the disorder self-energy $\Sigma_{\textrm{dis}}(\mathbf{k},\omega)$ that are absent from the self-consistent T-matrix approach and other common approximation schemes, which include regimes of nontrivial wavevector dependence and anomalous dependence upon the impurity concentration. Our study unravels puzzling, and so far largely inaccessible, manifestations of strong nonperturbative quantum interference effects in quantum materials and disordered phases of matter.


I. INTRODUCTION
Disorder plays an essential role in the wealth of phenomena observed in condensed matter. On the one hand, disorder influences the structural, optical and transport properties of metals and semiconductors, as well as the phase stability of unconventional states of matter, such as chiral supercondutors [1] and symmetry-protected topological insulators [2]. On the other hand, certain types of disorder can induce interesting behavior unseen in clean systems, ranging from the breakdown of the Fermi liquid description and quantum interference effects in mesoscopic conductors to the Kondo effect and novel quantum phases in strongly correlated systems [3][4][5][6][7][8].
Green's functions provide a powerful mathematical device to study the impact of disorder and correlations in many-body quantum systems. Of paramount importance in this context is the single-particle irreducible self-energy which, in effect, dresses bare Green's functions with the cooperative effects experienced by quasiparticles, and therefore determines the key features of the spectral function measured in angle-resolved photoemission experiments [9,10]. The disorder contribution to the self-energy in a homogeneous disordered system is defined in terms of disorder-free,Ĝ 0 (k, ω), and full, G(k, ω), single-particle Green's functions asΣ dis (k, ω) = G −1 0 (k, ω) − Ĝ (k, ω) −1 (here ... indicates configurational ensemble average) and contains detailed information on the effects of quenched disorder, such as the existence of impurity bound states, as well as the precise extent of quasiparticle renormalization induced by single impurity scattering events and quantum interference effects. Because the self-energy is intimately connected to vertex functions in diagrammatic theory [11], it is also an indispensable tool in the study of quantum transport phenomena; most notably as a means to obtain conserving approximations in self-consistent field-theoretical calculations [12].
The self-energy from electron-impurity interactions, as is presented in textbooks, is usually described in terms of a complex scalar,Σ dis (k, ω) = ∆(k, ω) ∓ i/[2τ (k, ω)], where the real part is responsible for the effective-mass renormalization and τ (k, ω) is the quasiparticle elastic scattering lifetime (here the signs ∓ hold for retarded/advanced Green's functions as required by analyticity). Closed analytical expressions for the singleparticle self-energy can be obtained by means of the diagrammatic technique based on perturbative expansions in the semiclassical parameter g ∝ (ωτ ) −1 [13]. The classical example is a Fermi gas with a dilute concentration of impurities, for which a lowest-order "rainbow" diagram calculation yields the relaxation time τ (ω) ∝ [n u 2 0 ν(ω)] −1 , with n the impurity density, u 0 the potential strength and ν(ω) the bare density of states. The weak-disorder approximation can be improved systematically by partial resummation of infinite series of diagrams (e.g., by means of the T -matrix or coherent potential approximations) and extended to multi-orbital Hamiltonians for both noninteracting and interacting cases; for a recent review see Ref. [14].
Despite its notable successes, the diagrammatic approach suffers from a major drawback: the topological complexity of scattering diagrams increases quickly with the perturbation order, which presents a formidable barrier to our understanding of nontrivial manifestations of disorder beyond the weak-coupling regime. Indeed, a range of intriguing phenomena triggered by strong (nonperturbative) disorder effects are currently the focus of intense investigation. Some examples include strong Anderson localization [15][16][17][18][19], rare region effects in three-dimensional topological semimetals [20][21][22][23] and frozen multifractality in chiral-symmmetric lattices [24][25][26][27][28], whose analysis has defied even the most advanced field-theoretic approaches. At the heart of such unusual phenomena is the nonperturbative accumulation of quantum coherent scattering processes, whose satisfactory description calls for the use of large-scale numerical approaches.
Here, we report a comprehensive numerical study of strong-coupling effects on the disorder self-energy in several electronic phases of matter. Our study is based on a new high-resolution real-space spectral method that gives access to previously unexplored features of the disorder self-energy, unveiling its rich internal matrix structure and full k-and ω-dependences. Our calculations, performed on very large systems, disclose several important features of the quasiparticle self-energy that are absent in the standard perturbative treatments, including surprisingly strong k-dependence generated by nonlocal correlations (quantum interference) and the emergence of offdiagonal self-energy elements near quantum criticality.
In what follows, we introduce the new disorder selfenergy framework, outline its favorable scaling properties and present its application to three distinct problems: (i) quantum criticality induced by chiral disorder in the BDI symmetry class; (ii) scattering resonances in a spin-orbit coupled ferromagnet; and (iii) disorder-enhanced p-wave superconductivity in graphene.

A. Chebyshev expansion of the spectral function
To set the stage, let us briefly review the polynomial expansion of the one-particle spectral function. Consider a general fermionic system on a d-dimensional lattice described by a Hamiltonian with a bounded spectrum,Ĥ. To enable a decomposition of the spectral function in terms of orthogonal polynomials, we first perform the following linear transformation where 1 is the identity operator defined on the Hilbert space of the lattice and E a(b) indicates the largest (smallest) eigenvalue ofĤ. Note that this procedure maps the eigenvalues of the Hamiltonian onto the canonical interval I = [ −1, 1]. Here, we employ Chebyshev polynomials of the first kind which are particularly well-suited to approximate nonperiodic functions over a finite interval on the real axis [29]. The spectral operatorÂ(ε) = δ(ε −ĥ) associated with the rescaled Hamiltonian in Eq. (1) can be decomposed into a Chebyshev series according tô where ω n (ε) = 2/[π(1 + δ n,0 )(1 − ε 2 ) 1/2 ] is the weight function entering in the orthogonality relations and ε ≡ (E − E + )/E − is the rescaled energy variable, with E the energy spectrum of the original Hamiltonian [30]. The favorable convergence properties of the generalpurpose Chebyshev expansion in Eq. (2) reflects its close relation to the Fourier series made manifest by the identity T n (ε) = cos(n arccos ε). We do not provide details on the underlying spectral theory but instead refer the reader to Boyd's book [29]. By virtue of the Chebyshev recurrence relations, the matrix coefficients in Eq. (2) can be obtained iteratively to any desired order. In a practical implementation, it suffices to evaluate the scalar Chebyshev expansion coefficients of the overlap A φψ (ε) = φ|Â(ε)|ψ . The choice of basis functions φ, ψ defines the target function of energy (see below). Because the recursive procedure [Eqs.
(3)-(4)] is highly stable, the reconstruction of the spectral function can be carried out in principle with very high energy resolution γ ≈ 1/M (in rescaled units), where M is the truncation order. Such features will be of key practical importance in the evaluation of the selfenergy operator, as we shall see briefly.
To illustrate the effectiveness of the spectral approach, we first consider the local density of states (LDOS) at a site i, defined as ν i (ε) = i|Â(ε)|i . The LDOS is the simplest spectral target function of energy derived from the spectral operator [note that it corresponds to a diagonal element ν i (ε) = A ii (ε) in the lattice basis], yet it contains rich information on electronic state hybridization and scattering processes. From Eq. (2), the M -order Chebyshev approximation to the LDOS is easily constructed as with µ in = i|T n (ĥ)|i [31].
To retrieve the Chebyshev coefficients {µ in }, the Eqs. (3)-(4) are iterated on the fly by exploiting the sparseness of minimal-basis (real-space) representations. The problem is thus solved efficiently by repeating the following (ii) large-scale evaluation of Chebyshev moments of the Hamiltonian matrix; and (iii) high-resolution spectral reconstruction of self-energy operator using an exact Chebyshev decomposition of the lattice Green's function. For problems that require a self-consistent mean-field approach, the procedure is repeated until convergence to the desired accuracy is achieved.
steps. (i) Starting with the initial vectors |ψ i 0 := |i and |ψ i 1 :=ĥ|ψ i 0 , act iteratively with the rescaled Hamiltonian using the Chebyshev rule |ψ i n+1 = 2ĥ|ψ i n − |ψ i n−1 ; and (ii) At each step compute the overlaps µ ni = ψ i 0 |ψ i n . From the knowledge of the Chebyshev moments {µ ni }, the LDOS in any desired energy range can be easily retrieved using Eq. (5). The overall computational cost is determined by the energy resolution desired for the spectral reconstruction of the LDOS. For typical sparse Hamiltonian matrices, the number of operations scales linearly with respect to both the number of sites and number of moments M , which makes the method particularly advantageous for single-electron problems [32]. The spectral approach has been successfully employed to reveal the electronic structure of a wide range of systems, including Anderson disorder models [33], graphene with atomic defects [34][35][36] and disordered superconductornormal metal interfaces [37]. Other recent applications include the calculation of time-dependent equilibrium Green's functions of superconductors [38] and dynamical structure factors in quantum spin chains [39]. (For a review of early work, see Ref. [30].)

B. High-resolution self-energy calculation
Having laid out the basic principles underlying the efficient reconstruction of the spectral function using Chebyshev polynomials, we now move on to tackle the nontrivial problem of determining the self-energy operator. The central objects of interest in this discussion are the re-tarded lattice Green's functions and disorder self-energŷ whereĤ 0 is the clean Hamiltonian and η plays the role of an energy resolution (see below). As customary, the real-space disorder is added to the Hamiltonian, H =Ĥ 0 +V dis , via random modifications of hopping amplitudes and on-site energies. Of particular interest is the wavevector (k) dependence and internal (orbital and spin) structure ofΣ η (ω). It is natural to ask whether the accurate LDOS polynomial scheme can be extended to the matrix inverse problem posed by Eq. (7). An immediate stumbling block is simply that the inversion procedure is prone to loss of accuracy, particularly for weak disorder, due to the parametrically small difference between the clean and disordered Green's functions. As shown below, such a hurdle can be overcome by devising a spectral algorithm that accurately reconstructs the Green's function projected onto the local basis elements, thereby effectively mapping the problem into an N × N matrix inversion, where N is the number of bands. A more subtle issue concerns the stability of the self-energy with respect to the broadening scheme of the Green's function of finite systems [40,41]. For many problems of interest, one must be able resolve the fine structure of the disordered Green's function and thus η must be comparable or smaller than the self-energy itself. In practical terms, η is bounded from below by the mean level spacing and a careful convergence analysis is required to obtain sensible thermodynamic results. We posit that the Chebyshev polynomial-based spectral approach is well suited to overcome the above ob-stacles, since it allows for real-space calculations with high accuracy and well-defined energy resolution. The standard strategy to mimic the broadening in Eq. (6) is to formally expand the Green's functionĜ (ω) = lim η→0 +Ĝ η (ω) in Chebyshev polynomials and regularize the resulting spectral series through convolution with a Lorentzian kernel [42,43]. Here we propose an alternative approach based on the direct expansion of the broadened Green's functionĜ η (ω) into Chebyshev polynomials of the first kind.
Let us start by defining the rescaled Green's function operator where we have used Eq. (1) and defined γ = η/E − . This rescaled Green's function then admits the following exact decomposition where a n = 2/(1 + δ n,0 ) and f (z) = i √ 1 − z 2 [44,45]. The main advantage of this variant of the familiar kernel polynomial method is that the energy levels are probed with known uniform resolution over the entire spectral range ε ∈ [−1, 1]. The energy resolution can be made as high as desired by a judicious truncation of Eq. (9); as a rule of thumb, the required number of Chebyshev iterations is M γ ≈ [1/γ]. This spectral scheme, combined with an efficient real-space implementation, will allow us to resolve the fine structure of the quasiparticle self-energy in large systems containing multi billions of orbitals, a task that has remained elusive thus far. Such a capability is key to exploring topological transitions and disordered systems at quantum criticality, where spectral convergence is already challenging at the level of considerably simpler average density of states [22,27,45].
With regards to the disorder averaging procedure [Eq. (7)], a brief discussion is in order. For the cases of interest here, the k-space Green's function is found to exhibit self-averaging behavior. Thus, the disorder self-energy of a single large sample is representative of the whole ensemble, greatly reducing the computational cost related to configurational averaging. In practical terms, the self-averaging property of the disorder self-energy is demonstrated numerically by analyzing the scaling of self-energy fluctuations with the system size. For the interested reader, we provide analytical proofs for two common classes of problems in Appendix V A. Our findings suggest that self-averaging is a universal property of the k-space disorder self-energy, which will be explored in future work.
To facilitate the evaluation of the disorder self-energy, let us introduce the orthogonal basis set of plane-wave states of wavevector k, {|k, α }, where α = 1...N labels the internal quantum numbers of the electronic system. We first evaluate the Green's function matrix elements, G αβ (k, z) = k, α|Ĝ (z) |k, β , with the desired spectral resolution. The projected Green's function Chebyshev moments, µ n αβ (k) = k, α|T n (ĥ)|k, β , are then computed with the recursive scheme described in Sec. II A. Specifically, we compute the N × N overlap matrix µ n αβ (k) = k, α|ψ βk n (n = 0, ..., M γ − 1) along the desired k paths, with |ψ n βk the n-th Chebyshev vector obtained by successive applications of Eqs. (3)-(4) with the initial vector |ψ 0 βk = |k, β . This will afford us with several computational advantages, in particular, the possibility to tackle very large systems (note that only 3 vectors of the Hilbert space dimension need to be stored for each k point). Next, the N × N k-space Green's function is reconstructed using the truncated expansion With a suitable choice of M γ , Eq. (11) yields numerically exact results for the broadened lattice Green's function within machine precision. Furthermore, the dependence of the Green's function with the resolution parameter (up to γ ≈ 1/M γ ) can be retrieved without having to recalculate the moments µ n αβ (k), which allows for assessing convergence on the fly. The momentum-space Green's function of the clean system can be calculated along the same lines (or simply through direct diagonalization of the corresponding Bloch Hamiltonian). The final step is a simple inversion of two N × N matrices, which reconstructs the momentum-space self-energy in a fully non-perturbative manner. This simple inversion of the lattice Green's function G η is justified because, for self-averaging systems, it is diagonal in k-space (see Appendix V A for a detailed discussion). The procedure may, of course, be repeated for any number of samples in order to recover the disorder-averaged Green's function in the traditional sense if self-averaging is not obvious [note that, in general, Σ η and G η in Eq. (12) should be understood as the disorder-averaged operators Σ η = Σ η and G η = G η .]. Our self-energy spectral framework is summarized in Fig. 1. Before closing this subsection, we briefly comment on the computational cost associated with the calculation of the self-energy. Because the recursive algorithm exhibits polynomial complexity [see Eq. (10) and text therein], one can extract the self-energy operator of considerably large systems with a small computational cost. For example, a single-orbital tight-binding model on a squarelattice with real first-neighbor hoping, 10 3 × 10 3 sites Significantly more complex problems can be tackled by optimizing the parallelization efficiency of matrix-vector multiplications using an adaptive real-space domain decomposition algorithm as implemented in the KITE package [32] (see Appendix V B for more details). Such a strategy was adopted in recent works reporting accurate studies of the effect of short-range disorder on the nodal density of states of topological semimetals [22,46].

A. Gade singularity of graphene: k-dependence
We now turn to the presentation of the results obtained with the spectral method introduced above. As a first case study, we consider graphene with vacancy defects [47][48][49]. Site dilution in a honeycomb layer provides an intriguing example of a random-hopping system which exhibits anomalous quantum critical behavior stemming from sublattice (chiral) symmetry [50,51]. Chiral-symmetric disorder induces critically delocalized states at the band center [so-called zero energy modes (ZEMs)] characterized by a Gade singularity in the density of states [24][25][26][27][28]. Moreover, quantum transport simulations indicate that dilute ZEMs can overcome Anderson localization in an infinite system, with conductivity pinned to the T = 0 ballistic value σ 0 = (4/π)e 2 /h irrespective of the vacancy concentration [45]. The puzzling behavior of ZEMs has been attributed to unusually strong nonperturbative quantum interference effects (beyond standard field-theoretic treatments), but direct evidence in k-space has remained elusive. In order to model the electronic structure of graphene with randomly distributed vacancies, we resort to a minimal tight-binding model for spinless fermions on the honeycomb latticê where c † i (c i ) adds (removes) an electron to the i-th site and i, j denotes nearest-neighbor pairs of sites. Furhermore, t ij = t if i and j are both undiluted sites, otherwise t ij = 0. The model in Eq. (13) possesses both timereversal (T 2 = +1) and particle-hole (C 2 = +1) symmetries, placing it in the chiral orthogonal (BDI) class of the topological classification [52]. Thus, random vacancies can be viewed as topological point defects which preserve the underlying non-spatial symmetries of the host crystal, most notably its chiral-sublattice symmetry, i.e. SĤS = −Ĥ. Here, S ≡ T C = σ z , with σ z the diagonal Pauli matrix defined in the A-B sublattice space.
Previous work has calculated the vacancy-induced selfenergy in the continuum limit of Eq. (13) by means of the self-consistent T -matrix (SCTM) approximation [53]. The SCTM can be evaluated analytically both near the band center and far from it, yielding a scalar self-energy of the form Im with c the vacancy concentration and Λ a suitable ultraviolet cutoff [53]. Although the SCTM provides a faithful description of the problem in the semiclassical regime (ω Γ), it cannot reproduce neither the Gade singularity in the density of states [27], nor the anomalous behavior in the conductivity of ZEMs [45], because it neglects quantum coherent multiple scattering. Moreover, the continuum model treats vacancies as δ-peak potentials centered at random positions, which results in a structureless (i.e. k-independent) self-energy operator at all energies. To overcome these limitations, we use our high-resolution spectral approach to map out the momentum-space self-energy of the lattice model. For definiteness, we focus on the case of compensated vacancies equally distributed on both sublattices. The Chebyshev-polynomial-based reconstruction of the selfenergy is carried out on very large systems with up to 10 9 sites and M = 2 16 = 65536 moments, giving us unprecedented access to sub-meV resolution over the entire (ω,k)-parameter space. The disorder selfenergy is projected onto the sublattice space according to Σ η αβ (k, ω) = k, α|Σ η (ω)|k, β , with α, β = A, B. Homogeneity implies Σ AA = Σ BB and Σ AB = Σ BA . Thus it suffices to find the AA and AB elements. Furthermore, particle-hole symmetry implies Σ AA (ω) = −Σ * AA (−ω) and Σ AB (ω) = Σ * AB (−ω). The self-averaging property combined with the very large lattice sizes provides the additional advantage that a single disorder landscape is required for our purposes (see Appendices V A and V C for additional details on the self-averaging property and spectral convergence, respectively).
We focus the subsequent discussions on the imaginary part of the self-energy operator, which encodes the quasiparticle lifetime. As a reference point, we calculate the Tmatrix and SCTM self-energies using the lattice Green's function of the clean model. The fully converged results, summarized in Fig. 2, contain a number of surprising findings. First, the disorder self-energy shows a strong kdependence in the vicinity of the band center, where the Gade singularity is located. Moreover, the point defects endow the self-energy with a nonzero off-diagonal component in that same region. Second, the concentration dependence of the self-energy exhibits anomalous scaling at the lowest energies, where a rich crossover between the quantum critical regime near ω = 0 and the pure semiclassical regime at high energies can be seen. All these features are missing from the SCTM approximation and, as argued below, provide fingerprints of the conjectured strong nonperturbative quantum interference effects induced by ZEMs in graphene. A close up of the self-energy matrix elements around the Gade singularity are shown in Figs. 2(a)-(d). The observed fine structure is confined to a narrow window of width δ ≈ 2tc 0.6 . The twin peaks in the diagonal elements Σ AA(BB) borne out by our high-resolution data become visibly sharper as one moves away from the K point along the path indicated in Fig. 2(e), with other paths showing similar behavior. Let us note that the distance between these peaks decreases, while their height increases, as we move away from the K point, such that the curves approach that of the T -matrix approximation. However, whereas the latter is structureless and diverges as Im Σ AA ∼ −1/ |ω| log |ω|, the numerically exact selfenergy is strongly k-dependent and bounded in the vicinity of the Dirac (K) point. The SCTM approximation effectively smears the bare T -matrix result (thus removing the divergence at ω = 0), however, neither approach captures the nonperturbative behavior seen at low energy. Strikingly, the ω = 0 self-energy at the Dirac point approaches zero as η → 0 at all concentrations, which implies a divergent elastic scattering time at the K point (see Appendix V C for a scaling analysis). We speculate that the exceedingly large quasiparticle lifetime protects ZEMs against backscattering, thus providing a new insight into the "mysterious" ZEM resilience observed in large-scale simulations of the dc conductivity [45]. We stress that the favorable scaling of our method is crucial to uncover the fine features of the self-energy, a task which requires both very large samples and fine resolution η δ.
We now briefly discuss the sublattice-coherence effects encoded in the off-diagonal elements of the k-space self-energy, Σ AB(BA) . Figure 2(d) shows that the offdiagonal component acquires the form of a symmetric Fano resonance with sizeable amplitude away from the Dirac point [e.g., ImΣ AB ≈ 0.3t for c = 0.003 and k = k * 0.76|K|x; see Fig. 2(e)]. Within the T -matrix or SCTM approach, the off-diagonal self-energy is identically zero, which shows that the AB-sublattice coherence captured by our real-space spectral method is inherently a nonperturbative effect resulting from coherent multiimpurity scattering events.
In order to unveil the extent of quantum coherent scattering effects, we extracted the vacancy concentration dependence of the self-energy over the entire spectrum. First, we note that the self-energy converges to its bare T -matrix form in the high electronic density regime, where chiral symmetry is absent and the system falls under the standard orthogonal Wigner-Dyson class; see  the energy Γ k < Γ T at which the k-dependency becomes very strong, which scales roughly as c 0.6 . For energies |ω| Γ k , the self-energy is k-independent and, for larger energies still, when |ω| Γ T , it matches the simple Tmatrix expression Σ(ω) ≈ −ic/ |ω| log |ω|. The linear dependence upon the concentration is the signature of the semiclassical regime, where single-impurity scattering events dominate. However, much more interesting is the vicinity of the BDI quantum critical point, where the perturbative picture breaks down [24][25][26][27][28]45]. Near the Dirac nodes (ω → 0), the self-energy becomes independent of the concentration. A close up of the concentration dependence at the Dirac point is shown in Fig. 2(f). This anomalous behavior is found to occur for energies that are within the energy window where the k dependence is the strongest, which indicates that both effects have origin in high order (multiple impurity) scattering processes. We note that the onset energy of this anomalous behaviour increases with the increase of vacancy concentration. This means that the manifestations of quantum criticality highlighted here can be pushed towards experimentally accessible energy scales in sufficiently disordered samples. Of equal interest is the intermediate energy regime [Fig. 2(g)], where considerable deviations from the semiclassical picture are also apparent. Here, the K-point self-energy is found to follow the scaling law Im Σ AA (0, ω) = c α f (ωc α ), see Fig. (2d), where f (x) is independent on c and α = 0.56 ± 0.02.

B. Spin-orbit effects in ferromagnets: SrRuO3
Next, we illustrate the versatility of our approach by computing the disorder self-energy matrix of a spinorbit-coupled ferromagnetic metal. The model system chosen for this study is the itinerant ferromagnet SrRuO 3 (SRO) [54], a well-studied oxide material that features, among other things, momentum-space monopoles of Berry curvature and interface-driven chiral spin textures [55][56][57][58][59]. Our focus here is on the spin-polarized twodimensional electron gases that are formed in SrRuO 3 embedded in a SrTiO 3 matrix [60]. Complex-oxide superlattices have attracted widespread interest because they provide a platform to engineer metallic states with unusual ferroelectric properties [61][62][63][64]. However, the effect of disorder upon their complex interfacial behavior remains largely unexplored.
To model an emergent spin-polarized electron gas in a SrRuO 3 /SrTiO 3 superlattice (which is known to be confined to the 4d orbitals of the Ru atoms in the SRO layers [65]), we employ a first-principles parameterized multiorbital TB model as developed in Ref. [58]. The details of the TB model, which accurately describes the predominant Ru t 2g -bands near the Fermi level and also includes spin-orbit coupling (SOC), are provided in Appendix V D 1. Our primary aim is to understand whether the self-energy operator acquires a nontrivial matrix structure. To this end, we investigate the quasiparticle self-energy generated by on-site disorder and point-like defects (vacancies). In the former, the random on-site energies within each unit cell in real space {ε i } are taken from a box distribution where W defines the disorder strength. With this choice, the on-site disorder is locally correlated since all the 4dorbitals inside a unit cell experience the same potential. Different choices are possible, but this simple prescription will be sufficient to illustrate the nontrivial role of disorder in this class of oxide materials. We first briefly describe the electronic structure of the spin-polarized gas formed at the interfacial layers of SRO superlattices. Figure 3(a) shows the average density of states and Fig. 3(b) shows the band structure along the path ΓXM Γ. In the absence of SOC, nodal loops are formed when the minority and majority bands intersect (see arrows). The majority and minority spin bands are hybridized when SOC is included, leading to a modulation of the equilibrium k-space spin-polarization density and an enhanced Berry curvature near the avoided anticrossings [55]. We now discuss the quasiparticle self-energy at the Γ point, Σ η (ω) ≡ Σ η (k = Γ, ω) -other k points behave similarly and hence are not further discussed here. The energy resolution in the Chebyshev polynomial expansion [Eq. (9)] is set to η = 1 meV. As it turns out, excellent spectral convergence is achieved after N 16384 iterations. Our results summarized in Figs. 3(c)-(d) disclose a rich self-energy structure. For conciseness, only the dominant matrix elements are shown (there are 18 nonzero matrix elements in total; see Appendix V D 2). These results counter conventional wisdom, which posits that the self-energy has essentially a scalar structure Σ η→0 + = −i Γ. The scalar component of the self-energy is dominant in relatively simple systems, such as graphene with point defects as discussed earlier (Sec. II A), which, appropriately far from the Gade singularity, displays an essentially scalar structure. However, the current example clearly illustrates that all symmetry-allowed matrix elements of the self-energy are generally nonzero, provided the existence of one-body interactions coupling different degrees of freedom inĤ 0 . In fact, the disorder selfenergy shares its matrix structure with the clean Hamiltonian. In particular, the SOC term inĤ 0 is responsible for the nonzero spin-flip components of the self-energy discussed below. The types of impurity potential and disorder statistics also play a crucial role in determining the self-energy matrix structure. Because we have considered a uniform on-site potential across all 4d-orbitals on a given impurity site [Eq. (14)], the impurity scattering effectively acts as a source of correlation between every orbital. The disorder correlation is thus responsible for the emergence of off-diagonal matrix elements at the lowest Born order in the diagrammatic expansion (such as Σ yz↑,xz↑ ), which would otherwise be forbidden. For further details, we direct the reader to Appendix V E.
It is interesting to contrast the results for random onsite disorder against a standard diagrammatic calculation performed at the bare Born approximation (BA) and self-consistent Born approximation (SCBA) levels. We find that for weak disorder (W 50 meV), the numerically exact, BA and SCBA self-energy are all in excellent agreement. For intermediate disorder strengths W ≈ 0.5 eV, the simple BA is no longer able to provide a satisfactory approximation to the self-energy. While the SCBA produces a better agreement, it is still unable to capture some of the finer details of the self-energy that our method captures (red arrows in Fig. 3c). This is to be expected, since, for strong disorder, quantum interference corrections are ubiquitous and cannot be captured by the SCBA. Interestingly, the combination of disorder correlation between the spins and a considerable (45 meV) SOC term induces a substantial spin-flip matrix element d xz,↑ d xy,↓ (see Appendix V D 2 for the other matrix elements).
For the second model of disorder, we use vacancies with concentration c and compare our results to the Tmatrix approximation [ Fig. 3(d)]. At low concentrations (c ∼ 0.1%), the T -matrix approximation is in complete agreement with our results, but starts to break down at higher concentrations (c ∼ 5%). To see this more clearly, we show the self-energy as a function of the concentration for fixed energies in Fig. 3(e). For low concentrations, the self-energy is proportional to the concentration of impurities in accord with the T -matrix result, but at higher concentrations, we start to see discrepancies which scale as ∼ c 1.4 near the peak, signaling the onset of nonperturbative disorder corrections. We note that such peaks cannot be attributed to van Hove singularities because they are absent at low defect concentration, and there is no correlation between the position of the peaks and the position of the singularities. We attribute them to resonances induced by multi-vacancy clusters, which only start to form at higher defect concentrations.
The renormalization of quasiparticles with a non-scalar self-energy, as predicted here, is expected to strongly impact the response of materials to external perturbations. For example, recent theoretical studies have alluded to robust spin-orbit scattering mechanisms underlying the extrinsic generation of spin Hall currents and currentinduced spin polarization, which can be traced back to the matrix structure of the disorder self-energy in spinorbit coupled materials [66][67][68]. On a fundamental level, the self-energy is connected to the four-point vertex functions of linear response theory through exact symmetry relations known as Ward identities [69][70][71] and thus the knowledge of all its matrix elements is essential to obtain physically sensible transport equations. Because our approach provides a systematic way to accurately evaluate the disorder self-energy of arbitrarily complex model Hamiltonians, regardless of the type and strength of disorder, it could provide new insights into the array of rich interfacial magnetic phenomena beyond the reach of diagrammatic calculations [72].
C. Disorder-enhanced p-wave superconductivity While our discussion so far has focused on lattice models with conventional quasiparticles, it is straightforward to generalize our approach to other condensed phases. As a final application, we employ our computational machinery to map out the mean-field phase diagram of a dirty superconductor. For definiteness, we focus on monolayer graphene, whose leading doping-dependent superconducting instabilities include chiral p-wave pairing states [73,74]. Chiral superconductivity has caused great excitement because it provides a platform to realize Majorana zero modes that are insensitive to local perturbations, and thus can be used to construct topological qubits [75]. Typically, disorder is detrimental for unconventional (non s-wave) superconductivity due to the breakdown of Anderson's theorem when the impurities violate the pairing symmetry thus acting as pair breakers [76][77][78][79]. However, there are exceptions to this rule (e.g., in d-wave cuprates, disorder is known to enhance the critical temperature through the appearance of superconducting islands around the impurities [80]). In this respect, the unusual electronic properties of graphene open up interesting possibilities. For example, it is known that the addition of scalar impurities in charge-neutral graphene has the counterintuitive effect of enabling conventional superconductivity for weak attractive interactions, while leading to a suppression of superconductivity in the strong attraction regime [81,82]. On the other hand, the phase diagram at finite doping, where superconductivity is expected to develop more easily due to a nonzero single-particle density of states, is far less explored. How is the doped graphene's ability to form superconducting states affected by impurity scattering? Can one tune the competition between different pairing states by tailoring the impurity potential (e.g., using adatoms adsorbed on particular lattice positions)? Here, we make a start on addressing these questions by computing the superconducting gaps at finite charge carrier density in the presence of disorder and competing order parameters. In the following, we demonstrate that p-wave superconductivity can be substantially enhanced by the presence of bond disorder (e.g. generated by gauge fields due to strain or adatoms at bridge sites).
Chebyshev-Bogoliubov-de Gennes formalism-The algorithm we use is a variant of the Chebyshev-Bogoliubovde Gennes algorithm proposed in Ref. [83], adapted here to reconstruct real-space Green's functions by means of the exact spectral decomposition [Eq. (9)] throughout the mean-field self-consistency cycle; see Fig. 1 (second panel). To speed up the evaluation of Chebyshev moments, we use a domain-decomposition technique as implemented in the open-source code KITE [32]. The proposed approach has two key features: (i) it is sufficiently powerful to handle systems with millions of sites, thus bypassing finite-size effects that severely restricted the accessible coupling strengths in previous studies; and (ii) it provides flexibility to treat disorder to different levels of accuracy, ranging from a self-consistent effectivemedium-type approximation all the way to a numerically exact treatment. This offers the possibility to capture important features of dirty superconductors missed by the commonly employed self-consistent T -matrix approximation, including multiple-impurity quantum interference phenomena and inhomogeneous pairing patterns.
The mean-field Hamiltonian reads asĤ =Ĥ 0 +Ĥ P , whereĤ 0 is the single-particle Hamiltonian of disordered graphene [Eq. (13)] and (15) is the pairing term. In the above expression, a † is (b † is ) adds an electron with spin s =↑, ↓ to site i on sublattice A (B) and g 0 (g 1 ) is the onsite (nearest-neighbour) interaction energy. The superconducting order parameters are ∆ 0,i = a i↓ a i↑ = b i↓ b i↑ and ∆ 1,ij = a i↓ b j↑ − a i↑ b j↓ and E 0 = −g 0 ∆ 2 0 − 3g 1 ∆ 2 1 is the condensation energy. The clean system achieves quantum criticality at zero chemical potential (µ = 0) and displays different superconducting phases depending on the interaction strength. There are regions of s-wave, p-wave and mixed symmetry. Away from half filling, every region is of mixed symmetry and the quasiparticle spectrum is gapped. In this regime, the order parameters are smooth across the whole diagram [73].
The order parameters are expressed as where f is the Fermi function, G is the retarded Green's function operator and i and j denote nearest neighbours. The exact choice of cutoff Λ depends on the origin of the pairing term. For conventional phonon-mediated superconductors, this is the Debye energy ω D , which restricts the integration to a thin shell around the Fermi order parameter a) b) c) d) level. On the other hand, unconventional superconductors, such as plasmon-mediated metal-coated superconducting graphene [73], generally have contributions from several energy regions. For the sake of simplicity, we let Λ → +∞, which captures the whole spectrum and is still able to accurately reproduce the clean phase diagram. In the simulations to be presented below, a set of random sites belonging to either sublattice are selected as the impurity sites, with a uniform concentration c = 5%. Around each impurity site, the nearest-neighbor hoppings are weakened by an amount ∆t.
Spatial dependence of order parameters-For a clean system, Eqs. (16)- (17) are easily diagonalizable, yielding a set of two self-consistent equations. When disorder is introduced, there will be four coupled self-consistent equations for each lattice site, which severely limits the system sizes accessible to exact diagonalization. Traditionally, dirty superconductors have been addressed by means of T -matrix and coherent potential approximation (CPA) schemes [84,85]. Here instead, we implement a different strategy that will allow us to keep the complexity to a minimum. Specifically, we restrict the order parameters to a constant uniform value in the bulk of the superconductor, while allowing them to vary in a circle of radius z around each impurity [ Fig. 4 (a)]. Since we expect the behavior of the order parameters to be similar in the vicinity of each impurity, we further restrict each order parameter to follow an identical spatial profile ∆ n = f n (r) (n = 0, 1) inside every circle and then find the function f n (r) which satisfies the self-consistent equations within these restrictions. This function is then computed using a stochastic evaluation of the matrix elements by means of the spectral approach described earlier (further details are given in Appendix V F). This approach is best suited for dilute point defects, where the regions rarely overlap. The choice of radius z regulates the approximation. A larger z allows us to capture more of the spatial dependency of the order parameters, at the cost of increased running time. The positions of the impurities are chosen randomly such that no two circles overlap. The spatial modulation of ∆ is expected to decay very rapidly [86], so we choose z to be around three lattice spacings. Figure 4(b) discloses a rich array of behaviors for different regions in the g 0 , g 1 plane. Most of this diagram for doped graphene can be understood within a virtual crystal approximation, as if the only effect of disorder is to renormalize the energy scales of the problem. The effective hopping is given by t = t + (c/z)∆t, with z = 3 the coordination number of the honeycomb lattice. Since ∆t is negative, the effective value of the interaction constants would increase to g 0 = g 0 t/ (t + c∆t/z) (likewise for g 1 ). Such a simple description is able to satisfactorily explain the behavior seen in regions of increased pwave and s-wave order parameters (see Appendix V F 3). However, there are two prominent features that cannot be captured by this heuristic argument. First, the threshold for superconductivity is substantially reduced when disorder is introduced, a behaviour typical of the presence of superconducting islands. Secondly, there is a line in the phase diagram across which the order parameters suffer an abrupt change that signals the onset of a crossover driven by disorder [white dashed line in Fig. 4(b)]. To better understand this transition, we plot in Fig. 4(c) the order parameters for fixed g 1 = −1.7t as a function of g 0 [horizontal gray line in Fig. 4(b)] for several values of the concentration. This is a regime of mixed symmetry, since both order parameters are nonzero. The discontinuity persists even at low (1%) concentrations of impurities, but the value of g 0 for which it happens increases with increasing concentration. We also plot two matrix elements of the self-energy as a function of energy at this value of g 1 for four different values of g 0 near the transition (Fig. 4d). The discontinuity in the order parameters is also reflected in these matrix elements. On both sides of the transition, we see the presence of a gap due to the finite value of both order parameters, but only the curves at the right of the transition have a van-Hovelike singularity. The existence of all off-diagonal matrix elements indicates that this disorder correlates different sublattices and spins, through the induced spatial inhomogeneity of the order parameters around the impurities. The self-energy can therefore provide valuable information about the local density of states around impurities through its connection to the LDOS [87]. Since these superconducting phases may be mediated by plasmons in a proximitized metal layer, the corresponding coupling parameters may also be controlled by changing the plasmonic properties of the metal or the distance between the metal and the graphene sheet. This opens up the possibility of sweeping the coupling parameters across the discontinuity line to look for the crossover, which may be identified experimentally through spectroscopic studies around impurities.

IV. CONCLUSION
We introduced a real-space numerical framework that gives access to the full wavevector and frequency dependence of the quasiparticle self-energy in arbitrarily complex disordered tight-binding models. For this purpose, we employed an exact Chebyshev decomposition of lattice Green's functions (Eq. 9) that provides full control over the resolution of the computation. The method was applied to 3 distinct problems (i.e., quantum criticality driven by chiral disorder in the honeycomb lattice, impurity resonances in a spin-orbit coupled ferromagnet and disorder-enhanced superconductivity in monolayer graphene), in unprecedented large systems, revealing a rich array of nonperturbative effects that challenge the standard perturbative picture of disordered systems.
For honeycomb lattices with a dilute concentration of vacancy defects, we uncovered nonzero off-diagonal components in the self-energy as well as a strong momentum dependency near the Gade singularity at zero energy. These previously inaccessible effects are ubiquitous even in the weak disorder limit, which shows that the quantum interference of electronic waves scattered by multiple defects plays a much deeper role than pre-viously thought. A striking result emerges in the long wavelength limit: the imaginary disorder self-energy at the Dirac (K) points vanishes ImΣ dis (K, 0) 0, within spectral resolution and accuracy, for any vacancy concentration. This suggests that the puzzling universal metallic conductivity σ Gade = (4/π)e 2 /h previously seen in large-scale quantum transport simulations of defected graphene systems [45] is a fundamental property of chiralsymmetry-protected zero energy modes with exceedingly large quasiparticle lifetime.
Secondly, our real-space self-energy framework is applied to 2D spin-polarized electron gases formed in SrRuO 3 /SrTiO 3 superlattices and compared against diagrammatic resummation schemes, including the T -matrix for vacancies and the self-consistent Born approximation for uncorrelated on-site disorder. Both are found to be in excellent agreement with the numerically exact self-energy at small disorder strength/impurity concentration, though they diverge significantly away from this regime. The most distinct effect of a large concentration of vacancies is the appearance of peaks in the self-energy whose height has a nonperturbative dependence on the concentration, which we attribute to scattering resonances from impurity complexes. A rich matrix structure of the self-energy is borne out by our study, illustrating how correlations in the disorder potential can manifest as off-diagonal matrix elements in the self-energy.
Finally, we applied the Chebyshev polynomial Green's function machinery to calculate the order parameters of superconducting monolayer graphene in the presence of dilute bond disorder. We found that for some regions of the phase diagram it is possible to enhance bulk s-wave, p-wave or both kinds of superconductivity by adjusting the amplitude of local fluctuations in the hopping parameters. For this purpose, a variation of the Chebyshev-Bogoliubov-de Gennes method was used to enable selfconsistent simulations of systems with millions of atomic orbitals. These results open up the intriguing possibility of tailoring superconducting phases in twisted bilayer graphene, a topic that is currently of much interest.
We briefly comment on possible extensions of the realspace spectral framework for the disorder self-energy that we introduced in this work. For conciseness, we restricted ourselves to orthogonal local basis sets, but this requirement can be easily relaxed at the cost of introducing an overlap matrix in the calculation of Chebyshev moments [viz. Eq. (2) and discussion therein]. The use of a nonorthogonal representation of the orbitals would open doors to accurate studies of disorder effects in complex problems and materials [88,89]. The method can also be easily extended to evaluate the self-energy resulting from other types of disorder such as local structures of defects, random hoppings and correlated disorder. Another interesting question for future study is whether the framework introduced here could shed new light on the behavior of mesoscopic systems without self-averaging properties.
Code availability statement.-The codes and scripts used in this study are available from the authors upon reasonable request.

A. Self averaging property
The self-averaging behavior of the disorder self-energy yields highly converged results in a computationally efficient manner and so must be justified. While a rigorous general proof is beyond the scope of this paper, we show below that under some rather general assumptions, the matrix elements G αβ (k, ω) = α, k|Ĝ(ω)|β, k (and hence the quasiparticle self-energy) satisfy the selfaveraging lemma where ... indicates disorder (configurational) averaging and D is the Hilbert space dimension of the lattice model that scales with system volume (a similar expression holds for the real part of the matrix elements). To simplify the discussion, we specialize to single-orbital models and thus omit the orbital index α, β hereafter. Let ξ k = k|Ĝ |k denote the matrix element that will be used to determine the self-energy. It is im- Specifically, we want to show that the imaginary parts of ξ k display self-averaging behavior, that is var ξ k ≡ (Im ξ k ) 2 − Im ξ k 2 ∝ D −1 . The argument is identical for the real part of ξ k . We consider two common classes of problems for lattice models defined with an arbitrary number of spatial dimensions: (i) systems characterized by perturbative (weak) disorder effects; and (ii) systems possessing exponentially localized single-particle states in their spectrum. Finally, we provide numerical evidence of our claim.

Weak disorder
If the diagrammatic expansion of the Green's function is convergent, then we can use an expansion in powers of V ≡V dis , the disorder potential: Im k|Ĝ 0 VĜ 0 n |k to evaluate the disorder average Im ξ k Im ξ k , keeping in mind that the term Im ξ k Im ξ k will remove all the terms in the diagrammatic expansion which do not connect both Green's functions. Defining G k 0 = k|Ĝ 0 |k for convenience and using R| k = D −1/2 e iR·k to express the disorder potential in real space, we obtain ... Figure 5. Feynman diagrams that contribute to the variance of Im G k up to fourth order inV dis .
We get a factor of 1/D from every disorder insertion V R and also a factor of D due to the sum over R. We assume that V R is an uncorrelated disorder potential with Gaussian statistics, i.e.
As explained below this assumption is not essential, but it aids in substantially simplifying the analysis. The configurational average introduces correlations between the disorder insertions as Kronecker deltas δ RR between different positions. Each δ RR effectively contributes with an additional factor of 1/D. Lastly, each loop in the diagrams (representing integrations over internal momenta) contributes with another factor of D. Figure 5 shows the diagrams that contribute to the variance up to fourth order inV dis . Counting all the powers of D, one can check that each term is associated with a factor of 1/D except for diagram (b). Instead, this diagram is proportional to (D − 1)/D, but the constant term gets cancelled precisely by Imξ k Imξ k and what is left is again proportional to 1/D. At higher orders in V dis , similar arguments can be made. If the upper branch of the diagrams is not connected to the lower branch, then it will get almost completely cancelled by Imξ k 2 , leaving only the 1/D contribution. If both branches are connected, the number of loops is not large enough to destroy the 1/D dependency.
While we have only strictly presented our argument for uncorrelated disorder, we argue that a generalization to correlated disorder should also be possible provided that the correlation length is finite. In such a scenario, averaging over disorder would introduce asymptotically decreasing functions of the distance between R and R in lieu of Kronecker deltas. In any case, a sum over the position (which would contribute with a factor of D as noted above) now contributes with a factor of order unity, effectively having the same effect as the Kronecker delta for the purposes of self averaging.

Localized states
Next, we analyze an important class of problems where diagrammatic methods break down [90]: strongly disordered systems with localized states in their spectrum. We assume that the value of ω is such that all the states in an energy window η around ω are localized, with a maximum localization length of ζ. We begin by expressing ξ k in terms of the eigenstates with energies {ε α } resolved in space where g R,k represents the contribution to ξ k from the sites around R. By assumption, these states are localized, so, for each R, only localized states with localization center within a distance 2ζ around R contribute. Let S be this region. This means that g R,k and g R ,k have appreciable correlation only if |R − R | < 2ζ. It is important to note that g R,k is independent of the system size, since the percentage of localized states is assumed to be an intensive property. This is a key assumption of this proof and fundamentally relies on the existence of a mobility edge. Note that g R,k is a random variable with a finite maximum absolute value because only a finite number N e of elements contribute to both the sum over R and the sum over α. Using the triangle inequality, For both sums, N e is the number of degrees of freedom inside a d-dimensional sphere of radius 2ζ. Therefore, the sum D −1 R g R,k can be seen as a sum of bounded random variables which are only correlated within a distance |r| < 2ζ of one another. Thus ξ k follows the central limit theorem and so var ξ k ∼ D −1/2 , hence proving the self-averaging property. We note the only assumptions made in this derivation were the locality of the localized wave functions and that only localized wave functions have relevant spectral weight in G k .

Numerical demonstration
Next, we demonstrate the self-averaging behavior numerically for the nontrivial example of a graphene system hosting a Gade singularity at the band center (ω = 0) generated by dilute point defects. To this end, we calculated the self-energy at the k = K point of very large lattices (D = 2L 2 = 10616832, 42467328, 169869312 and 679477248) for 100 realizations of disorder. This allowed us to obtain the standard deviation of this stochastic quantity (colored curves in Fig. 6), as a function of the linear system size. If self-averaging is taking place, then the standard deviation should decrease as ∼ L −1 , which is exactly what we obtain. The black dashed curves in Fig. 6 indeed have slope of −1 in a log-log scale, which indicates the correct scaling law.

B. Multi-scale domain decomposition and other computational details
Computations for this work were carried out with the open-source KITE code [32]. Periodic boundary conditions were employed in all calculations. KITE implements an efficient decomposition of the exact Green's function in terms of Chebyshev polynomials, which was used both for the calculation of the self-energy in all case studies as well as the local Green's function in the superconductor problem in Sec. III C. The real-space Green's functions are evaluated using the CPGF approach [Eq. (9)] and the calculation of α, k|G(ω)|β, k (or α, R|G(ω)|β, R for the superconductor problem) relies entirely on evaluating the Chebyshev moments α, k|T n (Ĥ)|β, k . The recursive nature of our spectral approach means that the complexity of the calculation of this matrix element scales linearly with the number of polynomials M .
Every object in the spectral approach is expressed in real space, exploiting the sparseness of the Hamiltonian matrix to improve the parallelization performance during the computation of the Chebyshev moments. In addition to being sparse, this HamiltonianĤ typically only connects sites that are close by neighbors. In the process of the matrix-vector multiplication |v =Ĥ|v , KITE divides the vector |v into equally-sized domains, which get assigned to different processing cores, and further subdivides these domains into tiles. Within each core, this multiplication is completed within each tile first, before moving on to the next tile to minimize cache misses when bringing the tile from memory. The size of the tile is adjusted according to the processor's cache size to maximize performance. This is at the core of the parallelization scheme in KITE [32]. Each processor performs the realspace matrix multiplication within its assigned domain, but it requires information about the other neighboring domains to correctly perform the multiplication around the borders. This is mitigated by keeping a copy of the neighboring domains' borders stored in memory for each processor. After the matrix-vector product has been calculated, the stored copy of the borders is updated before proceeding to the next iteration. This step is not parallelizable, but it scales with the area of the lattice rather than the volume, so the algorithm becomes more efficient with increasing lattice size. In Sec. II A, the k points chosen for the calculation of the self-energy belong to the (finite discrete) first Brillouin zone. Since the K point of the honeycomb lattice may be expressed as (b 1 − b 2 ) /3, in terms of the reciprocal lattice primitive vectors, the linear system sizes used were limited to multiples of 3. Failing to do so may induce noticeable errors due to the overlap with other momenta coming from the decomposition of k in terms of vectors belonging to the first Brillouin zone.

C. Spectral convergence
Convergence against several factors has been carefully assessed in all problems studied in this work, specifically: 1. The energy resolution η has to be as small as possible to accurately capture the singular nature of the Green's functions; Here we address points 1-4 with the help of Fig. 7, where we show the disorder self-energy of graphene calculated for a selected concentration of vacancies (0.3%), two different system sizes and several broadening factors η.
1. Due to the singular nature of this problem and the very fine resolutions required, full convergence with η is challenging to achieve, particularly at the band center. Sufficiently far from ω = 0, the orange curve in Fig. 7 appears to be converged. The inset shows the self-energy dependence as a function of η for states with ω = 0. This curve is well fitted by Im Σ AA (K, 0) = η 2/3 (orange line) which at η = 0 extrapolates to Im Σ AA (K, 0) = 0.
2. The curves no longer change when we increase further the number of polynomials, indicating that the Chebyshev series used to calculate the Green's functions is fully converged. The maximum number of polynomials was 131072.
4. The statistical fluctuations due to the disorder realizations are very small, of the order of 2 meV when L x = L y = 18432a and therefore do not influence our results (see previous section).
D. SRO

Multi-orbital model
The Hamiltonian used in Sec. III B consists of a sixorbital tight-binding model on a square lattice, which can be divided into four terms: represents the nearest-neighbor interactions, with t 1,x = t 2,y = t 2 , t 2,x = t 3,x = t 1,y = t 3,y = t 1 . The operator d † iaσ creates an electron in site i orbital a (yz = 1, xz = 2, xy = 3) and spin projection σ. The second term is the second-nearest-neighbor contribution to the Hamiltonian with g 1 = g 2 = t 3 , g 3 = t 4 and f 12 ij = f 21 ij = f if i and j are along the diagonal and f 12 ij = f 21 ij = −f if they are along the antidiagonal. The termsĤ 3 andĤ 4 encode the Zeeman interaction and spin-orbit coupling (SOC), respectively, with the expressionŝ where m is the amplitude of the Zeeman interaction, λ is the amplitude of the SOC, τ i is a Pauli matrix (τ 1 = σ x , τ 2 = σ y and τ 3 = σ z ) and ε abc is the Levi-Civita symbol. The SOC term was calculated by evaluating the matrix elementsL ·Ŝ in the angular momentum basis with = 2 restricted to the Cartesian set (i.e. {xy, xz, yz}).

E. Disorder operator correlations
In this section we show that the elements appearing in the self-energy matrix depend on the type of correlations appearing in the disorder potential. We find the dependency explicitly within the first Born approximation. In this paper, we use two kinds of disorder: Anderson disorder and dilute (short-range) disorder. We perform this analysis for Anderson disorder, but the proof is similar for the dilute case. Consider the following disorder oper-atorV ≡V dis in the Hamiltonian. The disorder operator V is diagonal in real space and is given bŷ where W iα is taken from the uniform distribution [see e.g., Eq. (14)] with width W and mean 0. Here, i indexes the unit cell and α indexes the degrees of freedom within the unit cell such as spin, orbital and/or sublattices. In momentum space, Now we assume that W iα is uncorrelated among different unit cells, but there might be correlations within any given unit cell. For example, if the α indices refer to spin up and down, one might require that W i↑ = W i↓ , which would impose W jα W iγ = δ ij W 2 12 for arbitrary α and γ. In contrast, if one is dealing with graphene and the α indices refer to the sublattice, one might require that the value of W in A be independent of that in B, yielding W jα W iγ = δ ij δ αγ Here, only the 1-point irreducible diagrams are to be kept. In the Born approximation, we obtain and it is clear that the matrix elements appearing in Σ depend on the correlation matrix M . If the disorder is completely uncorrelated, only the diagonal elements will survive. Despite this, higher-order terms beyond the Born approximation may contribute to the off-diagonal matrix elements of the self-energy.

Computational details
Here we provide additional details on the superconducting order parameter calculation. The starting point is the Chebyshev-Bogoliubov-de Gennes formalism [83], where the order parameters are obtained from Eqs. 16 and 17. In the clean case, the order parameters share the periodicity of the crystal. In the presence of disorder, we expect the order parameters to change appreciably in the vicinity of impurities and defects, but to remain relatively constant otherwise. In some cases, this modulation around an impurity may extend up to a few dozens of nanometers [91], but in our case, with dilute nonmagnetic impurities in graphene, the order parameter only changes appreciably on the order of a few unit cells [86].
For this task, we propose a new approach to calculate the order parameters which takes into account most of the spatial modulation while only requiring a small number of self-consistent equations be solved. Taking ∆ 0,i as an example (an identical procedure is used for ∆ 1,ij ), we start by defining a circle centered around each impurity site R n with radius z. Outside of these circles, the order parameter ∆ 0,i is restricted to be uniform: ∆ 0,i = ∆ b 0 . Inside every circle, the order parameter is assumed to behave identically, so for the sites r i such that |R n − r i | < z, the order parameter satisfies ∆ 0 (R n + r i ) = ρ 0,i for every R n and for some (yet to be determined) function ρ 0,i . In each self-consistent step, the value for the next ρ 0,i is defined as an average over identical sites within each circle ρ 0,i = N −1 imp Nimp n=1 ∆ 0 (R n + r i ), where N imp is the number of impurities. ∆ b 0 is calculated by averaging over the remaining sites. With this method, the number of self-consistent equations that have to be solved scales with the area of one circle ∼ z 2 instead of the number of lattice sites. We are able to greatly reduce the number of self-consistent equations that have to be solved while still capturing most of the spatial dependency of ∆ 0,i . The radius of the circles can be adjusted in order to better reflect the spatial dependency of the order parameters around the impurities.
Coming back to Eq. 16, the new self-consistent equation that has to be solved is R n + r i , ↑ |G (ω) |R n + r i , ↓ .
The notation has been slightly changed to reflect the vector addition of position vectors. While before only the site index was specified, now the position R n +r i is specified. This matrix element can be expressed in terms of ξ Rn |R n + r i , σ , thus casting the expression for ρ 0,i into where the bar denotes a random vector average. This expression can now be calculated with direct Chebyshev expansion of the Green operator in an efficient manner.

Order parameters
Using the mean-field model of the main text, we first calculated the bulk superconducting order parameters ∆ 0 and ∆ 1 in the clean system. In this case, the parameters are homogeneous and the sum in Eq. 18 runs over the whole lattice; see Fig. 9 (b) for ∆ 1 . The calculation was done at finite chemical potential µ = 0.4t. Both order parameters vary continuously over the g 0 , g 1 plane. A qualitatively similar picture exists for ∆ 0 .
Then, we performed the same calculation, but with the disorder specified in the main text (see Fig. 9 a)). Now the order parameters suffer a clear discontinuity which was not present before. The difference between these two graphs is the result presented in the main text, in a 2D color scheme [ Fig. 4(c)].

Renormalization
The main effect of the impurities is to renormalize the energy scales of the problem. Assuming the energy scales vary as g 0 = g 0 [1 + αt/ (t + ∆t)] for some α depending on the concentration, we get the left panel of Fig. 9(d).