Selective equal spin Andreev reflection at vortex core center in magnetic semiconductor-superconductor heterostructure

Sau, Lutchyn, Tewari and Das Sarma (SLTD) proposed a heterostructure consisting of a semiconducting thin film sandwiched between an s-wave superconductor and a magnetic insulator and showed possible Majorana zero mode. Here we study spin polarization of the vortex core states and spin selective Andreev reflection at the vortex center of the SLTD model. In the topological phase, the differential conductance at the vortex center contributed from the Andreev reflection, is spin selective and has a quantized value (dI/dV)Atopo=2e2/h\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${(dI/dV)}_{A}^{topo}=2{e}^{2}/h$$\end{document} at zero bias. In the topological trivial phase, (dI/dV)Atrivial\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${(dI/dV)}_{A}^{trivial}$$\end{document} at the lowest quasiparticle energy of the vortex core is spin selective due to the spin-orbit coupling (SOC). Unlike in the topological phase, (dI/dV)Atrivial\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${(dI/dV)}_{A}^{trivial}$$\end{document} is suppressed in the Giaever limit and vanishes exactly at zero bias due to the quantum destruction interference.

the MZM as well as other quasi-particle states. In the calculation of the differential tunneling conductance, the spin-polarized scanning tunneling microscope (STM) tip is modeled as a normal lead providing incident particles and receiving scattered particles, as depicted in Fig. 1. It is known that the total local differential tunneling conductance consists of the normal term proportional to the local density of states and an additional term arising from the Andreev reflection. However, in this work, we only consider the Andreev reflection part and focus on the similarity and difference for AR in topological phase and topological trivial phase.
The paper is organized as follows: we firstly describe the model Hamiltonian and present the spectra and the corresponding wave functions inside the vortex core in the topological phase. We then apply Fisher-Lee-Landauer-Bütikker formula to calculate the differential tunneling conductance. We then discuss the SOC induced SESAR in the topological trivial phase. Finally we will give a brief summary.

