Emergent vortex Majorana zero mode in iron-based superconductors

The vortex of iron-based superconductors is emerging as a promising platform for Majorana zero mode, owing to a magic integration among intrinsic vortex winding, non-trivial band topology, strong electron-electron correlations, high-Tc superconductivity and the simplification of single material. It overcomes many difficulties suffered in heterostructure-based Majorana platforms, including small topological gap, interfacial contamination, lattice imperfections, and etc. Isolated zero-bias peaks have been found in vortex of several iron-based superconductors. So far, studies from both experimental and theoretical aspects strongly indicate the realization of vortex Majorana zero mode, with a potential to be applied to topological quantum computation. By taking Fe(Te,Se) superconductor as an example, here we review original idea and research progress of Majorana zero modes in this new platform. After introducing the identifications of topological band structure and real zero modes in vortex, we summarize the physics behaviors of vortex Majorana zero modes systematically. Firstly, relying on the behavior of the zero mode wave function and evidence of quasiparticle poisoning, we analyze the mechanism of emergence of vortex Majorana zero modes. Secondly, assisted with some well-established theories, we elaborate the measurements on Majorana symmetry and topological nature of vortex Majorana zero modes. After that, we switch from quantum physics to quantum engineering, and analyze the performance of vortex Majorana zero mode under real circumstances, which may potentially benefit the exploration of practical applications in the future. This review follows the physics properties of vortex Majorana zero modes, especially emphasizes the link between phenomena and mechanisms. It provides a chance to bridge the gap between the well-established theories and the newly discovered iron home of Majoranas.


