Introduction

The quest to obtain a superconducting material with a high critical temperature is a task of utmost importance that has the potential to revolutionize all electric applications. Multiple ways to achieving room temperature superconductivity have been and are being explored1. Although the final goal still remains elusive, several strategies to increase the critical temperature appear promising.

One approach is to subject materials to extremely high pressure2. Hydrogen-rich compounds, such as hydrides, have demonstrated superconductivity at very high temperatures under pressure. A groundbreaking discovery was made in 2015 by Eremets and coworkers3, who reported superconductivity at Tc = 203K in H3S under a pressure of 155 GPa. Subsequent efforts achieved the critical temperature of Tc = 243 K for yttrium hydrides4, and Tc = 250 K in lanthanum hydrides5. The increase in the superconducting critical temperature is believed to be due to pressure-induced changes in the electronic bands and phonon spectrum of the material due to light hydrogen atoms.

Another promising direction is to explore superconductivity in layered structures. These include traditional high-Tc copper oxides (cuprates)6,7 and iron-based superconductors (iron pnictides)8. Recently, more exotic layered structures have come into focus, such as twisted bilayer graphene9 and interfaces between complex oxides10. Creating layered heterostructures can induce superconductivity in materials that are otherwise weakly or non-conducting. An example is the interface between lanthanum aluminate and strontium titanate, both of which are bulk insulators11. Interesting examples are also found among ultrathin 2D materials12, such as FeSe monolayers grown on STO substrates13, where a significant enhancement of the superconducting critical temperature has been reported. Superconductivity is known to be enhanced by superstructures9 and geometrical confinement14.

It is also well known that a higher critical temperature, Tc, can be achieved by increasing the density of electronic states (DOS) involved in Cooper pairing. In this context, compounds with a singular DOS near the Fermi surface have attracted significant interest. Notably, these include quasi-one-dimensional superconductors with DOS divergences due to dimensional effects15,16,17,18,19, systems resembling kagome materials20,21,22, where a nontrivial atomic lattice results in the formation of both van Hove singularities and a flat band, and twisted graphene layers, where van Hove singularities arise due to interlayer coupling23. A high critical temperature in systems with DOS divergences can be stabilized using a multiband structure, where interactions between the band condensates suppress superconductive fluctuations24,25.

Introducing disorder by adding impurities and defects can also enhance superconductivity. This may seem paradoxical, as disorder is commonly expected to suppress both conductivity and superconductivity. Strongly disordered superconductors demonstrate a conductor-insulator transition, accompanied by a pseudo-gap in the single-particle DOS and spatial inhomogeneity of the condensate26,27,28,29,30,31,32. However, in moderately disordered materials28,33,34, disorder can have the opposite effect, enhancing superconductivity and increasing the critical temperature35,36,37,38,39. This effect is related to the multi-fractal structure of electrons40, which increases the local DOS near the Fermi level18,41,42, thereby strengthening superconducting pairing. The potential for disorder to elevate the critical temperature has been explored in the context of cuprate superconductors and unconventional superconductivity in general43,44,45. It has been suggested that dopant atoms can enhance the coupling constants associated with Cooper pairing in these systems. Consequently, disorder can also be regarded as a tool to improve superconducting characteristics46,47,48,49,50,51,52,53. By manipulating the spatial distribution of the disordered potential, it is theoretically possible to achieve significantly higher temperatures at which Cooper pairing and the condensate state occur.

However, the presence of a condensate state alone does not necessarily indicate that the material is superconducting32,38. In disordered materials, the condensate can become spatially fragmented or clustered, with strong enhancement at some spatial points and suppression at others. For a supercurrent to be sustained, these regions of enhanced condensate must be interconnected to maintain long-range phase coherence. Disorder tends to disrupt this coherence, potentially isolating superconducting islands.

As a result, disordered materials exhibit two opposing tendencies: on one hand, disorder can facilitate the formation of the condensate state at higher temperatures, while on the other, the inhomogeneity of this state disrupts phase coherence, suppressing the supercurrent. Such systems have two critical temperatures: the local critical temperature \({T}_{c}^{\Delta }\), at which the condensate emerges, and the global one Tc, at which the supercurrent begins to flow38. A disordered material can have a very high \({T}_{c}^{\Delta }\) while at the same time a very low or zero Tc.

This work introduces a mechanism that enhances the superconducting critical temperature Tc in systems comprising multiple condensate components, one of which is strongly disordered. Such a situation arises naturally in multiband superconductors, for example, in MgB254,55,56, where substantial differences in disorder levels between bands can occur due to variations in orbital composition and impurity hybridization. Another example is provided by hybrid structures consisting of alternating layers of disordered and clean materials.

The enhancement of Tc originates from a synergistic interaction between the condensate components. The disordered condensate is localized and possesses a high intrinsic critical temperature \({T}_{c,1}^{\Delta }\), yet it is incapable of supporting a global supercurrent. In contrast, the second component is weakly disordered and able to sustain a supercurrent, although its critical temperature is significantly lower, \({T}_{c,2}\ll {T}_{c,1}^{\Delta }\). When these components are coupled, global superconductivity emerges with a critical temperature close to that of the disordered condensate, \({T}_{c}\simeq {T}_{c,1}^{\Delta }\). This mechanism—restoring global supercurrent at elevated temperatures via coupling to an auxiliary, nearly clean subsystem—is generic and applicable to a wide range of both natural and artificially engineered systems.

Results

The model

We consider a superconductor with two coupled subsystems, referred to as bands, each with its own charge carriers. This model can describe not only two-band superconductors but also layered materials where the subsystems are spatially separated. The system is modeled using a Hamiltonian with three distinct contributions:

$${{\mathcal{H}}}={H}_{1}+{H}_{2}+{H}_{12},$$
(1)

where Hν is the Hamiltonian for each of the bands ν = 1, 2, and H12 describes the interaction between the bands. Superconducting pairing in the bands is described using the lattice Hubbard model with on-site attraction,

$${H}_{\nu }={\sum}_{i,j,\sigma }{h}_{\nu ij}{c}_{\nu i\sigma }^{{\dagger} }{c}_{\nu j\sigma }-{g}_{\nu }{\sum}_{i}{n}_{\nu i\uparrow }{n}_{\nu i\downarrow },$$
(2)

