Spin-polarized edge currents and Majorana fermions in one- and two-dimensional topological superconductors
Kristofer Björnson, Sergey. S. Pershoguba, Alexander. V. Balatsky, Annica M. Black-Schaffer
SSpin-polarized edge currents and Majorana fermions in one- and two-dimensionaltopological superconductors
Kristofer Bj¨ornson , Sergey. S. Pershoguba , Alexander. V. Balatsky , , and Annica M. Black-Schaffer Department of Physics and Astronomy, Uppsala University, Box 516, S-751 20 Uppsala, Sweden Nordita, Center for Quantum Materials, KTH Royal Institute of Technology,and Stockholm University, Roslagstullsbacken 23, S-106 91 Stockholm, Sweden and Institute for Materials Science, Los Alamos National Laboratory, Los Alamos, NM 87545, USA
We investigate the persistent currents, spin-polarized local density of states, and spectral func-tions of topological superconductors constructed by placing ferromagnetic impurities on top of an s -wave superconductor with Rashba spin-orbit interaction. We solve self-consistently for the super-conducting order parameter and investigate both two-dimensional blocks and one-dimensional wiresof ferromagnetic impurities, with the magnetic moments pointing both perpendicular and parallelto the surface. We find that the topologically protected edge states of ferromagnetic blocks give riseto spin-polarized edge currents, but that the total persistent current flows in opposite direction towhat is expected from the dispersion relation of the edge states. We also show that the Majoranafermions at the end points of one-dimensional wires are spin-polarized, which can be directly relatedto the spin-polarization of the edge currents of two-dimensional blocks. Connections are also madeto the physics of the Yu-Shiba-Rusinov states for zero-dimensional impurities. PACS numbers: 74.90.+n, 03.65.Vf, 12.39.Dc
I. INTRODUCTION
During the past one and a half decade, topology hascome to play an important role in several different butrelated aspects of condensed matter physics. In 2001,Kitaev proposed that a one-dimensional (1D) spinless p -wave superconductor can provide a route towards topo-logical quantum computation, since it can host topologi-cally protected Majorana fermions. A few years later thesubject of topological insulators was born, with the es-sential feature being an insulating bulk and topologicallyprotected metallic interface states inside the bulk gap.
The close connection between Majorana fermions and theedge states in topological insulators was explicitly madeclear through the application of topological band theoryto superconductors.
An early proposal of such a topo-logical superconductor involved interfacing a topologicalinsulator with a conventional s -wave superconductor. This was followed by proposals to construct a relatedsetup using more conventional building blocks, such as s -wave superconductivity, magnetism, and Rashba spin-orbit interaction. This latter approach is particularlyinteresting because it can physically be realized in manydifferent ways and therefore leaves plenty of room for ex-perimental implementation.One of the most recent experimental investigations ofMajorana fermions involves putting ferromagnetic impu-rities on top of a superconductor to create a 1D topologi-cal superconductor wire, where surface effects give rise tothe needed Rashba spin-orbit interaction.
The samesetup can also be used to assemble structures such as2D ferromagnetic islands or blocks, which have attractedattention because of their topologically protected edgestates.
With this particular approach in mind, butalso knowing that the system can be implemented ina variety of ways, we here investigate both block and wire configurations of magnetic moments on top ofa Rashba spin-orbit coupled superconductor. To accu-rately model these systems, we solve self-consistently forthe superconducting order parameter.In particular, we investigate the persistent currents,local density of states (LDOS), and the spectral func-tions. We recently investigated the current patterns fora magnetic point impurity, providing insight into theproperties of the individual atoms in ferromagnetic do-mains and also connecting to the well-known physics ofYu-Shiba-Rusinov states of magnetic impurities in su-perconductors. Here, we instead start from the topo-logical band theory of a 2D ferromagnetic block and in-vestigate the resulting edge states, which are shown tobe spin-polarized. We then study 1D wires and connecttheir properties to the 2D blocks, by considering them asblocks with a width of a single site. Specifically, the end-point Majorana fermions are shown to be spin-polarizedwith a spin-polarization directly related to that of theedge states of the block. Through this approach we areable to fully elucidate the close connection between thetopological states in 2D and 1D, as well as the connectionto the Yu-Shiba-Rusinov states in 0D. In addition, wefind persistent currents around the ferromagnetic block,which surprisingly flows in opposite direction to what isexpected from the dispersion relation of the topologicallyprotected edge states. This apparent contradiction is re-solved as the gap-crossing edge states are indeed foundto give rise to a current compatible with their disper-sion relation, but lower energy quasiparticles are foundto produce even stronger counter-propagating currents, aphenomenon known in other superconductors. Thishas important experimental consequences, since physicalprobes such as superconducting quantum interference de-vices (SQUIDs), which measure the magnetic field gener-ated from the total persistent currents, will give opposite a r X i v : . [ c ond - m a t . s up r- c on ] S e p results to probes measuring transport properties involv-ing only states close to the Fermi level. II. MODEL
Motivated by the presence of Rashba spin-orbit inter-action, s -wave superconductivity, and magnetic impuri-ties in recent experimental setups, we consider thefollowing Hamiltonian H = H kin + H V z + H so + H sc , (1) H kin = − t (cid:88) (cid:104) i , j (cid:105) ,σ c † i σ c j σ − µ (cid:88) i ,σ c † i σ c i σ , H V z = − (cid:88) i ,σ,σ (cid:48) ( V z ( i )ˆ n · σ ) σσ (cid:48) c † i σ c i σ (cid:48) , H so = α (cid:88) ib (cid:16) e iθ b c † i + b ↓ c i ↑ + H . c . (cid:17) , H sc = (cid:88) i (cid:16) ∆ i c † i ↑ c † i ↓ + H . c . (cid:17) . Here i and j are site indices on a square lattice, b is avector pointing along the four types of nearest neighborbonds, θ b is the angular coordinate of b , and σ is thevector of Pauli matrices. The parameters of the modelare the strength of the hopping parameter t , chemicalpotential µ , Rashba spin-orbit interaction α , and su-perconducting pair potential V sc of the superconduct-ing substrate, as well as the effective Zeeman splittingcoming from the ferromagnetic impurities V z ( i ). The su-perconducting pair potential is used to self-consistentlydetermine the order parameter ∆ i = − V sc (cid:104) c i ↓ c i ↑ (cid:105) = − V sc (cid:80) E ν < v ∗ ν i ↓ u i ↑ , where u ν i σ and v ν i σ are the σ -spinelectron and hole components, respectively, of the ν theigenstate. We here set the hopping parameter t = 1,Rashba spin-orbit interaction α = 0 .
56, chemical po-tential µ = −
4, and superconducting pair potential V sc = 5 .
36. The specific values are not important, butour results are generally valid in a wide range of parame-ter space. The Zeeman term is set to zero everywhere,except for a block- or line-shaped region, where it is setto a finite V z , with the direction determined by the unitvector ˆ n . This allows us to study both 2D blocks and1D wires of ferromagnetic impurities embedded within aconventional superconductor. A. Currents
The most essential property of topologically nontrivialsystems is the appearance of gapless edge states. In thecase of quantum Hall systems, these give rise to persis-tent currents along the edges, while for topological in-sulators they instead result in persistent spin-currents.Likewise, topological superconductors have gapless edgestates, and it is therefore of interest to investigate the persistent currents in these systems. In fact, we findthat a study of these currents help us understand alsothe spin-polarization of Majorana fermions at the endpoints of 1D wires. We thus need not only expressionsfor calculating currents, but the spin-polarized currentsfor any given spin-polarization axis.To derive expressions for the currents, we consider thetime rate of change of the spin-density operator ˆ ρ i σ = c † i σ c i σ : d ˆ ρ i σ dt = i (cid:126) [ H , ˆ ρ i σ ] . (2)Because we are interested in arbitrary spin directions,the operators have to be understood to be written inthe basis of the particular spin direction of interest. If σ is in the direction ( θ, ϕ ) in spherical coordinates, theoperators are related to the operators in the basis of theHamiltonian according to c † i σ = cos (cid:18) θ (cid:19) c † i ↑ + sin (cid:18) θ (cid:19) e iϕ c † i ↓ . (3)Now let ¯ σ denote the opposite spin direction to σ anddefine ˆ 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) , (4)where a i ¯ σσ , a ib σσ , and a ib ¯ σσ are coefficients to be deter-mined later. It is then possible to show, as is done inAppendix B, that Eq. (2) can be written as d ˆ ρ i σ dt = ˆ S i σ − (cid:88) b (cid:16) ˆ J ib σσ + ˆ J ib σ ¯ σ (cid:17) , (5)where the sum runs over all bonds b for which thereis a term in the Hamiltonian (for Eq. (1) only nearestneighbor bonds). It is clear from Eq. (4) that ˆ S i σ is an on-site source operator which converts ¯ σ -spins into σ -spins.This operator is unrelated to currents and is thus of nofurther interest here. Next, ˆ J ib σσ is a current operator for σ -spins moving away from site i along bond b , and is fromEq. (5) seen to be responsible for decreasing the numberof σ -spins on site i . Finally, also ˆ J ib σ ¯ σ is responsible fortransferring σ -spins away from site i along bond b , butit flips the spin in the process.While ˆ J ib σσ allows us to investigate spin-polarized cur-rents, the total current is obtained by adding all contribu-tions to the current and we are therefore also interestedin the operator ˆ J ib = (cid:80) σσ (cid:48) ˆ J ib σσ (cid:48) . Here the summationover σ (cid:48) ensures that we add both spin-polarized currentsas well as ’spin-flipping’ currents, while the summationover σ ensures that we consider the outflow of both typesof spins. To be able to represent the current as a vectorfield, we also construct the on-site current operator byadding the currents along the bonds away from site i intothe directed-average vector operator ˆ J i = (cid:80) b b ˆ J ib ,and similarly for the spin-polarized and spin-flipping cur-rents.To fully construct the current operators we also needto evaluate the a -coefficients in Eq. (4). This is compli-cated by the fact that we are interested in spin-polarizedcurrents for arbitrary spin-polarization axes. In addi-tion, the final expressions for the currents are obtainedby evaluating expectation values (cid:104) ˆ J i (cid:105) and (cid:104) ˆ J i σσ (cid:48) (cid:105) , wherethe calculation also needs to be done in the particular ba-sis. A detailed derivation of the method used to calculatecurrents is provided in Appendix B. In the following wepresent the current in two complementary ways. First ofall, max i |(cid:104) ˆ J i (cid:105)| as a function of the magnetic moment al-lows us to identify the onset of persistent currents in thematerial. However, because the maximum is taken over the whole sample, it gives no information about wherethese currents appear. Therefore, for selected parame-ters we also show the current vector field. B. Spin-polarized LDOS
For 1D wires the topologically protected edge statesmanifest themselves in the form of Majorana fermions.For this reason we are interested in calculating the LDOSfor a cut along the wire. As we want to relate thespin-polarization of these Majorana fermions to the spin-polarization of the edge currents of 2D blocks, we arefurther interested in the spin-polarized LDOS for a givenspin-polarization axis. This is calculated using the ex-pression ρ ( i , θ, ϕ, E ) = (cid:88) ν (cid:2) u ∗ ν i ↑ u ∗ ν i ↓ (cid:3) (cid:20) cos ( θ ) cos( θ ) sin( θ ) e − iϕ cos( θ ) sin( θ ) e iϕ sin ( θ ) (cid:21) (cid:20) u ν i ↑ u ν i ↓ (cid:21) δ ( E − E ν ) . (6)The total LDOS is obtained by adding the spin-polarizedLDOS for two opposite spin-polarizations. C. Spectral function
The edge states of the 2D block also appear in thespectral function. Due to the limited size of the simu-lated system it is important to single out the contribu-tion to the spectral function coming from the edge states,otherwise bulk states pollute the edge spectrum. We dothis by introducing a state classification function C ( ν ),which is one for every state classified as an edge state,and zero otherwise. We here define a state to be an edgestate if more than 50% of its probability amplitude islocated inside a region, six sites wide and symmetricallypositioned across the boundary. With this definition, thecontribution to the spectral function coming from edgestates can be calculated as A ( k , E ) = (cid:88) ν,σ | u νk x σ | δ ( E − E ν ) C ( ν ) , (7)where u νk x σ is the Fourier transform of u νxσ , and u νxσ is the restriction of u ν i σ to the y -coordinate at which thebottom of the ferromagnetic block is located. III. DIMENSIONAL REDUCTION
For a homogeneous bulk it is possible to write the 2Dmodel of Eq. (1) as H ( k ) = (cid:20) H ( k ) ∆( k )∆ † ( k ) − H T ( − k ) (cid:21) , (8)where H ( k ) = (cid:15) ( k ) − V z σ z − L ( k ) · σ , (9)∆( k ) = i ∆ σ y , (10) (cid:15) ( k ) = − t [cos( k x ) + cos( k y )] − µ, (11) L ( k ) = α (sin( k y ) , − sin( k x ) , . (12)This Hamiltonian have been topologically classified andthe relevant condition for being in a topologically non-trivial phase is < | ∆ | < V z − (4 t + µ ) = V z , (13)where the last equality follows because µ = − t . More-over, making the restriction k y = 0 in Eq. (9) and defin-ing ˜ µ = µ + 2 t , we obtain H D ( k ) = − t cos( k x ) − ˜ µ − V z σ z + α sin( k x ) σ y , (14)which describes a 1D superconducting wire orientedalong the x -direction, recently intensively studied ine.g. Refs. [17,18,31]. It can easily be confirmed that thismodel, just like the 2D case, goes through a topologicalphase transition at | ∆ | = V z − (2 t + ˜ µ ) . This is oneexample of a dimensional reduction from two to one di-mension, and similar relations play an important role inthis work.Note, however, that if the 1D wire is embedded in a 2Dsurface, this type of dimensional reduction makes littlesense physically. First of all, setting k y = 0 correspondsto considering wave-functions that are completely delo-calized in the y -direction. Second, a 1D infinite wire isproduced by setting V z ( i ) = V z δ y , not by ignoring the y -direction. A proper dimensional reduction would ratherinvolve integrating out the y -dependence to arrive at anequation of similar form as Eq. (14). This would poten-tially renormalize the parameters, and the condition forbeing in the nontrivial phase should then be understoodin terms of these renormalized parameters. We explic-itly point this out here, because, while we have the strictcondition in Eq. (13) to rely on for identifying the topo-logical phase transition of a 2D geometry, we have nosuch condition for a 1D wire. Instead, we have to relyon a range of consequences, such as the appearance ofMajorana fermions and other indicators related to thecurrent, in order to identify the nontrivial phase of thewire. IV. RESULTSA. Block with magnetic moments perpendicular tosurface
We begin our investigation with a quadratic ferromag-netic block, with the magnetic moments perpendicular tothe surface (ˆ n = ˆ z ). In Fig. 1 the maximum value of thecurrent density and the size of the order parameter at thecenter of the ferromagnetic block is plotted as a functionof the Zeeman term, which quantify the strength of themagnetic impurity moments. The condition for being inthe topologically nontrivial phase, given in Eq. (13), isfulfilled in the shaded regions. It is clear that the maxi-mum current density is comparatively small in the trivialphase, raises up inside the topologically nontrivial phase,and then falls off as superconductivity is destroyed in theferromagnetic block. We also plot the spatial distributionof the persistent currents in Figs. 2-4 for the three pointsmarked 1-3 in Fig. 1. For small Zeeman splitting thecurrent flows anti-clockwise along the edge of the block,in agreement with a recent Ginzburg-Landau calculationfor a circular disk geometry. Next, in Fig. 3, the per-sistent current in the nontrivial phase is seen to not onlybe much larger, but it also flows in the opposite direc-tion to that in the trivial phase. We also note that thiscurrent is located mainly inside the ferromagnetic region.Below we show that these currents are related to the ap-pearance of edge states, but surprisingly propagates inthe direction opposite to what is expected from the edgestate spectrum. However, before we move on to a moredetailed study of the currents in the nontrivial phase,which is the main focus of this work, we also say a few x10 -2 z ' FIG. 1: (Color online.) Bottom: Order parameter in thecenter of the ferromagnetic block (blue line) as a functionof the Zeeman term V z in the direction ˆ n = ˆ z . Absolutevalue of the Zeeman term is indicated with dashed lines.Shaded regions indicate the topologically nontrivial phasewhere | V z | > | ∆ | >
0. Top: Maximum current density asa function of the Zeeman term. xy FIG. 2: (Color online.) Edge currents in the topologicallytrivial phase of a ferromagentic block magnetized along thedirection ˆ n = ˆ z , as indicated by grey arrows. The currentdensity is proportional to the arrow length, with the longestarrow corresponding to a current density of 2 . · − . Seealso (1) in Fig. 1. words about the currents in the third region.We see that the maximum current density remainsfairly large even after superconductivity is destroyed inthe ferromagnetic block region, and these currents arealso located around the edge of the block. Strong edgecurrents therefore remain even after the superconductinggap has been destroyed, even though it is not possible tomotivate their presence using topological arguments forthe block region. A few spikes also appear in the maxi-mum current density, and these are due to the formationof vortex-like currents inside the block region, as illus-trated in Fig. 4. The exact occurrence of these peaks asa function of the magnetization is dependent on the sizeof the block and these current patterns result from the xy FIG. 3: (Color online.) Edge currents in the topologicallynontrivial phase of a ferromagnetic block magnetized alongthe direction ˆ n = ˆ z . See also (2) in Fig. 1. xy FIG. 4: (Color online.) Currents at a current density peak inthe collapsed gap state of a ferromagnetic block magnetizedalong the ˆ n = ˆ z . See also (3) in Fig. 1. interplay of the Zeeman term, Rashba spin-orbit inter-action, and boundary conditions. A complete study ofthe currents in this regime is beyond the scope of thiswork, and we here only mention these effects to clarifythe properties of the non-zero current in the third region.
1. Reversed persistent current
Having covered the basic behavior of the currents, wenow investigate in detail the current in the nontrivialphase and relate it to the topologically protected edgestates. In Fig. 5 we plot the contribution to the spectralfunction from states classified as edge states for a cutalong the lower edge of the ferromagnetic block. The leftand the right figure corresponds to the nontrivial phasethat occurs for negative and positive Zeeman term, re- spectively. It is clear that the edge states have oppositedispersion, in agreement with the two phases having op-posite Chern numbers. We now note that the slope ofthe edge state in panel (2) of Fig. 5 implies that theedge state has a propagation direction to the right, since v ∼ dE/dk . However, this is in contradiction to the di-rection of the persistent current plotted in Fig. 3, whichflows to the left along the lower edge.The apparent contradiction is resolved once the con-tribution to the current from the edge states is singledout. In Fig. 6 we plot the contribution to the cur-rent from eigenstates with energy in the narrow interval − . < E < J ∼ (cid:82) E ( k ) 2. Spin-polarized current Next we study the spin-polarization of the currents. Inparticular, the x and y spin-polarization axes are of in-terest because the Rashba spin-orbit interaction couplesthese axes differently to the momentum. In Fig. 7 we plotthe total spin-polarized current as well as that carried byeigenstates in the energy interval − . < E < x -axis. The spin-polarized currentsfor y -up, x -down, and y -down spins are identical to theseunder a spatial rotation of π/ π , and 3 π/ 2, respectively.It is clear from Fig. 7(a) that strong spin-polarized cur-rents flow through the block, but the total current in thebulk cancels, since the spins polarized along the negative x -axis flow in opposite direction. However, in the edgeregions the two currents flow in the same direction and,together with the ’spin-flipping’ currents, these gives riseto the total edge currents in Fig. 3.We also note the partial depletion of x -up spin cur-rents towards the right end of the block in Fig. 7(a). Thepicture is radically different in Fig. 7(b), where only con-tributions to the currents from states close to the Fermilevel are included. Here, the x -up spins are located onthe right side, and with a similar argument to that inSec. IV A 1, we can understand the depletion of the total x -up spin current in Fig. 7(a) as a partial, but incom-plete, cancellation by the states close to the Fermi level.It is clear from both the total current in the bulk re-gion and the edge currents in Fig. 7(b) that x -up spinsfavor motion along the positive y -axis. This can be un-derstood as a consequence of the presence of Rashba spin-orbit interaction. The Rashba spin-orbit interaction cou-ples momentum and spin through k × σ , which leadsto spin-momentum locking for the currents. This spin-momentum locking and localization of different spins ondifferent edges are, as we will discuss next, also importantfor understanding 1D wire systems. B. Wire with magnetic moments perpendicular tosurface We now move on to an investigation of a 1D wire withthe magnetic moments still perpendicular to the surface,which corresponds to the 1D model in Eq. (14). As ex-plained in Sec. III, we no longer have a simple condi-tion, such as Eq. (13), to rely on to identify the topo-logical phase transition, because the system is embeddedin a 2D surface rather than being truly 1D. However, wefind that this present little problem since we have severalother means to identify the nontrivial phase. 1. Partial gap collapse and topological phase transition In Fig. 8 we plot the maximum current density and thevalue of the order parameter at the midpoint of a wire.First of all, we note that compared to Fig. 1 the horizon-tal axis has been rescaled and, compared to the block,a roughly twice as strong magnetization is required tosuppress superconductivity on the wires sites. A similarcalculation for a point impurity reveals that the magne-tization has to be about four times as large as as thatof the block geometry in order to destroy superconduc-tivity. This can be understood from the fact that eachedge site in a block is connected to a single nonmagneticsuperconducting site, while each site in the wire and forthe single impurity is connected to two and four suchsuperconducting sites, respectively.Next, we note that the partial collapse of the orderparameter in itself does not imply a topological phasetransition, even though we numerically always see these xy a xy FIG. 7: (Color online.) Spin-polarized current for spins along the positive x -axis for a ferromagnetic block magnetized along thedirection ˆ n = ˆ z . (a): Total spin-polarized current, (b): only contributions from states with energy in the interval − . < E < y -up, x -down, and y -down spins by a spatial rotation of π/ π ,and 3 π/ 2, respectively. x10 -3 70 V z FIG. 8: (Color online.) Similar to Fig. 1, but for a 1D wirealong the x -axis and magnetic moments along the directionˆ n = ˆ z . The topologically nontrivial phase is identified usingthe discontinuous onset of strong edge currents as well as theappearance of Majorana fermions. occurring at the same time. There are, however, goodreasons for why these can be expected to occur at thesame time. First of all, the topological phase transitionoccurs when ∆ and V z are related according to an expres-sion similar to Eq. (13), even though we can no longerrely on the last equality sign or the exact values of t and µ . The inequality is clearly valid when making the Zee-man term large enough, but also if the order parameteris small enough. A sudden drop in the order parameterwill therefore likely induce the topological phase transi-tion. Thus, assuming no a priori knowledge of the valueof ∆ at which the transition occurs, the likelihood of thetransition to occur is proportional to the size of the dropin the order parameter.An even stronger argument for why the two phenom-ena should occur in tandem can be made if we insteadconsider the drop in the order parameter to be induced by the topological phase transition. This relies on the factthat a topological phase transition is fundamentally re-lated to the interchange of energy levels across the Fermilevel, also known as a band inversion. Because the orderparameter is calculated using eigenstates below the Fermilevel, the value of the order parameter can be expectedto experience discontinuous changes each time a pairof eigenstates are interchanged across the Fermi level.Moreover, because the size of the order parameter inturn affects the eigenstates, this has enhancing feedbackeffects, which are only captured in fully self-consistentcalculations such as those performed here. Discontin-uous changes in the order parameter and the topologi-cal phase transition are therefore intimately tied to eachother. We note that this is closely related to the physicsof single magnetic impurities. The associated Yu-Shiba-Rusinov states, which in non-self-consistent calculationscross the Fermi level continuously, cross discontinuouslyin self-consistent models because of the feedback effectof the occupation numbers of these states on the orderparameter. 2. Signatures of nontrivial topology So far we have discussed the relation between the col-lapse of the gap and the topological phase transition,arguing that there are strong reasons to expect them tooccur at the same time. We have also mentioned thatthey always occur at the same time in our calculations,but so far not said what signatures we have used to iden-tify the topological phase transition. The first signatureis the strong peak in the maximum current density inFig. 8. Just like for the square block, the onset of thepeak is also associated with a reversal of the directionof the persistent current, as can be seen in Fig. 9. Most xy xy FIG. 9: (Color online.) Persistent currents around ferromag-netic wire with magnetic moments along the direction ˆ n = ˆ z .Left: Trivial phase. Right: Nontrivial phase. See also (4) and(5) in Fig. 8. Total x-up x-downE wire x wire x wire x FIG. 10: (Color online.) LDOS as a function of position alonga 1D wire when magnetic moments are perpendicular to thesurface (ˆ n = ˆ z ). Both total and spin-polarized LDOS areshown, using the x -axis as spin-polarization axis. The extentof the wire is indicated at the bottom. One Majorana fermionis clearly seen at each end of the wire. important though, this is also the point where Majoranafermions appear at the end points of the wire, as shownin Fig. 10. The Majorana fermion nature of these zero-energy states have been confirmed by checking that theyappear as locally non-degenerate Bogoliubov-de Gennesquasiparticles. 3. Majorana spin-polarization and block edge currents An interesting feature of the Majorana fermions seenin Fig. 10 are their spin-polarizations. The Majoranafermions at the right and left end of the wire are polar-ized along the x -up and x -down direction, respectively.That Majorana fermions have spin-polarization has al-ready been predicted, but we are here able to under-stand it in relation to the spin-polarized currents around kE y In fi nite Finite Single site FIG. 11: (Color online.) Schematic figures showing the edgestates along the rightmost edge of a ferromagnetic strip withvarying length. If the edge is made infinitely long the spec-trum is continuous, but for any finite length, the spectrumbecomes discrete. For a 1D wire placed along the x -axis, theMajorana fermion at its end point can be understood as thesingle state that remains when the length of the y -edge be-comes a single site. a 2D ferromagnetic block. In Fig. 7(b) it was shownthat the contribution to the x -up polarized current com-ing from states close to the Fermi level is located on theright edge. This is in complete agreement with the x -uppolarized Majorana fermion localized on the right end ofthe wire. This demonstrates the intimate connection be-tween the topologically protected edge states of a 2D and1D system.The connection between edge states and Majoranafermions can be made even clearer by considering theedge states plotted in Fig. 5. For simplicity of argu-ment, we first consider an infinite ferromagnetic striprather than a square block and let the infinite dimen-sion be along the y -direction. On each of the two edgesone spectral branch cuts the Fermi level, one of which isschematically depicted in the leftmost picture of Fig. 11.Next, consider the y -dimension wrapped around a cylin-der, such that it is still translationally invariant but withfinite length. This results in a discretization of the spec-tral function as shown in the middle picture of Fig. 11.By finally making the cylinder small enough such that the y -direction consists of a single site, a 1D wire with a sin-gle Majorana fermion at each edge results, as depicted tothe right in Fig. 11. This procedure closely resembles thedimensional reduction performed in Sec. III, and explic-itly shows the close relation between Majorana fermionsin 1D and edge states in 2D. However, this dimensionalreduction does not, just like that performed in Sec. III,correspond to the actual physical situation. This is im-portant to recognize because it might otherwise give thefalse impression that it is enough to discretize the edgespectrum to obtain a localized Majorana fermion. Afterall, the state at E = 0 in the middle panel of Fig. 11 isalso a Majorana fermion.The reason why the relation between the 2D block andthe 1D wire is slightly more complicated is that the two x -edges, which can be treated completely independentfor the infinite strip, is joined by two y -edges. The spec-tral branches of the different edges therefore join into onesingle branch with eigenstates delocalized over the wholeboundary. In particular, the Majorana fermions that oth-erwise would have existed at the right and left edge over-lap and hybridize. It is only once the y -dimension be-comes small enough, such that the two y -edges start hy-bridizing and thereby gaps out the spectrum along theseedges, that the two x -edges becomes independent of eachother. It is therefore clear why it is not enough to havea discrete spectrum, such as for a block with finite edgelength, to obtain localized Majorana fermions, and whya (quasi-)1D wire is required in spite of the close connec-tion between the two systems. C. Block with magnetic moments parallel tosurface Having investigated magnetic moments perpendicularto the surface, we now turn to a ferromagnetic block withthe magnetic moments pointing along the x -direction(ˆ n = ˆ x ). This system does not have a topologically non-trivial phase, but is of interest because of its relation tothe corresponding 1D wires.In Fig. 12 we plot the order parameter and maximumcurrent density. Interestingly, a strong peak is observedalso in this case, even though there is no topological phasetransition. However, the build up of the peak is not dis-continuous as in the previous cases. Instead, the maxi-mum current density is sizable already before supercon-ductivity starts to collapse in the ferromagnetic block.In fact, it is not only the maximum current density thatis large, but also the total persistent current. This isclear from Fig. 13, where it is seen that the current flowsthrough the whole sample, rather than being limited tothe boundary of the ferromagnetic block. A similar be-havior is predicted using a Ginzburg-Landau model inAppendix A.The current flowing through the bulk region of the fer-romagnetic block can be understood as a consequenceof x -up and x -down spins moving in opposite directionbecause of spin-orbit interaction, similar to what wasseen in Fig. 7. The addition of a Zeeman term in the x -direction creates an imbalance between the two spin-species and thereby also a total current through the sam-ple. The current continues to increase until the ferro-magnetic region enters the normal state. The peak inthe total current density is therefore not due to a suddenonset of new physics, but rather because the continuousbuild up of current comes to an end as superconductivityis destroyed. We also note that the region following thepartial collapse of the order parameter and up is com-plicated by the presence of vortex like currents insidethe ferromagnetic bulk region. See for example Fig. 14.For strong enough Zeeman term the currents becomes lo-calized on the left and right edges, similarly to what ispredicted using a Ginzburg-Landau model in AppendixA. x10 -3 70 V z FIG. 12: (Color online.) Same as in Fig. 1, but for magneticmoments along the direction ˆ n = ˆ x . xy FIG. 13: (Color online.) Persistent currents before the partialcollapse of the superconducting gap for a ferromagnetic blockmagnetized along the direction ˆ n = ˆ x . See also (6) in Fig. 12. xy FIG. 14: (Color online.) Persistent current pattern at maxi-mum current density peak at (7) in Fig. 12. Vortex-like solu-tions continue to appear in the ferromagnetic region also forsomewhat larger Zeeman terms, until the current eventuallybecomes localized on the left and right edges. Total z-up z-downE wire x wire x wire x FIG. 15: (Color online.) Same as in Fig. 10, but for ferrom-ganetic wire magnetized along the wire (ˆ n = ˆ x and wire alongˆ x ). Note that the Majorana fermions are spin-polarized alongthe z -axis rather than the x -axis. D. Wire with magnetic moments along wire We now consider a wire with the magnetic momentsalong the wire. This is the dimensional reduction of theblock in the previous section, where the width of theblock has been set to one along the y -direction. Thissystem can, just like the wire in Sec. IV B, be describedby a Hamiltonian of the form in Eq. (14), if the sub-stitution σ z → σ x is made. This substitution presentsno problem because the two Pauli matrices in Eq. (14)remain orthogonal, and the topological properties there-fore remain the same. It is in fact appropriate to thinkof the substitution as a rotation of spin indices aroundthe y -axis; σ x , σ z → − σ z , σ x . In Fig. 15 we plot theLDOS and spin-polarized LDOS using the z -axis as spin-polarization axis. It is clear that the Majorana fermionsspin-polarization that previously showed up in the x -component is now present in the z -component. More-over, direct comparison to Fig. 10 reveals that the x -upcomponent goes into the z -down components and viceversa, exactly as implied by the rotation of spin-indices.Finally, the maximum current density and the value ofthe order parameter at the wires midpoint are shown inFig. 16. The strong discontinuity in the maximum cur-rent density is correlated with the appearance of Majo-rana fermions. We also mention that the current patternfor the wire closely resembles that for the ferromagneticblock in Fig. 13. The current flows through the wireperpendicular to the direction of the wire, and the twopatterns are so similar that we refrain from showing it ina figure of its own. x10 -3 30 V z FIG. 16: (Color online.) Same as in Fig. 1, but for a ferro-magnetic wire magnetized along the wire (ˆ n = ˆ x and wirealong ˆ x ). The topologically nontrivial region has been identi-fied using the discontinuous onset of strong edge currents, aswell as the appearance of Majorana fermions. x10 -3 z FIG. 17: (Color online.) Same as in Fig. 1, but for wire within-plane magnetization perpendicular to wire (ˆ n = ˆ x and wirealong ˆ y ). E. Wire with in plane magnetic momentsperpendicular to wire The final setup we consider is a wire with in plane mag-netic moments perpendicular to the wire. Similarly asabove, this system can be obtained by making the widthof the ferromagnetic block in Sec. IV C equal to one inthe x -direction. However, because the dimension is re-duced in a different direction, the substitution of Paulimatrices needed to obtain an effective equation of theform in Eq. 14 is different. This time the substitution k x , σ y , σ z → k y , − σ x , σ x is required, and the two Paulimatrices for the spin-orbit and Zeeman term become thesame. The consequence is that the system has no topo-logically nontrivial phase, and no Majorana fermions ap-pear. The maximum current density and order parameterat the midpoint of the wire is presented in Fig. 17.1 V. SUMMARY We have studied ferromagnetic impurities on top of aconventional s -wave superconductor with Rashba spin-orbit interaction. We find a close relation between thetopological properties of 2D blocks and 1D wires of theferromagnetic impurities, as well as the physics of Yu-Shiba-Rusinov states of 0D point impurities. When the magnetic moments are directed perpendic-ular to the surface, both the 2D block and 1D wire hastopologically nontrivial phases. For a 2D block the non-trivial phase gives rise to spin-polarized edge currents.Interestingly, the total persistent current flows in op-posite direction to the current contribution from thetopologically protected gap-crossing edge states. Thismeans that different experimental probes will observedifferent current directions, depending on if they mea-sure the total persistent current or only that carried bylow-energy excitations. The Majorana fermions that ap-pear at the end points of 1D wires are spin-polarized,and their spin-polarization direction are directly relatedto the spin-polarized edge currents of the 2D block. Thespin-polarization of the Majorana fermions is along thedirection of the wire. When the magnetic moments are parallel to the sur-face the ferromagnetic block does not have a topologi-cally nontrivial phase. However, 1D wires can still havea topologically nontrivial phase, depending on if the wireis aligned with the magnetic moments or not. When thewire is parallel with the magnetic moments, a topolog-ically nontrivial phase exist and the corresponding Ma-jorana fermions are spin-polarized along the z -axis. Forboth the ferromagnetic block and the two wire directions,a substantial amount of current flows through the bulk,as long as superconductivity is present in the ferromag-netic region. VI. ACKNOWLEDGMENTS We thank G. Volovik, M. Eschrig, Y. Kedem, C. Tri-ola, and T. Ojanen for useful discussions. This workwas supported by the Swedish Research Council (Veten-skapsr˚adet), the G¨oran Gustafsson Foundation, and theSwedish Foundation for Strategic Research (SSF) (K.B.and A.B.-S.) and the European Research Council (ERC)DM-321031 and the US DOE BES E304 (S.S.P. andA.V.B.). A. Kitaev, Phys. Usp. , 131 (2001). C. L. Kane and E. J. Mele, Phys. Rev. Lett. , 226801(2005). B. A. Bernevig and S.-C. Zhang, Phys. Rev. Lett. ,106802 (2006). B. A. Bernevig, T. L. Hughes, and S.-C. Zhang, Science , 1757 (2006). M. Knig, S. Wiedmann, C. Br¨une, A. Roth, H. Buhmann,L. Molenkamp, X. L. Qi, and S.-C. Zhang, Science ,766 (2007). D. Hsieh, D. Qian, L. Wray, Y. Xia, Y. S. Hor, R. J. Cava,and M. Z. Hasan, Nature , 970 (2008). Y. Xia, D. Qian, D. Hsieh, L. Wray, A. Pal, H. Lin, A.Bansil, D. Grauer, Y. S. Hor, R. J. Cava, and M. Z. Hasan,Nat. Phys. , 398 (2009). H. Zhang, C.-X. Liu, X.-L. Qi, X. Dai, Z. Fang, and S.-C.Zhang, Nat. Phys. , 438 (2009). Y. L. Chen, J. G. Analytis, J.-H. Chu, Z. K. Liu, S.-K. Mo,X. L. Qi, H. J. Zhang, D. H. Lu, X. Dai, Z. Fang, S.-C.Zhang, I. R. Fisher, Z. Hussain, and Z.-X. Shen, Science , 178 (2009). M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. , 3045(2010). X. L. Qi and S.-C. Zhang, Rev. Mod. Phys. , 1057(2011). L. Fu and C. L. Kane, Phys. Rev. Lett. , 096407 (2008). M. Sato, Y. Takahashi, and S. Fujimoto, Phys. Rev. Lett. , 020401 (2009). M. Sato, Y. Takahashi, and S. Fujimoto, Phys. Rev. B ,134521 (2010). J. D. Sau, R. M. Lutchyn, S. Tewari, and S. Das Sarma,Phys. Rev. Lett. , 040502 (2010). J. Alicea, Phys. Rev. B , 125318 (2010). R. M. Lutchyn, J. D. Sau, and S. Das Sarma, Phys. Rev.Lett. , 077001 (2010). Y. Oreg, G. Refael, and F. von Oppen, Phys. Rev. Lett. , 177002 (2010). S. Nadj-Perge, I. K. Drozdov, J. Li, H. Chen, S. Jeon, J.Seo, A. H. MacDonald, B. A. Bernevig, and A. Yazdani,Science , 602 (2014). R. Pawlak, M. Kisiel, J. Klinovaja, T. Meier, S. Kawai, T.Glatzel, D. Loss, and E. Meyer, arXiv:1505.06078 (2015). M. Ruby, F. Pientka, Y. Peng, F. von Oppen, B. W. Hein-rich, and K. J. Franke, arXiv:1507.03104 (2015). J. R¨ontynen and T. Ojanen, Phys. Rev. Lett. , 236803(2015). Jian Li, Titus Neupert, Zhi Jun Wang, A. H. MacDon-ald, A. Yazdani, and B. Andrei Bernevig, arXiv:1501.00999(2015). V. Mourik, K. Zuo, S. M. Frolov, S. R. Plissard, E. P. A.M. Bakkers, and L. P. Kouwenhoven, Science , 1003(2012). L. P. Rokhinson, X. Liu, and J. K. Furdyna, Nat. Phys. ,795 (2012). A. Das, Y. Ronen, Y. Most, Y. Oreg, M. Heiblum, and H.Shtrikman, Nat. Phys. , 887 (2012). S. S. Pershoguba, K. Bj¨ornson, A. M. Black-Schaffer, A.V. Balatsky, arXiv:1505.01672 (2015). J. A. Sauls, Phys. Rev. B , 214509 (2011). H. Wu and J. Sauls, Phys. Rev. B , 184506 (2013). K. Bj¨ornson, A. M. Black-Schaffer, Phys. Rev. B ,024501 (2013). K. Bj¨ornson, A. M. Black-Schaffer, Phys. Rev. B ,134518 (2014). A. V. Balatsky, I. Vekhter, and J.-X. Zhu, Rev. Mod. Phys. , 373 (2006). D. Sticlet, C. Bena, and P. Simon, Phys. Rev. Lett. ,096802 (2012). Appendix A: Ginzburg Landau approach To complement the self-consistent calculations on alattice, we here also provide the results of a Ginzburg-Landau (GL) treatment. The GL method gives an ap-proximate description of the system, but is unable to cap-ture microscopic details such as spin-polarization, edgestates, density of states etc.. Nevertheless, it may beinteresting to compare the GL results with the exact nu-merical diagonalization discussed throughout the paper.We consider the following dimensionless GL free energyfunctional F = (cid:90) d r (cid:20) (cid:0) | ψ | − (cid:1) − g ψ S ( r )2 + a ( r ) | ψ | (A1)+ b | ∇ ψ | + ic ψ ∗ ∇ ψ − ∇ ψ ∗ ψ ) · ( ˆ z × S ) (cid:21) , Here, ψ ( r ) describes the superconducting order parame-ter in a 2D box, i.e. r = ( x, y ) and − L/ < x, y < L/ L = 1. The magnetic block of size d = 0 . L isplaced in the center with the magnetization along the x axis, i.e. S ( r ) = S ˆ x for − d/ < x, y < d/ 2. The fer-romagnetic vector S is analogous to the Zeeman term V z used throughout the main text. The first term inEq. (A1) is conventional and favors the superconduct-ing order at | ψ | = 1, whereas the second term describesthe diamagnetic Pauli energy disfavoring superconduct-ing order. The superconducting order is destroyed whenthe Pauli energy exceeds the the superconducting con-densation energy (the so-called Chandrasekhar-Clogstonlimit), which occurs for S > S p = 1 in Eq. (A1). Inorder to smoothly interpolate the Pauli term betweenthe superconducting and metallic phases we introduceda sigmoid function g ψ = 2 e v | ψ | + 1 , v = 5 . (A2)Note that g ψ ≈ ψ ≈ ψ (cid:54) = 0). To avoid pos-sible spurious numerical effects, we also add a boundaryterm with a large coefficient a ( r ) ∼ b = 0 . ξ ∼ √ b = 0 . 01. Thelast term, proportional to c = 0 . between the current and magne-tization S .We search for the minimum of the free energy (A1)numerically. We employ the gradient descent method inwhich the order parameter evolves in imaginary time τ as dψdτ = − δFδψ . (A3)Using Eq. (A3), the order parameter is updated itera-tively as ψ new = ψ old − dτ δF [ ψ old ] δψ . The time evolutionstops for the solution ψ ( r ) satisfying δF [ ψ ] δψ = 0, which,thus, minimizes the free energy.In Fig. 18 we plot the order parameter magnitude | ψ | ,phase Arg | ψ | , and the current j = ib ( ψ ∗ ∇ ψ − ∇ ψ ∗ ψ ) + c | ψ | ( ˆ z × S ) . (A4)The panels in the top (bottom) row correspond to thefield S = 0 . S = 2), which is below (above) the Paulilimit S p = 1. For the smaller field (top row), the super-conducting order | ψ | is uniform, whereas the supercon-ducting phase is non-uniform and is distributed in such away that the total current is continuous. In contrast,for the large field (bottom row), the superconductingorder is destroyed in the region covered with the ferro-magnetic block. The superconducting phase is stronglyuniform outside the ferromagnetic block. The currentis suppressed everywhere, except in small vicinity of theferromagnetic block.In general, the GL modeling agrees qualitatively withthose discussed in Section IV C: the current grows lin-early with the magnetic field until superconductivity isdestroyed. At the same time, we are not able to capturethe non-linear peak of the current shown in Fig. 12 aswell as the appearance of the vortex structure in Fig. 14. Appendix B: Derivation current expressions1. Introduction Here we derive the general expressions for calculatingthe currents in an arbitrary spin-basis. The imple-mentation of the numerical calculations are based onthe algorithm described in Sec. B 6. The method itselfis fairly straightforward. However, the derivation iscomplicated by the fact that the spin-polarized currentsare not necessarily calculated in the same basis as theHamiltonian. A short summary of the appendix follows: B 2: Notation Introduction of useful notation and derivations ofimportant identities. B 3: Generalized current operators The main problem is formulated. General sink-sourceand current operators are defined. In summary, thereare three important terms for any given spin: The sink-source term ˆ S s , the current ˆ J bσσ on bond b , and the spin-3 (a) (b) (c)(d) (e) (f)FIG. 18: (Color online.)(a, d) Order parameter magnitude | ψ | , (b, e) phase Arg ψ , and currents calculated using Ginzburg-Landau theory. The position of the ferromagnetic block is indicated with a red dashed rectangle. The ferromagnetic vector isin-plane S = S ˆ x . Top row corresponds to S = 0 . S = 2. The magnitude of the currents in panel (f) issix-fold magnified compared to panel (c) as indicated by the magnification ratio shown in these panels. flipping current ˆ J bσ ¯ σ on b , which together satisfies d ˆ ρ σ dt = ˆ S s − (cid:88) b (cid:16) ˆ J σσ + ˆ J σ ¯ σ (cid:17) . (B1)To derive these operators, we need to find a certain setof coefficients a λλ (cid:48) κκ (cid:48) s + b,η,s,η (cid:48) . B 4: Arbitrary basis General expressions are derived for how to calculatethe a λλ (cid:48) κκ (cid:48) s + b,η,s,η (cid:48) coefficients in an arbitrary basis from theknowledge of those same coefficients in the basis of theHamiltonian. B 5: Calculate sink-sources and currents Expressions for the expectation values of the sink-source and current operators are derived. B 6: Summary of method An algorithmic summary of the method is provided. B 7: Tabulating relevant a ’s The exact number of a λλ (cid:48) κκ (cid:48) s + b,η,s,η (cid:48) that needs to beevaluated in the basis of the Hamiltonian are deter-mined, and a useful way of tabulating them is introduced. B 8: Evaluating tables Tables are evaluated for the following types of termsin the Hamiltonian: chemical potential, kinetic, Zeeman, and Rashba spin-orbit interaction. The total table fora Hamiltonian containing all these terms is also provided. B 9: Note on superconductivity It is shown that the introduction of an s -wave super-conducting term in the Hamiltonian does not modifythe method described in the rest of this appendix. 2. Notation Here we introduce useful notation and show relationswhich will be important later on. This section is purelyalgebraic, and no physical motivation will be provided.The usefulness of the expressions will only become clearonce general currents are treated in Section B 3. a. Bond decomposition of commutators Consider a general quadratic Hamiltonian H = (cid:88) αλλ (cid:48) h αλλ (cid:48) , (B2)where h αλλ (cid:48) = (cid:88) ij f αiλjλ (cid:48) c † iλ c jλ (cid:48) . (B3)4Here i, j are site indices, λ, λ (cid:48) are spin indices, and α is an index distinguishing between different terms in theHamitlonian such as kinetic, chemical potential, Zeeman, etc..We now consider expressions of the form i (cid:126) (cid:2) h αλλ (cid:48) , c † sκ c sκ (cid:48) (cid:3) = i (cid:126) (cid:88) ij f αiλjλ (cid:48) (cid:16) δ js δ λ (cid:48) κ c † iλ c sκ (cid:48) − δ is δ λκ (cid:48) c † sκ c jλ (cid:48) (cid:17) = i (cid:126) (cid:88) b (cid:16) f αs + b,λ,sλ (cid:48) δ λ (cid:48) κ c † s + b,λ c sκ (cid:48) − f αsλ,s + b,λ (cid:48) δ λκ (cid:48) c † sκ c s + b,λ (cid:48) (cid:17) . (B4)The summation in the last expression is said to run overall bonds b connected to site s , where a bond between site s and s + b is said to exist if the Hamiltonian contains aproduct of creation and annihilation operators with thesetwo indices. In particular, the set of bonds has to be understood to include the bond between a site and itself ifany on-site terms are included in the Hamiltonian. UsingEq. (B4), it is possible to decompose the commutator forthe complete Hamiltonian as i (cid:126) (cid:2) H, c † sκ c sκ (cid:48) (cid:3) = i (cid:126) (cid:88) αλλ (cid:48) (cid:2) h αλλ (cid:48) , c † sκ c sκ (cid:48) (cid:3) = i (cid:126) (cid:88) bλ (cid:32)(cid:32)(cid:88) α f αs + b,λ,sκ (cid:33) c † s + b,λ c sκ (cid:48) − (cid:32)(cid:88) α f αsκ (cid:48) ,s + b,λ (cid:33) c † sκ c s + b,λ (cid:33) . (B5)Defining a λλ (cid:48) κκ (cid:48) i,η,j,η (cid:48) = δ λη δ κ (cid:48) η (cid:48) δ λ (cid:48) κ i (cid:126) (cid:88) α f αi,λ,j,λ (cid:48) ,b λλ (cid:48) κκ (cid:48) i,η,j,η (cid:48) = δ κη δ λ (cid:48) η (cid:48) δ λκ (cid:48) i (cid:126) (cid:88) α f αi,λ,j,λ (cid:48) , (B6) the expression can further be written as i (cid:126) (cid:2) H, c † sκ c sκ (cid:48) (cid:3) = (cid:88) bλλ (cid:48) ηη (cid:48) (cid:16) a λλ (cid:48) κκ (cid:48) s + b,η,s,η (cid:48) c † s + b,η c sη (cid:48) − b λλ (cid:48) κκ (cid:48) s,η,s + b,η (cid:48) c † sη c s + b,η (cid:48) (cid:17) = (cid:88) bλ (cid:16) a λκκκ (cid:48) s + b,λ,s,κ (cid:48) c † s + b,λ c sκ (cid:48) − b κ (cid:48) λκκ (cid:48) s,κ,s + b,λ c † sκ c s + b,λ (cid:17) . (B7) b. Useful expressions The Hamiltonian is Hermitian if H = H † , which usingEq. (B2) implies (cid:88) αλλ (cid:48) ij f αiλjλ (cid:48) c † iλ c jλ (cid:48) = (cid:88) αλλ (cid:48) ij f α ∗ iλjλ (cid:48) c † jλ (cid:48) c iλ . (B8) Identification of coefficients leads to the conclusion that (cid:80) α f α ∗ iλjλ (cid:48) = (cid:80) α f αjλ (cid:48) iλ , and using Eq. (B6) we have b λλ (cid:48) κκ (cid:48) iηjη (cid:48) = − a λ (cid:48) λκ (cid:48) κ ∗ jη (cid:48) iη . (B9)It is therefore possible to rewrite Eq. (B7) as5 i (cid:126) (cid:2) H, c † sκ c sκ (cid:48) (cid:3) = (cid:88) bλλ (cid:48) ηη (cid:48) (cid:16) a λλ (cid:48) κκ (cid:48) s + b,η,s,η (cid:48) c † s + b,η c sη (cid:48) + a λ (cid:48) λκ (cid:48) κ ∗ s + b,η (cid:48) ,s,η c † sη c s + b,η (cid:48) (cid:17) = (cid:88) bλ (cid:16) a λκκκ (cid:48) s + b,λ,s,κ (cid:48) c † s + b,λ c sκ (cid:48) + a λκ (cid:48) κ (cid:48) κ ∗ s + b,λ,s,κ c † sκ c s + b,λ (cid:17) , (B10)which eliminates the need to calculate b coefficients in-dependently. We also note from Eq. (B6) that a λλκκiηiη = b λλκκiηiη , which together with Eq. (B9) gives a λλκκiηiη = 0 . (B11)Using ¯ κ (cid:48) and ¯ η (cid:48) to denote the spins which are opposite to κ (cid:48) and η (cid:48) , respectively, we also find from Eq. (B6) that a λλ (cid:48) κκ (cid:48) i,η,j,η (cid:48) = a λλ (cid:48) κ ¯ κ (cid:48) i,η,j, ¯ η (cid:48) (B12) c. Obtaing a ’s through projection In Eq. (B6) explicit expressions for how to calculate a λλ (cid:48) κκ (cid:48) iηjη (cid:48) were provided. It will also be useful to have away of extracting coefficients from a general quadraticexpression, also when the explicit expression for thesecoefficients are not known. To this end we consider ageneral quadratic expression e = (cid:88) µν a µν c † µ c ν , (B13) where µ and ν are some arbitrary set of indices. It isclear that (cid:104) | c ρ ec † ρ (cid:48) | (cid:105) = a ρρ (cid:48) . (B14)Applying this to Eq. (B10) we can for b (cid:54) = 0 write (cid:104) | c s + b,η i (cid:126) (cid:2) H, c † sκ c sκ (cid:48) (cid:3) c † sη (cid:48) | (cid:105) = (cid:88) λλ (cid:48) a λλ (cid:48) κκ (cid:48) s + b,η,s,η (cid:48) = a ηκκκ (cid:48) s + b,η,s,η (cid:48) , (B15)while for b = 0 we have (cid:104) | c sη i (cid:126) (cid:2) H, c † sκ c sκ (cid:48) (cid:3) c † sη (cid:48) | (cid:105) = a ηκκκ (cid:48) s,η,s,η (cid:48) + a η (cid:48) κ (cid:48) κ (cid:48) κ ∗ s,η (cid:48) ,s,η , (B16)which can also be written as (cid:104) | c sη i (cid:126) (cid:2) H, c † sκ c sκ (cid:48) (cid:3) c † sη (cid:48) | (cid:105) = η (cid:54) = κ and η (cid:48) (cid:54) = κ (cid:48) ,a ηκκκ (cid:48) sηsη (cid:48) if η (cid:54) = κ and η (cid:48) = κ (cid:48) ,a η (cid:48) κ (cid:48) κ (cid:48) κ ∗ sη (cid:48) sη if η = κ and η (cid:48) (cid:54) = κ (cid:48) ,a ηκκκ (cid:48) s,η,s,η (cid:48) + a η (cid:48) κ (cid:48) κ (cid:48) κ ∗ s,η (cid:48) ,s,η if η = κ and η (cid:48) = κ (cid:48) . (B17) 3. Generalized current operators With the notation introduced in the previous section,we are now ready to formulate the problem in a waythat eventually will allow us to calculate spin-polarizedcurrents in any given basis. To this end we consider thelocal density operator for particles of spin σ on site s ,which can be written as ˆ ρ sσ = c † sσ c sσ . The time rate ofchange of this operator is d ˆ ρ sσ dt = i (cid:126) [ H, ˆ ρ sσ ] . (B18) The right hand side of this expression is of the form con-sidered in the previous section, and we therefore know6from Eq. (B10) that it can be written as d ˆ ρ sσ dt = (cid:88) bλ (cid:16) a λσσσs + b,λ,s,σ c † s + b,λ c sσ + a λσσσ ∗ s + b,λ,s,σ c † sσ c s + b,λ (cid:17) = (cid:0) a σσσσs,σ,s,σ + a σσσσ ∗ s,σ,s,σ (cid:1) c † sσ c sσ + a ¯ σσσσs, ¯ σ,s,σ c † s ¯ σ c sσ + a ¯ σσσσ ∗ s, ¯ σ,s,σ c † sσ c s ¯ σ + (cid:88) b (cid:54) =0 (cid:16) a σσσσs + b,σ,s,σ c † s + b,σ c sσ + a σσσσ ∗ s + b,σ,s,σ c † sσ c s + b,σ (cid:17) + (cid:88) b (cid:54) =0 (cid:16) a ¯ σσσσs + b, ¯ σ,s,σ c † s + b, ¯ σ c sσ + a ¯ σσσσ ∗ s + b, ¯ σ,s,σ c † sσ c s + b, ¯ σ (cid:17) , (B19)where ¯ σ is the opposite spin of σ . In the last step theon-site terms have been singled out of the expression,and the remaining sum over the bonds have been dividedinto one part free of spin-flips, and one involving spin-flips. Using Eq. (B11) we can further conclude that thefirst term is zero. The next two terms converts σ - and¯ σ -spins into each other on-site. When considering thedensity of a single spin species such as σ , these termstherefore acts as sink and source terms, respectively. Wetherefore define the sink-source operator on site s asˆ S s = a ¯ σσσσs, ¯ σ,s,σ c † s ¯ σ c sσ + a ¯ σσσσ ∗ s, ¯ σ,s,σ c † sσ c s ¯ σ . (B20)The third line describes processes where a σ spin jumpsbetween site s and site s + b . It is therefore the cur-rent into site s along bond b . We accordingly define thecurrent operator along bond b asˆ J bσσ = − (cid:16) a σσσσs + b,σ,s,σ c † s + b,σ c sσ + a σσσσ ∗ s + b,σ,s,σ c † sσ c s + b,σ (cid:17) , (B21)where it is implicitly understood that the operator is onlydefined for actual bonds between different sites. We havehere included an additional minus sign in order to definethe direction of the current to be along the bond awayfrom site s . The third term is similarly related to the inand out flow of σ -spins on site s . However, in contrast tothe ordinary currents the spins flip during the transition.We therefore define the spin-flipping current along bond b asˆ J bσ ¯ σ = − (cid:16) a ¯ σσσσs + b, ¯ σ,s,σ c † s + b, ¯ σ c sσ + a ¯ σσσσ ∗ s + b, ¯ σ,s,σ c † sσ c s + b, ¯ σ (cid:17) . (B22)With these definitions we can now write d ˆ ρ sσ dt = ˆ S s − (cid:88) b (cid:54) =0 (cid:16) ˆ J bσσ + ˆ J bσ ¯ σ (cid:17) . (B23)It is clear that in order to have a complete understand-ing of the currents in a given basis σ , it is sufficient toevaluate the three coefficients (cid:8) a ¯ σσσσs, ¯ σ,s,σ , a σσσσs + b,σ,s,σ , a ¯ σσσσs + b, ¯ σ,s,σ (cid:9) . (B24) These can all be written as a σ (cid:48) σσσs + b,σ (cid:48) ,s,σ . (B25)With knowledge of the eigenstates of the Hamiltonian,the expressions just derived allow us to calculate on-sitesink-sources, as well as currents on bonds, using S s = (cid:104) ˆ S s (cid:105) and J bσσ (cid:48) = (cid:104) ˆ J bσσ (cid:48) (cid:105) , respectively. However, it is alsouseful to know the current on a specific site, rather thanon the bonds. If b is the vector along bond b , from s to s + b , we also define the two types of currents at site s as J sσσ (cid:48) = 12 (cid:88) b (cid:54) =0 b J bσσ (cid:48) . (B26)This is a directed average of all currents flowing out fromsite s and forms a vector field, in contrast to the scalarcurrents defined on the bonds. 4. Arbitrary basis Although we in principle can proceed to calculate therelevant coefficients in Eq. (B24) using Eq. (B6), thisrequires the Hamiltonian to be expressed in the samebasis as σ . It would be a nuisance if each time we wantto calculate the current in a new basis, the Hamiltonianwould need to be rewritten in that basis, commutators bereevaluated, and finally also the expectation values (cid:104) ˆ S s (cid:105) and (cid:104) ˆ J bσσ (cid:48) (cid:105) . It is much better if the current in one basiscan be evaluated from the commutators and expectationvalues calculated in the basis of the original Hamiltonian.Consider therefore operators with arbitrary spin direc-tion c sσ = b ↑ c s ↑ + b ↓ c s ↓ ,c s ¯ σ = b ∗↓ c s ↑ − b ∗↑ c s ↓ ≡ ¯ b ↑ c s ↑ + ¯ b ↓ c s ↓ , (B27)where σ is some arbitrary spin, and ¯ σ is the oppositespin. Further, ↑ and ↓ here refers to up and down spinsin the basis of the Hamiltonian. If σ is along the direction( θ, ϕ ), in ordinary polar coordinates relative to the basisof the Hamiltonian, the coefficients can be chosen to be b ↑ = − ¯ b ∗↓ = cos( θ ,b ↓ =¯ b ∗↑ = cos( θ e − iϕ . (B28)and the relevant commutator can be written as i (cid:126) (cid:2) H, c † sσ c sσ (cid:3) = (cid:88) κκ (cid:48) b ∗ κ b κ (cid:48) i (cid:126) (cid:2) H, c † sκ c sκ (cid:48) (cid:3) . (B29)It is clear that the left hand side can be calculated for anyspin direction σ if the right hand side has been evaluatedin the basis of the Hamiltonian. However, we remem-ber that it is the expansion of the left hand side into a coefficients that allows us to calculate currents. Assume7therefore that the a coefficients in the expansion of thecommutators on the right hand side have been evaluated.We denote these coefficients by a capital A , to distinguishthem from the a coefficients in the expansion in the basisof the expression on the left hand side.We will now show how a coefficients in any basis can be calculated from a knowledge of the A coefficients. Tothis end, we remind ourselves that we in fact are onlyinterested in calculating the three a ’s in Eq. (B24). Fur-ther, using Eq. (B15) and (B17), these can be writtenas a ¯ σσσσs, ¯ σ,s,σ = (cid:104) | c s, ¯ σ d ˆ ρ σ dt c † sσ | (cid:105) = (cid:88) κκ (cid:48) κ (cid:48)(cid:48) κ (cid:48)(cid:48)(cid:48) b ∗ κ b κ (cid:48) ¯ b κ (cid:48)(cid:48) b ∗ κ (cid:48)(cid:48)(cid:48) (cid:104) | c s,κ (cid:48)(cid:48) i (cid:126) (cid:2) H, c † sκ c sκ (cid:48) (cid:3) c † sκ (cid:48)(cid:48)(cid:48) | (cid:105) ,a σσσσs + b,σ,s,σ = (cid:104) | c s + b,σ d ˆ ρ σ dt c † sσ | (cid:105) = (cid:88) κκ (cid:48) κ (cid:48)(cid:48) κ (cid:48)(cid:48)(cid:48) b ∗ κ b κ (cid:48) b κ (cid:48)(cid:48) b ∗ κ (cid:48)(cid:48)(cid:48) (cid:104) | c s + b,κ (cid:48)(cid:48) i (cid:126) (cid:2) H, c † sκ c sκ (cid:48) (cid:3) c † sκ (cid:48)(cid:48)(cid:48) | (cid:105) a ¯ σσσσs + b, ¯ σ,s,σ = (cid:104) | c s + b, ¯ σ d ˆ ρ σ dt c † sσ | (cid:105) = (cid:88) κκ (cid:48) κ (cid:48)(cid:48) κ (cid:48)(cid:48)(cid:48) b ∗ κ b κ (cid:48) ¯ b κ (cid:48)(cid:48) b ∗ κ (cid:48)(cid:48)(cid:48) (cid:104) | c s + b,κ (cid:48)(cid:48) i (cid:126) (cid:2) H, c † sκ c sκ (cid:48) (cid:3) c † sκ (cid:48)(cid:48)(cid:48) | (cid:105) . (B30)The two last expressions can be simplified straightfor-wardly by using Eq. (B15) and then summing over indicesappearing in δ -functions in Eq. (B6), which gives a σσσσs + b,σ,s,σ = (cid:88) κκ (cid:48) κ (cid:48)(cid:48) b ∗ κ | b κ (cid:48) | b κ (cid:48)(cid:48) A κ (cid:48)(cid:48) κκκ (cid:48) s + b,κ (cid:48)(cid:48) ,s,κ (cid:48) ,a ¯ σσσσs + b, ¯ σ,s,σ = (cid:88) κκ (cid:48) κ (cid:48)(cid:48) b ∗ κ | b κ (cid:48) | ¯ b κ (cid:48)(cid:48) A κ (cid:48)(cid:48) κκκ (cid:48) s + b,κ (cid:48)(cid:48) ,s,κ (cid:48) . (B31)The on-site term is slightly more complicated due to thespecial cases in Eq. (B17): a ¯ σσσσs, ¯ σ,s,σ = (cid:88) κκ (cid:48) κ (cid:48)(cid:48) κ (cid:48)(cid:48)(cid:48) b ∗ κ b κ (cid:48) ¯ b κ (cid:48)(cid:48) b ∗ κ (cid:48)(cid:48)(cid:48) (cid:16) A κ (cid:48)(cid:48) κκκ (cid:48) s,κ (cid:48)(cid:48) ,s,κ (cid:48)(cid:48)(cid:48) + A κ (cid:48)(cid:48)(cid:48) κ (cid:48) κ (cid:48) κ ∗ s,κ (cid:48)(cid:48)(cid:48) ,s,κ (cid:48)(cid:48) (cid:17) . (B32) 5. Sink-sources and currents Using Eqs. (B20)-(B22) together with Eq. (B31) andEq. (B32), we can now write the expectation values forthe sink-source and current operators as (cid:104) ˆ S s (cid:105) = (cid:88) κκ (cid:48) (cid:0) a ¯ σσσσs, ¯ σ,s,σ ¯ b ∗ κ b κ (cid:48) (cid:104) c † sκ c sκ (cid:48) (cid:105) H + a ¯ σσσσ ∗ s, ¯ σ,s,σ b ∗ κ ¯ b κ (cid:48) (cid:104) c † sκ c sκ (cid:48) (cid:105) H (cid:1) , (cid:104) ˆ J bσσ (cid:105) = − (cid:88) κκ (cid:48) b ∗ κ b κ (cid:48) (cid:16) a σσσσs + b,σ,s,σ (cid:104) c † s + b,κ c sκ (cid:48) (cid:105) H + a σσσσ ∗ s + b,σ,s,σ (cid:104) c † sκ c s + b,κ (cid:48) (cid:105) H (cid:17) , (cid:104) ˆ J bσ ¯ σ (cid:105) = − (cid:88) κκ (cid:48) (cid:16) a ¯ σσσσs + b, ¯ σ,s,σ ¯ b ∗ κ b κ (cid:48) (cid:104) c † s + b,κ c sκ (cid:48) (cid:105) H + a ¯ σσσσ ∗ s + b, ¯ σ,s,σ b ∗ κ ¯ b κ (cid:48) (cid:104) c † sκ c s + b,κ (cid:48) (cid:105) H (cid:17) . (B33)Here also the expectation values on the right hand sidehave been expanded in the basis of the Hamiltonian, which is indicated by the subscript H . Using that (cid:104) c † sκ c s + b,κ (cid:48) (cid:105) = (cid:104) c † s + b,κ (cid:48) c sκ (cid:105) ∗ , this can finally be written8as (cid:104) ˆ S s (cid:105) =2Re (cid:32) a ¯ σσσσs, ¯ σ,s,σ (cid:88) κκ (cid:48) ¯ b ∗ κ b κ (cid:48) (cid:104) c † sκ c sκ (cid:48) (cid:105) H (cid:33) , (cid:104) ˆ J bσσ (cid:105) = − (cid:32) a σσσσs + b,σ,s,σ (cid:88) κκ (cid:48) b ∗ κ b κ (cid:48) (cid:104) c † s + b,κ c sκ (cid:48) (cid:105) H (cid:33) , (cid:104) ˆ J bσ ¯ σ (cid:105) = − (cid:32) a ¯ σσσσs + b, ¯ σ,s,σ (cid:88) κκ (cid:48) ¯ b ∗ κ b κ (cid:48) (cid:104) c † s + b,κ c sκ (cid:48) (cid:105) H (cid:33) . (B34) 6. Summary of method Here we give a short algorithmic summary of themethod. We note that the two first steps only needs tobe performed once in the basis of the Hamiltonian. Re-peated application of the subsequent steps then allowsfor the sink-sources and currents in arbitrary bases to becalculated. a. Calculate expansion coefficients of commutators in thebasis of the Hamiltonian A λλ (cid:48) κκ (cid:48) s + b,η,s,η (cid:48) = δ λη δ λ (cid:48) κ δ κ (cid:48) η (cid:48) i (cid:126) (cid:88) α f αs + b,λ,s,λ (cid:48) . (B35)See Sec. B 8 for already evaluated terms for chemical po-tential, kinetic, Zeeman, and Rashba spin-orbit interac-tion. b. Calculate expectation values in the basis of theHamiltonian (cid:104) c † s + b,κ c sκ (cid:48) (cid:105) H . (B36) c. Choose spin direction Choose appropriate coefficients b ↑ and b ↓ for the spindirection σ of interest c sσ = b ↑ c s ↑ + b ↓ c s ↓ . (B37)For a spin pointing in the direction ( θ, ϕ ), use b ↑ = cos( θ ,b ↓ = sin( θ e − iϕ . (B38)The coefficients for the opposite spin follows as¯ b ↑ = b ∗↓ , ¯ b ↓ = − b ∗↑ . (B39) d. Calculate relevant a ’s in the chosen basis a ¯ σσσσs, ¯ σ,s,σ = (cid:88) κκ (cid:48) κ (cid:48)(cid:48) κ (cid:48)(cid:48)(cid:48) b ∗ κ b κ (cid:48) ¯ b κ (cid:48)(cid:48) b ∗ κ (cid:48)(cid:48)(cid:48) (cid:16) A κ (cid:48)(cid:48) κκκ (cid:48) s,κ (cid:48)(cid:48) ,s,κ (cid:48)(cid:48)(cid:48) + A κ (cid:48)(cid:48)(cid:48) κ (cid:48) κ (cid:48) κ ∗ s,κ (cid:48)(cid:48)(cid:48) ,s,κ (cid:48)(cid:48) (cid:17) ,a σσσσs + b,σ,s,σ = (cid:88) κκ (cid:48) κ (cid:48)(cid:48) b ∗ κ | b κ (cid:48) | b κ (cid:48)(cid:48) A κ (cid:48)(cid:48) κκκ (cid:48) s + b,κ (cid:48)(cid:48) ,s,κ (cid:48) ,a ¯ σσσσs + b, ¯ σ,s,σ = (cid:88) κκ (cid:48) κ (cid:48)(cid:48) b ∗ κ | b κ (cid:48) | ¯ b κ (cid:48)(cid:48) A κ (cid:48)(cid:48) κκκ (cid:48) s + b,κ (cid:48)(cid:48) ,s,κ (cid:48) . (B40) e. Calculate sink-sources and currents (cid:104) ˆ S s (cid:105) =2Re (cid:32) a ¯ σσσσs, ¯ σ,s,σ (cid:88) κκ (cid:48) ¯ b ∗ κ b κ (cid:48) (cid:104) c † sκ c sκ (cid:48) (cid:105) H (cid:33) , (cid:104) ˆ J bσσ (cid:105) = − (cid:32) a σσσσs + b,σ,s,σ (cid:88) κκ (cid:48) b ∗ κ b κ (cid:48) (cid:104) c † s + b,κ c sκ (cid:48) (cid:105) H (cid:33) , (cid:104) ˆ J bσ ¯ σ (cid:105) = − (cid:32) a ¯ σσσσs + b, ¯ σ,s,σ (cid:88) κκ (cid:48) ¯ b ∗ κ b κ (cid:48) (cid:104) c † s + b,κ c sκ (cid:48) (cid:105) H (cid:33) . (B41) f. Current vector field From the currents on the bonds, a current vector field,defined on site s , is finally calculated as J sσσ (cid:48) = 12 (cid:88) b (cid:54) =0 b (cid:104) ˆ J bσσ (cid:48) (cid:105) (B42) 7. Tabulating relevant a ’s Having summarized the method, it only remains toevaluate the a coefficients in the basis of the Hamilto-nian, which are denoted by a capital A . By counting thenumber of spin indices in Eq. (B6) to be six, it may seemlike we in general need to evaluate 2 = 64 a ’s per bondto be fully able to expand the commutators on the formin Eq. (B10). However, the three δ -functions ensure thatonly 2 = 8 of these are non-zero. For any given bond,we can therefore restrict ourselves to evaluate a λκκκ (cid:48) s + b,λ,s,κ (cid:48) .Further, Eq. (B12) allows us to reduce this number tofour, by fixing κ (cid:48) = ↑ . The right hand side of Eq. (B10)can then be completely determined for any site s from afull knowledge of the table( i, j ) A ↑↑↑↑ i ↑ j ↑ A ↑↓↓↑ i ↑ j ↑ A ↓↑↑↑ i ↓ j ↑ A ↓↓↓↑ i ↓ j ↑ (1,1)(1,2)... ... ... ... ...( n , n )9Here i, j are site indices, with n the number of sites. Theelements necessary to reconstruct Eq. (B10) for any givensite s can be found in the table by setting i = s + b and j = s . Note that also most terms in this table will bezero. From Eq. (B6) it is clear that an a will only be non-zero if the Hamiltonian contains some f αiλjλ (cid:48) connectingthe two sites. It is therefore sufficient to further limit theevaluation of the table to those pairs of sites that sharea bond.We also note that if we assume that each h α = (cid:80) λλ (cid:48) h αλλ (cid:48) is Hermitian, all results derived in this sectionapplies to the individual h α ’s. In particular, the entriesin the table for the full Hamiltonain H can be obtainedby adding the corresponding entries in the partial tablesof the individual h α ’s. 8. Evaluating tables a. Chemical potential The chemical potential can be written as h µ = (cid:88) iλ µ i c † iλ c iλ = (cid:88) ijλλ (cid:48) f µiλjλ c † iλ c iλ (cid:48) (B43)where f µiλjλ (cid:48) = δ ij δ λλ (cid:48) µ i . (B44)The table is therefore (all entries are to be multiplied by i/ (cid:126) ) ( i, j ) A ↑↑↑↑ i ↑ j ↑ A ↑↓↓↑ i ↑ j ↑ A ↓↑↑↑ i ↓ j ↑ A ↓↓↓↑ i ↓ j ↑ ( i, i ) µ i µ i Note: Although the chemical potential gives rise tosome non-zero terms, it does not give rise to a contribu-tion to any of the sink-source or current operators. Thatit does not contribute to the currents is clear from thefact that it only connects terms on the same site. To seethat it also does not contribute to the sink-source termwe note that the chemical potential term is rotationallyinvariant and can be considered to already be written inthe σ -basis for any given σ . That is, the A ’s in the tablecan be considered to be written directly in the σ -basis.Eq. (B20) reveals that only a λλ (cid:48) κκ (cid:48) s,η,s,η (cid:48) with λ (cid:54) = λ (cid:48) entersinto the sink-source term, all of which are zero in thetable above. b. Kinetic term The spin-degenerate kinetic energy can be written as h t = (cid:88) ijλ t ij c † iλ c jλ = (cid:88) ijλλ (cid:48) f tiλjλ (cid:48) c † iλ c jλ (cid:48) , (B45) where f tiλjλ (cid:48) = δ λλ (cid:48) t ij . (B46)The table is therefore (all entries are to be multiplied by i/ (cid:126) ) ( i, j ) A ↑↑↑↑ i ↑ j ↑ A ↑↓↓↑ i ↑ j ↑ A ↓↑↑↑ i ↓ j ↑ A ↓↓↓↑ i ↓ j ↑ ( i, j ) t ij t ij c. Zeeman term The Zeeman term can be written as h V z = (cid:88) iλ V z ( ˆn · σ ) λλ (cid:48) c † iλ c iλ = (cid:88) ijλλ (cid:48) f V z iλjλ (cid:48) c † iλ c jλ (cid:48) , (B47)where f V z iλjλ (cid:48) = δ ij δ λλ (cid:48) V z ( ˆn · σ ) λλ (cid:48) . (B48)The table is therefore (all entries are to be multiplied by i/ (cid:126) ) ( i, j ) A ↑↑↑↑ i ↑ j ↑ A ↑↓↓↑ i ↑ j ↑ A ↓↑↑↑ i ↓ j ↑ A ↓↓↓↑ i ↓ j ↑ ( i, i ) V z n z V z ( n x − in y ) V z ( n x + in x ) − V z n z d. Rashba spin-orbit interaction The Rashba spin-orbit interaction can be written as h SO = (cid:88) ijλλ (cid:48) α ij (cid:0) ( δ i,j − ˆ x − δ i,j +ˆ x ) ( − iσ y ) λλ (cid:48) + i ( δ i,j − ˆ y − δ i,j +ˆ y ) ( σ x ) λλ (cid:48) ) c † iλ c jλ (cid:48) = (cid:88) ijλλ (cid:48) f SOiλjλ (cid:48) c † iλ c jλ (cid:48) , (B49)where f SOiλjλ (cid:48) = α ij (cid:0) ( δ i,j − ˆ x − δ i,j +ˆ x ) ( − iσ y ) λλ (cid:48) + i ( δ i,j − ˆ y − δ i,j +ˆ y ) ( σ x ) λλ (cid:48) ) . (B50)The table is therefore (all entries are to be multiplied by i/ (cid:126) ) ( i, j ) A ↑↑↑↑ i ↑ j ↑ A ↑↓↓↑ i ↑ j ↑ A ↓↑↑↑ i ↓ j ↑ A ↓↓↓↑ i ↓ j ↑ ( i, i + ˆ x ) 0 − α i,i +ˆ x α i,i +ˆ x i, i − ˆ x ) 0 α i,i − ˆ x − α i,i − ˆ x i, i + ˆ y ) 0 iα i,i +ˆ y iα i,i +ˆ y i, i − ˆ y ) 0 − iα i,i − ˆ y − iα i,i − ˆ y e. Total For a total Hamiltonain of the form H = H µ + H t + H V z + H SO , (B51)the total table is obtanied as the sum of the individualtables (all entries are to be multiplied by i/ (cid:126) )( i, j ) A ↑↑↑↑ i ↑ j ↑ A ↑↓↓↑ i ↑ j ↑ A ↓↑↑↑ i ↓ j ↑ A ↓↓↓↑ i ↓ j ↑ ( s, s ) V z n z V z ( n x − in y ) V z ( n x + in y ) − V z n z ( s − ˆ x, s ) t s − ˆ x,s − α s − ˆ x,s α s − ˆ x,s t s − ˆ x,s ( s + ˆ x, s ) t s +ˆ x,s α s +ˆ x,s − α s +ˆ x,s t s +ˆ x,s ( s − ˆ y, s ) t s − ˆ y,s iα s − ˆ y,s iα s − ˆ y,s t s − ˆ y,s ( s + ˆ y, s ) t s +ˆ y,s − iα s +ˆ y,s − iα s +ˆ y,s t s +ˆ y,s We have here ignored the contributions from H µ , becauseas noted above all relevant terms are zero. We have alsoassumed that t ij only is non-zero for nearest neighbourterms. 9. Note on superconductivity All terms in the Hamiltonian treated so far arequadratic of the form c † iλ c jλ . However, even within amean-field treatment of ordinary s -wave superconduc-tivity, we will also encounter terms in the Hamiltonianwhich are of the form h SC = (cid:88) i ∆ i c † i ↑ c † i ↓ + ∆ ∗ c i ↓ c i ↑ , (B52)where ∆ i = − V sc (cid:104) c i ↓ c i ↑ (cid:105) , for some V sc . This term is aspin-singlet and therefore appears similarly in any givenspin-basis, such that we can assume that σ = ↑ . If theterm is included in the Hamiltonian, then it is clear from i (cid:126) (cid:104) h SC , c † s ↑ c s ↑ (cid:105) = i (cid:126) (cid:16) ∆ ∗ c s ↓ c s ↑ − ∆ s c † s ↑ c † s ↓ (cid:17) (B53) that Eq. (B23) has to be modified to read d ˆ ρ s ↑ dt = ˆ S s − (cid:88) b (cid:16) ˆ J b ↑↑ + ˆ J b ↑↓ (cid:17) + i (cid:126) (cid:16) ∆ ∗ s c s ↓ c s ↑ − ∆ s c † s ↑ c † s ↓ (cid:17) . (B54)Taking the expectation value of this expression and usingthat ∆ s = − V sc (cid:104) c s ↓ c s ↑ (cid:105) = (cid:16) − V sc (cid:104) c † s ↑ c † s ↓ (cid:105) (cid:17) ∗ = (∆ ∗ s ) ∗ wefind d (cid:104) ˆ ρ s ↑ (cid:105) dt = (cid:104) ˆ S s (cid:105) − (cid:88) b (cid:16) (cid:104) ˆ J b ↑↑ (cid:105) + (cid:104) ˆ J b ↑↓ (cid:105) (cid:17) . (B55)Rewriting this as d (cid:104) ˆ ρ s ↑ (cid:105) = (cid:32) (cid:104) ˆ S s (cid:105) − (cid:88) b (cid:16) (cid:104) ˆ J b ↑↑ (cid:105) + (cid:104) ˆ J b ↑↓ (cid:105) (cid:17)(cid:33) dt, (B56)it is possible to view any non-zero term as a possiblegenerator of density fluctuations. Although the super-conducting condensate can carry current, it is clear fromthe vanishing expectation value that it does not in itselfgenerate density fluctuations. In fact, it is clear fromEq. (B54) that the superconducting term is unable togenerate any currents at all, including currents whichpreserves the density, as it gives rise to a purely on-site contribution to d ˆ ρ s ↑ /dt . The superconducting termis therefore not responsible for generating any currents;the condensate only carries them once they exist. Thetreatment of only ordinary quadratic terms is thereforesufficient for calculating the current also in the presenceof ss