Discovery of higher-order topological insulators using the spin Hall conductivity as a topology signature

The discovery and realization of topological insulators, a phase of matter which hosts metallic boundary states when the d-dimension insulating bulk is confined to (d − 1)-dimensions, led to several potential applications. Recently, it was shown that protected topological states can manifest in (d − 2)-dimensions, such as hinge and corner states for three- and two-dimensional systems, respectively. These nontrivial materials are named higher-order topological insulators (HOTIs). Here we show a connection between spin Hall effect and HOTIs using a combination of ab initio calculations and tight-binding modeling. The model demonstrates how a non-zero bulk midgap spin Hall conductivity (SHC) emerges within the HOTI phase. Following this, we performed high-throughput density functional theory calculations to find unknown HOTIs, using the SHC as a criterion. We calculated the SHC of 693 insulators resulting in seven stable two-dimensional HOTIs. Our work guides novel experimental and theoretical advances towards higher-order topological insulator realization and applications.


INTRODUCTION
The discovery of symmetry protected topological phases of matter opened a new field in condensed matter physics. When a ddimensional topological material is confined to (d − 1)-dimension, symmetry protected gapless states, which are robust against perturbations, manifest at its boundaries 1-3 , named first-order topological insulator (TI). This fundamental characteristic of TIs is the so-called bulk-boundary correspondence principle 4 , which is mediated through a variety of symmetries, including time-reversal (TR) 1 , mirror planes 5 , non-symmorphic 6 , and particle-hole 7 . Recently, a novel category of nontrivial material that apparently violates the bulk-boundary correspondence was proposed, which is referred to as higher-order topological insulators (HOTIs) [8][9][10] . Specifically, HOTIs exhibit gapped states in the d− (bulk) and (d − 1)-dimensional boundaries, although exhibiting gapless protected states at the (d − 2) boundaries, e.g., hinge gapless states. Despite the abundance of first-order TIs in nature [11][12][13] , HOTIs are seemingly less common [14][15][16][17][18][19][20][21] , indicating a need to develop novel strategies or principles capable of discovering and design of such systems.
The metallic character of topological protected states historically established a connection between topological protection and quantized values of the conductivity, as first expressed by the TKNN number 22 relating the quantum Hall effect to the quantized Hall conductivity, i.e., σ xy = ne 2 /h. Analogously, Murakami et al. proposed a strategy to obtain finite spin Hall conductivity (SHC) in narrow-gap insulators such as HgS and PbSe 23 , which later were identified as first-order TIs 11 , suggesting for the first time a connection between the SHC and nontrivial topology. The SHC in first-order TIs is given by σ z xy ¼ C s e 2 =h, where C s is the so-called spin Chern number, defined as the first Chern number difference for each spin channel (C s = C ↑ − C ↓ ). Thus, while normal undoped insulators show zero midgap SHC 24 , TIs possess a perfect quantized SHC. Nevertheless, small deviations can occur due to the absence of conservation of the spin angular momentum z-component (S z ), owing to effects such as crystal field, hybridization, symmetry breaking terms 1 , or disorder 24 . Conversely, the relation between the SHC and HOTI materials has not been explored so far. The few known HOTIs is one of the bottlenecks for the understanding of this connection.
Here we demonstrate that the SHC is a signature of HOTI phases. For this purpose, we first employ a tight-binding model to elucidate the conditions leading to a HOTI phase and its connection with a large bulk-SHC inside the topological gap. With the connection between finite midgap SHC and HOTI phases established, we perform a systematic search for novel HOTIs. Except for time-reversal, no other symmetry is imposed in our method. We calculate the SHC of 693 insulators reported in the 2D materials database C2DB 25 , followed by a screening procedure consisting of the selection of only thermodynamically stable nonmagnetic insulating materials. Among the possible candidates, we discover seven stable two-dimensional HOTIs: BiSe, PbF, PbBr, PbCl, HgTe, and BiTe in both GaSe and CH prototype structures. Real-space projections of midgap states obtained from DFT calculations show symmetry-protected localization of topological states in nanoflakes of these systems. We compute the topological invariant for the representative candidate BiSe, validating our predictions.
For 2D HOTIs, the total SHC of the (d − 2)-finite system is expected to be zero since there cannot be any total currents in a system with open boundary conditions. However, in general, in the middle of any 0D flakes, the local SHC is equal to the bulk (periodic) value, e.g., for 0D TI flakes the local SHC is − e/(2π) 26 .
Here we focus on the bulk SHC, which corresponds to the local SHC value in the (d − 2) HOTI nanoflake, rather than the number of boundary topological protected states. The established numerical relation between higher-order topology and bulk SHC 1 Department of Physics, Fluminense Federal University, Niterói, Rio de Janeiro, Brazil. 2 Brazilian Nanotechnology National Laboratory (LNNano), CNPEM, Campinas, Brazil. 3 Center for Natural and Human Sciences, Federal University of ABC, Santo André, SP, Brazil. 4 Department of Physics and Department of Chemistry, University of North Texas, Denton, TX, USA. ✉ email: adalberto.fazzio@lnnano.cnpem.br advances the understanding of HOTIs and demonstrates an important design principle to predict these compounds.