Models and Results
Model Hamiltonian for device. We study a 2D semiconductor with a Rashba SOC, which is hybridized to an s-wave SC and under a Zeeman field 8 [see Fig. 2]. The system is described by SLYD model. The Cooper pairs in the semiconductor are induced through the proximity effect. It resembles an effective chiral p x + ip y topological SC at the interface. The model Hamiltonian reads, ( with m * , μ, α R and V z being the effective mass of electron (in the 2D thin film), chemical potential, strength of the Rashba SOC, and the Zeeman field, respectively. The Pauli matrices ( , , ) T is the proximity-induced on-site pairing gap function in the 2D semiconductor. According to the AZ classification 27 , the above Hamiltonian  D belongs to D class since only particle-hole symmetry Kτ x is preserved, where K is the complex conjugate and τ x is an operator describing particle-hole transformation.
The dispersions of the Hamiltonian in Eq. (2), corresponding to the helical chirality λ = ±1, are, Figure 1. Illustration of the system. A spin polarized STM tip detects the vortex in the heterostructure. Tightbinding lattice model is used in calculation.

Figure 2.
Illustration of the semiconductor-superconductor heterostructure studied in this paper. The semiconductor thin film is described by a 2D electron gas with a Rashba SOC. Superconductivity in the semiconductor is induced by proximity effect. A magnetic field is pointed down which induces a Zeeman energy V z in Eq. (2).
where we set  1 = for convenience. For V z ≠ 0 and |μ| ≤ V z , there is an energy gap 2V z at Γ point (k = 0). If the pairing potential Δ → = Δ r ( ) 0 is uniform and if the criterion μ > Δ + V z 2 0 2 2 is satisfied [8][9][10]28 , the system will open a gap at the dispersion's outer wings without closing the Zeeman gap at Γ point. In this case, the system is essentially the same as an effective 2D spinless p x + ip y topological SC with chiral MZM 6 , so we expect a localized MZM in the vortex core of the system.
To study the quasiparticle excitations of the Hamiltonian in Eq. (1) with a single vortex, we use Bogoliubov-de Gennes (BdG) equation, T is used here. Because of the p ar t i c l e -h ol e s y m m e t r y i n t h i s B d G e qu at i on ( 5 ) are eigenfunctions with eigenenergies E n and −E n , respectively. The Bogoliubov quasiparticle operator is defined as, n s ns s n s s The necessary condition for MZM is † γ γ = for zero-energy mode. There are three main practical approaches to solve BdG equation (5). The first one is to use corresponding 2D tight-binding model [29][30][31] of Hamiltonian in Eq. (1), and solve the problem in a lattice. The second one is to adopt a disc geometry to solve the BdG equation (5) and use orthogonal Bessel functions 8,[32][33][34][35][36] . And the third one is to adopt spherical geometry to utilize harmonic spherical function or associated Legendre polynomials 15, [37][38][39] . In this work, we shall use tight-binding model on 2D square lattice [see Fig. 1], where the Hamiltonian becomes, i i i i SC where t D and μ D are the nearest-neighboring hopping integral and the chemical potential, respectively, α R is the Rashba SOC strength, and r ij denotes the unit vector between site i and j. V z, D is the Zeeman energy due to the z-direction magnetic field. We assume that the pairing order parameter inside the vortex has form i 0 where r is the distance of the lattice site from the vortex core, and ϕ is the azimuthal angle of → r ( ), and ξ describes the size of the vortex. We diagonalize the BdG equation (5) for Hamiltonian in Eq. (7) to obtain the eigenvalues and the corresponding eigenvectors by using the Feast Eigenvalue Solver for large sparse matrix. In this paper, we consider a lattice size of N = 199 × 199 with open boundary condition, and the vortex is located at the center of the lattice. The Hamiltonian in Eq. (5) is a 4N × 4N hermitian matrix.
Spin polarized MZM in the Vortex Core. The energy spectra of the vortex states are plotted in Fig. 3 in topological phase region with parameters given in the figure caption. In the topological phase, we expect a MZM in the vortex core and a MZM at the edge in the infinitely large system. In a finite size system, the vortex core state and the edge state have a hybridization, leading to a pair of the MZMs with energies ± E 0 very close to zero. Our numerical calculations agree with this analysis and the calculated ]. Both satisfy the MZM condition † γ γ = in Eq. (6) to a high accuracy. Note that we only consider a single vortex in our model calculation, since the edge MZM will be easily destroyed in real system, where there are many vortices. Below we will only focus on the bound states in the vortex core. As shown in Fig. 4(a), one can see that the MZM's wave-function at the center of the vortex is fully polarized with spin-up: |u ↑ |≠0 and |u ↓ | = 0. As a comparison, the wave function of the first excited state is shown in Appendix A, which is spin-down at the center of the vortex. Transport calculation. To experimentally detect the MZMs, it is important to examine some unique transport features of the MZMs, such as quantized conductance 21  spin-polarized transport properties of the MZM in the vortex core in the model, which can be tested in STM/STS measurement. We consider the STM tip as a 1D normal lead, illustrated in Fig. 1, and the Hamiltonian  L for the semi-infinity lead is, here t L and μ L are lead's nearest-neighboring hopping coefficient and chemical potential respectively. V L is the potential of the magnetic field on the lead. The lead's vertex (site 0) contacts the device at site p, then where the lattice labels L0 and Dp are the connected (touched) points from N/S junction. Here we use t L = 1.2, μ L = 0 for the lead, and choose |V L | = 0.4 to polarize the spin, adjusting the spin polarized direction to the local wave function's of the vortex MZM. And, we set the coupling coefficient t c = 0.6 in the Giaever limit, which can simulate the barrier strength at the interface between normal lead and the SC device. Now we consider a single electron with energy E injected from the spin polarized normal lead (STM tip, shown in Fig. 1), and reflected by the SLTD device. We use the basis c c c c [ ] ↓ in our calculation. In this situation, the 4 × 4 submatrix of scattering matrix (S-matrix) at the probed point site p,   Base on the S-matrix, the differential tunneling conductance contributed from the Andreev reflection can be calculated by using Landauer-Bütikker formula 40 SESAR in topological phase. For the SESAR detecting system showed in Fig. 5, we calculate the differential tunneling conductance dI/dV(E) as a function of incident electron's energy in Eq. (14). Usually as for the STM experiment 15 , the measured conductance comes from two parts: normal conductance (proportional to local density of states) and Andreev reflection [44][45][46][47] . In this work, we only focus on the Andreev reflection part. Andreev reflection due to MZM is very different from that in usual SC. In the usual SC case, Andreev reflection is known to be weak and can be neglected in the Giaever limit. But Andreev reflection due to MZM is very different and its strength at zero energy remains 2e 2 /h as pointed out previously and also shown in our numerical results below in Fig. 5. The differential conductance dI/dV as functions of energy due to the Andreev reflection are plotted in Fig. 6 for a topological phase of the model. dI/dV is spin selective and shows a quantized zero-bias peak, i.e., dI/dV = 2e 2 /h at the center of the vortex core. We can analyse the S-matrix in Eq. (24), and see that the outgoing hole is spin up for spin-up incident electron. This is the reason why spin polarized STM experiment can see the unique signal of MZM. As we expect, the width of the zero-bias peak at the lattice site away from the core center becomes narrow and the spin polarization dependence becomes weak. SOC induced SESAR in topological trivial phase. In this subsection, we study the topological trivial phase. We use same models Eqs (7), (11), (12) and methods in previous sections, and keep parameters unchanged except setting, for simplicity, 2 which belongs to topological trivial region. We perform numerically calculation for the spin-polarized differential conductance, and the results are shown in Figs 7 and 8.
Interestingly, we find there is no SESAR signal at E = 0, namely dI/dV(0) = 0 exactly, as Figs 7 and 8. It indicates perfect quantum destructive interference 48,49 . And it can be explained by the particle-hole symmetry which makes the anomalous Green's function 26 vanishing in the vortex core center, because of E n  δ | | . However, it will not happen if there exists MZM whose energy is almost zero. Besides, in topological trivial case, for ground state wave function, the electronic component |u ↑ | is no longer equal to the hole part |v ↑ |. Since electron-hole reflection r he by spin-up particles u ṽ| When the SOC strength is increasing, as in Fig. 7, the SESAR signal becomes stronger and stronger. And it leads to our main conclusion that SESAR can be induced by SOC. This signal is different from the MZM-based SESAR, since there is no MZM in topological trivial phase. In Fig. 8, we vary the hopping element t c from 0.3 (dirty limit) to 1.2 (transparent limit), the conductance becomes more and more broad, consistent with the BTK theory 23 . In a short conclusion, the condition for the appearance of SOC-induced SESAR is large SOC strength and small tunneling barrier between STM tip and device.
To analytically investigate the effect of the SOC on the localized vortex core states, we start with a discussion of the Hamiltonian in continuum space on a spherical surface 26 . Note that the finite size effect can be reduced by increasing the radius of sphere. The Hamiltonian reads, }. The SC pairing function is Δ(θ) = Δ 0 tanh(R sin(θ)/ξ) with R as radius of sphere, and θ for polar angle, φ for azimuth angle. We firstly turn off SOC, i.e., α R = 0, and then turn on SOC as perturbation term. In the absence of SOC (α R = 0), the 4-by-4 BdG Hamiltonian can be decoupled into two 2-by-2 blocks. In each block, we notice that L z + σ z /2 or L z −σ z /2 provide a good quantum number, thus |m, i〉 (i = ± represents blocks) can be used to label the in-gap vortex states: (L z ±σ z /2)|m, ±〉 = (m±1/2)|m, ±〉. Therefore, the wave function for the localized vortex core states are |m, . The components for the corresponding wavefunction are, with P l m the associated Legendre polynomial. Note that the m = 0 channel is in the particle-hole relationship with the m = −1 channel, the m = 1 channel is in the particle-hole relationship with the m = −2 channel. And we also assume that the quasi-particle wave functions are all orthogonal and normalized.
Then, we switch on the SOC (α R ≠ 0), and assume it is small compared with η. For the following discussion, we only focus on the vortex core center. In the presence of SOC, K z = L z + (σ z − τ z )/2 provides a good quantum number 26 . Then, the SOC will mix the in-gap states, so that the corrected wave function should be eigenvector of K z . And we find that the corrected wave function up to the first order is exact since all spherical harmonics with l ≥ 1 vanish in the vortex core center. In other words, only u ↑/↓,0 and v ↑/↓,0 survives. Thus the corrected wave function is given by~ Thus the corrected wave function for |0, ±〉 and |−1, ±〉 are, u v ⁎ or u v ↓ ↓ ⁎ ) contributes to the Andreev reflection. From this point of view, we may expect that there is also equal spin Andreev reflection in the vortex core center, even when the system is topological trivial. However, as for 0, contain only nonzero electron component u ↓ or hole component v ↓ . They are also spin polarized down, but there is no Andreev reflection signal. Therefore, this phenomena can also be treated as SESAR.
Moreover, we should emphasize that such SESAR is totally induced by SOC, if the SOC strength is larger, ⁎ũ v R α will be larger, then the SESAR effect become more and more obvious, which is consistent with numerical results of square lattice in Fig. 7.

