Topological superconductivity with orbital effects in magnetic skyrmion based heterostructures
Maxime Garnier Andrej Mesaros Pascal Simon





[Corresponding author:  ]maxime.garnier1@u-psud.fr
Laboratoire de Physique des Solides, UMR 8502, CNRS,
Université Paris-Sud, Université Paris-Saclay, 91405 Orsay, France

Laboratoire de Physique des Solides, UMR 8502, CNRS,
Université Paris-Sud, Université Paris-Saclay, 91405 Orsay, France

Laboratoire de Physique des Solides, UMR 8502, CNRS,
Université Paris-Sud, Université Paris-Saclay, 91405 Orsay, France





Proximitizing magnetic textures and s-wave superconductors is becoming a platform for engineering topological superconductivity and Majorana fermions by the means of exchange processes. However, the consequences of orbital effects have not yet been fully taken into account. In this work, we investigate the magnetic skyrmion texture-induced orbital effects using a Ginzburg-Landau approach and clarify the conditions under which they can induce superconducting vortices. These orbital effects are then included in Bogoliubov-De-Gennes theory containing the exchange interaction, as well as superconducting vortices (when induced). We find that the topological phase is largely stable to all investigated effects, increasing the realistic promise of skyrmion-superconductor hybrid structures for realization of topological superconductivity.

Topological superconductivity with orbital effects in magnetic skyrmion based heterostructures
    Pascal Simon
    December 30, 2023
==============================================================================================







§ INTRODUCTION


The search for topological superconductors and  Majorana fermions has mostly relied on an approach combining conventional s-wave superconductivity, spin-orbit coupling and magnetism to generate p-wave pairing. Even though great advances have been made using this approach using semiconducting wires 
(see Ref. Lutchyn:2018 for a review) or magnetic atoms combined with 
superconducting substrates,<cit.> these experiments remain highly challenging so that it may be worth considering removing one of the ingredients. For example, one may consider removing the spin-orbit coupling since the exchange interaction between conduction electrons and a magnetic texture induces an effective spin-orbit interaction.<cit.> This seems to be a viable platform for inducing topological superconductivity.<cit.> In particular, particle-like topological spin configurations known as magnetic skyrmions<cit.> have recently been highlighted as prime candidates<cit.> with interesting prospects for manipulation due to their high degree of controllability.<cit.> However, in this context, the magnetic orbital effects have not yet been fully taken into account. 

Generally, systems in presence of supercurrents and spin-orbit coupling may support topologically non-trivial phases without a spin coupling of the Zeeman form.<cit.> This suggests that orbital effects are sufficient to create topological phases in magnet-superconductor heterostructures. In a more involved scenario, the Zeeman-form coupling (induced by the exchange interaction between skyrmion and superconductor) combined with the magneto-electric effect due to the spin-orbit coupling leads to the appearance of vortices,<cit.> which induce supercurrents and topological superconductivity.<cit.> If one now considers alternatives based on removing the explicit spin-orbit coupling, which is a key coupling ingredient in previous scenarios, one finds that topological superconductivity is also predicted with only the Zeeman-form coupling induced by the exchange interaction with the skyrmion.<cit.> However, even without explicit spin-orbit coupling, one should consider the effect of supercurrents, which may appear only the orbital effects generated by the skyrmion<cit.> (since the magneto-electric coupling vanishes together with the spin-orbit coupling). The orbital effects due to skyrmions cannot be generally neglected, as we argue below, and their effect on topological superconductivity has not been analyzed so far. As a consequence, it is not clear if vortices are to be expected and whether topological superconductivity persists when both exchange and orbital effects of the skyrmion are included in a superconductor without spin-orbit coupling.

In conventional type-II superconducting thin films, the effective penetration depth (or Pearl length) λ_ eff can be orders of magnitude larger than the film thickness d.<cit.> This implies that the screening currents are weak and may become negligible. Previous works<cit.> have used this observation to set to zero the magnetic vector potential in the superconductor. However, we argue that since the screening is weak, the magnetic field generated by the skyrmion penetrates the superconductor almost unaltered, and is thus not necessarily negligible. Therefore in contrast to previous works, we include the vector potential generated by the skyrmion as an orbital effect on the electrons in the superconductor. To fully understand the phase diagram, we consider the magnetic exchange as an independent effect on the electrons, since this term could be experimentally tuned by an insulating non-magnetic layer, which prevents the hopping of electrons between the skyrmion and the superconductor while not affecting the strength of the vector potential.

