Topological properties of multilayers and surface steps in the SnTe material class

TL;DR

This paper explores the topological properties of surface steps in SnTe multilayers, revealing how symmetries and invariants influence low-energy states and potential instabilities like magnetization, with implications for experimental observations.

Contribution

It identifies the topological invariants governing step states in SnTe multilayers and demonstrates their dependence on layer number and symmetry protections.

Findings

Step states occur when mirror- and spin-resolved Chern numbers differ across steps.

Particle-hole symmetry can protect Weyl points at surface steps.

Magnetic domain walls support low-energy bound states due to topological distinctions.

Abstract

Surfaces of multilayer semiconductors typically have regions of atomically flat terraces separated by atom-high steps. Here we investigate the properties of the low-energy states appearing at the surface atomic steps in \snte. We identify the important approximate symmetries and use them to construct relevant topological invariants. We calculate the dependence of mirror- and spin-resolved Chern numbers on the number of layers and show that the step states appear when these invariants are different on the two sides of the step. Moreover, we find that a particle-hole symmetry can protect one-dimensional Weyl points at the steps. Since the local density of states is large at the step the system is susceptible to different types of instabilities, and we consider an easy-axis magnetization as one realistic possibility. We show that magnetic domain walls support low-energy bound states because the regions with opposite magnetization are topologically distinct in the presence of non-symmorphic chiral and mirror symmetries, providing a possible explanation for the zero-bias conductance peak observed in the recent experiment [Mazur et al., Phys. Rev. B 100, 041408(R) (2019)]

Figures

Tables

Table 1. Table 1: Mirror Chern numbers C + subscript 𝐶 C_{+} for systems with odd number of layers N 𝑁 N ( C + = 0 subscript 𝐶 0 C_{+}=0 for even N 𝑁 N ). Parameters are m = 1.65 𝑚 1.65 m=1.65 eV, t 12 = 0.9 subscript 𝑡 12 0.9 t_{12}=0.9 eV, t 11 = 0.5 subscript 𝑡 11 0.5 t_{11}=0.5 eV , λ x = λ y = λ z = 0.3 subscript 𝜆 𝑥 subscript 𝜆 𝑦 subscript 𝜆 𝑧 0.3 \lambda_{x}=\lambda_{y}=\lambda_{z}=0.3 eV.

$N$	$C_{+}$
1	0
3	-2
5	2
7	2
9	-2
11	-2
13	2
15	2
17	-2
19	-2

Equations

$$\mathcal{H}(\mathbf{k})$$

$$H$$

$$\mathcal{H}(\mathbf{k})$$

$$h_x^{(1)}(k_x)=\left(\begin{array}{cccccccc}0& 1+e^{-i{k}_x}& 0& 0& 0& 0& 0& 0\\ 1+e^{i{k}_x}& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 1+e^{i{k}_x}& 0& 0& 0& 0\\ 0& 0& 1+e^{-i{k}_x}& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1+e^{-i{k}_x}\\ 0& 0& 0& 0& 0& 0& 1+e^{-i{k}_x}& 0\\ 0& 0& 0& 0& 0& 1+e^{i{k}_x}& 0& 0\\ 0& 0& 0& 0& 1+e^{i{k}_x}& 0& 0& 0\end{array}\right),$$

$$h_y^{(1)}(k_y)=\left(\begin{array}{cccccccc}0& 0& 0& 1+e^{-i{k}_y}& 0& 0& 0& 0\\ 0& 0& 1+e^{-i{k}_y}& 0& 0& 0& 0& 0\\ 0& 1+e^{i{k}_y}& 0& 0& 0& 0& 0& 0\\ 1+e^{i{k}_y}& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 1+e^{i{k}_y}& 0& 0\\ 0& 0& 0& 0& 1+e^{-i{k}_y}& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1+e^{-i{k}_y}\\ 0& 0& 0& 0& 0& 0& 1+e^{i{k}_y}& 0\end{array}\right),$$

$$h_z^{(1)}(k_z)=\left(\begin{array}{cccccccc}0& 0& 0& 0& 0& 1+e^{-i{k}_z}& 0& 0\\ 0& 0& 0& 0& 0& 0& 1+e^{-i{k}_z}& 0\\ 0& 0& 0& 0& 0& 0& 0& 1+e^{-i{k}_z}\\ 0& 0& 0& 0& 1+e^{-i{k}_z}& 0& 0& 0\\ 0& 0& 0& 1+e^{i{k}_z}& 0& 0& 0& 0\\ 1+e^{i{k}_z}& 0& 0& 0& 0& 0& 0& 0\\ 0& 1+e^{i{k}_z}& 0& 0& 0& 0& 0& 0\\ 0& 0& 1+e^{i{k}_z}& 0& 0& 0& 0& 0\end{array}\right),$$

$$h_{x,y}^{(2)}=\left(\begin{array}{cccccccc}0& 0& 1+e^{-i(k_x+k_y)}& 0& 0& 0& 0& 0\\ 0& 0& 0& e^{i{k}_x}+e^{-i{k}_y}& 0& 0& 0& 0\\ 1+e^{i(k_x+k_y)}& 0& 0& 0& 0& 0& 0& 0\\ 0& e^{-i{k}_x}+e^{i{k}_y}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& e^{-i{k}_x}+e^{i{k}_y}& 0\\ 0& 0& 0& 0& 0& 0& 0& 1+e^{-i(k_x+k_y)}\\ 0& 0& 0& 0& e^{i{k}_x}+e^{-i{k}_y}& 0& 0& 0\\ 0& 0& 0& 0& 0& 1+e^{i(k_x+k_y)}& 0& 0\end{array}\right),$$

$$h_{y,x}^{(2)}=\left(\begin{array}{cccccccc}0& 0& e^{-i{k}_x}+e^{-i{k}_y}& 0& 0& 0& 0& 0\\ 0& 0& 0& 1+e^{i(k_x-k_y)}& 0& 0& 0& 0\\ e^{i{k}_x}+e^{i{k}_y}& 0& 0& 0& 0& 0& 0& 0\\ 0& 1+e^{i(k_y-k_x)}& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 1+e^{i(k_y-k_x)}& 0\\ 0& 0& 0& 0& 0& 0& 0& e^{-i{k}_x}+e^{-i{k}_y}\\ 0& 0& 0& 0& 1+e^{i(k_x-k_y)}& 0& 0& 0\\ 0& 0& 0& 0& 0& e^{i{k}_x}+e^{i{k}_y}& 0& 0\end{array}\right),$$

$$h_{y,z}^{(2)}=\left(\begin{array}{cccccccc}0& 0& 0& 0& 1+e^{-i(k_y+k_z)}& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1+e^{-i(k_y+k_z)}\\ 0& 0& 0& 0& 0& 0& e^{i{k}_y}+e^{-i{k}_z}& 0\\ 0& 0& 0& 0& 0& e^{i{k}_y}+e^{-i{k}_z}& 0& 0\\ 1+e^{i(k_y+k_z)}& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& e^{-i{k}_y}+e^{i{k}_z}& 0& 0& 0& 0\\ 0& 0& e^{-i{k}_y}+e^{i{k}_z}& 0& 0& 0& 0& 0\\ 0& 1+e^{i(k_y+k_z)}& 0& 0& 0& 0& 0& 0\end{array}\right),$$

$$h_{z,y}^{(2)}=\left(\begin{array}{cccccccc}0& 0& 0& 0& e^{-i{k}_y}+e^{-i{k}_z}& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& e^{-i{k}_y}+e^{-i{k}_z}\\ 0& 0& 0& 0& 0& 0& 1+e^{i(k_y-k_z)}& 0\\ 0& 0& 0& 0& 0& 1+e^{i(k_y-k_z)}& 0& 0\\ e^{i{k}_y}+e^{i{k}_z}& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 1+e^{i(k_z-k_y)}& 0& 0& 0& 0\\ 0& 0& 1+e^{i(k_z-k_y)}& 0& 0& 0& 0& 0\\ 0& e^{i{k}_y}+e^{i{k}_z}& 0& 0& 0& 0& 0& 0\end{array}\right),$$

