Time-reversal-invariant $C_2$-symmetric higher-order topological superconductors

We propose a minimal lattice model for two-dimensional class DIII superconductors with $C_2$-protected higher-order topology. While this class of superconductors cannot be topologically characterized by symmetry eigenvalues at high symmetry momenta, we propose a simple Wannier-orbital-based real-space diagnosis to unambiguously capture the corresponding higher-order topology. We further identify and characterize a variety of conventional topological phases in our minimal model, including a weak topological superconductor and a nodal topological superconductor with chiral-symmetry protection. The disorder effect is also systematically studied to demonstrate the robustness of higher-order bulk-boundary correspondence. Our theory lays the groundwork for predicting and diagnosing $C_2$-protected higher-order topology in class DIII superconductors.

Introduction -Anyons are stable exotic quasiparticles with unconventional statistical braiding properties and serve as the cornerstone for topological quantum computation [1][2][3][4]. The most well-known anyonic quasiparticle is the zero-dimensional (0D) Majorana zero mode (MZM), which could in principle emerge as a vortex bound state of a two-dimensional (2D) p + ip chiral topological superconductor (TSC) [5], or as the end mode of a one-dimensional (1D) p-wave TSC [6,7]. Remarkable experimental progress has been made towards realizing Majorana physics in superconducting Rashba nanowires [8,9] and recently in iron-based superconductors [10], where promising signatures of MZMs such as zero-bias tunnel conductance peaks have been observed. However, a "smoking-gun" measurement of confirming MZM's anyonic nature is missing. Therefore, the experimental existence of MZMs is still not conclusive, which calls for more efforts on topological Majorana physics, perhaps in different systems.
Recently, there has been growing interest in understanding higher-order topology [11][12][13][14], where anomalous boundary physics could show up on the (d − n)dimensional boundary of a d-dimensional topological system with 1 < n ≤ d. In particular, a 2D higher-order TSC is defined to host MZMs on its 0D geometric corners (d = n = 2 here) [15][16][17][18][19][20], which is a novel promising Majorana platform. Such corner-localized MZMs are proposed to exist in iron-based superconductors [21][22][23] and doped WTe 2 [24]. However, in most of these proposals, the corner MZMs are relatively "fragile" and can be easily removed upon closing the edge energy gap. The robustness of the corner MZMs, however, can be enhanced against symmetry-preserving edge perturbations if we can introduce additional crystalline-symmetry protections [14,24]. Several groups have proposed various symmetry indicators to classify such symmetry-protected higher-order TSCs based on symmetry eigenvalues at high symmetry momenta [25][26][27][28][29]. In this context, class DIII superconductors with two-fold rotational symmetry C 2 are special in the sense that they always host the same C 2 symmetry eigenvalues at high symmetry mo-menta for all possible topologically distinct phases. Consequently, this class of superconductors does not admit an indicator-based classification [30]. While the absence of an indicator does not necessarily imply the absence of topological phases, understanding and characterizing higher-order topology for this symmetry class remains an important open question in spite of considerable recent activities in the topological classification of superconductors.
In this work, we propose a minimal example of a 2D class DIII higher-order topological superconductor with C 2 symmetry protection. Despite its simplicity, our minimal model hosts surprisingly rich topological physics, including a higher-order TSC phase with 0D cornerlocalized Majorana Kramers pairs (MKPs), a weak TSC phase with 1D Majorana edge band, and a chiralsymmetry-protected nodal TSC phase with edge Majorana flat band. The robustness of corner MKPs with respect to disorder effect is systematically studied, unambiguously demonstrating the higher-order bulk-boundary correspondence. Motivated by the recently proposed Majorana counting rule for class D superconductors [31], we develop a real-space Wannier-orbital-based counting rule to predict and diagnose various time-reversal-invariant topological phases emergent in our model. In particular, our real-space counting predicts that the presence of corner MKPs only depends on the information of both lattice geometry and Bogoliubov-de Gennes (BdG) Wannier orbitals, and not the microscopic details of Cooper pairings. We confirm this remarkable prediction through explicit calculations using different types of Cooper pairing.
Majorana counting rule in class DIII -We first define our real-space characterization of higher-order topology for C 2 -symmetric class DIII superconductors. It is convenient to map the electron and hole degrees of freedom into a pair of Majorana fermions α R,s and β R,s , where R is the real-space position and s is the spin index. The Majorana operators necessarily come in Kramers pairs due to the time-reversal symmetry. This Majorana representation allows us to generalize the concept of "Kitaev limit" in Ref. [31] to the time-reversal-invariant version. Specifically, a time-reversal-invariant Kitaev limit is achieved for a class DIII superconductor when every bulk MKP formed by electron and hole Kramers pairs is connected to exactly one other MKP via Majorana bonds. However, crucially distinct from class D, the MKP's spin degrees of freedom enable two different types of Majorana bondings: • equal-spin bonding: Majorana fermions with the same spin indices are bonded, which effectively leads to p-wave Cooper pairing; • opposite-spin bonding: Majorana fermions with opposite spin indices are bonded, which effectively leads to s-wave/d-wave Cooper pairing.
For every Kitaev limit, there exists a Kramers pair of maximally localized BdG Wannier orbitals sitting at the center of each bond [32]. Importantly, the position of such Wannier orbital does not rely on the explicit bonding type. We now define a set of counting numbers where n M i (n W i ) is the number of the MKPs (Wannier orbital pairs) at the maximal Wyckoff position q i for i = a, b, c, d. We denote q 1a = (0, 0), q 1b = ( 1 2 , 0), q 1c = ( 1 2 , 1 2 ) and q 1d = (0, 1 2 ). Then our real-space diagnosis states that higher-order topology shows up only when ∆ b,c,d ≡ 1 (mod 2) [32].
We emphasize that while this Majorana counting rule is derived in the Kitaev limit, we expect it to hold for a general Wannierizable superconductor. The procedure for reducing a non-Kitaev superconductor to its Kitaev limit generally involves smooth deformations such as minimizing the localization length of Wannier functions and adiabatic movements of Wannier orbitals. Clearly, our Majorana counting rule is topologically immune to these deformations as long as they are adiabatic and symmetric.
Model Hamiltonian -We now introduce our minimal model to demonstrate the real-space diagnosis for C 2protected time-reversal-invariant higher-order topology, where l = A/B is the sub-lattice index, a and b are lattice constants along x and y directions, and σ i (i = 0, x, y, z) are the Pauli matrices for spins. The subscript (×) denotes the equal (opposite)-spin bonding. The Hamiltonian can be transformed to fermion operators using α l s,R = (c l † s,R + c l s,R )/ The bonding scheme of Eq. 2 within a unit cell is depicted in Fig. 1(a). In the fermion representation, the momentumspace Bogoliubov-de Genes Hamiltonian is where τ and γ are Pauli matrices representing the particle-hole and sub-lattice degrees of freedom. It is easy to check that H(k) is invariant under time-reversal operator Θ = iσ y K, particle-hole operator P = τ x K, and rotation The gap closing condition of H(k) is determined by three independent parameters: Physically, t 0 characterizes the net bonding strength between onsite MKPs, while t x and t y characterize the bonding strength for nearest-neighbor MKPs along x and y directions, respectively. The Wannier orbital sits at the center of the strongest Majorana bond. For example, if t x > t 0,y , we can immediately tell the existence of Wannier orbital pairs sitting at both 1c and 1d. We note that the position of Wannier orbitals is independent of the spin texture of the bondings.
Meanwhile, we are able to analytically map out the phase diagram in terms of t x /t 0 and t y /t 0 , as shown in Fig. 1(b). Below, we ascertain the topological nature of each phase by smoothly deforming into the corresponding Kitaev limit and applying our Majorana counting rule to perform the topological diagnosis. While such a Kitaev reduction is not necessary for the topological diagnosis, this procedure makes it easier to directly "read out" the position information of Wannier orbitals, greatly facili-  tating the application of the counting rule.
Higher-order TSC -In the regime t x > t 0 + t y , we deform the Hamiltonian to the Kitaev limit where only t x > 0 and t 0 , t y → with t x . As shown in Fig. 1(b), while the two MKPs within one unit cell both originate from maximal Wyckoff positions 1a and 1d, the Wannier orbital pairs sit at 1b and 1c. Therefore, the counting numbers are ∆ d = −1 and ∆ b,c = 1, and thus our Majorana counting rule predicts this gapped phase as a higherorder TSC. Due to the presence of both C 2 and timereversal symmetry, we expect the symmetry-protected higher-order bulk topology to correspond to two C 2related zero-energy MKPs emerging on the boundary.
To confirm the above predictions, we now simulate this phase numerically in an open-boundary finite-size system. Without loss of generality, we choose t x = t x× and t y = t y× to be opposite-spin Majorana bonding (s-wavelike pairing) and t 0 = t 0 to act as the chemical potential. For a more realistic manifestation of the corner modes, we introduce a small perturbation along the vertical direction.
which should slightly delocalize possible topological bound states in our system. As shown in Fig. 2(a), the energy spectrum on a 25 × 25 lattice shows 4 degenerate Majorana modes at zero energy. By examining their spatial profiles, the four Majorana modes are divided into two spatially sep-arated MKPs which are localized exponentially around the two C 2 -related geometric corners. This unambiguously establishes the C 2 -protected higher-order topology in this phase, confirming the prediction from the Majorana counting rule.
On the other hand, this higher-order TSC phase manifests itself as an example that goes beyond the symmetry-eigenvalue-based topological characterization scheme. This is simply because a C 2 -invariant momentum is also a time-reversal-invariant momentum in momentum space. For class DIII systems, each Kramers pair at every high symmetry momentum carries (+i, −i) as their C 2 eigenvalues. Therefore, the symmetry information for such systems is essentially "featureless" and remains the same across possible topological phase transitions. Indeed, a rather complete symmetry indicator theory for BdG systems in all space groups has recently concluded the absence of momentum-space indicators for this symmetry class [30]. By contrast, our Majorana counting rule is based on a real-space description, thus fundamentally avoiding the difficulty caused by the featureless symmetry data.
Weak TSC -In the regime t y > t x + t 0 , we adiabatically approach the Kitaev limit by sending t x , t 0 → with t y . The spatial information of the bulk MKPs and the Wannier orbitals is shown in Fig. 1(b). We immediately see that ∆ d = −1 while ∆ b,c = 0, indicating the absence of higher-order topology. Indeed, in this limit, the system resembles a set of decoupled y−directional 1D time-reversal-invariant TSCs [33] stacking along x direction, manifesting itself as a weak topological superconductor protected by the x−directional translational symmetry. As each 1D TSC hosts two end MKPs, the collection of edge MKPs in our weak TSC phase forms two pairs of time-reversal-invariant Majorana flat bands that are individually localized at the upper and lower edges. To see this, we numerically plot the energy spectrum of a N x × N y lattice in Fig. 2(b). Specially, each upper and bottom edge hosts 2N x localized zero-energy states, which confirms the existence of the Majorana flat bands.
If we break the translational symmetry by doubling the unit cell along x direction, the weak topology is destroyed. However, such unit-cell doubling leads to an occupation of one MKP and zero Wannier orbital pair for all four maximal Wyckoff positions. This leads to ∆ b,c,d = −1, which satisfies the higher-order topological condition. In the Supplementary Material [32], we numerically confirm the existence of higher-order topology by adding x−directional dimer bonds to the weak phase, which is consistent with the counting rule. This is a demonstration of the versatility of our extended counting rule that can be used as an efficient real-space diagnosis and a powerful tool for model construction.
Nodal TSC -For |t x − t 0 | ≤ t y ≤ t x + t 0 , the system becomes gapless and hosts four 2D Dirac points. The where k 0 = cos −1 (t 2 y − t 2 0 − t 2 x )/(2t 0 t x ) , θ 0 = Arg(t 0 + it 0× ), θ x = Arg(t x + it x× ) and θ y = Arg(t y + it y× ). For class DIII superconductors, the chiral symmetry C, a combination of time-reversal and particle-hole symmetry, allows us to define a chiral winding number ν ∈ Z [34,35] to characterize the topological nature of the Dirac point [32]. Therefore, we numerically calculate the winding number, and find it to be non-zero for all four Dirac points. As shown in Fig. 3(a), Dirac points I and IV share ν = 1, while Dirac points II and III have ν = −1, which is consistent with the time-reversal symmetry requirement in the system.
A nontrivial winding number necessarily implies the existence of edge Majorana flat band at zero energy [34,35]. In Fig. 3(b), we numerically calculate the energy spectrum for the nodal SC phase in a ribbon geometry. As expected from the bulk-boundary correspondence, edge Majorana flat bands lying at zero energy are found to connect Dirac points with opposite ν, e.g. (I) with (II), and (III) with (IV). The existence of both bulk nontrivial winding number and robust edge Majorana flat bands together establish this gapless phase as a chiral-symmetry-protected nodal TSC phase.
Disorder Effect -We now study the robustness of the    higher-order TSC phase against local disorder. First, we note that there are two important energy scales in our system in the t 0 → 0 limit: (i) the bulk gap 2|t x − t y | and (ii) the edge gap 2|t y |. From the higherorder bulk-boundary correspondence, we expect that the corner MKPs could be destroyed only when the disorder strength exceeds the bulk energy gap, if the disorder globally preserves the protecting symmetries. Numerically, we introduce chemical potential fluctuations by randomizing the strength of on-site bonding t 0 following a Gaussian ensemble, which is centered at zero and with a standard deviation of ∆t. For each value of ∆t, we average the density of states over 100 random configurations. We now consider two different sets of bonding parameters such that the two corresponding systems share the same edge gaps but differ in their bulk gaps. As shown in Fig. 4(a) and (b), the critical disorder strength that suppresses corner Majorana modes is close to the value of the bulk gap (black dashed lines), and not the edge gap (white dashed lines). This clearly demonstrates the bulk origin and the consequent stability of higher-order topology in our system.
Conclusion -We propose a minimal model Hamiltonian for higher-order topology in class DIII systems, which contains rather rich topological phenomena such as higher-order TSC phase, weak TSC phase, and nodal TSC phase. While a momentum-space characterization becomes inapplicable, we succeed in deciphering the higher-order topology in real space with our Majorana counting rule, by identifying the spatial information of BdG Wannier orbitals. Our proposed counting rule should serve as an important diagnostic principle in the search for higher-order-topology-based Majorana platforms.
Acknowledgement -R.-X. Z. thanks Sheng-Jie Huang for helpful discussions. This work is supported by the Laboratory for Physical Sciences and Microsoft Corporation. R.-X. Z. is supported by a JQI Postdoctoral Fellowship.