where cνiσ are electron operators for band ν, lattice site i, and spin σ, \({n}_{\nu i\sigma }={c}_{\nu i\sigma }^{{\dagger} }{c}_{\nu i\sigma }\) is the electron number operator, gν > 0 is the coupling constant of the Cooper pairing, and the single-electron part of the Hamiltonian (2) is defined as

$${h}_{\nu ij}={t}_{\nu ij}+{V}_{\nu i}{\delta }_{ij},$$
(3)

where tνij are the hopping matrix elements that have a non-zero value of tν only between the nearest neighbors, and Vνi describes the disorder potential.

The Hamiltonian of the inter-band coupling reads as

$${H}_{12}=-{g}_{12}{\sum}_{i}\left({c}_{1i\uparrow }^{{\dagger} }{c}_{1i\downarrow }^{{\dagger} }{c}_{2i\uparrow }{c}_{2i\downarrow }+c.c.\right),$$
(4)

where g12 is the inter-band coupling constant. The model accounts for the inter-band transfer of Cooper pairs while neglecting the transfer of single electrons. Such inter-band tunneling is absent in clean systems, where the single-particle Hamiltonian is already diagonal in the basis of band states. However, it may become relevant in disordered materials, where strong disorder potentials facilitate scattering between bands, or in hybrid structures, where electrons can tunnel between distinct layers.

The model described by Eqs. (1)–(4) is not restricted to atomic-scale systems, where the index j denotes atomic positions within a crystal lattice. It can also be interpreted as a discretized version of the conventional BCS model for a continuous superconducting medium. In this context, the discretization lattice must be sufficiently fine to resolve all relevant characteristic length scales, such as the BCS coherence length. For convenience, we set the discretization spacing a = 1 as the unit of length in the system.

In the calculations, we assume the model is two-dimensional. It is directly applicable to layered structures; however, we believe that the results of the mean field calculations for this model can also describe 3D superconductors, for example, thin films.

The disordered potential due to impurities is modeled by the potential Vi = V at impurity sites, and Vi = 0 elsewhere57. The model is characterized by two key parameters: the disorder strength V and the impurity density Nd. We assume that band ν = 1 is strongly disordered. The other band, ν = 2, is weakly disordered so that the Anderson’s theorem is applicable and the critical temperature Tc,2 remains unaffected by disorder. For clarity, we assume that disorder is entirely absent in this band.

Impurity positions can be chosen either completely randomly or with spatial correlations. It is well known that spatial correlations in the disorder potential can significantly influence all superconducting properties, including the critical temperature38. The model above can also be used to describe materials with superstructures if one considers the limit of a fully correlated impurity potential distribution.

We note that the goal of this work is not to identify the optimal impurity configuration that maximizes the superconducting critical temperature. Instead, we analyze three representative scenarios: (i) a system with fully random impurity positions, (ii) a system with spatially correlated impurity positions, and (iii) a system featuring periodic 1D and 2D superstructures as the limiting case of the fully correlated disorder potential. In all cases, the model parameters are chosen such that the critical temperature of the strong band \({T}_{c,1}^{\Delta }\) is enhanced by a factor  ~2 compared to the clean, impurity-free case.

The analysis of the system is performed using the mean-field approximation, where the superconducting condensate is determined self-consistently by solving the Bogoliubov-de Gennes (BdG) equations together with the gap equation and the occupation number equation, as detailed in the Methods section. While the mean-field approximation is not strictly applicable to 2D systems due to the influence of fluctuations, the results are expected to remain valid. The validity of these results is supported by several factors: strong disorder, which effectively increases the system’s dimensionality, and the presence of two coupled band condensates. Furthermore, a direct mean-field solution to the BdG equation for 1D chains revealed that the results are equivalent to the corresponding 3D model when the interchain coupling is not very large58.

For the numerical solution, we use the tunneling constants t1,2 = 1 as the system’s energy unit. The pairing coupling constants are set to g1 = 2, g2 = 1, g12 = 0.25, while the band occupations are chosen as n1 = 0.5 and n2 = 0.8. These values correspond to electronic quasiparticles in both bands, with the Fermi level situated far from the center of the band. Note that these occupations correspond to the clean limit, and their relative values are altered by disorder.

The parameters are selected such that, in the absence of interband coupling, the critical temperature of band ν = 1 is significantly higher than that of band ν = 2, i.e., \({T}_{c,1}^{(0)}\gg {T}_{c,2}\). In the clean band ν = 2, all conventional definitions of the critical temperature coincide, yielding \({T}_{c,2}={T}_{c,2}^{(0)}={T}_{c,2}^{\Delta }\). Following standard terminology, band ν = 1 is referred to as the strong band, while band ν = 2 is designated as the weak band.

Despite relatively large values of the coupling constants, the system remains in the BCS regime. The band Fermi energies are εF,1 ≃ 2.5 for the strong band and εF,2 ≃ 3.5 for the weak band. The corresponding gap values at T = 0 are Δ1 ≃ 0.2 and Δ2 ≃ 0.1, ensuring that for both bands, the gap remains much smaller than the Fermi energy, i.e., Δν ≪ εF,ν. The BCS regime is further confirmed by comparing coherence lengths calculated using the standard BCS expression with the numerical results38. Both methods yield the same results: ξ1 ≃ 4 for the strong band and ξ2 ≃ 6 for the weak band, both of which exceed the average electron spacing l ~ 2 in both bands.

The disorder strength in band ν = 1 is set to V = 4, with a disorder density in the range 0.5≤Nd≤0.75. This level of disorder is sufficient to generate a highly inhomogeneous, fragmented (multifractal) condensate state. The fragmentation of the condensate disrupts its long-range phase coherence, as evidenced by the vanishing superfluid stiffness59,60, leading to a suppression of the actual superconducting critical temperature, Tc,1.

For the chosen parameters, the first band approaches the BCS-BEC crossover with kF1ξ1 ≃ 3. This regime exhibits both the much higher local critical temperature at which the condensate state appears and spatial localization, which hinders global superconductivity. In contrast, the parameters of the second band correspond to kF2ξ2 ≃ 6, which approaches the conventional BCS regime. The choice of parameters ensures the desired system design, where one band develops the localized condensate state with a high local Tc, while the other band behaves like a weakly coupled BCS superconductor. We note in passing that the BCS-BEC crossover regime is relevant to a great number of materials, including cuprates and iron-based superconductors, as well as multiband superconductors with one of the bands approaching the Lifshitz transition61,62,63.

