Sparse-Wire Diagnosis of Hinge-Localized States in Wannier Tight-Binding Hamiltonians

1. Introduction

The bulk–boundary correspondence of a conventional three-dimensional topological insulator connects a bulk invariant to gapless states on a two-dimensional surface. Higher-order topological insulators generalize this correspondence: their surfaces can be gapped while protected states occur on one-dimensional hinges or at zero-dimensional corners. Quantized electric multipole insulators established a particularly transparent route to this hierarchy [1,2]. Shortly thereafter, complementary classifications based on crystalline symmetries, inversion, and boundary mass patterns placed hinge states in a broad theoretical setting [3,4,5,6]. Experimental and first-principles studies of bismuth further showed that the concept is relevant to realistic electronic structures rather than only to solvable lattice models [7,8].

Modern material calculations commonly bridge density-functional theory and boundary-state analysis with maximally localized Wannier functions. The foundational localization and disentanglement procedures [9,10] and the Wannier90 implementation [11,12] produce compact real-space Hamiltonians that retain an accurate description of a selected energy window. Such Hamiltonians have enabled high-throughput databases [13] and specialized analysis packages. WannierTools, for example, calculates surface spectra, topological invariants, and transport-related properties [14]; Z2Pack tracks hybrid Wannier charge centres [15]; and Kwant provides a general framework for quantum transport and open-system tight-binding calculations [16]. Recent programs have expanded the ecosystem towards model fitting, magnetic exchange, excitons, and optoelectronics [17,18,19]. These packages evaluate boundary spectra for a geometry specified by the user; none screens the lattice axes automatically or reports corner-resolved spectral weights as a reproducible classification.

Nevertheless, obtaining a Wannier Hamiltonian does not make hinge-state screening automatic. A researcher must identify an appropriate bulk energy gap, decide which pair of surfaces should be opened, build a wire that preserves translation along a possible hinge, compute only a small spectral window of a potentially large matrix, and determine whether in-gap eigenstates occupy corners rather than the complete boundary. These tasks are simple in principle but vulnerable to inconsistent phase conventions, omitted Wigner–Seitz degeneracies, coarse Brillouin-zone sampling, inappropriate surface terminations, and finite-size localization unrelated to a higher-order bulk invariant.

HingeStateDetector was developed to make this intermediate step reproducible. The program accepts the two files most often available after a Wannier calculation, wannier90_hr.dat and POSCAR, and produces both human-readable plots and machine-readable evidence. It screens each primitive lattice axis as a possible wire direction, evaluates spectral and spatial criteria, checks Hermiticity, and records every numerical choice in JSON. Its automatic LIKELY_HINGE_STATE label identifies boundary-localized spectral evidence for subsequent convergence and bulk-topology tests.

This article describes the mathematical convention, implementation, output model, and verification of HingeStateDetector. The first example is an exact positive control derived from the three-dimensional chiral hinge model of Benalcazar, Bernevig, and Hughes [2]. The second is the experimentally established HOTI $\alpha$-Bi₄Br₄, for which theory, ARPES, and spatially resolved measurements support chain-parallel one-dimensional hinge states [20,21,22]. The third uses the official 30-orbital spinful Bi₂Se₃ example distributed with WannierTools. Bi₂Se₃ is a canonical strong topological insulator [23,24,25], so it is an instructive test of both realistic model size and interpretive discipline. The resulting workflow is summarized in Figure 1.

2. Model and Numerical Method

2.1. Real-Space Wannier Representation

A Wannier90 Hamiltonian is stored as matrices $H_{mn}\left(\mathbf{R}\right)$ between Wannier orbital n in a reference cell and orbital m in the cell displaced by integer lattice vector $\mathbf{R}=R_1{\mathbf{a}}_1+R_2{\mathbf{a}}_2+R_3{\mathbf{a}}_3$. With Wannier90’s Wigner–Seitz degeneracy $d_{\mathbf{R}}$, the Bloch Hamiltonian used here is

$$H_{mn}\left(\mathbf{k}\right)=\sum _{\mathbf{R}}{\displaystyle \frac{H_{mn}\left(\mathbf{R}\right)}{d_{\mathbf{R}}}}exp\left(2\pi i\mathbf{k}·\mathbf{R}\right),$$

