Odd-frequency superconductivity near a magnetic impurity in a conventional superconductor

Superconductor-ferromagnetic heterostructures have been suggested as one of the most promising alternatives of realizing odd-frequency superconductivity. In this work we consider the limit of shrinking the ferromagnetic region to the limit of a single impurity embedded in a conventional superconductor, which gives raise to localized Yu-Shiba-Rusinov (YSR) bound states with energies inside the superconducting gap. We demonstrate that all the sufficient ingredients for generating odd-frequency pairing are present at the vicinity of these impurities. We investigate the appearance of all possible pair amplitudes in accordance with the Berezinskii $SP^{\ast}OT^{\ast} = -1$ rule, being the symmetry under the exchange of spin, spatial, orbital (in our case $O=+1$) and time index, respectively. We study the spatial and frequency dependence of of the possible pairing amplitudes, analyzing their evolution with impurity strength and identifying a reciprocity between different symmetries related through impurity scattering. We show that the odd-frequency spin-triplet pairing amplitude dominates at the critical impurity strength, where the YSR states merge at the middle of the gap, while the even components cancel out close to the impurity. We also show that the spin-polarized local density of states exhibits the same spatial and frequency behavior as the odd-$\omega$ spin-triplet component at the critical impurity strength.


I. INTRODUCTION
Odd-frequency (odd-ω) superconducting (SC) pairing is a proposed unconventional dynamic SC state that is both non-local and odd in the relative time coordinate. 1 Berezinskii was the first to point out that the only requirement on the pair correlator is to be antisymmetric under a simultaneous exchange of spin (S), spatial (P * ) and time (T * ) labels, 1,2 which was later extended for multiorbital (multiband) systems 3 to include exchange of the orbital index (O). It is written concisely as SP * OT * = −1 and referred to as Berezinskii rule. By allowing for odd time (or, equivalently, frequency) dependence, the two possible symmetries (spinsinglet even-parity and spin-triplet odd-parity) were extended to two more (spin-singlet odd-parity and spintriplet even-parity). 4 The odd-ω spin-triplet and evenparity order was initially proposed in liquid 3 He and referred to as Berzinskii state. It was eventually ruled out in favor of the even-frequency, spin-triplet, odd-parity. 5,6 This idea drove the interest on realizing odd-ω states in solids. Later, it has been shown that phonon-mediated electron-electron interactions cannot stabilize an odd-ω SC order parameter (OP), 7 , but other mechanisms based on spin-fluctuation-mediated interactions dependent on the Cooper pairs spin may be favorable. Several systems, such as disordered systems, [8][9][10] heavy fermion and Kondo systems, [11][12][13][14] and, more recently, Dirac semimetals 15 were contemplated as theoretical possibilities of stabilizing an odd-ω SC OP. These materials were shown to exhibit exotic electromagnetic properties. 16,17 A different approach to inducing odd-ω pair correlations is through engineering heterostructures where a conventional (even-frequency, spin-singlet, s-wave) SC is proximitized to region that causes the breaking of some of its symmetries. This generates a corresponding odd-ω correlation component (see Sec. I. C. in Ref. 4). Historically the first and the most analyzed systems are the superconductor-ferromagnet (SF) heterojunctions where the odd-ω spin-triplet s-wave pair correlations survive disorder. 18,19 Issues concerning various geometries, magnetization profiles and the effects on different proximity and inverse proximity induced orders were considered (for a recent review see 20 ). More recently, the interest in these heterojunctions was revived in relation to their potential application in superconducting spintronic devices. [21][22][23][24] Conceptually, the sufficient ingredients for generating odd-ω pair correlations by proximity effect may be preserved while transitioning from the geometry of heterostrutures, through finite-sized ferromagnetic islands to the limit of an isolated magnetic impurity, as illustrated in Fig. 1. In this work we consider the extreme case of reducing the size of the ferromagnetic island, which corresponds to the controlled immersion of a single magnetic impurity atom in a clean superconductor, panel (c) in Fig. 1.
In this work, we consider the limiting situation of magnetic defects (Fig. 1(c)), namely an isolated magnetic impurity atom immersed in a bulk conventional SC. It is explicitly demonstrated that odd-ω pair correlations are generated, in accordance with the SP * OT * rule.Very recently, a complementary study has demonstrated the existence of odd-ω in non-magnetic potential impurities due to the renormalization of the SC OP. This corresponds to the limiting case of a metal-superconductor (NS) heterojunctions (where the N region has been reduced to an impurity), shown to to host odd-ω pair correlations much later than the already mentioned cases of SF heterojunctions when the spatial modulation of the SC OP is taken into account. [48][49][50] .
The simple geometry under consideration offers the advantage of exactly decomposing the spin, spatial and frequency dependence of the pair correlations, enabling a clear demonstration of the conversion of spin and time symmetry. Motivated by the experimental capabilities in measuring magnetic impurities in superconductors, 51,52 it is worthwhile to point out the features of the odd-ω Berezinskii pairing that are correlated with spectroscopic data. [53][54][55][56] We discuss the behavior of the spin-polarized (SP) local density of states (LDOS) as it can be related to the pair correlations.
The rest of the article is organized as follows: In Sec. II we describe the model Hamiltonian of the system formed by a Kondo-type impurity immersed in a BCS bulk SC. We introduce the Green's function and the Dyson equation in Sec. II A. The results are presented in Sec. III. In Sec. III A, we provide analytical expressions for the impurity Green's function. In Sec. III B, the reciprocity relations between the pair correlation components are introduced, discussing the parameters where odd-ω components dominate. In Sec. III C we show the SP LDOS close to the impurity and local magnetization DOS. In Sec. III D, III E the pair correlations are represented. Finally, analytic expressions and derivation details are given in the Appendices.

