Self-consistent calculation of the flux-flow conductivity in diffusive superconductors

In the framework of Keldysh-Usadel kinetic theory, we study the temperature dependence of flux-flow conductivity (FFC) in diffusive superconductors. By using self-consistent vortex solutions we find the exact values of dimensionless parameters that determine the diffusion-controlled FFC both in the limit of the low temperatures and close to the critical one. Taking into account the electron-phonon scattering we study the transition between flux-flow regimes controlled either by the diffusion or the inelastic relaxation of non-equilibrium quasiparticles. We demonstrate that the inelastic electron-phonon relaxation leads to the strong suppression of FFC as compared to the previous estimates making it possible to obtain the numerical agreement with experimental results.


I. INTRODUCTION
Vortex motion is an important process that determines resistive properties of type-II superconductors in the fluxflow regime. At magnetic fields B much weaker then the upper critical one, H c2 , the density of vortex lines is small and the total electric losses are given by the superposition of the individual vortex contributions. In this regime, the flux-flow resistivity ρ f is proportional to the density of vortex lines, ρ f ∼ B/H c2 , as described by the general expression suggested by Bardeen and Stephen 1 .
The inverse quantity σ f = 1/ρ f , the flux-flow conductivity (FFC) is therefore given by where σ n is the normal state conductivity and β is the numerical coefficient which is determined by the particular microscopic model. For superconducting materials with high rate of impurity scattering, the numerical value of β ≈ 0.9 at low temperatures has been reported by Gor'kov and Kopnin 2 (GK). This value of the dimensionless parameter has been obtained using the approximate vortex solution. Up to date the exact value of β has been unknown and it is reported in the present paper based on the fully-self consistent vortex structure calculations.
At elevated temperatures, two different regimes of the vortex motion has been considered, depending on the dominant mechanism of the relaxation 3 . One of them is the diffusion-controlled flux-flow, when the generation of non-equilibrium quasiparticles near the vortex line is balanced by their diffusion to the infinity. As temperature approaches T c this mechanism results in the divergent behaviour of FFC given by [3][4][5] : with the temperature-independent β 0 . Qualitatively this behaviour is explained by the vortex core size increase proportional to the Ginzburg-Landau coherence length where D is the diffusion coefficient. This dependence is in the qualitative agreement with experimental results 6 pointing to the significant increase of β as temperature approaches T c . However, quantitative agreement is lacking. Initially, the value of β 0 ≈ 1.1 has been reported 4 which by coincidence was in the good agreement with experiments 6 . However subsequently this result has been revised to β 0 ≈ 4.04 by Larkin and Ovchinnikov 3,5 (LO) which is several times larger than the measured values in various superconductors [7][8][9][10][11][12][13] .
When the temperature becomes sufficiently close to T c the relaxation is dominated by the inelastic electronphonon collisions. This regime is described by generalized time-dependent Ginzburg-Landau theory (GTDGL) yielding the FFC decreasing with temperature 3 where τ ph is the electron-phonon relaxation time. In the limit T → T c the gappless superconducting state is realized. In this case the decrease of β(T ) saturates at β = 1.45 14 .
The crossover between two regimes described by the Eqs. (2,3) occurs at the temperatures τ ph (T c − T ) ∼ 1 when the diffusion rate becomes of the order of electronphonon relaxation rate, Dξ −2 GL ∼ τ −1 ph . That yields an estimation of the maximal value max(β) ∼ T c τ ph obtained from Eqs. (2) and (3) at the upper and lower borders of their applicability respectively.
Although the estimations of max(β) obtained from Eqs. (2) and (3) agree by the order of magnitude, the temperature domains where these equations are valid do not overlap. Therefore, to find the behaviour of β in the transition interval from the diffusion-controlled to the GTDGL regime one needs to improve the accuracy of the calculation taking into account both mechanisms of relaxation. This is the problem we address in the present paper. We study the linear response FFC of the sparse vortex lattices in small magnetic field by solving numerically kinetic equations describing non-equilibrium states generated by moving isolated vortices. To find kinetic coefficients and driving terms we use vortex structures calculated self-consistently.
For the diffusion-controlled vortex motion, we calculate the temperature dependence β = β(T ) and compare it with the interpolation curve suggested in earlier works 3,6 . Taking into account the electron-phonon scattering, we demonstrate that it leads to the significant suppression of the FFC at intermediate temperatures, τ ph (T c −T ) ∼ 1, as compared to the estimations obtained from the Eqs. (2,3). Using the electron-phonon relaxation rate τ −1 ph as the fitting parameter we obtain numerically accurate fits to the experimentally measured temperature dependencies of FFC in Zr 3 Rh 13 and Nb-Ta 7,8 .
The structure of this paper is as follows. In Sec.II we introduce the Keldysh-Usadel description of the kinetic processes in dirty superconductors. Here the basic components of the kinetic theory are discussed including kinetic equations, self-consistency equations for the order parameter and the general expression for the viscous force acting on moving vortices. Sec. III introduces θparametrization of the theory. Calculated temperature dependencies of FFC are reported in Sec.IV for different regimes. The diffusion-controlled flux-flow is discussed in Sec.(IV A) and the influence of increasing electronphonon relaxation rate is studied in Sec.(IV B). The work summary is given in Sec. V.

