The superconducting radio-frequency (SRF) resonant cavities in particle accelerators enable high accelerating gradients with low power consumption. The best Nb cavities can have high quality factors Q ∼ 10¹⁰–10¹¹ and sustain accelerating fields up to 50 MV/m at T = 1.5–2 K and 0.6–2 GHz (Padamsee et al., 2018; Gurevich, 2023). The peak RF fields B ₀ ≃ 200 mT at the equatorial surface of Nb cavities can approach the thermodynamic critical field B _c ≈ 200 mT at which the screening current density flowing at the inner cavity surface is close to the depairing current density J _c ≃ B _c /μ ₀ λ—the maximum DC current density a superconductor can carry in the Meissner state (Tinkham, 2004), where λ is the penetration depth of the magnetic field. Thus, the breakdown fields of the best Nb cavities have nearly reached the DC superheating field B _sh ≃ B _c (Galaiko, 1966; Matricon and Saint-James, 1967; Christiansen, 1969; Chapman, 1995; Catelani and Sethna, 2008; Transtrum et al., 2011; Lin and Gurevich, 2012). The Q factors can be increased by material treatments such as high-temperature annealing followed by low-temperature baking which not only increase Q (B ₀) and the breakdown field but also reduce the deterioration of Q at high fields (Ciovati et al., 2010; Posen et al., 2020). High-temperature treatments combined with the infusion of nitrogen, titanium, or oxygen can produce an anomalous increase of Q (B ₀) with RF field amplitude B ₀ = μ ₀ H ₀ (Ciovati et al., 2016; Grassellino et al., 2017; Dhakal, 2020; Lechner et al., 2021). These advances raise the question about the fundamental limit of the breakdown fields of SRF cavities and the extent to which it can be pushed by surface nanostructuring and impurity management (Gurevich and Kubo, 2017; Gurevich, 2023).

Several ways of increasing the SRF breakdown fields by surface nanostructuring have been proposed. They include depositing high-T _c superconducting multilayers with thin dielectric interlayers (Gurevich, 2006; Gurevich, 2015; Kubo et al., 2014; Liarte et al., 2017; Kubo, 2021) or a dirty overlayer with a higher impurity concentration at the surface (Ngampruetikorn and Sauls, 2019). The DC superheating field of such structures has been evaluated using the London, Ginzburg–Landau (GL), and quasiclassical Usadel and Eilenberger equations in the limit of an infinite GL parameter κ = λ/ξ → ∞ in which the breakdown of the Meissner state at H ₀ = H _sh occurs once the current density at the surface reaches the depairing limit (Gurevich, 2006; Gurevich, 2015; Kubo et al., 2014; Liarte et al., 2017; Kubo, 2021; Ngampruetikorn and Sauls, 2019). Yet, it has been well-established that in a more realistic case of a finite κ, the breakdown of the Meissner state at H = H _sh occurs due to the exponential growth of periodic perturbations of the order parameter and the magnetic field with a wavelength λ c ∼ ( ξ 3 λ ) 1 / 4 along the surface, where ξ is the coherence length (Christiansen, 1969; Chapman, 1995; Transtrum et al., 2011). The effect of such periodic Turing instability (Cross and Hohenberg, 1993) on H _sh can be particularly important for Nb cavities with κ ∼ 1. Addressing the effect of finite κ (which in turn depends on the mean free path) on H _sh in superconductors with nanostructured surfaces is the goal of this work.

We present the results of numerical calculations of a DC superheating field for different superconducting geometries in materials with finite κ and determine optimal surface nanostructures that can withstand the maximum magnetic field in the vortex-free Meissner state. In particular, we consider a bulk superconductor with a thin impurity diffusion layer, a clean superconducting overlayer separated by an insulating layer from the bulk (e.g., Nb₃Sn-I-Nb₃Sn), a thin dirty superconducting layer on the top of the same superconductor (e.g., dirty Nb₃Sn-I-clean Nb₃Sn), and a thin high-T _c superconducting layer on the top of a low-T _c superconductor (e.g., Nb₃Sn-I-Nb). We calculate H _sh and determine an optimal layer thickness for each geometry by numerically solving the GL equations, taking into account both the non-linear screening of the applied magnetic field and the periodic instability of the Meissner state in a nanostructured superconductor.

