Stealth Majorana Zero Mode in a Trilayer Heterostructure MnTe/Bi2Te3/Fe(Te,Se)

TL;DR

This paper explains the origin of zero-energy states in a MnTe/Bi2Te3/Fe(Te,Se) heterostructure, identifying magnetic defects as sources of Majorana zero modes and proposing simpler structures for their control.

Contribution

It reveals how ferromagnetic Mn-Bi antisite defects induce Majorana zero modes and offers a simpler heterostructure design for their generation and modulation.

Findings

Defects induce Yu-Shiba-Rusinov states and Majorana zero modes.

Majorana modes respond differently to magnetic field and temperature.

Proposed simpler heterostructure for better control of Majorana modes.

Abstract

Recent experiment reported the robust zero-energy states with strange properties in a trilayer heterostructure MnTe/Bi2Te3/Fe(Te,Se). Here, we give comprehensive understandings about the magnetic and electronic properties of the heterostructure, and propose ferromagnetic Mn-Bi antisite defects are generated in the topmost sublayer of Bi2Te3 and hidden below the MnTe layer. We further reveal the defect can induce two types of quasiparticles. One is Yu-Shiba-Rusinov state from defect itself, and another is Majorana zero mode from the superconducting phase domain wall induced by the defect. The two types of quasiparticles have very different response to magnetic field, temperature etc. The coexistence and mutual cooperation of both can explain experimental observations. Furthermore, we propose more simple heterostructure with superiority to generate and finely modulate Majorana zero modes.

Tables

Table 1. Table 1: Relative energy of different magnetic structure. The unit is meV/Mn.

x-AFM1	x-AFM2	x-FM
0	0.23	17.3
z-AFM1	z-AFM2	z-FM
0.43	0.43	18.8

Equations

$$H_{tot}=H_{sur}+H_d+H_{\Delta }.$$

$$H_{sur}$$

$$H_d$$

$$H_r=H_{sur}+{\Delta }_s(\mathbf{r}){\tau }_x.$$

$$H_r(\mathbf{k}\to -i\nabla )\psi (r,\theta )=E\psi (r,\theta ),$$

$$H_{sur}=A{k}^2+({\Lambda }_0-B{k}^2){\sigma }_z{\kappa }_z+v_F(k_x{\sigma }_y-k_y{\sigma }_x).$$

$$H_{eff}=v_F(k_x{\sigma }_y-k_y{\sigma }_x)+{\Lambda }_0{\sigma }_z{\kappa }_z,$$

$$H_{eff}=\left(\begin{array}[c]{cc}h_{+}&0\\ 0&h_{-}\end{array}\right)$$

$$g{\mu }_B{B}_c={\Lambda }_0$$

$$\Psi (r,\theta )=-\Psi (r,\theta +2\pi ),$$

$$H_{BdG}(k\to -i{\nabla }_r)\Psi (r,\theta )=E\Psi (r,\theta ).$$

$$H_{BdG}=H_{BdG}^{+}\oplus H_{BdG}^{-}$$

$$H_{BdG}^{\pm }=(h_{\pm }-\mu ){\tau }_z+{\Lambda }_0{\sigma }_z+{\Delta }_0{\tau }_x$$

$$u_{\uparrow }(r)$$

$$u_{\downarrow }(r)$$

$$H_Z=g{\mu }_B{\sigma }_z{B}_z=m_z{\sigma }_z.$$

$$H_t={\epsilon }_0\sum _{i=1,2}{\gamma }_i^{†}{\gamma }_i+t(L_z)\sum _{⟨ij⟩}{\gamma }_i^{†}{\gamma }_j$$

$$H_{tb}=\sum _i(\mu +\sigma m_0)c_{i\sigma }^{†}{c}_{i\sigma }+\sum _{⟨ij⟩}(t_1-\sigma t_2)c_{i\sigma }^{†}{c}_{j\sigma }+i{\lambda }_R\sum _{⟨ij⟩}{c}_{i\sigma }^{†}({\mathbf{S}}^{\sigma {\sigma }^{\prime }}\times {\vec{d}}_{ij}{)}_z{c}_{j{\sigma }^{\prime }}$$

$$H_{tb}(k_x,k_y)=2t_1{h}_0(k)+\mu +\left[m_0-2t_2{h}_0(k)\right]{\sigma }_z-2{\lambda }_R{h}_{12}$$

$$H=H_{tb}+H_d+H_{\Delta }$$

$$H_d=-J(c_{i_0\uparrow }^{†}{c}_{i_0\uparrow }-c_{i_0\downarrow }^{†}{c}_{i_0\downarrow }),H_{\Delta }=\sum _i(\Delta c_{i\uparrow }^{†}{c}_{i\downarrow }^{†}+h.c.)$$

$$c_{i\sigma }=\sum _n^{{}^{\prime }}(u_{i\sigma }^n{\gamma }_n-\sigma v_{i\sigma }^{n\ast }{\gamma }_n^{†})$$

$$\sum _j{\hat{H}}_{ij}{\varphi }_j=E_n{\varphi }_i$$

$${\Delta }_i(T)=\frac{g}{2}\sum _n^{{}^{\prime }}(u_{i\uparrow }^n{v}_{i\downarrow }^{n\ast }+u_{i\downarrow }^n{v}_{i\uparrow }^{n\ast })\tanh (\frac{E_n}{k_{B}T})$$