II. KINETIC EQUATIONS AND THE FORCES ACTING ON THE MOVING VORTEX LINE
The quasiclassical Green's function (GF) is defined aš where g K is the (2×2 matrix) Keldysh component and g R(A) are the retarded (advanced) ones. The GFǧ = g(t 1 , t 2 , r) depends on times t 1,2 and a single spatial coordinate r. In dirty superconductorsǧ obeys the Keldysh-Usadel equation Hereτ 0,1,2,3 are Pauli matrices in Nambu space, D is the diffusion constant,Ĥ(r, t) = i∆−ieφτ 0 , where∆(t) = i|∆|τ 2 e −iϕτ3 is the gap operator, ϕ is the gap phase and φ is the electrostatic potential. In Eq.(5) the commutator is defined as [X, g] t = X(t 1 )g(t 1 , t 2 ) − g(t 1 , t 2 )X(t 2 ), similarly for anticommutator {X, g} t = X(t 1 )g(t 1 , t 2 ) + g(t 1 , t 2 )X(t 2 ). The symbolic product operator is given by t 2 ) and the covariant differential superoperator iŝ The collision integral in (5) is given by where the self energyΣ may contain contributions from different relaxation processes. Here we take into account only the electron-phonon scattering which plays an important role in the energy relaxation. The Keldysh-Usadel Eq. (5) is complemented by the normalization condition (ǧ •ǧ)(t 1 , t 2 ) =δ(t 1 − t 2 ) which allows to introduce parametrization of the Keldysh component in terms of the distribution function The deviation of f L from the equilibrium distribution is related to the effective temperature change, and f T is the charge imbalance on the quasiparticle branch.
Then, keeping the first order non-equilibrium terms we obtain the system of two coupled kinetic equation that determine both the transverse and longitudinal distribution function components f L,T = f L,T (ε, t) (the detailed derivation is given in the AppendixA) where the energy-dependent diffusion coefficients D T,L and the spectral charge current j e are given by In Eqs. (10,11,14) we use the covariant time derivative and spatial gradient defined by∂ t =τ 0 ∂ t + 2ieφτ 3 and ∇ = ∇ − ieA[τ 3 , ]. We omit the driving terms containing electric field which is justified in type-II superconductors with large Ginzburg-Landau parameters. In such systems the dominating driving terms are those containing order parameter gradients. The electron-phonon collision integral in the r.h.s. of kinetic equation (11) the components ofǏ are given by Eq.(7) with electronphonon self-energies 15 Here are the free phonon propagators, D RA = D R − D A and g RA =ĝ R −ĝ A . We parametrise the electron-phonon self-energy by dimensionless constant λ ph = (T c τ ph ) −1 where τ ph is electron-phonon relaxation time at T = T c . The force acting on the moving vortex line from the dissipative environment can be calculated according to the expression 3,16 where ν is the density of states andĝ nst is the nonstationary Green's function which can be obtained by the gradient expansion as followŝ Here f T denotes the deviation from the local equilibrium as discussed above. In Eq.(18) we neglect the contribution from the normal component of the charge current. This assumption is well justified for the small magnetic fields as compared to the upper critical one 17 .

