Robust and clean Majorana zero mode in the vortex core of high-temperature superconductor (Li0.84Fe0.16)OHFeSe

The Majorana fermion, which is its own anti-particle and obeys non-abelian statistics, plays a critical role in topological quantum computing. It can be realized as a bound state at zero energy, called a Majorana zero mode (MZM), in the vortex core of a topological superconductor, or at the ends of a nanowire when both superconductivity and strong spin orbital coupling are present. A MZM can be detected as a zero-bias conductance peak (ZBCP) in tunneling spectroscopy. However, in practice, clean and robust MZMs have not been realized in the vortices of a superconductor, due to contamination from impurity states or other closely-packed Caroli-de Gennes-Matricon (CdGM) states, which hampers further manipulations of Majorana fermions. Here using scanning tunneling spectroscopy, we show that a ZBCP well separated from the other discrete CdGM states exists ubiquitously in the cores of free vortices in the defect free regions of (Li0.84Fe0.16)OHFeSe, which has a superconducting transition temperature of 42 K. Moreover, a Dirac-cone-type surface state is observed by angle-resolved photoemission spectroscopy, and its topological nature is confirmed by band calculations. The observed ZBCP can be naturally attributed to a MZM arising from this chiral topological surface states of a bulk superconductor. (Li0.84Fe0.16)OHFeSe thus provides an ideal platform for studying MZMs and topological quantum computing.

The Majorana fermion, which is its own anti-particle and obeys non-abelian statistics, plays a critical role in topological quantum computing. It can be realized as a bound state at zero energy, called a Majorana zero mode (MZM), in the vortex core of a topological superconductor, or at the ends of a nanowire when both superconductivity and strong spin orbital coupling are present. A MZM can be detected as a zero-bias conductance peak (ZBCP) in tunneling spectroscopy. However, in practice, clean and robust MZMs have not been realized in the vortices of a superconductor, due to contamination from impurity states or other closely-packed Caroli-de Gennes-Matricon (CdGM) states, which hampers further manipulations of MZMs. Here using scanning tunneling spectroscopy, we show that a ZBCP well separated from the other discrete CdGM states exists ubiquitously in the cores of free vortices in the defect free regions of (Li0.84Fe0. 16)OHFeSe, which has a superconducting transition temperature of 42 K. Moreover, a Dirac-cone-type surface state is observed by angle-resolved photoemission spectroscopy, and its topological nature is confirmed by band calculations. The observed ZBCP can be naturally attributed to a MZM arising from this chiral topological surface states of a bulk superconductor. (Li0.84Fe0. 16)OHFeSe thus provides an ideal platform for studying MZMs and topological quantum computing.

I. INTRODUCTION
In recent years, many recipes have been proposed for realizing MZMs [1][2][3], as a critical step towards topological quantum computation [4,5]. For example, MZMs are predicted to exist at the ends of a semiconductor nanowire with strong spin-orbital coupling (SOC), when it is in proximity to a superconductor and under a sufficiently large Zeeman field [2][3][6][7][8]. A quantized ZBCP has been observed in the tunneling spectrum of a hybrid device between superconducting aluminum and an InSb nanowire [9], providing compelling evidence for a MZM. However, fundamental quantum computing operations such as braiding of MZMs in nanowires [10] are challenging and yet to be realized after a tremendous amount of ingenious technical endeavor [11,12].
Topological systems with additional dimensions provide a broader hunting ground in the search of MZM, e.g. in the vortex cores of a topological superconductor or superconducting heterostructures [1][2][3]. For conventional s-wave superconductors, confined quasi-particles in the vortex core give rise to CdGM bound states with E=μΔ 2 /E F , where μ is a half integer (±1/2, ±3/2…) [13]. However, for a chiral p-wave superconductor, bound states still exist but with integer μ (0, ±1, ±2…) -the E=0 state, or zero-energy mode, is a MZM [14]. Chiral p-wave superconductors are extremely scarce, but Fu and Kane have proposed that proximity effects from an s-wave superconductor on topological surface states would produce a two-dimensional system whose Hamiltonian effectively resembles a spinless p±ip superconductor, and thus it can host Majorana bound states in its vortices [15] (The spinless nature is due to the fact that topological surface state is spin non-degenerate). Based on this scenario, ZBCPs that may potentially correspond to MZMs have been found in the vortex cores of topological insulator / superconductor heterostructures (e.g. Bi2Te3/NbSe2 and Bi2Te3/FeTexSe1-x) [16][17][18], and topological surface states of bulk superconductors (e.g. CuxBi2Se3 and FeTexSe1-x) [19][20][21][22]. However, an unambiguous identification or isolation of a MZM is still lacking in these systems. For materials with a large Fermi energy (EF), such as Bi2Se3/NbSe2, CuxBi2Se3 and Bi2Te3/FeTexSe1-x, a mode at zero energy is buried under many densely packed CdGM states in the vortex core, their energies separated by only a few μeV. FeTexSe1-x has a small EF, so its ZBCP can be readily distinguished from the CdGM states. However, there are intrinsic impurity effects of interstitial Fe and heavy Te doping [23,24]. A zero-bias impurity state was observed on interstitial Fe which is insensitive to magnetic field [23], but whether this is related to a MZM is yet to be confirmed. A ZBCP can be observed in a fraction of vortices in certain FeTexSe1-x samples [21], while it cannot be observed in samples annealed to reduce impurities [24]. Therefore, to study the largely-unknown properties of MZMs, one needs to find a system that can host MZMs in a clean and robust manner, with weak pinning and scattering from impurities, which would lay the foundation for further operations such as braiding of MZMs.
To reach these goals, we need a superconducting system that has a large superconducting gap, a small EF, a stoichiometric lattice (at least in the sublattice responsible for superconductivity) with few defects, a short superconducting coherence length (ξ) to potentially reduce vortex pinning, and a topologically nontrivial electronic structure. The large gap usually also means high superconducting transition temperature (Tc) and large separations between bound states in vortices, which will make the operation temperature easier to reach in practice. As topological non-trivial band structures have been predicted to exist in many ironbased superconductors [25][26][27], (Li1-xFex)OHFeSe with a Tc as high as 42 K appears to be a promising candidate. It belongs to the heavily-electron-doped family of iron selenide that only have electron Fermi pockets, and its EF is merely 50~60 meV [28,29]. Unlike FeTexSe1-x, the superconducting FeSe layers in (Li1-xFex)OHFeSe are stoichiometric ( Fig. 1(a)). The remaining questions are whether it could exhibit topologically nontrivial behavior and if it hosts MZMs. To address these questions, we conduct density functional theory (DFT) and dynamical mean field theory (DMFT) calculations and ARPES measurements on (Li0.84Fe0.16)OHFeSe, and we find non-trivial band inversion and topological surface states. We studied the vortex states of (Li0.84Fe0.16)OHFeSe by low-temperature scanning tunneling microscopy (STM). In the core of impurity-free vortices, we indeed discovered discrete CdGM states which are clearly separated from a unique state located exactly at zero bias. Such a clean and robust zero energy mode exists in all the free vortex cores of (Li0.84Fe0.16)OHFeSe and matches all expectations of Fu and Kane's theory, providing an ideal platform to further explore the properties and applications of MZMs.

