Inversion-protected higher order topological superconductivity in monolayer WTe$_2$

TL;DR

This paper predicts that monolayer WTe$_2$ can host higher-order topological superconductivity with Majorana corner modes, driven by its inversion symmetry and boundary properties, offering a new platform for topological quantum computation.

Contribution

The study introduces a microscopic model showing Majorana corner modes in monolayer WTe$_2$, revealing a new class of inversion-protected higher-order topological superconductors.

Findings

Majorana corner modes found in spin-triplet phases.

In-plane magnetic field stabilizes unconventional pairing with Majorana modes.

Proposes monolayer WTe$_2$ as a platform for higher-order topological superconductivity.

Abstract

Monolayer WTe$_2$, a centrosymmetric transition metal dichacogenide, has recently been established as a quantum spin Hall insulator and found superconducting upon gating. Here we study the pairing symmetry and topological nature of superconducting WTe$_2$ with a microscopic model at mean-field level. Surprisingly, we find that the spin-triplet phases in our phase diagram all host Majorana modes localized on two opposite corners. Even when the conventional pairing is favored, we find that an intermediate in-plane magnetic field exceeding the Pauli limit stabilizes an unconventional equal-spin pairing aligning with the field, which also hosts Majorana corner modes. Motivated by our findings, we obtain a recipe for two-dimensional superconductors featuring "higher-order topology" from the boundary perspective: Generally a superconducting inversion-symmetric quantum spin Hall material whose normal-state Fermi surface is away from high-symmetry points, such as gated monolayer WTe$_2$, hosts Majorana corner modes if the superconductivity is parity-odd. We further point out that this higher-order phase is an inversion-protected topological crystalline superconductor and study the bulk-boundary correspondence. Finally, we discuss possible experiments for probing the Majorana corner modes. Our findings suggest superconducting monolayer WTe$_2$ is a playground for higher-order topological superconductivity, and possibly the first material realization for inversion-protected Majorana corner modes without utilizing proximity effect.

Tables

Table 1. Table 1: The parities of the irreducible representations under the 1 ​ T ′ 1 superscript 𝑇 ′ 1T^{\prime} lattice symmetry operations. The action of the symmetries on crystal momentum and internal indices and the used Nambu basis are shown in the text.

	${\eta }_{C_{2x}}$	${\eta }_{M_x}$	Examples
$A_g$	+	+	${\hat{s}}_0\otimes {\hat{\sigma }}_0\otimes {\hat{l}}_0$
$B_g$	-	-	${\hat{s}}_0\otimes {\hat{\sigma }}_z\otimes {\hat{l}}_x$
$A_u$	+	-	$k_x{\hat{s}}_x\otimes {\hat{\sigma }}_0\otimes {\hat{l}}_z$
$B_u$	-	+	$k_x{\hat{s}}_z\otimes {\hat{\sigma }}_0\otimes {\hat{l}}_z$

Table 2. Table 2: The inversion eigenvalues p k ~ , n subscript 𝑝 ~ k 𝑛 p_{\tilde{\textbf{k}},n} of all the occupied BdG BdG \rm BdG bands n = 1 , ⋯ , 8 𝑛 1 ⋯ 8 n=1,\cdots,8 , ordered with increasing energy, at each of the TRIMs k ~ = ( 0 , 0 ) ~ k 0 0 \tilde{\textbf{k}}=(0,0) , ( π , 0 ) 𝜋 0 (\pi,0) , ( 0 , π ) 0 𝜋 (0,\pi) and ( π , π ) 𝜋 𝜋 (\pi,\pi) . The indicator κ 𝜅 \kappa for inversion-protected topological superconductors is given by the sum of these eigenvalues divided by 4.

$n$	1	2	3	4	5	6	7	8
$\Gamma$	-	-	+	+	-	-	+	+
$X$	-	-	+	+	-	-	+	+
$Y$	+	+	+	+	+	+	+	+
$M$	-	-	+	+	-	-	+	+

Equations

$$H_0(\text{k})={\hat{s}}_0\otimes \left({\hat{h}}_0(\text{k})-\mu \right)+V_{soc}{\hat{s}}_z{\hat{\sigma }}_z{\hat{l}}_y.$$

$$H_{\mathrm{int}}=\sum _{\text{r}{\text{r}}^{\prime }}$$

$$=\sum _{\text{r}}$$

$$\displaystyle H^{\rm BdG}_{\textbf{k}}=\left(\begin{array}[]{cc}H_{0}(\textbf{k})&\Delta(\textbf{k})\\ \Delta^{\dagger}(\textbf{k})&-T^{\dagger}H_{0}^{\dagger}(\textbf{k})T\end{array}\right)$$

$$\kappa =\frac{1}{4}\sum _{\text{k}\in \mathrm{TRIM}}\sum _n{\xi }_{\text{k}n},$$

$$\nu =\kappa \mathrm{mod}2.$$

$$\kappa =2{\kappa }_N,{\nu }_N={\kappa }_N.$$

$$\displaystyle h_{0}(\textbf{k})=\left(\begin{array}[]{cccc}\varepsilon_{d}(\textbf{k})&0&t^{AB}_{d}g_{\textbf{k}}&t^{AB}_{0}f_{\textbf{k}}\\ 0&\varepsilon_{p}(\textbf{k})&-t^{AB}_{0}f_{\textbf{k}}&t^{AB}_{p}g_{\textbf{k}}\\ t^{AB}_{d}g^{\ast}_{\textbf{k}}&-t^{AB}_{0}f^{\ast}_{\textbf{k}}&\varepsilon_{d}(\textbf{k})&0\\ t^{AB}_{0}f^{\ast}_{\textbf{k}}&t^{AB}_{p}g^{\ast}_{\textbf{k}}&0&\varepsilon_{p}(\textbf{k})\end{array}\right)$$

$$H_{\mathrm{int}}=\sum _{\text{r}}$$

$$+$$

$$\displaystyle\Gamma=\left(\begin{array}[]{cc}\mathbb{V}(0)&\mathbb{V}(2k_{F})\\ \mathbb{V}(-2k_{F})&\mathbb{V}(0)\end{array}\right).$$

$$(-1{)}^{\nu }=\prod _{\text{k}\in \mathrm{TRIM}}\prod _{m=1}^{n_{\text{k}}/2}{\xi }_{\text{k}2m},$$

$$\nu =\sum _{\text{k}\in \mathrm{TRIM}}\frac{n_{\text{k}}^{-}}{2},$$

$$\kappa =\frac{1}{4}\sum _{\text{k}\in \mathrm{TRIM}}\sum _{m=1}^n{\xi }_{\text{k}m}=\sum _{\text{k}\in \mathrm{TRIM}}\frac{n}{4}-\sum _{\text{k}\in \mathrm{TRIM}}\frac{n_{\text{k}}^{-}}{2}.$$

$$(-1{)}^{\nu }=(-1{)}^{\kappa }$$

$$\nu =\kappa \mathrm{mod}2,$$

$$(-1{)}^{{\nu }_N}=\prod _{\text{k}\in \mathrm{TRIM}}\prod _{m=1}^{n_{\text{k}}/2}{\xi }_{\text{k}2m},$$

$${\nu }_N=\sum _{\text{k}\in \mathrm{TRIM}}\frac{n_{\text{k}}^{-}}{2},$$

$${\kappa }_N=\frac{1}{4}\sum _{\text{k}\in \mathrm{TRIM}}\sum _{m=1}^{n_{\text{k}}}{\xi }_{\text{k}m}=\sum _{\text{k}\in \mathrm{TRIM}}\frac{n_{\text{k}}}{4}-\sum _{\text{k}\in \mathrm{TRIM}}\frac{n_{\text{k}}^{-}}{2}.$$

$$(-1{)}^{{\nu }_N}=(-1{)}^{{\kappa }_N}$$

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