$$H_{BdG}^{+}\psi (r)=E\psi (r).$$

$$H_{BdG}^{+}=\left(\begin{array}[c]{cccc}(A-B)\nabla_{\nu,r}^{2}+\Lambda_{0}-\mu&-v_{F}(\partial_{r}+\frac{\nu+1}{r})&\Delta(r)&0\\ v_{F}(\partial_{r}-\frac{\nu}{r})&(A+B)\nabla_{\nu+1,r}^{2}-\Lambda_{0}-\mu&0&\Delta(r)\\ \Delta(r)&0&-(A-B)\nabla_{\nu+1,r}^{2}+\Lambda_{0}+\mu&v_{F}(\partial_{r}+\frac{\nu+1}{r})\\ 0&\Delta(r)&-v_{F}(\partial_{r}-\frac{\nu}{r})&-(A+B)\nabla_{\nu,r}^{2}-\Lambda_{0}+\mu\end{array}\right)\rule{7.22743pt}{7.22743pt}\rule{7.22743pt}{7.22743pt}$$

$$\psi(r)=\left(\begin{array}[c]{c}u_{1\uparrow}(r)\\ u_{1\downarrow}(r)\\ v_{1\downarrow}(r)\\ -v_{1\uparrow}(r)\end{array}\right).$$

$$u_{1\uparrow }(r)$$

$$u_{1\downarrow }(r)$$

$$l$$

$$\left[(A-B)({\partial }_r^2+\frac{1}{r}{\partial }_r-\frac{1}{4r^2})+{\Lambda }_0-\mu \right]u_{\uparrow }-\left[v_F({\partial }_r+\frac{1}{2r})+{\Delta }_0\right]u_{\downarrow }$$

$$\left[(A+B)({\partial }_r^2+\frac{1}{r}{\partial }_r-\frac{1}{4r^2})-{\Lambda }_0-\mu \right]u_{\downarrow }+\left[v_F({\partial }_r+\frac{1}{2r})+{\Delta }_0\right]u_{\uparrow }$$

$${\partial }_{r}u(r)$$

Full text

Since the discovery of high-temperature superconductivity in FeSe/SrTiO3FeSe-1 , the heterostructures involving Fe(Te,Se) have attracted intense interests to explore the new interfacial physical phenomena, such as spin density wave, interfacial phonon mode, topological physics, gating effect, competition between superconductivity and magnetism etcFeSe-2 ; FeSe-3 ; FeSe-4 ; FeSe-5 ; FeSe-6 ; FeSe-7 ; FeSe-8 ; FeSe-9 . Recently, an unexpected robust zero-energy state was experimentally observed in a trilayer heterostructure of MnTe/Bi2Te3/Fe(Te,Se)FeSe-10 , and the state shows some confusing properties, such as invisuality from surface topography but visuality from zero-bias dI/dV map, the robustness against external magnetic field, strange temperature evolution relation etc. Considering both Bi2Te3 and Fe(Te,Se) are two kinds of well-known topological materials, it is natural to connect the robustness of zero-energy state with the nontrivial topology of the system. However, the complexity of the trilayer heterostructures blurs and obscures such connection. Thus, the following urgent issues arise: is the robust zero-energy state a Majorana zero mode (MZM)? how is it generated? how to understand its strange properties? Clarification of these issues could establish such trilayer heterostructure as a new platform to study physics of MZMs, or supply an inspired strategy to construct more pragmatic system to realize MZMs.

In this work, we address the aforementioned issues through comprehensive understanding the magnetic and electronic properties of the trilayer heterostructures with density functional theory (DFT) calculations and extracting the underlying physics of zero-energy state with a simplified theoretical model. Our DFT calculations suggest single unit cell (UC) MnTe could probably generate weak antiferromagnetic (AFM) fluctuation rather than long-range AFM order. The consequent band structures are consistent with the experimental results from quasiparticle interference (QPI) measurement. Then, we propose there exists at least two different kinds of point defects in the trilayer heterostructure. One is from the vacancy of Te atom on the top surface of 1-UC MnTe/Bi2Te3/Fe(Te,Se). Another is the Mn-Bi antisite defect in form of Mn atom substituting Bi atom in top layer of Bi2Te3. The Mn-Bi antisite defect will generate a magnetic moment to couple the spin of the surface electrons of Bi2Te3 with ferromagnetic type. Such Mn-Bi antisite defect associated with ferromagnetic coupling not only induces the trivial Yu-Shiba-Rusinov (YSR) state but also drives a local quantum phase transition (QPT) of the Fe(Te,Se)-proximity-effect-enabled surface superconductivity of Bi2Te3. Namely, the superconducting order parameter of Bi2Te3 around the defect will acquire a sign-change beyond the QPT. Then, a degenerate pair of MZMs protected by an emergent mirror symmetry is induced and trapped by the phase domain wall near the Mn-Bi antisite defect. YSR state and MZM have very different response to the external magnetic field, temperature etc. Coexistence and mutual cooperation of both can explain all the experimentally observed confusing properties of the quasi-paricles lying in the superconducting gap. For the trilayer heterostructure with 2-UC MnTe, we find that the aforementioned picture is violated, and only trivial results can be obtained. Furthermore, we predict unexplored properties of the zero-energy state and propose more simple platforms to generate and modulate MZMs.

