Fate of disorder-induced inhomogeneities in strongly correlated d-wave superconductors

We analyze the complex interplay of the strong correlations and impurities in a high temperature superconductor and show that both the nature and degree of the inhomogeneities at zero temperature in the local order parameters change drastically from what are obtained in a simple Hartree-Fock-Bogoliubov theory. While both the strong electronic repulsions and disorder contribute to the nanoscale inhomogeneity in the population of charge-carriers, we find them to compete with each other leading to a relatively smooth variation of the local density. Our self-consistent calculations modify the spatial fluctuations in the pairing amplitude by suppressing all the double-occupancy within a Gutzwiller formalism and prohibit the formation of distinct superconducting-`islands'. In contrast, presence of such `islands' controls the outcome if strong correlations are neglected. The reorganization of the spatial structures in the Gutzwiller method makes these superconductors surprisingly insensitive to the impurities. This is illustrated by a very weak decay of superfluid stiffness, off-diagonal long range order and local density of states up to a large disorder strength. Exploring the origin of such a robustness we conclude that the underlying one-particle normal states reshape in a rich manner, such that the superconductor formed by pairing these states experiences a weaker but spatially correlated effective disorder. Such a route to superconductivity is evocative of Anderson's theorem. Our results capture the key experimental trends in the cuprates.


I. INTRODUCTION
The study of disordered high-temperature superconductors (HTSC) 1,2 is scientifically rewarding for multiple reasons. First, we learn to deal with strongly correlated electrons. Significance of strong repulsive correlations is paramount in these superconductors which make their parent (undoped) compound an antiferromagnetic Mott insulator. 3 Second, the superconductivity in the HTSC cuprates is believed to originate from the two-dimensional copper-oxide (CuO 2 ) planes, 4 and is rife with complex effects of enhanced quantum fluctuations from the reduced dimensionality. 5,6 Finally, these systems provide a hotbed for the complex interplay of electronic interactions and disorder, 7 both being strong. By tuning various parameters describing these systems including disorder, the cuprates can be made to undergo quantum phase transitions 8 to various non-superconducting phases. Each of such quantum phase transitions can vastly add to our knowledge of the condensed phases of matter.
While the physics of HTSC is far from being settled for the unusual normal state, the superconductivity in the clean systems is well described by a BCS-like ground state (GS) 9 at low temperatures (T ). Such a BCS-GS is identified with d-wave pairing amplitude, [10][11][12] and hence they are referred to as d-wave superconductors (dSC). The corresponding Bogoliubov quasiparticles result into a linear low-lying density of states (DOS). 13 The effect of impurities on this BCS-GS was developed within the conventional Abrikosov-Gorkov theory, 14 which infers that the non-magnetic impurities easily destroy the d-wave superconductivity. 15 Contrary to this wisdom, the cuprates are found significantly robust to disorder in the following senses: (1) The local density of states (LDOS) maintains the same low-lying V-shaped gap as in a pure dSC even in the presence of the disorder. [16][17][18] (2) Even though the gap-map develops nanoscale inhomogeneities, 19,20 the low-energy LDOS remains surprisingly homogeneous. 21 (3) The superfluid density undergoes only a modest reduction with the impurities. 22 (4) The superconducting gap in the nodal region gets suppressed by disorder while an enhancement of the antinodal gap is found 23 -a signature contributing to the nodal-antinodal dichotomy. 24,25 (5) Scanning tunnelling spectroscopy (STS) yielded an intriguing energyresolved nanoscale density modulation persisting on top of the inhomogeneous background at low T , as signalled by Fourier transformed local density of states (FT-LDOS). 21,26,27 While much of these surprises have been attributed to the strong electronic repulsions in these materials, a systematic calculation including both the correlations and disorder, 28 had always been challenging. The first step towards integrating strong correlations and disorder in a numerical calculation 29 has already uncovered interesting effects.
We incorporate both the strong interactions and disorder on the same footing for a two-dimensional dSC, and our main results are summarized as follows: (1) The nature and degree of the inhomogeneities in the order parameters change substantially from the scenario that neglects strong correlations. In particular, superconducting-'islands' cannot be discerned even at large disorder, though the nanoscale inhomogeneities develop. (2) Superfluid density, off-diagonal long range order and density of states show amazing insensitivity to the impurities. (3) The interplay of the strong correlations and bare disorder produces an effective background that is correlated in space. (4) Our results offer a simple description of the complex problem of disordered cuprates motivating a pairing mechanism similar to Anderson's theorem, 30 which traditionally describes only the s-wave superconductors (sSC).
Significant progress had been made in the past addressing separately the effects of disorder, 2 or strong correlations 31 on a dSC. The robustness of dSC to impurities had been partially addressed 32,33 within the Bogoliubov de-Gennes (BdG) mean-field theory that includes the inhomogeneities in the local order parameters, however, ignores any phase-space restrictions 34 from the strong correlations. We broadly refer to such calculations as inhomogeneous mean-field theory (IMT). The IMT produces an inhomogeneous pairing, 19,20 a gradual fall of the superfluid density with disorder, 32,35 a weaker filling of the low-lying DOS 32,33 and it helps visualizing the proposal for a 'swiss-cheese' 22,36 model. The IMT calculations had been extended to study the responses of a dSC to antiferromagnetism, 37 competing orders 38,39 and vortex-lattice. 40,41 However, the significance of the IMT results for HTSC is not obvious because of its mean-field treatment of the effective interactions, leaving the connection between the bare strong correlations and superconducting pairing artificial, and possibly dubious. The prime role of strong repulsive correlations in clean (pure) cuprates lies in prohibiting the double-occupancy of the charge carriers. A comprehensive theory for such a constrained GS is based on Gutzwiller's projected wave functions. [42][43][44][45] However, a simple framework, called Gutzwiller approximation (GA), 46 that restricts all double-occupancy by renormalizing the bare parameters of the Hamiltonian, had also been developed. The renormalization is easily rationalized, e.g., the electronic hopping is difficult due to the prohibition of double-occupancy, while the effects of exchange interactions are enhanced due to a large number of singly occupied sites. The resulting model is solved within the framework of a renormalized meanfield theory, 31,[47][48][49] which describes clean strongly correlated cuprates reasonably well. 50,51 Several attempts have been made recently to extend the GA analysis to the inhomogeneous superconductors. [52][53][54][55][56] We use the notation GIMT, for referring to such IMT calculations augmented with inhomogeneous GA. Incorporating the strong correlations, it has been shown 29,57,58 that the robustness of superconductivity goes far beyond what is obtained from a simple IMT. In particular, the site-averaged low-energy DOS remains nearly unaffected up to a rather large concentration (∼ 25%) of impurities 29 ! Such striking observation raises several questions: (1) How does the strong correlation make the impurities essentially disappear from the low-energy physics? (2) What happens to the disorder-induced inhomogeneities, whose presence is instrumental 32,59 for the inferences made by the simple IMT calculations? And finally, (3) How do the other observables, e.g., superfluid density or off-diagonal long range order (ODLRO) respond to the disorder? Answers to these questions are necessary for a better understanding of dirty superconductivity and we address them in this paper. Importance of inhomogeneous GA are currently being probed in a variety of contexts, e.g., on competing charge and spin stripes, 60 spin-density wave, 57 vortex-core physics, 61 antiphase superconducting domain, and on the moment formation around a single impurity. 62 The rest of the paper is organized as follows. In section II, we describe the model and parameters used for our study and also review the GA formalism for inhomogeneous systems. We discuss our results, together with a quantitative comparison between the IMT and GIMT schemes in section III, keeping our focus on the behavior of the physical observables. We further study the spatial structures in all the order parameters. We establish a connection between the nature of the inhomogeneities and the behavior of observables. In addition, we illustrate the importance of the underlying normal states in comprehending our findings. Finally, we summarize and conclude in section IV.