I. INTRODUCTION
Majorana zero mode (MZM) is a zero-energy quasiparticle with exotic properties. In a condensed matter system, it usually appears in a topological superconductor as a bound state [1][2][3][4][5]. MZMs have three fundamental properties. First of all, the hole and particle components of MZM are equal. The so-called "Majorana symmetry" requires self-conjugation of the generation and annihilation operators, demonstrating the charge neutrality of MZMs. MZM is a solid-state analogue to the Majorana fermions in the universe, of which the antiparticle is the particle itself. Similar to the Majorana fermions, MZMs can only exist individually. Two MZMs fuse into a normal fermion at a finite energy. The deviation from the zero energy increases exponentially with the two MZMs approaching to each other. Secondly, MZMs are topological in nature. The appearance of MZM is accompanied by nontrivial topological invariants. MZM can be viewed as a topological boundary state of an effective topological superconductor. Thirdly, MZMs obey non-Abelian anyonic statistics, owing to their fractional quantum dimension in a value of √ 2. MZM is regarded as "half" of a fermion in some literature. The braiding operations on these "halves" constitute a building block of fault- * dingh@iphy.ac.cn tolerant topological quantum computation [6][7][8][9]. These exotic properties propel the exploration of MZM on both fundamental physics research and potential practical application over the last decades. Majorana physics has become one of the most active topics in the frontier of condensed matter physics and quantum physics.
Owing to the particle-hole redundancy of superconductor, a superconducting quasiparticle is the quantum superposition of electron and hole. It provides an opportunity for emergence of MZM on some topological defects of a superconductor. MZM should be at the absolute zero energy, owing to the coexistence of Majorana symmetry and particle-hole symmetry (C † E = C −E ). Vortex is a typical topological defect. It appears spontaneously in type-II superconductors under a sufficient larger magnetic field [10][11][12]. In the vicinity of a vortex, the superconducting order parameter is spatially nonuniform. The amplitude of superconducting order parameter decreases when moving to the vortex center at which the amplitude vanishes, while the phase changes with the azimuth around the vortex. Note that the vorticity is defined as the ratio of the change of superconducting phase to the change of azimuth. Vortex can be regarded as a superconducting version of quantum well, in which superconducting quasiparticles form bound states. It can be described by Bogoliubov-de Genne (BdG) equation, H BdG = H 0 (r)τ z + ∆(r)e iθ(x,y) τ x (1) where H 0 (r) is the Fourier transformation of normal arXiv:2108.12850v1 [cond-mat.supr-con] 29 Aug 2021 states Hamiltonian H 0 (k) = ( k) 2 /(2m) − µ. ∆(r) is the amplitude of superconducting gap. θ(x, y) is the azimuth in the real space, describing the superconducting phase winding around the vortex. τ x and τ z are the Pauli matrices in particle-hole space. and m is the effective mass, and µ is the chemical potential. When the superconducting gap is nodeless, the solution of BdG equation manifests the well-defined vortex bound states. Owing to in-plane rotational symmetry, the energy of the vortex bound states can be expressed by the eigenvalue of angular momentum (ν). In conventional s-wave superconductors, the vortex bound state is known as Caroli-de Gennes-Matricon (CdGM) bound state [13][14][15]. Zero energy mode is forbidden in the level sequence, which can be expressed as E n ≈ ∆/k F ζ ≈ (n + 1/2)∆ 2 /E F , here n + 1/2 = ν, n is an integer, k F is the Fermi wavevector, E F is the Fermi energy, and ζ is the coherence length of the superconductor. Ultimately, the absence of zeroenergy vortex bound state is a requirement of quantum uncertainty principle. In a conventional superconductor, it is reflected in the antiperiodic boundary conditions of the quasiparticles, which is a natural result of the π-flux vortex.
There are two prerequisites in order to create MZM in a superconducting vortex. First, satisfaction of the zero mode condition. It requires to introduce an additional geometric phase with odd times of π into the BdG equation, so that, the energy level of vortex bound states turns to be E n ≈ n∆ 2 /E F . Second, stabilization of a single zero mode. It requires to elevate the spin degeneracy of the system. Following these clues, people first realized that the chiral p-wave superconductor (p x + ip y ) is a candidate which satisfies the requirements. The p x + ip y superconducting order parameter carries the intrinsic π phase, which counteracts the antiperiodic boundary condition of vortex, enables the zero-energy vortex bound state [16][17][18][19][20]. Theoretical analysis further manifests that the p x + ip y superconductors are intrinsic topological superconductors in the weak pairing phase [16]. In a spinless p x + ip y superconductor, a single MZM can be stabilized in a conventional vortex (φ = h/(2e)). But for spinful p x + ip y superconductors, two MZMs emerges in a conventional vortex, they are not protected and fuse to be a complex fermion immediately. Further studies resolve that the spin degree of freedom can be eliminated in a half-quantum vortex (φ = h/(4e)), in which the superconducting phase winding only couple to one of the spin components. Thus in the spinful case, a single MZM can be stabilized in a half-quantum vortex [17,18]. The idea of p x + ip y superconductors can be traced back to the studies of 5/2 fractional quantum Hall effect [21] in which a special ground state named as Paffian state (also known as Moore-Read state) was proposed. The Paffian state maps the spin-polarized electrons under a strong magnetic field into the spinless p x + ip y superfluid of composite fermions under the zero field. Majorana modes were proved to appear on topological defects (e.g. domain walls, vortex) of the Paffian state of 5/2 fractional quantum Hall [22]. Although the proposal of intrinsic topological superconductors is crystal clear in theory, both the p x + ip y superconductivity and the halfquantum vortices are extremely difficult to be realized in experiments [23]. It poses huge constrictions on the research of MZM at that time.
This deadlock was not broken until the advent of the Fu-Kane model in 2008. Fu-Kane model studies the quasiparticle behavior of a superconducting Dirac surface state in two dimension. It proves that a Dirac surface state proximitized by a conventional superconductor has equivalent low-energy quasiparticle as a spinless p x + ip y superconductor. MZM can emerge in a conventional vortex of a s-wave superconductor/topological insulator heterostructure. Fu-Kane Hamiltonian is where σ is the Pauli matrix in spin space. Under the spatial inversion transformation (r → −r), the spin of the Dirac surface state is rotated by 2π (give a minus sign) owing to spin-momentum locking (σ · k), thus the superconducting pairing c † k↑ c † −k↓ remain unchanged, compatible with the proximitized s-wave pairing. The Fu-Kane model satisfies the two prerequisites mentioned above for the creation of MZM. But unlike the p x + ip y superconductor, the extra phase, which is required to enable the zero mode, is provided by the spin-momentum locking of Dirac surface state, and the spin degeneracy of electrons is elevated by the spin-orbital coupling [24][25][26]. By taking the transformation c k = ψ ↑k + e iθ k ψ ↓k / √ 2, where θ k is the azimuth in k-space, Fu-Kane Hamiltonian becomes, Eq.(3) describes an isotropic electronic band taking spinless p x + ip y superconducting pairing. Thus it proves the equivalence between the Fu-Kane model and a spinless p x + ip y superconductor. A stable single MZM can emerge in a conventional vortex of the Fu-Kane model. Owing to the successful circumvention of p x + ip y superconductor, Fu-Kane model starts a new phase of Majorana research by facilitating the advances from theoretical hypothesis to experimental realization. Enlightened by the idea of Fu and Kane, various platforms were designed to realize Majorana modes in condensed matter systems, and some experimental evidence of MZM has been reported. These Majorana platforms can be categorized according to the physics mechanism as follows, 1) Zeeman gap assisted Rashba nanowires [31][32][33], including semiconductor nanowires [34][35][36][37] and Au nanowires [38][39][40]; 2) spin chains with strong exchange coupling, including magnetic atom chains [41][42][43][44] and magnetic carbon nanotubes [45]; 3) chiral topological superconductivity on a Yu-Shiba lattice [46][47][48][49][50][51]; 4) superconducting proximitized quantum anomalous Hall [52][53][54][55][56][57]; 5) Fu-Kane magnetic boundary states [58], including Majorana end typical band structure of iron-based superconductors, the orbital characters of each band are, α (dxz); β (dyz); γ (dxy); η (dxy); δ (dxz); ω (pz) ; (c) typical Fermi surfaces and superconducting order parameters of iron-based superconductors [27,28]; (d) mass renormalization among different compounds, indicating strong electron-electron correlations in iron-based superconductors [29]; (e) band structure near the Γ point of Fe(Te,Se) single crystals measured by ARPES [30]. mode in one dimension [59] and chiral Majorana modes in two dimension [60,61]; 6) helical Majorana modes in Fu-Kane Josephson junctions [58,62]; 7) MZM in Fu-Kane vortex [58,[63][64][65][66]. In these platforms, superconductivity is provided by an s-wave superconductor; and the realization of MZM relies on combing another topological material with the superconductor to form a heterostructure. Due to the low critical temperature (T c ) and small superconducting gap of the s-wave superconductor, the topological gap (∆ Top ≈ ∆ 2 /E F 0.1 meV), which separates the MZMs from the lowest excitations, is also very small. It requires a very low experimental temperature (T exp < 100 mK) to make a reliable observation of MZMs. It not only increases the experimental cost, but also makes MZM more vulnerable to other trivial quasiparticles, thus hinders the further studies of MZMs in experiments. In addition, the heterostructure induces uncontrollable disorders and fabrication complexity [67], and may induce other trivial mechanisms in the experiments which mimic Majorana-like signals, making the unequivocal confirmation of MZM difficult [37,57]. Therefore, a new leap in the Majorana research needs a paradigm innovation of material platform of MZMs.
In order to explore a better material platform of MZMs, here we review the difficulties suffered in the previous Majorana platforms systematically. Symptomatic solutions are found through the one-by-one analysis ( Fig. 1(a)). Firstly, the intrinsic p-wave superconductiv-ity is rare in nature and difficult to be realized. This difficulty is solved by the Fu-Kane model, by using the spin Berry phase of the topological surface state to replace the superconducting phase winding of intrinsic topological superconductor. This requires the coexistence of topological band structure and superconductivity. Instead of making a Fu-Kane heterostructure, we realized that in a single material such the coexistence requires a multiband structure to possess both topology and superconductivity. Secondly, the required experimental temperature is extremely low. A MZM is a superconducting quasiparticle, the survival temperature of the MZM is positively correlated to the value of T c . Incorporating a high-T c superconductor can effectively increase the survival temperature of the MZM. Furthermore, considering the necessary conditions of the existence of a welldefined vortex bound state, the superconductor should be full-gapped. So that the s± wave unconventional superconductors [68][69][70][71], and the nodeless d-wave unconventional superconductors [72][73][74][75] may be promising choices. Thirdly, the tiny topological gap, separating the MZM from the lowest quasiparticles, causes severe quasiparticle poisoning. In vortices, the topological gap is proportional to ∆ 2 /E F . A large topological gap requires a small Fermi energy, which can be realized in a material with strong electron-electron correlations. Fourthly, the heterostructure induces difficulties. The most direct way to remove those difficulties is to use a single material which integrates all of the three properties mentioned above. To sum up, a better Majorana platform could be potentially found in a material combining the properties of multiple bands, strong electron-electron correlations, high-T c superconductivity and topological band structure.
Our research experiences on both superconductors and topological matters inspired the idea that the iron-based superconductors could be such kind of candidates. First of all, the low-energy electronic states of iron-based superconductors [76,77] are mainly composed of the t 2g orbitals of Fe atom (d xz , d yz , d xy orbitals) and the p orbitals of chalcogenides. Near the Fermi level, three hole bands appear in the center of the Brillouin zone, and two electron bands appear at the corners of the Brillouin zone ( Fig. 1(b)). These bands pass through the Fermi level and form multiple Fermi pockets. It is observed experimentally that a full gap opens on the Fermi surfaces below T c [27]. Although a universally-acknowledged superconducting mechanism is not reached for the ironbased superconductors, a large number of theoretical and experimental results support the s±-wave pairing related to spin fluctuations ( Fig. 1(c)) [27,28,[78][79][80][81][82]. Secondly, many iron-based superconductors have strong electronelectron correlations [29], which leads to strong mass renormalization of their bands. For example, the mass renormalization factor of Fe(Te,Se) is close to that of cuprate superconductors (about 7) ( Fig. 1(d)). Therefore, the bulk bands of Fe(Te,Se) acquire very small Fermi energies (about 10 meV) ( Fig. 1(e)), which were clearly observed by angle-resolved photoelectron spectroscopy (ARPES) experiments [30].
Most importantly, some early experimental evidence of topological bands was found in Fe(Te,Se) single crystals that started the notion of the iron-based Majorana platform [87,98,99]. In 2014, by implementing in-situ surface dosing on FeTe 0.55 Se 0.45 single crystals, an electron band was found above the original Fermi level [98]. Later, the signature of band inversion was observed by synchrotron-based ARPES experiments with variable photon energy and polarization [87]. At the same time, a scanning tunneling microscopy/spectroscopy (STM/S) experiment found a robust zero-energy bound state on the interstitial Fe atom [99]. All these pieces of evidence suggest non-trivial topology in FeTe 0.55 Se 0.45 single crystals. These clues attracted the attentions on FeTe 0.55 Se 0.45 single crystals and inspired more systematic experiments to identify the superconducting topological surface states [100,101] and vortex MZM [102][103][104][105][106].
In this review, we systematically summarize the original ideas and research progress of the emergent vortex MZM in Fe(Te,Se) single crystals with rich details. We aim to bridge the gap between the classical Majorana theories and the emerging iron-based Majorana platform. The rest of the paper is organized as follows. In Section II, we introduce the topological band structure and the experimental observation of the surface states [87,100,101]. In Section III, the zero-bias conductance peak (ZBCP) was clearly observed in the vortices. By introducing a careful validation of vortex MZMs [102], we discuss the issue of "what" is the ZBCPs. In Section IV, we discuss "how" the vortex MZM emerges in a real material, e.g. Fe(Te,Se) single crystals. The wavefunction and quasiparticle poisoning behavior of the vortex MZM are carefully introduced. The effects of various realistic issues on emergence of vortex MZM are identified [102]. In Section V, we review the Majorana nature of the vortex MZM [105]. The experimental evidence of particle-hole equivalence is introduced in combination with the theory of Majorana induced resonant Andreev reflection. In Section VI, we review the topological nature of the vortex MZM [104]. Based on the energy eigenvalue and wavefunction distribution of the vortex bound states, we introduce the half-integer level shift of bound states, which is concomitant with emergence of vortex MZM. In Section VII, we switch our angle to focus the realistic details, which is an engineering concern. It has been observed that MZM is absent in some vortices on Fe(Te,Se). We analyze the possible mechanisms governing this phenomenon [104], which may benefit the exploration of practical applications in the future. In Section VIII, we present the conclusions and outlooks. The research progresses of the iron-based Majorana platform are briefly reviewed, and the main obstacles in the way to topological quantum computation are briefly analyzed.

II. TOPOLOGICAL BAND STRUCTURE
Iron-based superconductors can have multiple topological states, produced by the concurrence of strong interlayer coupling and spin-orbital coupling (SOC) [83,84,87,101]. Topological band inversion can be engineered by adjusting the k z dispersion of the p z orbital which belongs to the anions. Under appropriate conditions, the odd-parity p z orbital intersects with the even-parity d orbitals (d xz , d yz , d xy ) along Γ − Z. Those nontrivial band crosses form topological insulator state or topological Dirac semimetal state, contingent upon the same or different angular momentums of the crossing bands respectively. We note that the band inversion along k z is not unique to iron-based superconductors, it occurs ubiquitously in the 1T phase of transition metal dichalcogenides (such as PtTe 2 , PdTe 2 , and NiTe 2 ) [107][108][109]. In principle, topological band structure can be produced in all the iron-based superconductors under sufficient tuning. But some compounds are naturally so under ideal con- (a) first-principle calculation of band structure of FeSe (without SOC), the size of red circles represents the components of the pz orbital [87]; (b) crystal structure of Fe(Te,Se) [87]; (c) band inversion mechanism and orbital overlapping in Fe(Te,Se) [87]; (d) first-principle calculation of band structure of FeTe0.5Se0.5 (without SOC) [87]; (e) experimental band structure around Γ in FeTe0.55Se0.45 measured by ultra-high resolution laser ARPES [101]; (f) realistic topological band structure in FeTe0.55Se0.45 (with SOC). TDS stands for topological Dirac semimetal, and TI stands for topological insulator [101].
ditions, topological band inversions have been verified in Fe(Te,Se), LiFeAs and CaKFe 4 As 4 by first-principles calculations and ARPES measurements. Taking Fe(Te,Se) single crystals as an example, in this section we introduce the mechanism of band inversion and experimental progress on the studies of topological band structure.

A. Mechanism of band inversion
A first-principles calculation of FeSe single crystals is shown in Fig. 2(a) [87]. Three bands appear near the Fermi level in the absence of SOC. The corresponding irreducible representations are marked as Γ − 2 (p z /d xy antibonding orbital), Γ + 4 (d xy orbital), and Γ + 5 (d xz /d yz orbital), respectively. The odd-parity Γ − 2 band is always above the even-parity bands, so there is no topological band inversion in FeSe. Theoretically, the energy position of Γ − 2 band can be adjusted by the lattice constant. As shown in Fig. 2(c), the d xy orbital of Fe atoms lies on the Fe plane (shown as the cyan plane), its electron cloud spreads a little bit along the c-direction. On the contrary, the p z orbital of Se atoms is distributed vertically, with one end coupling with the d xy orbital on the Fe plane, and the other end crossing the van der Waals (vdW) gap and coupling with the p z orbital of the adjacent Se-Fe-Se layer. By changing the distance between the anion and the Fe plane, the relative strength of the intralayer p-d coupling and the interlayer p-p coupling can be effectively controlled. The competition between them determines the behavior of the Γ − 2 band. The tellurium substitution on the anionic site can increase the distance between the anion and the Fe plane, thus weakens the intralayer p-d coupling. So the Γ − 2 band moves to lower energy to adapt to this change. On the other hand, the tellurium substitution enhances the interlayer p-p coupling, which leads to a larger k z dispersion of the Γ − 2 band. As a consequence, the odd-parity Γ − 2 crosses the even-parity Γ + 4 and Γ + 5 along Γ-Z, producing the topological band inversion at Z (Fig. 2(d)).
With SOC, the Γ − 2 band evolves to be Γ − 6 , which is composed of the p z orbital; the Γ + 4 band evolves to be Γ + 7 , which is composed of the d xy orbital; and the Γ + 5 band evolves to be Γ + 6 and Γ + 7 , which are composed of the d xz orbital and the d yz orbital respectively. DFT calculations of FeTe 0.5 Se 0.5 single crystals show [87] that the d xy band (Γ + 7 ) is above the d yz band (Γ + 7 ), both of them cross the Fermi level around Γ. It leads to three intersections between the odd-parity p z band and the even-parity bands along Γ-Z. However, ultra-high reso- lution ARPES experiment observed only one bulk band (d yz ) crossing the Fermi level at Γ (Fig. 2(e)) [101]. This discrepancy invokes the investigation of the realistic details of the band inversion along Γ-Z. Owing to orbitaldependent mass renormalization, the band width of the d xy orbital is very small [29]. It pushes the d xy band sinking below the Fermi level, leads to the absence of the d xy Fermi surface. In addition, the d xy band intersects with the d yz band below the Fermi level, and a hybridization gap is clearly observed by ARPES around the intersection ( Fig. 2(e)) [101], which is an expected result owing to the same irreducible representation of the two bands (Γ + 7 ). Later on, theorists reproduced the realistic band structure (with SOC) of FeTe 0.55 Se 0.45 single crystals by using the k · p model ( Fig. 2(f)) [101]. Different from the DFT calculation, there are only two topological intersections along Γ-Z, one is a gapped intersection (p z and d xz ), in which a strong topological insulator state occurs, the other is a protected intersection (p z and d yz ), which is a topological Dirac semimetal state protected by the rotational symmetry. The strong topological insulator state is protected by the size of the hybridization gap. Tellurium substitution enhances the hybridization gap due to a larger SOC. It is beneficial for emergence of Dirac surface states.

B. Discovery of topological insulator state
Several pieces of evidence of Dirac surface states were observed by ARPES, confirming the topological insulator state in FeTe 0.55 Se 0.45 single crystals [100]. First, The linear-dispersed Dirac surface state was clearly observed in ARPES spectra. Second, the orbital characters were identified by the matrix element analysis in a polarization dependent ARPES [28,110]. Third, the signature of spin-momentum locking was directly observed through spin-resolved experiments [24][25][26]. Fourth, an ultra-low temperature ARPES measurement found an isotropic superconducting gap opening on the Dirac surface state below T c (14.5 K). Thus the requirements of Fu-Kane model [58] are fulfilled in a single material platform of FeTe 0.55 Se 0.45 , rather than a heterostructure as studied previously ( Fig. 1(a)). Furthermore, the E F of the Dirac surface state was found to be very small in FeTe 0.55 Se 0.45 , indicates a potential for a clear observation of vortex MZM.
ARPES is capable to distinguish both the momentum and the energy of electrons simultaneously, so it is widely used in band structure measurements on topological materials [111]. It is well known that, the spectral intensity of ARPES [112,113] is where f (ω) is the Fermi-Dirac function; is the spectral function of the material; is the ARPES matrix element, and determines the orbital selection rule of ARPES experiments. The ARPES signal exists only if M k f,i is nonzero. Since the emission photoelectron is a plane wave, ψ k f is an even function relative to the signal acquisition plane (namely the slit plane of the ARPES analyzer, see the pink plane in Fig. 3(a)). In order to get nonzero ARPES signal, the parity of photon polarization (A · p) and the orbital ψ k i should be same. In the case of normal emission, the wavefunction ψ k f is also an even function with respect to the plane vertical to the signal acquisition plane (the blue plane in Fig. 3(a)). The parity of A · p and ψ k i relative to the vertical plane also affects the selection rule. We summarize the orbital selection rule under different experimental conditions in the table of Fig. 3(a) [100].
The topological insulator state in Fe(Te,Se) single crystals is originated from the topological band inversion of p z and d xz bands. Those two orbitals compose the main orbital characters of the Dirac surface states due to the SOC. The expected orbital characters around Γ are shown in Fig. 3(d). The experimental geometry of the p-polarization contains both A y and A z , which is beneficial to the observation of p z orbital components, but suppresses the signal of the d xz orbital at Γ (Fig. 3(e)). Experimental results are shown in Fig. 3(b). A cone-like band structure can be clearly observed, but the top of the hole band is rather ambiguous which can only be distinguished in the second derivative. On the other hand, the experimental geometry of s-polarization only contains the A x polarized photons, which strongly suppresses the signal of the p z orbital near Γ, so that the topological surface states are almost invisible. But the s-polarization has a good selectivity to the d xz orbital, and indeed the hole-like d xz band was well resolved in the experiment (Fig. 3(c)). These results are fully consistent with the matrix-element analysis which resolves the orbital characters of each band. In summary, the hole-like bulk band below the SOC gap is mainly composed of the d xz orbital, the upper half branch of the Dirac surface states is mainly composed of the p z orbital, the lower half branch is composed of the p z and d xz orbitals, and the electronlike bulk band above the SOC gap can be inferred as the p z orbital ( Fig. 3(d)).
Spin-momentum locking is a fingerprint of the Dirac surface state. The Dirac electron obtains a π-flux spin Berry phase on the Fermi surface ( Fig. 4(a)), which is decisive evidence of topological insulator state. As shown in Fig. 4(a), a pair of spin-polarized energy distribution curves (Cut 1 and Cut 2) are measured at the selected Fermi points of the Dirac surface states ( Fig. 4(b) and Fig. 4(c)). The spin signals are reversed respect to the Dirac point, indicating that the spin polarization of the upper branch and lower branch of the Dirac surface states are different. In addition, the spins of Cut 1 and Cut 2 are opposite at the same energy. These results identify the spin-momentum locking of the observed lineardispersed bands, providing strong evidence of the Dirac surface state.
FeTe 0.55 Se 0.45 single crystals is regarded as a s±wave superconductor [79,114] (the momentum resolved s±-wave superconducting order parameters is shown in Fig. 1(c)). Since the Dirac surface state is mainly concentrated in the center of the Brillouin zone (k F ≈ 0.02Å −1 ), the superconductivity of the Dirac surface state, which is proximitized by the bulk states at different momenta, can only have a conventional s-wave. An early ARPES experiment observed two superconducting gaps (∆ 1 = 1.7 meV, ∆ 2 = 2.5 meV) near Γ, and one superconducting gap (∆ 3 = 4.2 meV) near M [115]. By incorporating the next-next-nearest-neighbor antiferromagnetic exchange coupling (J 3 ) in a strong coupling model [70], a theoretical gap equation, that is ∆ = |J 2 cos (k x ) cos (k y ) − J 3 [cos (2k x ) + cos (2k y )] /2|, was proposed to understand the observed three gaps in the experiment. It was found that ∆ 2 and ∆ 3 are scaled very well with the gap function. The extracted J 2 /J 3 values are close to the experimental results of neutron scattering [116]. However, the minimum gap ∆ 1 has the smallest k F among the three gaps, and greatly deviates from the expectations of the J 1 -J 2 -J 3 model [115]. These results suggest that ∆ 2 and ∆ 3 follow s±-wave pairing, while ∆ 1 may be with a different origin. Due to the improvement of ARPES resolution in the recent years, we achieved clearer observation of the superconducting gaps on FeTe 0.55 Se 0.45 [100]. New experiments show that the Dirac surface state opens a superconducting gap of 1.8 meV (∆ surface = 1.8 meV) at low temperature, with an isotropic distribution along the Fermi surface ( Fig. 4(e)). By measuring the temperature dependent energy distribution curves at the k F of the Dirac surface state, the gap closing temperature is resolved at about 14.5 K. ∆ surface , which is at a similar value as ∆ 1 in the earlier experiment, differs from the behavior of s±-wave, but conforms to the picture of the induced gap by k-proximity effect. It provides a reason- able understanding for the ill-scaling of ∆ 1 to J 1 -J 2 -J 3 model. Combined with different experimental measurements [30,79,[99][100][101][102][115][116][117][118][119], we summarize the superconducting gaps of FeTe 0.55 Se 0.45 single crystals as follows: near Γ, ∆ surface = 1.8 meV and ∆ dyz Bulk = 2.5 meV; near M, ∆ M Bulk = 4 meV. It is worth noting that the surface-sensitive experimental methods (e.g. ARPES and STM) are easier to observe the superconducting gap of the surface states [119].
The experimental parameters of the topological insulator state of FeTe 0.55 Se 0.45 are summarized in Fig. 4(f). Besides the basic features, two details are worth to be especially mentioned. First, the SOC gap of the topological insulator state is about 20 meV which is much smaller than that of other strong topological insulators (e.g. about 300 meV in Bi 2 Se 3 ) [24]. Secondly, the Dirac point of the topological surface states appears near the Fermi level, which eliminates the requirement of tuning the chemical potential. The Fermi energy of the Dirac surface state is only about 4.4 meV. The ratio of ∆/E F (about 0.5) is very large in FeTe 0.55 Se 0.45 single crystals, which is more than two orders of magnitude larger than that in Bi 2 Se 3 (about 10 −3 to 10 −2 ). This ideal parameter brings a hope for the realization of MZM with a large topological gap. In this sense, Fe(Te,Se) is a gift from nature [120].

C. Discovery of Dirac semimetal state
Topological Dirac semimetal is protected by the rotational symmetry of the principal axis [121]. As shown in Fig. 2(f), the Dirac semimetal state is composed by the p z and d yz orbitals in FeTe 0.55 Se 0.45 single crystals. Since the Dirac point is above the Fermi level, the band structure of the bulk Dirac cone cannot be directly measured by ARPES. However, the Dirac semimetal state can be identified in ARPES experiments through the anomalous spin-polarization on a bulk band, owing to the overlap Owing to the absence of the Fermi surface of the d xy band, the Dirac band may dominate the transport behavior on Fe(Te,Se). As shown in Fig. 5(g), the magnetoresistance shows semi-classical (quadratic) behavior below 6 T, but becomes linear between 6 and 40 T. Although the anomalous linear magnetoresistance can be explained by several conventional mechanisms, such as the averaging effect in polycrystalline materials with onedimensional Fermi surface [125] and heavily disordered systems [126], they are not consistent with the conditions of FeTe 0.55 Se 0.45 single crystals.
In the following we show that the magnetoresistance behavior can be well explained by the Dirac semimetal state. Linear magnetoresistance also exists under the quantum limit of Landau levels. The electron degeneracy on the Landau level increases with increasing magnetic field (the degeneracy density is g = eB/h). When the magnetic field is larger than the threshold of quantum limit, all the carriers occupy the lowest Landau level [125,127]. A strong magnetic field is not the only requirement for realization of the quantum limit, the energy gap between the lowest Landau level and the first excited state (∆E LL ) should be also larger than the Fermi energy (E F ) and thermal fluctuation (k B T ) [128]. For the conventional parabolic bands, the spectrum of Landau level is the type of harmonic oscillators, namely E n = eB/m * (n + 1/2) − µ, where m * is the effective mass and µ is the chemical potential. The quantum limit condition is, ∆E LL = eB/m * > k B T . In this case, the critical magnetic field of Fe(Te,Se) is estimated to be about 60 T (the parameters used in the estimation are T = 16 K, m * = 5m e ) [29,100,101], which is at odd with the experimental conditions (T = 16 K, B max = 40 T). For the Dirac bands, the spectrum of Landau level is E n = E D + sgn(n)ν F 2 eB|n|, where E D is the energy of the Dirac point, and ν F is the Fermi velocity. Since ∆E LL = ν F √ 2 eB the critical magnetic field is estimated to be about 2.8 T (the parameters are T = 16 K and k F = 0.1Å −1 ), that agrees with the experimental observations very well. Note that in an earlier magnetoresistance measurement the linear magnetoresistance occurs at 2 T [128], which is even closer to our estimation.
Later on, the Dirac surface states of Fe(Te,Se) were verified by other experimental groups [129,130]. The topological band structure of Fe(Te,Se) single crystals is summarized in Fig. 5(f), multiple topological states provide good opportunities for creating MZM in iron-based superconductors. It is worth noting that only the states around the Fermi level influence quasiparticle excitations, emergence of MZM does care about the circumstance of the Fermi surface. Iron-based superconductors have high mobility of normal carriers, which indicates that the chemical potential is difficult to be controlled by conventional methods such as the solid gating effect [131]. In Fe(Te,Se) single crystals, the topological insulator state appears at the Fermi level. As a result, in the following sections we will focus on the topological insulator state and introduce the vortex MZM emerges from it ( Fig. 6(a)).

III. VORTEX MAJORANA ZERO MODES
Since both the bulk and the surface bands have large value of ∆ 2 /E F , the topological gap between the vortex MZM and the lowest excitations is expected to be relatively large in FeTe 0.55 Se 0.45 single crystals, which is beneficial for the observation of a pure MZM. STM technique has excellent energy and spatial resolution, as well as multi-dimensional adjustability including temperature, magnetic field, and tunneling conditions [132,133]. In the case of single-particle tunneling in vacuum, the tunneling current between the tip and the sample is where V is the bias voltage; T is the temperature; x, y is the in-plane position of the sample; S is the distance between the tip and the sample; τ is the tunneling matrix element. When the tip wavefunction is swave, the matrix element can be written as τ (S, V, E) = exp −2 m e (φ t + φ s − 2E + eV )S/ , where φ t and φ s are the work functions of the tip and the sample respectively; ρ t is the tip DOS; ρ s is the spatially resolved local DOS of the samples; f (E, T ) is the Fermi-Dirac function.
When the bias voltage is small, τ is approximately irrelevant to the bias voltage V . When the tip position is fixed and the tip DOS is a constant, the local DOS of the sample can be directly measured by the differential conductance. It can be easily proved under the conditions mentioned above, that dI/dV (V, x, y) ∝ ρ s (V, x, y). STM is an ideal technique to validate the existence of vortex MZMs, which appear as a zero-energy peak on the density of states (DOS) [58,100], by discovering the ZBCPs on the dI/dV spectrum. FeTe 0.55 Se 0.45 single crystals was cleaved in-situ under ultra-high vacuum [102]. The fresh surface was studied by STM under ultra-low temperature. Since the Se-Fe-Se "sandwich" is weakly bonded through vdW interaction, it is cleaved in between the vdW gap, resulting in a Se/Te atomic bare surface. An atomic-resolved topography is shown in Fig. 6(b), in which Te atoms with a large atomic radius are shown as the bright spots, and Se atoms with a relatively small atomic radius as the dark spots. That indicates a good spatial resolution realized in the experiments, beneficial for the measurements on the wavefunction behavior of the vortex bound states. Under a magnetic field of 0.5 T perpendicular to the sample surface, a vortex lattice can be observed on the zero-bias conductance map (Fig. 6(c)), which is an important feature of the type-II superconductors [10]. The vortex size measured in STM experiment reflects the value of coherence length. It is different from other scanning probe microscopies such as scanning magnetic microscopes and scanning Hall microscopes where the vortex size is controlled by penetration depth (λ). At the vortex center, a very sharp ZBCP appears on the dI/dV spectrum, isolated in the middle of the superconducting gap. It is most likely a single zero-energy vortex bound state. At the position far away from the vortex center, the dI/dV spectrum recovers the shape of superconducting gap ( Fig. 6(d)). The observation of zero-energy vortex bound state is an important signature of the existence of MZM [58]. STM technique with higher spatial accuracy improves the reliability of the ZBCPs in Fe(Te,Se), which is an advantage comparing to transport methods [36,134].

A. Experimental validations
In order to verify the observation of a true MZM, the signal of ZBCP should be experimentally determined as a real single zero-energy vortex bound state. The exper-  [102]; (c) zero-bias conductance map which shows vortex lattice [102]; (d) a sharp zero-bias conductance peak measured at the center of a vortex. In order to make sure the observation of a zero-energy vortex bound state, three careful checks are listed as follows. Firstly, to make sure that the signal measured is indeed from vortex bound states [102], ZBC maps after and before applying a magnetic field was measured in (e). It shows the local environment of the vortex is clean and free of impurities [102]; (f) ZBCP is stable under different tunneling barriers. Secondly, to make sure that the observed ZBCP is truly a single peak [102]; (g) FWHM of ZBCP measured under different tunneling barriers; Thirdly, (h), (i) the observed ZBCP is truly a zero mode [102]. (h) simultaneous measurements of I(V ) curve and dI/dV curve on the center of a vortex core [102], (i) enlarged display of red box area in (h) [102].
imental tests were performed as follows.
First, the observed subgap state should be a real vortex bound state, rather than a quasiparticle peak caused by impurities, inhomogeneities or disorders. Before measuring the vortex, the material is characterized under the zero field. By performing the measurements of topography, zero-field dI/dV spectrum, and zero-field zero-bias conductance map, the good areas can be identified by the feature of ordered surface, the good spectral feature of the superconducting gap and a nearly vanishing conductance at the zero-bias, respectively. These areas were selected for the further studies under magnetic fields. For example, the vortex shown in the left panel of Fig. 6(e) supports a ZBCP under 2 T (Fig. 6(d)). As a check, a zero-field zero-bias conductance map was measured at the same spatial locations (right panel of Fig. 6(e)). The zero-bias DOS in the region almost vanishes without applying magnetic field, indicating no superconducting quasiparticle caused by other effects. Thus the DOS observed in Fig. 6(d) is totally caused by the vortex. For a further verification, a variable-tunnel-coupled dI/dV measurement was performed at the vortex ( Fig. 6(f)). Under the feedback regulation, the tip height of STM is adjusted accordingly as per the barrier conductance. The tunneling coupling can be manifested by the barrier conductance, which is defined as G N ≡ I t /V s . A large conductance corresponds to a small tip-sample distance and strong tunneling coupling. It was observed that the ZBCPs remain robust over two orders of magnitude in the barrier conductance (in unit of 2e 2 /h), which is consistent with the behavior of vortex bound states. On the contrary, for the quasiparticle states caused by impurities, inhomogeneities, or disorders, the energy of the peak changes with the tunneling coupling. For example, the impurity bound states share some similar conditions as the Andreev bound states in a quantum dot of hybrid nanowires. The tunneling coupling can influence the chemical potential, and can change the peak energy too [135].
Second, the observed subgap states should be a true single quantum state, rather than a mixture of multiple near-zero energy peaks. By fitting the spectra in Fig. 6(f), the full width at half maximum (FWHM) of the ZBCPs was extracted under different tunnel barriers. The measured FWHMs are less than 0.29 meV and basically remain unchanged with the change of the tunnel barrier ( Fig. 6(g)). Compared to the total energy broadening of the system which is about 0.28 meV [102], it shows that the measured width of ZBCP is mainly determined by the system resolution, strongly supports the observed ZBCP as a true single quantum state. The measurements displayed in Fig. 6(f) are within the weak tunneling regime. A scaling function (see Section V for details) is capable to fit the conductance values under different tunnel barriers [102]. It can be estimated that the intrinsic Majorana width, which is determined by the tunneling coupling (Γ t ), is approximately 2Γ t ≈ 170GN e 2 /h (µeV), much smaller than the thermal broadening (3.5k B T ) and the resolution of the STM. Therefore the FWHM of ZBCP is unchanged as a function of tunnel barrier ( Fig. 6(g)), manifesting the width evolution of a single quantum state under the weak tunneling regime. In addition, it is worth noting that under the strong tunneling regime the tunneling coupling increases faster than a linear function of G N , and the specific relation alters as a different thermal and dissipative broadening (see Section V A for details).
Third, the observed subgap states should be a true zero mode. The problem of zero-bias offset exists in STM. Generally, the zero voltage picked from the voltmeter is not the physical zero bias. In order to ensure that the measured ZBCP is a true zero mode, one needs to calibrate the zero drift carefully. The current and the voltage have a linear relationship at small bias voltages, thus the intersection of the I(V ) curves measured under different tunnel barriers is the true zero point of current and voltage. As shown in Fig. 6(h), the I(V ) curves at different tunnel barriers are measured at the vortex center. A dI/dV curve outputting from the lock-in amplifier is displayed simultaneously. The area marked by a red box is enlarged and displayed in Fig. 6(i). Obviously, the peak position of the dI/dV curve is exactly consistent with the intersection of the I(V ) curves with the error bar less than the sampling interval of the dI/dV curve (14 µeV).

B. Exam alternative explanations
The experimental validations demonstrate the observation of single isolated zero-energy vortex bound state in a material with the superconducting topological surface state, which fulfills the requirements of vortex MZM in Fu-Kane model. However, the ZBCP is a common feature appearing in the experiments of hybrid superconducting structures, that obscures the Majorana explanations. Although we performed experimental validations as mentioned above, for the sake of rigorousness, here we list and discuss alternative mechanisms that can lead to a ZBCP. These mechanisms caused difficulties in identifying a ZBCP as a MZM in other systems.
1) coherent Andreev reflection [136][137][138]. The most important mechanism in this category is the reflectionless tunneling. On a disordered superconductor-normal metal interface, scattering centers can cause the mirrorreflected electrons to shoot back to the sample again. The holes generated by the Andreev reflection can return along the incident path of the electrons. The phase conjugation of hole and electron induces an excess current at the zero bias. This effect can be destroyed by applying an external magnetic field which breaks the phase conjugation.
2) incoherent Andreev reflection [139,140]. A ZBCP can be induced by a cumulative effect. The possibility of incoherent Andreev reflection is very small in the weak tunneling regime, which is at odds with the observation of a strong ZBCP in vortex.
3) Kondo effect [141][142][143][144]. The zero-energy DOS can be induced by the spin-flip resonance when the conduction electrons couple to a localized impurity with a degenerate quantum ground state. The ZBCP of Kondo resonance splits under a magnetic field [144]. 4) Josephson current [145][146][147][148]. Cooper-pair tunneling induces a sharp ZBCP and a negative conductance at slightly lager bias voltages on the dI/dV spectrum. It can be well explained by the Josephson model with thermal fluctuations [147,148]. STM experiments use a metallic tip which cannot induce Josephson currents. 5) disorder-induced zero mode [149][150][151][152][153]. The most important mechanism in this category is the class-D weak antilocalization (WAL). When the size of the system (L) is longer than the mean free path (l), but shorter than the phase coherence length (l φ ), it is in the regime of quantum diffusive transport. In conventional materials, quantum diffusive transport can lead to weak localization (WL), which is proposed as the precursor of Anderson's localization in some literatures [154]. In topological materials, the additional geometric phase (π) reserves the interference condition of quantum diffusive transport, resulting in WAL and ZBCP [155]. Normally, a magnetic field destroys the electronic coherence, and eliminates the conventional WL and WAL. However, in the class-D superconductors (with time-reversal symmetry broken), the electronic coherence can be realized in an additional pathway, that is the electron and hole path of the Andreev reflection, it can induce WAL and ZBCPs even under a large magnetic field [150,151]. In principle, this mechanism can appear in Fe(Te,Se) single crystals, but is avoidable by selecting an good area without disorder, which is not difficult for the STM technique with high spatial resolution. 6) charge potential fluctuation induced zero-energy Andreev bound states [156][157][158][159]. This phenomenon commonly exists in quantum dots and the ends of hybrid nanowires [37,160]. In principle, it can also exist in a Fe(Te,Se) single crystals. The risk of this mechanism can be effectively avoided with the help of the zero-field test, a good area can be selected for measurements. If the ZBCP is caused by a smooth potential, it should already appear in the zero-field. 7) Yu-Shiba-Rusinov state of a magnetic impurity [161,162]. In experiments, a good area without magnetic impurity can be chosen for vortex measurement to avoid this risk (Section III A). 8) zero-energy surface Andreev bound state of unconventional superconductors [163][164][165][166][167][168][169][170][171][172][173][174]. The π phase interference between particle and hole trajectories could induce a zero-energy Andreev bound state. It occurs when the tunneling direction meets horizontal line node of the superconductor. For instance, a sharp ZBCP shows up in the tunneling spectrum of a d-wave superconductor measured along 110 . The gap structure of Fe(Te,Se) does not support this mechanism. 9) packed near-zero-energy vortex bound states. When the level spacing of the vortex bound states is very small, multiple near-zero-energy bound states are crowding near the zero energy. Under the condition of limited resolution, the multiple vortex bound states are convoluted to be a fake ZBCP [13][14][15]175]. This phenomenon is commonly observed in conventional superconductors (such as NbSe 2 ). In Section VI, we introduce the observation of the discrete levels of vortex bound states in Fe(Te,Se), which safely excludes this concern.