A. Electronic structure calculations and ARPES measurements
In Fig. 1(b), we present our DFT combined with DMFT calculations of the electronic structure of (Li0.75Fe0.25)OHFeSe (details of the calculation are described in Appendix and section S1 of Supplementary Materials). Along the Γ -Z direction, the three flat bands around the Fermi level are dominated by the Fe 3dxy, 3dxz and 3dyz orbitals whereas the dispersive band is composed of Se 4pz orbital and the 3dz 2 orbital of Fe in the (Li1-xFex)OH layer. This dispersive Se 4pz band has odd parity at the Γ and Z points. It crosses the Fermi level and the Fe 3d bands along Γ -Z, giving rise to a non-trivial band inversion and band topology. Here Fe atoms in the (Li1-xFex)OH layer play an important role as they hybridize strongly with the Se 4pz orbital and change the dispersion and position of the Se 4pz band. Without the Fe atom in the (Li1-xFex)OH layer, the band inversion along the Γ -Z direction would not have existed and the band topology would have been trivial (as demonstrated in Fig. S1(b)). The spin-orbital coupling (SOC) also plays an important role. Without SOC, the Fe 3dxz and 3dyz bands are degenerate at Γ point and are in the Γ5 + states, whereas the dispersive Se 4pz band is in the odd Γ2state at Γ point. With SOC, the doubly degenerate Γ5 + states are separated to a singlet lowlying Γ6 + state and an upper Γ7 + state, and the Γ2state changes to Γ6state [26]. Under the C4v symmetry inherent in the tetragonal crystal structure, the two Γ6 + and Γ6state-derived 6 bands with dominating Fe 3dxz and Se 4pz orbital characters hybridize with each other along Γ -Z direction, and open a ~2.5 meV gap around their crossing point (marked by the circle in Fig.  1(b)). The Z2 invariant is calculated to be 1 after this SOC gap opens [26], which suggests (Li0.75Fe0.25)OHFeSe is in a topologically non-trivial phase. Since one cannot make a transition between two states with different topologies without closing the band gap, namely, here between a topological nontrivial bulk with band inversion and the topological trivial "insulating" vacuum, the transition region, i.e. surface, would hosts gapless topological surface states (TSS) in the band gap. More specifically, the hybridization gap between the Γ6 + and Γ6states at the surface are removed by the discontinuity of the surface. Indeed, our calculations show that Dirac-cone-like surface states centered at appear on the (001) surface as shown in Fig. 1(c). The topological surface states have helical spin texture and can produce the MZM when it is in close proximity to s-wave superconductivity [15,26,27].
Guided by these calculations, we have conducted further ARPES measurements on high quality (Li0.84Fe0.16)OHFeSe films grown on LaAlO3 by hydrothermal epitaxy (Tc~42 K) [30,31] (see Appendix for details). Fig. 1(d) shows the ARPES spectrum of a cleaved film across the Brillouin zone (BZ) center as shown in Fig. 1(a). Multiple parabolic bands can be clearly resolved below the Fermi surface, and there is one flat band around -300 meV. Moreover, an electron band with its band bottom at ~ -50meV is observed at M ( Fig. 1(j, k)). These observed bands and the calculated ones in Fig. 1(b) have one-to-one correspondence, although qualitatively. The relative positions and bandwidths may differ due to correlation effects. We also note that certain bands, such as the other electron-like band at M, were not observed here due to known matrix element effect. A precise band calculation for heavily electron doped iron selenides is still so far challenging, however, qualitative features of the measured bands are reproduced by the calculation in Fig. 1(b), and the generic features here resemble several iron chalcogenide superconductors [32].
Intriguingly, in the region near the Fermi energy as plotted in Fig. 1(e), some finite spectral weight within the bulk band gap can be clearly observed. The second derivative of the region containing these spectral weight ( Fig. 1(f)) display an "X"-shaped band structure. Fig. 1(g) shows the momentum distribution curves (MDCs) of the data in Fig. 1(e) near EF. Through one-Lorentzian-peak fit (shown in Fig. S2(a)), we found the width of the spectral peak near k//=0 has a minimum at E=-20 meV ( Fig. 1(h)). This would indicates a position of band crossing point. Then we apply two-peak fit to MDCs with assuming each peak has a constant FWHM as that of the single peak at E=-20 meV (the scattering rate variations are expected to be negligible in such a small energy window). The results display a nearly-linear, Dirac-cone like band dispersion, as shown in Fig. 1(g), and its overall fitting quality is better than that of the one-peak fit (see Fig. S2(c) for comparison). In Fig. 1(i) we plot the E-k relations extracted from the two-peak fit, a linear fit is found to be plausible near the crossing point which gives a vF = 5.5(±1.5)×10 4 m/s and ED= 20(±2) meV (see Fig. S2(d) for more details). The kF is determined to be ~0.028 Å -1 through the MDC near EF.
We note the existence of topological surface states is determined by the band topology. The qualitative agreement between the measured and calculated band structures throughout BZ would suggest they have the same topology. Therefore, the predicted topological surface states provide the most natural explanation for such a nearly-linear, Dirac-cone-like band dispersion at Γ. Nevertheless, one should examine whether it is related to any other known electron pockets in the (Li1-xFex)OH and FeSe surface layers exposed by cleavage. On the FeSe surface layer, an electron-like pocket would only appear around Γ after heavily doping electrons through potassium dosing, and thus it is well above the Fermi energy in the case here [33]. There is possibly an electron pocket around Γ in the (Li1-xFex)OH surface layer, based on our previous quasi-particle interference (QPI) measurements of (Li0.8Fe0.2)OHFeSe bulk crystal [34]. However, its band bottom is at about 50 meV below EF with a much larger Fermi velocity, so it cannot be the Dirac cone observed here. Therefore, we conclude that the observed Dirac band dispersion can only be attributed to the topological surface state predicted in our calculations.
The FeSe surface of (Li1-xFex)OHFeSe is fully gapped [34,35], as shown in Fig. 2 below. If this topological surface state is on the FeSe surface, it thus should be gapped. A superconducting gap of about 10 meV can be observed on the electron-like band around , as shown in Fig. 1(j, k), and further illustrated in Fig. 1(l) (details described in caption). The absence of gap in our data for the topological surface state around may be attributed to the fact that this gap is small, while the resolution of ARPES is about 6 meV. Moreover, ARPES collects data over a millimeter-sized region and thus contains large contribution from the nonsuperconducting and disordered (Li1-xFex)OH surface, which has low energy spectral weight near .
B. STM characterization of the sample surface and superconducting gap STM measurement was conducted in a cryogenic STM with a base temperature of 0.4 K and had an energy resolution of 0.36 meV (see Appendix and Supplementary Materials section S3). Figure 2(a) shows typical topography of a cleaved (Li0.84Fe0.16)OHFeSe film. We observed both FeSe-terminated (left) and (Li1-xFex)OH-terminated (right) surfaces. The Fig. 2(a) inset shows the atomically-resolved FeSe lattice with a0 = 3.8 Å. Dimer-like defects at the Fe site were frequently observed (dashed circle), which are likely substitutional impurities. The typical tunneling spectrum (dI/dV) of FeSe surface is shown in Fig. 2(b) (blue curve). It displays a full superconducting gap with two pairs of coherence peaks (referred as Δ 1 and Δ 2 ) and a flat bottom, as observed in bulk crystals [34,35]. Since our band calculations and ARPES measurements revealed surface states at , it is possible that the double gaps actually originate from the bulk and the surface band, as discussed below. For the (Li1-xFex)OH surface, the tunneling spectrum typically shows metallic behavior with a weak dip off the Fermi level (green curve in Fig. 2 To determine the gap of the FeSe surface more precisely, we fit the tunneling spectrum with a Dynes function (see section S4 of Supplementary Materials for details). The observed gaps are broader than an isotropic gap with a similar size (grey curve in Fig. 2(c)). However, the line shape of the smaller gap (Δ 1 ) can be well fitted by an empirical anisotropic gap function: Δ 1 (k)=Δ 1 min +(Δ 1 max -Δ 1 min )|cos2θ| (red curve in Fig. 2(c)) after subtracting a sloping background.
The fit yields Δ 1 max = 8.2 meV, Δ 1 min = 5.7 meV and a small Dynes term Γ of 0.1 meV. Note that such a Δ 1 max value can also be estimated as half distance between the corresponding two coherence peaks, and so does Δ 2 max . However, the larger gap Δ 2 has additional broadening which makes the determination of gap anisotropy difficult. The green curve in Fig. 2(c) gives a simulated spectrum with Δ 2 max = 14.5 meV (see Supplementary Materials for details). We note that the size of Δ 1 (even Δ 1 max ) is obviously smaller than the ARPES measured gap (~10 meV) on the pocket, thus Δ 1 is unlikely from the pocket, and it can only be attributed to the gap of the topological surface state obtained through proximity from the bulk. On the other hand, Δ 2 max is larger than 10 meV and thus Δ 2 should correspond to the bulk gap on the pocket. As will be shown below, the gap size measured by ARPES reflects the average size of Δ 2 . QPI measurements were carried out on FeSe surface. Figure 2(d) presents a typical fast-Fourier-transformed (FFT) QPI image taken at E = 5 meV and T = 4.2 K -more data are presented in section S5 of Supplementary Materials. The dominant features are the ring-like patterns centered at (0, 0), (π, π) and (0, 2π), which arise from the scattering between electron pockets at [34,35]. Figure 2(e) summarizes the averaged FFT line cuts through q = (π, π), which displays an electron-like dispersion. Using q = 2k, a parabolic fit yields a band bottom at -57(±7) meV and kF = 0.21(±0.02) Å -1 , consistent with the reported band structure at of (Li1-xFex)OHFeSe [28,29]. In relation to the topological surface states at , we note that a spinhelical structure with spin-momentum locking would strongly suppress backscattering, causing them invisible to QPI.
Under a magnetic field of 10 T perpendicular to the surface, a zero-bias-conductance (ZBC) map is taken on a 36×36 nm 2 area of FeSe surface ( Fig. 2(f)), as shown in Fig. 2(g). Vortex cores are visible as bright circular regions. We found that many of the vortices are pinned by dimer-like defects, as indicated by the corresponding arrows in Figs. 2(f) and 2(g). The pinned vortices all have a "dark spot"-like feature near their center, which are induced by impurity states. However, we can still find un-pinned or "free" vortices which emerge in defect-free areas, such as the one marked by a dashed circle in Fig. 2(g).

C. Discrete bound states in the free (un-pinned) vortex cores
The tunneling data of four different free vortices are presented in Fig. 3. Figure 3(a) plots the dI/dV line-cut taken across Vortex 1 (inset image, along the arrow); the red curve is measured at the vortex center. It is remarkable that multiple (five) discrete peaks are observed near the core center. There is one peak located exactly at zero bias (a ZBCP), with the other peaks distributed symmetrically around it. The energy spacing between these 5 peaks is close to 1.5 meV. Figure 3(c) presents the spatial evolution of the spectra in a color plot. On leaving the core center, the intensities of these discrete peaks (marked by -2, -1, 0, 1, 2 here, also referred as E-2, E-1, E0, E1, E2 hereafter) decrease and vanish at ~2nm away. Intriguingly, when the discrete peaks fade out, a pair of much broader peaks (shaded regions in Fig. 3(a)) show up at higher energies. Those peaks keep shifting to higher energy on moving away from the core, giving an "X"-shaped pattern in Fig. 3

(c). Figures 3(b) and 3(d) present similar data taken on
another free vortex (Vortex 2). Multiple low-energy core states including a ZBCP were also observed, albeit with a slightly smaller energy spacing (~1 meV), and there are also "shifting" high-energy states away from the core center. The behavior of these states, except the one at zero bias, are expected for the CdGM vortex states in quantum limit [13,36].
Such a ZBCP accompanied by discrete CdGM states was repeatedly observed in free vortices. In regions with few defects, such as the one shown in Fig. 3(e), more than one free vortices (Vortex 3 and 4) could be found ( Fig. 3(f)). At zero field, a clean superconducting gap is observed over this region, which could vary by about 2 meV along the long cut ( Fig. 3(g)). Figures 3(h) and 3(i) show the dI/dV line-cuts of Vortex 3 and 4 respectively, focusing on the region close to the vortex center (±1.4 nm) and low energies (±5 meV). There are five discrete peaks with a ZBCP in most spectra, with the energy spacing of 0.8~1.0 meV, well above our energy resolution (0.36 meV).
The positions of E-2, E-1, E1, and E2 peaks vary in different vortices, and they even vary at different locations in the same vortex. As shown in Figs. 3(h) and 3(i), upon leaving the center, the positions of E-2, E-1, E1, and E2 peaks shift away from the dashed lines that represent the peak positions at the center. This is because the CdGM states are spaced by δE=Δ 2 /E F , and the superconducting gap varies in space ( Fig. 3(g)). In contrast, the energy position of the E0 peak is always at zero bias, independent of local gap, which suggests that it is most likely protected by some global properties, such as topology. In addition, since conventional CdGM states are not located at zero energy, but at half-integer multiples of Δ 2 /E F , the origin of the E0 peak is highly nontrivial.
The high-level core states are closely packed and form the broad "shifting" peaks. In theory the spacing between high-level states will decrease and make them undistinguishable, while the maximum intensity of core states will shift to high energy on leaving the core center [36][37][38], resulting in two "splitting" peaks. This behavior has been widely observed, such as in NbSe2 (Ref. [37]); however, individual core states are rarely seen due to the small δE=Δ 2 /E F for most conventional superconductors. Here we clearly observe both discrete low-level core states and the quasi-continuous high-level states, benefitting from the relatively large δE of (Li1-xFex)OHFeSe (as shown below) and sufficiently high resolution at 0.4 K (For comparison, see Supplementary Materials section S6 for the T=4.2 K data).
The well-separated low-level core states enable quantitative analysis on them. We first applied multiple-Gaussian-peak fitting to the summed low-energy spectra taken near the vortex center (over a ~ ±0.7 nm range to reduce the uncertainty from peak fluctuations). The results are shown in Figs. 4(a)-4(d) for Vortex 1~4, respectively. Here each core state peak corresponds to one Gaussian peak. The fitted peak energies are directly labeled in Figs Table SI in Supplementary Materials. From the fitting we can reveal several important facts: Firstly, the fitted energy of E0 for all the vortices are negligibly small, typically one order of magnitude smaller than the energy resolution (0.36 meV). The fluctuation range of E0 (error bars in Fig. 4(e)) always covers zero. Meanwhile, the FWHM of the E0 peak is in the range of 0.59~0.80 meV, larger than the energy resolution. However, one cannot attribute such an additional broadening (0.23~0.44 meV) to two unresolved conventional CdGM states sitting close to zero bias, assuming they have the energies of ± 1 2 Δ 2 /E F , because it would require a superconducting gap with a mean size of 3.6 ~ 5.0 meV (assuming EF = 57 meV), and the observed local gap size near these vortices are well above that (see Fig. S6 and Table SII, the size of Δ 1 =(Δ 1 max +Δ 1 min )/2 are all around 7 meV). Besides, the FWHM of E0 is obviously narrower than the FWHM of other CdGM states. As discussed below, the gap anisotropy may induce additional broadening but cannot well account for the finite width of E0. Therefore, within our experimental resolution, E0 is a zero-energy state with finite width, whose origin is yet to be explored. Secondly, the energies of E±2 and E±1 states are nearly symmetric with respect to the Fermi level (with different values for different vortices). However, the spacing between E2 and E1 is always slightly larger than that between E1 and E0. The ratio of (E 2 -E -2 )/(E 1 -E -1 ) falls between 2.1~2.5 for Vortex 1~4 (see Table SII). This is unexpected if we assume all five states are from one single electronic component, since the spacing between neighboring states will decrease on going to higher energy [36,38]. Therefore, these may actually arise from two bands with different gap size and/or band bottom position, as we did observe a double gap dI/dV spectrum and a topological surface state. Meanwhile, the local gap inhomogeneity and anisotropy should also affect the energy of core state for different vortices. Taking into account these effects, we measured the local gap size (Δ 1 max , Δ 1 min and Δ 2 max ) of the region where these vortices emerge (see Fig. S6 and Table SII for details). As Δ 1 is attributed to the topological surface state, and its approximate mean size is Δ 1 =(Δ 1 max +Δ 1 min )/2; considering the Dirac point is at 20 meV below the Fermi level (defined as E D ), then it could reasonably account for the E0 and E±2 states for most free vortices through the general formula E=μΔ 2 /E D and μ = 0, ±1. In Fig. 4(f) we plot for different vortices (red spots), a linear fit gives |E 2 | =(0.97Δ 1 ) 2 /E D . We note there is a more specific calculation on the vortex states of proximity induced superconductivity in topological surface state [39], which gives E=±0.83Δ 2 /(Δ 2 +E D 2 ) 1/2 for the two second lowest states. This formula can estimate E±2 states by taken Δ~1.06Δ 1 , which is also in between Δ 1 min and Δ 1 max (see Table SII).
For the topologically trivial bulk state around M with EF = 57 meV, its gap is Δ 2 . Its lowestorder core states should be at ± 1 2 Δ 2 2 /E F . The mean value of Δ 2 is difficult to determine but Δ 2 max can be estimated by the coherence peaks. We found a linear fit of Fig. 4(f)). The value of 0.72 Δ 2 max ≈ 10 meV is consistent with the ARPES gap size on the band around the M point (Figs. 1(j) and 1(k)), which should be close to Δ 2 . Although the two data points marked by dashed circles deviate from the fits, they still show a monotonic relation with gap. We speculate this could be due to some local variation of carrier concentration, which is difficult to determine precisely for each vortex. Despite this uncertainty, the most self-consistent understanding is that the E0 and E±2 are from the topological surface state while E±1 are from the trivial bulk band. Thirdly, the fitted FWHMs of the E±2 and E±1 states vary between 0.8 to 2 meV. We expect that this is at least partially due to the sizable anisotropy of Δ 1 and Δ 2 . The broadening of E±2 caused by gap anisotropy of Δ 1 (defined by α=(Δ 1 max -Δ 1 min )/(Δ 1 max +Δ 1 min )) can be estimated by yielding an additional broadening in the range of 0.7~1.5 meV for Vortex 1-4. Meanwhile, although the gap anisotropy of Δ 2 is unknown, if we take the same value with α and assume Δ 2 =0.72Δ 2 max , its contribution to the FWHM of E±1 is in the range of 0.3~0.6 meV. Then the resulting total broadenings are compatible with the fitted FWHMs. However, the estimated total broadening of E±2 (1.1~1.9 meV) is significantly larger than the FWHM of E0, indicative of the special properties of the latter. It is also worth noting that the core states have different but comparable weight, as reflected by the peak area summarized in Table SI. Moreover, the intensities of E±2 and E±1 are asymmetric. This is consistent with the particle-hole asymmetry expected for a superconductor with small EF [36,40,41]. Figure 4(g) shows the spatial dependence of the intensity of the core states (dI/dV values at the peak energies, extracted from Fig. 3(c) for Vortex 1). The E0 state is more concentrated at the vortex center (especially compared to E±2), consistent with the calculated behavior of the vortex state of a p+ip superconductor [42,43]. An exponential fit (red solid curve) gives a decay length of 1.4 nm, which is an estimation of coherence length. Theoretically, the low-level core states should have an additional oscillation with a period of λ F =2π/k F (characterized by Bessel functions [42,43]) besides the exponential decay. Here kF is 0.21 Å -1 for the bulk band and ~0.03 Å -1 for the surface band (from QPI and ARPES measurements), which gives λ F B = 3.0 nm and λ F S = 21.0 nm. These are both significantly larger than ξ, thus the λ F oscillation cannot be observed here.

