Field-dependent nonlinear surface resistance and its optimization by surface nano-structuring in superconductors

TL;DR

This paper develops a theory for the nonlinear surface resistance of dirty superconductors under strong RF fields, considering impurities and surface nanostructuring, to optimize superconducting RF performance.

Contribution

It introduces a comprehensive model accounting for impurities and surface layers, revealing how nanostructuring can minimize nonlinear RF losses in superconductors.

Findings

Minimum surface resistance occurs at specific RF field amplitudes.

Paramagnetic impurities can lower surface resistance over a range of fields.

Surface nanostructuring can further reduce surface resistance below ideal surfaces.

Abstract

We propose a theory of nonlinear surface resistance of a dirty superconductor in a strong radio-frequency (RF) field, taking into account magnetic and nonmagnetic impurities, finite quasiparticle lifetimes, and a thin proximity-coupled normal layer characteristic of the oxide surface of many materials. The Usadel equations were solved to obtain the quasiparticle density of states (DOS) and the low-frequency surface resistance $R_s$ as functions of the RF field amplitude $H_0$. It is shown that the interplay of the broadening of the DOS peaks and a decrease of a quasiparticle gap caused by the RF currents produces a minimum in $R_s(H_0)$ and an extended rise of the quality factor $Q(H_0)$ with the RF field. Paramagnetic impurities shift the minimum in $R_s(H_0)$ to lower fields and can reduce $R_s(H_0)$ in a wide range of $H_0$. Subgap states in the DOS can give rise to a residual surface resistance while reducing $R_s$ at higher temperatures. A proximity-coupled normal layer at the surface can shift the minimum in $R_s(H_0)$ to either low and high fields and can reduce $R_s$ below that of an ideal surface. The theory shows that the behavior of $R_s(H_0)$ changes as the temperature and the RF frequency are increased, and the field dependence of $Q(H_0)$ can be very sensitive to the materials processing. Our results suggest that the nonlinear RF losses can be minimized by tuning pairbreaking effects at the surface using impurity management or surface nanostructuring.

Equations

$$N(ϵ)=\text{mbox}Re\frac{N_s(ϵ+i\Gamma )}{\sqrt{(ϵ+i\Gamma {)}^2-{\Delta }^2}},ϵ>0,$$

$$\frac{D}{2}{\theta }^{\prime \prime }=s\sin \theta \cos \theta +{\tilde{\omega }}_n\sin \theta -\Delta (x)\cos \theta ,$$

$$s=s_0{e}^{-2x/\lambda }{\sin }^2\omega t+{\Gamma }_p,$$

$$\displaystyle s_{0}=\frac{\Delta}{\pi}\biggl{(}\frac{H_{0}}{H_{c}}\biggr{)}^{2},$$

$$\Delta (x)=2\pi Tg\sum _{{\omega }_n>0}^{\Omega }\sin \theta (x).$$

$${\theta }^{\prime }(-d)=0,\theta (\infty )={\theta }_s,$$

$${\sigma }_n{R}_B{\theta }_{-}^{\prime }=\sin ({\theta }_0-{\theta }_{-}),$$

$${\sigma }_n{\theta }_{-}^{\prime }={\sigma }_s{\theta }_0^{\prime },$$

$$\alpha =\frac{N_n}{N_s}\frac{d}{{\xi }_s},$$

$$\beta =\frac{4e^2}{\hslash }{R}_B{N}_n\Delta d,$$

$$\frac{D}{2}{\theta }^{\prime \prime }=s\sinh \theta \cosh \theta -i\widetilde{ϵ}\sinh \theta +i\Delta (x)\cosh \theta ,$$

$$n(ϵ)=\frac{N(ϵ)}{N_s}=\mathrm{Re}{G}^R(ϵ).$$

$${\Delta }_s=\Delta -\Gamma -\frac{\pi s}{4},T\ll T_c,$$

$$R_s$$

$${\sigma }_1=\frac{{\sigma }_n}{\pi }{\int }_0^{\pi /\omega }dt{\int }_{-\infty }^{\infty }dϵ[f(ϵ,t)-f(ϵ+\omega ,t)]M(ϵ,\omega ,x,t),$$

$$M=\mathrm{Re}{G}^R(ϵ)\mathrm{Re}{G}^R(ϵ+\omega )+\mathrm{Re}{F}^R(ϵ)\mathrm{Re}{F}^R(ϵ+\omega ).$$

$${\sigma }_1=\frac{{\sigma }_n}{\pi }{\int }_0^{\pi /\omega }dt{\int }_{-\infty }^{\infty }\frac{(1-e^{-\omega /T})M[ϵ,\omega ,s(x,t)]dϵ}{(1+e^{-ϵ/T})(e^{ϵ/T}+e^{-\omega /T})}.$$

$$R_{\rm MB}=\frac{\mu_{0}^{2}\omega^{2}\lambda^{3}\Delta}{\rho_{n}T}\ln\biggl{(}\frac{CT}{\omega}\biggr{)}e^{-\Delta/T},$$

$$\displaystyle\tau_{r}=\tau_{1}\biggl{(}\frac{T_{c}}{T}\biggr{)}^{1/2}e^{\Delta/T},\hskip 28.45274pt\tau_{s}=\tau_{2}\biggl{(}\frac{T_{c}}{T}\biggr{)}^{7/2},$$

$$\delta f_{ϵ}(x,t)=\frac{h_{ϵ}(x,t)}{{\cosh }^2(ϵ/2T)},$$

$$\cos v_{ϵ}\cosh u_{ϵ}\frac{\partial h_{ϵ}}{\partial t}=D\nabla \cdot ({\cos }^2{v}_{ϵ}\nabla h_{ϵ})+$$

$$\frac{1}{2T}\cos v_{ϵ}\sinh u_{ϵ}\left(\frac{\pi }{4}-\sin v_{ϵ}\cosh u_{ϵ}\right)\frac{\partial s}{\partial t}-I_{ph},$$

$$I_{ph}=2\pi b\cosh \frac{ϵ}{2T}\cos v_{ϵ}\cdot$$

$${\int }_{-\infty }^{\infty }d{ϵ}^{\prime }(h_{ϵ}-h_{{ϵ}^{\prime }})\frac{(ϵ-{ϵ}^{\prime }{)}^2\cos v_{{ϵ}^{\prime }}\cosh (u_{ϵ}-u_{{ϵ}^{\prime }})}{\cosh ({ϵ}^{\prime }/2T)\sinh (|ϵ-{ϵ}^{\prime }|/2T)}.$$

$$\frac{\partial h_{ϵ}}{\partial t}+{\gamma }_{ϵ}(t)h_{ϵ}=F_{ϵ}(t),$$

$${\gamma }_{ϵ}=2\pi b\frac{\cosh (ϵ/2T)}{\cosh u_{ϵ}}{\int }_{-\infty }^{\infty }d{ϵ}^{\prime }\frac{(ϵ-{ϵ}^{\prime }{)}^2\cos v_{{ϵ}^{\prime }}\cosh (u_{ϵ}-u_{{ϵ}^{\prime }})}{\cosh ({ϵ}^{\prime }/2T)\sinh (|ϵ-{ϵ}^{\prime }|/2T)},$$

$$F_{ϵ}(t)=\frac{1}{2T}\tanh u_{ϵ}\left(\frac{\pi }{4}-\sin v_{ϵ}\cosh u_{ϵ}\right)\frac{\partial s}{\partial t}.$$

$$h_{ϵ}=\frac{1}{2T{\gamma }_{ϵ}}\tanh u_{ϵ}\left[\frac{\pi }{4}-\sin v_{ϵ}\cosh u_{ϵ}\right]\frac{\partial s}{\partial t}.$$

$$h_{ϵ}=\frac{1}{2T}{\int }_0^{s(t)}\tanh u_{ϵ}\left[\frac{\pi }{4}-\sin v_{ϵ}\cosh u_{ϵ}\right]ds.$$

$$is\sinh \theta \cosh \theta +ϵ\sinh \theta -{\Delta }_s\cosh \theta =0.$$

$$\displaystyle\sin v=\frac{1}{2s\cosh u}\Bigl{(}-\Delta_{s}+\sqrt{\Delta_{s}^{2}-s^{2}\sinh^{2}2u}\Bigr{)},$$

$$\displaystyle\sinh 2u=\frac{1}{s}\Bigl{[}(q+\epsilon s\Delta_{s})^{1/3}-(q-\epsilon s\Delta_{s})^{1/3}\Bigr{]},$$

$$\displaystyle q=\bigl{[}\epsilon^{2}s^{2}\Delta_{s}^{2}+\frac{1}{27}(\epsilon^{2}-\Delta_{s}^{2}+s^{2})^{3}\bigr{]}^{1/2}.$$

$${ϵ}_g=({\Delta }_s^{2/3}-s^{2/3}{)}^{\frac{3}{2}},$$

$${\Delta }_s={\Delta }_0-\frac{\pi }{4}s.$$

Full text

I Introduction

The physics of electromagnetic response of superconductors, and the fundamental limits dissipation in the Meissner state at low temperatures and frequencies has recently attracted much interest. The issue of ultra-low-dissipation is particularly important for microresonators for quantum computing or radio-frequency (rf) superconducting cavities for particle accelerators  Zmuidzinas ; Devoret ; 2002_Hein ; 2012_Gurevich_review ; Padamsee_book . At temperatures $T$ well below the critical temperature $T_{c}$ and frequencies $\omega$ smaller than the gap frequency $2\Delta/\hbar$, $s$-wave superconductors have very small surface resistance $R_{s}\propto\exp(-\Delta/T)$  1958_Mattis_Bardeen . Indeed, the Nb cavities typically have $R_{s}\sim 10~{}\rm n\Omega$ at $2\,{\rm K}$ and $1\,{\rm GHz}$, which translates into huge quality factors $Q\propto 1/R_{s}\sim 10^{10}$-$10^{11}$. The surface resistance depends on the amplitude $H_{0}$ of the rf magnetic field $H(t)=H_{0}\sin\omega t$ and can be significantly altered by the materials treatments. For instance, $R_{s}$ of electropolished Nb cavities  Padamsee_book at $2\,{\rm K}$ and $1\,{\rm GHz}$ increases with the rf field amplitude, consistent with the reduction of a quasiparticle gap and the superfluid density by the rf pairbreaking currents  parment ; bardeen ; kopnin . This manifests itself in a field-dependent London penetration depth $\lambda(H)$ and the nonlinear Meissner effect  nme1 ; nme2 ; nme3 ; nme4 ; nme5 ; nme6 . Yet the Nb cavities infused with nitrogen  2013_Grassellino , titanium  2013_Dhakal or other impurities  cornell ; jlab can exhibit a striking field-induced reduction of $R_{s}(H_{0})$ by factors of 2-4 as $H_{0}$ increases from [math] to $\lesssim 0.5H_{c}$, where $H_{c}$ is the thermodynamic critical field.

