## Abstract

The iron-based superconductor FeTe* _{x}*Se

_{1−x}is one of the material candidates hosting Majorana vortex modes residing in the vortex cores. It has been observed by recent scanning tunneling spectroscopy measurement that the fraction of vortex cores having zero-bias peaks decreases with increasing magnetic field on the surface of FeTe

*Se*

_{x}_{1−x}. The hybridization of two Majorana vortex modes cannot simply explain this phenomenon. We construct a three-dimensional tight-binding model simulating the physics of over a hundred Majorana vortex modes in FeTe

*Se*

_{x}_{1−x}. Our simulation shows that the Majorana hybridization and disordered vortex distribution can explain the decreasing fraction of the zero-bias peaks observed in the experiment; the statistics of the energy peaks off zero energy in our Majorana simulation are in agreement with the experiment. These agreements lead to an important indication of scalable Majorana vortex modes in FeTe

*Se*

_{x}_{1−x}. Thus, FeTe

*Se*

_{x}_{1−x}can be one promising platform having scalable Majorana qubits for quantum computing.

## INTRODUCTION

Majorana zero modes, being their own antiparticle in many-body systems, have attracted great interest for quantum computing (*1*, *2*). Two of the promising platforms to host these exotic modes are a one-dimensional (1D) superconducting proximitized semiconductor nanowire (*3*–*5*) and an *s*-wave superconducting surface with a Dirac cone (*6*). The observation of the quantized Majorana conductance (*7*) in the nanowire and the distinct zero-bias peaks (ZBPs) (*8*–*10*) in the vortex cores bring excitement and encouragement in this field. The Majorana zero modes may form qubits storing quantum information; the key property of the distant Majorana modes for quantum computing is the quantum information in the Majorana qubits can be protected by topology. On the other hand, according to the DiVincenzo criteria (*11*), achieving quantum computing requires a platform hosting scalable qubits. Seeking for scalable qubit systems is one of the major obstacles to pursuing quantum computing.

It has been observed that iron-based superconductors (SCs) (FeTe* _{x}*Se

_{1−x}; 0.55 ≤

*x*≤ 0.6) have multiple ZBPs of tunneling in the vortex cores (

*9*). Here, we incorporate these experimental data of the scanning tunneling microscopy (STM) and the simulation of the Majorana surface physics to show the evidence that an

*s*-wave superconducting surface with a single Dirac cone, as an effective 2D

*p*±

*ip*SC, can host multiple Majorana zero modes (

*6*), which can form scalable Majorana qubits. This iron-based SC can potentially be an ideal platform for Majorana quantum computing (

*2*) and interacting Majorana systems (

*12*).

Since in the presence of the vortex the SC order parameter phase winding aligns the spin rotation texture in the surface Dirac cone, a Majorana vortex mode (MVM) with zero energy (Majorana zero mode) emerges in the vortex core. The first candidate of this effective *p* ± *ip* SC is the heterostructure of a topological insulator (Bi_{2}Te_{3}) and an *s*-wave SC (NbSe_{2}). Recent STM measurement observed strong zero-bias conductance peaks at the vortex cores of the heterostructure (*13*, *14*), consistent with MVM interpretation in tunnel spectroscopy (*15*). However, the measurement resolution is not high enough to rule out other low-energy modes, which are the Caroli–de Gennes–Matricon (CdGM) modes (*16*). The two different modes residing in the vortex cores had been indistinguishable until the topology of the type II iron-based topological SC (FeTe* _{x}*Se

_{1−x}) was revealed (

*17*). This material has bulk Fermi surfaces and Dirac cone–type spin-helical surface states, for which SC pairing can be induced by the bulk

*s*-wave superconductivity (

*18*). The Fermi energies of the bulk and surface states are small enough to lead to the energy of the CdGM modes being much greater than the STM resolution. Hence, while the tunneling peaks away from zero-bias voltage have been observed and indicated the CdGM modes (

*19*), the tunneling observation of the ZBPs in the vortex cores (

*8*) has provided the first clue of MVM existence. The further hint of the MVM existence is supported by the recent observation of the tunneling conductance plateau that is near the quantized value (

*20*,

*21*).

Recently, the high-energy-resolution STM (*22*) probed hundreds of vortex cores on the (001) surface of FeTe* _{x}*Se

_{1−x}in the different magnetic fields. The experimental data show that the fraction

*R*

_{ZBP}of the vortex cores having the ZBPs of tunneling spectra monotonically decreases as the magnetic field increases (

*9*). Another experimental group (

*23*) independently has observed the similar decreasing behavior of the ZBP rate, although the rate