D. The core states of impurity-pinned vortices
A large number of vortices are pinned by dimer-like defects (Figs. 2(f)-2(g)), as the dimerlike defects at the Fe site can locally suppress superconductivity by inducing impurity states. Figures 5(a) and 5(b) show the tunneling spectra measured at T = 0.4 K on two typical dimerlike defects under B = 0 T and 10 T (or 11 T). Multiple sharp impurity state peaks can be observed at B = 0 T (blue curves), and the peak positions can vary for different defects (see Fig.  S8 for additional data). In some cases, a zero-bias peak can show up at the defect site in zero field (Fig. 5(b)). At high field, vortices pinned by these two defects were observed in the ZBC map (the insets of Figs. 5(a) and 5(b)). The spectra on the defect sites now display signatures of split impurity states, such as the double peaks marked by arrows in Fig. 5(a) (see also Fig.  S8). The zero-bias peak in Fig. 5(b) is also split away from zero bias. This suggests the dimerlike defects are mostly magnetic [44]. Figures 5(c) and 5(d) show the tunneling spectra along the arrows in the insets of Figs. 5(a) and 5(b), respectively. There is no zero-bias peak at the vortex center. Upon moving away from the defect site the impurity states decay, while the broad high-energy core state peaks appear (shaded regions in Figs. 5(c)-5(d)). These observations apparently suggest a competition between CdGM and impurity states. We note that the pinned vortices always have a "hole" in their center in the ZBC map, which also indicates that the low-energy core states are strongly suppressed by the magnetic impurity (see also Figs. S5(d-f)). Nevertheless, we note that in some case the zero-bias vortex state recovers a certain distance away from the defect site, as highlighted by the green curve in Fig. 5(d).