The behavior of $R_{s}(H_{0})$ is determined by multiple mechanisms including interplay of temporal oscillations of the density of states (DOS) and the kinetics of nonequilibrium quasiparticles under the strong rf field. The field dependence of $R_{s}(T,H_{0})$ is also sensitive to the ratios $\omega/T$ and $T/\Delta$, as well as the electronic structure and compositional inhomogeneities at the surface. It was shown 2014_Gurevich_PRL that the well-known effect of broadening of the DOS gap peaks by the pairbreaking current  2003_Anthore_PRL ; 1964_Maki ; 1965_Fulde ; 1969_Maki can result in a pronounced minimum in $R_{s}(H_{0})$, in agreement with experiment  2013_Dhakal ; 2014_Ciovati_APL . Such microwave reduction of the surface resistance  2014_Gurevich_PRL ; 2017_Gurevich_SUST is a manifestation of a general effect by which $R_{s}$ can be reduced by engineering an optimum broadening of the DOS peaks at the surface using pairbreaking mechanisms. These mechanisms can be due to the rf Meissner currents, magnetic impurities, local reduction of the pairing constant, or a proximity-coupled normal layer which models nonstoichiomentry and metallic suboxide layers at the surface  2017_Gurevich_Kubo .

Sparse magnetic impurities in Nb are among realistic materials features  2005_Casalbuoni ; 2011_Proslier which broaden the DOS peaks 1966_Fulde_Maki ; 1969_Maki ; 2006_Balatskii ; 2011_Fominov ; 2012_Kharitonov and can reduce the low-field $R_{s}$ by $\sim 50\%$ despite a small reduction of $T_{c}$  2017_Gurevich_Kubo . The DOS broadening and the appearance of subgap states at quasiparticle energies $|\epsilon|<\Delta$ have been revealed by numerous tunneling experiments  2003_Zasadzinski ; 2013_Dhakal . Such DOS has been commonly described by the phenomenological Dynes model which incorporates a constant quasiparticle lifetime $\hbar/\Gamma$ by the replacement $\epsilon\to\epsilon+i\Gamma$  1978_Dynes ; 1984_Dynes :

[TABLE]

where $N_{s}$ is the density of states in the normal state. Numerous STM experiments have shown that the DOS broadening can be significant, particularly in thin films and bilayers blstm1 ; blstm2 ; stmPb . The subgap states have been attributed to inelastic scattering of quasiparticles on phonons inelast ; kopnin , Coulomb correlations coulomb , anisotropy of the Fermi surface anis , inhomogeneities of the BCS pairing constant larkin , magnetic impurities 2006_Balatskii , spatial correlations in impurity scattering 2006_Balatskii ; meyer , or diffusive surface scattering. arnold

The Dynes model with a constant $\Gamma$ has not been derived from a microscopic theory (see, e.g., Ref. feigel, for an overview of different mechanisms), yet Eq. (1) gives an insight into how the DOS broadening could affect $R_{s}(T)$ at different temperatures. For instance, at $(\omega,\Gamma)\ll T\ll\Delta$, the surface resistance is mostly determined by thermally-activated quasiparticles with $\epsilon\approx\Delta$. In this case the DOS broadening reduces $R_{s}(T)\propto e^{-\Delta/T}\ln(T/\Gamma)$  2017_Gurevich_SUST . However, the effect of subgap states reverses at very low temperatures $T<(\Gamma,\omega)$ for which $R_{s}(T)$ is dominated by quasiparticles with $\epsilon\ll\Delta$ if $\Gamma(\epsilon)$ does not vanish at $\epsilon\to 0$. In this case a finite DOS at the Fermi level $N(0)=N_{s}\Gamma/\Delta$ increases $R_{s}$ as compared to the BCS model, giving rise to a residual surface resistance $R_{i}\propto(\Gamma/\Delta)^{2}$ at $T\to 0$  2017_Gurevich_SUST ; 2017_Gurevich_Kubo . This dependence of $R_{s}$ on $\Gamma$ can be used to minimize $R_{s}(T)$ in a particular temperature region by engineering an optimum DOS at the surface. While ways of affecting $\Gamma$ in the bulk are not well understood, a better way of engineering the optimum DOS is to use the materials treatment and nanostructuring of the surface of superconducting materials which usually have a thin layer of weakened superconductivity. For instance, the Nb surface is covered by a layer of dielectric ${\rm Nb_{2}O_{5}}$ oxide followed by a few nm thick layer of normal metallic suboxides. Other materials such as ${\rm Nb_{3}Sn}$, ${\rm MgB_{2}}$, or iron-based superconductors can exhibit a significant surface nonstoichiometry, which can be modeled by a thin normal (N) layer coupled with the bulk by the proximity effect or a superconducting (S′) layer separated by a thin insulating (I) layer from the S substrate mlag1 ; mlag2 ; mltk1 ; mltk2 ; mlc . Modification of realistic surface structures can cause profound changes in the low-energy DOS which determines the surface resistance. In particular, it has been shown  1989_Golubov ; 1996_54_Belzig ; 2017_Gurevich_Kubo that a proximity-coupled N layer causes a disturbance of DOS at $\epsilon\approx\Delta$ which extends into the S region over distances much greater than the coherence length. As a result, a thin N layer with a moderately transparent NS interface can reduce $R_{s}$ at weak RF fields by $\sim 15\%$ relative to the ideal surface 2017_Gurevich_Kubo . Here both the thickness of N layer or the transparency of the N-S interface can be tuned by the materials processing and heat treatment.

Nonlinear dc screening for an ideal surface  nme1 ; nme2 ; nme3 ; nme4 ; nme5 ; nme6 and proximity-coupled N-S sandwiches, including the superconductivity breakdown under a strong dc magnetic field   1996_53_Belzig ; 1997_Fauchere ; 1989_Belzig ; 2003_Galaktionov have been thoroughly investigated in the literature. Yet the surface resistance, which is rather sensitive to the details of DOS at $|\epsilon|<\Delta$, paribreaking effects and nonequilibrium kinetics of quasiparticles, has been understood to a much lesser extent. Recently we calculated $R_{s}$ for a proximity-coupled N layer under a weak rf field  2017_Gurevich_Kubo but a theory of the field-dependent nonlinear surface resistance in the Meissner state affected by different pairbreaking effects at the surface is lacking. Such a theory should take into account the microwave reduction of $R_{s}$ due to current pairbreaking  2014_Gurevich_PRL ; 2017_Gurevich_SUST , the nonlinear rf response of a proximity-coupled N layer and nonequilibrium kinetics of quasiparticles due to collision with phonons and impurities under strong rf field kopnin ; kramer . In this work we address some of these issues and calculate the field-dependent surface resistance for an imperfect surface using the Usadel equations for dirty superconductors  1970_Usadel ; 1999_Belzig_review ; 2004_Golubov_review . Our results show how the dependence of $R_{s}(H_{0})$ on the rf amplitude is affected by multiple realistic materials features and suggest ways by which $R_{s}(H_{0})$ could be optimized by tuning the concentration of paramagnetic impurities or properties of the N layer. These results may be useful for improving the quality factors of the resonant cavities for particle accelerators and microresonators for quantum information processing and photon detectors.

The paper is organized as follows. In Sec. II, the geometry of the problem, the Usadel equations, and the boundary conditions are formulated. We obtain formulas for a low-frequency nonlinear surface resistance and evaluate a contribution of nonequilibrium effects. In Sec. III, we calculate $R_{s}(H_{0})$ for a superconductor with paramagnetic impurities and show that they can produce a significant minimum in $R_{s}(H_{0})$ as a function of $H_{0}$. In Sec. IV, we investigate the effect of the Dynes parameter $\Gamma$ on $R_{s}$ and show that a finite $\Gamma$ results in a residual resistance at $T\ll T_{c}$ and affects a nonlinear field dependence of $R_{s}(H_{0})$ in the way similar to that of paramagnetic impurities at intermediate temperatures. In Sec. V, we consider the effect of SIS′ surface nanostructuring which can shift the minimum in $R_{s}(H_{0})$ to higher fields. In Sec. VI, we calculate $R_{s}(T,H_{0})$ for a superconductor covered with a thin, proximity-coupled N layer and investigate the field dependence of $R_{s}(H_{0})$ and the position of the minimum of $R_{s}(H_{0})$ as functions of temperature, thickness of the N layer and the contact resistance $R_{B}$. In Sec. VII, we discuss implications of our results for engineering an optimum DOS at the surface to reduce the rf losses in superconductors under strong electromagnetic fields.

II Green’s functions and surface resistance

II.1 Usadel equation

Consider the geometry shown in Fig. 1 which represents a superconductor with an ideal surface or a superconductor S covered with a proximity-coupled $N$ layer or a layer of another superconductor $S^{\prime}$ at $-d\leq x<0$. For a dirty superconductor, the equilibrium normal and anomalous Green’s functions $G=\cos\theta$ and $F=\sin\theta$ satisfy the thermodynamic Usadel equation:

[TABLE]

