1 Introduction↩︎

The potential of quantum computation has sparked significant research interest in recent decades. However, achieving practical quantum computation necessitates the implementation of fault-tolerant techniques. A promising approach is topological quantum computation (TQC)  [1]–[3]. The cornerstone of TQC lies in topological excitations known as non-Abelian anyons, which emerge in many-body systems. These non-Abelian anyons form qubits that are topologically protected, rendering them resilient against local noise and enabling extremely low error rates that meet the threshold requirements for fault-tolerant quantum computation. In TQC, quantum gates are implemented through the braiding of non-Abelian anyons, referred to as braiding gates. Braiding results are obtained by measuring the signals associated with information after the fusion of non-Abelian anyons—a process involving coupling and hybridization of these anyons. This approach offers a promising route towards realizing fault-tolerant quantum computation and paves the way for more robust and reliable applications in quantum computing.

In recent years, the Majorana zero mode (MZM) has emerged as a prominent non-Abelian anyon for realization in solid-state systems [4]–[10]. Signatures of MZMs have been observed in one-dimensional (1D) systems, such as nanowires and ferromagnetic atomic chains [11]–[16]. However, since MZMs cannot be braided within a single 1D structure, the creation of a 1D network becomes crucial for implementing braiding operations [17], [18]. In topological semiconductor-superconductor heterostructure networks, braiding can be achieved by locally controlling electron density using electrode gates [19]–[22].

Compared to 1D or quasi-1D systems, two-dimensional (2D) topological superconductors (TSCs) offer a more natural platform for the braiding of anyons. Various proposals have been put forward to realize MZMs in 2D systems, including fractional quantum Hall effects [23], intrinsic \(p\)-wave superconductivity [24], topological insulators proximitized by s-wave superconductors [25], and 2D semiconductor/superconductor hybrid structures in proximity to magnetic insulators or subjected to external magnetic fields [26], [27]. In particular, superconductor vortex is a natural host for MZM and other sub-gap states [16], [28]. MZM signatures within vortices have been observed in hybrid 2D systems [29], [30] and iron-based superconductors over the past decade [31]–[38], opening up new avenues for topological quantum computation with vortices.

The observation of 2D MZMs has stimulated extensive investigations into corresponding studies on braiding methods [39]–[48]. For example, various approaches have been proposed including nanoscale tips or tip-shaped pinning sites [40], [41], electrically/magnetically controlled braiding techniques [42]–[44], [49], and optically mediated methods [50], [51]. However, the majority of these methods fail to satisfy the essential criteria required for successful implementation of MZM braiding protocols. In particular, gate-controlled techniques relying on electrostatics are incompatible with 2D semiconductor-superconductor hybrids or iron-based superconductors due to the screening effect of metal electrodes, which prevents effective electron density tuning. Meanwhile, tip-based approaches face integration challenges.

In this study, we propose a novel braiding scheme to address the challenges of 2D MZMs using spintronics devices. Similar to the control elements in Field-Programmable Topological Array (FPTA) [52], magnets with nanometer scales can be individually configured and controlled by electrical current pulses. We implement spin-transfer-torque (STT) devices, which have been extensively studied in the field of spintronics [53], to enable independent, simultaneous, and localized programmability of various properties including electromagnetic fields, synthetic spin-orbit fields, and superconductivity.

Leveraging STT devices for braiding offers significant advantages in terms of efficiency, control speed, and applicability compared to existing proposals. Previous studies by Fatin et al. in Ref. [54] and Zhou et al. in Ref. [55] have demonstrated the theoretical feasibility of manipulating and braiding MZMs using magnetic tunneling junctions (MTJs) in 2DEG-superconductor hybrids. Inspired by these groundbreaking works, we propose the utilization of STT-matrix as a novel approach to braid MZMs within vortex-type TSCs. It is noteworthy that our work differs from Ref. [54] in addressing different material systems, employing different numerical methods, expanding to time-dependent simulations and braiding dynamics.

The article is structured as follows: Section 2 provides an overview of the fundamental concepts and system structures of MZM braiding devices. In Section 3, we present a detailed description of the model and method based on time-dependent Ginzburg-Landau (TDGL) theory and time-dependent Bogoliubov-de Gennes (TDBdG) equations, with particular emphasis on integrating the STT-matrix into the TDGL equations. In Section 4, we simulate and demonstrate the manipulation of MZMs by the STT-matrix, including their ability to drive single vortices (Subsec. 4.1) and analyze braiding dynamics (Subsec. 4.2). Additionally, we propose a fusion method for MZMs in the STT setup (4.3). Finally, conclusions are drawn in Section 5.

2 Fundamental concept and architectural framework↩︎

Our proposals are formulated based on three fundamental principles. Firstly, within TSCs, a vortex hosts an MZM. As mentioned in the first section, \(p\)-wave superconductor vortices have been proposed as one of the most promising platforms for hosting MZMs, and experimental evidence of MZMs in vortices has been reported. Secondly, a vortex can be individually manipulated by a local magnetic field. Previous studies have demonstrated the creation and manipulation of vortices using field coils on superconducting quantum interference devices [56]–[58]. Thirdly, the local magnetic field can be provided by the stray field of STT device matrix [54], [55]. Combining these principles implies that braiding the MZMs is equivalent to braiding the vortices in \(p\)-wave superconductors, which is further equivalent to configuring the STT-matrix.

An STT device typically comprises an MTJ multiple layer structure, consisting of two ferromagnetic layers separated by a thin non-magnetic spacer layer. One of the ferromagnetic layers, referred to as the “reference layer," maintains a fixed magnetization direction during operation. The other ferromagnetic layer, known as the”free layer", can have its magnetization direction altered by the spin-polarized current flowing through the structure. The change in magnetization direction of the free layer corresponds to the alteration of the STT state, which is achieved through the interaction between the angular momentum of polarized current and magnetic moments of the free layer [59]. Depending on the direction of current, either alignment or anti-alignment can be established between the magnetization directions of the free layer and reference layer. Consequently, programmability of STT states arises from distinct net stray fields exhibited during alignment or anti-alignment. The size of the STT device can now be scaled down to a few tens of nanometers, enabling the integration of billions of these devices into a single chip. As a result, the STT-matrix can provide exceptional precision and scalability for vortex manipulation.