We emphasize that the specific choice of parameters is intended solely for illustrative purposes, as the qualitative results do not depend on their exact numerical values. The conclusions of this work remain valid as long as the key elements of the model are preserved: it comprises two coupled components, with the condensate in the strong band being localized or lacking long-range coherence due to a strong disorder potential — the same disorder that enhances its local critical temperature \({T}_{c,1}^{\Delta }\) — while the weak band remains in the clean limit with a much lower intrinsic critical temperature.

The calculations are performed on a finite system of size N × N with N = 32, which is sufficiently large relative to the characteristic length scales of the model for the chosen parameters. To further minimize finite-size effects, periodic boundary conditions are imposed.

Random disorder

First, we investigate scenario (i), where impurity positions are completely random and uncorrelated. The impurity density is set to Nd = 0.5. The results of these calculations are summarized in Fig. 1, where panel (a) displays the temperature dependence of the sample-averaged gaps:

$${\bar{\Delta }}_{\nu }=\frac{1}{{N}^{2}}{\sum}_{i}| {\Delta }_{\nu i}| .$$
(5)

This panel compares the results for the cases where the strong band is disordered (filled circles) and clean (empty circles). Disorder leads to a substantial increase in the critical temperature \({T}_{c}^{\Delta }\), at which the condensate emerges. It is related to the fact that disorder reduces the coherence length, thereby enhancing the local density of single-particle states (DOS) and effectively increasing the superconducting pairing strength. The coherence length decreases significantly to ξ1 ~ 2, driving superconductivity in that band toward the BCS-BEC crossover regime with an enhanced maximal gap \({\Delta }_{1}^{\,{\mbox{max}}\,}\simeq 0.5\) and the critical temperature \({T}_{c,1}^{\Delta }\simeq 1.8{T}_{c,1}^{(0)}\). The intrinsic critical temperature of the weak band is much lower \({T}_{c,2}\ll {T}_{c,1}^{\Delta }\). However, inter-band coupling induces a non-zero gap in the weak band at temperatures close to Tc [cf Fig. 1 (a)] due to the proximity effect.

Fig. 1: Uncorrelated disorder.
figure 1

a Temperature dependence of the sample-averaged gap function \({\bar{\Delta }}_{\nu }\) for the strong band ν = 1 (red) and the weak band ν = 2 (blue), showing results for the cases when strong band is disordered (filled circles) and in its clean limit (empty circles). Low temperature \({T}_{c,2}^{(0)}\) is the critical temperature of the second band in the absence of inter-band coupling. Color density plot with the spatial distribution of the gap function Δ1i (b) and Δ2i (c) at temperature \(T=0.9\,{T}_{c}^{(0)}\). d Profile of the gap function for strong band ν = 1 (red) and weak band ν = 2 (blue), calculated along the dashed lines shown in b and c with \({\Delta }_{1}^{max}=0.365\) and \({\Delta }_{2}^{max}=0.085\). The minimum values of the gap function in the strong and weak bands are indicated by the red and blue dashed lines, respectively. e Histograms of the absolute value of the gap function for the strong band ν = 1 (red) and the weak band ν = 2 (blue). f Homogeneous contribution \({\Delta }_{2}^{(0)}\) to the gap in weak band ν = 2, compared between numerical calculations (circles) and estimates using Eq. (7). g Superfluid stiffness \({D}_{s,2}^{(0)}\) for weak band ν = 2 calculated for the cases when strong band ν = 1 is disordered (filled circles) and in the clean limit (empty circles). The dashed lines represent estimations \({D}_{s,2}^{(0)}\propto | {\Delta }_{2}^{(0)}{| }^{2}\) (orange) and \({D}_{s,2}^{(0)}\propto | {\bar{\Delta }}_{2}{| }^{2}\) (red).

Full size image

As mentioned in the Introduction, the condensate existence does not necessarily imply that the system can conduct supercurrent. Strong disorder leads to a highly inhomogeneous, fragmented condensate state, where the long-range coherence is destroyed. The inhomogeneity is illustrated by the spatial profiles of the band gaps ∣Δ1i∣ and ∣Δ2i∣, shown in Fig. 1b and c, respectively, that are calculated at the critical temperature \({T}_{c}^{(0)}\) of the clean system.

The profiles in Fig. 1b and c look similar, reflecting the fact that the condensate in the weak band takes place due to the proximity effect. However, the contrast between the maximum and minimum values for the two bands differs significantly. For the strong band in Fig. 1b, the contrast reaches a few orders of magnitude, while for the weak band in Fig. 1c, it is much smaller. This is further seen in the profiles of ∣Δ1i∣ and ∣Δ2i∣ shown in Fig. 1d. Here, the profiles are calculated along the dashed lines in Fig. 1b and c, respectively.

The profiles in Fig. 1d demonstrate that for the weak band, the relative difference between the maximal and minimal gap values is much smaller, because in this case the gap has a notable minimal value. In contrast, the minimal value in the strong disordered band is practically zero. Comparing the profiles for both bands, one concludes that the gap function in the weak band can be represented approximately as

$${\Delta }_{2i}={\Delta }_{2}^{(0)}+\gamma {\Delta }_{1i},$$
(6)

where \({\Delta }_{2}^{(0)}\) is the minimal (homogeneous) contribution, and γ is a proportionality constant. Indeed, Fig. 1(d) shows that, with the exception of the minimal values, the profiles appear rather similar, and the intervals where the gap has its minimal and maximal values coincide for both bands.

The homogeneous part of the gap reveals itself in the statistical distribution P(Δ) of the gap values, shown in Fig. 1e. Notice that for both bands, the distribution has a two-peak structure. It reflects the condensate fragmentation (clustering), where the gap is large inside condensate clusters and small between them. The peaks are much closer one to another for the weak band due to its lower contrast. The homogeneous contribution in the weak band is reflected in the minimal non-zero value of the lower peak. In contrast, for the strong band, this minimal value of the lower peak is practically zero.