Here $D$ is the electron diffusivity, the prime denotes differentiation with respect to $x$, $\omega_{n}=\pi T(2n+1)$ are the Matsubara frequencies, $\tilde{\omega}_{n}=\omega_{n}+\Gamma$, and $\Gamma$ accounts for a finite quasiparticle lifetime in the Dynes model. The pairbreaking parameter $s=DQ^{2}/2+\Gamma_{p}$ contains two contribution: the first term comes from the effect of pairbreaking currents and the second term accounts for spin-flip scattering on magnetic impurities. Here ${\bf Q}=\nabla\theta+2\pi{\bf A}/\phi_{0}$ is the gauge invariant phase gradient, A is the vector potential, $\phi_{0}$ is the flux quantum, and $\Gamma_{p}$ is the spin-flip parameter due to magnetic impurities 2006_Balatskii ; 2011_Fominov ; 2012_Kharitonov . For the planar geometry shown in Fig. 1, $(\phi_{0}/2\pi)Q(x,t)=-\mu_{0}\lambda H_{0}e^{-x/\lambda}\sin\omega t$ at $x\geq 0$ and $Q(x,t)=Q(0,t)$ at $x<0$, where $\lambda$ is the bulk London penetration depth. We assume that $H_{0}$ is well below the superheating field hs1 ; hs2 ; hs3 so a weak dependence of $\lambda$ on $H_{0}$ due to the nonlinear Meissner effect nme1 ; nme2 ; nme3 ; nme4 ; nme5 ; nme6 is negligible. We also neglect the field attenuation in the N surface layer with the thickness $d\ll\xi\ll\lambda$. Hence,

[TABLE]

where $H_{c}=(N_{s}/\mu_{0})^{1/2}\Delta$ is the thermodynamic critical field, $\Delta$ is the pair potential at $T=0$, and $\Gamma_{p}$ is assumed uniform. Here $\Delta(x,T,s)$ satisfies the BCS self-consistency equation,

[TABLE]

Here $g$ is the pairing constant in the S region, and the summation over $\omega_{n}$ is cut off at the Debye frequency $\Omega$.

Eqs. (2)-(5) are supplemented by the following boundary conditions at $x=-d$ and $x=\infty$:

[TABLE]

Here $\theta_{s}$ defines the bulk Green’s function satisfying Eq. (2) with $\theta^{\prime\prime}\to 0$. We also use the standard boundary conditions at the N-S interface: 1988_Kuprianov_Lukichev

[TABLE]

Here $\theta_{-}\equiv\theta|_{x=-0}$, $\theta_{0}\equiv\theta|_{x=+0}$, $R_{B}$ is the N-S contact resistance, and $\sigma_{n}$ and $\sigma_{s}$ are the normal-state conductivities in the N and S regions, respectively. It is convenient to define the dimensionless parameters:

[TABLE]

where $N_{n}$ and $N_{s}$ are the normal densities of states in N and S regions, and $\xi_{n}=\sqrt{D_{n}/2\Delta}$ and $\xi_{s}=\sqrt{D_{s}/2\Delta}$ are the respective coherence lengths. More general boundary conditions for quasiclassical Green’s functions are given in Refs. 1999_Nazarov, ; 2015_Eschrig, . Yet Eqs. (7) and (8) can be used if the NS interface has a small transmission coefficient $t\sim d/\beta\xi_{0}\ll 1$ namely $d/\xi_{0}\ll\beta<\infty$, where $\xi_{0}$ is a coherence length in the clean limit 2017_Gurevich_Kubo . For a thin N layer with $d\ll\xi_{s}<\xi_{0}$, the condition $t\ll 1$ includes the essential cases with both $\beta>1$ and $\beta\ll 1$.

Calculation of the electromagnetic response requires retarded Green’s functions $G^{R}=\cosh\theta$ and $F^{R}=\sinh\theta$. In the case of dc currents and magnetic fields, $\theta$ satisfies the real-frequency Usadel equation:

[TABLE]

where $\tilde{\epsilon}=\epsilon+i\Gamma$. The quasiparticle DOS is given by

[TABLE]

For a uniform superconductor with $s=0$, $G^{R}=\tilde{\epsilon}/\sqrt{\tilde{\epsilon}^{2}-\Delta^{2}}$ and Eq. (12) reduces to Eq. (1).

The pair potential $\Delta_{s}$ for an ideal surface in the case of weak pairbreaking ($s\ll\Delta,\Gamma\ll\Delta$) was calculated in Appendix A:

[TABLE]

where $\Delta=2\Omega\exp(-1/g)$.

II.2 Surface resistance

We consider here type-II superconductors $(\lambda\gg\xi)$ for which the nonlocal BCS electromagnetic response 1958_Mattis_Bardeen simplifies to the local relation $\textbf{J}_{\omega}=\sigma(\omega)\textbf{E}_{\omega}$ between the Fourier components of the current density $\textbf{J}_{\omega}$ and the electric field $\textbf{E}_{\omega}$. Here $\sigma(\omega)=\sigma_{1}-i\sigma_{2}$ is a complex conductivity, where $\sigma_{2}(\omega)=1/\mu_{0}\omega\lambda^{2}$ accounts for the Meissner effect, and the quasiparticle conductivity $\sigma_{1}(\omega)$ determines the rf dissipation. To express $R_{s}$ in terms of $\sigma_{1}$, we use $E_{\omega}(x)=-i\omega A_{\omega}(x)$ and $J_{\omega}(x)=-i\omega\sigma(x,\omega)A_{\omega}(x)$ and calculate the power dissipation per unit surface $R_{s}H_{0}^{2}/2=(1/2)\int_{0}^{\infty}{\rm Re}[EJ^{*}]dx$ by integrating the local power density $\sigma_{1}E^{2}(x)$. For the case shown in Fig. 1(a,b), this yields

[TABLE]

where $d\ll\lambda$. We used Eq. (14) previosly to calculate a low-field $R_{s}$ in a superconductor with a proximity-coupled N layer at the surface  2017_Gurevich_Kubo .

Extension of Eq. (14) to high rf fields requires taking into account nonlinearities of the electromagnetic response in both $\sigma_{1}$ and $\sigma_{2}=(\mu_{0}\omega\lambda^{2})^{-1}$. The dependence of $\lambda$ on $H_{0}$ usually referred to as the nonlinear Meissner effect nme1 ; nme2 ; nme3 ; nme4 ; nme5 ; nme5 , is rather weak at $T\ll T_{c}$, $\omega\ll\Delta$, and $H_{0}\lesssim 0.5H_{c}$, so it will be neglected in this work. By contrast $\sigma_{1}$ controlled by the current-induced DOS broadening and nonequilibrium kinetics of quasiparticles is far more sensitive to a low-frequency rf field than $\sigma_{2}$ determined by the net superfluid density. Calculation of the field-dependent $\sigma_{1}$ in a dirty limit can be done using the time-dependent Usadel equation and kinetic equations for quasiparticles under strong low-frequency rf field kopnin ; kramer . In this case Eq. (14) can describe the field-dependent surface resistance, if the quasiparticle conductivity $\sigma_{1}(H_{0})$ is averaged over the rf period, taking into account temporal oscillations of the DOS and the distribution function of quasiparticles caused by the current pair-breaking parameter $s(x,t)$. The nonlinear conductivity $\sigma_{1}(H_{0})$ derived in Appendix B is given by 2014_Gurevich_PRL :

[TABLE]

where $f(\epsilon,t)$ is a distribution function of quasiparticles and $M[\epsilon,\omega,x,s(x,t)]$ is a spectral function:

[TABLE]

For weak or low-frequency rf fields, $f(\epsilon)$ tends to the equilibrium Fermi distribution $f_{0}(\epsilon)=(e^{\epsilon/T}+1)^{-1}$, and Eq. (15) takes the form:

[TABLE]

This formula determines the local nonlinear conductivity in a type-II superconductor where both the magnitude of the rf field and $\sigma_{1}(x,H_{0})$ vary slowly over $\xi$. At low fields $H_{0}\ll(\omega/\Delta)^{3/4}H_{c}$ and frequencies $\omega\ll T$, the surface resistance $R_{s}=\mu_{0}^{2}\omega^{2}\lambda^{3}\sigma_{1}/2$ can be calculated from Eqs. (16) and (17) using $G^{R}=\epsilon/\sqrt{\epsilon^{2}-\Delta^{2}}$, $F^{R}=\Delta/\sqrt{\epsilon^{2}-\Delta^{2}}$ at $\Gamma=0$ and $d=0$. This yields Zmuidzinas ; 2017_Gurevich_SUST ; 1958_Mattis_Bardeen :

[TABLE]

where $\rho_{n}=1/\sigma_{n}$, and $C\approx 9/2$. The logarithmic factor in Eq. (18) results from two close square root singularities at $\epsilon=\pm\Delta$ and $\epsilon=\pm(\Delta+\omega)$ in $M(\epsilon,\omega)$. These singularities characteristic of the idealized BCS model disappear as the DOS peaks are broaden by pair-breaking current 2014_Gurevich_PRL , bulk subgap states 2017_Gurevich_SUST or realistic surface features 2017_Gurevich_Kubo . In turn, strong rf fields can drive quasiparticles out of equilibrium, so $f(\epsilon)$ in Eq. (15) should in general be calculated by solving a kinetic equation kopnin ; kramer .

II.3 Nonequilibrium effects

Deviation of $f(\epsilon,x,t)$ from $f_{0}(\epsilon)$ is determined by an absorbed rf power and a rate of the power transfer from quasiparticles to phonons. At $T\ll\Delta$ the power transfer bottleneck is provided by scattering of quasiparticles on phonons which determine an inelastic scattering time $\tau_{s}$ and a recombination time of Cooper pairs $\tau_{r}$, where $\tau_{s}$ and $\tau_{r}$ depend on $T$ and $\epsilon$  kopnin ; kaplan . For instance, $\tau_{r}$ and $\tau_{s}$ at $\epsilon=\Delta$ in the absence of the rf field are given by kaplan

[TABLE]