RESULTS
SHC as a signature of HOTIs Initially, we demonstrate how a HOTI shows a non-zero midgap SHC. We start from an eight-band tight-binding Hamiltonian, proposed by Schindler et al. in ref. 14 . The model is adapted for a two-dimensional material and is given by where we consider two full spin orbitals μ and ν on each site. The site index is represented by i and c y i (c i ) is the creation (annihilation) fermionic operator. The Pauli matrices ρ, τ, and σ stand for the orbital μ, ν, and spin, respectively. The term d e = ±1 for e = x, y, while M, t, ξ, and Δ are the mass, hopping, spin-orbit coupling (SOC) strength, and mixing term, respectively. This model realizes a 2D HOTI that preserves TR-symmetry (T ¼ Àis y K) and four-fold rotational symmetry (C 4 ¼ τ z e Àiπsz=4 ) 27 . Setting the parameters M = 2t = 2ξ = 1 and Δ = 0.25, the bulk (2D) band structure with a 2 eV bandgap is obtained as seen in Fig. 1a. In a torus geometry, a gapped band structure is obtained (see Supplementary information Fig. 1), typical of HOTIs since its nontrivial topology is manifested in (d − 2). Figure 1b shows the 0D energy spectrum, i.e., square nanoflake geometry (open boundary conditions). The inset highlights the zero energy modes (red circles) which are localized at the flake corners (d − 2 boundary states), as shown in Fig. 1c. In Fig. 1a side plot, we show the calculated HOTI model SHC σ z xy % 0:75 e 2 =h, showing that SOC introduces a non-zero SHC into HOTI phases. To deepen our analysis, we calculate the SHC landscape in a given range of model parameters, 0 ≤ t ≤ 2, 0 ≤ Δ ≤ 1, and 0 < M ≤ 10, keeping the SOC strength fixed at ξ = 0.5. The resulting map in Fig. 1d shows a clear inverse relation between the SHC and Δ. To clarify their relationship we analyze cross sections of the SHC landscape in Fig. 1e-g. For Δ = 0 (Fig. 1g), we obtain a TI with a quantized SHC in the blue region. The phase boundary between the TI and NI region depends on the |M|/t ratio and is depicted as a black line. For Δ = 0.25 the HOTI phase is obtained, as can be seen in the red region, where the SHC is no longer quantized, which is evidenced by its pale color. Such reduction in SHC is more pronounced as Δ increases, according to Fig. 1g, where Δ = 0.75. As previously stated, the Δ term is responsible for driving a TI into the HOTI phase. This analysis of the model Hamiltonian evidences the gradual suppression of the SHC as the system goes further into the HOTI phase, i.e., when Δ increases.

