Correlated disorder as a way towards robust superconductivity

Ordinary superconductors are widely assumed insensitive to small concentrations of random nonmagnetic impurities, whereas strong disorder suppresses superconductivity and even makes superconductor-insulator transition occur. In between these limiting cases, a most fascinating regime can take place where disorder enhances superconductivity. Hitherto, almost all theoretical studies have been conducted under the assumption that disorder is completely independent and random. In real materials, however, positions of impurities and defects tend to correlate with each other. This work shows that these correlations have a strong impact on superconductivity making it more robust and less sensitive to the disorder potential. Superconducting properties can therefore be controlled not only by the overall density of impurities and defects, but by their spatial correlations as well. Strong random disorder is known to suppress superconductivity. In this work, authors showed that when impurities are correlate with each other, superconductivity is more robust and thus its properties can be controlled by spatial correlations of impurities and defects.

F ounding theoretical studies of Abrikosov and Gorkov 1 of superconductivity in a disordered material were followed up by a general physical consideration by Anderson 2 according to which superconductivity with the s pairing symmetry is insensitive to weak nonmagnetic disorder. Subsequent experimental studies, however, demonstrated that superconductivity is suppressed in strongly disordered samples. A very strong disorder can lead to a metal-insulator transition in the normal state, to the appearance of a pseudogap in spectrum, to larger spatial fluctuations of superconductive pairing, and increased Δ/T c ratio [3][4][5][6][7][8][9][10] . Theoretical investigations have largely supported the experimental results [11][12][13] .
The recent advent of the quasi-2D materials era has opened up new horizons in physics of disordered systems including superconductors. Many phenomena important for electrical transport, such as renormalization of the electron-electron interactions, Anderson localization, and phase fluctuations, are boosted in lowdimensional structures, resulting in suppression of superconductivity. The detrimental effects of lowered dimension in combination with disorder were observed in amorphous as well as highly crystalline 2D superconductors [14][15][16][17][18][19] .
The disorder can also result in a remarkable enhancement in superconductivity figures of merit. The disorder potential controls the interplay between the superconductive pairing and longrange phase coherence. The stronger disorder increases spatial inhomogeneity, which enhances the local pairing correlations and superconducting gap, comparing with the clean system 20 . Recent studies have revealed an unprecedented disorder-induced increase in the superconductivity pairing instability in quasi-1D single crystals made of MoSe chains weakly bound by Na atoms 21 . An anomalous superconductivity enhancement due to large structural disorder was also reported for TaS 2 monolayers 22 . Disorder-related effects are assumed responsible for a large increase of the critical temperature in the recently discovered superconducting NbSe 2 monolayers 23 . Theoretical analysis attributes the enhancement to the disorder-induced multi-fractal structure of the electronic wave functions 24,25 , which was observed by numerically solving microscopic theory equations in low-dimensional samples 26 . The calculations predict the cluster formation and survival of the spectral gap even at extremely large disorder 27 .
On the other side, strong disorder enhances phase fluctuations, thereby reducing the superfluid density (stiffness), and suppressing superconductivity on the global scale. With the changing disorder strength the system can, therefore, have an optimal degree of inhomogeneity, where the superconductivity enhancement and the critical temperature T c are maximal 28 . Notice, that a similar mechanism underlies superconductivity enhancement in samples where inhomogeneities are related to the quantum size effects [29][30][31] .
Introducing random impurities and defects in otherwise ordered materials can thus be regarded as a tool to control superconducting characteristics of materials, a design element to engineer superconductors with desired functionalities. However, one of the main problem with this approach is to identify key parameters of the impurity distribution which affect the superconductive state. Solving it requires detailed numerical simulations that account for the disorder distribution in real materials.
Most theoretical investigations of disordered superconductors are based on the models with the spatially uncorrelated disorder, which are analyzed using the perturbation expansion methods. In real systems, however, the disorder is almost never completely random. The inhomogeneities are often arranged in a certain structure, characterized by the long-range spatial correlations. There is a growing appreciation of the fact that such correlations can change essential properties of disordered systems qualitatively [32][33][34] .
Prominent examples of this type include opening resonant transmission channels 35 , modifying metal-insulator transitions 36 , and changing mobility edges 37 . Effects related to the disorder correlations have been experimentally observed in gases of cold atoms 38,39 and photonic systems 40 . Also, correlated deviations from the ideal periodicity of a crystal structure plays a significant role for functional materials 41 that utilize ionic conductivity 42 and their ferroelectric 43,44 , thermoelectric 45 , and photoelectric properties 46 . Achieving desired functionalities by manipulating patterns of the structural disorder is one of the goals of contemporary material science studying disordered crystals 33,47 .
In recent experiments with dirty superconducting films 23 , a visual analysis of the disorder distribution reveals spatial longrange correlations. It is clear, the superconducting state will be much more complex when the disorder is spatially correlated. The correlations introduce an additional characteristic lengthscale that competes with those defined by the BCS coherence length ξ and the Fermi wave length λ F . It has been previously shown that superconductivity can be strongly boosted in both of the opposite limiting cases: when the external potential is fully random and uncorrelated 26,27 , and when it has a fixed ordered structure without the random component 20,48,49 . However, very little if anything is known about a most relevant situation, where an inherently random distribution of material inhomogeneities acquires a certain degree of spatial correlations. The goal of this work is to initiate studies of an interplay of "order and disorder" in superconductors, with a strategic goal to find the optimal balance for the most robust superconducting state.
Here, we investigate the influence of the long-range correlations of the disorder potential on the superconducting state at zero temperature. The analysis is done by calculating superconducting properties of 2D sample using a disordered model with spatial correlations. Our analysis shows that long-range correlations can considerably change statistical properties of the Cooper pairing and enhance multi-fractal features of the order parameter. As the result, global superconductivity becomes more robust with respect to the disorder strength.