*R*

_{ZBP}is much lower than that of the previous group operating the higher-energy-solution STM. Although as illustrated in Fig. 1 the Majorana hybridization might be a straightforward explanation, this picture cannot directly lead to the monotonically decreasing rate

*R*

_{ZBP}of the observed ZBPs without considering other factors. The reason is that the energy splitting of the two Majorana hybridization exhibits an oscillation with the exponential decay as the intervortex distance increases (

*24*,

*25*). In this regard, the ZBP rate

*R*

_{ZBP}should have an oscillating behavior, which has not been observed in FeTe

*Se*

_{x}_{1−x}.

Instead of the focus on the physics of individual vortex, the main goal of this manuscript was to investigate whether in all the observed vortex cores the global features of the ZBPs are linked to the nature of multiple MVMs. We build a 3D Bogoliubov–de Gennes (BdG) tight-binding (TB) model capturing the effective low-energy physics (*E* < 0.2 meV) of over a hundred MVMs on the (001) surface of FeTe* _{x}*Se

_{1−x}; the TB model includes only one surface Dirac cone with an

*s*-wave superconducting pairing, which leads to the presence of MVMs. In this regard, we simplify the TB model in the second quantization form of

*Se*

_{x}_{1−x}(

*8*,

*18*). Moreover, we insert Abrikosov vortices by adding phase windings in the order parameter and vector potential (

*12*,

*26*) in the hopping terms (see the explicit form of the TB model in the Supplementary Materials). It is known that the self-consistent theory can determine the values of the SC parameters. Circumventing the subtle self-consistent electronic structure of multiple vortices, we directly adapt the values of the physical parameters from the experimental data (

*8*,

*18*,

*27*,

*28*) for the TB model.

Our key finding is that the Majorana hybridization and the disordered vortex distribution result in the decreasing ZBP rate *R*_{ZBP} with the increasing magnetic field. Moreover, by analyzing the experimental data (*9*), we find the statistics of the observed energy peaks agrees with the one from the simulation, including the physics of the Majorana hybridization. That is, the agreement between the analyzed experimental data and the Majorana simulation suggests the existence of the scalable MVMs. In addition, we study the interplay of the vortex location and the ZBPs in the vortices to show that the subtle relation between the energy spectra and the Majorana couplings is consistent with our new experimental observation. To confirm multiple MVMs in the iron-based SCs, we propose spin-selected STM for the measurement of the universal Majorana spin signature. With further experimental confirmation in the future, the surface of this iron-based SC can be a platform hosting scalable MVMs for quantum computing.

## RESULTS

### Majorana physics on the surface of FeTe_{x}Se_{1−x}

_{x}

We start with one single vortex through the iron-based SC FeTe* _{x}*Se

_{1−x}. The iron-based SC has a bulk-band inversion gap so that a Dirac cone protected by time-reversal symmetry is present on its surface. Furthermore, the recent angle-resolved photoemission spectroscopy (ARPES) measurement (

*18*) has shown that an

*s*-wave superconducting gap is induced in the surface Dirac cone. Hence, this superconducting Dirac cone surface state leads to the nontrivial topology (

*6*). Since our focus is on only the low-energy physics of FeTe

_{0.55}Se

_{0.45}, the effective physics of the topological SC can be described by a 2D BdG Hamiltonian of the surface Dirac cone with an

*s*-wave superconducting pairing

*v*

_{F}is the Fermi velocity, and Δ

_{0}is the superconducting gap. To reveal the topological superconductivity, we introduce a vortex at

**r**= 0 so that the order parameter has the additional winding phase Δ = Δ

_{0}e

^{iθ}.

Hence, with broken translation symmetry, since the momentum *k* is not a good quantum number, the surface Hamiltonian has to be rewritten in real space*29*, *30*)*J _{i}*(

*k*

_{F}

*r*) is the Bessel function of the first kind, the coherence length ξ =

*v*

_{F}/Δ

_{0}, and the Fermi momentum of the surface Dirac cone

*k*

_{F}= μ/

*v*

_{F}. The coherence length ξ represents the spatial size of the MVM, and the Fermi wavelength 1/

*k*

_{F}represents the oscillation length of the MVM due to the oscillation of the Bessel functions. We note that the Majorana coherence length is physically different from the superconducting coherence length, which represents the length scale of the vanishing superconducting gap in a vortex, although these two lengths are close in the Bardeen-Cooper-Schrieffer (BCS) theory.

