# Quasiparticle relaxation in superconducting nanostructures

###### Abstract

We examine energy relaxation of non-equilibrium quasiparticles in “dirty” superconductors with the electron mean free path much shorter than the superconducting coherence length. Relaxation of low-energy non-equilibrium quasiparticles is dominated by phonon emission. We derive the corresponding collision integral and find the quasiparticle relaxation rate. The latter is sensitive to the breaking of time reversal symmetry (TRS) by a magnetic field (or magnetic impurities). As a concrete application of the developed theory, we address quasiparticle trapping by a vortex and a current-biased constriction. We show that trapping of hot quasiparticles may predominantly occur at distances from the vortex core, or the constriction, significantly exceeding the superconducting coherence length.

## I Introduction

Current interest to the dynamics of Bogoliubov quasiparticles in superconductors is motivated in no small part by the efforts aimed at building a quantum computer. The actively explored “cavity QED” architecture relies on quantum coherence of qubits built of conventional superconductors MHDRJS2013 (). Realization of the topological quantum computing requires coherence of devices made of proximitized semiconductor quantum wires brought into the p-wave superconducting state by applied magnetic field Karzig2017 (). In any of the concepts, the presence of quasiparticles is detrimental to the coherence. The Q-factors of the parts comprising the cavity-QED qubit are reduced by quasiparticles. They also are able to “poison” the Majorana states, which are central for the topological quantum computing.

The development of the qubit technology has advanced also the ability to monitor the quasiparticles population and dynamics. Time-resolved measurements performed with the transmon Chen2014 () and fluxonium quibits Pop2014 (); Vool2014 () allowed the experimentalists to measure the rates of quasiparticle trapping by a single vortex in a superconducting strip, to identify minute dissipative currents of quasiparticles across a Josephson junction (thus resolving a longstanding “-problem” Barone ()), and to monitor the spontaneous temporal variations of the quasiparticle density.

Measurements Pop2014 (); Vool2014 () did confirm that at low temperatures (less than of Al) the quasiparticle density, albeit low, far exceeds the equilibrium value. Furthermore, statistics of temporal variations of the density substantially differs from thermal noise. Sources of excess quasiparticles remain unknown, and planting quasiparticle traps Rajauria2009 (); Rajauria2012 (); Riwar2017 () remains a viable way of improving the device performance. Trap is a spatial region with a suppressed value of superconducting gap. Suppression may be achieved, e.g., through the proximity effect, or through local violation of time-reversal invariance (as it naturally happens in and around the core of a vortex). Energy loss in the trap (mostly due to phonon emission) prevents a quasiparticle from exiting into the region with the nominal gap value.

The importance of the quasiparticle energy relaxation in device applications, and the newly acquired ability of precise measurements Chen2014 (); Pop2014 (); Vool2014 () of the quasiparticles dynamics prompts us to revisit the kinetic theory of quasiparticles interacting with phonons in a disordered superconductor. We derive the corresponding collision integrals and relaxation rates which then may be used in sophisticated phenomenological models of quasiparticle diffusion and trapping Riwar2017 (); Ullom2012 ().

In considering the electron-phonon interaction in disordered metals, we follow the seminal works of Tsuneto Tsuneto () and Schmid Schmid (), who established the correct form of the electron-phonon interaction in the limit of short electron mean free path, (here is the phonon wave vector). We incorporate the electron-phonon interaction in the general framework of Keldysh non-linear sigma-model. It allows us to consider on equal footing normal metals and superconductors, and it becomes especially convenient for describing the effect of breaking the time-reversal symmetry (TRS).

The kernel of the collision integral for electrons in normal metal which we find in the unified technique, agrees with the earlier results Reizer (); Yudson () obtained diagrammatically; this kernel depends only on the energy transferred in the collision to a phonon. Considering the Bogoliubov quasiparticles, we are able to cast the result for the collision integral in the conventional terms of the quasiparticle energy distribution functions. In the presence of TRS, the corresponding kernel factorizes on two terms: the normal-state kernel and a combination of the Bogoliubov transformation parameters. Factorization takes place also if TRS is broken; in that case, the second factor is determined by the proper solution of the Usadel equation. In either of the two cases, the second factor depends separately on the initial and final energy of a quasiparticle.