The aim of the present work is twofold. First, we investigate the orbital effects of the magnetic field created by the skyrmion on a conventional type-II superconductor without spin-orbit coupling and clarify the conditions of existence of vortices in a Ginzburg-Landau framework. Second, using the Bogoliubov-de-Gennes formalism we address the robustness of the topological phase induced by the exchange interaction to the orbital effects and the possible presence of vortices. We find that the exchange-induced topological phase is robust to the inclusion of orbital effects, even in the presence of vortices. This contributes to making the skyrmion-superconductor heterostructure an even more promising platform for the realization of topological superconductivity.

The article is organized as follows. In Sec. <ref>, we present the system under study, the hypotheses made on the superconductor, derive an expression for the magnetic field generated by a single skyrmion and set up our Ginzburg-Landau analysis. We establish the phase diagram of the superconductor in Sec. <ref> and study the implications on the topological superconducting phase in Sec. <ref>.





§ SETUP



The system under investigation is composed of a conventional type-II  superconducting film without spin-orbit coupling (such as Al) in proximity to an insulating magnetic thin film harboring a skyrmion (such as Cu_2 O Se O_3) <cit.> as shown in fig: setup. We denote by d the thickness of the superconducting layer and by h the thickness of the magnet.

We focus on insulating magnets so as to minimize feedback mechanisms on the magnet due to the superconductor. <cit.> Even though magnetic skyrmions exist in two types, namely (fig: setup) and Bloch (as in Cu_2 O Se O_3),<cit.> this distinction does not modify the physics described here as will be clarified.

Since in-plane critical fields of superconducting films are usually larger than that of bulk superconductors<cit.>, in analyzing the magnetostatics of our system we focus on the magnetic field component perpendicular to the superconductor (z axis). The orbital coupling to the superconductor is hence through the vector potential which satisfies ∇  A( r, z) = B_z( r)ẑ, where r labels the position in the plane. We take B_z to be the magnetic induction created by a lone skyrmion, thereby neglecting any feedback effect of the superconductor on the magnetic material. This assumption is consistent with neglecting the screening currents and their fields. In the limit where the skyrmion is confined to a plane, the magnetic induction B_z( r, z) it creates becomes equal to its magnetization component m_z( r) at least near the plane (see Ref. QinWang_SKL_2018 and Supplementary Material (SM) <ref>). More precisely, we find that  B_z decays away from the plane on a lengthscale given by the radius  of the skyrmion, and for z smaller than this lengthscale we can define A by ∇  A = μ_0   m_z( r)ẑ where μ_0 is the vacuum magnetic permeability. The radius  is defined by the magnetization profile m( r) of the skyrmion texture in polar coordinates r = (r, θ) 

    m( r)  =  M
    [ sin f(r) cos(θ + γ); sin f(r) sin(θ + γ);            cos f(r) ]

The skyrmion is characterized by the radial winding number p ∈ℕ that counts the number of spin flips as one moves away from the core r = 0 along the radial direction. f(r) is the radial profile of the skyrmion that we choose to be f(r) = k_s r for r ≤, where we have introduced k_s = π/λ_s with the spin-flip length λ_s = /p. The global angular offset γ called helicity<cit.> allows to describe both (γ = 0, π) and Bloch (γ = ±π/2) skyrmions. The norm M of m( r) defines the saturation magnetization of the magnet.

Focusing on the (γ = 0) case, the magnetization eq: skyrmion_def can be written in cylindrical coordinates

    m( r) = M (sin(k_s r)  u_r + cos(k_s r) u_z )

A suitable vector potential is

    A( r) =
    μ_0 M/k_s(cos(k_s r)k_s r + sin(k_s r))   u_θ

which indeed gives B_z = μ_0 M cos(k_s r). Note that although this expression was derived for a skyrmion, the magnetic induction generated by a Bloch skyrmion is qualitatively the same (see Ref. QinWang_SKL_2018 and SM <ref>).