II. MODEL AND METHOD
We consider a clean metal with spin-degenerate and isotropic dispersion ξ k ≈ v F (|k| − k F ) in D spatial dimensions whose Fermi surface is a sphere of radius k F and the band velocity at the Fermi level (Fermi velocity) is v F . The metal is an intrinsic conventional SC, with a spatially uniform 57 singlet s-wave SC OP ∆. The Hamil-tonian of the system is given by H = H BCS +H imp , where the metal is a BCS mean-field Hamiltonian: Here, c kσ (c † kσ ) is the second-quantized electron annihilation (creation) operator that destroys (creates) an electron with wavevector k and spin projection σ ∈ {↑, ↓}.
Inside the metal, a point-like magnetic impurity at position x = 0 is immersed. The impurity Hamiltonian is a Kondo-type Hamiltonian described by where U m is the impurity strength and n is a unit vector ((n · n) = 1) that determines the polarization of the magnetic impurity. M(x) is the real-space spin polarization due to the conduction electrons: with i = (x, y, z) being a Cartesian component label.
Here, c † σ (x) is a creation operator that creates an electron at position x with spin polarization σ. It is related to the momentum-space via a Fourier transform: and analogously for c σ (x). The minus sign in the exponent is in accordance with the choice for the translation operator T (a) = exp [−i (a · P)] (we choose units in which = k B = 1) so that T (a) c † σ (x) T † (a) = c † σ (x + a). The normalization L −D/2 ensures that the real-space second-quantized operators satisfy the anticommutation relation: We define the Nambu spinor: which satisfies the commutation relations with the spin operator: where the Pauli matricesτ µ andσ µ operate in the Nambu and the spin sub-spaces, respectively. We define an imaginary-time-ordered Green's function matrix: whereˇdenotes a 4 × 4 matrix in the Nambu-spin space. The imaginary-time Heisenberg operators evolve with the total Hamiltonian of the system H according to: In case of a time-independent Hamiltonian, the Green's function Eq. (8) is only a function of the time difference τ − τ and may be expanded in fermionic Matsubara frequencies: Then, through an analytic continuation i ω m → ω, the Green's function becomes a function of complex frequency.
In the matrix structure adopted in Eq. Eq. (6), the SP * OT * = −1 rule has the form: or, equivalently, as a function of Matsubara frequency: 58 (11b) Therefore, Eq. (11b) implies for the singlet G eh s , and triplet components G eh t : the following symmetry properties under P * and T * : This in accordance with the fact that spin-singlet pair correlations are odd under exchange of spin indices (S = −1), while spin-triplet are even (S = +1).