A. Strong correlations and Gutzwiller approximation for an inhomogeneous d-wave superconductor
We describe the phases of the strongly correlated cuprates by the Hubbard model, which is a minimal lattice model representing repulsively interacting band electrons, Here, the first term indicates electrons' hopping on a twodimensional square lattice, with c † iσ (c iσ ) as the creation (annihilation) operator of the electron on the site i and with spin σ. We take t ij = −t, when i and j are nearest neighbors, denoted as ij , and t ij = t ′ , when i and j are next-nearest neighbors, with the notation of ij . We choose t ij = 0 for all other pairs of i and j. U is the onsite interaction energy between electrons occupying the same site and mimics a Coulomb-type repulsion where inter-site contributions are neglected in favor of screening. The local density operator with spin σ at site j is denoted as n jσ = c † jσ c jσ . In the strongly correlated limit U ≫ t, the Schrieffer-Wolff transformation on H Hubb yields an effective "t − J" model: 63 where all terms up to O(t 2 /U ) are kept, and J ij = 4t 2 ij /U . Here, we used the notation J ij = J for ij , and J ij = J ′ for ij , andc iσ = c iσ (1 − n iσ ) is the annihilation operator in the 'projected space' prohibiting double-occupancy at the site i. The terms appearing in the second and third line of equation (2) represent the 'three-site terms', which along with J ′ -term are often neglected in the standard t − J model, even though they are of the same order as the J-term. In fact, these terms are capable of generating additional broken symmetry phases in the Cooper-channel, as we discuss in section II B. We used the quotes in "t − J" to emphasize that it has contributions beyond the standard t − J model. The notations ijm and ijm are for the 3 sites with (i, j) and (j, m) being nearest neighbors and next-nearest neighbors to each other respectively.
We introduce disorder by redefining H "t−J" to H "t−J" + iσ (V i − µ)n iσ , where V i is the (non-magnetic) impurity potential at site i and µ is the chemical potential that fixes the average density of electrons in the system to a desired value.
There are different ways of choosing V i 's, depending on the specifics of the experimental realizations. It is, however, important to quantify the degree of disorder by a small number of parameters defining the distribution of V i , for an appropriate comparison with the experiments.
For our calculations, we used three models of disorder: (i) Box-disorder : V i 's are drawn from a 'box' distribution, such that, V i ∈ [−V /2, V /2] uniformly, (ii) Concentration disorder (or conc-disorder) : V i = V 0 on a given n imp fraction of the sites, which are randomly chosen, and V i = 0 for all other sites. (iii) Out-of-plane disorder, as detailed in the supplementary material. While the box-disorder is quantified by a single parameter V , the conc-disorder needs two parameters for its characterization: V 0 and n imp . We will focus, in particular, on the box-disorder which is better suited to address several conceptual issues as we will see in section III C. All data are averaged over 10-15 independent realizations for a given strength of disorder.
Strong electronic repulsions restrict the Hilbert space of H "t−J" by eliminating all double-occupancies. One simple method to implement such restriction is the GA, in which the higher order off-site correlations are neglected. 46 Following Ref. 54, we consider the ground-state wave function in the projected space |ψ = P |ψ 0 , where |ψ 0 is the groundstate wave function in the Hilbert space that allows doubleoccupancy and P is the projection operator, P = i P i , with Here, γ i are the local fugacity factors obtained by demanding conservation of the local electron densities, ρ iσ = n iσ = n iσ 0 . In GA, the effects of projection on impure dSC are absorbed by local Gutzwiller renormalization factors. 53 The expectation value of a general operator A in different ground-state wave functions can be written as Â 0 = ψ 0 |Â|ψ 0 / ψ 0 |ψ 0 and Â = ψ|Â|ψ / ψ|ψ , and GA amounts to the following simple relations: Here g t ij , g 31 ijm , g 32 ijm and g xy ij are the Gutzwiller renormalization factors (GRF) which depend on an appropriate combinations of local doping x i , x j , x m (here x i = 1 − ρ i ). The GRFs for our system (without any magnetic order) are given by, We emphasize that since the GRFs are different for S i .S j (renormalized by g xy ij ) and n i n j (not renormalized) terms in H "t−J" , a flat renormalization of the exchange coupling, J, would not suffice.
Upon complete removal of the double-occupancy by GA, the total energy of the system can be calculated in terms of the relevant local order parameters: ρ i = σ n iσ ≡ σ n iσ 0 , ∆ ij ≡ c j↓ c i↑ 0 + c i↓ c j↑ 0 , and τ ij ≡ c † i↓ c j↓ 0 ≡ c † i↑ c j↑ 0 . We minimize ψ 0 |H "t−J" |ψ 0 , with respect to |ψ 0 with additional constraint ψ 0 |ψ 0 = 1 and i ρ i = ρ (ρ being the chosen total density), leading to the minimization of the functional W , where, employing the variational relation (δW )/(δ ψ 0 |) = 0. Here β and µ are the Lagrangian multipliers for the constraints mentioned above. Upon minimizing W , we obtain the following BdG mean-field Hamiltonian, 47,57,58,60 satisfying the eigenvalue equation: H MF |ψ 0 = β|ψ 0 . Upon expanding the derivatives in equation (8), we obtain the explicit form of the mean-field Hamiltonian, We denote the Cooper-channel order parameters ∆ ij on the bond connecting sites i and j as ∆ δ i where j = i + δ, with δ = ±x or ± y,δ = ±(x ± y), and the primed-summation means: The coefficients G 1 to G 5 in equation (9) are built with the GRFs defined in equation (6) and their exact expressions are included in the A. µ HS i and W FS iδ are the Hartree-and Fockshift terms arising from the mean-field decomposition of interactions and are also given in the A. We found that the major contributions of µ HS i come from the derivatives of the GRF for the parameters describing cuprates.