(1)

where the components of $\mathbf{k}$ are fractional reciprocal coordinates. The parser reads the number of Wannier functions, the number and ordered list of lattice vectors, all degeneracies, and the complex matrix elements. Repeated entries are accumulated. Matrix elements smaller than a user-controlled cutoff can be removed, although the default preserves the supplied model.

The POSCAR provides the lattice matrix and is used to assign physical coordinates to repeated cells. Cartesian and direct atomic coordinates are accepted, including a selective-dynamics line, but atomic coordinates do not enter the Hamiltonian construction. This distinction is essential: hr.dat contains an orbital-indexed Hamiltonian but not the spatial coordinate of each Wannier centre, and POSCAR contains atomic positions but not an unambiguous orbital-to-atom map. Consequently, the current implementation resolves probability at the repeated-cell level. An orbital-resolved assignment would require additional information such as wannier90_centres.xyz or a user-supplied orbital map.

Before any boundary calculation, the program samples Equation (1) and evaluates

$${ϵ}_v^{max}=\underset{\mathbf{k}}{max}{ϵ}_v\left(\mathbf{k}\right),{ϵ}_{v+1}^{min}=\underset{\mathbf{k}}{min}{ϵ}_{v+1}\left(\mathbf{k}\right),$$

(2)

for all adjacent band indices v. A sampled global gap exists when ${ϵ}_{v+1}^{min}-{ϵ}_v^{max}>g_{min}$. Among these gaps, the largest interval containing the requested reference energy is preferred. If none contains it, the code selects the nearest interval using a gap-normalized distance. This rule is useful because hr.dat does not encode the Fermi energy or the electron count. It also makes the limitation explicit: a coarse bulk mesh can miss a narrow indirect overlap, so gap convergence must precede a strong physical claim.

2.2. Mixed-Boundary Wire Hamiltonian

For a trial periodic axis p, the other two axes u and v are truncated to $N_u$ and $N_v$ complete unit cells. A basis state is indexed by $(i_u,i_v,m)$, and the matrix dimension is

$$D=N_u{N}_v{N}_{\mathrm{W}},$$

(3)

where $N_{\mathrm{W}}$ is the number of Wannier orbitals. A hopping with displacement $(R_u,R_v,R_p)$ contributes only if both $i_u+R_u$ and $i_v+R_v$ remain within the finite cross section. Translation along p is retained through the phase factor $exp\left(2\pi i{k}_p{R}_p\right)$. Thus,

$$\begin{array}{c}\langle j_u,j_v,m\left|H_p\left(k_p\right)\right|i_u,i_v,n\rangle \hfill \\ =\sum _{\mathbf{R}}{\displaystyle \frac{H_{mn}\left(\mathbf{R}\right)}{d_{\mathbf{R}}}}{e}^{2\pi i{k}_p{R}_p}\\ \hfill \times {\delta }_{j_u,i_u+R_u}{\delta }_{j_v,i_v+R_v}.\end{array}$$

(4)

The implementation assembles Equation (4) directly as a SciPy sparse matrix. Sparse diagonalization is central for realistic models because the full eigenspectrum scales cubically with D and stores many states far from the target gap. Instead, HingeStateDetector applies ARPACK through SciPy’s shift–invert interface and requests a user-defined number of eigenpairs around the bulk-gap midpoint. NumPy and SciPy supply the array and eigensolver infrastructure [26,27]; Matplotlib produces static figures [28].

The real-space hopping convention, boundary clipping, and phase factor are tested together by evaluating

$${\delta }_H\left(k_p\right)={\displaystyle \frac{\parallel H_p\left(k_p\right)-H_p^{†}\left(k_p\right)\parallel }{max[\parallel H_p\left(k_p\right)\parallel ,{ϵ}_{\mathrm{mach}}]}}.$$

(5)

The matrix is numerically symmetrized before diagonalization. If the maximum ${\delta }_H$ exceeds the requested tolerance, a positive-looking spectrum is downgraded to an inconclusive result rather than silently accepted.

2.3. Spatial Weights and Candidate Criterion

For a normalized eigenvector ${\psi }_{\alpha k_p}$, the probability in cross-sectional cell $(i_u,i_v)$ is