Because these two key lengths characterize the MVM, experimentally determining these characteristic lengths is essential to confirm the existence of the MVMs. For FeTe* _{x}*Se

_{1−x}, we adapt SC gap Δ

_{0}= 1.8 meV, Fermi velocity

*v*

_{F}= 25 nm·meV, and Fermi momentum

*k*

_{F}= 0.2 nm

^{−1}based on the recent experimental measurements of ARPES (

*8*,

*18*); hence, the values of the two characteristic lengths are given by the coherence length ξ = 13.9 nm and the Fermi wavelength 1/

*k*

_{F}= 5 nm. Using these characteristic lengths, we build a 3D TB for the simulation, and fig. S1 (A and B) shows the wave function of the single MVM in the continuum model (

*3*) and the TB model are in perfect agreement.

### Hybridization of two MVMs

The hybridization of the two MVMs plays an important role to affect the presence of the ZBPs in the vortex cores. By considering two vortices with distance *R _{ij}*, we have two MVMs residing in the two vortex cores, respectively, as illustrated in Fig. 2A. The energy hybridization of the two MVMs is captured by a distance

*R*–dependent term and a phase term separately affected by all of the vortices (

_{ij}*t*∝

_{ij}*t*

_{R}t_{phase}). For

*k*

_{F}

*r*≫ 1, because of the asymptotic form of the Bessel function, the distance-dependent hybridization energy of the two MVMs in the continuum model is given by (

*24*,

*25*)

Although this hybridization approximation holds in the large *k*_{F}*r* limit, as shown in Fig. 2B, this is in good agreement with the TB model, which naturally includes the Majorana hybridization (see the Supplementary Materials for the details of the TB model). It is known that the phase factor in the continuum model is given by θ = π/4; however, in reality, θ depends on the details of the system. For example, for the TB model of the two vortices, θ = 1.40 in our simulation for FeTe_{0.55}Se_{0.45}. In the case of the iron-based SC, the Majorana hybridization exhibits the oscillation with the length period π/*k*_{F} = 15.7 nm. On the one hand, in the STM experiment (*9*), the intervortex distance is varied from 46.5 to 19.2 nm as the magnetic field is increased from 1 to 6 T. Since the oscillation length of the Majorana hybridization is within the experimental range, the ZBP rate *R*_{ZBP} controlled by the hybridization should be somehow oscillating as the magnetic field increases. However, the oscillation has never been observed. The main goal of this manuscript is to resolve this puzzle by considering additional factors in the Majorana physics.

The hybridization energy is also affected by the phases from all the vortices on the surface. We introduce the third vortex away from those two hybridized MVMs as an example in Fig. 2C; the idea can be easily extended to multiple vortices. While the third vortex contributes additional superconducting gap phases to the first two MVMs (*31*), *h*/2*e* magnetic flux with London penetration depth [λ ∼ 500 nm (*28*, *32*)] goes through each vortex, and the vector potential of the magnetic flux from each vortex affects the strength of the Majorana hybridization (*12*, *26*). Apart from the oscillation and exponential decay (*5*), the Majorana hybridization energy is proportional to an additional term (*12*, *26*)**A** is the vector potential stemming from the magnetic flux and *33*), when all of the intervortex distances are much less than the London penetration depth λ, the magnetic flux is diluted so that the vector potential **A** can be neglected (see Eq. A11 in the Supplemental Materials). Hence, ω_{12} = (π + Δ_{3})/2, where the π phase stems from the first two vortices and Δ_{3} is the phase difference at the two coupling vortices from the third vortex as shown in Fig. 2C. When the third vortex is moved away from the first two vortices so that the distances are much greater than the London penetration depth (*R*_{13}, *R*_{23} ≫ λ), Δ_{3} becomes small, and the magnetic flux from the third vortex screens out the contribution of ω_{12} from the third vortex

That is, the far-away vortex no longer affects the hybridization of the two Majoranas. This gives rise to the minimal simulation size of the TB model, which must be greater than the London penetration depth.

### Majorana vortex triangular lattice

We extend the system from few vortices to triangular vortex lattice on the surface of FeTe* _{x}*Se

_{1−x}. The vector potential

**A**from the penetrating magnetic field cannot be neglected, since the magnetic flux is accumulated in the presence of the multiple vortices in the wide spatial range. The simple continuum model (

*2*) of the single vortex is unable to capture the surface physics of the vortex lattice including the magnetic flux. Therefore, our TB model provides a straightforward way to reveal the Majorana lattice physics, and the mechanism of the TB model naturally includes the physics of the Majorana hybridization. There are two important characteristic lengths for the simulation, Majorana coherence length

