Introduction

The unification of nontrivial spin texture and superconductivity via advanced interface engineering is a futuristic approach to create and manipulate non-Abelian Majorana bound states (MBS) for their controlled usage in fault-tolerant topological quantum computing1,2,3,4,5. The nanoscale control of magnetism not only relaxes the need for a specific form of Rashba spin-orbit coupling, but also motivates for a magnetic field-free platform for the braiding of the MBS6,7,8,9,10,11,12. Despite numerous successes in the search for the MBS in one-dimensional geometries, the associated limitations such as the intrinsic instabilities of one-dimensional systems, the need for fine tuning of parameters, and the technological obstacles in physical implementation, suggest to look for a two-dimensional platform13,14. The discovery of topological superconductivity in phase-controlled planar Josephson junctions is, therefore, a major step towards the realization of a two-dimensional array of MBS for designing scalable braiding protocols15,16,17,18. The Josephson junction geometry provides additional control to tune the MBS by changing the shape of the junction, strain, and unconventional spin-orbit coupling19,20,21,22,23. A time-reversal invariant topological superconductivity can also be induced by placing the Josephson junction on top of a strong topological insulator24. Previous works on the Josephson junction-based platforms, however, revealed the requirements of a strong intrinsic Rashba spin-orbit coupling and π-phase biasing of the Josephson junction. These constraints pose serious challenges in the detection and manipulation of the MBS under realistic conditions. Chiral magnetism in proximity to an s-wave superconductor generates exotic effects including the appearance of the Majorana modes25,26,27,28,29,30,31,32; however, the location and stability of the Majorana states in these platforms are difficult to anticipate due to their nonlocalized nature.

In our considered geometry, the planar Josephson junction, composed of a two-dimensional electron gas and an s-wave superconductor, is placed on top of a Néel-type skyrmion crystal (SkX) in such a way that the two-dimensional electron gas experiences the spatially varying magnetic field from the bottom SkX and it is also proximitized to the electron pairing from the top superconductors, as described in Fig. 1a. The interplay between the SkX spin texture and the proximity-induced superconductivity leads to topological superconductivity near the middle quasi-one-dimensional channel of the Josephson junction with localized MBS at its two ends. The advantages of using the SkX are: (i) the chiral magnetism generates a robust fictitious gauge field, that can also be visualized as a spin-orbit coupling, and a local Zeeman field which remove the stringent criteria of a strong Rashba-type spin-orbit coupling and, therefore, essentially expands the region of parameter space to realize the MBS, (ii) the existence of the MBS can be further controlled externally by tuning the skyrmion radius, and (iii) the current SkX-based Josephson junction is not required to be phase biased and the MBS can be found robustly at a phase difference φ = 0, between the two superconducting regions, unlike the usual planar Josephson junctions that are required to be phase biased with φ = π, to minimize the critical magnetic field for the topological transition and to maximize the chemical-potential range within which the MBS appear15. Using the zero-energy feature of the quasiparticle states with a topological energy gap, sharp localization of these states, charge-neutrality condition, two order parameters, viz. Majorana polarization and curvature of the density of states, we confirm the existence of the MBS in our set up. The tunable phase difference and the skyrmion radius together provide broad, flexible control of the MBS which is indispensable to achieve the long-sought-after goal of the braiding.

Fig. 1: Device geometry and a skyrmion crystal.
figure 1

a Planar Josephson junction on top of a skyrmion crystal. The two-dimensional electron gas (2DEG) exhibits both proximity-induced superconductivity from the top superconductor (SC) layers and spatially varying magnetism from the bottom skyrmion crystal that is created in the ferromagnet due to the competition between exchange interactions in the ferromagnet (FM) and the heavy metal or heavy insulator (HM/HI), with a field or anisotropy. The zero-energy Majorana bound states (shown as yellow bubbles) are localized at the two ends of the quasi-one-dimensional metallic channel. b The skyrmion crystal spin texture, spontaneously generated in a Monte Carlo simulation using a 100 × 100 × 6 lattice with ferromagnetic exchange interaction strength J = 1, Dzyaloshinskii-Moriya interaction strength D = 0.3J, magnetic field Hz = 0.1J, spin amplitude S = 1, and easy-plane anisotropy A = 0.01J. The colorbar in b denotes the z component of the magnetization mz.