$$h_{x,z}^{(2)}=\left(\begin{array}{cccccccc}0& 0& 0& 0& 0& 0& 1+e^{-i(k_x+k_z)}& 0\\ 0& 0& 0& 0& 0& e^{i{k}_x}+e^{-i{k}_z}& 0& 0\\ 0& 0& 0& 0& e^{i{k}_x}+e^{-i{k}_z}& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1+e^{-i(k_x+k_z)}\\ 0& 0& e^{-i{k}_x}+e^{i{k}_z}& 0& 0& 0& 0& 0\\ 0& e^{-i{k}_x}+e^{i{k}_z}& 0& 0& 0& 0& 0& 0\\ 1+e^{i(k_x+k_z)}& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 1+e^{i(k_x+k_z)}& 0& 0& 0& 0\end{array}\right),$$

$$h_{z,x}^{(2)}=\left(\begin{array}{cccccccc}0& 0& 0& 0& 0& 0& e^{-i{k}_x}+e^{-i{k}_z}& 0\\ 0& 0& 0& 0& 0& 1+e^{i(k_x-k_z)}& 0& 0\\ 0& 0& 0& 0& 1+e^{i(k_x-k_z)}& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& e^{-i{k}_x}+e^{-i{k}_z}\\ 0& 0& 1+e^{i(k_z-k_x)}& 0& 0& 0& 0& 0\\ 0& 1+e^{i(k_z-k_x)}& 0& 0& 0& 0& 0& 0\\ e^{i{k}_x}+e^{i{k}_z}& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& e^{i{k}_x}+e^{i{k}_z}& 0& 0& 0& 0\end{array}\right).$$

$$m_x=\left(\begin{array}{cccccccc}0& 1& 0& 0& 0& 0& 0& 0\\ 1& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0& 0\\ 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1\\ 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0\end{array}\right),$$

$$m_z=\left(\begin{array}{cccccccc}0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 1\\ 0& 0& 1& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0& 0& 0\\ 1& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0& 0\end{array}\right)$$

$$r_z=\left(\begin{array}{cccccccc}0& 0& 1& 0& 0& 0& 0& 0\\ 0& 0& 0& 1& 0& 0& 0& 0\\ 1& 0& 0& 0& 0& 0& 0& 0\\ 0& 1& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 0& 0& 1\\ 0& 0& 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 0& 0& 1& 0& 0\end{array}\right).$$

$${\tilde{\mathcal{H}}}_N(k_x,k_y)=U(k_x{)}^{†}{\mathcal{H}}_N(k_x,k_y)U(k_x),\tilde{\mathcal{T}}=\mathcal{K}U(k_x{)}^T({\sigma }_y\otimes \text{mathbbm}{1}_3\otimes \text{mathbbm}{1}_8)U(k_x).$$

$$\chi =U(k_x+2\pi {)}^{†}U(k_x),\chi {\tilde{\mathcal{H}}}_N(k_x,k_y){\chi }^{†}={\tilde{\mathcal{H}}}_N(k_x+2\pi ,k_y).$$

$${\tilde{\mathcal{T}}}_2=\tilde{\mathcal{T}}\chi ,{\tilde{\mathcal{T}}}_2{\tilde{\mathcal{H}}}_N(k_x,k_y){\tilde{\mathcal{T}}}_2^{-1}={\tilde{\mathcal{H}}}_N(-k_x-2\pi ,-k_y),$$

$${\vec{A}}_n=i⟨n,\vec{k}|{\vec{\nabla }}_{\vec{k}}|n,\vec{k}⟩,$$

$${\Omega }_n={\vec{\nabla }}_{\vec{k}}\times \overrightarrow{A_n}.$$

$$\mathcal{C}=\frac{1}{2\pi }\underset{BZ}{\int }{d}^{2}k\sum _{n\le n_F}{\Omega }_n.$$

$$\mathcal{C}=\frac{1}{\pi }\underset{BZ}{\int }{d}^{2}k\sum _{\genfrac{}{}{0ex}{}{n\le n_F}{n^{\prime }>n_F}}\mathfrak{Im}\frac{⟨n,\vec{k}|{\partial }_{k_x}\mathcal{H}|n^{\prime },\vec{k}⟩⟨n^{\prime },\vec{k}|{\partial }_{k_y}\mathcal{H}|n,\vec{k}⟩}{{\left(E_{\vec{k}}^{(n)}-E_{\vec{k}}^{(n^{\prime })}\right)}^2},$$

$${\mathcal{H}}_{\mathrm{step}}(k)=\Delta {\sigma }_0{\tau }_x+vk{\sigma }_z{\tau }_0+\tilde{v}k{\sigma }_0{\tau }_z+\Lambda {\sigma }_y{\tau }_z+\tilde{\Lambda }k{\sigma }_x{\tau }_y,$$

$$S_z{\mathcal{H}}_{\mathrm{step}}(k)S_z={\mathcal{H}}_{\mathrm{step}}(k),S_z={\sigma }_z{\tau }_0\mathrm{for}\mathrm{\Lambda }=\tilde{\mathrm{\Lambda }}=0,$$

$$C{\mathcal{H}}_{\mathrm{step}}(k)C=-{\mathcal{H}}_{\mathrm{step}}(k),C={\sigma }_y{\tau }_y$$

$$M_x{\mathcal{H}}_{\mathrm{step}}(-k)M_x={\mathcal{H}}_{\mathrm{step}}(k),M_x={\sigma }_x{\tau }_x.$$

$${\sigma }_y{\tau }_x{\mathcal{H}}_{\mathrm{step}}^{\ast }(-k){\sigma }_y{\tau }_x={\mathcal{H}}_{\mathrm{step}}(k).$$

$${\mathcal{H}}_{\mathrm{step}}^Z=h_z{\sigma }_z{\tau }_0+c{h}_x{\sigma }_x{\tau }_0,$$

$$U^{\dagger}\big{[}{\cal H}_{\rm step}(k)+{\cal H}^{Z}_{\rm step}\big{]}U=\begin{pmatrix}0&D(k)\\ D^{\dagger}(k)&0\end{pmatrix},\ U=\begin{pmatrix}\tau_{z}&-\tau_{z}\\ \tau_{x}&\tau_{x}\end{pmatrix},\ D(k)=\begin{pmatrix}(-i\tilde{\Lambda}-v-\tilde{v})k-h_{z}&\Delta+ch_{x}-i\Lambda\\ \Delta-ch_{x}-i\Lambda&(-i\tilde{\Lambda}-v+\tilde{v})k-h_{z}\end{pmatrix}.$$

$$\mathrm{Im}\{\det [D(k\to 0)]\}>0\mathrm{and}\mathrm{Im}\{\det [\mathrm{D}(\mathrm{k}\to \pm \infty )]\}<0$$

Full text

Sn1-xPbxTe1-ySey  systems have attracted interest due to realization of 3D topological crystalline insulator phase Fu (2011); Hsieh et al. (2012); Dziawa et al. (2012); Tanaka et al. (2012); Xu et al. (2012), prediction of 2D topological phases Liu and Fu (2015); Safaei et al. (2015); Liu et al. (2014); Bassanezi et al. (2018); Araújo et al. (2018) and appearance of low-energy states at the defects Tang and Fu (2014); Sessi et al. (2016). Robust 1D modes were observed at the surface steps separating regions of even and odd number of layers Sessi et al. (2016) and interpreted as topological flat bands using a model obeying a chiral symmetry Rechciński and Buczko . In a more accurate description step modes can have a band width but the local density of states (LDOS) is large so that the system is susceptible to formation of correlated states Tang and Fu (2014); Volovik (2018); Ojajärvi et al. (2018). Recent experiments indicate that an order parameter emerges at low temperatures and it is accompanied with an appearance of a robust zero-bias peak (ZBCP) in the tunneling conductance Mazur et al. ; Das et al. (2016). The temperature and magnetic field dependence of the energy gap are consistent with superconductivity and under such circumstances the ZBCP is often interpreted as an indication of Majorana zero modes Mazur et al. ; Das et al. (2016); Wang et al. (2016); Aggarwal et al. (2016), which are intensively searched non-Abelian quasiparticles Lutchyn et al. (2018). Thus, this finding calls for a critical study of different mechanisms which may explain the appearance of the ZBCP.

