# Mixed higher-order topology and nodal and nodeless flat band topological
phases in a superconducting multiorbital model

Rodrigo Arouca1 Tanay Nag2,1<EMAIL_ADDRESS>Annica M.
Black-Schaffer1 1Department of Physics and Astronomy, Uppsala University,
75120 Uppsala, Sweden  2Department of Physics, BITS Pilani-Hyderabad Campus,
Telangana 500078, India

###### Abstract

We investigate the topological phases that appear in an orbital version of the
Benalcazar-Bernevig-Hughes (BBH) model in the presence of conventional spin-
singlet $s$-wave superconductivity and with the possibility of tuning an in-
plane magnetic field. We chart out the phase diagram by considering different
boundary conditions, with the topology of the individual phases further
examined by considering both the Wannier and entanglement spectra, as well as
the Majorana polarization. For weak to moderate values of magnetic field and
superconducting pairing amplitude, we find a second-order topological
superconductor phase with eight zero-energy corner modes. Further increasing
field or pairing, half of the corner states can be turned into zero-energy
edge-localized modes, thus forming a type of hybrid-order phase. Then, we find
two different putative first-order topological phases, a nodal and a nodeless
phase, both with zero-energy flat bands localized along mirror-symmetric open
edges. For the nodal phase, the flat bands are localized between the nodes in
reciprocal space, while in the nodeless phase, with its a full bulk gap, the
zero-energy boundary flat band spans the whole Brillouin zone.

##  I Introduction

The study of topological materials is an extremely active area of research in
condensed matter physics. They present phases of matter that are not
characterized by spontaneous symmetry breaking but rather by topological
invariants. In the Altland-Zirnbauer classification [1, 2], time-reversal,
particle-hole, and chiral symmetries classify ten possible topological
classes, indicating the kind of invariant and the branches of the symmetry-
protected boundary states. The number of possible symmetry-protected
topological classes of free fermions has further been increased by including
crystalline symmetries to this original classification [3, 4, 5, 6, 7].

In addition to new topological classes, crystalline symmetries also allow for
the presence of higher-order topological phases [8, 9, 10, 11], where the
topological invariant computed in the bulk is not related to modes appearing
on the whole surface of the material but rather on a smaller set. As an
example, the original model of a higher-order topological insulator, the two-
dimensional (2D) Benalcazar-Benervig-Hughes (BBH) model [8, 9], hosts a
second-order topological phase with protected zero-energy modes appearing at
the corners of the system. The last years have seen a profusion of work on
higher-order topological insulators [12, 13, 14, 15, 16, 17, 18, 19, 20, 21,
22, 23], with experimental realization in materials [11, 24] and a variety of
metamaterials [25, 26, 27, 28, 29, 30, 31, 32, 33].

A particularly interesting class of higher-order topological systems is
composed of higher-order topological superconductors (HOTSC), both for static
and driven Hamiltonians [34, 35, 36, 37, 38, 39, 40, 41, 42]. The corner or
hinge states in HOTSC appear at zero energy and can, as such, host Majorana
zero modes (MZMs) [43, 44, 45, 46, 47, 48, 47, 49, 50], where their spatial
separation can be controlled and degree of degeneracy can also be tuned under
certain circumstances. As such, they must present an equal amount of particle
and hole components and are thus their own antiparticles, which are promising
for applications in quantum computing due to their non-Abelian properties
[51]. Topologically protected boundary states can also appear in nodal
superconductors [52, 53, 54, 55, 56, 57, 58, 59]. Topological nodal
superconductor is a non-trivial phase that is not contained in the Altland-
Zirnbauer classification since the bulk of the system is gapless at the nodal
points [54, 55, 57]. Nevertheless, the presence of the nodal points is
protected by symmetry, and there exists a bulk-boundary correspondence between
the topology of the bulk nodal points and boundary-localized flat bands
located between the nodes. It would be interesting to uncover systems where
multiple different superconducting topological phases are readily realized,
including higher-order topology or various nodal states. This would not only
give realizations of the individual topological phases, and possible
intriguing combinations thereof, but also offer tunability between the
different phases and their characteristic properties as well as the
possibility to uncover previously unknown phases.

Figure 1: (a) Schematic realization of the superconducting BBH model: BBH
system with conventional spin-singlet $s$-wave superconducting pairing
$\Delta$ induced by proximity effect from a substrate and an in-plane magnetic
field $B_{x}$. Different edges are indicated by Roman numerals. Dashed area
represents the region where the entanglement spectrum is computed. (b) Summary
of phase diagram: Solid lines represent analytical expressions for the phase
boundaries, while different colors represent different topological phases. See
main text for further definitions.

In this work, we set out to study a simple and versatile superconducting
model, and we find a plethora of different topological phases. Aiming at least
for higher-order topology, we choose to investigate the orbital version of the
BBH model in the deep topological limit with conventional spin-singlet
$s$-wave superconductivity induced by proximity effect from a substrate and
using an applied in-plane magnetic field as an additional easily accessible
tunable parameter, all illustrated in Fig. 1(a). The in-plane magnetic field
$B_{x}$ breaks the $C_{4}$ symmetry responsible for protecting the corner
states in the BBH model, while the proximity-induced superconducting order
parameter $\Delta$ transforms these states in electron-hole excitations.
Moreover, we find that the superconducting term in this orbital model is
represented by an unusual matrix structure, giving rise to a multitude of
different topological phases. Using the Wannier spectrum [7, 9, 11], the
entanglement spectrum [60, 61, 62, 63, 64] for a quarter of the lattice
indicated as in Fig. 1(a), and the Majorana polarization [65, 66], we
completely characterize the topological phases and obtain the rich phase
diagram in Fig. 1(b).

To briefly summarize the phase diagram, for $B_{x}=0$ and $\Delta\approx 0$
the result is the superconducting version of the second-order topological
phase of the BBH model characterized by eight corner states and shown as the
HOTSC phase (yellow) in Fig. 1(b). The pairing makes the corner states of the
BBH model transform into Andreev bound states (ABS), built up from two
degenerate MZMs, located at each corner. The presence of both finite pairing
and an in-plane magnetic field makes even more interesting topological phases
appear. For larger values of both $\Delta$ and $B_{x}$ we first find another
higher-order topological phase, a type of hybrid ordered phase (green). In
this phase, some (four) of the original zero-energy corner states stay,
becoming isolated MZMs, while the rest turn into MZMs localized along in the
$y$-direction, edges II and IV, in Fig. 1(a). This presents an intriguing mix,
or hybrid, of a second-order and a dipolar topological phase, but where the
number edge localized states do not grow with system size as usually expected.
Further increasing $B_{x}$ and $\Delta$, we find two other interesting
topological phases. These are both first-order topological phases with flat
bands at zero energy. One of the phases is more traditional, a nodal flat band
phase (blue). Here the bulk of the system is gapless with symmetry-protected
nodes and the zero-energy flat bands localized on the $y$-edge, occur in the
region of momentum space between the nodes, as expected for a topological
nodal superconductor. However, in the other topological flat band phase the
bulk stays completely gapped, but we nevertheless find flat bands at zero
energy, now spanning across the whole Brillouin zone and localized at the
edges along the $x$-direction, which we name a nodeless flat band phase (red).
The association of these flat bands to a bulk invariant remains unclear, but
we find that they clearly appear in a quantized Wannier spectrum. Both flat
band phases present four zero-energy states per momenta, two at each edge,
with a large Majorana polarization, such that we consider them to be MZMs.
These results show how a seemingly simple model can generate a plethora of
topological phases, accessible by tuning physical parameters. Importantly, the
results further establish the existence of both a hybrid-order phase and a
fully gapped (i.e. nodeless) topological phase with zero-energy flat bands
boundary states.

The rest of this work is structured as follows. In Sec. II we briefly
summarize the main properties of the BBH model and discuss the form of the
Hamiltonian, especially the superconducting pairing and the magnetic field
term, which we add to the BBH model. In Sec. III we introduce the three
topological invariants we use in this work, explaining what aspects of
topological phases can be characterized by each invariant. In Sec. IV we
report our main results. First, we consider gap closings as a function of
$\Delta$ and $B_{x}$, which indicate the phase boundaries. For each phase, we
analyze the behavior of the energy spectra for different boundary conditions
as well as the different topological invariants, in order to characterize all
phases. In Sec. V we discuss the dual [67] and surface [35, 59] Hamiltonians
for this model, which help us understand the presence of nodes starting with
$s$-wave pairing and the higher-order phases of the model. Finally, in Sec. VI
we conclude and summarize our main results.

##  II Model

The paradigmatic model of a higher-order topological insulator, the BBH model
[9, 8], is described by the Bloch Hamiltonian

$h_{\rm BBH}(\mathbf{k})=\mathbf{d}(\mathbf{k})\cdot\bm{\gamma},$ (1)

where $d_{1}=t\sin(k_{y}a)$, $d_{2}=\lambda+t\cos(k_{y}a)$,
$d_{3}=t\sin(k_{x}a)$, $d_{4}=\lambda+t\cos(k_{x}a)$,
$\gamma_{i}=\sigma_{2}\otimes s_{i}$ for $i=1,2,3$, and
$\gamma_{4}=\sigma_{1}\otimes s_{0}$. The model is defined on a square lattice
with a lattice parameter $a$. The original BBH model can be considered as the
collection of spinless fermions on a square lattice with four sublattices,
being an extension of the two-dimensional (2D) Su-Schrieffer-Heeger model [68]
with a $\pi$-flux. Here, we instead interpret $\sigma_{i}$ and $s_{i}$ as
Pauli matrices ($\sigma_{0}=s_{0}=\mathbbm{1}_{2}$) associated with orbital
and spin degree of freedom. This Hamiltonian then has spin-orbit coupling
between orbital and spin on the same site (terms proportional to $\lambda$)
and between different sites (terms proportional to $t$). This choice is
motivated for considering the proximity effect from a conventional
superconductor, where the pairing is $s$-wave (onsite) intra-orbital
(proportional to $\sigma_{0}$) spin-singlet (proportional to $s_{y}$), which
is not possible for spinless fermions. We here further consider the deep
topological limit of Eq. (1), by setting $\lambda=0$. This both generates
nontrivial topology in the normal state and removes the need for on-site spin-
orbit coupling.

Since it is known that a conventional superconductor with spin-orbit coupling
can host topological superconductivity in an in-plane magnetic field [67, 44,
45], we also add an in-plane external magnetic field $B_{x}$ along $x$,
accounted for by a Zeeman term proportional to $s_{x}$. An interesting
question then is to investigate the interplay between superconductivity and
magnetic field with the intrinsic higher-order topology of the BBH model. The
Bogoliubov-de Gennes (BdG) Hamiltonian for this system in the particle-hole
basis becomes

$h_{\rm
BdG}(\mathbf{k})=\left(\begin{matrix}\mathbf{h}(\mathbf{k})&-i\mathbf{\Delta}\\\
i\mathbf{\Delta}&-\mathbf{h}^{*}(-\mathbf{k})\end{matrix}\right)=\,\mathbf{D}(\mathbf{k})\cdot\mathbf{\Gamma},$
(2)