III. DISCUSSION AND CONCLUSION
Our data represent the cleanest and most robust zero-bias mode observed in vortices so far. The robustness is manifested in the following aspects: (i) It is always present at zero bias in free vortex cores, regardless of variations in the underlying superconducting gap.
(ii) It can survive to high magnetic field due to the short coherence length. The ξ of (Li1-xFex)OHFeSe is significantly smaller than that of FeTexSe1-x (~3 nm) [21] and Bi2Te3/NbSe2 (~30 nm) [16], which also means a higher probability to observe multiple free vortices and a lower probability of pinning by impurities during vortex motion. These are critical for braiding and fabricating qubits.
(iii) The large superconducting gap and high Tc makes the system robust against temperature fluctuations and relatively easy to study in experiments.
The cleanness is manifested in three aspects: (i) The ZBCP is unambiguously separated from other low-lying CdGM states and free from impurity effects in free vortices.
(ii) The FeSe layer is intrinsically stoichiometric, so that improved sample quality will allow more free vortices in larger defect-free regions, enabling further manipulation of the MZMs.
(iii) The width of the ZBCP is the narrowest of all the core states. The proximity effect of bulk superconductivity induces topological superconductivity in chiral topological surface states [1][2][3]15], which naturally generates a MZM in the vortex core and explains our observation. Nevertheless, in section S10 of the Supplementary Materials, we use a two-band model to simulate the vortex states with a variety of possible relevant pairing functions (fully gapped to be consistent with our experiment) of the bulk superconducting state, such as pure s-wave, d+id', and nodeless d-wave, in both the presence and absence of SOC. None of these scenario would lead to a robust ZBCP. Other origins may lead to a ZBCP in tunneling spectrum such as Kondo effect, SIS tunneling can also be excluded, since there is no impurity in free vortex and the tip used here is non-superconducting. Therefore, so far MZM is the most likely origin of the observed clean and robust ZBCP.
To summarize, our experimental and theoretical findings compellingly demonstrate that (Li0.84Fe0.16)OHFeSe, a heavily electron doped FeSe-based superconductor, is topologically nontrivial. We present the cleanest and most robust zero-energy modes in vortices so far, which enable us to obtain the important properties of MZMs with unprecedented accuracy and reliability, including their spatial distribution, scattering rate, response to magnetic impurities, and tunneling quantum efficiency as compared to conventional CdGM states, which would put strong constraints on theory. Our work thus presents an ideal and practical platform to further study the properties of MZMs, explore their manipulation such as braiding, and construct MZM-based qubits, which opens a new, clear route to rapid progress in both the fundamental understanding and potential applications.