where $\tau_{1}$ and $\tau_{2}$ are materials constants. For Nb with $T_{c}=9.2\,{\rm K}$, $\Delta=1.9T_{c}$, $\tau_{1}=3\cdot 10^{-12}\,{\rm s}$, and $\tau_{2}=8\cdot 10^{-11}\,{\rm s}$  kaplan , Eq. (19) gives $\tau_{r}\sim 4\cdot 10^{-8}\,{\rm s}$ and $\tau_{s}\simeq 2\cdot 10^{-8}\,{\rm s}$ at $T=2\,{\rm K}$. As $T$ increases to 4K, the time constants $\tau_{r}\simeq 4\cdot 10^{-10}\,{\rm s}$ and $\tau_{s}\simeq 2\cdot 10^{-9}\,{\rm s}$ diminish but remain much longer than the relaxation time constant of the superconducting condensate, $\hbar/\Delta\sim 10^{-12}$ s. At $T\ll T_{c}$ and $\omega\ll\Delta/\hbar$ an exponentially small density of quasiparticles has practically no effect on the dynamics of the condensate which reacts nearly instantaneously to the time-dependent currents. As a result, the spectral function $M$ in Eq. (16) is determined by the quasi-static Green’s functions. However, $\sigma_{1}$ can be very sensitive to a slow dynamics of nonequilibrium quasiparticles which control the distribution function $f$. Depending on the magnitude of the rf field, and the relation between $\omega$ and the relaxation time constants, quasiparticles can either follow the temporal variations of $M(t)$ if $(\tau_{r},\tau_{s})\omega\gg 1$ or relax quickly to the equilibrium state if $(\tau_{r},\tau_{s})\omega\ll 1$. Both $\tau_{r}$ and $\tau_{s}$ increase strongly as $T$ decreases, so nonequilibrium effects become more pronounced at $T\ll T_{c}$.

We evaluate non-equilibrium effects in $f$ using the kinetic equation for the matrix distribution function in a dirty superconductor in a low-frequency magnetic field kopnin ; kramer ; lo . In this case the quasiparticle electron-hole symmetry and the charge neutrality are preserved and the matrix kinetic equation reduces to only one equation for a correction to the odd in $\epsilon$ distribution function $\delta f_{\epsilon}(t)$. It is convenient to define

[TABLE]

where the amplitude $h_{\epsilon}(x,t)$ quantifies deviation of $f$ from equilibrium (nonequilibrium effects are negligible if $h_{\epsilon}\ll 1$). Using Eq. (13) and the parameterization $G^{R}=\cosh(u+iv)$, $F^{R}=\sinh(u+iv)$, the kinetic equation obtained in Ref. kramer, can be recast in the form:

[TABLE]

where $I_{ph}$ is the electron-phonon collision integral:

[TABLE]

Here phonons are assumed to be at equilibrium, and the low-energy phonon spectral function $\mu(\omega)=\alpha^{2}(\omega)F(\omega)=b\omega^{2}$ is used, where the constant $b$ was calculated in Ref. kaplan, for different materials.

We neglect the diffusion term in Eq. (21) assuming that the rf currents vary slowly over $\xi$ and the diffusion length of quasiparticles. We also neglect the contribution of $h_{\epsilon^{\prime}}$ to the collision integral (42) for qualitative evaluation of $h_{\epsilon}$ (see Ref. kramer, ), so Eq. (21) simplifies to:

[TABLE]

Here $\gamma_{\epsilon}$ is the rate of electron-phonon inelastic collisions, and $F_{\epsilon}(t)$ is a driving term. Because $u_{\epsilon}(t)$ and $v_{\epsilon}(t)$ depend on $s(t)$, both $\gamma_{\epsilon}(t)$ and $F_{\epsilon}(t)$ are functions of time.

At low frequencies $\omega\ll\gamma_{\epsilon}$, the term $\partial h/\partial t$ in Eq. (23) can be neglected, and

[TABLE]

At high frequencies $\omega\gg\gamma_{\epsilon}$, the electron-phonon collisions term in Eq. (23) can be neglected, and

[TABLE]

The functions $h_{\epsilon}(t)$ given by Eqs. (26) and (27) become of the same order of magnitude at $\omega\sim\gamma_{\epsilon}$.

Shown in Figs. 2 and 3 are the nonequilibrium amplitudes $h_{\epsilon}$ calculated from Eqs. (26) and (27). Here $\theta=u+iv$ was calculated by numerically solving the uniform quasistatic Usadel equation (11). Hence, it follows that: 1. The magnitude of $h_{\epsilon}$ increases as $T$ decreases, consistent with the above qualitative analysis. 2. The function $h_{\epsilon}$ has broadened peaks at $\epsilon\simeq\pm\Delta$ and decreases rapidly at $|\epsilon|>\Delta$ as $|\epsilon|$ increases. 3. Nonequilibrium subgap states which appear in $h_{\epsilon}$ become more pronounced as the rf field increases and $T$ decreases. 4. As the subgap parameter $\Gamma$ increases, the amplitude of the gap peaks in $h_{\epsilon}$ diminishes.

As follows from Figs. 2 and 3, the nonequilibrium effects become negligible $(h_{\epsilon}\ll 1)$ as $\omega$ and $H_{0}$ decrease or $\Gamma$ increases. For instance, at $H_{0}\lesssim 0.5H_{c}$ and $\omega\simeq(2-3)\cdot 10^{-3}\Delta$ (which corresponds to 1-2 GHz for Nb), the condition $h_{\Delta}\lesssim 1$ is satisfied at 4K at $\Gamma/\Delta=0.05$ and at 2K at $\Gamma/\Delta=0.1$. The reduction of $h_{\epsilon}$ as $\Gamma$ increases reflects acceleration of the energy relaxation of quasiparticles due to a finite DOS at $|\epsilon|<\Delta$ and the smearing out the DOS gap singularity. In turn, faster energy relaxation and diminishing nonequilibrium effects are facilitated by a proximity coupled N layer at the surface which provides a source of low-energy quasiparticles with $\epsilon<\Delta$. This layer plays a role of a ”quasiparticle trap” which can reduce the density of nonequilibrium quasiparticles in a superconductor qpt1 ; qpt2 . Given the wide range of the materials parameters, $H_{0}$, $\omega$ and $T$ where the nonequilibrium effects are not of prime importance, we assume here the Fermi-Dirac distribution in Eq. (15) to identify the material features which can be tuned to optimize $R_{s}(H_{0})$.

III An ideal surface with magnetic impurities in the bulk

Consider the effect of pairbreaking magnetic impurities on the nonlinear surface resistance at $\Gamma=0$. Here $G^{R}=\cosh\theta$ and $F^{R}=\sinh\theta$ in the spectral function $M$ are calculated from a quasi-static Usadel equation (11) for $\theta=u+iv$, where the low-frequency rf currents $(\omega\ll\Delta)$ vary slowly over $\xi$ so $\theta^{\prime\prime}$ can be neglected:

[TABLE]

Separating real and imaginary parts of Eq. (28) gives a quadratic equation for $\sin v$ and a cubic equation for $\sinh 2u$ with the solutions 2014_Gurevich_PRL :

[TABLE]

The density of states $n={\rm Re}G^{R}=\cosh u\cos v$ calculated from Eqs. (29)-(31) is shown in Fig. 4 for different values of the pairbreaking parameter $s$. Here the DOS vanishes at energies $|\epsilon|<\epsilon_{g}$, where  1969_Maki :

[TABLE]

Because of current and paramagnetic pairbreaking effects the quasiparticle gap $\epsilon_{g}$ is smaller than $\Delta$, vanishing at $s>0.43\Delta_{s}$. If $\Gamma_{p}<0.22$ the gap $\epsilon_{g}$ remains finite even at the superheating field $H_{\rm sh}\simeq 0.84H_{c}$, consistent with the previous result obtained for $\Gamma_{p}=0$  nme6 .

III.0.1 Surface resistance

We calculate the nonlinear surface resistance $R_{s}(H_{0})$ using Eqs. (14), (16), (17) for an ideal surface $(d=0)$. In this case the $x$-integration in Eq. (14) can be changed to integration over $s$, giving

[TABLE]

where $\tau=\omega t$. Shown in Fig. 5 is $R_{s}(H_{0})$ calculated from Eqs. (34) and (35) for different values of $\Gamma_{p}$, where $R_{s}(H_{0})$ is normalized to $R_{\rm MB}$ at $H_{0}=0$ and $\Gamma_{p}=0$. The dashed curve corresponds to the case of no magnetic impurities $(\Gamma_{p}=0)$ considered in Ref. 2014_Gurevich_PRL, .

One of the striking features of $R_{s}(H_{0})$ shown in Fig. 5 is a minimum in $R_{s}(H_{0})$ resulting from the interplay of the current-induced broadening of the DOS peaks which reduces $R_{s}$ and the reduction of $\epsilon_{g}$ which increases $R_{s}$  2014_Gurevich_PRL . This effect of the rf field suppression of $R_{s}(H_{0})$ can be qualitatively understood from Eq. (18) in which the logarithmic term results from two close square root singularities in $M(\epsilon)$ for the idealized DOS of the BCS model at $\omega\ll T$. This logarithmic singularity at $\omega=0$ is suppressed by either magnetic or current pairbreaking or the Dynes DOS broadening. As was pointed out above, any weak broadening of the DOS peaks reduces $R_{s}$. For instance, the pairbreaking effects take over if they broaden the DOS peaks to the width $\delta\epsilon\sim s^{2/3}\Delta^{1/3}$ exceeding $\omega$, so that the gap singularity in $M(\epsilon)$ is cut off at the field-dependent $\delta\epsilon$ rather than $\omega$. Thus, $R_{s}$ at $s>\omega^{3/2}\Delta^{-1/2}$ can be roughly evaluated by replacing $\omega\to\delta\epsilon$ in Eq. (18), giving a logarithmic decrease of $R_{s}$ with $H_{0}$:

[TABLE]

The decrease of $R_{s}(H_{0})$ with $H_{0}$ eventually stops at the field $H_{m}$ defined by $\delta\epsilon(H_{m})\sim T$. Hence,

[TABLE]

At $H_{0}=H_{m}$ and $\omega\lesssim\Gamma_{p}^{2/3}\Delta^{1/3}$ the surface resistance reaches a minimum which is by the factor $\sim\ln(T/\Delta^{1/3}\Gamma_{p}^{2/3})$ smaller than $R_{s}(0)$. As $\Gamma_{p}$ increases, the minimum in $R_{s}(H_{m})$ gets shallower and shifts to smaller fields. The minimum in $R_{s}(H_{0})$ disappears above a critical concentration of magnetic impurities defining the threshold pairbreaking parameter $\Gamma_{p}\simeq T^{3/2}\Delta^{-1/2}$ which decreases with $T$. This qualitative consideration is consistent with the numerical calculations of $R_{s}(H_{0})$ done at $T/\Delta=0.11$ and $\omega/\Delta=0.004$. For this case shown in Fig. 5, the minimum disappears at $\Gamma_{p}\gtrsim 0.01\Delta$.