The additional (compared to the normal state) energy dependence of the kernel affects the dependence of the quasiparticle relaxation rate on its energy. These rates, in turn, determine the effectiveness of trapping. As an example of application of the developed theory, we consider trapping of a quasiparticle by an isolated vortex and a current-biased constriction. In both cases there is a pattern of super-currents, slowly decaying as a function of distance, , from the vortex core, or from a constriction with -dimensional superconducting leads. These super-currents lead to a weak breaking of TRS and thus suppression of the energy gap and modification of the energy dependence of the density of states (DOS). Such a suppression allows quasiparticles to be trapped already very far from the vortex core or the constriction. For low enough phonon temperature and relatively “hot” quasiparticles this peripheral shallow trapping proves to be more efficient than the deep trapping by the core of the vortex, or the constriction. We discuss possible relation of the theory to experiments Chen2014 (); Siddiqi2014 ().

The paper is organized as follows: in Section II we review the theory of electron-phonon interactions in disordered normal metals. We derive the corresponding Keldysh non-linear sigma-model and use it to obtain the electron-phonon collision integral in the dirty limit. In Section III we generalize the sigma-model on superconductors, including those with broken TRS, and in Section IV derive kinetic equation for the quasiparticle distribution. Section V is devoted to applications of the theory to trapping by vortex and current-biased constriction as well as discussion of the existing experiments. We summarize with a brief discussion in Section VI. Two Appendices present an alternative derivation of the sigma-model and summarize results for ultrasound attenuation.

## Ii Electron-phonon interactions in disordered normal metals

### ii.1 Interaction vertex

Theory of electron-phonon interactions in normal disordered metals had a long and, at times, controversial history. Early considerations were based on the Fröhlich Hamiltonian Frohlich (), which assumes screened Coulomb interactions between electron density and induced lattice charge, , created by phonon displacement . Here is the uniform lattice charge density. Due to global neutrality it is exactly equal to the electron density:

(1) |

where normal metal DOS is and . While perfectly legitimate in the clean case, the Fröhlich Hamiltonian misses an important piece of the physics in the “dirty” limit , where is phonon wavenumber and is electron elastic mean free path.

As was first realized by Pippard Pippard (), phonons not only deform the lattice, but also displace impurities, transforming formerly static impurity potential into the dynamic one, . Colloquially, this leads to the electron density being dragged along with the lattice displacement and providing a perfect compensation for the induced lattice charge . In other words, the displaced impurity potential provides fast elastic relaxation of the electron distribution around the Fermi surface locally deformed by phonons. These ideas were put on a quantitative basis by Tsuneto Tsuneto () and Schmid Schmid (), who showed that in the limit the Fröhlich Hamiltonian should be substituted by another effective electron-phonon interaction vertex:

(2) |

where and are electrons creation and annihilation operators and is the traceless tensor

(3) |

Notice that, in view of Eq. (1), the last term here represents the Fröhlich coupling . Upon averaging over the Fermi surface it is exactly compensated by the first term in Eq. (3). The remaining coupling is of a quadrupole nature, as seen from Eqs. (2), (3). This leads to a significantly weaker electron-phonon coupling, than the one inferred from the Fröhlich term Altshuler1978 (). A number of subsequent studies Reizer (); Yudson (); Shtyk () reaffirmed validity of the Schmid coupling (2), (3) from various perspectives.

The most straightforward way to derive Eqs. (2), (3) Tsuneto (); Shtyk-thesis () is by performing a unitary transformation, which yields a Hamiltonian in the co-moving reference frame, where the impurity potential is static. We shall nor repeat this derivation here. Instead, we accept Eqs. (2), (3) as a starting point and derive an effective non-linear sigma model which incorporates electron-phonon interaction in the Schmid form. In Appendix A we provide an alternative derivation of the sigma-model, which proceeds in the laboratory reference frame and deals with a dynamic random potential . We show that it brings the same effective sigma model, justifying the use of the effective electron-phonon vertex in the Schmid form (2), (3).