where $h(\mathbf{k})=h_{\rm BBH}(\mathbf{k})+B_{x}\,\sigma_{0}\otimes s_{x}$
is the total normal-state Hamiltonian. For the second equality we use the
vectors $\mathbf{D}=(d_{1},d_{2},d_{3},d_{4},B_{x},\Delta)$ and
$\mathbf{\Gamma}$ to write the Hamiltonian in terms of matrices in the
particle-hole, orbital, and spin degrees of freedom. Using $\tau_{i}$ as the
Pauli matrices in the particle-hole degrees of freedom, we construct the
matrices $\Gamma_{i}=\tau_{3}\otimes\gamma_{i}$,
$\Gamma_{4}=\tau_{3}\otimes\gamma_{4}$,
$\Gamma_{5}=\tau_{2}\otimes\sigma_{0}\otimes s_{2}$, and
$\Gamma_{6}=\tau_{3}\otimes\sigma_{0}\otimes s_{1}$. Before investigating the
different phases of this system in detail, we next discuss the topological
invariants that we use to characterize them.

##  III Topological characterization

As hinted by the phase diagram in Fig. 1(b), the model Eq. (1) presents a
variety of topological phases with a mix of first- and second-order as well as
nodal topological phases. In addition, since we are dealing with a spinful
topological superconductor, we can have different symmetry-protected boundary
states. Therefore, we do a thorough analysis of the topology using three
different invariants, which can identify different aspects of the topological
phases.

We use the Wannier spectrum [8, 7] to investigate the presence of a non-
trivial polarization in the lattice. Since the presence of higher-order
topological phases is not completely characterized by the Wannier spectra, we
also use a real-space invariant, the entanglement spectrum [60, 69], to verify
whether the boundary modes are of higher-order topological origin. Finally, to
understand whether these boundary states may be MZMs or ABS, we use the
Majorana polarization [65, 66], which indicates how much of a combination of
electron and hole a state in a superconductor is. In addition to these
invariants, we study the energy spectra using different boundary conditions,
the local density of states (LDOS) at zero frequency, and wavefunctions of the
states in the middle of the spectrum to further characterize and verify the
phases. Below, we briefly review these topological invariants and explain the
signatures of each topological phase.

###  III.1 Wannier spectrum

In the modern theory of polarization, the charge polarization is directly
related to the Berry phase [70, 71, 72, 7] being, therefore, a topological
property of a material. However, the direct numerical calculation of these
quantities is not so practical due to an overall ill-definition of the phase
of wavefunctions [7]. A very convenient alternative is based on the use of
Wilson loops [73, 8, 7]. For a 2D system with periodic boundary conditions
along $x$ and open boundary conditions along $y$, the Wilson loop components
are defined by [7] 111The Wilson loop, in general, needs to be defined in
terms of a reference point $\mathbf{k}_{i}$ [8] from which the loop is made.
However, for the topological invariants considered here, the Wannier spectrum
is the same for all reference points, and we thus choose $k_{i}=-\pi$ for both
$k_{x}$ and $k_{y}$.

$\left(W\right)_{mn}^{x}=\braket{u_{m}(\pi)}{\prod\limits_{k_{x}}^{\pi\leftarrow-\pi}\mathcal{P}(k_{x})}{u_{n}(-\pi)},$
(3)

where $\mathcal{P}$ is a projector

$\mathcal{P}(k_{x})=\sum\limits_{m}\ket{u_{m}(k_{x})}\bra{u_{m}(k_{x})}$ (4)

over the occupied eigenstates $\ket{u_{m}(k_{x})}$ of the (semi-periodic)
Hamiltonian with momentum $k_{x}$.

$W^{x}$ is a unitary matrix which can be associated with the so-called Wannier
Hamiltonian $H_{W}^{x}$

$W^{x}=e^{i2\pi H_{W}^{x}}.$ (5)

The eigenvalues $\nu^{x}$ of $H_{W}^{x}$ are the Wannier spectrum of the
system with open boundaries along $y$. A gauge transformation, associated with
the change of phases of the wavefunctions, can change $\nu^{x}$ by integer
values, making these quantities, in general, defined mod 1. We choose a gauge
where $\nu^{x}$ take values between $0$ and $1$ [7, 8]. The presence of
crystalline symmetries imposes some constraints in the Wannier spectrum [7,
8]. For instance, symmetry upon reflection along the $x$-axis makes $\nu^{x}$
come in $(\nu,1-\nu)$ pairs [8]. In this way, $\nu^{x}=0.5$ are reflection
invariant and indicate the presence of boundary modes protected by this
symmetry, created by a non-trivial polarization along $x$[7, 8]. The same
holds, mutatis mutandis, to periodic boundary conditions along $y$ with
boundaries open along $x$, obtaining a corresponding $\nu^{y}$. Thus, when
either $\nu^{x}$ or $\nu^{y}$ presents modes at $0.5$, or half-quantized
modes, we obtain a dipolar phase. For a higher-order topological phase, there
are instead half-quantized midgap states in both $\nu^{x}$ and $\nu^{y}$ [8].

An illustrative example is the BBH model in Eq. (1). This model presents
reflection symmetry along both $x$ and $y$, which restricts $\nu^{x/y}$ to
appear in pairs. For $\left|\lambda/t\right|<1$, the system is in a second-
order, or quadrupolar, topological phase characterized by four corner modes at
zero energy, corresponding four $0.5$ eigenvalues in both $\nu^{x}$ and
$\nu^{y}$ [8]. For reference, for $\left|\lambda/t\right|>1$, the model is in
a trivial phase, with no midgap states in both the energy and Wannier spectra
[8]. If we further allow $t$ or $\lambda$ to be different along $x$ and $y$,
we can obtain a phase with polarization just along one of the directions
showing first-order topology.

Before moving on, we discuss why we are not using the nested Wilson loop
spectrum [8] to characterize the higher-order topological phases. The nested
Wilson loop is computed by using the eigenvectors of the Wannier Hamiltonian
$H_{W}$, Eq. (5), but for fully periodic boundary conditions in the general
expression of the Wilson loop in Eq. (3). It is often taken as a clear
indicator of higher-order topology, presenting midgap states in its spectrum
when there are symmetry-protected higher-order modes [8]. For example, in the
BBH model, Eq. (1), the nested Wilson loop is a phase, which is equal to zero
in the trivial phase and $\pi$ in the quadrupolar phase [8]. However, for the
nested Wilson loop spectrum to present quantized values, one needs inversion
symmetry. But our full system in Eq. (2) breaks mirror symmetry along $x$:
$M_{x}h_{\rm BdG}(k_{x},k_{y})(M_{x})^{-1}\neq h_{\rm BdG}(-k_{x},k_{y})$ with
$M_{x}=\tau_{3}\otimes\sigma_{1}\otimes s_{3}$, while it preserves mirror
symmetry along $y$: $M_{y}h_{\rm BdG}(k_{x},k_{y})(M_{y})^{-1}=h_{\rm
BdG}(k_{x},-k_{y})$ with $M_{y}=\tau_{3}\otimes\sigma_{1}\otimes s_{1}$.
Consequently, inversion symmetry, generated by $\mathcal{I}=M_{x}M_{y}$, is
broken as $\mathcal{I}h_{\rm BdG}(k_{x},k_{y})\mathcal{I}^{-1}\neq h_{\rm
BdG}(-k_{x},-k_{y})$. Therefore, we cannot use the nested Wilson loop as a
bulk invariant to diagnose our topological phases. We instead revert to the
Wannier spectra along $x$ and $y$.

###  III.2 Entanglement spectrum

Since we cannot use the nested Wilson loop spectrum, but only the Wannier
spectrum along $x$ or $y$, an alternative invariant that can characterize
higher-order topological phases is useful. One such invariant has recently
turned out to be the entanglement spectrum [60, 61, 62, 74, 63, 64, 75]. In
the same way that a non-trivial polarization in the lattice can be determined
using the Wannier spectrum, it can also be diagnosed by the entanglement
spectrum [60, 61, 62, 74, 63, 64, 75]. To obtain the entanglement spectrum, we
compute the correlation matrix in the occupied state $\ket{\Omega}$

$C_{\mathbf{r},\tau,\sigma,s;\mathbf{r}^{\prime},\tau^{\prime},\sigma^{\prime},s^{\prime}}=\braket{\Omega}{c^{\dagger}_{\mathbf{r},\tau,\sigma,s}c_{\mathbf{r}^{\prime},\tau^{\prime},\sigma^{\prime},s^{\prime}}}{\Omega},$
(6)

where $\ket{\Omega}$ represents the (many-body) fermionic ground state and
$c^{\dagger}_{\mathbf{r},\tau,\sigma,s}$ creates a particle ($\tau=1$) or hole
($\tau=-1$) in orbital $\sigma$ with spin $s$ at position $\mathbf{r}=(x,y)$.
The entanglement spectrum $\xi$ consists of the eigenstates of the correlation
function constrained to some finite region in real space. One can intuitively
understand the cut(s) needed to create such as finite region as creating
artificial boundaries in the system, such that the presence of boundary states
may appear in properties of the entanglement spectrum. For instance, for a
system with inversion symmetry, the entanglement spectrum of a lattice cut in
half displays modes at $0.5$ in the topological phase [76], analogously to
what happens to the Wannier spectrum. A cut that preserves the symmetries that
protect the corner modes can also be used to diagnose the presence of a higher
order topological phase [77, 69]. Since these modes are protected by $C_{4}$
symmetry, we cut the system in half along both $x$ and $y$, obtaining a
quarter of the original lattice, as indicated by the dashed area in Fig. 1(a).
For this cut, the number of $\xi$ modes between $1/4$ and $3/4$ is equal to
the number of $C_{4}$-protected midgap states [77, 69].

###  III.3 Majorana polarization

The Majorana polarization $P$ is defined as [65, 66, 78]

$P_{m}(x,y)=\sum\limits_{\sigma,s}2\psi_{x,y,\tau=1,\sigma,s;m}\psi_{x,y,\tau=-1,\sigma,s;m},$
(7)

where $\psi_{\mathbf{r},\tau=1,\sigma,s;m}$
($\psi_{\mathbf{r},\tau=-1,\sigma,s;m}$) is the particle (hole) component of
the wavefunction of the $m$-th eigenstate at position $\mathbf{r}=(x,y)$ with
orbital $\sigma$ and spin $s$. It is a tool to characterize how much particle-
hole symmetric a state is. In particular, the comparison between the absolute
value of $P$ with the probability density

$\Psi_{m}(x,y)=\sum\limits_{\tau,\sigma,s}|\psi_{x,y,\tau,\sigma,s;m}|^{2},$
(8)

which measures the spatial properties of the wavefunction, quantifies how much
of the wavefunction is particle-hole symmetric. In addition, the phase of $P$
is proportional to the relative phase between the electron and hole components
of the wavefunction. As an example, for the MZMs of the paradigmatic 1D Kitaev
chain, the Majorana polarization is $P=e^{\pm\pi i}$ [66], with the different
signs corresponding to different edges of the system. In contrast, an ABS can
present any value of $P$.

For systems that present only one isolated state per boundary, the Majorana
polarization becomes an unambiguous indicator of a MZM. However, if there are
many MZMs per boundary, since they are degenerate in energy, one may obtain
different values of $P$ for different linear combinations of the states at
zero energy. Thus, even if we numerically find a high value of the Majorana
polarization compared to the probability density, it can be unclear whether
two such putative MZMs can actually recombine into a complex fermion. For such
recombination to not occur, different spin and orbital degrees of freedom
generally have to be in play. Such ambiguousness is the case for some
topological phases in our system and we can thus not use $P$ as a unique
indicator in these cases. Nevertheless, we still investigate the midgap states
in terms of the Majorana polarization to provide an additional tool whenever
it is distinctive. Since we always obtain a real $P$, we use only its real
value to check the sign across the lattice.