We show that Sn1-xPbxTe1-ySey multilayers are a paradigmatic system for realization of topological phases due to emergent symmetries of the low-energy theory, and ZBCP can appear in the absence of superconductivity. The important 2D topological invariants are the mirror-resolved Chern number $C_{\pm}$ (due to structural mirror symmetry) and spin-resolved Chern number $C_{\uparrow(\downarrow)}$ (due to approximate spin-rotation symmetry). For odd number of layers $N$ the mirror symmetry is a point-group operation whereas for even $N$ it is a non-symmorphic (NS) symmetry (Fig. 1), so that adding one layer can change the topology of the system Bassanezi et al. (2018); Araújo et al. (2018). We calculate the dependence of the Chern numbers on $N$ and show that the step states appear when these invariants are different on the two sides of the step. The theory and experiment Sessi et al. (2016); Rechciński and Buczko attribute step states only to odd-height steps, but we predict that also even-height steps can exhibit step-states consistent with other experiment Iaia et al. .

We discuss the conditions under which a spontaneous symmetry breaking gives rise to an energy gap at the step and study an easy-axis magnetic order as one possibility. We show that magnetic domain walls (DWs) support low-energy bound states because the regions with opposite magnetization are topologically distinct in the presence of NS chiral and mirror symmetries. Due to the appearance of DWs the Fermi level is pinned to the energy of the DW states for a range of electron density providing an explanation for the ZBCP observed in the experiment Mazur et al. . The observed temperature and magnetic field dependencies are consistent with our theory.

Our starting point is a $p$-orbital tight-binding Hamiltonian describing bulk topological crystalline insulator in Sn1-xPbxTe1-ySey-material class Hsieh et al. (2012)

[TABLE]

where we have chosen a cubic unit cell with internal sites at the corners labeled by $i=1,\dots,8$ sup , $\varepsilon_{\alpha\beta}$ is a Levi-Civita symbol, $L_{\alpha}=-i\varepsilon_{\alpha\beta\gamma}$ are the $3\times 3$ angular momentum $L=1$ matrices, $\sigma_{\alpha}$ are Pauli matrices acting in the spin space, $\Sigma$ is a diagonal $8\times 8$ matrix with entries $s_{i}=\pm 1$ at the two sublattices [(Sn,Pb)/(Te,Se) atoms], and $h_{\alpha}^{(1)}(k_{\alpha})$ and $h_{\alpha,\beta}^{(2)}(k_{\alpha},k_{\beta})$ are $8\times 8$ matrices describing hopping between the nearest-neighbor and next-nearest-neighbor lattice sites in the directions $\hat{\alpha}$ and $\hat{\alpha}+\epsilon_{\alpha\beta}\hat{\beta}$, respectively sup . We allow the possibility to tune the spin-orbit coupling terms $\lambda_{\alpha}$ ($\alpha=x,y,z$) to be different from each other although in the real material $\lambda_{\alpha}=\lambda$. When not otherwise stated we use $m=1.65$ eV, $t_{12}=0.9$ eV, $t_{11}=0.5$ eV and $\lambda=0.3$ eV foo (a).

We focus on Hamiltonian ${\cal H}_{N}(\vec{k})$ for $N$ layers in the $z$-direction. The mirror symmetries for odd/even N can be written as $M_{z}^{o/e}(k_{x})=\sigma_{z}\otimes(2L_{z}^{2}-1)\otimes m_{z}^{o/e}(k_{x})$, where $m_{z}^{o}$ is a point-group reflection and $m_{z}^{e}(k_{x})$ is a (momentum-dependent) NS operation consisting of reflection and a shift by a half lattice vector [Figs. 1(a,b)]. To calculate $C_{\pm}$ we split ${\cal H}_{N}(\vec{k})$ into two blocks in the $M_{z}^{o/e}(k_{x})$ eigenspace. Since $M_{z}^{o/e}(k_{x})$ anticommutes with time-reversal symmetry (TRS) operator ${\cal T}={\cal K}\sigma_{y}\otimes\mathbbm{1}_{3}\otimes\mathbbm{1}_{8}$, these blocks carry opposite Chern numbers $C_{\pm}$ Chiu et al. (2016). We find that $C_{+}$ oscillates between $+2$ and $-2$ for $N=2n+1$ (we exclude $N=1$) and $C_{\pm}=0$ for $N=2n$ ($n\in\mathbb{N}$) sup ; foo (b). We point out that for even number of layers the non-symmorphic nature of the mirror symmetry guarantees that the Chern number calculated for the blocks must always be equal to zero. Therefore, one can equivalently conclude that the mirror-resolved Chern number does not exist as a topological invariant for even number of layers. The edge state spectra for $N=3$ and $N=4$ are shown in Figs. 1(c,d). As predicted by $C_{\pm}=\mp 2$ we see two pairs of gapless edge modes in the case $N=3$, but surprisingly we find four pairs of edge modes if $N=4$. These edge modes are consistent with $C_{\pm}=0$ because there exists small gaps which do not vanish by increasing system size.

To understand the existence of edge modes in the case of $N=4$ we notice that the momentum in the $z$-direction is quantized and the low-energy degrees of freedom are associated with a motion within the $(x,y)$-plane foo (c). Therefore, the components of the spin-orbit coupling $\lambda_{\alpha}$ ($\alpha=x,y,z$) contribute differently to the spectrum, and the dominant effect comes from $\lambda_{z}\sigma_{z}L_{z}$. Hence, turning off $\lambda_{x}$ and $\lambda_{y}$ is a good approximation [cf. Figs. 1(d) and 2(b)] and this leads to a spin rotation symmetry with respect to the $z$-axis. Thus, by block-diagonalizing ${\cal H}_{N}(\vec{k})$ we can calculate $C_{\uparrow}=-C_{\downarrow}$ as a function of $N$ sup . The results are shown in Fig. 2(a) and suggest that $C_{\uparrow}$ grows linearly with $N$ and takes values $C_{\uparrow}=2+4m$ ($m\in\mathbb{Z}$) for odd $N$ and $C_{\uparrow}=4m$ ($m\in\mathbb{Z}$) for even $N$. The careful study whether this behavior persist for arbitrary thickness goes beyond the scope of the current manuscript, but importantly this numerical evidence is already enough that we obtain a topological description of step modes for reasonably large systems, including all the system considered in this manuscript. In particular, the result for $N=4$ is consistent with number of edge modes in Figs. 1(d), 2(b). For $\lambda_{x}=\lambda_{y}=0$ the tiny gaps originally present in the spectrum vanish completely.

These Chern numbers provide an interpretation for the appearance of the step modes because they appear whenever $C_{\uparrow}$ is different on the two sides of the step [Figs. 1(e)-(h)] sup . For a step separating even $N$ and odd $N$ $\Delta C_{\uparrow}=2+4m$ ($m\in\mathbb{Z}$) and therefore at least two pairs of helical step modes foo (d) exist at these steps in agreement with Refs. Sessi et al. (2016); Rechciński and Buczko [Fig 1(e),(g)]. These step modes are weakly gapped because the spin-rotation symmetry is only present as an approximate symmetry. However, we can use $C_{+}$ to show that these gaps vanish in the limit $N\to\infty$. Namely, by using a symmetrized step-construction shown in Fig. 1(f),(h) and fixing $N$ so that $C_{+}=\pm 2$ on the different sides of the step, we find that in a system containing steps on both surfaces there exists $|\Delta C_{+}|=4$ pairs of gapless edge modes protected by the mirror symmetry. In the limit $N\to\infty$ the step modes at the different surfaces are completely decoupled, and therefore each step supports two pairs of gapless edge modes. Interestingly, we find that (depending on $N$) also even-height steps can exhibit $|\Delta C_{+}|=4$ or $|\Delta C_{\uparrow}|=4|m|$ ($m\in\mathbb{Z}$) pairs of step modes consistent with the experiment Iaia et al. .

In the case of a one atom-high step and finite $N$ the step modes are weakly gapped even though $C_{\pm}$ are different on two sides of the step because the step breaks the mirror symmetry and hybridizes the mirror blocks. However, the step modes can still be exactly gapless for certain widths $N_{y}$ in $y$-direction if $\lambda_{x}=0$ [Fig. 2(c)]. This effect comes from the existence of an effective particle-hole symmetry which together with a mirror symmetry gives rise to an antiunitary chiral symmetry ${\cal S}={\cal K}\sigma_{y}\otimes(2L_{x}^{2}-1)\otimes(i\Sigma m_{x})$, where $m_{x}$ is a mirror reflection with respect to the $x$-plane interchanging the sublattices sup . Due to this symmetry the Hamiltonian ${\cal H}_{N,N_{y}}(k_{x})$ supports a $\mathbb{Z}_{2}$ Pfaffian-invariant, which protects a 1D Weyl point if it changes sign as a function of $k_{x}$ Zhao et al. (2016); Brzezicki and Cuoco (2017). As demonstrated in Fig. 2(c) this kind Weyl points can be realized at the steps if $\lambda_{x}=0$.