B. Iterative self-consistency: BdG method
We diagonalize H MF using the BdG transformations, 64 c iσ = n (γ n,σ u i,n − σγ † n,−σ v * i,n ), where γ † n,σ and γ n,σ are the creation and annihilation operators of the Bogoliubov quasiparticles. The resulting eigen-system is then solved selfconsistently for all the local order parameters. Inclusion of J ′ and the three-site terms in the Hamiltonian, allows additional Cooper-channel orders in H MF , namely; s xy and d xy , in addition to the standard d-wave (d x 2 −y 2 ) and the extended s-wave (s xs ) which take the following forms in the clean systems: where, ∆ sxy and ∆ (0) dxy are the uniform values of the pairing amplitudes with respective symmetries. We denote ∆ d x 2 −y 2 by ∆ d for the rest of the paper for notational simplification. Starting with initial guesses for all the local order parameters, we numerically find out the eigenvalues and eigenfunctions of the BdG Hamiltonian. With these, we recalculate the order parameters and compare them with the initial guesses. If the two do not match on all the sites, the whole process is iterated with a new choice of the order parameters, until self-consistency is achieved.

C. Model and parameters
We investigated the physics of H MF within the framework of GIMT for a wide range of parameters. We will describe our results for T = 0 in the next section for U = 12 and t ′ = 1/4 in the units of t for all our GIMT calculations, which are believed to describe a typical cuprate. 65 All other energies are also expressed in the units of t. Typical size of the unit cell for our calculations had been 30 × 30, and we focus on the results for the average density of electrons, ρ = 0.8 which coincides with near optimal doping in the cuprates. In general, HTSC materials show a low temperature symmetry breaking primarily in the Cooper-channel at the optimal doping and hence we do not consider other orders. [66][67][68][69][70] For the box-disorder we tune the disorder strength up to V = 3, which is sufficient to destroy dSC in the IMT calculations, as we will see in section III B. The chosen strength of the conc-disorder is V 0 = 1 for most of our calculations setting them to the 'Born scattering limit'. While we report some results even for V 0 ∼ 5 − 10 to address the effect of strong scatterers, we avoid the unitary limit, 57 e.g., in the case of substitutionary Zinc impurities. 71 They arguably produce a local antiferromagnetic surrounding, 57 and can lead to effects beyond the scope of our study. In order to carry out an appropriate comparison, a value of U = 3.69 is chosen for the IMT calculations which yields the same homogeneous (V = 0) d-wave gap from the GIMT scheme with the cuprate parameters. An alternative IMT scheme could also be set up and is discussed in the supplementary material. The H MF for plain IMT is recovered by setting the GRF in equation (6) to unity, however, suffers from a technical problem. The three-site terms deplete the d x 2 −y 2 order heavily leading to a cancellation of much of the d-wave pairing amplitudes. The GIMT calculations, however, remains free from such a deficit. Further, the s xy and d xy orders from the J ′ term were found insignificant in the GIMT calculations, as will be illustrated in figure 1(c). Thus, we present bulk of our results using just the t−t ′ −J model for the GIMT calculations with box-disorder and for IMT calculations with both the box-and the conc-disorder, dropping the three-site terms and J ′ term in favor of simplicity. Our GIMT results with conc-disorder allow all the Cooper-channel orders as it uses the full "t − J" model.

III. RESULTS
We begin this section showing our results quantifying the nature and degree of the inhomogeneities in the local order parameters in both the GIMT and plain IMT calculations.