###  III.4 Spectral characterization

In addition to the topological invariants discussed above, it is useful to
consider how the energy spectrum and the wavefunctions of the system behave in
every phase. Symmetry-protected topological states appear at zero energy at
the boundaries for a gapped or nodal bulk. Since, in our case, we have surface
modes localized both on the edges and corners of the lattice, we extract the
energy spectra for different boundary conditions. For the bulk spectrum we
apply fully periodic boundary conditions and generally sample the Brillouin
zone taking paths connecting the high-symmetry points of the square lattice
$\Gamma=(0,0)$, $X=(\pi/a,0)$, $Y=(0,\pi/a)$, and $M=(\pi/a,\pi/a)$. For edge
or corner states we apply open boundary conditions in the appropriate
directions. For the midgap states we plot both the wavefunction
$\Psi_{m}(x,y)$, to verify in which part of the boundary these states are
localized, and the Majorana polarization to verify whether these states are
MZMs or not. Since the Hamiltonian is particle-hole symmetric, there is a
correspondence between states with infinitesimally small positive energy
($\epsilon=0^{+}$) and those with infinitesimally small negative energy
($\epsilon=0^{-}$). Therefore, when showing the localization and Majorana
polarization properties of the midgap modes of the HOTSC and hybrid phases,
we, for convenience, focus on the states with energy $\epsilon=0^{-}$. To
obtain complete information on the localization of all low-lying states in the
system, we also show the local density of states (LDOS) at zero frequency
$\omega=0$.

##  IV Topological phases

In this section we present our main analysis of the topological phases of the
model in Eq. (1) using the topological invariants discussed in Sect.III. We
focus on the deep topological limit of the BBH model, setting $\lambda=0$ for
simplicity, which allow us to obtain analytical expressions of the phase
boundaries, but we remark that the topological phases are present for
$|\lambda/t|<1$. We refer to Appendix A for results on finite $\lambda$. For
$\lambda=0$, the topological phase diagram is displayed in Fig. 1(b). Here, we
start by detailing how we obtain the phase boundaries, followed by a detailed
description of each of the phases.

###  IV.1 Phase boundaries

A quantum phase transition is accompanied by the closing of the gap of the
system [79]. Therefore, the topological phase boundaries are obtained by
considering the energy gap, $\delta$ between the highest valence band and the
lowest conduction band around zero energy. In Fig. 2 we plot $\delta$ as a
function of the superconducting order parameter $\Delta$ and magnetic field
$B_{x}$ for different boundary conditions. Starting with fully periodic
boundary conditions in Fig. 2(a), we obtain that the bulk gap closes along
several different lines. For $B_{x}=\sqrt{2t^{2}+\Delta(\Delta-2t)}$ (magenta
dashed line) the bulk spectrum closes at momenta $(k_{x},k_{y})=(0,\pi/2)$ and
$(\pi,\pi/2)$ and for $\Delta=\sqrt{2t^{2}+B_{x}(B_{x}-2t)}$ (yellow dashed
line) it closes at momentum $(\pi/2,\pi/2)$. These two lines separates the
nodal flat bands phase (blue region in Fig. 1(b)) from the nodeless flat bands
phase (red region in Fig. 1(b)) and hybrid phase. Further, the bulk spectrum
also closes at $B_{x}=\sqrt{2t^{2}+\Delta(\Delta+2t)}$ (green dashed line) at
momentum $(k_{x},k_{y})=(\pi,\pi/2)$ and at
$\Delta=\sqrt{2t^{2}+B_{x}(B_{x}+2t)}$ (cyan dashed line) at momentum
$(\pi/2,\pi/2)$, which separates the trivial phase from the nodal flat band
phase. These results also establish that the HOTSC (yellow in Fig. 1(b)),
hybrid (green), trivial (white), and nodeless flat band (red) phases are all
fully gapped in the bulk.

To complement the results for fully periodic boundary conditions, we also
analyze how the gap closes for open boundary conditions along $y$ in Fig. 2(b)
and along $x$ in Fig. 2(c). This analysis of different boundary conditions
brings three important pieces of information. First, we notice that there is a
new gap-closing line at $\Delta=\sqrt{2}t-B_{x}$ (white dashed line), which
indicates the phase boundary between the HOTSC and hybrid phases in Fig. 1(b).
This line just appears for $x$-open boundary conditions, see Fig. 2(c),
illustrating how the gap only closes along the $y$-direction. Already this
result points to the higher order topology of its neighboring phases. Second,
we notice that the nodal flat bands phase (blue in Fig 1(b)) actually
seemingly hosts a small but finite gap, seen as red arcs in Figs. 2(a,c), but
it does become completely gapless in Fig. 2(b). In fact, for all values of
$B_{x}$ and $\Delta$ in this phase, the gap closes at $k_{y}=\pi/2$ and
different $k_{x}$. Therefore, the spectrum is actually gapless in all three
figures and the red arcs are just finite-size effects due to a necessarily
finite $k_{x}$-point sampling in Figs. 2(a,c). This verifies that the nodal
flat bands phase is a bulk nodal phase. Finally, we find that one of the bulk
gapped regions, the nodeless flat band phase (red in Fig. 1(b)) is also gapped
for $x$-open boundary conditions but notably not for $y$-open boundary
conditions, as seen in Fig. 2(b). This indicates zero-energy states localized
to edges along the $x$-direction.

Figure 2: Energy gap $\delta$ between the highest valence band and lowest
conduction band around zero energy, derived from Eq. (2), as a function of
superconducting pair amplitude $\Delta/t$ and in-plane magnetic field$B_{x}/t$
for (a) fully-periodic, (b) $x$-periodic (open in $y$) and (c) $y$-periodic
(open in $x$) boundary conditions. Parameters used: $\lambda=0$, system size
$32$ unit cells in each direction.

In the next subsections, we detail the properties of each of the non-trivial
topological phases, while the trivial phase is discussed in Appendix B. We
characterize the general properties considering both the energy spectra
extracted above and the topological invariants discussed in Sec. III, for
representative values of $B_{x}$ and $\Delta$ in all phases.

###  IV.2 HOTSC phase

For small $B_{x}$ and $\Delta$, we find a phase that we call the HOTSC phase,
indicated by the yellow region in Fig. 1(b). The main features of this phase
are summarized in Fig. 3. First, considering the energy $\epsilon$ spectrum
for fully periodic in Fig. 3(a), $x$-periodic in Fig. 3(b), $y$-periodic in
Fig. 3(c), and fully open in Fig. 3(d) boundary conditions, we realize that
the bulk of the system is gapped, while there are eight states at zero energy
present only for fully open boundary conditions, Fig. 3(d). The system is thus
gapped under $x$\- or $y$-periodic boundary conditions in this phase. Plotting
the LDOS at zero frequency/energy in Fig. 3(e), we see that these states are
corner states, explaining why they appear just for open boundary conditions.
To further understand the properties of these zero energy modes, we examine
the Wannier spectrum $\nu$ along both $x$ in Fig. 3(f) and along $y$ in Fig.
3(g). The half-integer values of both $\nu_{x,y}$ indicate the higher-order
character of this phase. We note, however, that within increasing value of
$\Delta$ and $B$, but still in the HOTSC phase, the midgap values of $\nu^{x}$
deviate significantly from $0.5$ due to the breaking of mirror symmetry along
$x$. Nevertheless, we still have robust corner modes with properties similar
to the ones discussed below. The higher-order aspect is corroborated by the
entanglement spectrum $\xi$ in Fig. 3(h), which shows distinct isolated half-
quantized modes inside the midgap region.

Figure 3: Main features of the HOTSC phase. Energy spectrum $\epsilon$ for
fully periodic (a), $x$-periodic (b), $y$-periodic (c), and open (d) boundary
conditions. Inset in (d) a zoom-in on the modes in the middle of the spectrum
with energy $\epsilon=0$. LDOS at zero energy (e). Wannier spectrum along $x$
(f) and $y$ (g). Entanglement spectrum $\xi$ (h). Blue dashed lines mark the
values of 0.25 (1/4) and 0.75 (3/4). Parameters used: $B_{x}=0.2t$,
$\Delta=0.1t$, $\lambda=0$, system size $32$ unit cells for directions with
open boundary conditions and $100$ $k-$points.

Finding a total of eight corner modes, two at each corner is expected
considering the BBH model at $B_{x}=\Delta=0$ [80, 38]. Due to the finite
$\Delta$, however, the corner modes are now a combination of particles and
holes as they are Bogoliubov quasiparticles in a superconductor. In Figs.
4(a-d) and (e-h) we analyze, respectively, the zero-energy Boguliubov
quasiparticles by plotting the wavefunction $\Psi$ from in Eq. (8) and
Majorana polarization $P$ from Eq. (7)) of the four states at zero (negative)
energy. As seen, these modes are fully localized at the corners, as shown in
both $\Psi$ and $P$. Moreover, since $P$ is much smaller than $\Psi$, not
being even visible in this scale, and since there are two states per corner,
we designate these states as ABSs.

Figure 4: HOTSC boundary states features. (a-d) Wavefunction probability
density $\Psi$ and (e-h) Majorana polarization $P$ for the four states with
$\epsilon=0^{-}$ in the HOTSC phase. Parameters same as in Fig. 3.

###  IV.3 Hybrid phase

For increasing values of $\Delta$ and $B_{x}$, the system enters into the
hybrid phase, green region in Fig. 1(b). While present only in a narrow region
of the phase diagram, it displays interesting features.

We analyze the general properties of this phase in Fig. 5. We find that it has
energy spectrum similar to the HOTSC phase: both the bulk, Fig. 5(a), and the
edges, Figs. 5(b,c), are gapped while there are eight midgap states at zero
energy for open boundary conditions, Fig. 5(d). These midgap states are still
located at the corners, as shown in the LDOS at zero energy, Fig. 5(e).
However, when inspecting the Wannier spectrum, Figs. 5(f,g), we see that
$\nu^{y}$ is half-quantized, while $\nu^{x}$ exhibits a gap around $0.5$.
Further, the entanglement spectrum in Fig. 5(h) shows a continuous array of
midgap eigenvalues that are present symmetrically around $\xi=0.5$. Therefore,
this hybrid-order phase can be considered topologically distinct from the
previous HOTSC phase. We remark that since $\nu^{x}$ is not quantized for
larger values of $\Delta$ and $B$ in the HOTSC phase, the topological phase
transition between the two phases is due to the change of localization of some
of the zero energy states, as shown below.

Figure 5: Main features of the hybrid phase. Same panels as in Fig. 3.
Parameters used: $B_{x}=0.9t$, $\Delta=0.8t$, $\lambda=0$, system size $32$
unit cells for directions with open boundary conditions and $100$ $k-$points.

By going back to Figs. 2(a,b), we see that there is no change in the energy
gap from the phases HOTSC and hybrid for systems with periodic conditions
along $x$, indicating that the change from the phase HOTSC to the hybrid phase
is due to properties related to the I and III edges, which are along the
$y$-direction. As a consequence, only when periodic boundary conditions are
applied along the $y$-direction, the gap closing, indicating a topological
phase transition, is noticed between the HOTSC and hybrid phases.