III. θ-PARAMETRIZATION
In general, the normalization condition allows one to parametrize GF by complex variables θ andφ. For the axially symmetric vortices the latter coincides with the vortex phaseφ = ϕ. In this case we havê The complex parameter θ = θ(r), depending only on the distance to the vortex centre r is given by the solution of the Usadel equation where ∇ 2 r = ∂ 2 r + r −1 ∂ r , see Appendix C. The boundary conditions for Eq.(22) read where ∆ 0 = |∆(r = ∞)|. Electron-phonon scattering with characteristic time τ in Eqs. (22,23) regularizes spectral functions near the gap edge singularity. At low temperatures electron-phonon scattering does not affect the calculation results. In the vicinity of T c , its value is important since the inelastic relaxation dominates the dissipation. To describe the effects of electron-phonon scattering on the relaxation we calculate τ self-consistently within the relaxation-time approximation described in Appendix B. In this approach To determine the gap profile, we use a stationary selfconsistency equation written in the form Here the summation runs over Matsubara frequencies ω n = (2n + 1)πT , n = 0 . . . ∞, and the angle θ M n (r) parametrizes imaginary-frequency GF obtained by the transformation θ → −iθ M n from the Eq. (20,21) in the upper and lower half-planes, respectively. To obtain θ M n (r) we solve the Usadel equation (22) with θ → −iθ M n and ε → iω n . We assume that the condition ω n τ ≫ 1 is always satisfied and neglect relaxation time correction while solving Eq. 22 for the imaginary frequencies. The boundary conditions read as θ M n (0) = 0 and The driving terms in kinetic equations (10,11) are given by the time-derivatives of the order parameter which for the steady vortex motion can be written as where v L is the vortex velocity. For the axially-symmetric vortex, this form of driving terms allows for the separation of the variables using the ansatz Here the amplitudesf L,T =f L,T (r) are given by the ordinary differential equations which can be written in the compact form as follows where ϑ = Reθ and η = Imθ. For detailed derivation see Appendix C. The last two terms in Eq.(31) describe scattering-out and scattering-in contributions to the inelastic relaxation of the non-equilibrium longitudinal imbalance. The integrals are given by Kinetic equations (30,31) are solved numerically within the interval 0 ≤ r ≤ r c , where r c is the cell radius. For the regime of diffusion-controlled dissipation we choose the interval large enough so that the result is not sensitive to r c . When discussing the crossover to the inelastic relaxation-driven dissipation we set the interval to be larger than the inelastic relaxation length, √ Dτ , which determines the decay off L at large distances. We use the following boundary conditions Here the condition at r = 0 in Eq.(34) follows from the regularity of the solutions at the origin, while condition at r c provides the disappearance of charge imbalance and the absence of the heat flow into the bulk. The viscous friction force acting on individual moving vortex can be written as F env = −̺v L . We present viscosity coefficient in the form ̺ = π ν(α + γ) separating the contributions of the driving terms related to the gap modulus and phase gradients, see Appendix C. In general the flux-flow conductivity can be expressed through the vortex viscosity as follows 4 where φ 0 is the magnetic flux quantum. Taking into account the normal-state Drude conductivity, σ n = 2e 2 νD, we write the FFC in the form (1) with The upper critical field H c2 is determined by the Maki where ψ is digamma function. The low-temperature limit gives H c2 = φ 0 T c /(2Dγ 0 ), where γ 0 = 1.781. Close to the critical temperature one obtains H c2 = φ 0 /(2πξ 2 GL ) and ξ GL = πD /8(T c − T ) is the Ginzburg-Landau correlation length.

A. Diffusion-controlled flux flow
When the temperature is sufficiently far from the critical one the electron-phonon relaxation terms in the kinetic equation (31) can be neglected being much smaller than the diffusion one. Qualitatively this approximation means that the non-equilibrium quasiparticles generated near the vortex can drift to the infinity at the rate exceeding the one of inelastic relaxation. This regime is called the diffusion-controlled flux-flow and it is realized in the temperature domain (T c − T )τ ph ≫ 1. Below we analyse this scenario separately for different temperature intervals.