As verified by the experiment, the 1-UC MnTe shares the same crystal constant with 1-UC Bi2Te3. In DFT calculations, we treat 1-UC Bi2Te3 as a substrate and neglect the influence of bare Fe(Te,Se) substrate, since Fe(Te,Se) weakly modulates the band structures of Bi2Te3 except supplying proximity superconductivityBiTe-FeTeSe-1 ; BiTe-FeTeSe-2 . Let the MnTe layer relax freely until the force and energy criterion is satisfied. The final structure is shown in Fig. 1 (a). Te-Mn-Te sandwich structure is as same as it in $\alpha$-MnTe crystalalpha-MnTe-1 ; alpha-MnTe-2 . Take magnetic order of bulk $\alpha$-MnTe as a reference, various magnetic orders for 1-UC MnTe are considered. Our DFT results are consistent with the calculations in Ref. FeSe-10 . Furthermore, we calculate the band structures of 1-UC MnTe/Bi2Te3 by assuming several different kinds of long-range AFM or FM orders (See SMsSM for details). Fig. 1 (c) shows that the bands crossing Fermi level around $\Gamma$ point split into Rashba-type, which is inconsistent with the parabolic-type extracted by QPI measurementFeSe-10 . Thus, we propose there only exist weak AFM fluctuation in 1-UC MnTe, even though the DFT predict 1-UC MnTe with AFM-1 order has the lowest energy. We argue such divergence is not rare in low-dimensional systems. For example, DFT predicts checkerboard AFM order in monolayer FeSe/SrTiO3FeSe-mag-theo-1 ; FeSe-mag-theo-2 ; FeSe-mag-theo-3 ; FeSe-mag-theo-4 , but no experiment observes such order FeSe-mag-exp-1 ; FeSe-mag-exp-2 ; FeSe-mag-exp-3 ; FeSe-mag-exp-4 ; FeSe-mag-exp-5 ; FeSe-mag-exp-6 . The observation of gapless Dirac-cone-like surface state in MnBi2Te4 also hints the missing long-rang magnetic order in topmost MnTe layerMnBiTe-1 ; MnBiTe-2 ; MnBiTe-3 ; MnBiTe-4 . Our calculations also find tiny magnetic anisotropy in 1-UC MnTe, which further indicates the missing long-rang AFM order in 1-UC MnTe according to Mermin-Wagner theoremM-W-theo . Furthermore, the experiment observed the enhancement of superconductivity in 1-UC MnTe trilayer heterostructure. We propose such enhancement comes from the gain effect of AFM fluctuation in 1-UC MnTe accoring to spin-fluctuation-mediated superconductivity theory. Here, we do not extend the discussion about enhancement of superconductivity and will address it elsewhere. Once the long-range AFM order is missing, the calculated band structure of 1-UC MnTe/Bi2Te3 shown in Fig. 1 (b) has a simple parabolic-type band crossing Fermi level, which is consistent with result of QPI measurementFeSe-10 . Thus, we conclude that the first effect of 1-UC MnTe in trilayer structure is to introduce AFM fluctuation to enhance superconductivity and the second is to change the Fermi level by charge doping. Fig. 1 (d) shows only trivial band structure survives in the case with 2-UC MnTe. We do not further discuss it.

Remarkably, two types of Mn-Bi and Bi-Te antisite defects were experimentally verified in MnBi2Te4Mn-impurity . It is natural to expect similar antisite defects could exist on the interface between 1-UC MnTe and 1-UC Bi2Te3 in such heterostructure. In particular, the heterostructure is bulit with the molecular beam epitaxy (MBE) method. Fig. 2 (a) shows structures of two types of typical defects, i.e., the defect on the topmost sublayer 1-UC MnTe formed by Te vacancy and the Mn-Bi antisite defect with a Mn atom substituting a Bi atom on the topmost sublayer 1-UC Bi2Te3. Fig. 2 (b) shows the calculated differential charge density distribution around the Mn-Bi antisite defect. It is obvious the charge transfer mainly happens between the bottommost sublayer of 1-UC MnTe and the topmost layer of 1-UC Bi2Te3. Then, charge density distribution is nearly not affected on the topmost surface of heterostructure. It indicates the STM scanning topographic image is not easy to find and fix positions of the Mn-Bi antisite defects. Fig. 2 (c)-(d) shows the DFT simulated STM image, where only Te-vacancy defect can be observed. If the observed zero-energy states have close relationship with the Mn-Bi antisite defects, the randomness of them can be easily understood now. Our calculations also find that the Mn-Bi antisite defect spontaneously generates a local magnetic moment of 4.8 $\mu_{B}$ due to the magnetism of Mn atom. In the following, we will show that the ferromagnetic Mn-Bi antisite defects play a key role to form the observed zero-energy states, which are indeed the MZMs. Since Mn-Bi antisite defects in the inner layer of trilayer heterostructure cannot be detected by the outer-layer charge distribution but by the defect-bound superconducting quasi-particle tunnel spectrum, we call the zero-energy states the stealth MZMs. Thus, we conclude the third effect of 1-UC MnTe in trilayer structure is to introduce ferromagnetic Mn-Bi antisite defects in the topmost sublayer 1-UC Bi2Te3 in the trilayer heterostructure.