$$P_{\alpha k_p}(i_u,i_v)={\displaystyle \sum _{m=1}^{N_{\mathrm{W}}}}{\left|{\psi }_{\alpha k_p}(i_u,i_v,m)\right|}^2.$$

(6)

The cross section is partitioned as illustrated in Figure 2. Corner masks contain square regions of ${\ell }_c$ complete cells at the four corners. The edge mask contains the remaining cells within ${\ell }_c$ layers of the boundary, and the bulk mask contains the interior. The corresponding weights $W_c$, $W_e$, and $W_b$ are sums of Equation (6); they obey $W_c+W_e+W_b=1$ up to rounding.

A state is a candidate if its energy lies inside the sampled bulk gap after a fractional margin has been removed, $W_c$ exceeds the corner threshold, $W_c>W_e$, and $W_b$ is below the bulk threshold. A periodic direction receives the LIKELY_HINGE_STATE label if at least a specified fraction of sampled $k_p$ points contains one or more candidates. These thresholds are operational screening parameters, and all values are written to the output report so that the classification can be reproduced and varied systematically.

The cell-resolved LDOS provides a second, visually direct representation:

$$\rho (i_u,i_v,E)={\displaystyle \frac{1}{N_k}}\sum _{k_p,\alpha }{P}_{\alpha k_p}(i_u,i_v){\displaystyle \frac{\eta /\pi }{{(E-E_{\alpha k_p})}^2+{\eta }^2}}.$$

(7)

Here $\eta$ is the Lorentzian half-width. HingeStateDetector stores the complete three-dimensional array $\rho (i_u,i_v,E)$ rather than only the plotted subset. In all LDOS figures, marker area follows the bounded linear mapping $s=s_{min}+(s_{max}-s_{min})(\rho /{\rho }_{max})$, while color is independently mapped with the jet palette. The two visual channels make high-LDOS points substantially larger than the background, while the NPZ array remains the quantitative source.

3. Computational Implementation and Workflow

Table 1 summarizes the computational scope. The principal executable, hinge_detector.py, is intentionally self-contained. Data classes represent the parsed structure, hopping model, selected gap, and result from one orientation. The parser and dense bulk Hamiltonian are separated from the sparse WireBuilder; the latter precomputes boundary-compatible index maps for every hopping block and reuses them over the one-dimensional momentum mesh. Independent functions implement region masks, state weights, LDOS accumulation, plotting, and JSON serialization. This separation makes it possible to test the marker scaling and topological/trivial demonstration models without running a large material calculation.

The minimum invocation is

python3 hinge_detector.py \

  --hr wannier90_hr.dat --poscar POSCAR \

  --output hinge_results --periodic-axis auto \

  --size 8 8 --verbose

Automatic orientation mode executes three wire calculations, one for each primitive lattice direction. A recommended production workflow begins with a $6\times 6$ or $8\times 8$ cross section and modest k mesh, identifies the most promising direction, then repeats that direction with increasing cross section, bulk mesh, one-dimensional momentum mesh, and requested eigenpair count. The companion run_convergence.py automates this two-stage procedure and writes a consolidated convergence.json.

Each run produces a summary.json containing the overall status, selected axis, sampled gap, model dimensions, complete command parameters, orientation reports, runtime, and a textual interpretation. Each orientation directory contains report.json, hinge_data.npz, a corner-weighted band plot, a two-dimensional LDOS slice, and a static three-dimensional LDOS point cloud. If Plotly is installed, an interactive HTML point cloud is also written. Exit status 0 denotes at least one likely orientation, status 1 denotes a completed negative screen, status 2 means that no sampled global bulk gap was found, and status 3 denotes invalid input or a numerical error. In batch workflows, status 1 is therefore a scientific result rather than a crashed job. The source code, tests, input archives, run parameters, and numerical arrays underlying the reported figures are available at the following DOI: 10.5281/zenodo.21778263.

Table 2 lists the most consequential controls. The defaults are chosen for exploratory screening and should be converged for quantitative material calculations. In particular, the physical localization length need not equal one cell. Increasing corner-layers is appropriate for wider wires, but the mask should occupy a decreasing fraction of the cross-sectional area during a convergence study. Likewise, num-eigs must be large enough to include all bands inside and near the selected gap. A too-small request can make a complex surface manifold look artificially sparse.