Low-temperature limit
At low temperatures, the sizeable quasiparticle density exists only inside vortex cores where the superconducting order parameter is suppressed. In this case, it is sufficient to consider only zero-energy GF for which parameter θ is purely imaginary, ϑ = 0. The dissipation is dominated by the charge relaxation processes described by the distribution function f T . The effective temperature change described by the distortion of f L can be neglected. Then Eq. (30) can be written as follows As a result, the coefficients that determine FFC in Eq.(36) are given by Previously, the value of β ≈ 0.9 has been reported by GK 2 . The calculation was based on the approximate vortex solution taken from Ref. 19. This vortex structure was obtained by solving iteratively the self-consistency Eq. (26). Each iteration step was performed as follows. First for a given vortex profile the GFs at each Matsubara frequency were determined by solving the Usadel equation. Then these GFs were substituted to the selfconsistency equation in order to calculate the updated order parameter distribution. The iteration procedure used in Ref. 19 started from the gap function which is known also as the Clem ansatz 20 . However, instead of taking the sufficient number of iterations to reach selfconsistency only a single iteration step was performed in Ref. 19. In this way the approximate vortex profile was obtained which was used later to calculate the flux-flow conductivity at low fields 2 .
To get the correct order parameter distribution we have performed sufficient number of iterations so that to ensure that the order parameter changes with each update become negligible. With the help of the fully self-consistent vortex structure obtained in this way we found that β = 0.77. Previously reported value 0.9 is overestimated by 17%. The disparity between initial gap distribution, the one obtained after the first iteration and the exact gap function together with corresponding values of β are shown in Fig. 1.

High-temperature limit
For elevated temperatures but still in the diffusivecontrolled limit (T c − T )τ ph ≫ 1, the non-equilibrium states are dominated mostly by the change in the number of quasiparticles determined by the f L mode while the charge imbalance f T yields the subdominant contribution. In this regime the FFC was calculated within the local density approximation for the spectral functions 3,5 . This approximation results in the expression (2) with β 0 = 4.04, see Appendix D for the calculation details.
The local density approximation is well justified in the limit T → T c . However, to stay in the diffusion-controlled regime the temperature cannot be taken infinitesimally close to the critical one. Thus it is interesting to improve the accuracy of β calculation for small but finite values of T c − T . For this purpose we find the order parameter solving the self-consistency equation (26) numerically. After that we considered Eq.(22) for the spectral functions, where small parameter 1/τ regularizes the gap edge singularities. We fixed the value λ ph ∼ 10 −6 so that relaxation time τ appears to be sufficiently large and diffusion-controlled FFC remains unaffected up to the temperatures 1 − T /T c ∼ λ ph . By starting with initial distributions for ϑ and η, we calculated relaxation time τ according to Eq. (25) and then solved numerically Eq.(22) to get new functions ϑ and η. By repeating this procedure iteratively, we found spectral functions with sufficient accuracy. By using these solutions, we calculated relaxation rate ν out for non-equilibrium longitudinal imbalance, Eq. (32), and solved the kinetic equations (30)-(31) by omitting scattering-in term j in . Note that the relaxation term in the kinetic equation (31) allows to apply the zero boundary conditions in the bulk. Fig. 2 demonstrates the exactly calculated ϑ, η and the distribution functionf L compared to those obtained within the local-density approximation. With these functions we calculated integrals α and γ, see expression (C7). As a result, we obtained the divergent behaviour (2) with the dimensionless parameter β 0 ≈ 3.7. Therefore, localdensity approximation overestimates β 0 by 9%.