Majorana Bonding and Wannier orbitals
In this section, we discuss the relation between Majorana bonding and the position of resulting Wannier orbital pair. We assume two Majorana Krammer pairs (α s , α −s ) and (β s , β −s ) localized at R α and R β respectively. The interaction between the two pairs in the basis (α s , α −s , β s , β −s ) is where A = −A * due to the anti-commutation of Majorana fermions. T has vanishing diagonal matrix elements because time-reversal invariance forbids coupling between time-reversed partners. The Hamiltonian T thus has four solutions, corresponding to two Dirac fermions where u 1 , u 2 , v 1 , v 2 are two-component vectors. The position relative to the bond center is given by the matrix with U being unitary. As a result, we can define a transformation so that M ν 1 = ν 1 but {M, X} = 0. Thus (ν * 1 ) T Xν 1 = 0 and similarly, (ν * 3 ) T Xν 3 = 0. As a result, we have shown that the Wannier orbital pair is always localized at the middle of the bond, regardless of the bond type.

Majorana counting rule for class DIII systems
We provide a derivation of the Majorana counting rule for class DIII systems presented and utilized in the main text. This derivation is similar to that of the counting rule for class D systems [31]. In the Kitaev limit, we are able to "count" N , the number of dangling MKPs on the C 2 -symmetric boundary. In particular, the system is higher-order topological if N ≡ 1 (mod 2). N consists of two parts: (i) a positive contribution from the dangling Wannier orbital pairs N W ; (ii) a negative contribution N M from the dangling MKPs that break C 2 symmetry for a finite boundary and need to be removed. Notably, we only consider both N W and N M contributions from the ones at maximal Wyckoff positions q 1i for i = a, b, c, d. Thus, we define n M i and n W i as the number of MKPs and Wannier orbitals at q 1i within one unit cell, respectively. For a finite size system with L x × L y lattice sites, we find that which leads to Here we have defined the counting number ∆ i = n W i − n M i . On the other hand, we want to avoid possible weak topology, which leads to the constraint ∆ b ≡ ∆ c ≡ ∆ d (mod 2). Together with Eq. 11, we arrive at the counting rule presented in the main text, that the system is higher-order topological (i.e. N is odd) only if ∆ b,c,d ≡ 1 (mod 2). We calculate the energy spectrum and spatial density of zero-energy modes to demonstrate the higher-order topology of the system in (d). The inset shows a zoom-in plot of the corner MKPs.