A. Distribution of the order parameters
We plot along the y-axis of figures 1(a) and 1(b) the number of times a value of ρ occurs anywhere in the system against ρ itself along the x-axis from the GIMT and IMT findings. The resulting normalized distribution P (ρ) widens with V , however, the disorder dependence of P (ρ) is markedly different in the two calculations. For example, at V = 2.5, P GIMT (ρ) is finite only for ρ ∈ [0.6, 0.95], where as P IMT (ρ) is non-zero over a much wider range of ρ = 0.1 to ρ = 1.45. The Gutzwiller projection ensures ρ i ≤ 1. It is interesting to note that P GIMT (ρ) develops an asymmetry with increasing V with a larger weight towards larger values of density. Such asymmetry is absent in P IMT (ρ). The explanation for this feature, will be discussed in section III C. Results from the conc-disorder (which includes three-site terms and J ′ terms in GIMT) is broadly similar to the box-disorder with respect to the inhomogeneities, and are presented as the insets in the figures 1(a) and 1(b). As expected, P (ρ) develops finite weight for depleted values of ρ on the disorder sites (at ρ ∼ 0.7 for GIMT and ρ ∼ 0.4 for IMT calculations), in addition to the large contribution around its average, ρ ≈ 0.8. It is interesting that the electron population on the impurities is much larger in GIMT than in plain IMT. This ensures that the electronic correlations prohibit strong inhomogeneities. Further, it reduces the effective disorder strength, if we were to define this strength locally in terms of depleting population of the electrons in our context of repulsive scatterers. The reduction of the effective disorder has serious consequences on the physical properties of the superconductor and will be addressed in section III D.
The local density plays a crucial role in deciding the nature of the local d-wave order parameter. The double occupant sites (and empty sites) in the IMT calculation de- plete the d x 2 −y 2 -wave pairing amplitude locally, defined by . This is because these sites have extreme values of local density (namely, ρ i = 0 or 2) and do not participate in the particle-hole mixing, which is essential for a non-zero value to ∆ d (i). As a result, the dwave pairing amplitude from plain IMT reduces with increasing V (box-disorder) or with n imp (conc-disorder). The corresponding distribution, P IMT (∆ d ), as shown in figure 1(d), broadens with increasing V with large weight progressively shifting towards a smaller value. Such an evolution finally leads to a rather skewed In fact, this same physical picture 59,72 yielded a very similar disorder-dependence of P IMT (∆ s ) for a disordered sSC. On the other hand, the d-wave pairing amplitude in the GIMT calculation, which is free of doubly-occupied sites (and hence, free of near-empty sites too, satisfying the constraint of maintaining ρ = 0.8), stays around its homogeneous value (see figure 1(c)). Of course, there are inhomogeneities in the GIMT results, which give the increasing width in P GIMT (∆ d ) with V . However, the fluctuations in the pairing amplitude are weak (at least for moderate V ) and symmetric about ∆ (0) d yielding a P GIMT (∆ d ) nearly Gaussian in shape for all V . While P (ρ) indicates that the degree of inhomogeneity is different in the GIMT and IMT schemes, the results for P (∆ d ) also establish that the nature of inhomogeneity is substantially different in the two calculations. Such a disorder evolution of P GIMT (∆ d ) is strikingly different from the predictions of simple Abrikosov-Gorkov theory, or its extension in the form of self-consistent T-matrix approximation, 15,73 in which ∆ d (assumed homogeneous) depletes monotonically with V . It is interesting to note that the GIMT calculation yields ∆ d (i) ≥ ∆ (0) d on certain locations, particularly at large V , whereas the IMT value of ∆ d (i) is more or less limited by ∆ (0) d . The similarity of our P GIMT (∆ d ) and that found from the nanoscale inhomogeneity 19,20 in Bi 2 Sr 2 CaCu 2 O 8+x (BSCCO) is heartening, while the nature of P IMT (∆ d ) is unfound in cuprates. We draw the attention of the readers to the fact that in plotting P (∆ d ) in figure 1, we introduced an additional prefactor J(3g xy ij + 1)/8 such that the GIMT and IMT calculations produce same ∆ d at V = 0 (g xy ij = 1 for IMT). We also present the GIMT distribution of one sub-dominant Cooper-channel order for conc-disorder, ∆ dxy , as the right inset of figure 1(c). The tiny width of P (∆ dxy ) and also the smallness of max{∆ dxy (i)} at all V illustrate the insignificance of ∆ dxy compared to ∆ d . Their redundancy is further verified by repeating the GIMT calculations where ∆ dxy , ∆ sxy and ∆ sxs are suppressed by hand, without showing any change in the results. The disorder dependence of P (∆ sxy ) and P (∆ sxs ) are essentially same as that of P (∆ dxy ). Upon establishing their insignificance we drop all the sub-dominant orders for the box-disorder calculation in favor of simplicity and clarity.

B. Disorder dependences of observables
The qualitative as well as quantitative differences in the GIMT and plain IMT results for P (ρ) and P (∆ d ) raises the question: How do the inhomogeneities, or the lack of it in the presence of strong correlations, affect the disorder dependence of various physical observables? We address this issue in figures 2 and 3 where we show the evolution of ODLRO, superfluid density and average DOS with V . The results bring out a striking contrast in the V -dependence of these quantities with or without the strong correlations.

Off-diagonal long range order (ODLRO)
First, we discuss the superconducting ODLRO defined as: is the singlet Cooper-pair creation operator on the links connecting the neighboring sites at i and i + δ. F δ,δ ′ (i − j) is independent of δ and δ ′ in the limit of large |i − j|. Since F δ,δ ′ (i − j) can be interpreted as simultaneous hopping of a singlet cooper-pair on a link, the GRF corresponding to this process is g t i,j g t i+δ,j+δ ′ . The evolution of ∆ OP (normalized by its value ∆ (0) OP at V = 0) with box-disorder is presented in figure 2(a). With increasing V , ∆ OP from the IMT results decreases monotonically, initially with a gradual fall, followed by a rapid decrease beyond V ≥ 1, causing it to almost vanish by V ∼ 2.5. Such a behavior is consistent with the V -dependence of P IMT (∆ d ) in figure 1(d). The strong correlations make the superconductor robust against disorder. There is hardly any significant fall of the GIMT value of ∆ OP at V = 2.5, a strength by which the transition to a nonsuperconducting state occurs within the IMT scheme. The Gaussian-type shape of P GIMT (∆ d ) with minimally depleted centroid for all V is crucial to ascertain that ∆ OP sees little changes in the whole range of V within the GIMT calculation. The comparison of the GIMT and IMT results for ∆ OP from conc-disorder, as presented in the inset of figure 2(a), shows a trend similar to the box-disorder results. A numerically inexpensive scheme for approximating ∆ OP is further discussed in B.
The length scale of decay of F δ,δ ′ (i − j) with |i − j| defines the coherence length (ξ) of a superconductor. The evolution of F δ,δ ′ (i − j) for several V , as presented in the lower inset of figure 2(a) from the GIMT method, implies that there is hardly any V -dependence of ξ. Given that ξ is already small (only a few lattice spacings) in the clean cuprates, it is not expected to show any significant reduction with disorder. A similar insensitivity of ξ on V is also found in the IMT results.

Superfluid stiffness
The defining characteristic of a superconductor lies in the Meissner 74 effect, which is quantified by the stiffness of the BCS-GS wave function to an externally applied phase twist. This stiffness translates into its finite superfluid stiffness, which is proportional to the superfluid density. Within the Kubo formalism the superfluid stiffness, denoted as D s , is measured by, where k x is the kinetic energy along the x-direction and Λ xx is the long wavelength limit of transverse (static) current-current correlation function. 75 It is calculated by Fourier transforming the impurity averaged Matsubara Green's function; In order to obtain a good qresolution of Λ xx from a finite simulation box, we used an effectively larger system employing 'repeated zone scheme' (RZS) (see C) for a single realization of disorder. The final q y → 0 extrapolation of the disorder averaged Λ xx obtained from single unit cell was guided by Λ RZS xx (q y ). The disorder dependence of ∆ OP and D s are quite similar, while the IMT causes them to almost vanish by V ∼ 2.5, there is hardly any appreciable reduction of them in the GIMT results.
The inset of figure 2(b) shows the individual disorder dependence of −k x and Λ xx (q y → 0). The −k x from the GIMT method are smaller than the IMT values due to a complete removal of double-occupancy. The kinetic energy, which is the diamagnetic contribution to D s , decreases slowly with disorder and the rate of decrease is similar for both the IMT and GIMT calculations, because it is proportional to the total density. For V = 0, Λ xx is zero in the absence of any quantum fluctuations making D s equal to −k x . With increasing V , Λ xx grows in the IMT scheme and almost equals −k x by V ∼ 2.5 and continue to reduce their differences for larger V . However, Λ xx remains insensitive to V within the GIMT scheme making D s robust to impurities.
Such a strong insensitivity of D s to V is surprising though consistent with the broad experimental trends. 22 There are recent data on the temperature and doping dependences of D s 76 for the cuprates. We hope similar experiments will be extended to careful disorder dependence for a justified comparison with our results.