The paper is organized as follows. The GL equations and methods of numerical detection of H _sh and the wavelength λ _c of a critical perturbation causing the instability of the Meissner state are presented in Section 2. The results of numerical calculations of H _sh for impurity diffusion layers and various S-I-S structures are given in Section 3 and Section 4, respectively. Section 5 contains discussion of the results, and Section 6 provides the conclusion with a summary. Computational details are given in Section Method and other technical details are given in Supplementary Appendices A, B.

We first consider a semi-infinite uniform superconductor in a magnetic field H ₀ applied along the z-axis, parallel to the planar surface. In this case, the induced supercurrents flow in the xy plane and GL equations for the complex order parameter ψ = Δe ^iφ , and two components of the vector potential A _x and A _y can be reduced to two coupled partial differential equations for the amplitude Δ(x, y, t) and the z-component of the magnetic field H (x, y, t). As shown in Supplementary Appendix A, these equations can be written in the following dimensionless form: f ̇ − f + f 3 − ∇ 2 f + κ 2 f 3 ∂ x h 2 + ∂ y h 2 = 0 , (1) ∇ ⋅ ∇ h f 2 = h κ 2 . (2) Here, f (x, y) = Δ(x, y)/Δ₀, Δ₀(T) is the equilibrium order parameter in the bulk, h ( x , y ) = H ( x , y ) / 2 H c , and all lengths are in units of the coherence length ξ and κ = λ/ξ. Despite the presence of the time derivative f ̇ in Eqs (1, 2), they are essentially the quasi-static GL equations, but not the true time-dependent Ginzburg–Landau (TDGL) equations (Watts-Tobin et al., 1981; Sheikhzada and Gurevich, 2020) which describe a non-equilibrium superconductor at T ≈ T _c . Here, f ̇ is added just to detect the instability of the Meissner state in numerical simulations upon slow ramping the applied magnetic field. This procedure allows us to find the field H ₀ = H _sh above which the GL equations no longer have stationary solutions. Another way of numerical calculation of H _sh is based on finding the applied field at which the linearized Eqs (1, 2) have zero eigenmode with f ̇ = 0 (Christiansen, 1969; Chapman, 1995; Transtrum et al., 2011; Liarte et al., 2017), as summarized in Supplementary Appendix B. It turns out that the direct solving Eqs (1, 2) with the ad hoc term f ̇ is much faster than solving the eigenmode problem. For slow magnetic ramp rates H ̇ 0 used in our simulations, the resulting H _sh calculated by these two methods only differ by ≃ 1%, as it is shown in the next Sections. Equations (1, 2) were solved with the following boundary conditions: h 0 , y = h 0 t , h x , 0 = h x , L y , f x , 0 = f x , L y , f L x , y = 1 , h L x , y = 0 , (3) where h 0 = H 0 / 2 H c and the lengths L _x and L _y of the simulation box L _x × L _y were chosen to be ≃ (50–150)ξ depending on κ. The details of the numerical calculations are given in the Supplementary Method.

Shown in Figures 1A,B are f (x, y) calculated at κ = 10 and the applied fields H ₀ slightly below and above H _sh . At H ₀ < H _sh , the order parameter f(x) is reduced at the surface by the flowing screening currents. At H ₀ > H _sh , the stationary f(x) becomes unstable with respect to spontaneously growing periodic perturbations δf (x, y, t) along the surface, as shown in Figure 1B. This Turing instability (Cross and Hohenberg, 1993) occurs with respect to a small disturbance δf (x, y, t) ∝ δf(x)e ^iky+Γt , where the increment Γ(H ₀, k) depends on the wave vector k of spatial oscillations of f (x, y) along the surface as shown in Figure 2A. Below the superheating field, Γ(k) is negative so perturbations with all k decay exponentially and the Meissner state is stable. At the superheating field, Γ(k) first vanishes at a critical wave number k = k _c at which Γ(k) is maximum. At H ₀ = H _sh + 0, the increment Γ(k) becomes positive at k = k _c , making the Meissner state unstable with respect to a growing critical perturbation δ f ( x , y , t ) ∝ δ f ( x ) e i k c y + Γ ( k c ) t with the wavelength λ _c = 2π/k _c , while all other perturbations with k ≠ k _c decay exponentially. We calculated H _sh by slowly ramping the applied field and detecting the onset of the exponential growth of f (x, y, t) with time as described in Method. The critical wavelength λ _c was evaluated from the maximum peak in the spatial Fourier transform of δf (0, y), as shown in Figure 2B. This instability is a precursor of the penetration of the vortex structure with the initial period λ c ∼ ( ξ 3 λ ) 1 / 4 smaller than the stationary vortex spacing ∼ λ ξ at H ₀ ≃ H _sh and κ > 1. The aforementioned direct method for the calculation of H _sh is based on the numerical detection of the field threshold above which the stationary Meissner state does not exist. Here, the time scales of the transition to the vortex state at H ₀ > H _sh are irrelevant, provided that H _sh is calculated at low enough magnetic ramp rates at which H _sh is independent of h ̇ 0 . We calculated H _sh at h ̇ 0 ≃ 1 0 − 5 and verified that H _sh is indeed practically independent of h ̇ 0 , which is also consistent with the calculations of H _sh using both the TDGL and full non-equilibrium equations for dirty superconductors (Sheikhzada and Gurevich, 2020).