Kitaev limit of the gapped phases
The time-reversal-invariant 2D Hamiltonian with l = A/B as the sub-lattice index is given by Without loss of generality, we take We can then define a local transformation to decouple the spinful system into two time-reversal-related copies: where O[θ] = cos θσ 0 + i sin θσ 2 is the rotation matrix, |α and |β denote the two-component spinors for brevity, and θ 0 = Arg(t x + it x× ), θ x = Arg(t x + it x× ), θ y = Arg(t y + it y× ). Under such a transformation, the Hamiltonian becomes a direct sum of two effectively "spinless" Hamiltonian H = h h where h =i m,n,l t 0 β l m,n α l m,n + t x β l m,n α l m+1,n + i m,n t y β A m,n β B m,n − t y α B m,n α A m,n+1 , and t 0 = t 2 0 + t 2 0× , t x = t 2 x + t 2 x× , t y = t 2 y + t 2 y× . The bulk gap closing of h only depends on three parameters t 0,x,y , instead of the previous six parameters in the original Hamiltonian. In particular, the system becomes nodal when |t 0 − t x | ≤ t y ≤ t 0 + t x .
In Fig. 5, we show the real-space Kitaev limit of each gapped phases of and the corresponding Wannier orbital pair configuration. The Kitaev limit is obtained by keeping the strongest bond, and deforming the weaker ones to zero. However, it should be emphasized that the unit cell can be enlarged by adding dimerizing bonds even though these bonds are weak, which could change the values of the counting numbers and could lead to a change of topology. For example, we consider adding an additional dimer bond t x along the x-direction to gap out the weak topological phase in (c), as shown in Fig. 5 (d), which is defined as In Fig. 5 (e) with t = 0.1, we find two corner-localized MKPs that signals the higher-order topology. This agrees with the prediction of the counting rule that ∆ b,c,d = −1.

Nodal superconductor
When |t 0 − t x | < t y < t 0 + t x , the system hosts four nodal points at where k 0 = cos −1 (t 2 y − t 2 0 − t 2 x )/(2t 0 t x ) . The P and Θ naturally define a chiral symmetry C = iτ x ⊗ σ y ⊗ γ 0 . Then where D denotes a diagonal matrix. The winding number around a Dirac point is given by [34] ν = ∇ k Det[A(k)]dk.
When θ x = θ 0 and θ y = 0 and π/2, there are four separate Dirac points with non-zero winding numbers shown in the main text. Due to the anti-commutation {Θ, U C } = 0, we have N T (−k) = −T † N (k)T ⇒ det A(k) = det A(−k).
As a result, the winding numbers around two opposite-momentum Dirac points are the same.