### ii.2 Non-linear sigma model

We now perform the standard Andreev1999 (); Kamenev2011 () averaging over the static disorder and introduce the non-local field to split emerging four-fermion term. The resulting action, including electron-phonon coupling, Eqs. (2), (3), is now quadratic in the fermionic fields which may be integrated out in the usual way, leading to

(4) |

where the inverse bare electron Green function is given by .

From this point on, one proceeds along the standard root of deriving Keldysh non-linear sigma-model Andreev1999 (); Kamenev2011 (). To this end one passes to the Keldysh structure, by splitting the contour on forward and backward branches and performing Keldysh rotation. Upon this procedure the fields acquire the matrix structure, e.g. , where denotes classical and quantum Keldysh components and are the two vertex matrices in the Keldysh space.

One then realizes that the soft diffusive modes of the action are described by the manifold and therefore one can write , where is the Green function in coinciding spatial points,

(5) |

and is a distribution function. Rotation matrices, , belong to an appropriate symmetry group. One then introduces dressed Green function and rewrites the action (4) as

(6) |

Finally, one expands the logarithm here to the lowest non-vanishing orders. This way one obtains the standard non-linear sigma-model action (first neglecting electron-phonon -term):

(7) |

where is the diffusion constant and is the dimensionality of electron system. We focus now on the phonon-induced term. It is easy to see that the first order in term vanishes due to the fact that the integral over the Fermi surface . It is this point, where the Schmid coupling, Eqs. (2), (3), is qualitatively different from the Frölich one (the latter would bring the first order term). Going to the second order in , one finds:

(8) |

We use now

(9) |

along with

(10) | |||

to find

(11) |

where

(12) |

The local vertex (11) is the leading term describing interaction of phonons with the electronic degrees of freedom in disordered metals, in limit. The naive deformation potential term is absent due to the perfect screening manifested in the traceless form of the electron-phonon vertex (3). See also Appendix A for more discussion of this issue. The second order cross-term between the two terms in the logarithm in Eq. (6) leads to . It is of the order of the leading term (11) and thus should not be kept within the accuracy of the adopted approximations.

The effective electron-phonon sigma model, Eqs. (7) and (11), should be supplemented with the standard phonon action. In the Keldysh technique it is given by

(13) |

where is the material mass density, labels longitudinal and transversal polarizations encoded by the projectors

(14) |

and is the acoustic phonons dispersion with the speed of sound . Indexes and Pauli matrix act in the Keldysh space. We will need an imaginary part of the corresponding retarded propagator:

(15) | |||

The corresponding Keldysh component is given by the fluctuation-dissipation relation: , where is the bosonic distribution function.

The effective action, Eqs. (7), (11), and (13) with the vertices defined in Eqs. (12) and (14) serves as the starting point for investigating the kinetics of electrons and phonons. We relegate the evaluation of the ultrasonic attenuation to Appendix B, where we reaffirm the known results Schmid (); Abrikosov (); Shtyk (); Shtyk-thesis () obtained by different techniques, and proceed to study the electron kinetics.

### ii.3 Electron-phonon collision integral

To derive collision integral for electron-phonon interactions one first integrates over the phonon displacements with the help of Eq. (15) to obtain the collision action from Eq. (11):

(16) |

where and summation over are understood. One now looks for the stationary point equation for the action . It’s Keldysh (1,2) component constitutes the kinetic equation Andreev1999 (); Kamenev2011 () for the distribution function in Eq. (5),

(17) |

where the collision integral is given by

(18) |

The variational derivative here ought to be restricted to the sigma-model target space, . A way to insure this is to use parameterization and expand the action (16) to the linear order in . Here ’s are infinitesimal generators of the symmetry transformations. Because of the local nature of the vertex in Eq. (16), the integration may be substituted by the stationary point: . This way one obtains:

