Superheating field in superconductors with nanostructured surfaces

1 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¹¹ 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.

2 GL equations and numerical detection of H _*sh* and k _*c*

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: 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: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.

f ̇ − f + f 3 − ∇ 2 f + κ 2 f 3 ∂ x h 2 + ∂ y h 2 = 0 , (1)

∇ ⋅ ∇ h f 2 = h κ 2 . (2)

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)

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). 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.

H s h H c ≈ 5 3 + 0.545 κ , (4)

λ k c ≈ 0.956 κ 3 / 4 . (5)

3 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 ξ 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: 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.

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)

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.

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₃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.

4.1 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₃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ξ ₂, where ξ ₂ 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/λ ₂ at d < 9ξ ₂ in a thin overlayer to k _*c* ≈ 7.2/λ ₂ at d > 9ξ ₂, corresponding to the instability of the Meissner state in a semi-infinite superconductor. The calculated k _*c* at d > 9ξ ₂ is approximately 10% smaller than k _*c* ≈ 8/λ ₂ given by the asymptotic Eq. (5) at κ ₂ = 17.

4.2 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 ₂ are in units of their respective bulk values of ξ ₂ and Δ₂. In turn, the GL equations in the overlayer are as follows: Equations (8–11) were solved for a dirty Nb₃Sn overlayer on a cleaner Nb₃Sn with a mean free path l = 2 nm, λ 1 ≈ λ 2 ( ξ 2 / l ) 1 / 2 ≈ 135 nm, ξ 1 ≈ ( l ξ 2 ) 1 / 2 ≈ 3 nm, κ ₁ = 45 in the overlayer, and κ ₂ = 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ξ ₂. 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 λ ₂ in the overlayer with a larger λ ₁ (Kubo et al., 2014; Gurevich, 2015). Here, the cusp-like dependence of H _*sh*(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 κ 2 3 / 2 ≈ 8 at d < d _*m* and k c λ 2 = 0.956 κ 1 1 / 4 κ 2 1 / 2 ≈ 10.2.

f ̇ 2 = ∇ 2 f 2 + f 2 − f 2 3 − κ 2 2 f 2 3 ∂ x h 2 2 + ∂ y h 2 2 , (8)

∇ ⋅ ∇ h 2 f 2 2 = h 2 κ 2 2 , (9)

f ̇ 1 = f 1 − f 1 3 + ξ 1 2 ξ 2 2 ∇ 2 f 1 − ξ 1 2 ξ 2 2 κ 1 2 f 1 3 ∂ x h 1 2 + ∂ y h 1 2 , (10)

∇ ⋅ ∇ h 1 f 1 2 = h 1 ξ 2 2 ξ 1 2 κ 1 2 . (11)

4.3 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 ₂ and the field h ₂ in the substrate are described in Eqs (8, 9), and the GL equations for f ₁ and h ₁ in the overlayer are given by where T _*c*1 and T _*c*2 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.

f ̇ 1 = ζ f 1 − f 1 3 + s ∇ 2 f 1 − κ ̃ 2 f 1 3 ∂ x h 1 2 + ∂ y h 1 2 , (12)

∇ ⋅ ∇ h 1 f 1 2 = λ 2 2 h 1 λ 1 2 ζ κ 2 2 , (13)

ζ = 1 − T / T c 1 1 − T / T c 2 , s = ξ 1 2 ξ 2 2 ζ , κ ̃ 2 = ξ 1 2 λ 1 4 ξ 2 2 λ 2 4 κ 2 2 ζ 3 , (14)

We solved the GL equations for a clean Nb₃Sn overlayer on a bulk Nb using κ ₂ = 50/22 and κ ₁ = 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ξ ₂. 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₃Sn overlayer, while H _*sh* at d > d _*m* is determined by the superheating field of Nb₃Sn enhanced at d ≈ d _*m* by the current counterflow caused by the Nb substrate. The corresponding dependence of k _*c*(d) on the overlayer thickness is shown in Figure 13. The jumpwise change of k _*c*(d) reflects the switch from the instability of the Meissner state at the inner surface of the low-T _*c* substrate at d < d _*m* to the instability at the outer surface in the high-T _*c* overlayer at d > d _*m*, which is similar to that shown in Figure 11. For the parameters used in the simulations, such Nb₃Sn-I-Nb structures with d ≈ d _*m* can boost the superheating field up to ∼ 2.2 times higher than the bulk H _*sh*2 of Nb (Gurevich, 2006; Gurevich, 2015) and approximately 5.3% higher than the bulk H _*sh*1 of Nb₃Sn.

5 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₃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). where H _*sh*1 and H _*sh*2 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 _*sh*1 at d ≫ λ ₁. The maximum H _*sh*(d _*m*) is given by: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.

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)

H s h d m = H sh 1 2 + 1 − λ 2 2 λ 1 2 H sh 2 2 1 / 2 . (18)

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 _*sh*1 = H _*sh*2), Eqs (15–17) give d _*m* = 0 and H _*sh*(d) = H _*sh*1 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 _*sh*1 and H _*sh*2 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 _*c*2, 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 _*sh*2(T) of the bare substrate. If T → T _*c*2, the penetration depth λ 2 ( T ) ∝ ( T c 2 − T ) − 1 / 2 diverges and H _*sh*2(T) ∝ T _*c*2 − T vanishes, while λ ₁ and H _*sh*1 remain nearly independent of T. This case is characteristic of Nb₃Sn-I-Nb for which T _*c*1 ≃ 2T _*c*2, the crossover thickness d _*m*(T) increases with T and diverges logarithmically at T → T _*c*2. 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 _*sh*2(T) ∝ T _*c*2 − T at T ≈ T _*c*2, 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 _*c*2 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 _*sh*1(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.

H s h T ≃ λ 2 / λ 1 sinh d / λ 1 H sh 2 T ∝ T c 2 − T , (19)

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 ≲ λ ₁.

6 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 ₀ = 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.

Author contributions

WP performed all numerical simulations, analyzed the results, and wrote the first draft of the manuscript. AG initiated and supervised the project and wrote the revised manuscript. All authors contributed to the article and approved the submitted version.

Funding

This work was supported by DOE under grant DE-SC 100387–020 (ODU) and by the Virginia Military Institute (VMI).

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors, and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/femat.2023.1246016/full#supplementary-material

Supplementary material Supplementary material at the source

Supplementary material Supplementary material at the source