Results and discussion
The calculations are done using the model of the structural disorder where the spatial correlations of the disorder potential in the inverse space obey the power law in the limit of small q. The power exponent α determines the correlation degree. When α = 0 the random potential V i at different lattice sites is fully uncorrelated. The spatial correlations of the random potential in the real space are characterized by the correlation length ξ V calculated in the Methods section.
Spatial profile of the disorder. We start the analysis with showing in the row (a) of Fig. 1 typical spatial profiles of a disorder potential V i obtained for different values of the correlation degree α = 0, 1, 2, 3. For the fully uncorrelated disorder with α = 0 the color density plot representation for the profile appears as fully randomly distributed array of "pixels" or grains. When α increases the pixels are smoothed out gradually forming textures or structures of larger dimensions, on the scale ξ V (see Fig. 2a). It is to be noted that positions, shapes and orientation of those textures are still random, such that the ensemble-averaged system remains homogeneous (if the number of disorder realization is large enough).
Spatial profile of the order parameter. Using the above configurations of the disorder potential we solve the Bogoliubovde Color density plot of the spatial profile of the random potential V i (row a) and of the order parameter Δ i (rows (b-e)), calculated assuming the disorder strength of V ¼ 1:0; 1:5; 2:0; 2:5 and the disorder correlation degree α = 0, 1, 2, 3 (panel columns). Blue and white colors denote, respectively, superconducting S and normal N domains. Each panel has its own scaling chosen for a better contrast and visibility of the regions with weak superconductivity: dark blue corresponds to the gap value Δ i > 0:1Δ i , where Δ i is the average gap value. First, we consider the case of uncorrelated disorder investigated extensively in earlier works 11,12,27 . It is well known that when the disorder is weak the disorder parameter is homogeneous. This trivial case is not shown in Fig. 1. However, the panel obtained at relatively small strength V ¼ 1 and α = 0 demonstrates that the order parameter is close to being homogeneous, and the normal phase occurs only in a restricted number of small-size domains.
When the strength V increases, the superconducting state becomes inhomogeneous. Blue superconducting (S) clusters with larger values of the order parameter are mixed with white normal state (N) regions. The clustering in the order parameter profile appears not strongly correlated with the disorder potential: one sees in Fig. 1 that the disorder profile at α = 0 shows no textures with sizes comparable to clusters of order parameter V > 1. The clustering is a manifestation of the multi-fractal properties of the BdG solutions. The degree of correlations between the potential and the order parameter can be described by the statistical correlator where V i ¼ jV i þ U i À μ i j, and which characterizes the suppression of superconductivity by the random potential. Results for the calculations are shown in Fig. 2c. At α = 0 and V ¼ 1 one obtains C ΔV ' À0:5 which corresponds to a relatively weak negative correlation between the gap Δ i and potential V i .
The order parameter profile, shown in Fig. 1 for α ≠ 0, demonstrates how it changes with the disorder potential strength V and the correlation degree α. First, one notes that the order parameter becomes more homogeneous at larger α: sizes of normal clusters decrease making the system "more superconductive". At large values of α, a typical size of the order parameter cluster (N or S) is comparable to that of textures of the disorder potential. This implies the order parameter and disorder potential become more correlated, which is demonstrated in Fig. 2c.
The increased correlations is a consequence of an intuitively obvious fact that superconductivity is most likely to exist in valleys of minimal values of V i . For weakly correlated disorder, superconductive clusters, determined by the BCS coherence length ξ, are much larger than disorder potential textures, determined by ξ V , reducing correlations between the order parameter and disorder potential. However, when disorder is strongly correlated, sizes of the clusters and textures are comparable and the correlation increases.
This visual analysis of results in Fig. 1 demonstrates two opposite tendencies. A strong disorder with larger V suppresses superconductivity, which is accompanied by the increased inhomogeneity and clusterization with the larger area of N islands. In contrast, a larger degree of correlation α enhances superconductivity and smears out the clusters giving rise to a more homogeneous profile.
Order parameter distribution. Statistical properties of the order parameter Δ i are visualized by its absolute value distribution, shown in Fig. 3a-c for a few values of V and α. The distribution is plotted as a function of the relative order parameter Δ= Δ where Δ ¼ jΔ i j is the averaged order parameter. In a superconductor with uncorrelated disorder, this distribution is well described by the log-normal function 11,12,27,50 PðΔÞ  The distribution changes itself when the disorder is correlated and α ≠ 0. In the weak disorder in Fig. 3a, the distribution is still well approximated by Eq. (3) for all considered values of α. The maximum is, however, shifted to larger Δ reflecting the already noted trend [ Fig. 1] that increasing disorder correlations and amplitude act in the opposite directions. This trend can be seen also for strong disorder. Figure 3b, c demonstrate the maximum is shifted to larger Δ= Δ when α increases. However, when the disorder has large correlations and is strong, i.e., both α and V are large, P(Δ) deviates from the lognormal form (3) at small Δ. For α = 3, the deviation is visible at V ¼ 1:5 [ Fig. 3b], while at V ¼ 2 its is clearly seen also for α = 2 [ Fig. 3c and the inset]. Notice that for large α and V, P(Δ) deviates qualitatively from Eq. (3). It approaches a finite value in the limit Δ → 0 and have the opposite slope comparison to the log-normal distribution.
The changes in the distribution P(Δ) with V and α are quantitatively characterized by calculating the width σ and Δ max position of its peak, shown in Fig. 3d and e, respectively, as functions of α at several values of V. The results reveal the same opposite trends. The width σ increases at larger V and decreases at larger α, and the peak position Δ max decreases at large V and increases at large α. Finally, Fig. 3f illustrates the superconductivity suppression and the appearance of N-state clusters in the order parameter profile by plotting the probability to have (close to) zero order parameter P(Δ → 0) (here we use the criterion Δ < 0:1 Δ). The figure shows this probability increases with raising V and decreases when α increases.
This analysis of the order parameter distribution and its defining characteristics quantify visual changes seen in the order parameter profile in Fig. 1.
Superconductive correlations. We now turn to spatial correlations of the order parameter. As mentioned in the introduction a weakly disordered system is naturally characterized by the competition between the BCS coherence length and the disorder correlation length. However, for the strong disorder the superconducting state clusterizes becoming inhomogeneous. In this case, the system is characterized by spatial correlations of the order parameter with the correlation function defined as where u ðnÞ i and v ðnÞ i are eigenfunctions of the HF-BdG Hamiltonian. This function describes correlations on the longer scale. In a homogeneous superconductor state it becomes constant equal to j Δj 2 at large distances, reflecting the off-diagonal long-range order 51 in the mean-field approximation. The disorder destroys the long-range order and the global superconducting correlations. The corresponding characteristic length is found as In a clustered superconducting state in a disordered sample ξ Δ can be interpreted as a measure of correlations between different clusters. It differs from the much shorter BCS coherence length, which defines the superconducting pairing scale, associated with the correlation function which has a meaning of a squared absolute value of the Cooper pair wave function |Ψ(i, j)| 2 , averaged over disorder realizations. The BCS coherence length ξ is calculated substituting this expression in place of S Δ in Eq. (5) 27 . For a clean system this expression yields a standard result for the BCS coherence length ξ 0 . For a dirty superconductor with uncorrelated disorder this expression gives a result close to the perturbation theory [see Supplementary Note 3]. The spatial dependence of S Δ (r), calculated at V ¼ 1 for several values of α, is illustrated in Fig. 4a. As expected, the correlation drops at large distances due to disorder. However, it saturates at large distances leaving a residual correlation across the entire sample. The value of this residue increases notably with raising degree of the disorder correlation α.
A quantitative measure of the long-range correlations described by S Δ (r) can be extracted from the small-q asymptotic of its Fourier transform. It is shown in Fig. 4b in the log-log scale. At q → 0, it gives a linear dependence, so that its asymptotic is also the power law S Δ (q) ∝ q −β . By fitting the linear dependence we obtain the power exponent β as a function of α, which is shown in the inset in Fig. 4b. β is small when α ≲ 1, but it rises sharply at α ≳ 2. The rising β means that the long-range superconducting correlations grow rapidly with the correlations in the disordered potential. In principle, this is consistent with the assumption that α = 2 is the border value separating two different regimes in the chosen model for disorder correlations. The numerical dependence is well fitted by the expression β = β 0 + (α/4) γ with β 0 = 1.25 and γ = 2.75.
Correlation lengths ξ and ξ Δ in Fig. 5 illustrate their dependence on V and α. Comparing them with Fig. 2a confirms the expected relation ξ V < ξ < ξ Δ , which holds for all considered parameters. For uncorrelated vanishingly small disorder V ! 0, ξ yields its clean system limit (σ 0 ≃ 9), whereas ξ Δ approaches the effective sample dimension which is L= ffiffi ffi 6 p ' 16. When V increases, both ξ and ξ Δ reduce monotonically. The rising correlation degree α results in the opposite trend making both lengths monotonically increase. Here too, the disorder strength and correlation degree act on the correlation lengths in the opposite way.
Superfluid stiffness. In a strongly disordered material, the definition of superconductivity using characteristics of the local microscopic state is inadequate. In a weakly disordered sample the advance of superconductivity is unequivocally connected to the non-zero order parameter. However, in the case of strong disorder the superconductive state is strongly non-homogeneous [ Fig. 1] forming weakly connected superconductive clusters with randomized phases, and the definition of superconductivity via the local or average order parameter is no longer possible. The onset of superconductivity in a finite strongly disordered sample can be characterized by the (Meissner) superfluid stiffness D 0 s 52,53 , which is a coefficient in the expression for the "phase rigidity" term in the energy of the superconducting condensate where θ is the phase of the superconducting order parameter. From the macroscopic point of view this is an elastic energy associated with the twisted phase of the superconductive condensate. A non-zero stiffness D 0 s means a spatial variation of the phase requires additional energy and thus the phase is rigidly fixed ("stiff"). It is the rigidity of this phase that endows a superconductor with the ability to sustain a super-current.
Within the mean-field theory the stiffness is determined by the linear response to an externally applied static vector potential 52 , where Λ xx is the long wavelength limit of the transverse currentcurrent correlator averaged over both the superconductive state and the disorder realizations, with j x q; τ ð Þ being the paramagnetic current and the diamagnetic component ÀK x is the disorder-averaged kinetic energy along the x direction. The calculation equations for D 0 s within the BdG approach are presented in the Supplementary Note 2.
The mean-field theory calculations can be further improved by taking into account phase fluctuations within the effective XY model and the self-consistent Harmonic approximation. This gives the renormalized stiffness as 12 where the mean square fluctuation of the nearest neighbor phase difference is obtained as with κ being the mean-field compressibility κ = dn/dμ and ϵ q ¼ 2½2 À cosðq x Þ À cosðq y Þ the single-particle spectrum. Figure 6a plots the stiffness as a function of the disorder strength V for a few values of correlation α. The figure gives both full D s (full circles) and the mean-field result D 0 s (open circles) for comparison. As is expected intuitively, the fluctuations reduce the stiffness and suppress superconductivity. Notice, that Eq. (11) does not account for all contributions of the quantum fluctuations. However, the role of the fluctuations is expected to decrease when the disorder becomes correlated and α grows [see detailed calculations in Supplementary Note 4]. This implies that the approximate result in Eq. (11) is sufficient to qualitatively describe the α-dependence of the stiffness.
In agreement with previous works 12,27,54 for the noncorrelated disorder (α = 0), the stiffness declines when the disorder strength V grows [ Fig. 6a]. Here we obtain the similar dependence also for correlated disorder (α ≠ 0). However, increasing correlation degree α results in the larger stiffness with a consequence that the critical disorder strength, where the stiffness becomes zero, increases at larger α.
It must be noted, the calculation of stiffness shows the global superconductivity ceases at much smaller values of V then anticipated from the order parameter profile and its distribution. For example, the stiffness becomes zero at V ≳ 0:75 at α = 0. In contrast, the order parameter is non-zero in almost the entire sample at V ¼ 1:0 and α = 0 [ Fig. 1]. It is further confirmed by the distribution P(Δ) [Fig. 3a] and its value at Δ ≃ 0 [ Fig. 3f], both showing the order parameter is zero practically nowhere in the sample. One sees, globally non-zero order parameter is not sufficient to ensure global superconductivity.
The relations between the increased stiffness and larger disorder correlation are also traced when looking at the distribution in Fig. 3a-c. The stiffness is known to decrease due to phase twists at places where it costs less energy, i.e. where the order parameter and the condensate density n s are minimal. When the correlation degree α increases, the distributions in Fig. 3a-c shift towards larger Δ's, leading to the larger stiffness.
It is worth to consider diamagnetic ÀK x and paramagnetic Λ xx contributions to the stiffness separately [ Fig. 6b]. ÀK x is a slowly decreasing function of V. It does not depend on α because it is proportional to the total electron density n e .
The paramagnetic contribution Λ xx behaves differently. For the clean system with V ¼ 0, the total current in the system is proportional to its momentum. Therefore, the respective current operator commutes with the Hamiltonian and the currentcurrent correlator contributing to Λ xx ðq x ¼ 0; q y ! 0; ω ¼ 0Þ vanishes. In a disordered system, the operator of the total current no longer commutes with the Hamiltonian, and Λ xx ≠ 0. The value of Λ xx depends on the degree of the non-commutativity between the current and Hamiltonian operators, controlled by disorder strength and correlations [ Fig. 6b].
A summary of the stiffness calculations is presented in Fig. 6c as a phase diagram in the α À V plane. The region of non-zero stiffness is marked by blue (S), and the region of zero-stiffness is denoted by red (N). Solid and open points are obtained using the fluctuation corrected (D s ) and bare (D ð0Þ s ) stiffness. The phase diagram in Fig. 6c highlights a general tendency: the correlations inhibit the destructive influence of disorder on superconductivity. The dashed lines, connecting the points, serve as guides to an eye, but are expected to follow the real crossover lines, separating S and N domains, closely. Notice, they increase monotonically with rising α. Their slope increases at α ≳ 2, the special value for the model.