$$+[m_0+m_1(\cos k_x+\cos k_y)]{\hat{\tau }}_z\otimes {\hat{s}}_0\otimes {\hat{\rho }}_z$$

$$+v\sin k_x{\hat{\tau }}_0\otimes {\hat{s}}_z\otimes {\hat{\rho }}_x+v\sin k_y{\hat{\tau }}_z\otimes {\hat{s}}_0\otimes {\hat{\rho }}_y$$

$$+\Delta \sin k_x{\hat{\tau }}_x\otimes {\hat{s}}_z\otimes {\hat{\rho }}_0+\Delta \sin k_y{\hat{\tau }}_y\otimes {\hat{s}}_0\otimes {\hat{\rho }}_z$$

$$\kappa =4(0):-(m_0-{ϵ}_0)<2m_1<m_0-{ϵ}_0$$

$$\kappa =3:-(m_0+{ϵ}_0)<2m_1<-(m_0-{ϵ}_0)$$

$$\kappa =2:2m_1<-(m_0+{ϵ}_0),$$

$${\mathcal{H}}_{\Gamma }^{\uparrow /\downarrow }(r,\theta )$$

$$\displaystyle=\left(\begin{array}[]{cccc}m_{+}(r)&0&\pm\Delta e^{\mp i\theta}k_{r}&0\\ 0&m_{-}(r)&0&\pm\Delta e^{\pm i\theta}k_{r}\\ \pm\Delta e^{\pm i\theta}k_{r}&0&-m_{+}(r)&0\\ 0&\pm\Delta e^{\mp i\theta}k_{r}&0&-m_{-}(r)\end{array}\right),$$

$$m_{+}(r)<0,m_{-}(r)>0\text{mbox}forr<R_a$$

$$m_{+}(r)>0,m_{-}(r)>0\text{mbox}for{R}_a<r<R_b$$

$$m_{+}(r)>0,m_{-}(r)<0\text{mbox}forr>R_b.$$

$$\displaystyle\psi^{\uparrow/\downarrow}_{23}(r,\theta)=e^{-\frac{1}{\Delta}|r-R_{a}|}e^{il\theta}\left(\begin{array}[]{c}e^{\mp i\frac{\theta}{2}}\\ 0\\ \pm ie^{\pm i\frac{\theta}{2}}\\ 0\end{array}\right),$$

Full text

Introduction— Extensive experimental and theoretical effort has been devoted to transition metal dichalcogenides (TMD), a family of materials with chemical formula $\rm{MX_{2}}$ (M $=$ transition metal, X $=$ S, Se, Te) known to host a rich variety of intriguing ground states, such as topological insulators and semimetalsQian et al. (2014); Tang et al. (2017); Fei et al. (2017); Wu et al. (2018); Li et al. (2017); Deng et al. (2016), charge density wavesCastro Neto (2001); Wilson et al. (1974, 2001); Sipos et al. (2008); Xi et al. (2015a); Ritschel et al. (2015); Li et al. (2015), and various types of possibly unconventional superconductivityYuan et al. (2014); Hsu et al. (2017); Ye et al. (2012); Lu et al. (2015); Shi et al. (2015); Xi et al. (2015b); Li et al. (2015); Sajadi et al. (2018); Fatemi et al. (2018). Moreover, tuning among these phases is possible by widely accessible experimental knobs, for example changing the thickness, pressureKusmartseva et al. (2009); Sipos et al. (2008); Kang et al. (2015); Pan et al. (2015); Chi et al. (2018), electrostatic gatingYe et al. (2012); Li et al. (2015); Yu et al. (2015), and recently even the twist angle between monolayersWu et al. (2018, 2018). Recently, a centrosymmetric member of the TMD family, monolayer $\rm WTe_{2}$, has been establishedQian et al. (2014); Fei et al. (2017); Tang et al. (2017); Jia et al. (2017); Peng et al. (2017); Wu et al. (2018) as a quantum spin Hall (QSH) insulatorKane and Mele (2005a); Bernevig and Zhang (2006). Remarkably, in this same material, superconductivity at temperatures around 1K was soon after reported under tunable electrostatic gatingSajadi et al. (2018); Fatemi et al. (2018). We are thus motivated to understand the nature of this superconductivity given the prevailing expectation that inducing superconductivity in already topological materials is a promising route for achieving topological superconductors.

Theoretically a known necessary condition for two-dimensional (2D) time-reversal topological superconductors requires negative pairing potentials on an odd number of Fermi surfaces that enclose time-reversal invariant momenta (TRIMs)Qi et al. (2010); Zhang et al. (2013). The presence of the inversion symmetry, however, enforces two-fold degeneracy of the Fermi surfaces and thus sets up a “no-go” theorem that precludes such superconductors from being topological. Nonetheless, recent developments suggest that inversion can unexpectedly enrich the topological structure of a systemKhalaf et al. (2018); Khalaf (2018); Ono et al. (2018), and enable new topological crystalline superconductors (TCsc) that are completely beyond the previous paradigmQi et al. (2010); Zhang et al. (2013). In particular, there exists a type of inversion-protected TCsc in dimension $d$ that has no Majorana boundary modes in $d-1$ dimension, yet is still topologically distinct from a trivial superconductorKhalaf (2018); Ono et al. (2018). This suggests the possibility that such inversion-protected TCsc belongs to the so-called “higher order topological phases”Benalcazar et al. (2017); Wang et al. (2018a); Yan et al. (2018); Wang et al. (2018b); Schindler et al. (2018); Langbehn et al. (2017); Shapourian et al. (2018); Wang et al. (2018b); Bultinck et al. (2018); Xu et al. (2019); Zhang et al. (2019, 2019), and may host Majorana boundary modes in $d-2$ or lower dimension.

Here, we propose a surprisingly simple recipe for this exotic inversion-protected TCsc: (1) the normal state is an inversion-symmetric QSH material with Fermi pockets away from TRIMs, and (2) the superconductivity is parity-odd. Given that gated monolayer $\rm WTe_{2}$ readily satisfies criteria (1), unconventional superconductivity with odd parity becomes the last piece of the puzzle for an inversion-protected TCsc that could host exotic Majorana corner modes.

In fact, in $\rm WTe_{2}$ there is ample reason to suspect that electron correlations might be strong, and odd-parity superconductivity is therefore plausible. First is the fact that the reported superconductivityFatemi et al. (2018); Sajadi et al. (2018) occurs at a low carrier density, while ab initio calculations do not reproduce the low-energy normal state band structure found by angle-resolved photoemission spectroscopy (ARPES)Tang et al. (2017) and scanning tunnelling microscopy (STM)Tang et al. (2017) studies unless one goes beyond the generalized-gradient approximation Qian et al. (2014); Zheng et al. (2016); Tang et al. (2017). Moreover, the reported in-plane upper critical field $H_{c2}^{\parallel}$ is 2.5-4.5 times higher than the Pauli limit $H_{p}$Fatemi et al. (2018); Sajadi et al. (2018). While an $H_{c2}^{\parallel}$ higher than the Bardeen-Cooper-Schrieffer theory prediction in centrosymmetric materials can occur when the normal state has a high spin-orbit scattering rateKlemm et al. (1975) or when the g factor deviates from twoSajadi et al. (2018), another possible origin is a spin-triplet (and thus odd-parity) paired state with spin aligning in the field direction.

In this work, we report the pairing symmetry and topological nature of the newly discovered superconductivity in gated monolayer WTe2. First, we solve the linearized gap equations to obtain a superconducting phase diagram in terms of microscopic interactions. By investigating the boundary modes in different phases, we find Majorana corner modes in odd-parity phases and surprisingly, also in the field-induced equal-spin phase emerging upon the suppression of conventional pairings. Then, we obtain a general recipe from the boundary perspective for achieving such 2D superconductivity with corner Majoranas. Finally, we point out that such higher-order phase is an inversion-protected topological crystalline superconductor that can be characterized by a bulk invariant we propose, and address the bulk-boundary correspondence. Our recipe provides a new route towards materializing a novel topological phase of matter, as well as realizing Majorana zero modes, which is the first step for topological quantum computation.