IV. MAJORANA MECHANISM OF EMERGENCE
The simplification of Fu-Kane model avoids secondary factors that interfere the elaboration of Majorana physics, and demonstrates emergence of zero-dimensional vortex MZM from a two-dimensional superconducting Dirac surface states [58]. However, single crystals such as Fe(Te,Se) are not as perfect as the theory. The three-dimensionality of the real materials has influence on MZM. Under a vertical magnetic field, there are not only zero-dimensional vortices on the surface, but also one-dimensional vortex lines throughout the bulk [12]. The zero-dimensional vortex on the surface is actually the ends of the one-dimensional vortex lines. Moreover, FeTe 0.55 Se 0.45 single crystal is not a perfect topological insulator. The coexistence of Dirac surface states and bulk bands (Fig. 2(e)) is an inherent requirement of the superconducting self-proximity effect. There are multiple Fermi surfaces formed by the Dirac surface states and other bulk bands [100]. Furthermore, FeTe 0.55 Se 0.45 single crystals suffers multiple imperfections including disorders, inhomogeneities, and defects, which induce quantum dissipation of the quasiparticle bound states. Thus a one-by-one experimental investigation on those factors can help our understanding of emergence of the vortex MZM, by clarifying the promoting and destructive issues of formation of MZM.
In this section, with the help of wavefunction behavior and quasiparticle poisoning of vortex MZM [102], we introduce the mechanism of emergence of vortex MZM in realistic circumstance of Fe(Te,Se). It shows that MZM is produced by the Dirac surface states, and damaged by the material-imperfection-induced quasiparticle poisoning. We also compare the two-dimensional Fu-Kane model [58] with a three-dimensional vortex line model in the last part of this section [176][177][178]. The behavior of vortex MZM can be well described by Fu-Kane model at sufficient low temperatures.