4. Verification with a Chiral Hinge Model

4.1. Positive-Control Construction

The positive control implements the four-orbital chiral hinge Hamiltonian associated with Equation (132) of Benalcazar, Bernevig, and Hughes [2]. The chosen parameters place each constant-$k_3$ two-dimensional subsystem in the appropriate quadrupole regime over part of the Brillouin zone, producing chiral branches at the four hinges of a wire open along ${\mathbf{a}}_1$ and ${\mathbf{a}}_2$. A generator writes the model in the standard Wannier90 text format together with a cubic POSCAR; the detector therefore exercises the same input path used for a first-principles Hamiltonian. The archived metadata identify the publication and DOI from which the model was transcribed.

The bulk mesh identifies a gap centred at zero. With a $6\times 6$ cross section, 31 momenta, and 24 eigenpairs, automatic orientation screening produces the result shown in Figure 3. Directions ${\mathbf{a}}_1$ and ${\mathbf{a}}_2$ are inconclusive: some states approach a single-cell corner weight of 0.468, but none satisfies the combined energy and spatial criteria. For ${\mathbf{a}}_3$ periodic, 35.5% of momenta contain candidates, the maximum corner weight is unity, and the mean candidate weights are $(W_c,W_e,W_b)=(0.859,0.127,0.013)$. The maximum relative Hermiticity error is $5.8\times {10}^{-17}$.

The denser $12\times 12$ calculation in Figure 4 resolves four branches crossing at zero energy. Their red coloring corresponds to nearly unit corner weight around the crossing, while states merging into the bulk continua evolve continuously towards smaller corner weight. The LDOS at zero energy is concentrated at all four cross-sectional corners. Because LDOS marker area increases with intensity, the contrast does not depend only on subtle color differences. This simultaneous spectral flow and four-corner localization is stronger evidence than either an in-gap energy or an LDOS peak alone.

4.2. Finite-Size and Performance Behavior

Cross sections from $6\times 6$ to $16\times 16$ remain positive, as summarized in Table 3 and Figure 5. The mean candidate corner weight approaches approximately 0.916, while the mean interior weight falls below 0.01 for the two largest calculations. The four eigenvalues nearest zero at $k_3=0$ have a span of $1.76\times {10}^{-16}$ for $N=6$ and $7.40\times {10}^{-32}$ for $N=12$ and 16. The exceedingly small values reflect the exact structure and symmetry of this benchmark rather than a general promise of machine-zero splitting in materials.

The cost is not a clean power law because the calculations deliberately used different values of $N_k$ and requested eigenpairs, and because sparse LU fill-in depends on connectivity. The $16\times 16$ run requested 64 eigenpairs at 101 momenta and took approximately 2822 s, whereas the $12\times 12$ calculation took 28 s. These measurements support the staged strategy: orientation screening should use small wires and a narrow spectral window, followed by targeted convergence only for promising directions. For a material with 30 Wannier orbitals, Equation (3) gives $D=1920$ already at $8\times 8$; a blind three-axis large-size scan would therefore be wasteful.

5. Real-Material Positive Control: $\alpha$-Bi₄Br₄

The van der Waals material $\alpha$-Bi₄Br₄ provides the real-material positive control required to test HingeStateDetector beyond an ideal lattice model. Its HOTI identity is established independently of the present calculation. First-principles theory identifies a rotational-symmetry-protected three-dimensional crystalline phase with one-dimensional hinge states along the quasi-1D chains [20]. ARPES measurements subsequently showed that the stacking of quantum-spin-Hall building blocks realizes a three-dimensional HOTI [21]. More recent spatially and spin-resolved spectroscopy found gapped surface states together with in-gap states connecting the surface valence and conduction bands, providing additional hinge-state evidence [22]. Together, these studies establish $\alpha$-Bi₄Br₄ as a real-material HOTI positive control independently of the present detector calculation.

