## Correlated disorder as a way towards robust superconductivity

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

$${S}_{V}(q)\propto {q}^{-\alpha }$$

(1)

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 Bogoliubov – de Gennes (BdG) equations and obtain the order parameter profile. The results for different values of the disorder strength \({{{{{{{\mathcal{V}}}}}}}}=1,1.5,2,2.5\) are shown in the rows (b–e) in Fig. 1, where columns correspond to the selected values of *α*. The blue color indicates domains of the superconductive phase of Δ_{i} ≠ 0 and the white color marks domains of the normal phase with Δ_{i} ≃ 0. The deep blue color indicates \({{{\Delta }}}_{i} \, > \, 0.1\overline{{{{\Delta }}}_{i}}\), where \(\overline{{{{\Delta }}}_{i}}\) is the full-ensemble and lattice average of the order parameter. Note, that each panel in the plot has its own scaling for a better contrast, to highlight details of domains of weak superconductivity.

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 \({{{{{{{\mathcal{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 \({{{{{{{\mathcal{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 \({{{{{{{\mathcal{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

$${C}_{{{\Delta }}{{{{{{{\mathcal{V}}}}}}}}}=\frac{\overline{| {{{\Delta }}}_{i}| {{{{{{{{\mathcal{V}}}}}}}}}_{i}}-\overline{| {{{\Delta }}}_{i}| }\,\overline{{{{{{{{{\mathcal{V}}}}}}}}}_{i}}}{{\sigma }_{{{\Delta }}}{{{{{{{\mathcal{V}}}}}}}}},\quad {\sigma }_{{{\Delta }}}=\sqrt{\overline{| {{{\Delta }}}_{i}{| }^{2}}-{\overline{| {{{\Delta }}}_{i}| }}^{2}}$$

(2)

where \({{{{{{{{\mathcal{V}}}}}}}}}_{i}=| {V}_{i}+{U}_{i}-{\mu }_{i}|\), and which characterizes the suppression of superconductivity by the random potential. Results for the calculations are shown in Fig. 2c. At *α* = 0 and \({{{{{{{\mathcal{V}}}}}}}}=1\) one obtains \({C}_{{{\Delta }}{{{{{{{\mathcal{V}}}}}}}}}\simeq -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 \({{{{{{{\mathcal{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 \({{{{{{{{\mathcal{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 \({{{{{{{\mathcal{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 \({{{{{{{\mathcal{V}}}}}}}}\) and *α*. The distribution is plotted as a function of the relative order parameter \({{\Delta }}/\bar{{{\Delta }}}\) where \(\bar{{{\Delta }}}=\overline{| {{{\Delta }}}_{i}| }\) 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({{\Delta }})=\frac{\bar{{{\Delta }}}}{{{\Delta }}\sqrt{2\pi }\sigma }\exp \left\{-\frac{1}{2{\sigma }^{2}}{\left[\ln \left({{\Delta }}/\bar{{{\Delta }}}\right)-\lambda \right]}^{2}\right\},$$

(3)

where *σ* and *λ* are constants, regarded as fitting parameters. Dashed lines in Fig. 3a–c show the best fit using Eq. (3).

At *α* = 0 Eq. (3) yields an excellent fit for all values of \({{{{{{{\mathcal{V}}}}}}}}\). For weak disorder with \({{{{{{{\mathcal{V}}}}}}}}=1\) [blue line in Fig. 3a], *P*(Δ) is maximal at \({{\Delta }}\simeq 0.5\bar{{{\Delta }}}\). In the limit \({{\Delta }}/\bar{{{\Delta }}}\to 0\), the distribution *P*(Δ) approaches zero [Fig. 3a] which implies the absence of N domains. When the disorder amplitude increases to \({{{{{{{\mathcal{V}}}}}}}}=1.5\), the maximum of the distribution shifts to smaller values of Δ [Fig. 3b], approaching the limit Δ → 0. This corresponds to profiles in Fig. 1, obtained at \({{{{{{{\mathcal{V}}}}}}}}=1.5\) and *α* = 0, where N-state clusters appear.

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 \({{\Delta }}/\bar{{{\Delta }}}\) when *α* increases.

However, when the disorder has large correlations and is strong, i.e., both *α* and \({{{{{{{\mathcal{V}}}}}}}}\) are large, *P*(Δ) deviates from the log-normal form (3) at small Δ. For *α* = 3, the deviation is visible at \({{{{{{{\mathcal{V}}}}}}}}=1.5\) [Fig. 3b], while at \({{{{{{{\mathcal{V}}}}}}}}=2\) its is clearly seen also for *α* = 2 [Fig. 3c and the inset]. Notice that for large *α* and \({{{{{{{\mathcal{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 \({{{{{{{\mathcal{V}}}}}}}}\) and *α* are quantitatively characterized by calculating the width *σ* and \({{{\Delta }}}_{\max }\) position of its peak, shown in Fig. 3d and e, respectively, as functions of *α* at several values of \({{{{{{{\mathcal{V}}}}}}}}\). The results reveal the same opposite trends. The width *σ* increases at larger \({{{{{{{\mathcal{V}}}}}}}}\) and decreases at larger *α*, and the peak position \({{{\Delta }}}_{\max }\) decreases at large \({{{{{{{\mathcal{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 \({{\Delta }} \, < \, 0.1\bar{{{\Delta }}}\)). The figure shows this probability increases with raising \({{{{{{{\mathcal{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

$${S}_{{{\Delta }}}\left({{{{{{{{\bf{r}}}}}}}}}_{i},{{{{{{{{\bf{r}}}}}}}}}_{j}\right)={\langle {{{\Delta }}}_{i}^{* }{{{\Delta }}}_{j}\rangle }_{s}=\mathop{\sum}\limits_{nm}{\left\langle {u}_{i}^{(m)* }{u}_{j}^{(n)}{v}_{i}^{(m)* }{v}_{j}^{(n)}\right\rangle }_{s},$$

(4)

where \({u}_{i}^{(n)}\) and \({v}_{i}^{(n)}\) 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 \(| \bar{{{\Delta }}}{| }^{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

$${\xi }_{{{\Delta }}}^{2}=\frac{{\sum }_{i,j}{S}_{{{\Delta }}}({{{{{{{{\bf{r}}}}}}}}}_{i},{{{{{{{{\bf{r}}}}}}}}}_{j}){({{{{{{{{\bf{r}}}}}}}}}_{i}-{{{{{{{{\bf{r}}}}}}}}}_{j})}^{2}}{{\sum }_{i,j}{S}_{{{\Delta }}}({{{{{{{{\bf{r}}}}}}}}}_{i},{{{{{{{{\bf{r}}}}}}}}}_{j})}.$$

(5)

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

$$S({{{{{{{{\bf{r}}}}}}}}}_{i},{{{{{{{{\bf{r}}}}}}}}}_{j})={\left\langle {\left|\mathop{\sum}\limits_{n}{u}_{i}^{(n)}{v}_{j}^{(n)}\right|}^{2}\right\rangle }_{s}=\mathop{\sum}\limits_{nm}{\left\langle {u}_{i}^{(m)* }{u}_{i}^{(n)}{v}_{j}^{(m)* }{v}_{j}^{(n)}\right\rangle }_{s},$$

(6)

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 \({{{{{{{\mathcal{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 \({{{{{{{\mathcal{V}}}}}}}}\) and *α*. Comparing them with Fig. 2a confirms the expected relation *ξ*_{V} < *ξ* < *ξ*_{Δ}, which holds for all considered parameters. For uncorrelated vanishingly small disorder \({{{{{{{\mathcal{V}}}}}}}}\to 0\), *ξ* yields its clean system limit (*σ*_{0} ≃ 9), whereas *ξ*_{Δ} approaches the effective sample dimension which is \(L/\sqrt{6}\simeq 16\). When \({{{{{{{\mathcal{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}_{s}^{0}\)^{52,53}, which is a coefficient in the expression for the “phase rigidity” term in the energy of the superconducting condensate

$$E\left[\theta \right]=\frac{{D}_{s}^{0}}{2}\int d{{{{{{{\bf{r}}}}}}}}{\left|\nabla \theta \right|}^{2},$$

(7)

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}_{s}^{0}\) 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},

$${D}_{s}^{0}=-{{{\Lambda }}}_{xx}\left({q}_{x}=0,{q}_{y}\to 0,\omega =0\right)+\left\langle -{K}_{x}\right\rangle ,$$

(8)

where Λ_{xx} is the long wavelength limit of the transverse current-current correlator averaged over both the superconductive state and the disorder realizations,

$${{{\Lambda }}}_{xx}\left({{{{{{{\bf{q}}}}}}}},\omega \right)=\frac{1}{N}\int\limits_{0}^{\infty }d\tau {e}^{i\omega \tau }\left\langle {j}_{x}\left({{{{{{{\bf{q}}}}}}}},\tau \right)\,{j}_{x}\left(-{{{{{{{\bf{q}}}}}}}},0\right)\right\rangle$$

(9)

with \({j}_{x}\left({{{{{{{\bf{q}}}}}}}},\tau \right)\) being the paramagnetic current and the diamagnetic component \(\left\langle -{K}_{x}\right\rangle\) is the disorder-averaged kinetic energy along the *x* direction. The calculation equations for \({D}_{s}^{0}\) 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}

$${D}_{s}={D}_{s}^{0}{e}^{-{\left\langle {\theta }_{ij}^{2}\right\rangle }_{0}/2},$$

(10)

where the mean square fluctuation of the nearest neighbor phase difference is obtained as

$${\left\langle {\theta }_{ij}^{2}\right\rangle }_{0}=\frac{2}{N\xi }\mathop{\sum}\limits_{q < \pi }{\left(\frac{{\epsilon }_{{{{{{{{\bf{q}}}}}}}}}}{{D}_{s}\kappa }\right)}^{\frac{1}{2}},$$

(11)

with *κ* being the mean-field compressibility *κ* = *d**n*/*d**μ* and \({\epsilon }_{{{{{{{{\bf{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 \({{{{{{{\mathcal{V}}}}}}}}\) for a few values of correlation *α*. The figure gives both full *D*_{s} (full circles) and the mean-field result \({D}_{s}^{0}\) (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 non-correlated disorder (*α* = 0), the stiffness declines when the disorder strength \({{{{{{{\mathcal{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 \({{{{{{{\mathcal{V}}}}}}}}\) then anticipated from the order parameter profile and its distribution. For example, the stiffness becomes zero at \({{{{{{{\mathcal{V}}}}}}}}\,\gtrsim\, 0.75\) at *α* = 0. In contrast, the order parameter is non-zero in almost the entire sample at \({{{{{{{\mathcal{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 \(\left\langle -{K}_{x}\right\rangle\) and paramagnetic Λ_{xx} contributions to the stiffness separately [Fig. 6b]. \(\left\langle -{K}_{x}\right\rangle\) is a slowly decreasing function of \({{{{{{{\mathcal{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 \({{{{{{{\mathcal{V}}}}}}}}=0\), the total current in the system is proportional to its momentum. Therefore, the respective current operator commutes with the Hamiltonian and the current-current correlator contributing to \({{{\Lambda }}}_{xx}({q}_{x}=0,{q}_{y}\to 0,\omega =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 \(\alpha -{{{{{{{\mathcal{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}_{s}^{(0)}\)) 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.