The appearance of the homogeneous part \({\Delta }_{2}^{(0)}\) is counterintuitive since the proportionality Δ2i ∝ Δ1i is expected due to the proximity effect. Instead, the gap of the weak band demonstrates the homogeneous part for all considered parameters of the model. It is not a coincidence and can be derived directly from the gap equation (see Methods section). To do this one assumes no back action effects – the condensate in the strong band is affected by the weak band negligibly. With this assumption, one indeed obtains Eq. (6) where the homogeneous contribution is given by

$${\Delta }_{2}^{(0)}=\frac{{g}_{12}}{{g}_{1}}\frac{{g}_{2}{c}_{1}\ln \left({c}_{2}/T\right)}{1-{g}_{2}{c}_{1}\ln \left({c}_{2}/T\right)}{\bar{\Delta }}_{1},$$
(7)

c1,2 are model-dependent constants, and \({\bar{\Delta }}_{1}\) is the average gap of the strong band.

The temperature dependence of the homogeneous contribution \({\Delta }_{2}^{(0)}\) is governed by the logarithmic factor in Eq. (7) and by the average \({\bar{\Delta }}_{1}\). Consequently, \({\Delta }_{2}^{(0)}\) vanishes at \({T}_{c}^{\Delta }\). In the vicinity of \({T}_{c}^{\Delta }\), Eq. (7) provides the asymptotic expression \({\Delta }_{2}^{(0)}\propto {\bar{\Delta }}_{1}/T\). This estimate is plotted in Fig. 1f demonstrating an excellent agreement with the numerical results at temperatures \(T\gtrsim 0.5{T}_{c}^{(0)}\).

The homogeneous contribution to the gap means that the condensate in the weak band is not fragmented as its clusters remain connected. Consequently, the long-range coherence is preserved, and the system remains globally superconductive up to temperatures \(T\to {T}_{c}^{\Delta }\). This conclusion is supported by calculating the superfluid stiffness \({D}_{s,2}^{(0)}\) for this band (for details of the superfluid stiffness calculations at finite temperatures, see Methods section). Its logarithmic temperature dependence is shown in Fig. 1g. The stiffness disappears at \({T}_{c}\simeq {T}_{c}^{\Delta }\). For comparison, empty circles show the results for the clean limit. In the latter case, the stiffness disappears at much lower temperature \({T}_{c}^{(0)}\) [cf Fig. 1a], confirming the rise of the critical temperature in the disordered systems.

Notice, that the calculation of \({D}_{s,2}^{(0)}\) does not account for thermal fluctuations, and, therefore, the results should be interpreted with caution, particularly when a system is strongly disordered. Nevertheless, this quantity provides a useful indication of whether the system is superconductive. Plotting it alongside \(| {\bar{\Delta }}_{2}{| }^{2}\) and \(| {\Delta }_{2}^{(0)}{| }^{2}\) in Fig. 1g reveals a special role of the homogeneous part of the condensate. In a clean superconductor, the superfluid stiffness is proportional to the condensate density. For a strongly disordered system, this relation is more complex59,60. However, Fig. 1g shows that the temperature dependence of \(| {\Delta }_{2}^{(0)}{| }^{2}\) is very close to that of \({D}_{s,2}^{(0)}\), suggesting that \({D}_{s,2}^{(0)}\propto | {\Delta }_{2}^{(0)}{| }^{2}\) in the wide temperature range. In contrast, the temperature dependence of the sample-averaged condensate density \(| {\bar{\Delta }}_{2}{| }^{2}\) yields a much poorer approximation for the superfluid stiffness. One concludes that the homogeneous condensate component is the most important factor contributing to the superfluid stiffness and thus to superconducting current.

Correlated disorder

We now turn to scenario (ii), where impurity positions exhibit spatial correlations. Such correlations are present in many real materials, including various classes of superconductors48,64,65. Previous studies have demonstrated that spatial correlations in the disorder landscape can significantly affect superconducting properties, influencing both the local emergence of the condensate and the global phase coherence. To generate an impurity distribution with spatial correlations, we develop an algorithm based on the approach outlined in earlier works38,66,67,68. The details of this method are provided in the Methods section.

Results for a correlated impurity distribution with density Nd = 0.7 are shown in Fig. 2a1–a4, where panels (a1) and (a2) display the band gaps ∣Δ1i∣ and ∣Δ2i∣, respectively. Comparing these with the distributions for uncorrelated disorder in Fig. 1b, c reveals that the correlated disorder is characterized by significantly larger condensate clusters, separated by similarly large intervals of suppressed condensate.

Fig. 2: Correlated disorder.
figure 2

a1–a4 Superconductor with correlated disorder. Color density plot with the spatial distribution of the gap function Δ1i (a1) and Δ2i (a2) for the strong and the weak band, respectively, at temperature \(T=0.9\,{T}_{c}^{(0)}\). a3 temperature dependence of the sample-averaged gap function \({\bar{\Delta }}_{\nu }\) for the strong band ν = 1 (red) and the weak band ν = 2 (blue), showing results for the cases when strong band is disordered (filled circles) and in its clean limit (empty circles). a4 Homogeneous contribution \({\Delta }_{2}^{(0)}\) to the gap in weak band ν = 2, compared between numerical calculations (circles) and estimates using Eq. (7). b1–b4 Superconductor with 1D superstructure. Color density plot with the spatial distribution of the gap function Δ1i (b1) and Δ2i (b2) at temperature \(T=0.9\,{T}_{c}^{(0)}\). b3 Temperature dependence of the sample-averaged gap function \({\bar{\Delta }}_{\nu }\) for the strong band ν = 1 (red) and the weak band ν = 2 (blue), showing results for the cases when strong band is disordered (filled circles) and in its clean limit (empty circles). b4 Homogeneous contribution \({\Delta }_{2}^{(0)}\) to the gap in weak band ν = 2, compared between numerical calculations (circles) and estimates using Eq. (7). c1–c4 Superconductor with 2D superstructure. Color density plot with the spatial distribution of the gap function Δ1i (c1) and Δ2i (c2) at temperature \(T=0.9\,{T}_{c}^{(0)}\). c3 Temperature dependence of the sample-averaged gap function \({\bar{\Delta }}_{\nu }\) for the strong band ν = 1 (red) and the weak band ν = 2 (blue), showing results for the cases when strong band is disordered (filled circles) and in its clean limit (empty circles). c4 Homogeneous contribution \({\Delta }_{2}^{(0)}\) to the gap in weak band ν = 2, compared between numerical calculations (circles) and estimates using Eq. (7).

