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

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 with bulk energy gaps, 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].
Topologically protected boundary states can also appear in nodal superconductors [53][54][55][56][57][58][59][60][61][62][63].Topological nodal superconductivity is a non-trivial phase that is however not contained in the Altland-Zirnbauer classification since the bulk of the system is gapless at the nodal points [55,56,58].Nevertheless, the presence of the nodal points is protected by symmetry, and there also 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 and various nodal states.This would both provide realizations of different individual topological phases and possible intriguing combinations thereof, possibly even uncovering previously unknown phases, and also offer tunability between the different phases and their characteristic properties.To be precise, being motivated by both recent studies on the BBH model in the normal, i.e. non-superconducting, state and already existing HOTSC [39,[64][65][66], we in this work seek the uncover the different superconducting phases in the BBH model in the presence of a tunable magnetic field.We are concerned with mapping out the full phase diagram, and, in particular, we focus on unexpected topological superconducting phases generated by the intricate interplay between the higher-order topology of the normal state and superconductivity and magnetic field.
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 swave 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 proximityinduced superconducting order parameter ∆ 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 [67][68][69][70][71] for a quarter of the lattice indicated as in Fig. 1(a), and the Majorana polarization [72,73], 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 ∆ ≈ 0 the result is a superconducting version of the second-order topological phase of the BBH model, characterized by eight corner states and displayed as the HOTSC phase (yellow) in Fig. 1(b).The pairing makes the corner states of the BBH model transform into Andreev bound states, 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 and unexpected topological phases appear.For larger values of both ∆ 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 on the edges 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 of edge localized states do not grow with system size, as usually expected.Further increasing B x and ∆, we find two other also interesting topological phases.These are both first-order topological phases with flat bands at zero energy as their boundary modes.One of this phases is a traditional nodal phase: a nodal flat band phase (blue), where 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, just as expected for a topological nodal superconductor.In contrast, in the other topological flat band phase the bulk stays completely gapped, but we nevertheless still find flat bands at zero energy, now spanning across the whole Brillouin zone and localized at the edges along the x-direction.We name this latter, unexpected, phase 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.
Our results show how a seemingly simple model, a normal state with spin-orbit coupling together with proximity-induced conventional superconductivity and in-plane magnetic field, can generate a plethora of different topological phases.The different kind of surface states can further be accessible by tuning physical parameters, in particular the magnetic field that can be externally controlled.Importantly, our results 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, in addition to already known topological phases.Our findings thus increase the catalogue of emergent topological superconducting phases and establish their microscopic origin.
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 indicators 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 ∆ 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 [74] and surface [39,60] 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 [8,9], is described by the Bloch Hamiltonian where The model is defined on a square lattice with a lattice parameter a.In addition to chiral (Π = σ 3 ⊗ s 0 ) symmetry1 , this model is also reflection symmetric along x (M x = iσ 1 ⊗s 3 ) and y (M y = iσ 1 ⊗s 1 ), besides presenting symmetries [8].For λ < t, this model presents topological corner modes protected by both chiral and C 4 symmetries.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 [75] with a π-flux.Here, we instead interpret σ i and s i as Pauli matrices (σ 0 = s 0 = 1 2 ) associated with orbital and spin degree of freedom, respectively.This Hamiltonian then has spin-orbit coupling between orbital and spin on the same site (terms proportional to λ) and between different sites (terms proportional to t).This choice is motivated when considering the proximity effect from a conventional superconductor, where the pairing is s-wave (onsite) intra-orbital (proportional to σ 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 λ = 0.This both generates nontrivial topology in the normal state and removes the 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 [45,46,74], 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 resulting Bogoliubov-de Gennes (BdG) Hamiltonian for this system in the particle-hole basis becomes where Note that by this definition, ∆ = ∆σ 0 ⊗ s 2 which represents an intra-orbital spin-singlet pairing, as expected by proximity effect from an external conventional superconductor.This model has chiral Π = τ 1 ⊗ σ 0 ⊗ s 0 and particle-hole symmetry P = τ 1 ⊗ σ 0 ⊗ s 0 K.The discussion about crystalline symmetries is present in Sect.III A. We remark that in the deep topological limit of the BBH model, which is the parameter range for λ/t that we focus on in our work, only the corners are gapless in the normal state.Therefore, in a more realistic model, in which superconductivity is introduced self-consistently in the topological BBH system, pairing terms would technically only be present around the corners (or edges) of the system.However, since our main focus in this work is to understand the variety of topological phases that arise due to the interplay of the BBH model with superconductivity, we choose to add s-wave pairing without selfconsistent treatment.Our approach follows similar recent treatments in the literature of higher-order topological superconductors, see, for instance, Refs.[39,64,65].
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 already hinted by the phase diagram in Fig. 1(b), the model in Eq. ( 2) 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 first use the Wannier spectrum [7,8] 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 indicator, the entanglement spectrum [67? ], 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 [72,73], which indicates how much of a combination of electron and hole a state in a superconductor is.In addition to these indicators, we study the energy spectra using different boundary conditions and the local density of states (LDOS) at zero energy to further characterize and especially verify the various phases and the phase transitions inbetween them.Below, we briefly review these tools, all already used in different previous settings and models, in order to provide a comprehensive background to be able to explain the signatures of each topological phase.

A. Wannier spectrum
In this subsection, to make this work more selfcontained, we pedagogically review the use of the Wannier spectrum to characterize first-and higher-order topological phases as discussed, for instance, in Refs.[7,76,77].In the modern theory of polarization, the charge polarization is directly related to the Berry phase [7,[76][77][78] 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 [7,8,79].For a 2D system with periodic boundary conditions along x and open boundary conditions along y, the Wilson loop components are defined by [7] where P is a projector over the occupied eigenstates |u m (k x )⟩ of the (semiperiodic) Hamiltonian with momentum k x .W x is a unitary matrix which can be associated with the so-called Wannier Hamiltonian H The eigenvalues ν x of H x W 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 ν x by integer values, making these quantities, in general, defined mod 1.We choose a gauge where ν 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 ν x come in (ν, 1−ν) pairs [8].In this way, ν x = 0.5 are reflection invariant and indicate the presence of boundary modes protected by this symmetry, created by a nontrivial polarization along x [7,8].The same holds, mutatis mutandis, to periodic boundary conditions along y with boundaries open along x, obtaining a corresponding ν y .Thus, when either ν x or ν 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 ν x and ν y [8].
An illustrative example is the BBH model in Eq. ( 1).This model presents reflection symmetry along both x and y, which restricts ν x/y to appear in pairs.For |λ/t| < 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 ν x and ν y [8].For reference, for |λ/t| > 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 λ 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 higherorder 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 π 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 . Therefore, we cannot use the nested Wilson loop as a bulk invariant to diagnose our topological phases.In fact, we computed the nested Wilson loop for this model, which presents non-quantized values for any finite ∆, reinforcing that it is not a good invariant.We instead revert to the Wannier spectra along x and y.

B. Entanglement spectrum
Since we cannot use the nested Wilson loop spectrum, but only the Wannier spectrum along x or y, an alternative tool to characterize higher-order topological phases is useful 4 .One such tool has recently turned out to be the entanglement spectrum [67,[69][70][71][83][84][85].In the same way that a non-trivial polarization in the lattice can be determined using the Wannier spectrum, it can also be diagnozed by the entanglement spectrum [67,[69][70][71][83][84][85].To obtain the entanglement spectrum, we compute the correlation matrix in the occupied state |Ω⟩ where |Ω⟩ represents the (many-body) fermionic ground state and c † r,τ,σ,s creates a particle (τ = 1) or hole (τ = −1) in orbital σ with spin s at position r = (x, y).The entanglement spectrum ξ 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 the 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 [68], 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 [70? ].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).

C. Majorana polarization
To determine the Majorana nature of the boundary states the Majorana polarization P is useful, defined as [72,73,86] where ψ r,τ =1,σ,s;m (ψ r,τ =−1,σ,s;m ) is the particle (hole) component of the wavefunction of the m-th eigenstate at position r = (x, y) with orbital σ and spin s.This is a tool to characterize how much particle-hole symmetric a state is.In particular, the quantity x,y P m (x, y) x,y,τ,σ,s compares the value P with the usual probability density of a state m.Therefore, C quantifies how much of the wavefunction is particle-hole symmetric [73].
For systems that present only one isolated state per boundary, the Majorana polarization becomes an unambiguous indicator of a MZM.We note that, however, if there are many putative 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 still be unclear whether two such putative MZMs can actually recombine into a complex fermion.For such recombination to be able 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 completely 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 its sign across the lattice.

D. Spectral characterization
In addition to the topological invariants and indicators 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 normally 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 several 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 Γ = (0, 0), X = (π/a, 0), Y = (0, π/a), and M = (π/a, π/a).For phases with edge or corner states we also apply open boundary conditions in the appropriate directions.For all midgap states we also plot the sum of the Majorana polarization for states at zero energy and compute C in Eq. ( 8) to verify whether these states are MZMs or not.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 energy or frequency ω = 0.

IV. TOPOLOGICAL PHASES
In this section, we present our main analysis of the topological phases of the model in Eq. ( 2) using the topological invariants discussed in Sec.III.We focus on the deep topological limit of the BBH model, setting λ = 0 for simplicity, which allow us to obtain analytical expressions of the phase boundaries, but we remark that the topological phases are present for |λ/t| < 1.We refer to Appendix A for results on finite λ.For λ = 0, the topo-logical 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.

A. Phase boundaries
A quantum phase transition is accompanied by the closing of the gap of the system [87].Therefore, the topological phase boundaries are obtained by considering the energy gap, δ between the highest valence band and the lowest conduction band around zero energy.In Fig. 2 we plot δ as a function of the superconducting order parameter ∆ and magnetic field B x for different boundary conditions.For fully periodic boundary conditions, we can additionally obtain the critical lines analytically using the Bloch Hamiltonian, which we represent in dashed lines in Fig. 2(a), overlayed on the value of δ for each ∆ and B numerically computed using a real space Hamiltonian.For B x = 2t 2 + ∆(∆ − 2t) (magenta dashed line) the bulk spectrum closes at momenta (k x , k y ) = (0, π/2) and (π, π/2) and for ∆ = 2t 2 + B x (B x − 2t) (yellow dashed line) it closes at momentum (π/2, π/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 = 2t 2 + ∆(∆ + 2t) (green dashed line) at momentum (k x , k y ) = (π, π/2) and at ∆ = 2t 2 + B x (B x + 2t) (cyan dashed line) at momentum (π/2, π/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 additional pieces of information.First, we notice that there is a new gap-closing line at ∆ = √ 2t − 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 at this phase transition.Second, we notice that the nodal flat bands phase (blue in Fig 1(b)) seemingly hosts a small but finite gap, illustrated by the faint red arc-like features in Figs.2(a,c), but it does become completely gapless in Fig. 2(b).An analysis of the spectrum of this system for semi-periodic boundary conditions in Sect.IV D show that, in fact, for all values of B x and ∆ in this phase, the gap closes at k y = π/2 and different k x .Therefore, the spectrum is actually gapless in all three figures, and the red arcs are just finite-size effects 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 yopen boundary conditions, as seen in Fig. 2(b).This indicates zero-energy states localized to edges along the x-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 ∆ in all phases.

B. HOTSC phase
For small B x and ∆, 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 ϵ spectrum for fully periodic in Fig. 3(a), x-periodic in Fig. 3(b), yperiodic in Fig. 3(c), and fully open in Fig. 3(d) boundary conditions, we realize that the bulk 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 both x-and y-periodic boundary conditions in this phase.Plotting the LDOS at zero 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 ν along both x in Fig. 3(f) and along y in Fig. 3(g).The half-integer values of both ν x,y indicate the higher-order character of this phase5 .The higher-order aspect is corroborated by the entanglement spectrum ξ in Fig. 3(h), which shows distinct isolated half-quantized modes inside the midgap region.
Finding a total of eight corner modes, two at each corner, can be thought of as expected since the original BBH model at B x = ∆ = 0 host similar corner modes in this characteristic quadrupolar phase [39,88].To analyze these eight corner modes in our superconducting system, we consider the total (i.e. the sum) Majorana polarization for the states at zero energy in Fig. 4. We see that it is also localized in the corners of the system with a pattern that changes signs in the different corners, which is an important signature of MZMs [73].However, when computing C (defined in Eq. ( 8) and using the grey shaded area in Fig. 4), we obtain a much smaller than one (maximum of around 0.25, with an overall increase with increasing ∆).This indicates that, while these are zero-energy states, they should be classified as degenerate zero-energy Andreev bound states and not individual MZMs.In summary, the HOTSC phase can be viewed as a superconducting extension of the quadrupolar phase of the BBH model at least regarding its surface states, but where the corner states are now a combination of particles and holes as they are Bogoliubov quasiparticles, but not still not MZMs.

C. Hybrid phase
With increasing values of ∆ and B x , the system enters into what we name the hybrid phase, the 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, with a similar set of data as for the HOTSC phase earlier.We find that it has energy spectrum similar to the HOTSC phase: both the bulk, Fig. 5 Fig. 5(e).However, when inspecting the Wannier spectrum, Figs.5(f,g), we see that ν y is still half-quantized, while ν x exhibits a gap around 0.5.Further, the entanglement spectrum in Fig. 5(h) now shows a continuous array of midgap eigenvalues that are present symmetrically around ξ = 0.5.Therefore, this hybrid-order phase can be considered to be topologically distinct from the previous HOTSC phase.We remark that since ν x is also not fully quantized for larger values of ∆ and B in the HOTSC phase, the topological phase transition between the HOTSC and the hybrid phases is due to the change of localization of some of the zero energy states, as we show below.
In fact, 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, 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.
Although not very visible6 in the total zero-energy LDOS in Fig. 5 (e), we find that four, i.e. half, of the zeroenergy states in the hybrid phase now have a significant weight not just at the corners but also partially along the y-direction.To illustrate these localization properties better we plot the Majorana polarization divided up into two sets in Figs.6(a,b), respectively.Here it is now clear that four of the states are still corner states, just as in the HOTSC phase, but four other states are now substantially delocalized along the y-direction.We also find that both sign changing between different corners and C = 1, provides strong indications that these states are now MZMs.We thus conclude that with all eight zeroenergy modes now possible to spatially separate they can  In summary, we see that in this hybrid phase hosts a combination of edge and corner modes.Since there is no bulk gap closing (see Fig. 2 (a)) between the HOTSC and hybrid phase, we expect the change between the two phases to be caused by a change in the properties of the edge, which is similar to what happens in some other extrinsic higher-order phases [89].In fact, when considering the edge theory for this Hamiltonian in Sec.V B, we see that the magnetic field can change the mass profile at the edges, delocalizing some of the corner modes.The delocalized states are now instead localized at the edge and are still characterized by a quantized ν y .These states appear similar to the ones in the phases with polarization p p x = 0 and p y = 0.5 (or vice-versa) in the original BBH model when the hopping along x and y is different [8,18], named a dipolar phase in Ref. [18], which is an example of a boundary obstructed atomic insulator [90].Due to the intriguing combination of corner and edge boundary modes, we name this a hybrid-order phase, as it has an inherent mix of different topologies.However, we note that our use of the word hybrid should not be confused with the situation where more standard first-and secondorder phases are appearing at the same time, also recently called a hybrid phase [91][92][93][94][95]. Instead, our hybrid phase is a standard second-order phase appearing jointly with a dipolar phase.A clear difference is that in our hybrid phase the number of zero-energy boundary states remains constant (four corner modes and four edge modes), independent on system size, while any phase with a standard first-order character sees the number of edge modes grow with system size.We remark that the only invariant we use that can distinguish this phase from the HOTSC phase is the nested entanglement spectrum, which shows a continuous array of values in contrast to the HOTSC phase where there are just modes at 0.5.

D. Nodal flat bands phase
Increasing either ∆ or B x , we arrive at yet other phase, a 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 ̸ = 0, π.The bulk gap thus vanishes away from the high-symmetry line and that is why the gap closing is 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).
Further evidence that the zero-energy states are localized along the II and IV edges is found in the LDOS at zero energy 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, since the absence of corner modes clearly discards these phase as any higherorder phase.We find this fully consistent with the topological invariants: an absence of midgap modes in ν x in Fig. 7(f), while a number of 0.5 modes in ν y appears in Fig. 7(g).We further find that the number of zero-energy states and 0.5 modes in ν y scales with the number of unit cells along x, which further corroborates that this is an edge phenomenon along the x-direction.Here ν y is thus an invariant that characterizes the presence of these zeroenergy states, which is related by a bulk-boundary correspondence to the nodes in the bulk [2,56,58].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 as for e.g. the HOTSC phase.Instead, we find a discontinuous array of mid-gap eigenvalues symmetrically placed around 0.5.
To better understand the overall properties of the zeroenergy flat bands, 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, finding four states with zero energy.After diagonalization, we then have a wavefunction ψ(k x = 0, y).To obtain a complete real space wavefunction, we multiply which is equivalent to a partial Fourier transform.This is the wavefunction we use to compute the Majorana polarization P in Figs. 8. We find that P changes sign between the two edges for all modes, and importantly we find C = 1 when summing over the gray region, indicating that these edge states have Majorana properties.However, 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 still recombine into complex fermionic modes.Still, because these modes both present a strong Majorana polarization and we have only two possible particle-hole symmetric pairs that we can build with both spin and orbital degrees of freedom, we choose 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 two per edge and one per corner, as discussed in Sec.IV C.

E. Nodeless flat bands phase
Finally, moving closer to the diagonal of the phase diagram, with ∆/t ≈ B x /t > 1, we find a nodeless flat band phase, 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 midgap 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 ν x is completely gapped in Fig. 9(f), ν y in Fig. 9(g) has a finite number of 0.5 modes, which grows with the number of unit cells along x.This ν 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, gathering reasonably close 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.
We investigate the Majorana polarization of the states at zero energy in this phase in Fig. 10.We here again 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 here present without any corresponding nodal points in the bulk.Instead we have a flat zero-energy edge state entirely separated from the bulk gap spectrum.As such, they are a very different kind of surface state compared to those of the nodal flat band phase.We currently actually do not know how they can be fully classified in terms of symmetry-protected topological phases.Nevertheless, we still know that they are topologically protected by symmetry since the number of states at zero energy is determined by the number of 0.5 states in ν y , which is a topological invariant.Finally, we note that these flat bands states present C = 1, which makes us designate them 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 feasible, 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.

A. 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 D, 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 spinorbit 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 λ = 0): where we obtain a dual Hamiltonian h D , following the nomenclature of Ref. [74], with the form where we omit the outer product symbol ⊗ 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)τ 1 ⊗ σ 0 ⊗ s 3 , with extended s-wave intra-orbital spin-singlet symmetry given by −t cos(k x a)τ 1 ⊗ σ 0 ⊗ s 1 , and with s-wave oddinterorbital spin-triplet symmetry given by B x τ 1 ⊗σ 2 ⊗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 ̸ = 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 ∆ 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.

B. 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 [96], 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,39], 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 BdG in Eq. ( 2) and obtain h Γ (k) = tk y aΓ 1 +tΓ 2 +tk x aΓ 3 +tΓ 4 +∆Γ 5 +B x Γ 6 .(11) Here k = (k x , k y ) represents the continuum momentum in relation to the Γ point.
We can obtain the low-energy descriptions of any corner states of h Γ by substituting k x/y → −i∂ 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 Γ localized on the I and III edges are where N x is a normalization factor, |c l,m | 2 = 1, and the spinors χ l,m are given by and where l = ±1, m = ±1, and s = ±1 are eigenvalues of τ 3 , σ 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 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 τ 3 , while they are a combination of different spin states with the effective spin degrees of freedom represented by σ 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 where N y is a normalization factor, we have 1, and the spinors χ II/IV lm are given by and where again l = ±1, m = ±1, and s = ±1 are eigenvalues of τ 3 , σ 3 and s 3 , respectively.Again projecting on the solutions of Eqs. ( 18) and ( 19), we obtain the effective Dirac Hamiltonians 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 τ 3 ⊗ σ 0 , while on the II and IV edges it is proportional to τ 3 ⊗ σ 3 .Since eigenvalues of σ 0 and σ 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 τ 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 σ 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. CONCLUNDING REMARKS
In this work, we investigate the topological phases of an orbital BBH model proximitized to a conventional spinsinglet 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 resulting phase diagram by considering different boundary conditions and investigate the topology of each phase by calculating both the Wannier and entanglement spectra, as well as the Majorana polarization.
At weak superconducting pairing ∆ 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, at least when considering the surface states, although the latter now also have particle and hole character since the system is superconducting, but still not being MZMs.Beyond this expected HOTSC phase, we also find several other, much more unexpected, topological phases.First, an unusual hybrid phase (green in Fig. 1(b)) presents an atypical mix between a second-order and a dipolar topological phase.Here four zero-energy corner states from the HOTSC phase are preserved as expected for a second-order topological phase, 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 this hybrid phase, but the number of zero-energy states remains fixed, at eight states.Second, two additional first-order phases present symmetryprotected 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 surprisingly 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, importantly demonstrate the existence of an unusual hybrid mixing of topological phases and shows that a zero-energy edge flat band can exists also for a nodeless, i.e. fully gapped, bulk, such that the edge states are not continuously connected to the bulk bands.
In terms of experimental feasibility of the BBH superconductor model developed in this work, we can refer to the atom-optics setup for the realization of lattice tightbinding topological models [97].Apart from metamaterials [34,35,98], higher-order topological insulator phases have also been experimentally observed in van der Waals stacking of bismuth-halide [24], as well as Bi x Sb 1−x alloys [99].The s-wave superconductor can be, in principle, placed in proximity to the above materials with appropriate substrates, such that HOTSC phases are realized in the above material.It is worth noting that superconductivity can be induced in the surface states of Bi 2 Se 3 , HgTe etc. via the proximity effect [100][101][102], with a superconducting gap even around ∆ 0 ∼ 0.5 meV [102], despite the fully gapped bulk.Given the above experimental developments in the solid-state topology, we believe that the HOTSC phases obtained in our work are at least within future experimental reach.
We further note that we study the HOTSC model starting from a higher-order topological insulator limit of the underlying BBH model.It would also be interesting to study a BBH superconductor model starting from a metallic bulk, where by tuning the chemical potential different topological phases may be uncovered in addition to those found in this work.One can hence investigate the three-parameter interplay between chemical potential, superconducting gap function, and magnetic field for engineering of HOTSC phases.Moreover, the possibility to alternate between topological states with different spatial localizations by tuning the magnetic field may lead to interesting applications where these states are useful, including topological quantum computation.
Note added: During the final preparation of this work we became aware of the recent work Ref. [103], 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 three dimensions (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.
In the main text, we focus on the properties of Eq. ( 2) for λ = 0. Here, we discuss the effect of having a finite λ in the phases discussed in the main text.For large λ > t, the system is just a trivial superconductor.For smaller but finite λ, the phase boundaries of the system change considerably but the topological phases largely persist.In Fig. 11 we set λ = 0.5t and we repeat the plot of the energy gap δ as a function of B x and ∆ 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.

Dipolar phase
With λ finite, but smaller than t, we find a different kind of first-order topological phase appearing for higher values of either ∆ 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 also shown in the zero-energy LDOS in Fig. 12(e).The ν x invariant along x is not quantized, see Fig. 12(f), while the ν y presents two halfquantized values, see Fig. 12(g), indicating that there are four symmetry-protected modes.The entanglement spectrum shows a symmetric distribution of eigenvalues excluding 0.5 within the mid gap region, as depicted in Fig. 12(h).
In Fig. 13, we show the sum of the Majorana polarization for these states at zero energy.We see that they are localized on the edges of the system and also present C = 1, indicating that they are MZMs.These states are very similar to the edge localized, or dipolar, states that appear in the hybrid phase, see Fig. 6(a).Moreover, this phase resembles (a superconducting version of) the phase with p x = 0 and p y = 0.5 in the BBH model [8] called dipolar phase in Ref. [18], such that we also call it a dipolar phase.We note that the P of the zero-energy boundary states appropriately also resembles an electric dipole.For the I and III edges, we look for localized solutions on x.We split h Γ = h 0x +h ky , where h 0x = tΓ 2 +tk x aΓ 3 + tΓ 4 is the part which determine the zero states localized along x, while we treat h ky = tk y aΓ 1 + ∆Γ 5 + B x Γ 6 as its perturbation, reasonable given the fact that B x , ∆ ≪ t.Next we perform the substitution k x → −i∂ x and look for solutions of h 0x ψ(x, y) = [tΓ 2 − it∂ x aΓ 3 + tΓ 4 ] ψ(x, y) = 0. (C2) Using the ansatz ϕ α (x) = exp(q x x/a) and that the χ α are linearly independent, we get the matrix equations [tτ 3 σ 2 s 2 − itq x τ 3 σ 2 s 3 + tτ 3 σ 1 s 0 ] χ α = 0, (C3) where we omit the external product for convenience.
Multiplying the whole equation by τ 3 ⊗ σ 2 ⊗ s 2 and dividing by t, we obtain (τ 0 σ 0 s 0 + q x τ 0 σ 0 s 1 − iτ 0 σ 3 s 2 ) χ α = 0, (C4) which is diagonal in τ and σ, but not in s.The above equation has non-trivial solutions when the determinant of the term between parenthesis is zero: where m = ±1 is the eigenstate of σ 3 .Note that q x does not depend on the eigenvalue l of τ 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 = − √ 2, while for the III edge, q x = √ 2.

FIG. 1 .
FIG. 1.(a) Schematic realization of the superconducting BBH model: BBH system with conventional spin-singlet s-wave superconducting pairing ∆ induced by proximity effect from a substrate and an in-plane magnetic field Bx.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.

FIG. 2 .
FIG. 2. Energy gap δ between the highest valence band and lowest conduction band around zero energy, calculated from Eq. (2), as a function of superconducting pair amplitude ∆/t and in-plane magnetic field Bx/t for (a) fully-periodic, (b) xperiodic (open in y) and (c) y-periodic (open in x) boundary conditions.Parameters used: λ = 0, system size 32 unit cells in each direction.Dashed lines represent analytical results, see main text for a description.

FIG. 3 .
FIG. 3. Main features of the HOTSC phase.Energy spectrum ϵ 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 ϵ = 0. LDOS at zero energy (e).Wannier spectrum along x (f) and y (g).Entanglement spectrum ξ (h).Parameters used: Bx = 0.2t, ∆ = 0.1t, λ = 0, system size 32 unit cells for directions with open boundary conditions and 100 k−points.
FIG. 5. Main features of the hybrid phase.Energy spectrum ϵ for fully periodic (a), x-periodic (b), y-periodic (c), and open (d) boundary conditions.Inset in (d) is a zoom-in on the modes in the middle of the spectrum with energy ϵ = 0. LDOS at zero energy (e).Wannier spectrum along x (f) and y (g).Entanglement spectrum ξ (h).Parameters used: Bx = 0.9t, ∆ = 0.8t, λ = 0, system size 32 unit cells for directions with open boundary conditions and 100 k−points.

FIG. 7 .
FIG. 7. Main features of the nodal flat band phase.Energy spectrum ϵ 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 ϵ = 0. LDOS at zero energy (e).Wannier spectrum along x (f) and y (g).Entanglement spectrum ξ (h).Parameters used: Bx = 1.5t, ∆ = 0.9t, λ = 0, system size 32 unit cells for directions with open boundary conditions and 100 k−points.

FIG. 8 .
FIG. 8. Sum of the Majorana polarizations P for the four states with ϵ = 0 for kx = 0 in the nodal flat band phase, using a partial Fourier transform, as explained in the main text.The gray dashed area indicates the region where C is calculated.For all four states, C = 1, indicating that they are MZMs.Parameters same as Fig. 7.

FIG. 9 .
FIG. 9. Main features of the phase with nodeless flat bands.Energy spectrum ϵ 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 ϵ = 0. LDOS at zero energy (e).Wannier spectrum along x (f) and y (g).Entanglement spectrum ξ (h).Parameters used: Bx = 1.8t, ∆ = 2t, λ = 0, system size 32 unit cells for directions with open boundary conditions and 100 k−points.

FIG. 10 .
FIG. 10.Sum of the Majorana polarizations P for the four states with ϵ = 0 for kx = 0 in the nodeless flat band phase, using a partial Fourier transform, as explained in the main text.The gray dashed area indicates the region where C is calculated.For all four states, C = 1, indicating that they are MZMs.Parameters same as Fig. 9.

FIG. 11 .
FIG. 11.Energy gap δ derived from Eq. (2) as a function of ∆/t and Bx/t for (a) fully-periodic, (b) x-periodic (open in y) and (c) y−periodic (open in x) boundary conditions.Parameters used: λ = 0.5t, system size 32 unit cells in each direction.

FIG. 12 .
FIG. 12. Main features of the dipolar phase.Energy spectrum ϵ 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 ϵ = 0. LDOS at zero energy (e).Wannier spectrum along x (f) and y (g).Entanglement spectrum ξ (h).. Parameters used: Bx = 1.5t, ∆ = 0.2t, λ = 0.5, system size 32 unit cells for directions with open boundary conditions and 100 k−points.
s 1 is the total normalstate Hamiltonian.For the second equality we use the vectors D = (d 1 , d 2 , d 3 , d 4 , ∆, B x ) and Γ to write the Hamiltonian in terms of matrices in the particle-hole, orbital, and spin degrees of freedom.Using τ i as the Pauli matrices in the particle-hole degrees of freedom, we construct the matrices Γ