With the above prepared information, we can construct an effective model to answer the question how to generate the robust zero-energy state. The effective Hamiltonian can be expressed as follows,

[TABLE]

Here,

[TABLE]

Here, $H_{sur}$ is the low-energy effective Hamiltonian to describe the surface bands of 1-UC MnTe/Bi2Te3 thin filmTI-model , as indicated by the red-dashed rectangle in Fig. 1 (b). Note that the effect of structural inversion asymmetry is neglected, because the Rashba-like splitting induced by such effect is inconsistent with the QPI measurementFeSe-10 . $H_{d}$ describes the coupling between the ferromagnetic Mn-Bi antisite defect and surface electron. $H_{\Delta}$ is the proximity-effect-induced superconducting pair term from Fe(Te,Se). $v_{F}$ is Fermi velocity. $\mu$ is chemical potential. $\Lambda_{0}$ measures gap of Dirac cone opened by finite-thickness effect. $A$ and $B$ are coefficients of the second-order terms. $J(\mathbf{r})=J\delta(\mathbf{r})$ describing the ferromagnetic coupling with a constant $J$. $\mathbf{S}_{def}$ is the moment of defect with only $z$-component. $\Delta(\mathbf{r})$ is the superconducting order parameter from proximity effect. Note that $\Delta(\mathbf{r})$ is a constant without the defect. $\sigma$, $\kappa$ and $\tau$ are Pauli matrix to span spin, psudo-valley and particle-hole subspaces, respectively. $H_{tot}$ can be numerically solved in the lattice caseprevious-work (See SMsSM for details). The influence of ferromagnetic defect can be understood as follows. As $J(\mathbf{r}=0)$ increases from zero, $\Delta(\mathbf{r}=0)$ at the impurity site is depressed and goes to zero at a critical $J_{c}$, beyond which, $\Delta(\mathbf{r}=0)$ has a sign-change RMP-1 ; QPT-1 ; QPT-2 ; Order-1 ; Order-2 ; Order-3 . It means there exists a QPT at critical $J_{c}$. Meanwhile, the hole-type and electron-type branches of YSR statesYSR-1 ; YSR-2 ; YSR-3 undergo approaching-crossing-separating processes as $J(\mathbf{r}=0)$ increases from zero. The crucial physical quantity $J_{c}$ can be estimated by the calculations. We find $J_{c}\sim 4.2v_{F}k_{F}$ with $v_{F}k_{F}\sim 4.5$mV measuring the energy scale of surface Dirac cone of Bi2Te3. Then, $J_{c}\sim 19$mV, which can be further reduced by the atomic spin-orbit coupling.

Besides the impurity-induced trivial YSR states, we now turn to other quasi-particle related to the zero-energy state. Note that beyond QPT, i.e., $J>J_{c}$, there exists a phase domain wall for $\Delta(\mathbf{r})$ as shown in Fig. 3 (a). Namely, $\Delta(\mathbf{r})<0$ for $r<a_{0}$ and $\Delta(\mathbf{r})>0$ for $r>a_{0}$ with $a_{0}$ a characteristic length close to lattice constant of Bi2Te3. Then, in the case of $J>J_{c}$, $H_{tot}$ in Eq. (1) can be further reduced to the following form

[TABLE]

Here, $\Delta_{s}(\mathbf{r})=-\Delta_{1}$ for $r<a_{0}$ and $\Delta_{s}(\mathbf{r})=\Delta_{2}$ for $r>a_{0}$ with $\Delta_{1/2}>0$ and $0<\Delta_{1}<\Delta_{2}$. The phase of $\Delta_{1/2}$ is uniform and is omitted. In continuum limit, the eigen-equation of $H_{r}$ is

[TABLE]

which can be solved under boundary conditions with a 0-$\pi$ disk junction geometry as shown in Fig. 3 (a). The 0-$\pi$ disk junction can be obtained by bending a 0-$\pi$ line junction. Such deformation enables $\psi(r,\theta)$ to gain a phase factor $e^{i\kappa_{0}\tau_{0}\sigma_{z}\theta/2}$previous-work1 . Thus, the wave function $\psi(r,\theta)$ must obey anti-periodic boundary condition $\psi(r,\theta+2\pi)=-\psi(r,\theta)$previous-work ; previous-work1 ; previous-work2 . We have numerically proven such boundary conditionprevious-work . The details to solve Eq. (5) are shown in SMsSM . Without considering the second order of $k$ terms in $H_{sur}$ of Eq. (2), a pair of MZMs can be obtained. The wave function of the first MZM takes the form $\psi_{1}(r,\theta)=[e^{-i\theta/2}u_{\uparrow}(r),e^{i\theta/2}u_{\downarrow}(r),e^{-i\theta/2}v_{\downarrow}(r),-e^{i\theta/2}v_{\uparrow}(r),0,0,0,0]^{T}$ with $u_{\sigma}(r)=-v_{\sigma}(r)$. $u_{\sigma}(r)=a_{\sigma}J_{\mp 1/2}(k_{F}r)e^{r/\xi_{1}}$ for $r<a_{0}$ and $u_{\sigma}(r)=b_{\sigma}J_{\mp 1/2}(k_{F}r)e^{-r/\xi_{2}}$ for $r>a_{0}$. $J_{\mp 1/2}(k_{F}r)$ is Bessel functions with $\mp 1/2$ for spin up and down, respectively. The Fermi wave vector is $k_{F}=\sqrt{\mu^{2}-\Lambda_{0}^{2}}/v_{F}$. Decay length $\xi_{1/2}=v_{F}/\Delta_{1/2}$. The wave function of the second MZM is $\psi_{2}(r,\theta)=\mathcal{T}\psi_{1}(r,\theta)$ with $\mathcal{T}=i\kappa_{x}\sigma_{y}\tau_{0}\mathcal{K}$ the emergent time-reversal symmetry operator. If the second order of $k$ terms in $H_{sur}$ of Eq. (2) is taken into account, the unchanged in-plane spin texture of the surface state still gives effective topological superconductivityhelical-SC , and a pair of MZMs is robust. In such a case, the analytical solutions are not available, but the numerical calculations give the clear solutions of MZMs, as shown in Fig. 3 (b) (See SMsSM for details).