Conclusions
This work studies superconductivity in systems with spatially correlated disorder, by solving a full set of microscopic HF-BdG equations for a finite 2D sample. The long-range correlations of the disordered potential are taken into account using the model where the Fourier components of the potential have random phases. Only the case of zero temperature is considered.
The calculations reveal a very notable role the disorder correlations play in shaping superconductivity, affecting its properties both on the local and global scales. It was shown that the superconductivity is enhanced by increasing the degree of disorder correlations. It is manifested in the changed character of spatial profile, statistical distribution and spatial correlations of the order parameter. The statistical distribution demonstrate that multifractal features of the HF-BdG eigenstates, reported earlier for uncorrelated disorder, hold also when the disorder becomes correlated.
The correlations in the disorder potential enlarge both the BCS pairing correlations and the spatial correlations of the order parameter. The sizes and configuration of clusters in the inhomogeneous superconducting state are shown to depend strongly on the correlations degree. The global superconductivity of the sample, quantified by the superfluid stiffness, is also boosted by the disorder correlations.
Although our results are obtained for a relatively simple model of a non-structured disorder, they clearly demonstrate that spatial correlations in the disorder potential are a very important factor affecting superconductivity. It can change the whole range of superconductive properties and cannot be ignored when discussing the disorder-related effects in such materials. It is also clear, that deliberate manipulations of the disorder correlations, when possible, can be used to manipulate superconducting characteristics. This adds another design parameter to create superconductors with desired functional properties.