For conventional superconductors, the Ginzburg-Landau (GL) free energy functional F is

    F = B_c^2/μ_0∫ d^3 r[-|Δ̃|^2 + 1/2|Δ̃|^4 + ξ^2|D̂Δ̃|^2 + B^22B_c^2]

where Δ̃ is the superconducting order parameter normalized by its thermodynamic value in the absence of fields or gradients, B_c is the thermodynamic critical field, D̂ = -i∇ +2e/ħA is the covariant derivative with 2e the charge of Cooper pairs, and ξ is the superconducting coherence length.<cit.>
Finally, B = ∇ A is the magnetic induction inside the superconductor.  The free energy eq: GL_F is measured with respect to the normal-state free energy. As we are interested in the behavior of the superconductor in an external magnetic field, the correct thermodynamic potential to consider is the Gibbs free energy G =  F - ∫ d^3 r  H( r)· B( r) where H is the magnetic field.
As mentioned above, the typical electromagnetic response of a type-II superconducting film in a homogeneous magnetic field occurs on a lengthscale λ_ eff much larger than any other lengthscale in the problem. As a consequence, the inhomogeneous response on lengthscales of the order of the skyrmion's spin-flip length λ_s cannot be inferred easily. As outlined above we neglect the screening currents so that we can approximate B by μ_0 H inside the superconductor.<cit.>
With this assumption, G simply reduces to F and is given by eq: GL_F where B is now the skyrmion-generated induction given below eq: SK_vect_pot.

We further suppose that all quantities are independent on the z coordinate, which is valid as long as the superconducting thin film thickness d is smaller than the skyrmion radius , so that in this regime the thickness d can be factored out of the free energy.

The final ingredient of our model is a superconducting vortex. To establish its presence or absence we use a standard vortex ansatz<cit.>  Δ̃( r) = Δ̃_α(r)  e^i αθ where Δ̃_α(r) =tanh^|α|(r/ξ) and α∈ℤ is the phase winding of the vortex. Note that even if this ansatz corresponds to an Abrikosov vortex (in  a bulk sample), we expect that the details of the ansatz don't matter much as long as the order parameter amplitude decays on a lengthscale ξ and vanishes at the vortex core. 
Including the vector potential due to the skyrmion eq: SK_vect_pot and the vortex ansatz, the total free energy G_ tot is:


    G_ tot = B_c^2/μ_02π d/k_s^2∫dr̃ r̃[-|Δ̃_α|^2 + 1/2|Δ̃_α|^4 +  k_s^2 ξ^2 (∂_r̃|Δ̃_α|)^2  + k_s^2 ξ^2ν^2/r̃^2|Δ̃_α|^2(α/ν + cos(r̃) + r̃sin(r̃))^2 +p^2π^2/8(μ_0 M/B_c)^2
        
      ]


where we have defined the dimensionless coordinate r̃ = k_s r and the parameter ν≡2πμ_0 M/k_s^2 ϕ_0 where ϕ_0 = h/2e is the superconducting flux quantum. The ν parameter can be interpreted as the ratio |ϕ_S| / (2ϕ_0) of the flux of the skyrmion |ϕ_S| through a disc of radius λ_s (the single spin flip distance) to the normal-metal flux quantum 2ϕ_0.




§ PHASE DIAGRAM OF THE SUPERCONDUCTOR WITHOUT EXCHANGE EFFECTS



The free energy, eq: adim_G, depends on three parameters:


  (i) the ratio of the characteristic skyrmion flux and the flux quantum ν = |ϕ_S|/2 ϕ_0,

  (ii) the ratio of the superconducting coherence length to the spin-flip length of the skyrmion k_s ξ=p πξ/,

  (iii) the ratio μ_0 M / B_c.

Given some values of these parameters, our strategy is to numerically compute the free energy for different windings α (including α=0 for absence of vortex) and find the value that gives the lowest free energy. If the free energy is negative, the normal state is realized instead of the superconducting one.
Furthermore, as the magnetic energy term (last term in eq: adim_G) doesn't depend on α, we start by neglecting it. 


In fig: PhDiag we show the phase diagram obtained for a p=4 skyrmion  where we have introduced a short-distance cutoff k_s l for the dimensionless variable r̃ to deal with the logarithmic divergence in the α = 0 case. (This cutoff has only minor consequences on our results, see SM <ref>, and its value is chosen to be comparable to the lattice spacing in the Bogoliubov-de-Gennes Hamiltonian of Sec. <ref>).