A. Dyson equation
The Green's function defined in Eq. (10) in the presence of impurity potential Eq. (2) satisfies the following Here,Ǧ cl (x, x ; ω) is the Green's function matrix for the clean system which is invariant under spatial translations:Ǧ andǦ cl (k; ω) is assumed to correspond to a conventional s-wave spin-singlet superconductor with a SC OP ∆ and an isotropic normal dispersion ξ k invariant under spin rotations. This means thatǦ The integral over k in Eq. (14) is evaluated in the wide band limit in Appendix A and the result is expressible in terms of the functions φ c/s D (r, ω). Due to the translation invariance and isotropy, these functions only depend on r = |x − x |. The ω-dependence enters through while the imaginary part is odd. The real-space Green's function is given by: where N F is the DOS per spin component and its dependence on the model parameters is given by Eq. (A2). Note that (φ c ) D (0; ω) = 1, and (φ s ) D (0; ω) = 0. Finally, the solution of Eq. (12a) is of the form: so that all the effects of the magnetic impurity potential are expressed in terms of the impurity Green's functioň G imp (x, x ; ω). This function lacks translation symmetry due to the presence of a scattering center. Nevertheless, it has an analytic expression that may be obtained within the T -matrix approach. Details of its evaluation are given in Appendix B. The T -matrix has a denominator that is a quadratic function of ω, giving two poles at ±ε 0 , corresponding to the energy of the YSR states. They are given by: 25 with the dimensionless impurity strength:

III. RESULTS
In this section, we present results that are directly related to the pair correlations around the impurity and experimentally measurable single-particle properties.