We then compare some of our numerical results with the known analytical approximations for H _sh and k _c at κ ≫ 1, given as follows (Christiansen, 1969; Chapman, 1995; Transtrum et al., 2011; Liarte et al., 2017). H s h H c ≈ 5 3 + 0.545 κ , (4) λ k c ≈ 0.956 κ 3 / 4 . (5) At κ = 10–20, Eqs (4, 5) give the instability wavelength λ c = 6.57 ( ξ 3 λ ) 1 / 4 ≃ ( 1.17 − 0.21 ) λ and H _sh approximately (23–16)% higher than H _sh = 0.745H _c in the limit of κ → ∞ in which λ _c → 0 and the breakdown of the Meissner state at H = H _sh occurs once the current density at the surface reaches the depairing limit (Christiansen, 1969). Thus, even at κ = 10–20, characteristic of a dirty Nb or a clean stoichiometric Nb₃Sn (Orlando et al., 1979; Posen and Hall, 2017), the periodic instability along the surface occurs on the scale of the order of the field penetration depth, so the self-consistent GL calculation of H _sh is required.

We consider a dirty layer at the surface with a higher impurity concentration, as shown in Figure 3A. In our simulations, such a layer was modeled by a spatially varying coherence length and penetration depths ξ 2 ( x ) / ξ ∞ 2 = λ ∞ 2 / λ 2 ( x ) = 1 − α ⁡ exp ( − x / l d ) , as shown in Figure 3B. Here, ξ _∞ and λ _∞ are the corresponding bulk values far away from the surface, l _d is the thickness of the diffusion layer, and the parameter α < 1 quantifies the reduction of ξ(0) = (1 − α)^1/2 ξ _∞ and the enhancement of λ(0) = (1 − α)^−1/2 λ _∞ at the surface. The ratio ξ 2 ( x ) / ξ ∞ 2 = λ ∞ 2 / λ 2 ( x ) is controlled by the impurity function χ(ℏv _F /2πT _c l(x)) (Werthamer, 1969) with an inhomogeneous mean free path l(x), which is defined in Supplementary Appendix A. The resulting GL equations take the following form: f ̇ = f − f 3 + ∇ ⋅ S γ ∇ f − κ 2 S γ f 3 ∂ x h 2 + ∂ y h 2 , (6) ∇ ⋅ ∇ h S γ f 2 − h κ 2 = 0 , (7) where κ = λ _∞ /ξ _∞ , S γ = ξ 2 ( x ) / ξ ∞ 2 = 1 − α ⁡ exp ( − x / l d ) , and the lengths are in units of ξ _∞ . The boundary conditions are the same as in Eq. (3). Different impurity profiles were investigated by changing α and l _d at κ = 2 and κ = 10, respectively, representing a cleaner and dirtier Nb.

The calculated dependencies of H _sh (l _d ) on the diffusion layer thickness at different α for κ = 2 and κ = 10 are shown in Figures 4A, B, respectively. One can see that H _sh (l _d ) first increases with l _d , reaches a maximum, and then decreases with l _d approaching a lower value of H _sh at l _d ≫ ξ _∞ . At κ = 2, H _sh (l _d ) is maximum at l _d /ξ _∞ = 0.8, 0.9, 1.5 for α = 0.2, 0.5, 0.8. Similarly, at κ = 10, H _sh (l _d ) is maximum at l _d /ξ _∞ = 4, 5, 10. Here, the diffusion layer can increase H _sh by ≈ 9 % at κ = 2 and by ≈ 14 % at κ = 10 as compared to a superconductor with an ideal surface. A qualitatively similar non-monotonic dependence of H _sh on l _d was also obtained by solving the quasiclassical Eilenberger equations in the entire temperature range 0 < T < T _c (Ngampruetikorn and Sauls, 2019). The maximum in H _sh (d) results from a current counterflow induced in the dirty surface layer by a cleaner substrate with a smaller λ _∞ < λ(0) (Kubo et al., 2014; Gurevich, 2015), the magnitude of the peak increases as the diffusion layer gets dirtier. The curves H _sh (l _d ) cross over at larger l _d for which H _sh (∞) is determined by the surface GL parameter κ(0) = κ _∞ /(1 − α). As a result, H _sh (∞) decreases as the material gets dirtier, in agreement with Eq. (4).