Now, we have two types of quasi-particles in superconducting state of the 1-UC MnTe/Bi2Te3/Fe(Te,Se), as shown in Fig. 3 (b). One is trivial YSR state from the ferromagnetic Mn-Bi antisite defect itself, another is MZM bounded by the superconducting phase domain wall induced by the defect. These two types of quasi-particles have very different response to the external magnetic field, temperature etc. However, the coexistence and mutual cooperation of two types quasi-particles can explain all the observed phenomena. The first one is to apply an external magnetic field. The energy levels of YSR must shift through the Zeeman coupling as shown in Fig. 3 (d). This behavior was observed by the experiment for the non-zero-energy bound states. However, a pair of MZMs is not split under the out-of-plane magnetic field, and the relevant forms of wave function of MZM solutions of Eq. (5) are unchanged except $k_{F}$ being renormalized to$\sqrt{\mu^{2}-(\Lambda_{0}\pm\mu_{B}B_{z})^{2}}$. $\mu_{B}$ and $B_{z}$ are Bohr magneton and $z$-direction of magnetic field strength. Our calculations show that $\mu\sim 100$meV and $\Lambda_{0}\sim 45$meV$\sim 375$T. The experimental applied magnetic field $<10$T has negligible effect to the MZMs. Thus, the pair of MZMs are robust against out-of-plane magnetic field as shown in Fig. 3 (c). Such robustness is actually protected by an hidden mirror symmetry, i.e., $\mathcal{M}=i\kappa_{x}\sigma_{y}\tau_{y}$ with mirror plane along the phase domain wall. obviously, in-plane $y$-direction magnetic field can break the mirror symmetry and split the pair of MZMs. We suggest further experiment to verify this prediction. Another issue is the temperature effect. The zero-energy state can only be observed below a temperature which is quite lower than $T_{c}$. This behavior can also be understood by our self-consistent calculations, from which, one can find that the QPT is quenched at quite lower temperature than $T_{c}$ (See SMsSM for details). The third interesting property is the behavior of the tunneling spectra under different tunneling transmission by tuning tunneling barrier. The experiment of trilayer heterostructure observed a saturated value of 0.22 in unit of $2e^{2}/h$ which is much smaller than vortex case with nearly a saturated value of 1 in unit of $2e^{2}/h$. The smaller value 0.22 strongly indicates the observed zero-energy state in the trilayer heterostructure experiment comes from the topmost sublayer of 1-UC Bi2Te3 as we propose, and the MnTe layer inevitably lower the tunneling transmission. Namely, the MZMs are stealth as we have pointed out.

Based on uncovering and understanding the underlying physics for the zero-energy states, we improve and propose new heterostructure to realize the MZMs, as shown in Fig. 3 (e). In the left plot of Fig. 3 (e), the MnTe layer is not necessary, and substrate Fe(Te,Se) can also be replaced by other superconductor such as NbSe2. The only matter is to introduce a FM impurity or defect on the topmost layer of Bi2Te3 thin film. The superiority of heterostructure in Fig. 3 (e) is the coupling between FM impurity or defect and top surface of Bi2Te3 thin film can be finely tuned. It would provide a method to manipulate the MZMs. The heterostructure can further be simplified by adopting a single superconducting topolocial material to replace the Bi2Te3/ superconductor heterostructure, as shown in the right plot of Fig. 3 (e).

Different scenarios can be argued to generate zero-energy state. The first one is the spontaneous vortex state. In this case, The MZMs can be taken as end states of one-dimensional vortex line. When the line is very short like the thin film case, the two end states can couple and splitsplitting . In our case, the pair of MZMs lies in the phase domain wall around the defect in the topmost sublayer of 1-UC Bi2Te3. Thus, there is no such coupling, and they are robust. The second one is the possible magnetic skyrmion. In such a case, the theory has proven only the skyrmion with topological charge $Q=2$ can bind a MZMsky-MZM-1 ; sky-MZM-2 ; sky-MZM-3 . Usually, $Q=2$ skyrmion is not stable and split into two $Q=1$ skyrmions. Besides, the DFT calculations give the AFM-type not FM-type exchange coupling, and the latter is necessary to form the familiar skyrmions. Furthermore, the size, shape and density of the skyrmions are fragile to the change of external magnetic fieldsky-split-1 ; sky-split-2 ; sky-split-3 ; sky-split-4 . Thus, the skyrmion scenario is not preferable.