The phase diagram shows that there exists a superconducting phase without vortex (red color) for a relatively large range of superconducting coherence lengths and skyrmion fluxes ν≲ 1.5. For sufficiently strong skyrmion flux and small ξ/, there exist phases with a vortex.
We have also checked that including the magnetic energy term and varying the size of the skyrmion don't affect qualitatively our results, see SM <ref> and <ref>. 
We thus conclude that vortices are not always expected when proximitizing a conventional superconductor without spin-orbit coupling by a magnetic skyrmion.



§ IMPLICATIONS FOR TOPOLOGICAL SUPERCONDUCTIVITY




We now turn to the consequences of the orbital effects on the topological superconducting phase induced by the exchange interaction.<cit.> As shown in the previous section, the inclusion of the orbital effects due to the skyrmion has two consequences: the first is that electrons see a magnetic vector potential and the second is that vortices may be present.




 §.§ Bogoliubov-de-Gennes setup


In the Bogoliubov-de-Gennes (BdG) formalism, the total Hamiltonian H describing the electrons can be written H = 1/2∫ d r Ψ^†( r)  ( r)  Ψ( r) where $̋ is the BdG Hamiltonian. Throughout the article, we work in the Nambu basisΨ^†(r) = (ψ_^†(r) , ψ_^†(r),   ψ_(r) , -ψ_(r) )whereψ_σ^†(r) is the field operator creating an electron with spin projectionσ= ,  at positionr = (r, θ)in two dimensions (2D). 
Following the minimal coupling procedurep̂ →p̂ +e A, the 2D BdG Hamiltonian in the presence of both orbital and exchange effects reads

    ( r)  = (1/2m(p̂ + e Aτ_z)^2 - μ) τ_z + J	 σ· m( r)  + Δ_0  τ_x

whereτ_iandσ_i,i = x,  y,  z, are Pauli matrices acting in particle-hole and spin space, respectively,p̂ = -iħ∇is the momentum operator,mis the effective mass,μthe chemical potential,the bares-wave pairing,Jis the exchange interaction, and orbital effects are due to the vector potential  A. The vectorm(r)is the skyrmion texture as parametrized in eq: skyrmion_def and we set the helicityγ= 0since it can be unitarily removed from the Hamiltonian.<cit.>  Hereafter, the exchange interactionJincludes the saturation magnetizationMof the magnet. We emphasize that the strengths of orbital effects and exchange interactions can be tuned independently. Experimentally, the couplingJcan be reduced by an insulating layer between the magnet and the superconductor.

The BdG Hamiltonian eq: BdG_Skyrmion_full has a generalized rotation symmetry and total angular momentum operator which, in absence of any vortices, readsJ_z = L_z + 1/2ħσ_zwhere the orbital angular momentum readsL_z = -i ħ∂_θ. The eigenvalues ofJ_zprovide a quantum numberm_J(in units ofħ) labeling independent blocks of the BdG Hamiltonian. Angular momenta are half-odd-integer as required by the singlevaluedness of the wavefunction. We discretize the Hamiltonian according tor →r_j = jawith lattice constantachosen as the length unit (a ≡1) and the hopping parametert = ħ^2/(2ma^2)is chosen as the energy unit (t ≡1, see SM <ref> and the Methods section of Ref. Garnier_SkTopoSC for additional technical details). For completeness, the lengthy expression of the discretized HamiltonianH_m_J^(1)corresponding to eq: BdG_Skyrmion_full is given in SM <ref>. In all computations we use hard wall boundary conditions at the skyrmion's edge.




 §.§ Orbital effects without exchange and vortices

Neglecting the exchange interaction (J = 0in eq: BdG_Skyrmion_full), the relevant angular momentum operator isJ_z = L_z. The discretizedm_L-dependent HamiltonianH_m_L^(2)is given in SM <ref>, withm_L ∈ℤthe eigenvalue ofL_z(in units ofħ). The vector potential contributes two terms: (1) a space-dependent renormalization of the chemical potential  theA^2term, and (2) a term∝1that depends on both space and angular momentum.