Since we have now established the topological origin of the step modes in the non-interacting system, we turn our attention to the correlation effects (e.g. spin, charge, orbital or superconducting order), which are inevitably present due to the large LDOS in the limit $N\gg 1$ Tang and Fu (2014); Volovik (2018); Ojajärvi et al. (2018); sup . Our aim is to show that there exists a mechanism for the appearance of the ZBCP in the absence of superconductivity (without analyzing the competition between different types of order sup ). We require that the order parameter opens an energy gap, which means that it breaks the symmetries associated with $C_{\uparrow(\downarrow)}$ and $C_{\pm}$. Thus, we assume that there exists a magnetic instability in the vicinity of the steps (due to magnetic impurities or electron-electron interactions sup ) giving rise to a Zeeman field $H_{Z}=\mathbf{h}\cdot\vec{\sigma}$. Because the step modes are approximately spin-polarized along the $z$-direction the directions of $\mathbf{h}$ within the $(x,y)$-plane are efficient in opening an energy gap, and due to spin-orbit coupling the gap depends on the direction of $\mathbf{h}$ within the $(x,y)$-plane sup . Therefore, the system realizes an easy-axis ferromagnet and the topological defects are DWs foo (e). In the following we consider $\mathbf{h}=(h_{x},0,0)$.

To study the topological DW states we determine the low-energy theory for a single step, the emergent symmetries and the topological invariants. Although $h_{x}$ breaks TRS the spectrum of a system with two steps still exhibits Kramers degeneracy at $k_{x}=\pi$ [Fig. 1(g)] due to a remaining NS TRS ${\cal T}^{\prime}(k_{x})={\cal K}\mathbbm{1}_{2}\otimes(2L_{y}^{2}-1)\otimes g(k_{x})r_{z}$, where $g(k_{x})$ is a diagonal matrix with entries $e^{\pm ikx/2}$ and $r_{z}$ denotes $\pi$-rotation with respect to $z$ axis sup . ${\cal T}^{\prime}(k_{x})$ squares to $+1$ ($-1$) at $k_{x}=0$ ($k_{x}=\pi$) which yields Kramer degeneracy only at $k_{x}=\pi$. We assign half of the states at $k_{x}=\pi$ to each step by selecting from each Kramer’s doublet the state with larger projection on each step. Our low-energy theory is obtained by expanding the Hamiltonian around $k_{x}=\pi$ in one of the subspaces of the projected states sup . By assuming $\lambda_{y}=0$ the system supports a $k$-dependent mirror symmetry $M_{x}(k_{x})=\sigma_{x}\otimes(2L_{x}^{2}-1)\otimes g(k_{x})$ and NS chiral symmetry $S(k_{x})=i\sigma_{y}\otimes\mathbbm{1}_{3}\otimes\Sigma\,m_{x}g(k_{x})$ [Fig. 3(a)] sup .

We find that in the presence of these symmetries there exists three topologically distinct phases shown in Fig. 3(b). The trivial phase for $|h_{x}|<h_{c}$ is separated from the two non-trivial ones by the energy gap closings at $h_{x}=\pm h_{c}$. The non-trivial phases are characterized by a NS chiral $\mathbb{Z}_{2}$ invariant $\nu=1$ sup ; Shiozaki et al. (2015), and the phases at $h_{x}<-h_{c}$ and $h_{x}>h_{c}$ are topologically distinct because the band-inversions occur in the different mirror sectors $M_{x}(\pi)=\pm 1$ [Fig. 3(b)]. It is not a priori known whether the DWs between the topologically distinct phases in this symmetry class support DW states Shiozaki et al. (2015). However, our calculations show that sharp interfaces between trivial and non-trivial phases (two non-trivial phases) support one (two) low-energy bound state(s) per DW [Figs. 3(c-f)]. The number of low-energy states in each case is consistent with the number of zero-energy DW states expected in the case of smooth DWs with slowly varying spin textures sup . Although these states resemble the zero-energy DW states considered in the context of Dirac equation Jackiw and Rebbi (1976) and Su-Schrieffer-Heeger (SSH) model Su et al. (1979, 1980); Heeger et al. (1988), there is an important difference because they are realized in a model belonging to a different symmetry class. Namely, the appearance of the DW breaks the symmetries, and therefore the energies of these states in the case of sharp DWs remain non-zero even if the DWs are well-separated sup . Moreover, the energies depend on the tight-binding parameters and in this sense they resemble the topological DW states in systems with more complicated unit cells consisting of three or more atoms Nurul Huda et al. (2018). Nevertheless, for realistic system parameters the DW states appear close to the zero energy sup .

Because the DW states have non-zero energy, in high-resolution tunneling spectroscopy one would observe two splitted peaks in the conductance. However, already small broadening of the energy levels leads to a single ZBCP [Fig. 3(g),(h)] sup . Moreover, for sufficiently large density of DWs the bound states hybridize and form a band inside the energy gap leading to a single ZBCP where the height of the peak depends on the density of DWs [Fig. 3(e),(h)]. The ZBCP is robust against variations of the density because in analogy to the SSH model Heeger et al. (1988) we expect that the DWs are the lowest energy charged excitations in the system, and therefore small density of excess electrons (excess holes) is accommodated in the system by increasing the number of DWs, so that up to a critical variation of the density the Fermi level is pinned to the energy of the DW states. Similar situation occurs in quantum Hall ferromagnets where the lowest energy charged excitations are skyrmions Girvin (1999) which appear due to excess electrons and have been observed experimentally Barrett et al. (1995). Also the parametric dependencies of the ZBCP and the energy gap are consistent with the experiment Mazur et al. . The increase of temperature suppresses the order parameter and the energy gap. The external magnetic field breaks the degeneracy of states with opposite magnetization leading to a confinement between the DWs similarly as a symmetry breaking term in the SSH model Heeger et al. (1988), so that the number of DWs and the magnitude of the ZBCP decrease. Furthermore, we find that by increasing the Zeeman field the energy gap of the system decreases sup . Therefore, all the observations can be explained without requiring the existence of superconductivity. Finally, the observation that magnetic dopants enhance the ZBCP and the energy gap Mazur et al. makes it more plausible that the effect originates from magnetic instability instead of superconductivity. The systematic analysis of the correlated states which are consistent with the observations Mazur et al. and the exploration of the possible common origin of the zero-bias anomalies in various topological semiconductors and semimetals Mazur et al. ; Das et al. (2016); Wang et al. (2016); Aggarwal et al. (2016) are interesting directions for future research sup .

Supplementary material for ”Topological properties of multilayers and surface steps in the SnTe material class”.

I Unit cell, Hamiltonian and symmetries

Our starting point is a $p$-orbital tight-binding Hamiltonian describing bulk topological crystalline insulator in Sn1-xPbxTe1-ySey-material class Hsieh et al. (2012)

[TABLE]

where $\hat{c}_{j\alpha}(\mathbf{r})$ are vectors of fermionic operators corresponding to $p_{x}$-, $p_{y}$- and $p_{z}$-orbitals and the indices denote the sublattice $j\in\{1,2\}$ [(Sn,Pb)/(Te,Se) atoms], spin $\alpha$ and lattice site $\mathbf{r}$. Here $\hat{\sigma}_{\alpha,\beta}$ is a vector of Pauli matrices, $\hat{d}_{{\bf r}{\bf r^{\prime}}}$ are unit vectors pointing from ${\bf r}$ to ${\bf r^{\prime}}$ and the next-nearest-neighbor hoppings satisfy $t_{11}=-t_{22}$. As discussed in the main text the 3D bulk Hamiltonian can be represented as

[TABLE]