A. Wavefunction of vortex MZM
A MZM wavefunction occupies a certain spatial volume which is determined by the coherence length of the superconducting Dirac fermions. ZBCPs can be observed within the certain spatial range around the vortex center. The conductance of ZBCP can be calculated by modular square of MZM wavefunction [58,179], that is where where r is the distance from the vortex center; J i (x) is a Bessel function; ∆(r) is taken as a step function [102]. Thus the modular square of MZM wavefunction is, Note that the truncation threshold of ∆(r) does not affect the results. C is a normalized constant; ∆ 0 , E F and ξ represent the superconducting gap, the Fermi energy and the coherence length of the Dirac surface states, respectively. The dI/dV spectra were carefully studied along the line crossing a vortex (Fig. 7(a)). It was found that a subgap ZBCP remains at the zero energy as measuring at various spatial positions, and the intensity of ZBCP decreases when moving away from the vortex center (Figs. 7(b)-(d)). These observations are fully consistent with the theoretical expectation of MZM. The spatial distribution of the zero-bias conductance of ZBCP was extracted in Fig. 7(e). In order to determine the parameters of the Dirac surface states, we compared a wide-range dI/dV spectrum with the ARPES spectrum ( Fig. 7(f)). The features at 6.1 meV and -14.9 meV of the dI/dV spectrum are in good agreement with the bottom of the conduction band and the top of the valence band shown in the ARPES measurement respectively. A linear depression, which is the DOS feature of a Dirac point [180], appears at -4.4 meV of the dI/dV spectrum, which is coincident with the position of the Dirac point observed in ARPES. Based on these results, the experimental parameters of the Dirac surface states can be determined as follows, ∆ 0 = 1.8 meV, E F = 4.4 meV, ξ = ν F /∆ 0 = 123Å. By substituting these values into the MZM wavefunction (Eq.(11)), the Fu-Kane model fully reproduces the spatial distribution of ZBCP observed in our STM experiment (Fig. 7(g)). It attributes the ZBCP to an underlying Dirac surface state, thus supporting the appearance of Majorana zero mode in the vortex.
Although the vortex line penetrates throughout the three-dimensional bulk, the two-dimensional approximation of Fu-Kane model is still valid for the analysis of behaviors of the vortex MZM under sufficient lowtemperatures. In this situation, the three-dimensionality is not important for the Majorana modes. A MZM can be regarded as a topological quasiparticle appearing on the point defect (vortex) of a two-dimensional topological superconductor. The influences of the one-dimensional vortex line on emergence of vortex MZM will be discussed later (Section IV C and Section VII). The appearance of one-dimensional dispersive vortex bound states can affect the evolution of vortex MZM on the surface, manifesting the three-dimensionality of the material.

B. Evidence of quasiparticle poisoning
In perfect s-wave superconductors, conduction electrons condense to Cooper pairs completely at zero kelvin, leads to a hard gap without any subgap quasiparticles [181]. When temperature is raised, thermally-excited quasiparticles gradually occupy the subgap states. At a finite temperature, the equilibrium quasiparticle density is where ρ N (0) is the normal state DOS at the Fermi level. The quasiparticle lifetime can be estimated as follows, However, multiple experiments have found that the quasiparticle density far exceeds the equilibrium value (Eq.(12)), implying a reduction of quasiparticle lifetime [182,183]. The non-equilibrium occupation of superconducting quasiparticles applies severe disturbance on experiments, but its origin has yet to be well understood. As a result, it was called "quasiparticle poisoning" [182][183][184][185][186][187][188][189][190]. The degree of quasiparticle poisoning can be characterized by the "hardness" of the superconducting gap. A hard gap has steep gap edges and low subgap DOS, indicating negligible amounts of unpaired quasiparticles.
Quasiparticle poisoning is a dissipation process that the coherent quantum states couple to a fermionic bath.
It reduces the lifetime and weakens the signals of the quantum states in the DOS spectrum [184][185][186][187]. The conservation of fermion parity is one of the prerequisites of the topological quantum computation bases on Majorana qubits. Quasiparticle poisoning changes fermion parity randomly, causes information loss of Majorana qubits. It constrains the speed of gate operations of MZMs which should be safely within the lifetime determined by the poisoning rate [184][185][186][187][188][189][190]. Moreover, from experimental view, the quasiparticle poisoning reduces the zero-bias conductance and increase the energy broadening of MZM, causing the disappearance of MZM signal when the experimental resolution is limited. Quasiparticle poisoning is one of the most destructive factors to MZM and topological quantum computation. In ironbased superconductors, the evidence of quasiparticle poisoning was carefully studied. In addition to the basic quasiparticle poisoning, experiments show that the bulk bands may provide additional quasiparticle poisoning at high temperatures, as introduced below.
First, the spatial variation of basic quasiparticle poisoning is studied at the STM base temperature (0.55 K).
As shown in Fig. 8(a), the dI/dV spectra of the vortex center (red curve) and vortex edge (black curve) are measured on three vortices. A sharp ZBCP appears at the vortex center, and a gap feature appears at the edge. The basic quasiparticle poisoning can be quantified by the subgap DOS extracted from the spectra of vortex edge, and specifically, the quasiparticle background was defined as the integral of dI/dV over the range from -1 to 1 meV. By comparing the FWHM of ZBCP at the vortex center and the quasiparticle background at the vortex edge, we found that stronger poisoning correlates with larger energy broadening of MZM, indicates a shorter quasiparticle lifetime (the lower right of Fig. 8(a)). It makes the ZBCP broadening of some vortex MZMs much larger than the energy resolution of the system. At the positions with weak poisoning, the broadening of MZM is close to the resolution limit of the system. Obviously, the degree of quasiparticle poisoning of vortex MZM is spatially nonuniform, which is related to the intrinsic inhomogeneity of FeTe 0.55 Se 0.45 single crystals (Section VI C and Section VII for details).
Next, we introduce the additional quasiparticle poisoning effect which is investigated by temperature dependent measurements on several selected vortices. In order to minimize the influence of vortex creep [191,192], the temperature was ramping up very slowly which ensures that the selected vortices stay at the same positions during the temperature-variable experiments. The left and right panels of Fig. 8(b) show the vortex map (upper) and the dI/dV (r, V ) across the vortex (lower) at 0.55 K and 4.20 K respectively. Surprisingly, the spectral feature of vortex MZM disappears completely at 4.2 K, which is much lower than the gap-closing temperature of the Dirac surface states. (as shown in Section II B, the gapclosing temperature of the Dirac surface states is about 14.5 K). The early destruction of vortex MZM cannot be explained by a simplest model (e.g. Fu-Kane model). As depicted in Eq.(11), temperature modulates the intensity of vortex MZM through the superconducting gap, which basically remains unchanged below T c /2 [181]. Thus the dI/dV spectra of vortex MZM are expected to follow the thermal convolution below 8 K for FeTe 0.55 Se 0.45 single crystals. As shown in Fig. 8(c), the temperature evolution of the ZBCPs was measured at the center of a same vortex. The width of the ZBCPs is jointly determined by the tunneling coupling, instrumental resolution, thermal broadening and the poisoning. The first two issues do not change with temperature, providing that the experimental setups remain the same. If the quasiparticle poisoning is also unchanged, the dI/dV spectrum below T c /2 can be numerically reproduced by convolving the base-temperature dI/dV spectrum with the derivative of the Fermi-Dirac function. The numerical results after convolution are shown as the gray curves in Fig. 8(c). It shows that the suppression of the ZBCP exceeds the degree determined by the thermal effect at high temperatures. This observation was reproduced independently on nine vortices, indicates that additional destructive fac-tors may emerge at high temperatures, responsible to the early disappearance of vortex MZM (Fig. 8(b)).
In order to make a quantitative investigation, we define ZBCP amplitude as the conductance difference between peak and valley. A nonzero ZBCP amplitude is required for an observable MZM in experiment. Temperature evolutions of ZBCP amplitude are extracted on three vortices. As shown in Fig. 8(d), the extracted ZBCP amplitude of those vortices goes zero at about 3 K. This value is compatible to the minigap (δ) of the one-dimensional vortex [99], that δ = ∆ 2 /E F Bulk ≈ (k B T ) T = 3 K, indicates that the excessive suppression of vortex MZM is correlated to the bulk bands. The quasiparticle density increases exponentially with the temperature, and these thermally-excited quasiparticles can disturb the vortex MZM on the surface. The bulk quasiparticles are easier to be thermally excited in FeTe 0.55 Se 0.45 single crystals, owing to a smaller bulk ∆ 2 /E F comparing to the Dirac surface state. It supplies a new fermionic bath for additional quasiparticle poisoning, causing the excessive suppression of the vortex MZM with the rising temperature (the right of Fig. 8(f)). When k B T ≈ δ, the one-dimensional vortex line is somewhat conducting, and the thermally-excited bulk quasiparticles strongly mix with the vortex MZM, and pull the MZMs more inward into the bulk, or even annihilate them by meeting their partners on each end of the vortex line. In addition, at sufficient low temperatures (k B T δ), the one-dimensional vortex line become a good insulator for superconducting quasiparticles, and the vortex MZM can be approximated as emerging on a point defect of an isolated two-dimensional topological superconductor. The three-dimensionality of the material can be neglected in this situation.
Temperature-dependent measurements have been done on nine vortices. By fitting the ZBCP amplitudes with C/T (see a theory in Section V), C parameters can be extracted on each vortex, which reflects the survival temperature of MZMs. Two examples of this fitting are shown in the left panel of Fig. 8(e). The extracted C parameters were summarized in the right panel of Fig. 8(e) as a function of the amplitude at the base temperature (0.55 K). The data converge very well, demonstrating a positive correlation between the C parameters and the ZBCP amplitude at the base temperature. Owing to the inhomogeneity, basic quasiparticle poisoning is nonuniform ( Fig. 8(a)). It produces different degree of poisoning at different positions, which is manifested by different ZBCP amplitudes at the base temperature. Reducing quasiparticle poisoning is helpful to retain MZMs at higher temperatures.