APPENDIX: EXPERIMENTS AND METHODS
The high-quality single-crystalline superconducting films of (Li0.84Fe0.16)OHFeSe were grown on a LaAlO3 substrate by a matrix-assisted hydrothermal epitaxial method, as described in Refs. [30,31]. The full-widths at half-maximum (FWHM) of their X-ray rocking curves are 0.1~0.12 degrees, indicative of their high quality. The thicknesses of different films vary from 100 to 400 nm. STM measurements were conducted in a UNISOKU cryogenic STM at T = 0.4 K or 4.2 K. The sample was cleaved at 78 K in ultrahigh vacuum with a base pressure of 5 × 10 -11 Torr and immediately transferred into the STM module. Pt-Ir tips were used after being treated on a clean Au (111) substrate. dI/dV spectra were collected by a standard lock-in technique with a modulation frequency of 741 Hz and a typical modulation amplitude ΔV of 0.1 mV at 0.4 K and 1.0 mV at 4.2 K. ARPES data were taken under an ultra-high vacuum of 5×10 -11 mbar, with a Fermi Instruments discharge lamp (21.2 eV He-Iα light) and a Scienta DA30 electron analyzer. The energy resolution is 6 meV and the angular resolution is 0.3°.
We used fully charge self-consistent density functional theory combined with dynamical mean-field theory (DFT+DMFT) [45] to calculate the electronic structure of (Li, Fe)OHFeSe in the paramagnetic state. The DFT part is based on the linearized augmented plane wave method as implemented in WIEN2K [46]. A Hubbard U of 5.0 eV and Hund's coupling J = 0.8 eV were used in the calculations, consistent with previous calculations [47]. The DMFT quantum impurity problem was solved using the continuous time quantum Monte Carlo method [48] at temperature T = 116K. We used the experimentally determined crystal structure including the internal atomic positions in our calculations [49]. The surface states were calculated through the iterative Green's function method [50], taking into account the renormalization and shifting of the DFT bands due to strong electronic correlation effects. More specific details of the calculation are described in section S1 of supplementary materials. Red markers indicate the Dirac-cone like dispersion extracted from the two-Lorentzian-peak fitting (blue and dashed curves). (h) The FWHM obtained from the one-peak fit (shown in Fig. S2(a)), as a function of energy. Error bars represent 95% confidence bounds. (i) E-k dispersion extracted from the two-peak fit in panel (g), and the linear fit around the crossing point. Error bars represent 95% confidence bounds of the fitted peak position.

