Scalable Majorana vortex modes in iron-based superconductors

TL;DR

This paper models Majorana vortex modes in FeTe$_x$Se$_{1-x}$ superconductors, explaining experimental observations and suggesting the material's potential for scalable quantum computing with Majorana qubits.

Contribution

The study constructs a realistic 3D tight-binding model that explains experimental ZBP fraction decrease and energy peak statistics, indicating scalable Majorana vortex modes in the superconductor.

Findings

Simulation matches experimental ZBP decrease with magnetic field.

Energy peak statistics agree with Majorana physics.

FeTe$_x$Se$_{1-x}$ is a promising platform for Majorana qubits.

Abstract

A vortex in an s-wave superconductor with a surface Dirac cone can trap a Majorana bound state with zero energy leading to a zero-bias peak (ZBP) of tunneling conductance. The iron-based superconductor FeTe$_x$Se$_{1-x}$ is one of the material candidates hosting these Majorana vortex modes. It has been observed by recent scanning tunneling spectroscopy measurement that the fraction of vortex cores possessing ZBPs decreases with increasing magnetic field on the surface of this iron-based superconductor. We construct a three-dimensional tight-binding model simulating the physics of over a hundred Majorana vortex modes in FeTe$_x$Se$_{1-x}$ with realistic physical parameters. Our simulation shows that the Majorana hybridization and disordered vortex distribution can explain the decreasing fraction of the ZBPs observed in the experiment. Furthermore, we find the statistics of the energy peaks off zero energy in our simulation with the Majorana physics in agreement with the analyzed peak statistics in the vortex cores from the experiment. This agreement and the explanation of the decreasing ZBP fraction lead to an important indication of scalable Majorana vortex modes in the iron-based superconductor. Thus, FeTe$_x$Se$_{1-x}$ can be one promising platform possessing scalable Majorana qubits for quantum computing. In addition, we further show the interplay of the ZBP presence and the vortex locations qualitatively agrees with our additional experimental observation and predict the universal spin signature of the hybridized multiple Majorana vortex modes.

Equations

$$\displaystyle\hat{H}_{\rm{BdG}}=\hat{H}_{\rm{TI}}+\sum_{\bm{r}}\Big{[}\Delta({\bm{r}})C_{{\bm{r}}\downarrow}^{\dagger}C_{{\bm{r}}\uparrow}^{\dagger}-\Delta({\bm{r}})C_{{\bm{r}}\uparrow}^{\dagger}C_{{\bm{r}}\downarrow}^{\dagger}+h.c.\Big{]}.$$

$$H_{\mathrm{BdG}}^{\mathrm{surf}}(\mathbf{k})=\left(\begin{array}{cccc}-\mu & {\nu }_F(k_x-i{k}_y)& 0& -\Delta \\ {\nu }_F(k_x+i{k}_y)& -\mu & \Delta & 0\\ 0& {\Delta }^{\ast }& \mu & {\nu }_F(k_x+i{k}_y)\\ -{\Delta }^{\ast }& 0& {\nu }_F(k_x-i{k}_y)& \mu \end{array}\right),$$

$$H_{\mathrm{BdG}}^{\mathrm{surf}}(\mathbf{r})=\left(\begin{array}{cccc}-\mu & -i{\nu }_F{e}^{-i\theta }({\partial }_r-\frac{i}{r}{\partial }_{\theta })& 0& -{\Delta }_0{e}^{i\theta }\\ -i{\nu }_F{e}^{i\theta }({\partial }_r+\frac{i}{r}{\partial }_{\theta })& -\mu & {\Delta }_0{e}^{i\theta }& 0\\ 0& {\Delta }_0{e}^{-i\theta }& \mu & -i{\nu }_F{e}^{-i\theta }({\partial }_r+\frac{i}{r}{\partial }_{\theta })\\ -{\Delta }_0{e}^{-i\theta }& 0& -i{\nu }_F{e}^{-i\theta }({\partial }_r-\frac{i}{r}{\partial }_{\theta })& \mu \end{array}\right),$$

$$|\mathrm{MVM}⟩=e^{-r/\xi }(J_0(k_{F}r),J_1(k_{F}r)e^{i\theta },J_0(k_{F}r),J_1(k_{F}r)e^{-i\theta }{)}^T,$$

$$t_{ij}\propto t_R=\frac{\cos (k_F{R}_{ij}+\theta )}{\sqrt{R_{ij}}}{e}^{-R_{ij}/\xi }.$$

$$t_{12}\propto t_{\rm{phase}}=\sin\omega_{12},\ {\rm{where}}\ \omega_{12}=\int_{\bm{r}_{1}}^{\bm{r}_{2}}\Big{(}\frac{1}{2}\nabla\theta-\frac{e}{\hbar}\bm{A}\Big{)}\cdot d\bm{l},$$

$$(\frac{1}{2}\nabla \theta -\frac{e}{\hslash }\mathbf{A}{)}_3\approx 0.$$

$$H_{\mathrm{TI}}=$$

$$+{\nu }_F\sin k_{z}a{\rho }_x{\sigma }_z-\mu {\rho }_0{\sigma }_0,$$

$$H_{\mathrm{BdG}}=$$

$$+{\nu }_F\sin k_{z}a{\tau }_0{\rho }_x{\sigma }_z-\mu {\tau }_z{\rho }_0{\sigma }_0+{\Delta }_0{\tau }_y{\rho }_0{\sigma }_y,$$

$${\hat{H}}_{\mathrm{BdG}}$$

$$+\sum _{\mathbf{r},\mathbf{\delta }}\left(\begin{array}{cc}{C}_{\mathbf{r}}^{†}& C_{\mathbf{r}}\end{array}\right)\left(\begin{array}{cc}{H}_{\mathrm{nn}}(\mathbf{\delta })& 0\\ 0& -H_{\mathrm{nn}}^{\ast }(\mathbf{\delta })\end{array}\right)\left(\begin{array}{c}{C}_{\mathbf{r}+\mathbf{\delta }}\\ C_{\mathbf{r}+\mathbf{\delta }}^{†}\end{array}\right),$$

$${\hat{H}}_{\mathrm{BdG}}$$

$$+\sum _{\mathbf{r},\mathbf{\delta }}\left(\begin{array}{cc}{C}_{\mathbf{r}}^{†}& C_{\mathbf{r}}\end{array}\right)\left(\begin{array}{cc}{H}_{\mathrm{nn}}(\mathbf{\delta })e^{-\frac{ie}{\hslash }{\int }_{\mathbf{r}}^{\mathbf{r}+\mathbf{\delta }}\mathbf{A}(\mathbf{r})\cdot d\mathbf{l}}& 0\\ 0& -H_{\mathrm{nn}}^{\ast }(\mathbf{\delta })e^{\frac{ie}{\hslash }{\int }_{\mathbf{r}}^{\mathbf{r}+\mathbf{\delta }}\mathbf{A}(\mathbf{r})\cdot d\mathbf{l}}\end{array}\right)\left(\begin{array}{c}{C}_{\mathbf{r}+\mathbf{\delta }}\\ C_{\mathbf{r}+\mathbf{\delta }}^{†}\end{array}\right),$$