Density of states (DOS)
We now turn towards the discussion of the disorderdependence of the site-averaged DOS, given by: where E n 's are the BdG eigenvalues. We focus on N (ω) for ω up to the coherence peak energy. While N (ω) consists only of δ-function peaks at T = 0, it is necessary to obtain a smooth trace such that the relevant features are truly identified and is free of any artefact of smoothing. In order to achieve this, we generate a denser spectrum by using RZS (see C). We broaden each of the δ-function in N (ω) by an amount comparable to the average spacing of the E n 's, and the final N (ω) is presented in figures 3(a) and 3(b). After establishing that the additional Cooper-channel orders, e.g., ∆ sxs , ∆ sxy , ∆ dxy have insignificant contribution to N (ω), just as in all the other observables, we restrict the DOS calculations implemented with RZS only to t − t ′ − J model for all classes of disorder.
While the average N (ω) calculated within plain IMT shows significant gap filling in figure 3(b) for large V , results from the GIMT scheme in figure 3(a) show that the low-lying DOS remains unaffected with impurities even up to V = 3. Strong correlations thus appear to prohibit the formation of the lowenergy states in the GIMT results. While the coherence peak height gets reduced with V in both the methods, the rate of this reduction is significantly weaker in the GIMT. Our results for conc-disorder are similar to Ref. 29. The strong correlations thus seem to protect the V-shape of the N (|ω| ≤ ∆ (0) d ), a feature established in the STS-experiments. 21,27 In order to develop a deeper understanding of the robustness of low-energy DOS in GIMT, we calculated the LDOS for two regions with large and small ∆ d . We average the LDOS for these regions with a weak ∆ d (i) ∈ [0.18, 0.21] and with a strong ∆ d (i) ∈ [0.25, 0.28] in figure 3(c) and find that LDOS in the two regions behave quite differently even for N (|ω| ≤ ∆ (0) d ). We will also show in section III C that the regions of large (small) ∆ d (i) corresponds also to the regions of large (small) ρ i within the GIMT calculations, so that the inhomogeneity in the local density is also reflected in the same LDOS. While the anti-correlation between the coherence peak heights and the local gap values is clearly identified similar to experiments, 21 we did not find the expected homogeneity in LDOS for |ω| < ∆

Fourier transformed local density of states (FT-LDOS)
Another observable, that received attention in the recent times, is the energy resolved FT-LDOS, defined by N (q, ω) = i e −iq.ri N i (ω), where N i (ω) is the LDOS at site i, obtained from equation (18) without performing the summation over i. FT-LDOS is extracted from the STS data indicating the existence of a periodic and energy resolved modulation of the local density. An interpretation of the STS data comes from the simple 'octet'-model, 77,78 based on the scattering of the d-wave quasiparticles from a single impurity. The resulting theory yields an octet of wave-vectors q connecting the constant quasiparticle-energy contours for which the scattering probability is maximal. The Fermi surface obtained from the dispersion of these octet-members, agrees well with its direct measurement from the angle-resolved photo-emission spectroscopy. 79 The extension of the octet-analysis to manyimpurities is a subtle issue. Without invoking strong correlations, results from many-impurity calculations found sensitive dependence on the details of the disorder 80 for a reasonable fit to the STS data, and the octet peaks are often masked by a large noise. 81 Further, non-dispersive FT-LDOS peaks have been attributed to competing orders 38 for under-doped samples. However, to the best of our knowledge, the role of strong correlations on the FT-LDOS have never been addressed in a calculation.
We present the FT-LDOS from the GIMT calculations in figure S3. Each density-plot shows the power spectrum P (q, ω) = |N (q, ω)| 2 /N for a given ω on the first Bril-louin zone (BZ): q x , q y ∈ [−π, π], disorder averaged over 10 independent realizations of conc-disorder with n imp = 0.09. We apply RZS for a smooth data from a larger system (see C). Octet peaks can easily be recognized, such as q 1 ∼ (0, ±π/L 1 ) or (±π/L 1 , 0). The dispersion of this peak results into a continuous change in L 1 ∼ 1.36 for ω = 0.05 to L 1 ∼ 7.48 for ω = 0.2. Another (dispersing) peak, likely the q 7 in the octet terminology, 77 can also be discerned along the diagonal direction of the BZ (better resolved for negative ω). Our results recover much of the experimental trends, e.g., the shape and location of the FT-LDOS peaks, their dispersion (except may be q 5 ), among other things. However, the asymmetry of the FT-LDOS profiles for ±ω is found to be much stronger than the experimental data. The qualitative features of figure S3 remain unaltered with a change of n imp to 0.05 from 0.09. A similar calculation of FT-LDOS within the IMT scheme (shown in the supplementary material) was found to capture the key points of the figure S3, however, appropriate choice of the model-parameters for the IMT calculations is crucial for a justified comparison.
Before closing the discussions on the disorder dependence of the observables, we note that the degree of inhomogeneity in ∆ d reduces upon including strong electronic repulsions (See figure 1(c)), at least for the moderate disorder (0.5 ≤ V ≤ 2.0). However, the behavior of the observables remain very different compared to the self-consistent T-matrix predictions 15 based on the assumption of homogeneous but renormalized ∆ d for any V . The inadequacy of self-consistent T-matrix approximation for describing disordered cuprates is ω = 0.05 (π, π) (-π, -π) q 1 q 5 GIMT ω = 0.15 (π, π) (-π, -π) q 1 q 4 q 5 ω = 0.2 (π, π) (-π, -π) discussed in Ref. 33, which is a weak coupling theory built on simple metallic normal state. The strong interactions amplify such differences further due to modifications of the normal states, as we will see in section III E.