Continuing to analyze the wavefunctions in Figs. 6(a-d) and Majorana
polarization in Figs. 6(e-h), we find that half, Figs. 6(a,b), of the zero-
energy states with $\epsilon=0^{-}$ have a significant weight not just at the
corners but also partially along the $y$-direction, even better visualized in
the Majorana polarization, Figs. 6(e,f), while the other half of the
$\epsilon=0^{-}$ states, Figs. 6(c,d), are still localized at the corners of
the lattice, as in the HOTSC phase. As a consequence, we find a combination of
corner modes, as appropriate for a second-order topological phase and clearly
a remnant of the HOTSC second-order phase, and more extended edge-localized
modes. The latter is also expected for a first-order topological phase, but
here we note that we only have a finite number of modes localized at the edge,
four to be precise, which do not grow with system size, and thus the modes are
not representable of a standard first-order topological phase. Instead, the
edge modes are those of a dipolar phase [12, 81], which can be seen as created
as the edge of a single 1D system extending in the $y$-direction. Due to this
combination of boundary modes, we name this a hybrid-order phase, being an
interesting mix of different topologies. However, note that our use of the
word hybrid should not be confused with the situation where more standard
first- and second-order phases are appearing at the same time, also called a
hybrid phase [82, 83, 84, 85, 86]. Here, the hybrid phase is instead a
standard second-order phase appearing jointly with a dipolar phase.

Moreover, with all eight zero-energy modes now spatially distinct and also
presenting similar values for $\Psi$ and $\left|P\right|$ as seen in Fig. 6,
we find them to be fully particle-hole symmetric. Therefore, the increase of
$\Delta$ and $B_{x}$ moving from the HOTSC to the hybrid phase turn the
surviving four corner modes in the HOTSC phase corner-localized, single, MZMs
and the other four into edge-extended MZMs.

Figure 6: Hybrid phase boundary states features. Same panels as in Fig. 4 for
the four states with $\epsilon=0^{-}$ in the hybrid phase. Parameters same as
in Fig. 5.

###  IV.4 Nodal flat bands phase

Increasing either $\Delta$ or $B_{x}$, we arrive at the nodal flat band phase,
represented in blue in Fig. 1(b). This phase is characterized by bulk nodal
points, i.e. there exists momenta where the bulk energy gap (superconducting
gap) closes. Between these bulk nodal points, we find flat bands that appear
localized to the boundaries of the system. We investigate the general
properties of this phase in Fig. 7. Considering the energy dispersion along
the high-symmetry points in Fig. 7(a), the bulk of the system seems to be
gapped, in contradiction to Fig. 2(a). However, analyzing the energy spectrum
for the system with open boundary conditions along $y$, in Fig. 7(b) and $x$
in Fig. 7(c), we understand that this happens because the bulk gap is only
zero at specific points $k_{x}a\neq 0,\pi$. The bulk gap vanishes away from
the high-symmetry line and that is why gap-closing in not visible in Fig.
7(a). Additionally, there exist zero-energy flat bands along $k_{x}$, which
connect these bulk nodal points, while all bands are dispersive along $k_{y}$.
This suggests that the flat bands are boundary states localized along the
edges II and IV in Fig. 1(a). These flat bands cause the energy spectrum with
fully open boundary conditions in Fig. 7(d) to host a macroscopic degeneracy
at zero energy, as shown in the inset of Fig. 7(d).

Figure 7: Main features of the nodal flat band phase. Same panels as in Fig.
3. Parameters used: $B_{x}=1.5t$, $\Delta=0.9t$, $\lambda=0$, system size $32$
unit cells for directions with open boundary conditions and $100$ $k-$points.

Further evidence that the zero-energy states are localized along the II and IV
edges is shown in the LDOS at zero frequency in Fig. 7(e), which shows weight
just at these edges. This indicates that these are the boundary states of a
first-order topological phase, also since the absence of corner modes clearly
discards these phase as any higher-order phase. We find this fully consistent
with the topological invariants: an absence of midgap modes in $\nu^{x}$ in
Fig. 7(f), while a number of 0.5 modes in $\nu^{y}$ appears in Fig. 7(g). We
further find that the number of zero-energy states and 0.5 modes in $\nu^{y}$
scales with the number of unit cells along $x$, which further corroborates
that this is an edge phenomenon along the $x$-direction. Here $\nu^{y}$ is
thus an invariant that characterizes the presence of these zero-energy states,
which are related by a bulk-boundary correspondence to the nodes in the bulk
[2, 55, 57]. Another indication that this is not a second-order phase is the
entanglement spectrum in Fig. 7, which does not show sharply quantized modes
at 0.5. We instead find a discontinuous array of mid-gap eigenvalues
symmetrically placed around 0.5.

Figure 8: Nodal flat band phase boundary states features. Same panels as in
Fig. 4 for the four states with $\epsilon=0$ for $k_{x}=0$ in the nodal flat
band phase, using a partial Fourier transform, as explained in the main text.
Parameters same as in Fig. 7.

To better understand the overall properties of the zero-energy flat bands
occurring between the bulk nodal points, we plot in Fig. 8 the wavefunction
$\Psi$ in Figs. 8 (a-d) and the corresponding Majorana polarization $P$ in
Figs. 8 (e-h), for four of the zero-energy modes. It is here most convenient
to consider a system with open boundary conditions along $y$ and periodic
boundary conditions along $x$ (corresponding to the spectrum of Fig. 7(b)) for
$k_{x}=0$. After diagonalizing we then have a wavefunction
$\tilde{\psi}(k_{x}=0,y)$ and to obtain a complete real space wavefunction, we
multiply

$\psi_{k_{x}=0}(x,y)=1/\sqrt{L}\exp(ik_{x}x)\tilde{\psi}(k_{x}=0,y).$

This is the wavefunction we use to compute the probability density and
Majorana polarization in Figs. 8 and we find four-zero energy modes for
$k_{x}=0$. We also find four-zero energy modes for each other value of
momentum. The Majorana polarization $P$ has the same maximum value as the
wavefunction probability $\Psi$ and also changes sign between the two edges
for all modes, indicating that they resemble that of MZMs. Still, since there
is a degeneracy of two zero-energy modes per edge per momentum in this phase,
we cannot for sure classify these states as MZMs without considering if a
linear combination of them cannot recombine into complex fermionic modes.
However, since these modes present a strong Majorana polarization and we have
two possible particle-hole symmetric pairs that we can build with spin and
orbital degrees of freedom, we choose to still call these states flat band
MZMs. We note that these flat band MZMs are extended along the edges, and
their number increases with system size, in contrast to the MZMs obtained for
the hybrid phase that are always just one per edge and one per corner, as
discussed in Sec. IV.3.

###  IV.5 Nodeless flat bands phase

Finally, moving closer to the diagonal of the phase diagram, with
$\Delta/t\approx B_{x}/t>1$, we find the nodeless flat band phase in red in
Fig. 1(b). For this phase, the bulk system is gapped, as seen in Fig. 9(a), as
is the semi-infinite spectrum with the open boundary along $x$, illustrated in
Fig. 9 (c). In Fig. 9(b), we instead find zero-energy flat bands spanning the
whole Brillouin zone along $k_{x}$. Fully open boundary conditions also
generate zero-energy mid-gap states in Fig. 9(d). Plotting the zero-energy
LDOS we find that these flat bands are localized on the II and IV edges in
Fig. 9(e). Further, although $\nu^{x}$ is completely gapped in Fig. 9(f),
$\nu^{y}$ in Fig. 9(g) has a microscopic number of $0.5$ modes, which grows
with the number of unit cells along $x$. This $\nu^{y}$ profile confirms that
the zero-energy states are an edge phenomenon along the $x$-direction.
Finally, the entanglement spectrum in Fig. 9(h) shows a discontinuous array of
mid-gap eigenvalues around 0.5 while they gather more closely to 0.5. This is
qualitatively similar to the nodal flat band case, however, the discontinuous
profile is more asymmetric around 0.5 in this case compared to the previous
case.

Figure 9: Main features of the phase with nodeless flat bands. Same panels as
in Fig. 3. Parameters used: $B_{x}=1.8t$, $\Delta=2t$, $\lambda=0$, system
size $32$ unit cells for directions with open boundary conditions and $100$
$k-$points. Figure 10: Nodeless flat band boundary states features Same panels
as in Fig. 4 for the four states with $\epsilon=0$ for $k_{x}=0$ on the
nodeless flat band phase, using a partial Fourier transform, as explained in
the main text. Parameters same as in Fig. 9.

To investigate the spatial properties of the states at zero energy in this
phase, shown in Fig. 10, we focus on the semi-periodic system with $k_{x}=0$
(corresponding to the spectrum of Fig. 9 (b)) as done for the nodal phase.
Although the modes present a similar spatial profile, we remark that these
flat bands are present without corresponding nodal points in the bulk. As
such, they are a very different kind of surface state compared to those of the
nodal flat band phase. We currently do not know how they can be fully
classified in terms of symmetry-protected topological phases. Nevertheless,
they are still topologically protected by symmetry since the number of states
at zero energy is determined by the number of $0.5$ states in $\nu^{y}$, which
is a topological property. Finally, we note that these flat bands states again
present a Majorana polarization $P$ of the same amplitude as the of the
wavefunction density $\Psi$, Figs. 10(e-h), which makes them designated as
flat band MZMs in a similar way to in the nodal flat band case.

##  V Analytical results

After having described in detail the properties of all the different
topological phases numerically, we also present some results feasible to
achieve analytically in order to enhance the overall understanding. While a
complete analytical treatment at present seems not possible, we can still
obtain selective analytical results. Here, we first present analytical
calculations supporting the existence of a nodal phase despite an $s$-wave
superconducting order parameter. For this purpose, we define a dual
Hamiltonian. In addition to a dual Hamiltonian, we also extract a low-energy
Hamiltonian and are able to show analytically that corner states exist in its
regime of validity.

###  V.1 Nodal superconductivity

One of the phases, the nodal flat band phase, depicted in blue in Fig. 1(b)
and discussed in detail in Sec. IV.4, hosts bulk nodes. This is despite the
model only containing on-site, or in $k$-space isotropic $s$-wave, pairing,
which are usually associated with a fully gap, while nodes are usually thought
of as requiring a $k$-dependent order parameter. The presence of nodes in the
nodal flat band phase can be understood as a consequence of the spin-orbit
coupling in the normal state Hamiltonian, as such terms can induce an
effective momentum-dependent pairing. The equivalence between the spin-orbit
coupling and nodal superconductivity can be seen explicitly by performing the
unitary transformation of the Hamiltonian of Eq. (2) (for $\lambda=0$):

$h_{D}(\mathbf{k})=Uh_{\rm
BdG}(\mathbf{k})U^{\dagger};\,\,U=\begin{pmatrix}\sigma_{0}\otimes
s_{0}&i\sigma_{2}\otimes s_{2}\\\ i\sigma_{2}\otimes s_{2}&\sigma_{0}\otimes
s_{0}\end{pmatrix},$ (9)

where we obtain a dual Hamiltonian $h_{D}$, following the nomenclature of Ref.
[67], with the form

$\displaystyle h_{D}(\mathbf{k})$ $\displaystyle=$
$\displaystyle\tau_{1}\big{[}-t\sin(k_{y}a)\sigma_{0}s_{3}-t\cos(k_{x}a)\sigma_{0}s_{1}$
(10) $\displaystyle+$ $\displaystyle
B_{x}\sigma_{2}s_{3}\big{]}+\Delta\tau_{3}\sigma_{1}s_{0}-t\cos(k_{y})\tau_{0}\sigma_{2}s_{2},$

where we omit the outer product symbol $\otimes$ for convenience.