Model— Monolayer WTe2 is stable in the $1T^{\prime}$ structure, which is a buckled honeycomb lattice that is distorted into a rectangular lattice consisting of in-plane and buckled zigzag chains of W and Te atoms, respectively, see Fig. 1(a). This lattice is nonsymmorphic, with a two-fold screw rotation $C_{2x}$ and a glide mirror symmetry $M_{x}$111$M_{x}$ can be transformed back to the conventional mirror symmetry by shifting the mirror plane away from inversion center., each with a half-unit-cell translation along the chain direction $\hat{x}$. The lattice also has inversion symmetry $I_{0}$, resulting from the product of the two symmetries.

To study the dominant pairing channels in gated monolayer WTe2, we start from a minimal tight-binding model previously obtained by other authors from a low-energy fit to ab initio calculationsMuechler et al. (2016); Ok et al. (2018)222Ref. Lau et al., 2018 showed that more spin-orbit coupling terms are required for a better fit to the experimentally data. We nonetheless expect our results to change only quantitatively.. The Hamiltonian is written in a basis of spin s and four Wannier orbitals. These Wannier orbitals are labeled by the sublattice $\sigma=A,B$ they are on and by whether they transform as $d_{x^{2}-y^{2}}$ or $p_{x}$ orbitals ($l=d,p$). The $l=d,p$ orbitals are derived from W and Te atoms, respectively. Each degree of freedom is denoted by the corresponding Pauli matrices, $\hat{s}$, $\hat{\sigma}$, and $\hat{l}$, respectively. The full normal-state Hamiltonian is

[TABLE]

Here, the $s_{z}$-preserving intrinsic spin-orbit coupling $V_{soc}$Ok et al. (2018) is the lowest order term in k that obeys time-reversal, screw rotation, and glide mirror symmetries, while the spin-degenerate part $\hat{h}_{0}(\textbf{k})$Muechler et al. (2016); Ok et al. (2018) is a $4\times 4$ matrix in the basis of $\hat{\sigma}\otimes\hat{l}$ [Supplementary Material (SM) Sec. I]. As a zeroth-order approximation to the gating effects, we set the overall chemical potential $\mu=0.5$. The resulting two electron pockets are centered along the $\Gamma-X$ line [Fig. 2(a)], as observed by ARPESTang et al. (2017).

We consider short-ranged density-density interactions that preserve the lattice symmetries up to nearest-neighbor unit cells [Fig. 1(b)]:

[TABLE]

where $\sigma^{(^{\prime})}$, $l^{(^{\prime})}$, and $a=1,2$ indices are summed over, $n_{s\sigma l}(\textbf{r})$ is the density with spin $s$ and orbital $l$ locating at sublattice $\sigma$ in the unit cell centered at r, and $n_{\sigma l}(\textbf{r})=\sum_{s}n_{s\sigma l}(\textbf{r})$. Here $U^{l}$ denotes the on-site interactions for orbital $l$, $V_{1}^{ll^{\prime}}$ and $V_{2}^{ll^{\prime}}$ denote the nearest- and next nearest-neighbor interactions respectively on the zigzag chains with intra (inter)-orbital characters for $l^{\prime}=l(\bar{l})$ [Fig. 1(b)], and $\boldsymbol{\delta}_{a}$ denotes corresponding lattice vectors [SM Sec. I]. For simplicity, in the following we consider the case where $U^{l}=U$, and $V_{1}^{ll^{\prime}}=V_{2}^{ll^{\prime}}=V$.

Method and phase diagram— To analyze the dominant pairing channel for given interactions $U$ and $V$, we first classify the symmetries of possible pairing gaps. The normal state preserves two nonsymmorphic symmetries $C_{2x}=e^{ik_{x}a_{x}/2}(-i\hat{s}_{x}\otimes\hat{\sigma}_{x}\otimes\hat{l}_{0}$), $k_{y}\rightarrow-k_{y}$, and $M_{x}=e^{ik_{x}a_{x}/2}(-i\hat{s}_{x}\otimes\hat{\sigma}_{0}\otimes\hat{l}_{z}$), $k_{x}\rightarrow-k_{x}$. The mean-field Bogoliubov-de Gennes (BdG) Hamiltonian

[TABLE]

therefore obeys $g^{BdG}_{\textbf{k}}H^{\rm BdG}_{\textbf{k}}(g^{BdG}_{\textbf{k}})^{\dagger}=H^{\rm BdG}_{g\textbf{k}}$, where $T=is_{y}\mathcal{K}$, $\textbf{k}\rightarrow-\textbf{k}$ is the time-reversal operation with $\mathcal{K}$ the complex conjugation, and $g^{BdG}_{\textbf{k}}=\rm diag$ $[g_{\textbf{k}},\eta_{g}g_{\textbf{k}}]$ describes how the two symmetries $g=C_{2x},M_{x}$ act on the Nambu basis $[c_{\textbf{k}\uparrow},c_{\textbf{k}\downarrow},c^{\dagger}_{\textbf{k}\downarrow},-c^{\dagger}_{\textbf{k}\uparrow}]$. Thus, the pairing gaps transform as $g_{\textbf{k}}\Delta_{\textbf{k}}g^{\dagger}_{\textbf{k}}=\eta_{g}\Delta_{\textbf{k}}$, and we can classify all possible pairing gaps into four irreducible representations $A_{g}$, $B_{g}$, $A_{u}$, and $B_{u}$ according to their parities $\eta_{g}=\pm 1$ under the symmetry transformations $g$, see Table 1.

Next, we determine which irreducible representation has the highest $T_{c}$ by solving the linearized gap equationSigrist (2005) $\Delta_{\alpha^{\prime}\beta^{\prime}}(\textbf{k}^{\prime})=-\sum_{\textbf{k}^{\prime\prime}\textbf{k}}\Gamma_{\alpha^{\prime}\beta^{\prime},\beta^{\prime\prime}\alpha^{\prime\prime}}(\textbf{k}^{\prime},\textbf{k}^{\prime\prime})\times\chi_{\beta^{\prime\prime}\alpha^{\prime\prime},\alpha\beta}(\textbf{k}^{\prime\prime},\textbf{k},T)\Delta_{\alpha\beta}(\textbf{k})$, where Greek indices contain all the internal indices $(s,\sigma,l)$, and repeated indices are summed over. Here, the interaction $\Gamma_{\alpha^{\prime}\beta^{\prime},\beta\alpha}(\textbf{k}^{\prime},\textbf{k})$ is the Fourier transform of $\Gamma_{\alpha^{\prime}\beta^{\prime},\beta\alpha}(\textbf{r},\textbf{r}^{\prime})$ in Eq. I, and $\chi_{\beta^{\prime\prime}\alpha^{\prime\prime},\alpha\beta}(\textbf{k}^{\prime\prime},\textbf{k},T)$ is the non-interacting static pairing susceptibility at temperature $T$. Solving the lineared gap equation amounts to solving the eigenvalue problem of the effective interaction projected onto the Fermi surface $\tilde{\Gamma}(\textbf{p}^{\prime},\textbf{p})=-\sqrt{P}_{\textbf{p}^{\prime}}\Gamma(\textbf{p}^{\prime},\textbf{p})\sqrt{P}_{\textbf{p}}$, where $\textbf{p}^{(\prime)}$ is the incoming (outgoing) momentum on the Fermi surface, and $P_{\textbf{p}}=(\sum_{n=1,2}|\textbf{p},n\rangle\langle\textbf{p},n|)\otimes(\sum_{n=1,2}|-\textbf{p},n\rangle\langle-\textbf{p},n|)$ projects an electron-pair state to the two degenerate non-interacting bands $n$ on the Fermi surface at momenta p and $-\textbf{p}$. The eigenvector $\psi(\textbf{p})$ of $\tilde{V}$ with the most negative eigenvalue $\lambda$ is the solution to the linearized gap equation with the highest $T_{c}\propto\exp(-1/|\lambda|)$. We can then determine how $\psi(\textbf{p})$ behaves under symmetries $C_{2x}$ and $M_{x}$ under different interactions and obtain the superconducting phase diagram of $H=H_{0}+H_{\rm int}$.