where $\varepsilon_{\alpha\beta}$ is a Levi-Civita symbol, $L_{\alpha}=-i\varepsilon_{\alpha\beta\gamma}$ are the $3\times 3$ angular momentum $L=1$ matrices and we allow the possibility to tune the spin-orbit coupling terms $\lambda_{\alpha}$ ($\alpha=x,y,z$) to be different from each other although in the real material $\lambda_{\alpha}=\lambda$. The matrices $h_{\alpha}^{(1)}(k_{\alpha})$ and $h_{\alpha,\beta}^{(2)}(k_{\alpha},k_{\beta})$ describe hopping between sites of 8-site unit cells (see Fig. 1) being either nearest or next-nearest neighbors. These are given by,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

$\Sigma$ is a diagonal matrix describing mass modulation in a unit cell and has diagonal entries $(-1,1,-1,1,-1,1,-1,1)$. For $\lambda_{\alpha}=0$ the system, which is translationally invariant in at least $x$ direction, has a NS chiral symmetry of the form $S_{\alpha}(k_{x})=i\sigma_{\alpha}\otimes\mathbbm{1}_{3}\otimes\Sigma\,m_{x}g(k_{x})$, where $m_{x}$ is a mirror-$x$ reflection of the unit cell given by

[TABLE]

and $g(k_{x})$ is a diagonal gauge matrix with the entries $(e^{i\frac{k_{x}}{2}},e^{-i\frac{k_{x}}{2}},e^{-i\frac{k_{x}}{2}},e^{i\frac{k_{x}}{2}},e^{i\frac{k_{x}}{2}},e^{i\frac{k_{x}}{2}},e^{-i\frac{k_{x}}{2}},e^{-i\frac{k_{x}}{2}})$. The plane of the mirror $m_{x}$ cuts the unit cell into halves and interchanges the sublattices and this is why it appears in the chiral symmetry operator. Its relation to the Hamiltonian is $\{{\cal H}({\mathbf{k}}),S_{\alpha}(k_{x})\}=0$ for every $k$-point. In the main text we always use $S(k_{x})\equiv S_{y}(k_{x})$ so the $y$ subscript is omitted for brevity. We also point out that the sign of $\lambda$ is not important for our results because we can construct a chiral symmetry $S_{0}(k_{x})$ by using the identity matrix $\sigma_{0}$ in the expression given above, such that the Hamiltonians with positive and negative spin-orbit coupling are related to each other as ${\cal H}_{\lambda}(\mathbf{k})S_{0}(k_{x})=-S_{0}(k_{x}){\cal H}_{-\lambda}(\mathbf{k})$. The same gauge $g(k_{x})$ matrix is used for constructing the $x$-plane mirror reflection in the $k$-space $M_{x}(k_{x})=\sigma_{x}\otimes(2L_{x}^{2}-1)\otimes g(k_{x})$. It satisfies the usual relation with the Hamiltonian $M_{x}(k_{x}){\cal H}(k_{x},k_{y},k_{z})M_{x}^{\dagger}(k_{x})={\cal H}(-k_{x},k_{y},k_{z})$. Analogically the $z$-plane mirror reflection $M_{z}(k_{x})$ can be defined as $M_{z}(k_{x})=\sigma_{z}\otimes(2L_{z}^{2}-1)\otimes g(-k_{x})m_{z}$, where $m_{z}$ transforms the sites in the unit cell with a matrix

[TABLE]

and we get $M_{z}(k_{x}){\cal H}(k_{x},k_{y},k_{z})M_{z}^{\dagger}(k_{x})={\cal H}(k_{x},k_{y},-k_{z})$. Finally, there exists also a symmorphic $\pi$-rotation symmetry $R_{z}$ in the $xy$ plane $R_{z}{\cal H}(k_{x},k_{y},k_{z})R_{z}^{\dagger}={\cal H}(-k_{x},-k_{y},k_{z})$, where $R_{z}=\sigma_{z}\otimes(2L_{z}^{2}-1)\otimes r_{z}$ and

[TABLE]

By combining $R_{z}$ with $M_{x}(k_{x})$ we can obtain a $y$-plane mirror reflection $M_{y}(k_{x})=R_{z}M_{x}(k_{x})$ satisfying $M_{y}(k_{x}){\cal H}(k_{x},k_{y},k_{z})M_{y}(k_{x})^{\dagger}={\cal H}(k_{x},-k_{y},k_{z})$, and by combining $R_{z}$ with $M_{z}(k_{x})$ we get an inversion symmetry $I(k_{x})=R_{z}M_{z}(k_{x})$ satisfying $I(k_{x}){\cal H}({\mathbf{k}})I^{\dagger}(k_{x})={\cal H}(-{\mathbf{k}})$. The inversion symmetry can be combined with time-reversal symmetry ${\cal T}={\cal K}\sigma_{y}\otimes\mathbbm{1}_{3}\otimes\mathbbm{1}_{8}$ where ${\cal K}$ is complex conjugation. This way we obtain an operator ${\cal Q}(k_{x})={\cal T}I(k_{x})$ satisfying ${\cal Q}(k_{x}){\cal H}({\mathbf{k}}){\cal Q}^{-1}(k_{x})={\cal H}({\mathbf{k}})$. Since ${\cal Q}^{2}(k_{x})=-1$ for all $k_{x}$, every band is Kramers degenerate within the whole Brillouin zone (BZ).

When the system is not fully 3D but consists of a finite number of atomic layers $N$ stacked in the $z$ direction the spatial part of the symmetry operators may become richer if they permute the layers. When $N$ is even the supercell can be constructed by stacking cubic cells vertically but in case of odd $N$ we have to stack $(N-1)/2$ of whole cubic cells and one half extra. This means that for even $N$ the $z$-plane mirror reflection is a NS symmetry $M_{z}^{e}(k_{x})$, where the intracell reflection $m_{z}$ is otherwise the same as in the case of bulk $M_{z}(k_{x})$ but it also reverts the order of the unit cells in $z$ direction. On the other hand, for odd $N$ the reflection plane coincides with the central layer of the system, and this leads to a symmorphic $M_{z}^{o}$ version of the $M_{z}(k_{x})$ operator (with no $k_{x}$ dependence).

In both even and odd case $M_{z}^{o,e}(k_{x})$ anticommutes with ${\cal T}$, i.e. ${\cal T}M_{z}^{o,e}(k_{x})+M_{z}^{o,e}(-k_{x}){\cal T}=0$, and commutes with the multilayer 2D Hamiltonian $[{\cal H}_{N}(k_{x},k_{y}),M_{z}^{o,e}(k_{x})]=0$. This implies that the double degeneracy of every band is resolved in the $M_{z}^{o,e}(k_{x})$ invariant subspaces. It is due to the peculiarity of the NS symmetry that in the case of even $N$ the $M_{z}^{e}(k_{x})$ subspaces keep internal TRS. To understand this Brzezicki and Cuoco (2017) we can first order the eigenvectors of $M_{z}^{e}(k_{x})$ based on the eigenvalues and write them in the columns of matrix $U(k_{x})$. Then, in the basis given by $U(k_{x})$ the mirror symmetry $M_{z}^{e}(k_{x})$ is diagonal, and we denote the Hamiltonian and TRS as

[TABLE]

The Hamiltonian commutes with $M_{z}^{e}(k_{x})$ so it has block-diagonal form but TRS is block-off-diagonal due to anticommutation with $M_{z}^{e}(k_{x})$. Apparently there is no TRS inside the diagonal blocks of $\tilde{{\cal H}}_{N}$. One can however notice that $U(k_{x})$ is not $2\pi$- but $4\pi$-periodic due to nonsymmorphicity so one can define a non-trivial unitary matrix $\chi$ that produces a $2\pi$ shift of the momentum in the $U(k_{x})$ basis

[TABLE]

Now, knowing that eigensubspaces of $M_{z}^{e}(k_{x})$ are interchanged after a shift of $2\pi$ it is easy to notice that $\chi$ has a purely block-off-diagonal structure, just like $\widetilde{{\cal T}}$. Therefore, by combining these two operators we obtain an intra-block TRS, i.e.,

[TABLE]

which is a time-reversal with respect to the $(k_{x},k_{y})=(-\pi,0)$ point. This relation is enough to guarantee that the mirror-resolved Chern numbers vanish in the case of even number of layers. Another way of resolving the global double degeneracy is to have a spin-rotation symmetry around some axis which happens when only one component of SOC is present, i.e., having only $\lambda_{z}$ non-zero yields a symmetry given by $S_{z}=\sigma_{z}\otimes\mathbbm{1}_{3}\otimes\mathbbm{1}_{8}$. The spin subspaces do not depend on spatial symmetries so the spin-resolved Chern numbers are allowed for any $N$.