As follows from Fig. 5, moderate magnetic pairbreaking decreases $R_{s}(H_{0})$ at low fields while increasing $R_{s}(H_{0})$ at higher fields. Here $\Gamma_{p}\simeq 0.01\Delta$ corresponds to an optimum DOS broadening due to magnetic pairbreaking. At higher fields the combined magnetic and current pairbreaking parameter $s$ exceeds the optimal value and $R_{s}(H_{0})$ increases sharply with $H_{0}$ so the position of the minimum in $R_{s}(\Gamma_{p})$ depends on $H_{0}$. Shown in Fig 6 are $R_{s}(\Gamma_{p})$ curves calculated for various $H_{0}$ and $T$. Here the optimal $\Gamma_{p}$ proportional to the density of magnetic impurities decreases as $H_{0}$ increases but increases as $T$ increases. For $\Gamma_{p}=\hbar v_{F}/2\ell_{s}\sim 0.01\Delta$, the optimum spin-flip scattering mean free path $\ell_{s}\sim 10^{2}\xi_{0}$ is much larger than $\xi_{0}$ in the absence of low-energy bound states on magnetic impurities 2017_Gurevich_Kubo .

The evolution of $R_{s}(H_{0})$ as temperature increases is shown in Fig. 7. Here the minimum in $R_{s}(H_{0})$ becomes shallower and shifts to lower fields as $T$ decreases. This behavior is consistent with the above qualitative discussion around Eqs. (36) and (37) that at low temperatures $T\sim\delta\epsilon\sim\Delta(s/\Delta)^{2/3}$, the effect of DOS broadening on $R_{s}(H_{0})$ weakens and $R_{s}(H_{0})$ is dominated by the reduction of the quasiparticle gap $\epsilon_{g}(H_{0})$ with $H_{0}$. The evolution of $R_{s}(H_{0})$ with the frequency $\omega$ is similar to the effect of $\Gamma_{p}$ and $T$ considered above. For instance, the results of our calculations presented in Fig. 8 show that the minimum in $R_{s}(H_{0})$ becomes shallower and shifts to lower fields as $\omega$ increases, and disappears at $\omega\gtrsim T$. It should be noted that this behavior of $R_{s}(H_{0},T,\omega,\Gamma_{p})$ obtained for the equilibrium distribution function of quasiparticle can change at low temperatures where the nonequilibrium effects become important. For instance, the condition of frozen density of quasiparticles caused by slow electron-phonon relaxation $\gamma_{\epsilon}\ll\omega$ at low $T$ can significantly deepen the minimum in $R_{s}(H_{0})$ as compared to the equilibrium distribution function 2014_Gurevich_PRL ; 2017_Gurevich_SUST

IV An Ideal surface with bulk subgap states

In this section we address the effect of bulk subgap states on the dependencies of the nonlinear surface resistance on the rf field, frequency and temperature.

IV.1 Green’s functions and density of states

Retarded Green’s functions are obtained by solving the real-frequency Usadel equation (28) with $\epsilon\to\epsilon+i\Gamma$ and $\Delta_{s}$ given by Eq. (13). This uniform Usadel equation for $\theta=u+iv$ yields two coupled equations for $v$ and $u$:

[TABLE]

The density of states $n=\cosh u\cos v$ calculated from Eqs. (38) and (39) is shown in Fig. 9. Here the gap peaks in the DOS are smeared even at $s=0$ and subgap states appear at $\epsilon=0$. The density of these subgap states increases with $s$. For instance, at $\epsilon=0$, the solution of the uniform Usadel equation in the first order in $(\Gamma,s)/\Delta\ll 1$ gives:

[TABLE]

where Eq. (13) was used.

IV.2 Surface resistance

The surface resistance at $\Gamma>0$ is calculated by generalization of Eq. (34):

[TABLE]

Here $\tau=\omega t$, and the spectral function $M$ is defined by Eq. (16), where $G^{R}$ and $F^{R}$ are obtained from the solutions of Eqs. (38) and (39).

Shown in in Fig. 10 is $R_{s}(H_{0})$ calculated for different values of $\Gamma$. Here the effect of subgap states on the field dependence of $R_{s}(H_{0})$ appears similar to the effect of magnetic impurities shown in Fig. 5. The dashed curve in Fig. 10 corresponds to $R_{s}(H_{0})$ for an ideal surface with no magnetic impurities and bulk subgap states 2014_Gurevich_PRL . The minimum in $R_{s}(H_{0})$ comes from the interplay of the broadening of the DOS peaks and a reduction of the quasiparticle gap by current. Here the DOS peaks are broaden by both rf currents and a finite lifetime of quasiparticle, so the minimum in $R_{s}(H_{0})$ shifts to lower fields as $\Gamma$ increases. If $\Gamma$ exceeds a critical value $\Gamma_{c}\sim T^{3/2}\Delta^{-1/2}$ (see Eq. (37)), the optimum DOS peak width occurs at $H_{0}=0$, and $R_{s}(H_{0})$ increases monotonically with $H_{0}$. Yet the effect of subgap states on the temperature dependences of $R_{s}(H_{0},T)$ turns out to be rather different from that of magnetic impurities.

Shown in Fig. 11 are the Arrhenius plots of $R_{s}$ versus $1/T$ calculated for different values of $s$ and $\Gamma$. At higher temperatures, $\ln R_{s}$ follows the linear dependence of $1/T$ expected from Eq. (18) of the BCS model. Here moderate broadening of the DOS peaks due to current pairbreaking and subgap states reduce $R_{s}$. At lower temperatures, $\ln R_{s}$ deviates from the linear dependence and levels off at a finite residual surface resistance, which comes from a finite density of subgap states at $\epsilon=0$ 2017_Gurevich_Kubo ; 2017_Gurevich_SUST . Because $n(\epsilon=0)$ depends on $s$ (see Eq. (40)), the residual surface resistance also depends on the rf field.

Figure 12 shows the effect of temperature on the field dependence of $R_{s}(H_{0})$. Here the minimum in $R_{s}(H_{0})$ becomes shallower and shifts to lower fields as $T$ decreases, similar to that is shown in Fig. 7. At the lowest temperature in Fig.  12, the monotonic increase of $R_{s}(H_{0})$ with $H_{0}$ mostly reflects the field dependence of the residual surface resistance. The effect of the rf frequency on the field dependence of $R_{s}$ shown in Fig. 13 is similar to that in Fig. 8 for magnetic impurities.

V S*′*-I-S surface nanostructuring

In this section we consider a way of tuning the nonlinear surface resistance by a higher - $T_{c}$ or $H_{c}$ superconducting (S′) layer separated by a thin dielectric (I) layer from the bulk superconductor (S) as shown in Fig. 1b. Here the I layer is assummed thick enough to suppress the Josephson coupling of the S and S′ layers which remain coupled through the electromagnetic interaction of screening currents. The use of S′-I-S multilayers was proposed to increase the superheating field and the onset of dissipative penetration of vortices in superconductors under strong rf fields mlag1 . It has been shown mlag2 ; mltk1 ; mltk2 ; mlc that, if the penetration depth $\lambda^{\prime}$ of the S′ layer exceeds $\lambda$ of the S-substrate, the S′-I-S structure can withstand the fields exceeding the superheating fields of both S and S′ regions at the optimal thickness of the S′ layer $d\sim\lambda^{\prime}$. This results from a current counterflow mlag2 ; mltk1 induced in the S′ layer by the S-substrate with $\lambda<\lambda^{\prime}$. Here we address manifestations of this counterflow effect in the field-dependence of $R_{s}(H_{0})$ for a S′-I-S structure.

V.1 Density of states

The retarded Green’s functions and $\Delta(x)$ are calculated from the Usadel equations in the $S$ and $S^{\prime}$ regions with the respective pair-breaking parameters $s(x)$ obtained by solving the London equations with the boundary conditions $\partial_{x}A(-d)=\mu_{0}H(t)$ and the continuity of $A(x)$ and $\partial_{x}A(s)$ at the $S-S^{\prime}$ interface at $x=0$. Using the solution of the London equation mlag2 ; mltk1 , we obtain:

[TABLE]

The amplitudes $s_{0}$ and $s_{0}^{\prime}$ are given by

[TABLE]

Here the prime corresponds to the respective parameters of the $S^{\prime}$ layer, and the factor $C$ describes the field attenuation by the $S^{\prime}$ layer at the $S$-$S^{\prime}$ interface.

Shown in Fig. 14 is the DOS calculated for a ${\rm Nb_{3}Sn}$ layer on the bulk ${\rm Nb}$. Here the DOS in the $S^{\prime}$ layer is less broadened as compared to the $S$ region, even though it is the $S^{\prime}$ layer which is directly exposed to the external field. Such reduction of the pairbreaking effect in the $S^{\prime}$ layer results from a counterflow induced by the $S$ substrate if $\lambda^{\prime}>\lambda$ mlag2 ; mltk1 . This condition is satisfied for Nb3Sn-Nb multilayers with typical values of $\lambda^{\prime}\simeq(2-3)\lambda$.

V.2 Surface resistance

The surface resistance of the $S^{\prime}$-$I$-$S$ structure is:

[TABLE]

Here the first term in the brackets describes the contribution of the $S$ region in the same way as in Eq. (41), where $I(s,\omega)$ is defined by Eq. (42), and $C^{2}$ accounts for the field attenuation by the $S^{\prime}$ layer. If $S$ and $S^{\prime}$ are made of the same material, Eq. (47) reduces to Eq. (41).