fig: spec_orbonly below shows the numerically obtained BdG spectum contrasting the cases of weak and strong orbital effects,ν=1(fig: spec_orbonlya) andν= 5(fig: spec_orbonlyb). In the latter case, we neglect the predicted appearance of vortices (see fig: PhDiag) to better isolate the orbital effects. The electronic local density of states (LDoS) along ther = 0 →r = radial line in theν=5case is shown in fig: orbexcha.


We observe that the bare superconducting gapis preserved when orbital effects are added. This is somewhat expected since spin degeneracy is not lifted in the Hamiltonian. In theν= 5case, a few localized states come down in energy but stick to the top (resp. bottom) of the gap form_L < 0(resp.m_L >0). We find that these states are largely localized (not necessarily near the core of the skyrmion) resembling an analogue of Caroli-de-Gennes-Matricon states.<cit.> In fact, similar in-gap spectrum due to purely orbital effects was found in a self-consistent calculation<cit.> of the magnetic field due to a preformed normal metal dot in a superconductor (without vortex).





 §.§ Orbital effects and exchange interaction

When restoring the exchange interaction, the 2D BdG Hamiltonian in eq: BdG_Skyrmion_full has total angular momentumJ_z = L_z + 1/2ħσ_zand the full radial tight-binding HamiltonianH_m_J^(1)is given in SM <ref>.

To assess the consequences of orbital effects, we first return to the model with only exchange interaction (ν= 0), which leads to topological superconductivity,<cit.> see fig: spec_orbonlyc. For completess, let us recall briefly the properties of this topological superconductor. First, the effective gap is ofp-wave origin with an amplitude given byΔ_eff = π/λ_s Δ_0/J √(J + μ) where all energies are in units oftandλ_sis expressed in units ofa. This evaluates toΔ_eff ≈7% which is consistent with the numerical data. Within the gap, there are two types of states: namely a nearly flat band located at the edge of the skyrmion together with dispersing states located near the core. These states are attributed to impurity-like states induced by the discretized magnetic texture. In the absence of orbital effects, the nearly-flat band is in fact slightly chiral and can be assigned a topological character thereby forming a chiral Majorana edge mode.<cit.>

We now turn to the case of combined exchange and orbital effects forν= 1. fig: orbexchb shows that the topological superconductor described above essentially survives the inclusion of orbital effects. Specifically, the strongest effect on the spectrum is for the dispersing in-gap “impurity states”. In contrast, the momentumm_J^*at which the bulk gap closes still matches the value predicted for purely exchange coupling (see Ref. Garnier_SkTopoSC for details).
Regarding the chiral edge mode, interestingly its velocity increases and changes sign. This can be phenomenologically explained by the chiral symmetry interpretation of Ref. Garnier_SkTopoSC. In absence of orbital effects, the smallness of the velocity was attributed to an only weakly broken chiral symmetry given byS = σ_y  τ_y. When including orbital effects (ν≠0), different chiral-symmetry-breaking (CSB) terms appear,  the one proportional to𝕀. This term spatially decays asr^-1, so it should affect chirality stronger than the CSB term in the pure exchange model, which decays asr^-2. 




 §.§ Adding superconducting vortices

The interplay between exchange effects and vortices in a skyrmion proximitized by ans-wave superconductor was studied in Ref. Rex_SK_Vortex_2019 in absence of orbital effects. Here we analyse the full problem by including orbital effects together with the exchange interaction and a superconducting vortex.

Using the notation of the previous section, a superconducting vortex of winding numberαis represented by modifying the pairing Hamiltonian_̋SC = τ_xin eq: BdG_Skyrmion_full to

    _̋ v = Δ_α(r)   e^iαθτ_+ + Δ_α(r)   e^-iαθτ_-

whereτ_±= (τ_x ±i τ_y)/2.[We made sure of the consistency between the GL approach and the BdG one by requiring that the superconducting order parameter is an expectation value of the form ⟨ψψ⟩, which ensures that α has the same meaning in both approaches.] We neglect the spatial variation ofΔ_α(r)and assume that the amplitude of the order parameter is constant, independent ofαand equal to. In this situation, the total angular momentum readsJ_z = L_z + 1/2ħσ_z - α/2ħτ_zand we still denote its eigenvalue bym_J. The momentum is quantized according tom_J ∈ℤ(resp.m_J ∈ℤ + 1/2) ifαodd (resp.αeven). The total radial HamiltonianH_m_J^(3)is given in SM <ref>.