C. Spatial distribution of local orders and their relation to inhomogeneity
Having seen the dramatic robustness of the observables to impurities in the GIMT scheme, it appears natural to look for the relationship between such insensitivity and the inhomogeneity compared to plain IMT. We develop an important insight by studying the spatial distribution of the different order parameters on the lattice, as we discuss below.
We present the GIMT and IMT results for ∆ d (i) on the left and right columns of figure 5(a) respectively. Similarly, the local density, ρ i , from the GIMT and IMT calculations is shown in the left and right columns of figure 5(b). All results are shown for a particular realization of box-disorder. There are important points to note here: Firstly, the spatial inhomogeneity in density from the GIMT calculation does not evolve with V (not even in its magnitude) as seen from figures 5(b1), (b3) and (b5). While this observation reinforces the finding of figure 1(a), the additional information here is on the spatial structures, as explained below. The profiles of ρ i in plain IMT, shown in figures 5(b2), (b4) and (b6), are similar to that from the GIMT in figures 5(b1), (b3) and (b5) respectively, though its magnitude differs significantly (See also the figure 1(b)). Secondly, the evolution of spatial profile in ∆ d (i) from GIMT with V is also very weak as shown in figures 5(a1), (a3) and (a5), compared to the same calculated in plain IMT scheme as given in figure 5(a2), (a4) and (a6). Thirdly, an impressive spatial correlation is found between the GIMT profiles of ρ i and ∆ d (i) for all V implying that both are large (or small) approximately on the similar region of space (further illustrated in figure 7(b)), whereas such spatial correlations are completely absent in the IMT results. The IMT data in figure 5(a6) illustrates the formation of 'islands' with large ∆ d (i), separated by regions where ∆ d (i) ≈ 0. We verify that the sites with large ∆ d are indeed the sites where |V i − µ HS i | ≈ 0 for large V , and correspond to the spatial regions allowing maximum particle-hole mixing in plain IMT.
Clearly, the origin of the inhomogeneities leading to the island formation in ∆ d (i) within the IMT scheme is essentially the same as in the case of an sSC. 59,72 The lack of the evolution of the GIMT spatial structures in ∆ d (i), on the other hand, illustrates that there are no well defined 'islands'. Primarily, those sites lying in the deepest valleys of the disorder potential become nearly doubly occupied in plain IMT. As a result, no significant spatial correlation between ∆ d (i) and ρ i are expected at large V , which is clear from our data. Strong electronic correlations in the GIMT, however, ensure that the local density never goes beyond half-filling anywhere in the system resulting into the observation that the sites with large ρ i (close to half filling) also satisfy |V i − µ HS i | ≈ 0 approximately (at large V ) for ρ = 0.8. It is thus not surprising that we see a strong spatial correlation between local density and pairing amplitude in the GIMT results on the left columns of figures 5(a) and 5(b) (See also figures 7(b1) and (b2)), quite  in contrast with the IMT findings. This is further illustrated in the supplementary material for out-of-plane disorder where the contrast is even more prominent.
Above observation raises a conceptual question: Inhomogeneities in ∆ d (i) within the IMT calculations arise with a length-scale of ξ, while the fluctuations in ρ follows (uncorrelated) disorder with the natural length-scale k −1 F , the Fermi wavelength. This makes a clear separation of the two length-scales. But, such a distinction is difficult as the spatial profile of one follows the other in the GIMT results! This turns out not so crucial for the cuprates, in which both these scales are of the order of a few lattice spacing. In fact, a recent extraction 82 of the length scales associated with disorder seems to validate the GIMT picture.

D. Spatially correlated effective disorder
The mean-field Hamiltonian of equation (9) can be cast in the following effective form: where the explicit expressions for t eff and V eff can be extracted comparing equation (9) and (19), both of which become inhomogeneous in the presence of the disorder. The inhomogeneities in t eff from the IMT calculation are very weak though, and that in V eff is largely similar to the bare disorder. Obviously, the physics of the strong correlations within GIMT is also contained in the equation (19) through local GRF and their derivatives with respect to ρ i . The physical origin for the differences between the IMT and GIMT results can be motivated using the following argument. A strong repulsive interaction smears out any local accumulation of the charge carriers. Such a smearing plays an important role in GIMT, but is missing from plain IMT scheme. As a result, the spatial distribution of electrons ought to be different between the two calculations. Such re-organization of the electronic density must modify all other local order parameters in a self-consistent manner and hence should also alter all the physical observables. The question remains: How is this correlation driven physics incorporated in the H MF ? A significant insight is obtained by studying the spatial fluctuations of t eff and V eff individually, which is presented in figure 6. We plot t eff (i) = δ t eff (i, δ) and V eff (i) on the lattice both from the IMT and the GIMT calculations. Both the t eff (i) and the V eff (i) show qualitatively similar spatial variations in plain IMT ( figure 6(a) and (b)). On the contrary, a strong spatial anti-correlation develops in the fluctuating parts of t eff and V eff within the GIMT scheme as demonstrated in figures 7(b3) and (b4) using scatter plots. Comparison of their spatial profiles in figures 6(c) and (d) shows that t eff (i) is small precisely in those regions where V eff (i) is large and vice versa. Such a competing behavior 83 of the two is a result of the strong repulsive interactions, and we elaborate this point in the following.
According to equation (19), t eff (i) is the probability of hopping of an electron onto the site i from the neighbors. A site with V eff (i) > V eff (i + δ) will support ρ i < ρ i+δ leading to a larger g t i,i+δ (compared to when V eff were same for i and i + δ). This, in turn, enhances t eff (i), increasing the final ρ i compared to what V eff (i) alone would have predicted in the first place! So high 'hills' of disorder profile does not remain as sparsely populated by electrons in GIMT as expected in the IMT calculation. Exactly opposite mechanism smears out the charge accumulation in GIMT from sites with deep 'valleys' of disorder potential.
On the other hand, t eff (i) is largely independent of ρ i within plain IMT framework, though it develops a very weak densitydependence through the Fock-shifts. In fact, t eff (i) becomes weaker in regions with high-hills and deep-valleys of disorder profile within plain IMT. Such sites tend to get decoupled from hopping in the absence of the 'feedback-loop' in terms of g t i,i+δ , quite in contrast to the GIMT scenario discussed above. This picture is substantiated through our results in figures 6(e) and (f), where we demonstrate the competing effects of t eff (i) and V eff (i) on ρ i , if the two were to work individually. It is thus no surprise that the GIMT local density (Left column, figure 5(b), incorporating both t eff and V eff ) becomes less inhomogeneous than the IMT finding.
Finally, we note that one of the important effects of the strong correlation lies in the significant reduction 58 of the magnitude of V eff (i) from the bare V i . We demonstrate this by presenting the distribution of P (V eff ) in figures 8(a) and 8(b) both from the GIMT and plain IMT calculations. The homogeneous component of V eff (the contribution at V = 0) from the Hartree-shift is already subtracted in this analysis. We see that while in the IMT calculations the range of V eff (i) actually increases a bit from the bare value of ±V /2, strong electronic correlations indeed cause its drastic reduction as seen from the resulting P GIMT (V eff ). Interestingly, P GIMT (V eff ) becomes more asymmetric with increasing V with larger weight accumulating for negative V eff (i). In contrast, the P IMT (V eff ) maintains its symmetry around zero just as in the bare P (V i ) by construction. Nevertheless, we found that the asymmetric P GIMT (V eff ) still produces V eff = 0 for all V . In fact, it is this asymmetry that translates into an oppositely asymmetric P GIMT (ρ i ) in figure 1(a). We comprehend this intriguing feature in the following way: For a non-interacting problem at large V , the local density ranges from zero to doubleoccupancy on the highest hill and the deepest valley, respectively. The GA, however, restricts ρ i ≤ 1, breaking the symmetry of the non-interacting problem. A self-consistent interplay of the strong correlations and the disorder now ensures that the effective disorder has more 'attractive' points such that the maximal ρ i is limited by half-filling everywhere, still maintaining the global density at the desired value (ρ = 0.8, in our case). An evolution of spatial correlation in V eff (i) with disorder, producing larger patches of attractive points is further illustrated using colour-density plots in figure 7(a), indicating a growing spatial correlation of such attractive points. Let us emphasize that we verified the crucial physical insights, emerging from the study of the spatial distribution of the local observables, survive in all the realizations of disorder we considered.