For reproducibility, we use the NIST JARVIS-WTBH entry JVASP-32759 [13]. Its reduced formula is BiBr and its space group is $C2/m$ (No. 12); the supplied POSCAR contains eight Bi and eight Br atoms, corresponding to two Bi₄Br₄ formula units. The official archive contains the mutually matched POSCAR, wannier90_hr.dat, Wannier90 input and output, and the JARVIS validation record. The Hamiltonian has 128 spinor Wannier orbitals constructed from Bi and Br $s,p$ projections and 85 nonzero lattice-vector blocks. We use its mesh Fermi energy, 2.95712475 eV, as the reference energy. The JARVIS Hamiltonian is an independent high-throughput Wannierization of the same material and need not be numerically identical to the models used in Refs. [20,21,22].

An exploratory $3\times 3$ calculation scanned all primitive lattice axes with 11 momenta and 16 eigenpairs. The best direction was lattice axis 2, which is the quasi-one-dimensional chain and expected hinge direction. We then retained this direction for $4\times 4$ and $5\times 5$ wires using 21 momenta and 24 shift–invert eigenpairs. Increasing the bulk mesh from $4^3$ to $5^3$ changed the sampled gap only from 0.169929 to 0.167591 eV, with the final interval $[2.92857,3.09616]$ eV. The $5\times 5$ sparse Hamiltonian has dimension 3200 and remained numerically Hermitian to $6.28\times {10}^{-17}$.

Figure 6 summarizes the resulting positive control. At $5\times 5$, 52.4% of the sampled chain momenta contain at least one candidate, the maximum corner weight is 0.718, and the mean candidate weights are $(W_c,W_e,W_b)=(0.626,0.359,0.0149)$. The LDOS at 3.069 eV is concentrated at the two diagonally opposite termination corners selected by this complete unit-cell cut. Marker areas are linearly proportional to normalized LDOS and are capped so that neighboring cells remain visible. The interior weight stays below 0.026 throughout the $3\times 3$–$5\times 5$ sequence, while the corner weight remains larger than the remaining edge weight. These results reproduce the expected combination of a chain-parallel in-gap spectrum and cross-sectional hinge localization for a material whose higher-order topology is known independently.

The complete-cell boundary of the JARVIS model need not be atomically identical to the cleaved surfaces studied experimentally, and the finite wire contains additional surface-derived levels. Nevertheless, the selected direction, in-gap spectrum, and cross-sectional localization agree with the established hinge phenomenology of $\alpha$-Bi₄Br₄, providing a real-material validation of the detector.

6. Real-Material Stress Test: Bi₂Se₃

Bi₂Se₃ is a useful real-material input for two reasons. Its strong topological-insulator band structure is established by first-principles theory, photoemission, and low-energy modeling [23,24,25], while its Wannier model is publicly distributed as an official WannierTools example [14]. The archived POSCAR and wannier90_hr.dat contain a rhombohedral cell and a 30-orbital spinful basis built from Bi and Se p orbitals. The Hamiltonian has 547 nonzero lattice-vector blocks. We use the tutorial reference energy 4.4195 eV and retain the files without re-fitting or modifying the hoppings. The source URLs and raw input archive are preserved with the manuscript.

A $4^3$ bulk mesh for initial screening gave a sampled gap of 0.315 eV. The first two primitive-axis wires were inconclusive, while the third was positive. We then calculated $6\times 6$ and $8\times 8$ cross sections along this direction. For the $8\times 8$ wire, an $8^3$ bulk mesh gives $[4.26337,4.57866]$ eV. Thirty-one momenta and 32 shift–invert eigenpairs produce the spectrum in Figure 7. The mean candidate weights are $(0.657,0.175,0.168)$ and the maximum corner weight is 0.701. The maximum wire Hermiticity error is $6.4\times {10}^{-17}$.

This case tests the specificity of the geometric classifier. The candidate bands remain near the lower bulk-gap edge instead of displaying the unambiguous gap-spanning spectral flow of the BBH control. The representative LDOS is dominated by one corner, and the interior weight increases from 0.131 at $6\times 6$ to 0.168 at $8\times 8$ for the fixed one-cell mask. The threshold-based classifier has therefore found geometrically localized wire states without the spectral signature of the positive control. Because Bi₂Se₃ is canonically described as a first-order strong topological insulator, surface-state quantization, termination-specific potentials, and unit-cell grouping are plausible sources of this signal.