*v*

_{F}/Δ= 13.9 nm and London penetration depth λ ∼ 500 nm. For the simulation system, not only must the surface size be much greater than the Majorana coherence length, but also the virtual surface size must be much greater than the London penetration depth (see Supplementary Appendix B for details). Namely, the virtual size means that virtual vortices, which are present outside the physical simulation system, with the magnetic flux and their SC phase windings added in the physical system. Since λ ≫

*v*

_{F}/Δ for FeTe

*Se*

_{x}_{1−x}, the real surface size for the simulation can be much smaller than the virtual surface size. By using this simulation setup, the MVMs in the central region of the system do not encounter the finite size effect, since the coupling between any two MVMs in the center and on the boundary is suppressed because of the exponential decay of the Majorana wave function (

*4*). Furthermore, in Eq. 7, the contribution of the phase factor ω

_{ij}from vortices outside the virtual system has been screened out.

The distributions of the triangular vortex lattice with an intervortex distance of 52.5 and 31.5 nm are shown in Fig. 3 (A and B). Each vortex core in the central region of the surface has the identical local density of states (LDOSs) with fixed *R*_{v} as an effective lattice translational symmetry due to the screening effect from the magnetic flux. In the absence of the magnetic flux, the translation symmetry in the central region is completely broken as shown in fig. S1, since the virtual surface size for the simulation is always not big enough to reach the screening boundary as the effective London penetration depth λ → ∞. Hence, introducing the magnetic flux plays an important role to depict the real experimental systems with finite London penetration depth.

Choosing the effective temperature (*T*) for the simulation is also essential to capture the observed LDOS of the vortex cores. Here, we assume that the tunneling spectra measured by STM are identical to the LDOS from the TB model simulation. In the experiment, the average full width at half maximum (FWHM) of the LDOS peak in the low field (*B* = 1 T, *R*_{v} = 46.5 nm) is around 80 μeV as shown in fig. S4. Hence, we choose *T* = 20 μeV so that, as shown in Fig. 3C (*R*_{v} = 52.5 nm), the FWHM of the ZBPs in the TB model is given by 76 μeV, which is reasonably close to the experimental value (see fig. S4) (*9*). Since each vortex core has an identical LDOS with the fixed *R*_{v}, Fig. 3D shows the LDOS as a function of the intervortex distance *R*_{v}. The ZBPs persist for *R*_{v} ≥ 48 nm since the temperature broadening erases the effect of the small energy splitting. The FWHM of the ZBPs is always close to 80 μeV for different *R*_{v}’s in this region.

For *R*_{v} < 48 nm, the oscillation of the peak location shares a similarity with the hybridization energy (*5*) of the two MVMs as a function of the intervortex distance as shown in Fig. 2C. When the MVMs get closer, the LDOS starts to break the particle-hole symmetry due to the Majorana hybridization and wave function overlapping (see Supplementary Appendix C for the detailed reason). Figure 3C (*R*_{v} = 31.5 nm) shows different peak heights with the opposite energies. On the other hand, another ZBP appears at *R*_{v} ≈ 36 nm.

The lattice translational symmetry leads to either 100 or 0% ZBP rate *R*_{ZBP} at the vortex cores. Therefore, the triangular vortex lattice cannot explain the monotonically decreasing fractional ZBP rate. This leads to further investigation of the vortex distribution on the surface.

### Disordered Majorana vortex distribution

In reality, because of the imperfection and impurity of the crystal, the distribution of the vortices in FeTe* _{x}*Se

_{1−x}is not a perfect triangular lattice anymore. The experimental observation (

*9*) shows in the Fourier space the vortex distributions resemble distorted triangular lattice in the low magnetic field and glass distribution in the high field. For the simulation, we set up these two types of vortex distributions in the TB model for different average intervortex distances (

*9*)). We examine whether on the top surface the LDOS in each vortex core within the central region (100 × 100 nm) has a ZBP. First, we examine the Majorana physics in the low field (

*B*= 1 T) corresponding to the average intervortex distance

*R*

_{ZBP}= 85.0% is fractional, while red circles and yellow diamonds indicate vortex cores with or without ZBPs in Fig. 4A. The LDOS of the selected vortex cores with or without ZBPs is shown in Fig. 4 (C and D, respectively). On the other hand, in the high field

*B*= 3.33 T, which corresponds to short intervortex distance (

*34*) in Fourier space, as shown in Fig. 4 (E and F) in the Fourier and real spaces, respectively. The ZBP rate