First principles calculations
We performed a screening over the 2D materials database C2DB 25 . The C2DB is constructed based on prototypes of known 2D materials, i.e., Graphene, MoS 2 , BN, and others. With these prototypes, a combinatorial approach created around 4000 2D compounds for which the properties were obtained via density functional theory (DFT) calculations. The authors established criteria to categorize the materials as low, medium, and high stability. Thermodynamic stability is evaluated by computing each candidate's convex Hull. The dynamic stability is based on the phonon spectra of experimentally synthesized 2D materials; see ref. 25 for details. We performed a screening over the C2DB database to select materials with appropriate characteristics. Initially, we selected thermodynamic and dynamic stable materials. After this initial screening, we narrowed our search over nonmagnetic insulators, resulting in 693 2D materials.
Following the initial screening we performed fully relativistic density functional theory [28][29][30] calculations using the QUANTUM ESPRESSO package 31 . We calculated the spin Hall conductivity (SHC) of all selected 2D materials and validated our results with the C2DB topological classification reported in ref. 32 . Then, we further selected all materials with large SHC and not classified as topological insulators or topological crystalline insulators. Finally, for the predicted HOTIs we calculated the electronic structure of triangular nanoflakes, corresponding to open boundary conditions, using the VASP software code 33 . To avoid spurious interactions, 15 Å of vacuum is used in all calculations.
The SHC was calculated using local effective Hamiltonians constructed via the pseudo-atomic orbital (PAO) projection method 34,35 , implemented in the PAOFLOW code 36 . The method consists in projecting the plane wave Kohn-Sham states, composed of several thousand of basis functions, into a compact subspace spanned by PAO orbitals. The PAO Hamiltonian allows the efficient and accurate computation of the SHC 37 , topological properties 38,39 , dynamical properties 40 , and others. A comparison between the DFT and PAO electronic structures is presented in the Supplementary information Fig. 6. The SHC was calculated via linear response theory using the Kubo formula, where s is the quantization axis (spin polarization), j is the applied electrical field orientation, and i gives the spin current direction ( σ s ij ). The Fermi-Dirac distribution is represented by f ð k ! Þ and Ω s n;ij ð k ! Þ is the spin Berry curvature given by where m and n are band indices, ψ and E are the PAO Hamiltonian wavefunction and eigenvalues, respectively. The spin current operator (j s i ) is defined as the anticommutator of the Pauling matrix (σ s ) and velocity operators ({σ s , v i }). For the SHC the Brillouin zone sampling was increased to 72.0/Å −1 . The use of nondegenerate perturbation theory in the Kubo formula can lead to unphysical midgap SHC in trivial insulators 41 . To circumvent this issue PAOFLOW uses degenerate perturbation theory, see ref. 37 for details. In the case of spin rotation symmetry breaking, the spin quantization axis is not well defined, leading to a reduction of the midgap SHC in TIs 1,42 . Nevertheless, it does not induce unphysical midgap SHC in trivial insulators. In our calculations, the only nonzero midgap SHCs are due to TIs or HOTIs (see the SHC histogram for trivial and HOTIs in Supplementary Fig. 7). The Hamiltonian term responsible for the higher-order topology in inversion symmetry compounds also induces the fractional SHC. This suggests that the fractional SHC is an indicator of higher-order topology in compounds that are not topological insulators. Our approach is applicable to HOTIs having significant SOC, see Supplementary Fig. 5. Analog to topological crystalline insulators, since SOC is not explicitly required for HOTIs, the recently predicted graphdiyne 18 and graphyne 16 having negligible SOC would not be identified.
BiSe electronic structure Next, we focus on the electronic structure of a representative candidate, BiSe on the GaSe prototype (BiSe-GaSe). In Fig. 2a the BiSe-GaSe crystal structure is displayed. It features inversion, TR, and three-fold rotational symmetry. The inset presents its Brillouin zone and TR invariant momenta (TRIM) points. The fully relativistic band structure, Fig. 2b, reveals a direct bandgap of 0.46 eV at the Γ point. Its SHC clearly exhibits a constant value of σ z xy % 0:5 e 2 =h in the bandgap region. Additionally, a distinguishing characteristic of a HOTI is the appearance of (d − 2) (in this case, corner) states inside the bulk bandgap whenever the (d − 2) geometry preserves the underlying symmetries. As such, since BiSe is protected by C 3 , we constructed a triangular nanoflake with edges of approximately 50 Å, shown in Fig. 2d. Along with the nanoflake geometry, the real-space wave function projection of the six eigenstates closest to Fermi level is shown, evidencing its localization at the corners. The nanoflake eigenvalues are displayed in Fig. 2c, where the corner states are highlighted in red. The expected degeneracy is lifted owing to the interaction between corner states, as previously reported for other system 16 . A detailed discussion of interaction and localization of corner states is presented in the Supplementary information (Supplementary Figs. 2, 3, 4, and 5). Indeed, using the local SHC, Rauch and Töpler have shown that bulk and boundary states give different contributions to the SHC 26 .
Topological invariant. To definitively demonstrate BiSe higherorder topological properties we introduce an adaptation for two dimensions of the invariant defined in ref. 14 for 3D HOTIs. The symmetry-protected topological phases are driven by band inversions at K-points that preserve the symmetries protecting the bulk topology 43 . Such band inversions are characterized by a topological invariant 44,45 . For instance, in inversion symmetry (IS) insulators protected by the TR-symmetry, band inversions are characterized by the parity eigenvalues product (ξ i = ±1) for all occupied bands at all TRIM points (the Z 2 topological invariant) 24 . Similarly, in IS-HOTIs protected by C 3 and TR symmetries, band inversions at high-symmetry K-points preserving those symmetries can be divided into two groups of states according to the C 3 eigenvalue (γ C3 ¼ Àπ or ±π/3), see Fig. 2b. For each band inversion, there is an invariant given by the products of the IS eigenvalues for all occupied bands, i.e., ν 45 . Fig. 2 Two-dimensional HOTI material BiSe (in the GaSe prototype). a 2D bulk crystal structure top and side views. The inset shows the Brillouin zone high-symmetry points. b Bulk band-structure along the high symmetry lines and the spin Hall conductivity (σ z xy ) as a function of the energy. The dashed blue line represents the Fermi level. The gray shaded area encompass the bulk band gap. Green and purple dots represent states with C 3 eigenvalues γ C3 ¼ Àπ or ±π/3, respectively. Filled and unfilled dots represent + and − parity eigenvalues (ξ i = ±1) for each occupied band at the TRIM K-points. c Energy eigenvalues of a triangular flake with edges of approximately 50 Å. The localized corner modes are highlighted in red. d Flake structure with corner states charge density highlighted in blue.