C. Picture vortex MZM: two-versus three-dimensional model
Although the low-temperature behavior of MZM can be well described by the two-dimensional Fu-Kane model, the three-dimensionality of the real materials, especially the coexistence of Dirac surface states and bulk bands, does impact emergence of the vortex MZM on the sample surface. A more realistic three-dimensional model is highly required, not only for extending the Fu-Kane model, but also for the precise understanding of vortex MZMs in real materials.
After the proposal of Fu-Kane model, people soon noticed that the realistic properties of a material may complicate the fate of MZM. In most topological insulators, their bulk is a good metal rather than an in-sulator [24], which raises a question "Does the vortex MZM still exist in this situation?" [176,193]. To answer the question, a three-dimensional model was built by taking the electron-doped Bi 2 Se 2 as an example [176][177][178]. This model studied the topological phase transition of the one-dimensional vortex lines along the c-axis.
Here we summarize the evolution of dispersive bound states along the vortex line in the three-dimensional model (Fig. 9); 1) when the chemical potential is below the bulk band gap (µ < |m Γ |), the vortex line is completely insulating, and there is no subgap bound state; 2) when the chemical potential rises and reaches the edge of bulk bands (µ = |m Γ |), some dispersive vortex bound states appear along the vortex line, and there is a very small quasiparticle gap between different energy levels (δ ≈ ∆ 2 /E F Bulk ). In this case, the vortex line can be insulating only at an extremely low temperatures (lower than δ); 3) as the chemical potential continues to increase, the quasiparticle gap continues to decrease, and the Berry phase of bulk band varies accordingly; when µ = µ c , the Berry phase of bulk band reaches π, the dispersive bound states become gapless, and causes a topological vortex phase transition; 4) when µ > µ c , quasiparticle gap reopens, identifying a new topological phase.
The critical points, where the topological vortex phase transition occurs, can be directly determined by calculating the Berry phase of the bulk bands or the evolution of the bound states in a vortex line [176]. However, the identification of the topological regime, whether before or after the critical points, relies on the presence or absence of the Dirac surface state. In order to make this point clearly, we next discuss the vortex bound states for the cases of superconducting topological insulators and normal insulators, and analyze their behaviors in the views of two-and three-dimensional models. 1) when the normal state is a topological insulator. In the Fu-Kane model, a vortex MZM exists at any chemical potential, changing the chemical potential only changes the quasiparticle gap of MZMs, the vortex MZM never disappears. In the three-dimensional model, the vortex MZM can only exist in a limited energy window near the Dirac point. The vortex line is a one-dimensional topological superconductor belonging to the class-D [194][195][196][197]. When the chemical potential is located near the Dirac point, a vortex MZM appears as a boundary state of the topological superconductor. With increasing the chemical potential, a one-dimensional bound state appears in the vortex line as soon as a bulk band is involved. As the c-axis localization length of MZM is inversely proportional to δ [176], the appearance of the dispersive bound state drags the wavefunction of vortex MZM inward into the bulk. At the critical chemical potential where the topological phase transition occurs, the vortex MZM disappears. 2) when the normal state is a normal insulator. In this case, the two-dimensional model is vacuum, there is no MZM at any chemical potential. In the three-dimensional model, no MZMs exist near the zero chemical potential, but with the change of the chemical potential, a topological vortex phase transition occurs at the critical point where the Berry phase of bulk band crossing π. After the phase transition, a vortex MZM appears on the sample surface [176].
The concept of topological vortex phase transition plays a vital role on the understanding of the threedimensionality in a real material. Recently, the threedimensional model has been extended to more specific circumstances, including 1) the vortex phase transitions in Fe(Te,Se) single crystals [198]; 2) the vortex phase transitions when the topological bands are coupled with trivial bands. In this case, the topological regime deforms, but still includes the Dirac points [199]; 3) the topological vortex phase transitions induced by Zeeman coupling [200]; 4) the topological vortex phase transitions of the weak topological insulator [201] and the Dirac semimetal [202,203].
In this section, we explain in details about the mechanism of emergence of vortex MZM in a three-dimensional material Fe(Te,Se). We show that the Dirac surface state is the direct contributor to emergence of vortex MZM. At sufficient low temperatures, the behavior of vortex MZM can be well described by the Fu-Kane model, and the three-dimensionality of Fe(Te,Se) can be ignored on the context of vortex MZM. On the other hand, the quasiparticle poisoning effect damages the MZM, particularly, the bulk bands introduce additional quasiparticle poisoning at high temperatures, which accelerates the disappearance of the vortex MZM. The bulk bands also affect emergence of the vortex MZM through the topological vortex phase transition. In a material with inhomogeneity, the quantum critical point may be reached on some areas, leading to the coexistence of two topological phases in a same piece of sample. After a thorough introduction of the mechanism of emergence, the physics properties of vortex MZM will be discussed in the following sections.

V. MAJORANA NATURE OF VORTEX ZERO MODE
The self-conjugation of the creation and the annihilation operators (γ † = γ) is the essential nature of MZM. A MZM can be expressed in the form of Bogoliubov quasiparticle by the operators of complex fermion (c † = c), The wavefunctions of electron component and hole component are required to be equal (u * = v), which is a direct representation of the Majorana nature of MZM [1][2][3][4][5][6][7][8]. In this section, we introduce both the theoretical and experimental investigations of the intrinsic quantum conductance of Majorana modes (2e 2 /h). As proved in Law-Lee-Ng theory [204], the quantized Majorana conductance is a consequence of resonant Andreev reflection which is induced by the particle-hole equivalence of Majorana modes. The quantized Majorana conductance is immune to variations of the tunneling barrier under ideal conditions. In Reference [105], stable conductance plateaus were observed on the vortex zero mode of the FeTe 0.55 Se 0.45 single crystals by means of variabletunnel-coupled STS. It was found that the plateau behavior is unique to the MZM and is absent at other states, including the finite-energy CdGM states, the zero-field zero-bias conductance, and the continuum state outside the superconducting gap. Remarkably, a nearly 2e 2 /hquantized conductance plateau was observed in one vortex. These results indicate the appearance of Majorana-induced resonant Andreev reflection (MIRAR), imply the Majorana nature of vortex MZM.

A. Theory of Majorana conductance
MIRAR is a superconducting version of classical resonant tunneling [205][206][207]. A double-barrier quantum well can be constructed by band engineering in a semiconductor heterostructure (Fig. 10(a)). Quasibound states appear in the quantum well with quantization in energy. Under the condition of symmetric tunnel barrier, the transmission coefficient of electrons tunneling across the quantum well is, where τ qs , E n are the lifetime and energy of the quasibound states, respectively. Perfect transmission occurs when the incident electron is aligned with the energy level of the quasibound state. This phenomenon is similar with the resonant transmission in an optical Fabry-Perot cavity. The barriers of the double quantum well like two mirrors of the electron wave. When the electron energy is aligned with E n , the reflected electron waves are eliminated by constructive interference, leads to a resonant transmission coefficient of 1. Classical resonant tunneling can also occur between two identical tips ( Fig. 10(c)), when the tunneling amplitudes between the two tips are equal (t 1 = t 2 ). A schematic of classical resonant tunneling is shown in Fig. 10(b). It is worth noting that the relationship between the tunneling coupling strength (Γ t ) and the tunneling amplitude (t) is Γ t = 2πρ 0 |t| 2 , where ρ 0 refers to the DOS related to the coupling. Resonant Andreev reflection can be easily derived from classical resonant tunneling by replacing the outgoing electron with a hole (the electron barrier of the outgoing part is also replaced with a hole barrier) ( Fig. 10(d)) [204]. A single tip can act as both the lead of incident electron and reflective hole, which guarantees t e = t h in Andreev reflection. In order to realize the resonant condition (Γ e t = Γ h t ), the wavefunctions of the electron and hole components should be equal (ρ e = ρ h ). This requirement is naturally satisfied by Majorana symmetry (Fig. 10(e)). MZM-mediated Andreev reflection is a resonant process. Under the ideal conditions, Majorana conductance is irrelevant to the strength of tunneling coupling. On the contrary, conventional Andreev bound states (such as the subgap CdGM state) can not guarantee the particle-hole equivalence, thus the conditions of resonant Andreev reflection are not satisfied (Fig. 10(f)).
Law-Lee-Ng theory calculate the MIRAR induced quantum conductance by taking chiral Majorana modes as an example [60,61]. As shown in Fig. 10(g), a superconductor island is put on top of a topological insulator, while the area outside is covered by a magnetic insulator.
Effective spinless p x + ip y superconductivity appears under the island, outside the island a Zeeman gap opens at the Dirac points (E Z = gµ B M ), where g is the Landeg factor, µ B is the Bohr magnetic moment, and M is the effective magnetization. When E z > 2 ∆ 2 + µ 2 , chiral Majorana modes appear on the boundary of the superconducting island. The quantized energy level of the chiral Majorana modes is The corresponding quantized momentum is where m is an integer; L is the perimeter of the superconducting island; ν m is the Fermi velocity of the chiral Majorana mode; n is the number of vortices appearing in the superconducting island. In Eq. (16) and Eq. (17), the first term is derived from periodic boundary conditions; the second term is the contribution of the spin Berry phase of the Dirac surface state; the third term is derived from the vortex. The chiral Majorana mode is tunnel-coupled with a single lead at point a ( Fig. 10(g)) with the coupling amplitude of t. The total Hamiltonian of the system is