Summary
In this work, we have studied selective equal spin Andreev reflection at vortex core center in magnetic semiconductor-superconductor heterostructure described by the model proposed by Sau, Lutchyn, Tewari and Das Sarma. We solve the BdG equation for a single vortex of the model in 2D square lattice. In the topological phase, the Majorana zero mode is localized at the vortex core and its spin component at the center is completely parallel to the external magnetic field, which leads to spin selective Andreev reflection. In the topological trivial phase, there is no Majorana zero mode inside the vortex. However, the spin-orbit coupling induces a spin selective Andreev reflection at the bias of the lowest quasiparticle energy. The Majorana zero mode induced spin selective ScIenTIfIc RepoRtS | (2018) 8:7853 | DOI:10.1038/s41598-018-26184-z Andreev reflection is robust and gives a quantized value of differential conductance 2e 2 /h, which is independent of the tunneling barrier. The usual vortex quasiparticle induced spin selective Andreev reflection gives a vanishing value of the differential conductance at zero bias due to quantum destructive interference and is sensitive to the barrier in the tunneling.

Methods
Transport Methods for dI/dV of N/S Junction. Then, to calculate the differential conductance in Eq. (14), what we need is the S-matrix for N/S junction. From the Fisher-Lee Formula 50 , the whole 4N × 4N S-matrix can be calculated R R † , Σ R is the self-energy induced by the lead. It's a 4N × 4N matrix but only nonzero for the probed site p on the device, that we will show later, as well as Γ. And the total (retarded) Green's function reads Only the 4 × 4 sub-matrix G R (p, p) contributes to the calculation because of the sparsity of Γ. Infinitesimal positive number η D = 10 −5 is adopted. To calculate the inverse of large sparse matrix, we use Intel MKL PARDISO Solver in Fortran code.
As we know, the lead will contribute a self-energy 43 to the device, which could be described by a 4N × 4N matrix Σ R , η L is a infinitesimal positive number, we used η L = 10 −5 . Although the lead's Hamiltonian has infinite dimension in matrix form, same as the Green's function, the only submatrix we need in the calculation is the Green's function on the lead's vertex site 0, because of the sparsity of the coupling matrix τ. And the surface Green's function G (0, 0) L R can be calculated by decimation method through Eq. (27). In order to get the surface Green's function, let's write Eq. (27) in block matrix form (2 ) ( ) and the effective interaction between adjacent sites A = diag{−t L ,−t L , t L , t L }. Note that the only value that we need is g G (0, 0) If we discard all even number for n and m, and take the transform  Repeat these steps in Eq. (32), we will update the coefficient and abandon the Green's function between nearby sites consistently.
After sufficient number of iterations, the coefficient A becomes the effective interaction between pretty far sites which must be a sub-matrix comprised of small values. From d′′g 11 −A′′g′′ 21 = 1, we finally obtain the surface Green's function for lead Apply sparse matrix to Transport calculation. In the calculation of Eq. (25), it's difficult to take the inverse directly, since Hamiltonian of SLTD model device  D in dominator is a big matrix of size 4N × 4N. Traditionally, the methods of recursive Green's function could be applied here to calculate the Green's function on the contact point p, G R (p, p), for S-matrix calculating Eq. (24). There is an alternative approach using sparse matrix. In this regard, we store  D as sparse matrix, which is suitable for tight-binding lattice model. Then the j-th column of total Green's function, defined as G j R , can be solved by the linear equation and G R (p, p) can be drawn from relevant columns of total Green's function.
Intel MKL PARDISO Solver could be used to solve the linear equation of sparse matrix Eq. (36). We also tried methods of recursive Green's function, by calculating Green's function of one lattice's column gradually. Same result of the differential conductance [ Fig. 6(a)] was obtained, but it cost more time in calculation. So taking inverse of sparse matrix directly shows its efficiency advantage.
ScIenTIfIc RepoRtS | (2018) 8:7853 | DOI:10.1038/s41598-018-26184-z to the spin property of this first vortex excitation, we see it is spin polarized down at vortex core, i.e., u ↓ ≠ 0 and u ↑ = v ↑ = v ↓ = 0. It is consistent with the calculation in ref. 15,34 Appendix B: Normal Andreev Reflection. For normal Andreev reflection, incident electrons are reflected by SC device as holes with opposite spin direction. Though electrons with energy less then Δ reflect as holes, this process will be suppressed by sufficiently small t c (high barrier) unlike SESAR induced by MZM. For comparison, we calculated the Andreev reflection coefficient in normal Andreev reflection case. The Andreev reflection coefficient T A here is defined by A h e he † As showed in Fig. 10, for tiny coupling coefficient t c , the Andreev reflection coefficient T A vanished for all energy E lower than SC gap Δ. Since t c increasing in certain range, the T A that the reflection in SC gap contributed rises, and finally form a plateau valued 2 for two channels spin up and down. These behaves of Andreev reflection coefficient are consistent with the BTK theory.