We first use Eq. (47) to calculate $R_{s}(H_{0})$ for a dirtier $S^{\prime}$ layer deposited onto the same but a cleaner material $S$ with $\Delta=\Delta^{\prime}$. For instance, Fig. 15 shows the results for a dirty Nb layer on a cleaner Nb substrate. Here the dirtier layer always increases $R_{s}(H_{0})$ at low fields consistent with the Mattis-Bardeen Eq. (18). However, as $H_{0}$ increases $R_{s}(H_{0})$ at $d\sim(0.25-0.5)\lambda^{\prime}$ becomes smaller than $R_{s}(H_{0})$ at $d=0$. Mitigation of the high-field increase of $R_{s}(H_{0})$ by a dirty surface layer reflects the counterflow effect which reduces the surface current in the $S^{\prime}$ layer, while the $S^{\prime}$ layer partially screens the bulk $S$ region. As a result, the minimum in $R_{s}(H_{0})$ shifts to a higher field at $d\simeq d_{m}$, where $d_{m}$ is an optimal thickness at which the $S$-$I$-$S^{\prime}$ structure can screen the magnetic field exceeding the bulk superheating fields of both $S$ and $S^{\prime}$   mlag2 ; mltk1 . The minimum in $R_{s}(H_{0})$ moves back to $H_{0}\sim 0.4H_{c}$ if $d\gtrsim 2\lambda^{\prime}$.

Next we calculate $R_{s}(H_{0})$ for a $S^{\prime}$-$I$-$S$ structure made of different materials with $\Delta^{\prime}>\Delta$ in which case the behavior of $R_{s}(H_{0})$ changes significantly as compared to the above case of $\Delta=\Delta^{\prime}$. For instance, Fig. 16 shows $R_{s}(H_{0})$ calculated for a ${\rm Nb_{3}Sn}$-${\rm Nb}$ multilayer. Because $\Delta^{\prime}_{\rm{Nb_{3}Sn}}\simeq 2\Delta_{\rm{Nb}}$, the surface resistance of Nb3Sn is smaller than $R_{s}$ for Nb at $T\ll T_{c}$. As a result, $R_{s}$ reduces as $d$ increases, in contrast to the case for a Nb-Nb multilayer. However, thick ${\rm Nb_{3}Sn}$ films with $d\gg\lambda^{\prime}$ are prone to the dissipative vortex penetration at lower fields $H_{0}\sim H_{c1}^{\rm{Nb_{3}Sn}}\ll H_{c1}^{\rm{Nb}}$, resulting in a significant increases of $R_{s}$. Thus, we are interested in optimum thicknesses $d\lesssim\lambda^{\prime}$ for which the field onset of vortex penetration is greatly increased mlag1 while the screening field limit is pushed up beyond the superheating fields of both $S$ and $S^{\prime}$ superconductors mlag2 ; mltk1 . As an illustration, Fig. 16 shows $R_{s}(H_{0})$ calculated for $\Gamma^{\prime}=0.05\Delta$ and $\Gamma=0.01\Delta$. In this case the minimum in $R_{s}(H_{0})$ in a thick $S^{\prime}$ layer with $d\gg\lambda^{\prime}$ disappears due to the large $\Gamma^{\prime}$. However, a thin $S^{\prime}$ coating with $d\sim\lambda^{\prime}$ not only preserves the minimum in $R_{s}(H_{0})$ but also pushes it to higher fields exceeding $H_{c}$ of the $S$ substate. These results suggest that a thin $\rm{Nb_{3}Sn}$ coating can be used to produce high gradient resonant cavities with no field degradation of the quality factors $Q(H_{0})$ which can exceed $Q^{\rm{Nb}}(H_{0})$ up to fields close to the superheating field of ${\rm Nb_{3}Sn}$.

VI Proximity coupled normal layer at the surface

In this section we consider the effect of a thin N layer on $R_{s}(H_{0})$, generalizing our approach developed for the calculations of $R_{s}$ at low fields 2017_Gurevich_Kubo to the case of strong current pairbreaking in type-II superconductors with $\lambda\gg\xi$ and proximity-coupled systems 1970_Usadel ; 1999_Belzig_review ; 2004_Golubov_review . For the geometry shown in Fig. 1(b), $\theta(x)$ is nearly constant across a thin N layer with $d_{N}\ll\xi_{N}$, so the problem can be reduced to solving the Usadel equation only in the S region and taking into account the effect of the N layer by a self-consistent boundary condition at $x=0$. In what follows we outline the essential steps of this approach and give the final formulas necessary for the calculations of $R_{s}(H_{0})$, relegating all technical details to Appendix C.

VI.0.1 Retarded Green’s functions and density of states

The real-frequency Usadel equation for the functions $u$ and $v$ in the S-region can be written in the form

[TABLE]

where

[TABLE]

Here $u_{s}$ and $v_{s}$ are the solutions of the uniform Usadel Eq. (38) and (39). As shown in Appendix C, Eqs. (48) and (49) are supplemented by the following boundary condition on the S-side at $x=+0$:

[TABLE]

According to the boundary condition (7), the complex function $\theta(x)=u(x)+iv(x)$ has a jump at $x=0$ for any nonzero contact resistance $R_{B}$. Here $\theta(-0)=u_{-}+iv_{-}$ on the N side of the N-S interface is related to $\theta(x+0)=u_{0}+iv_{0}$ by two coupled equations:

[TABLE]

where $r=D_{S}/D_{N}$, and the parameter $\Psi$ in Eq. (52) and (53) which accounts for a short-range disturbance of the pair potential caused by the N layer in the S region was calculated in Appendix C:

[TABLE]

Here $\Delta(x)$ is given by Eq. (13), and $\theta_{s}$ is determined by the uniform thermodynamic Usadel equation:

[TABLE]

The closed set of Eqs. (48)-(58) determine self-consistently the Green functions and the pair potential. We solved these equations numerically to calculate the spectral function $M(\epsilon,s)$ and the DOS which are the key input characteristics in the calculations of $R_{s}(H_{0})$. For instance, Fig. 17 shows the DOS at (a) the N and (b) S side of the interface at $x=0$. For a small contact resistance $R_{B}$, that is, a nearly transparent interface with $\beta\ll 1$ (solid curves), the N and S regions are strongly coupled and the behaviors of both $N_{N}(\epsilon)$ and $N_{S}(\epsilon)$ are similar to that of $N(\epsilon)$ shown in Fig. 9 for an ideal surface. As $R_{B}$ and $\beta$ increases (dashed curves), the proximity coupling of the N and S regions weakens, and minigap states in the N layer appear. These states are further broaden and shifted to lower energies by current.

Figure 18 shows the evolution of the DOS in the S region with the distance from the NS interface calculated for: (a) $s/\Delta=0.005$ and (b) $s/\Delta=0.1$. As $x$ increases, the effect of the N layer on the DOS weakens and $N(\epsilon)$ approaches the bulk DOS in the S region. A spatial extent of the DOS disturbance by the proximity effect is given by $L\sim\xi_{s}(\Delta/\delta\epsilon)^{1/4}$, where $\delta\epsilon$ is the DOS peak width  2017_Gurevich_Kubo . As $s$ increases $\delta\epsilon$ increases and the proximity induced DOS disturbance becomes more short-ranged.

VI.0.2 Surface resistance

The surface resistance is given by Eqs. (47):

[TABLE]

Here the first term in the brackets comes from the contribution of the N layer. The integral term describes the contribution of the S region, where $I(x)$ defined by Eq. (42) is calculated by solving Eqs. (48)-(58) taking into account both the proximity effect and the spatial variation of the current pairbreaking parameter $s(x)$.

Shown in Fig. 19 is $R_{s}(H_{0})$ calculated from Eq. (59) for different thicknesses of the N layer, where the dashed curves correspond to $R_{s}(H_{0})$ for an ideal surface with $\Gamma=0.005\Delta$. In all cases the minimum in $R_{s}(H_{0})$ is due to the interplay of the current-induced DOS broadening and the reduction of the quasiparticle gap. As a result, the N layer shifts the minimum in $R_{s}(H_{0})$ to lower fields as $d$ increases because the proximity-induced DOS broadening becomes more pronounced as the N layer gets thicker, so that the optimum width of DOS peaks is achieved at smaller $H_{0}$. For strongly-coupled N layers represented in Fig. 19(a), the minimum in $R_{s}(H_{0})$ disappears as the N layer thickness exceeds a critical value $d>d_{c}=\alpha\xi_{s}N_{s}/N_{n}$ which is still much smaller than $\xi_{s}$ at $\alpha\simeq 0.1$.

As follows from Fig. 19(a), for thin strongly-coupled N layers with $\alpha\lesssim 0.05$, there is a crossover in the field dependence of $R_{s}(H_{0})$: the low-field $R_{s}(H_{0})$ is smaller than the corresponding $R_{s}$ for an ideal surface but this relation reverses as $H_{0}$ increases. This reduction of $R_{s}$ at small $d$ comes from the long range proximity effect in which the length $L\sim\xi_{s}(\Delta/\delta\epsilon)^{1/4}$ of the DOS disturbance produced by the N layer extends deep into the S region providing an optimum DOS broadening  2017_Gurevich_Kubo . The resulting reduction of $R_{s}(H_{0}=0)$ in the S region overwhelms the increase of $R_{s}$ due to dissipation in the N layer if $\delta\epsilon\sim\Gamma\ll\Delta$. However, as $H_{0}$ increases, the current-induced DOS broadening takes over, shortening the range of DOS disturbance $L$ as shown in Fig. 18 and increasing $R_{s}$ in the S region. It should be noted that this behavior of $R_{s}(H_{0})$ is only characteristic of strongly-coupled N layers with $\beta\ll 1$. For weakly coupled N layers with large contact resistances and $\beta>1$, the field crossover in $R_{s}(H_{0})$ disappears, as shown in Fig. 19(b).