The direct numerical calculation of H _sh involves detecting the instability of the Meissner state with respect to an infinitesimal perturbation δ f ( x , y ) = δ f ( x ) e i k c y + Γ t with a finite wave number k _c and an increment Γ(H ₀) changing the sign from negative at H ₀ < H _sh to positive at H ₀ > H _sh . Figure 5 shows the fast Fourier transform of a snapshot of δf (0, y, t) along y calculated at H ₀ = H _sh + 0 at α = 0.5, κ = 10, and different ratios l _d /ξ _∞ . One can see that δf (0, y) has several harmonics even at the slow field ramp rate h ̇ 0 = 5 ⋅ 1 0 − 5 used in our simulations. Such multi-mode temporal oscillations of δf (0, y) can result from a nonlinear mode coupling above the Turing instability threshold (Cross and Hohenberg, 1993), as well as a finite size of the computational box. In this case, the critical wave number k _c would correspond to the highest peak in the Fourier spectrum of δf _k (0, t). Yet, Figure 5 reveals two uneven peaks whose heights change differently as the ratio l _d /ξ _∞ is varied. For instance, at l _d /ξ _∞ = 5, the critical wave number k _c is determined by the higher left peak observed in Figure 5, but as l _d /ξ _∞ is increased to 6, the right peak becomes higher than the left one, so k _c changes jumpwise at l _d /ξ _∞ ≈ 5.5. The peak shifts toward higher λ _∞ k values, providing a constant k _c in this range of l _d /ξ _∞ . The switching of k _c between two values as l _d /ξ _∞ is increased can be seen in “Fourier Transform. mp4” in Supplementary Video S1. To see the extent to which this ambiguity in k _c may affect H _sh , we have also calculated k _c and H _sh from the sign change of the second variation of the free energy δ ² F caused by small perturbations of δf (x, y) and δh (x, y). In this method (Christiansen, 1969; Chapman, 1995; Transtrum et al., 2011), H _sh is determined by the conditions δ ² F (k _c , H _sh ) = 0 and ∂δ ² F/∂k _c = 0. Figures 6A, B show λ _∞ k _c as a function of l _d /ξ _∞ at κ = 2 and κ = 10 and different α computed from the second variation δ ² F, as described in Supplementary Appendix B. One can see that the peaks in k _c shown in Figure 5 are in the range of λk _c ≈ 5–7 qualitatively matching λk _c (l _d ) ≈ 5 at α = 0.5, as shown in Figure 6. Yet, the H _sh (l _d /ξ _∞ ) curves calculated by these two methods at α = 0.5 turned out to be very similar (the difference is approximately 1%), as shown by the dashed lines in Figure 4.

Using the direct simulation method outlined in Section 2, we have calculated H _sh (d) and the critical wave number k _c (d) for various S-I-S structures: the S layer of thickness d is separated by an insulator from the S substrate of the same material, a dirty S layer on the top of a cleaner superconductor (e.g., dirty Nb-I-clean Nb), and a thin high T _c overlayer on the top of a low T _c superconductor (e.g., Nb₃Sn-I-Nb). Here, the I layer is assumed to be thick enough to fully suppress the Josephson coupling between the S overlayer and the bulk, but thinner than the S overlayer. The screening of the applied field in an S-I-S multilayer is shown in Figure 7.