S2. Additional details on the extraction of the surface state dispersion
To explore the band dispersion of surface states, we applied both one and two Lorentzianpeak fittings to the MDCs crossing Γ, as shown in Figs. S2(a) and (b). The MDCs have been normalized by their intensities near k//=0 before the fitting. We found that the peak width (FWHM) in one-peak fitting displays a minimum at E=-20 meV, as shown in Fig. 1(h), which would indicate a position of band crossing. In two-peak fitting (Fig. S2(b)), we assume the scattering rate variations are negligible in a small energy window, so that each peak has a constant FWHM as that of the single peak width at E=-20 meV. The resulting peak positions then clearly displays a nearly linear, Dirac-cone like dispersion (marked by red arrows). To compare the "goodness" of one-and two-peak fits, in Fig. S2(c) we plot the sum of squares due to error (χ 2 ) of these fittings. One sees that at the energies away from the crossing point (~ -20 meV), two-peak fit gives lower χ 2 than that of the one-peak fit, which suggests a dispersive band is more appropriate to account for the spectral weight at Γ.
Then the extracted E-k dispersion with error bars is plotted in Fig. S2(d). At the energies close to crossing point (±5 meV), we found that a linear fit to dispersion gives lower χ 2 with respect to (two) quadratic fit (as labelled in the figure), suggesting that it is more likely from a Dirac-cone. However, at higher energies the dispersion observably deviates from linear fit and became more quadratic. Since the calculated SOC gap of (Li, Fe)OHFeSe is actually small (~2.5meV), such deviation may be due to the influence of bulk state, despite that in theory the dispersion right at the Dirac point should be always linear. Due to limited energy/momentum resolution here, a more precise determination of the whole surface state dispersion will require further study. Here the linear fit near the crossing point gives ED= 20(±2) meV and vF = 5.5(±1.5)×10 4 m/s. The kF is determined to be ~0.028 Å -1 directly through the MDC near EF. The goodness of one and two-peak fits for MDCs at various energy, reflected by the sum of squares due to error (χ 2 ). At the energies away from the band crossing point (~ -20 meV), the two-peak fit is better than one-peak fit. (d) The E-k dispersion extracted from the two-peak fit (error bars represent the 95% confidence bounds of the fitting parameter). Red and green curves are linear and two quadratic fits to the dispersion near the band crossing point, respectively. Linear fit has a smaller χ 2 than that of the quadratic fit.