Predicted 2D HOTIs
Finally, the detailed analysis is repeated for the other discovered HOTIs. The predicted materials present the structural prototypes GaSe and CH, which correspond to AB stacking centrosymmetric bilayers preserving the TR-and C 3 symmetries. The compounds are formed by elements with high intrinsic SOC, such as Bi and Pb. This results in i) small bulk bandgaps, i.e., 260, 406, 58, 398, 57, and 94 meV for BiTe, PbBr, BiTe-CH, PbCl, PbF, and HgTe, respectively; and ii) large bulk SHC. In Fig. 3 we show the crystal structure, bulk band structure and SHC, flake structure with corner states charge density, and corner states eigenvalues for the PbBr in the CH prototype. As expected, we obtain a large midgap SHC of ≈ − 0.5e 2 / h and midgap states localized at the triangular flake corners. The other three materials are shown in the Supplementary information, see Supplementary Figs. 8-12. Despite not being quantum spin Hall insulators, the midgap SHC is a constant fraction of the quantum conductivity e 2 /h. These constant values of SHC are different for compounds formed by Bi and Pb atoms. Nevertheless, the magnitude of the SHC is strongly dependent on the structure. For instance, the same BiTe composition in the structural prototypes GaSe and CH have SHC of 0.75 and 0.5 e 2 /h, respectively.

DISCUSSION
We have shown the relation between the spin Hall effect and the higher-order topological insulator (HOTI) phase. Our tight-binding model predicts a non-zero bulk midgap spin Hall conductivity (SHC). This finding allowed us to search for novel HOTIs. We combined density functional theory calculations with local effective Hamiltonians to calculate the SHC of 693 insulators. We found seven stable two-dimensional HOTIs candidates: BiSe, BiTe (in two different structures), PbF, PbBr, PbCl, and HgTe. All discovered HOTIs display metallic states in the bulk bandgap for (d − 2)-dimensional (0D) structures preserving the C 3 symmetry, being spatially localized at the corners. This confirms the existence of the HOTI states protected by the C 3 symmetry and evidence the relationship between the d-dimensional bulk SHC and the topologically protected states in (d − 2) dimensions. Despite the success of our method, predicting 7 new 2D HOTIs, one can possibly still find HOTIs in the 2D materials we investigated. These materials can be formed by lighter elements with smaller SOC and thus will result in a small SHC. While we establish numerically the relation between the SHC and the higher-order topological phase, the development of analytical expressions for the SHC in HOTI models and implementations of HOTI invariants suitable for highthroughput calculations are both desirable. For instance, the recently proposed fractional corner anomaly HOTI invariant 47 requires the density of states of both the d − 2 and d − 1 systems, posing a challenge for current high-throughput strategies. Our work advances the understanding of higher-order topological phases showing for the first time its connection with the spin Hall effect.
Note added − During the peer review process we became aware of the predicted 3D semi-conducting HOTI α-Bi 4 Br 4 21 . Our preliminary results indicate the validity of the SHC criteria also for 3D HOTIs.

DFT calculations
The exchange-correlation (XC) is described within the generalized gradient approximation (GGA) via the Perdew-Burke-Ernzerhof (PBE) functional 48 . Ionic potentials were described by projector augmented-wave (PAW) 49 pseudopotentials available in the pslibrary database 50 . Both the GGA exchange-correlation functional and the PAW ionic pseudopotentials are the same as used in the C2DB database construction. The wavefunctions and charge density cutoff energy were 40% larger than those recommended. The reciprocal space sampling was performed with a K-point density of 6.0/Å −1 for structural optimization and 12.0/Å −1 for selfconsistent calculations. The lattice parameters reported in the database were used, and a full structural optimization was performed until Hellman-Feynman forces were smaller than 0.01 eV/Å.

DATA AVAILABILITY
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request. The codes used in this work can be downloaded at QUANTUM ESPRESSO https://www.quantum-espresso.org/and PAOFLOW http://aflowlib.org/src/paoflow/.
Received: 18 September 2020; Accepted: 11 March 2021; Fig. 3 PbBr in the CH prototype. a Crystalline structure, b bulk band structure and spin Hall conductivity, c flake structure with corner states charge density, and d corresponding 0D eigenvalues.
M. Costa et al.