The GL calculations of the DC superheating field at T ≈ T _c self-consistently take into account the essential non-linear field screening and the periodic instability of the Meissner state in the entire range of the GL parameters which can be tuned by the impurities. This approach shows that the thicknesses of the impurity diffusion layer or S-I-S layers can be optimized to increase H _sh above the superheating fields of individual components. For instance, optimizing the diffusion length can enhance H _sh by ≃ 5–20% at κ = 10 and by ≃ 2–9% at κ = 2. An optimized dirty Nb₃Sn overlayer deposited onto the Nb₃Sn field by ≃ 10% as compared to H _sh of a clean Nb₃Sn. This effect manifests itself in a non-monotonic dependence of H _sh on the dirty layer thickness due to the current counterflow induced at the surface by a superconducting substrate with a shorter field penetration depth. Such behavior of H _sh (d) is consistent with the previous calculations of H _sh based on the London (Gurevich, 2006; Gurevich, 2015; Kubo et al., 2014) or Usadel (Kubo, 2021) and Eilenberger (Ngampruetikorn and Sauls, 2019) theories at κ → ∞. To see the extent to which the London model is consistent with the GL results, we consider H _sh (d) calculated for an S-I-S multilayer in the London limit (Gurevich, 2015). H s h d = cosh d / λ 1 + λ 2 / λ 1 sinh d / λ 1 H sh 2 , d < d m (15) H s h d = λ 1 + λ 2 ⁡ tanh d / λ 1 λ 1 ⁡ tanh d / λ 1 + λ 2 H sh 1 , d > d m (16) d m = λ 1 ⁡ ln λ 1 λ 1 + λ 2 H sh 1 H sh 2 + H sh 1 2 H sh 2 2 + 1 − λ 2 2 λ 1 2 , (17) where H _sh1 and H _sh2 are the bulk superheating fields of the overlayer and the substrate, respectively. Equation (15) describes H _sh (d) of S-I-S structures with thin overlayers (d < d _m ), where the Meissner state first breaks down at the surface of the substrate. Here, the high-H _c overlayer partly screens the substrate, allowing it to stay in the Meissner state at a higher applied field H ₀ = H _sh (d) than the bare substrate. If d > d _m , the Meissner state first breaks down at the outer surface of the overlayer so that H _sh (d) → H _sh1 at d ≫ λ ₁. The maximum H _sh (d _m ) is given by: H s h d m = H sh 1 2 + 1 − λ 2 2 λ 1 2 H sh 2 2 1 / 2 . (18) The maximum H s h ( d m ) = ( H sh 1 2 + H sh 2 2 ) 1 / 2 in an S-I-S multilayer occurs if ( λ 2 / λ 1 ) 2 ≪ 1 .

Equations (15–17) do not take into account the reduction of the superfluid density by current, non-linear field screening and the periodic instability of the Meissner state at a finite κ. The London model does not account for the size effect of reducing H _sh (d) if the overlayer thickness is smaller than the decay length ∼ λ 1 ξ 1 of the critical perturbation, as was discussed in Section 4.1. Indeed, for identical materials of the substrate and overlayer (λ ₁ = λ ₂, H _sh1 = H _sh2), Eqs (15–17) give d _m = 0 and H _sh (d) = H _sh1 independent of d, which is inconsistent with the reduction of H _sh (d) at d ≲ 10ξ ₂, as shown in Figure 8. Yet, the London model captures the non-monotonic thickness dependence H _sh (d) calculated from the GL theory if the overlayer has different properties than those of the substrate, and the input parameters H _sh1 and H _sh2 in Eqs (15–17) are exact bulk superheating fields for given values of κ and H _c , respectively (Gurevich, 2015). For instance, Figure 10 compares the GL and the London H _sh (d) calculated for an Nb₃Sn(dirty)-I-Nb₃Sn structure with a dirty overlayer for which the London model works reasonably well. For the Nb₃Sn-I-Nb multilayers considered in Section 4.3, we observed a rather good agreement between H _sh (d) calculated from the GL theory and Eqs (15–17), as shown in Figure 12. Such surprising accuracy of Eqs (15–17) was also observed by Kubo (2021) in the Usadel simulations of dirty S-I-S multilayers in the entire temperature range of 0 < T < T _c at κ → ∞.