In conclusion, we reveal the underlying physics of the robustness of zero-energy state in a trilyer heterostructure MnTe/Bi2Te3/Fe(Te,Se), and find that there are two types of quasi-particles. The first one is the YSR state from the FM Mn-Bi antisite defect itself and the second one is MZM from the superconducting phase domain induced by the defect. We show that consideration of both types of quasi-partcles can explain all the experimental observations. Furthermore, we propose some more simple platforms with superiority to generate and finely tune MZMs. Our studies provide new strategies to explore Majorana physics.

I First-principles calculations

I.1 Calculating methods

First-principles calculations were performed by density functional theory (DFT) using the Vienna ab initio simulation package (VASP) vasp1 ; vasp2 . The plane-wave basis with an energy cutoff of 360 eV was adopted. The electron-ion interactions were modeled by the projector augmented wave potential (PAW) PAW and the exchange-correlation functional was approximated by the Perdew-Burke-Ernzerhof-type (PBE) generalized gradient approximation (GGA) PBE . The strong correlated interaction of d-orbital electrons in Mn is taken into account by GGA+U method, where the Hubbard U value is adopted as 5 eV. And the spin-orbital coupling (SOC) is taken into account during all the self-consistent calculations.

I.2 Electronic and magnetic structure of the trilayer heterostructure

As detected by the experiment, the 1UC MnTe shares the same crystal constant with 1QL Bi2Te3. Here we treat 1QL Bi2Te3 as a substrate and let the MnTe layer relax freely until the force and energy criterion is satisfied. The final structure is shown in Fig. 5(b), the Te-Mn-Te sandwich structure is as same as it in $\alpha$-MnTe crystal. In bulk $\alpha$-MnTe, the magnetic order is in-plane intra-layer FM and inter-layer AFM, which is usually called A-AFM. However when it is reduced to a single layer, the ground magnetic state tends to be in-plane AFM. As shown in Table. 1,

the relative energy of different magnetic structures are calculated, form where it is easy to find that the collinear anti-ferromagnetic (x-AFM1, as shown in 4(b)) structure has the lowest energy, however, the energy difference is very small between two different AFM state or directions. Compared with some typical magnetic 2D materials whose magnetic anisotropy energy (MAE) usually reaches $\mu$eV, the MAE of 1UC-MnTe is quite small. Additionally, the band structures of 1UC-MnTe/Bi2Te3 with different magnetic orders have been calculated and the results are displayed in Fig. 4. It is obvious that with a long-range AFM magnetic order, the band will obtain a giant Rashba splitting, which could be easily detected by experiments, however that is negated by experiments.

So here we treat the MnTe layer as spin fluctuation instead of any long-range magnetic order. Under this setting, the band structure is obtained and shown in Fig. 5(c) where the red dots label the weight of the surface atomic layer’s contribution. The surface state of MnTe/Bi2Te3 is a typical massive Dirac cone.

For 2UC-MnTe structure, due to the inter-layer interaction it tends to form a A-AFM magnetic structure. Adopting this order, the band structure of 2UC-MnTe/Bi2Te3 is obtained and shown in Fig. 6 and Fig. 7, both of them is relatively trivial and that should the reason why Majorana cannot be induced. Also the band structure of 2UC-MnTe without long-range magnetic order has been considered and the results are shown in the main text.

We have investigated several typical defects that could exist in this heterostructure, which includes Mn atom vacancy, Mn-Bi antisite and Bi-Mn substitution. The results show that the Mn-Bi antisite can easily introduce a strong magnetic moment near the surface. The exchange interaction takes effect between the local magnetic moment and electrons of surface states. That is quite similar with the interstitial impurity in iron-based superconductor Fe(Te,Se), the Mn/Fe atom is located at the center of Te atoms. The magnetic moment of Mn impurity is obtained as 4.8 $\mu_{B}$.

II Theoretical Model

II.1 The effective Hamiltonian

At the ultra-thin limitation of a topological insulator, for example 1QL Bi2Te3 massive-Dirac ; massive-Dirac-2 , its surface state can be derived as

[TABLE]

Here, $A$ and $B$ are two constants, $\sigma$ and $\kappa$ are Pauli matrices which spans the spin and psudo-valley space, respectively. $\Lambda_{0}$ is the gap induced by the hybridization between surface states. To better understand the picture of this model, we can neglect the k-quadratic term first, and then Eq. (6) could be reduced to an effective Hamiltonian as:

[TABLE]

or in a matrix form as:

[TABLE]