Full size image

Figure 2a3 gives the temperature dependence of the average gaps \({\bar{\Delta }}_{1}\) and \({\bar{\Delta }}_{2}\), which determine the critical temperature \({T}_{c}^{\Delta }\). For comparison, Fig. 2a3 plots the temperature dependence for the clean system, which gives \({T}_{c}=1.9{T}_{c}^{(0)}\).

As in the case of uncorrelated disorder, Fig. 2a1 and a2 demonstrate a substantial difference in the contrasts of the band gap profiles, which is related to the same mechanism. The disorder gives rise to the condensate fragmentation in the strong band, while in the weak band the gap acquires the homogeneous contribution \({\Delta }_{2}^{(0)}\). Its temperature dependence, extracted from the numerical solution, is plotted in Fig. 2a4. The figure also shows the estimate by Eq. (7), which demonstrates a very good quantitative agreement with the numerical result.

Superstructures

Finally, we examine the extreme case (iii) of a fully correlated potential that forms superstructures. Of particular interest are periodic superstructures, where the critical temperature of the superconducting transition can be very high due to quantum confinement resonances41,69,70,71. However, unlike earlier works, we focus on superstructures where the condensate state is localized and global superconductivity is suppressed. In this case, the condensate is fragmented, and in some parts of the sample, it is nearly completely absent. We consider two prototype periodic structures: 1D stripes with the density of potential sites Nd = 0.5 and a 2D chessboard with Nd = 0.75, as illustrated in Fig. 2b1–b4 and c1–c4.

Calculated profiles of ∣Δ1i∣ and ∣Δ2i∣ are depicted in Fig. 2b1 and b2 for the stripe patterns, and in Fig. 2c1 and c2 for the chessboard configurations, respectively. In these figures, the stripes and squares correspond to regions of maximum and minimum gap values. As observed in the disordered models discussed above, the gap in the stronger band reveals a significantly larger contrast between its maximum and minimum values, indicating that the condensate density is nearly entirely suppressed between the clusters. The fragmentation gives rise to an increase in the critical temperature, with \({T}_{c}\simeq 1.6{T}_{c}^{(0)}\) for the stripes and \({T}_{c}\simeq 2.2{T}_{c}^{(0)}\) for the chessboard.

As in disordered samples, the gap in the weak band comprises a homogeneous contribution, \({\Delta }_{2}^{(0)}\), which persists between the condensate clusters. Its appearance helps the system condensate to sustain the long-range coherence. The temperature dependence of \({\Delta }_{2}^{(0)}\) is shown in Fig. 2b4 and c4, alongside estimates by Eq. (7). As before, the estimate yields an excellent agreement with the numerical results. The homogeneous component only vanishes when the temperature reaches Tc, confirming that the critical temperature for the superconducting transition Tc coincides with that of the condensate \({T}_{c}^{\Delta }\).

Discussion

The above examples describe a mechanism to achieve high-temperature superconductivity by manipulating the spatial distribution of the external potential. The latter creates spatially fragmented localized condensate islands, with a significantly higher critical temperature. The enhancement of the critical temperature is related to the localization and multi-fractal structure of single-particle states, which enhance the Cooper pairing. On the other side, this gives rise to clustering of the condensate state, suppressing the long-range coherence of the condensate, which in turn hinders the superconductive current. As a result, the system develops a condensate state without actually becoming superconductive.

However, global coherence and supercurrent are restored when the disordered material is coupled to an auxiliary “weak” subsystem, where the condensate state is induced by the tunneling of Cooper pairs (proximity effect). In contrast to the original material, the condensate in the weak subsystem comprises a homogeneous part. It connects the condensate clusters, restoring the long-range coherence of the condensate state and enabling superconducting current through the entire sample. This contribution vanishes only with the condensate itself. Consequently, the superconducting critical temperature coincides with that of the condensate state in the strong subsystem.

This mechanism for the superconductivity enhancement is based, eventually, on the proximity effect. Localized Cooper pairs of the material tunnel to an auxiliary subsystem where they can propagate freely, inducing a supercurrent. However, the proximity effect in this system differs from the standard situation, appearing when a superconductor is connected to a normal conductor. Superconducting correlations in the weak subsystem are expected to decay when moving away from a superconducting cluster. In contrast, here induced superconductivity correlations approach a constant value. This appears similar to the phenomenon of Bose-condensation. Bosonic Cooper pairs, created in the strong subsystem where the condensate forms localized clusters, are injected into the clean subsystem, where the pair ground state is homogeneous.

The considered mechanism also has parallels with the phenomenon of suppressing fluctuations in two-band superconductors. The latter takes place where one of the bands has a high mean-field critical temperature, but the superconductivity is disrupted by large fluctuations. However, even a small coupling to the other band suppresses the fluctuations, keeping the superconducting critical temperature close to its high mean-field value25. In this work, superconductivity is compromised by spatial fluctuations of the condensate, but the latter are mitigated by coupling to the auxiliary clean band.

Mechanisms for suppressing both thermal and spatial fluctuations can work alongside. This may take place when spatial fragmentation of the condensate is accompanied by strong thermal fluctuations. Then coupling to the weak auxiliary subsystem restores the long-range coherence, suppressing the fluctuations as well.

As noted above, direct tunneling between different bands — neglected in the present model — may become relevant in strongly disordered or hybrid systems. However, our conclusions remain valid unless the tunneling is sufficiently strong to fully hybridize the single-particle states of the model. In this limiting case, the resulting behavior depends on the distribution of disorder across the newly hybridized bands. If one hybridized band remains strongly disordered while the other stays relatively clean, the system effectively preserves the key features of the original model, and the main conclusions of this work continue to hold.