(19) | |||

The phonon matrix element in the collision integral (II.3) in normal metal is found to depend only on the energy difference, , where

(20) |

and the superscripts stand for longitudinal and transverse modes, respectively. Expressing the fermion and boson distributions in terms of the respective occupation numbers, and , we may bring to the standard “in minus out” form,

(21) |

where and in equilibrium.

Employing Eqs. (12), (15) and (20) one finds:

(22) |

where is the area of unit sphere and is the effective phonon dimensionality. The coefficients are given ; (in these instances is the dimensionality of the electron system).

These results for the normal metal were derived in Ref. Schmid (); Reizer (); Yudson () using diagrammatic techniques. Here we reproduced them through the sigma-model technique, which is much more suitable for treating the superconducting case, considered below. We notice that of Eq. (22) is factor smaller than that of Ref. Altshuler1978 () for . The latter was obtained with the Fröhlich coupling (i.e. disregarding impurities shifting with the lattice deformations). The two approaches give comparable results for the longitudinal phonons at where they both match with the clean limit expectation . The transversal phonons give the dominant contribution to the collision integral in the disordered limit , since typically . However, in the opposite – clean limit, the transversal matrix element is , Yudson () and is less important than the longitudinal one.

For comparison, the electron-electron collision integral may be written in the form of Eq. (21), with the bosonic occupation number and a different matrix element given by: , Altshuler-Aronov (). As a result the ratio of electron-electron and phonon matrix elements in normal metals is

(23) |

where is the ion mass and we used that and . Therefore electron-electron relaxation in normal metals dominates for the energy transfer .

## Iii Disordered superconductors

### iii.1 Sigma-model

The non-linear sigma-model is readily extended to disordered superconductors Skvortsov (); Kamenev2011 (). It is written in terms of the local pair correlation function , where is the four component spinor in the Nambu and Keldysh subspaces. As a result is a matrix, as well as the matrix in the time, , space. It satisfies the non-linear condition . Its dynamics is governed by the action:

(24) |

where is the order parameter matrix. To discuss broken TRS later on, we have also included a vector potential through the long derivative:

(25) |

Hereafter are Pauli matrices in the Nambu space and . Here the operation involves trace in Nambu-Keldysh space, as well as trace in time (or equivalently energy) space and the spatial integration.

The electron-phonon interactions are given by Eq. (11) (with factor to compensate for the Nambu doubling of the degrees of freedom), where the displacement field is proportional to matrix in the Nambu space. The corresponding collision action, obtained by integrating out the phonon degrees of freedom, is given by Eq. (16) (again with factor ). Its variation over leads to the collision integral in the form of Eq. (II.3). The major difference of the superconducting case is that the matrices in Eq. (II.3) are rotated in Nambu space, as explained below.

Taking variation of the effective action (24), (11) with respect to the as explained after Eq. (18), one obtains the saddle point Usadel equation Usadel (); Kamenev2011 ()

(26) |

We look for a solution of this equation in the standard form respecting causality:

(27) |

with retarded, advanced and Keldysh components being matrices in the Nambu subspace. The non-linear constraint is resolved as

(28) |

where is a distribution matrix in the Nambu space, which may be written as LarkinOvchinnikov () . Here are longitudinal (odd with respect to energy permutation) and transverse (even in energy permutation) components of the quasiparticle distribution function. These two are responsible for the transport of energy and charge correspondingly. (The conventional distribution functions are obtained by Wigner transformation with and .) Since the transversal component usually decays fast to zero, we shall primarily focus only on the long-lived longitudinal component of the non-equilibrium quasiparticle distribution and often omit the superscript for brevity . In thermal equilibrium , while .

The non–linear constraints , Eq. (28), may be explicitly resolved in the Nambu space by the angular parametrization BelzigWilhelm (); Kamenev2011 ():

(29) | |||

where and are complex, coordinate– and energy–dependent angles. Here

(30) |

notice that in presence of the phase the matrix is not a complex conjugate of . The full saddle point -matrix (27), (28) then acquires the form