We therefore interpret this output as a boundary-localization signal rather than a higher-order assignment. Establishing such an assignment would require surface and termination calculations, larger-wire scaling, and an independent crystalline invariant or boundary-mass analysis. The example simultaneously verifies the parser and sparse solver for a realistic spin-orbit-coupled Hamiltonian.

7. Robustness, Applicability, and Limitations

The two mandatory files make HingeStateDetector broadly reusable. The hr.dat file must follow the standard Wannier90 real-space format, including degeneracies, and the POSCAR lattice must correspond to the same primitive or conventional cell used for the Hamiltonian. The target energy must lie near the intended insulating gap, and the system must have a sampled global gap. Metals and semimetals are outside the current classifier’s assumptions.

Only wires parallel to a supplied lattice vector are generated automatically. Many physically relevant hinges are defined by non-primitive Miller planes, surface reconstructions, cleavage choices, or enlarged surface cells. A general transformation to an integer supercell basis and a termination builder would be needed for these cases. Even for an axis-aligned wire, cutting complete cells selects one particular termination. Adding or removing a surface layer can alter trivial boundary bands without changing the bulk Hamiltonian, so termination robustness is a physical test rather than a software detail.

Cell-level localization is another deliberate compromise. Summing over all Wannier functions in a repeated cell is invariant under unitary rotations among those orbitals and is sufficient for localization at well-separated wire corners. It cannot reveal which atom, orbital, or Wannier centre carries the state, and it can exaggerate “corner” localization when the cross section contains only a few cells. For example, the four one-cell corner masks occupy 44% of a $3\times 3$ cross section before edges are separately counted. Production calculations therefore require larger cross sections and should monitor both the absolute localization length and the fractional mask area.

The bulk-gap estimator is based on a regular mesh rather than an optimizer over the Brillouin zone. This choice is transparent and parallelizable, but a small indirect overlap can be missed. Increasing bulk-mesh, checking high-symmetry paths, and comparing the Wannier bands against the parent first-principles calculation are necessary safeguards. The latter calculation may originate from VASP [29], Quantum ESPRESSO [30], or another electronic-structure package.

The classifier intentionally combines several imperfect indicators. Corner weight rejects states distributed along complete surfaces, the edge comparison rejects states spread around the perimeter, the interior threshold penalizes bulk-like states, and the momentum-fraction test suppresses isolated accidental levels. No fixed threshold is basis-, size-, and material-independent. The machine-readable arrays make it possible to replace these defaults with continuous diagnostics, localization-length fits, inverse participation ratios, clustering, or supervised ranking without repeating the expensive eigensolutions.

The boundary diagnostics complement the bulk invariants and surface observables available from programs such as Z2Pack and WannierTools [14,15]. Combining these layers yields a reproducible topological-materials workflow, as illustrated by the analytically known BBH benchmark and the independently established Bi₄Br₄ positive control.

8. Conclusions

HingeStateDetector converts a common pair of electronic-structure outputs into a reproducible search for hinge-localized wire states. The implementation handles Wannier90 degeneracies, selects a sampled bulk gap near a user-specified energy, assembles three mixed-boundary sparse Hamiltonian families, solves only the relevant spectral window, and quantifies corner, edge, and interior probability. Static and interactive LDOS products complement the band spectrum, while JSON and NPZ outputs preserve the complete numerical evidence.

The chiral BBH benchmark reproduces the expected axis-selective, four-corner, gap-spanning hinge modes and remains positive through a $16\times 16$ size sequence. The 128-orbital JARVIS-WTBH model of the experimentally established HOTI $\alpha$-Bi₄Br₄ selects the chain direction and retains corner-dominated candidates through a $3\times 3$–$5\times 5$ sequence. The official WannierTools Bi₂Se₃ model demonstrates that the same code also handles a realistic 30-orbital spin-orbit Hamiltonian and provides a specificity stress test for the geometric classifier.

Future development should add arbitrary integer surface transformations, explicit surface termination editing, use of Wannier-centre coordinates, orbital- and atom-resolved LDOS, iterative Green-function alternatives to finite wires, symmetry annotations, and interfaces to bulk-invariant calculations. In its present form, HingeStateDetector provides a transparent and reproducible boundary-localization layer that identifies hinge-like spectral weight and retains the information needed for subsequent topological analysis.