Non-Hermitian higher-order topological superconductors in two-dimension: statics and dynamics

Being motivated by intriguing phenomena such as the breakdown of conventional bulk boundary correspondence and emergence of skin modes in the context of non-Hermitian (NH) topological insulators, we here propose a NH second-order topological superconductor (SOTSC) model that hosts Majorana zero modes (MZMs). Employing the non-Bloch form of NH Hamiltonian, we topologically characterize the above modes by biorthogonal nested polarization and resolve the apparent breakdown of the bulk boundary correspondence. Unlike the Hermitian SOTSC, we notice that the MZMs inhabit only one corner out of four in the two-dimensional NH SOTSC. We extend the static MZMs into the realm of Floquet drive. We find anomalous $\pi$-mode following low-frequency mass-kick in addition to the regular $0$-mode that is usually engineered in a high-frequency regime. We further characterize the regular $0$-mode with biorthogonal Floquet nested polarization. Our proposal is not limited to the $d$-wave superconductivity only and can be realized in the experiment with strongly correlated optical lattice platforms.

The realm of topological quantum matter is transcended from the Hermitian system to the non-Hermitian (NH) system due to the practical realization of TI phases in meta-materials [67][68][69][70] where energy conservation no longer holds [71,72]. The NH description has a wide range of applications, including systems with source and drain [73,74], in contact with the environment [75][76][77], and involving quasiparticles of finite lifetime [78][79][80]. Apart from the complex eigenenergies and nonorthogonal eigenstates, the NH Hamiltonians uncover a plethora of intriguing phenomena in TI [72,[81][82][83][84] that do not have any Hermitian analog. For instance, NH Hamiltonian becomes non-diagonalizable at exceptional points (EPs) where eigenstates, corresponding to degenerate bands, coalesce [85,86]; line and point are two different types of gaps in these systems that can be adiabatically transformed into a Hermitian and NH systems, * arnob@iopb.res.in † tanay.nag@physics.uu.se respectively [82]; the conventional Bloch wave functions do not precisely indicate the topological phase transitions under the open-boundary conditions (OBCs) leading to the breakdown of the BBC [87][88][89][90][91][92][93]; consequently, the non-Bloch-wave behavior results in the skin effect where the bulk states accumulate at the boundary [87][88][89]94], and the structure of topological invariants become intricate [81,[95][96][97]. The EPs are studied in the context of Floquet NH Weyl semimetals [98,99]. While much has been explored on the HOTI phases in the context of NH systems [100][101][102][103][104][105][106][107][108][109], HOTSC counterpart, along with its dynamic signature, is yet to be examined. Note that NH 1D nanowire with s-wave pairing and p-wave SC chain are studied for the Majorana zero modes (MZMs) [110][111][112][113][114][115][116]. We, therefore, seek the answers to the following questions that have not been addressed so far in the context of proximity induced HOTSC with nonhermiticity -(a) How does the BBC change as compared to the Hermitian case? (b) Can one use the concept of biorthogonal nested-Wilson-loop to characterize the MZMs there similar to that for HOT electronic modes [101]? (c) How can one engineer the anomalous FHOTSC phase for the NH case?
Considering the NH TI in the proximity to a d-wave superconductor, we illustrate the generation of the NH second-order topological superconductor (SOTSC). The breakdown of BBC is resolved with the non-Bloch nature of the NH Hamiltonian, where phase boundaries, obtained under different boundary conditions, become concurrent with each other (see Fig. 1). The SOTSC phase is characterized by the non-Bloch nested polarization. We demonstrate the NH skin effect where MZMs and bulk modes both display substantial corner localization (see Fig. 2). We further engineer the regular 0and anomalous π-mode employing the mass-drive in highand low-frequency regimes, respectively (see Figs. 3 and 4). We characterize the regular dynamic 0-mode by the non-Bloch Floquet nested polarization.
provided |E(k EP )| = 0 which is in a complete contrast to the Hermitian system with E(k) = 0 at gapless point. A close inspection of Eq. (S1) suggests that four-fold degenerate energy bands yield |E(0, 0)| = 0 [|E(π, π)| = 0] for  break-down of conventional BBC due to the non-Bloch nature of the NH Hamiltonian [87,[91][92][93]. This apparent ambiguity in BBC affects the calculation of topological invariants, which we investigate below. Fig. 2 (a) demonstrates the complex-energy profile Im[E] vs Re[E] of Hamiltonian (S1) in real space for m 0 <|m 0 |. We find the line gap for the NH system irrespective of the boundary conditions as the complexenergy bands do not cross a reference line in the complexenergy plane. The origin, marked with red dot in Fig. 2 (a), indicates the MZMs under OBC that are further shown the by the eight mid-gap states in Re[E]-m (state index) plot (see Fig. 2 (b)). Analyzing the local density of states (LDOS) of the above MZMs, we find sharp localization only at one corner out of the four corners [100] (see Fig. 2 (c)). This is a consequence of the mirror-rotation symmetries M xy or M xȳ even though MZMs spatially coincide. There might be additional protection from the bulk modes coming due to the emergent short-range nature of the superconducting gap [118]. The MZMs are localized over more than a single corner when M xy or M xȳ is broken [117]. The MZMs are also found to be robust against onsite disorder that respects mirror-rotation and chiral symmetries (see supplemental material [117]). In addition, we remarkably find that the LDOS of the bulk modes also exhibits substantial corner localization as depicted in the insets of Fig. 2(c) [87][88][89]. The above features, reflecting the non-Block nature of the system, are referred to as the NH skin effect [87,88]. This is in contrast to the Hermitian case where only the zero-energy modes can populate four corners of the 2D square lattice [13,43,47].
Topological characterization.-To this end, in order to compute the topological invariant from H(k) characterizing the SOTSC phase under OBC, we exploit the non-Bloch nature. We need to use the complex wavevectors to describe open-boundary eigenstates such that [81]. Upon replacing k x,y → k x,y −iγ x,y /λ x,y , the renormalized topological mass m 0 acquires the following form in the limit k x,y → 0 and γ x,y → 0 (2) Note that |m 0 | = |m 0 − m 0 | denotes phase boundary of the SOTSC phase as obtained from Fig. 1 from the non-Bloch NH Hamiltonian H (k ) [119,120].
being the the number of discrete points considered along i-th direction and e i being the unit vector along the said direction. Notice that, the bi-orthogonalization guarantees the following n |Ψ R n (k ) Ψ L n (k )| = I and Ψ L n (k )|Ψ R n (k ) = δ mn ; where, n runs over all the energy levels irrespective of their occupations. The first-order polarization ν x,µ (k y ) is obtained from the eigenvalue equation for W x,k as follows Note that unlike the Hermitian case, here W x,k is no longer unitary resulting in ν x,µ (k y ) to be a complex number [101]. Importantly, |ν R x,µ (k ) ( ν L x,µ (k )|) designates bi-orthogonalized right (left) eigenvector of W x,k associated with µ = 1, · · · , 4-th eigenvalue. For a (second-order topological) SOT system, the real part of first-order polarization exhibits a finite gap in spectra such that it can be divided into two sectors as ±ν x (k y ) where each sector is two-fold degenerate. Such a structure of Wannier centres in the non-Bloch case might be relied on the mirror symmetry of the underlying Hermitian Hamiltonian [10,66]. In order to characterize the SOT phase, we calculate the polarization along the perpendicular y-direction by projecting onto each ±ν x branch. This allows us to employ the nested Wilson loop as follows [10,66] x,µ2 (k ) n with F y,k mn = Ψ L m (k + ∆ y e y )|Ψ R n (k ) . The indices µ 1,2 ∈ ±ν x run over the projected eigenvectors of W x,k only. We evaluate W ±νx y,k for a given value of k x that is the base point while calculating W x,k (3).
The nested polarization ν ±νx y,µ (k x ) can be extracted by solving the eigenvalue equation for W ±νx The average nested Wannier sector polarization ν ±νx y,µ for the µ -th branch, characterizing the 2D SOTSC, is given by We explore the SOT phase diagram by investigating Fig. 1(d)). The blue (brown) region indicates the SOTSC and trivial phase. The green line in Fig. 1 (d), separating the above two phases, represents the phase boundarym 0 = 2 + γ 2 as demonstrated in Eq. (2). On the other hand, the phase boundaries, obtained from bulk Hamiltonian (S1), are found to be m +,± 0 = 2 ± √ 2γ that are depicted by yellow lines in Fig. 1 (d). Therefore, the topological invariant, computed using the non-Bloch Hamiltonian H (k ), can accurately predict the MZMs as obtained from the real space Hamiltonian under OBC (see Fig. 1 (c)). This correspondence for very higher values of γ no longer remains appropriate due to the possible break down of Eq. (2). Even though M x,y are broken, ν ±νx y,µ ( ν ±νy x,µ ) yields half-integer quantization provided M xy (M xȳ ) is preserved. Note that based on mirror rotation and sublattice symmetries, the NH SOTSC can be shown to exhibit integer quantization in winding number similar to NH SOTI [100] (see supplemental material [117]).
Floquet generation of NH SOTSC.-Having studied the static NH SOTSC, we seek the answer to engineer dynamic NH SOTSC out of trivial phase by periodically kicking the on-site mass term of the Hamiltonian H(k) (Eq. (S1)) as [57,64] Here, m 1 and T represent the strength of the drive and the time-period, respectively. The Floquet operator is formulated to be where TO denotes the time ordering. Notice that m 0 > |t x + t y + γ 2 x + γ 2 y | such that the underlying static NH Hamiltonian H(k) remains in the trivially gapped phase. Having constructed the Floquet operator U (k, T ), we resort to OBCs and diagonalize the Floquet operator to obtain the quasi-energy spectrum for the system. We depict the real part of the quasi-energy µ m as function of the state index m in Fig. 3 (a) where frequency of the drive is higher than the band-width of the system. The existence of eight MZMs is a signature of the NH Floquet SOTSC phase. The LDOS for the MZMs displays substantial localization only at one corner in Fig. 3(b).  Fig. 1 (c).
Insets show the NH skin effect where the bulk modes at finite energy also have a fair amount of corner localization.
In order to topologically characterize the above MZMs, we again make use of the non-Bloch form. Instead of the static Hamiltonian, we derive the high-frequency effective Floquet Hamiltonian, in the limit T → 0 and m 1 → 0 to analyze the situation with Γ 21 = −σ y s z , Γ 31 = σ x , Γ 41 = −τ y σ z . Upon substitution of k → k +iβ, the modified mass term in H Flq (k ) reads as Evaluating the effective Floquet nested Wannier sector polarization ν F,±νx y,µ numerically from non-Bloch Floquet operator U (k , T ) [48,66], we obtain the Floquet phase diagram in m 1 -γ plane as shown in Fig. 3 (c). The non-Bloch Floquet operator can be considered as the dynamic analog of the non-Bloch NH Hamiltonian H (k ). In particular, we use bi-orthogonalized |Ψ R F (k ) ( Ψ L F (k )|), representing the occupied right (left) quasistates of U (k , T ) with quasi-energy −π/T < Re[µ] < 0, to construct the Wilson loops W F x,k for the driven case. eigenvector of W F x,k . Interestingly, this is similar to the static phase diagram where the phase boundary is accurately explained by Eq. (11).
We further analyze the problem for lower frequency regime to look for anomalous Floquet modes at quasienergy Re[µ] = ±π [57,66]. We depict one such scenario for Ω = 2π/T = 5.0 in Fig. 4 (a) where eight anomalous π-modes appear simultaneously with regular eight 0-modes. The corresponding LDOS for 0-mode and π-mode are shown in Fig. 4(b) and (c), respectively. Interestingly, the 0-mode and the π-mode populate different corners of the system. As a signature of NH skin effect, we show the LDOS for two bulk states in the inset I1 and I2 of Fig. 4(b). The localization profile of the zero-energy states and bulk states are unique to the NH system that can not be explored in its Hermitian counterpart.
Discussions.-The number of MZMs can be tuned in our case by the application of magnetic field similar to the Hermitian SOTSC phase [64]. The long-range hopping provides another route to enhance the number of MZMs that can in principle be applicable for the non-Hermitian case as well [121,122]. Interestingly, Floquet driving delivers an alternative handle to generate long-range hopping effectively out of the short-range NH model such that the number of MZMs are varied (see supplemental material [117]). Interestingly, Hermitian and non-Hermitian phases belong to the Dirac and non-Hermitian Dirac universality classes [123,124]. In the case of HOT phases, one expects different critical exponents with re-spect to the usual Dirac model. The breakdown of BBC and skin effect are intimately related to such non-Hermitian Dirac universality class. The edge theory, computed from Hermitian HOT model, is modified due to the non-Hermiticity with the possible non-Bloch form. Given the experimental realization of spin-orbit coupling [125,126], non-Hermiticity [127,128] and theoretical proposals on topological superfluidity [129,130] in optical lattice, we believe that the cold atom systems might be a suitable platform for the potential experimental realization of our findings [74,131,132]. However, we note that the superconductivity might be hard to achieve in the NH setting.
Summary and conclusions.-In this article, we consider 2D NH TI, proximized with d-wave superconductivity, to investigate the emergence of NH SOTSC phase. From the analysis of EPs on the bulk NH Hamiltonian under PBC, one can estimate the gapped and gapless phase in terms of the real energies (see Fig. 1). By contrast, the MZMs, obtained from the real space NH Hamiltonian under OBC, do not immediately vanish inside the bulk gapless region (see Fig. 2). This apparent breakdown of the BBC can be explained by the non-Bloch nature of the NH Hamiltonian that further results in the MZMs residing at only one corner while the bulk modes populate the boundaries. While the later is dubbed as NH skin effect. We propose the nested polarization for topologically characterize the MZMs upon exploiting the non-Bloch form of the complex wave vectors. This resolves the anomaly between the phase boundaries, obtained from OBC and PBC, in the topological phase diagram. Finally, we adopt a mass-kick drive to illustrate the Floquet generation of NH SOTSC out of the trivial phase and characterize it using non-Bloch Floquet nested Wannier sector polarization (see Fig. 3). In addition, we demonstrate the emergence of anomalous π-mode following such drive when the frequency is lowered (see Fig. 4). The mirror symmetries M x,y play crucial role in characterizing the anomalous π-modes [48,66]. Therefore, such characterization in the absence of mirror symmetries is a future problem.
Acknowledgments.-A.K.G. acknowledges SAMKHYA: High-Performance Computing Facility provided by Institute of Physics, Bhubaneswar, for numerical computations. We thank Arijit Saha for useful discussions. In this section, we show that how the first-and second-order Wannier centres behave under the above spatial symmetries. The first-order Wannier centres ν x (k y ) and second-order Wannier centres ν νx y behave in the following way [10]-mirror-rotation along x-direction M x : ν x (k y ) → −ν x (k y ) [see Fig. S2], and ν νx y (k x ) → ν −νx y (−k x ); M x causes the first-order branches to appear in pairs, mirror-rotation along y-direction. M y : ν x (k y ) → ν x (−k y ), and ν νx y (k x ) → −ν νx y (k x ); M y defines the shape of the first-order branches. The four-fold rotation C 4 and mirror rotations M xy , M xȳ interchange the branches, C 4 : ν x (k y ) → −ν y (k x ), and ν νx y (k x ) → ν −νy x (−k y ), M xy : ν x (k y ) → ν y (k x ), and ν νx y (k x ) → ν νy x (k y ) M xȳ : ν x (k y ) → −ν y (−k x ), and ν νx y (k x ) → −ν −νy x (−k y ). In our case, the first-order branches do not follow the above relations as M x,y symmetries are broken. However, the second-order polarization i.e., nested polarization mod( ν ±νx y,µ , 1.0) and mod( ν ±νy x,µ , 1.0) both are found to be 0.5 following the mirror rotation symmetries generated by M xy and M xȳ . Moreover, under these mirror rotation symmetries, one can show that ν νy . We further notice that when M xy is broken, ν ±νx y,µ and ν ±νy x,µ do not produce a quantized value of 0/0.5. Therefore, these spatial symmetries can predict the quantization of the second-order polarization that determines the second-order topology. In this regard, we would like to comment that the winding number, based on the mirror rotation symmetry M xy and the sublattice symmetry S, is shown to topologically characterize the NH higher-order topological insulator (HOTI) [100]. We believe that it is possible to define such a winding number yielding quantized values in the SOTSC phase as our NH model preserves mirror rotation and sublattice symmetries.
Note that winding number, usually characterizing a first-order topological phase provided Hamiltonian preserves sublattice symmetry [133], can also identify a second-order topological phase while appropriately defined with other symmetry constraints. The winding number for a Hamiltonian physically means how many times the Hamiltonian encircles a given point when a parameter is cyclically varied. On the other hand, half integer eigenvalues of nested Wilson loop refers to the fact that Wannier centers, located at a high-symmetry point with respect to the mirror symmetries, lie half-way between the lattice sites for higher-order topological phase [10]. Therefore, both the above invariant captures the non-trivial nature of the underlying wave-functions in the Brillouin zone through either Wannier might be intimately connected to winding number as stated above. However, the detailed mathematical proof is an open question yet, to the best of our knowledge, that we leave for future studies.

S4. s-WAVE NH SOTSC
In the main text, we discuss d-wave proximized to engineer NH SOTSC phase. We here illustrate that NH higherorder topological superconductor (HOTSC) can be contemplated for s-wave superconductivity with pairing amplitude ∆ s as follows [28,63] H(k) = N·Γ; where, N = {λ x sin k x +iγ x , λ y sin k y +iγ y , m 0 −t x cos k x −t y cos k y , ∆ s , Λ(cos k x − cos k y )}, Γ = {τ z σ x s z , τ z σ y s 0 , τ z σ z s 0 , τ x σ 0 s 0 , τ 0 σ x s x }. The last term proportional to Λ represents C 4 symmetry breaking Wilson-Dirac mass term. We depict the corresponding complex-energy bands in complex energy plane, eigenvalue spectra highlighting the mid-gap states, and the LDOS associated to the MZMs and bulk modes in Fig. S3 (a), (b), and (c), respectively. The NH SOTSC phase with ∆ s can also be obtained using in-plane magnetic field instead of C 4 symmetry breaking [30,63]. The topological characterization of these phases we leave for future studies.  In this section, we investigate the effect of on-site disorder to investigate the stability of the MZMs. We consider an onsite disorder potential of the form V (i, j) = i,j V ij Γ 3 . Here, V ij is randomly distributed in the range V ij ∈ − w 2 , w 2 ; while w accounts for the strength of the disorder potential. The disordered NH Hamiltonian H (i, j) = H(i, j) + V (i, j)Γ 3 respects the chiral symmetry while investigating the extended Hermitian Hamilto-nianH : SH S −1 = −H with S = µ 0 τ y σ 0 s 0 . We also note that on-site disorder preserves M xy . In the presence of the onsite disorder, the real-space Hamiltonian H is given as Here, Ψ i,j is a 8 × 1 matrix consisting of the annihilation operator in the particle-hole, orbital, and spin subspace at a lattice site (i, j). And i, j x ( i, j y ) represents the nearest neighbor hopping along x (y) direction. We consider 500 disorder configurations and depict the eigenstates E m as a function of the state index m and the LDOS of the MZMs in Fig. S4 for different disorder strengths. The MZMs remain localized at a given corner due to the mirror symmetry preserving nature of the on-site disorder. We can conclude from the eigenvalue spectra and the LDOS plot that the MZMs are robust and against the onsite disorder.

S6. TUNING THE NUMBER OF MCMS DYNAMICALLY
In this section, we discuss how Floquet engineering permits us to generate more in-gap states (both 0-and π-gap) by tuning the different driving parameters suitably. The generation of multiple in-gap states is attributed to higher-order hoppings generated in the dynamical system. The mass-kick protocol, introduced in Eq. (9) of the main text, can be used to create multiple Majorana corner modes (MCMs). We depict the real part of the quasienergy spectra as a function of the driving frequency Ω in Fig. S5 (a). In order to understand the generation of the MCMs more clearly, we exemplify two cases in Fig. S5 (b) and (c) for Ω = 4.6 and 1.9, respectively. In Fig. S5 (b), the number of 0-MCMs is twice that of the static case, whereas the number of π-MCMs remains the same. While in Fig. S5 (c), the number of π-MCMs is thrice compared to the case in the main text (see Fig. 4 (a) of the main text).