E. Underlying normal state: Anderson's theorem for cuprates?
The intriguing evolution of the spatially correlated effective disorder in the GIMT calculations immediately tells us that they must also modify the properties of the underlying nor-mal state (NS). Here, the NS is defined as the solution of the same H MF of equation (19), but suppressing all the Cooperchannel order parameters ∆ δ i , for all i and δ. We caution the readers that our usage of the term 'normal state' might lead to a confusion, because it is frequently used to identify the T > T c phase of a HTSC, where the d-wave order is already destroyed by the thermal fluctuations. Our NS, on the other hand, is a true GS at T = 0, had there been no pairing. As the locally renormalized GRF leads to spatially correlated disorder, the NS GIMT must differ from the NS IMT , where the later is essentially the solution of the Anderson model 84 of noninteracting electrons in disordered background. The NS GIMT plays a crucial role for comprehending our GIMT results hinting towards a physical picture for the dirty cuprates, as we discuss below.
The robustness of the observables to the impurities, as demonstrated in section III B, is actually reminiscent of Anderson's theorem. 30 While the theorem does not apply to a weak-coupling d-wave superconductor, 14 our results open up such a possibility for strongly correlated cuprates. In order to support this scenario, we first argue that our GIMT results are consistent with Anderson's mechanism provided we extend the scope of pairing between all the 'correlated' singleparticle normal states (NS GIMT ). To establish this, we solve for H Normal MF (equation (19) with all ∆ δ i terms suppressed) to obtain the NS GIMT parameters t eff (i) and V eff (i) in a selfconsistent manner in the first step. These parameters carry information on the local orders in and Hartree-and Fockchannels. In the next step, we set up a BdG-type matrix H BdG in terms of these NS GIMT parameters in the diagonal  (21), with P (∆ d ) taken directly from the GIMT results of figure 1(c). The strong resemblance of the resulting N (ω), particularly for small |ω|, with that in figure 3(a) demonstrates the irrelevance of the spatial inhomogeneity. blocks, and with ∆ δ i in the off-diagonal blocks, where ∆ δ i are obtained from the corresponding GIMT solution of the full H MF , including the dSC order. ThisH BdG is then diagonalized only once in the final step, and we recalculate the new (non self-consistent) local pairing amplitudes∆ i δ from this sole diagonalization. We find an impressive agreement between ∆ i δ and∆ i δ for all i and δ. Such a non self-consistent calculation is equivalent to an Anderson-type idea (though the pairing in this case is not limited to the time-reversed states alone). This interpretation works because the effective attraction does not seem to renormalize the spatially fluctuating NS orders in the diagonal blocks ofH BdG . If instead, we had any significant mismatch between ∆ i δ and∆ i δ , a feedback loop would have been required for self-consistency which could, in general, update NS order parameters as well, implying that the pairing modifies single-particle states. Such a conclusion goes against Anderson's philosophy, according to which pairing from effective attraction is independent of disorder. This is because, the impurities are already consumed in producing the corresponding "exact eigenstates". Note that in the case of the sSC, the normal states are essentially the exact eigenstates of the effective single-particle Hamiltonian, as there are no strong correlations.
It is not only the pairing between the normal states which is responsible for the insensitivity of superconductivity to the impurities. The extent of localization of these states is also crucial for the inferences made by Anderson's theorem as we discuss in the following. While such a pairing theory for a sSC can be formulated in principle for any strength of disorder tailoring the disorder-induced inhomogeneities in it, the robustness of s-wave superconductivity breaks down at large disorder (V ≥ t) due to strong localization. 59 It is interesting to note that the insensitivity of observables for dirty cuprates actually extends even beyond V > 2.5 (See figures 2 and 3), a strength by which traditional Anderson's theorem (for sSC) weakens drastically. In order to settle this apparent conflict, we study the differences between the NS from the GIMT and IMT schemes, focusing on their localization properties. We quantify the extent of localization by calculating the participation ratio (PR) for the effective one-particle NS near µ, defined as: where ψ n (i) is the normalized eigenfunctions of the NS for the eigenvalue ε n ≈ µ. For a better statistics, we actually average PR over all those n's such that |ε n − µ| ≤ 0.03. The PR measures the fraction of sites of the system contributing to the nth wave function, and is a good estimate of the localization area in two-dimensions. The behavior of PR with the disorder strength V in figure 8(c) shows that the one-particle NS GIMT (with ε n ≈ µ) remains fairly extended up to V = 1.75 on a showing LDOS on the impurity site i0 and also on i0 + δ. Such a gap-filling heals quickly over a short distance, and LDOS follows the profile at V = 0 on the sites beyond the next-nearest neighbor from i0. N (i0 + δ, ω ≃ 0) was found to increase with V0. scale of the system size. In contrast, the NS IMT , pairing of which leads to the Anderson's theorem for sSC, gets localized relatively fast! The normal state DOS, N NS (ω ≈ µ), that incorporates the GA was found to be insensitive to V (though shows a strong sample to sample fluctuations), 86 consistent with the extended nature of the NS GIMT . The robustness of observables for dirty superconductors follows as a consequence of traditional pairing philosophy of Anderson, as long as the exact eigenstates (or NS in our terminology) remains fairly extended in the scale of the system size. A similar picture emerges for the dirty cuprates from the above discussions: (a) The GIMT findings can be thought of as a result of pairing between correlated normal states. (b) Such normal states remain delocalized up to a large disorder.
An important consequence of the conventional Anderson's theorem 30 is the irrelevance of any detailed spatial information of the system. Such irrelevance for our case is seen by calculatingÑ (ω) using the following prescription: 85 where, N BCS is the density of states of a homogeneous d-wave BCS superconductor and is given by 13 Here Before closing the current discussion, we must mention that the above arguments supporting the robustness of cuprate superconductors as being similar to the conventional Anderson's theorem, are all indirect ones, though they add up to a significant claim! It will be interesting to substantiate these claims also by direct demonstration. Such a formalism is beyond the scope of this study and will be addressed elsewhere.
It is important to address the fate of the cuprates as the disorder strength is increased beyond what is reported here. It is easy to see that the GIMT scheme, which removes all doubleoccupancy by enforcing ρ i ≤ 1 (for all i), would not work on the sites with V i ≥ U . However this is not a deficit of the GA. It would only require a modification of allowing occasional ρ i > 1 based on relative strength of U and V i . Such calculations, however, are beyond the scope of our study and will be addressed elsewhere. Nevertheless, we found that the GIMT self-consistency cracks up technically even for V > 3.0 (boxdisorder). In fact, making it work for large V requires significant improvements in the iterative self-consistency in terms of applying accelerated convergence methods. Our usage of Broyden 87 and modified Broyden 88 methods was crucial in pushing the self-consistency to larger disorder strengths.
It turns out that conc-disorder helps to push the selfconsistent calculations towards a higher V 0 , because we could tune n imp independently to keep iterative self-consistency under control. The average DOS calculated with V 0 = 5 and n imp = 0.03, 0.05 are presented in figure 9(a) which shows the first evidence of a gap-filling within the GIMT scheme. 89 Figures 9(b) and 9(c) present the LDOS from the GIMT cal-culation around a single impurity, where the strength of the sole impurity is V 0 = 7 and V 0 = 10, an energy scale larger than t, t ′ , J and J ′ . Results for both V 0 = 7 and V 0 = 10 show that growing V 0 produces larger N (r, ω ≃ 0), where r represents the nearest neighbors of the impurity site. A simple picture emerges from repeating the calculations for several combinations of V 0 and n imp : The gap filling begins when V 0 grows beyond the bandwidth of the homogeneous, Gutzwiller-renormalized normal state, i.e., V 0 > 8g t t, where g t = (1 − ρ)/(1 − ρ/2). This observation motivates the formalism of an alternate IMT scheme elaborated further in the supplementary material.