Methods
Microscopic equations for superconductivity. To describe superconductivity with the s-wave coupling we employ the BCS Hamiltonian for the tight-binding model for the single-particle stateŝ whereĉ i denotes the electron operator at site i of the lattice, t ij is the tunneling amplitude between sites i and j assumed non-zero t ij = − t only for the nearest neighbors, g > 0 is the on-site BCS coupling constant, and V i is the impurity potential. Here we consider a 2D lattice with the vector indices denoting sites i ¼ ði x ; i y Þ. The mean field Bogoliubovde Gennes (BdG) equations for the Hamiltonian (12) write aŝ where μ is the system chemical potential, u ! and v ! are vectors with components v i and u i , the operatorsĤ andΔ are defined by their matrix elements with δ ij being the Kroneker symbol. The order parameter Δ i and Hartree potential U i are found from the self-consistency equations where 〈…〉 are quantum mechanical averages. Notice that the Hartree selfconsistency Eq. (15b) modifies the effective potential for electrons in an inhomogeneous system. It plays a crucial role when the disorder is strong and cannot be neglected as in weakly disordered systems 55 . The model neglects the long-range repulsive Coulomb interactions between electrons, as is also done in many previous studies of the role of disorder in superconductive systems 26,27 . Although the influence of the interactions can become notable in disordered systems with localized carriers, it is expected to decrease when the disorder becomes correlated. We also do not consider the d-wave pairing that can be strongly enhanced by combining the Coulomb repulsion with a fully random potential for electrons 56 .
Correlated disorder model. Disordered models with the power-law correlation asymptotic have long been proposed 57 , in particular, to investigate the correlationinduced changes in the Anderson localization 58,59 , and the percolation threshold 60 . Those models, sometimes referred to as speckle, have been applied to model, e.g., of cold gases in randomized optical lattices 38,59,[61][62][63] . The power-law correlations in a spatial distribution of linear or planar crystal dislocations have been recently confirmed experimentally 64 , and the asymptotic (1) agreed well with those observations. The power-law dependence (1) is realized by assuming a model with the random potential V i on a lattice site i given by where r i is the lattice position, q j = (2πj x /N, 2πj y /N) is a discrete inverse space vector, j x,y = 1, …, N and q j = |q j |. Finally, ϕ j are random phases, uniformly distributed in the interval 0; 2π ½ Þ. The spatial correlations of the random potential are defined as with brackets 〈…〉 s denoting the statistical ensemble average. Notice, that for the large system this quantity depends only on the difference r i − r j , but this is violated in practical calculations for a finite system with a limited number of disorder potential realizations. In the inverse space this correlator is obtained by evaluating a sum One can easily see the model (16) Fig. 2b for a few selected values of α (S V (q) depends only on q for a homogeneous system). For all values of α, the log-log dependence is linear with the slope − α in full agreement with Eq. (1). Notice, that V q / ffiffiffiffiffiffiffiffiffiffiffi S V ðqÞ p , where V q is the Fourier transform of the potential defined in Eq. (16) but without the random phases.
Spatial correlations in the real space can be characterized by the correlation length ξ V calculated as The calculated ξ V is shown in Fig. 2a. It increases monotonically as a function of α, indicating a larger correlation length at larger α. At α = 0 this length is ξ < 1, meaning the potential is fully uncorrelated at the most close lattice sites. Finally, it is convenient to shift the potential as to ensure its average value is zero. Here the numerical averaging denoted by the overline is done both over lattice sites and the statistical distribution. This constant shift in the potential excludes the change in the chemical potential. The amplitude of the disorder potential is defined by the quantity referred to as disorder strength.
The model in Eq. (16) was employed to describe real stochastic processes with the long-range correlations including nucleotide sequences in DNA molecules 65,66 , and other systems where disorder correlations affects electronic transport, plasma fluctuations 67 , and patterns in surface growth 68 .
Numerical parameters. Equations (13) are solved numerically on a square lattice with dimensions A = N × N with the periodic boundary conditions. We take N = 40 which is an optimal trade-off between the numerical load and our goal to describe a bulk system, where results are already not sensitive to the sample size. The chemical potential μ is chosen so that the average electron density is set to n e = 0.875. We note that a particular value of n e is not important for as long as we sufficiently far from the half filling of n e = 1 with additional symmetries.
The analysis is done in the limit of nominally large coupling constant of g = 1, and the limit of very large Debye window of ℏω D = 5, which implied the coupling involves all single-particle states. (All energy quantities in the problem are expressed in units of the hopping integral t.) Note, however, that a typical average value of the gap function Δ ' 0:04 is much less than the Fermi energy E F ≃ 3.76.
The ratio of the two quantities Δ=E F ' 0:01, which implies that we are in the weak coupling regime of BCS superconductivity. Finally, the physical quantities of interest are statistically averaged over N s = 50 independent disorder realizations that has the same V and α.

Data availability
The data-sets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.