Platform for braiding Majorana modes with magnetic skyrmions

After a decade of intense theoretical and experimental efforts, demonstrating braiding of Majorana modes remains an unsolved problem in condensed matter physics due to platform-speciﬁc challenges. In this work, we propose topological superconductor–magnetic multilayer heterostructures with on-chip microwave cavity readout as a platform for initializing, braiding, and reading out Majorana modes. Stray ﬁelds from a skyrmion in the magnetic layers can nucleate a vortex in the topological superconducting layer. Such a vortex is known to host Majorana bound states at its core. Through analytical calculations within London and Thiele formalisms, and through micromagnetic simulations, we show that our nucleation and braiding scheme can be effectively realized with a variety of existing options for magnetic and superconducting layers. Furthermore, we show that the coupling of the Majorana bound states to the electric ﬁeld of a resonator leads to an experimentally observable parity-dependent dispersive shift of the resonator frequency. Our work paves the way for realizing Majorana braiding in the near future.


I. INTRODUCTION
Demonstration of non-Abelian exchange statistics is one of the most active areas of condensed matter research and yet experimental realization of braiding of Majorana modes remains elusive [1,2].Most efforts so far have been focused on superconductor-semiconductor nanowire hybrids, where Majorana bound states (MBS) are expected to form at the ends of a wire or at boundaries between topologically trivial and nontrivial regions [3][4][5][6].Recently, it became clear that abrupt interfaces may also host topologically trivial Andreev states with experimental signatures similar to MBS [7,8], which makes demonstrating braiding in nanowire-based platforms challenging.Phase-controlled long Josephson junctions (JJs) open a much wider phase space to realize MBS with a promise to solve some problems of the nanowire platform, such as enabling zero-field operation to avoid detrimental flux focusing for in-plane fields [9,10].However, MBS in long JJs suffer from the same problems as in the original Fu-Kane proposal for topological insulator-superconductor JJs, such as poor control of flux motion along the junction and presence of sharp interfaces in the vicinity of MBS-carrying vortices which may host Andreev states and trap quasiparticles.For instance, MBS spectroscopy in both HgTe and InAs-based JJs Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license.Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.
shows a soft gap [11], despite a hard SC gap in an underlying InAs/Al heterostructure.
In the search for alternate platforms to realize Majorana braiding, spectroscopic signatures of MBS have been recently reported in scanning tunneling microscopy studies of vortex cores in iron-based topological superconductors (TSCs) [12].Notably, a hard gap surrounding the zero-bias peak at a relatively high temperature of 0.55 K, and a 5 K separation gap from trivial Caroli-de Gennes-Matricon (CdGM) states were observed [13,14].Moreover, vortices in a TSC can be field-coupled to a skyrmion in an electrically separated magnetic multilayer (MML) [15,16], which can be used to manipulate the vortex.This allows for physical separation of the manipulation layer from the layer wherein MBS reside, eliminating the problem of abrupt interfaces faced by nanowire hybrids and JJs.Finally, recent advances in the field of spintronics provide a flexible toolbox to design MMLs in which skyrmions of various sizes can be stabilized in zero external magnetic field and at low temperatures [16][17][18].Under the right conditions, stray fields from these skyrmions alone can nucleate vortices in the adjacent superconducting layer.In this paper, we propose TSC-MML heterostructures hosting skyrmion-vortex pairs (SVPs) as a viable platform to realize Majorana braiding.By patterning the MML into a track and by driving skyrmions in the MML with local spin-orbit torques (SOTs), we show that the SVPs can be effectively moved along the track, thereby facilitating braiding of MBS bound to vortices.
The notion of coupling skyrmions (Sk) and superconducting vortices (Vx) through magnetic fields has been studied before [15,[19][20][21][22][23][24][25].Menezes et al. [26] performed numerical simulations to study the motion of a skyrmionvortex pair when the vortex is dragged via supercurrents and  [27] proposed an analytical model for the motion of such a pair where a skyrmion and a vortex are coupled via exchange fields.Nothhelfer et al. [28] proposed using SVPs for Majorana braiding in a superconductor (nontopological) exchange coupled with a ferromagnet, where an SVP can host an MBS under favorable energy conditions.However, experimental evidence of an MBS-stabilizing topological phase has not yet been established in such heterostructures and no promising superconductor-ferromagnet combination has been proposed with existing materials.Moreover, practical constraints governing the dynamics of an SVP in the context of Majorana braiding such as maximum Sk-Vx force, estimated braiding time, and power dissipated during the braiding process remain largely unexplored in the literature.And most crucially, no in situ nondemolition experimental technique has been proposed to measure MBS in these heterostructures.In this paper, we propose replacing the conventional s-wave superconductor in such superconductor-ferromagnet heterostructures with a topological superconductor, which inherently guarantees MBS at its vortex cores [29] with ample experimental evidence supporting that prediction [12][13][14].Additionally, through micromagnetic simulations and analytical calculations within London and Thiele formalisms, we study the dynamics of a SVP subjected to external spin torques.We demonstrate that the SVP moves without dissociation up to speeds necessary to complete Majorana braiding (in the adiabatic limit) within the estimated quasiparticle poisoning time.We further eliminate the problem of in situ MBS measurements by proposing a novel on-chip microwave readout technique.By coupling the electric field of the microwave cavity to dipole-moments of transitions from Majorana modes to CdGM modes, we show that a topological nondemolition dispersive readout of the MBS parity can be realized.Moreover, we show that our platform can be used to make the first experimental observations of quasiparticle poisoning times in topological superconducting vortices.
The paper is organized as follows: In Sec.II we present a schematic and describe our platform.In Sec.III we present the conditions for initializing a skyrmion-vortex pair and discuss its equilibrium properties.In particular, we characterize the skyrmion-vortex binding strength.In Sec.IV we discuss the dynamics of a SVP in the context of braiding.Crucially, towards the end of the section, we discuss the conditions that enable adiabatic approximation for the braiding operation.Then in Sec.V, we present details of our microwave readout technique.Finally, we discuss the scope of our platform in Sec.VI.

II. PLATFORM DESCRIPTION
Our setup consists of a thin TSC layer that hosts vortices grown on top of a MML that hosts skyrmions, as shown in Fig. 1(a).A thin insulating layer (of the order a few nm) separates the magnetic and superconducting layers, ensuring electrical separation between the two layers and preventing any significant leakage currents across the layers. 1 Vortices in a TSC are expected to host MBS at their cores [12][13][14].Stray fields from a skyrmion in the MML can nucleate such a vortex in the TSC, forming a bound skyrmion-vortex pair under favorable energy conditions (see Sec. III).This phenomenon has been recently experimentally demonstrated in Ref. [16], where stray fields from Néel skyrmions in Ir/Fe/Co/Ni magnetic multilayers nucleated vortices in a bare niobium superconducting film.
The MML consists of alternating magnetic and heavy metal (HM) layers, as shown in Fig. 1(b).The size of a skyrmion in a MML is determined by a delicate balance between exchange, magnetostatic, anisotropy and Dzyaloshinskii-Moriya interaction (DMI) energies [34,35]and the balance is highly tunable, thanks to advances in spintronics [17,18,36].Given a TSC, this tunability allows us to find a variety of magnetic materials and skyrmion sizes that can satisfy the vortex nucleation condition [to be detailed in Eq. ( 1)].In the Appendix, we provide a specific example of FeTe x Se 1−x topological superconductor coupled with Ir/Fe/Co/Ni magnetic multilayers as one promising realization of our platform.
Due to large intrinsic spin-orbit coupling, a charge current through the heavy-metal layers of a MML exerts spin-orbit 1 Here we estimate the effect of leakage currents J l into the TSC layer.To setup a braiding current of J = 2 × 10 8 A/cm 2 in our MML, we must apply a voltage of V = ρJL ≈ 1 V across the MML, where ρ is the effective resistivity of the MML and L ∼ 50r sk is the length of the braiding track.At this bias voltage, ratio of Lorentz force F l = J l φ 0 d s to F max as defined in Eq. ( 6) is estimated as F l /F max = 10 −15 for 5 nm and F l /F max = 10 −10 for 1.2 nm insulating thickness [33].As such, leakage currents into TSC have negligible effect on SVPs.torques (SOTs) on the magnetic moments in the MML, which have been shown to drive skyrmions along magnetic tracks [37,38].In our platform, to realize Majorana braiding we propose to pattern the MML into a track as shown in Fig. 1(a) and use local spin-orbit torques to move skyrmions along each leg of the track.If skyrmions are braided on the MML track, and if skyrmion-vortex binding force is stronger than total pinning force on the SVPs, then the MBS-hosting vortices in TSC will closely follow the motion of skyrmions, resulting in the braiding of MBS.Here we assume that the braiding operation is performed at speeds much lower than the Fermi velocity of the TSC, so that we can invoke adiabatic principle and ensure that the finite speed of the braiding operation does not perturb the MBS.This condition is discussed more rigorously in Sec.IV.We also note in that section that there is an upper threshold speed with which a SVP can be moved.By using experimentally relevant parameters for TSC and MML in the Appendix, we show that our Majorana braiding scheme can be realized with existing materials.
We propose a nondemolition microwave measurement technique for the readout of the quantum information encoded in a pair of vortex Majorana bound states (MBS).A similar method has been proposed for the parity readout in topological Josephson junctions [39][40][41][42][43] and in Coulomb-blockaded Majorana islands [44].Dipole moments of transitions from MBS to CdGM levels couple dispersively to electric fields in a microwave cavity, producing a parity-dependent dispersive shift in the cavity's resonance frequency.Thus by probing the change in the cavity's natural frequency, the parity of the MBS can be inferred.Virtual transitions from the Majorana manifold to the excited CdGM manifold induced due to coupling to the cavity electric field are truly parity conserving, making our readout scheme a so-called topological quantum nondemolition technique [42,43].The readout scheme is explained in greater detail in Sec.V.
As discussed above, in our platform we consider coupling between a thin superconducting layer and magnetic multilayers.We note that in thin superconducting films, vortices are characterized by the Pearl penetration depth, given by = λ 2 /d s , where λ is the London penetration depth and d s is the thickness of the TSC film.Typically, these penetration depths are much larger than skyrmion radii r sk in MMLs of interest.Furthermore, interfacial DMI in MML stabilizes a Néel skyrmion as opposed to a Bloch skyrmion.So hereinafter, we only study coupling between a Néel skyrmion and a Pearl vortex in the limit r sk .

III. INITIALIZATION AND SKYRMION-VORTEX PAIR IN EQUILIBRIUM
Figure 1(c) illustrates the process flow of our initialization scheme.Skyrmions can be generated individually in MMLs by locally modifying magnetic anisotropy through an artificially created defect center and applying a current through adjacent heavy-metal layers [45].Such defect centers have been experimentally observed to act as skyrmion creation sites [46].When the TSC-MML heterostructure is cooled below the superconducting transition temperature (SC T C ), stray fields from a skyrmion in the MML will nucleate a vortex and an antivortex in the superconducting layer if the nucleation leads to a lowering in the overall free energy of the system [15].An analytical expression has been obtained for the nucleation condition in Ref. [47] ignoring contributions of dipolar and Zeeman energies to total magnetic energy: a Néel skyrmion nucleates a vortex directly on top of it if Here, d m is the effective thickness, M 0 is the saturation magnetization, A is the exchange stiffness, and K is the perpendicular anisotropy constant of the MML; α K and α A are positive constants that depend on the skyrmion's spatial profile (see the Appendix), r sk is the radius of the skyrmion in the presence of a Pearl vortex, 2 φ 0 is the magnetic flux quantum, and (ξ ) is the Pearl depth (coherence length) of the TSC.Although a complete solution of the nucleation condition must include contributions from dipolar and Zeeman energies to total energy of a MML, such a calculation can only be done numerically and Eq. ( 1) can still be used as an approximate estimate.
For the choice of materials listed in the Appendix, the left side of the equation exceeds the right side by 400%, strongly suggesting the nucleation of a vortex for every skyrmion in the MML.Furthermore, skyrmions in Ir/Fe/Co/Ni heterostructures have also been experimentally shown to nucleate vortices in niobium superconducting films [16].
We proceed to characterize the strength of a skyrmion (Sk)-vortex (Vx) binding force because it plays a crucial role in determining the feasibility of moving the skyrmion and the vortex as a single object.The spatial magnetic profile of a Néel skyrmion is given by where ζ = ±1 is the chirality and θ (r) is the angle of the skyrmion.For r sk , the interaction energy between a vortex and a skyrmion is given by [47]: where d = d m /r sk , J n is the nth-order Bessel function of the first kind, and R = r/r sk is the normalized horizontal displacement r between the centers of the skyrmion and the vortex.m z,θ (q) contains information about the skyrmion's spatial profile and is given by [47] where θ (x) is determined by the skyrmion ansatz.
We now derive an expression for the skyrmion-vortex restoring force by differentiating Eq. ( 2) with respect to r: ( For small horizontal displacements r r sk between the centers of the skyrmion and the vortex, we can approximate the 2 The radius of a skyrmion is not expected to change significantly in the presence of a vortex [47].We verified this claim with micromagnetic simulations.For the materials in the Appendix, when vortex fields are applied to a bare skyrmion, its radius increases by less than 10%.So, for numerical calculations in this paper, we use the bare skyrmion radius for r sk .

Sk-Vx energy as
with an effective spring constant Figures 2(a) and 2(b) show the binding energy and restoring force between a vortex and skyrmions of varying thickness for the materials listed in the Appendix.Here we used the domain wall ansatz for the skyrmion with θ (x) = 2 tan −1 [ sinh(r sk /δ) sinh(r sk x/δ) ], where r sk /δ is the ratio of skyrmion radius to its domain wall width and x is the distance from the center of the skyrmion normalized by r sk .As seen in Fig. 2(b), the restoring force between a skyrmion and a vortex increases with increasing separation between their centers until it reaches a maximum value F max and then decreases with further increase in separation.We note that F max occurs when the Sk-Vx separation is equal to the radius of the skyrmion, i.e., when R = 1 in Eq. ( 3): As the size of the skyrmion increases, the maximum binding force F max of the SVP increases.For a given skyrmion size, increasing the skyrmion thickness increases the attractive force until the thickness reaches the size of the skyrmion.Further increase in MML thickness does not lead to an appreciable increase in stray fields outside the MML layer and, as a result, the Sk-Vx force saturates.
It is important to note that stray fields from a skyrmion nucleate both a vortex and an antivortex (Avx) in the superconducting layer [15,[48][49][50].While the skyrmion attracts the vortex, it repels the antivortex.Equations ( 2) and (3) remain valid for Sk-Avx interaction, but switch signs.The equilibrium position of the antivortex is at the location where the repulsive skyrmion-antivortex force F Sk−Avx is balanced by the attractive vortex-antivortex force F V x−Avx [51,52].Figure 2(c) shows F V x−Avx against F Sk−Avx for the platform in the Appendix.We see that, for thicker magnets, the location of the antivortex is far away from that of the vortex, where the Avx can be pinned with artificially implanted pinning centers [53,54].For thin magnetic films, where the antivortex is expected to be nucleated right outside the skyrmion radius, we can leverage the Berezinskii-Kosterlitz-Thouless (BKT) transition to negate F V x−AVx for Vx-Avx distances r < [55][56][57][58].Namely, when a Pearl superconducting film is cooled to a temperature below T C but above T BKT , vortices and antivortices dissociate to gain entropy, which minimizes the overall free energy of the system [59].While the attractive force between a vortex and an antivortex is nullified, a skyrmion in the MML still attracts the vortex and pushes the antivortex towards the edge of the sample, where it can be pinned.Therefore we assume that the antivortices are located far away and neglect their presence in our braiding and readout schemes.

IV. BRAIDING
Majorana braiding statistics can be probed by braiding a pair of MBS [1] which involves swapping positions of the two vortices hosting the MBS.We propose to pattern the MML into interconnected Y junctions, as shown in Fig. 3 (also see the corresponding skyrmion braiding Supplemental video [32]) to enable that swapping.Ohmic contacts in HM layers across each leg of the Y-junctions enable independent application of charge currents along each leg of the track.These charge currents in-turn apply spin-orbit torques on the adjacent magnetic layers and enable skyrmions to be moved independently along each leg of the track.As long as skyrmion and vortex move as a collective object, braiding In step II, γ 1 and γ 2 are brought close to rf lines by applying charge currents from C to A and D to B, respectively.γ 1 and γ 2 are then initialized by performing a dispersive readout of their parity (see Sec. V).Similarly, γ 3 and γ 4 are initialized after applying charge currents along P to R and Q to S, respectively.In step III, γ 2 is moved aside to make room for γ 4 by applying currents from B to X followed by applying currents from X to C. In step IV, γ 4 is braided with γ 2 by applying currents along S to X and X to B. Finally, in step V, the braiding process is completed by bringing γ 2 to S by applying currents from A to X and from X to S. Parities (i.e., fusion outcomes) of γ 1 and γ 4 , and γ 3 and γ 2 are then measured in step VI.Fusion outcomes in each pair of MBS indicate the presence or absence of a fermion corresponding to a parity of ±1 [30,31].(b) Initial position of the skyrmions labeled A and B in the micromagnetic simulation for skyrmion braiding (see the Appendix).[(c)-(h)] Positions of the two skyrmions at the given times as the braiding progresses.Charge current j = 2 × 10 12 A/m 2 was applied.This skyrmion braiding is further illustrated in the accompanying Supplemental video [32].
of skyrmions in the MML leads to braiding of MBS hosting vortices in the superconducting layer.Below we study the dynamics of a SVP subjected to spin torques for braiding.We calculate all external forces acting on the SVP in the process and discuss the limits in which the skyrmion and the vortex move as a collective object.
For a charge current J in the HM layer, the dynamics in the magnetic layer is given by the modified Landau-Lifshitz-Gilbert (LLG) equation [60,61]: (7) where we have included damping-like term from the SOT and neglected the field-like term because it does not induce motion of Néel skyrmions for our geometry. 3Here, γ is the gyromagnetic ratio, α is the Gilbert damping parameter, and H eff is the effective field from dipole, exchange, anisotropy, and DMI interactions.p = sgn( SH ) Ĵ × n is the direction of polarization of the spin current, where SH is the spin Hall angle, Ĵ is the direction of charge current in the HM layer, and n is the unit vector normal to the MML.η = h SH /2eM 0 d m quantifies the strength of the torque, h is the reduced Planck's constant and e is the charge of an electron.
Assuming skyrmion and vortex move as a collective object, semiclassical equations of motion for the centers of mass of the skyrmion and the vortex can be written using collective coordinate approach as done in Ref. [27]: and where R sk (R vx ), m sk (m vx ), and q sk (q vx ) are the position, mass, and chirality of the skyrmion (vortex).k is the effective spring constant of the Sk-Vx system, given in Eq. ( 5).F SOT = π 2 γ ηr sk sJ × n is the force on a skyrmion due to spin torques in Thiele formalism, where s = M 0 d m /γ is the spin density [64,65].The third term on the right side of Eq. ( 8) gives Magnus force on the skyrmion, with G sk = 4π sq sk ẑ, and the fourth term characterizes a dissipative force due to Gilbert damping.Similarly, the second term on the right side of Eq. ( 9) gives the Magnus force on the vortex with G vx = 2π sn vx q vx ẑ, with n vx being the superfluid density of the TSC, and the third term characterizes viscous force with friction coefficient α vx .U sk, pin (U vx, pin ) gives the pinning potential landscape for the skyrmion (vortex).The last term in Eq. ( 9) represents the restoring force on a vortex due to its separation from a skyrmion and is valid when | R sk − R vx |< r sk .Here, k is the effective spring constant characterizing the Sk-Vx force, as given by Eq. ( 4).We consider steady-state solutions of the equations of motion assuming that the skyrmion and the vortex are bound.We discuss conditions for the dissociation of a SVP later.For a given external current J, velocity v of a SVP in steady state is obtained by setting Rsk = Rvx = 0 and Ṙsk = Ṙvx = Ṙ in Eqs. ( 8) and ( 9): In general, the SVP moves at an angle ϕ relative to F SOT due to Magnus forces on the skyrmion and the vortex, with Armed with the above equations, we extract some key parameters that determine the feasibility of our braiding scheme.First, if F SOT from external currents is unable to overcome the maximum pinning force on either the skyrmion (F pin,sk ) or the vortex (F pin,vx ), the SVP will remain stationary.This gives us a lower threshold J − on the external current which is obtained by weighing F SOT against the pinning forces: Second, once the SVP is in motion, drag and Magnus forces that act on the skyrmion and the vortex are proportionate to their velocity.If the net external force on a vortex in motion is larger than the maximum force with which a skyrmion can pull it (F max ), then the skyrmion and the vortex dissociate and no longer move as a collective object.This sets an upper bound v + on the SVP speed which can be obtained by balancing F max with the net force from Magnus and drag forces on the vortex.This maximum speed plays a key role in determining whether our braiding and readout scheme can be completed within the quasiparticle poisoning time: An upper bound on the SVP speed implies an upper bound J + on the external current which can be obtained by putting v + in Eq. ( 10): Another parameter of critical importance, the distance of closest approach between two skyrmion-vortex pairs (r min ) plays a crucial role in achieving significant overlap of the MBS wave functions centered at the vortex cores and is given by balancing Sk-Vx attractive force by Vx-Vx repulsive force: Finally, the power P dissipated in heavy metal layers due to Joule heating from charge currents has to be effectively balanced by the cooling rate of the dilution refrigerator: where n hm is the number of heavy-metal layers, L (W ) is the length (width) of the active segment of the MML track, t hm is the thickness of each heavy-metal layer, and ρ hm is the resistivity of a heavy-metal layer.By applying a current J − < J < J + locally in a desired section of the MML track, each SVP can be individually addressed (see the Supplemental skyrmion braiding video [32]).For the materials listed in the Appendix, the maximum speed v + with which a SVP can be moved is over 1000 m/s.At this top speed, SVPs can cover the braiding distance [the sum of the lengths of the track in steps I-VI of Fig. 3(a)] of 50r sk in about 0.15 ns, but the process generates substantial Joule heating.At a reduced speed of 0.25 m/s, SVPs cover that distance in 7 µs generating 30 µW of heat during the process, which is within the cooling power of modern dilution refrigerators.SVPs can be braided at faster speeds if the dilution fridges can provide higher cooling power or if the resistivity of heavy metal layers in the MML can be lowered.Although quasiparticle poisoning times in superconducting vortices have not been measured yet, estimates in similar systems range from hundreds of microseconds to seconds [66][67][68].Our braiding time falls well within such estimates for quasiparticle poisoning times, indicating the viability of our platform.Furthermore, the ability to easily tune braiding times in our platform by varying magnitude of currents in heavy metal layers can be used to investigate the effects of quasiparticle poisoning on the braiding protocol.
As will be shown in Sec.V, Vx-Vx distances <10ξ should be sufficient to perform a dispersive readout of MBS parity in adjacent vortices.For the materials listed in the Appendix, the distance of closest approach between two vortices is r min = 40 nm.The shape of the MML track further limits how close two vortices can be brought together [see step II in Fig. 3(a)].With the geometry of the track taken into account, Vx-Vx distance less than 10ξ can still be easily achieved, enough to induce a detectable shift in the cavity's resonance frequency during the dispersive readout.
We now discuss the conditions for invoking adiabatic approximation for the braiding operation.The dragging of MBS by skyrmions during braiding of SVPs can be treated as an external time-dependent perturbation on top of the static superconducting manifold [28].If the rate at which these SVPs are moved is much slower than the rate governing electron dynamics within the superconducting ground state, we can invoke the adiabatic approximation argument.This translates to the condition that the braiding time be much larger than the timescale governing electron transitions in the time-independent superconducting manifold: L/v h/δE or L/v hE F / 2 , where v is the velocity of a SVP, h is the reduced Planck constant, E F is the Fermi energy, is the superconducting gap, and δE is the energy separation between the MBS manifold and the excited CdGM manifold in the spectrum of stationary TSC vortex.Using numbers from our platform of L ∼ 10r sk , v ≈ 1 m/s, δE ∼ 2 /E F ≈ 0.9 meV, the left side of the inequality (L/v) exceeds the right side ( h/δE ) by five orders of magnitude!This guarantees that the effect of time dependence of skyrmion's fields on MBS due to motion of a SVP can be safely neglected despite all the complexities of the combined TSC-MML system.Furthermore, clever universally applicable braiding protocols have been proposed, such as the so-called "bang bang" protocol which further lowers diabatic errors during the braiding process by two to three orders compared with the naive average velocity protocol [69].
Figures 3(b)-3(h) show the results of micromagnetic simulation of braiding skyrmions in a smaller section of the MML (for computational reasons) for the example platform.The details of the simulation are given in the Appendix.The simulation results demonstrate the effectiveness of using local SOT to move individual skyrmions and realize braiding.Finally, as discussed in this section, due to the strong skyrmion-vortex binding force, vortices hosting MBS in the TSC will braid alongside the skyrmions.

V. READOUT
Quantum information is encoded in the charge parity of MBS hosted in a pair of vortices which we propose to readout with a dispersive measurement technique.Figure 4(a) summarizes our readout scheme: for a two-vortex system, the top figure shows single-particle energy levels and the bottom figure shows many-body energy levels as the vortices are brought close to each other.Energy separation E M in the figure between the ground state and the Majorana occupancy state increases with decreasing distance between the vortices due to hybridization of MBS at the vortex cores.In the dispersive limit, a microwave cavity electric field can drive virtual transitions from the ground-state Majorana manifold to the excited CdGM manifold (only one CdGM level, labeled 1, is considered in the figure).The transitions allowed by selection rules, labeled ω −M,1 and ω M,1 are shown in the many-body spectrum.Each of these virtual transitions causes a statedependent dispersive shift in the cavity's natural frequency and the parity of the vortex pair can be inferred from relative change in cavity frequency, as explained in the caption for Fig. 4(a).Note that each of the allowed transitions is truly parity conserving since microwaves cannot change the number of fermions.Since parity states are true eigenstates (as opposed to approximate eigenstates) of the measurement operation, our readout scheme can be dubbed as a topological quantum nondemolition technique [42,43].We now proceed to calculate dipole coupling strengths of the allowed transitions to cavity electric field and the corresponding dispersive shift.
In BCS mean-field theory, the coupling to electric field can be described by the Hamiltonian where E(t ) = E 0 cos ωt is the microwave-induced timedependent electric field which is approximately uniform over the scale of the vortices. 4The electric field couples to the dipole operator d of the electronic states in the vortices.We have written it in terms of the electron field operator in Nambu spinor notation, ˆ = (ψ ↑ , ψ ↓ , ψ † ↓ , −ψ † ↑ ) T .The Pauli matrix n represent oscillations in the wave functions of a clean system.In a disordered (real) system the oscillations are expected to be smeared out.The inset shows the probability density for the MBS hosted by a vortex pair 400 nm apart.The simulation was done for an effective 2D model (a 1000 × 600 nm 2 rectangle) of a 3D topological insulator surface, see Refs.[70][71][72].We used ξ = 15 nm, vortex radius r = ξ , and E F = 1.125 in the vortex.
τ z acts on the particle-hole indices.We use effective singleband model to evaluate eigenstates, wave functions and field operators of the topological superconductor surface [29].At low energies, field operators for the TSC can be expanded in terms of its eigenstates as where γ1 and γ2 are the Majorana operators for vortices 1 and 2, and ˆ ( †) 1 is the annihilation (creation) operator for the lowest CdGM state.The corresponding four-component spinor wave functions multiply the operators in Eq. (18).We now proceed to derive expressions for dispersive shift produced in cavity's resonance frequency due to the aforementioned electric dipole coupling.
At low frequencies much below the level spacing δE of the vortex quasiparticle bound states, ω δE / h, the cavity microwave field does not excite the quasiparticle states of the vortices.We shall also assume that these quasiparticle states are not occupied, for example due to quasiparticle poisoning.Under these conditions, the vortex pair stays in its ground state manifold consisting of the two states of unoccupied or occupied nonlocal MBS.With sufficiently weak microwave driving we can use dispersive readout to measure the charge parity σ z = i γ1 γ2 [43,73].The dispersive Hamiltonian of the resonator-vortex pair system reads [43] where â, â † are the harmonic-oscillator annihilation and creation operators for the resonator.The MBS parity-dependent dispersive frequency shift is where we denote and ω is the resonator bare frequency, E 0 is the electric field amplitude, and δE is the energy gap separating the MBS from the first-excited CdGM mode.We ignore here the exponentially small energy splitting between the MBS [29], which would give subleading corrections to χ ; we will see that χ itself will be exponentially small in the vortex separation (due to the parity-sensitive transition dipole matrix elements Evaluating the dipole transition matrix elements d 1,±M microscopically is somewhat involved since proper screening by the superconducting condensate needs to be carefully accounted for and is beyond the BCS mean-field theory [74][75][76][77][78]. Nevertheless, to estimate d 1,±M we can use Eq. ( 17) by replacing r ≈ l ẑ in it, with l ≈ a B being the effective distance to the image charge in the superconductor and ẑ the surface normal vector [74].Here a B denotes the Bohr radius.We evaluate the dimensionless matrix elements of the effective dipole "charge" d • ẑ/l by using a numerical simulation of the Majorana and CdGM states in a double vortex system depicted in Fig. 4(b).The numerical simulations will be detailed in a future presentation [72].
In Fig. 4(b) we plot the parity-sensitive coupling g 2 n that largely determines the dispersive shift χ , Eq. ( 20).
The coupling decays exponentially in the distance between the Majorana pair since, at large distances, the local coupling to microwaves, Eq. ( 17), cannot access the nonlocal parity information.Nevertheless, we find that even a relatively distant vortex pair can provide a parity-dependent shift g 2 n ∼ 10 −2 (elE 0 ) 2 .Since the relevant dipole moment is normal to the superconductor surface, we can couple to the dipole by using a microwave resonator above the surface, producing a large perpendicular electric field.With a resonator zero-point voltage V 0 ≈ 100 µV at a ≈10 nm distance from the vortices, we obtain elE 0 ≈ 1 µeV × (l/Å) ≈ 2.4 × 10 2 h MHz × (l/Å).(We estimate that such high zero-point voltages can be achieved in high-inductance resonators [79].)Taking a low-lying CdGM state with δE ≈ 10 µeV , we obtain χ/2π ∼ 20 MHz × (l/Å) 2 where l ≈ a B 1 Å is the typical dipole size [74].We thus see that the MBS vortex parity measurement is well within standard circuit QED measurement capabilities [73].We note that the above estimate does not include the resonant enhancement, the second factor in Eq. ( 20), which may further substantially increase the frequency shift.Finally, we note that the dipole operator d also has a nonzero diagonal matrix element d M in the Majorana state [77], leading to a term E 0 • d M σ z (â + â † ) in Eq. ( 19).This term in principle allows one to perform longitudinal readout of the MBS parity.However, making longitudinal readout practical may require parametric modulation of the coupling, in our case d M , which may be difficult [42,73].

VI. CONCLUSION
Measuring braiding statistics is the ultimate method to conclusively verify the existence of non-Abelian excitations.We proposed a unified platform to initialize, braid, and readout Majorana modes, avoiding abrupt topological-trivial interfaces at each stage.We derived general expressions for braiding speeds with spin currents, distance of closest approach between two Majorana modes, and the resultant dispersive shift in cavity resonance frequency.We showed that our setup can be readily realized with existing options for TSC and MML materials.

APPENDIX: EXAMPLE PLATFORM
While our platform can be realized with a variety of existing combinations of MML and TSC, we choose Ir/Fe/Co/Ni and FeTe 0.55 Se 0.45 for the respective roles [12,16].For the TSC, we used λ = 500 nm, ξ = 15 nm and chose a thickness of d s = 50 nm.This gives a Pearl length of = 5 µm.We propose to use 10 layers of Ir/Fe/Co/Ni heterostructure with 0.5 nm of Fe and 0.5 nm of Co per layer giving an effective magnet thickness of d m = 10 nm.
We performed micromagnetic simulations using MUMAX [80,81] with parameters similar to those in Ref. [16]: M 0 = 1450 emu/cc, A = 13.9 pJ/m, K = 2000 kJ/m 3 , D = 3 mJ/m 2 , and B ext = 180 mT, where D is the Dzyaloshinskii-Moriya interaction (DMI) constant and B ext is the external magnetic field.We modified the values of uniaxial anisotropy constant and DMI constant from those in Ref. [16] to stabilize skyrmions of large size.We also used external magnetic field in the simulation to stabilize isolated skyrmions.In the experimental setup, this magnetic field can be applied via exchange interaction in magnetic multilayers.Thanks to versatile spintronics tools [17,18,36,82], these modifications are quite feasible to realize.Using a grid size of 2 × 2 × 2 nm 3 , isolated skyrmions of size r sk = 35 nm were stabilized in the simulation.
To demonstrate skyrmion braiding, the MML was patterned into a Y-shaped track and a spin current corresponding to a charge current of j = 2 × 10 12 A/m 2 in HM layers was locally applied in each active leg of the MML track.To save computational resources, vacuum layers (representing heavymetal layers) have not been included in MUMAX simulations.The inclusion of vacuum layers does not affect the conclusions of this paper in any way.Here, we set spin Hall angle, SH = 0.15, and the ratio of field-like SOT term to dampinglike SOT term to 0.1.
For quantitative analysis of Eqs. ( 12)-( 16), we assume clean limit and neglect pinning forces on skyrmion and vortex.We include drag and Magnus forces on skyrmion, and drag force on vortex.We neglect the Magnus force on the vortex beause it is expected to be small compared with its drag force [83].We model the drag force on the vortex using friction coefficient α vx = η vx d s , where η vx is the Bardeen-Stephen viscous drag coefficient given by η vx = φ 2 0 /(2πξ 2 ρ n ) [84,85].Here, ρ n is the normal-state resistivity of the TSC.Following parameters were used for calculations: Gilbert damping parameter for the MML α = 0.1 and ρ n = 10 −7 m, giving a η vx = 3 × 10 −8 N s/m 2 .We propose that each leg of the MML track in Fig. 1(a) be 10r sk wide and 10r sk long.This implies that skyrmions have to cover a total distance of 50r sk for the entirety of the braiding process in Fig. 3(a).We used the resistivity of iridium, ρ hm = 5 × 10 −7 m to calculate the Joule heat dissipated during the braiding process.
FIG. 1.(a) Schematic of the Majorana braiding platform.Magnetic multilayer (MML) is patterned into a track and is separated from TSC by a thin insulating layer.Green lines represent on-chip microwave resonators for a dispersive parity readout setup.The left inset shows a magnified view of an SVP and the right inset shows the role of each layer.(b) Expanded view of the composition of an MML.(c) Process flow diagram for our Majorana braiding scheme.Here, T c is superconducting transition temperature and T BKT is Berezinskii-Kosterlitz-Thouless transition temperature for the TSC

FIG. 3 .
FIG. 3. (a) Schematic of our braiding process: manipulations of four skyrmions in the MML track are shown.MBS at the centers of vortices bound to each of these skyrmions are labeled γ 1 -γ 4 .Ohmic contacts in HM layers of the MML are shown in brown and rf readout lines are shown in green.II-VI show the steps involved in braiding γ 2 and γ 4 .In step II, γ 1 and γ 2 are brought close to rf lines by applying charge currents from C to A and D to B, respectively.γ 1 and γ 2 are then initialized by performing a dispersive readout of their parity (see Sec. V).Similarly, γ 3 and γ 4 are initialized after applying charge currents along P to R and Q to S, respectively.In step III, γ 2 is moved aside to make room for γ 4 by applying currents from B to X followed by applying currents from X to C. In step IV, γ 4 is braided with γ 2 by applying currents along S to X and X to B. Finally, in step V, the braiding process is completed by bringing γ 2 to S by applying currents from A to X and from X to S. Parities (i.e., fusion outcomes) of γ 1 and γ 4 , and γ 3 and γ 2 are then measured in step VI.Fusion outcomes in each pair of MBS indicate the presence or absence of a fermion corresponding to a parity of ±1[30,31].(b) Initial position of the skyrmions labeled A and B in the micromagnetic simulation for skyrmion braiding (see the Appendix).[(c)-(h)] Positions of the two skyrmions at the given times as the braiding progresses.Charge current j = 2 × 10 12 A/m 2 was applied.This skyrmion braiding is further illustrated in the accompanying Supplemental video[32].

FIG. 4 .
FIG. 4. (a) Schematic of our readout process.When two vortices are brought close, microwave transitions can be dispersively driven from the MBS to the excited hybridized CdGM levels (only level 1 is shown).Parity of the Majorana mode can be inferred from the difference in the cavity frequency shift produced by ω −M,1 and ω M,1 transitions [see Eq. (20)].The allowed fermion parity-conserving transitions are shown in both single-particle and many-particle representations.In the latter, dashed and solid lines denote states in the two fermion parity sectors.The transition of frequency ω −M,1 (blue arrows) corresponds to breaking a Cooper pair and exciting the MBS and CdGM levels (MBS being initially unoccupied).When the MBS is occupied, the transition of frequency ω M,1 (red arrow) excites the MBS quasiparticle into the CdGM level.The dipole transition matrix elements are different for the two processes, enabling parity readout.(b) MBS-parity sensitive dipole transition strength versus vortex pair separation.We denote by g 2n = (|E 0 • d n,−M | 2 − |E 0 • d n,M | 2 )the dipole transition strength between the Majorana level and the nth CdGM level.We plot the dimensionless strength normalized by U = e|E 0 |l.As expected from MBS hybridization, g n decays approximately exponentially in the distance between the two vortices.Oscillations in g 2 n represent oscillations in the wave functions of a clean system.In a disordered (real) system the oscillations are expected to be smeared out.The inset shows the probability density for the MBS hosted by a vortex pair 400 nm apart.The simulation was done for an effective 2D model (a 1000 × 600 nm 2 rectangle) of a 3D topological insulator surface, see Refs.[70][71][72].We used ξ = 15 nm, vortex radius r = ξ , and E F = 1.125 in the vortex.