fig: vortex_resultsa, b show the excitation spectrum and the electronic LDoS for a vortex withα= -1without orbital effects, while fig: vortex_resultsc, d show the same with orbital effects ofν= 2( the phase diagram in fig: PhDiag).
fig: vortex_resultsa,b show that the features of the topological superconductor of the exchange model are almost unaffected by the presence of the vortex (without orbital effects), in accord with Ref. Rex_SK_Vortex_2019.

fig: vortex_resultsc,d show that the inclusion of orbital coupling leads to similar effects as in the case of the absence of vortex (previous subsection), namely, (1) it mixes the in-gap impurity states (similar impurity states were found due to purely orbital effects in presence of vortex on a metallic dot in a superconductor<cit.>), and (2) it changes the slope of the topological chiral mode.
For these specific parameters, note that the chiral symmetry breaking due to the vortex (proportional toα/r^2, see eq: TB_Ham_vortex in SM <ref>) adds constructively to the symmetry breaking by the orbital term (proportional toν/r, see eq: TB_Ham_full_supp in SM <ref>), because the change in the slope of the chiral mode is significantly higher when the vortex is included together with the orbital terms.




§ DISCUSSION


Our analysis shows that the exchange-induced topological phase is robust to the inclusion of skyrmion-generated orbital effects as well as superconducting vortices. These effects have to be taken into account since the magnetic field generated by an isolated skyrmion can reach themTrange<cit.>. We have shown using a Ginzburg-Landau approach that proximitizing a skyrmion with a superconductor does not necessarily lead to the formation of vortices. However, even if the superconductor does not develop vortices, the electrons still experience a magnetic vector potential whose effects on the topological superconducting phase were not fully understood in materials without spin-orbit coupling. Our results demonstrate that the inclusion of orbital effects does not invalidate the previously established understanding of the topological superconductor and contributes to making the skyrmion-superconductor hybrid structure a promising platform for the realization of topological superconductivity. Note however that we have neglected the Zeeman effect that would effectively render the exchange interaction anisotropic (see eq: blat). Nevertheless  this renormalization is far below the bare exchange interaction strength, so would not matter even if the exchange strength was reduced by a non-magnetic insulating layer between the skyrmion and the superconductor. 

Because isolated skyrmions with arbitrary winding numbers and helicity  are theoretically more likely to induce topological superconductivity,<cit.> a possible interesting direction would be 
to extend the present magnetostatic calculations for such skyrmions.
More generally, our work calls for a fully self-consistent calculation of both the magnetic and the superconducting order to confirm all the features of the system.




§ ACKNOWLEDGMENTS

The authors acknowledge useful conversation with Freek Massee, Marco Aprili, Stanislas Rohart and Vardan Kaladzhyan. M. G. thanks Marc Gabay for stimulating and insightful discussions.


 Supplementary Material: Magnetic-skyrmion-induced orbital effects in superconductors

Maxime Garnier, Andrej Mesaros and Pascal Simon  
Laboratoire de Physique des Solides, UMR 8502, CNRS, 

Université Paris-Sud, Université Paris-Saclay, 91405 Orsay, France


§ MAGNETOSTATICS OF THE SKYRMION LATTICE

In this appendix we support our hypothesis that the perpendicular inductionB_zcreated by a single skyrmion is proportional to thezcomponent of the magnetization. We do so by focusing on a skyrmion lattice following Ref. QinWang_SKL_2018.
A triangular skyrmion lattice can be approximated using a superposition of three helical spin orders with wavevectorsQ_i=1,  2,  3yielding the so-called triple-Qparametrization<cit.>. In three-dimensional space equiped with the unit vectorsu_x,y,z, consider an infinitely thin magnetic film that lies in thez = 0plane. The three wavevectorsQ_ihave the same norm|Q_i|= Qfor alliand they make a2π/3angle with respect to each other. Our precise choice isQ_1 = Q u_x,Q_2 = Q (-u_x + √(3)u_y)/2andQ_3 = Q (-u_x - √(3)u_y)/2.

