\documentclass[aps,prl,onecolumn]{revtex4}
\usepackage{bm,color,amsmath,amssymb,mathrsfs,latexsym,graphicx,psfrag,appendix}

% simplifies using the up and down arrows to denote spin
\newcommand{\up}{\uparrow}
\newcommand{\dw}{\downarrow}

% Roman functions for real and imaginary parts
\newcommand{\re}{\mathrm{Re}}
\newcommand{\im}{\mathrm{Im}}

\newcommand{\Hamil}{{\cal H}}
\newcommand{\bk}{{\mathbf k}}
\newcommand{\bq}{{\mathbf q}}
\newcommand{\br}{{\mathbf{r}}}

\def\be{\begin{equation}}
\def\ee{\end{equation}}
\def\bea{\begin{eqnarray}}
\def\eea{\end{eqnarray}}
\def\bD{{\mathbf{D}}}
\def\tr{\textrm{tr}}

\begin{document}

\appendix
\section{Proof that $|\tr(\mathcal{M})|$ equals the number of MZM's}
Using $\{\mathcal{H},\mathcal{M}\}=0$, we know that for any given eigenstate $|u\rangle$ of $M$ with eigenvalue $\pm1$, $\mathcal{H}|u\rangle$, if not a null vector, must be an eigenstate of $M$ with the eigenvalue $\mp1$. Denote by $N_\pm$ the number of eigenvalues of $M$ that have eigenvalues $+1$ and $-1$ respectively, and assume that $N_+>N_-$. The corresponding eigenstates are denoted by $|u_i^\pm\rangle$, forming subspaces $\Psi_\pm$. Obviously, since any nonzero $H|u_i^+\rangle$ for $i=1,...,N_+$ is a state in $\Psi_-$, there must be at least $N_+-N_-$ of them that are null vectors, i.e., $N_+-N_-$ eigenstates of $\mathcal{H}$ of zero eigenvalue. The same argument proceeds for the case $N_->N_+$, leading to the result that there must be at least $N_--N_+$ eigenstates of $\mathcal{H}$ of zero eigenvalue.

\section{MZM's in a general TCI having mirror Chern number $C_m$ with induced superconductivity}

In the main text, we explicitly show that on the $(001)$-surface of mirror TCI SnTe, there are two MZM's protected by $M_T$. Here we extend to all TCI's having (i) mirror Chern number $C_m$, (ii) TRS, (iii) Cooper pairing that preserves the above two symmetries and (iv) same sign of pairing amplitude on all pieces of the Fermi surface. All these make the system have both $M$ and $T$, so that when there is vortex, their combination $M_T$ is preserved, while they are no longer symmetries separately.

Given a mirror plane, a 3D BZ may have one (e.g., FCC) or two (e.g., simple cubic) mirror symmetric subspaces. For a surface termination that preserves the mirror symmetry, i.e., perpendicular to the mirror plane(s), there are one or two mirror symmetric lines in its SBZ. On a mirror symmetric line, bands with opposite mirror eigenvalues ($\pm{i}$ in spinful systems) can cross each other, leaving a set of Dirac points in the SBZ. On a mirror symmetric line in SBZ, TRS sends a right going mode to a left going mode while changing the mirror eigenvalue, making the number of right/left going modes with eigenvalue $+i$ exactly equals the number of left/right going modes with eigenvalue $-i$.

We now make an assumption that there is only one, instead of two, mirror symmetric plane in 3D BZ, or one mirror symmetric line in the SBZ. We will relax the assumption at a later stage. The mirror Chern number is defined as $C_m=n^+_R-n^+_L+n^-_L-n^-_R$, where $n^{+/-}_{L/R}$ is the number of left/right going modes with eigenvalue $+i/-i$. We can adiabatically tune the system such that (i) if $C_m>0$, all right going modes have $M=+i$ and all left going ones have $-i$ or (ii) if $C_m<0$, all right going modes have $M=-i$ and all left going ones have $+i$. Along the mirror symmetric line in SBZ, the Dirac points are symmetric about the origin. A Dirac point can be located either at a TRS momentum or a generic point, and we denote it by $D_0$ and $D$, respectively.

The pairing amplitude at $\bk$ is defined as
\bea
\delta(\bk)=\langle{\Omega}|\hat{\Delta}\psi(\bk)\hat{T}\psi(\bk)\hat{T}^{-1}|\Omega\rangle,
\eea
where $\hat\Delta\equiv\tilde\Delta_{\alpha\beta}c^\dag_\alpha(\bk)c^\dag_\beta(-\bk)+h.c.$ is the pairing operator, $\psi(\bk)$ is the eigenstate annihilation operator in the \emph{normal state} at $\bk$, and $|\Omega\rangle$ is the Fermi liquid ground state (non-superconducting). This amplitude is well defined for any non-degenerate $\bk$-point. Using TRS\cite{Qi2010}, it can be proved that $\delta(\bk)=\delta(-\bk)\in{real}$. To the lowest order of pairing strength, the Bougliubov excitation at $\bk$ has energy $\sqrt{(\epsilon(\bk)-\mu)^2+\delta(\bk)^2}$. This means $\delta(\bk)$ cannot change sign on any connected FS in a gapped superconductor. At the same time, if there are multiple disjoint FS's, then two FS's may have opposite signs, unless they are related to each other by TRS. Sign changes between disjoint FS's usually implies unconventional pairing mechanisms. The 1111-family of iron-based superconductors belongs to this special class. For simplicity, let us for now assume that all FS's have the same sign of pairing amplitude, and the generalization to the generic case is made in Sec.\ref{sec:total_number}.

\subsection{Surface Hamiltonians in the normal state}
We first write down the effective theory ($k\cdot{p}$-model) for each $D_0$ and $D$ without superconductivity. For $D_0$, the little group is generated by $M$ and $T$, represented by $i\sigma_y$ and $K(i\sigma_y)$. So the effective theory reads:
\bea
h_0(\bq)=-\mu\sigma_0+v_1\sigma_yq_x+(v_2\sigma_x+v_3\sigma_z)q_x+O(|q|^2).
\eea
At $D$, the little group is generated by only $M=i\sigma_y$, so
\bea
h(\bq)=(v_0q_x-\mu)\sigma_0+v_1\sigma_yq_x+(v_2\sigma_x+v_3\sigma_z)q_x+O(|q|^2).
\eea

\subsection{Bulk superconductivity for a Dirac cone centered at TRIM}

Next we consider a homogeneous onsite BCS pairing. Since $D_0$ is a TRS momentum, the Cooper pairing is intra-pocket. The BdG Hamiltonian in general looks
\bea\label{eq:BdG_H0}
H_0(\bq)=\left(\begin{matrix}
h_0(\bq) & \tilde\Delta\\
\tilde\Delta^\dag & -h^T_0(-\bq)
\end{matrix}\right).
\eea
The form of $\tilde\Delta$ is determined by the symmetry constraints, namely,
\bea
(i\sigma_y)^T\tilde\Delta(i\sigma_y)&=&\tilde\Delta\textrm{ from mirror symmetry and}\\
\nonumber
(i\sigma_y)^T\tilde\Delta^\ast(i\sigma_y)&=&\tilde\Delta\textrm{ from TRS}.
\eea
These constraints require
\bea
\tilde\Delta=\Delta_2\sigma_0+i\Delta_1\sigma_y,
\eea
where $\Delta_{1,2}\in{Real}$. Diagonalizing the homogeneous Hamiltonian exactly, we obtain
\bea\label{eq:h0_dispersion}
E(\bq)=\pm\sqrt{v_1^2q_x^2+(v^2_2+v_3^2)q_y^2+\Delta_1^2+\Delta_2^2+\mu^2\pm2\sqrt{q_x^2v_1^2\mu^2+q_y^2(v_2^2+v_3^2)(\Delta_2^2+\mu^2)}}.
\eea
It is straightforward to show that the dispersion in Eq.(\ref{eq:h0_dispersion}) is always gapped as far as $\Delta_1^2+\Delta_2^2\neq0$. This means the surface state around $D_0$ is always gapped, against any parameter change. On this gapped FS, we can calculate the sign of the pairing amplitude: $\textrm{sign}(\delta(\bq))=\textrm{sign}(\Delta_1)$.

\subsection{Bulk superconductivity for two Dirac cones centered at $\pm\bD$}

For a Dirac cone around $D$, there must be another around $-D$ due to TRS. The gauge at $D$ is chosen such that
\bea
\hat{M}f_{+,\tau}(\bq)\hat{M}^{-1}=(i\sigma_y)_{\tau\tau'}f_{+,\tau'}(M\bq),
\eea
where $f_{+,\tau}(\bq)$ is the electron annihilation operator at $\bD+\bq$. We can thus fix the gauge at $-\bD-\bq$ by TRS symmetry, i.e., $f_{-,\tau}(\bq)\equiv\hat{T}f_{+,\tau}(-\bq)\hat{T}^{-1}$. Since $[\hat{M},\hat{T}]=0$, at $-\bD$ ($\bq=0$), we also have the representation of the mirror symmetry as $M=i\sigma_y$. In this gauge, the $k\cdot{p}$-model for the cone centered at $-\bD$ is
\bea
h_-(\bq)=-(v_0q_x-\mu)\sigma_0+v_1q_x\sigma_y-(v_2\sigma_x+v_3\sigma_z)q_y.
\eea
We have assumed that Cooper pairs have zero momentum, so the pairing is only between $f_+(\bq)$ and $f_-(-\bq)$, leading to the following BdG Hamiltonian:
\bea\label{eq:BdG_H}
H(\bq)=\left(\begin{matrix}
h_+(\bq) & 0 & 0 & \tilde\Delta\\
0 & h_-(\bq) & -\tilde\Delta^T & 0\\
0 & -\tilde\Delta^\ast & -h^T_+(-\bq) & 0\\
\tilde\Delta^\dag & 0 & 0 & -h^T_-(-\bq)
\end{matrix}\right).
\eea
Eq.(\ref{eq:BdG_H}) is in a block-diagonalized form, and the pairing matrix $\tilde\Delta$ is constrained by
\bea
(i\sigma_y)^T\tilde\Delta(i\sigma_y)&=&\tilde\Delta,\\
\nonumber
\tilde\Delta^\dag&=&\tilde\Delta.
\eea
The only possible form of $\tilde\Delta$ is then $\tilde\Delta=\Delta_1\sigma_0+\Delta_2\sigma_y$, where $\Delta_{1,2}\in{real}$. Eq.(\ref{eq:BdG_H}) can also be exactly diagonalized, giving the dispersion:
\bea\nonumber
E^2(\bq)=v_1^2q_x^2+(v_2^2+v_3^2)q_y^2+(\mu\pm{}v_0q_x)^2+\Delta_1^2+\Delta_2^2
\pm2\sqrt{[\pm{}q_xv_1(\pm{}v_0q_x-\mu)+\Delta_1\Delta_2]^2+q_y^2(v_2^2+v_3^2)[(v_0q_x\pm\mu)^2+\Delta_2^2]}.
\eea
Tedious algebra shows that when $|\Delta_1|>|\Delta_2|\ge0$, the spectrum is fully gapped and when $|\Delta_2|>|\Delta_1|\ge0$, we have four band crossings at zero energy, namely: 
\bea\nonumber
(q_x,q_y)=(\pm\frac{\Delta_1\mu}{v_0\Delta_1-v_1\Delta_2},\pm\sqrt{\frac{(\Delta_2^2-\Delta_1^2)[(v_0\Delta_1-v_1\Delta_2)^2+\mu^2v_1^2]}{|v_0\Delta_1-v_1\Delta_2|(v_2^2+v_3^2)}}).
\eea
This surface superconducting phase has four nodes in four quadrants respectively, and two nodes related by mirror symmetry have opposite winding numbers of $d$-vectors leaving the total winding number zero. But this is not our current interest which is fully gapped surface superconductivity. Therefore, we take $|\Delta_1|>|\Delta_2|$, and the dispersion is always gapped against any change in parameters. We can also calculate the sign of the pairing amplitude, obtaining $\textrm{sign}(\delta(\bq))=\textrm{sign}(\Delta_1)$.

\subsection{Superconducting vortex bound states for a Dirac cone centered at TRIM}

We are now ready to consider a vortex in the Cooper pairing. In principle, a vortex in real space will induce Cooper pairing between Dirac cones that are not opposite to each other in the momentum space, but we will for now ignore this effect, which means that $r_0\delta{D}\gg1$, where $r_0$ is the size of the vortex and $\delta{D}$ is the distance in $k$-space between nearest Dirac cones. After we calculate, under this working assumption, the number of symmetry protected MZM's, we then can argue that so far as the bulk gap is open, this symmetry protected quantum number remains the same under any perturbation that we ignore during the calculation, although the explicit wavefunctions of the MZM's generically vary.

First we study the surface states around $\bD_0$. The vortex is introduced by replacing $\Delta_{1,2}$ in Eq.(\ref{eq:BdG_H0}) by $\Delta_{1,2}(r)e^{-i\theta}$, where $\Delta_{1,2}(0)=0$ and $\Delta(\infty)=\Delta_{1,2}>0$.
Since $H_0$ is gapped for any parameter set, we can always adiabatically tune the parameters to $v_1=-v_2\equiv{v}$, $v_0=v_3=\mu=0$ and $\Delta_2=0$. The number of MZM's for this set of parameters is the same as the number of MZM's for the original parameters, since the surface remains fully gapped during the change. Solving the Schrodinger equation
\bea
\left(\begin{matrix}
0 & -ve^{-i\theta}(\partial_r-\frac{i\partial_\theta}{r}) & 0 & \Delta_1(r)e^{-i\theta}\\
ve^{i\theta}(\partial_r+\frac{i\partial_\theta}{r}) & 0 & -\Delta_1(r)e^{-i\theta} & 0\\
0 & -\Delta_1(r)e^{i\theta} & 0 & ve^{i\theta}(\partial_r+\frac{i\partial_\theta}{r})\\
\Delta_1(r)e^{i\theta} & 0 & -ve^{-i\theta}(\partial_r-\frac{i\partial_\theta}{r}) & 0
\end{matrix}\right)
\left(\begin{matrix}
\psi_1\\
\psi_2\\
\psi_1^\ast\\
\psi_2^\ast
\end{matrix}\right)=0,
\eea
we obtain one MZM mode
\bea
\gamma&=&i\int{d^2\br}[f^\dag_\dw(\br)-f_\dw(\br)]e^{-\int_0^r|\frac{\Delta_1(r')}{v}|dr'},\;\textrm{if}\;v\Delta_1>0;\\
\gamma&=&\int{d^2\br}[f^\dag_\dw(\br)+f_\dw(\br)]e^{-\int_0^r|\frac{\Delta_1(r')}{v}|dr'},\;\textrm{if}\;v\Delta_1<0.
\eea
Using $T=i\sigma_yK$ and $M=i\sigma_y$, it is easy to verify that 
\bea\label{eq:TRIM}
\hat{M}_T\gamma\hat{M}^{-1}_T=\textrm{sign}(v_1\Delta_1)\gamma.
\eea

\subsection{Superconducting vortex bound states for Dirac cones centered at $\pm\bD$}

For the Cooper pairing between the surface electrons around $\bD$ and $-\bD$, the Schrodinger equation can be block diagonalized into two equations, namely
\bea
\left(\begin{matrix}
0 & -ve^{-i\theta}(\partial_r-\frac{i\partial_\theta}{r}) & \Delta_1(r)e^{-i\theta} & 0\\
ve^{i\theta}(\partial_r+\frac{i\partial_\theta}{r}) & 0 & 0 & \Delta_1(r)e^{-i\theta}\\
\Delta_1(r)e^{i\theta} & 0 & 0 & ve^{i\theta}(\partial_r-\frac{i\partial_\theta}{r})\\
0 & \Delta_1(r)e^{i\theta} & -ve^{-i\theta}(\partial_r-\frac{i\partial_\theta}{r}) & 0
\end{matrix}\right)
\left(\begin{matrix}
\psi_1\\
\psi_2\\
\psi_3\\
\psi_4
\end{matrix}\right)=0,
\eea
and
\bea
\left(\begin{matrix}
0 & -ve^{i\theta}(\partial_r+\frac{i\partial_\theta}{r}) & -\Delta_1(r)e^{-i\theta} & 0\\
ve^{-i\theta}(\partial_r-\frac{i\partial_\theta}{r}) & 0 & 0 & -\Delta_1(r)e^{-i\theta}\\
-\Delta_1(r)e^{i\theta} & 0 & 0 & ve^{i\theta}(\partial_r+\frac{i\partial_\theta}{r})\\
0 & -\Delta_1(r)e^{i\theta} & -ve^{-i\theta}(\partial_r-\frac{i\partial_\theta}{r}) & 0
\end{matrix}\right)
\left(\begin{matrix}
\psi_3^\ast\\
\psi_4^\ast\\
\psi_1^\ast\\
\phi_2^\ast
\end{matrix}\right)=0,
\eea
where we have used that MZM must be invariant under PHS. Also, since $M_T$ is a symmetry, the combination $P\times{M}_T$ is also a symmetry that anti-commutes with the Hamiltonian. A MZM must also be an eigenstate of $P\times{M}_T$ with eigenvalue being either $+1$ or $-1$. This places constraints on $\psi_{1,2,3,4}$, namely, $\psi_1=\pm\psi_4$ and $\psi_2=\mp\psi_3$. Solving the equations analytically, we have two MZM solutions:
\bea
\gamma_1&=&\int{d}r^2[f^\dag_{+,\dw}(\br)-f_{-,\up}(\br)+h.c.]e^{-\int_0^r|\frac{\Delta_1(r')}{v}|dr'},\\
\nonumber
\gamma_2&=&i\int{d}r^2[f^\dag_{+,\dw}(\br)-f_{-,\up}(\br)-h.c.]e^{-\int_0^r|\frac{\Delta_1(r')}{v}|dr'},
\eea
if $\textrm{sign}(v\Delta_1)>0$ and
\bea
\gamma_1&=&\int{d}r^2[f^\dag_{+,\dw}(\br)+f_{-,\up}(\br)+h.c.]e^{-\int_0^r|\frac{\Delta_1(r')}{v}|dr'},\\
\nonumber
\gamma_2&=&i\int{d}r^2[f^\dag_{+,\dw}(\br)+f_{-,\up}(\br)-h.c.]e^{-\int_0^r|\frac{\Delta_1(r')}{v}|dr'},
\eea
if $\textrm{sign}(v\Delta_1)<0$. It is straightforward to verify:
\bea\label{eq:nonTRIM}
\hat{M}_T\gamma_{1,2}\hat{M}^{-1}_T=\textrm{sign}(v\Delta_1)\gamma_{1,2}.
\eea

\subsection{Total number of MZM's}
\label{sec:total_number}
From the analysis, we see that there is one MZM from a Dirac cone centered at a TRIM, and two MZM's from a pair of Dirac cones centered at $\pm\bD$. Due to our assumption that there is no sign change between pairing amplitudes on different FS's, $\Delta_1$ has the same sign for all Dirac cones, therefore, every MZM transforms under $M_T$ as
\bea\label{eq:temp1}
\hat{M}_T\gamma_{i=1,...,|C_m|}\hat{M}^{-1}_T=\textrm{sign}(C_m\Delta_1).
\eea
In the basis of $\gamma_i$'s, the symmetry operation $M_T$ is simply $\mathcal{M}=\textrm{sign}(C_m\Delta_1)I_{|C_m|\times|C_m|}$, so $|\tr(\mathcal{M})|=|C_m|$. Therefore, there are exactly $|C_m|$ MZM's protected by $M_T$.

If there are two mirror symmetric planes in the 3D BZ, then there are two mirror Chern numbers denoted by $C_{m0}$ and $C_{m\pi}$. There are two mirror symmetric lines in the SBZ, along which there are $|C_{m0}|$ and $|C_{m\pi}|$ Dirac points, respectively. An analysis same as the above and the result in Eq.(\ref{eq:temp1}) apply to Dirac cones on both lines. The total number of protected MZM's is hence
\bea
N_{MZM}=|C_{m0}+C_{m\pi}|.
\eea

It is easy to generalize to cases where the different surface Fermi pockets have different signs of pairing amplitude. There are in total $|C_m|$ Fermi pockets on the surface, and the pairing signs are symmetric about the origin due to TRS. The total number of MZM is simply the sum of the signs of all the pockets, namely
\bea
N_{MZM}=|\sum_{m}\textrm{sign}(\Delta_m)|,
\eea
where $\textrm{sign}(\Delta_m)$ is the pairing sign of the $m$-th Fermi pocket.


\section{The two `semi'-Majorana modes in pure samples}

Out of the four MZM's in Eq.(10) of the main text, two can be gapped out by the following perturbation
\bea\label{eq:perturbation}
\delta{H}&=&i\lambda(\gamma_1\gamma_2+\gamma_2\gamma_3+\gamma_3\gamma_4+\gamma_4\gamma_1)\\
\nonumber
&=&i\lambda(\gamma_1-\gamma_3)(\gamma_2-\gamma_4).
\eea
From Eq.(10) in the main text, we can see that $\gamma_1-\gamma_3$ is composed of electrons from Dirac cones centered at $\bD_{1,3}$, and $\gamma_2-\gamma_4$ is composed of electrons from Dirac cones centered at $\bD_{2,4}$. Since $\bD_{1,3}\pm\bD_{2,4}\sim(\pi,\pi)$, all terms in Eq.(\ref{eq:perturbation}) carry large momentum scattering or large momentum pairing. Large momentum scattering is suppressed if the impurity potential is smooth; large momentum pairing is suppressed if the vortex size is much larger than the lattice constant. When both effects are suppressed, the coefficient $\lambda$ in Eq.(\ref{eq:perturbation}) must be very small. Since Eq.(\ref{eq:perturbation}) is the only allowed hybridization term, there are four, instead of two MZM's.
\end{document}