(31) |

where is the Keldysh space matrix and

(32) |

The expectation value of the order parameter satisfy the self-consistency equation, obtained by variation of over the quantum component . This leads to (we assume ):

(33) |

where is the BCS interaction constant and is the Debye frequency cutoff.

In the absence of the vector potential, i.e. with unbroken TRS, substituting Eqs. (29) into the retarded and advanced components of the Usadel equation (26) one finds for the Nambu angle:

(34) |

For one thus finds that is real and

(35) |

Within the energy gap, , the angle is , with real . For all energies the following symmetry relation holds . The local DOS is expressed through the Nambu angle as

(36) |

where is the step function.

### iii.2 Superconductors with broken TRS

In many cases of the practical interest the vector potential (and hence the phase ) changes slowly on the scale of the superconducting coherence length. In these cases one may disregard the gradient terms in the action (24) and write it as:

(37) |

where is the energy scale associated with the local breaking of TRS. For a vortex , where is distance from the core, and thus , where is the coherence length. For a thin film of width in a parallel magnetic field one finds .

Taking retarded and advanced components of the Usadel equation (26) without gradients (or equivalently, substituting the saddle point ansatz (27) – (III.1) into the action (37) and taking variation over the complex Nambu angles , ), one finds the saddle point condition:

(38) |

It’s solution is depicted in Fig. 1 and admits an important symmetry relation:

(39) |

The local DOS is expressed through the Nambu angle as

(40) |

it is shown in Fig. 2. Within the energy gap, , DOS is zero, i.e. and thus the angle is , with real . This brings . The right hand side of the latter condition reaches maximum at . Substituting this back into Eq. (38) one finds for the energy gap Larkin1965 ()

(41) |

where the last approximate relation holds for . The gap closes at . Immediately above the gap, , DOS takes the form:

(42) |

At it reaches its maximum value and merges with the BCS result at , see Fig. 2.

At the self-consistency relation (33) takes the form:

(43) |

where according to Eq. (38) and the last integral runs along the contour depicted in Fig. 1. Performing the elementary integration one finds Larkin1965 ()

(44) | |||

where is the order parameter at . Since at , the self-consistency condition looses a non-trivial solution at . On the other hand, the gap closes at . Therefore in the narrow range the order parameter is finite, while where is no gap in DOS, Fig. 3. This is the phenomenon of gapless superconductivity.

For one finds from Eq. (44) . The linear in suppression of the order parameter may be also found from Ginzburg-Landau equation for . Notice that this suppression of the order parameter is parametrically weaker than suppression of the gap, Eq. (41). Therefore for the weak breaking of TRS, , one may drop the distinction between and .

## Iv Kinetics of quasiparticles

### iv.1 Kinetic equation

The kinetic equations are given by the Keldysh component of the Usadel equation (26). Employing Wigner representation and projecting onto and Nambu components, one obtains equations for the longitudinal and the transversal distribution functions:

(45) | |||

(46) |

where local DOS is given by Eq. (40) and other parameters are defined as, LarkinOvchinnikov (); BelzigWilhelm (); Kamenev2011 ():

(47) |

(48) |

(49) | |||||

The mass, , exists only in the absence of TRS. It goes to zero at large energy as , but acquires a large value near the gap. Such a mass provides a rapid decay of the transversal component of the distribution function to zero. We thus focus here only on the slow longitudinal relaxation.

The corresponding collision integral is given by Eq. (II.3) (with factor to compensate for Nambu doubling of the degrees of freedom), where one should use the Nambu-rotated -matrices, Eq. (31). This yields, e.g.:

(50) | |||

where and . As a result, the kinetic equation for the quasiparticles occupation number acquires a form:

(51) | |||

where superconducting phonon matrix element is:

(52) | |||||

where , the normal state matrix element is given by Eqs. (20), (22) and DOS is given by Eqs. (40), (42). Motivated by standard TRS notations, we introduced

(53) |

which is only defined for , see Fig. 4. Employing Eqs. (38)–(42), one may show that