In the dual Hamiltonian in Eq. (10), the superconducting pairing appears in
multiple terms. There is pairing with $p$-wave intra-orbital spin-triplet
pairing symmetry given by $-t\sin(k_{y}a)\tau_{1}\otimes\sigma_{0}\otimes
s_{3}$, with extended $s$-wave intra-orbital spin-singlet symmetry given by
$-t\cos(k_{x}a)\tau_{1}\otimes\sigma_{0}\otimes s_{1}$, and with $s$-wave odd-
interorbital spin-triplet symmetry given by
$B_{x}\tau_{1}\otimes\sigma_{2}\otimes s_{3}$. The first two terms always
present nodes at some $k_{x}$ and $k_{y}$ coordinates and result in the nodal
profile of the superconducting order parameter. If such nodes in the
superconducting order parameter also overlap with the normal-state Fermi
surface, then the system will have a nodal energy gap. In the presence of the
third term $B_{x}\neq 0$, the nodal structure caused by the first two terms
can still be retained with nodes then appearing at different values of $k_{x}$
and $k_{x}$ in appropriate parameter regimes. This qualitatively explains why
a nodal state at all can be possible and also shows how it is intricately
linked to the $\Delta$ and $B_{x}$ parameters. We further note that in the
pairing terms, the dependence on $k_{x}$ comes from a hopping term
$\cos(k_{x}a)$, while $\sin(k_{y}a)$ appears due to the spin-orbit coupling of
the normal state. This difference, together with the fact that the magnetic
field is applied along $x$, explains the asymmetry seen in the properties of
the system with periodic boundary conditions along $x$ or $y$. Finally, the
existence of a $p$-wave term is also an underlying reason for the existence of
MZMs in several of the different phases.

###  V.2 Corner states

In addition to the dual Hamiltonian, we can also analytically study a low-
energy continuum Hamiltonian. This is a technique commonly used to connect the
presence of localized boundary states with a non-trivial mass profile, as
first constructed by Jackiw and Rebbi[87], and where the mass profile tells us
some properties of the topological phase. For instance, for a higher-order
topological system, the mass profile needs to change sign between two adjacent
edges [8, 35], with the corner in between hosting the localized boundary
state. To obtain the most simple low-energy Hamiltonian, we keep just terms
that are first-order in momentum in $h_{\rm BdG}$ in Eq. (1) and obtain

$h_{\Gamma}(\mathbf{k})=tk_{y}a\Gamma_{1}+t\Gamma_{2}+tk_{x}a\Gamma_{3}+t\Gamma_{4}+\Delta\Gamma_{5}+B_{x}\Gamma_{6}.$
(11)

Here $\mathbf{k}=(k_{x},k_{y})$ represents the continuum momentum in relation
to the $\Gamma$ point.

We can obtain the low-energy descriptions of any corner states of $h_{\Gamma}$
by substituting $k_{x/y}\rightarrow-i\partial_{x/y}$ and looking for localized
solutions. The complete calculation is reported in detail in Appendix C, and
here we comment on the results. Considering appropriate boundary conditions,
the solutions of $h_{\Gamma}$ localized on the I and III edges are

$\psi^{\textrm{I}/\textrm{III}}(x,y)=\sum\limits_{l,m=\pm
1}c_{l,m}\mathcal{N}_{x}e^{\mp\sqrt{2}\frac{x}{a}}e^{ik_{y}y}\chi_{l,m}^{\textrm{I}/\textrm{III}},$
(12)

where $\mathcal{N}_{x}$ is a normalization factor,
$\sum\left|c_{l,m}\right|^{2}=1$, and the spinors
$\chi_{l,m}^{\textrm{I}/\textrm{III}}$ are given by

$\chi_{l,m}^{\textrm{I}}=\ket{l,m}\frac{\ket{s=1}+\left(m-\sqrt{2}\right)\ket{s=-1}}{\sqrt{4-2m\sqrt{2}}},$
(13)

and

$\chi_{l,m}^{\textrm{III}}=\ket{l,m}\frac{\ket{s=1}-\left(m+\sqrt{2}\right)\ket{s=-1}}{\sqrt{4+2m\sqrt{2}}},$
(14)

where $l=\pm 1$, $m=\pm 1$, and $s=\pm 1$ are eigenvalues of $\tau_{3}$,
$\sigma_{3}$ and $s_{3}$, respectively. We further obtain the edge
Hamiltonians by projecting the low-energy Hamiltonian in the subspace of Eqs.
(13) and (14), yielding

$\displaystyle h^{\textrm{I}}(k_{y})$ $\displaystyle=$ $\displaystyle-
tk_{y}a\,\tau_{3}\otimes\sigma_{2}-\frac{B_{x}}{\sqrt{2}}\,\tau_{3}\otimes\sigma_{0},$
(15) $\displaystyle h^{\textrm{III}}(k_{y})$ $\displaystyle=$ $\displaystyle-
tk_{y}a\,\tau_{3}\otimes\sigma_{2}-\frac{B_{x}}{\sqrt{2}}\,\tau_{3}\otimes\sigma_{0}.$
(16)

These are Dirac-like Hamiltonians with a mass term proportional to $B_{x}$.
Note that both solutions in Eqs. (13) and (14) are polarized in the effective
particle-hole space, represented by $\tau_{3}$, while they are a combination
of different spin states with the effective spin degrees of freedom
represented by $\sigma_{3}$. Note here that since the edge Hamiltonians are
projected Hamiltonians, the original notion of particle-hole or spin is
replaced by their effective notions. Turning to the localized solutions on the
II and IV edges, we obtain similarly as on the other edges

$\psi^{\textrm{II}/\textrm{IV}}(x,y)=\sum\limits_{l,m=\pm
1}\tilde{c}_{l,m}\mathcal{N}_{y}e^{ik_{x}x}e^{\mp\sqrt{2}\frac{y}{a}}\chi_{l,m}^{\textrm{II}/\textrm{IV}},$
(17)

where $\mathcal{N}_{y}$ is a normalization factor, we have
$\sum\left|\tilde{c}_{l,m}\right|^{2}=1$, and the spinors
$\chi_{lm}^{\textrm{II}/\textrm{IV}}$ are given by

$\chi_{l,m}^{\textrm{II}}=\ket{l,m}\frac{\ket{s=1}-m\left(1+\sqrt{2}\right)\ket{s=-1}}{\sqrt{4+2\sqrt{2}}},$
(18)

and

$\chi_{l,m}^{\textrm{IV}}=\ket{l,m}\frac{\ket{s=1}+m\left(1-\sqrt{2}\right)\ket{s=-1}}{\sqrt{4+2\sqrt{2}}},$
(19)

where again $l=\pm 1$, $m=\pm 1$, and $s=\pm 1$ are eigenvalues of $\tau_{3}$,
$\sigma_{3}$ and $s_{3}$, respectively.

Again projecting on the solutions of Eqs. (18) and (19), we obtain the
effective Dirac Hamiltonians

$\displaystyle h^{\textrm{II}}(k_{x})$ $\displaystyle=$ $\displaystyle
tk_{x}a\,\tau_{3}\otimes\sigma_{2}-\frac{1}{\sqrt{2}}B_{x}\,\tau_{3}\otimes\sigma_{3},$
(20) $\displaystyle h^{\textrm{IV}}(k_{x})$ $\displaystyle=$ $\displaystyle
tk_{x}a\,\tau_{3}\otimes\sigma_{2}-\frac{1}{\sqrt{2}}B_{x}\,\tau_{3}\otimes\sigma_{3}.$
(21)

Having extracted all four edge Hamiltonians, Eqs. (15), (16), (20), and (21),
we note that all four edges have a mass term proportional to $B_{x}$. On the I
and III edges, it is proportional to $\tau_{3}\otimes\sigma_{0}$, while on the
II and IV edges it is proportional to $\tau_{3}\otimes\sigma_{3}$. Since
eigenvalues of $\sigma_{0}$ and $\sigma_{3}$ can be different in terms of
their signs, there exists a sign change in the mass term for adjacent edges.
For a given subspace defined by $\tau_{3}$, we can hence obtain a positive
mass $B_{x}$ on edges II and IV, while having a negative mass $-B_{x}$ on
edges I and III, when $\sigma_{3}$ is projected on the effective spin-down
states. This leads to the formation of corner modes between adjacent edges. On
the other hand, for the effective spin-up projections, all the edges have
negative mass terms given by $-B_{x}$, and hence corner modes are not expected
to appear. The effective spin-polarized nature of the corner modes can be
attributed to the sign-polarized Majorana polarization for a given corner
mode. This analysis qualitatively explains the emergence of corner states for
our model, Eq. (2), in both the HOTSC and hybrid phases.

##  VI Conclusions

In this work, we investigate the topological phases of an orbital BBH model
proximitized to a conventional spin-singlet $s$-wave superconductor and in the
presence of an in-plane magnetic field. The interpretation of the BBH model in
terms of orbital and spin degrees of freedom makes the superconducting pairing
have a matrix structure not present in the original dimerized lattice of the
BBH model and allows for an intriguingly rich topological phase diagram,
despite using only a conventional superconductor. We map out the phase diagram
by considering different boundary conditions and investigate the topology of
the individual phases by calculating both the Wannier and entanglement
spectra, as well as the Majorana polarization.

At weak superconducting pairing $\Delta$ and magnetic field $B_{x}$, we find a
HOTSC phase (yellow in Fig. 1(b)) with eight zero-energy corner modes. This
phase can be seen as analogous to the standard second-order topological phase
in the BBH model but now with particle and hole degrees of freedom since the
system is superconducting. Beyond this expected HOTSC phase, we also find
several other, unexpected, topological phases. First, an unusual hybrid phase
(green in Fig. 1(b)) presents an atypical mix between the second-order and a
first-order topological phase. Here four zero-energy corner states from the
HOTSC phase are preserved, while a dipolar phase contributes another four
zero-energy edge states, all experiencing MZM character. Notably, the number
of edge states does not grow with system size in the hybrid phase despite its
mixed-order character. Second, two additional first-order phases present
symmetry-protected zero-energy flat bands on opposite edges, with either nodal
or nodeless bulk dispersion. The nodal flat bands phase (blue in Fig. 1(b))
presents flat bands MZMs that are straightforwardly connected by a bulk-
boundary correspondence to the bulk nodes. However, in the nodeless flat bands
phase (red in Fig. 1(b)) the bulk remains fully gapped and the flat bands MZMs
now spans the whole edge Brillouin zone and are protected by a quantized
Wannier spectrum. These results not only establish the rich phase diagram of a
superconducting BBH model, but also demonstrate the existence of both unusual
hybrid mixing of topological phases and nodeless flat bands.

Note added: During the final preparation of this work we became aware of the
recent work Ref. [88], where zero-energy flat bands are also found to appear
as topological boundary states in a nodeless superconductor. However, in that
case the flat band is protected by a bulk invariant in 3D. Our flat band
instead appears in a 2D system and is protected by a 1D edge invariant. It
remains to be studied if deeper similarities between the two exist.

###### Acknowledgements.

We are grateful for discussions with A. K. Ghosh, A. Saha, S. Mondal, B. Roy,
V. Juricic, O. A. Awoga, D. Chakraborty, and A. Bhattacharya. We acknowledge
financial support from the Knut and Alice Wallenberg Foundation. Part of the
simulations were enabled by resources provided by the National Academic
Infrastructure for Supercomputing in Sweden (NAISS) and the Swedish National
Infrastructure for Computing (SNIC) at UPPMAX, partially funded by the Swedish
Research Council through grant agreements no. 2022-06725 and no. 2018-05973.