The STT-matrix, as shown in Fig. 1 (a), functions as the control layer within the system. Each individual STT can be independently controlled using conventional addressing techniques [60]. By altering the current direction, specific STTs can be selectively activated (ON) or deactivated (OFF), enabling precise alignment or anti-alignment of the free layer with the reference layer for manipulation and propulsion of vortices within the TSC layer. In Figs 1 (b)-(c), we demonstrate MZMs within TSC vortices controlled by STTs (see formula in Sec. 3.3). When only one STT is activated, a vortex is induced by this particular STT, resulting in the emergence of a point-like MZM within the vortex. Since MZMs always appear in pairs, another MZM manifests as a 1D edge mode [Fig. 1 (b) and (d)], akin to edge MZMs observed near a s-wave superconductor and ferromagnetic insulator in semiconductor systems [61]. Conversely, activating two STT devices generates a pair of vortices hosting point-like MZMs [Fig. 1 (c) and (e)].

3 Model and Method↩︎

In this section, our objective is to establish the comprehensive framework model as described in the preceding section and devise a methodology for simulating MZMs braiding.

Firstly, vortices are induced by STT-matrix, which is modeled via a Finite Element Analysis method (see Appendix 8 and Ref. [62]). The vortex dynamics are governed by the TDGL equations, which are a powerful tool for studying the dynamics of vortex systems in type-II superconductors with mixed phases [63]. While the TDGL theory was originally formulated for dirty or gapless superconductors near the critical temperature \(T_c\) [64], its generic phenomenological features enable its validity across a broader range of conditions. Specifically, for iron-based superconductors, where the penetration length \(\lambda\) is typically several orders of magnitude larger than the coherence length \(\xi\) [65]–[67], the dirty limit construction of TDGL is well-justified. To incorporate the effect of the STT-matrix, we modify the TDGL equations by including local magnetic fields associated with it. This modification enables us to simulate and analyze individual vortex dynamics effectively [68], [69].

The model of TSCs is an effective 2D representation. The MZMs are obtained by solving the static model of the TSCs. By incorporating vortex dynamics derived from TDGL equations, the braiding dynamics of MZMs can be solved using TDBdG equations.

It should be noted that although we utilize a 2D model as an illustrative example in our simulations, we expect that our approach can also effectively capture the behavior of quasi-2D TSCs or the surface states of bulk topological materials.

3.1 TDGL equations↩︎

The TDGL equations describe a system of nonlinear equations governing the behavior of the superconductor order parameter \(\Psi\) and electromagnetic potential \((\Phi, \mathbf{A})\). These equations can be expressed as [68], [70]:

\[\begin{align} \label{eqTDGLDimension} \frac{\hbar^2}{2mD}(\frac{\partial}{\partial t}+i\frac{q}{\hbar}\Phi)\Psi & = & -\frac{1}{2m}(\frac{\hbar}{i}\nabla-q\mathbf{A})^2\Psi+\alpha\Psi\nonumber\\ & & -\beta|\Psi|^2\Psi,\\ \sigma(\frac{\partial \mathbf{A}}{\partial t}+\nabla\Phi) & = & \frac{q\hbar}{2mi}(\Psi^*\nabla\Psi-\Psi\nabla\Psi^*)\\ & & -\frac{q^2}{m}|\Psi|^2 \mathbf{A}-\frac{1}{\mu_0}\nabla\times(\nabla\times \mathbf{A}-\mathbf{B}).\nonumber \end{align}\tag{1}\] where \(\hbar\) represents the reduced Planck’s constant, \(m\) denotes the mass of a Cooper pair, and \(q=2e\) signifies the charge carried by a Cooper pair (where \(e\) is the elementary charge carried by an electron). The parameters \(\alpha\) and \(\beta\) correspond to the phenomenological aspects of Gibbs free energy in the GL theory. Specifically, \(\alpha(T)\) is a temperature-dependent parameter given by \(\alpha(0)(1-T/T_c)\), where \(T_c\) represents the critical temperature of the superconductor. On the other hand, \(\beta\) remains as a constant value. Additionally, \(D\), \(\sigma\), and \(\mu_0\) are material-specific parameters: \(D\) refers to diffusion coefficient; \(\sigma\) corresponds to conductivity; and finally, \(\mu_0\) stands for permeability of free space. Furthermore, we consider an external magnetic field denoted as \(\mathbf{B}=\mathbf{B}(x,y,t)\) which varies with time and can be controlled using the STT-matrix.

Related to GL equations, there are two important characteristic lengths: the London penetration depth, denoted as \(\lambda=\sqrt{{m\beta}/(q^2 \mu_0 \alpha)}\), and the GL coherence length, denoted as \(\xi={\hbar}/{\sqrt{2m\alpha}}\). The ratio of these two lengths, referred to as the Ginzburg-Landau parameter, \(\kappa={\lambda}/{\xi}\), serves as an indicator for distinguishing between type I superconductors (\(\kappa<1/\sqrt{2}\)) and type II superconductors (\(\kappa>1/\sqrt{2}\)).

The dimensionless coordinates and time are rescaled by \(\lambda\) and \(t_0=\xi^2/D\), respectively. Consequently, the dimensionless derivatives can be expressed as \[\begin{align} \tilde{\nabla}=\lambda \nabla,\quad \partial_\tau=t_0 \partial_t. \end{align}\] The dimensionless quantities can be defined by introducing the corresponding unit, \[\begin{align} \begin{aligned} \Psi_0 &= \sqrt{\frac{\alpha}{\beta}}, & \sigma_0 &= \frac{1}{\mu_0 D \kappa^2}, & \Phi_0 &= \frac{\hbar D \kappa}{q\xi^2}, \\ \mathbf{A_0} &= \frac{\hbar}{q\xi}, & \mathbf{B_{0}} &= \frac{\hbar}{q \kappa \xi^2}, & \end{aligned} \end{align}\]