We emphasize that the proposed mechanism for enhancing superconductivity is likely applicable to both naturally occurring and artificially engineered or modified materials. It can be effective in heterostructures, hybrid systems, and compounds with layered atomic architectures featuring spatially separated condensate components. Moreover, it is relevant for materials with surface superstructures, where localized surface states play the role of the strong band, while the bulk behaves as the weak one. A notable example is bulk graphite with impurity-induced surface superstructures. According to recent reports72, this system exhibits superconductivity at nearly room temperature. Additionally, the recent direct detection of local carrier pairing exhibiting an unexpectedly large gap value Δ ≈ 2 eV in bismuth-based high-temperature superconductors73 positions these materials as promising options for realizing the mechanism proposed in this study to enhance the critical temperature.

Methods

Bogoliubov-de Gennes equations

The condensate state is described using the mean-field approximation for the model (1)–(4). The corresponding mean-field Hamiltonian for each band ν = 1, 2 reads as

$${{{\mathcal{H}}}}_{\nu }={\sum}_{ij\sigma }{h}_{\nu ij}{c}_{\nu i\sigma }^{{\dagger} }{c}_{\nu j\sigma }+{\sum}_{i}{\Delta }_{\nu i}{c}_{\nu i\uparrow }^{{\dagger} }{c}_{\nu i\downarrow }^{{\dagger} }+h.c.,$$
(8)

where the gap functions Δνi are found self-consistently by solving a system of gap equations

$${\Delta }_{1i}={g}_{1}\langle {c}_{1i\uparrow }{c}_{1i\downarrow }\rangle +{g}_{12}\langle {c}_{2i\uparrow }{c}_{2i\downarrow }\rangle ,$$
(9a)
$${\Delta }_{2i}={g}_{2}\langle {c}_{2i\uparrow }{c}_{2i\downarrow }\rangle +{g}_{12}\langle {c}_{1i\uparrow }{c}_{1i\downarrow }\rangle$$
(9b)

which connects the gaps with the anomalous averages 〈cνi↑cνi↓〉. Notice that the model does not account for the cross-band averages. The Hamiltonians (8) are diagonalized separately by solving Bogoliubov-de Gennes equations (BdG) for eigenfunctions \({u}_{\nu i}^{(n)}\), \({v}_{\nu i}^{(n)}\), and eigenvalues Eνn,

$${\sum}_{j}{\tilde{h}}_{\nu ij}{u}_{\nu j}^{(n)}+{\Delta }_{\nu i}{v}_{\nu i}^{(n)}={E}_{\nu n}{u}_{\nu i}^{(n)},$$
(10a)
$$-{\sum}_{j}{\tilde{h}}_{\nu ij}{v}_{\nu j}^{(n)}+{\Delta }_{\nu i}^{* }{u}_{\nu i}^{(n)}={E}_{\nu n}{v}_{\nu i}^{(n)},$$
(10b)

where matrix elements of the single-particle Hamiltonian are

$${\tilde{h}}_{\nu ij}={h}_{\nu ij}+({U}_{\nu i}-\mu ){\delta }_{ij}.$$
(11)

Here, the chemical potential μ and potential Uνi are found from another set of self-consistency equations,

$${N}_{e}={\sum}_{\nu ,i}{n}_{\nu i},\quad {U}_{\nu i}={V}_{\nu i}-\frac{{g}_{\nu }}{2}{n}_{\nu i},$$
(12)

where \({n}_{\nu i}=\langle {c}_{\nu i\sigma }^{{\dagger} }{c}_{\nu i\sigma }\rangle\) and Vνi is the disorder potential. Using solutions of the BdG equations, the anomalous averages are calculated as

$$\langle {c}_{\nu i\uparrow }{c}_{\nu i\downarrow }\rangle ={\sum}_{n}{u}_{\nu i}^{(n)}{v}_{\nu i}^{(n)* }(1-2{f}_{\nu n}),$$
(13)

and the normal averages are

$${n}_{\nu i}=2{\sum}_{n}\left[{\left\vert {u}_{\nu i}^{(n)}\right\vert }^{2}{f}_{\nu n}+{\left\vert {v}_{\nu i}^{(n)}\right\vert }^{2}(1-{f}_{\nu n})\right].$$
(14)

with fνn = f(Eνn/T) being the Fermi distribution.

Homogeneous contribution to the gap

To obtain an estimation for the homogeneous contribution to Δ2, we use the fact that ∣Δ2∣ < ∣Δ1∣ and neglect the back action of Δ2 on Δ1. This implies that we solve the corresponding gap equation for Δ1 separately, and then use the result to find Δ2 from the remaining equations. The corresponding gap equation for the second band becomes

$${\Delta }_{2i}={g}_{2}{F}_{2i}+\frac{{g}_{12}}{{g}_{1}}{\Delta }_{1i},$$
(15)

where

$${F}_{2i}={\sum}_{n}{u}_{2i}^{(n)}{v}_{2i}^{* (n)}\left(1-2{f}_{2n}\right)$$
(16)

with \({u}_{2i}^{(n)}\) and \({v}_{2i}^{* (n)}\) being eigenvectors of the BdG equation for the second band. We now look for these eigenvectors in the form of the Fourier series,

$$\left(\begin{array}{l}{u}_{2i}\\ {v}_{2i}\end{array}\right)={\sum}_{k}\left(\begin{array}{l}{u}_{{{\bf{k}}}}\\ {v}_{{{\bf{k}}}}\end{array}\right){e}^{i{{{\bf{kr}}}}_{i}},$$
(17)

where the sum over k involves the normalizing factor 1/N2. In the inverse space, the BdG equations read as

$${\xi }_{{{\bf{k}}}}{u}_{{{\bf{k}}}}+{\sum}_{{{{\bf{k}}}}^{{\prime} }}{\Delta }_{2}({{\bf{k}}}-{{{\bf{k}}}}^{{\prime} }){v}_{{{{\bf{k}}}}^{{\prime} }}={E}_{{{\bf{k}}}}{u}_{{{\bf{k}}}},$$
(18a)
$$-{\xi }_{{{\bf{k}}}}{v}_{{{\bf{k}}}}+{\sum}_{{{{\bf{k}}}}^{{\prime} }}{\Delta }_{2}{({{\bf{k}}}-{{{\bf{k}}}}^{{\prime} })}^{* }{u}_{{{{\bf{k}}}}^{{\prime} }}={E}_{{{\bf{k}}}}{v}_{{{\bf{k}}}},$$
(18b)