The lead term is
where ν F is the Fermi velocity of the electrons and ψ σ (x) is the fermion field. The Majorana term is The coupling term is where η(a) is the Majorana field at a, f and g are complex numbers with a modulus of 1. By a transformation on fermion field as follows, the Hamiltonian of the system can be simplified as It is apparently that the Majorana modes only couple to one spin (ψ 1 ) in the process of Andreev reflection, while the other spin (ψ 2 ) is completely decoupled [204,[208][209][210][211]. The Andreev process mediated by Majorana mode is equal-spin. The incident electrons and the reflected holes are denoted as (ψ 1,k (−)) and (ψ † 1,−k (+)), respectively. The electron-hole scattering can be described in the form of S-matrix, S-matrix can be calculated by H LLN , that Resonant Andreev tunneling (|s he | 2 = 1) occurs when θ(k, n)/2 = mπ (π is an integer). In this condition, k is just the quantized Majorana momentum (k = (2m − 1 − n)π/L). In other words, resonant Andreev reflection occurs when the incident electrons are aligned with the quantized energy level of the chiral Majorana mode. On resonance, the probability of the particle-hole reflection is 100% and irrelevant to t.
The tunneling current on resonance can be calculated by S-matrix, where Γ t = 2t 2 ν m /L is the tunneling coupling strength.
The differential conductance is The 2e 2 /h-quantized resonant conductance of Majorana modes is derived by the Law-Lee-Ng theory (Fig. 10(h)). It is worth noting that the quantum conductance is a universal property of Majorana quasiparticles, rooted on the Majorana nature of particles-holes equivalence, and irrelevant to the details used in the model, including the forms of Majorana mode (edge mode [204] or zero mode [212][213][214]), the theoretical technique (S-matrix [204,214] or Green function [212]), and the coupling strength (single particle tunneling [204,212] or quantum ballistic transport [213]). Therefore, it is a direct manifestation of essential properties of Majorana symmetry. Law-Lee-Ng theory investigates the Majorana conductance in a dissipationless system at absolute zero temperature. With those ideal conditions, the zero-bias conductance of MZM should be at the quantum value (2e 2 /h). A quantized plateau should appear in a variable-tunnelcoupled measurement of Majorana modes. However, the conductance of Majorana modes is much smaller than 2e 2 /h in many experiments [36], indicates the influence of imperfect conditions on the conductance behavior of Majorana modes [212][213][214][215][216][217][218]. In the following, we consider quasiparticle poisoning effect and thermal effect, to investigate their influence on the behavior of the zero-bias conductance.
First, we analysis the Majorana conductance in a dissipationless system at a finite temperature. The temperature evolution of Majorana conductance is, The zero-bias conductance of MZM satisfies a k B T /Γ t scaling behavior (Fig. 10(i)) [215][216][217], that is G s = is the scaling function. In a dissipationless system the quantized conductance can be observed only when the tunneling coupling is much larger than the thermal broadening.
Next, we analysis the Majorana conductance in a dissipative system at a finite temperature. At the zero kelvin, the conductance of MZM with quasiparticle poisoning is where Γ p = 2πρ p |t p | 2 is the coupling strength between the fermionic bath and the MZM [186]. The quasiparticle poisoning effect weakens the zero-bias conductance of MZM by a fraction of Γ t /(Γ t + Γ p ). Similarly, the MZM conductance at a finite temperature can be obtained by convolving the derivative of Fermi-Dirac function (df FD /dE) with the zero-temperature conductance. In a dissipative system, the zero-bias conductance no longer follows the k B T /Γ t scaling behavior. We calculate the evolution of zero-bias conductance under different dissipation strength (inset of Fig. 10(i)). Dissipation causes the decrease in G s , but G s still approaches to the quantum conductance when Γ t Γ p .
In summary, the width of MZM is determined by the temperature broadening (3.5k B T ), the dissipation and the tunneling coupling broadening (2(Γ p + Γ t )). A quantized Majorana conductance plateau can be experimentally observed only when Γ t is much larger than k B T and Γ p .

B. Conductance plateau of vortex zero modes
The variable-tunnel-coupled experiment takes the advantages of the feedback regulation of STM, which controls the tip height by the setpoint of tunneling current (I t ) and bias voltage (V s ). It can be used to change the tunneling coupling continuously in an experiment. The conductance of the tunnel barrier (G N ≡ I t /V s ) is regarded as a measure of tunneling coupling (Γ t ) between the tip and the sample (G N is positively correlated with Γ t ). Quantized Majorana conductance can be hopefully explored by achieving a large value of G N , in which Γ t is much larger than k B T and Γ p .
By positioning the STM tip directly above the vortex MZM in a FeTe 0.55 Se 0.45 single crystal ( Fig. 11(a)), the evolution of the zero-bias conductance of MZM was measured continuously with increasing the tunneling coupling strength [105]. As shown in Fig. 11(b), the zerobias conductance of MZM tends to be saturated when the tunneling coupling is strong enough (G N ≈ 0.3G 0 , G 0 = 2e 2 /h). As the tunneling coupling is further increased, the zero-bias conductance of MZM shows a plateau behavior, while the conductance of the electron continuum outside the superconducting gap keeps increasing. It implies that tunneling coupling between the tip and the vortex MZM may be unconventional. In order to exclude other possible trivial mechanisms of the zero-bias conductance plateau, the conductance behavior of the CdGM bound states at finite energies ( Fig. 11(h)), the continuum outside the superconducting gap under the zero field (Fig. 11(i)), and the zero-field zero-bias conductance ( Fig. 11(j)) were measured with the change of tunneling coupling strength. Those measurement show that no conductance plateau appears on those cases. The conductance measurement under the zero field eliminates the possibility of quantum ballistic transport [219][220][221][222][223][224]. In addition, the wavefunction of CdGM states is electron-hole inequivalence (Fig. 10(f)), which is at odds with the requirements of the resonant Andreev reflection. A repeatable experiment show that the conductance plateau behavior is unique to the vortex MZM.
In Ref. [105], variable-tunnel-coupled STS measurements were performed on 60 topological vortices, among which 29 vortices became unstable in the process of approaching the tip, possibly owing to vortex creep. On the other 31 vortices, the zero-bias conductance of the MZM shows a plateau behavior after the tunneling coupling exceeding a threshold. The plateau conductance (G P ) of the 31 vortex MZMs are summarized in Fig. 11(g). It is obviously that in most cases the zero-bias conduc- tance reaches a plateau at a non-universal value (0.2G 0 to 1.0G 0 ), the central value of the G P distribution is about 0.6G 0 . A typical case of non-quantized conductance plateau is shown in Fig. 11(c) and Fig. 11(d). However, among the 31 measurements, in one vortex, the zero-bias conductance saturates on a nearly quantized value when the barrier conductance G N = 0.7G 0 (Fig. 11(e) and Fig. 11(f)). It is consistent with the theoretical expectation of MIRAR, and may present the Majorana nature of MZM.
The non-universal plateau conductance is a common behavior of vortex MZMs in FeTe 0.55 Se 0.45 single crystals. Further experiments [105] found that, artificially increasing the instrumental broadening (increasing the lock-in excitation voltage) can reduce the plateau conductance of vortex MZM. Furthermore, by comparing the behavior of vortex MZM at different positions, it was found that the basic quasiparticle poisoning suppresses the plateau conductance. We note that the nonuniversal plateau conductance cannot be described by the theory explained in Section V A. Although the conductance plateau behavior is unique to the vortex MZM, which strongly implies the appearance of MIRAR, there is no direct experimental evidence requiring that the nonuniversal conductance plateau must be induced by Majorana modes. Further theoretical and experimental stud-ies are needed to reach a fully understanding of the nonuniversal plateau conductance.
The discovery of zero-bias conductance plateau of vortex MZMs, especially the nearly quantized case (Fig. 11(e) and Fig. 11(f)), indicates a direct measurement of the Majorana nature on vortex MZM, which supports a Majorana origin of the observed ZBCP in the vortex of FeTe 0.55 Se 0.45 single crystals.

VI. TOPOLOGICAL NATURE OF VORTEX ZERO MODE
MZM is characterized by non-trivial topological invariants [194][195][196][197]. For intrinsic topological superconductors (e.g. p x + ip y superconductor [16] and p-wave Kitaev chains [225]), the topological properties of their quasiparticle spectra can be demonstrated by the mapping from S 2 (k-space) to S 2 (spinor space). For connate topological superconductors (e.g. Fe(Te,Se) single crystals), the vortex MZM can be regarded as the end mode of an one-dimensional topological superconductor (the vortex line) in the view of three-dimensional model, or as the bound state on the topological defects (the vortex) in the view of two-dimensional model (Section IV C). However, emergence of MZM in a vortex requires neither the topologically non-trivial superconductivity, nor the topologically non-trivial band structure [226], but needs to integrate all the ingredients, i.e. band topology, superconductivity and vorticity, to reach a non-trivial global topological invariant. For a system with N Fermi surfaces, the topological invariant of a vortex MZM can be presented as [227,228], where w i is the winding of the quasiparticle spectrum on the i-th Fermi surface, m i is the vorticity of the i-th Fermi surface. When Z 2 = 1, a vortex MZM can emerge. Note that Z 2 is an emergent topological invariant, which indicates the opportunity for emergence of a vortex MZM in trivial materials [226,227,229]. For instance, the fractional vortex of multi-component superconductors (e.g. spin triplet superconductors, pair density wave, and nematic superconductors) can satisfy the requirements of emergence of vortex MZM [227,230].
Since a vortex has an intrinsic winding number, the nontrivial Z 2 in the Fu-Kane hamiltonian is mainly attributed to the appearance of Dirac surface states [58]. It inspired an experimental method for identifying the topological nature of vortex MZM by investigating the influence of Dirac surface states on the behaviors of vortex bound states. Dirac surface states possess an intrinsic spin angular momentum, which leads to a half-integer level shift on the vortex bound states level sequence. It not only produces a zero-energy Majorana mode, but concomitantly changes the DOS spatial pattern sequence of vortex bound states. These characteristics directly reflect the topological nature of vortex MZM. In this section, we discuss the behaviors of vortex bound states with an underlying superconducting Dirac surface state, and introduce the experimental observations of the topological nature of vortex MZM.

A. Vortex bound state with topology
In conventional s-wave superconductors, the total angular momenta (ν) of vortex bound states are halfintegers (ν = ±1/2, ±3/2, ±5/2···) owing to the underlying parabolic band structure. When the level spacing is not too large, the vortex bound states follow the half-integer level sequence (Fig. 12(a)), i.e., E ν ≈ ν∆ 2 /E F , with no zero-energy modes. It is a common behavior of conventional s-wave superconductors. However, a pronounced ZBCP was observed at the vortex center of a conventional superconductors NbSe 2 in an early STM experiment. When moving away from the vortex center, the ZBCP splits into two symmetrical peaks, displaying a spatially-dispersive distribution ( Fig. 12(b)) [14,175,231]. These experimental observations are apparently at odds with the theoretical prediction which requires non-zero and discrete solutions in energies. However, in nodal superconductors (such as d-wave superconductors), the subgap quasiparticles mix with the continuum outside the gap, leading to a spatially-dispersive distribution of the vortex DOS [232][233][234][235]. Although it seemly matches the experimental results, it is apparently inconsistent with the fully-gapped superconductivity of NbSe 2 .
The anomalous behavior of vortex bound states is caused by the realistic condition that the experimental temperature fails to reach the quantum limit of vortex bound states [15]. When the thermal broadening is smaller than the level spacing of the vortex bound states [236], the quantum limit condition is reached. The threshold temperature of quantum limit (T QL ) is estimated as, T QL = T c ∆/E F . When the experimental temperature exceeds T QL , thermal broadening dominates the behavior of vortex bound states, so their energy levels overlap with each other, leading to the spatiallydispersive distribution. Discrete vortex bound states can be only observed in the quantum limit. The realization of quantum limit is a prerequisite for investigating the influence of Dirac surface states on the behaviors of vortex bound states. Conventional superconductors usually have a large Fermi energy (E F ≈ 2 to 10 eV) and a small superconducting gap (∆ ≈ 1 meV). It leads to an extremely small T QL , which is difficult to be reached in experiments. The realization of quantum limit prefers a material with a high critical temperature and a small Fermi energy. As discussed in Section I, several compounds of iron-based superconductors naturally meet these requirements ( Fig. 1(a)) [237][238][239][240][241][242].
A superconducting vortex is a magnetic field induced There are only parabolic bulk bands involved [104]; (b) the quantum limit is difficult to reach in conventional s-wave superconductors so that a large ZBCP observed in the center of vortex core is generally due to multiple overlapping of densely packed non-zero peaks [231]; (c) integer quantized level sequences of the bound state in Fu-Kane model. The intrinsic spin Berry phase carried by the Dirac surface states induces the half-integer level shift [104]; (d) the zero-doping limit is defined as the situation that the chemical potential approaches the energy of the Dirac point. In this case, a vortex MZM is the only allowed subgap bound state [104]; (e) theoretically calculated angular momentum resolved wavefunction of the BdG eigenstate, the blue and green curves are spin-down and spin-up components, respectively [210]. Inset: calculated spin-integrated 2D local density of states of three lowest levels of the bound states in the cases of (c) and (a), respectively [104]; (f) theoretically calculated eigenvalue of BdG Hamiltonian near the zero chemical potential limit.
topological defect in a type-II superconductor. The amplitude of the superconducting gap gradually decreases to zero along the radial direction towards the vortex center. The phase of superconducting gap winds around the vortex which is accordant with the supercurrent circulating at the boundary (London penetration depth). Generally, a superconducting vortex can be written as where ∆ 0 is the gap amplitude at the region far away from the vortex; h(r) is the radial distribution function of the gap amplitude, it can be written as h(r) = tanh(r/ξ) or h(r) = r/ r 2 + ξ 2 ; m, the vorticity, which is an integer and represents how many turns that the superconducting phase changes as it goes around a vor-tex per circle. φ(x, y) is the spatial distribution of the superconducting phase. For example, a single vortex can be expressed as φ(x, y) = arctan(y/x), and the vortex-antivortex pair can be expressed as φ(x, y) = arctan[2ay/(x 2 + y 2 − a 2 )], in which the vortex and the antivortex locate at (a, 0) and (−a, 0), respectively [243]. Note that the definition of the vortex and the antivortex depends on the positive and negative value of the vortex winding numbers, respectively. The vortex winding number is Z v = c arg[∆(r)]ds, where c is a closed loop around the vortex center [244].
The vortex bound states induced by a superconducting Dirac surface state can be derived by incorporating the surface Dirac cone in the BdG equation with a single vortex [58,205,208], Due to the rotational symmetry of the vortex, the vertical component (K z ) of the total angular momentum is a good quantum number, namely [K z , H BdG ] = 0, the energy eigenvalues of vortex bound states can be expressed by the eigenvalues (ν) of K z . When ∆/E F is not too large, the level sequence of vortex bound states can be approximately expressed as where m is the vorticity; sgn(µ) is the sign of the chemical potential of the Dirac surface states; ν is the vertical component of the total angular momentum. It has been proved that with the incorporation of Dirac electrons, the angular momentum is, l z is the orbital angular momentum which is an integer. S z is the spin angular momentum, which is +1 (-1) for the spin-up (down) component [18,210,211]. Different from the case of conventional s-wave superconductors, the angular momentum (ν) of the vortex bound states shown in Eq.(34) can be any integers (ν = 0, ±1, ±2, ±3, ). Accordingly, the energy levels of vortex bound states inherit the integer quantization, in which ν = 0 level is the vortex MZM ( Fig. 12(c)). This exotic behavior is caused by the additional S z introduced by the Dirac surface states [104]. As a contrast, the angular momentum in a conventional s-wave superconductor is This Dirac-surface-state-induced half-integer level shift of vortex bound states is a direct expression of the topological nature of the vortex MZM. The half-integer level shift of the vortex bound states is accompanied by a change of the DOS pattern sequence on each level. As we introduced in Section IV A, the theoretical wavefunction of a vortex MZM is proportional to Bessel function J i (x), of which the order is determined by the orbital angular momentum of vortex bound states, namely i = l z . According to the properties of Bessel function, the maximum of the l z = 0 component appears at the vortex center. With increasing l z , the DOS maximum gradually moves away from the vortex center, leading to a low-intensity valley. Combining with the analysis shown in Eq.(33)-Eq. (35), we obtain the behaviors of angular momentum (ν, l z , S z ) and wavefunction on each level. In Fig. 12(e), we show a theoretical simulation of the DOS distribution of vortex bound states under the conditions of m = −1 and µ > 0, where the purple (green) curves represent the spin-down (up) component [210]. By summing up all the components on each energy level, we obtain the DOS pattern sequence of vortex bound states, which is observable by the constant-energy conductance map of a STM measurement. It can be found that the lowest energy level (ν = 0) and one of the two second lowest levels (ν = +1 or -1) are of "solid circle pattern" while the rest are of "hollow ring patter" (the top row of Fig. 12(e)). On the contrary, in the case of conventional s-wave superconductors, only one of the two lowest energy levels (ν= +1/2 or -1/2) has "solid circle pattern" (the bottom row of Fig. 12(e)). This difference in the pattern sequence is an important feature of the half-integer level shift induced by the Dirac surface state [104].
By implementing the wavefunction and angular momentum analysis, we derived some conclusions of the influence of Dirac surface states on the vortex bound states, as listed below: 1) reversing the direction of magnetic field changes the sign of vorticity, consequently changing the positive and negative correspondence between ν and E ν , but it does not change the sign of the energy of the second solid-circle-pattern level [106]; 2) the sign of the chemical potential determines the sign of the energy of the second solid-circle-pattern level. Specifically, the second solid-circle-pattern level appears in the unoccupied side (positive energy) when µ > 0, and vice versa; 3) at the vortex center, the spin of the zero-energy level is always parallel to the magnetic field, while the spin of the second solid-circle-pattern level is always anti-parallel to the magnetic field. It is the spin-resolved property of the vortex MZM [66,[208][209][210][211].
In the last part of this section, we focus on the behavior of vortex bound states under the zero-doping limit. As mentioned previously, the integer quantization is robust in term of the angular momentum, but not robust in term of energy. When ∆/E F is very large, the energy of the second lowest vortex bound states will be very close to the superconducting gap edge, and this would induce quantum confinement on the higher levels that locate between the gap edge and the second lowest level. Thus the energies of the vortex bound states are no longer integer quantized.
Theoretically, when the chemical potential approaches to the Dirac point (E F → 0), which is referred as the zero-doping limit, the MZM is the only allowed subgap bound state (Fig. 12(d)) [245,246]. The energy eigenvalue of the BdG equation has analytical solutions under some special circumstances. For example, when h(r) = r/ r 2 + ξ 2 , the energies of vortex bound states in a single vortex (|m| = 1) can be derived as, [247]. By using Eq.(36), we calculate the level spectrum of vortex bound states when m = −1 and µ → 0 + as shown in Fig. 12(f), finding that only an isolated vortex MZM exists within the superconducting gap, while other non-zero vortex bound states are pushed to the gap edge. The topological gap of vortex MZM is the largest under the zero-doping limit.
It is worth noting that when µ is exact zero, an additional pseudo-chiral symmetry appears in the Fu-Kane hamiltonian, it changes the topological classification of vortices from class-D to class-BDI, and the topological invariant from Z 2 to Z [196]. Majorana hybridization of the vortex MZM is forbidden in class-BDI [248][249][250], resulting in failure of Majorana topological qubits. However, the pseudo-chiral symmetry guarantees the degeneracy of multiple vortex MZMs. A Majorana flat band appears in a lattice formed by the vortex MZMs, and the four-body Majorana interactions play an important role in the Majorana flat bands [251,252]. It provides a rare opportunity to generate novel quantum phenomena [253], including SYK model and Majorana fractional quantum Hall effect. The MZM interaction is beyond our scope here. The readers interested in this topic may refer to Ref. [253].

B. Observation of integer-quantized topological vortex
The realization of quantized vortex bound states in Fe(Te,Se) [239] enables the measurements of Diracsurface-state-induced integer quantization, which reflects the topological nature of vortex MZM. In Section IV and Section V, we only focus on the properties of the vortex MZM located at the zero energy, but behaviors of other vortex bound states that accompany the presence of the MZM were missed previously [102,105]. With a deep understanding of the topological nature of vortex MZM, the characteristics of the global level sequence of vortex bound states have been further studied. In this section, we introduce the relevant experimental results.
On the basis of previous work, the vortices of Fe(Te,Se) were further studied in experiments. By paying more attention to the global behavior of vortex bound states rather than focusing only on the zero modes, we found that both the vortex MZM and other finite-energy vortex bound states exist in the topological vortices. At extremely low temperatures, they all show discrete features, indicating realization of the quantum limit in FeTe 0.55 Se 0.45 single crystals ( Fig. 13(a)). The energy/spatial positions of the vortex bound states (vor-tex1) were extracted in Fig. 13(b) which clearly shows the level quantization. In Fig. 13(c), the energy of each levels (E L ) was normalized by the level spacing (∆E 1 ). An integer-quantized level sequence can be clearly observed. The energy eigenvalues of the vortex bound states can be numerically calculated as a function of angular momentum, which was performed on vortex1, and the results are fully consistent with experiment ( Fig. 13(c)). Furthermore, a statistical analysis was done among all the topological vortices measured. Figure 13(i) is a histogram that shows the distribution of normalized level energies of 35 vortices, and integer quantization was further supported by the statistical analysis.
The integer-quantized level sequence is attributed to an additional half-integer angular momentum introduced by the Dirac surface states (Section VI A). In addition to the integer quantization, the Dirac surface states also produce double "solid circle pattern" levels in the DOS spectra of the vortex bound states. In order to ensure that the integer level sequence observed in Fig. 13(a)- Fig. 13(c) does come from the underlying Dirac surface state, rather than a coincidence, the constant energy conductance maps were measured experimentally at the energy of the 0-, +1-and +2-levels of the vortex bound states (Fig. 13(j)). It shows that the first two levels have solid-circle DOS pattern, and the +2-level is hollowring-like. Considering the positive chemical potential of the Dirac surface states in FeTe 0.55 Se 0.45 , the second "solid-circle" level should appear on the unoccupied side (positive bias) [205]. The experimental results are fully consistent with the theoretical expectations ( Fig. 12(e)). The observation of the double "solid circle" spatial pattern indicates that the observed integer-quantized bound states in FeTe 0.55 Se 0.45 single crystals are emerged from the superconducting Dirac surface states. We note that the clear observation of the hollow-ring-like pattern at ν = +2 level is attributed to the tiny k F (about 0.02Å) of the Dirac surface state (Fig. 3). The ring radium of the levels with l z = 0 is proportional to 1/k F (Fig. 12(e)). A smaller k F value results in a larger spatial oscillation period, which is easier to be observed on the constant energy conductance map of STM. However, k F of the bulk states of Fe(Te,Se) is about 0.1Å, corresponding to a small spatial oscillation period in the order of 1 nm. It is difficult to survive under spatial inhomogeneity.
In some vortices, a prominent MZM was isolated in the middle of the superconducting gap, there seemly no other bound states except the MZM, obviously violates the integer quantized level sequence. In Fig. 13(d)- Fig. 13(f), we reanalyze a topological vortex belong to this case (vor-tex11 displayed in Fig. 7). It can be found that the prominent MZM is not the only subgap state, there are three discrete vortex bound states near the gap edge. It is obvious that the vortex bound states do not conform to the integer quantization ( Fig. 13(f)). A numerical simulation, following the same model used in Fig. 13(c) but with a smaller chemical potential, fully reproduced experimental results, indicates vortex11 is near the zero doping limit. Furthermore, the zero-bias conductance of MZM was extracted on the three vortices ( Fig. 13(g)). It shows that the spatial extension of the MZM near the zero-doping limit is wider. We calculated the analytical Majorana wavefunction (|u| 2 (r)) under different E F value (Fig. 13(h)), and found that smaller E F values correspond to larger spatial distribution of vortex MZM [58,210]. It is fully consistent with the experi- (c) comparison between experimentally observed and theoretically calculated level energy in topological vortex1; (d) same as (a), but measured on vortex11, which is close to the zero chemical potential limit; (e) overlapping display of dI/dV spectra selected from (d); (f) same as (c), but shows the case of vortex11; (g) comparison of observed MZM line profile in topological vortices under integer quantization (open circles) and near the zero chemical potential limit (dark stars); (h) calculated MZM wavefunction under different chemical potential by Fu-Kane model; (i) histogram of averaged level energies normalized by the first level spacing, i.e., the ratio EL/∆E. The statistical analysis is performed among all the 35 topological vortices which show integer quantized CdGMs levels; (j) experimentally observed spatial pattern of the lowest three levels of bound state in a topological vortex [104]. ments, indicates that the break down of integer quantization in vortex11 is a result of near zero-doping limit.

C. Inhomogeneity enables ordinary vortex
Even though topology dictates the existence of two types of discrete bound state spectra, i.e., integer quantization in a topological vortex and half-integer quantization in an ordinary vortex, in ordinary circumstances a given material belongs to just one of the classes. This restricts a single sample to one type of spectrum and thus forbids a direct comparison. However, the intrinsic inhomogeneity of FeTe 0.55 Se 0.45 single crystals [255][256][257][258][259][260][261][262], while playing a destructive role on emergence of vortex MZM (Section IV), enriches the sample properties, so that a direct comparison between topological and ordinary vortices can be realized in a same region, provides a rare opportunity to identify the half-integer level shift.
Intrinsic inhomogeneity is identified in FeTe 0.55 Se 0.45 single crystals on multiple aspects, including chemical composition, chemical potential, disorder/scattering potential, superconducting gap, and etc. Some influences of them have been demonstrated in the previous sections, that 1) spatial fluctuations of the superconducting gap and the chemical potential can induce the near-zerodoping limit behavior in some topological vortices (Section VI B), 2) stronger spatial fluctuations of the chem-FIG. 14. Inhomogeneity of a material enables coexistence of ordinary and topological vortices. (a) dI/dV (r, V ) line-cut intensity plot measured on ordinary vortex8. Half-odd-integer quantized bound states are clearly observed [104]; (b) comparison between experimentally observed and theoretically calculated level energy in ordinary vortex8 [104]; (c) histogram of averaged level energies normalized by the first level spacing, i.e., the ratio EL/∆E. The statistical analysis is performed among all the 26 ordinary vortices which show the half-odd-integer quantized CdGM levels [104]; (d) surface disorder transforms a strong topological insulator to a normal insulator. The scattering potentials gradually increase from left to right [254]; (e) concentration of the dopants could drive a strong topological insulator to be a normal insulator or a weak topological insulator in Fe(Te,Se). The bands in green (red) represent pz (dxz/dyz) orbital with odd (even) parity [201].
ical potential can drive a topological phase transition of the vortex line (Section IV C), 3) inhomogeneous disorder/scattering potential may contribute to nonuniform quasiparticle poisoning, induces different zero-bias conductance of MZM on different vortices (Section IV B). In those cases, inhomogeneity directly affects the behavior of low-energy superconducting quasiparticles, and we call them "weak inhomogeneity" in this review. In contrast, the "strong inhomogeneity" not only affects lowlying quasiparticles directly, but also indirectly modify the underlying band structure.
In a FeTe 0.55 Se 0.45 single crystal, strong inhomogeneity can destroy strong topological insulator states in some portions. A single piece of sample can be regarded as a topological crystal embedded by many non-topological nanocrystals. The Dirac surface state is not fully seated on the bare surface of the sample, but wiggles along the new boundaries between topological and non-topological crystals ( Fig. 16(a)). It leads to the disappearance of the Dirac surface state on some non-topological regions where the half-integer-quantized ordinary vortices can appear. Experimental studies found that topological or ordinary vortices usually appear in groups, It supports the picture of disappearance of the Dirac surface states on some surface regions of FeTe 0.55 Se 0.45 (Fig. 15). Two possible mechanisms of "strong inhomogeneity" which can locally destroy the strong topological insulator states are described below.
1) disorder induced non-magnetic scattering drives the transition to a normal insulator state. Generally, the strong topological insulator states are protected by the time-reversal symmetry, and it cannot be destroyed by any non-magnetic impurities [24]. However, the topological protection is only valid under the condition of weak impurity scattering. When the non-magnetic scattering potential is comparable to the inverted SOC gap of the bulk bands, the strong topological insulator state can be destroyed. A theoretical calculation confirms the disappearance of the Dirac surface state under the condition of strong scattering potential (Fig. 14(d)). The SOC gap of the bulk bands in FeTe 0.55 Se 0.45 is estimated to be about 20 meV (Fig. 4(f)), which is much smaller than that of the classical topological insulator Bi 2 Se 3 (about 300 meV). Although the topological protection of the Dirac surface state is robust and universal within the topological phase, the topological phase itself is not robust owing to the small SOC gap of FeTe 0.55 Se 0.45 . The strong inhomogeneity of the material breaks the topological band structure at some regions, in which ordinary vortices appear [254,263,264].
2) The composition fluctuation of the anionic atomic drives the transition to normal insulator or weak topological insulator state [265,266] (Fig. 14(e)). As shown in Section II A, there is no topological bands in Se-rich Fe(Te,Se) single crystals. With proper doping concentration of Te atoms, a topological band inversion occurs at Z, the material enters a strong topological insulator state. However, it is worth noting that in an overdopped sample, the topological band inversion also occurs at Γ, so the material becomes a weak topological insulator, and the Dirac surface state only exists on the side surface [265,266]. Therefore, the overdoping of both Se and Te atoms may eliminate the Dirac surface states on the (001) surface, resulting in the appearance of ordinary vortex in some regions [201,267,268].
Ordinary vortex was observed on the surface of FeTe 0.55 Se 0.45 single crystals (Fig. 14(a)-Fig. 14(c)). After implementing a standard analysis (used in Section VI B) on 26 ordinary vortices, we identified a good half-integer level sequence of the vortex bound states. This result can also be fully reproduced by numerical simulation. Owing to the strong inhomogeneity of FeTe 0.55 Se 0.45 , these ordinary vortices emerge at some regions on the (001) surface without Dirac electrons, coexisting with topological vortices in a same piece of sample. A direct comparison between the two classes of vortices under different topologies displays the half-integer level shift of vortex bound state ( Fig. 12(a) & (c)), and reveals the topological nature of vortex MZM.