and we have \[\begin{align} \begin{aligned} \psi &=\Psi/\Psi_0,~ & \tilde{\sigma} &=\sigma/\sigma_0, & \phi &=\Phi/\Phi_0,\\ \tilde{\mathbf{A}} &=\mathbf{A}/\mathbf{A_0}, & \mathbf{b} &=\mathbf{B}/\mathbf{B_{0}}. & \end{aligned} \end{align}\]

After performing algebraic calculations and applying a gauge transformation (see 6 for detailed procedures), we can obtain a dimensionless form as follows: \[\begin{align} \partial_\tau \psi & = & -(\frac{i}{\kappa}\nabla+\mathbf{A})^2\psi+\psi-|\psi|^2\psi,\\ \sigma \partial_\tau \mathbf{A} & = & \frac{1}{2i\kappa}(\psi^*\nabla\psi-\psi\nabla\psi^*)-|\psi|^2 \mathbf{A}\nonumber\\ & & -\nabla\times(\nabla\times \mathbf{A}-\mathbf{b}), \end{align}\] where the redefined variables \(\tilde{\nabla}\), \(\tilde{\sigma}\) and \(\tilde{\mathbf{A}}\) are now denoted as \(\nabla\), \(\sigma\) and \(\mathbf{A}\).

Extended boundary conditions are employed, specifically designed for numerical simulations [68], which involve enlarging the superconductor with a non-superconducting and insulating material. Mathematically, the TDGL equations are modified by introducing a region-dependent function. Consequently, the dimensionless boundary conditions can be expressed as \[\begin{align} \psi|_{\partial \Omega} & = & 0,\\ \nabla\times \mathbf{A}|_{\partial \Omega} & = & \mathbf{b}|_{\partial \Omega}, \end{align}\] in which \(\partial \Omega\) represents the external boundary, and \(\mathbf{b}\) is the dimensionless magnetic field. The initial condition for the TDGL equations is \[\begin{align} \psi|_{\tau=0}=1,\mathrm{\, \in\, superconductor\, region }. \end{align}\]

Given the boundary and initial conditions, the quantities \(\psi=\psi(t)\) and \(\mathbf{A}=\mathbf{A}(t)\) can be determined by solving the TDGL equations. Here, our primary focus lies on the observable quantity, namely the Cooper pair density \(\rho\), which can be expressed as \(\rho=|\Psi|^2=|\psi|^2\Psi_0^2\).

3.2 Model of the STT-matrix↩︎

In the previous subsection, we present the formalism of TDGL theory. By imposing appropriate boundary and initial conditions, numerical methods can be employed to obtain the dynamics of vortices. Unlike prior studies on vortex lattices formed in uniform external magnetic fields, our focus here lies on the vortex lattice generated by the STT-matrix, in order to achieve individual vortex manipulation. To model the STT-matrix, we abstract it as a series of localized magnetic fields acting upon the superconductor. Through integration with the STT-matrix, the TDGL equations enable simulation of vortex dynamics driven by STTs.

The localized magnetic field generated by an activated STT device can be approximately represented by a Gaussian magnetic field, which corresponds to the external magnetic field \(b(x,y)\) in the TDGL equations, \[\begin{align} \label{eqbSTT} b_{i,j}(x,y)=b_{0}e^{-[(x-X_{i,j})^2+(y-Y_{i,j})^2]/R_v^2}, \end{align}\tag{2}\] where \(b_{0}\) is the maximum magnetic field, \((X_{i,j}, Y_{i,j})\) is the central coordinate of the STT in the \(i^{\rm th}\) row and the \(j^{\rm th}\) column. In the case of a Gaussian-shaped magnetic field, \(R_v\) represents the characteristic radius of an individual STT’s magnetic field. It should be noted that while selecting a Gaussian-shaped magnetic field for each STT is not crucial, any model that can capture the localized nature of the STT’s magnetic field can be employed for vortex control.

The distance between the nearest STTs in the STT-matrix, denoted as \(L_{grid}\), is defined as \(L_{grid} = X_{i,j}-X_{i-1,j}=Y_{i,j}-Y_{i,j-1}\) for a square STT-matrix, as illustrated in Fig. 1 (a). The net magnetic field from the entire STT-matrix, denoted as \(b_{\rm net}\) at time \(\tau\), is the aggregate of all individual STT magnetic fields, \[\begin{align} b_{\rm net}(\tau)=\sum_{i,j} s_{i,j} (\tau) b_{i,j}(x,y). \end{align}\] The time dependence arises from the status function \(s_{i,j} (\tau)\), which is a binary function representing the ON or OFF state of the \(i,j\)-th STT, \[\begin{align} s_{i,j} (\tau) = \begin{cases} 1, & \mathrm{ON}\\ 0, & \mathrm{OFF} \end{cases}. \end{align}\] The status function \(s_{i,j} (\tau)\) serves as the crucial link between the control layer and the magnetic field texture, effectively demonstrating the programmable nature of the STT-matrix.

3.3 Model of TSCs and the TDBdG equations↩︎

To describe the TSCs, we employ the 2D effective model, which offers the advantage of low computational cost while still capturing the essential properties of TSCs [44], [71], [72].

The Hamiltonian is given by \[\begin{align} H_{\rm TSC}=c_{x,y}^\dagger \mathcal{H}_0 c_{x,y} +\Delta(x,y) c_{x,y}^\dagger c_{x,y}^\dagger + \mathrm{H.c.}, \end{align}\] where \(c_{x,y}/c_{x,y}^\dagger\) represent the annihilation/creation operators for particles at position \((x, y)\), \(c_{x,y}=(c_{x,y,\uparrow},c_{x,y,\downarrow})^T\), and \(\Delta(x, y)\) denotes the superconducting pair potential. Here, \(\mathcal{H}_0\) corresponds to the Hamiltonian governing surface states in a 3D topological insulator: \[\begin{align} \mathcal{H}_0=\hbar v_F (\boldsymbol{\sigma}\times\mathbf{k})\cdot \mathbf{e}_z+\frac{W\lambda}{2}\sigma_z|\mathbf{k}|^2, \end{align}\] with \(v_F\) representing Fermi velocity. We introduce a mass term \(W\) to address fermion doubling issues and restore a finite gap in Brillouin zone boundaries [72]. In these equations, \(\boldsymbol{\sigma}=(\sigma_x,\sigma_y,\sigma_z)\) refers to Pauli matrices acting on spin space and \(\mathbf{k}=(k_x,k_y)\) represents momenta on the surface.