When a step along $x$ is introduced in the system, see Fig. 2(a), the dimension of the $k$-space is further reduced to 1D and we denote the Hamiltonian as ${\cal H}_{\{N,N^{\prime}\}}(k_{x})$ where $N$ and $N^{\prime}$ denote the number of layers at the either side of the step edge. Note that apart from the effective Hamiltonian obtained by projecting the Hamiltonian to a single step we always considered systems which are periodic in the $y$ direction, so that the steps come in pairs and the $y$-plane mirror $M_{y}(k_{x})$ is conserved. Then for $\lambda_{x}=0$ we can define an antiunitary chiral symmetry ${\cal S}={\cal T}M_{x}(k_{x})S_{x}(k_{x})$ which is symmorphic and satisfies ${\cal S}^{2}=+1$. It allows for transformation that makes the Hamiltonian purely imaginary. Its explicit form is ${\cal S}={\cal K}\sigma_{y}\otimes(2L_{x}^{2}-1)\otimes(i\Sigma m_{x})$ and the relation with the Hamiltonian is usual anticommutation although due to antiunitarity it does not yield two off-diagonal blocks of ${\cal H}_{\{N,N^{\prime}\}}(k_{x})$. A part of ${\cal S}$ symmetry is ${\cal T}S_{x}(k_{x})$ which is referred as an effective particle-hole symmetry in the main text. What is peculiar about ${\cal S}$ is that it is compatible with a Zeeman field term in any direction. Finally, if we assume $\mathbf{h}=(h_{x},0,0)$ and $\lambda_{y}=0$ the unitary symmetries that survive are mirror $M_{x}(k_{x})$ and chirality $S_{y}(k_{x})\equiv S(k_{x})$. The $y$-plane mirror $M_{y}(k_{x})$ is not a symmetry anymore because it anticommutes with the Zeeman field term but together with ${\cal T}$ it forms a NS time-reversal symmetry ${\cal T}^{\prime}(k_{x})={\cal T}M_{y}(k_{x})$ that can be expressed as ${\cal T}^{\prime}(k_{x})={\cal K}\mathbbm{1}_{2}\otimes(2L_{y}^{2}-1)\otimes g(k_{x})r_{z}$ and it satisfies ${\cal T}^{\prime}(k_{x}){\cal H}_{\{N,N^{\prime}\}}(k_{x}){\cal T}^{\prime-1}(k_{x})={\cal H}_{\{N,N^{\prime}\}}(-k_{x})$. It is another peculiarity of a NS symmetry that ${\cal T}^{\prime 2}(0)=1$ but ${\cal T}^{\prime 2}(\pi)=-1$ so that the Kramers degeneracy appears only at $k_{x}=\pi$ Brzezicki and Cuoco (2017).

II Calculation of the mirror- and spin-resolved Chern numbers

Chern number for a given band can be calculated from the definition involving Berry connection and Berry curvature. This can be quite non-trivial when the bands are degenerate Gradhand et al. (2012) but for the simple case of non-degenerate bands one defines the Berry connection for the $n$-th band as

[TABLE]

and the Berry curvature as

[TABLE]

The integral of the Berry curvature over the Brillouin zone summed over all occupied bands gives an integer called Chern number

[TABLE]

This definition of Chern number is however not convenient for numerical calculation because it requires a smooth gauge of vectors $\left|n,\vec{k}\right\rangle$ which for large matrices are known only numerically. This difficulty can be overcome by special gauge-invariant prescriptions for Berry curvature Fukui et al. (2005); Deng et al. (2012) but since they yield Berry curvature band by band they become singular whenever there is a band crossing. If we do not need to know the curvatures of individual bands, and this is the case here, the most practical way is to express the Chern number as

[TABLE]

where in the sum $n$ runs over the occupied bands and $n^{\prime}$ over unoccupied ones. The above equation is so-called Kubo formula that describes Hall conductivity in the quantum Hall systems Nagaosa et al. (2010). In our case we do not insert the full Hamiltonian into the Kubo formula because it would give zero Chern number due to TRS. Instead we insert one of the diagonal blocks of ${\cal H}$ in the eigenbasis of $M^{o,e}_{z}(k_{x})$ for the mirror-resolved Chern number and $S_{z}$ for the spin-resolved Chern number. The resulting mirror-resolved Chern numbers $C_{+}$ are given in the Table 1. We note that they keep on oscillating between values $\pm 2$ for $N>1$ as the number of layers grows. Assuming that this behavior persists for arbitrary $N$ we predict that the step states appear generically in the limit of large $N$ when every surface step can be well approximated by a symmetrized step [Fig. 1(f) of the main text].

III Approximate spin-rotation symmetry versus topology

In the main text we argue that taking $\lambda_{x}=\lambda_{y}=0$ is a good approximation for the low-energy theory of a system where the layers are stacked along the $z$-direction. Namely, because of the confinement in the $z$ direction the $p_{x}$- and $p_{y}$-orbitals are active and thus $\lambda_{z}$ is responsible for the dominating spin-orbit coupling term. One can however ask whether this is still a good approximation for a very thick multilayer that seems to approach the three-dimensional case. In such a case the non-trivial topology is guaranteed by the mirror Chern number in the diagonal mirror planes, for instance in the $(110)$ plane. If all components of the spin-orbit coupling are present this invariant for the parameters used in the main text takes values ${\cal C}_{\pm}=\pm 2$ and leads to four Dirac points around X point at the $(001)$ surface of the system. Now we can ask what happens if we assume conserved spin-rotation symmetry. Is it possible that this will not affect the topology of a 3D system? The answer is positive if we choose the axis of the spin-rotation symmetry to be perpendicular to the $(110)$ plane when we calculate the mirror Chern number i.e. we obtain the same mirror Chern number as in the case where all three components of $\lambda_{i}$ are included.

On the other hand, the bulk-boundary correspondence exists for surfaces that are perpendicular to the $(110)$ plane so that the electrons are necessarily confined to a plane which is not compatible with the spin-rotation axis used in the calculation of the mirror Chern number. Nevertheless, we find that all calculations lead to correct approximate results as long as the axis of the spin-rotation symmetry is consistently chosen to be perpendicular to the plane where the electrons are confined to be. In Fig. 3 we show that the surface Dirac cones are correctly reproduced in the low energy spectrum of a system with $N=200$ layers stacked in the $(001)$ direction for $\lambda_{z}=0.3$ eV and $\lambda_{x}=\lambda_{y}=0$.

IV Dependence of the step mode dispersion on the number of layers

As shown in Fig. 4, for increasing number of layers $N$ the system tends to develop Dirac cones with closing of the bulk gap that are consistent with the $(001)$ surface states of the TCI. The dispersion of the step modes connecting these cones becomes more and more flat with increasing $N$.

V Completeness of the topological description of the step modes

For large number of layers $N$ one expects that number of the step modes appearing at an atomic step should depend only on the height of the step rather than on $N$ itself. This is because adding or substracting a layer far away from the surface should not affect the step states. Having this in mind it is possible to argue that the step states cannot be entirely determined by the mirror- and spin-resolved Chern numbers of the multilayer systems. By looking the spin-resolved Chern numbers given in Fig. 2 of the main text we see that it is possible that two different even-layer systems have the same spin-resolved Chern numbers (e.g. $N=20$ and $N=22$) and the same mirror-resolved Chern numbers (always zero), and we expect that such situation can occur for arbitrarily large $N$. We can now make an interface of these systems so that there is a unit height step on top and bottom surfaces. Such kind of system must support step modes because using the arguments discussed in the main text we know that for large $N$ the unit height steps support step modes and the surfaces are only weakly coupled in this limit. Thus it remains an open question whether one can define additional topological invariant that would describe completely all possible interfaces.

VI Dependence of the energy gap on the Zeeman field

The behavior of the energy gap $E_{g}$ of a system with step states is non-trivial when the Zeeman field is switched on. In Fig. 5(a) we show the dependence of $E_{g}$ on $h_{x}$. Consistently with the results presented in the main text above some critical value of $h_{x}$ the gap opens up to maximal value around $h_{x}=0.06t_{12}$ but then it decreases and tends to close again. In our explanation the closing and reopening of the energy gap will not be seen experimentally because of the spontaneous magnetization present due to magnetic impurities/magnetic instability of step mode electrons. The dependence on the direction of the field in the $xy$ and $xz$ planes is shown in Figs. 5(b-c). We expect that in the case of magnetic instability there exists an easy-axis along $x$-direction because the magnetization in $x$-direction maximizes the gap and leads to a minimum of the free energy.