where $h_{\pm}=v_{F}(k_{x}\sigma_{y}-k_{y}\sigma_{x})\pm\Lambda_{0}\sigma_{z}$. It is obvious that both $h_{+}$ and $h_{-}$ describes a massive Dirac cone which is gapped by an external magnetic field, however, the difference is the direction of this magnetic field is opposite. So this hybrid interaction can be equivalently regarded as the existence of a pseudo-magnetic field. As shown in Fig. 8, we can assume that two Dirac cones are placed in two opposite pseudo-magnetic field (pMF) and thus gain an energy gap whose amplitude is $\pm\Lambda_{0}$. And to preserve the time reversal (TR) symmetry, the two massive Dirac cones must locate in two opposite momentum point, i.e. $k$ and $-k$ respectively. For example, strain induced pMF in graphene has opposite directions when applied to the Dirac cone located at $K$ and $-K$ point pMF-1 . While the massive Dirac cones in ultra-thin topological insulator are degenerate at $\Gamma$ point.

So under the TR operation, the two massive Dirac cones will map to each other. Therefore, the TR operator should be $\hat{T}=i\sigma_{y}\kappa_{x}K$ on the basis $(c_{1\uparrow},c_{1\downarrow},c_{2\uparrow},c_{2\downarrow})$ where the subscript $1(2)$ labels two psudo-valleys and $\uparrow(\downarrow)$ labels the spin.

At the existence of an external magnetic field, the two massive Dirac cones correspondingly have different response, depending on the external magnetic field is parallel or antiparallel to its pMF. For the parallel condition, obviously, the gap is enlarged by the external field; while for the antiparallel condition, the gap will minish to close and then reopen, this process is actually a topological phase transition. We can estimate the critical magnetic field of the topological phase field by this equation

[TABLE]

where $g$ is the Lande factor of electrons, $\mu_{B}$ is the Bohr magneton. In 1QL Bi2Te3, $\Lambda_{0}\simeq 180meV$ which corresponds to a critical magnetic field $B_{c}\sim 1500\ T$; In 1UC MnTe/Bi2Te3, $\Lambda_{0}\simeq 45meV$ which corresponds to a critical magnetic field $B_{c}\sim 375\ T$. In above two systems, the critical magnetic field is extremely high where the system have already entered the Landau-level range. So when the external magnetic field is not so large, the surface state of 1QL Bi2Te3 or 1UC-MnTe/Bi2Te3 is very stable. Of course, the degenerate of two massive Dirac cones will be lifted by the external magnetic field, but the splitting energy is small. Additionally, TR symmetry is broken by the external magnetic field which results the TR partners $h_{+}$ and $h_{-}$ decouple with each other.

II.2 Emergence of Majorana Zero-energy Mode

The sub-surface defect can induce a strong local magnetic moment near the surface, which will be seen in the next section. The exchange interaction will trigger a quantum phase transition which makes a $\pi$-phase area around the defect when the surface enters superconductivity. At the boundary of 0 and $\pi$-phase area, there exists the bound state whose wave function obeys the anti-periodic condition anti-peri

[TABLE]

and satisfies the BdG equation

[TABLE]

Here, the BdG Hamiltonian can be expressed as

[TABLE]

with

[TABLE]

the intra-surface pairing is considered and the basis is $\Psi^{{\dagger}}=(c_{1k\uparrow}^{{\dagger}},c_{1k\downarrow}^{{\dagger}},c_{1-k\downarrow},-c_{1-k\uparrow},c_{2k\uparrow}^{{\dagger}},c_{2k\downarrow}^{{\dagger}},c_{2-k\downarrow},-c_{2-k\uparrow})$. That means we can separately solve the BdG equation for these two independent parts. The solution of $h_{+}$ or $h_{-}$ is equivalent to the condition where an external magnetic field exists, which we have solved in our previous work. Here the results are displayed

[TABLE]

where $k_{F}=\frac{\sqrt{\mu^{2}-\Lambda_{0}^{2}}}{v_{F}}$, $\xi_{0}=\frac{\Delta_{0}}{v_{F}}$, $c_{1}$ and $c_{2}$ are two normalized constants.

At the presence of an out-of-plane external magnetic field, a Zeeman term is added to the Hamiltonian

[TABLE]

Usually, $m_{z}\ll\Lambda_{0}$, thus the solution is lightly influenced.

In real materials, sometimes the quadratic term cannot be ignored. Actually it only has impact on out-of-plane spin texture, and the in-plane helical spin texture is not changed so that the effective p-wave superconductivity still can be induced, which is verified by experiments. And our numerical calculations will prove that the existence of Majorana zero-energy mode is not disturbed by these quadratic term at all.

In our previous work, we have proven that the 0-$\pi$ disk junction is equivalent to the vortex condition. So the inference is naturally raised that Majorana mode also exists at the core of a vortex in this heterostructure. However, experimental observation shows that there is no zero-energy state in the vortex. That can be explained by the tunneling between two Majorana fermions which locates at the top and bottom surface, respectively. In vortex condition, the two Majorana fermions can be seen as end states of a 1D vortex line, when the line is long enough two end states will not overlap with each other, however when the line is short, equivalently the topological insulator is too thin, one Majorana fermions can tunnel to another. A simple effective model is applied to describe this process, that is

[TABLE]

where $\gamma_{i,j}$ is the annihilation operator of the Majorana fermion, $\varepsilon_{0}$ is the intrinsic energy of the Majorana fermion which is set as 0, and $t(L_{z})$ is the hopping integral between two Majorna fermions which is a function of the thickness of topological insulators. As derived from Ref. splitting , $t_{L_{z}}$ exponentially decays with the thickness $L_{z}$ increasing. So when the thickness is large enough, the energy spectrum of the vortex core is a robust zero-energy peak; On the contrary, if the thickness of TI is too thin, the zero-energy peak will split into two peaks locating at $\pm t$, respectively. That coincides well with the experimental observation where two symmetric peaks appear in STS spectrum when the vortex is detected.