where

$${\Delta }_{2}({{\bf{k}}})={\sum}_{i}{\Delta }_{2i}\,{e}^{-i{{\bf{k}}}{{{\bf{r}}}}_{i}},$$
(19)

and \({\xi }_{{{\bf{k}}}}=-2{t}_{2}(\cos {k}_{x}a+\cos {k}_{y}a)-{\mu }_{2}\) with a being the lattice constant, t2 the hopping integral, and μ2 the chemical potential. Notice that the sum over k is restricted by the Debye frequency ωD. In this equation, we use the fact that g1 is small and neglect the Hartree potential in Eq. (12).

We now utilize the fact that the characteristic size of condensate clusters or the localization length ξ is much larger than the de Broglie wavelength 1/kF of electrons taking part in the Cooper pairing, with ℏkF being the Fermi momentum. Then, one can neglect the slow spatial dependence of the gap to obtain the homogeneous contribution to Δ2:

$${\Delta }_{2}({{\bf{k}}})\simeq \left({g}_{2}{\bar{F}}_{2}+\frac{{g}_{12}}{{g}_{1}}{\bar{\Delta }}_{1}\right)\delta ({{\bf{k}}}),$$
(20)

where δ(k) is the Dirac function, \({\bar{\Delta }}_{1}\) is given in Eq. (5), and

$${\bar{F}}_{2}=\frac{1}{{N}^{2}}{\sum}_{i}{F}_{2i}.$$
(21)

Note that this approximation is eventually equivalent to the Anderson approximation, where the quasiparticle states are looked after in the form of ui = uϕi and vi = vϕi, where ϕi are single-electron states. For the clean system, single-electron states are simply plane waves \({\phi }_{i}\propto \exp (i{{\bf{k}}}{{\bf{r}}})\).

When one inserts this approximation for the gap into the BdG equation (18), it can then be solved explicitly, yielding a well-known result,

$${v}_{{{\bf{k}}}}^{2}=\frac{1}{2}\left(1-\frac{{\xi }_{{{\bf{k}}}}}{{E}_{{{\bf{k}}}}}\right),\quad {u}_{{{\bf{k}}}}^{2}=\frac{1}{2}\left(1+\frac{{\xi }_{{{\bf{k}}}}}{{E}_{{{\bf{k}}}}}\right),$$
(22)

with the quasi-particle energy being

$${E}_{{{\bf{k}}}}=\sqrt{{\xi }_{{{\bf{k}}}}^{2}+{\left({g}_{2}{\bar{F}}_{2}+\frac{{g}_{12}}{{g}_{1}}{\bar{\Delta }}_{1}\right)}^{2}}.$$
(23)

We now substitute this result into the gap equation (15), which yields the gap in the form of the two contributions,

$${\Delta }_{2i}={\Delta }_{2}^{(0)}+\frac{{g}_{12}}{{g}_{1}}{\Delta }_{1i},$$
(24)

where the homogeneous part is given by

$${\Delta }_{2}^{(0)}={g}_{2}{\sum}_{{{\bf{k}}}}\left({g}_{2}{\bar{F}}_{2}+\frac{{g}_{12}}{{g}_{1}}{\bar{\Delta }}_{1}\right)\frac{1-2{f}_{{{\bf{k}}}}}{2{E}_{{{\bf{k}}}}}.$$
(25)

We note that Eq. (24) reproduces our numerical result in Eq. (6) that the gap splits into a homogeneous part \({\Delta }_{2}^{(0)}\) and the one that is proportional to Δ1i.

To find the value of \({\Delta }_{2}^{(0)}\) explicitly, we note that following Eq. (15), this gap contribution is related to the anomalous correlation function as \({\Delta }_{2}^{(0)}={g}_{2}{F}_{2i}\), which is thus constant being equal to its sample average \({\bar{F}}_{2}\). Then, using the substitution \({g}_{2}{\bar{F}}_{2}={\Delta }_{2}^{(0)}\) in Eq. (25) one obtains the homogeneous part of the gap in the form

$${\Delta }_{2}^{(0)}=\frac{{g}_{12}}{{g}_{1}}\frac{{g}_{2}\alpha }{1-{g}_{2}\alpha }\,{\bar{\Delta }}_{1},\quad \alpha ={\sum}_{{{\bf{k}}}}\frac{\left(1-2{f}_{{{\bf{k}}}}\right)}{2{E}_{{{\bf{k}}}}}.$$
(26)

We note that this is still the equation for \({\Delta }_{2}^{(0)}\) as Ek is also the function of this quantity. However, \({\Delta }_{2}^{(0)}\) can be easily estimated when T ~ Tc, and the gap in Eq. (23) can be neglected. In this case, one obtains

$$\alpha \simeq {c}_{1}\ln \left(\frac{{c}_{2}}{T}\right),$$
(27)

where c1 ∝ N2(0) and c2 ∝ ℏωD, with N2(0) being the number of states at the Fermi level. The temperature dependence given by Eqs. (26) and (27) is shown in Fig. 1(f).

Spatially correlated impurity distribution

To describe impurity distribution that is spatially correlated, forming a granulated (clusterized) disorder potential, we employ the model of spatially correlated random functions defined on the lattice as38

$${Z}_{i}=\frac{1}{{N}^{2}}{\sum }_{{j}_{x},{j}_{y}=1}^{N/2}{q}_{j}^{-\beta /2}\cos \left({{{\bf{q}}}}_{j}{{{\bf{r}}}}_{i}+{\phi }_{j}\right),$$
(28)

where ri is a point on the lattice, qj = (2πjx/N, 2πjy/N) is a discrete inverse space vector, jx,y = 1, …, N, qj = ∣qj∣, ϕj are random phases uniformly distributed in the interval \(\left[0,2\pi \right)\), and β is a constant that defines asymptotic dependence of the spatial correlations.

To generate a random impurity distribution used to obtain Fig. 2a1 and a2 we use random Zi from Eq. (28) to construct another random function as

$${Y}_{i}={\left[1-\frac{| {Z}_{i}| }{\max \left(| {Z}_{i}| \right)}\right]}^{2},$$
(29)

which lies in the interval [0, 1]. Next we obtain the normalized quantity