Intermediate temperatures
For the temperatures within the broad range between limiting cases considered above, the contributions of both the f L and f T modes are generically of the same order of magnitude. Therefore in order to calculate the FFC it is necessary to solve system of coupled kinetic equations (30)-(31). This can be done only numerically, and exact temperature dependence β = β(T ) in the diffusioncontrolled regime has never been calculated before. Previously only the interpolation curve between GK and LO results has been suggested 3 . Below we compare this interpolation curve with the result of an exact numerical calculation which is done in the same way as discussed above in Sec. IV A 2 by repeating all steps at different temperatures.
Shown by the red curve in Fig. 3 is the obtained temperature dependence β = β(T ) which is qualitatively similar to the interpolation curve (black line). Both dependencies feature the gradual increase from Bardeen-Stephen limit, β ∼ 1, at small temperatures to the large values of β at high temperatures due to decrease in diffusion relaxation rate. However, calculated dependence β(T ) is significantly lower compared to the interpolation curve suggested previously in Ref. 3.

B. Suppression of FFC by inelastic relaxation
Inelastic electron-phonon scattering provides an additional relaxation mechanism which affects FFC. This relaxation channel plays an important role at temperatures close to the critical one when the spatial gradients of the distribution functions become small due to increase in the correlation length and superconducting energy gap is suppressed.
The crossover between diffusion-controlled and inelastic relaxation-controlled branches of the β(T ) dependence occurs at the temperatures τ ph (T c − T ) ∼ 1 where none of these approximations can be applied. The behaviour of β(T ) in this region of parameters has not been studied before. To analyse an interplay between two relaxation regimes we calculate numerically FFC for different inelastic scattering rates determined by the value of parameter λ ph . We apply the same numeric procedure as was discussed in Sec. IV A 2. Fig. 4 shows the result of the calculation. Inelastic electron-phonon scattering suppresses the maximal value of FFC and smears the crossover from solely diffusiondriven to inelastic relaxation controlled regimes. Such a behaviour is caused by suppressed generation of nonequilibrium quasiparticles due to the presence of electronphonon relaxation channel so that non-equilibrium longitudinal imbalance becomes weaker. This follows from kinetic equation (31) where electron-phonon relaxation tends to suppress the source term determined by the density gradient.
To demonstrate the consistency of our numerics we first estimated the effect of scattering-in contribution to collision integral in kinetic equation (31). To do this we solved kinetic equations without scattering-in term for the energies in the interval [−20∆ 0 . . . 20∆ 0 ] and then calculated collision integral (33). By using its value we solved kinetic equations again and obtained corrected distribution functions together with more accurate values of FFC shown in the left panel of Fig. 4 by dashed curves. The effect of scattering-in term is rather small.
Secondly, we calculated FFC as temperature approaches the critical one. In this limit the GTDGL theory 21 becomes of relevance. By neglecting electric field and charge imbalance mode we found that the remaining contribution to FFC gradually approaches the Tinkham term of the GTDGL model, see Fig. 4. At that, the crossover towards electron-phonon relaxation controlled regime and the GTDGL theory takes place very close to the critical temperature where electron-phonon relaxation rate is at least ten times larger than the one for the diffusion. The opposite condition D/ξ GL = 10T c λ ph is satisfied at the temperature T /T c = 0.6 for λ ph = 0.1 and below this limit FFC is well approximated by diffusion mechanism only, see Fig. 4. This suggests that the temperature interval where FFC is characterized by coexistence of diffusion driven and inelastic scattering controlled mechanisms of relaxation can be quite wide and none of estimations (2,3) can give adequate description in this region.
The modification of FFC caused by electron-phonon scattering allows to obtain good numerical agreement with experimental data by using the inelastic relaxation rate λ ph = (τ ph T c ) −1 as the fitting parameter. In real superconducting systems, FFC is strongly affected by electron-phonon relaxation so that applicability of Eqs.
(2) and (3) appears to be very limited. In this case, the overall temperature behaviour of FFC can be found only numerically due to multi-component mechanism of the non-equilibrium quasiparticle relaxation during the motion of the vortexes.
In Fig. 4 we compare numerically calculated curves with experimental data for Na-Ta system 7,8 and amorphous superconductor Zr 3 Rh 13 . For the former case, good fit is achieved for the value λ ph = 0.1 which corresponds to the electron-phonon relaxation time about 10 −11 s. This value coincide by the order of magnitude with ones reported previously for niobium 22,23 . For the latter system, electron-phonon relaxation time is found to be about 10 −10 s.