The effects of contact resistance on $R_{s}(H_{0})$ is shown in Fig. 20. In the case of $D_{n}=0.5D_{s}$ represented by Fig. 20(a) the proximity effect reduces $R_{s}$ at low fields at $\beta\lesssim 0.1$ but the optimal $\beta$ at high RF fields becomes much larger. For instance, the deepest minimum in $R_{s}(H_{0})$ and the lowest high field $R_{s}(H_{0})$ occur at intermediate coupling $\beta\simeq 0.5$. For larger $\beta\gtrsim 1$, the minigap reduction in the N layer becomes dominant mechanism which increases $R_{s}(H_{0})$. A rather different evolution $R_{s}(H_{0})$ with $\beta$ occurs for a highly conductive N layer with $D_{n}=10D_{s}$ represented in Fig.  20(b). In this case the N layer gives the main contribution to $R_{s}$ but the effect of the contact resistance on $R_{s}$ is rather intricate: $R_{s}$ first increases with $\beta$ reaching the maximum low-field values at $\beta\simeq 7$ and then decreases as $\beta$ further increases. Here $R_{s}(H_{0}))$ at $\beta\gtrsim 3$ exceeds $R_{s}$ of a completely decoupled N layer at $\beta\to\infty$. Such unexpected behavior can be understood by analysing the Green functions in the N and S regions shown in Fig. 21. At $\beta=7$ the DOS in the N layer is nearly equal to that of of normal state but ${\rm Re}F$ does not vanish. Due to the contribution of ${\rm Re}F$ in the spectral function $M$ in Eq. (16), the N layer still has proximity induced pair correlations but $R_{s}$ becomes larger than $R_{s}$ of a fully decoupled N layer with $\beta\to\infty$. As the field increases, ${\rm Re}F$ decreases, so $R_{s}$ approaches $R_{s}$ of the normal state.

Shown in Fig. 22 is the effect of diffusivity ratio $D_{n}/D_{s}$ on $R_{s}(H_{0})$ calculated for different values of $\beta$. As the ratio $D_{n}/D_{s}$ increases, the minimum first shifts to lower fields because the effect of current pair-breaking in the N layer becomes more pronounced so that $R_{s}(H_{0})$ turns to increase at lower fields. However, this trend reverses at larger $\beta$ as it is evident from Fig. 22(c), where the minimum disappears and the field dependence of $R_{s}(H_{0})$ is determined by the behavior of the Green functions discussed above in relation to $R_{s}(H_{0})$ for a nearly decoupled N layer shown in Figs. 20(b) and Fig. 21.

Now we turn to the effects of temperature and the rf frequency on the field dependence of $R_{s}(H_{0})$. Fig. 23 shows that the evolution of $R_{s}(H_{0})$ with $T$ can change as the interface resistance increases and the N layer becomes more decoupled. For instance, at $\beta=1$, the minimum in $R_{s}(H_{0})$ becomes shallower and shifts to lower fields as $T$ decreases. This behavior of $R_{s}(H_{0},T)$ is controlled by the same mechanism discussed above in relation to Figs. 7 and 12. However, this evolution of $R_{s}(H_{0})$ with $T$ reverses for $\beta=7$ corresponding to the solid curves in Fig. 23. In this case the minimum becomes more pronounced and shifts to higher fields as $T$ decreases. Here the field dependence of $R_{s}(H_{0})$ is mostly determined by the effect of the N layer similar to that is shown in Fig. 22(c), whereas the contribution of the S region diminishes at $T$ decreases. The effect of frequency on $R_{s}(H_{0})$ shown in Fig. 24 is qualitatively similar to that was discussed above in relation to Figs. 8 and 13. In both cases our calculations do not take into account nonequilibrium effects which become more pronounced at lower temperatures and higher frequencies.

VII Discussion

In this work we address mechanisms of nonlinear RF losses and show that $R_{s}(H_{0})$ can be reduced by optimizing the DOS using a multitude of pairbreaking mechanisms, including the N-S or S-I-S′ surface nanostructuring, impurity management or the rf Meissner currents. Here the nonmonotonic field dependence of $R_{s}(H_{0})$ results from the broadening of the DOS gap peaks by the pair-breaking rf currents  2014_Gurevich_PRL . Combining the rf current pairbreaking with pairbreaking caused by the material features 2017_Gurevich_SUST ; 2017_Gurevich_Kubo opens up opportunities to tune the field dependence of $R_{s}(H_{0})$ by materials processing and surface treatment. For instance, magnetic impurities shift the minimum in $R_{s}(H_{0})$ to lower fields and can reduce $R_{s}(H_{0})$ in a wide range of $H_{0}$ (Figs. 5). Bulk subgap states can give rise to a residual surface resistance while reducing $R_{s}$ at higher temperatures (see Figs. 10 and 11). A S-I-S′ multilayer with the optimal thickness $d\sim\lambda^{\prime}$ can shift the minimum in $R_{s}(H_{0})$ to higher fields exceeding the superheating field of the S-substrate. A proximity-coupled N layer not only affects the field position and the depth of the minimum in $R_{s}(H_{0})$ but can also reduce the surface resistance below $R_{s}$ for an ideal surface without the N layer (see Figs. 19, 20, and 22).

These results can explain why the field dependence of the quality factor $Q(H_{0})\propto 1/R_{s}(H_{0})$ can be very sensitive to the materials processing. For instance, magnetic impurities in Nb have been attributed to oxygen vacancies in ${\rm Nb_{2}O_{5}}$, while other impurities can segregate at the surface during chemical polishing and electro polishing  2005_Casalbuoni ; 2011_Proslier . It has been found that low-temperature baking increases the Curie constant, which in turn increases the spin flip parameter $\Gamma_{p}$ and affects $R_{s}(H_{0})$ curves as shown in Fig. 5. Tunneling experiments have shown that the DOS is strongly affected by infusion of Ti at high-temperatures, which sharpens the DOS peaks and results in a pronounced minimum in $R_{s}(H_{0})$ 2013_Dhakal . These observations are consistent with the results of this work in which the behavior of $R_{s}(H_{0})$ is controlled by the DOS at the surface and the values of $\Gamma$ and $\Gamma_{p}$ which manifest themselves in different $R_{s}(H_{0})$ curves shown in Fig. 10.

A common feature of many superconducting materials is a thin proximity-coupled N layer on the surface, for example, metallic suboxides underneath the dielectric ${\rm Nb_{2}O_{5}}$ pentooxide layer on the Nb surface. Chemical or heat treatments change the properties of the N layer, particularly the thickness, impurity concentration, and the interface contact resistance 1992_Robertazzi ; 1993_Ekin . As a result, the materials processing used for superconducting resonator cavities  2014_Padamsee_review ; 2012_Aontoine_review , such as electropolishing 1989_Saito and heat treatments at high temperatures 2013_Dhakal ; 2013_Grassellino and low temperatures 2017_Grassellino ; 2018_Dhakal ; 2018_Grassellino , can significantly affect $R_{s}(H_{0})$. These observations are qualitatively consistent with our calculations which show that various $R_{s}(H_{0})$ curves similar to those observed on the Nb cavities can occur, depending on the N layer properties and the value of $R_{B}$, as illustrated by Figs. 19, 20, and 22. Yet a quantitative comparison of our theory with experiment would require independent measurements of multiple parameters characterizing a particular material, for example, the thickness and conductivity of the N layer, $R_{B}$, $\Gamma_{N}$ and $\Gamma_{S}$, as well as the way these parameters change after different materials treatments. This information can be extracted from different surface characterization techniques, for example, scanning tunneling microscopy or X-ray photoelectron spectroscopy. Such surface measurements combined with the rf measurements of $Q(H_{0})$ and the theory analysis could clarify whether it is a reduction of the thickness of the N layer or contact resistance or the Dynes parameters in the bulk which plays a major role in the development of the field-induced reduction of $R_{s}(H_{0})$ after the Ti or N infusion or other heat treatments in Nb.

Tuning of $R_{s}(H_{0})$ can also be done by depositing a S′ layer with higher-$H_{c}$ and $\lambda^{\prime}>\lambda$ on the S surface, which can withstand magnetic fields of the order of the superheating field $H_{s}^{\prime}$ while significantly increasing the field onset of vortex penetration mlag1 . Moreover, the current counterflow caused by the S substrate in the S′ layer with $\lambda^{\prime}>\lambda$ allows the S-I-S′ structure to remain in the Meissner state at $H_{0}$ exceeding both superheating fields $H_{s}$ and $H_{s}^{\prime}$ at an optimum thickness $d\sim\lambda^{\prime}$  mlag2 ; mltk1 ; mltk2 . This counterflow effect can shift the minimum in $R_{s}(H_{0})$ to higher fields exceeding $H_{s}$ of the S-region at $d\simeq\lambda^{\prime}$, as shown in Fig. 16 for a Nb3Sn-I-Nb structure. Also, the Nb covered with a dirtier Nb layer, which models the Nb surface after low-temperature baking  muonSR , mitigates the high-field increase of $R_{s}(H_{0})$  2017_Gurevich_SUST ; hs3 as shown in Fig. 15. This conclusion is also qualitatively consistent with experiments on the Nb cavities (see, e.g., reviews 2012_Aontoine_review, ; 2014_Padamsee_review, or Refs. 2018_CEH, ; 2018_Grassellino, for recent developments).

For a superconductor with an ideal surface, the minimum in $R_{s}(H_{0})$ shown in Figs. 7 and 12 becomes shallower as $T$ decreases. On the other hand, for a weakly-coupled N layer with a large contact resistance $(\beta\gtrsim 1)$, the minimum in $R_{s}(H_{0})$ becomes more pronounced as the temperature decreases (see Fig. 23) because the N layer contribution to $R_{s}(H_{0})$ becomes dominant at low $T$. These conclusions follow from our calculations in which the equilibrium distribution function of quasiparticles $f_{0}(\epsilon)$ was assumed. As shown in Sec. II C, this condition is satisfied at low frequencies $\hbar\omega\ll k_{B}T$, the subgap states and minigaps in the N layer accelerate the energy relaxation of quasiparticles reducing the deviation of $f(\epsilon)$ from the Fermi distribution. Yet the calculations of $R_{s}(H_{0})$ at low temperatures and higher frequencies would require the non-equilibrium quasiparticle distribution function obtained by solving the kinetic equations for strong RF fields kopnin ; kramer ; lo .

We considered different ways of minimization of $R_{s}(H_{0})$ by optimizing the DOS at the surface. Yet it should be noted that the optimum parameters actually depend on the rf amplitude. For instance, Fig. 6 shows that the optimum spin flip parameter $\Gamma_{p}$ is different for different $H_{0}$, although the optimum concentration of paramagnetic impurities is roughly determined by the spin-flip scattering mean free path $\ell_{s}\sim 10^{2}\xi_{0}$, at which $R_{s}$ can be significantly reduced in a wide range of $H_{0}$. The subgap states affect $R_{s}$ in a way similar to that of paramagnetic impurities at higher temperatures. While the physics and materials mechanisms behind the Dynes parameter $\Gamma$ are not well understood, tunneling spectroscopy can give valuable insights into how $\Gamma$ is affected by various materials processing. The observed values of $\Gamma$ can then be used as input parameters in the theory of this work to identify the optimum materials treatments which minimize the nonlinear surface resistance.

