Skyrmion Control of Majorana States in Planar Josephson Junctions

Planar Josephson junctions provide a versatile platform, alternative to the nanowire-based geometry, for the generation of the Majorana bound states, due to the additional phase tunability of the topological superconductivity. The proximity induction of chiral magnetism and superconductivity in a two-dimensional electron gas showed remarkable promises to manipulate topological superconductivity. Here, we consider a Josephson junction involving a skyrmion crystal and show that the chiral magnetism of the skyrmions can create and control the Majorana bound states without the requirement of an intrinsic Rashba spin-orbit coupling. Interestingly, the Majorana bound states in our geometry are realized robustly at zero phase difference at the junction. The skyrmion radius, being externally tunable by a magnetic field or a magnetic anisotropy, brings a unique control feature for the Majorana bound states.

The unification of non-trivial 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 computing [1][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 MBS [6][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 platform [13,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 protocols [15][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 coupling [19][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 insulator [24]. 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 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 spatiallyvarying 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) usual planar Josephson junctions are required to be phase biased with a phase difference ϕ = π, between the two superconducting regions, to minimize the critical magnetic field for the topological transition and to maximize the chemical-potential range within which the MBS appear [15]; the current SkX-based Josephson junction is not required to be phase biased and the MBS can be found robustly at ϕ = 0. 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. 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 inter-face 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 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, H z 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 SkX [1,2,33]. Besides the engineered interfaces, the SkX naturally appears in a wide variety of materials [36][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 artificiallycreated without the need for any external magnetic field by nanopatternization [40].
The spin texture B i 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 Bogoliubovde 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 twodimensional 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 B i , is described by the Hamiltonian where t = 2 /(2m * a 2 ) 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 = −U i c i↑ c i↓ is calculated self-consistently using the onsite attractive interaction strength U i of the induced superconducting states in the two-dimensional electron gas. U i = 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 interface [41,42], without any spin texture. The g factor and the effective mass are set to g = 50 and m * = 0.017m 0 for InSb [43,44]. The lattice grid spacing used is a = 10 nm [45] with which the hopping energy becomes t = 22.44 meV. The amplitude of the spin texture B i is set to B 0 = 0.3 T, compatible with the saturation magnetization M s = 1.7 × 10 6 A/m for CoFe [10,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. 46). We present results for a planar Josephson junction with length L y = 2 µm, transverse length of the SC leads L x = 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 larization, defined as [5,6] where u n i↑ and v n i↑ are the Bogoliubov-de Gennes quasiparticle and quasihole amplitudes, respectively, corresponding to the n th eigenstate, spin ↑, and site i. 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 spinselective Andreev reflection technique [49]. 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 and δ(E −E n ) 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, |P M,1 | and ∂ 2 D ∂E 2 , may provide additional insight in the experimental detection of the MBS, besides the conventional zero-bias conductance peak [8] which often leads to ambiguity due to other possible zero-bias states in a superconductor [9].
As shown in Fig. 2c, ∂ 2 D ∂E 2 takes finite values in the same ranges of µ as that of the Majorana polarization |P M,1 |. The Majorana oscillations, in the form of delta-function-like peaks, is also noticeable in ∂ 2 D ∂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 ρ i LDOS = σ (|u iσ | 2 + |v iσ | 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 ρ i LDOS 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 corresponding to the lowest positive-energy eigenstate at µ = 0.5 meV. The profiles of ρ i LDOS and ρ i CDOS , at µ = −0.5 meV where the Josephson junction is in the topologically-trivial superconducting regime, are shown in Fig. 2f and Fig. 2g, respectively. In this case, both the quasiparticle state and the charge density are distributed near the middle of the quasione-dimensional channel. Interestingly, a comparison of Fig. 2e and Fig. 2g, implies an order-of-magnitude suppression in ρ i CDOS 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 heterointerfaces [53][54][55]. In our Monte Carlo simulations, the skyrmion size was varied by tuning the magnetic field and the DMI, as shown in Figs. 3a-c. The Bogoliubov-de Gennes quasiparticle spectra at different skyrmion sizes, shown in Figs. 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 stronglylocalized MBS in the skyrmion-tuned planar Josephson junctions can, therefore, have advantageous over other platforms for MBS realization in fault-tolerant topological quantum computing.
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 superconductivity [3,4], 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. platforms hosting the MBS, is the phase difference ϕ between the two superconducting regions of a Josephson junction. The theoretical prediction [15] and the subsequent experimental discoveries [16,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. 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 |P 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 |P M,1 |, however, survives up to ϕ ≈ π. At ϕ = π, |P 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 Z 2 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, 10 10 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 (1).

Self-consistent Bogoliubov-de Gennes formalism.
The Hamiltonian (2), which is quadratic in the fermionic operatorsĉiσ, can be solved by exact diagonalization via a unitary transformationĉiσ = n u n iσγn + v n * iσγ † n , wherê γ † n (γn) is a fermionic creation (annihilation) operator of the quasiparticle/quasiphole state in the n th energy eigenstate.
The quasiparticle amplitudes u n iσ and the quasihole amplitudes u n iσ are determined by solving the Bogoliubov-de Gennes equations: j Hijψ n j = nψ i n , where ψ n i = [u n i↑ , u n i↓ , v n i↑ , v n i↓ ] T is the basis wave function and n is the n th energy eigenvalue. The Hamiltonian Hij is expressed in the following matrix form where where Ui = U inside the two superconducting regions and zero in the middle quasi-one-dimensional metallic channel, f ( n) = 1/(1 + e n/kB 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 Si = S(sin θi cos φi, sin θi sin φi, cos θi), where S is the spin amplitude and the spin angles θi and φi are defined as 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 The Hamiltonian (2) in the main text, in the transformed basis, becomes where the new hopping integral is given by α 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,ky ,σ = j e iky y di,j,σ. The resulting Hamiltonian H(x, ky) is then used to compute the topological invariant Q, the Z2 topological index associated with the broken chiral symmetry, given by [58] Q = sgn Pf{H(x, ky = π)} Pf{H(x, ky = 0)} where 'Pf' denotes the Pfaffian. A topologically nontrivial phase is determined by Q = −1, while Q = 1 represents the trivial phase.

Supplementary Note 1: Majorana oscillation in the Josephson junctions
The oscillation of the Majorana bound states (MBS) is a wellknown phenomenon that appears due to the finite size effect. In a one-dimensional geometry, e.g. a Rashba spin-orbit coupled nanowire with a Zeeman exchange coupling, the probability density of the MBS, which are localized dominantly at the two ends, decay exponentially with distance towards the middle of the nanowire. The overlap of these two MBS wave functions give rise to a finite splitting in energy. This split energy gap between the MBS pair oscillates with varying a parameter, such as the chemical potential or the Zeeman energy. In our considered planar Josephson junction geometry, the finite length of the quasi-one-dimensional metallic channel in the middle, therefore, naturally leads to the oscillations of the zero-energy MBS with varying chemical potential µ. Furthermore, the finite width of the quasi-one-dimensional channel provides extra room for delocalization of the MBS at the two ends, contributing additively to the Majorana oscillation. The oscillation amplitude, however, decreases with increasing the length of the Josephson junction. In Fig. S1, we show the quasiparticle spectra with varying µ at different values of the phase difference ϕ between the two superconducting regions of the Josephson junction. Evidently, the oscillation amplitude increases with increasing ϕ from 0 to π. At ϕ = π, the oscillation becomes large and the MBS completely vanish. The oscillation can also be visible in the Majorana polarization |PM,1| of the lowest positive eigenstate. In Fig. S2, we show the variation in |PM,1| for the lowest positive eigenstate with µ. At ϕ = 0, the oscillation is small in amplitude and appears on the top of a significant value of the Majorana polarization (|PM,1|≈1). With increasing ϕ from 0 to π, the value of |PM,1| decreases and the oscillation amplitude increases simultaneously. These results explain why the MBS disappear with increasing the phase difference ϕ from 0 to π, as discussed in the main text.

Supplementary Note 2: Effect of Rashba spin-orbit coupling
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 superconductivity in a one-dimensional [3] or two-dimensional geometry [4], leading to the desired topological superconductivity. To explore the mutual influence of the Rashba spin-orbit coupling and the Skyrmion crystal (SkX)-generated spin-orbit coupling on the emergence of the zero-energy MBS, we consider the Hamiltonian HRSOC=−iα/(2a) ij ,σ,σ (σ ×rij) z σσ c † iσ c jσ for the Rashba spin-orbit coupling, where α is the strength of the Rashba spin-orbit coupling, and obtain the quasiparticle spectrum based on the total Hamiltonian H BdG + HRSOC. In range of µ start to become delocalized, possibly because a very large Rashba spin-orbit coupling interferes destructively with the synthetic spin-orbit coupling. But remarkably, with larger values of α, new MBS start to appear at larger values of the chemical potential, e.g. within the range 2.2 meV µ 2.6 meV. Although the mutual effects of the intrinsic Rashba spin-orbit coupling and the SkX-generated spin-orbit coupling are rather complex, within some ranges of µ, they cooperate to generate robust MBS with a larger topological energy gap.
Supplementary Note 3: Skyrmion crystal below the quasi-1D channel We studied the possibility of obtaining the MBS with the SkX, lying only below the non-superconducting region of the planar Josephson junction. As shown by the spectra in is no signature of the formation of robust MBS with strongly localized nature. The overall influence of the skyrmion spin texture on the quasiparticle spectrum is lesser than the case when the SkX is placed below the entire Josephson junction. We would further like to stress on the fact that the planar Josephson junction operates in a slightly different mechanism than the 1D nanowire setting. In the 1D nanowire set up, the topological superconductivity is realized within the nanowire, while in the planar Josephson junction, the MBS appear in the middle region, a region that is void of superconducting pairing. Therefore, the influence of the skyrmion texture on the superconducting regions is also important to realize the strongly localized MBS in the non-superconducting channel.

Supplementary Note 4: Order parameters for the Majorana bound states
From the solution of the Bogoliubov-de Gennes equations of the planar Josephson junction with the realistic SkX texture, we identify the MBS in the quasiparticle spectrum and track the transition to the topological superconducting phase by computing two quantities: (i) the generalized Majorana polarization PM, proposed in Refs. [5,6], and (ii) the second derivative of the total density of states at zero energy, ∂ 2 D ∂E 2 , also called as the curvature of the density of states. The second derivative of the local density of states at a boundary site and at zero energy was used as an order parameter in Ref. [7] to track the topological superconducting transition. These quantities may provide additional insight to detect the MBS in experiments, besides the conventional zero-bias conductance peak [8,12] which often leads to ambiguity due to other possible zero-bias states in a superconductor [9]. To test the above two quantities in our set up, we first consider a chain with Rashba spin-orbit coupling and a uniform Zeeman exchange coupling, attached on the top and at the middle of a quasi-one-dimensional s-wave superconductor. The BdG quasiparticle spectrum of such a system, shown in Fig. S5a, depicts the emergence of the zero-energy MBS within a range of the chemical potential µ. The µ-variation of the two above-discussed quantities, |PM,1| and ∂ 2 D ∂E 2 , are plotted in, respectively, Fig. S5b and c. Interestingly, both these quantities reveal sharp jumps on entering the topological superconducting phase, similar to an order parameter in an ordinary phase transition. The Majorana oscillations, although not noticeable in the quasiparticle spectrum, is visible in both |PM,1| and ∂ 2 D ∂E 2 . In the infinite-length limit, the Majorana oscillations vanish and both these quantities reveal a quantized feature, as shown previously in Ref. 5. The sharp jumps, therefore, attest that these two quantities can reliably be used to track the topological superconducting transition in our computational framework. The superconducting potential is U = 2 meV and the pairing amplitude ∆i is treated self-consistently. The colorbar represents the Majorana polarization |PM,n| of the eigenstate with energy n. b The µ-variation of |PM,1|. c The µ-variation of ∂ 2 D ∂E 2 , the curvature of the density of states at zero energy.