Currents induced by magnetic impurities in superconductors with spin-orbit coupling
Sergey S. Pershoguba, Kristofer Björnson, Annica M. Black-Schaffer, Alexander V. Balatsky
CCurrents Induced by Magnetic Impurities in Superconductorswith Spin-Orbit Coupling
Sergey S. Pershoguba , Kristofer Bj¨ornson , Annica M. Black-Schaffer , and Alexander V. Balatsky , Nordita, Center for Quantum Materials, KTH Royal Institute of Technology,and Stockholm University, Roslagstullsbacken 23, S-106 91 Stockholm, Sweden Department of Physics and Astronomy, Uppsala University, Box 516, S-751 20 Uppsala, Sweden and Institute for Materials Science, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA (Dated: October 8, 2018)We show that superconducting currents are generated around magnetic impurities and ferro-magnetic islands proximity coupled to superconductors with finite spin-orbit coupling. Using theGinzburg-Landau theory, T-matrix calculation, as well as self-consistent numerical simulation on alattice, we find a strong dependence of the current on the direction and magnitude of the magneticmoment. We establish that in the case of point magnetic impurities, the current is carried by theinduced Yu-Shiba-Rusinov (YSR) subgap states. In the vicinity of the phase transition, where theYSR states cross at zero energy, the current increases dramatically. Furthermore, we show that thecurrents are orthogonal to the local spin polarization and, thus, can be probed by measuring thespin-polarized local density of states.
PACS numbers: 71.55.Ak, 73.23.Ra, 74.81.Bd
Superconductor-ferromagnet heterostructures were re-cently proposed as a viable platform for realizing topo-logical superconductivity (TS) [1–3], which can host Ma-jorana fermion quasiparticles at vortex cores and bound-aries [4–6]. The Majorana fermions obey non-Abelianstatistics and may be utilized for topological quantumcomputation [7–9]. The key ingredients driving thesesystems into the topologically nontrivial regime are thespin-orbit coupling (SOC) and magnetism. Recently, thesearch for experimental realizations of TS has also led toengineering the Yu-Shiba-Rusinov (YSR) [10–12] bandsinduced by magnetic atoms on the surface of a supercon-ductor [13–25]. Following this recipe, zero-energy peaksin the tunneling spectrum were recently measured at theends of a one-dimensional (1D) chain of magnetic atoms[26]. Such a tunneling spectrum could be the evidence ofthe Majorana edge states, although alternative explana-tions are also possible [27].The interplay of SOC and magnetism has another re-markable consequence. Consider a two-dimensional (2D)surface of a three-dimensional (3D) material. The effec-tive Hamiltonian of the surface h ( p ) = p m + λ ( σ × p ) z contains the Rashba SOC due to the absence of the inver-sion symmetry at the surface. Then, the velocity opera-tor v = dh ( p ) d p = p m + λ ˆ z × σ contains a spin-dependentterm that gives an extra contribution to the current j extra = λ ˆ z × (cid:104) σ (cid:105) . (1)A ferromagnet proximity coupled to the superconductorwould render a finite spin polarization (cid:104) σ (cid:105) (cid:54) = 0 and, thus,generate a current as schematically shown in Fig. 1(a).The phenomenon of driving a current with magnetismis known as the magnetoelectric effect. This effect mayvanish in metals due to dissipation but survives in su-perconductors lacking inversion symmetry [28–32]. Themagnetoelectric effect was also recently discussed in apure 1D model of TS [33]. (a)(b) (c)FIG. 1. (color online) (a) Schematic representation of the non-local currents (red arrows) induced by a ferromagnetic (FM)island on the surface of a superconductor with the RashbaSOC. GL solutions for the current around a circular ferro-magnetic island (gray area) with S = S ˆ z (b) and S = S ˆ x (c). In this work, we show that the magnetoelectric currentis universally generated around single magnetic impuri-ties and ferromagnetic islands, which have been recentlystudied in the context of TS [13–26]. More specifically, wefirst derive the extra terms in the Ginzburg-Landau (GL)free energy corresponding to Eq. (1). For a small ferro-magnetic island on a superconductor with SOC, we find astrong dependence of the current on the relative orienta-tion of the ferromagnetic moment. The current circulatesaround the ferromagnetic island and is short ranged forthe ferromagnetic moment normal to the surface. Onthe other hand, the current has a dipolar power law de-cay for the ferromagnetic moment parallel to the surface. a r X i v : . [ c ond - m a t . s t r- e l ] S e p Next, we discuss the current generated around a pointmagnetic impurity and show that the current is carriedby the impurity-induced YSR states. We also performa self-consistent numerical calculation and find a strongnonmonotonic dependence of the current on the strengthof the ferromagnetic moment. The current strongly peaksat the phase transition, where the YSR states cross zeroenergy E = 0. We further demonstrate that the cur-rent can be mapped by measuring the spin-polarized lo-cal density of states (SP LDOS), which, thus, providesa probe of both the current and the phase transition.Our findings are, therefore, highly relevant for the ongo-ing search of the Majorana bound states in ferromagneticchains [13–26]. Ginzburg-Landau treatment.-
We start by consideringa ferromagnetic island deposited on a 2D surface of a con-ventional s -wave superconductor with the Rashba SOCas illustrated in Fig. 1(a) and described by the Hamilto-nian H = 12 (cid:90) d r Ψ † ( r ) [ h ( p ) τ z + ∆ τ x − S ( r ) · σ ] Ψ( r ) ,h ( p ) = p m + λ ( σ × p ) z − µ, p = − i ( ∇ x , ∇ y ) . (2)Here Ψ = ( ψ ↑ , ψ ↓ , ψ †↓ , − ψ †↑ ) T is a four component spinor, σ and τ are the Pauli matrices acting in the spin andparticle-hole Nambu space, ∆ is the superconductinggap, and we set e = (cid:126) = 1. The ferromagnet and itscoupling to the superconductor are described by the spa-tially dependent vector S = ( S x , S y , S z ). An intuitiveand qualitatively correct picture of the currents can bederived using the GL free energy F = (cid:90) d r (cid:104) n s m A + α ( ˆ z × S ) · A + β ( ∇ S z × ˆ z ) · A (cid:105) , A = ( A x , A y ) = A + ∇ θ , (3)which is valid at length scales larger than the supercon-ducting coherence length ξ sc . In the first term propor-tional to the superfluid density n s , vector A encapsu-lates both the superconducting phase θ and the vectorpotential A . The second and third terms describe thecoupling between the Rashba SOC and magnetism andare derived in the appendix. For example, in the limit p F λ (cid:29) mλ (cid:29) ∆, the coefficients are [34] α = mλ/ π , β = m λ / πp F , and, thus, only present at finite SOC,i.e. when λ (cid:54) = 0. The term proportional to α , known asthe magnetoelectric term [28–30, 32, 35, 36], is allowedonly in the absence of inversion symmetry.Within the above framework, we now discuss the cur-rents induced by a ferromagnetic island of a uniformdisc geometry, which we model as S ( r ) = S θ H ( R − r ),where θ H ( z ) is the Heaviside theta function. We find thecurrent from Eq. (3) as j = δFδ A (cid:12)(cid:12)(cid:12)(cid:12) A =0 = n s m ∇ θ + α ( ˆ z × S ) + β ( ∇ S z × ˆ z ) . (4) First consider the an out-of-plane ferromagnetic moment S = S ˆ z and θ constant. Then the current is givenby the last term in Eq. (4). The current is localizednear the boundary as j ( r ) = − βS (ˆ r × ˆ z ) δ ( r − R ) andcirculates around the ferromagnetic island as shown inFig. 1(b). Since the GL equations are valid at r > ξ sc ,the δ function in the current solution is artificially broad-ened to a scale of the superconducting coherence length ξ sc for visualization purposes. For an in-plane moment S = S ˆ x , both of the first two terms in Eq. (4) arenonzero. The contribution given by the term α ( ˆ z × S )is constant over the region covered by the ferromagneticregion and discontinuous at the boundary. However, thefirst term n s m ∇ θ fixes this discontinuity. Indeed, the vari-ation of the free energy over θ gives the continuity equa-tion: 0 = ∇ · j = n s m ∇ θ + α ∇ · ( ˆ z × S ). The lastexpression is the 2D Poisson equation with a source termthat we solve for θ and plot the currents in Fig. 1(c), seeappendix for more details. The current is constant overthe region covered by the ferromagnet, i.e. j ( r ) = d R for r < R , and has a dipolar profile outside of it, i.e. j ( r ) = r ( d · r ) r − d r for r > R . Here the effective “dipole”moment is defined d = αR ( ˆ z × S ). We also note that ifthe coefficients α and β are large, vortex solutions for thesuperconducting phase θ are favored by the free energyexpression Eq. (3). Microscopic calculation.-
To complement the aboveGL analysis, we also study microscopically the cur-rents generated around a single point magnetic impurity,i.e. we set S = S δ ( r ) in the Hamiltonian Eq. (2). In con-trast to the GL approach, the Green’s function method,used below, allows to study effects to infinite order in S and also at distances smaller than the superconduct-ing coherence length, i.e. for r (cid:28) ξ sc . We evaluate theGreen’s function of the superconductor in the T-matrixapproximation G rr (cid:48) ( ω ) = g rr (cid:48) ( ω ) + g r ( ω ) T ( ω ) g r (cid:48) ( ω ) , (5) T ( ω ) = − S · σ S · σ g ( ω ) . (6)The Green’s function of a clean superconductor in realspace at r is (for r (cid:28) ξ sc ) g r ( ω ) = − π ω + ∆ τ x √ ∆ − ω [ f ( r ) + i ( σ × r ) f ( r )] , (7)where f ( r ) = (cid:2) ρ + J ( p + F r ) + ρ − J ( p − F r ) (cid:3) and f ( r ) = r (cid:2) ρ + J ( p + F r ) − ρ − J ( p − F r ) (cid:3) , J and J are Bessel func-tions, p ± F ≈ p F ∓ λm are the Fermi momenta of the spin-polarized Rashba bands, and ρ ± ≈ ρ p ± F p F are the corre-sponding density of states ( ρ = m/π ). Equation (7) iscalculated with the assumption µ (cid:29) ∆ >
0. The second,spin-dependent, term in Eq. (7) is a consequence of theRashba SOC and vanishes if λ = 0. The poles of theT-matrix give the energies of the impurity-induced YSR (a) (b) (c) (d)FIG. 2. (color online) (a),(c): Currents around a point magnetic impurity calculated using the T-matrix approach. Theimpurity is located at r = 0, and the direction of the magnetic moment S is shown in the inset. The vectors in green indicatethe direction of the in-plane spin polarization determined from SP LDOS. The current and the spin polarization are orthogonal,which is consistent with Eq. (1). In order to enhance the figures, the currents are plotted away from the impurity, i.e. for p F r > .
5. (b),(d): SP LDOS calculated at r . The insets show the SP-LDOS in the vicinity of the positive YSR state.Red and blue lines correspond to SP-LDOS calculated for opposite directions ˆ x and − ˆ x , whereas the dashed green line is theresulting spin polarization. Panels (a) and (b) correspond to the out-of-plane magnetic moment, whereas panels (c) and (d)correspond to the in-plane magnetic moment. subgap states [10–12, 37] E ± YSR = ± ∆ (cid:34) − (cid:18) πρS (cid:19) (cid:35)(cid:44) (cid:34) (cid:18) πρS (cid:19) (cid:35) , (8)which are unaffected by the Rashba SOC [38]. The ener-gies of the YSR states, however, depend on the ferromag-netic vector magnitude S . For a critical value S = 2 /πρ ,the energies of the YSR states reach E = 0, and thesystem undergoes a quantum phase transition as the twoYSR states cross [37, 39]. For simplicity, let us temporaryfix S = 2 / √ πρ , which corresponds to E ± YSR = ± ∆ / j ( r ) = lim r (cid:48) → r δ → +0 (cid:90) dω e iωδ πi (9) × Tr (cid:20) τ z (cid:18) i ∇ (cid:48) − ∇ m + λ ˆ z × σ (cid:19) G rr (cid:48) ( ω ) (cid:21) . In addition to the usual gradient term [40] in the paren-thesis, there is also a spin-dependent contribution due tothe Rashba SOC. We evaluate the current in Eq. (9) usingthe Green’s function in Eq. (5) and plot it in Figs. 2(a)and (c) for the cases of out-of-plane S = S ˆ z and in-plane S = S ˆ x moments, respectively. We note that onlythe pole in the T-matrix corresponding to the YSR stategives rise to nonzero currents. Both panels show con-centric patterns of current centered around the impurity.In the case of the out-of-plane moment, the current cir-culates around the impurity. In contrast, in the caseof the in-plane moment orientation, the current pointspredominantly in the y direction. The currents shownin Fig. 2(a) and (c) for point magnetic impurities agreequalitatively with the patterns obtained for the circu-lar island within the GL theory and shown in Figs. 1(b)and (c). However, in contrast with Fig. 1, the currentsin Fig. 2 display fine Friedel oscillations on the scale of r ∼ /p F . Note that the current in panel (c) is not con-tinuous. This can be understood by using the analogywith the Ginzburg-Landau current (4). For the in-planevector S , the current consists of the bare term α ( ˆ z × S ),as well as the condensate term n s m ∇ θ . These two distinctcontributions to the current are discontinuous, however,their sum is continuous. Since, the T-matrix calculationis not self-consistent, it does not take into account thereaction of the condensate that would fix the discontinu-ity. We discuss a fully self-consistent calculation, whichdemonstrates the continuity of the currents, in the nextsection, as well as in appendix.According to Eq. (1) the current and the spin polar-ization are coupled. Thus, we expect a nonzero in-planespin polarization even away from the impurity site thatsustains the nonlocal currents shown in Figs. 2(a) and(c). So, we evaluate the SP LDOS using the Green’sfunction as ρ j r ( ω ) = − π Im Tr (cid:20) τ z σ j G rr ( ω + iδ ) (cid:21) , (10)where j = x, y, z denotes the polarization axis. From theSP LDOS we also define the energy-dependent local spinpolarization as σ j r ( ω ) = ρ j r ( ω ) − ρ − j r ( ω ) . (11)In Figs. 2 (b) and (d), we plot both the SP LDOS and thespin polarization at the point r with solid and dashedlines, respectively. First, consider the out-of-plane mag-netic moment S = S ˆ z in panel (b). The SP LDOSpeaks at the superconducting coherence peak, i.e. at ω = ∆, as well as at the subgap YSR state energy,i.e. at ω = E YSR = ∆ /
2. The SP LDOS correspond-ing to the opposite directions j = ± x , shown with redand blue lines, are notably different at the YSR state.Therefore, the YSR state has a finite spin polarizationalong the x axis, shown with a dashed green line. Thisfeature is a consequence of the spin structure of the (a) (b) (c) (d)FIG. 3. (color online) Self-consistent numerical calculation of the currents induced by a point magnetic impurity with theout-of-plane moment S = S ˆ z . (a) Bogolyubov-de Gennes spectrum (top) and magnitude of current j (bottom) as a functionof S . (b)-(d) Spatial profile of the currents plotted at discrete points on the lattice for increasing magnitude of S as indicatedin panel (a). The current reaches the maximum value j c and switches direction for S = S c , i.e. where the YSR states crosszero energy. Note that in order to enhance visibility of the current in panels (b) and (d), the arrows representing the currentare magnified tenfold as indicated by the magnification ratio in the top-left corner of each panel. Green’s function Eq. (7) and vanishes in the absence ofthe Rashba SOC. Now, consider panel (d) correspond-ing to an in-plane moment S = S ˆ x . The YSR statein this case has a dominating spin polarization in the+ x direction with only a small admixture of the oppo-site spin. In Figs. 2(a) and (c), we plot the direction ofthe in-plane spin polarization for the positive YSR state σ r ( E +YSR ) = [ σ x r ( E +YSR ) σ y r ( E +YSR ) ] at the point r = 0as well as r = 7 ˆ x /p F and r = 7 ˆ y /p F . Note that thespin polarization of the negative YSR state is opposite,i.e. σ r ( E − YSR ) = − σ r ( E +YSR ). The current and spin po-larization are consistently orthogonal, which agrees withEq. (1). So, it is possible to map the current generatedby magnetic impurities and ferromagnetic islands usingspin-polarized scanning tunneling microscopy (SP STM). Self-consistent numerical modeling.-
The T-matrixapproximation discussed above predicts currents whichare qualitatively consistent with the GL results. How-ever, the T-matrix approach does not capture the influ-ence of the magnetic impurity on the superconductingorder parameter. It is known that the superconductingorder is strongly renormalized and may even change sign[37, 39, 41] in the vicinity of the magnetic impurity. In or-der to take this into account, we also perform a fully self-consistent numerical simulation[42] of the point magneticimpurity on a lattice [43–45] and show the results for anout-of-plane magnetic moment S = S ˆ z in Fig. 3. Panels(b)-(d) show the current for increasing values of the fer-romagnetic moment S . Note that the Friedel oscillationsare not fully visible here since the calculation is done fora coherence length such that ξ sc < /p F . In panel (a)we show the Bogolyubov-de Gennes spectrum (top) andthe magnitude of the current (bottom) as a function of S . For small S (b), the current circles around the impu-rity, which is consistent with both previous Figs. 1 and2. With further increase of S , the current grows and ul-timately undergoes a first-order discontinuous transitionat a critical value of magnetic vector S = S c . There, thecurrent abruptly reverses direction, as shown in Fig. 3(c)and (d), and reaches its maximal magnitude. This isaccompanied by the YSR states crossing at zero energy, and the superconducting order parameter reversing signat the impurity site. We note that the YSR states alsohave a first-order avoided crossing at zero energy [39] asshown in Fig. 3(a). With further increase of S , supercon-ductivity is suppressed and the currents diminish in thevicinity of the impurity. More details on the numericalsimulation can be found in the appendix. Concluding remarks.-
We have shown that super-conducting currents are generated by ferromagnetic is-lands and single magnetic impurities in 2D supercon-ductors with spin-orbit coupling. The currents originatefrom the magnetoelectric effect and are a direct conse-quence of combining SOC and magnetism. The discussedcurrents are unavoidable in ferromagnet-superconductorheterostructures, which have been proposed as a plat-form for topological superconductivity with the Majo-rana boundary states [1–3, 13–26]. We find a strong de-pendence of both the spatial pattern and magnitude ofthe currents on the direction of the ferromagnetic mo-ment. The currents are localized on the scale of the co-herence length in the case of the out-of-plane local mag-netic moment, whereas the currents have a dipolar powerlaw decay in the case of the in-plane magnetic moment.The presence of these non-local currents may induce long-range interactions between local magnetic moments on asuperconductor [46], which could qualitatively change thebehavior of the Majorana modes in such systems. Fur-thermore, by analyzing the currents in detail, we findthat they are carried by the subgap YSR states inducedby point magnetic impurities. The YSR states are spin-polarized, and the current is orthogonal to the local spinpolarization. Moreover, the current magnitude peakssharply at the phase transition, where the YSR statescross at zero energy. Thus, by using SP STM it shouldbe possible to map out the currents as well as detect thephase transition, which is paramount for finding TS andthe Majorana modes.We thank G. Volovik, M. Eschrig, Y. Kedem, andC. Triola for useful discussions. This work was supportedby the European Research Council (ERC) DM-321031and the US DOE BES E304 (S.S.P. and A.V.B.) and theSwedish Research Council (Vetenskapsr˚adet), the G¨oranGustafsson Foundation, and the Swedish Foundation for Strategic Research (SSF) (K.B. and A.B.-S.). [1] R. M. Lutchyn, J. D. Sau, and S. Das Sarma, “Ma-jorana Fermions and a Topological Phase Transition inSemiconductor-Superconductor Heterostructures,” Phys.Rev. Lett. , 077001 (2010).[2] Y. Oreg, G. Refael, and F. von Oppen, “Helical Liquidsand Majorana Bound States in Quantum Wires,” Phys.Rev. Lett. , 177002 (2010).[3] J. D. Sau, R. M. Lutchyn, S. Tewari, and S. Das Sarma,“Generic new platform for topological quantum compu-tation using semiconductor heterostructures,” Phys. Rev.Lett. , 040502 (2010).[4] A. Yu. Kitaev, “Unpaired Majorana fermions in quantumwires,” Phys. Usp. , 131 (2001).[5] J. Alicea, “New directions in the pursuit of Majoranafermions in solid state systems,” Rep. Prog. Phys. ,076501 (2012).[6] C. W. J. Beenakker, “Search for Majorana fermions insuperconductors,” Annu. Rev. Condens. Matter Phys. ,113 (2013).[7] N. Read and D. Green, “Paired states of fermions intwo dimensions with breaking of parity and time-reversalsymmetries and the fractional quantum Hall effect,”Phys. Rev. B , 10267 (2000).[8] D. A. Ivanov, “Non-Abelian Statistics of Half-QuantumVortices in p-Wave Superconductors,” Phys. Rev. Lett. , 268 (2001).[9] C. Nayak, S. H. Simon, A. Stern, M. Freedman, andS. Das Sarma, “Non-Abelian anyons and topologicalquantum computation,” Rev. Mod. Phys. , 1083(2008).[10] L. Yu, “Bound state in superconductors with paramag-netic impurities,” Acta Phys. Sin. , 75 (1965).[11] H. Shiba, “Classical Spins in Superconductors,” Prog.Theor. Phys. , 435 (1968).[12] A. I. Rusinov, “Superconductivity near a paramagneticimpurity,” JETP Lett. , 85 (1969).[13] T. P. Choy, J. M. Edge, A. R. Akhmerov, and C. W. J.Beenakker, “Majorana fermions emerging from magneticnanoparticles on a superconductor without spin-orbitcoupling,” Phys. Rev. B , 195442 (2011).[14] S. Nadj-Perge, I. K. Drozdov, B. A. Bernevig, andA. Yazdani, “Proposal for realizing Majorana fermionsin chains of magnetic atoms on a superconductor,” Phys.Rev. B , 020407 (2013).[15] J. Klinovaja, P. Stano, A. Yazdani, and D. Loss,“Topological Superconductivity and Majorana Fermionsin RKKY Systems,” Phys. Rev. Lett. , 186805 (2013).[16] M. M. Vazifeh and M. Franz, “Self-Organized TopologicalState with Majorana Fermions,” Phys. Rev. Lett. ,206802 (2013).[17] B. Braunecker and P. Simon, “Interplay between Classi-cal Magnetic Moments and Superconductivity in Quan-tum One-Dimensional Conductors: Toward a Self-Sustained Topological Majorana Phase,” Phys. Rev.Lett. , 147202 (2013).[18] F. Pientka, L. I. Glazman, and F. von Oppen, “Topo-logical superconducting phase in helical Shiba chains,” Phys. Rev. B , 155420 (2013).[19] S. Nakosai, Y. Tanaka, and N. Nagaosa, “Two-dimensional p-wave superconducting states with mag-netic moments on a conventional s-wave superconduc-tor,” Phys. Rev. B , 180503 (2013).[20] K. P¨oyh¨onen, A. Weststr¨om, J. R¨ontynen, and T. Oja-nen, “Majorana states in helical Shiba chains and lad-ders,” Phys. Rev. B , 115109 (2014).[21] Y. Kim, M. Cheng, B. Bauer, R. M. Lutchyn, andS. Das Sarma, “Helical order in one-dimensional mag-netic atom chains and possible emergence of Majoranabound states,” Phys. Rev. B , 060401 (2014).[22] I. Reis, D. J. J. Marchand, and M. Franz, “Self-organizedtopological state in a magnetic chain on the surface of asuperconductor,” Phys. Rev. B , 085124 (2014).[23] P. M. R. Brydon, S. Das Sarma, H.-Y. Hui, and J. D.Sau, “Topological Yu-Shiba-Rusinov chain from spin-orbit coupling,” Phys. Rev. B , 064505 (2015).[24] J. R¨ontynen and T. Ojanen, “Topological superconduc-tivity and high Chern numbers in 2D ferromagnetic Shibalattices,” arXiv:1412.5834.[25] J. Li, T. Neupert, Z. J. Wang, A. H. MacDonald,A. Yazdani, and B. A. Bernevig, “A novel platformfor two-dimensional chiral topological superconductiv-ity,” arXiv:1501.00999v1.[26] S. Nadj-Perge, I. K. Drozdov, J. Li, H. Chen, S. Jeon,J. Seo, A. H. MacDonald, B. A. Bernevig, and A. Yaz-dani, “Observation of Majorana fermions in ferromag-netic atomic chains on a superconductor,” Science ,602 (2014).[27] J. D. Sau and P. M. R. Brydon, “Bound states of a ferro-magnetic wire in a superconductor,” arXiv:1501.03149.[28] L. S. Levitov, Yu. V. Nazarov, and G. M. Eliashberg,“Magnetostatics of superconductors without an inversioncenter,” JETP Lett. , 445 (1985).[29] V. M. Edelstein, “Characteristics of the Cooper pairing intwo-dimensional noncentrosymmetric electron systems,”Sov. Phys. JETP , 1244 (1989).[30] V. M. Edelstein, “Magnetoelectric Effect in Polar Super-conductors,” Phys. Rev. Lett. , 2004 (1995).[31] S. K. Yip, “Two-dimensional superconductivity withstrong spin-orbit interaction,” Phys. Rev. B , 144508(2002).[32] E. Bauer and M. Sigrist, eds., Non-CentrosymmetricSuperconductors , Vol. 847 (Springer Berlin Heidelberg,2012).[33] T. Ojanen, “Magnetoelectric Effects in Superconduct-ing Nanowires with Rashba Spin-Orbit Coupling,” Phys.Rev. Lett. , 226804 (2012).[34] Note that we do the calculation at T = 0, where themagnetoelectric coefficients have a discontinuous jumpbetween the normal and superconducting phases. In con-trast at finite temperature T (cid:54) = 0, the magnetoelectriccoefficients interpolate smoothly as ∝ ∆ between thetwo phases.[35] K. V. Samokhin, “Magnetic properties of superconduc-tors with strong spin-orbit coupling,” Phys. Rev. B , , 373 (2006).[38] Y. Kim, J. Zhang, E. Rossi, and R. M. Lutchyn,“Impurity-induced bound states in superconductors withspin-orbit coupling,” arXiv:1410.4558.[39] M. I. Salkola, A. V. Balatsky, and J. R. Schrieffer, “Spec-tral properties of quasiparticle excitations induced bymagnetic moments in superconductors,” Phys. Rev. B , 12648 (1997).[40] A. A. Abrikosov, L. P. Gor’kov, and I. Y Dzyaloshinskii, Quantum Field Theoretical Methods in Statistical Physics (Pergamon, New York, 1965).[41] T. Meng, J. Klinovaja, S. Hoffman, P. Simon, andD. Loss, “Superconducting Gap Renormalization aroundtwo Magnetic Impurities: From Shiba to Andreev BoundStates,” arXiv:1501.07901.[42] The numerical simulation is done on a square latticewith nearest-neighbor hopping t = 1, spin-orbit coupling λ = 0 . t , chemical potential µ = − t . The supercon-ducting gap is determined self-consistently using a pairpotential v sc = 5 . t . Panels (b), (c) and (d) correspondto a magnetic impurity with S = 1 . t, . t , and 7 . t , re-spectively. For numerical reasons limiting the lattice size,the pair potential is chosen such that the superconductorcoherence length is of the order of the lattice constant.[43] A. M. Black-Schaffer and S. Doniach, “Self-consistent so-lution for proximity effect and Josephson current in bal-listic graphene SNS Josephson junctions,” Phys. Rev. B , 024504 (2008).[44] K. Bj¨ornson and A. M. Black-Schaffer, “Vortexstates and Majorana fermions in spin-orbit coupledsemiconductor-superconductor hybrid structures,” Phys.Rev. B , 024501 (2013).[45] K. Bj¨ornson and A. M. Black-Schaffer, “Probing vor-tex Majorana fermions and topology in semiconduc-tor/superconductor heterostructures,” Phys. Rev. B ,214514 (2015).[46] N. Y. Yao, L. I. Glazman, E. A. Demler, M. D.Lukin, and J. D. Sau, “Enhanced Antiferromagnetic Ex-change between Magnetic Impurities in a Superconduct-ing Host,” Phys. Rev. Lett. , 087202 (2014). Appendix A: Derivation of the extra terms in theGinzburg-Landau theory
The non-local coupling of the vector A = ( A x , A y ) = A + ∇ θ and ferromagnetic polarization S in 2D momen-tum q = ( q x , q y ) space is given by the extra term in thefree energy F extra = − (cid:90) q S a ( q ) K ab ( q ) A b ( − q ) . (A1) The tensor K ab is obtained by integrating out thefermions K ab ( q ) = 12 (cid:90) ω, p Tr (cid:104) σ a g p + q ( iω ) v b g p − q ( iω ) (cid:105) (A2)where the integration variables are given below the inte-gral sign for brevity, i.e. (cid:82) ω, p = (cid:82) dωd p (2 π ) . The integrationover continuous frequency ω corresponds to zero temper-ature T = 0. The factor 1 / g p ( iω ) = 1 iω − [ ξ ( p ) − λ ( σ × p ) z ] τ z − ∆ τ x , (A3) v = p m + λ ( ˆ z × σ ) . (A4)We expand the tensor K ab ( q ) up to the first order in qK ab ( q ) = K ab (0) + q c ∂ q c K ab (0) + O ( q ) . (A5)and find the tensors K ab (0) and ∂ c K ab (0) in the followingtwo subsections K ab (0) = (cid:15) zba α, α = λm π , (A6) ∂ c K ab (0) = iδ az (cid:15) zcb β, β = m λ πp F . (A7)Here, the coefficients α and β were evaluated in the limitof the small superconducting gap ∆ and large Fermi mo-mentum p F = √ mµ , i.e. λp F (cid:29) mλ (cid:29) ∆. It mayalso be useful to instead express the coefficients usingthe density of states ρ = m/π and Rashba momentum p R = mλ as α = ρ λ , β = ρ λ p R p F . (A8)We can now substitute Eqs. (A6) and (A7) in Eq. (A5)and rewrite Eq. (A1) in real space r = ( x, y ) as F extra = (cid:90) r α ( ˆ z × S ) · A + β ( ∇ S z × ˆ z ) · A . (A9)The term α was also recently derived in Ref. [36]. Weare not aware of a previous derivation of the term β .
1. Derivation of the tensor K ab (0) . Let us calculate the first term in this expansion, i.e. K ab ( ) = 12 (cid:90) ω, p Tr [ σ a g p ( p ) v b g p ( p )] . (A10)We define the operatorsΠ ± = 12 [1 ± ( σ × ˆ p ) z ] , ˆ p = p p , (A11)which project onto the spin eigenstates of the Rashbacoupling λ ( σ × p ) z corresponding to the two eigenstates ± λp . We expand the identity operator I via the projec-tion operators (A11) as I = Π + + Π − and simplify theGreen’s function g p ( iω ) = Π + iω + ξ + ( p ) τ + ∆ τ ( iω ) − E ( p ) +Π − iω + ξ − ( p ) τ + ∆ τ ( iω ) − E − ( p ) , (A12)where we defined the energies due to the Rashba splitting ξ ± ( p ) = ξ ( p ) ± λp, E ± ( p ) = (cid:113) ξ ± ( p ) + ∆ . (A13)We substitute Eq. (A12) in Eq. (A10) and obtain fourterms. Only the terms containing both Π + and Π − havedistinct poles, which produce a non-vanishing result uponthe frequency integration. Thus we obtain K ab (0) = 12 (cid:90) p Tr σ { σ a Π + v b Π − + σ a Π − v b Π + }× (cid:90) ω Tr τ { [ iω + ξ + ( p ) τ + ∆ τ ][ iω + ξ − ( p ) τ + ∆ τ ] } [( iω ) − E ( p )][( iω ) − E − ( p )] . (A14)where the first trace is over the spin space, whereas sec-ond trace is over the Nambu space. The trace in the sec-ond line of Eq. (A14) gives 2[ − ω + ξ + ( p ) ξ − ( p ) + ∆ ].Using the angular integration in momentum space, wesimplify the trace over spin matrices in the first line ofEq. (A14) to λ(cid:15) zab . Then, we integrate over frequencyand obtain K ab (0) = (cid:15) zba λ (cid:90) p E + ( p ) E − ( p ) − ξ + ( p ) ξ − ( p ) − ∆ E + ( p ) E − ( p )[ E + ( p ) + E − ( p )] . (A15)We can evaluate the integral in various limiting cases K ab (0) = (cid:15) zba mλ π (cid:26) λ p F / , µ (cid:29) ∆ (cid:29) λp F (cid:29) mλ , , λp F (cid:29) mλ (cid:29) ∆ , (A16)where p F = √ mµ is the Fermi momentum. The firstline of Eq. (A16) can also be obtained from Eq. (7) ofRef. [36].
2. Derivation of the tensor ∂ q c K ab (0) . We expand Eq. (A2) and obtain the first order coeffi-cient ∂ q c K ab = 14 (cid:90) ω, p Tr { σ a [ ∂ p c g p ( iω ) v b g p ( iω ) − g p ( iω ) v b ∂ p c g p ( iω )] } . We integrate by parts to shift the position of the derivate ∂ p c , use the identity ∂ p c g = g τ z v c g (which follows fromEqs. (A3) and (A4)) and obtain ∂ q c K ab = 12 (cid:90) ω, p Tr { σ a ∂ p c g p ( iω ) v b g p ( iω ) } = 12 (cid:90) ω, p Tr { σ a g p ( iω ) τ z v c g p ( iω ) v b g p ( iω ) } Below, we omit additional terms containing powers ofthe Green’s function g n under the trace because theyvanish upon the frequency integration (since poles lie onthe same side of the imaginary plane). We substitute theexpressions for the velocity Eq. (A4), omit the terms thatvanish under the angular integration in the momentumspace and simplify the equation to ∂ q c K ab (0) = λ (cid:15) z ˜ cc (cid:15) z ˜ bb (cid:90) ω, p Tr (cid:8) σ a g p ( iω ) τ z σ ˜ c g p ( iω ) σ ˜ b g p ( iω ) (cid:9) (A17)Here, we dropped the term containing p b p c Tr { σ a g p ( iω ) τ z g p ( iω ) } which vanishes upon thefrequency integration. We substitute the expansion ofthe Green’s functions Eq. (A12) in Eq. (A17) and tracesover the Nambu and spin matrices decouple such that (cid:88) s ,s ,s = ± Tr σ (cid:8) σ a Π s σ ˜ c Π s σ ˜ b Π s (cid:9) [( iω ) − E s ( p )][( iω ) − E s ( p )][( iω ) − E s ( p )] × Tr τ { [ iω + ξ s ( p ) τ z + ∆ τ x ] τ z [ iω + ξ s ( p ) τ z + ∆ τ x ] × [ iω + ξ s ( p ) τ z + ∆ τ x ] } . First we evaluate the trace over the spin matricesTr σ (cid:8) σ a Π s σ ˜ c Π s σ ˜ b Π s (cid:9) → (1 − s s ) i (cid:15) a ˜ c ˜ b , where the angular integration in the momentum spacewas invoked to simplify the expression. The remainingtrace over the Nambu matrices can also be evaluated asTr τ { [ iω + ξ s ( p ) τ z + ∆ τ x ] τ z [ iω + ξ s ( p ) τ z + ∆ τ x ] × [ iω + ξ s ( p ) τ z + ∆ τ x ] } = − ω [ ξ s ( p ) + ξ s ( p ) + ξ s ( p )]+ 2 ξ s ( p ) ξ s ( p ) ξ s ( p ) + 2∆ [ ξ s ( p ) + ξ s ( p ) − ξ s ( p )] . We integrate over frequency, substitute all terms inEq. (A17) and obtain ∂ q c K ab (0) = (A18) iδ az (cid:15) zcb (cid:90) p λ ∆ ξ ( p ) p (cid:20) ξ ( p ) λpE ( p ) + ξ ( p ) λpE − ( p ) + 1 E + ( p ) − E − ( p ) (cid:21) Similar to Eq. (A16), we evaluate the integral in the limitof small ∆ and obtain ∂ q c K ab (0) = iδ az (cid:15) zcb m λ πp F , λp F (cid:29) mλ (cid:29) ∆ . (A19) Appendix B: Ginzburg-Landau solution of currentsinduced by a ferromagnetic disc
In this section, we provide details of the calculation ofthe current j = δFδ A (cid:12)(cid:12)(cid:12)(cid:12) A =0 = n s m ∇ θ + α ( ˆ z × S ) + β ( ˆ z × ∇ S z ) . (B1)around a ferromagnetic region with disc geometry S ( r ) = S θ H ( R − r ) , (B2)where the index in θ H denotes the Heaviside theta func-tion to contrast it with the phase of the condensate θ .Let us first consider the case where the spin is out-of-plane, i.e. S = S ˆ z . Then, taking into account Eq. (B2),Eq. (B1) becomes j ( r ) = n s m ∇ θ + β [ ˆ z × ∇ S z ( r )] (cid:12)(cid:12)(cid:12) θ =0 (B3)= βS (ˆ r × ˆ z ) δ ( r − R ) . (B4)The current is thus localized around the boundary of theferromagnetic region as shown in Fig. 1(b). Note thatEq. (B4) corresponds to vanishing superconducting phase θ = 0, which is valid for small constant β . For larger val-ues of β , the vortex configuration of the superconductingphase, i.e. θ ( r ) = atan( y/x ), minimizes the free energyand the full expression Eq. (B3) for the current must beused.Now let us consider the case of an in-plane ferromag-netic vector, i.e. S = S ˆ x . Then, the current Eq. (B1)becomes j ( r ) = n s m ∇ θ + α [ ˆ z × S ( r )] . (B5)Notice here that the Euler-Lagrange equation for the su-perconducting phase θ gives the continuity equation forthe current0 = ∇ · j = n s m ∆ θ + α ∇ · [ ˆ z × S ( r )] . (B6)We rewrite Eq. (B6) as n s mα ∆ θ = ˆ z · [ ∇ × S ( r )] = − ˆ z · (ˆ r × S ) δ ( r − R ), use the Green’s function G L ( r ) = π ln( r ) for the 2D Laplace operator, whichsatisfies ∆ G L ( r ) = δ ( r ), and find the solution θ ( r ) = − mα(cid:15) zab S b n s (cid:82) d r (cid:48) G L ( r − r (cid:48) )ˆ r (cid:48) a δ ( r (cid:48) − R ), whichafter integration gives a simple result θ ( r ) = mα n s r r · ( S × ˆ z ) (cid:2) R + r − (cid:12)(cid:12) R − r (cid:12)(cid:12)(cid:3) (B7)= mαn s r · ( S × ˆ z ) (cid:26) , r < R, R r , r > R. (B8)Using Eqs. (B7) and (B2) we calculate the current (B5) j ( r ) = (cid:26) d R , r < R, r ( d · r ) r − d r , r > R, (B9) (a)(b)FIG. 4. (Color online.) Current around a point magneticimpurity in a square tight-binding model obtained using self-consistent calculations. The lattice constant is a and the cur-rent vector is indicated on each tight-binding site with theimpurity marked with a black dot. The direction of the im-purity magnetization is shown in the insets. where the “dipole” moment is defined as d = αR ( ˆ z × S ) / B = H + 4 π M induced by the ferromagnet of magnetization M . Inthe magnetostatics problem, the constant magnetization M is equivalent to α [ ˆ z × S ( r )], and the magnetic fieldstrength H is the term n s m ∇ θ in Eq. (B5). The diver-genceless magnetic induction B is equivalent to the cur-rent j . Appendix C: T-matrix calculation1. Green’s function in the real space
We calculate the Green’s function in the real space as g r ( ω ) = (cid:90) p e i p · r ω − [ ξ ( p ) + λ ( σ × p ) z ] τ z − ∆ τ x (C1)We substitute the expansion of the Green’s function(A12) in Eq. (C1) g r ( ω ) = 12 (cid:90) p (cid:20) i ( σ × ˆ p ) ω − ξ + ( p ) τ z − ∆ τ x + 1 − i ( σ × ˆ p ) ω − ξ − ( p ) τ z − ∆ τ x (cid:21) , where ξ ± ( p ) = ξ ( p ) ± λp . The angular integration inthe momentum space transforms the Rashba term in themomentum space ( σ × ˆ p ) into a corresponding term inreal space and produces the Bessel functions J and J g r ( ω ) = 12 (cid:90) p (cid:20) J ( pr ) + i ( σ × ˆ r ) J ( pr ) ω − ξ + ( p ) τ z − ∆ τ x + J ( pr ) − i ( σ × ˆ r ) J ( pr ) ω − ξ − ( p ) τ z − ∆ τ x (cid:21) , For small distance r (cid:28) ξ , we can substitute the momentain the Bessel function with their average values p = p ± F ,which satisfy equations ξ ± (cid:0) p ± F (cid:1) = 0. Then numeratorsdo not depend on the integration variable p , and, thus,integration gives g r ( ω ) = − π ω + ∆ τ x √ ∆ − ω [ f ( r ) + i ( σ × r ) f ( r )] , (C2)where f ( r ) = 12 (cid:2) ρ + J ( p + F r ) + ρ − J ( p − F r ) (cid:3) ,f ( r ) = 12 r (cid:2) ρ + J ( p + F r ) − ρ − J ( p − F r ) (cid:3) . In Eq. (C2), the density of states for each branch of theRashba spectrum has been defined as ρ ± = ρ ± mλ/p ± F , ρ = mπ . (C3)For example at r = 0, the Green’s function has a conven-tional form g ( ω ) = − πρ ω + ∆ τ x √ ∆ − ω , (C4)where ρ = ρ + + ρ − .
2. T-matrix and spin polarization
Using Eq. (C4), we can expand the T-matrix in Eq. (5) T ( ω ) = − S · σ − πρ S · σ ω +∆ τ x √ ∆ − ω , which after insertion of projectors (1 ± τ x ) / (cid:88) s = ± sτ x − S · σ − πρ S · σ ω + s ∆ √ ∆ − ω . We multiply the numerator and denominator of the frac-tion by 1 + πρ S · σ ω + s ∆ √ ∆ − ω to eliminate the σ matrices inthe denominator (cid:88) s = ± sτ x − S · σ − ( πρS ) ω + s ∆) ∆ − ω (cid:18) πρ S · σ ω + s ∆ √ ∆ − ω (cid:19) and after a few transformations rewrite the expression inthe final form T ( ω ) = (cid:88) s = ± sτ x s ∆ − ω (cid:104) ( πρS ) (cid:105) ( ω − E s YSR ) × (cid:18) S · σ + πρS ω + s ∆ √ ∆ − ω (cid:19) , (C5)where the YSR energies E ± YSR are given in Eq. (8). Wesubstitute Eq. (C5) in the expression for the current (9)and find that the last term in the parenthesis on the sec-ond line vanishes after taking the trace. The remainingterm proportional to S · σ has a pole at the YSR energyand gives a non-zero contribution to the current. Notethat, in general, we expect the contribution to the currentboth from the localized subgap states and the delocalizedsupragap states.Using the expansion of the T-matrix in Eq. (C5), wecalculate the spin-polarized LDOS (10) in the vicinity ofthe positive YSR state ρ ± xr ˆ x ( ω ) = [ f ( r ) ± rf ( r )] S ∆ δ ( ω − E +YSR )2 (cid:20) (cid:16) πρS (cid:17) (cid:21) (C6)for, e.g., perpendicular local moment S = S ˆ z . The firstterm in square brackets is responsible for distinct SP-LDOS in the opposite directions ± ˆ x , whereas the sec-ond fraction in Eq. (C6) determines the overall strengthof the YSR state. The terms f ( r ) and rf ( r ), definedin Eq. (C2), are of the same order sufficiently far fromthe impurity, and the YSR state acquires strong in-planespin polarization σ xr ˆ x = ρ + xr ˆ x − ρ − xr ˆ x ≈ ρ + xr ˆ x according toEq. (C6). Such a large spin polarization should be pos-sible to detect experimentally using SP-STM.0 Appendix D: Numerical simulation
We consider the following tight-binding model on asquare lattice (see e.g. Ref. [44, 45]): H = H kin + H so + H sc + H imp , (D1) H kin = − t (cid:88) (cid:104) i , j (cid:105) ,σ c † i σ c j σ − µ (cid:88) i ,σ c † i σ c i σ ,H so = iλ (cid:88) ib c † i + b σ ˆ z · ( σ × b ) σσ (cid:48) c i σ (cid:48) + H . c .,H sc = (cid:88) i ∆ i c † i ↑ c † i ↓ + H . c ..H imp = − c † σ ( S · σ ) σσ (cid:48) c σ (cid:48) , where i and j are site indices and b is a unit vector point-ing along one of the four types of bonds on the squarelattice. The parameters of the model are the strength ofthe nearest neighbor hopping t , the chemical potential µ ,and the ferromagnetic vector S . Also, a superconductingpair potential v sc is used to self-consistently determinethe superconducting order parameter through∆ i = − v sc (cid:104) c i ↓ c i ↑ (cid:105) = − v sc (cid:88) E ν < v ∗ ν i ↓ u i ↑ . (D2)We solve Eqs. (D1) and (D2) self-consistently and cal-culate currents around point magnetic impurities usingthe expressions derived below in Appendix D 1. Figure 3corresponds to the following parameters: t = 1, µ = − λ = 0 .
56, and v sc = 5 .
36. These values give a genericband structure of a lightly hole-doped Rashba SOC semi-conductor in proximity to a conventional s -wave super-conductor. Panels (b), (c) and (d) correspond to a mag-netic impurity with S = 1 . , .
72, and 7 .
36 respectively.Figure 4 is plotted for the parameters t = 1, µ = − λ = 0 . v sc = 5 .
36, and S = 2 .
56, which is just below S c .
1. Calculating current
Numerical expressions for the current are derived byconsidering the time rate of change of the density opera- tor ˆ ρ i = (cid:80) i σ c † i σ c i σ (see e.g. Ref. [43]): d ˆ ρ i dt = i (cid:126) [ H , ˆ ρ i ] . (D3)Let σ and ¯ σ be opposite spins and defineˆ S i σ = a i ¯ σσ c † i σ c i ¯ σ + a i ∗ ¯ σσ c † i ¯ σ c i σ , ˆ J ib σσ = − (cid:16) a ib σσ c † i + b σ c i σ + a ib ∗ σσ c † i σ c i + b σ (cid:17) , ˆ J ib σ ¯ σ = − (cid:16) a ib ¯ σσ c † i + b ¯ σ c i σ + a ib ∗ ¯ σσ c † i σ c i + b ¯ σ (cid:17) , (D4)where b runs over the vectors from site i to its nearestneighbors and a i ¯ σσ , a ib σσ , and a ib ¯ σσ are coefficients to be de-termined. It can be shown that Eq. (D3) can be writtenas d ˆ ρ i dt = (cid:88) σ ˆ S i σ − (cid:88) b σ (cid:16) ˆ J ib σσ + J ib σ ¯ σ (cid:17) . (D5)It is clear from the definition of Eq. (D4) that S i σ is anon-site source operator which converts ¯ σ spins into σ ,while J ib σσ is a current operator for σ -spins away fromsite i along bond b . Similarly J ib σ ¯ σ is a current operatorfor σ -spins away from site i along bond b , but for spinswhich are flipped in the transfer process. Summing thetwo types of currents for both spin species, we arrive atthe expression ˆ J ib = (cid:88) σσ (cid:48) ˆ J ib σσ (cid:48) , (D6)for the total operator for currents away from site i alongbond b . The current operator on site i can now be definedas ˆ J i = (cid:88) b b ˆ J ib . (D7)Evaluation of the commutator in Eq. (D3) reveals thatthe relevant coefficients in Eq. (D4) are a ib σσ = − it,a ib ¯ σσ = − ˆ z · ( σ × b ) ¯ σσ λ . (D8)Using Eq. (D4), (D7), and (D8) the current can finallybe calculated as (cid:68) ˆ J i (cid:69) = (cid:88) b σ,E ν < b (cid:8) it (cid:0) v ∗ i + b σ u i σ − v ∗ i σ u i + b σ (cid:1) + λ (cid:2) ˆ z · ( σ × b ) ¯ σσ v ∗ i + b ¯ σ u i σ + ˆ z · ( σ × b ) ∗ ¯ σσ v ∗ i σ u i + b ¯ σ (cid:3)(cid:9) ..