Expressed in the Nambu basis \(C_{x,y}=(c_{x,y},c_{x,y}^\dagger)^T\), the TSC Hamiltonian can be written compactly as \[\begin{align} H_{\rm TSC}=C_{x,y}^\dagger \mathcal{H}_{\rm TSC} C_{x,y}, \end{align}\] where the BdG Hamiltonian \(\mathcal{H}_{\rm TSC}\) is given by \[\begin{align} \mathcal{H}_{\rm TSC}=h_1 \nu_z + h_2 +\Delta_R \nu_x - \Delta_I \nu_y, \end{align}\] with \(h_1 = i \hbar v_F (\sigma_y \partial_x-\sigma_x \partial_y)\), and \(h_2 = -\frac{W\lambda}{2}\sigma_z (\partial_x^2+\partial_y^2)\). Here, \(\Delta\) represents a complex pairing potential with real and imaginary parts denoted as \(\Delta_R\) and \(\Delta_I\), respectively. The Pauli matrices \(\nu_{x,y,z}\) act on the particle-hole space.

The dimensionless form of the BdG Hamiltonian is required to coherently integrate into the dimensionless TDGL equations, which can be deduced using the same procedure as in the previous section. The resulting final dimensionless Hamiltonian is given by: \[\begin{align} \label{eqBdGHamiltonian} \tilde{\mathcal{H}}_{\rm TSC} & = & i \tilde{v}_F (\sigma_y \tilde{\partial}_x - \sigma_x \tilde{\partial}_y) \nu_z - \tilde{W}\sigma_z (\tilde{\partial}_x^2\nonumber\\ & & + \tilde{\partial}_y^2) +\tilde{\Delta}_R\nu_x - \tilde{\Delta}_I\nu_y, \end{align}\tag{3}\] where the dimensionless quantities are defined as follows: \[\begin{align} \begin{aligned} \tilde{\mathcal{H}}_{\rm TSC} &=\frac{H_{\rm TSC}}{\Delta_0}, & \tilde{v}_F &=\frac{\hbar v_F}{\lambda\Delta_0}, \\ \tilde{W} &=\frac{W}{2\lambda\Delta_0}, & \tilde{\Delta}_{R,I} &=\frac{{\Delta}_{R,I}}{\Delta_0}. & \end{aligned} \end{align}\] Here, \(\Delta_0\) represents the amplitude of superconducting pair potential without vortex.

With the dimensionless \(\tilde{\mathcal{H}}_{\rm TSC}\), the TDBdG Schrödinger equation can be also written in a dimensionless form \[\begin{align} \label{eqTDBdG} \tilde{\mathcal{H}}_{\rm TSC}(\tau)\phi_{\rm TSC}(\tau)=id\partial_\tau \phi_{\rm TSC}(\tau), \end{align}\tag{4}\]