VII Dependence of the domain wall state energies on the system parameters

The scaling properties of the energy gap of the step modes $E_{g}$ and the energy of the DW states $E_{0}$ in case of sharp DWs are complicated and depend on the value of the Zeeman field. Here we consider the field in the $x$ direction. In Fig. 6 we show the values of $h_{x}=h_{x}^{\rm max}$ that maximize $E_{g}$ for systems of 2/3 and 6/7 layers as a function of $N_{y}$. We observe saturation of both $h_{x}^{\rm max}$ and $E_{g}$ in the limit $N_{y}\to\infty$. We now consider the case of optimal energy gap and study the energies of the DW states for magnetic domains with opposite magnetization $\pm h_{x}$ as a function of the system size. The results are presented in Fig. 7 for various system sizes and we consider $E_{0}/E_{g}$ as the figure of merit. In these figures $N_{y}$ is taken in a range of roughly $N_{x}<N_{y}<5/2N_{x}$. We see that $E_{0}/E_{g}$ is typically between $0.1$ and $0.2$ and there is only a weak dependence on the system size.

We also find that the energies $E_{g}$ and $E_{0}$ are surprisingly immune to the changes of the Zeeman field direction. In Fig. 8 we show the low energy spectrum as function of the field direction for a fixed magnitude of the Zeeman field.

VIII Low-energy theory for the step modes and domain walls

By solving the eigenstates corresponding to one of the steps at $k_{x}=\pi$ and by expanding the Hamiltonian for momentum $k_{x}=\pi+k$, the low-energy theory for the step modes in the absence of Zeeman field can be written as

[TABLE]

where $\Delta$ is the energy of the step modes at $k=0$ (in the absence of $\lambda_{x}$), $v$ and $\tilde{v}$ describe the velocities of the step modes ($|v|>|\tilde{v}|$), and $\Lambda$ and $\tilde{\Lambda}$ are small terms arising due to the spin-orbit coupling $\lambda_{x}$. We find that in the lowest order approximation $\lambda_{y}$ does not have any effect in the low-energy theory. If $\lambda_{x}=0$ the model obeys spin-rotation symmetry around $z$-axis

[TABLE]

and supports two pairs of gapless counterpropagating edge modes as described in the main text. Due to $\lambda_{x}\neq 0$ ($\Lambda\neq 0$, $\tilde{\Lambda}\neq 0$) the spin-rotation symmetry is weakly broken and the step modes are gapped.

Within this low-energy model the non-symmorphic symmetries become usual symmorphic symmetries, which are obtained by projecting them to the states at $k_{x}=\pi$. The non-symmorphic chiral symmetry defines a chiral symmetry

[TABLE]

and the $k$-dependent mirror symmetry leads to a relation

[TABLE]

Additionally the system obeys time-reversal symmetry

[TABLE]

This is the most general linearized Hamiltonian in the presence of these symmetries.

The Zeeman field component $h_{x}$ preserves both $M_{x}$ and $C$, whereas $h_{z}$ breaks $M_{x}$ but preserves $C$. Therefore, we consider the effect of both of these components on the topology of the system. In the low-energy theory around $k_{x}=\pi$ they give rise to a Hamiltonian

[TABLE]

where the factor $c$ takes into account the renormalization of the effect of $h_{x}$ within the low-energy theory and $h_{z}$ is practically unmodified due to the fact that the step modes are spin-polarized as a good approximation. Since both of these terms obey the chiral symmetry [Eq. (11)], we can utilize the chiral symmetry to block-off diagonalize the Hamiltonian

[TABLE]

It is easy to see ${\rm Im}\{\det[D(k)]\}$ is independent of $h_{x}$ and $h_{z}$ at $k=0$ and $k\to\pm\infty$. Moreover, the sign of ${\rm Im}\{\det[D(k)]\}$ is different at $k=0$ and $k\to\pm\infty$ [see Fig. 9], and in the following we consider

[TABLE]

as a constraint which is always required to be satisfied. Within the low-energy theory this requirement is robustly satisfied even if we introduce an arbitrary small perturbation to the Hamiltonian. Moreover, our considerations are valid beyond the low-energy theory because if we equate the $k\to\pm\infty$ limit with $k\to 2\pi$ this constraint is always satisfied in the full theory where all the higher order terms in the expansion of ${\cal H}_{\rm step}(k)$ are correctly taken into account. This follows from the fact that the chiral symmetry $C$ originates from non-symmorphic chiral symmetry, which guarantees that a relation ${\rm Im}\{\det[D(0)]\}=-{\rm Im}\{\det[D(2\pi)]\}$ should always be satisfied Shiozaki et al. (2015). Importantly, in the presence of the constraint (16) we can identify topologically distinct phases by plotting the trajectories of the $Z(k)=\det[D(k)]$ in the complex plane [see Fig. 9].

If we first set $h_{z}=0$ we reproduce the phase diagram discussed in the main text. Here the non-symmorphic chiral $\mathbb{Z}_{2}$ invariant $\nu$ can be determined by plotting the trajectory of $Z(k)$ from $k=0$ to $k=\infty$ in analogy to the trajectory plotted from $k=0$ to $k=2\pi$ in the case of lattice Hamiltonian obeying non-symmorphic chiral symmetry Shiozaki et al. (2015). In the case of $\nu=0$ the trajectory is on the left side of the origin and in the case $\nu=1$ on the right side of the origin [see Fig. 9]. Therefore, these two phases cannot be smoothly deformed into each other without going through the origin (i.e. closing the energy gap). Moreover, the trajectory of $Z(k)$ from $k=-\infty$ to $k=0$ reproduces the same trajectory so that the behavior of the Hamiltonian at the momenta $k\to\pm\infty$ can identified with the behavior of the lattice Hamiltonian at momentum $k=2\pi$, and therefore we expect that the full Hamiltonian where all the higher order terms in the expansion of ${\cal H}_{\rm step}(k)$ are correctly taken into account would also reproduce the same behavior. Therefore, although we have identified the $\nu=0$ and $\nu=1$ phases using the low-energy theory we expect that they are topologically distinct phases also in the case of the full lattice Hamiltonian (in analogy to the topologically distinct phases discussed in Ref. Shiozaki et al. (2015)). As discussed in the main text the band-inversions at $h_{x}=\pm h_{c}$ occur at different mirror sectors so that we can additionally label the non-trivial phases as $\nu=1_{\pm}$ as long as $h_{z}=0$.

The $z$-component of the Zeeman field $h_{z}$ breaks the mirror symmetry but preserves the chiral symmetry. Therefore for $h_{z}\neq 0$ the trajectories $k=0$ to $k=\infty$ and $k=-\infty$ to $k=0$ are no longer identical. However, as long as both of these trajectories remain on the same side of the origin, the Hamiltonians are topologically equivalent with the corresponding $\nu=0$ (trajectory on the left side) or $\nu=1$ (trajectory on the right side) Hamiltonians discussed above. Interestingly, in the presence of $h_{z}\neq 0$ there exists also a third topologically distinct phase within the low-energy theory, which we have labelled $\nu=1/2$. In this case the trajectories from $k=-\infty$ to $k=0$ and $k=0$ to $k=\infty$ are on the different sides of the origin (so that in a sense half of the trajectory is on the left and half on the right). This phase does not have lattice analogue Shiozaki et al. (2015) but it is still a topologically distrinct phase within the low-energy theory, because there necessarily exists a gap closing when deforming the Hamiltonian belonging to this topological class into a Hamiltonian belonging to either $\nu=0$ or $\nu=1$ topological classes in the presence of the constraint (16).