Appendix A Pair potential

For a superconductor with an ideal surface, the solution of the thermodynamic Usadel equation (2) in the first order in $s/\Delta\ll 1$ is given by:

[TABLE]

Substituting Eq. (60) into Eq. (5) yields the equation for $\Delta$ in the case of weak pairbreaking:

[TABLE]

At $T\ll T_{c}$ we replace the summation in Eq. (61) with integration over $\omega_{n}$, express $1/g=\ln(2\Omega/\Delta_{0})$ in terms of the gap $\Delta_{0}$ at $s=T=0$ and obtain:

[TABLE]

At $s=0$ Eq. (62) yields:

[TABLE]

Hence, the subgap states reduce the superfluid density and fully suppress superconductivity at $\Gamma>\Delta_{0}/2$. For weak DOS broadening and current pairbreaking considered in this work, the solution of Eq. (62) in the first order in $s$ and $\Gamma$ is given by:

[TABLE]

Appendix B Nonlinear conductivity

Here we present the derivation of the time-averaged $\sigma_{1}(H)$ following the approach of Ref. (2014_Gurevich_PRL, ) based on the time-dependent Usadel equations for the $4\times 4$ quasiclassical Greens function $\check{G}(\textbf{r},t,t^{\prime})$ lo ; kopnin ; 1999_Belzig_review :

[TABLE]

where $\hat{G}^{R}$ and $\hat{G}^{A}$ are the retarded and advanced Green functions, $\check{G}\circ\check{G}=\check{1}$, $\hat{G}^{K}=\hat{G}^{R}\circ\hat{f}-\hat{f}\circ\hat{G}^{R}$ is the Keldysh function expressed in terms of a distribution function of quasiparticles, $\hat{f}(\epsilon,t)$, the hat denotes matrices in the Nambu space, $\hat{\Pi}=\nabla+i\pi\mathbf{A}\hat{\sigma}_{z}/\phi_{0}$, $\check{\sigma}_{z}$ is a diagonal matrix composed of the Pauli matrices $\hat{\sigma}_{z}$, and $\circ$ means time convolution. For a dirty type-II superconductor with $\lambda\gg\xi$, the relation between the current density and the vector potential is local but time-dispersive:

[TABLE]

where $\hat{G}_{z}=\hat{\sigma}_{z}\cdot\hat{G}$. In Eq. (77) the gradient terms $\nabla G$ are neglected, and the matrix $\hat{f}$ reduces to a single distribution function since the transverse electromagnetic field does not cause the electron-hole imbalance.

We use the standard mixed Fourier representation

[TABLE]

where $t_{0}=(t+t^{\prime})/2$. Expanding the time convolution in small derivatives with respect to $t^{\prime}$ in Eq. (77) for slow variations of $\textbf{A}(\textbf{r},t^{\prime})$ and neglecting linear in $\omega$ terms yields:

[TABLE]

The time-averaged nonlinear conductivity $\sigma_{1}=\bar{\sigma}$ is defined in terms of the mean dissipated power $q=\bar{\sigma}E_{0}^{2}/2$ induced by the ac electric field $\mathbf{E}=-c^{-1}\partial_{t}\mathbf{A=E}_{0}\sin\omega t$ for which $\mathbf{A}(t)=(c\mathbf{E}_{0}/\omega)\cos\omega t$. Here

[TABLE]

where

[TABLE]

Changing variables $t=t_{0}+t_{1}/2$ and $t^{\prime}=t_{0}-t_{1}/2,$ yields

[TABLE]

In the dirty limit $P_{1}(\epsilon,\epsilon^{\prime},t_{0})$ and $P_{2}(\epsilon,\epsilon^{\prime},t_{0})$ are even functions of $t_{0}$, so only $\sin\omega t_{1}$ contributes:

[TABLE]

Using $\mbox{Tr}\hat{G}_{z}\cdot\hat{G}_{z}=2G\cdot G+F\cdot F^{{\dagger}}+F^{{\dagger}}\cdot F$ gives

[TABLE]

Because $\mbox{Re}(G,F)^{R}=-\mbox{Re}(G,F)^{A}$, and $\mbox{Im}(G,F)^{R}=\mbox{Im}(G,F)^{A}$, we can present Eqs. (84) and (85) in the form

[TABLE]

In the parameterization $G^{R,A}_{\epsilon}=\pm\cosh u\cos v+i\sinh u\sin v$ and $F^{R,A}_{\epsilon}=\pm\sinh u\cos v+i\cosh u\sin v$, Eq. (86) becomes

[TABLE]

If $s=0$, we have $\cos v_{\epsilon}=1$, $\cosh u_{\epsilon}=$ $\epsilon/\sqrt{\epsilon^{2}-\Delta^{2}}$, $\sinh u_{\epsilon}=\Delta/\sqrt{\epsilon^{2}-\Delta^{2}}$ at $|\epsilon|>\Delta$. If $|\epsilon|<\Delta$ the integrand in Eqs. (86) and (87) vanishes and $\bar{\sigma}$ reproduces the Mattis-Bardeen result for $\omega<\Delta$:

[TABLE]

At $\exp(-\Delta/T)\ll 1$ the main contribution to this integral comes from a narrow range $\epsilon-\Delta\sim T\ll\Delta$, where one can take $z=\epsilon-\Delta$ and $\epsilon^{2}-\Delta^{2}\approx 2\Delta z$. Then

[TABLE]

Hence

[TABLE]

Here a modified Bessel function reduces to $K_{0}(\omega/2T)\simeq\ln(4T/\omega)-\gamma_{E}$ at $\omega\ll 2T$, where $\gamma_{E}=0.577$ is the Euler constant. Then Eq. (90) gives $R_{s}=\mu_{0}\omega^{2}\lambda^{3}\sigma_{1}/2$ which reproduces Eq. (18) with $C=4e^{-\gamma_{E}}\approx 9/2$.

Appendix C Proximity coupled normal layer

To calculate $\theta$ in the N region, we solve the Usadel equation with $\Delta=0$:

[TABLE]

Here $r=D_{\rm S}/D_{\rm N}$ and $\tilde{\omega}_{{\rm N}n}=\omega_{n}+\Gamma_{\rm N}$, and we use the dimensionless form with dimensionless lengths $x\to x/\xi_{s}$ and energies $\omega_{n}\to\omega_{n}/\Delta$, $s\to s/\Delta$, $\Gamma\to\Gamma/\Delta$.

If $d\ll\xi_{n}$ the function $\theta(x)$ varies weakly across the N layer, so the solution of Eq. (91) satisfying the boundary condition $\theta^{\prime}(-d)=0$ can be approximated by

[TABLE]

The boundary condition (7) at $x=0$ takes the form

[TABLE]

where $\alpha$ and $\beta$ are defined by Eqs. (9) and (10).

To obtain the relation between $\theta_{-}$ at the N side of the boundary and $\theta_{0}^{\prime}$ at the S side of the boundary, we solve the Usadel equation in the S region ($x\geq 0$):

[TABLE]

Here $\tilde{\omega}_{{\rm S}n}=\omega_{n}+\Gamma_{\rm S}$. In the vicinity of the NS interface, the proximity effect causes a weak short-range disturbance $\Delta(x)=\Delta_{s}+\delta\Delta(x)$. we can approximate by 2017_Gurevich_Kubo :

[TABLE]

with $\Psi$ to be determined self-consistently. Integrating Eq. (95) from [math] to $+0$ and using $\theta^{\prime}(-0)=\alpha\eta$, we obtain the boundary condition for Eq. (95) at $x=+0$:

[TABLE]

A thin N layer produces a weak disturbance of the pair potential, $\Psi\ll 1$  2017_Gurevich_Kubo , so Eq. (95) can be linearized in $\delta\theta(x)=\theta(x)-\theta_{s}$, giving $\delta\theta^{\prime\prime}=-k^{2}\delta\theta$ where $k$ is defined by Eq. (57). Hence,

[TABLE]

Substituting Eqs. (96) and (95) into Eq. (5) linearized with respect to $\delta\theta_{0}$ and integrating it over $x$ yields the self-consistency Eq. (56) for $\Psi$.

Consider now the retarded Green’s functions. The solution of the real frequency Usadel equation in the N region is similar to that of Eq. (92):

[TABLE]

where $\tilde{\epsilon}_{\rm N}=\epsilon+i\Gamma_{\rm N}$. In the S region we have:

[TABLE]

where $\tilde{\epsilon}_{\rm S}=\epsilon+i\Gamma_{\rm S}$, and $\Psi$ is given by Eq. (56). Substituting $\theta=u+iv$ into Eq. (100) and separating real and imaginary parts at $x=+0$ yields Eqs. (48) and (49).

The equation connecting the boundary values $\theta_{-}$ and $\theta_{0}$ readily follows from Eqs. (93) and (94):

[TABLE]

Taking here $\theta_{0}=u_{0}+iv_{0}$ and $\theta_{-}=u_{-}+iv_{-}$ and separating real and imaginary parts yields Eqs. (54) and (55). Another boundary condition is obtained by integrating Eq. (100) from $x=-0$ to $x=+0$, which yields $r\theta^{\prime}_{-}-\theta_{0}^{\prime}=i\Psi\cosh\theta_{0}$. Substituting here $\theta^{\prime}_{-}=\alpha(s\sinh\theta_{-}\cosh\theta_{-}-ir\tilde{\epsilon}_{n}\sinh\theta_{-})$ from Eq. (C) gives

[TABLE]

Separating here real and imaginary parts results in Eqs. (52) and (53).

Acknowledgements.

The work of T. K. was supported by Japan Society for the Promotion of Science (JSPS) KAKENHI Grant Number JP17H04839 and JP17KK0100. The work of A. G. was supported by NSF under Grant No. PHY-1416051 and DOE under Grant No DE-SC100387-020.