$${X}_{i}=V\frac{{Y}_{i}}{\max \left(| {Y}_{i}| \right)},$$
(30)

with V being the disorder strength, and use it to construct the impurity distribution of the potential as

$${V}_{i}=\left\{\begin{array}{l}V,\quad {X}_{i}\ge \epsilon ;\quad \\ 0,\quad {X}_{i} < \epsilon .\quad \end{array}\right.$$
(31)

Parameter ϵ ∈ [0, V] is the threshold that controls the size of superconducting clusters. If ϵ = V, the potential is absent, and if ϵ = 0 then Vi = V for all i. In the calculations for Fig. 2a1 and a2 we set β = − 4 and ϵ = 0.5V.

Superfluid stiffness

The superfluid stiffness is obtained using a standard expression59,60

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

Here, the kinetic energy along the direction x is given as

$$\langle -{K}_{x}\rangle =\frac{4t}{N}{\sum}_{i,n}\left[{v}_{i}^{(n)}{v}_{i+{e}_{{{\bf{x}}}}}^{(n)}\left(1-{f}_{n}\right)+{u}_{i}^{(n)}{u}_{i+{e}_{{{\bf{x}}}}}^{(n)}{f}_{n}\right].$$
(33)

The transverse current-current correlator is obtained using the expression

$${\Lambda }_{xx}\left({q}_{x}=0,{q}_{y}\to 0,\omega =0\right)=2{t}^{2}$$
(34)
$$\begin{array}{l}\times \sum\limits_{m,n}\sum\limits_{i,j}\left[\left({A}_{ij}^{(m,n)}+{C}_{ij}^{(m,n)}\right)\frac{{f}_{n}-{f}_{m}}{{E}_{n}-{E}_{m}}\right.\\ +\left.\left({B}_{ij}^{(m,n)}+{D}_{ij}^{(m,n)}\right)\frac{{f}_{n}-{f}_{m}-1}{{E}_{n}+{E}_{m}}\right],\end{array}$$
(35)

with

$$\begin{array}{rcl}{A}_{ij}^{(m,n)}&=&{u}_{i}^{(n)}{u}_{j+{e}_{{{\bf{x}}}}}^{(n)}{u}_{j}^{(m)}{u}_{i+{e}_{{{\bf{x}}}}}^{(m)}-{u}_{j+{e}_{{{\bf{x}}}}}^{(n)}{v}_{i+{e}_{{{\bf{x}}}}}^{(n)}{u}_{j}^{(m)}{v}_{i}^{(m)}\\ &&-{u}_{i+{e}_{{{\bf{x}}}}}^{(n)}{u}_{j+{e}_{{{\bf{x}}}}}^{(n)}{u}_{j}^{(m)}{u}_{i}^{(m)}+{u}_{j+{e}_{{{\bf{x}}}}}^{(n)}{v}_{i}^{(n)}{u}_{j}^{(m)}{v}_{i+{e}_{{{\bf{x}}}}}^{(m)}\\ &&-{u}_{i}^{(n)}{u}_{j}^{(n)}{u}_{j+{e}_{{{\bf{x}}}}}^{(m)}{u}_{i+{e}_{{{\bf{x}}}}}^{(m)}+{u}_{j}^{(n)}{v}_{i+{e}_{{{\bf{x}}}}}^{(n)}{u}_{j+{e}_{{{\bf{x}}}}}^{(m)}{v}_{i}^{(m)}\\ &&+{u}_{i+{e}_{{{\bf{x}}}}}^{(n)}{u}_{j}^{(n)}{u}_{j+{e}_{{{\bf{x}}}}}^{(m)}{u}_{i}^{(m)}-{u}_{j}^{(n)}{v}_{i}^{(n)}{u}_{j+{e}_{{{\bf{x}}}}}^{(m)}{v}_{i+{e}_{{{\bf{x}}}}}^{(m)},\end{array}$$
$$\begin{array}{rcl}{B}_{ij}^{(m,n)}&=&{u}_{i}^{(n)}{u}_{j+{e}_{{{\bf{x}}}}}^{(n)}{v}_{j}^{(m)}{v}_{i+{e}_{{{\bf{x}}}}}^{(m)}+{u}_{j+{e}_{{{\bf{x}}}}}^{(n)}{v}_{i+{e}_{{{\bf{x}}}}}^{(n)}{v}_{j}^{(m)}{u}_{i}^{(m)}\\ &&-{u}_{i+{e}_{{{\bf{x}}}}}^{(n)}{u}_{j+{e}_{{{\bf{x}}}}}^{(n)}{v}_{j}^{(m)}{v}_{i}^{(m)}-{u}_{j+{e}_{{{\bf{x}}}}}^{(n)}{v}_{i}^{(n)}{v}_{j}^{(m)}{u}_{i+{e}_{{{\bf{x}}}}}^{(m)}\\ &&-{u}_{i}^{(n)}{u}_{j}^{(n)}{v}_{j+{e}_{{{\bf{x}}}}}^{(m)}{v}_{i+{e}_{{{\bf{x}}}}}^{(m)}-{u}_{j}^{(n)}{v}_{i+{e}_{{{\bf{x}}}}}^{(n)}{v}_{j+{e}_{{{\bf{x}}}}}^{(m)}{u}_{i}^{(m)}\\ &&+{u}_{i+{e}_{{{\bf{x}}}}}^{(n)}{u}_{j}^{(n)}{v}_{j+{e}_{{{\bf{x}}}}}^{(m)}{v}_{i}^{(m)}+{u}_{j}^{(n)}{v}_{i}^{(n)}{v}_{j+{e}_{{{\bf{x}}}}}^{(m)}{u}_{i+{e}_{{{\bf{x}}}}}^{(m)},\end{array}$$

where \({u}_{i}^{(n)}\) and \({v}_{i}^{(n)}\) are solutions of the BdG equations, En is the quasi-particle energy, and fn = f(En/T) is the Fermi distribution. In the summation we take only the solutions with En > 0. Finally, tensors \({C}_{ij}^{(m,n)}\) and \({D}_{ij}^{(m,n)}\) are determined, respectively, from \({A}_{ij}^{(m,n)}\) and \({B}_{ij}^{(m,n)}\) via the replacement \({u}_{i}^{(n)}\iff {v}_{i}^{(n)}\).