Spontaneous generation of vortices by a nonuniform Zeeman field in a two-dimensional Rashba coupled superconductor
aa r X i v : . [ c ond - m a t . s up r- c on ] J u l Spontaneous generation of vortices by a nonuniform Zeeman field in atwo-dimensional Rashba coupled superconductor
A. G. Mal’shukov
Institute of Spectroscopy, Russian Academy of Sciences, Troitsk, Moscow, 108840, RussiaMoscow Institute of Physics and Technology, Institutsky per.9, Dolgoprudny, 141700 Russia andNational Research University Higher School of Economics, Myasnitskaya str. 20, Moscow, 101000 Russia
The interplay of the electron exchange interaction and spin-orbit coupling results in spontaneoussupercurrents near magnetic insulator islands, which are placed on the top of a two-dimensional (2D)superconductor, and whose magnetization is parallel to 2D electron gas. It is shown that in contrastto the well studied situation, where such an effect involves only topologically trivial spatial variationsof the superconducting order parameter, one should take into account supercurrent vortices. Thelatter are spontaneously generated around the island’s boundary of an arbitrary shape and result inscreening of the Zeeman field. This problem has been considered for electrons subject to a strongRashba spin-orbit coupling, including Dirac systems as well. In the latter case vortices can carryMajorana zero modes.
Introduction - The Zeeman interaction of Cooperpair electrons with strong magnetic or exchange fieldsproduces a depairing effect. It destroys the supercon-ducting state [1, 2], or makes the superconducting orderparameter nonuniform [3, 4]. The spin-orbit coupling(SOC) of electrons substantially modifies the effects ofthe Zeeman field. One of the most striking manifesta-tion of the interplay between this field and SOC is themagnetoelectric effect [5, 6], which takes place even ata weak Zeeman interaction. It results in a spontaneoussupercurrent, which was predicted [7–9] to circulate neara magnetic insulator island deposited on the top of a 2Dsuperconductor or a thin film whose electrons are subjectto the Rashba [10] SOC. At the same time, the Zeemanfield is oriented parallel to the 2D electron gas. In spa-tially uniform systems this effect results in a helix spatialvariations of the order parameter [5, 11–16]. In this casethe electric currents, which are created by the magne-toelectric effect and by the order-parameter phase gra-dient, compensate each other, so that the net current iszero. In contrast, for the varying in space Zeeman fieldsuch a compensation is absent, which results in a finitesupercurrent [7–9]. Basically, the same magnetoelectriceffect leads to the anomalous Josephson effect in the socalled φ junctions [17–24] where the Zeeman field andSOC inside the junction induce the spontaneous Joseph-son current, even in the absence of the external phasedifference between superconducting contacts.In thermal equilibrium the spatial variations of the or-der parameter must minimize the free energy of the elec-tronic system which interacts with a Zeeman field. Sofar, the theory was mostly restricted to situations wherethe order parameter has a trivial topology. At the sametime, it is important to understand the role of topologi-cally nontrivial spatial configurations of the order param-eter, such as vortices. As known, vortices may be inducedin type II superconductors under the orbital effect of anexternal magnetic field. Therefore, it is natural to ask,whether the exchange interaction of electron spins with a proximized magnetic insulator could induce vortices.For example, the creation of single edge vortices whichcan carry Majorana zero modes was considered [25] in aspecial case of a Zeeman field produced by a long thinmagnetic wire on top of a 2D Dirac superconductor. Inthis connection it is important to consider magnetic is-lands of various shapes and to find out general condi-tions for the formation of vortices. Moreover, a theory isneeded not only for topological superconductors, but alsofor spin-orbit coupled superconductors having the trivialtopology. Since the number of vortices and their posi-tions strongly depend on the strength and orientation ofthe Zeeman field, their interaction with magnetic excita-tions opens new possibilities to control hybrid system in-volving superconducting and magnetic components, withgreat potential for practical applications.It will be shown below that a chain of vortices may begenerated along the border of a magnetic island, as wellas at a domain wall. At the same time, supercurrents,which are created by such chains, screen the currents thatare produced by the direct magnetoelectric effect and insome cases completely compensate them in the bulk ofan island and outside it, irrespective of its shape.The parameters of such a vortex state are calculatedby minimizing the free energy functional which, in turn,will be obtained using the semiclassical theory of Green’sfunctions [26–29] for disordered superconductors. Twomodels will be considered: a 2D superconductor withthe Dirac electron band and a 2D superconductor withthe parabolic band and the strong Rashba SOC. It is as-sumed that the electronic band splitting, which is asso-ciated with this SOC, is much larger than the supercon-ducting gap, while the latter is larger, or comparable tothe Zeeman field. Within the chosen model the s-wavesuperconductivity in the 2D gas has either an intrinsicorigin, or is induced by the proximity effect. The modelis supported by recent experimental works on 2D super-conductors with strong spin-orbit effects. For example,there are evidences that the Dirac surface band migrates FIG. 1: (Color online) A chain of supercurrent vortices isspontaneously formed around a Zeeman-field island, whichhas a form of a long strip (top), or an arbitrary form (bottomright). It also appears at a domain wall (bottom left). Themagnetization direction is depicted by the arrow from the surface of a three-dimensional topological in-sulator (TI) towards the top of a thin superconductingover-layer and acquires a superconducting gap [30, 31].Self-formation of a thin superconducting Rashba-coupledfilm at the interface of the Pb superconductor and TI wasobserved [32], a 2D superconducting film of Tl and Pbwas grown on a semiconductor substrate [33, 34]. Thesymmetry of the superconducting state in these systemsis not known. Therefore, their relevance to the consid-ered model has yet to be found out.
Free energy - The analysis is based on the free en-ergy functional F which depends on the order-parameterof the superconductor ∆ = | ∆ | exp iθ . This functionalis derived from the one-particle Hamiltonian. For a 2Dsuperconductor the latter has the form H = τ ˆ k m + τ α ( σ x ˆ k y − σ y ˆ k x ) + σ Z − τ µ + τ V imp + Re[∆] τ − Im[∆] τ , (1)where µ is the chemical potential, Z = ( Z x , Z y ) is theZeeman field produced by the exchange interaction withspins of the magnetic island, ˆk = − i∂/∂ r − ( e/c ) τ A , and σ is a vector composed of Pauli matrices σ j ( j = x, y, z ). A is the vector-potential of the magnetic field which isinduced by the supercurrent. The Pauli matrices τ i , i = 1 , ,
3, operate in the Nambu space. The secondterm in Eq.(1) represents the Rashba SOC. In the caseof Dirac electrons, such as electrons on the surface ofa three dimensional (3D) TI, the first parabolic termshould be ignored. In the following, the sufficientlystrong Rashba coupling α will be assumed, so that thesplitting ∆ so = 2 αk F of the electron energy is muchlarger than | ∆ | , where k F is the Fermi wave vector. Theterm V imp in Eq.(1) describes the random impurity po-tential. It will be assumed that the mean elastic scatter-ing rate 1 /τ imp of electrons on impurities is much larger than | ∆ | and Z . At the same time 1 /τ imp ≪ ∆ so .The free energy will be calculated within the semi-classical approximation which is valid when µ ≫| ∆ | , Z, /τ imp and | ∇ Z | ≪ k F Z . The Eilenberger-Usadelformalizm [26–29] may be employed in this parameterrange to derive semiclassical equations for the Greenfunction. Furthermore, due to the relatively large ∆ so these equations may be projected onto two spin-splittedbands [35], or on a single band in the case of TI [36–38].It is important that the low temperature k B T ≪ | ∆ | is assumed. Therefore, the free energy is different fromthe Ginzburg-Landau functional, which is valid near thecritical temperature where | ∆ | ≪ k B T . However, onecan use the fact that | ∆ | strongly varies only near vortexcores, where it turns to zero, while it stays approximatelyconstant in the most of the system. In contrast, the su-percurrent which is associated with the phase gradient ∇ θ decreases slowly outside the core. Especially, in thinfilms the vortex size is large because the screening effectof the induced magnetic field is weak [39]. As long asthe size of the vortex is much larger than the supercon-ductor’s coherence length ξ , a contribution of the core inthe vortex energy may be neglected. Therefore, in thefollowing | ∆ | will be fixed equal to the spatially uniformbulk value ∆ >
0. Such an approximation implies aweak depairing effect of the Zeeman field. According toRef. [35], in a strongly Rashba coupled system this effectcan be ignored at Z ≪ p ∆ /τ imp .At fixed | ∆ | = ∆ and k B T ≪ ∆ one may calculatethe free energy by expanding semiclassical Green’s func-tions over ∇ θ and Z , up to quadratic terms. As shownin the Supplemental material, the free energy is given by F = π N F D ∆ Z d r (cid:16) ( ∇ θ − ec A + κ F ) + βF (cid:17) , (2)where F x = − Z y /v F and F y = 2 Z x /v F . N F and D arethe state density and diffusion constant of 2D electrons,respectively. κ and β are dimensionless parameters. Ina 2D Rashba system D = τ imp ( µ/m + α ) [35], N F = m/ π , β = 2(1 − r ) / (1 + r ) and κ = 2 r/ (1 + r ),where r = α/ p µ/m + α . At the same time, in a Diracsystem D = 2 α τ imp , N F = µ/ πα , β = 0 and κ = 1. Simple examples - Let us assume that Z is finiteinside a strip which is infinite in the x -direction and hasthe width w . Z is directed parallel to the y -axis, as shownin Fig.1. Therefore, F y = 0 and F x >
0. If w ≪ λ s ,where λ s is the effective magnetic screening length, onemay ignore A in Eq.(S3). Note that in thin films thislength is much larger than the London penetration depth[39]. By varying Eq.(S3) with respect to θ we arrive tothe equation ∇ θ + κ ∇ F = 0 (3)In the considered geometry a solution of this equationis θ = 0. The phase is finite only near remote edges ofthe magnetic strip at | x | → ∞ where ∇ F = 0 [7]. Asa result, the integrand in Eq.(S3) becomes ( κ + β ) F .The first term originates from the kinetic energy of thesupercurrent which is induced by the magnetoelectric ef-fect, while the second term represents a direct depairingeffect of the Zeeman field.However, the above solution of Eq.(3) is not unique.One may consider vortex solutions, where θ winds up aninteger of 2 π around singular points. Usually, vorticesincrease the superconductor’s energy. In contrast, in theconsidered situation they can reduce it, thus leading toa spontaneous formation of vortices around the island.Let us consider two regular chains of vortices with theperiod b ≫ ξ , which are placed along the boundariesof the strip at y = − w/ y = w/
2, respectively,where ξ ≪ w . It is assumed that counterclockwise andclockwise vortices with the vorticity 1 are placed on theopposite boundaries. In this case θ in Eq.(S3) is given bythe sum of phases θ i = φ ( r − r i ), where φ ( r − r i ) is theangle between r − r i and the x -axis for the i -th vortex.Since Z is small, its effect on the coordinate dependenceof θ i is neglected. It is convenient to represent ∇ θ as ∇ θ = ∇ θ + d , (4)where ∇ θ is the average of ∇ θ over x , while d is thedeviation from the average. As shown in the Supple-mentary material, ∇ θ y = 0, while ∇ θ x = 2 π/b at − w/ < y < w/
2, and ∇ θ x = 0 elsewhere. By denotingthe first term in Eq.(S3) as F kin and substituting thereEq.(4) we arrive to F kin = π N F D ∆ (cid:18) Lw ( κF − π/b ) + Z d rd (cid:19) . (5)It is seen that at κF = 2 π/b the first term vanishes, whileat 2 πw ≫ b the second term, as shown in the Supplemen-tary material, plays the role of the line energy, because d is not small only within a shell ∼ b near the island’sboundary. As a result, in the leading order with respectto b/w one may take into account only the first termin Eq.(5). Hence, the formation of the vortex ”coat” atthe boundary of the magnetic island becomes energet-ically favorable. The vortex lattice period is given by b = 2 π/κF = πv F /κZ . With decreasing w the secondterm in Eq.(5) must be taken into account. Fig.2 demon-strates that in this case the lattice period b deviates fromits asymptotic ( w → ∞ ) value b = πv F /κZ . Arbitrary island’s shape - The above theory may beextended to the case of an arbitrary island’s shape. Let usassume that the island is much smaller than the magneticscreening length and ignore A in Eq.(S3). In Eq.(3) thephase may be written in the form θ = θ v + θ F , where θ F is associated with the magnetoelectric effect, while θ v originates from vortices. The former turns to 0 at F = 0.In contrast, θ v does not vanish with F and satisfies theequation ∇ θ v = 0. Let us denote ∇ θ + κ F = J . From FIG. 2: (Color online) Dependence of the inverse distancebetween vortices on the strip width, at various strengths ofthe Zeeman field. From top to bottom: ξκZ/v F = 0 . , . . this definition one obtains ∇ × J = ∇ × ( κ F + ∇ θ v ) . (6)Note, that the second term on the right is not zero, be-cause θ v is singular in the vortex center. In contrast,the nonsingular ∇ × ( ∇ θ F ) = 0. Therefore, it does notenter in Eq.(6). The current conservation (Eq.(3)) gives ∇ J = 0. Therefore, J may be represented in the form J = ∇ × B , where B is parallel to the z -axis. By substi-tuting the so expressed J into Eq.(6) we arrive to ∇ B z = − ( ρ F + ρ v ) , (7)where ρ F = κ [ ∇ × F ] z , ρ v = [ ∇ × ∇ θ v ] z and the sub-script z denotes the vector’s z -component. ρ F plays therole of the ”magnetic charge”, which is responsible forthe magnetoelectric effect. At the same time, ρ v may becalled the ”vortex charge” density. Integration of ρ v overa small area enclosing a vortex gives 2 πl , where the inte-ger l is the vorticity. Therefore, ρ v = 2 π P i l i δ ( r − r i ),where r i is the position of the i -th vortex. If the distancebetween vortices is much smaller than the size of an is-land one can average ρ v over r i . It is important that onecan always choose the averaged vortex spatial distribu-tion ρ v in such a way that ρ F + ρ v = 0 in Eq.(7). Asa result B z = 0, as well as the averaged J = 0. Hence,by ignoring fluctuations of ρ v , the first term in Eq.(S3)which can be written as J , turns to zero. Therefore, thefree energy is minimized by an appropriate choice of thevortex distribution. If the Zeeman field is uniform insidean island, the magnetic and vortex charges are concen-trated on its boundary. With decreasing island’s size theline energy given by the second term in Eq.(5) becomesimportant and the employed above macroscopic approachfails. In this case a few vortices can help to reduce thefree energy [25]. Domain wall - If the Zeeman field is not uniform in-side an island, the magnetic charge ρ F = κ [ ∇ × F ] z =2 κ ( ∇ · Z ) /v F is finite there. As an example, let us con-sider a domain wall in Fig.1, where Z is directed in thepositive y -direction for y > y < ρ F = κ [ ∇ × F ] z = − κF x δ ( y ). Hence, in Eq.(7), ρ F + ρ v turns to zero,if a row of counterclockwise vortices with the density κF x /π is placed along the domain wall. Similar to pre-viously discussed examples, such a vortex distributionminimizes the free energy. If the size of the island in the y -direction is finite, the magnetic charges with the linedensity ρ F = κF x / π reside on both boundaries togetherwith two compensating chains of vortices. Magnetic screening - Usually, the magnetic screen-ing is produced by supercurrents which induce the mag-netic field directed in the opposite direction to the ex-ternal field. As a result, the magnetic field can notpenetrate deep into superconductors. In 3D samples itdecreases exponentially together with the screening cur-rent. However, in thin films such a screening effect isweak. It leads only to a r − reduction of the current at r & λ s [39]. A similar situation takes also place in thecase when the supercurrent and vortices are induced bythe considered here magnetoelectric effect of the Zeeman(exchange) field. The screening effect must be taken intoaccount for magnetic islands whose size is comparablewith λ s . The magnetic field is represented in Eq.(S3) by A . This vector-potential satisfies the Maxwell equation ∇ A = − πc δ ( z ) j . (8)For a disordered superconductor the 2D electric cur-rent density is given by j = eπ | ∆ | N F D ( ∇ θ + κ F − e/c ) A ). By substituting this current in Eq.(8) thevector-potential may be expressed in terms of J = ∇ θ + κ F . As a result, instead of J − e/c ) A , that enters inthe free energy Eq.(S3), we obtain the expression whoseFourier transform is given by J q − e/c ) A q = J q qq + k s , (9)where k s = 4 π e | ∆ | N F D/c . It is seen that at the large q we arrive to an unscreened expression ( A = 0), whileat q . k s a crossover occurs from 1 /r to 1 /r spatialdependence of the vortex current [39]. A typical k − s islarger than 1 µ m.Let us take into account in Eq.(S3) only ∇ θ and A which are averaged over vortex positions. The to-tal energy F tot must include also the energy density( rot A ) / π of the induced magnetic field. By expressing A in terms of J and by integrating the magnetic energyover z we arrive to F tot = π N F D ∆ X q | J q | q + | q × J q | k s q ( q + k s ) . (10) It follows from this equation that F tot turns to zero if J q = 0. It is the same condition which regulates the for-mation of vortices in the absence of the screening effect.At the same time, the line energy, which is given by thesecond term in Eq.(5), was ignored in Eq.(S5). It is rea-sonable to expect, however, that the screening effect cannot modify the surface energy, if the distance d betweenvortices is much smaller than k − s , because the surfaceterm in Eq.(5) is determined by a thin shell ∼ d nearthe boundary. Anyway, the screening effect is not impor-tant for large enough islands where the surface energyconstitutes only a small fraction. Conclusion - In conclusion, the spatial configurationof the superconducting order parameter in a 2D spin-orbit coupled superconductor with an inhomogeneousZeeman interaction has been considered. The Zeemanfield was assumed to be directed parallel the 2D elec-tron gas. Such a field may be produced by the exchangeinteraction of 2D electrons with polarized electrons ofa magnetic insulator film which is placed on top of thesuperconductor. It is shown that a chain of vortices isspontaneously formed on the island’s boundary, as well asalong a domain wall, even at a weak Zeeman interaction,which may be much less than the magnitude of the orderparameter. The distance between vortices increases for aweaker Zeeman field and decreases down to the coherencelength when this field becomes comparable to the orderparameter. This situation is quite different from the spe-cial island’s geometry which was considered in Ref.[25],where the Zeeman field must be sufficiently strong to cre-ate a pair of vortices with localized Majorana zero modesat the edges of a long magnetic wire. It is shown that theformation of the vortex ”coat” takes place for sufficientlylarge islands having an arbitrary shape. It is importantthat such a coat suppresses the supercurrent which is in-duced by the Zeeman interaction inside the island andoutside it, except for a thin shell along the boundary, ora domain wall.This work has been focused on an analysis of equi-librium properties of the considered hybrid system. Onthe other hand, one would expect many interesting non-stationary phenomena associated with interaction of vor-tices with collective spin excitations of a magnetic insula-tor film. Furthermore, in the case when the superconduc-tor has a nontrivial topology Majorana bound states maylocalize at vortices and there are various possibilities tocontrol their positions and mutual distance by changingthe magnetization direction.
Acknowledgements - The work was partly supportedby the Russian Academy of Sciences program ”Actualproblems of low-temperature physics.” [1] B. S. Chandrasekhar, Appl. Phys. Lett. , 7 (1962) [2] A. M. Clogston, Phys. Rev. Lett. , 266 (1962)[3] P. Fulde and R. A. Ferrel, Phys. Rev. , A550 (1964)[4] A. I. Larkin and Yu. N. Ovchinnikov, Zh. Eksp. Teor. Fiz. , 1136 (1964) [Sov. Phys. JETP , 762 (1965)][5] V. M. Edelstein, Sov. Phys. JETP , 1244 (1989)[6] S. K. Yip, Phys. Rev. B , 144508 (2002)[7] A.G. Mal’shukov, Phys. Rev. B , 054511 (2016).[8] S. S. Pershoguba, K. Bj¨ornson, A. M. Black-Schaffer, andA. V. Balatsky, Phys. Rev. Lett. , 116602 (2015).[9] K. M. D. Hals, Phys. Rev. B , 134504 (2017)[10] Yu. A. Bychkov and E. I. Rashba, J. Phys. C , 6039(1984).[11] V. P. Mineev and K. V. Samokhin, Zh. Eksp. Teor. Fiz. , 747 (1994) [Sov. Phys. JETP 78, 401 (1994)][12] V. Barzykin and L. P. Gor’kov, Phys. Rev. Lett. ,227002 (2002)[13] D. F. Agterberg, Physica C , 13 (2003)[14] R. P. Kaur, D. F. Agterberg, and M. Sigrist, Phys. Rev.Lett. , 137002 (2005)[15] D.F. Agterberg and R.P. Kaur, Phys. Rev. B , 064511(2007)[16] O. Dimitrova and M.V. Feigel’man, Phys. Rev. B ,014522 (2007)[17] I. V. Krive, A. M. Kadigrobov, R. I. Shekhter and M.Jonson, Phys. Rev. B , 214516 (2005)[18] A. Reynoso, G.Usaj, C.A. Balseiro, D. Feinberg,M.Avignon, Phys. Rev. Lett. , 107001 (2008)[19] A. Zazunov, R. Egger, T. Martin, and T. Jonckheere,Phys.Rev. Lett. , 147004 (2009)[20] A. G. Mal’shukov, S. Sadjina, and A. Brataas, Phys. Rev.B , 060502 (2010)[21] J.-F. Liu and K. Chan, Phys. Rev. B , 125305 (2010)[22] T. Yokoyama, M. Eto, Y. V. Nazarov, Phys. Rev. B ,195407 (2014)[23] F. Konschelle, I. V. Tokatly and F. S. Bergeret, Phys.Rev. B ,125443 (2015)[24] A. Assouline, C. Feuillet-Palma, N. Bergeal, T. Zhang, A.Mottaghizadeh1, A. Zimmers, E. Lhuillier, M. Marangolo,M. Eddrief, P. Atkinson, M. Aprili, H. Aubin, NatureCommunications 10, 126 (2019)[25] A. G. Mal’shukov, Phys. Rev. B , 134514 (2020)[26] A. I. Larkin, and Y. N. Ovchinnikov, Zh. Eksp. Teor. Fiz. , 2262 (1968) [Sov. Phys. JETP , 1200 (1965)].[27] G. Eilenberger, Z.Phys. , 195 (1968)[28] J. Rammer, H. Smith, Rev. Mod. Phys. , 323 (1985)[29] J. M. Serene, D. Rainer, Phys. Rep. , 221 (1983)[30] N. Sedlmayr, E. W. Goodwin, M. Gottschalk, I. M. Day-ton, C. Zhang, E. Huemiller, R. Loloee, T. C. Chas-apis, M. Salehi, N. Koirala, M. G. Kanatzidis, S. Oh,D. J. Van Harlingen, A. Levchenko, and S. H. Tessmer,arXiv:180512330[31] C. X. Trang, N. Shimamura, K. Nakayama, S. Souma,K. Sugawara, I. Watanabe, K. Yamauchi, T. Oguchi, K.Segawa, T. Takahashi, Yoichi Ando, and T. Sato, Nat.Communication , 159 (2020)[32] M. Bai, F. Yang, M. Luysberg, J. Feng, A. Bliesener,G. Lippertz, A. A. Taskin, J. Mayer, and Y. Ando,arXiv:1910.08331 (2019)[33] A. V. Matetskiy, S. Ichinokura, L. V. Bondarenko, A. Y.Tupchaya, D. V. Gruznev, A. V. Zotov, A. A. Saranin, R.Hobara, A. Takayama, and S. Hasegawa, Phys. Rev. Lett. , 147003 (2015)[34] T. Sekihara, R. Masutomi, and T. Okamoto, Phys. Rev.Lett. , 057005 (2013) [35] Manuel Houzet and Julia S. Meyer, Phys. Rev. B ,014509 (2015).[36] A. Zyuzin, M. Alidoust, and D. Loss, Phys. Rev. B ,214502 (2016).[37] I. V. Bobkova, A. M. Bobkov, A. A. Zyuzin, andM.Alidoust, Phys. Rev. B , 134506 (2016)[38] Henning G. Hugdal, Jacob Linder, and Sol H. Jacobsen,Phys. Rev. B 95, 235403 (2017)[39] J. Pearl, Appl. Phys. Lett, , 65 (1964). Supplemental Materials
S1. CALCULATION OF THE FREE ENERGY
The free energy in Eq.(2) of the main text is given by the second-order expansion with respect of ∇ θ, A , Z and thefirst order expansion in their cross-products. It is convenient to transform Hamiltonian Eq.(1) with the help of theunitary operator exp( iτ θ/ H = τ m (ˆ k + τ A ) + τ α [ σ x (ˆ k y + τ A y ) − σ y (ˆ k x + τ A x )] + σ Z − τ µ + τ V imp + ∆ τ , (S1)where A = ∇ θ/ − ( e/c ) A and ˆ k = − i∂/∂ r . ∆ is chosen real and positive. It depends strongly on coordinates closeto the vortex cores. However, sufficiently far from them, at the distance much larger than the coherence length thisdependence is weak. At small A and Z one may neglect their effect on ∆ and take the latter equal to ∆ , which isthe unperturbed order parameter. ˜ H can be represented in the form ˜ H = ˜ H + δH , where ˜ H is the unperturbedHamiltonian and the interaction Hamiltonian is given by δH = 12 m (ˆ k A + A ˆ k + τ A ) + α ( σ x A y − σ y A x ) + Z σ (S2)Let us define ˜ H λ = ˜ H + λδH . Then, the correction to the free energy can be represented as [1] δ F = − k B T X ω n Z dλ Z d r Tr[ δHG λ,ω n ( r , r ′ )] r = r ′ , (S3)where the trace is taken over Nambu and spin variables and G λ,ω n ( r , r ′ ) is the imaginary-time Green function whichis calculated with the Hamiltonian ˜ H λ . This function should be calculated up to the first-order with respect to δH .First we consider the Rashba spin-orbit coupling (SOC) in a topologically trivial 2D system. It will be assumedthat the splitting of electron energy bands ∆ so , which is caused by the SOC, is much larger than ∆ and the elasticscattering rate 1 /τ . In this case it is convenient to transform ˜ H to the helical basis. In this basis, instead of spinprojection, the quantum states are labeled by the helicity index ν = ± . Accordingly, there are two bands with theenergies ξ ν k = k / m + ναk − µ . Each band crosses the Fermi energy ( ξ ν k = 0) at the corresponding wave-number k νF .The Fermi velocity v F = p µ/m + α is the same in both bands, while densities of states N νF = ( m/ π )(1 − να/v F )are different. Therefore, there are different electron scattering rates 1 /τ ± = 2 πN ± F | V imp | , where | V imp | is expressedin terms of the Born scattering amplitude for random uncorrelated impurities [1].The calculation of δ F in Eq. (S3) can be simplified, if the semiclassical Green function is used instead of G λ,ω n ( r , r ′ ).The semiclassical approximation is valid when the Fermi wave length ∼ k − F is much less then other relevant parametershaving the dimension of length. In the considered case of the strong SOC there are two helical bands. Therefore,semiclassical Green functions are defined for each band, while their nondiagonal terms, which correspond to mixingbetween bands, may be ignored at the large ∆ so = | ξ + k − ξ − k | . The semiclassical functions are defined in the followingway: g νλ,ω n ( n k , r ) = i τ π Z dξ ν k G λ,ω n ( k , r ) , (S4)where n k = k F /k F and G λ,ω n ( k , r ) is obtained from G λ,ω n ( r , r ′ ) by Fourier transform with respect to r − r ′ , and bysetting ( r + r ′ ) / → r . Note, that the semiclassical function depends only on the direction of k , which in turn is fixedon the Fermi line. Semiclassical equations for g νλ,ω n ( r ) in the presence of the strong Rashba SOC and magnetic fieldhave been derived in Ref.[35] of the main text. These results will be employed for the calculation of the free energyEq. (S3). In terms of quasiclassical functions Eq. (S3) can be written in the form δ F = − i πk B T X ω n Z d n k Z d r Z dλ Tr[ τ ( v F n k A + Z σ ) g λ,ω n ( n k , r )] , (S5)where g λ,ω n ( n k , r ) = N + F g + λ,ω n ( r )(1 + σ × n k ) + N − F g − λ,ω n ( r )(1 − σ × n k ) . (S6)It is easy to see that due to the angular averaging only the odd in k part of g contributes in the integral Eq. (S5).In the case of the strong impurity scattering the Green function is almost isotropic. Therefore, it can be representedas g νλ,ω n ( n k , r ) = g νλ,ω n ( r ) + n k g νλ,ω n ( r ) [2]. The second term, which represents the anisotropic correction, is small, sothat | g | ≪
1. From Eqs.(S5-S6) one thus obtains δ F = − iv F πk B T X ω n Z d r Z dλ Tr τ [ τ ( N + F g + λ,ω n + N − F g − λ,ω n )2 A + τ ( N + F g + λ,ω n − N − F g − λ,ω n ) F ] , (S7)where F x = − Z y /v F and F y = 2 Z x /v F . As shown in Ref.[35] of the main text, the functions g ± λ,ω n can be expressedin terms of the isotropic function ˜ g λ,ω n ( r ) = (1 / g + λ,ω n ( r ) + g − λ,ω n ( r )). The corresponding equations have the form N + F g + λ,ω n + N − F g − λ,ω n = − v F N F τ ˜ g λ,ω n (cid:18) τ ∇ + + 1 τ − ∇ − (cid:19) ˜ g λ,ω n ,N + F g + λ,ω n − N − F g − λ,ω n = − v F N F τ ˜ g λ,ω n (cid:18) τ ∇ + ˜ g λ,ω n − τ − ∇ − ˜ g λ,ω n + iλτ + τ − F [ τ , ˜ g λ,ω n ] (cid:19) , (S8)where N F = ( N + F + N − F ) /
2, so that τ ± /τ = N F /N ± F . The gradients ∇ ± are given by ∇ ± ∗ = ∇ ∗ + λi A ± F )[ τ , ∗ ] . (S9)As can be seen from Eq. (S7) it is sufficient to calculate the linear in A and Z correction to ˜ g . The equation forthis function has been derived in Ref.[35]. An analysis of this equation shows that the linear correction is absent, ifone takes into account Eq. (3) of the main text. Therefore, one may set ˜ g equal to the zero-order Green function˜ g λ,ω n = ( τ ω n + τ ∆ ) / p ω n + ∆ and ignore the gradient in the right-hand side of Eq. (S9). By substituting sosimplified Eq. (S9) into Eq. (S8) and Eq. (S7) we arrive to the free energy given by Eq. (2) of the main text. Thesum over ω n was calculated at the low temperature k B T ≪ ∆ .In the case of a Dirac system with a single Dirac point there is only a single helical band, which should be takeninto account at sufficiently large µ . In this respect, the above algebra is simplified, while the corresponding equationfor the semiclassical Green function has been calculated in Refs.[36,37] of the main text. The free energy is given byEq. (2). However, in contrast to the Rashba superconductor the second term in Eq. (2) is absent. S2. VORTEX COAT ON A ZEEMAN ISLAND BOUNDARY
The phase θ of the order parameter is induced by two regular chains of vortices on two opposite edges of a strip at y = − w/ y = w/
2, respectively, where w is the strip width. It is assumed that counterclockwise and clockwisevortices with coordinates r + i = ( bi + s, − w/ ,
0) and r − i = ( bi, w/ , b is the chain period, i = 1 , , ... , and s denotes the relative shift of the chains in the x -direction. In this case ∇ θ inEq.(2) of the main text may be written as ∇ θ = X i (cid:18) − ˆz × ( r − r + i )( r − r + i ) + ˆz × ( r − r − i )( r − r − i ) (cid:19) , (S10)where ˆz is a unit vector in the z -direction. ∇ θ can be averaged over x . Let us denote the so averaged ∇ θ as ∇ θ .From Eq.(S10) it is easy to calculate it as ∇ θ = − πb ˆx θ ( x + w/ θ ( w/ − x ), where ˆx a unit vector in the x -direction.Note, that this average is finite only inside the strip. Therefore, by choosing the appropriate lattice period b , one maycompensate the Zeeman term in Eq.(2) of the main text. However, it is necessary to take into account the deviation δ ∇ θ of ∇ θ from its average value. By taking the sum over i in Eq.(S10) this deviation can be written in the form δ ∇ θ = πb X m =0 e ig m x (cid:16) e −| g m || y − w/ | V + m − e −| g m || y + w/ | e − ig m s V − m (cid:17) , (S11)where g m = 2 πm/b , V ± x =sign( y ∓ w/
2) and V ± y = i sign( g m ). It is seen that δ ∇ θ decreases fast when the distancefrom the edges increases. At 2 πw ≫ b it contributes in the second term of Eq.(7) only within a thin shell ∼ b/ π nearthe edges. This term plays a role of the the surface energy. Let us denote it as F s . By substituting Eq.(S11) in thesecond term of Eq.(7) we obtain F s = π N F D ∆ Lg (cid:20) gξ + ln (cid:0) − e − wg + isg (cid:1) + ln (cid:0) − e − wg − isg (cid:1)(cid:21) , (S12)where g = 2 π/b . Since gξ ≪
1, the superconductor coherence length ξ plays the role of the logarithmic cutoff nearvortex cores. If the width of the strip is much larger than the distance between vortices the first term in Eq.(S12)dominates and the surface energy does not depend on the strip width. By minimizing the free energy Eq.(7) oneobtains b , as a function of w and F . The result is shown in Fig.2 of the main text. [1] A. A. Abrikosov, L. P. Gor’kov, and I. E. Dzyaloshinskii, Methods of Quantum Field Theory in Statistical Physics (Dover,New York, 1975)[2] K. D. Usadel, Phys. Rev. Lett.25