$$\Delta (\mathbf{r})={\Delta }_0\tanh (r/{\xi }_0)e^{i\theta },$$

$$B_z(\mathbf{r})-{\lambda }^2{\nabla }^2{B}_z(\mathbf{r})={\Phi }_0\delta (\mathbf{r}),$$

$$B_z(r)=\frac{{\Phi }_0}{2\pi {\lambda }^2}{K}_0(r/\lambda ),$$

$$\phi(r)=2\pi\int_{0}^{r}B_{z}(r^{\prime})r^{\prime}dr^{\prime}=\big{(}1-\frac{r}{\lambda}K_{1}(r/\lambda)\big{)}\Phi_{0}.$$

$$\bm{A}({\bm{r}})=\Phi_{0}(\frac{y}{r^{2}}\hat{x}-\frac{x}{r^{2}}\hat{y})(1-\frac{r}{\lambda}K_{1}\big{(}\frac{r}{\lambda})\big{)}$$

$$\Delta (\mathbf{r})={\Delta }_0\prod _j\tanh (|\mathbf{r}-{\mathbf{r}}_j|/{\xi }_0)\frac{x-x_j+(y-y_j)i}{|\mathbf{r}-{\mathbf{r}}_j|}$$

$$\bm{A}({\bm{r}})=\Phi_{0}\sum_{j}\frac{(y-y_{j})\hat{x}+(x-x_{j})\hat{y}}{|{\bm{r}}-{\bm{r}}_{j}|}\Big{(}\frac{1}{|{\bm{r}}-{\bm{r}}_{j}|}-\frac{1}{\lambda}K_{1}(\frac{|{\bm{r}}-{\bm{r}}_{j}|}{\lambda})\Big{)}$$

$$\frac{e}{\hslash }{\int }_{\mathbf{r}}^{\mathbf{r}+\mathbf{\delta }}\mathbf{A}(\mathbf{r})\cdot d\mathbf{l}\approx \sum _j{\theta }_j(1-\frac{|\bar{\mathbf{r}}-{\mathbf{r}}_j|}{\lambda }{K}_1(\frac{|\bar{r}-r_j|}{\lambda })),$$

$$|{\varphi }_{\mathbf{r}}^j⟩=({\mathbf{u}}_{\mathbf{r}}^j{\mathbf{\nu }}_{\mathbf{r}}^j),$$

$$N({\bm{r}},E)=\sum_{E_{j}<0}\Big{(}|\bm{u}^{j}_{\bm{r}}|^{2}\delta(\beta(E-E_{j}))+|\bm{\nu}^{j}_{\bm{r}}|^{2}\delta(\beta(E+E_{j}))\Big{)},$$

$$|{\varphi }_1(\mathbf{r})⟩$$

$$|{\varphi }_2(\mathbf{r})⟩$$

$$N_i(\mathbf{r},E)=2|{\mathbf{u}}_i(\mathbf{r}){|}^2\delta (\beta E).$$

$$H_h=it{\gamma }_1{\gamma }_2.$$

$$(\gamma_{1}-i\gamma_{2}){|\rm{BCS}\rangle}=\big{(}\bar{\bm{u}}({\bm{r}})c^{\dagger}_{\bm{r}}+\bar{\bm{\nu}}({\bm{r}})c_{\bm{r}}\big{)}|\rm{BCS}\rangle,$$

$$\bar{N}(\mathbf{r},E)=$$

$$u_{i}({\bm{r}})=e^{-|{\bm{r}}-{\bm{r}}_{i}|/\xi+i\phi_{i}/2}\big{(}J_{0}(k_{F}|{\bm{r}}-{\bm{r}}_{i}|),e^{i\theta_{i}}J_{1}(k_{F}|{\bm{r}}-{\bm{r}}_{i}|)).$$

$$\bar{u}({\mathbf{r}}_1)=$$

Full text

Introduction

Majorana zero modes, being their own antiparticle in many-body systems, have attracted great attention for quantum computing Kitaev2001 ; RMP_braiding . Two of the promising platforms to host these exotic modes are a one-dimensional superconducting proximitized semiconductor nanowire Gil_Majorana_wire ; Sau_semiconductor_heterostructures ; Roman_SC_semi and an s-wave superconducting surface with a Dirac cone Fu2008Superconducting . The observation of the quantized Majorana conductance Zhang:2018aa in the nanowire and the distinct zero bias peaks (ZBP) Wang2018Evidence ; 2018arXiv181208995M ; PhysRevX.8.041056 in the vortex cores brings 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 the topology. On the other hand, according to the DiVincenzo criteria DiVincenzo:2000aa , achieving quantum computing requires a platform hosting scalable qubits. Seeking for scalable qubit systems is one of the major obstacles to pursue quantum computing.