For a  skyrmion lattice, the magnetizationm(r, z)withr = (x, y)reads

    m_ lat( r, z )  = m_0  δ( z )  u_z 
        + A ∑_i = 1^3 [ cos( Q_i· r)  u_z + sin( Q_i· r) u_i  ]  δ( z )

where we have defined the unit vectorsu_i = Q_i/Q. We also define the skyrmion radiusas= 2π/(Q√(3))(see caption of Fig. A <ref>). eq: SkL_approx is an approximation in the sense that it is not a proper micromagnetic solution since it is not normalized.

In terms of the thicknesshand the saturation magnetizationM_latof the magnetic film, the parameterAreadsA= M_lat h. In this approach, we define the skyrmion radiusas= 2π/(Q√(3)). An example of such a  skyrmion lattice is presented in Fig. A <ref>.

One way to compute the magnetic field[Strictly speaking B is the magnetic induction but as we focus on the region outside the magnetized medium, the magnetic induction B and the magnetic field H are related by B = μ_0   H so that we call B the magnetic field.]B_lat (r, z )created by the magnetization distribution is to compute the magnetic vector potentialA_lat(r, z )sinceB_lat = ∇ A_lat. Then,A_latis found by solving the Poisson equation∇^2A_lat(r, z ) = -μ_0 J_m(r, z ) in the Coulomb gauge∇ ·A_lat = 0, whereJ_m(r, z ) =  ∇ m_lat (r, z )is the Amperean current density that contains both “bulk” and “surface” contributions. Taking the curl of the solutionA_lat(r, z )of this equation yields 

    B_ lat( r, z )     = μ_0 A/2 e^-Q| z |∑_i = 1^3 [ Q (1 -   sgn( z) )cos( Q_i· r) u_z .
       .+ (Q   sgn( z) -Q + 2δ( z ))sin( Q_i· r)  u_i   ]

This field displays the “single-sided flux” phenomenon or Halbach effect<cit.> meaning that the magnetic field is only present on one side of the plane. The apparent discontinuity in the perpendicular component of the magnetic induction is an artifact of the model and can be regularized by taking into account the finite thickness of the magnetic film so that this effect is indeed physical<cit.>. The perpendicular decay length of the magnetic field is given byQ^-1and can be expressed in terms of the skyrmion radiusasQ^-1 = √(3)/(2 π) ≈0.3   .

With the idea to proximitize the skyrmion lattice by a superconductor, we are only interested in the magnetic field near the plane|Q z |→0. In the case discussed here, the relevant limit isQ z →0^-.
In this limit, thezcomponent of the magnetic field reads

    B^z_ lat( r, 0^- )     =    μ_0 A Q ∑_i = 1^3   cos( Q_i· r)
 close to the planeB^z_latis proportional to thezcomponent of the magnetization. Focusing onz ≪Q^-1,  we argue that this result applies to our case of the isolated skyrmion as long as the thicknessdof the superconductor is much smaller than the skyrmion radius since∝Q^-1. Additionally, we haveA = M_lat h, we find that the amplitude of the magnetic field near the surface of the magnet is given byμ_0 M_lat Q h. Note that this result also holds in the case of a Bloch skyrmion lattice with the only difference that the magnetic field is evenly shared between the two sides of the plane<cit.>. 



§ INFLUENCE OF PARAMETERS ON THE PHASE DIAGRAM

The phase diagram displayed in fig: PhDiag in the main text was computed for ap = 4skyrmion with a cutoffk_s l = 10^-1and we neglected the magnetic energy term proportional toμ_0 M/B_cwhereB_cis the thermodynamic critical field of the superconductor. We now discuss their effects.



 §.§ Effect of the short-distance cutoff

We expect that reducing the cutoff will make the superconducting state with windingα= 0energetically defavorable since is will give more weight to the divergence asr →0. This is indeed what we obtain in Fig. B <ref> by changingk_s l = 10^-1tok_s l = 10^-12.

The normal-superconductor transition line is almost unaffected by this change but theα= 0phase still exists for weak enough skyrmions.



 §.§ Effect of the magnetic energy term