S3. Calibration of STM energy resolution at T = 0.4 K and bias voltage offset
The energy resolution of an STM is limited by thermal and electrical noise broadening. The total broadening can be estimated as 3.5kBTeff, where Teff is the effective electron temperature.
To check Teff, we measured the superconducting gap of a Pb/Si(111) film at T = 0.4 K, as shown in Fig. S3(a). A standard BCS fit gives Δ = 1.39 meV, Teff = 1.18 K and a small Dynes term Γ = 0.005 meV which accounts for a finite quasi-particle lifetime. The total broadening (defined as the energy resolution in the text) of the STM is then given by 3.5kBTeff = 0.36 meV.
The STM bias voltage applied to the sample usually has a small offset. This offset can be calibrated by measuring I-V curves at different setpoints (tunneling resistance), because all the I-V curves should intersect at a single point where V = 0 and I = 0. Figure S3(b) shows such a calibration performed on the metallic (Li,Fe)OH surface which has a more linear I-V curve, giving an accurate determination of the bias offset of -0.093(±0.001) meV. All tunneling spectra presented in this paper have been corrected for this bias offset.

S4. Fitting of the superconducting gap of (Li1-xFex)OHFeSe
The double gapped tunneling spectrum of (Li1-xFex)OHFeSe cannot be fitted by an isotropic gap (or two), as shown by grey dashed curve in Fig. 2(c). Taking into account the gap anisotropy, the smaller gap Δ 1 can be well fitted by a gap function of: Δ 1 (k) = Δ 1 min + (Δ 1 max -Δ 1 min )|cos(2θ k )|, where Δ 1 max and Δ 1 min are the maximum and minimum size of Δ 1 in k-space.
The superconducting DOS is given by the Dynes formula [S10]: The total tunneling conductance is given by: is the derivative of the Fermi-Dirac function at Teff = 1.18 K. Before fitting, a linear background is subtracted from the original dI/dV to reduce the line-shape asymmetry (blue dots in Fig. 2(c)). The fitting yields Δ 1 max = 8.2 meV, Δ 1 min = 5.7 meV, and Γ = 0.1 meV. The larger gap Δ 2 is much broader than Δ 1 and even an anisotropic gap function cannot give a satisfying fit. The green dashed curve in Fig. 2(c) is a simulated curve using the combined gap function Δ(k)=Δ 1 min +(Δ 1 max -Δ 1 min )|cos(2θ k )|+A[Δ 2 min +(Δ 2 max -Δ 2 min )|cos(2θ k )|] , where A is the relative weight between the two gaps. The parameters are chosen as Δ 1 max = 8.2 meV, Δ 1 min = 5.7 meV, Δ 2 max = 14.5 meV, Δ 2 min = 8.5 meV, Γ = 0.1 meV and A = 0.8. It is seen that Δ 2 has additional broadening with comparing to simulated curve, which indicates that the anisotropy of Δ 2 cannot be precisely derived from fitting. However Δ 2 max can still be estimated as half the distance between two outer coherence peaks, as the simulated curve always peaks at the maximum value of the gap.

S5. Additional QPI Results
Additional dI/dV maps and FFT images are shown in Fig. S4. S4. dI/dV maps taken on a 36 × 36 nm 2 area of the FeSe surface, and the corresponding FFT images (four-fold symmetrized). The energy for each map is labeled in the image. Each map has 300 × 300 pixels.
S6. The vortex core states of free and pinned vortices measured at T = 4.2 K Figs. S5(a)-S5(c) show dI/dV spectra taken across the same free vortex in the Fig. 3(a) inset, but measured at T = 4.2 K. At this elevated temperature, only a single broad peak near zero bias can be observed at the core center, and it "splits" upon moving away from the core center. This is due to enhanced thermal broadening at 4.2 K, which smears out the discrete low-level core states. As the most likely consequence, the zero bias mode was not resolved in our early STS work on Li1-xFexOHFeSe (Ref. [34]).
Figs. S5(d)-S5(f) show dI/dV spectra taken across a pinned vortex, measured at T = 4.2 K. One can clearly see from the color plot that the low-energy core states are gapped out by the impurity at the core center.