##  Appendix A Finite $\lambda$

In the main text, we focus on the properties of Eq. (2) for $\lambda=0$. Here,
we discuss the effect of having a finite $\lambda$ in the phases discussed in
the main text. For large $\lambda>t$, the system is just a trivial
superconductor. For smaller but finite $\lambda$, the phase boundaries of the
system change considerably but the topological phases largely persist. In Fig.
11 we set $\lambda=0.5t$ and we repeat the plot of the energy gap $\delta$ as
a function of $B_{x}$ and $\Delta$ from Fig. 2. We still retain the phases
discussed in the main text, although the phase boundaries change considerably
and now cannot be computed analytically. In addition, we obtain a new phase,
which we call the dipolar phase, described next.

Figure 11: Energy gap $\delta$ derived from Eq. (2) as a function of
$\Delta/t$ and $B_{x}/t$ for (a) fully-periodic, (b) $x$-periodic (open in
$y$) and (c) $y-$periodic (open in $x$) boundary conditions. Parameters used:
$\lambda=0.5t$, system size $32$ unit cells in each direction.

###  A.1 Dipolar phase

With $\lambda$ finite, but smaller than $t$, we find a different kind of
first-order topological phase appearing for higher values of either $\Delta$
and $B_{x}$, when one is much larger than the other. This phase, as shown in
Fig. 12(a-d), has a gapped bulk but with some localized midgap states on the I
and III edges, as shown in the zero-energy LDOS in Fig. 12(e). The $\nu^{x}$
invariant along $x$ is not quantized, see Fig. 12(f), while the $\nu^{y}$
presents two half-quantized values, see Fig. 12(g), indicating that there are
four symmetry-protected modes. The entanglement spectrum shows a symmetric
discontinuous distribution of eigenvalues excluding 0.5 within the mid gap
region as depicted in Fig. 12(h).

Figure 12: Main features of the dipolar phase. Same panels as in Fig. 3.
Parameters used: $B_{x}=1.5t$, $\Delta=0.2t$, $\lambda=0.5$, system size $32$
unit cells for directions with open boundary conditions and $100$ $k-$points.
I am not sure whether we are interchanging between $\nu_{x}$ and $\nu_{y}$.

In Fig. 13, we plot the wavefunction probability $\Psi$ in Fig. 13(a-d) and
the Majorana polarization $P$ in Fig. 13(e-h), for the four states with
negative energy in the inset of Fig. 12(d), just as for all other phases
investigated earlier. The states with finite energy, seen Figs. 13(a,b), are
identically localized on the I and III edges, respectively, and their $P$
oscillate around zero along the edge for both edges, see Figs. 13(e,f). In
contrast, the states with zero energies are localized more in the corner of
these edges, see Figs. 13(c,d), and $P$ change sign only once along the edges
on the lattice, see Figs. 13(g,h). These zero-energy states, therefore,
resemble that of the edge states in the hybrid phase discussed in Section
IV.3, while the corner states of the hybrid phase are notably absent here. Due
to the similarity with the dipolar edge states in the hybrid phase, we call
this a dipolar phase. We note that the $P$ of the zero-energy boundary states
appropriately also resembles an electric dipole.

Figure 13: Dipolar phase boundary states features. Same panels as Fig. 4 for
the four states with smallest (negative energy), but where only the last two
states (c,d,g,h) have energy $\epsilon=0^{-}$, in the dipolar phase.
Parameters same as in Fig. 12.

##  Appendix B Trivial phase

For high values of only $\Delta$ or $B_{x}$, we find a topologically trivial
phase. For $\lambda=0$, this occurs for
$\Delta>\sqrt{2t^{2}+B_{x}\left(2t+B_{x}\right)}$ or
$\Delta<-t+\sqrt{B_{x}^{2}-t^{2}}$, corresponding to the white region in Fig.
1(b). Although this phase is topologically trivial, for completeness, we
display the main features of this phase in Fig. 14. The system is gapped for
all boundary conditions, see Figs. 14(a-d), which is supported by the LDOS at
zero frequency, see Fig 14(e), which shows no occupation. Further, all
topological invariants indicate a trivial system, see Figs. 14(f-h).

Figure 14: Main features of the trivial phase. Same panels as in Fig. 3.
Parameters used: $B_{x}=2t$, $\Delta=0.5t$, $\lambda=0$, system size $32$ unit
cells for directions with open boundary conditions and $100$ $k-$points.

##  Appendix C Derivation of edge Hamiltonians

Here we provide the detailed derivation of the edge Hamiltonian expressions
used in the main text. We continue in the deep topological limit of the BBH
model with $\lambda=0$. Since the continuum Hamiltonian Eq. (11) is composed
of linear operators in $x$, $y$, and matrices acting on the particle-hole
$\tau$, orbital $\sigma$, and spin $s$ degrees of freedom, we can write a
generic wavefunction as

$\psi(x,y)=\sum\limits_{\alpha}c_{\alpha}\phi^{\alpha}(x)\varphi^{\alpha}(y)\chi_{\alpha},$
(22)

where $\phi^{\alpha}$ ($\varphi^{\alpha}$) is a complex spatial function of
$x$ ($y$), and $\chi_{\alpha}$ are eight component spinors labelled by
$\alpha$. Using this form, we can obtain solutions that are localized on
specific edges (or corners) by checking how the Hamiltonian acts individually
in each part.

For the I and III edges, we look for localized solutions on $x$. We split
$h_{\Gamma}=h_{0x}+h_{ky}$, where
$h_{0x}=t\Gamma_{2}+tk_{x}a\Gamma_{3}+t\Gamma_{4}$ is the part which determine
the zero states localized along $x$, while we treat
$h_{ky}=tk_{y}a\Gamma_{1}+\Delta\Gamma_{5}+B_{x}\Gamma_{6}$ as its
perturbation, reasonable given the fact that $B_{x},\Delta\ll t$. Next we
perform the substitution $k_{x}\rightarrow-i\partial_{x}$ and look for
solutions of

$h_{0x}\psi(x,y)=\left[t\Gamma_{2}-it\partial_{x}a\Gamma_{3}+t\Gamma_{4}\right]\psi(x,y)=0.$
(23)

Using the ansatz $\phi_{\alpha}(x)=\exp(q_{x}x/a)$ and that the
$\chi_{\alpha}$ are linearly independent, we get the matrix equations

$\left[t\tau_{3}\sigma_{2}s_{2}-itq_{x}\tau_{3}\sigma_{2}s_{3}+t\tau_{3}\sigma_{1}s_{0}\right]\chi_{\alpha}=0,$
(24)

where we omit the external product for convenience. Multiplying the whole
equation by $\tau_{3}\otimes\sigma_{2}\otimes s_{2}$ and dividing by $t$, we
obtain

$\left(\tau_{0}\sigma_{0}s_{0}+q_{x}\tau_{0}\sigma_{0}s_{1}-i\tau_{0}\sigma_{3}s_{2}\right)\chi_{\alpha}=0,$
(25)

which is diagonal in $\tau$ and $\sigma$, but not in $s$. The above equation
has non-trivial solutions when the determinant of the term between parenthesis
is zero:

$\left|\begin{matrix}1&q_{x}-m\\\
q_{x}+m&1\end{matrix}\right|=2-q_{x}^{2}=0\rightarrow q_{x}=\pm\sqrt{2}.$ (26)

where $m=\pm 1$ is the eigenstate of $\sigma_{3}$. Note that $q_{x}$ does not
depend on the eigenvalue $l$ of $\tau_{3}$. A solution that is localized on
the I edge should have negative $q_{x}$, while one localized on the III edge
should have a positive $q_{x}$. Therefore, for the I edge, $q_{x}=-\sqrt{2}$,
while for the III edge, $q_{x}=\sqrt{2}$.

Focusing first on the I edge, we write the localized solutions as

$\chi_{l,m}=\ket{l,m}\otimes\left(a\ket{s_{3}=1}+b\ket{s_{3}=-1}\right),$

such that $a$ and $b$ are given by

$\displaystyle\begin{pmatrix}1&-\sqrt{2}-m\\\
-\sqrt{2}+m&1\end{pmatrix}\begin{pmatrix}a\\\ b\end{pmatrix}=0$
$\displaystyle\therefore b=\left(m-\sqrt{2}\right)a,\quad
a=\sqrt{\frac{1}{4-2m\sqrt{2}}},$ (27)

where in the last equality, we used that $a^{2}+b^{2}=1$. Next we add the
other terms in the Hamiltonian of Eq. (11), to obtain the energy of the modes
localized on the I edge by computing $\chi_{l,m}^{\dagger}h_{kyee}\chi_{l,m}$.
For that, we notice that

$\displaystyle\left[\bra{s_{3}=1}+\left(m-\sqrt{2}\right)\bra{s_{3}=-1}\right]$
$\displaystyle s_{1}$
$\displaystyle\left[\ket{s_{3}=1}+\left(m^{\prime}-\sqrt{2}\right)\ket{s_{3}=-1}\right]=m+m^{\prime}-2\sqrt{2}$
(28)
$\displaystyle\left[\bra{s_{3}=1}+\left(m-\sqrt{2}\right)\bra{s_{3}=-1}\right]$
$\displaystyle s_{2}$
$\displaystyle\left[\ket{s_{3}=1}+\left(m^{\prime}-\sqrt{2}\right)\ket{s_{3}=-1}\right]=i(m-m^{\prime}),$
(29)

such that we get

$\displaystyle
tk_{y}a\chi_{lm}^{\dagger}\Gamma_{1}\chi_{l^{\prime}m^{\prime}}$
$\displaystyle=$ $\displaystyle-
tk_{y}a\left(\tau_{3}\right)_{l,l^{\prime}}\left(\sigma_{2}\right)_{m,m^{\prime}}$
(31)
$\displaystyle\Delta\chi_{lm}^{\dagger}\Gamma_{5}\chi_{l^{\prime}m^{\prime}}$
$\displaystyle=$ $\displaystyle 0$ (32) $\displaystyle
B_{x}\chi_{lm}^{\dagger}\Gamma_{6}\chi_{l^{\prime}m^{\prime}}$
$\displaystyle=$
$\displaystyle-\frac{B_{x}}{\sqrt{2}}\left(\tau_{3}\right)_{l,l^{\prime}}\left(\sigma_{0}\right)_{m,m^{\prime}},$
(33)

where we use that the $\Gamma_{6}$ term vanishes because it is proportional to
$(m-m^{\prime})\left(\sigma_{0}\right)_{m,m^{\prime}}=0$. Combining all terms,
we arrive at the edge Hamiltonian

$h^{\textrm{I}}(k_{y})=-tk_{y}a\tau_{3}\otimes\sigma_{2}-\frac{B_{x}}{\sqrt{2}}\tau_{3}\otimes\sigma_{0}.$
(34)

For the III edge, we instead need to use the solution $q_{x}=\sqrt{2}$.
Performing an equivalent procedure to above, we obtain

$\chi_{lm}=\ket{\tau_{3}=l}\otimes\ket{\sigma_{3}=m}\otimes\frac{\ket{s_{3}=1}-\left(m+\sqrt{2}\right)\ket{s_{3}=-1}}{\sqrt{4+2m\sqrt{2}}},$

and

$h^{\textrm{III}}(k_{y})=-tk_{y}a\tau_{3}\otimes\sigma_{2}-\frac{B_{x}}{\sqrt{2}}\tau_{3}\otimes\sigma_{0}.$
(35)