In Fig. 2(b) we present this phase diagram as a function of $U$ and $V$. We find that while on-site attractions favor the even-parity “trivial” representation $A_{g}$ as expected, the odd-parity representations $A_{u}$ and $B_{u}^{(^{\prime})}$ (the superscript denotes different pair spin textures) dominate over a large portion of the phase diagram where the nearest-neighbor attraction $V$ dominates. In particular, the degenerate $A_{u}$ and $B_{u}^{\prime}$ gaps at repulsive $U$ are equal-spin triplet in the out-of-plane direction ($|\uparrow\uparrow\mp\downarrow\downarrow\rangle$), and the $B_{u}$ gap at attractive $U$ has $s_{z}=0$ ($|\uparrow\downarrow+\downarrow\uparrow\rangle$). This $SU(2)$-symmetry breaking is due to the intrinsic spin-orbit coupling $V_{soc}$.

We can understand qualitatively the competition between even- and odd-parity pairings from their real-space gap structures. To this end, we write down the mean-field Hamiltonian in Eq. 5 in real space and solve the self-consistency equations $\Delta_{\alpha^{\prime}\beta^{\prime}}(\textbf{r},\textbf{r}^{\prime})=-\sum_{\alpha\beta}\Gamma_{\alpha^{\prime}\beta^{\prime},\beta\alpha}(\textbf{r},\textbf{r}^{\prime})\langle c_{\textbf{r}^{\prime}\beta}c_{\textbf{r}\alpha}\rangle$ by iteration. We consider the short-ranged interactions $H_{\rm int}$, and show results for representative points for even- and odd-parity pairings [see blue stars in Fig. 2(b)]. We find the dominant component in the even-parity $A_{g}$ gap to be the on-site pairings as expected, while the dominant contribution to the odd-parity $B_{u}$ gap comes from the next-nearest-neighbor d-orbital pairing along the chains in the $\hat{\textbf{x}}$ direction [see the bonds with $\Delta_{1}$ in Fig. 2(c) and (d)]. It is then clear that attractions $U^{d}$ and $V_{2}^{dd}$ in Eq. I are the main terms responsible for $A_{g}$ and $B_{u}$ pairings respectively. While on-site attractions $U^{d}$ are uniform in momentum space and promote even-parity pairing, attractive $V_{2}^{dd}$ terms enhance scattering processes with large momentum-transfer $2k_{F}$ across the two pockets, which promotes odd-parity pairing [SM Sec. II]; hence the balance between even- and odd-parity pairings as shown in the phase diagram.

*Corner Majoranas in *WTe2— To understand the topological properties of these phases, we examine the boundary modes of different paired states in the phase diagram. While the spin-singlet $A_{g}$ phase is topologically trivial as expected, we find that spin-triplet phases exhibit exotic boundary modes. Specifically, in our model for superconducting WTe2 given by $H_{0}$ and the self-consistently obtained $B_{u}$ pairing 333For numerical convenience we take the self-consistent $B_{u}$ symmetry solution and multiply by 10, so that the resulting superconducting gaps are always much larger than any finite-size gaps of the normal bulk or edge states for tractable lattice sizes. , we numerically demonstrate the existence of zero-energy corner-localized states on an open-boundary geometry [Fig. 3]. We further verify that with an increasing system size $L$, these corner states tend exponentially toward zero energy [SM Sec. III], which unambiguously demonstrate the existence of Majorana Kramers pairs localized at two opposite corners. We also find similar Majorana corner modes the other spin-triplet phase $A_{u}$ [SM Sec. III].

Even if the realistic WTe2 lies in the even-parity pairing $A_{g}$ regime in Fig. 2(b), we find that an intermediate in-plane magnetic field can surprisingly drive a first-order phase transition and stabilize a new equal-spin phase $B_{u}^{\prime\prime}$444The superscript denotes a different spin orientation from that of $B_{u}^{(^{\prime})}$. aligning with the applied field near the Pauli limit [Fig. 4(a)]. This is consistent with the in-plane critical field exceeding the Pauli limit reported by recent experimentsFatemi et al. (2018); Sajadi et al. (2018). Importantly, this field-induced $B_{u}^{\prime\prime}$ phase also exhibits two majoranas localized near opposite corners [Fig. 4(b)(c)]. Due to the broken time-reversal symmetry, these two corner modes are single Majoranas instead of Majorana Kramers pairs [SM Sec. III]. We therefore emphasize that even if the superconductivity in the realistic WTe2 belongs to the even-parity $A_{g}$ representation, it is stil possible to obtain single Majorana corner modes by applying an in-plane field.

*A recipe for 2D higher-order superconductors—

*We point out that these 2D higher-order superconducting phases can in fact be achieved by a general recipe555There could be more recipes to achieve such a phase.: Our studies on WTe2 suggests that corner majoranas might occur generically from the combination of a gated QSH state with odd parity superconductivity. This recipe is most intuitive from the boundary perspective. Consider such a QSH normal state at a doping level where it still exhibits counter-propagating modes well-localized on the edge. In the absence of pairing, the corresponding BdG Hamiltonian has two electron-like and two hole-like zero-energy eigenstates with edge-localized wavefunctions. When we introduce an odd-parity pairing potential, which changes sign in real space when projected onto opposite edges and inevitably vanishes at the domain walls, the electron- and hole-like edge states will mix and acquire finite energies except at the two inversion-related points where the projected pairing vanishes. The resulting “leftover” zero-energy modes, whose point-like wavefunctions will likely to be trapped at corners for realistic samples, therefore lead to two Majorana Kramer’s pairs localized on two opposite corners. Although the bulk-boundary correspondence is not rigorously proven, we analytically show that corner majoranas naturally exist in a minimal model we construct for superconductors built from our recipe [SM Sec. VII].

Bulk invariant perspective of the recipe— This 2D higher-order topological superconducting state is in fact a type of TCsc protected by inversion symmetry. Based on studies of various symmetry-protected topological phasesFu and Kane (2007); Khalaf et al. (2018); Ono et al. (2018), we conjecture that the bulk topology in inverion-protected TCscs could be inferred from the inversion eigenvalues of occupied BdG bands at TRIMs. With these BdG parity data, we define a symmetry indicator as the bulk invariant for a 2D inversion-protected TCsc in the presence of time-reversal symmetry

[TABLE]

inspired by indicators proposed for 3D systemsKhalaf et al. (2018); Ono et al. (2018). Here $\xi_{\textbf{k},n}$ are the parity eigenvalues of the occupied BdG bands at TRIMs kFu and Kane (2007), and this indicator is stable to adding trivial normal bands for restricted cases where the normal state is half-filled. Application of this formula thus requires extending the “normal” inversion operator $I_{0}$ to Nambu space. For odd-parity superconductors, which are defined by superconducting gaps satisfying $I_{0}\Delta_{\textbf{k}}I_{0}^{-1}=-\Delta_{\textbf{k}}$, the operator I=\rm diag$$(I_{0},-I_{0}) defines the inversion operation for BdG Hamiltonians. For even-parity superconductors, this inversion operator $I$ has no minus sign in the hole part, so $\kappa$ is always 0Ono et al. (2018). By identifying trivial BdG parity data as those from “atomic superconductors”, which are constructed by placing zero-dimensional electron- or hole-like bogoliubons at Wycoff positions, we can see that our indicator is stable upon mod 4. This indicates that the classification of 2D inversion-protected TCsc is $Z_{4}$.

To identify which of the four states features corner Majoranas, we relate our index $\kappa$ to the well-known $Z_{2}$ index $\nu$ for 2D time-reversal superconductors[SM Sec. IV]:

[TABLE]

It is thus clear that $\kappa=0,2$ phases do not have edge Majoranas while $\kappa=1,3$ phases do. Nonetheless, the $\kappa=2$ phase is topologically distinct from the trivial $\kappa=0$ phase, hinting that the former has corner Majoranas.

In fact, the phases hosting corner Majoranas in WTe2 has $\kappa=2$, which we explicitly verified using $H^{\rm BdG}$ with self-consistently obtained $\Delta$ in Fig. 2(d) [SM Sec. V]. Not only for this particular example, here we show that general 2D higher-order superconductors constructed from our recipe have $\kappa=2$. To see this, we relate the $Z_{4}$ indicator $\kappa$ for a time-reversal parity-odd BdG system to the $Z_{2}$ topological index $\nu_{N}$Kane and Mele (2005b) for its normal state:

[TABLE]

Here, $\kappa_{N}=0,1$ is the $Z_{2}$ indicator defined analogously as in Eq. 6 but for normal-state HamiltoniansOno et al. (2018). Importantly, the latter relation holds for a metallic state only when the numbers of occupied bands are all the same at all TRIMs [SM Sec. VI]. Now, we follow our recipe and take the normal state to be a gated QSH state whose FS does not circle any TRIM, just as gated WTe2. In this case, Eq. 8 holds and we have $\kappa_{N}=\nu_{N}=1$. Upon introducing odd-parity pairing, the resulting superconductor therefore has $\kappa=2$.

Discussion— For the odd-parity paired states we find in WTe2, which we find to be inversion-protected higher-order TCsc, we expect that the Majorana corner modes cannot be removed without closing the bulk gap if inversion is preserved. When the inversion symmetry is broken, while the Majoranas are no longer protected by the 2D bulk topology, they are still protected by the gaps on the 1D edges. In this case the paired state becomes the so-called “extrinsic” higher-order topological superconductorGeier et al. (2018). We thus expect these Majorana corner modes can in principle be probed by STM or transport measurements.