VII. FROM QUANTUM PHYSICS TO QUANTUM ENGINEERING
In a real material, the three-dimensionality and intrinsic inhomogeneity induce variable behaviors and complex distributions of vortices, including two distinct classes of vortices and disappearance of MZM in some vortices. As every coin has two sides, the imperfection of samples plays different roles depending on the topics concerned. On the perspective of physics study, it is a gift of nature, which enables a direct comparison of the vortex bound states under two different topologies. But from the view of practical applications, it becomes a nightmare for the exploration of topological quantum computation. In the previous sections, we describe in details the basic properties and experimental observations of the vortex MZM in iron-based superconductors. In combination with the simplified theoretical models, that only taking the perfect cases into account, the analysis always captures the key nature of exotic phenomena (Section IV). This is the methods of physics research. However, real materials are not as simple as the theoretical models, the imperfections in materials and environments introduce other factors, which interfere with the target physical phenomena. Those problems must be solved in developing practical applications such as topological quantum computation. This is the task of quantum engineering. Therefore, in order to realize the ultimate goal of topological quantum computation, it is of profound significance to experimentally study the morphology of vortex MZM on the sample surface and theoretically explore the physics mechanisms that affect presence or absence of vortex MZM.
Here we introduce the distribution characteristics of vortex MZMs on FeTe 0.55 Se 0.45 [104]. Early experiments were performed at a higher temperature (about 450 mK), under such a condition, the percentage of the vortices that host MZM is less than 20%, which is difficult to conduct statistical studies [102]. The investigation of quasiparticle poisoning (Section IV B) shows that the signal of vortex MZMs is stronger at lower temperatures. It stimulated new experiments conducted under more extreme conditions. As shown in Fig. 15, the vortex morphology was studied on three regions which were randomly selected and far away from each other. All the vortices appearing on the three regions were carefully studied by spatial dependent dI/dV measurements of the vortex bound states (in total 76 vortices were measured). By performing standard data analysis as used in Fig. 13 and Fig. 14, their types were identified on each vortex. As marked in Figs. 15(a), (c), (e), the yellow circles represent the topological vortices with vortex MZM, the blue circles represent the ordinary vortices without vortex MZM. The statistical results of the three regions are listed in Figs. 15(b), (d), (f). We derived three conclusions from these data: 1) the occurrence probability of vortex MZM has large spatial fluctuation. It shows that the probability varies from 37% to 75% depending on the regions. Since the three regions are randomly selected, the large-scale inhomogeneity should play an important role on the occurrence probability of vortex MZM, thus emphasizes the importance of sample quality. 2) most of the vortices (about 76%, among all three regions) show observable quantization sequences of vortex bound states, either integer quantization in topological vortices or half-odd-integer quantization in ordinary vortices. 3) ordinary or topological vortices appear in groups. In Fig. 15, the vortices within the same class are encircled by green dotted lines. It supports the picture that the Dirac surface states disappear on some surface areas, while remaining intact in the others.
The disappearance of MZM in some vortices arose con- troversies on the field of iron-based Majorana platform. In order to understand these phenomena, here we systematically summarize the possible microscopic mechanisms which can influence the presence or absence of vortex MZM [104]. First of all, strong inhomogeneity eliminates the Dirac surface states on some areas of the (001) surface. Topological regions (with Dirac surface states) and conventional regions (without Dirac surface states) can coexist on sample surface, in which topological vortices (with vortex MZM) and ordinary vortices (without vortex MZM) appear, respectively.
Second, the following effects affect presence or absence of vortex MZM at a given spatial position.
1) vortex topological phase transition changes the topological invariant [176] (see a detailed discussion of the influence of topological phase transition of vortex lines in Section IV C). For the vortices in a topological region, the occurrence of topological phase transition destroys the vortex MZM, while in a conventional region, the phase transition creates the vortex MZM.
2) low-lying quasiparticles suppress the vortex MZM. First, the vortex MZM is protected by the bulk mini gap (δ) of the one-dimensional vortex lines. The vertical decay length of the Majorana wavefunction is proportional to 1/δ. On approaching to the critical point of quantum phase transition, the Majorana wavefunction goes deeply into the bulk of the material, thus become invisible in STM measurementa. On the quantum critical point, the one-dimensional vortex line becomes gapless. The two MZMs on the top and bottom surfaces annihilate with each other. Second, the quasiparticle poisoning is not uniform owing to the intrinsic inhomogeneity of the material. It results in different intensities of vortex MZM measured at the same experimental temperature. Third, additional quasiparticle poisoning exists at higher temperatures (Section IV B), which induces the rapid smearing of vortex MZMs with raising temperatures.
These mechanisms are summarized into a phase diagram of vortex MZM when the magnetic field is week (Fig. 16), in which the red dot line is an indication of the phase region which can be covered by current experiments.
Finally, it is worth pointing out that under high magnetic fields the in-plane Majorana hybridization plays a vital role on controlling presence or absence of a vortex MZM. It was observed that the occurrence probability of the vortex MZM decreases with the increasing magnetic field [269,270], this phenomenon was later explained by a theoretical simulation of Majorana hybridization in a disordered vortex lattice [271]. In the conventional regions, the corresponding bulk states can be normal insulators or weak topological insulators. Consequently, the Dirac surface state moves deeper into the bulk and goes around the conventional region, as indicated by the gray boundary inside the crystal. In other topological regions (gray color), where the Dirac surface states remain intact, the corresponding bulk states are still in the strong topological insulating phase; (b) schematic phase diagram of vortex MZMs appearing in topological regions (topological vortices). The gradient blue areas in (b) and (c) indicate the phase sector that MZMs can be detected by STM/S experiments. In the dark blue sector, the Majorana wavefunction is more localized on the sample surface, while in brighter positions, the Majorana wavefunction strongly hybridizes with the bulk quasiparticles and moves deeper beneath the surface, leading to a weak ZBCP signal measured by STM/S. The vertical axis demonstrates the evolution of MZMs as a function of effective temperature which can be represented by the extrinsic broadening of observed ZBCPs. The horizontal axis demonstrates the MZMs evolution as a function of quantum parameters, e.g., the chemical potential (µ) measured from the Dirac point. The black dots with an arrow indicate the quantum critical points in which a vortex phase transition happens. Across the critical point, the vortex line turns to be topologically trivial and MZMs disappear in the topological region. The red dashed line indicates the achievable region in experiments; (c) a schematic phase diagram of vortex MZMs appearing in conventional regions (ordinary vortices). There are no MZMs in our measurements in those vortices. The observable MZMs can only exist above the critical points when the vortex phase transition turns the trivial vortex line into a 1D topological superconductor in the conventional region [104].