A. Analytic expression for the impurity Green's function
We begin our analysis deriving the expression for the impurity Green's function using the T-matrix formalism, described Sec. B. From Eqs. (B1a), (B2) and (16a), we have the following explicit expression for the impurity correction to the Green's function matrix: where is a denominator with poles at the YSR bound-state energies ±ε 0 , Eq. (18), andĝ s (x, x ; ω) andĝ t (x, x ; ω) denote the singlet and triplet component of the Green function. We note that the triplet component is aligned with the magnetization direction of the impurity. The singlet and triplet contributions can be decomposed in a basis of 4 functions with definite parity under the exchange x ↔ x : The first 3 functions are even, while the last is odd under coordinate exchange. All of the functions, according to Eqs. (A5), are even in frequency. Analytic expressions and the spatial behavior of these functions are investigated in more detail in Appendix C for different D dimensions. Using these basis functions, we may writeĝ α (x, 2Jm ω =τ 3∆ . (23d)

B. Reciprocal relations between singlet and triplet pair correlations
The SP * OT * = −1 condition for single-band systems, where O = +1, goes over to SP * T * = −1. Thus, there are 4 possibilities of allowed combined symmetries (S, P * , T * ). The enumeration of these possibilities, together with their labeling and the corresponding contributions from Eqs. (16) and (23) are given in Table I. The eh (ω), the situation is reversed. Having in mind the general matrix structure ofĝ t a andĝ s a given in Eqs. (23), we see that for any pair of coordinates x, x , there is a general reciprocity relation between the singlet and triplet components with different parity in frequency ω. We point-out, however, that the basis functions Eqs. (22) also carry (even) frequency dependence, so, the total frequency dependence is coordinate dependent.
Another reciprocity between the electron-hole components follows by noticing that the same ratio holds: g s a (ω) = r I (ω)ĝ t a (ω), (a = 1, 4), for two sets of basis functions with opposite parity under coordinate exchange. We note that at ω → 0 the ratio vanishes, implying thatĝ s a (ω = 0) = 0 (a = 1, 4). In the same way, at the critical impurity strength, J m = 1, g t a (ω) = 0 (a = 1, 4). Taking the electron-hole components, we may rewrite the ratio in Eq. (24a) as: Thus, there is a reciprocity between exchange of the parity with respect to coordinate exchange (a = 1 and a = 4) and ω, while keeping the same spin symmetry. The ratio of these pair correlations is independent of the impurity strength J m . A similar relation can be derived for the remaining two basis functions with same parity under coordinate exchangeĝ s a (ω) = r II (ω)ĝ t a (ω), (a = 2, 3), We note that [g α 3 ] eh (ω) = 0, which may be interpreted as a selection rule, which guarantees that there is no electron-hole component generated by impurity scattering if neither the spin, nor the parity symmetry is changed.

C. Spin-polarized local density of states
We begin plotting local single-particle properties, such as the SP LDOS. The LDOS ν(x, ω), and the LMDOS m(x, ω) are evaluated as: The magnetization has the same direction as the impurity spin polarization. The spin-up (σ = +1) and spin-down (σ = −1) components of the SP LDOS are evaluated as ν σ = (ν + σ |m|)/2, and can be experimentally measured with a magnetic tip STM in spectroscopic mode. From them, working backwards, one may determine m. When representing intensity plots of a quantity f (SP LDOS, LMDOS, real or imaginary part of a particular pair correlation) containing narrow high-intensity peaks around the YSR bound-state energies or the gap edges, and having a smooth background, we use a scaling factor s according to the sigmoid scaling function: This is an odd-valued (s(−f ) = −s(f )), smooth singleparameter function, where the parameter f 0 roughly corresponds to a region of values for f for which the scaling is linear. This scaling has the advantage that it collapses the whole real line on a finite segment [−1, 1], while keeping a linear resolution for low magnitudes, which is particularly convenient in plots with narrow isolated peaks and a varying background. In Fig. 2, we display the spin-up LDOS (top row, panels (a), (d), (g)), the spin-down LDOS (middle row, panels (b), (e), (h)) and the LMDOS (bottom row, panels (c), (f), (i)) as function of the frequency ω/∆ and the position |x|/ξ away from the impurity.
For illustration purposes, we take a short SC coherence length ξ/λ F = 3, which implies a value for the SC OP in units of the Fermi energy (ε F ), ∆/ε F = λ F /(2 π ξ) = 0.053. Also, we evaluate the Green's function at ω → ω + iη where η/∆ = 0.0005 in Eqs. (26) to avoid divergencies. Each column in Fig. 2  explicitly verified by comparing the middle column panels with analogous plots for J m = 2.0. For a clean system, J m = 0 (left column), there is no spin polarization and the LDOS is homogeneous with BCS coherence peaks at ±|∆|. As the impurity strength is increased (middle column), two sub-gap SP peaks appear at ±ε 0 : YSR bound states. Below the critical impurity strength (J m < 1), the state at −ε 0 has the same spin polarization as the impurity, while the one at ε 0 has the opposite polarization. The spin polarization of these peaks does not oscillate with position. At the same time, the above-gap continuum shows a small spin polarization that has the same direction as the impurity for both positive and negative frequencies, but oscillates with position according to the oscillations of the basis functions ( Fig. 9 of Appendix C).
At critical impurity strength, J m = 1 (right column), the two sub-gap peaks merge at zero-frequency. Remarkably, they do not cancel, but, instead display oscillation with position. As it will be shown in the following, similar features are also observed in the spin-triplet odd-ω local pair correlation. Finally, for J m > 1, the YSR states cross and their polarization is reversed. Guided by the pair correlations of the conventional SC OP, which are on-site, we first consider local (x = x ) pair correlations in the vicinity of the impurity. For this choice, the symmetry under exchange of coordinates is necessarily even P * = +1, so there are only 2 possible components consistent with the SP * T * = −1 rule, namely even-ω spin-singlet G eh (−,+,+) (x, x, ω) and odd-ω spin-triplet G eh (+,+,−) (x, x, ω). The conventional correlation belongs to the contribution G eh (−,+,+) , displayed in Fig. 3. In the clean limit (panels (a) and (b)), the system is uniform and there is no x dependence. The correlation is real and even-ω for sub-gap frequencies (−|∆| < ω < |∆|), while it becomes imaginary and odd for above-gap frequencies, as expected from a conventional superconductor. As the impurity strength is increased from J m = 0 (panels (c), (d), (e), and (f)), the features above the bulk gap persist. At intermediate impurity strengths, the build-up of YSR states at ±ε 0 (for J m = 0.5 ε 0 = 0.6 |∆|, in panels (c) and (d)) leads to additional features inside the superconducting gap. The real part of the (−, +, +) component changes sign on both sides of the peak while the imaginary part has the same sign. This feature persists for the real and the imaginary part of any correlation function around the YSR, Figs. 4-8. Further we note that the parity of the real and imaginary parts of any correlation are opposite. 58 The real part, being even in frequency, has the same sign between the YSR bound state energies (|ω| < ε 0 ) and changes sign as the frequency crosses the YSR bound state energies, vanishing at the gap edges (ε 0 < |ω| < ∆). The imaginary part, being odd in frequency, has sharp peaks around the YSR bound-state energies with opposite sign and spatial modulation determined by that of Fig. 9. As the strength of the impurity approaches the critical value, J m = 1, the two peaks with opposite sign merge at zero energy and cancel out. Therefore, the conventional superconducting correlations, (−, +, +), are suppressed for the sub-gap frequency region and in the proximity to the impurity. For J m > 1 the sub-gap features appear again, showing the same symmetry as in the J m = 0.5 case, with a global sign in the imaginary part.
The odd-ω on-site component G eh (+,+,−) is displayed in Fig. 4. Here we point out only the differences with Fig. 3.
We note that the correlation is induced by the effect of the impurity (only from G t imp eh ) and, therefore, completely vanishes for a clean superconductor (J m = 0). At finite J m , we observe the appearance of signatures outside and inside the gap. This component is triplet with spin polarization in the direction of the impurity, as illustrated by Eq. (20). As shown in Fig. 4, the real part is odd (and the imaginary part is even) in ω. Differently from the previously commented case, as the peaks merge (panel Fig. 4(f)), they do not cancel, but, enhance instead. The same enhancement persists also in the real part. The spatial modulation is identical to the one of Fig. 3, which is in agreement with the reciprocity of the two matrices having the same spatial dependence. One final feature is the small contribution of the continuum states exhibiting the same spatial modulation than the YSR.
It is important to note that at the critical impurity strength the conventional (−, +, +) pair correlation is completely suppressed in the vicinity of the impurity, as shown in Fig. 3 (e) and (f). Therefore only local component that survives is the (+, +, −), which corresponds to the Berezinskii pairing, 4 (panels (e), and (f)). Another possibility to consider is the non-local pair correlations. Because the impurity site is preferential, we are interested in fixing one of the coordinates at the impurity site x = 0. In Figs. 5-8, we show the electronhole separation, when a hole is created at the impurity site and an electron is annihilated at a site x. Fig. 5 displays the pair correlation with the same symmetry as the one in Fig. 3. In the absence of an impurity (panels (a) and (b)), it displays a characteristic periodicity with a period ξ. The spatial modulation is according to the different behavior of the basis functions for this choice of coordinates (Fig. 10). Unlike the local correlations, where the imaginary part is zero for sub-gap frequencies, in this case there is a finite component that decays with separation over a length scale determined by 1/ Re (κ(ω)), introduced in relation to Eqs. (16). As expected, the frequency dependence, including the parity and the behavior around the YSR bound state energies is the same as in the local case, since it is determined by the basis functions. Thus, at critical impurity strength, J m = 1, (Fig. 5(e), (f)) the imaginary and real parts vanish in accordance with the exact cancellation of the local correlation at the impurity site.

E. Non-local pair correlations spectra
An analogous comparison may be made for the nonlocal (Fig. 6) and local (Fig. 4) correlations with symmetry (+, +, −). It exhibits the same spatial dependence as Fig. 5. As in Fig. 4, the real part of the correlation is odd in frequency while the imaginary one is even. Simi-larly to the local case, the (+, +, −) correlation does not vanish for J m = 1, but increased instead. It means that, close to the critical impurity strength, the odd correlations dominate.
In addition to the two components discussed before, there are two more correlation components appearing in the non-local case. According to Table I,  As shown, both of the odd-parity non-local pair correlations require a non-zero impurity scattering, as illustrated by the vanishing of Fig. 7(a), (b) and Fig. 8(a), (b). For J m = 1, both components exhibit a finite value, where the symmetry of the real part is set by the T * value. For J m = 1, the (+,-,+) component vanishes while the (-,-,-), odd-ω gets enhanced. By comparing the behavior at critical impurity strength on Fig. 7(e), (f), and Fig. 5(e), (f), we conclude that there exist another type of reciprocity. In this case the parity and spin symmetry are exchanged, while keeping the frequency symmetry unchanged, as hinted in Sec. III B. We note that the the odd-parity pair correlations have spatial dependence determined by the only odd basis function ψ 4 , which has a different spatial behavior than the even ones (Fig. 10), with no change in sign in both the real and imaginary part for sub-gap frequencies.

IV. CONCLUSIONS
To conclude, we have demonstrated that the Yu-Shiba-Rusinov (YSR) bound states formed at the vicinity of magnetic impurities in conventional superconductors harbor spin-triplet, odd-ω pair correlations, which may be referred to as Berezinskii-YSR (BYSR) states. In case of spin-rotation invariant bulk SC, the spin polarization of the BYSR state is collinear with the impurity spin. We have demonstrated that the odd-frequency pairing dominates close to the impurity at the critical coupling strength, where YSR states merge at the middle of the gap. This leads to a regime where odd-frequency correlations can be probed experimentally. We note that the local odd-frequency component exhibits the same frequency and spatial dependence than the local magnetization density of states, which can be measured using spinpolarized scanning tunneling spectroscopy. By solving exactly for the Green's function, we were able to identify a complete reciprocity as a function of the impurity coupling strength under exchange of the S and T * symmetries and with the same parity coordinate, P * . Finally, we have also demonstrated the presence of all the allowed non-local pair correlations consistent with the SP * OT * = −1 rule and summarized in Table I The real space Green's function is obtained by inverting Eq. (15a) and performing the momentum integral according to Eq. (14). When performing the momentum integral the following steps are taken: Here, dΩk stands for integration over the directionsk with Ω D being the solid angle in D dimensions. The angle θ is the angle between r and k. In step Eq. (A1a), the summation over momenta was converted to an integration over a continuous variable not restricted to the First Brillouin Zone, using D-dimensional spherical coordinates. This approximation is valid as long as the Fermi wavevector is far away from the Brillouin Zone boundaries. In step A1b, the linear approximation of the dispersion relation ξ k ≈ v F (k − k F ) was employed. Here, the Fermi momentum k F is defined as ξ kF = 0 and the Fermi velocity is the derivative of the dispersion relation evaluated at k F , v F = dξ k /dk| k=k F . This approximation is valid as long as the relevant energy scales δε that enter in the functionĜ cl (such as ∆ and ω), as well as in the exponential (v F /r) satisfy the constraint |δε| |m * | v 2 F /2. Here, m * is the band mass (m * ) −1 = d 2 ξ k /dk 2 | k=kF . For a parabolic band, the right-hand side of this inequality is simply the Fermi energy ε F , i.e. the energy difference between the bottom of the dispersion relation at k = 0 and the Fermi level at k = k F . For a linear band ε F = v F k F . So, the linear approximation is valid as long as the energy scales of the problem are much lower than the Fermi energy. In step Eq. (A1c), in accordance with the above assumptions, the Fermi energy is taken to infinity. To make the integrals convergent, any high powers of ξ were neglected as well, essentially approximating the DOS per spin component as being flat with energy, and equal to its value at the Fermi level: where λ F = 2π/k F is the Fermi wavelength. This approximation is in agreement with the above quasi-continuum approximation as long as k F π/a, with a being the lattice constant. All these approximations are referred to as the wide-band approximation.
Also, having in mind the non-local pair correlations presented in Sec. III E, we plot in Fig. 10 the basis functions Eqs. (22) for x = 0 as a function of r = x for the same choice of frequencies and parameters as in Fig. 9. Similarly to the local case, the imaginary part vanishes for sub-gap frequencies and different from zero only above the superconducting gap. We note that the overall envelope decay is much slower, because these functions behave as the first powers of φ c/s 3 (r; ω) given in Table II.