Results

Theoretical set up

For the generation of the SkX, we consider a heterointerface of a thin-layer ferromagnet and a heavy compound (metal or insulator). The advantage of the heavy compound is that it helps to generate a large Dzyaloshinskii-Moriya interaction (DMI) at the interface between the ferromagnet and the heavy compound. The cooperation between the DMI and the ferromagnetic exchange interaction of the ferromagnet produces a triangular SkX, in the presence of a magnetic field or an anisotropy. Our Monte Carlo simulations reveal that columns of skyrmions, arranged in a triangular array, appear spontaneously within a six-layer ferromagnet, although the DMI exists predominantly at the interface between the ferromagnet and the heavy compound, as shown in Fig. 1b. We perform simulated annealing using the Metropolis energy-minimization algorithm, formulated with the following Hamiltonian

$${\mathcal{H}}\ = \ -J\mathop{\sum}\limits_{\langle ij\rangle }{{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j}-D\mathop{\sum}\limits_{\langle ab\rangle }(\hat{z}\times {\hat{r}}_{ab})\cdot ({{\bf{S}}}_{a}\times {{\bf{S}}}_{b})\\ \, -{H}_{z}\mathop{\sum}\limits_{i}{S}_{zi}-A\mathop{\sum}\limits_{a}| {S}_{za}{| }^{2},$$
(1)

where J is the nearest-neighbor ferromagnetic exchange interaction strength in the ferromagnet, D is the DMI strength at the bottom ferromagnet layer that interfaces with the heavy compound, Hz is the perpendicular magnetic field, A is the easy-plane magnetic anisotropy at the bottom ferromagnet layer, i,j are the site indices in the entire ferromagnet, and a,b are the two-dimensional site indices at the bottom ferromagnet layer. The DMI, present dominantly at the interface between the ferromagnet and the heavy compound, generates a Néel-type SkX33,34,35. Besides the engineered interfaces, the SkX naturally appears in a wide variety of materials36,37,38,39 that can also be utilized in the proposed device geometry, instead of the combination of the ferromagnet and the heavy compound. Also, the SkX can be artificially created without the need for any external magnetic field by nanopatternization40.

The spin texture Bi on the top layer of the ferromagnet, obtained from the Monte Carlo simulations, is used to obtain the low-energy spectrum of the planar Josephson junction by solving self-consistently the Bogoliubov-de Gennes equations. Since the SkX lies underneath the two-dimensional electron gas without a finite separation between them, it is reasonable to assume that the deviation in the magnetic fringing field in the two-dimensional electron gas from the original SkX texture is negligible. The proximity-induced superconductivity in the two-dimensional electron gas, which is subject to the SkX spin texture Bi, is described by the Hamiltonian

$${{\mathcal{H}}}_{{\rm{BdG}}}\ = \ -t\mathop{\sum}\limits_{\langle ij\rangle ,\sigma }({c}_{i\sigma }^{\dagger }{c}_{j\sigma }+H.c.)+\mathop{\sum}\limits_{i,\sigma }(4t-\mu ){c}_{i\sigma }^{\dagger }{c}_{i\sigma }\\ \, -\frac{1}{2}g{\mu }_{B}\mathop{\sum}\limits_{i,\sigma }{({{\bf{B}}}_{i}\cdot {\boldsymbol{\sigma }})}_{\sigma {\sigma }^{\prime}}{c}_{i\sigma }^{\dagger }{c}_{i{\sigma }^{\prime}}+\mathop{\sum}\limits_{i}({{{\Delta }}}_{i}{c}_{i\uparrow }^{\dagger }{c}_{i\downarrow }^{\dagger }+H.c.),$$
(2)

where t = 2/(2m*a2) is the hopping energy, m* is the effective mass of electrons, a is the unit spacing of the lattice grid, μ is the chemical potential, and Δi is the induced local s-wave pairing amplitude on the two sides of the Josephson junction that are attached to the top Al layer. The pairing amplitude Δi = −Uicici〉 is calculated self-consistently using the onsite attractive interaction strength Ui of the induced superconducting states in the two-dimensional electron gas. Ui = U in the two-dimensional electron gas below the Al superconductors and zero in the middle metallic channel. The value U = 2 meV is determined by setting Δi = 0.2 meV, the estimated proximity-induced gap magnitude for a two-dimensional electron gas with an SC interface41,42, without any spin texture. The g factor and the effective mass are set to g = 50 and m* = 0.017m0 for InSb43,44. The lattice grid spacing used is a = 10 nm (a is the lattice grid spacing used to discretize the kinetic energy term \(\frac{{p}^{2}}{2m}\) within finite-difference approximation. It is immaterial as long as the lattice lengths remain fixed) with which the hopping energy becomes t = 22.44 meV. The amplitude of the spin texture Bi is set to B0 = 0.3 T, compatible with the saturation magnetization Ms = 1.7 × 106 A/m for CoFe10,11. HM/ferromagnet interfaces with Pt, Pd, Ag, Ir, and Au as the HM have been developed together with Co, Fe, and their alloys (see e.g., ref. 45). We present results for a planar Josephson junction with length Ly = 2 μm, transverse length of the SC leads Lx = 200 nm, and width of the quasi-one-dimensional metallic channel W = 50 nm.

Emergence of the MBS

The low-energy spectrum, shown in Fig. 2a, reveals that there exist multiple ranges of the chemical potential within which the zero-energy MBS appear. To determine the Majorana character of the quasiparticle states, we compute the Majorana polarization, defined as46,47

$${{\mathcal{P}}}_{{{\mathcal{M}}}_{,n}}=2\mathop{\sum}\limits_{i}{u}_{i\downarrow }^{n}{v}_{i\downarrow }^{n* }-{u}_{i\uparrow }^{n}{v}_{i\uparrow }^{n* },$$
(3)

where \({u}_{i\uparrow }^{n}\) and \({v}_{i\uparrow }^{n}\) are the Bogoliubov-de Gennes quasiparticle and quasihole amplitudes, respectively, corresponding to the nth eigenstate, spin , and site i. As evident from Fig. 2a, \(| {{\mathcal{P}}}_{{{\mathcal{M}}}_{,n}}| \approx 1\) indicates the occurrence of a pair of robust MBS with a finite topological energy gap. The Majorana polarization \(| {{\mathcal{P}}}_{{\mathcal{M}},1}|\) of the first positive eigenstate, plotted with μ in Fig. 2b, acquires finite values within the range of μ, in which the MBS emerge. The delta-function-like peaks in \(| {{\mathcal{P}}}_{{\mathcal{M}},1}|\) are the signatures of the Majorana oscillations, which is also clearly seen in the low-energy spectrum in Fig. 2a, originating due to the overlap of the MBS wave functions at the two ends of the finite-length quasi-one-dimensional channel. The Majorana oscillations in \(| {{\mathcal{P}}}_{{\mathcal{M}},1}|\) have also been confirmed from the calculations of a one-dimensional wire (for results in the wire geometry, see Supplementary Note 4). The Majorana polarization, with a modification in the expression used in Eq. (3), was proposed to be probed in this planar Josephson junction geometry using the spin-selective Andreev reflection technique48. To further characterize the evolution of the topological superconductivity with changing a parameter, such as μ, we look at the curvature of the density of states at zero energy \(\frac{{\partial }^{2}D}{\partial {E}^{2}}\), where D(E) is defined as49

$$D(E)=\mathop{\sum}\limits_{i,n,\sigma }(| {u}_{i\sigma }^{n}{| }^{2}+| {v}_{i\sigma }^{n}{| }^{2})\delta (E-{E}_{n}),$$
(4)

and δ(E − En) is modeled using a Gaussian with broadening 0.001 meV (t). The second derivative is computed using the second-order finite-difference method. These two quantities, \(| {{\mathcal{P}}}_{{\mathcal{M}},1}|\) and \(\frac{{\partial }^{2}D}{\partial {E}^{2}}\), may provide additional insight in the experimental detection of the MBS, besides the conventional zero-bias conductance peak50, which often leads to ambiguity due to other possible zero-bias states in a superconductor51.

Fig. 2: Emergence of the Majorana bound states with changing chemical potential.
figure 2

a The quasiparticle spectrum of the planar Josephson junction at phase difference φ = 0 with varying chemical potential (μ), showing the emergence of the zero-energy Majorana bound states. The colorbar represents the Majorana polarization \(| {{\mathcal{P}}}_{{\mathcal{M}},n}|\) that displays the Majorana character of the quasiparticle states. The skyrmion crystal with a skyrmion diameter Dsk = 10a was obtained using a magnetic field Hz = 0.95J and a Dzyaloshinskii-Moriya interaction strength D = 1.6J in the Monte Carlo calculations. b, c The Majorana polarization \(| {{\mathcal{P}}}_{{\mathcal{M}},1}|\) of the first positive eigenstate and the curvature of the density of states at zero energy \(\frac{{\partial }^{2}D}{\partial {E}^{2}}\), with varying μ, showing the μ range within which the Majorana bound states appear. The delta-function-like peaks are associated with the oscillations of the MBS with changing μ. d, e The profiles of the local density of states (DOS) and the charge DOS at μ = 0.5 meV. The green rectangle in d indicates the quasi-one-dimensional metallic channel of the planar Josephson junction at the ends of which the Majorana bound states appear. f, g The profiles of the local DOS and the charge DOS at μ = −0.5 meV, in the nontopological regime.

As shown in Fig. 2c, \(\frac{{\partial }^{2}D}{\partial {E}^{2}}\) takes finite values in the same ranges of μ as that of the Majorana polarization \(| {{\mathcal{P}}}_{{\mathcal{M}},1}|\). The Majorana oscillations, in the form of delta-function-like peaks, is also noticeable in \(\frac{{\partial }^{2}D}{\partial {E}^{2}}\), albeit with changes in the sign. To visualize the location of the zero-energy MBS, we show, in Fig. 2d, the profile of the local density of states \({\rho }_{{}_{{\rm{LDOS}}}}^{i}\ =\ {\sum }_{\sigma }(| {u}_{i\sigma }{| }^{2}+| {v}_{i\sigma }{| }^{2})\), corresponding to the lowest positive-energy eigenstate at μ = 0.5 meV where the Josephson junction is in the topological superconducting regime. The sharp peaks in \({\rho }_{{}_{{\rm{LDOS}}}}^{i}\) indicate that the MBS are localized predominantly near the two ends of the quasi-one-dimensional channel. Figure 2e shows the profile of the charge density of states \({\rho }_{{}_{{\rm{CDOS}}}}^{i}\ =\ {\sum }_{\sigma }(| {u}_{i\sigma }{| }^{2}-| {v}_{i\sigma }{| }^{2})\) corresponding to the lowest positive-energy eigenstate at μ = 0.5 meV. The profiles of \({\rho }_{{}_{{\rm{LDOS}}}}^{i}\) and \({\rho }_{{}_{{\rm{CDOS}}}}^{i}\), at μ = −0.5 meV where the Josephson junction is in the topologically trivial superconducting regime, are shown in Fig. 2f, g, respectively. In this case, both the quasiparticle state and the charge density are distributed near the middle of the quasi-one-dimensional channel. Interestingly, a comparison of Fig. 2e, g, implies an order-of-magnitude suppression in \({\rho }_{{}_{{\rm{CDOS}}}}^{i}\), which is reminiscent of the local charge-neutrality signature of the MBS and is another confirmation of the Majorana character of this state. The above results establish that the spin-orbit coupling, generated by the SkX, alone can lead to the emergence of the MBS in the planar Josephson junction devices.

Skyrmion control

The skyrmion size in a SkX is tunable, with remarkable precession, by an external magnetic field, magnetic anisotropy and advanced symmetry protocol at heterointerfaces52,53,54. In our Monte Carlo simulations, the skyrmion size was varied by tuning the magnetic field and the DMI, as shown in Fig. 3a-c. The Bogoliubov-de Gennes quasiparticle spectra at different skyrmion sizes, shown in Fig. 3d-f, imply that the presence of the zero-energy MBS at a given chemical potential can be turned ON or OFF by only changing the skyrmion properties. The minimum diameter of the skyrmions, realized in our Monte Carlo simulations, is 10 lattice spacings for which the MBS appear in the discussed planar Josephson junction setting. With increasing the skyrmion size, we find that the range of the chemical potential within which the MBS appear decreases effectively, however, the oscillation amplitude of the MBS is suppressed gradually. Therefore, the skyrmions offer a unique ability to manipulate the localization length of the MBS in the planar Josephson junctions. The strongly localized MBS in the skyrmion-tuned planar Josephson junctions can, therefore, have advantageous over other platforms for MBS realization in fault-tolerant topological quantum computing.

Fig. 3: Skyrmion control of the Majorana bound states.
figure 3

ac The skyrmion crystals of different skyrmion diameters a Dsk = 12a, b Dsk = 14a, c Dsk = 16a, obtained in the Monte Carlo calculations using a a magnetic field Hz = 0.8J, a Dzyaloshinskii-Moriya interaction strength D = 1.4J, b Hz = 0.45J, D = J, and c Hz = 0.23J, D = 0.6J. df The corresponding quasiparticle spectra of the planar Josephson junction at the phase difference φ = 0 with varying chemical potential, obtained by solving the Bogoliubov-de Gennes equations with the above skyrmion crystal spin configurations. The colorbar in df represents the Majorana polarization \(| {{\mathcal{P}}}_{{\mathcal{M}},n}|\) of the quasiparticle states.

The broken inversion symmetry at the interface between the two-dimensional electron gas and the superconductor, often leads to a sizable intrinsic Rashba spin-orbit coupling, which is usually considered as the primary mechanism for modifying the pairing symmetry of the induced superconductivity55,56, leading to the desired topological superconductivity. We find that the MBS remain robust in the presence of the intrinsic Rashba spin-orbit coupling (for details on the effect of Rashba spin-orbit coupling, see Supplementary Note 2). By placing the SkX texture only below the middle quasi-1D channel, we didn’t find robust MBS formation (for details, see Supplementary Note 3) and hence it suggests the significance of placing the SkX texture below the entire Josephson junction.

Phase control

Another important control parameter, that sets the planar Josephson junctions apart from other platforms hosting the MBS, is the phase difference φ between the two superconducting regions of a Josephson junction. The theoretical prediction15 and the subsequent experimental discoveries16,17 suggest that the Josephson junction needs to be biased by a phase difference φ = π to minimize the critical Zeeman field, required for inducing the topological superconductivity. Remarkably, in the current Josephson junction set up with the SkX, the topological superconductivity is induced at φ = 0, as depicted by the quasiparticle spectrum with varying φ in Fig. 4a. With increasing φ, the MBS move gradually from zero to higher energies, indicating an enhancement in the localization length of the MBS. The MBS appear again at zero energy above φ ≈ 3π/2.

Fig. 4: Phase control of the Majorana bound states.
figure 4

a Quasiparticle spectrum of the planar Josephson junction at a chemical potential μ = 0.5 meV with changing phase difference φ between the two superconducting regions. The skyrmion diameter of the skyrmion crystal was Dsk = 10a. b Colormap of the Majorana polarization of the first positive-energy state in the plane spanned by the phase difference φ and the chemical potential μ. The Majorana bound states appear in the parameter regimes with \(| {{\mathcal{P}}}_{{\mathcal{M}},1}| \approx\)1, φ = 0 being the most favorable scenario in this set up. The cyan and magenta dots in b represent the topological invariant, respectively, Q = − 1 (topologically nontrivial) and Q = +1 (topologically trivial).

This dephasing effect of the MBS can be understood from the Majorana oscillations—as we find that the oscillation increases with increasing φ in the range 0 < φ ≤ π (for details, see Supplementary Note 1). The finite length of the quasi-one-dimensional metallic channel gives rise to the oscillations of the zero-energy MBS with varying chemical potential μ. Furthermore, the finite width of the metallic channel provides extra room for delocalization of the MBS at the two ends, contributing additively to the Majorana oscillation. In the previous works on this geometry, a magnetic field is applied in the plane of the planar Josephson junction, which locks the phases of the MBS to φ = π. In our set up, there is no magnetic field applied to a particular direction, but a chiral magnetism exists throughout the junction and, therefore, a π phase biasing is not required in this case. The middle metallic region of the Josephson junction can be perceived as a quasi-one-dimensional void region surrounded by the superconducting two-dimensional electron gas. Additional phase difference between the two superconducting sides, therefore, only causes disruption to the induced topological superconductivity. This phenomenon generically takes place at several values of the chemical potential, as shown in Fig. 4b, where we plot the Majorana polarization \(| {{\mathcal{P}}}_{{\mathcal{M}},1}|\) of the lowest positive-energy eigenstate in the parameter space spanned by the phase difference φ and the chemical potential μ. For the chosen range of μ values, the Majorana polarization decreases substantially within the range π/2 φ 3π/2. The Majorana oscillation in \(| {{\mathcal{P}}}_{{\mathcal{M}},1}|\), however, survives up to φ ≈ π. At φ = π, \(| {{\mathcal{P}}}_{{\mathcal{M}},1}|\) vanishes completely, indicating the disappearance of the MBS. Therefore, φ = 0 is the most favorable condition to realize the MBS in our Josephson junction set up and the phase difference can be further tuned to control the presence of the MBS. To check the consistency of our assignment of MBS, we also compute the Z2 topological index Q at several points of the above phase diagram using an effective one-dimensional Hamiltonian of the planar Josephson junction with the SkX. In Fig. 4b, we show the topological invariant which confirms the phase diagram.

Conclusion

The skyrmions bring outstanding control functionalities to the planar Josephson junctions for the creation and manipulation of the zero-energy MBS and their localization properties. The SkXs, being realized in an abundance of magnetic materials and also artificially created in patterned magnetic materials, offer a feasible approach for advanced manipulation of the zero-energy MBS. The proposed planar Josephson junction, combined with a SkX, has the major advantages that there is no need for a strong intrinsic Rashba-type spin-orbit coupling and phase-biasing constraint for the realization of the zero-energy MBS. The enhanced tunability of the MBS in the proposed two-dimensional platform opens up opportunities for designing feasible MBS braiding protocols for the fault-tolerant topological quantum computing, and investigating Majorana spectroscopy using the multi-terminal superconducting quantum interference devices.

Methods

Monte Carlo simulations

The SkX spin configurations were obtained using a Lx × Ly × Lz lattice with periodic boundary conditions along the x and y directions and open boundary conditions along the z direction. A bias-free sampling method, that provides a full and uniform coverage of the phase space spanned by the spin angles, was used for generating the completely random spin configurations. The calculation was started at a high temperature T = 10J with a random spin configuration and the temperature was lowered slowly down to a low value T = 0.001J in 2000 steps. At each temperature step, 1010 Monte Carlo spin updates were performed. In each spin update step, a new spin direction was chosen randomly within a small cone spanned around the initial spin direction. The new spin configuration was accepted or rejected according to the Metropolis energy-minimization algorithm by comparing the total energies of the previous and the new trial spin configurations, calculated using the Hamiltonian Eq. (1).

Self-consistent Bogoliubov-de Gennes formalism

The Hamiltonian Eq. (2), which is quadratic in the fermionic operators \({\hat{c}}_{i\sigma }\), can be solved by exact diagonalization via a unitary transformation \({\hat{c}}_{i\sigma }={\sum }_{n}{u}_{i\sigma }^{n}{\hat{\gamma }}_{n}+{v}_{i\sigma }^{n* }{\hat{\gamma }}_{n}^{\dagger }\), where \({\hat{\gamma }}_{n}^{\dagger }\) (\({\hat{\gamma }}_{n}\)) is a fermionic creation (annihilation) operator of the quasiparticle/quasiphole state in the nth energy eigenstate. The quasiparticle amplitudes \({u}_{i\sigma }^{n}\) and the quasihole amplitudes \({u}_{i\sigma }^{n}\) are determined by solving the Bogoliubov-de Gennes equations: \({\sum }_{j}{{\mathcal{H}}}_{ij}{\psi }_{j}^{n}={\epsilon }_{n}{\psi }_{n}^{i}\), where \({\psi }_{i}^{n}={[{u}_{i\uparrow }^{n},{u}_{i\downarrow }^{n},{v}_{i\uparrow }^{n},{v}_{i\downarrow }^{n}]}^{T}\) is the basis wave function and ϵn is the nth energy eigenvalue. The Hamiltonian \({{\mathcal{H}}}_{ij}\) is expressed in the following matrix form

$${{\mathcal{H}}}_{ij}=\left(\begin{array}{llll}{{\mathcal{H}}}_{\uparrow \uparrow }&{{\mathcal{H}}}_{\uparrow \downarrow }&0&{{{\Delta }}}_{i}\\ {{\mathcal{H}}}_{\downarrow \uparrow }&{{\mathcal{H}}}_{\downarrow \downarrow }&-{{{\Delta }}}_{i}&0\\ 0&-{{{\Delta }}}_{i}^{* }&-{{\mathcal{H}}}_{\uparrow \uparrow }^{* }&-{{\mathcal{H}}}_{\uparrow \downarrow }^{* }\\ {{{\Delta }}}_{i}^{* }&0&-{{\mathcal{H}}}_{\downarrow \uparrow }^{* }&-{{\mathcal{H}}}_{\downarrow \downarrow }^{* }\end{array}\right),$$
(5)

where \({{\mathcal{H}}}_{\uparrow \uparrow ,\downarrow \downarrow }=-t(1-{\delta }_{ij})+(4t-\mu ){\delta }_{ij}-1/2\sigma g{\mu }_{B}{B}_{z}\), where σ = ± for , , and \({{\mathcal{H}}}_{\uparrow \downarrow }=-1/2g{\mu }_{B}({B}_{x}+i{B}_{y})\). The s-wave pairing amplitude Δi = −Uicici〉 is computed using \({{{\Delta }}}_{i}=-{U}_{i}{\sum }_{n}\left[\right.{u}_{i\uparrow }^{n}{v}_{i\downarrow }^{n* }(1-f({\epsilon }_{n}))+{u}_{i\downarrow }^{n}{v}_{i\uparrow }^{n* }f({\epsilon }_{n})\left]\right.\), where Ui = U inside the two superconducting regions and zero in the middle quasi-one-dimensional metallic channel, \(f({\epsilon }_{n})=1/(1+{e}^{{\epsilon }_{n}/{k}_{B}T})\) is the Fermi-Dirac distribution function at temperature T. The self-consistency iterations were performed until the pairing amplitudes Δi were converged at every lattice sites.

Calculation of topological invariant

A triangular SkX can be described by the following spin structure

$${{\bf{S}}}_{i}=S(\sin {\theta }_{i}\cos {\phi }_{i},\,\sin {\theta }_{i}\sin {\phi }_{i},\,\cos {\theta }_{i}),$$
(6)

where S is the spin amplitude and the spin angles θi and ϕi are defined as

$${\theta }_{i}=\pi \,\min \left(\frac{| {\bf{r}}-{{\bf{R}}}_{i}| }{R},1\right),{\phi }_{i}=\arctan \left(\frac{y-{y}_{i}}{x-{x}_{i}}\right),$$
(7)

Ri = (xi, yi) is the skyrmion center near position r = (x, y) and R is the skyrmion radius. To obtain an effective Hamiltonian for the planar JJ, we perform the following unitary transformation that rotates the local spin Si with respect to the SkX plane normal, the z direction

$$\left(\begin{array}{l}{c}_{i\uparrow }\\ {c}_{i\downarrow }\end{array}\right)={\hat{{{\Gamma }}}}_{i}\left(\begin{array}{l}{d}_{i\uparrow }\\ {d}_{i\downarrow }\end{array}\right)$$
(8)

where

$${\hat{{{\Gamma }}}}_{i}=\exp \left(-i\frac{{\theta }_{i}}{2}\frac{\hat{z}\times {{\bf{S}}}_{i}}{| \hat{z}\times {{\bf{S}}}_{i}| }\cdot \sigma \right).$$
(9)

The Hamiltonian (2) in the main text, in the transformed basis, becomes

$${\mathcal{H}}\ = \, -\mathop{\sum}\limits_{\langle ij\rangle ,\sigma ,{\sigma }^{\prime}}{t}_{\sigma {\sigma }^{\prime}}^{\prime}{d}_{i\sigma }^{\dagger }{d}_{j{\sigma }^{\prime}}+\mathop{\sum}\limits_{i,\alpha }\left(4t-\mu -\frac{1}{2}g{\mu }_{B}{B}_{0}{\sigma }_{\sigma \sigma }^{z}\right){d}_{i\sigma }^{\dagger }{d}_{i\sigma }\\ \, +\mathop{\sum}\limits_{i}\left({{{\Delta }}}_{i}{d}_{i\uparrow }^{\dagger }{d}_{i\downarrow }^{\dagger }+{\rm{H.c.}}\right),$$
(10)

where the new hopping integral is given by

$${t}^{\prime}=\left(\begin{array}{ll}-t&-{\alpha }^{* }\\ \alpha &-{t}^{* }\end{array}\right),$$
(11)

α being the strength of the generated SOC.

To obtain an effective one-dimensional Hamiltonian, we consider an infinitely long junction (i.e., Ly → ) and perform a partial Fourier transform \({d}_{i,{k}_{y},\sigma }={\sum }_{j}{e}^{i{k}_{y}y}{d}_{i,j,\sigma }\). The resulting Hamiltonian \({\mathcal{H}}(x,{k}_{y})\) is then used to compute the topological invariant Q, the Z2 topological index associated with the broken chiral symmetry, given by57

$$Q={\rm{sgn}}\left[\frac{{\rm{Pf}}\{{\mathcal{H}}(x,{k}_{y}=\pi )\}}{{\rm{Pf}}\{{\mathcal{H}}(x,{k}_{y}=0)\}}\right]$$
(12)

where ‘Pf’ denotes the Pfaffian. A topologically nontrivial phase is determined by Q = −1, while Q = 1 represents the trivial phase.