III Numerical calculations

III.1 Quantum phase transition induced by magnetic impurity

Now we consider a tight-binding model in a hexagonal lattice, whose Hamiltonian is

[TABLE]

where $\sigma$ is spin index, $\langle\rangle$ denotes the nearest neighbor site. Taken a Fourier transformation, (18) becomes as

[TABLE]

where $h_{11}=\cos k_{x}+2\cos\frac{k_{x}}{2}\cos\frac{\sqrt{3}k_{y}}{2}$, $h_{12}=\sin\frac{k_{x}}{2}\cos\frac{\sqrt{3}k_{y}}{2}\sigma_{y}-\sin\frac{\sqrt{3}k_{y}}{2}\cos\frac{k_{x}}{2}\sigma_{y}$. Expanding (19) near $\Gamma$ point, we obtain the effective Hamiltonian $h_{+}$. considering the system enters superconductivity and a magnetic impurity located at site $i_{0}$, the total Hamiltonian can be expressed as :

[TABLE]

where

[TABLE]

describes the exchange and proximity-induced superconductivity in the system. We can perform the Bogoliubov transformation

[TABLE]

where ′ denotes summation over the positive eigenvalues, and numerically solve the equations

[TABLE]

in the Nambu spinor representation $\phi_{i}=(u_{i\uparrow}^{n},u_{i\downarrow}^{n},v_{i\uparrow}^{n},v_{i\downarrow}^{n})^{T}$. The order parameter $\Delta_{i}$ should be self-consistently determined as

[TABLE]

To get the numerical results, we have adopted a $25\times 25$ lattice. The parameters are set as $t_{1}=1$ meV, $t_{2}=0.1\ t_{1}$, $\mu=0.8\ t_{1}$, $\Lambda_{0}=0.5\ t_{1}$, $\Delta_{0}=2\ t_{1}$, $\lambda_{R}=3.5\ t_{1}$. The results are shown in Fig. 9, where the quantum phase transition happens when $J$ crosses over $J_{c}$, meanwhile two branches of YSR states meet at zero energy and the order parameter at the impurity mutates, as shown in Fig. 9(a)(b). The consequence of QPT is the emergence of phase domain wall around the impurity as shown in Fig. 9(c). As the temperature increases, the QPT could be quenched far below the superconducting critical temperature Tc, as shown in Fig. 9(d).

III.2 Majorana zero mode

In Section. II we have proven that by adjusting the order of basis, the BdG Hamiltonian can be block diagonal. Thus we only need two solve one part of them. The radial form of the first-part BdG equation is

[TABLE]

Here $H_{BdG}^{+}$ and $\psi(r)$ has the following form

[TABLE]

[TABLE]

where $\nabla_{\nu,r}^{2}=\partial_{r}^{2}+\frac{1}{r}\partial_{r}-\frac{\nu^{2}}{r^{2}}$. And applying the Majorana condition

[TABLE]

we obtain the 1D radial differential equation list as

[TABLE]

It is hard to get the analytic solution. Here we will discuss some special conditions where analytic solutions can be given : (a) $B=0$, it’s equivalent to a Rashba system gapped by an external magnetic field, and the requirement of a Majorana mode yields $\mu^{2}+\Delta_{0}^{2}>\Lambda_{0}^{2}$; (b) $A=B=0$, the solution is listed as (14) and (15);

Next we will perform the difference method diff-method to get the energy spectrum and check the existence of Majorana zero-energy mode. Note that

[TABLE]

we can divide the radius into discrete points and apply this transformation

[TABLE]

Under the new basis

[TABLE]

where $\phi(r_{i})=\left[w_{1\uparrow}(r_{i}),w_{1\downarrow}(r_{i}),z_{2\downarrow}(r_{i}),-z_{2\uparrow}(r_{i})\right]$ and $r_{i}=r_{0}+ih$, Eq. 25 can be written as

[TABLE]

where

[TABLE]

and

[TABLE]

with $\nabla^{\prime}_{\nu,r}=\frac{\frac{1}{4}-\nu^{2}}{r^{2}}+\frac{1}{2h^{2}}$. The energy spectrum is obtained as Fig. 10, where the Majorana zero-energy mode appears in the zero-angular-momentum channel, i.e. $l=0$ ($\nu$=-1/2).

Note that at $l\neq 0$ ($\nu\neq-1/2$), there exist some near-zero-energy modes as shown in Fig. 10(b). We see they are just some trivial edge modes originating from the finite-size-disk geometric configuration. Firstly, the energy of these edge modes is much larger than the Majorana mode whose energy is about 10*-13* meV in our numerical calculations; Secondly, the wave functions of these edge modes are confined near the edge of the disk, in reality, the radius should be infinite so these edge modes cann’t emerge in experiments. However, the wave function of the Majorana mode concentrate on the core of the disk, i.e. around the impurity and satisfy the Majorana condition (LABEL:c1)-(30), as shown in Fig. 11(a)-(d).