We now reinstate the magnetic energy term we have neglected so far. As it isα-independent, it's only effect will be to move the normal-superconducting transition line depending on the ratioμ_0 M/B_c. However, in our parametrization the parameterνis not independent fromM. In order to estimate the effects of this term, we setμ_0 M/B_c = 1for allν. This is a crude approximation that largely overestimates the strength of this term on a large area of the phase diagram. The phase diagram is shown in Fig. B <ref>
This shows that the inclusion of the magnetic term does not qualitatively affect our results. It further shows that the inhomogeneity of the skyrmion yields a stable superconducting phase even ifμ_0 M > B_c.



 §.§ Influence of the skyrmion size

The results for ap = 10skyrmion are displayed in Fig. B <ref>. These results show that the conclusions of the main text hold qualitatively except that winding number different from -1 have now disappeared from the phase diagram.


§ RADIAL TIGHT-BINDING HAMILTONIANS

In this section we give the expression of the second-quantized discretized radial Hamiltonians used to obtain the results of the main text. In all Hamiltonians we have introducedt = ħ^2/(2m a^2),k_s = k_s awherea ≡1is the radial lattice spacing chosen as the distance unit. We use the notationC^†_j =  (c^†_(j a), c^†_(j a),c_(j a), -c_(j a) )andc_σ(j a)to denote the discretized versions of the spinorΨ(r)and the field operatorψ_σ(r). We have chosenj = 0.5, …, L + 0.5to avoid the singularity at the origin.



 §.§ Orbital effects and exchange interaction

This case corresponds to eq: BdG_Skyrmion_full of the main text. The total angular momentum isJ_z = L_z + 1/2ħσ_zwith eigenvaluem_J(in units ofħ). The discretized Hamiltonian reads

    H_m_J^(1)   = ∑_j -t   C^†_j+1τ_z  C_j  +    +  C^†_j [2t - μ - t/4 j^2(1 - 4 m_J^2 - q^2 ) - t/j^2q m_J σ_z   +t/4 k_s^2 ν^2 (cos( k_s j)/ k_s j + sin( k_s j))^2]τ_z  C_j  
        +  C^†_j [t  k_s ν/j(cos( k_s j)/ k_s j + sin( k_s j))(m_J 𝕀 - q/2σ_z) + J cos( k_s j) σ_z + J sin( k_s j) σ_x + Δ_0 τ_x]  C_j


 §.§ Orbital effects without exchange and vortices

This case corresponds to settingJ = 0in eq: BdG_Skyrmion_full of the main text. or equivalentlyJ = q = 0in eq: TB_Ham_full_supp. The total angular momentum isJ_z = L_zwith eigenvaluem_L(in units ofħ). The discretized Hamiltonian is

    H_m_L^(2) = ∑_j -t   C^†_j+1τ_z  C_j     +    +  C^†_j [2t - μ - t/4 j^2(1 - 4 m_L^2) +t/4 k_s^2 ν^2 (cos( k_s j)/ k_s j + sin( k_s j))^2]τ_z  C_j  
        +  C^†_j [t  k_s ν/j(cos( k_s j)/ k_s j + sin( k_s j))m_L 𝕀 + Δ_0 τ_x]  C_j


 §.§ Orbital effects, exchange interaction and superconducting vortex

Replacing the superconducting pairing term in eq: BdG_Skyrmion_full  by eq: vortex_param (both in the main text), the total angular momentum isJ_z = L_z + 1/2ħσ_z - α/2ħτ_zwith eigenvaluem_J(in units ofħ) whereαis the vortex winding. The discretized Hamiltonian reads

    H_m_J^(3)   = ∑_j = 1^L -t   C^†_j+1τ_z  C_j  +    +  C^†_j [2t - μ - t/4 j^2(1 - 4 m_J^2 - q^2 - α^2)  +t/4 k_s^2 ν^2 (cos( k_s j)/ k_s j + sin( k_s j))^2]τ_z  C_j  
        +  C^†_j [t  k_s ν/j(cos( k_s j)/ k_s j + sin( k_s j))(m_J 𝕀 - q/2σ_z + α/2τ_z) + t/j^2(- q m_J σ_z + α m_J τ_z - α q/4σ_zτ_z)τ_z
    .
       . + J cos( k_s j) σ_z + J sin( k_s j) σ_x + Δ_α(r)  τ_x]  C_j