Acknowledgment—This work is supported by Microsoft and Laboratory for Physical Sciences. R.-X.Z. is supported by a JQI Postdoctoral Fellowship. JS was supported by the NSF-DMR1555135 (CAREER). The authors acknowledge the University of Maryland super-computing resources (http://hpcc.umd.edu) made available for conducting the research reported in this paper. This research was supported in part (through helpful discussions at KITP) by the National Science Foundation under Grant No. NSF PHY-1748958. Note added– After posting this work we became aware of Ref. Skurativska et al., 2019, which mainly discussed the formulation of symmetry indicators for inversion-protected TCsc in any $d$ dimension. Their $d=2$ case agrees with our conjecture in Eq. 6 for the cases we focus on.

I I.   Normal state and interacting Hamiltonians

For the normal state, we adapt the first-principles-derived model reported in Ref. Ok et al., 2018 for the normal-state Hamiltonian $H_{0}$ [see Eq. 1 in the main text] of monolayer WTe2. The spin-degenerate part $\hat{h}_{0}(\textbf{k})$ in $H_{0}$ is given by

[TABLE]

in the basis of $\hat{\sigma}\otimes\hat{l}$. Here the momentum dependence is contained in the functions $\varepsilon_{l}(\textbf{k})=\mu_{l}+2t_{l}\cos(k_{x})+2t^{\prime}_{l}\cos(2k_{x})$, $f_{\textbf{k}}=1-e^{-ik_{x}}$, $g_{\textbf{k}}=(1+e^{-ik_{x}})e^{ik_{y}}$. Following Ref. Ok et al., 2018, we fix the tight-binding parameters (in eV) as: $\mu_{d}=0.4935,\mu_{p}=-1.3265,t_{d}=-0.28,t^{\prime}_{d}=0.075,t_{p}=0.93,t^{\prime}_{p}=0.075,t_{0}^{AB}=1.02,t_{d}^{AB}=0.52,t_{p}^{AB}=0.40,V_{soc}=0.115$.

For the interacting part, we consider short-ranged density-density interactions with the following explicit form:

[TABLE]

where $\sigma$ and $l$ indices are summed over, $n_{s\sigma l}(\textbf{r})$ is the density with spin $s$ and orbital $l$ locating at sublattice $\sigma$ in the unit cell centered at r, $n_{\sigma l}(\textbf{r})=\sum_{s}n_{s\sigma l}(\textbf{r})$, and $\boldsymbol{\delta}_{p}=0$, $\boldsymbol{\delta}_{d}=\hat{y}$. As mentioned in the main text, $U^{l}$ denotes the on-site interactions for orbital $l$, and $V_{1}^{l\tilde{l}}$ ($V_{2}^{l\tilde{l}}$) denotes the nearest-neighbor (next nearest-neighbor) interactions on the zigzag chains with intra- or inter-orbital characters for $\tilde{l}=l$ and $\tilde{l}=\bar{l}$, respectively.

II II.   Competition between even- and odd-parity paired states

Given the understanding about the real space configurations discussed in the main text, we now attempt to understand the phase diagram in Fig. 2(b) in the main text within a simplified “two-patch” scheme, where we ignore the intra-pocket momentum-dependence and consider effectively two points located at $\textbf{k}=(\pm k_{F},0)$ instead of two pockets. Such simplification works well in small-pocket limit, and is similar in spirit to the Eliashberg formalism, where the momentum-dependence is assumed to be uniform within a pocket. Here we consider only the interactions responsible for the most dominant components in the self-consistency solutions with $A_{g}$ and $B_{u}$ symmetries [see the bonds with $\Delta_{1}$ in Fig. 2(c) and (d) in the main text], i.e. $U^{d}$ and $V_{2}^{dd}$ terms in Eq. 2 in the main text. Such interaction has a simple momentum dependence of $\mathbb{V}(\textbf{q})=U+V\cos(q_{x})$, where we set $U^{d}=U$, $V_{2}^{dd}=V$, and q the momentum transfer. Thus the interaction matrix $\Gamma$ in the linearized gap equation on p.3 in the main text can be simplified into a $2\times 2$ matrix in the basis of incoming and outgoing momenta running over $\textbf{k}^{(^{\prime})}=(\pm k_{F},0)$:

[TABLE]

The two eigenvalues $\frac{\mathbb{V}(0)\pm\mathbb{V}(2k_{F})}{2}$ correspond to the parity-even and odd eigenvectors $(1,1)$ and $(1,-1)$ respectively, and the eigenvector with the more negative eigenvalue corresponds to the dominant pairing gap. Evidently, repulsive $V$ forbids the parity-odd $A_{u}/B_{u}$ whereas repulsive $U$ forbids the parity-even $A_{g}$, given that $\cos(2k_{F})<0$. This is true for a wide range of $k_{F}$ values, including that from the model in Eq. 1 in the main text. As for the cases where $U$ and $V$ are both attractive, the balance between $A_{g}$ and $B_{u}$ is tilted by the sign of the interaction with large momentum transfer $\mathbb{V}(2k_{F})$. To be precise, when the $2k_{F}$ component contributed by $V$ dominates over the momentum-independent $U$ ($\mathbb{V}(2k_{F})=U+V\cos(2k_{F})>0$), the odd-parity $B_{u}$ is favored over the even-parity $A_{g}$. For the lightly gated WTe2 model $H_{0}$ with $\mu=0.5$, the two Fermi pockets are located at $\textbf{k}=(\pm k_{F},0)$ with $k_{F}\sim 1.2$.

III III.   Majorana corner modes in various spin-triplet phases

In this section, we show more numerical evidences for Majorana corner modes in different spin-triplet phases we find in WTe2 phase diagram, the field-free phases $B_{u}$ and $A_{u}$, and the field-induced equal-spin phase $Bu^{{}^{\prime\prime}}$, as examples for our recipe (as these phases all clearly satisfy our recipe).

We first show that the corner-localized Kramers’ doublets we find for the $B_{u}$ phase [see Fig. 3 in the main text, also duplicated here in Fig. 5(a)(c)] indeed have zero energies by studying the finite-size scaling of the low-lying eigenvalues of the BdG Hamiltonian. Since we are interested only in the spectrum near zero, here we can partially diagonalize $H^{\rm BdG}$ using Lanczos techniques on extremely large lattices with open boundary conditions along $\hat{x}$ and $\hat{y}$ directions. It is clear from Fig. 5(b) that with an increasing lateral dimension $L$ of lattices with open boundary conditions along $\hat{x}$ and $\hat{y}$ directions, the corner-localized bound states tend exponentially toward zero energy with an increasing gap to the lowest-lying bulk quasiparticle excitations. In Fig. 5(d), we show the probability distribution of the lowest supragap state $|\psi_{2}|^{2}$ extends over the bulk. This is in contrast to the lowest eigenstate $|\psi_{0}|^{2}$ we showed in Fig. 3(b) in the main text [duplicated here in Fig. 5(c)], which sharply localized to the upper-left and bottom-right corners.

We now turn to the other spin-triplet phase $A_{u}$ in the superconducting phase diagram for WTe2 in Fig. 2b in the main text. Since the normal state is a gated QSH material with Fermi pockets all away from TRIMs [see Fig. 2a in the main text] and the pairing is spin-triplet, phase $A_{u}$ satisfies our recipe and should host Majorana Kramers’ pairs localized on opposite corners, just as phase $B_{u}$. Indeed, in Fig. 6(a) we show that the BdG spectrum indeed exhibits two zero-energy Kramers’ pairs, and the finite-size scaling of low-lying eigen-energies in Fig. 6(b) further confirms that the doublets’ energy $\epsilon_{0}\rightarrow 0$. Moreover, just as the $B_{u}$ case, the density distribution of these zero-energy states $|\psi_{0}|^{2}$ are sharply localized to two opposite corners, and that of the lowest supragap state $|\psi_{2}|^{2}$ extends over the bulk [see Fig. 6(c) and (d) respectively].

As we discussed in the main text, an in-plane magnetic field can drive a phase transition from the spin-singlet phase $A_{g}$ to the equal-spin phase $Bu_{{}^{\prime\prime}}$ [see Fig. 4 in the main text]. Here, we further show the log-scale density profile for the single Majorana modes we find in phase $Bu_{{}^{\prime\prime}}$ [see Fig. 7(c)] to demonstrate their localization properties. For a localized state $\psi$, we expect $\left|\psi(\vec{r})\right|^{2}\leq C\exp\left(\left|\vec{r}-\vec{r}_{0}\right|/\xi\right)$, where $\vec{r}_{0}$ is the localization center, $\xi$ a characteristic localization length, and $C$ some constant. Therefore, from Fig. 7(c) we can better see that the wavefunction of a Majorana mode is exponentially localized to the two opposite corners from the roughly linear change in the color scale.

Finally, in Fig. 8 we present another set of results for the Majorana corner modes in phase $Bu_{{}^{\prime\prime}}$ on a rectangular geometry with different dimensions to show the robustness of our findings. In contrast to the time-reversal symmetric phases $A_{u}$ and $B_{u}$, the time-reversal broken BdG spectrum here shows two single zero-energy modes instead of two Kramers’ doublets [see Fig. 8(a)]. In Fig. 8(b), we further show the finite-size scaling for the lowest two eigen-energies to confirm that the BdG spectrum is gapped and that Majorana modes indeed have zero energies. We then show that the Majorana modes are localized near two opposite corners [see Fig. 8(c)], and that the lowest gapped states are bulk states [see Fig. 8(d)]. It is interesting to find that the latter contains a remnant of the QSH edge states in the absense of pairing.

IV IV.   Relating the Z4 indicator and Z2 index for 2D superconductors

In this section, we explain how to obtain Eq. 5 in the main text. Consider a time-reversal and inversion-symmetric superconductor in 2D. Due to the presence of inversion symmetry, the $Z_{2}$ topological index $\nu$ for 2D time-reversal superconductors (class DIII) can be written in the following way:

[TABLE]

where $\xi_{\textbf{k}m}$ is the parity of the filled BdG band $m$ at a time-reversal invariant momentum (TRIM). Here, the parity is the eigenvalue of a BdG band for the BdG inversion operator $I$ defined below Eq. 4 in the main text, $n_{\textbf{k}}$ is the total number of filled bands at k, and $n_{\textbf{k}}=n$ is independent of k for superconductors. Only half of the bands enter the product since the system is two-fold degenerate. When $(-1)^{\nu}=1$ (0), the system is topological (trivial). Similar to the normal state case studied previouslyKhalaf et al. (2018), we can also express this index directly as

[TABLE]

where $n^{\pm}_{\textbf{k}}$ is the number of occupied BdG bands with parity $\pm 1$.

For cases where the normal state is half-filled, one can also define a Z4 symmetry indicator $\kappa$ for 2D inversion-symmetric superconductors (in the presence of time-reversal symmetry) the way we do in Eq. 4 in the main text. To relate $\kappa$ to $\nu$, we further write it as

[TABLE]

Here we have used the fact that $\sum_{m=1}^{n}\xi_{\textbf{k}m}=n^{+}_{\textbf{k}}-n^{-}_{\textbf{k}}=n-2n^{-}_{\textbf{k}}$. Now it is clear that

[TABLE]

since there are four TRIMs in 2D and that $n$ is even due to the two-fold degeneracy. We therefore conclude that

[TABLE]

which is the Eq. 5 in the main text.

V V.   Parity eigenvalues for BdG bands

To calculate the symmetry indicator $\kappa$ for the BdG Hamiltonian of monolayer WTe2 with the odd-parity $B_{u}$ pairing, we need to calculate the parity eigenvalues of the filled BdG bands according to the definition of $\kappa$ [see Eq. 4 in the main text]. We numerically obtain these eigenvalues at each of the high-symmetry points for our model $H_{0}$ with the self-consistently obtained $B_{u}$ pairing, as shown in Table 2.

VI VI.   Relating the Z2 indicator and Z2 index for topological ‘metals’

In this section, we explain why the second equation in Eq. 5 in the main text holds for topological ‘metals’ that are ‘effectively gapped’. We first review the insulator case presented in Ref. Khalaf et al., 2018. Consider an insulator with both time-reversal and inversion symmetries. Very similar to the DIII superconductor case in section IV., the Z2 index $\nu_{N}$ for time-reversal topological insulators can be calculated simply byFu and Kane (2007)

[TABLE]

where $\xi_{\textbf{k}m}$ is the parity of the filled normal band $m$ at a time-reversal invariant momentum (TRIM). Here, the parity is the eigenvalue of a normal band for the ‘normal’ inversion operator $I_{0}$ defined in the main text, $n_{\textbf{k}}$ is the total number of filled normal bands at k, and $n_{\textbf{k}}=n$ is independent of k for insulators. Only half of the bands enter the product since the system is two-fold degenerate. When $(-1)^{\nu_{N}}=1$ (0), the system is topological (trivial). As shown in Ref. Khalaf et al., 2018, we can also express this index directly as

[TABLE]

where $n^{\pm}_{\textbf{k}}$ is the number of occupied normal bands with parity $\pm 1$.

Due to the presence of inversion symmetry, one can also calculate the Z2 symmetry indicator $\kappa_{N}$ for inversion-protected topological crystalline insulators given by

[TABLE]

Here we have used the fact that $\sum_{m=1}^{n_{\textbf{k}}}\xi_{\textbf{k}m}=n^{+}_{\textbf{k}}-n^{-}_{\textbf{k}}=n_{\textbf{k}}-2n^{-}_{\textbf{k}}$. Now it is clear that

[TABLE]

if $\alpha\equiv\sum_{\textbf{k}\in\rm{TRIM}}\frac{n_{\textbf{k}}}{4}$ is even. Since both indices are $Z_{2}$, we have $\nu_{N}=\kappa_{N}$.

For insulators, $\alpha=n$(number of TRIMs)$/4$, where $n$ is even due to the two-fold degeneracy. Eq. 26 in the Supplementary Material is thus not guaranteed for 1D insulators, but holds for both 3D and 2D insulators, where the numbers of TRIMs are 8 and 4, respectively. As for metals with Fermi pockets enclosing any TRIM, $n_{\textbf{k}}\neq n$ is different for different TRIMs in general, and thus $\alpha$ is not guaranteed even. Nonetheless, for metals whose Fermi pockets are away from TRIMs, $n_{\textbf{k}}=n$ is still k-independent, and thus Eq. 26 in the Supplementary Material still holds in 2D and 3D.

VII VII.   Bulk-boundary correspondence for $\kappa=2$ inversion-protected TCsc

As mentioned in the main text, the fact that a $\kappa=2$ inversion-protected TCsc in 2D has no Majorana edge modes but is still topologically distinct from the trivial $\kappa=0$ phase suggests the possibility of Majorana corner modes. Such bulk-boundary correspondence for a higher-order topological superconductor, however, has not been rigorously proven to the best of our knowledge. Thus to gain more understanding about this bulk-boundary correspondence, in this section we will first present a minimal model for a $\kappa=2$ TCsc in 2D, which can be tuned across phase boundaries to $\kappa=3$ and 4 (trivial) phases. Then within this model, we will show analytically how zero-dimensional Majorana modes arise on the boundary of the $\kappa=2$ phase when placed against the trivial phase.

VII.1 A.   Minimal model for an inversion-protected TCsc

Our eight-band minimal model $\mathcal{H}=\sum_{\textbf{k}}\mathcal{H}(\textbf{k})$

[TABLE]

consists of a regularized Bernevig-Hughes-Zhang (BHZ) like model for a QSH state with odd-parity pairing. Here $\hat{\tau}$, $\hat{s}$, and $\hat{\rho}$ are Pauli matrices for particle and hole, spin $s=\uparrow,\downarrow$, and orbital $\rho=s,p_{-}$. This model $\mathcal{H}$ is invariant under time reversal operation $\Theta=is_{y}\mathcal{K}$, $\textbf{k}\rightarrow-\textbf{k}$, particle-hole transformation $\mathcal{P}=\hat{\tau}_{x}\mathcal{K}$, $\textbf{k}\rightarrow-\textbf{k}$, and the inversion operation defined for odd-parity superconductors I $=\hat{\tau}_{z}\otimes\hat{\rho}_{z}$, $\textbf{k}\rightarrow-\textbf{k}$.

This minimal model exhibits two bulk topological transitions in the parameter space of $\epsilon_{0}$, $m_{0}$, and $m_{1}$: one is from the inversion symmetry indicator $\kappa=2$ to $\kappa=3$, and the other from $\kappa=3$ to $\kappa=4$. In particular, we can tune through different phases by tuning $m_{1}$ at a fixed $m_{0}$ and $\epsilon_{0}\geq 0$ as follows

[TABLE]

where all band inversions occur at $\Gamma$. Taking the topologically trivial $\kappa=4$ phase as the reference point, the spectrum first undergoes a single band inversion to enter the $\kappa=3$ phase, then follows another band inversion to enter the $\kappa=2$ phase.

VII.2 B.   Edge modes in a three-domain geometry

To understand what kind of boundary modes a $\kappa=2$ phase can host when placed against a trivial phase, we first study the edge modes in a geometry of concentric rings with three domains along the radial direction $r$: $\kappa=2$ phase for $r<R_{a}$, $\kappa=3$ phase for $R_{a}<r<R_{b}$, and $\kappa=4$ phase for $r>R_{b}$ [see Fig. 9 (a) in Supplementary]. Since $\kappa$ is defined modulo 4, the phase $\kappa=4\equiv 0$ on the outside is the trivial phase. Given that all band inversions in $\mathcal{H}$ happen at $\Gamma$, it is more convenient to realize this geometry by working with the $k\cdot p$ model $\mathcal{H}_{\Gamma}(\textbf{k})$ around $\Gamma$ written in the polar coordinate $(r,\theta)$. We focus on the $v\ll\Delta$ limit to avoid the possible existence of edge modes, in which case the existence of the submerged corner modes becomes ambiguous. Together with the fact that $s=\uparrow,\downarrow$ is a good quantum number, we arrive at the Hamiltonian for each spin species in the $\hat{\tau}\otimes\hat{\rho}$ basis

[TABLE]

where $k_{r}=-i\partial_{r}$ and $m_{\pm}(r)\equiv\epsilon_{0}(r)\pm(m_{0}(r)+2m_{1}(r))$. Here we have set $1/r\rightarrow 0$ since we focus on the asymptotic form of the edge modes (i.e. consider $|r-R_{a,b}|\gg 0$). The three-domain geometry can then be realized by setting

[TABLE]

Since the $\kappa=3$ superconducting phase in the middle domain has a nontrivial $Z_{2}$ index $\nu=1$, we expect to find zero-energy eigenstates $\psi^{s}_{23}(r,\theta)$ and $\psi^{s^{\prime}}_{34}(r,\theta)$ localizing along the domain walls at $r=R_{a}$ and $r=R_{b}$, respectively. By taking $m_{+}(r)=$sgn$(r-R_{a})$, $m_{-}(r)=-$sgn$(r-R_{b})$, and $\Delta>0$ for simplicity, we find one pair of helical edge modes per domain wall [see Supplementary Material Sec. IV A]:

[TABLE]

and

[TABLE]

where $l$ is the orbital angular momentum taking half integers. Given that $\psi^{s}_{23/34}(r,\theta)$ obeys Majorana condition up to an overall phase and that the spin-up modes $\psi^{\uparrow}_{23/34}$ and spin-down modes $\psi^{\downarrow}_{23/34}$ propagate along the domain walls in opposite directions [see Supplementary Material Sec. IV B], we have arrived at one pair of helical Majorana edge modes per domain wall [see Fig. 9 (a) in Supplementary Material].

VII.3 C.   Majorana corner modes in a $\kappa=2$ phase

We are now ready to study the boundary modes between a $\kappa=2$ and a $\kappa=4$ phase. In the following we will shrink the $\kappa=3$ domain by bringing $R_{a}\rightarrow R_{b}$, and see if there exists any symmetry-allowed perturbation that gaps out the edge modes $\psi^{s}_{23}$ and $\psi^{s^{\prime}}_{34}$ at the two domain walls.

We first write down the rotational invariant perturbations $\mathcal{H}^{\prime}_{\rm rot}(r,\theta)$ up to linear oder, and project them onto the edge modes from the two domain walls [see Supplementary Material Sec. IV C]. We find that $\mathcal{H}^{\prime}_{\rm rot}(r,\theta)$ only couples edge modes propagating in the same direction, i.e. $\psi^{s}_{23}$ and $\psi^{s}_{34}$. The two pairs of helical edge modes therefore remain gapless in the presence of rotational invariant perturbations.

Next we consider the rotational-breaking perturbations. Here we focus on the lowest order terms $\mathcal{H}^{\prime}$, which have no spatial dependence. After projecting all the symmetry-allowed terms onto the edge modes $\psi_{23}(r,\theta)$ and $\psi_{34}(r,\theta)$, we find that there are only two non-vanishing terms that are hermitian and couple components in the edge modes that are counter-propagating

[TABLE]

Both of these terms are pairing terms that are odd under inversion, i.e. $[\mathcal{H}^{\prime},\textit{I}~{}]=0$, where I is the inversion operator for odd-parity superconductors defined earlier. Moreover, their corresponding amplitudes after projection $\tilde{\mathcal{H}}^{\prime}_{a,b}(\theta)=\int dr\psi_{23}^{\dagger}(r,\theta)\mathcal{H}^{\prime}_{a,b}\psi_{34}(r,\theta)$ have an angular dependence of $\tilde{\mathcal{H}}^{\prime}_{a}(\theta)\propto i\sin\theta$ and $\tilde{\mathcal{H}}^{\prime}_{b}(\theta)\propto i\cos\theta$, respectively. The full Hamiltonian of the edge modes therefore has the form

[TABLE]

where $\alpha$ and $\beta$ are real numbers, and $\hat{\lambda}$, $\hat{s}$ are Pauli matrices for edge modes $\psi^{s}_{23}$, $\psi^{s}_{34}$, and spin $s=\uparrow,\downarrow$, respectively. Since the back-scattering term $\alpha\cos\theta+\beta\sin\theta$ has opposite signs at any $\theta$ and $\theta+\pi$, it is bound to vanish at $\theta_{0}=\tan^{-1}(\frac{\beta}{\alpha})+\frac{\pi}{2}$ and $\theta_{0}+\pi$. This means that even at the lowest order, rotational breaking perturbations will gap the edge modes $\psi_{23}$ and $\psi_{34}$ in an odd-parity way, and a $\kappa=2$ TCsc will host at least two ‘leftover’ zero-dimensional zero-energy Kramer’s pairs located at $\theta_{0}$ and $\theta_{0}+\pi$ in an open geometry [see Fig. 9 (b) in Supplementary Material]. The specific value of $\theta_{0}$ is given by the microscopics, and these two zero-dimensional Majorana Kramers pairs are often trapped at the opposite corners of the considered geometry. Importantly, the two Majorana pairs can annihilate each other only when inversion symmetry is broken. Thus within our minimal model for an inversion-protected TCsc with $\kappa=2$, we have shown how zero-dimensional Majorana Kramers pairs emerge on the boundaries, and hence an inversion-protected higher-order topological superconductor.

VIII VIII.   Details regarding derivations in section V

VIII.1 A.   Edge modes in the three-domain geometry

In this subsection, we show how we obtain the Majorana edge modes $\psi^{s}_{23}(r,\theta)$ and $\psi^{s^{\prime}}_{34}(r,\theta)$ in the three-domain geometry. Consider the domain wall between the $\kappa=2$ and $\kappa=3$ phases. The Hailtonian $\mathcal{H}_{\Gamma}^{s}$ in Eq. 33 in Supplementary Material section III. becomes block diagonal after taking the limit $v\rightarrow 0$, and it is clear from Eq. 33 in Supplementary Material section III. that only the block with $m_{+}(r)$ terms experiences a sign-changing mass term and is thus expected to trap a zero-energy eigenstate $\psi^{\uparrow}_{23}(r,\theta)$ localized at $r=R_{a}$. In the following, we solve for the asymptotic for m of $\psi^{\uparrow}_{23}(r,\theta)$. This $2\times 2$ block (represented by Pauli matrix $\tau$) with $m_{+}(r)$ term for the spin-up sector is given by

[TABLE]

The zero-energy eigenstate is then given by $h_{+}^{\uparrow}\psi^{\uparrow}_{23}=0$, but we can obtain the edge state more conveniently by solving the zero-energy mode for a rotated Hamiltonian $\tilde{h}_{+}^{\uparrow}\phi=0$, where $\tilde{h}_{+}^{\uparrow}=\hat{\tau}_{x}\hat{U}_{\theta}h_{+}^{\uparrow}\hat{U}^{\dagger}_{\theta}$, and $\hat{U}_{\theta}=e^{i\hat{\tau}_{z}\theta/2}$ is a unitary transformation. After some algebra, we find that

[TABLE]

Now $\tilde{h}_{+}^{\uparrow}$ is $\theta$-independent, and we can thus write its zero-energy eigenstate as $\phi(r)=f(r)\xi$, where the spinor $\xi$ obeys $\hat{\tau}_{y}\xi=a\xi$. For simplicity, we take $m_{+}(r)=$sgn$(r-R_{a})$ and $\Delta>0$. Then by solving the differential equation for $f(r)$ and requiring that $f(r)$ localized at the boundary $r=R_{a}$, we find that $a=1$ and $\phi(r)=e^{-1/\Delta|r-R_{a}|}(1,0,i,0)^{T}$. Together with the fact that the orbital part of the angular momentum has the form $L_{z}=-i\partial_{\theta}$, we can then obtain the zero-mode for $h_{+}^{\uparrow}(r,\theta)$ with angular dependence

[TABLE]

Note that the Majorana condition $\Xi\psi^{\uparrow}_{23}=\psi^{\uparrow}_{23}$ is satisfied up to an overall phase. $\psi^{\downarrow}_{23}(r,\theta)$ and $\psi^{\uparrow/\downarrow}_{34}(r,\theta)$ can be obtained in a similar way.

VIII.2 B.   Propagating directions of the edge modes

To determine the propagating directions of these edge modes, we need to take back the terms containing $k_{\theta}=-i\partial_{\theta}$ in the $k\cdot p$ Hamiltonian $\mathcal{H}^{\uparrow/\downarrow}_{\Gamma}(\textbf{k})$ and consider them as perturbations. These are the terms that emerge when we write $k_{\pm}=k_{r}+ik_{\theta}$. In the Hamiltonian $\mathcal{H}^{\uparrow/\downarrow}_{\Gamma}(r,\theta)$ in Eq. 33 in Supplementary Material, we did not include these terms since we expect them to be small in the limit of large $r$. Such perturbation has the form

[TABLE]

Recall that the edge modes are given by

[TABLE]

and

[TABLE]

where $l$ is the orbital angular momentum taking half integers. We thus find the energy correction arising from rotational motion in $\theta$ to be

[TABLE]

where $j=23,34$ denotes the two domain walls. We therefore conclude that the spin-up and spin-down modes are right- and left-movers respectively for both domain walls, which amounts to one pair of helical edge modes per domain wall [see Fig. 9 (a) in Supplementary Material].

VIII.3 C.   Rotational invariant perturbations

In this subsection, we write down the general form of rotational invariant perturbations $H^{\prime}_{\rm rot}(\textbf{k})$ and its projection on to the edge modes $\psi^{s}_{23}$ and $\psi^{s^{\prime}}_{34}$. The rotational operator is given by $C_{\theta}=e^{-iJ_{z}\theta/2}$, $k_{\pm}\rightarrow e^{\mp i\theta}k_{\pm}$, where the angular momentum has the form $J_{z}=\hat{\tau}_{z}\otimes\hat{s_{z}}\otimes\hat{\sigma_{z}}/2$, and $k_{\pm}=k_{x}\pm ik_{y}$. The perturbations that preserve time-reversal symmetry $\Theta$, particle-hole symmetry $\mathcal{P}$, inversion symmetry for odd-parity pairing I, and also obey $C_{\theta}H^{\prime}_{\rm rot}(\textbf{k})C_{\theta}^{\dagger}=H^{\prime}_{\rm rot}(\textbf{k})$ have the general form of

[TABLE]

Here, $H^{\prime}_{\rm rot}(\textbf{k})$ is written in the basis of $\hat{\tau}\otimes\hat{s}\otimes\hat{\rho}$, and $m_{1/2}$, $D_{0}$, $A$, $B$, and $C$ are free parameters.

Since the edge modes with the same spin (opposite spins) $\psi^{s}_{23}$ and $\psi^{s}_{34}$ ($\psi^{\bar{s}}_{34}$) propagate in the same direction (opposite directions), $H^{\prime}_{\rm rot}$ has to couple edge modes with opposite spins in order to create a gap. However, it is obvious from Eq. 78 in Supplementary Material that there exists no spin-flipping terms that are allowed by the symmetries considered. The edge modes $\psi^{s}_{23}$ and $\psi^{s^{\prime}}_{34}$ thus remain gapless in the presence of rotational invariant perturbations.