VIII. CONCLUSION AND OUTLOOK
In this article, we have present a systematic and comprehensive review on vortex MZMs in Fe(Te,Se) single crystals, ranging from the origin idea to research progresses, from classical theories to new experiments, from band structure to quasiparticles, and from fundamental physics to realistic details. We aim to bridge the gap between the well-established Majorana theories and the emerging "iron home" of Majoranas, to help the readers to thoroughly understand and reasonably evaluate the emergent vortex MZMs in the iron-based Majorana platform.
Since Fe(Te,Se) was first discovered as a carrier of Majorana modes, the topological properties of iron-based superconductors soon become a hot topic in condensed matter physics. Over the past few years, a large number of theories and experiments have emerged. Here we try to summarize them as follows, 1) independent verifications of appearance of vortex MZM in Fe(Te,Se) single crystals [269][270][271][272]. 2) new developments in theories of the vortex topological phase transition in iron-based superconductors [200][201][202][203]. 3) discovery of more topological compounds of iron-based superconductors. Topological band structure has been found to be universal in iron-based superconductors [101]. Vortex MZMs were observed in (Li,Fe)OHFeSe [273,274] and CaKFe 4 As 4 [106]; 4) exploration of new compounds of iron-based superconductors which may support vortex MZM at higher temperatures. Experimental studies identified a topological band structure in a high-T c superconductor Fe(Te,Se) monolayer [267,268], indicating that the monolayer of iron-based superconductor may be a high-temperature Majorana platform above the liquid helium temperature [75]. 5) developments of defect states, both in theory and experiment. Evidence of the Majorana mode has been reported on point-like impurities [99,[275][276][277][278], step edges [279], domain walls [272], and atomic line defects [280][281][282]. 6) evidence of time-reversal symmetry breaking [283][284][285]. 7) possibilities of intrinsic topological superconductivity [286,287]. 8) studies of heterostructures that combine a Fe(Te,Se) single crystal with a topological material [288][289][290][291][292]. 9) improvement of sample quality [293]. 10) design of Majorana Kramer pairs which preserve the time-reversal symmetry. Majorana Kramer pairs are predicted to appear on the boundary of some iron-based superconductors, providing the s± superconducting pairing applied on the topological surface state [294,295]. 11) design of high-order MZMs [296][297][298][299]. 12) exploration of Majorana braiding [300][301][302]. 13) Majorana research on other connate topological superconductors [303][304][305]. The booming developments of this field not only open up a new horizon for Majorana physics but also strongly support the research on high-temperature superconductivity.
The vortex of iron-based superconductors is emerging as one of the most reliable platform for Majorana zero modes. On the basis of deep understanding of their properties, hybridizing, braiding, fusing MZMs and reading out the quantum information of Majorana qubits become important research directions in the future. The realization of these goals requires joint efforts on the aspects of theory, material, and technique. A practical braiding strategy should be theoretically designed under the real situation of vortex MZM [306]; The sample quality should be optimized to improve the surviving temperature of MZM, and achieve homogeneous electronic environments in the bulk and on the surface; A controllable technique for vortex manipulation should be explored, which is capable for fast braiding within the quasiparticle lifetime. Those advances could remove the barriers for construction of topological qubits and merge the two large fields (Fig. 17), that are the condensed matter physics and the quantum computation, in a small piece of crystals of iron-based superconductors. This preprint is selected in CNKI Journal Translation Project https://jtp.cnki.net/bilingual/Navi/Detail? pykm=WLXB&year=2020&issue=11, of which the original version was published in Chinese (Acta Physica Sinica 69, 110301 (2020)).