where \(d=\frac{\hbar D}{\xi^2 \Delta_0}=\frac{\hbar}{\Delta_0 t_0}\) is a dimensionless constant that can be absorbed by rescaling the time variable to \(\tau'=\tau/d\). After this rescaling, the equation becomes \[\begin{align} \tilde{\mathcal{H}}_{\rm TSC}(\tau')\phi_{\rm TSC}(\tau')=i\partial_{\tau'} \phi_{\rm TSC}(\tau'). \end{align}\] This formulation represents the TDBdG formula used for numerical calculations of MZM braiding dynamics.

To incorporate vortices into the TDBdG Hamiltonian, it is necessary to modify the superconducting pair potential near a vortex pinned by the \(i,j\)-th STT to \[\begin{align} \Delta_{i,j}(\boldsymbol{r},\tau)=\Delta_0 \tanh\frac{|\boldsymbol{r}-\boldsymbol{r_{i,j}}(\tau)|}{\xi}\mathrm{exp}(-i\mathrm{Arg}(\boldsymbol{r}-\boldsymbol{r_{i,j}}(\tau))) \end{align}\] where \(\boldsymbol{r_{i,j}}(\tau)=(x_{i,j}(\tau),y_{i,j}(\tau))\) represents the central coordinates of the vortex pinned by the \(i,j\)-th STT. Renormalize \(\Delta_0\) to unity, the overall superconducting pair potential \(\Delta(\boldsymbol{r},\tau)\) is given by \[\begin{align} \label{eqDelta} \Delta(\boldsymbol{r},\tau)=\prod_{i,j} \Delta_{i,j}(\boldsymbol{r},\tau), \end{align}\tag{5}\] where the product is taken over all vortices.

3.4 Simulation Procedure↩︎

After introducing the basic components, a comprehensive simulation of the MZM braiding dynamics can be conducted.

The initial step involves designing the braiding paths and determining the ON/OFF operation time sequences of the STTs based on specific braiding process. Subsequently, utilizing the defined braiding paths and STT time sequences, the TDGL equations yield the Cooper pair density \(\rho(\tau)\), thereby enabling analysis of vortex dynamics. The trajectories of vortices \(\{\mathbf{r}_{i,j}(\tau) = (x_{i,j}(\tau), y_{i,j}(\tau)) \mid i,j=1,2,3,\ldots \}\) are then extracted from \(\rho(\tau)\) for all \(\tau\) and passed to TDBdG equations through Eq. (5 ). By solving static BdG equations, we identify states of MZMs. Finally, by selecting an initial MZM state, we can obtain corresponding braiding dynamics along the designated paths.

4 MZMs braiding↩︎

4.1 Single Vortex Dynamics↩︎

To demonstrate the vortex-manipulation capability of STTs, we first simulate the motion of a single vortex under the dragging of STTs.

As depicted in Fig. 2 (a), a vortex initially immobilized by device STT-A (indicated by the left dashed circle) undergoes drift towards the position of STT-B (right dashed circle) upon altering the ON/OFF states of both STT-A and STT-B, as shown in Fig. 2 (b). Eventually, the vortex reaches and settles into a position near STT-B, representing another stable state. We refer to this process as an elementary move. The TDGL simulation also indicates that the vortex dissipation time (after the corresponding STT being turned off) is much longer than the braiding time scale.

We investigated the dynamics of vortex elementary motion with varying distances \(L_s\) between the two STTs. The displacements of the vortex, denoted as \(D(\tau)\), were first extracted from the simulation results obtained by solving the TDGL equations. Subsequently, we determined the velocities of the vortex, represented as \(v(\tau)(=dD(\tau)/d\tau)\). As illustrated in Fig. 2 (c), for small \(L_s\), immediate capture of the vortex by STT-B occurred after its release by STT-A due to its proximity within the effective region of stray field generated by STT-B, resulting in a rapid acceleration in velocity. Conversely, for large \(L_s\), initially outside reach of STT-B’s stray field, a transient wandering phase was observed before eventual capture by STT-B.

Notably, the maximum velocity (\(v_{max}\)) reaches saturation when \(L_s\) exceeds a certain threshold value. This behavior is further elucidated in Fig. 2 (d), where data points corresponding to different values of \(L_s\) are collected to demonstrate two distinct patterns with fitted lines. The maximum velocity of vortex motion is intricately linked to both the speed of braiding operations and adiabaticity. Besides, we will also see in Subsec. 4.2 that large accelerations of vortex (small \(L_s\)) leads to a high braiding error.

As mentioned earlier, for a type-II superconductor, the GL parameter \(\kappa\) should satisfy \(\kappa>1/\sqrt{2}\). In Fig. 2, we set the value of \(\kappa\) to 4. It is worth noting that our model remains valid across a wide range of \(\kappa\), encompassing both the clean superconductor regime (small values of \(\kappa\)) and the iron-based superconductor regime (large values of \(\kappa\)).

4.2 Braiding↩︎

The previous subsection demonstrated the STT’s ability to drive the vortex, and it is natural to extend this capability in an STT-matrix. By properly arranging the status of the STTs, it becomes possible to freely move the vortex along the trajectory defined by of the STT lattice. The collection of these statuses is referred to as a time sequence, which consists of status functions for all STTs \(\{s_i(\tau)\}\) defined in section 3. This implies that any braiding operation can be implemented by utilizing the interaction between STTs and vortices. To illustrate this point, it would be beneficial to begin with an introduction of non-Abelian anyon braiding.

In the TQC proposal, the computational process of manipulating non-Abelian anyons is represented by their space-time trajectories, known as worldlines. These worldlines can be formally described using the braiding group framework. According to this theory, a general braiding operation \(\mathbb{B}\) can be decomposed into a combination of elementary operations \(\mathbb{B}=b_ib_j \cdots b_k\), where \(i,j,\cdots k\in\pm\{1,\ldots,n-1\}\) (with \(n\) being the number of non-Abelian anyons). The elementary operation \(b_i\) (\(i>0\)), referred to as a generator of the braiding group, represents a counterclockwise exchange between the \(i\)-th and \((i+1)\)-th non-Abelian anyons, while the corresponding operator \(b_{-i}\) is the inverse of \(b_i\) representing a clockwise braiding. This decomposition ensures that any general braiding operation can be achieved by performing all necessary generators, which can be facilitated by digitally controlling the STT-matrix.

In practice, the braiding operation requires two essential preparatory steps. Firstly, it is necessary to establish a mapping between the braiding worldlines and the corresponding braiding paths in the STT lattice. Secondly, one must design appropriate time sequences for the STTs along these braiding paths, as illustrated in Figs. 3 (a)-(b).

To demonstrate the braiding of MZMs, we calculate the dynamics of MZM braiding using the TDBdG equations. The system under consideration consists of four vortices with even interval, with each vortex hosting an individual MZM denoted as \(\gamma_{i=1,2,3,4}\). For illustrative purposes, we choose the braiding operation \(\mathbb{B}_X=b_2b_2\) (i.e., exchange the positions of \(\gamma_2\) and \(\gamma_3\) twice, which is an effective \(X\)-gate), which is selected for its simplicity and ability to yield a distinct signature in the braiding dynamics.

In the language of TQC, four MZMs can be encoded into one logical qubit defined in terms of fermionic modes, specifically the even subspace. For example, we have \[\begin{align} \{|\bar{0}\rangle=|00\rangle, |\bar{1}\rangle=|11\rangle=c_1^\dagger c_2^\dagger |00\rangle\}, \end{align}\] where \(c_1=\frac{1}{2}(\gamma_1+i\gamma_2)\), \(c_2=\frac{1}{2}(\gamma_3+i\gamma_4)\) represent the fermionic operators. The parity of a pair of MZMs can be defined as the fermion number associated to the fermionic operator. Here, \(|\bar{0}\rangle\) and \(|\bar{1}\rangle\) denote the logical qubit. By performing a braiding operation \(\mathbb{B}_X\), MZMs \(\gamma_2\) and \(\gamma_3\) acquire a \(\pi\)-phase due to their non-Abelian statistics [44], [45], [73], resulting in \(\gamma_2\to-\gamma_2\) and \(\gamma_3\to-\gamma_3\), while leaving \(\gamma_{1}\) and \(\gamma_{4}\) unchanged. Consequently, the associated logical qubit undergo a transformation from \(|\bar{0}\rangle \to |\bar{1}\rangle\) and \(|\bar{1}\rangle \to |\bar{0}\rangle\), which essentially corresponds to a quantum NOT gate operation on the logical qubit representation. This becomes evident by noting that after braiding, the fermionic operators transform as follows: \[\begin{align} c_1=\frac{\gamma_1+i\gamma_2}{2}\to \frac{\gamma_1-i\gamma_2}{2}=c_1^\dagger, \\ c_2=\frac{\gamma_3+i\gamma_4}{2}\to \frac{-\gamma_3+i\gamma_4}{2}=-c_2^\dagger. \end{align}\]

The nontrivial transformation of the operators can also be elucidated in the picture of single-particle spectrum within the MZM subspace. The states in the MZM subspace can be obtained by solving the static BdG Hamiltonian Eq. (3 ). Theoretically the MZMs are the states appear inside the superconductor gap with zero energy. However, due to the overlaps of their wavefunctions, the energies split to values close to zero. The particle-hole symmetry ensures the simultaneous existence of \(|\pm E\rangle\) states. The splitting of energy depends on the distance between the MZMs. Therefore, by singling out the two pairs of particle-hole symmetric subgap eigenstates, we have identified and collected a set of states constituting the MZM subspace, denoted as \(\{|\pm E_{i}\rangle (i=1,2)\}\).

After the braiding operation \(\mathbb{B}_{X}\), the states undergoes a transition: \({ |\pm E_i \rangle} \rightarrow {| \mp E_i \rangle}, (i = 1 , 2)\). This nontrivial transition can be understood as follows. Suppose we start from the many-body state \(|\bar{0}\rangle=|00\rangle\), which annihilated by the operators \(c_1\) and \(c_2\) (\(c_1|\bar{0}\rangle=c_2|\bar{0}\rangle=0\)). This corresponds to the occupied single-particle states being \(|-E_1\rangle\) and \(|-E_2\rangle\), as they each satisfy the annihilation conditions \(c_1|-E_1\rangle=0\) and \(c_2|-E_2\rangle=0\). Take \(|-E_1\rangle\) as an example. It should satisfy the constraints before and after the braiding, in order to deduce the final state. Concretely, the initial state is annihilated by \(c_1\) (\(c_1|-E_1\rangle=0\)), but after the braiding, the operator becomes \(c_1\to c_1^\dagger\), thus requiring the final state be annihilated by \(c_1^\dagger\) (\(c_1^\dagger|f\rangle=0\)), which means the final state is \(|f\rangle=|+E_1\rangle\). Collecting each case of these states, it is implied that the braiding operation \(\mathbb{B}_X\) causes a state transition as \(|\pm E_i\rangle\to|\mp E_i\rangle,i=1,2\).

To characterize the fidelity of braiding operation, we focus on two physical quantities. The first physical quantity is the state transition probability, defined as \[\begin{align} P_{\Phi}(\tau)=|\langle\phi(\tau)|\Phi(0)\rangle|^2, \end{align}\] where \(|\phi(\tau)\rangle\) represents the braided state and \(|\Phi(0)\rangle\) denotes the initial MZM state labeled as \(|\pm E_i\rangle, i=1,2\) at time \(\tau=0\). This probability quantifies the likelihood of a braided state projecting onto the encoded space constituted by MZMs. The second physical quantaty is the leakage out of the MZM subspace during operation, which can be denoted as \(1-P_{\rm MZM}\), where \(P_{\rm MZM}=\sum_{\Phi} P_{\Phi}(\tau), \Phi \in \{|\pm E_i \rangle | i=1,2\}\).

By applying the TDBdG equations Eq. (4 ) and choosing the initial state as \(\phi(\tau=0)=|-E_1\rangle\), the entire braiding process of the braided state \(\phi(\tau)\) can be calculated. The state transition probabilities throughout the entire \(\mathbb{B}_X\) braiding process are depicted in Figs. 3 (c) and (d). It is evident that the process initiates with a high \(P_{|-E_1\rangle}\) and ends with a high \(P_{|+E_1\rangle}\). Ideally, \(P_{|+E_1\rangle}\) should reach unity after \(\mathbb{B}_X\), while \(P_{\rm MZM}\) should remain at unity for the entirety of the operation.

However, the braiding fidelity can be influenced by two factors. Firstly, a finite coupling strength between MZMs, denoted as \(E_M\), can transform MZMs into fermionic modes and cause loss of their non-Abelian statistics. The finite value for \(E_M\) primarily results in a reduction in \(P_{|+E_1\rangle}\). Secondly, the nonadiabatic effects arising from vortex dynamics significantly impact the system. The velocity nonuniformity induced by vortex motion can be regarded as a perturbation to MZMs, particularly during the processes of vortex releasing and repinning. When \(L_s\) is small, the rapid changes in vortex dynamics [as illustrated in Fig. 2 (c)] lead to pronounced fluctuations in velocity. This results in substantial leakage from the MZM computational subspace into the quasi-continuum, thereby reducing both \(P_{\rm MZM}\) and \(P_{|+E_1\rangle}\).

To investigate the impact of finite \(E_M\) and nonadiabaticity on the braiding process in the STT-matrix scheme, we examined the dependence of \(P_{|+E_1\rangle}\) and \(P_{\rm MZM}\) on \(L_v\) and \(L_s\). Here, \(L_v\) represents the interval between \(\gamma_i\) and \(\gamma_{i+1}\), serving as an indicator for \(E_M\). As depicted in Fig. 3 (e), when \(L_v\) is small, there is a significant reduction in \(P_{|+E_1\rangle}\) while maintaining a high value for \(P_{\rm MZM}\). Conversely, when \(L_v\) becomes large enough to neglect the effect of \(E_M\), both \(P_{|+E_1\rangle}\) and \(P_{\rm MZM}\) approach unity. With further increase in \(L_v\), gradual drops are observed in both \(P_{|+E_1\rangle}\) and \(P_{\rm MZM}\) due to an increase in elementary moves. The value of \(L_s=3\), remains constant for all data points presented in Fig. 3 (e). In Fig. 3 (f), we fix \(L_v=18\), while varying \(L_s\). It can be clearly seen that both \(P_{|+E_1\rangle }\) and \(P_{\text{MZM}}\) increase with increasing \(L_s\). This behavior arises from the fact that smaller \(L_s\) will result in larger accelerations experienced by vortices and cause larger leakage (recalling the analysis in Sec. 4.1). Additionally, smaller values of \(L_s\), require more elementary moves.

4.3 Fusion↩︎

The fusion of anyons [1]–[3] is a crucial step in bringing them together to decode information in the computational basis, which plays an integral role in the read-out process of TQC. Moreover, the utilization of fusion alone (without braiding) can also facilitate measurement-based TQC [74].

In the STT-matrix braiding scheme, the fusion operation can be easily performed by bringing together the vortices hosting MZMs that are intended to be fused, as illustrated in Figs. 4 (a) and (b). The controlled movement of these two MZMs can be achieved through the application of STTs (see Sec. 4.1), ultimately leading them to approach each other. Subsequently, the MZMs are pinned at a pair of neighboring STTs set in an ON state.

Remarkably, we observed that even if one of the two STTs is switched OFF, its associated vortex will still migrate towards the remaining ON-STT. In essence, this implies that two vortices can symmetrically localize around a single ON-STT. However, due to their repulsive interaction [63], these vortices remain separated by a distance of around \(\lambda\).

After the STT operations, the fusion of the two MZMs within the neighboring vortices occurs, resulting in an observable finite energy splitting denoted as \(E_M\). The dependence of \(E_M\) on the distance between MZMs, referred to as \(L_M\), is illustrated in Fig. 4 (c). Due to the relatively large Fermi velocity in metal, the MZM wavefunction does not show oscillations along with \(L_M\) [75].

The quantum information stored in the fused MZMs can be extracted through charge flipping or other related physical effects, as proposed in previous works on read-out methods for 2D-TSC systems [44], [76]. These methods are compatible with the STT-matrix scheme; however, their implementation is beyond the scope of this study.

5 Conclusion↩︎

In conclusion, we investigate the braiding of MZMs in 2D TSC vortices using an STT-matrix platform. To validate our proposal, we develop a theoretical simulation method based on the TDGL theory and TDBdG equations. By incorporating the STT-matrix into the TDGL equations through modeling local magnetic stray fields, we simulate individual MZM manipulation within vortices driven by STTs. Subsequently, this method is employed to study vortex dynamics induced by STTs and further explore MZM braiding dynamics using TDBdG equations. We demonstrate that manipulating the STT-matrix enables achieving MZM movement, braiding, and fusion operations. Additionally, we show that nonadiabaticity and finite coupling of MZMs significantly impact their braiding dynamics, which can be adjusted in our proposed scheme. Our approach is general and independent of concrete TSC models, such as the Fu-Kane model, semiconductors with Rashba spin-orbit interaction, and iron-based superconductor model. We believe the studies on the MZMs braiding in the vortex states of TSC will advance TQC.

Acknowledgment. This work was supported by the Innovation Program for Quantum Science and Technology (Grant No. 2021ZD0302401), the Hunan Provincial Science Foundation for Distinguished Young Scholars (Grant No.2021JJ10043), and the Open Research Fund from State Key Laboratory of High Performance Computing of China (HPCL) (Grant No.201901-09 and No.2023-KJWHPCL-07.

G.Y.H and J.B.F contribute equally to this work.

6 Derivation of TDGL equations↩︎

The derivation of dimensionless TDGL equations is presented in detail in this section. We commence with the dimensional TDGL equations [68], [70]. \[\begin{align} \label{eqTDGLDimension2} \frac{\hbar^2}{2mD}(\frac{\partial}{\partial t}+i\frac{q}{\hbar}\Phi)\Psi & = & -\frac{1}{2m}(\frac{\hbar}{i}\nabla-q\mathbf{A})^2\Psi+\alpha\Psi\nonumber\\ & & -\beta|\Psi|^2\Psi,\\ \sigma(\frac{\partial \mathbf{A}}{\partial t}+\nabla\Phi) & = & \frac{q\hbar}{2mi}(\Psi^*\nabla\Psi-\Psi\nabla\Psi^*)\\ & & -\frac{q^2}{m}|\Psi|^2 \mathbf{A}-\frac{1}{\mu_0}\nabla\times(\nabla\times \mathbf{A}-\mathbf{B}).\nonumber \end{align}\tag{6}\] Parameter definitions can be found in the main texts.

We define: \[\begin{align} \begin{aligned} \partial_t &=\frac{D}{\xi^2}\partial_\tau, & \nabla &=\frac{\tilde{\nabla}}{\lambda}, & \Phi &=\frac{\hbar D \kappa}{q\xi^2} \phi, & \Psi &=\sqrt{\frac{\alpha}{\beta}}\psi,\\ \sigma &=\frac{\tilde{\sigma}}{\mu_0 D \kappa^2}, & \mathbf{A} &=\frac{\hbar}{q\xi}\tilde{\mathbf{A}}, & \mathbf{B} &=\frac{\hbar}{q \kappa \xi^2} \mathbf{b}. & \end{aligned} \end{align}\] Substituting the parameters above into Eq. 1 (or Eq. 6 ), we get a dimensionless form: \[\begin{align} (\frac{\partial}{\partial \tau}+i\kappa\phi)\psi & = & -(\frac{i}{\kappa}\tilde{\nabla}+\tilde{\mathbf{A}})^2\psi+\psi-|\psi|^2\psi,\\ \sigma(\frac{\partial \tilde{\mathbf{A}}}{\partial \tau}+\tilde{\nabla}\phi) & = & \frac{1}{2i\kappa}(\psi^*\tilde{\nabla}\psi-\psi\tilde{\nabla}\psi^*)-|\psi|^2 \tilde{\mathbf{A}}\nonumber\\ & & -\tilde{\nabla}\times(\tilde{\nabla}\times \tilde{\mathbf{A}}-\mathbf{b}) \end{align}\]

To further simplify the TDGL equations, we perform the gauge transformation on the electromagnetic potential. \[\begin{align} \begin{array}{rcl} \psi' &=& \psi e^{i\kappa \chi}, \\ \tilde{\mathbf{A}}' &=& \tilde{\mathbf{A}}+\tilde{\nabla} \chi, \\ \phi' &=& \phi-\partial_\tau \chi, \end{array} \end{align}\] where \(\chi=\chi(x,y,z,t)\) is an arbitrary real function. Fixing the gauge by choosing \(\phi=\partial_\tau \chi\), the electric potential \(\phi'\) will vanish globally.

Since TDGL equations are gauge invariant, they now have the form: \[\begin{align} \partial_\tau \psi & = & -(\frac{i}{\kappa}\tilde{\nabla}+\tilde{\mathbf{A}})^2\psi+\psi-|\psi|^2\psi,\\ \sigma \partial_\tau \tilde{\mathbf{A}} & = & \frac{1}{2i\kappa}(\psi^*\tilde{\nabla}\psi-\psi\tilde{\nabla}\psi^*)-|\psi|^2 \tilde{\mathbf{A}}\nonumber\\ & & -\tilde{\nabla}\times(\tilde{\nabla}\times \tilde{\mathbf{A}}-\mathbf{b}). \end{align}\] Redefine \[\begin{align} \begin{array}{rcl} \tilde{\nabla} &\to& \nabla, \\ \tilde{\mathbf{A}} &\to& \mathbf{A}, \\ \tilde{\sigma} &\to& \sigma, \end{array} \end{align}\] and get the final form of TDGL equations in the main text: \[\begin{align} \partial_\tau \psi & = & -(\frac{i}{\kappa}\nabla+\mathbf{A})^2\psi+\psi-|\psi|^2\psi\\ \sigma \partial_\tau \mathbf{A} & = & \frac{1}{2i\kappa}(\psi^*\nabla\psi-\psi\nabla\psi^*)-|\psi|^2 \mathbf{A}\nonumber\\ & & -\nabla\times(\nabla\times \mathbf{A}-\mathbf{b}) \end{align}\]

7 Experimental Feasibility↩︎

As commercially mature spintronic components, STT devices provide a robust foundation for the proposed MZM braiding scheme. STT technology has advanced rapidly since the first experimental demonstrations of STT-induced magnetization switching in metal multilayers [77] and magnetic tunneling junctions (MTJs) [78]. Its commercialization in non-volatile memories (MRAMs) [79] and scalability to gigabit integration on 22 nm/28 nm nodes [80], [81] underscores the practical potential of an STT-matrix for MZM control.

Regarding key device parameters and material compatibility, the proposed scheme is compatible with experimentally validated 2D-TSCs. Integration of the STT-matrix with 2D-TSC materials can be achieved via flip-chip packaging technology  [82], [83] or polymer-assisted 2D material transfer methods [84]. Validated 2D-TSC materials include iron-based superconductors [31]–[38] and topological insulator/s-wave superconductor hybrids [29], [30]. Such materials host vortex-bound MZMs and support the vortex dynamics modeled here. The diameter of STT-matrix elements can be scaled down to sub-10 nm [85], matching the superconducting coherence length of iron-based TSCs [31], [86]. This dimensional compatibility ensures localized stray fields can precisely pin individual vortices. It should be noted that current STT material fabrication processes are primarily optimized for fast memory devices; therefore, they require optimization to generate moderate magnetic fields suitable for MZM control.

For dynamic performance, the magnetization flip time of STT free layers has reached 200 ps\(\sim\)​10 ns [87]. This guarantees that the single vortex manipulation time is significantly shorter than the theoretical coherence time of the topological qubits [88]. Furthermore, thermal effects are negligible: STT switching energy is \(\sim 10\) fJ per operation at room temperature [89], [90], increasing only slightly at cryogenic temperatures while remaining manageable. Given that dilution refrigerators provide \(\sim 10\;\mu\)W of cooling power at mK temperatures, the STT-matrix can operate at gigahertz frequencies without disrupting superconductivity.

8 FEA Modeling of STT devices↩︎

To characterize the stray field of the STT device, we conducted a Finite Element Analysis (FEA) simulation for a single STT.

In this FEA method, the STT device is modeled as a simplified cylindrical structure consisting of three layers: a 30-nm-thick reference layer on top, followed by a 2-nm-thick tunneling layer separating it from the bottom 20-nm-thick free layer. The radius of the cylindrical structure is 90 nm. It should be noted that real STT devices often have complex multiple-layer structures; however, our FEA model captures the essential local distribution of stray fields.

To simulate the distribution of stray field in the ON and OFF states, we manually manipulate the magnetization direction of the free layer to align it either parallel or antiparallel with the reference layer. The linecuts of the simulated \(B_z\)-component of the stray field along a horizontal axis located 40 nm below the STT is illustrated in Fig. 5 (a). It can be seen that for the ON state, there is a significant presence of stray field at the STT center, which can be approximately fitted by a Gaussian function in the TDGL simulation. Conversely, in the OFF state, magnetic field strength is considerably weaker and can be regarded as zero during the TDGL simulation.

From the linecuts, we observe a deviation in the negative part of \(B_z\) from the Gaussian fitted lines adjacent to the peak. This deviation arises due to the closed nature of magnetic field lines in the STT device. However, this deviation does not cause any other consequences, such as the formation of a vortex-antivortex pair, since the reversed magnitude of \(B_z\) is too small. To demonstrate this, we simulate a vortex by accounting the negative part of \(B_z\) using FEA data in TDGL equations. The resulting vortex shown in Fig. 5 (b) exhibits no difference compared to the results obtained using Gaussian magnetic fields shown in the main text.

9 An example of STT-matrix device↩︎

Here, we demonstrate the fabrication and measurements of a very-first-stage STT device. As the SEM image shown in Fig. 6 (a), an STT-Matrix with eight STT-units are fabricated. The STT-arrays are based on a multilayer MTJ material stack Fig. 6 (b). Top electrodes and a (global) bottom electrode are connected to the STT-array. In fact, these STT-devices are similar to those used in non-volatile memories. We characterized the STT-matrix using a write current pulse with various amplitude \(I_{w1},I_{w2},\ldots\) and a following read current \(I_{\rm read}\) of 10 \(\mu\)A in amplitude Fig. 6 (c). The write and read pulse durations are 1 ms and 300 ms, respectively. As shown in Fig. 6 (d), the resistance of the measured STT is bistable, and can be switched when the write pulse goes beyond a critical point. This is the so called STT-effect. The switching of the STT resistance indicates a switch of the free-layer magnetization. As just a demonstration, the measurements are performed at room temperature. Parameters like device compatibility at low temperature, the strength of stray-fields and the effect of write pulse Joule heat etc. need to be further investigated.

References↩︎