*R*

_{ZBP}= 14.1% becomes much smaller by comparing with the low magnetic field (Fig. 4A). Note that the LDOS of the selected vortex cores with or without ZBPs in Fig. 4 (G and H) has multiple peaks. Since the energies (≥0.96 meV) of the CdGM modes are much higher than the energy peaks (≤0.2 meV) of the LDOS, the multiple peaks stem from the hybridization of the multiple MVMs instead of CdGM modes (see Supplementary Appendix A for the CdGM energies.)

We consider the two different types of vortex distributions for the TB model simulation. Our key result presenting in Fig. 5A shows that the overall ZBP rate *R*_{ZBP} in the TB model qualitatively decreases when the average intervortex distance becomes shorter with the increasing magnetic field. We note that some of the physical parameters (e.g., Δ) may vary as the magnetic field is changed. Since the upper critical magnetic field is close to 50 T (*35*), we assume that the magnetic field strength (1 to 6 T) does not affect the physical parameters, except for the average intervortex distance. In the low magnetic field (* _{x}*Se

_{1−x}. We further examine whether inhomogeneous superconducting gaps or inhomogeneous chemical potentials can affect the LDOS of the vortex cores. The inhomogeneity of the superconducting gaps does not alter the LDOS in the vortex cores, whereas nonuniform chemical potentials can affect the LDOS (see fig. S12 for the details). However, because Se and Te are isovalent, the spatial inhomogeneity of the chemical potential should be very small (

*36*). Only the vortex distribution is the main physical factor affecting the LDOS of the vortex cores.

### Statistics of energy peaks

Examining the statistics of the first peak in each vortex core closest to zero energy in the experimental data (*9*), we find another important clue that supports the surface physics stemming from the hybridization of multiple MVMs. By excluding the vortex cores having ZBPs, most of the first peaks in the remaining vortex cores as shown in Fig. 5B (blue bars) are very close to zero energy. We also use the simulation result for the Majorana surface physics to plot the number distribution of the first peaks in Fig. 5C. The experimental data and simulation share a similar peak distribution showing the number of the peaks decreases as the energy moves away from zero energy. In this regard, these peaks slightly away from zero energy can be explained by the hybridization of the multiple MVMs. This peak statistics is another evidence supporting the existence of the multiple MVMs on the surface of FeTe* _{x}*Se

_{1−x}.

On the other hand, we examine the vortex cores with ZBPs as the second energy peaks in the experimental data (*9*). The distribution of the second peaks in those vortices closest to zero energy is always empty near zero energy with 60-μeV gap. The second peaks can arise from the CdGM mode with the lowest nonzero energy (*37*) or the hybridization of the MVMs. Figure 4G demonstrates that the second peak stems only from the Majorana hybridization in the TB model simulation since the energies of the CdGM modes from the surface superconductivity are greater than 0.96 meV. However, it is difficult to distinguish the two different origins by merely looking at the energy peaks in the experiment. Although it requires further examination, we suggest that the wide spreading of the second peaks in Fig. 5B can result from many factors, such as the coupling of the multiple CdGM modes nearby from the bulk superconductivity, the different strengths of the Majorana hybridization, and the inhomogeneous spatial distribution of the SC gap and chemical potential, since the energy of the CdGM mode is roughly proportional to Δ^{2}/μ.

The emergence of the second peaks closest to zero energy can further explain the observation of the low ZBP rates (*8*) when the vortex cores in FeTe* _{x}*Se

_{1−x}were probed by lower-energy-resolution STM regardless of the strength of the magnetic field. Consider one ZBP and a side peak closest to zero energy. We note that a CdGM mode can lead to the second peak and breaks particle-hole symmetry of LDOS (

*38*). The high-energy-resolution STM can distinguish the two peaks as illustrated in Fig. 5D, whereas in the low-energy-resolution STM, the two energy peaks merge to one peak away from zero energy as illustrated in Fig. 5E. In this regard, the ZBP cannot be observed by the low-energy-resolution STM when the energy of the second peak is roughly within the energy resolution. Thus, this is the reason that the low-energy-resolution STM has low ZBP rates.

### ZBPs and vortex locations

Since the relation between the presence of the ZBPs and the Majorana hybridization on the surface is subtle, it has been studied (*9*) that the presence of the ZBPs in the vortex cores does not exhibit simple correlations with spatial distributions of atoms, zero-magnetic-field in-gap bound states, and superconducting gaps. On the other hand, using time dependence measurements, we observe that the creeping vortex, which slightly moved on the surface of FeTe_{x}Se_{1−x}, has a significant change of the LDOS as shown in fig. S8. Furthermore, our additional STM measurement shows that at *B* = 4 T, in some cases, the LDOS of the vortex is not affected by other moving vortices when its surrounding vortex distribution is unchanged, as shown in fig. S9. To understand how the LDOS in the vortices is affected by the Majorana hybridization physics, we study how a single moving vortex alters the LDOS in each vortex core in the field of view (FOV). We simulate the TB model of the iron-based SC in the high magnetic field (4 T, *R*_{ZBP} from 47.6 to 31.7%, while the second movement (Fig. 6C) of the central vortex does not lead to the marked change of the ZBP rate. This significant change in the first movement demonstrates that the ZBP rate *R*_{ZBP} fluctuates in the wide range with the same average intervortex distance as shown in Fig. 5A.

The LDOS for each vortex core varies as the central vortex moves. The moving vortex, which initially has a ZBP (Fig. 6D, α), does not have a ZBP after the movements. Unexpectedly, this movement can also affect the LDOS of the vortex 120 nm away from the moving vortex (cf. the coherence length ξ = 13.9 nm), as Fig. 6D shows that the ZBP appears after the first movement. It is known that the coupling between the moving MVM and the distant MVM is exponentially suppressed; however, the unexpected LDOS change stems from that the coupling between the distant MVM (δ) and its nearby MVMs is slightly altered by the phase change (ω_{ij} in Eq. 6) from the central vortex (α) movement, since the distance between those vortices (α, ϵ) is less than the London penetration depth (~500 nm). The vortices near the moving vortex have significant changes of the LDOS as shown in Fig. 6D(α, β, γ, δ), since in this short length scale, the couplings of the moving MVM and its nearby MVMs are varied by the movement, and the coupling change between the other MVMs stems from the phase change (ω_{ij}) from the moving vortex. In the lower magnetic fields, the moving vortex affects the LDOS of its nearby vortices as shown in figs. S6 and S7. However, the LDOS of the vortices away from the moving vortex is not markedly altered, and there are slight changes in the ZBP rates.

Now, we are back to our new experimental data of the moving vortices in figs. S8 and S9. First, as shown in fig. S8, the LDOS of the creeping vortex exhibits a marked change after the creeping vortex moved 2 nm. This change can be explained by the change of the Majorana hybridization with the surrounding MVMs, because we have shown that the 7-nm vortex movement can significantly change the LDOS in Fig. 6D (β, γ). Furthermore, since the spatial range of the FOV probed by STM for the creeping vortex in fig. S8 is limited in a 15 nm by 15 nm square, it is possible that the surrounding vortices, which were not observed by STM, also moved at the same time. The location change of the surrounding vortices can affect the LDOS of the creeping vortex. Second, as shown in fig. S9, the LDOS of the observed vortex is unchanged after the movement of the vortices, which are away from the observed one. We have shown in Fig. 6D(ϵ, ζ) that the LDOSs of the vortices away from the moving vortex have only slight changes. Since in the experimental data, the smallest distance between the observed vortex and the moving vortices is around 50 nm, which is much longer than the Majorana coherence length, the LDOS of the observed vortex can be unchanged. Thus, the simulation of the vortex movement can qualitatively explain the observed LDOS in the vortices before and after the vortex movement.

### Spin LDOS of multiple MVMs

We have shown the multiple clues of the recent STM measurements in agreement with the results of the TB model hosting hundreds of MVMs on the surface. However, the observation of the tunneling peaks in the vortices is not enough to conclude the Majorana nature on the surface of FeTe* _{x}*Se

_{1−x}. It is crucial to predict and to observe the universal signature of the multiple MVMs even when the Majorana hybridization is strong and the ZBPs are removed.

In the literature (*29*), it has been proposed that the spin-polarized wave function of the isolated MVM oscillates spatially as the signature of the Majorana observation. On the one hand, practically the oscillation behavior might be difficult to be observed since the exponential decay of the Majorana wave function flattens the spin oscillation, so the current STM might not have enough sensitivity to detect the oscillation away from the vortex cores. On the other hand, since the Majorana coherence length ξ = 13.9 nm is comparable with the average intervortex distance *4*) of each individual MVM is kept in the original form. Thus, before experimentally looking for the spin signature of the multiple MVMs in the future, we first study the TB model of FeTe* _{x}*Se

_{1−x}to confirm that the spin oscillation is not altered by the Majorana hybridization and to check the requirement to observe the spin oscillation by the spin-selected STM probe.

The spin oscillation of the isolated MVM wave function (*4*) stems from the Bessel functions *J*_{0}(*k*_{F}*r*) and *J*_{1}(*k*_{F}*r*) for spin up and down, respectively, where *r* is the distance away from the vortex core. Therefore, because of the properties of the Bessel functions, the first peaks appear at *r* = 0 and *r* = 1.84/*k*_{F} for spin up and down, respectively, and the other oscillation peaks are exponentially suppressed. In the low magnetic field, the TB model shows similar spin spatial peaks for the vortices with and without ZBPs as shown in fig. S10. That is, for the isolated and hybridized MVMs, the LDOS of spin up exhibits peaks appearing at the vortex centers, and the LDOS of spin down vanishes at the vortex centers and has peaks away from the centers for *r* = 1.84/*k*_{F}. Although we focus on the spatial locations of the peaks, the energy peaks might be away from zero energy due to the Majorana hybridization. However, since, as shown in fig. S10, the peak height of spin down, which is suppressed by the exponential decay, is roughly six times smaller than the one of spin up, the spin down peak might be difficult to be detected by STM without high sensitivity. Alternatively, we integrate the LDOS within ±0.15-meV energy range, which is greater than the current STM energy resolution. Figure 7 shows the FOV for the two different spin distributions in the weak (1 T; Fig. 7, A and B) and strong (4 T; Fig. 7, C and D) magnetic fields. It is similar with the conventional STM probe that the spin up peaks in Fig. 7 (A and C) indicate the locations of the vortex centers. Differently, as shown in Fig. 7 (B and D), the FOV of spin down shows that each vortex core is encircled by a donut-shaped peak, since the LDOS spatial peaks of spin down appear away from the vortex cores. In the high magnetic field, because of the high density of the vortices, the donuts are connected to other neighbors. Although the MVMs strongly hybridize, preserving the MVM wave function leads to the observation of the donuts. This observation will be a smoking gun to show the FeTe* _{x}*Se

_{1−x}surface having the multiple MVMs.

## DISCUSSION

The decreasing fraction of vortices having ZBPs on the surface of FeTe* _{x}*Se

_{1−x}can be successfully explained by the hybridization of the MVMs with disordered vortex distributions. In other words, the recent observation of the decreasing ZBP rate supports the existence of multiple MVMs. While the oscillation of the Majorana coupling as a function of the intervortex distance has been erased by the distorted triangular lattice and glass of the vortex distributions, the Majorana coupling of the exponential decay as a function of the intervortex distance survives. Since the Majorana coupling increases with the increasing magnetic field, the stronger hybridization moves more MVMs away from zero energy and then reduces the ZBP rate. In other words, the reducing ZBP rate is the first evidence showing the hybridizations of multiple MVMs. Furthermore, by analyzing the energy peak statistics from the STM measurement, we find the numbers of the first LDOS peaks closest to zero energy, in the vortex cores not having ZBPs, have a pyramid-shaped distribution except for zero energy. This distribution, which is consistent with the TB simulation, can link to the Majorana hybridization moving the peaks away from zero energy.

The interplay between the LDOS and the detailed couplings of the MVMs is fascinatingly complicated. The slight movement (~14-nm coherence length) of the single vortex can significantly affect the LDOS of the vortex cores near the moving vortex, while in the high magnetic field, the vortex movement can alter the absence or presence of some ZBPs in vortex cores away from the moving vortex (distance less than the London penetration depth). It has been shown that manipulating the Majorana movement can manifest quantum logic gates (*39*, *40*) for quantum computing. Encouraged by the recent success controlling the movement of a single superconducting vortex (*41*), we believe that FeTe* _{x}*Se

_{1−x}can be an ideal platform to achieve Majorana quantum computing.

Since the spin textures of the MVMs survive in the presence of the Majorana hybridization, this universal spin feature of multiple hybridized MVMs probed by the spin-selected STM will be a smoking gun of scalable Majorana qubits in the iron-based SC platform. We predict that the multiple MVMs exhibit the spatial donut-shaped distribution of spin down surrounding the vortex cores. After the experimental confirmation, to have scalable Majorana qubits, we apply low magnetic field in FeTe* _{x}*Se

_{1−x}so that each MVM is well separated to stay at zero energy and that the quantum information stored in the Majorana qubits can be safely protected by topology. Having a platform of the scalable Majorana qubits fulfills one key of the DiVincenzo criteria for quantum computing. This platform can potentially become a primary step to build scalable quantum gates formed by MVMs and to read out Majorana qubits in the future.

## METHODS

The Supplementary Materials provide additional results from STM measurements and a detailed and systematic 3D TB model construction describing the physics of multiple MVMs on the FeTe* _{x}*Se

_{1−x}surface. On the one hand, the experimental methodology was already written in recent work (

*9*). On the other hand, as our numerical analysis is lengthy, here we point out the keys of the approach.

We distilled the main physical features of FeTe* _{x}*Se

_{1−x}to capture the Majorana nature by only considering the

*s*-wave pairing surface Dirac cone, since on the superconducting surface Dirac cone, the two characteristic lengths, the Majorana coherence length ξ =

*v*

_{F}/Δ

_{0}and the Fermi length 1/

*k*

_{F}, determine the main Majorana physics. First, by adopting the physical values of the SC gap Δ

_{0}= 1.8 meV, the Fermi velocity of the Dirac cone

*v*

_{F}= 25 nm·meV,

*k*

_{F}= 0.2 nm

^{−1}, and its chemical potential μ =

*v*

_{F}

*k*

_{F}= 5 meV, we build a TB model describing a single MVM on the (001) surface. When the simulation is extended to include hundreds of MVMs on the same surface, the magnetic fluxes through the vortices significantly affect the couplings of the MVMs. We introduce the vector potential in the TB model describing

*h*/2

*e*magnetic flux through each vortex with the London penetration depth λ = 500 nm. By performing the Lanczos algorithm, we found all of the energy states within the ±0.25-meV energy range. Using these energy states, we computed the LDOS in the vortex cores with effective temperature

*T*= 20 μeV, which leads to the FWHM close to the experimental observation. The criterion to determine ZBPs in the vortex cores is that the energy peak location of the LDOS is within the ±20-μeV energy range.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/9/eaay0443/DC1

Appendix A. 3D TB model of FeTe* _{x}*Se

_{1−x}.

Appendix B. Simulation details.

Appendix C. Particle-hole symmetry broken by Majorana hybridization.

Appendix D. Simulation for a moving vortex.

Appendix E. STM measurement of creeping vortex and vortex relocation.

Appendix F. Simulation for LDOS line profiles of ZBP and non-ZBP vortices.

Appendix G. Inhomogeneous superconducting gaps and chemical potentials.

Appendix H. ZBP heights.

Fig. S1. The line profile of the Majorana wave function.

Fig. S2. The system size of the 3D TB model for the simulation.

Fig. S3. The TB simulation without magnet flux.

Fig. S4. FWHM of the tunneling peaks in the STM experiment (*9*).

Fig. S5. The angles and phases for a point near vortex.

Fig. S6. The LDOS changes with *B* = 1 T before and after the central vortex α moves 14 nm down.

Fig. S7. The LDOS changes with *B* = 2 T before and after the central vortex α moves 14 nm down.

Fig. S8. Time dependence measurements on a creeping vortex in *B* = 4 T with averaged intervortex distance

Fig. S9. The tunneling spectra before and after the vortex relocation in the same FOV.

Fig. S10. (Spin) LDOS line profiles of the two selected vortices.

Fig. S11. Inhomogeneous superconducting gaps for triangular vortex lattice with intervortex distances.

Fig. S12. Inhomogeneous chemical potentials for triangular vortex lattice with intervortex distances.

Fig. S13. ZBP heights from the simulation and the STM experiment in different magnetic fields.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We thank S. Das Sarma, M. Franz, L. Kong, D. Liu, T. Liu, X. Liu, P. A. Lee, T. Posske, D. Wang, H.-H. Wen, R. Wiesendanger, and H. Zheng for the helpful discussions and comments.

**Funding:**C.-K.C. and F.-C.Z. are supported by the Strategic Priority Research Program of the Chinese Academy of Sciences (grant XDB28000000). T.M. and T.H. are supported by CREST project JPMJCR16F2 from Japan Science and Technology Agency. T.M. is also supported by a Grants-in-Aid for Scientific Research (KAKENHI) (no. 19H01843) by the Japan Society for the Promotion of Science (JSPS). F.-C.Z. is also supported by Beijing Municipal Science and Technology Commission (no. Z181100004218001), NSF of China (grant no.11674278), and National Basic Research Program of China (no. 2014CB921203).

**Author contributions:**C.-K.C. and Y.H. built and performed the TB model. T.M. and T.H. analyzed the experimental data. C.-K.C. wrote the manuscript. All authors discussed the results and contributed to the manuscript.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).