This topological phase diagram enables a full description of the properties of the magnetic DWs where the magnetization $(h_{x}(x),0,h_{z}(x))$ varies slowly as a function of $x$. For such kind of smooth magnetic DWs the number of gap closings (zero energy DW states) is determined by the number of topological phase transitions occurring along the trajectory $(h_{x}(x),0,h_{z}(x))$ as a function of $x$. For example if one considers a DW between $\nu=1_{+}$ and $\nu=1_{-}$ phases any smooth DW will lead to at least two gap closings. The low-energy theory does not allow us to calculate the DW states appearing in the case of sharp DWs where the magnetization $(h_{x}(x),0,h_{z}(x))$ varies quickly as a function of $x$. However, the numerical calculations described the main text and in the supplementary material are consistent with the idea that the number of low-energy bound states appearing in the case of sharp DWs is the same as the number of gap closing points in the case of smooth DWs.

IX Calculation of the tunneling conductance

The tunneling current $I$ due to applied voltage $V_{dc}$ is described by Bruus and Flensberg (2004)

[TABLE]

where $n_{F}$ is Fermi function and $T(E)$ the transmission probability given by

[TABLE]

Here $A_{L\nu}(E)$ and $A_{R\mu}(E)$ are the spectral functions of states $\nu$ and $\mu$ localized on the left $L$ and right $R$ side of the tunneling barrier, and $T_{\nu\mu}$ is the tunneling matrix element. In the case of soft point-contact spectroscopy we may assume that one of the systems has a constant density of states so that

[TABLE]

and thus in zero-temperature limit

[TABLE]

In realistic systems the spectral functions have a broadening $\Gamma$ due to electron-phonon and electron-electron interactions, disorder and coupling of the tip to the system. The simplest approximation in which these effects can be accounted for is Lorenzian form

[TABLE]

where $\epsilon_{\mu}$ is the energy of state $\mu$. This way we obtain

[TABLE]

with $G_{0}={eI_{0}}/{\pi\Gamma}$, and this formula has been used for plotting the tunneling conductances in the main text. Moreover, we define $G_{\rm peak}=G(0)$ and $\Gamma_{0}=\sqrt{3}E_{0}$, and in the main text we use $\Gamma=\Gamma_{0}$ in all figures. Because the DW states are not exactly at zero-energy one would observe a double-peak structure in the tunneling conductance in the case of small broadening. However, for broadening $\Gamma\geq\Gamma_{0}$ one observes a single zero-bias peak (see Fig. 10).

X Domain wall states in the presence of randomly arranged magnetic impurities

In the case when the magnetism at the step is due to the magnetic impurities one can ask whether the disorder in the arrangement of these impurities will affect the low-energy states at the domain walls. To show that these states are not affected by the disorder we implement a system with $10\%$ concentration of magnetic impurities distributed randomly in space with magnetic moments pointing along $x$. Then we create a domain wall by choosing opposite $h_{x}$ in two halves of the system with magnitude such that $0.1|h_{x}|>h_{c}$ to compensate the dilution of magnetic centers. The results for the energy spectrum, differential conductance and LDOS [Figs. 11(a-c)] are very similar to the ones obtained without disorder (Fig. 3 in the main text).

Recent point-contact experiments done on the SnTe class of materials show not only zero-bias conductance peaks at low temperature, but also a soft gap feature at higher temperatures Das et al. (2016); Mazur et al. . Within our theoretical framework we expect that the main effect of temperature will be melting of the magnetic order. To simulate this we assume that the randomly distributed magnetic impurities also have random direction of the magnetic moments in the (x,y)-plane. This leads to the energy spectrum and differential conductance shown in Fig. 12. As we can see the soft gap feature is correctly reproduced.

XI Competition between different types of correlated states at the steps and possible experimental signatures

We emphasize that within the framework of our theoretical investigations the correlated phase realized in the experiments Das et al. (2016); Mazur et al. does not necessarily need to be an easy-axis ferromagnetic. Rather we use the easy-axis ferromagnet as an example of non-superconducting state which is consistent with the parametric dependencies of the observed tunneling conductance – demonstrating that the experimental observations can be explained without requiring the existence of superconductivity. In this section we briefly discuss the possible options for the correlated states realized at the steps and the experimental methods for distinguishing them from each other.

The appearance of correlated states at the steps due to the approximately flat bands have analogies to several other systems supporting flat bands, such as symmetry-broken states (ferromagnets, exciton condensates and more complicated order) in quantum Hall systems Girvin (1999); Moon et al. (1995); Nomura and MacDonald (2006); Alicea and Fisher (2006); Yang et al. (2006); Kharitonov (2012a, b); Young et al. (2013), competition between flat-band surface superconductivity and magnetism in semimetals Kopnin et al. (2011); Volovik (2018); Ojajärvi et al. (2018); Pamuk et al. (2017); Löthman and Black-Schaffer (2017) and the correlated states in graphene Moiré superlattices Cao et al. (2018a, b); Sharpe et al. (2019); Lu et al. (2019); Liu et al. . The generic feature of these systems is that they are often characterized by a competition between different types of order parameters which are almost degenerate with each other Nomura and MacDonald (2006); Alicea and Fisher (2006); Yang et al. (2006) and the ground state depends on the details of microscopic model Kharitonov (2012a); Ojajärvi et al. (2018). Therefore, it is practically impossible to predict what type of correlated state is realized in the flat-band systems by using purely theoretical methods. Nevertheless, the combination of theoretical and experimental investigations has lead to impressive progress in our understanding of these systems Kharitonov (2012a, b); Young et al. (2013). There the idea is essentially to systematically analyze the properties of the possible candidate states and to propose experimental tests to find out which of them is realized. In the case of step states, by relying on the analogy to the previously studied systems we expect that both easy-axis ferromagnet (flat-band ferromagnetic state caused by the Coulomb interactions) and flat-band superconductivity (e.g. due to phonon-mediated attractive interactions) are realistic possibilities Girvin (1999); Moon et al. (1995); Ojajärvi et al. (2018). However, since there exists four approximately flat bands also more complicated order might appear in analogy to the SU(4)-symmetry-broken states studied in graphene in the quantum Hall regime Nomura and MacDonald (2006); Alicea and Fisher (2006); Yang et al. (2006); Kharitonov (2012a, b); Young et al. (2013).

Therefore we need to find out which of the candidate states are consistent with the experimental observations. As we discussed above the intrinsic properties of the non-superconducting symmetry-broken states at the step defects can lead to parametric dependencies of the tunneling conductance which can mimic all the features which are often attributed to the unconventional superconductivity and Majorana modes Das et al. (2016); Mazur et al. ; Wang et al. (2016); Aggarwal et al. (2016). This means that the current experimental evidence is not conclusive and new experiments are necessary. Four-probe STM and spin-polarized STM measurements can determine whether the correlated state at the steps is related to superconductivity or magnetic order. Moreover, our theory predicts that the topological nature of the step modes is independent of the thickness of the sample. Therefore, the correlated states can be studied also in thin films, where the dispersion of the step modes can be controlled with the thickness and the density with the help of gate voltages. This makes it possible to study the competition between different phases Ojajärvi et al. (2018); Löthman and Black-Schaffer (2017) in a controllable way.

There are several indications that our approach might capture the essential physics related to the experiments Das et al. (2016); Mazur et al. ; Wang et al. (2016); Aggarwal et al. (2016). In particular, it provides a possible explanation why the zero-bias anomalies appear in topologically nontrivial multilayer systems (topological step modes appear only if the systems on the different sides of the step are topologically distinct) and a microscopic mechanism for the appearance of symmetry-broken states in materials where no spontaneous symmetry-breaking happens in the bulk (steps support approximately flat bands with a large density of states). Moreover, our theory explains why these zero-bias anomalies are observed in materials with a sublattice degree of freedom (this is the origin of the non-symmorphic chiral symmetry). Since the explanation does not require the existence of superconductivity it is also consistent with the lack of direct signatures of superconductivity Mazur et al. and it provides a possible explanation for the observation that the increase of the concentration of magnetic dopants enhances the effect Mazur et al. (magnetic instability instead of superconducting one). From the viewpoint of magnetic instability the magnetic impurities give rise to another mechanism (in addition to the Coulomb interactions between electrons) favouring magnetic order and therefore if there exists competition between different phases they help to stabilize the magnetic order. On the other hand, magnetic impurities are practically always destructive for the superconducting order both in conventional and unconventional superconductors Balatsky et al. (2006); Alloul et al. (2009); Mackenzie et al. (1998), and also from the random matrix theory viewpoint the disorder in this symmetry class of superconductors would lead to subgap states at arbitrarily low energies and therefore to the destruction of the superconducting gap.