For the II and IV edges, the solutions should instead be localized on $y$.
Therefore, we take $h_{0y}=t\Gamma_{2}+tk_{y}a\Gamma_{1}+t\Gamma_{4}$ to look
for localized solutions and consider the perturbation to be
$h_{ky}=tk_{x}a\Gamma_{3}+\Delta\Gamma_{5}+B_{x}\Gamma_{6}$, making
$k_{y}\rightarrow-i\partial_{y}$, and look for solutions of

$h_{0y}\psi(x,y)=\left[-it\partial_{y}a\Gamma_{1}+t\Gamma_{2}+t\Gamma_{4}\right]\psi(x,y)=0.$
(36)

Again, we use an ansatz $\varphi_{\alpha}(y)=\exp(q_{y}y/a)$ to obtain the
matrix equation

$\left[q_{y}\tau_{0}\sigma_{0}s_{0}-\tau_{0}\sigma_{0}s_{3}+\tau_{0}\sigma_{3}s_{1}\right]\chi_{\alpha}=0,$
(37)

after multiplying $h_{0y}$ by $i/t\,\tau_{3}\sigma_{2}s_{1}$. In this way,
$q_{y}$ is determined by

$\left|\begin{matrix}q_{y}-1&m\\\
m&q_{y}+1\end{matrix}\right|=q_{y}^{2}-2=0\rightarrow q_{y}=\pm\sqrt{2}.$ (38)

The solution localized on the II edge is the one with $q_{y}=-\sqrt{2}$, and
the one on the IV edge has $q_{y}=\sqrt{2}$. For the II edge, the solution is
given by

$\chi_{l,m}=\ket{l,m}\otimes\frac{\ket{s_{3}=1}-m\left(1+\sqrt{2}\right)\ket{s_{3}=-1}}{\sqrt{4+2\sqrt{2}}},$

and

$h^{\textrm{II}}(k_{x})=tk_{x}a\tau_{3}\otimes\sigma_{2}-\frac{1}{\sqrt{2}}B_{x}\,\tau_{3}\otimes\sigma_{3},$
(39)

while for the IV edge, the solution is given by

$\chi_{l,m}=\ket{l,m}\otimes\frac{\ket{s_{3}=1}+m\left(1-\sqrt{2}\right)\ket{s_{3}=-1}}{\sqrt{4+2\sqrt{2}}},$

and

$h^{\textrm{IV}}(k_{x})=tk_{x}a\tau_{3}\otimes\sigma_{2}-\frac{1}{\sqrt{2}}B_{x}\,\tau_{3}\otimes\sigma_{3}.$
(40)

