Superheating Field in Superconductors with Nanostructured Superheating Field in Superconductors with Nanostructured Surfaces Surfaces

We report calculations of a DC superheating ﬁ eld H sh in superconductors with nanostructured surfaces. Numerical simulations of the Ginzburg – Landau (GL) equations were performed for a superconductor with an inhomogeneous impurity concentration, a thin superconducting layer on top of another superconductor, and superconductor – insulator – superconductor (S-I-S) multilayers. The superheating ﬁ eld was calculated taking into account the instability of the Meissner state with a non-zero wavelength along the surface, which is essential for the realistic values of the GL parameter κ . Simulations were performed for the material parameters of Nb and Nb 3 Sn at different values of κ and the mean free paths. We show that the impurity concentration pro ﬁ le at the surface and thicknesses of S-I-S multilayers can be optimized to enhance H sh above the bulk superheating ﬁ elds of both Nb and Nb 3 Sn. For example, an S-I-S structure with a 90-nm-thick Nb 3 Sn layer on Nb can boost the superheating ﬁ eld up to ≈ 500 mT, while protecting the superconducting radio-frequency (SRF) cavity from dendritic thermomagnetic avalanches caused by local penetration of vortices.


Introduction
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 -10 11 and sustain accelerating fields up to 50 MV/m at T = 1.5-2K and 0.6-2 GHz (Padamsee et al., 2018;Gurevich, 2023).The peak RF fields B 0 ≃ 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 / μ 0 λ-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 0 ) 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 0 ) with RF field amplitude B 0 = μ 0 H 0 (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 0 = 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 3 Sn-I-Nb 3 Sn), a thin dirty superconducting layer on the top of the same superconductor (e.g., dirty Nb 3 Sn-I-clean Nb 3 Sn), and a thin high-T c superconducting layer on the top of a low-T c superconductor (e.g., Nb 3 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.

GL equations and numerical detection of H sh and k c
We first consider a semi-infinite uniform superconductor in a magnetic field H 0 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: Here, f (x, y) = Δ(x, y)/Δ 0 , Δ 0 (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 y/~ instability of the Meissner state in numerical simulations upon slow ramping the applied magnetic field.This procedure allows us to find the field H 0 = 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: 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 0 slightly below and above H sh .At H 0 < H sh , the order parameter f(x) is reduced at the surface by the flowing screening currents.At H 0 > 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 0 , 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 0 = 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 ikcy+Γ(kc)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 0 ≃ 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 0 > 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 ≃ 10 −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).
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 3 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.

Superconductor with an impurity diffusion layer
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   (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: where , 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 ikcy+Γt with a finite wave number k c and an increment Γ(H 0 ) changing the sign from negative at H 0 < H sh to positive at H 0 > H sh .Figure 5 shows the fast Fourier transform of a snapshot of δf (0, y, t) along y calculated at H 0 = 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 • 10 −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 δ 2 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 δ 2 F (k c , H sh ) = 0 and ∂δ 2 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 δ 2 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.

S-I-S structures
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 3 Sn-I-Nb).Here, the I layer is assumed to be thick enough to fully suppress

S overlayer on the top of the S-substrate
The GL Eqs (1, 2) for the S overlayer separated by the I layer from the substrate made of the same superconductor were solved in both S-domains with the boundary conditions given by Eq. (3) supplemented by the conditions of continuity of h (d + 0, y) = h (d − 0, y), parallel electric field E y (d + 0, y) = E y (d − 0, y), and zero current ∂ y h (d + 0) = ∂ y h (d − 0) = 0 through the I layer.Figure 8 shows that H sh is a function of the thickness of the S overlayer d calculated at κ = 17 representing Nb 3 Sn.Here, a very thin S overlayer reduces H sh (d), which then gradually increases with d, reaching a higher bulk value of H sh at d > 9ξ 2 , where ξ 2 is the coherence length in the S-substrate.The reduction of H sh (d) at small d results from the I layer blocking the perpendicular currents produced by the critical perturbation and reducing its decay length in the bulk from ~ λξ to d.In turn, the critical wave number k c (d) along the surface shown in Figure 9 increases jumpwise from k c ≈ 1.8/λ 2 at d < 9ξ 2 in a thin overlayer to k c ≈ 7.2/λ 2 at d > 9ξ 2 , corresponding to the instability of the Meissner state in a semiinfinite superconductor.The calculated k c at d > 9ξ 2 is approximately 10% smaller than k c ≈ 8/λ 2 given by the asymptotic Eq. ( 5) at κ 2 = 17.

Dirty S overlayer on a cleaner S-substrate
A dirty S overlayer with a higher concentration of non-magnetic impurities on a cleaner S substrate of the same material is considered, assuming that both have the same T c unaffected by non-magnetic impurity scattering (Tinkham, 2004).Superconductivity in the bulk is described by the following GL equations where index 2 corresponds to the substrate parameters in which the lengths and f 2 are in units of their respective bulk values of ξ 2 and Δ 2 .
In turn, the GL equations in the overlayer are as follows: Equations (8-11) were solved for a dirty Nb 3 Sn overlayer on a cleaner Nb 3 Sn with a mean free path l = 2 nm, λ 1 ≈ λ 2 (ξ 2 /l) 1/2 ≈ 135   nm, ξ 1 ≈ (lξ 2 ) 1/2 ≈ 3 nm, κ 1 = 45 in the overlayer, and κ 2 = 17 in the bulk.Figure 10 shows the calculated dependence of H sh on the overlayer thickness which has a maximum at the optimum thickness d m ≈ 9ξ 2 .Such an optimal dirty overlayer can increase H sh by approximately 10% as compared to the bulk H sh .The behavior of H sh (d) at a finite κ turns out to be similar to that calculated using the London and GL theories in the limit of κ → ∞ in which the enhancement of H sh at d ≃ d m results from the current counterflow induced by the substrate with a shorter λ 2 in the overlayer with a larger λ 1 (Kubo et al., 2014;Gurevich, 2015).
Here, the cusp-like dependence of H sh

High-T c superconducting overlayer
Finally, we consider an S-I-S structure comprising a high-T c layer on the top of a lower-T c substrate.The order parameter f 2 and the field h 2 in the substrate are described in Eqs (8, 9), and the GL equations for f 1 and h 1 in the overlayer are given by where T c1 and T c2 are the critical temperatures of the overlayer and the substrate, respectively, and the order parameter and lengths are normalized to the respective parameters of the substrate.Equations (12-14) are supplemented by the boundary conditions given by Eq.
(3) and the conditions of field continuity and zero current through the I layer.We solved the GL equations for a clean Nb 3 Sn overlayer on a bulk Nb using κ 2 = 50/22 and κ 1 = 17 (Orlando et al., 1979;Posen and Hall, 2017).The calculated superheating field H sh (d) shown in Figure 12 has a maximum at d m ≈ 4ξ 2 .This behavior of H sh (d) is similar to that of H sh (d) considered in the previous section: H sh (d) at d < d m is limited by the instability of the Meissner state in the Nb substrate partly screened by the Nb 3 Sn overlayer, while

FIGURE 11
The critical wave number k c (d) calculated for the Nb 3 Sn(dirty)-I-Nb 3 Sn structure by solving the quasi-static GL equations directly as described in Section 2.  (Gurevich, 2006;Gurevich, 2015) and approximately 5.3% higher than the bulk H sh1 of Nb 3 Sn.

Discussion
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 3 Sn overlayer deposited onto the Nb 3 Sn field by ≃ 10% as compared to H sh of a clean Nb 3 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).
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 (λ 1 = λ 2 , H sh1 = H sh2 ),  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ξ 2 , 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 3 Sn(dirty)-I-Nb 3 Sn structure with a dirty overlayer for which the London model works reasonably well.For the Nb 3 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 c2 − T) −1/2 diverges and H sh2 (T) ∝ T c2 − T vanishes, while λ 1 and H sh1 remain nearly independent of T. This case is characteristic of Nb 3 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: Hence, H sh (T) can be significantly higher than H sh2 (T) ∝ T c2 − T at T ≈ T c2 , particularly if d > λ 1 .As an illustration, Figure 14 shows H sh (T) calculated from Eqs (15-17) for different ratios d/λ 1 and the parameters of Nb 3 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/λ 1 < 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 −3 and the TDGL relaxation time of the order parameter τ Δ = πZ/8k B (T c − T) (Gurevich, 2023).For Nb and Nb 3 Sn at T ≈ T Nb c 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 λ 1 (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 ≲ λ 1 .

Conclusion
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 0 = 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.T/Tc2 0.9
as shown in Figure3B.Here, ξ ∞ and λ ∞ are the corresponding

Frontiers
FIGURE 2 (A) Qualitative dependence of the instability increment Γ(H 0 , k) on the wave vector of perturbation k at different applied fields H 0 .(B) Snapshot of δf(y) at x =0, and H 0 = H sh +0 calculated at κ =10.
FIGURE 3 (A) Impurity diffusion layer at the surface shown by the dark gray contrast.(B) Variations of normalized coherence length and penetration depth across a dirty layer with α =0.5.

FIGURE 7
FIGURE 7Superconductor-insulator-superconductor structure.The vertical black line represents the insulating layer, and the red line shows the screened magnetic field profile H(x).

FIGURE 8
FIGURE 8Superheating field H sh (d) calculated for the Nb 3 Sn-I-Nb 3 Sn structure.

FIGURE 9 Frontiers
FIGURE 9Critical wave number k c (d) calculated for the Nb 3 Sn-I-Nb 3 Sn structure by solving the quasistatic GL equations directly as described in Section 2.
(d) is controlled by the instability of the Meissner state in the substrate at d < d m and by the instability of the Meissner state in the overlayer at d > d m , the overlayer partly screening the substrate and allowing it to withstand external fields higher than the bulk H sh .The corresponding critical wave number k c (d) is shown in Figure 11.The jumpwise change of k c (d) reflects the switch from the instability of the Meissner state at the inner surface of the substrate at d < d m to the instability at the outer surface in the overlayer at d > d m .Such a jump in k c also follows Eq. (5) which gives k c λ 2 0.956κ 3/2 2 ≈ 8 at d < d m and k c λ 2 0.956κ 1/4 1 κ 1/2 2 ≈ 10.2.

FIGURE 10
FIGURE 10 Superheating field H sh (d) calculated for the Nb 3 Sn(dirty)-I-Nb 3 Sn structure.The red dashed line shows H sh (d) calculated from the London Eqs (15-17) with H sh1 =0.855H c and H sh2 =0.91H c taken from the asymptotic limits of H sh (d) at d ≫ λ 1 and d =0, respectively.

FIGURE 12 Frontiers
FIGURE 12 Superheating field H sh (d) calculated for the Nb 3 Sn-I-Nb structure.The red dashed line shows H sh (d) calculated from the London Eqs (15-17) with H sh1 =2.28H c and H sh2 =1.08H c taken from the asymptotic limits of H sh (d) at d ≫ λ 1 and d =0, respectively.
) 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 0 = 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 ≫ λ 1 .The maximum H sh (d m ) is given by:

FIGURE 13
FIGURE 13The critical wave number k c (d) calculated for the Nb 3 Sn-I-Nb structure by solving the quasistatic GL equations directly as described in Section 2.

FIGURE 14 Frontiers
FIGURE 14 Temperature dependencies of H sh (T) calculated from Eqs (15-17) for different ratios d/λ 1 and the superconducting parameters of Nb 3 Sn-I-Nb specified in the text.The sharp change in the behavior of H sh (T) upon decreasing T at d/λ 1 =2 occurs 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 smaller d/λ 1 , such a transition in H sh (T) takes place at lower T outside the GL temperature range shown in the figure.Here, the blue line with d/λ 1 =0 represents H sh2 (T) of the bare substrate.