V. SUMMARY
To summarise, we have calculated the FFC in diffusive superconductors for small magnetic fields and arbitrary temperatures taking into account the electronphonon relaxation and using the self-consistent vortex solutions. At first, we have obtained the exact value of the dimensionless parameter β = 0.77 which determines the FFC in the low temperature limit T → 0. Second, we calculated the overall temperature dependence of β in the domain of the diffusion-controlled flux flow, that is at τ ph (T c − T ) ≫ 1. Significant deviations from the previously reported interpolation curve are obtained. Finally, we studied the crossover between the diffusioncontrolled and generalized TDGL regimes which occurs at τ ph (T c − T ) ∼ 1. The maximal value of β obtained in this region is much smaller than expected from the estimations based on the Eqs. (2) and (3) at the border of their applicability. Consequently, we obtained significant suppression of FFC near T c by changing the electronphonon relaxation rate and achieved better agreement with the experimental data.

VI. ACKNOWLEGEMENTS
The work was supported by the Estonian Ministry of Education and Research (grant PUTJD141) and the Academy of Finland.
The diagonal elements of matrix equation (A1) give equations forĝ R/A which have same form as (A1) witĥ The non-diagonal element reads as To obtain last relation we substituted parametrization (8) and used the associative property of differential su-peroperator∂ r (g 1 •g 2 ) =∂ r g 1 •g 2 +g 1 •∂ r g 2 . To get rid of the last two terms we subtract the spectral components of the Eq.(A1) to obtain finally the equation In our consideration, collision integral describes electronphonon scattering channel only. This term is responsible for establishing the equilibrium in the system. To proceed we introduce the mixed representation in the time-energy domain as followsĝ(t 1 , t 2 ) = ∞ −∞ĝ (ε, t)e −iε(t1−t2) dε 2π , where t = (t 1 + t 2 )/2. To keep the gauge invariance we introduce the modified GFĝ new =Ŵ (t 1 , t)ĝ(t 1 , t 2 )Ŵ (t, t 2 ) where the link operator is defined in the text. This transformation removes the scalar potential term from the kinetic equations and adds the chemical potential shift eφ. We absorb this shift by substituting f T → f T + eφ∂ ε f 0 , where f 0 (ε) = tanh[ε/(2T )] is equilibrium distribution, so that f T hereafter denotes the deviation from the local equilibrium. Then keeping the first order terms in frequency, we get the gradient approximation where E = −∇φ − ∂ t A is electric field.
Here we assume the first order in deviation from equilibrium so that equilibrium distribution f 0 is substituted in last term in (A5). With the same accuracy we obtain In (A6) we keep only terms which contribute to the kinetic equations.
In the mixed representation the kinetic Eq.(A4) has the following gauge-invariant form Here we took into account only first-order terms in the deviation from equilibrium and introduced gaugecovariant time derivative∂ t =τ 0 ∂ t + 2ieφτ 3 .
To obtain the equations (10) and (11) in the main text we trace Eq. (A7) with Nambu matricesτ 0 andτ 3 respectively. Here we took into account that Tr(ĝ Rτ 3ĝ A ) = 0 because of the relationĝ A = −τ 3ĝ R+τ 3 and the general form of the equilibrium spectral functionĝ R = g 3τ3 +g 2τ2 e −iϕτ3 . Then we neglect the driving terms with electric field and electron-phonon relaxation of the charge imbalance to get the Eq.(10). We keep the electronphonon collision integral in the Eq.(11) which plays an important role in vortex dynamics being the only energy relaxation channel.
After integration, this can be written as F env = −̺v L , where the viscosity coefficient is given by ̺ = π ν(α + γ) and Here we have took into account thatf T and η are even, whilef L and ϑ are odd functions of energy ε.
Appendix D: Derivation of the LO result Following LO 5 , analytical result for diffusion-driven FFC can be obtained by noticing that near T c the diffusion terms in the Usadel equation (22) are much smaller than the gap field. As a result, local density approximation can be implemented, where ϑ and η are determined by their homogeneous expressions with bulk gap substituted by local value of gap field.