## References

  * Altland and Zirnbauer [1997] A. Altland and M. R. Zirnbauer, Nonstandard symmetry classes in mesoscopic normal-superconducting hybrid structures, Phys. Rev. B 55, 1142 (1997).
  * Schnyder  _et al._ [2008] A. P. Schnyder, S. Ryu, A. Furusaki, and A. W. W. Ludwig, Classification of topological insulators and superconductors in three spatial dimensions, Phys. Rev. B 78, 195125 (2008).
  * Fu [2011] L. Fu, Topological crystalline insulators, Phys. Rev. Lett. 106, 106802 (2011).
  * Slager  _et al._ [2013] R.-J. Slager, A. Mesaros, V. Juričić, and J. Zaanen, The space group classification of topological band-insulators, Nat. Phys. 9, 98 (2013).
  * Kruthoff  _et al._ [2017] J. Kruthoff, J. de Boer, J. van Wezel, C. L. Kane, and R.-J. Slager, Topological classification of crystalline insulators through band structure combinatorics, Phys. Rev. X 7, 041069 (2017).
  * Po  _et al._ [2017] H. C. Po, A. Vishwanath, and H. Watanabe, Symmetry-based indicators of band topology in the 230 space groups, Nat. Commun. 8, 50 (2017).
  * Neupert and Schindler [2018] T. Neupert and F. Schindler, Topological crystalline insulators, Topological Matter: Lectures from the Topological Matter School 2017 , 31 (2018).
  * Benalcazar  _et al._ [2017a] W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Electric multipole moments, topological multipole moment pumping, and chiral hinge states in crystalline insulators, Phys. Rev. B 96, 245115 (2017a).
  * Benalcazar  _et al._ [2017b] W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Quantized electric multipole insulators, Science 357, 61 (2017b).
  * Schindler  _et al._ [2018a] F. Schindler, A. M. Cook, M. G. Vergniory, Z. Wang, S. S. Parkin, B. A. Bernevig, and T. Neupert, Higher-order topological insulators, Science advances 4, eaat0346 (2018a).
  * Schindler  _et al._ [2018b] F. Schindler, Z. Wang, M. G. Vergniory, A. M. Cook, A. Murani, S. Sengupta, A. Y. Kasumov, R. Deblock, S. Jeon, I. K. Drozdov, _et al._ , Higher-order topology in bismuth, Nat. Phys. 14, 918 (2018b).
  * Ezawa [2018] M. Ezawa, Higher-order topological electric circuits and topological corner resonance on the breathing kagome and pyrochlore lattices, Phys. Rev. B 98, 201402 (2018).
  * van Miert and Ortix [2018] G. van Miert and C. Ortix, Higher-order topological insulators protected by inversion and rotoinversion symmetries, Phys. Rev. B 98, 081110 (2018).
  * Benalcazar  _et al._ [2019] W. A. Benalcazar, T. Li, and T. L. Hughes, Quantization of fractional corner charge in ${C}_{n}$-symmetric higher-order topological crystalline insulators, Phys. Rev. B 99, 245151 (2019).
  * Călugăru  _et al._ [2019a] D. Călugăru, V. Juričić, and B. Roy, Higher-order topological phases: A general principle of construction, Phys. Rev. B 99, 041301 (2019a).
  * Nag  _et al._ [2019] T. Nag, V. Juričić, and B. Roy, Out of equilibrium higher-order topological insulator: Floquet engineering and quench dynamics, Phys. Rev. Res. 1, 032045 (2019).
  * Benalcazar and Cerjan [2020] W. A. Benalcazar and A. Cerjan, Bound states in the continuum of higher-order topological insulators, Phys. Rev. B 101, 161116 (2020).
  * Agarwala  _et al._ [2020] A. Agarwala, V. Juričić, and B. Roy, Higher-order topological insulators in amorphous solids, Phys. Rev. Res. 2, 012067 (2020).
  * Arouca  _et al._ [2020] R. Arouca, S. N. Kempkes, and C. Morais Smith, Thermodynamics of a higher-order topological insulator, Phys. Rev. Res. 2, 023097 (2020).
  * Saha  _et al._ [2023] S. Saha, T. Nag, and S. Mandal, Multiple higher-order topological phases with even and odd pairs of zero-energy corner modes in a c3symmetry broken model, Europhysics Lett. 142, 56002 (2023).
  * Nag  _et al._ [2021] T. Nag, V. Juričić, and B. Roy, Hierarchy of higher-order floquet topological phases in three dimensions, Phys. Rev. B 103, 115308 (2021).
  * Ghosh  _et al._ [2022a] A. K. Ghosh, T. Nag, and A. Saha, Systematic generation of the cascade of anomalous dynamical first- and higher-order modes in floquet topological insulators, Phys. Rev. B 105, 115418 (2022a).
  * Călugăru  _et al._ [2019b] D. Călugăru, V. Juričić, and B. Roy, Higher-order topological phases: A general principle of construction, Phys. Rev. B 99, 041301 (2019b).
  * Noguchi  _et al._ [2021] R. Noguchi, M. Kobayashi, Z. Jiang, K. Kuroda, T. Takahashi, Z. Xu, D. Lee, M. Hirayama, M. Ochi, T. Shirasawa, _et al._ , Evidence for a higher-order topological insulator in a three-dimensional material built from van der waals stacking of bismuth-halide chains, Nat. Mater. 20, 473 (2021).
  * Imhof  _et al._ [2018] S. Imhof, C. Berger, F. Bayer, J. Brehm, L. W. Molenkamp, T. Kiessling, F. Schindler, C. H. Lee, M. Greiter, T. Neupert, _et al._ , Topolectrical-circuit realization of topological corner modes, Nat. Phys. 14, 925 (2018).
  * Serra-Garcia  _et al._ [2019] M. Serra-Garcia, R. Süsstrunk, and S. D. Huber, Observation of quadrupole transitions and edge mode topology in an lc circuit network, Phys. Rev. B 99, 020304 (2019).
  * Ni  _et al._ [2019] X. Ni, M. Weiner, A. Alù, and A. B. Khanikaev, Observation of higher-order topological acoustic states protected by generalized chiral symmetry, Nat Mater. 18, 113 (2019).
  * Xue  _et al._ [2019] H. Xue, Y. Yang, F. Gao, Y. Chong, and B. Zhang, Acoustic higher-order topological insulator on a kagome lattice, Nat. Mater. 18, 108 (2019).
  * El Hassan  _et al._ [2019] A. El Hassan, F. K. Kunst, A. Moritz, G. Andler, E. J. Bergholtz, and M. Bourennane, Corner states of light in photonic waveguides, Nat.Photonics 13, 697 (2019).
  * Kempkes  _et al._ [2019] S. Kempkes, M. Slot, J. van den Broeke, P. Capiod, W. Benalcazar, D. Vanmaekelbergh, D. Bercioux, I. Swart, and C. M. Smith, Robust zero-energy modes in an electronic higher-order topological insulator, Nat. Mater. 18, 1292 (2019).
  * Cerjan  _et al._ [2020] A. Cerjan, M. Jürgensen, W. A. Benalcazar, S. Mukherjee, and M. C. Rechtsman, Observation of a higher-order topological bound state in the continuum, Phys. Rev. Lett. 125, 213901 (2020).
  * Zheng  _et al._ [2022] S. Zheng, X. Man, Z.-L. Kong, Z.-K. Lin, G. Duan, N. Chen, D. Yu, J.-H. Jiang, and B. Xia, Observation of fractal higher-order topological states in acoustic metamaterials, Science Bulletin 67, 2069 (2022).
  * Li  _et al._ [2022] J. Li, Q. Mo, J.-H. Jiang, and Z. Yang, Higher-order topological phase in an acoustic fractal lattice, arXiv preprint arXiv:2205.05298 (2022).
  * Wang  _et al._ [2018] Y. Wang, M. Lin, and T. L. Hughes, Weak-pairing higher order topological superconductors, Phys. Rev. B 98, 165144 (2018).
  * Ghosh  _et al._ [2021a] A. K. Ghosh, T. Nag, and A. Saha, Floquet generation of a second-order topological superconductor, Phys. Rev. B 103, 045424 (2021a).
  * Vasiloiu  _et al._ [2022] L. M. Vasiloiu, A. Tiwari, and J. H. Bardarson, Dephasing-enhanced majorana zero modes in two-dimensional and three-dimensional higher-order topological superconductors, Phys. Rev. B 106, L060307 (2022).
  * Wang  _et al._ [2023a] Y. Wang, G. Rai, S. Haas, and A. Jagannathan, Edge and corner superconductivity in a two-dimensional topological model, Phys. Rev. B 107, 104507 (2023a).
  * Ghosh  _et al._ [2021b] A. K. Ghosh, T. Nag, and A. Saha, Hierarchy of higher-order topological superconductors in three dimensions, Phys. Rev. B 104, 134508 (2021b).
  * Ghosh  _et al._ [2021c] A. K. Ghosh, T. Nag, and A. Saha, Floquet generation of a second-order topological superconductor, Phys. Rev. B 103, 045424 (2021c).
  * Ghosh and Nag [2022] A. K. Ghosh and T. Nag, Non-hermitian higher-order topological superconductors in two dimensions: Statics and dynamics, Phys. Rev. B 106, L140303 (2022).
  * Ghosh  _et al._ [2022b] A. K. Ghosh, T. Nag, and A. Saha, Dynamical construction of quadrupolar and octupolar topological superconductors, Phys. Rev. B 105, 155406 (2022b).
  * Roy [2020] B. Roy, Higher-order topological superconductors in $\mathcal{P}$-, $\mathcal{T}$-odd quadrupolar dirac materials, Phys. Rev. B 101, 220506 (2020).
  * Kitaev [2001] A. Y. Kitaev, Unpaired majorana fermions in quantum wires, Physics-Uspekhi 44, 131 (2001).
  * Lutchyn  _et al._ [2010] R. M. Lutchyn, J. D. Sau, and S. Das Sarma, Majorana fermions and a topological phase transition in semiconductor-superconductor heterostructures, Phys. Rev. Lett. 105, 077001 (2010).
  * Oreg  _et al._ [2010] Y. Oreg, G. Refael, and F. von Oppen, Helical liquids and majorana bound states in quantum wires, Phys. Rev. Lett. 105, 177002 (2010).
  * Mourik  _et al._ [2012] V. Mourik, K. Zuo, S. M. Frolov, S. Plissard, E. P. Bakkers, and L. P. Kouwenhoven, Signatures of majorana fermions in hybrid superconductor-semiconductor nanowire devices, Science 336, 1003 (2012).
  * Nadj-Perge  _et al._ [2013] S. Nadj-Perge, I. K. Drozdov, B. A. Bernevig, and A. Yazdani, Proposal for realizing majorana fermions in chains of magnetic atoms on a superconductor, Phys. Rev. B 88, 020407 (2013).
  * Nadj-Perge  _et al._ [2014] S. Nadj-Perge, I. K. Drozdov, J. Li, H. Chen, S. Jeon, J. Seo, A. H. MacDonald, B. A. Bernevig, and A. Yazdani, Observation of majorana fermions in ferromagnetic atomic chains on a superconductor, Science 346, 602 (2014).
  * Baldo  _et al._ [2023] L. Baldo, L. G. G. V. D. D. Silva, A. M. Black-Schaffer, and J. Cayao, Zero-frequency supercurrent susceptibility signatures of trivial and topological zero-energy states in nanowire junctions, Superconductor Science and Technology 36, 034003 (2023).
  * Awoga  _et al._ [2023] O. A. Awoga, M. Leijnse, A. M. Black-Schaffer, and J. Cayao, Mitigating disorder-induced zero-energy states in weakly coupled superconductor-semiconductor hybrid systems, Phys. Rev. B 107, 184519 (2023).
  * Nayak  _et al._ [2008] C. Nayak, S. H. Simon, A. Stern, M. Freedman, and S. Das Sarma, Non-abelian anyons and topological quantum computation, Rev. Mod. Phys. 80, 1083 (2008).
  * Löfwander  _et al._ [2001] T. Löfwander, V. S. Shumeiko, and G. Wendin, Andreev bound states in high-tc superconducting junctions, Superconductor Science and Technology 14, R53 (2001).
  * Ryu and Hatsugai [2002] S. Ryu and Y. Hatsugai, Topological origin of zero-energy edge states in particle-hole symmetric systems, Phys. Rev. Lett. 89, 077002 (2002).
  * Schnyder  _et al._ [2012] A. P. Schnyder, P. M. R. Brydon, and C. Timm, Types of topological surface states in nodal noncentrosymmetric superconductors, Phys. Rev. B 85, 024522 (2012).
  * Matsuura  _et al._ [2013] S. Matsuura, P.-Y. Chang, A. P. Schnyder, and S. Ryu, Protected boundary states in gapless topological phases, N. J. Phys. 15, 065001 (2013).
  * Potter and Lee [2014] A. C. Potter and P. A. Lee, Edge ferromagnetism from majorana flat bands: Application to split tunneling-conductance peaks in high-${T}_{c}$ cuprate superconductors, Phys. Rev. Lett. 112, 117002 (2014).
  * Schnyder and Brydon [2015] A. P. Schnyder and P. M. Brydon, Topological surface states in nodal superconductors, J. Phys.: Condens. Matter 27, 243201 (2015).
  * Chakraborty  _et al._ [2022] D. Chakraborty, T. Löfwander, M. Fogelström, and A. M. Black-Schaffer, Disorder-robust phase crystal in high-temperature superconductors stabilized by strong correlations, npj Quantum Materials 7, 44 (2022).
  * Bonetti  _et al._ [2023] P. M. Bonetti, D. Chakraborty, X. Wu, and A. P. Schnyder, van hove, rashba, and hubbard meet to form first-order and higher-order topological superconductors, arXiv preprint arXiv:2304.07100 (2023).
  * Li and Haldane [2008] H. Li and F. D. M. Haldane, Entanglement spectrum as a generalization of entanglement entropy: Identification of topological order in non-abelian fractional quantum hall effect states, Phys. Rev. Lett. 101, 010504 (2008).
  * Hughes  _et al._ [2011a] T. L. Hughes, E. Prodan, and B. A. Bernevig, Inversion-symmetric topological insulators, Phys. Rev. B 83, 245132 (2011a).
  * Fukui and Hatsugai [2018] T. Fukui and Y. Hatsugai, Entanglement polarization for the topological quadrupole phase, Phys. Rev. B 98, 035147 (2018).
  * You  _et al._ [2020a] Y. You, J. Bibo, and F. Pollmann, Higher-order entanglement and many-body invariants for higher-order topological phases, Phys. Rev. Res. 2, 033192 (2020a).
  * Zhu  _et al._ [2020a] P. Zhu, K. Loehr, and T. L. Hughes, Identifying ${C}_{n}$-symmetric higher-order topology and fractional corner charge using entanglement spectra, Phys. Rev. B 101, 115140 (2020a).
  * Sticlet  _et al._ [2012] D. Sticlet, C. Bena, and P. Simon, Spin and majorana polarization in topological superconducting wires, Phys. Rev. Lett. 108, 096802 (2012).
  * Bena [2017] C. Bena, Testing the formation of majorana states using majorana polarization, Comptes Rendus Physique 18, 349 (2017).
  * Sato  _et al._ [2009] M. Sato, Y. Takahashi, and S. Fujimoto, Non-abelian topological order in $s$-wave superfluids of ultracold fermionic atoms, Phys. Rev. Lett. 103, 020401 (2009).
  * Su  _et al._ [1979] W. P. Su, J. R. Schrieffer, and A. J. Heeger, Solitons in polyacetylene, Phys. Rev. Lett. 42, 1698 (1979).
  * Zhu  _et al._ [2020b] P. Zhu, K. Loehr, and T. L. Hughes, Identifying ${C}_{n}$-symmetric higher-order topology and fractional corner charge using entanglement spectra, Phys. Rev. B 101, 115140 (2020b).
  * King-Smith and Vanderbilt [1993] R. D. King-Smith and D. Vanderbilt, Theory of polarization of crystalline solids, Phys. Rev. B 47, 1651 (1993).
  * Resta [1998] R. Resta, Quantum-mechanical position operator in extended systems, Phys. Rev. Lett. 80, 1800 (1998).
  * Spaldin [2012] N. A. Spaldin, A beginner’s guide to the modern theory of polarization, Journal of Solid State Chemistry 195, 2 (2012).
  * Wilson [1974] K. G. Wilson, Confinement of quarks, Phys. Rev. D 10, 2445 (1974).
  * Herviou  _et al._ [2019] L. Herviou, N. Regnault, and J. H. Bardarson, Entanglement spectrum and symmetries in non-hermitian fermionic non-interacting models, SciPost Physical 7, 069 (2019).
  * Ortega-Taberner  _et al._ [2022] C. Ortega-Taberner, L. Rødland, and M. Hermanns, Polarization and entanglement spectrum in non-hermitian systems, Phys. Rev. B 105, 075103 (2022).
  * Hughes  _et al._ [2011b] T. L. Hughes, E. Prodan, and B. A. Bernevig, Inversion-symmetric topological insulators, Phys. Rev. B 83, 245132 (2011b).
  * You  _et al._ [2020b] Y. You, J. Bibo, and F. Pollmann, Higher-order entanglement and many-body invariants for higher-order topological phases, Phys. Rev. Res. 2, 033192 (2020b).
  * Maiellaro  _et al._ [2023] A. Maiellaro, J. Settino, C. Guarcello, F. Romeo, and R. Citro, Hallmarks of orbital-flavored majorana states in josephson junctions based on oxide nanochannels, Phys. Rev. B 107, L201405 (2023).
  * Vojta [2003] M. Vojta, Quantum phase transitions, Rep. Prog. Phys. 66, 2069 (2003).
  * Ma  _et al._ [2022] H. Ma, Z. Zhang, P.-H. Fu, J. Wu, and X.-L. Yu, Electronic and topological properties of extended two-dimensional su-schrieffer-heeger models and realization of flat edge bands, Phys. Rev. B 106, 245109 (2022).
  * Wang  _et al._ [2021] B. Wang, X. Zhou, H. Lin, and A. Bansil, Higher-order topological insulator phase in a modified haldane model, Phys. Rev. B 104, L121108 (2021).
  * Ghorashi  _et al._ [2020] S. A. A. Ghorashi, T. Li, and T. L. Hughes, Higher-order weyl semimetals, Phys. Rev. Lett. 125, 266804 (2020).
  * Zhang  _et al._ [2020] X. Zhang, Z.-K. Lin, H.-X. Wang, Z. Xiong, Y. Tian, M.-H. Lu, Y.-F. Chen, and J.-H. Jiang, Symmetry-protected hierarchy of anomalous multipole topological band gaps in nonsymmorphic metacrystals, Nature Communications 11, 65 (2020).
  * Yang  _et al._ [2021] Y. Yang, J. Lu, M. Yan, X. Huang, W. Deng, and Z. Liu, Hybrid-order topological insulators in a phononic crystal, Phys. Rev. Lett. 126, 156801 (2021).
  * Wang  _et al._ [2023b] Z. Wang, X. Wang, Z. Hu, D. Bongiovanni, D. JukiÄ‡, L. Tang, D. Song, R. Morandotti, Z. Chen, and H. Buljan, Sub-symmetry-protected topological states, Nature Physics 19, 992 (2023b).
  * Hossain  _et al._ [2024] M. S. Hossain, F. Schindler, R. Islam, Z. Muhammad, Y.-X. Jiang, Z.-J. Cheng, Q. Zhang, T. Hou, H. Chen, M. Litskevich, _et al._ , Discovery of a hybrid topological quantum state in an elemental solid, arXiv preprint arXiv:2401.04845 (2024).
  * Jackiw and Rebbi [1976] R. Jackiw and C. Rebbi, Solitons with fermion number ½, Phys. Rev. D 13, 3398 (1976).
  * Lapp  _et al._ [2023] C. J. Lapp, J. M. Link, and C. Timm, Complete zero-energy flat bands of surface states in fully gapped chiral noncentrosymmetric superconductors, arXiv preprint arXiv:2310.14800 (2023).