IV. CONCLUSION
In conclusion, we have studied the qualitative and quantitative effects of the disorder-induced inhomogeneities on a strongly correlated dSC. We established that such fluctuations in the spatial orders are fundamentally different from when strong correlations are neglected. Such differences in GIMT do not allow formation of superconducting-'islands', though produces nanoscale inhomogeneities. The surprising insensitivity of ∆ OP , D s and N (ω) to impurities is reminiscent of Anderson's theorem, allowing a simple understanding of the complex physics of dirty cuprates. Even though the spatial inhomogeneity in the order parameters are somewhat weaker from plain IMT, the outcome of our GIMT results remain substantially different from the prediction of the Abrikosov-Gorkov theory or self-consistent T-matrix approximation. The additional effects of strong correlations in the GIMT scheme make the underlying normal states substantially different from the conventional metallic ones. The GIMT normal states differ from the localized phase of the Anderson model because it generates spatially correlated effective disorder. The superconducting GS within the GIMT scheme is interpreted as a result of pairing between these correlated normal states. While a large part of the differences between the GIMT and IMT can be attributed to the changes in the effective disorder, they cannot complete the story just by themselves. A complex and correlated interplay of the effective disorder and pairing modifies the scenario in a self-consistent manner. Our results from the GIMT method help comprehending much of the experimental trends. Extension of the GIMT scheme to variationally include inhomogeneous double-occupancy (consistent with local disorder) is an important future direction. It will be interesting to consider the effect of quantum phase fluctuations. While they play an important role on a dSC, 90 our GIMT results of the lesser fluctuation in the pairing amplitude call for revisiting the effects of the quantum phase fluctuations riding on top of the GIMT state in a dSC. We speculate that our results will motivate further research to address surprises from the dirty superconductors.
(J → J ′ , δ →δ) for n = 1, 2, 4 and (δ ′ ) is included in the superscript, wherever applicable. The Fock-shifts in equation (9) are given by, where τ δ i (≡ τ ij ) lives on the bonds, here j = i + δ. The Hartree-shift in equation (9) gets its major contribution from the derivatives of GRF with respect to local density:  FIG. S1. The spatial profile of ρi and ∆ d (i), for out-of-plane 'smooth' disorder with nimp = 0.03 and V OP = 1. The left column presents the GIMT results for ∆ d (i) (top) and ρi (bottom), and the IMT results for the same appear on the right panels. While the main features are qualitatively similar to those from the box-disorder in figure 5(a) and figure 5(b) in the main paper, the extended nature of disorder in this case produces smooth variations of the order parameters. The similarity between the spatial structures of ρi and ∆ d (i) in the GIMT results is evident! The two structures, however, are quite different in the IMT results which can be understood by looking into the regions with larger local density (ρi > 1, magenta in the bottom-right panel), which weaken ∆ d in that neighborhood. Resulting SC-'islands' are restricted only to the annular regions separating large and small ρi. In contrast, these impurities produce an impression of 'islands' in the GIMT calculations, but both for ρi and ∆ d (i).
e.g.,Ṽ /t = V /t andJ/t = J/t. This automatically explains the reduction of the effective disorder in the GIMT calculation, because,Ṽ ∼ V × (t/t)! In fact, the observation of gap filling in N (ω) for V 0 > 8g t t (see section III E in the main paper) motivates this alternative IMT scheme, which already incorporates the reduction of V eff (i). Similar idea is extended to our "t − J" model of equation (2) of the main paper. such that all the order parameters be identical to those from the GIMT at V = 0. We find the gap-filling to be weaker compared to the previous IMT findings, See figure 3(b) in the main paper. Significant differences with the GIMT results of figure 3(a) of the main paper, however, persist. The inset compares ∆OP from the two IMT schemes and the GIMT scheme.