It has been observed that the iron-based superconductors (${\rm FeTe}_{x}{\rm Se}_{1-x}$, $0.55\leq x\leq 0.6$) possess multiple ZBPs of tunneling in the vortex cores 2018arXiv181208995M . Here, we incorporate this experimental data of the scanning tunneling microscope (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\pm ip$ superconductor (SC), can host multiple Majorana zero modes Fu2008Superconducting , which can form scalable Majorana qubits. This iron-based superconductor can potentially be an ideal platform for Majorana quantum computing RMP_braiding and interacting Majorana systems Chiu_strong_inter_Majorana .

Since in the presence of the vortex the superconductor 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\pm ip$ superconductor is the heterostructure of a topological insulator (${\rm Bi}_{2}{\rm Te}_{3}$) and an s-wave superconductor (${\rm NbSe}_{2}$). Recent scanning tunneling spectroscopy measurement observed strong zero-bias conductance peaks at the vortex cores of the heterostructure Xu2015Experimental ; Sun2016Majorana , consistent with MVM interpretation in tunnel spectroscopy Law2009Majorana . However, the measurement resolution is not high enough to rule out other low-energy modes, Caroli-de Gennes-Matricon (CdGM) modes CAROLI1964307 . The two different modes residing in the vortex cores had been indistinguishable until revealing the topology of the type-II iron-based topological superconductor (${\rm FeTe}_{x}{\rm Se}_{1-x}$) 10.1093/nsr/nwy142 . This material possesses bulk Fermi surfaces and Dirac-cone-type spin-helical surface states, for which SC pairing can be induced by the bulk $s$-wave superconductivity Zhang2018Observation . 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 Chen:2018aa , the tunneling observation of the ZBPs in the vortex cores Wang2018Evidence 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 2019arXiv190406124Z ; Chen_2019 .

Recently, the high-energy resolution STM doi:10.1063/1.5049619 probed hundreds of vortex cores on the ($001$) surface of ${\rm FeTe}_{x}{\rm Se}_{1-x}$ in the different magnetic fields. The experimental data show that the fraction $R_{\rm{ZBP}}$ of vortex cores possessing the ZBPs of tunneling spectra monotonically decreases as the magnetic field increases 2018arXiv181208995M . Another experimental group 2019arXiv190901686C independently has observed the similar decreasing behavior of the ZBP rate, although the rate $R_{\rm{ZBP}}$ is much lower than the previous group operating the higher energy solution of the 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_{\rm{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 PhysRevLett.103.107001 ; :2010tb . In this regard, the ZBP rate $R_{\rm{ZBP}}$ should have an oscillating behavior, which has not been observed in ${\rm FeTe}_{x}{\rm Se}_{1-x}$.

Instead of the focus on the physics of individual vortex, the main goal of this manuscript is 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 tight-binding (TB) model capturing the effective low-energy physics ($E<0.2$meV) of over a hundred MVMs on the (001) surface of ${\rm FeTe}_{x}{\rm 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

[TABLE]

and we construct the Hamiltonian of the topological insulator $\hat{H}_{\rm{TI}}$ to possess a surface Dirac cone being identical to the one in ${\rm FeTe}_{x}{\rm Se}_{1-x}$ Zhang2018Observation ; Wang2018Evidence . Moreover, we insert Abrikosov vortices by adding phase windings in the order parameter and vector potential Liu_2015 ; Chiu_strong_inter_Majorana in the hopping terms (see the explicit form of the TB model in Supplementary Materials). It is known that the self-consistent theory can determine the values of the superconductor parameters. Circumventing the subtle self-consistent electronic structure of multiple vortices, we directly adapt the values of the physical parameters from the experimental data Zhang2018Observation ; Wang2018Evidence ; PhysRevB.95.224505 ; PhysRevB.81.180503 for the TB model.

Our key finding is that the Majorana hybridization and the disordered vortex distribution result in the decreasing ZBP rate $R_{\rm{ZBP}}$ with the increasing magnetic field. Moreover, by analyzing the experimental data 2018arXiv181208995M 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 superconductors, 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 superconductor can be a platform hosting scalable MVMs for quantum computing.

Results

Majorana Physics on the surface of ${\rm FeTe}_{x}{\rm Se}_{1-x}$

We start with one single vortex through the iron-based superconductor FeTexSe1-x. The iron-based superconductor possesses a bulk-band inversion gap so that a Dirac cone protected by time-reversal symmetry is present on its surface. Furthermore, the recent ARPES measurement Zhang2018Observation 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 non-trivial topology Fu2008Superconducting . Since our focus is on only the low-energy physics of FeTe0.55Se0.45, the effective physics of the topological superconductor can be described by a 2D Bogoliubov-de-Gennes (BdG) Hamiltonian of the surface Dirac cone with an s-wave superconducting pairing

[TABLE]

in the basis of $C_{\bm{k}}=(c_{\uparrow{\bm{k}}},c_{\downarrow{\bm{k}}},c_{\uparrow{\bm{k}}}^{\dagger},c_{\downarrow{\bm{k}}}^{\dagger})$, where $\mu$ is the chemical potential, $\nu_{F}$ the Fermi velocity, and $\Delta_{0}$ the superconducting gap. To reveal the topological superconductivity, we introduce a vortex at ${\bm{r}}=0$ so that the order parameter has the additional winding phase $\Delta=\Delta_{0}e^{i\theta}$.

Hence, with broken translation symmetry, since the momentum $\bm{k}$ is not a good quantum number, the surface Hamiltonian has to be rewritten in real space

[TABLE]

in the basis of $C_{\bm{r}}=(c_{\uparrow{\bm{r}}},c_{\downarrow{\bm{r}}},c_{\uparrow{\bm{r}}}^{\dagger},c_{\downarrow{\bm{r}}}^{\dagger})$. By solving the eigenvalue problem, the wavefunction of the MVM with zero energy trapped in the vortex is in the form of PhysRevLett.115.177001 ; PhysRevB.94.224501

[TABLE]

where $J_{i}(k_{F}r)$ is the Bessel function of the first kind, the coherence length $\xi=\nu_{F}/\Delta_{0}$, and the Fermi momentum of the surface Dirac cone $k_{F}=\mu/\nu_{F}$. The coherence length $\xi$ 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 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 ${\rm FeTe}_{x}{\rm Se}_{1-x}$, we adapt superconductor gap $\Delta_{0}=1.8$meV, Fermi velocity $\nu_{F}=25$nm$\cdot$meV, and Fermi momentum $k_{F}=0.2$nm*-1* based on the recent experimental measurements of angle-resolved photoemission spectroscopy Zhang2018Observation ; Wang2018Evidence ; hence, the values of the two characteristic lengths are given by the coherence length $\xi=13.9$nm and the Fermi wavelength $1/k_{F}=5$nm. Using these characteristic lengths, we build a three-dimensional TB for the simulation and Fig. S1(a,b) shows the wavefunction of the single MVM in the continuum model (3) and the TB model are in perfect agreement.

Hybridization of two Majorana vortex modes

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. 2(a). The energy hybridization of the two MVMs is captured by a distance $R_{ij}$ dependent term and a phase term affected by all of the vortices separately ($t_{ij}\propto t_{R}t_{\rm{phase}}$). For $k_{F}r\gg 1$, due to the asymptotic form of the Bessel function, the distance-dependent hybridization energy of the two MVMs in the continuum model is given by PhysRevLett.103.107001 ; :2010tb

[TABLE]

Although this hybridization approximation holds in the large $k_{F}r$ limit, as shown in Fig. 2(b) this is in good agreement with the TB model, which naturally includes the Majorana hybridization. (See Supplementary Materials for the details of the TB model.) It is known that the phase factor in the continuum model is given by $\theta=\pi/4$; however, in reality, $\theta$ depends on the details of the system. For example, for the TB model of the two vortices $\theta=1.40$ in our simulation for FeTe0.55Se0.45. In the case of the iron-based superconductor, the Majorana hybridization exhibits the oscillation with the length period $\pi/k_{F}=15.7$nm. On the one hand, in the STM experiment 2018arXiv181208995M , the intervortex distance is varied from $46.5$nm to $19.2$nm as the magnetic field is increased from $1$T to $6$T. Since the oscillation length of the Majorana hybridization is within the experimental range, the ZBP rate $R_{\rm{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 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. 2(d); the idea can be easily extended to multiple vortices. While the third vortex contributes additional superconducting gap phases to the first two MVMs PhysRevLett.111.136401 , $h/2e$ magnetic flux with London penetration depth ($\lambda\sim 500$nm PhysRevB.81.180503 ; PhysRevB.82.184506 ) goes through each vortex and the vector potential of the magnetic flux from each vortex affects the strength of the Majorana hybridization Chiu_strong_inter_Majorana ; Liu_2015 . Apart from the oscillation and exponential decay (5), the Majorana hybridization energy is proportional to an additional term Liu_2015 ; Chiu_strong_inter_Majorana

[TABLE]

where $\bm{A}$ is the vector potential stemming from the magnetic flux and $\int_{\bm{r}_{i}}^{\bm{r}_{j}}\nabla\theta\cdot d\bm{l}=\theta_{j}-\theta_{i}$ is the phase difference between the superconducting gap at the first two hybridized MVMs. This coupling relation is directly included in our TB model with the magnetic flux. For few vortices residing (three in Fig. 2(d)) in a small superconductor SERRIERGARCIA2017109 when all of the intervortex distances are much less than the London penetration depth $\lambda$, the magnetic flux is diluted so that the vector potential $\bm{A}$ can be neglected (see Eq. 18). Hence, $\omega_{12}=(\pi+\Delta_{3})/2$, where $\pi$-phase stems from the first two vortices and $\Delta_{3}$ is the phase difference at the two coupling vortices from the third vortex as shown in Fig.  2(d). 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}\gg\lambda$), $\Delta_{3}$ becomes small and the magnetic flux from the third vortex screens out the contribution of $\omega_{12}$ from the third vortex

[TABLE]

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 FeTexSe1-x. The vector potential $\bm{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 $\Delta/\nu=13.9$nm and London penetration depth $\lambda\sim 500$nm. For the simulation system, not only the surface size must be much greater than the Majorana coherence length but also the virtual surface size must be much greater than the London penetration depth. (See sec. B in Supplementary Materials for the 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 are added in the physical system. Since $\lambda\gg\Delta/\nu$ for FeTexSe1-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 due to the exponential decay of the Majorana wavefunction (4). Furthermore, by Eq. 7 the contribution of the phase factor $\omega_{ij}$ from vortices outside the virtual system has been screened out.

The distributions of the triangular vortex lattice with intervortex distance $52.5$nm and $31.5$nm are shown in Fig. 3(a,b). Each vortex core in the central region of the surface possesses the identical local density of states (LDOS) 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. S3 since the virtual surface size for the simulation always is not big enough to reach the screening boundary as the effective London penetration depth $\lambda\rightarrow\infty$. Hence, introducing the magnetic flux plays an important role to depict the real experimental systems with finite London penetration depth.

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

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

The lattice translational symmetry leads to either $100\%$ or $0\%$ ZBP rate $R_{\rm{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, due to the imperfection and impurity of the crystal, the distribution of the vortices in ${\rm FeTe}_{x}{\rm Se}_{1-x}$ is not perfect triangular lattice anymore. The experimental observation 2018arXiv181208995M 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 ($18\leq\bar{R}_{v}\leq 48$nm cf. the experimental range $19\sim 46.5$nm 2018arXiv181208995M ). We examine if on the top surface LDOS in each vortex core within the central region ($100\times 100$nm) possesses a ZBP. First, we examine the Majorana physics in the low field ($B=1T$) corresponding to the average intervortex distance $\bar{R}_{v}=46.5$nm in the experiment. Fig. 4(a,b) show the vortex distribution of the distorted triangular lattice in real and Fourier spaces respectively. Due to the translational symmetry breaking, the ZBP rate $R_{\rm{ZBP}}=85.0\%$ is fractional while red circles and yellow diamonds indicate vortex cores with/without ZBPs in Fig. 4(a). The LDOS of the selected vortex cores with/without ZBPs is shown in Fig. 4(c,d) respectively. On the other hand, in the high field $B=3.33T$, which corresponds to short intervortex distance ($\bar{R}_{v}=25.5$nm), the vortices distribute like glass, which has a ring-shaped distribution RAMAN:1927aa in Fourier space, as shown in Fig. 4(e,f) in the Fourier and real spaces respectively. The ZBP rate $R_{\rm{ZBP}}=14.1\%$ becomes much smaller by comparing with the low magnetic field (a). It is worth to notice that the LDOS of the selected vortex cores with/without ZBPs in Fig. 4(g,h) possess multiple peaks. Since the energies ($\geq 0.96$meV) of the CdGM modes are much higher than the energy peaks ($\leq 0.2$meV) of the LDOS, the multiple peaks stem from the hybridization of the multiple MVMs, instead of CdGM modes. (See sec. A in Supplementary Materials for the CdGM energies.)

We consider the two different types of vortex distributions for the TB model simulation. Our key result presenting in Fig. 5(a) shows overall the ZBP rate $R_{\rm{ZBP}}$ in the TB model qualitatively decreases when the intervortex distance becomes shorter with the increasing magnetic field. We note that some of the physical parameters (e.g. $\Delta$) may vary as the magnetic field is changed. Since the upper critical magnetic field is close to $50$T PhysRevB.81.184511 , we assume that the magnetic field strength ($1\sim 6$T) does not affect the physical parameters, except for the average intervortex distance. In the low magnetic field ($\bar{R}_{v}\geq 33$nm) the ZBP rates of the two different vortex distributions share the similar decay consistent with the experimental observations, while the vortex distribution transits from the distorted lattice to glass distribution in the high magnetic field region. Since the ZBP statistics from the simulation does not exhibit any oscillation as a function of the intervortex distance, losing the perfect triangular lattice distribution of the vortices can be a primary factor smearing the oscillation feature of the Majorana coupling. The stronger the average strength of the Majorana coupling is, the shorter the average intervortex distance is. The stronger Majorana coupling increases the chance of the MVMs away from zero energy. Therefore, the disordered vortex distribution and the hybridization of the MVMs are the keys leading to the monotonic decreasing ZBP rate. In other words, the observation of the decreasing ZBP rate provides a clue to the multiple Majorana hybridizations on the surface of ${\rm FeTe}_{x}{\rm Se}_{1-x}$. We further examine if 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 non-uniform chemical potentials can affect the LDOS (See sec. S12 in Supplementary Materials for the details). However, because Se and Te are isovalent, the spatial inhomogeneity of the chemical potential should be very small PhysRevB.98.144519 . 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 2018arXiv181208995M , we find another important clue that supports the surface physics stemming from the hybridization of multiple MVMs. By excluding the vortex cores possessing ZBPs, most of the first peaks in the remaining vortex cores as shown in fig. 5(b)(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. 5(c). 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 ${\rm FeTe}_{x}{\rm Se}_{1-x}$.

On the other hand, we examine the vortex cores with ZBPs as the second energy peaks in the experimental data 2018arXiv181208995M . The distribution of the second peaks in those vortices closest to zero energy is always empty near zero energy with $60\mu$eV gap. The second peaks can arise from the CdGM mode with the lowest non-zero energy 2019arXiv190102293K or the hybridization of the MVMs. Fig. 4(g) 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. Although it requires further examination, we suggest that the wide spreading of the second peaks in Fig. 5(b) can result from many factors — 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 superconductor gap and chemical potential since the energy of the CdGM mode is roughly proportional to $\Delta^{2}/\mu$.

The emergence of the second peaks closest to zero energy can further explain the observation of the low ZBP rates Wang2018Evidence when the vortex cores in ${\rm FeTe}_{x}{\rm Se}_{1-x}$ were probed by lower-energy-resolution STM regardless the strength of the magnetic field. Consider one ZBP and another peak closest to zero energy. We note that a CdGM mode can lead to the second peak and breaks particle-hole symmetry of LDOS PhysRevLett.80.2921 . The high-energy-resolution STM can clearly distinguish the two peaks as illustrated in Fig. 5(d), whereas in the low-energy-resolution STM the two energy peaks merge to one peak away from zero energy as illustrated in Fig. 5(e). 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 2018arXiv181208995M 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 ${\rm FeTe}_{x}{\rm Se}_{1-x}$, has a significant change of the LDOS as shown in Fig. S8. Furthermore, our additional STM measurement shows that at $B=4T$ 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 superconductor in the high magnetic field (4T, $\bar{R}_{v}=24$nm). We start with the vortex distribution in Fig. 6(a) and the vortex in the center of the FOV moves down $7$nm twice. The first movement (Fig. 6(b)) dramatically changes the distribution of the ZBPs in the vortex cores and the ZBP rate $R_{\rm{ZBP}}$ from $47.6\%$ to $31.7\%$, while the second movement (Fig. 6(c)) of the central vortex does not lead to the dramatic change of the ZBP rate. This significant change in the first movement demonstrates that the ZBP rate $R_{\rm{ZBP}}$ fluctuates in the wide range with the same average intervortex distance as shown in Fig. 5(a).

LDOS for each vortex core varies as the central vortex moves. The moving vortex, which initially possesses a ZBP (Fig. 6(d)$\alpha$), does not have a ZBP after the movements. Surprisingly, this movement can also affect the LDOS of the vortex $120$nm away from the moving vortex (cf. the coherence length $\xi=13.9$nm) as Fig. 6(d)$\epsilon$ shows 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 surprising LDOS change stems from that the coupling between the distant MVM ($\delta$) and its nearby MVMs is slightly altered by the phase change ($\omega_{ij}$ in Eq.6) from the central vortex ($\alpha$) movement, since the distance between those vortices ($\alpha,\ \epsilon$) is less than the London penetration depth ($\sim 500$nm). The vortices near the moving vortex have significant changes of the LDOS as shown in Fig. 6(d)($\alpha,\beta,\gamma,\delta$), since in this short length scale the coupling 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 ($\omega_{ij}$) from the moving vortex. In the lower magnetic fields, the moving vortex affects the LDOS of its nearby vortices as shown in Fig. S6,S7. However, the LDOS of the vortices away from the moving vortex is not altered dramatically and there are slight changes in the ZBP rates.

Now we are back to our new experimental data of the moving vortices in Fig. S8,S9. First, as shown in Fig. S8 the LDOS of the creeping vortex exhibits a dramatic 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. 6(d)($\beta,\gamma$). Furthermore, since the spatial range of the FOV probed by STM for the creeping vortex in Fig. S8 is limited in $15$nm$\times 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. 6(d)($\epsilon,\zeta$) that the LDOS 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 Majorana vortex modes

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 ${\rm FeTe}_{x}{\rm 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 PhysRevLett.115.177001 it has been proposed that the spin-polarized wavefunction 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 wavefunction 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 $\xi=13.9$nm is comparable with the average intervortex distance $\bar{R}_{v}$ from $19$nm to $46.5$nm, in the presence of the Majorana hybridization it is uncertain that the wavefunction (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 ${\rm FeTe}_{x}{\rm Se}_{1-x}$ to confirm the spin oscillation is not altered by the Majorana hybridization and check the requirement to observe the spin oscillation by the spin-selected STM probe.

The spin oscillation of the isolated MVM wavefunction (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, due to the properties of the Bessel functions, the first peaks appear at $r=0$ and $r=1.84/k_{F}$ for spin up and down 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 as well as 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 6 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 $\pm 0.15$meV energy range, which is greater than the current STM energy resolution. Fig. 7 shows the FOV for the two different spin distributions in the weak ($1$T:(a,b)) and strong ($4$T:(c,d)) magnetic fields. It is similar with the conventional STM probe that the spin up peaks in (a,c) indicate the locations of the vortex centers. Differently, as shown in Fig. 7(b,d) the FOV of spin down shows 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, due to the high density of the vortices, the donuts are connected to other neighbors. Although the MVMs strongly hybridize, preserving the MVM wavefunction leads to the observation of the donuts. This observation will be a smoking gun to show the ${\rm FeTe}_{x}{\rm Se}_{1-x}$ surface possessing the multiple MVMs.

**Discussion

**The decreasing fraction of vortices possessing ZBPs on the surface of ${\rm FeTe}_{x}{\rm 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 possessing ZBPs, has 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 ($\sim 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 2019arXiv190106083L ; Alicea_Majorana_wire for quantum computing. Encouraged by the recent success controlling the movement of a single superconducting vortex Ge:2016aa , we believe that ${\rm FeTe}_{x}{\rm 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 superconductor 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 ${\rm FeTe}_{x}{\rm Se}_{1-x}$ so that each MVM is well separated to stay at zero energy and 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 provides additional results from STM measurements as well as a detailed and systematic three-dimensional TB model construction describing the physics of multiple MVMs on the ${\rm FeTe}_{x}{\rm Se}_{1-x}$ surface. On the one hand, the experimental methodology was already written in the recent work 2018arXiv181208995M . On the other hand, as our numerical analysis is lengthy, here we point out the keys of the approach.

We distill the main physical features of ${\rm FeTe}_{x}{\rm 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 $\xi=\nu_{F}/\Delta_{0}$ and the Fermi length $1/k_{F}$, determine the main Majorana physics. First, by adopting the physical values of the SC gap $\Delta_{0}=1.8$meV, the Fermi velocity of the Dirac cone $\nu_{F}=25$nm$\cdot$meV, $k_{F}=0.2$nm*-1*, its chemical potential $\mu=\nu_{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/2e$ magnetic flux through each vortex with the London penetration depth $\lambda=500$nm. By performing the Lanczos algorithm, we find all of the energy states within $\pm 0.25$meV energy range. Using these energy states, we compute the LDOS in the vortex cores with effective temperature $T=20\mu$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 $\pm 20\mu$eV energy range.

Acknowledgments: The authors 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) (numbers 19H01843) by the Japan Society for the Promotion of Science (JSPS). F.-C.Z. is also supported by Beijing Municipal Science $\&$ Technology Commission (No. Z181100004218001), National Science Foundation of China (Grant No.11674278), National Basic Research Program of China (No.2014CB921203). Author contributions: C.-K.C and Y.H. built and performed the tight-binding 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.

**Supplementary materials

**Appendix A: 3D tight-binding model of FeTexSe1-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

Appendix A 3D tight-binding model of ${\rm FeTe}_{x}{\rm Se}_{1-x}$

Many physical mechanisms are present on the surface of ${\rm FeTe}_{x}{\rm Se}_{1-x}$ probed by the recent experiments Zhang2018Observation ; Wang2018Evidence ; 2018arXiv181208995M . Since the focus is on the quantum states with energy within $\pm 0.2$meV, only Majorana physics emerges in this energy range. The reason is that the superconductor gap of the bulk and surface is much greater than $0.2$meV Wang2018Evidence and the low energy modes (CdGM modes) located in the vortex cores possess energy higher than $0.2$meV 2019arXiv190102293K . Hence, the low-energy physics already excludes the CdGM modes and the quasiparticles with energy higher than the superconductor gaps. Only MVMs emerge in the low-energy physics. The non-trivial topological superconductivity of ${\rm FeTe}_{x}{\rm Se}_{1-x}$, which hosts MVMs, stems from the surface Dirac cone. Eq. 3 in the manuscript shows that the phase winding of the vortex aligns the spin rotation texture of the surface Dirac cone so that the vortex can host an MVM. Although the $d-$wave orbitals are involved in the basis of the surface Dirac cone 10.1093/nsr/nwy142 , the emerging of the MVMs results from the spin texture of the Dirac cone, regardless of the choice of the orbitals. To capture the salient physics of multiple MVMs here, we start with the model of a non-trivial 3D topological insulator (TI) hosting a surface Dirac cone Zhang:2009aa , which is identical to the one in ${\rm FeTe}_{x}{\rm Se}_{1-x}$. The TI Hamiltonian in momentum space is written as

[TABLE]

where $M({\bm{k}})=m-\cos k_{x}a-\cos k_{y}a-\cos k_{z}a$, $\nu_{F}=25$nm$\cdot$meV, $k_{F}=0.2$nm*-1*, and chemical potential $\mu=\nu_{F}/k_{F}=5$meV based on the experimental measurement Zhang2018Observation ; Wang2018Evidence . The Pauli matrices $\sigma_{\alpha}$ and $\rho_{\beta}$ represent spin and orbital degrees of freedom respectively. Although for ${\rm FeTe}_{x}{\rm Se}_{1-x}$ its space group is $P4/nmm$, and the lattice constants have been known $a=b=3.8\AA$ $c=6.1\AA$ PhysRevB.79.054503 , we use cubic lattice and choose the lattice constant $a=2$nm in all of the directions to simplify the problem and focus on the physics of the surface Dirac cone, since the value of the lattice constant does not affect the low-energy Majorana physics.

By introducing s-wave superconductivity, the TI Hamiltonian is extended to the BdG Hamiltonian

[TABLE]

where the value of the s-wave surface superconducting gap ($\Delta_{0}\approx 1.8$meV) has been measured in the experiments Zhang2018Observation ; Wang2018Evidence ; 2018arXiv181208995M ; Chen:2018aa . Now we have the estimated values of the two Majorana characteristic lengths — Majorana coherence length $\xi=\nu_{F}/\Delta=13.9$nm and the Fermi wavelength $1/k_{F}=\nu_{F}/\mu=5$nm. Hence, the Fermi velocity $\nu_{F}$ of the surface Dirac cone is an essential quantity to determine the two characteristic lengths of the MVM. The primary platform we consider is the $(001)$ surface with vortices without translation symmetries; in this regard, the Bogoliubov-de Gennes (BdG) Hamiltonian can be rewritten in real space

[TABLE]

where the on-site part $H_{0}(m)=\nu_{F}m\rho_{z}\sigma_{0}-\mu\rho_{o}\sigma_{0}$, $\bm{\delta}=\pm a\hat{x},\pm a\hat{y},\pm a\hat{z}$ indicates the three directions of the hopping, and the nearest neighbor hopping part $H_{\rm{nn}}(\pm a\hat{n})=-\nu_{F}\rho_{z}\sigma_{0}/2\pm i\nu_{F}\rho_{x}\sigma_{n}/2$ for $n=x,y,z$. To reduce the number of the layers in the $z$ direction for the numerical simulation, we specifically choose $m=2$ so that the surface Dirac cones with small momentum can be almost located on the top and bottom layers. On these layers, the neutral point of the Dirac cone is located at $(0,0)$ in the $(001)$ surface BZ, which is consistent with ${\rm FeTe}_{x}{\rm Se}_{1-x}$.

Since ${\rm FeTe}_{x}{\rm Se}_{1-x}$ is a type-II superconductor, multiple Abrikosov vortices can be present with magnetic flux through the superconductor. In the presence of the vortices, the BdG Hamiltonian is in the form of

[TABLE]

where the superconductor gap $\Delta({\bm{r}})$ has an additional phase stemming from the vortices and the nearest neighbor hopping is modified by the Peierls substitution integral with the vector potential $\bm{A}({\bm{r}})$ describing the magnetic flux. In order to perform the numerical simulation, it is necessary to know the explicit expressions of $\Delta({\bm{r}})$ and $\bm{A}({\bm{r}})$.

First, we consider a single vortex located at ${\bm{r}}=(x,y)=(0,0)$ along the $z$ direction and then the superconducting gap is given by

[TABLE]

where $\xi_{0}$ is the superconducting coherence length. According to the experimental results (refer to Fig. 2(g,h) in 2018arXiv181208995M ), the superconducting gap appears roughly $5$nm away from the vortex cores, so $\xi_{0}=4.6$nm is chosen for the simulation. We note that the Majorana coherence length $\xi=13.9$nm, which is the length scale of the MVM, and the superconducting coherence length $\xi_{0}$ represent the different physical entities. To obtain the vector potential $\bm{A}({\bm{r}})$, we start with the magnetic field in the $z$ direction, which concentrates in the vortex core and can be described by the London equation deGennesbook

[TABLE]

where $\Phi_{0}=h/2e$ indicates the $\pi$-flux of the total magnetic flux on the entire surface and $\lambda$ is the London penetration depth. For ${\rm FeTe}_{x}{\rm Se}_{1-x}$, the London penetration length $\lambda\approx 500$nm PhysRevB.95.224505 ; PhysRevB.81.180503 . Since the London penetration length is much greater than the two characteristic lengths of the MVM ($\xi=13.9$nm and $1/k_{F}=5$nm), the magnetic field is diluted in the characteristic length scale. It is reasonable to neglect the magnetic flux effect for a single MVM in the continuum model in the manuscript. However, the magnetic flux plays a pivotal role in the physics of multiple MVMs; we still include the magnetic flux effect in the TB model with a single vortex.

By solving the London equation, the magnetic field in the $z$ direction is given by

[TABLE]

where $K_{i}(r/\lambda)$ is the modified Bessel function of the second kind. Furthermore, the total magnetic flux through the disk with the center at $(0,0)$ and the radius $r$ has the expression of

[TABLE]

As $r\rightarrow\infty$, the asymptotic behavior of the modified Bessel function leads to $\phi\rightarrow\Phi_{0}$ indicating the $\pi$-flux forming the single Abrikosov vortex. By using the relation $\phi(r)=\oint\bm{A}(\bm{r})\cdot d\bm{l}$ in the closed integral path with the fixed radius $r$, the vector potential can be explicitly written as

[TABLE]

Now all of the parameters and functions in the BdG Hamiltonian (11) with a single vortex are given; the wavefunction of the MVM trapped in the vortex on the surface of ${\rm FeTe}_{x}{\rm Se}_{1-x}$ can be solved in this TB model. Fig. S1(a,b) shows the TB model (11) and the continuum model (3) are in great agreement although the superconducting gap in the continuum model, which does not vanish near the vortex core (no $\tanh(r/\xi_{0})$ dependent), is slightly different from the superconducting gap in the tight binding model. We note that the MZM wavefunctions are almost identical in the presence and absence of the magnetic flux since the spatial range $r$ is much less than the London penetration depth Chiu:2011fk . Furthermore, the first three lowest energies of the CdGM modes trapped in the vortex in the TB are given by $\pm 0.96,\ \pm 1.30,\ \pm 1.43$meV as shown in Fig. S1(c). The energy levels of the CdGM modes are in order of $\Delta^{2}/\mu$, it is not necessary to have equal energy spacing, when the CdGM energy is close to the superconductor gap.

Since the surface of the type-II superconductor practically hosts more than one Abrikosov vortex, the physics of a single vortex is not enough to capture the mechanism of the topological superconductivity on the entire surface. Therefore, we consider multiple Abrikosov vortices located at ${\bm{r}}_{j}=(x_{j},y_{j})$ and then the superconducting gap is related to the locations of the vortices

[TABLE]

The vector potential for the multiple vortices is given by

[TABLE]

Thus, the complete BdG Hamiltonian of the TB model for multiple vortices, which are given, can be solved in principle.

Appendix B Simulation details

The major problem for the simulation is the choice of the system size ($L_{x},L_{y},L_{z}$) in the unit of the lattice constant. We expect to have $L_{z}$ as small as possible without losing the Majorana physics so that the size of the (001) surface can be extended to include more vortices. The two MVMs connected through the vortex along the $z$ direction can hybridize. This hybridization is expected to be strongly suppressed so the ZBP for an MVM can be seen on the surface. Here, we choose $m=2$ in the TB model (11) so that the surface Dirac cone has a rapidly exponential decay along the $z$ direction PhysRevB.86.155146 . Fig. S2(a) shows the energy hybridization $\Delta E$ of the only two MVMs on the top and bottom layers decays exponentially as a function of $L_{z}$. Since for $L_{z}=3$ the energy hybridization $\Delta E=3\mu$eV is within $\pm 20\mu$eV energy range of the ZBP criterion, we safely choose $L_{z}=3$ for the simulation.

Choosing the size ($L_{x},L_{y}$) of the $(001)$ surface is subtle. Ideally, the lengths of the surface edges should be much greater than the Majorana coherence length ($\xi=13.9$nm) and the London penetration depth ($\lambda=500$nm). That is, the surface size for the simulation has to be at least roughly greater than $1000\times 1000$ in the unit of the lattice constant ($a=2$nm). It will cost tremendous computer resources to solve the TB model in this size. Instead, we judiciously circumvent this scaling problem. We consider two types of surface sizes for the simulation — real size and virtual size. The real size indicates the range including almost all of the MVM couplings; the couplings outside this range, which exponentially decay, can be neglected. In other words, the range of the summations in the BdG Hamiltonian (11) is given by $|r^{x}|<L_{x}^{r}/2$ and $|r^{y}|<L_{y}^{r}/2$, where $(r^{x},r^{y})={\bm{r}}$. Hence, as shown in Fig. S2(b) we choose the real size $L_{x}^{r}\times L_{y}^{r}=464\times 480$nm, which is much greater than the Majorana coherence length ($13.9$nm). On the other hand, the virtual size indicates the range including almost all of the phase effect from the vortices. That is, the coupling of the two nearby MVMs is affected by phases from far away vortices. When the distance between the coupled MVMs and the far-away vortices is much greater than the London penetration depth, the phase effect is screened out by Eq. 7. Therefore, the virtual size must be greater than the real size and the London penetration depth so that the phase effect for the Majorana coupling in the real size range is included in the virtual size range. As shown in Fig. S2(b), we choose virtual size $L_{x}^{v}\times L_{y}^{v}=3712\times 3840$nm. In other words, the range of the summation in the order parameter (17) and the vector potential (18) are given by $|r^{x}_{j}|<L_{x}^{v}/2$ and $|r^{y}_{j}|<L_{y}^{v}/2$, where $(r^{x}_{j},r^{y}_{j})={\bm{r}}_{j}$.

Another time-consuming process for the simulation is to compute the integral $\frac{e}{\hbar}\int_{{\bm{r}}}^{{\bm{r}}+\bm{\delta}}\bm{A}({\bm{r}})\cdot d\bm{l}$ for each hopping term in the BdG Hamiltonian (11) because without the vector potential the TB simulation is incorrect as shown in Fig. S3. Since the distance $\delta$ of the two nearest neighbor atoms is small, we use this approximation

[TABLE]

where $\bar{{\bm{r}}}={\bm{r}}+\delta/2$, $\theta_{j}$ is the angle between line ${\bm{r}}_{j},{\bm{r}}$ and line ${\bm{r}}_{j},{\bm{r}}+\bm{\delta}$.

By using the workstation with 256GB memory and Xeon Gold 6130 processor, we perform the Lanczos algorithm to find the energy states roughly within $\pm 0.25$meV energy range in the BdG Hamiltonian (11). Consequently, we have the states with energy $E_{j}$ close to zero in the form of

[TABLE]

where $\bm{u}^{j}_{\bm{r}},\bm{\nu}^{j}_{\bm{r}}$ are 4-component vector. These states can reveal the hybridization of multiple MVMs. To connect the physical observable, by using these states we can compute the LDOS on the top surface of ${\rm FeTe}_{x}{\rm Se}_{1-x}$

[TABLE]

where $\delta(\beta E)=1/\cosh^{2}(\beta E)$. We choose $1/\beta=k_{B}T=20\mu$eV for the simulation leading to $76\mu$eV FWHM of the tunneling peaks. This value of the FWHM is close to the averaged value of the FWHM (80$\mu$eV) in the experiment 2018arXiv181208995M in the low magnetic field as shown in Fig. S4.

Appendix C Particle-hole symmetry broken by Majorana hybridization

For an isolated MVM with zero energy, its LDOS exhibits perfect particle-hole symmetry. On the other hand, Fig. 3(c) shows that when the intervortex distance of the Majorana triangular lattice becomes shorter, the LDOS of the vortex exhibits the broken particle-hole symmetry. The reason is that the Majorana hybridization and wavefunction overlapping break particle-hole symmetry. To capture the relation between the Majorana physics and particle-hole symmetry, we start with two isolated MVMs with zero energy located at ${\bm{r}}_{1}$ and ${\bm{r}}_{2}$

[TABLE]

where $\gamma_{1}({\bm{r}})=\bm{u}_{1}({\bm{r}})c^{\dagger}_{\bm{r}}+\bm{u}_{1}^{*}({\bm{r}})c_{\bm{r}}$, $\gamma_{2}({\bm{r}})=\bm{u}_{2}({\bm{r}})c^{\dagger}_{\bm{r}}+\bm{u}_{2}^{*}({\bm{r}})c_{\bm{r}}$, and $|\rm{BCS}\rangle$ is the BCS wavefunction. Here we use $\bm{\nu}_{i}^{*}=\bm{u}_{i}$ since a MVM with zero energy is its own antiparticle. By Eq. 21, each MVM has LDOS in the form of

[TABLE]

By definition, the peak is located at zero energy and the LDOS obeys $N_{i}({\bm{r}},E)=N_{i}({\bm{r}},-E)$ indicating preserving particle-hole symmetry. As these two MVMs couple weakly, the Majorana wavefunctions are kept in the original forms, and the hybridization Hamiltonian can be generally written in the second quantization form

[TABLE]

By defining the BCS wavefunction satisfying $(\gamma_{1}+i\gamma_{2})|\rm{BCS}\rangle=0$, we have an eigenstate $(\gamma_{1}-i\gamma_{2})|\rm{BCS}\rangle$ with energy $E=-t$. The explicit form of the wavefunction is written as

[TABLE]

where $\bar{\bm{u}}({\bm{r}})=\bm{u}_{1}({\bm{r}})+i\bm{u}_{2}({\bm{r}})$ and $\bar{\bm{\nu}}({\bm{r}})=\bm{u}_{1}^{*}({\bm{r}})+i\bm{u}_{2}^{*}({\bm{r}})$. The LDOS of the Majorana hybridized state is given by

[TABLE]

This is the general form of the LDOS. Now we consider the specific case that the hybridization of the two MVMs residing at two vortex cores (${\bm{r}}={\bm{r}}_{1},\ {\bm{r}}_{2}$) respectively. In the Hamiltonian (3), the two phases $(\theta_{i},\phi_{i})$ of the order parameter stem from the phase winding of vortex $\#i$ and other vortices respectively so that $\Delta=\Delta_{0}e^{i(\theta_{i}+\phi_{i})}$. By definition, we have $\theta_{i}+\phi_{i}=\theta_{j}+\phi_{j}$. Here we assume $\phi_{i}(r)$ is a constant for MVM $\#i$ since for vortex $\#i$ $\theta_{i}$ is used in the polar coordinate to solve the eigenvalue problem. With the additional phase $\phi_{i}$, the wavefunction (4) of the MVM $\#i$ has a modification

[TABLE]

Consider the LDOS at one of the vortex core, say ${\bm{r}}={\bm{r}}_{1}$. We simplify the problem by neglecting the magnetic flux effect ($\bm{A}({\bm{r}})$); the distribution of the two vortices gives the specific values of the three phases $\phi_{1}=-\eta,\ \theta_{2}=-\eta,\ \phi_{2}=\pi-\eta$ as illustrated in Fig. S5, where $\eta$ is an arbitrary angle. We have the explicit form of the hybridized wavefunction

[TABLE]

where $d\equiv|{\bm{r}}_{1}-{\bm{r}}_{2}|$ is the distance between the two MVMs. Here we use $J_{0}(0)=1,\ J_{1}(0)=0$. The LDOS at ${\bm{r}}={\bm{r}}_{1}$ can be separated to two parts

[TABLE]

where

[TABLE]

While the even function $N_{\rm{e}}(E)$ of $E$ preserves particle-hole symmetry, the odd function $N_{\rm{o}}(E)$ of $E$ breaks particle-hole symmetry. Due to the exponential decay $e^{-d/\xi}$, the distance between the two MVMs controls the strength of the particle-hole symmetry breaking. This idea can be extended to the relation between the intervortex distance and the particle-hole symmetry breaking in the Majorana lattice system. Hence, the short intervortex distance results in broken particle-hole symmetry in the vortex core as shown in Fig. 3(c).

Appendix D Simulation for a moving vortex

Fig. 6 in the main text shows the LDOS changes of the vortex movement only in the high magnetic field ($\approx 4$T). Here we further show the vortex movement in two lower magnetic fields ($1$T and $2$T) as shown in Fig. S6,S7. The relation between the Majorana hybridization and the moving vortex is still complicated. For example, as shown in Fig. S7 the vortex core $\delta$, which originally does not possess a ZBP, has a ZBP after vortex $\alpha$ moves away from vortex $\delta$. However, the overall behavior of the LDOS is similar to the high magnetic field. That is, the LDOS of the vortices near the moving vortex has significant changes. For vortices away from the moving vortex, the LDOS exhibits no change (Fig. S6$\delta$) or slight change (Fig. S7$\epsilon$).

Appendix E STM measurement of creeping vortex and vortex relocation

Fig. S8 shows the LDOS of the creeping vortex has a dramatic change after this vortex moved $2$nm. The range of the FOV is $15\times 15$nm, which is comparable to the Majorana coherence length ($\xi=13.9$nm). Therefore, it is possible that the surrounding vortices, which were not seen by the STM, also moved to affect the LDOS of the creeping vortex. Fig. S9 shows the LDOS of the vortex core was not affected by the vortex relocation, since the shortest distance between the probed vortex and the visible relocated vortices is roughly $50$nm ($>\xi$).

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

Fig. S10 shows (spin) LDOS line profiles of the two selected vortices in Fig. 4. As shown in panel (a), the spatial length ($5$nm) of the ZBP is consistent with the STM experiment observation of ${\rm FeTe}_{x}{\rm Se}_{1-x}$ Wang2018Evidence ; 2018arXiv181208995M . In addition, this observable length is identical to the spatial length of the splitting energy peaks as shown in panel (b). Panel (c,d) show that the spin up line profiles with and without ZBPs are almost the same with the line profile including both spins in panel (a,b). The spin down line profiles in panel (e,f) show that the spatial peaks appear away from the vortex cores regardless of the splitting of the energy peaks. The heights of the spin up peaks are roughly 6 times smaller than the ones of the spin up.

Appendix G Inhomogeneous superconducting gaps and chemical potentials

It has been observed that the (001) surface of ${\rm FeTe}_{x}{\rm Se}_{1-x}$ has inhomogeneous superconducting gap distribution (see supplementary materials in 2018arXiv181208995M ). We investigate whether this inhomogeneity affects the Majorana hybridization by examining LDOS of vortex cores in triangular vortex lattice. We use the spatial distribution of the superconducting gap based on the supplementary materials in 2018arXiv181208995M for the simulation. Fig. S11 shows in the presence of the inhomogeneous superconducting gaps each vortex in the same magnetic field has almost identical LDOS. Hence, the inhomogeneous gaps do not seriously affect the Majorana hybridization on the surface of ${\rm FeTe}_{x}{\rm Se}_{1-x}$. On the other hand, the inhomogeneous chemical potentials significantly affect the Majorana hybridization. As shown in Fig. S12 this inhomogeneity leads to different LDOS for each vortex core in the low and high magnetic fields.

Appendix H ZBP heights

As the magnetic field increases, the intervortex distance decreases so that the Majorana hybridization enhances. It is not surprising that the averaged ZBP height decreases with the increasing magnetic field. As shown in Fig. S13(a,b) the averaged ZBP heights of the LDOS in the vortex cores from the simulation exhibit this decreasing behavior. The experimental data only in region A have this similar behavior of the tunneling conductance as shown in Fig. S13(c). However, in Fig. S13(d) the data in region B do not exhibit deceasing ZBP heights. The reason of the inconsistency might be that in different magnetic fields the tunnelings of the STM were not in the same conditions. It is important to further experimentally check the decreasing averaged ZBP heights in ${\rm FeTe}_{x}{\rm Se}_{1-x}$.