In the GL region T ≈ T _c2, Eqs (15–17) predict a significant change in the temperature dependence of H _sh (T) of the S-I-S multilayer with a higher-T _c overlayer as compared to H _sh2(T) of the bare substrate. If T → T _c2, the penetration depth λ 2 ( T ) ∝ ( T c 2 − T ) − 1 / 2 diverges and H _sh2(T) ∝ T _c2 − T vanishes, while λ ₁ and H _sh1 remain nearly independent of T. This case is characteristic of Nb₃Sn-I-Nb for which T _c1 ≃ 2T _c2, the crossover thickness d _m (T) increases with T and diverges logarithmically at T → T _c2. In turn, H _sh (d, T) obtained by Eq. (15) is limited by the small superheating field of the substrate partially screened by the high-T _c overlayer: H s h T ≃ λ 2 / λ 1 sinh d / λ 1 H sh 2 T ∝ T c 2 − T , (19) Hence, H _sh (T) can be significantly higher than H _sh2(T) ∝ T _c2 − T at T ≈ T _c2, particularly if d > λ ₁. As an illustration, Figure 14 shows H _sh (T) calculated from Eqs (15–17) for different ratios d/λ ₁ and the parameters of Nb₃Sn-I-Nb specified in Section 4.3. One can see both the square root temperature dependence given by Eq. (19) at T ≈ T _c2 and a sharp change in H _sh (T) upon decreasing T as d _m (T) becomes shorter than d and H _sh (T) crosses over to a nearly constant H _sh1(T) of the overlayer. For d/λ ₁ < 2, such a transition in H _sh (T) happens at lower T outside the GL temperature range shown in the figure.

The relation between the static H _sh calculated here and the dynamic superheating field H _sd (T, ω) representing the fundamental field limit of superconductivity breakdown in SRF cavities depends on the rf frequency ω, temperature, and the material purity (Gurevich, 2023). The calculation of H _sd for S-I-S structures generally requires solving complex equations of non-equilibrium superconductivity, which in some cases can be reduced to TDGL equations at T ≈ T _c (Watts-Tobin et al., 1981). The dynamic superheating field of an alloyed superconductor with an ideal surface at T ≈ T _c was calculated from the microscopic theory (Sheikhzada and Gurevich, 2020), where it was shown that H _sd (T) approaches the static H _sh (T) at low frequencies ω ≪ ω _c but can be by a factor 2 larger than H _sh at ω ≫ ω _c . Here, the crossover frequency ω c ∼ min ( τ ϵ − 1 , τ Δ − 1 ) is set by the inelastic electron–phonon scattering time τ _ϵ (T) ∝ T ⁻³ and the TDGL relaxation time of the order parameter τ _Δ = πℏ/8k _B (T _c − T) (Gurevich, 2023). For Nb and Nb₃Sn at T ≈ T c N b = 9.2 K, both τ _ϵ (9K) ∼ 10^–11 s and τ _Δ ∼ 10^–11 s at T _c − T = 0.2 K are much shorter than the rf period at 1 GHz. In this case, the superconducting and quasiparticle screening currents follow practically instantaneously the driving rf field, and the quasistatic H _sh (T) considered here is applicable. The dynamic superheating field at lower temperatures T ≃ 2 K at GHz frequencies has not yet been calculated from a microscopic theory.

In this work, H _sh was calculated for S-I-S structures with ideal surfaces and interfaces without topographical and material defects or weakly coupled grain boundaries in the overlayer and the substrate. Topographical and other surface defects can locally reduce the field onset of the dissipative penetration of vortices and reduce the global H _sh , as was shown by TDGL simulations (Vodolazov, 2000; Pack et al., 2020; Wang et al., 2022). Likewise, H _sh (T) can be reduced by weakly coupled grain boundaries causing premature proliferation of mixed Abrikosov–Josephson vortices or phase slips (Sheikhzada and Gurevich, 2017). The I interlayer in S-I-S coating can mitigate these detrimental effects by the following: 1. increasing the cavity breakdown field by thin high-H _c overlayers and 2. confining vortices penetrating at surface defects in a thin overlayer and blocking flux penetration in the cavity wall, where it can trigger thermo-magnetic avalanches, causing global superconductivity breakdown (Gurevich, 2015; Gurevich, 2023). The S-I-S coating can provide these two goals if the overlayer thickness does not exceed λ ₁ (Gurevich, 2006). In this work, we calculated the upper limit of H _sh and showed how the S-I-S geometry can be optimized to increase H _sh at d ≲ λ ₁.

Our numerical GL calculations of the DC superheating field in superconductors with nanostructured surfaces cover the entire range of 1 < κ < ∞ and account for both the non-linear Meissner screening and the instability with a finite wave number k _c at H ₀ = H _sh . We showed that there are optimum thicknesses of the impurity diffusion layer and the superconducting overlayer which maximize H _sh . These results suggest the possible ways of increasing the breakdown fields by surface nanostructuring and can help understand the ways of optimizing SRF cavities to achieve higher accelerating gradients.