S7. Local superconducting gap of the regions where free vortices emerge
We measured the local superconducting gap in the regions near Vortex 1~4, to exclude any impurity state in these regions and obtain the local gap size. Figure S6(a) summarizes the dI/dV spectra measured near the center position of Vortex 1~4 (at T = 4.2 K). All the spectra display a clean superconducting gap with double coherence peaks and no in-gap states. To extract the values of Δ 1 max and Δ 1 min , we fit the "inner" gap of each spectrum using the same method described in section S4, as shown in Figs. S6(b)-S6(e). The values of Δ 2 max are obtained by measuring the half distance between two outer coherence peaks in each spectrum. The value of Δ 1 max , Δ 1 min and Δ 2 max for Vortex 1~4 are summarized in Table SII. The values of Δ 1 max , Δ 1 min and Δ 2 max for Vortex 1~4 are listed in Table SII. S8. Quantitative characterization of the vortex core states As described in the main text, we fit the summed spectra near the centers of Vortex 1~4 with multiple Gaussian peaks. The fitting curves are shown in Figs. 4(a)-4(d) and the fit parameters are summarized Table SI, including the peak position, FWHM and peak area. Uncertainties are listed following the main value. We note the energies of the core states exhibit spatial fluctuations, which is likely due to variation of local gap size (as shown in Figs. 3(h)-3(i)). To characterize this uncertainty, we performed similar Gaussian-peak fitting to each individual spectrum taken near the core center of Vortex 1~4 (typically within ±0.7 nm of the core center). The yielded peak positions at different sites reflect their spatial fluctuation. The total range of peak fluctuation are shown in Table SI (in parentheses), and correspond to the error bars in Fig. 4(e) (the uncertainties in the fit peak positions are negligibly small compared to the fluctuation range). In Fig. S7 we plot all the fitted peak positions measured at different sites near the core center.
The bulk electronic structure of (Li1-xFex)OHFeSe has two electron pockets at M. Before they are folded into the reduced 2Fe/cell Brillouin zone, they can be viewed as a pocket around X and another around Y in the 1Fe/cell Brillouin zone. Upon folding, the two pockets intersect -see Fig. S9(a). Following [S13] we adopt a k·p model and compactify it on a lattice. This is sufficient to describe the low energy quasiparticles. The normal-state single-particle Hamiltonian is, in momentum space, h k =ϵ k +d k σ 3 +g k s σ 1 Henceforth, , , are Pauli matrices acting on the two effective orbitals, and , , are Pauli matrices acting on spins. Agterburg et al proposed, in the continuum limit [S13]: where m is the effective mass, and α and β are coefficients. The above form of h k takes proper account of the symmetry of the effective orbitals at the M point near the Fermi level. We rotate the coordinate frame (by 45 o about z) and spin axis independently, so that with modified coefficients. Notice that rotation of the spin axis by constant Euler angles does not alter the singlet pairing. The advantage of the rotated h k is an emergent superficial C4v symmetry, with the two orbitals behaving effectively as xz and yz. (Notice that the resulting Fermi pockets are elongated along x and y in the new frame.) We then compactify the model on a lattice, with ϵ k →-2t( cos k x + cos k y )-μ, d k →2t'( cos k x -cos k y ), g k →2λ( sin k y , sin k x ) where t' accounts for the Fermi pocket anisotropy, and λ measures the strength of the spinorbital coupling (SOC). We set t = 1, μ = -3.5t (or E f = 0.5t) to have shallow electron pockets mimicking the experimental situation. The small parameters t' and λ will be specified in due course. A further advantage of the rotated hk arises after the compactification: the symmetric hopping integrals, the d-wave like anisotropy and the SOC can all be defined on nearestneighbor bonds.
In the superconducting state, the vortex bound states would be the usual Caroli-de Gennes-Matricon states for s-wave pairing. Here we first consider a less trivial case, so-called nodeless d-wave pairing. In such a case, there is a full gap on both pockets, but with a sign change from X to Y (see Fig. S9(a)). One may worry that upon hybridization, the energy gap on the reconstructed bands, illustrated in Fig. S9(b), would have to be nodal. However, this does not have to be so if the hybridization is from SOC, as shown in Ref. [S13]. We now ask how the nodeless d-wave would impact on the vortex bound states. We write the pairing part of the Hamiltonian as, in momentum space (for the uniform case), where Δ 1 is the onsite part, with opposite phase on the two bands before reconstruction, see Fig. S9(a) for illustration, Δ 2 is the amplitude of d-wave pairing on nearest bonds, and is2 is the spin antisymmetric tensor accounting for singlet pairing (note the definitions of Δ 1 and Δ 2 here are different from those in the main text). Note that σ3 transforms as d-wave under rotation, hence both components in the gap function behave as d-wave. The resulting pairing gap is nodeless if |Δ 1 |>4|Δ 2 |. Since Δ 2 leads to gap variation on the Fermi surface, which is weak in experiments (about 10% from ARPES), we shall ignore Δ 2 for the moment (and we verified that including this part does not alter the results qualitatively). Writing Δ k in real space with the non-uniform pairing in a vortex state, and using Eq. S1 for Δ 1 , we calculate the local density of states (LDOS) along a line cut approaching the vortex core. We set t ' = 0.0125t and λ = 0.025t. Experimentally, the pocket anisotropy is tiny, and the energy scale of SOC is small, of the order of 2 meV versus Ef ~ 50 meV. Our parameters are in rough correspondence. We set ξ = 6, but a comparison to the experimental coherence length would be difficult since the model is not defined on the material lattice. As for the resulting pairing gap Δ 0 (the mean value of the gap on the Fermi surfaces of the two reconstructed bands), we use a relatively larger value to improve the resolution in LDOS, without affecting the qualitative behavior of the bound states, as argued before. In Fig. S9(c) we use Δ 0 = 0.2t. Interestingly, we see a clear ZBCP at the vortex center. Experimentally, the band structure parameters should not change significantly, but the gap size varies slightly from place to place (see Table SII). Nonetheless, the ZBCP is very robust in our STS data. To see whether this is the case in the nodeless d-wave scenario, we set Δ = 0.26t = 1.3Δ 0 -the resulting LDOS is shown in Fig. S9(d) (red curve). Here the conductance peak at the vortex core is shifted away from zero energy. Therefore, the zero mode in Fig. S9(c) is accidental, in contrast to the robustness in the experiment. In fact, we are able to obtain the bound state energy level analytically in the limit of t ' = 0, E∼lω 0 ±g, l = 0, ±1, ±2, ⋯ where g = |g k f | is the energy of SOC on the Fermi surface. The integer angular momentum is a nontrivial result of the spin-momentum texture caused by SOC. Clearly, zero modes may appear but only accidentally when g is an integer multiple of ω 0 , and such zero modes are not Majorana, since the quasiparticle and its anti-particle have different angular momentum and hence cannot be identical. We remark that the stability of the above energy level requires g > 0.
If we set g = 0 formally, we would have MZMs at l = 0 but they are coupled and gapped out by corrections beyond the quasi-classical approximation. This is best seen by ignoring SOC in the first place: The two bands are independent, each manifesting trivial s-wave pairing, even though there is a relative negative sign between the pairing in the two bands.
Other pairing symmetries: For completeness here we include the vortex states for other pairing symmetries in Fig. S10. Figure S10(a) is for pure s-wave in the two-band model as above with the same SOC. There is no ZBCP, similar to the pure s-wave case without SOC. Figure S10(b) is for chiral d-wave in a simple one-band model. We should mention that on a square lattice the two d-wave pairing functions are both one-dimensional irreducible representations of the point group, hence in general they do not mix. We include the chiral dwave here because it is fully gapped and topologically nontrivial. There are gapless edge states in this case. Nonetheless, there is no zero mode in the vortex core. This can actually be understood by solving the vortex problem in the quasi-classical approximation. The result is that the energy levels are similar to that in the pure s-wave case for all chiral pairings of even internal angular momentum L ( L = 0, ±2 for s-and chiral d-wave), although the bound state wave function depends on. Finally, we present the result for a nodal d-wave in Fig. S10(c) for comparison. It is known that in this case the LDOS develops a single broad bump, a kind of resonance, around zero energy and near the vortex core. None of the above cases is consistent with our experimental data.