Boundary-obstructed topological high-T$_c$ superconductivity in iron pnictides
Xianxin Wu Wladimir A. Benalcazar Yinxiang Li Ronny Thomale Chao-Xing Liu Jiangping Hu
cond-mat.supr-con cond-mat.mes-hall cond-mat.str-el
















 






c_k,σ^†    c_k,σ
ç_̧ķ↑̧   c_-k↓
c_k↑^†   c_-k↓^†
c_p↑ c_-p+q↓
c_p↑^† c_-p+q↓^†

#1⟨#1⟩



    1/2Phys. Rev. A Phys. Rev. B Phys. Rev. D Rev. Mod. Phys.  Phys. Rev. Lett. 
        addendum
           
           
           Acknowledgements
               
               
               
               
               
               
        
        
        
        
        
        
        
        
        
        Department of Physics, the Pennsylvania State University, University Park, PA, 16802, USABeijing National Laboratory for Condensed Matter Physics,
        and Institute of Physics, Chinese Academy of Sciences, Beijing 100190, Chinaxianxinwu@gmail.comDepartment of Physics, the Pennsylvania State University, University Park, PA, 16802, USATin Ka-Ping College of Science, University of Shanghai for Science and Technology, Shanghai, 200093, ChinaInstitut für Theoretische Physik und Astrophysik,
          Julius-Maximilians-Universität Würzburg, 97074 Würzburg,
          Germanycxl56@psu.eduDepartment of Physics, the Pennsylvania State University, University Park, PA, 16802, USAjphu@iphy.ac.cnBeijing National Laboratory for Condensed Matter Physics,
        and Institute of Physics, Chinese Academy of Sciences, Beijing 100190, ChinaCAS Center of Excellence in Topological Quantum Computation and Kavli Institute of Theoretical Sciences,
        	University of Chinese Academy of Sciences, Beijing 100190, China
        
        Non-trivial topology and unconventional pairing are two central guiding principles in the contemporary search for and analysis of superconducting materials and heterostructure compounds. Previously, a topological superconductor has been predominantly conceived to result from a topologically non-trivial band subject to intrinsic or external superconducting proximity effect. Here, we propose a new class of topological superconductors which are uniquely induced by unconventional pairing. They exhibit a boundary-obstructed higher-order topological character and, depending on their dimensionality,  feature unprecedently robust Majorana bound states or hinge modes protected by chiral symmetry. We predict the 112-family of iron pnictides, such as Ca_1-xLa_xFeAs_2, to be a highly suited material candidate for our proposal, which can be tested by edge spectroscopy. Because of the boundary-obstruction, the topologically non-trivial feature of the 112 pnictides does not reveal itself for a bulk-only torus band analysis without boundaries, and as such had evaded previous investigations. Our proposal not only opens a new arena for highly stable Majorana modes in high-temperature superconductors, but also provides the smoking gun evidence for extended s-wave order in the iron pnictides.
        
        Boundary-obstructed topological high-T_c superconductivity in iron pnictides
            Jiangping Hu
            December 30, 2023
        ============================================================================
        
        
        
        
        
        Iron-based high temperature (high T_c) superconductors have recently appeared as an exciting platform to realize topological superconductivity at high temperatures<cit.>. Due to an intrinsic superconducting proximity effect, the surfaces in these materials can host Majorana zero modes (MZMs), evidence of which has been observed in the vortex core of Fe(Te,Se) and (Li_1-xFe_x)OHFeSe crystals<cit.>. That is, a band inversion between the Fe d and the ligand p orbitals is found, which then culminates with the superconducting proximity imposed by the particle-particle instability at the Fermi level. While this band inversion appears to be rather generic for the pnictides, current experimental evidence suggests that the topological superconducting phase necessitates significant tuning of Fermi level and chemical composition.
        
        
        In unconventional high T_c superconductors, the pairing symmetry is as essential as it can be difficult to directly identify it. In cuprates, only several years after their discovery, a d-wave pairing symmetry has been unambiguously proved by detecting the π phase shift  in corner junction interferometer experiments<cit.>. In iron-based superconductors, the pairing symmetry has been subject of a long-lasting debate ever since their discovery a decade ago. An s_±-wave pairing, possessing a sign-reversed gap on hole and electron pockets in momentums space,  has been  proposed for iron-based superconductors with some indirect evidences in neutron scattering and scanning tunneling microscopy <cit.>. So far, however, no decisive experiment has been proposed to distinguish s_±-wave from a sign-unchanged s-wave pairing because both states share the same A_1g symmetry character. The s_±-wave pairing was suggested to realize DIII topological superconductors<cit.>. As we will show in this work, s_±-wave and s-wave, despite their identical symmetry character, can give rise to different topological phases, and thus, we claim that the topological aspects in iron-based high T_c superconductors can shine a light on this outstanding problem<cit.>.
        
        Higher-order topological phases<cit.> are a new family of phases of matter with the defining property of hosting fractional charges or topological states at corners or hinges of the material. In 2D insulators, these phases manifest fractional corner charges protected by crystalline symmetries and are directly related to the positions of the Wannier centers of the occupied bands<cit.> or to the topology of its Wannier bands <cit.>. In superconductors, despite the absence of a Wannier description, particle-hole and/or chiral symmetries can also protect the existence of corner-localized MZMs. Recently, several proposals have been put forward for the realization of higher-order topological superconductors (HOTSC) in 2D or 3D <cit.>,
        some of which include certain iron-based superconductor compounds<cit.>.
        Setting aside that the complete characetrization of HOTSC by topological invariants is still missing, there is, in addition, a significant disconnect between the toy models studied in those works and actual material candidates.
        
        In order to understand the new class of superconductors we propose in this work, we develop a novel fused perspective on the current fields of unconventional superconductivity and higher-order topological states of matter. Specifically, we show that s_±-wave symmetry pairing, together with the topological properties in the 112-family of iron pnictides, drives the material into a HOTSC that hosts a Kramers pair of MZMs at each corner of each unit-layer. Different from previous proposals, the MZMs are particularly robust, as they do not depend on crystal symmetries, but are protected by chiral symmetry.  The demonstration of high-order topology in this family of compounds will provide a "smoking-gun" evidence for s_±-wave pairing. Remarkably, the topological phase we find for this material has the property of being adiabatically connected to the trivial phase in the absence of boundaries, but, once boundaries are introduced, an edge-localized obstruction topologically separates it from the trivial phase. This property, referred to as boundary topological obstruction, was originally identified in the minimal model for a higher-order topological insulator (HOTI) hosting a quantized quadrupole moment protected by reflection symmetries <cit.>. In insulators, boundary topological obstructions have recently been explained in terms of the Wannier centers of the occupied bands <cit.>. Here, we generalize the concept of boundary topological obstructions to superconductors – for which a Wannier description is absent – and identify the 112-family of iron pnictides as the first intrinsic material realization of boundary-obstructed HOTSC. The existence of MZMs in boundary-obstructed HOTSCs is a clear signature of its nontrivial topology. Upon a phase transition into a trivial phase, the bulk remains gapped, only the edges become gapless, providing one-dimensional channels for the MZMs to hybridize as they disappear into the trivial phase.
        
        
        The detection of MZMs in this material would be a decisive evidence for s_±-wave paring in iron-based superconductors for the following reason: The 112-family of iron pnictides, including Ca_1-xLa_xFeAs_2<cit.> and (Ca,Pr)FeAs_2<cit.>, with T_c up to 47 K<cit.>, are intrinsic topological insulator/high T_c superconductor heterostructures<cit.> with a staggered intercalation between zigzag chainlike As1 layers with the quantum spin Hall state and the superconducting Fe-As layers along the c axis, as shown in Fig.<ref>(a). The edge Dirac cones from the As1 layers at two orthogonal (100) and (010) edges are in proximity to projections of bulk pockets around Γ and M from adjacent FeAs layers, respectively. The s_±-wave pairing with opposite gap functions on pockets around Γ and M points will create the Majorana Kramers pairs at corners, as demonstrated in Fig.<ref>(c).
        
        In what follows, we first investigate s_± pairing in Ca_1-xLa_xFeAs_2 and relate the bulk spectrum to that of its edges, which gives rise to the mechanism that realizes the corner MZMs. We then demonstrate that under s_± pairing, the topological phase in this material is boundary-obstructed by showing that across the topological phase transition (TPT) between the topological and trivial phases, the energy gap only closes at the boundary of a slab configuration, while it remains open in the bulk. Finally, since the two phases are protected by chiral symmetry, we propose a new quantity, the edge winding number, as the invariant that captures the corresponding topological obstruction, and show that this invariant jumps by an integer across a TPT.
        
        s_± pairing in Ca_1-xLa_xFeAs_2– We take CaFeAs_2 as a typical example in 112-type. Besides the hole pockets around Γ and electron pockets around M from FeAs layers, there are additional Fermi surfaces in CaFeAs_2 attributed to the zig-zag As1 layers (Fig.<ref>(c)). As the correlation effect is relatively weak in As atoms, the pairing state of CaFeAs_2 is expected to be dominantly determined by the Fermi surfaces from FeAs layers. Here we adopt a five-band model whose band structure fits well with those in DFT calculations [see Sec. I in the supplementary material (SM)]. Fig.<ref>(a) displays the spin susceptibility. The peak around (π,0) is attributed to the Fermi surface nesting between the hole and electron pockets. Consequently, the dominant pairing is s_±-wave from the effective repulsive electron-electron interactions mediated by spin fluctuations<cit.>. Fig.<ref>(b) shows the typical gap function of s_± pairing from random phase approximation (RPA) calculations for CaFeAs_2, revealing a sign change in superconducting gaps between hole and electron pockets (see Sec. II in SM). The sign change and gap size can also be well described by a simple form factor cos k_x×cos k_y in one-Fe unit cell originating from the next-nearest neighbor antiferromagnetic exchange coupling<cit.>. Unlike the d-wave state in cuprates, in which a sign change occurs in orthogonal directions, s_± pairing with a sign change in momentum space is extremely difficult to be detected via Josephson interferometry. Although some evidence for the existence of sign change in pairing states have been provided in  inelastic neutron scattering and scanning tunneling microscopy measurements<cit.>, a decisive signature for s_±-pairing is still missing in iron pnictides. In the following, we show that s_± pairing gives rise to a HOTSC phase in CaFeAs_2 with Majorana corner states, which can be regarded as the direct evidence for s_±-pairing.
        
        Effective Model for As1 layers – We start with the tight-binding model for the As1 layers. A two-dimensional four-band model on a square lattice (Fig.<ref>(b)) has been derived to capture the band structure attributed to p_x and p_y orbitals of two As1 atoms in one unit cell.  We introduce the operator ψ^†_kσ=(c^†_Axσ(k),c^†_Ayσ(k),c^†_Bxσ(k),c^†_Byσ(k)), where c^†_αησ(k) is a Fermionic creation operator with  σ,  η and α being spin, orbital and sublattice indices respectively. The tight-binding Hamiltonian reads
        
            ℋ_0=∑_kσψ^†_kσh(k)ψ_kσ.
        
        The matrix elements in the Hamiltonian h(k) are provided in Sec. III in the SM. A band inversion occurs at the X point and it generates two gapless Dirac cones in the bulk dispersion along X-M path without spin-orbit coupling (SOC), protected by the screw axis along y<cit.>. Two small pockets appear around the X point deriving from these cones, demonstrated by orange circles in Fig.<ref>(c) and supported by ARPES experiments<cit.>. Further including SOC opens a gap in the Dirac cones and the As1 layers becomes ℤ_2 topologically nontrivial, leading to an intrinsic topological insulator/high T_c superconductor heterostructure in CaFeAs_2. Around the X point, the bulk dispersion of As1 layers can be described by an effective Hamiltonian ℋ^ X_eff=∑_kψ̃^†_kh_eff(𝐤)ψ̃_k with basis ψ̃^†_k=(c^†_Xk,-↑,c^†_Xk,+↑,c^†_Xk,-↓,c^†_Xk,+↓), where "+/-" denotes the eigenvalue of C_2z for eigenstates at the X point and
        
            h_eff(k)=ϵ_0(𝐤) + M(𝐤)σ_z-A_1k_xs_0σ_2 + A_2k_ys_3σ_1.
        
        Here ϵ_0(𝐤)=C_0+C_1k^2_x+C_2k^2_y and M(k)=M_0-B_1k^2_x-B_2k^2_y, and σ and s are Pauli matrices in orbital and spin degrees of freedom. The parameters in the model are given in Sec. IV of SM. The edge states can be obtained by solving a semi-infinite system along x or y direction within tight-binding model. For the (100) edge, a distorted Dirac cone appears around the Γ̅ point located above the Fermi level, as shown in Fig.<ref>(c), while, for the (010) edge, a Dirac cone occurs around the X̅ point and is embedded in the bulk conduction bands, as shown in Fig.<ref>(d). We emphasize that the Fermi level for both edges only crosses the lower part of the edge Dirac cone within a relatively large electron doping region from the substitution of La/Pr for Ca in experiments.
        
         Below T_c, the As1 layers become superconducting through the proximity effect to the adjacent FeAs layers. We model the superconducting pairing on As1 layers the same way as in FeAs layers and consider a spin-singlet intra-orbital pairing within the same sublattice. The corresponding pairing Hamiltonian reads
         
            ℋ_SC   =   ∑_ανσkσ[Δ_0+2Δ_1(cosk_x+cosk_y)]
                  c^†_ανσ(k)c^†_ανσ̅(-k)
            +h.c.,
        
        where Δ_0 and Δ_1 are the on-site pairing and pairing between the next nearest neighbor sites, respectively. The Δ_1 term gives rise to the well-known s_±-wave pairing in iron based superconductors. Owing to the absence of spin-flip SOC terms, the Bogoliubov-de Gennes (BdG) tight-binding Hamiltonian ℋ_BdG=ℋ_0+ℋ_SC can be written as two block diagonal parts ℋ_BdG^↑↓ and ℋ_BdG^↓↑ (see Sec. V in SM). In each block, time reversal and particle-hole symmetries are absent but chiral symmetry is preserved, thus, each block belongs to class AIII in the ten fold classification.
        
         Boundary-obstruction and phase transitions – To study the edge properties of the As1 layers, we consider a slab configuration for the above BdG Hamiltonian ℋ^↑↓_BdG with an open boundary along either x or y directions. With the on-site pairing term Δ_0, the edge states from As1 layers open a gap, as shown in Fig.<ref> (a1) and (b1). Now we investigate the effect of s_±-wave pairing on these edge states. At the (100) edge, the gap of edge states around Γ̅ monotonously increases with increasing Δ_1, as shown in Fig.<ref> (a1) to (a3). For the (010) edge, however, the gap around X̅ exhibits a rather different behavior. Upon Δ_1 increasing to 0.051 eV, the gap closes and a pair of gapless modes with linear dispersion appears [Fig. <ref> (b2)]. Further increasing Δ_1 reopens the gap, suggesting a TPT that separates a Δ_0-dominated phase from a Δ_1-dominated phase. Note that throughout this transition there is no gap closing in the bulk states (see Sec. VIII in SM). In fact, we show in Sec. VII of the SM that all the symmetry indicator invariants due to the C_2z and reflection symmetries in the lattice structure of ℋ_BdG identically vanish due to time-reversal symmetry, which is a necessary condition for the existence of boundary-obstructed phases.
        
        To calculate the effective pairing at the edges, we first analytically obtain the wavefunctions of the edge states by solving ℋ_eff with open boundary conditions and then projecting the bulk pairing on the edge states<cit.>. The obtained pairings at the Dirac points on the (100) and (010) edges are
        
            Δ^(100)_eff   =   Δ_0+2Δ_1M_0/B_1,
            Δ^(010)_eff   =   Δ_0-2Δ_1M_0/B_2.
        
         We find that Δ_0 provides the same pairing at the two edges while the s_±-wave pairing Δ_1 provides an opposite pairing sign. This can be heuristically understood from Fig.<ref>(c). The edge Dirac states appear around Γ̅ and X̅ for (100) and (010) edges; the momentum independent term proportional to Δ_0 induces the same positive superconducting gap at both edges; the s_±-wave pairing, on the other hand, induces superconducting gaps with opposite signs as the corresponding Dirac cones are in proximity to the bulk superconducting gap around Γ and M at the (100) and (010) edges, respectively (see Fig.<ref>(c)). At a fixed non-zero value of Δ_0, with increasing Δ_1, Δ^(100)_eff monotonously increases while Δ^(010)_eff first decreases to zero, followed by increasing in amplitude albeit with opposite sign, which is consistent with the aforementioned numerical calculations. The finite chemical potential at the edges with respect to Dirac points should be taken into consideration but their relative small values have a negligible effect on the effective pairing (see Sec. VI in SM).
        
        To characterize the topological nature of the TPT, let us focus on the blocks ℋ_BdG^↑↓, and ℋ_BdG^↓↑, each of which belongs to class AIII. Hamiltonians in this class are topologically characterized by the 1D winding number<cit.>, defined by ν_1=i/2π∫_BZdk Tr[q^-1_k∂_kq_k]. Here the unitary q_k matrix is the off-diagonal part of the so-called Q matrix, given by
        Q_k=1-2P_k=(
        [     0   q_k; q^†_k     0 ]) within the eigenbasis of chiral symmetry, where P_k is the projection operator of the Hamiltonian for a slab model ℋ^↑↓_BdG,slab with N lattice sites. We first consider the winding number on a slab configuration with open boundary along the (010) direction. The total winding number ν_1 is zero across the TPT (see Sec. VIII in SM).  Motivated by the fact that the bands close at the (010) edges of the slab during the TPT, we define a site-resolved winding number by projecting the total winding number ν_1 into the lattice site basis as
        
            ν^i_1=i/2π∑_γ∫_BZdk[q^-1_k∂_kq_k]_iγ,iγ,
        
        such that ν_1=∑_i=1^Nν^i_1. Here, i is the index for lattice site and γ denotes the sublattice or orbital index (details are given Sec. VIII in SM). This site-resolved winding number resembles the site-resolved polarization defined in Ref. benalcazar2017quadPRB to calculate the edge-localized dipole moments in insulators, including the boundary-obstructed quadrupole topological insulator <cit.>. The dependence of ν^i_1 on the site i for the HOTSC and trivial phases is displayed by the blue and red curves in Fig.<ref>(a). For both curves, one can see that the contribution to the winding number mainly comes from edges, with the opposite edges having opposite contributions. We also notice that the profile of ν^i_1 near one edge (see the zoom-in Fig.<ref>(a)) has a substantial difference between the HOTSC and trivial phases. We define the edge winding number ν^T_1 and ν^B_1 for the top and bottom edges, respectively, by summing ν^i_1 for the upper half part (i=1,...,N/2) and lower half part (i=N/2+1,...,N) of the slab and examine the winding number change |Δν^T/B_1| between two phases for the upper and lower half parts as a function of the lattice size N. As shown in Fig.<ref>(b), we find that the winding number change |Δν^T_1| approaches 1 in the thermodynamic limit (N→∞). Thus, the edge winding number ν^T/B_1 characterizes the different topologies across the TPT in our new class of superconductor.
        
        Majorana corner states – In the Δ_1 dominated phase, the (100) and (010) edges have the opposite superconducting gaps, belonging to topologically distinct phases. As a consequence, MZMs are expected to occur at the corners where they meet. We adopted the Hamiltonian ℋ_BdG with both open boundary conditions along x and y directions, and performed calculations with several lattice sizes. There are eight mid-gap states and their energies are zero up to finite-size effects. We examine the energy splitting due to the hybridization of the MZMs at different corners by plotting the energies E_M of mid-gaps states as a function of lattice size in Fig.<ref>(c). The linear relationship between log_10(E_m) and lattice size suggests that E_M goes to zero in the thermodynamic limit, compatible with the existence of zero-energy states exponentially localized at the corners of the lattice, as shown in Fig.<ref>(d). At each corner, there are two zero-energy states forming a Majorana Kramers pair. The appearance of MZMs at the corners supports our prediction of a HOTSC phase with boundary-obstruction in CaFeAs_2.
        
        Discussion –
        As the 112-family of iron pnictides is 3D and the interlayer coupling is relatively weak, the helical Majorana states, localized at the hinges between (100) and (010) surfaces, have a weak dispersion along k_z. MZMs can also appear at corners/hinges between non-orthogonal surfaces, as long as the corresponding surface states have the opposite effective superconducting gap in proximity to the bulk s_± pairing. Therefore, the s_± pairing is directly manifested by the appearance of MZMs at hinges, protected by chiral symmetry and irrespective of crystalline symmetries. We contrast our proposal with that of Ref. ZhangF2013. In it, a heteroustructure couples Rashba semiconductors with nodeless iron-based superconductors. Unlike our proposal, the realization of such a structure would require a large Rashba SOC. It is also worth noting that the toy model in Ref. WangQY2018, realizes a bulk-and boundary-obstructed HOTSC with and without the fourfold rotational symmetry, respectively. The 112-family of iron pnictides to some extent can be considered as a realization of the toy model in Ref. WangQY2018 with boundary-obstructed HOTSC. The MZMs we find are robust against impurities and disorders in real material scenarios that always break crystalline symmetries. In realistic materials, the symmetries are reduced compared with the adopted model. Under this scenario, SOC terms could in principle appear that flip the spin sectors, however, the Majorana Kramers pairs will survive as they are protected by time reversal symmetry. The detection of hinge MZMs in CaFeAs_2 hence provides a “smoking-gun" evidence for the s_± pairing in iron pnictides.  As the Fermi level only crosses the lower part of edge Dirac cones for edges in As1 layers within a wide electron-doping range, the MZMs can survive in electron-doped compounds Ca_1-xLa_xFeAs_2, providing a new high-temperature platform for MZMs without fine-tuning.
        
        In particular, due to the weak interlayer coupling in 112-family of iron pnictides, the cleavage occurs between Ca/La layers and As1 layers or FeAs layers, generating step edges on (001) surfaces. At the ends of step edges, which can be viewed as corners of As1 layers, MZMs can give rise to a zero-bias peak in the  scanning tunneling microscopy measurements (see Sec. VIII in SM) and their localized nature can be manifested in spatial zero-energy mapping (see Fig.<ref>(d)). Hinge MZMs are expected to induce a zero-bias anomaly in transport measurements with junctions in contact with the corresponding hinge<cit.>. The Majorana Kramers pair can split under an in-plane magnetic field. As the Sb doping can enhance both T_c and the quantum spin Hall gap of the As1 layer<cit.>, Ca_1-xLa_xFe(As_1-ySb_y)_2 can be a good choice for samples in experimental measurements.
        
        Conclusion – We propose the 112 family of iron pnicitides as the first material realization of boundary-obstructed topological superconductivity, owing to their intrinsic s_± pairing and effective topological insulator/high-T_c superconductor heterostructure profile. The edge topological obstruction, indepedent of crystalline symmetries and uniquely characterized by the edge winding number, provides a robust platform for the realization of MZMs, which en passant also constitutes decisive evidence for  s_±-wave pairing in the 112 pnictides.
        
        
        
         Acknowledgements:
        We thank  J. Yu and N. Hao for helpful discussions. W.A.B. thanks the support of the Eberly Postdoctoral Fellowship at the Pennsylvania State University. C.X.Liu acknowledges the support of the Office of Naval Research (Grant No. N00014-18-1-2793), DOE grant (DE-SC0019064) and Kaufman New Initiative research grant of the Pittsburgh Foundation.  J. Hu is supported by the Ministry of Science and Technology of China 973 program (Grant No. 2017YFA0303100), National Science Foundation of China (Grant No. NSFC-11888101), and the Strategic Priority Research Program of CAS (Grant No. XDB28000000). The work in Würzburg is funded by the Deutsche
        Forschungsgemeinschaft (DFG, German Research Foundation) through
        Project-ID 258499086 - SFB 1170 and through the Würzburg-Dresden Cluster of Excellence on Complexity and Topology in Quantum Matter –ct.qmat Project-ID 39085490 - EXC 2147.
        
        
        
        
        
        
        
        
        
        
        
        § BAND STRUCTURE AND TIGHT-BINDING MODEL FOR FEAS LAYERS IN CAFEAS_2
        
        We adopted to a 2D five-band model capture the low-energy electronic structure of CaFeAs_2. The hopping parameters are given in Table <ref> and the obtained band structure is displayed in Fig.<ref>(a), which fit well with DFT. In
        general, a ten-band model should be used to describe the low-energy
        physics of CaFeAs_2 due to the absence of the glide mirror symmetry and four-fold rotational symmetry. But this five-band model can already provide a good description for the Fermi surfaces and therefore we adopt it in the following RPA calculations. It is noted that the largest hole pockets around Γ point are attributed to Ca d and As1 p_z orbital and the small pockets around X point are attributed to As1 p_x,y orbitals<cit.>. As the correlation effect in As and Ca atoms is relatively weak, we neglect the pockets from them in our model and focus on the Fermi surfaces derived from the FeAs layers.
        
        
        
        
        § RPA CALCULATIONS
        
        In this section, we explain the formalism of the multiorbital RPA approach<cit.>, adopted in the main text. The adopted onsite Coulomb interaction terms in CaFeAs_2 are,
        
            H_int=U∑_i,αn_iα↑n_iα↓+U'∑_i,α<βn_iαn_iβ    +J_H∑_i,α<β,σσ^'c^†_iασc^†_iβσ^'c_iασ^'c_iβσ
                     +J'∑_i,α≠βc^†_iα↑c^†_iα↓c_iβ↓c_iβ↑,
        
        where n_i,α=n_i,α,↑+n_i,α,↓. U, U', J and J' represent the intra- and inter-orbital repulsion, the Hund's rule and pair-hopping terms. In the following calculations, we use Kanamori relations U=U'+2J and J=J' as requried by the lattice symmetry.
        
        The multi-orbital susceptibility is defined as,
        
            χ_l_1l_2l_3l_4(q,τ)   =   1/N∑_kk'⟨ T_τ c^†_l_3σ(k+q,τ)c_l_4σ(k,τ)c^†_l_2σ'(k'-q,0)c_l_1σ'(k',0) ⟩ .
        
        In momentum-frequency space, the multi-orbital bare susceptibility is given by
        
            χ^0_l_1l_2l_3l_4(q,iω_n)=-1/N∑_kμνa^l_4_μ(k)a^l_2*_μ(k)a^l_1_ν(k+q)a^l_3*_ν(k+q)n_F(E_μ(k))-n_F(E_ν(k+q))/iω_n+E_μ(k)-E_ν(k+q),
        
        where μ and ν are the band indices, n_F is the usual Fermi distribution, l_i(i=1,2,3,4) are the orbital indices, a^l_i_μ(k) is the l_i orbital component of the eigenvector for band μ resulting from the diagonalization of the tight-binding Hamiltonian H_0 and E_μ(k) is the corresponding eigenvalue. With interactions, the RPA spin and charge susceptibilities are given by
        
            χ^RPA_s(q)=χ^0(q)[1-U̅^sχ^0(q)]^-1,
            χ^RPA_c(q)=χ^0(q)[1+U̅^cχ^0(q)]^-1,
        
        where U̅^s (U̅^c) is the spin (charge) interaction matrix,
        
            U̅^s_l_1l_2l_3l_4(q)   =   
            U        l_1=l_2=l_3=l_4,
            
            U'         l_1=l_3≠ l_2=l_4,
            
            J       l_1=l_2≠ l_3=l_4,
            
            J'       l_1=l_4≠ l_2=l_3,
            
            U̅^c_l_1l_2l_3l_4(q)   =   
            U       l_1=l_2=l_3=l_4,
            
            -U'+2J        l_1=l_3≠ l_2=l_4,
            
            2U'-J        l_1=l_2≠ l_3=l_4,
            
            J'       l_1=l_4≠ l_2=l_3,
        
        In the main text, we plot the largest eigenvalues of the susceptibility matrix χ^RPA_s,l_1l_1l_2l_2(q,0). Within RPA approximation, the effective Cooper scattering interaction on Fermi surfaces is,
        
            Γ_ij(k,k')=   ∑_l_1l_2l_3l_4a^l_2,∗_v_i(k)a^l_3,∗_v_i(-k)Re[Γ_l_1l_2l_3l_4(k,k',ω=0)]a^l_1_v_j(k')a^l_4_v_j(-k'),
        
        where the momenta k and k' is restricted to different FSs with k∈ C_i and k'∈ C_j. The  orbital vertex function Γ_l_1l_2l_3l_4 in spin singlet and triplet channels <cit.> are
        
            Γ^S_l_1l_2l_3l_4(k,k',ω)=   [3/2U̅^sχ^RPA_s(k-k',ω)U̅^s+1/2U̅^s-1/2U̅^cχ^RPA_c(k-k',ω)U̅^c+1/2U̅^c]_l_1l_2l_3l_4,
            Γ^T_l_1l_2l_3l_4(k,k',ω)=   [-1/2U̅^sχ^RPA_s(k-k',ω)U̅^s+1/2U̅^s-1/2U̅^cχ^RPA_c(k-k',ω)U̅^c+1/2U̅^c]_l_1l_2l_3l_4,
        
        where χ^RPA_s and χ^RPA_c are the RPA spin and charge susceptibility, respectively. The pairing strength functional for a specific pairing state is given by,
        
            λ[g(k)]=-∑_ij∮_C_idk_/v_F(k)∮_C_jdk'_/v_F(k')g(k)Γ_ij(k,k')g(k')/(2π)^2∑_i∮_C_idk_/v_F(k)[g(k)]^2,
        
        where v_F(k)=|∇_kE_i(k)| is the Fermi velocity on a given Fermi surface sheet C_i. The pairing vertex function in spin singlet and triplet channels are symmetric and antisymmetric parts of the interaction, that is, Γ^S/T_ij(k,k')=1/2[Γ_ij(k,k')±Γ_ij(k,-k')]. Fig.<ref>(b) displays the pairing eigen value as a function of interaction parameters for CaFeAs_2 in RPA calculations with J/U=0.1. We find that the s_±-wave pairing is always dominant, which is attributed to the Fermi surfaces nesting between electron and hole pockets.
        
        
        
        
        § TIGHT-BINDING MODEL FOR AS1 LAYERS
        
         The crystal structure of CaFeAs_2 is shown Fig.1(a) in the main text and the As layers are shown in Fig.1(b).  A two-dimensional four-band model has been derived  to capture  the band structure attributed to p_x and p_y orbitals of two As-1 in a unit cell. As a unit cell includes  two As1 atoms, we divide the As1 lattices  into two sublattices.  We introduce the operator ψ^†_kσ=(c^†_Axσ(k),c^†_Ayσ(k),c^†_Bxσ(k),c^†_Byσ(k)), where c^†_αησ(k) is a Fermionic creation operator with  σ,  η and α being spin, orbital and sublattice indices respectively. The tight-binding Hamiltonian can be written as:
         
            ℋ_TB=∑_kσψ^†_kσh(k)ψ_kσ.
        
          The matrix elements in the Hermitian h(k) matrix are given by
        
            h_11=h_33=ϵ_X+2t^11_1cosk_x +2t^11_2cosk_y, 
                  h_13=(2t^13_1e^i(x_0-1)k_x+2t^13_2e^ix_0k_x)cos(k_y/2),
                  h_14=h_23=-(2it^14_1e^i(x_0-1)k_x+2it^14_2e^ix_0k_x)sin(k_y/2),
                  h_22=h_44=ϵ_Y+2t^11_2cosk_x +2t^11_1cosk_y,
            
                  h_24=(2t^24_1e^i(x_0-1)k_x+2t^24_2e^ix_0k_x)cos(k_y/2),
        
        with x_0=0.418 being the difference between the x components of B and A sublattices.  The corresponding tight binding parameters are specified in unit of eV as,
        
            ϵ_X=-0.30,   ϵ_Y=-0.109,    t^11_1=-0.149,  
                  t^11_2=0.128,    t^13_1=0.89,    t^13_2=0.649, 
                  t^14_1=1.169,   t^14_2=-1.740,    t^24_1=0.567  
                   t^24_2=1.213.
        
        Now we consider the atomic SOC term in the As1 atoms and the SOC Hamiltonian can be written as
        
            H_so   =   iλ/2∑_ασkσ c^†_α xσ(k) c_α yσ(k)+h.c.
               =   iλ/2∑_αk[ c^†_α x↑(k) c_α y↑(k)-c^†_α x↓(k) c_α y↓(k)]+h.c.
               =   iλ/2∑_k(ψ^†_k↑h_so,↑ψ_k↑+ψ^†_k↓h_so,↓ψ_k↓)
        
        where the spin-orbit coupling strength for As atoms is λ=0.19 eV and h_so,↑=-h_so,↓=τ_0⊗-iσ_2. τ, σ, s are the Pauli matrices in sublattice, orbital and spin space. The total Hamiltonian is H_0=H_TB+H_so. Due to the conservation of s_z, this Hamiltonian can be block diagonal H_0=H_0↑+H_0↓, with h_0↑/↓(k)=h(k)+iλ/2 h_so,↑/↓. The obtained band structure with SOC is displayed in Fig.<ref>(a), in agreement with DFT calculations in Fig.<ref>(b)<cit.>.
        
        
        
        § EFFECTIVE MODEL AROUND X POINT FOR AS1 LAYERS IN CAFEAS_2
        
        For a point group operation, the wavefunction at k transforms as,
        
            P̂_gψ_k,α(r)   =   ∑_βψ_gk,β(r)D_βα[g],
        
        where D[g] is the representation matrix for the point group operation g.
        If we consider gk=G+k, then ψ_gk,β(r)=V^†_βα(G)ψ_k,α(r)=e^-iG·τ_βψ_k,β(r). Therefore, we have,
        
            P̂_gψ_k,α(r)   =   ∑_βψ_k,β(r)e^-iG·τ_βD_βα[g].
        
        For As1 lattice, the point group at X is C_2v. Around X point, the effective Hamiltonian can be written as ℋ^X_eff=∑_kψ̃^†_kh_eff(𝐤)ψ̃_k with ψ̃^†_k=(c^†_Xk,-↑,c^†_Xk,+↑,c^†_Xk,-↓,c^†_Xk,+↓), where "+/-" denotes the eigenvalue of C_2z for eigenstates at the X point and k is relative to X(π,0). The representation matrices for the two generators in C_2v are: D(C_2z)=is_3σ_3 and D(M_xz)=is_2. Here σ and s are Pauli matrices in orbital and spin spaces. The time reversal symmetry operator is T̂=-is_2𝒦. The corresponding Hamiltonian matrix reads,
        
            h_eff(𝐤)     =   ϵ_0(𝐤)+ ([            M(𝐤)  iA_1k_x+A_2k_y               0               0; -iA_1k_x+A_2k_y           -M(𝐤)               0               0;               0               0            M(𝐤)  iA_1k_x-A_2k_y;               0               0 -iA_1k_x-A_2k_y           -M(𝐤);                 ]),
               =   ϵ_0(𝐤) + M(𝐤)σ_z-A_1k_xs_0σ_2 + A_2k_ys_3σ_1
        
        where ϵ_k=C_0+C_1k^2_x+C_2k^2_y and M(k)=M_0-B_1k^2_x-B_2k^2_y. This Hamiltonian satisfies all the symmetry operations: ĝH(k)ĝ^-1=H(gk) with ĝ=C_2z,M_xz,T̂. Fitting the above effective model to the TB band structure around X point, the resulting band structure is displayed in Fig.<ref> and the corresponding parameters are given in Table <ref>. These bands show good agreement around X point.
        
        
        
        § BDG HAMILTONIAN FOR BULK STATES
        
        As the As1 atoms are located just above the Fe site, we consider the s_± pairing for As1 layer (intrasublattice pairing) from the proximity effect from the adjacent FeAs layers. In the superconducting phase, the BdG Hamiltonian, for the basis Ψ^†_k=(ψ^†_k↑,ψ^T_-k↓, ψ^†_k↓,ψ^T_-k↑), is given by
        
            ℋ_BdG=1/2∑_kΨ^†_kh_BdG(𝐤)Ψ_k, 
                  h_BdG=([ h_B,↑↓(k)         0;         0 h_B,↓↑(k);           ])
                        =([     h_0↑(k)-μ       Δ_↑↓(k)             0             0;     Δ_↑↓^†(k) -h^*_0↓(-k)+μ             0             0;             0             0     h_0↓(k)-μ       Δ_↓↑(k);             0             0     Δ_↓↑^†(k) -h^*_0↑(-k)+μ;               ]).
        
        Here we consider the intrasublattice and intraorbital singlet pairing and Δ_↑↓(k)=-Δ_↓↑(k)=[Δ_0+2Δ_1(cosk_x+cosk_y)]τ_0σ_0. h_B,↑↓ and h_B,↓↑ are related by time reversal symmetry or particle hole symmetry. However, for h_B,↑↓, time reversal symmetry and particle hole symmetry are broken but chiral symmetry is preserved. Therefore, it belongs to AIII class and the topological invariant is ℤ in 1D.
        
        Around X point, with basis Ψ̃^†_k=(ψ̃^†_k,ψ̃_-k) the BdG Hamiltonian can written as ℋ_BdG=1/2∑_kΨ̃^†_kh̃_BdG(𝐤)Ψ̃_k. h̃_B,↑↓(k) can be explicitly written as,
        
            h̃_B,↑↓(k)=([   ϵ_0(k)+M(𝐤)-μ  iA_1k_x+A_2k_y          Δ_X(k)               0; -iA_1k_x+A_2k_y   ϵ_0(k)-M(𝐤)-μ               0          Δ_X(k);          Δ_X(k)               0  -ϵ_0(k)-M(𝐤)+μ -iA_1k_x+A_2k_y;               0          Δ_X(k)  iA_1k_x+A_2k_y  -ϵ_0(k)+M(𝐤)+μ;                 ]),
        
         with Δ_X(k)=Δ_0+Δ_1(k^2_x-k^2_y).
        
         
        
        § EFFECTIVE PAIRING FOR THE EDGE STATES
        
         To calculate the effective pairing on edges of As1 layers, we first
        analytically obtain the wavefunctions of edge states by solving H_eff with open boundary conditions and then project the
        bulk pairing on edge states (omit the ϵ_0(𝐤) term in the effective model)<cit.>. The obtained pairings on (100) and (010) edge with finite chemical potential are,
        
            Δ^(100)_eff   =   Δ_0+2Δ_1M_0/B_1-2Δ_1(1+B_2/B_1)μ_1^2/A^2_2,
            Δ^(010)_eff   =   Δ_0-2Δ_1M_0/B_2+2Δ_1(1+B_1/B_2)μ_2^2/A^2_1,
        
        where μ_1/2 is the chemical potential on edge states and is defined relative to the corresponding Dirac points. According to the Fig.2(c) and (d) in the main text and Table II, the chemical potential is in the range from 0.1 to 0.15 eV and thus |μ_1/2/A_2/1|^2≪ 1. Therefore the finite chemical potential has a negligible effect on the effective pairing on edge states. This will not has a qualitative change on the analysis of pairing on edge states in the main text.
        
        
        
        
        
        
        
        § SYMMETRIES OF THE HAMILTONIAN AND ABSENCE OF BULK OBSTRUCTIONS
        
        
        Bulk-obstructed topological phases have gapped Hamiltonians adiabatically disconnected from each other. The adiabatic disconnection manifest in the spectrum by the fact that, when transitioning from one phase to another, the Hamiltonian goes through a gapless point to achieve a band inversion. In the presence of crystalline symmetries, such band inversions typically occur at high-symmetry points of the Brillouin zone and can be tracked by calculating symmetry indicators topological invariants. In fact, whenever a topological phase has nontrivial symmetry indicators, it will necessarily have to close its bulk bands when deformed into the trivial phase. Therefore, a boundary-obstructed phase necessarily has trivial symmetry indicators.
        
        In what follows, we show that the Hamiltonian for the As1 layers in Eq.<ref>
        
        
        
        
        
        
        
        
        
        
        has vanishing symmetry indicator invariants and thus it is not constrained to be bulk-obstructed.
        
        For completeness, let us us begin by pointing out the local symmetries. This Hamiltonian obeys particle-hole, time-reversal, and chiral symmetries, which are, respectively,
        
            Ξ̂ h_BdG() Ξ̂^-1   = - h_BdG(-),   Ξ̂ = (s_x ⊗ρ_x ⊗τ_0 ⊗σ_0) 𝒦
            Θ̂ h_BdG() Θ̂^-1   = h_BdG(-),   Θ̂=(is_y ⊗ρ_z ⊗τ_0 ⊗σ_0)𝒦
            Π̂ h_BdG() Π̂^-1   = -h_BdG(),   Π̂=s_z ⊗ρ_y ⊗τ_0 ⊗σ_0,
        
        where 𝒦 is complex conjugation. The degrees of freedom in which the Pauli matrices s, ρ, τ and σ operate are: spin, Nambu basis, sublattice, and orbital, respectively.
        
        
        
         §.§ Crystalline symmetries
        
        In addition to the local symmetries, our material has C_2z rotation and refrlection symmetries, as shown in Fig. <ref>. We now show that for both C_2z and reflection symmetries, the symmetry indicator invariants vanish due to time reversal symmetry for all phases.
        
        
        
          §.§.§ C_2 symmetry
        
        This Hamiltonian has additional C_2z rotational symmetry,
        
            r̂_2 h_BdG() r̂_2^-1 = h_BdG(-),   r̂_2=-is_z ⊗ρ_0 ⊗τ_x ⊗σ_0.
        
        The rotation operator obeys r̂_2^2=-1 and [r̂_2,Θ̂]=0.
        At the high symmetry points of the Brillouin zone (HSP), Π=Γ,  X, Y, M, we have Π=-Π, and h_BdG() commutes with both Θ̂ and r̂_2. In particular, let us calculate the C_2z eigenvalues of the energy eigenstates |u_n(Π)⟩,
        
            r̂_2 |u_n(Π)⟩ = r_p(Π) |u_n(Π)⟩.
        
        Here r_p(Π), for p=1 or 2 is the rotation eigenvalue at Π, and can be any of the two values  r_1(Π)=+i or r_2(Π)=-i. We then define the rotation symmetry indicator invariants
        
            [Π_p]=#r_p(Π)-#r_p(Γ)  ∈ℤ,
        
        where #r_p(Π) is the number of eigenstates in the negative energy bands with eigenvalue r_p(Π), for Π= X, Y, and M. Due to the fact that the number of negative energy bands is constant across the BZ, we have
        
            ∑_p=1,2#r_p(Π) = ∑_p=1,2#r_p(Γ) → [Π_1] + [Π_2]=0.
        
        Now, due to TRS, we have
        
            r̂_2 Θ̂|u_n()⟩ = Θ̂r̂_2 |u_n()⟩ = r^*_p Θ̂|u_n()⟩,
        
        which means that the rotation eigenvalues have to come in complex conjugate pairs. For the rotation invariants, Eq. <ref>, this implies that
        
            #r_1(Π)=#r_2(Π)    #r_1(Γ)=#r_2(Γ)
        
        or, in terms of the invariants,
        
            [Π_1] TRS= [Π_2].
        
        This, along with the restriction in Eq. <ref>, results in vanishing rotation invariants at all HSPs,
        
            [Π_p] TRS=0.
        
        Thus, we see that C_2 symmetry does not force the phases of h_BdG( k) to be bulk-obstructed.
        
        
        
          §.§.§ Reflection symmetry
        
        The Hamiltonian also has reflection symmetries,
        
            m̂_2 h_BdG(k_x,k_y) m̂_2^-1 = h_BdG(k_x,-k_y),   m̂_2= i s_x ⊗ρ_z ⊗τ_0 ⊗σ_z,
        
        The reflection operator obeys m̂_x^2=-1 and [m̂_x, Θ̂]=0. We can proceed as for C_2z symmetry to define the symmetry indicator invariants. However, the same constraints as in Eq. <ref> and Eq. <ref> apply, leading to trivial invariants at all HSPs.
        
        
        
        
        
        § PHASE TRANSITION ON THE EDGE AND WINDING NUMBER
        
        For a fixed Δ_0, with increasing Δ_1, there is no gap closing for the bulk state and (100) edge states but a phase transition for the (010) edge states, as displayed in Fig.<ref>. It indicates the resulting high-order topology is derived from the edge state not from bulk states, indicating that the system has boundary topological obstruction. In insulators, boundary topological obstructions have recently been explained in terms of the Wannier centers of the occupied bands. Following this, we calculate the corresponding Wannier centers and Fig.<ref> and <ref> show the Wannier center spectra as a function of k_y/k_x with a increasing Δ_1 and a fixed Δ_0=0.005 eV (the phase transition point is at Δ_1=0.051 eV). A small nonzero Δ_1 will immediately open a gap around k_x/k_y=0 and the Wannier center spectra donot exhibit significant changes across the phase transition on (010) edge states, suggesting that they cannot be utilized  to characterize the edge topological phase transition. Therefore, we need to introduce other topological invariants to identify the observed topological phase transition.
        
        
        
         §.§ winding number
        
        
        For a gapped system, we can define the Projection operator at k point as,
        
            P_k=∑^N_occ_n=1|u_n(k)⟩⟨ u_n(k)|,
        
        where |u_n(k)⟩ is the eigen states for BdG Hamiltonian and N_occ is the number of the occupied bands. It is convenient to introduce the Q matrix by,
        
            Q_k=1-2P_k,
        
        where Q matrix satisfies Q^†_k=Q_k and Q^2_k=1. The eigenvalues of Q_k are ±1. For systems with chiral symmetry, the Q matrix can be written as a block off-diagonal form in the basis of chiral operator,
        
            Q_k=(
            [     0   q_k; q^†_k     0 ]),
        
        where q_k satisfies q_kq^†_k=q^†_kq_k=1. Hence, the q matrix defines a map from the base space (Brillouin zone) to the space of unitary matrix U(N). The relevant homotopy group for projectors is π_d[U(N)]=𝒵 for odd d and N≥(d+1)/2<cit.>. The winding number in 1D and 3D are defined as,
        
            ν_1   =   i/2π∫_BZdkTr[q^-1∂_kq],
            ν_3   =   1/24π^2∫_BZd^3kϵ^ijlTr[(q^-1∂_k_iq)(q^-1∂_k_jq)(q^-1∂_k_lq)].
        
        In the weak pairing limit, q matrix can be written as,
        
            q_k=∑_ne^iθ_nk|n,k⟩⟨ n,k|,
        
        with e^iθ_nk=(ϵ_nk+iδ_nk)/|ϵ_nk+iδ_nk| and δ_nk=⟨ n,k|U_TΔ^†_k|n,k⟩ is the gap function<cit.> (the time-reversal symmetry T̂=U_T K). Here |n,k⟩ is the eigen state for normal Hamilton without superconductivity (not the BdG Hamiltonian). To the leading order, near the Fermi surface we have,
        
            e^iθ_nk=v_F(k_⊥-k_F)+iδ_nk_F/√(v^2_F(k_⊥-k_F)^2+δ^2_nk_F)
        
        In the limit δ_nk_F→ 0, we have θ_nk→ -π sgn(δ_kk)η(k_F-k_⊥) with η(x) being the step function. The derivative of θ is given by,
        
            ∂_k_⊥θ_nk=-π sgn(δ_nk)δ(k_⊥-k_F).
        
        Generally, we have ∇θ_nk=-πv_nksgn(δ_nk)δ(ϵ_nk).
         Omitting the interband pairing, the widing number in 1D can be further written as<cit.>,
        
            ν_1   =   i/2π∫_BZdkTr[q^-1∂_kq]=
            i/2π∑_n∫_BZdke^-iθ_nk∂_ke^iθ_nk=-∑_n1/2π∫^2π_0dk ∂_k θ_nk
               =   1/2∑_ϵ_nk=0sgn[∂_kϵ_nk]sgn[δ_nk].
        
        Here ν_1 is a real number. From the definition, only the bands crossing the Fermi level can contribute the winding number. For gapped bulk system with gapless localized edge states, the dominant contribution in the slab model comes from the edge states. In the following, we define the site-resolved winding number to identify the topological phase transition on the edges.
        
        
        
         §.§ spatial dependent winding number
        
        For surfaces or edges in the slab model, we introduce the spatial dependent winding number to identify the phase transition on the surfaces or edges. Generally we introduce Ψ̅^†_k=(ψ^†_k,11,ψ^†_k,12,...,ψ^†_k,1m,...,ψ^†_k,iγ,...,ψ^†_k,Nγ,ψ_-k,1γ̅,...,ψ_-k,iγ̅,...,ψ_-k,Nγ̅), and the BdG Hamiltonian for the slab can be written as,
        
            ℋ^slab_BdG=∑_kΨ̅^†_kh̅_BdG(𝐤)Ψ̅_k, 
                  h̅_BdG(k)=([      h̅_σ(k)         Δ̅_k;       Δ̅^†_k -h̅^*_σ̅(-k);              ]),
        
        where i denotes the lattice site (sublattice index), N is the lattice size, γ denotes the other indices(orbital or spin) and γ̅ represent the index with flipped spin. The chiral symmetry operator is S=U_S=τ_y ⊗1_Nm and it can be expressed as Sh̅_BdG(k)S^-1=-h̅_BdG(k), which results h̅_σ(k)=h̅^*_σ̅(-k) and Δ̅^†_k=Δ̅_k. Under a unitary transformation U, the chiral symmetry matrix become diagonal and the Hamiltonian matrix and Q matrix become off-diagonal,
        
            U^† U_S U   =   (
            [  1  0;  0 -1 ])
            
            U   =   1/√(2)(
            [   1   1;  i1 -i1 ]),
            h̃̅̃_BdG(k)=U^†h̅_BdG(k)U   =   (
            [               0 h̅_σ(k)-iΔ̅^†_k;   h̅_σ(k)+iΔ̅_k               0 ]),
            
            Q^σσ̅_k   =   (
            [         0   q^σσ̅_k; q^σσ̅,†_k         0 ]).
        
        Here q^σσ̅_k is Nm× Nm matrix and N is lattice size and m is the dimension of tight-binding Hamiltonian. The new basis is  Ψ̃̅̃^†_k=Ψ̅^†_kU=(ψ^†_k,11+iψ_-k,11̅,...,ψ^†_k,1m+iψ_-k,1m̅,...,ψ^†_k,iγ+iψ_-k,iγ̅,...,ψ^†_k,Nγ+iψ_-k,Nγ̅,ψ_-k,1γ̅-iψ^†_k,1γ,...,ψ_-k,iγ̅-iψ^†_k,iγ,...,ψ_-k,Nγ̅-iψ^†_k,Nγ) and each has a definite lattice site index. The total winding number can be written as,
        
            ν_1   =   i/2π∫_BZdkTr[q^-1_k∂_kq_k]
               =   i/2π∑_iγ∫_BZdkTr[q^-1_k∂_kq_k]_iγ,iγ
               =   ∑_i=1^Nν^i_1,
        
        with site-resolved winding number ν^i_1=i/2π∑_γ∫_BZdk[q^-1_k∂_kq_k]_iγ,iγ.
        The integral can be numerically evaluated on a discretized Brillouin zone.
        
        For (010) slab of As1 layers, the total winding number for is always zero as shown in Fig.<ref>(a), therefore it cannot be used to identify the phase transition. We calculated the corresponding spatial dependent winding numbers and Fig.<ref>(b) shows the typical profiles in two topological phases. The profiles in two phases are significant different with each other. To further quantitatively characterize the edge phase transition, we separate the top and bottom edges and define the winding number for them ν^T_1 and ν^B_1 by summing the winding number of the upper half and lower half edges, respectively.
        Although ν^T,B_1 are generally not an integer, the its change is an integer during the topological phase transition. For the (010) edge, with increasing Δ_1, the sign change of superconducting pairing on edge states will change winding number of top edge and bottom edge by -1 and 1 in the thermodynamic limit (demonstrated in the main text), respectively. Meantime, the winding number of left and right edge for (100) edge donot change. Therefore, the left/right and top/bottom edges are in topologically distinct phases and when they meet at corners, Majorana modes will appear. The representative local density of states at corners and edges are displayed in Fig.<ref>(b), where Majorana modes will induce a zero-bias peak at corners.
