BerryEasy: a GPU enabled python package for diagnosis of nth-order and spin-resolved topology in the presence of fields and effects

Multiple software packages currently exist for the computation of bulk topological invariants in both idealized tight-binding models and realistic Wannier tight-binding models derived from density functional theory. Currently, only one package is capable of computing nested Wilson loops and spin-resolved Wilson loops. These state-of-the-art techniques are vital for accurate analysis of band topology. In this paper we introduce BerryEasy, a python package harnessing the speed of graphical processing units to allow for efficient topological analysis of supercells in the presence of disorder and impurities. Moreover, the BerryEasy package has built-in functionality to accommodate use of realistic many-band tight-binding models derived from first-principles.

Crucial to this progress has been the development of powerful community codes for the construction and analysis of tight-binding models.
These programs include Z2Pack [6], WannierTools [7], Wannier90 [8], PythTB [9], Kwant [10], WannierBerri [11], and Pybinding [12] among others.While sharing many features, each package has individual strengths.Until recently a common issue was that none of these packages supported built-in functionality for computing state-of-the-art topological diagnostics including the nested Wilson loop [13,14] and the spinresolved Wilson loop [15,16].Both computations are critical for a comprehensive analysis of band topology.Recently, an auxiliary code was constructed to work in conjunction with the PythTB package for the computation of these quantities [16].However, as mentioned, each code has separate strengths.PythTB is a terrific option for investigating band topology, however attempting to account for the presence of disorder/impurities in supercells can be cumbersome.Furthermore, it is not an ideal program for investigating transport or edge spectral density of many-band models.These are situations in which the WannierTools package has excelled due to its speed advantage through use of Fortran.However, as a Fortran based program it is challenging to manipulate, and does not offer access to the spinresolved Wilson loop or nested Wilson loop at the moment.This situation has motivated the development of a new auxiliary package, BerryEasy, which works in tandem with the PyBinding software for computation of topological properties.PyBinding is designed with a python front end and C++ back end, offering a balance between speed and ease of use.It is thus ideal for investigating the effects of external fields, disorder and impurities on transport and density of states, but to this point it has lacked functionality for diagnosing band topology.In this paper, we detail how BerryEasy offers a simple interface for the computation of the Wilson loop [6,[17][18][19][20][21][22][23][24][25][26], nested Wilson loop [13,14], and spin-resolved Wilson loop [15,16] for models defined in PyBinding.Given its popularity, we note that BerryEasy can additionally be interfaced with Kwant, although not all functionalities available via PyBinding are available with Kwant at this point.For clarity we focus on interfacing of BerryEasy with PyBinding in this work.For documentation and a tutorial on how to interface BerryEasy with Kwant please see the accompanying webpage [27].Importantly, BerryEasy can be run using a CPU or GPU.Operating the program on a GPU is made possible by CuPy [28] and decreases the computation time dramatically.This is significant as computation of topological invariants is generally expensive, particularly in supercells, due to the need for exact diagonalization of a discretized Hamiltonian at many points in reciprocal space.In the BerryEasy workflow, the Hamiltonian is rapidly built by Py-Binding's C++ back-end and the subsequent analysis by BerryEasy on GPUs enjoys significant speed advantages.As a result it is possible to efficiently investigate the fate of band-topology upon introduction of disorder, impurities, external fields and other physical situations for which finite-size effects and disorder averages must be considered.Finally, we incorporate functionality for defining a PyBinding lattice using a Wannier tight-binding model created using the Wannier90 software package.This extends all functionalities in idealized models to realistic sys-arXiv:2312.13051v2[cond-mat.mtrl-sci]13 Jun 2024 tems.
The remainder of this work is organized as follows.In Sec.(II), Wannier center charges and their use to determine topological properties is defined.The numerical approach used by BerryEasy for their determination is also presented.In Sec.(III), nested Wannier center charges are defined along with their computation in BerryEasy.In Sec.(IV), spinresolved Wannier center charges are defined along with details of their computation in BerryEasy.Finally, before concluding in Sec.(VI), details of how to analyze a realistic tight-binding model produced by Wannier90 are provided in Sec.(V).In each section multiple examples are provided and in select cases these examples are accompanied by code snippets.These code snippets are offered to illustrate the ease with which BerryEasy can be implemented and do not represent the full code required to produce the results shown in the main body.Tutorials containing all necessary code to reproduce the data and figures seen in this work are available on the accompanying webpage [27].

II. WANNIER CENTER CHARGES
Analysis of Wannier center charges (WCCs) as calculated via Wilson loops has proven to be a fundamental building block in diagnosis of band topology [6,[17][18][19][20][21][22][23][24][25][26].For clarity, we provide a brief overview of the formalism for computation of WCCs in this section and provide examples for implementation of the computation in the BerryEasy package.Mathematically, the Wannier center charge for band n of a one-dimensional system is defined as [6], (1) where a is the lattice constant, |u n (k)⟩ is the Bloch wavefunction corresponding to band n and A n (k) is the Berry gauge connection.We note that this computation requires that periodic boundary conditions be enforced and path-ordering must be explicitly enforced if the Berry gauge connection is non-Abelian.When the Berry gauge connection is Abelian, computation of the Wilson loop is equivalent to computation of the one-dimensional winding number which, in the Altland-Zirnbauer table [29][30][31][32][33][34], is Z classified for class AIII insulators, such as the famous Su-Schreiffer-Heeger model [35].If the Berry gauge connection is non-Abelian, as is the case for spinful insulators in class AII, the Wilson loop is instead Z 2 valued.

A. Numerical Implementation:
In the BerryEasy package Wannier center charges are computed by discretizing the integral in eq. ( 1) and explicitly imposing path-ordering.This is accomplished via the formula, where F n,k+N ∆k = ⟨u nk+∆k |u nk ⟩, ∆k = 2π/(aN ) where a is the lattice constant, and |u nk ⟩ is the Bloch function of band n at reciprocal space coordinate k.The accuracy of the integral is determined by the parameter, N .This is an integer value which determines the number of discretized segments used in computing the integral and is set by the user in BerryEasy.We emphasize that periodic boundary conditions must be imposed along the path of integration, i.e. ψ(k) = ψ(k + 2π/a).
When considering supercells, including those which incorporate disorder and impurities, the above implementation remains valid provided we work in the language of twisted boundary conditions (TBCs).For a one-dimensional supercell of size L with a unit lattice constant, TBCs are implemented as ψ(x) = e iθx ψ(x + L) where 0 ≤ θ x < 2π.Following Ref. [36], the single particle wavefunctions can then be Fourier transformed such that they are a function of the discrete momenta, k = 2πn x /L + θ x /L, where 0 ≤ n x < L.
In BerryEasy, when interfacing with PyBinding we fix θ x = 0 by setting perioidc boundary conditions with a period equivalent to the supercell size and compute eq. ( 2) as a function of k = 2πn x /L.As such, when working with a supercell of size L, the reciprocal lattice vectors must be properly defined as 2π → 2π/L.Upon making this definition, the same numerical algorithms are implemented to rapidly compute the Wannier charge centers for the supercell.

B. Hybrid Wannier center charges
Importantly, the WCC formalism has been extended to allow for the determination of bulk topological invarants in higher physical dimensions through computation of hyrbid WCCs.This is simply the computation of WCCs as a function of a transverse momenta.For example, in a two-dimensional system if we wish to compute the WCCs along k x as a function of transverse momenta k y , we modify eq.(1) to FIG. 1: Wannier center charge (WCC) spectra for occupied valence band of eq. ( 4).Spectral flow of WCCs indicates a unit Chern number.
the form, We emphasize that in the above formula, both directions must support periodic boundary conditions.Indeed, to utilize WCCs for topological classification periodic boundary conditions must be imposed along all directions in the d-dimensional system we wish to classify.This is distinct from stating that pe-riodic boundary conditions must be imposed along all directions.For example, if we wish to analyze a two-dimensional x − y plane embedded in a threedimensional system, the third direction (z) boundary conditions need not be periodic.By analyzing spectral flow of the WCCs as a function of the transverse momenta, the presence or absence of topological invariants can be determined [6].
As an illustrative example, we consider a model of a Chern insulator on a square lattice [37][38][39][40].The Bloch Hamiltonian takes the form, ) where t has units of energy, σ j=1,2,3 correspond to the three Pauli matrices and the lattice constant has been set to unity.This model supports a non-trivial Chern number, C = 1 which can be computed as, where Alternatively, the first Chern number can be computed using hybrid WCCs as, where the sum is over the occupied bands.Upon defining the Bloch Hamiltonian in the PyBinding package, we execute the calculation in the formalism of BerryEasy as, import BerryEasy a s be vec=lambda t1 , t 2 : [ t1 , t2 , 0 ] #i n t e g a t i o n a l o n g kx a s f u n c t .o f ky ds =100 #d i s c r e t i z i n g i n t e g r a t i o n i n t o 100 s t e p s ds2 =100 #d We clarify that in the above code vec = lambda t1, t2 : [t1, t2, 0] is a Python lambda function defined by the user for determining the twodimensional plane in which to perform computation of WCCs.Namely, t1 specifies the direction of integration in units of the reciprocal lattice vectors and t2 specifies the direction of the transverse momenta which will be varied to search for spectral flow of the WCCs.Both directions must be fixed to support periodic boundary conditions within PyBinding.In addition, the rvec term specifies the reciprocal lattice vectors and must be a 3 × 3 array.This is the case even if the system is one or two-dimensional.In one dimensional systems the second and third row are discarded by BerryEasy and in two-dimensional systems the third row is discarded.For clarity the above sample code shows only the final step in computing the WCCs via BerryEasy.
In addition, the tight-binding model, eq. ( 4), must be defined in the formalism of PyBinding.In the It is now important to consider a spinful system sup-porting time-reversal symmetry T with T 2 = −1.In this case the system belongs to class AII and supports a Z 2 topological invariant under the ten-fold classification scheme [29].Importantly, it was shown that a direct connection can be made between the Z 2 index [41,42] and the WCC spectra in Ref. [20].We briefly summarize that in a time-reversal symmetric system each Kramers pair is composed of two eigenstates which admit equal and opposite Chern numbers, As a result, a non-trivial Z 2 index is indicated by the presence of a WCC spectra simultaneously detailing a smooth interpolation from x(0) = 0(a) to x = a(0).
As an example of this behavior we will modify the Bloch Hamiltonian in eq. ( 4) to the form, where σ 0,1,2,3 (τ 0,1,2,3 ) are the 2 × 2 identity matrix and three Pauli matrices respectively, operating on the spin (orbital) indices.This is the celebrated Bernevig-Hughes-Zhang model [43,44]  Again, details of the code for defining eq. ( 7) as a PyBinding instance are available in the tutorial on the accompanying webpage [27].The results in Fig.
(2) demonstrate the resulting fully connected WCC spectra indicating a non-trivial Z 2 index.This formalism can be easily extended to three-dimensional systems as computation of the weak and strong Z 2 indices is accomplished via computation of the Z 2 index in each high-symmetry plane [41].

C. Example: Topological Insulator with vacancies
As discussed in the introduction, a benefit of the PyBinding package is the capability to include fields and defects with ease.Accounting for defects is critical to the investigation of any realistic experimental setup but computing topological invariants in the presence of defects such as vacancies can be exceedingly challenging.In this example, we consider the BHZ model of a two-dimensional quantum spin-Hall insulator with a non-trivial Z 2 index [43,44].We will create a supercell of 21 × 21 unit cells.We then remove lattice sites contained within a radius r cen- tered at the origin.As r is increased the Z 2 index will be computed.For a Jupyter notebook containing full details of this example please consult Ref. [27].
Vacancies can be included in lattice tight-binding models simply in the PyBinding package through the use of the site-state modifier function, the complete code to construct the model is given on the webpage [27].As stated we will consider a 21 × 21 supercell of the model given in eq.(7).A plot of the supercell and bulk density of states for varying values of r is shown in We note as r is increased the density of states within the bulk gap begins to become populated.This can be understood in a straightforward manner: increasing r causes the vacancies to appear as a boundary with edge states.The limited size of the boundary means the hybridization of these edge states maintains an energetic gap, protecting the bulk topology.
To test this, the Z 2 index can be computed via the Wannier center charge spectra of the supercell in the presence of the vacancies.The results of this computation are available in Figs.

III. NESTED WCC SPECTRA
There are currently limited community codes which accommodate computation of the nested Wilson loop.The nested Wilson loop or nested WCC spectra was introduced in Refs.[13,14,45] as a method for computing the WCC spectra of the Wannier Hamiltonian.As an example, the Wannier Hamiltonian for a two-dimensional system takes the form, H W 1 (k 2 ), computed as, where P indicates path-ordering and A 1 (k 2 ) is the Berry gauge connection for occupied bands.We therefore note that the Wannier Hamiltonian can be constructed numerically in BerryEasy by replacing the right hand side of eq. ( 8) with the dicretized formulation of the integral in eq.(2).It is also important to emphasize that this computation requires periodic boundary conditions are imposed along k 1 and k 2 .While claims that the Wannier Hamiltonian is equivalent to the surface Hamiltonian have been made to justify the use of the nested Wilson loop in diagnosis of higher-order topological insulators where the surfaces can take the form of lower dimensional topological insulators, this remains unresolved with counter examples appearing in the literature [46].Neverthe-less, the nested Wilson loop has proven incredibly useful and powerful in identifying higher-order topological insulators (HOTIs).While the nested Wilson loop formalism is provided in detail in Ref. [13], here we provide an example for its implementation using the BerryEasy package.For clarity we will utilize the well studied chiral HOTI Bloch Hamiltonian [14] which takes the form, If ∆ 0 is set to zero the model reduces to that of a strong topological insulator.Furthermore, setting ∆ 0 = 0 and k z = 0 this model reduces to the BHZ two-dimensional QSH insulator investigated previously.However, by setting ∆ 0 = 1, the surface states are gapped on the (100) and (010) surfaces.As a result, the WCC spectra is gapped, as seen in

IV. SPIN-RESOLVED WCC SPECTRA
While first introduced by Prodan[15] over a decade previously, the spin-resolved Chern number (C s ) has found renewed importance in the diagnosis of band topology for spinful higher-order and fragile topological insulators as the crystalline symmetry preserving perturbations which serve to gap the edge states and bring about higher-order topology often violate the spin-rotation symmetry, forcing the use of this method to diagnosis the bulk invariant.The spin-resolved WCC spectra is detailed at length in Ref. [15,16], where a python package for its computation within the PythTB package is provided.This is a terrific community code.Nevertheless, we have found a need for implementation of the spin-resolved Wilson loop in PyBinding given the ease with which realistic fields and perturbations can be included.As an example, we again utilize the Bloch Hamiltonian given in eq. ( 9).As seen in Fig. (4(a)), the WCC spectra is gapped.Computation of the spin-resolved WCC spectra requires defining the projected spin operator (PSO), P (k)ŝP (k), where P (k) is the projector onto occupied bands and ŝ is a chosen spinquantization axis.In the absence of spin-orbit coupling the eigenvalues of the PSO are fixed as +1 and −1, corresponding to the spin-up and spin-down eigenstates respectively.Since we have introduced spin-orbit coupling in eq. ( 9), the eigenvalues can adiabatically deviate from ±1.Nevertheless, a gap in the eigenvalue spectra of the PSO remains when selecting ŝ = s z = σ 3 ⊗ τ 3 , allowing for calculation of the WCC spectra for the negative eigenstates of the PSO, x− (k y ), as detailed in Lin et.al [16].Such The study of disorder in topological quantum matter is of supreme importance yet it presents a computational challenge [40,[47][48][49][50].In order to demonstrate the benefits of the BerryEasy program and its integration with the PyBinding package, we consider the case of a disordered quadrupolar insulator [13].Disorder is generally difficult to implement in alternative packages, however, PyBinding's built in onsite modifiers combined with the BerryEasy package allow for an efficient route to establishing bulk topological invariants in disordered systems.We will consider a spinful version of the celebrated Benalcazar-Bernevig-Hughes model [13], the Bloch Hamiltonian is given as, This Hamiltonian supports chiral symmetry, generated as S −1 HS = −H, where S = σ 0 ⊗ τ 2 .In Refs.[51][52][53], the robustness of the bulk topology to chiral symmetry preserving disorder was studied in depth.It was shown that the bulk spin-Chern number, as defined using the method of Prodan [15] fixing ŝ = σ 0 ⊗ τ 3 in the PSO, was robust prior to the closing of the average bulk energetic gap.Following Ref. [52], this disorder will be implemented as , where w j is selected randomly from a uniform distribution, [−W/2, W/2].The bulk average density of states can be calculated rapidly for a 100×100 system using the built-in Kernel Poly- In order to determine the spin-Chern number in the presence of disorder we employ the coupling matrix approach.This is implemented by first creating a supercell of size 21 × 21 lattice cells and applying the chiral-symmetry preserving disorder.We then impose periodic boundary conditions such that the 21 × 21 system represents the fundamental unit cell and the reciprocal lattice vectors are altered as 2π → 2π/(21a), where a is the lattice constant which we set to unity.As we are working with a supercell, the reciprocal lattice vectors and the number of valence bands must be modified to account for the increase in the dimensions of the Bloch Hamiltonian.
We then perform the spin-Chern number computation.This is accomplished utilizing the W SpinLine function in BerryEasy, which allows for explicit construction of closed lines in momentum space along which the Wilson loop is computed.These closed paths are constructed by discretizing the Brillouin zone into a grid of 20 × 20 plaquettes for which the integrated Berry flux is computed.This represents a spin-resolved generalization of the coupling matrix method introduced in Ref. [36].The spin-Chern number is then averaged over 20 disorder configurations.The results shown in  The effects of disorder in solid-state systems admitting non-trivial band topology continues to be an area of active research.In general, the bulk topological invariant is robust to the introduction of disorder which leaves the bulk mobility gap intact.However, measurement of the bulk topological invariant in disordered systems is known to be a computationally demanding endeavor due to the need to simultaneously avoid finite size effects and average over many disorder configurations.
In this example, we have showcased the ability of BerryEasy package to compute the spin-resolved Chern number, a capability offered by only one other community code, and to perform the computation rapidly in a large supercell accounting for disorder utilizing the power of GPUs.In Fig. ( 5)(e), we plot the time for computing integrated Berry flux for the occupied states of eq. ( 10) through a single plaquette in an L × L supercell of the BerryEasy package when utilizing a Tesla V100 GPU vs Intel-Xeon CPU as offered in Google Colab.
Comparison with PythTB: To further demonstrate the speed advantage of BerryEasy, we plot the total time required to compute the Chern number for eq.( 4) on a supercell of linear system size L using BerryEasy and PythTB to an identical accuracy in Fig.
( 5)(f).The results demonstrate the enormous advantage of BerryEasy, particularly when run on a GPU.This speed advantage is of vital importance to achieve accurate results as finite size effects can be minimized and a greater number of disorder configurations can be considered without requiring extended computational time.The code utilized for the results in Fig. ( 5) is available online [27].

GPU environment requirements:
The GPU enabled version of BerryEasy relies on usage of the CuPy package [28].As a result, requirements for GPU environment are the same as those for use of CuPy.Namely, CuPy requires a NVIDIA CUDA GPU with CUDA Toolkit V11.2 or higher.

V. WANNIER90 TIGHT-BINDING MODELS
Finally, progress in the diagnosis of band topology has been accelerated by the ability to directly utilize the tools listed above in realistic tight-binding models derived from density functional theory using the Wannier90 software package [8].A primary advantage of the PythTB, Z2Pack and WannierTools software packages is their ability to directly interface with the output of Wannier90 for construction of a tight-binding model from which topological quantities can be computed using the built in functionalities.The BerryEasy package provides a built-in version of wanPB [54] to create PyBinding instances from the output of Wannier90, thereby extending the functionality of the tools detailed above to realistic many-band models generated from density functional theory.
As an example, we examine a known Z 2 topological insulator in two-dimensions, bilayer-bismuth [55][56][57], also known as β-bismuthene.Besides specifying the correct number of valence bands and the correct reciprocal lattice vectors, the necessary code is virtually unchanged from that displayed for eq.( 7).The results are then plotted showing that the hybrid WCC spectra is gapless in Fig. (6(b)).As such the system can be classified as a strong topological insulator.While such a computation can be performed using any of the software packages given previously (this is demonstrated by comparing the results to those obtained by the Z2Pack software package to show that they are identical in Fig. (6(c))), we remark that the ability to implement perturbations such as disorder or impurities now remains possible even in the context of complicated Wannier tight-binding models through the interface with PyBinding.
Having determined that β-bismuthene supports a non-trivial Z 2 index, it is important to determine the magnitude of the spin-Chern number as in the presence of a non-trivial Z 2 index, it has been shown that ν 0 = |C s |Mod2.As a result, for β-bismuthene, previous works have predicted C s = ±1 or C s = ±3 [56], but been unable to distinguish between the two classifications.Here we show that the spin-Chern number can be directly computed for β-bismuthene using BerryEasy.To do so it is first necessary to determine the preferred spin-direction for which the PSO is fully gapped.This is accomplished within the BerryEasy package using the spin spectrum func- tion.This function takes as parameters the kpoint at which to compute the spectra of the PSO, the PyBinding instance, the bands from which to form the PSO, and the spin operator.Through a trial and error process we have identified the proper spin-operator which produces a gapped PSO for the ground-state of β-bismuthene.This is accomplished by discretizing the Brillouin zone into a grid of 400 k-points, collecting the eigenvalues of the PSO at each point on the grid, and plotting the eigenvalues to identify the presence of a spectral gap.The results in Fig. (7(a)) clearly demonstrate that the PSO is gapped for the chosen spin operator.
Having identified the proper spin-operator, we compute the spin-resolved WCC spectra in a manner identical to that performed for eq.(9).To determine the spin-Chern number we then invoke eq. ( 6), plotting n x−,n (k 2 ) in Fig. (7(b)).The results clearly demonstrate that the system supports C s = +1, resolving the ambiguity posed in earlier studies of the system.

VI. SUMMARY
In summary, multiple excellent community codes currently exist for the computation of topological quantities in both Wannier tight-binding models and idealized toy models.However, each package is ideally suited to different use cases.The BerryEasy package fills what we view as a current gap in this spectra of offerings by providing access to state-ofthe-art diagnostics for toy models and realistic models while simultaneously offering the possibility to account for application of fields and effects by interfacing with the PyBinding package.As a result, it is our hope that this package will prove useful both for those researchers investigating new forms of topolog-ical order in clean systems as well as those interested in understanding the complex effects of disorder and application of external fields in realistic systems being actively studied in experimental settings.While BerryEasy is capable of interfacing with Kwant [10], future versions of the BerryEasy package are currently being developed to allow for all functions offered via PyBinding.Future versions are further being developed to accommodate tight-binding models defined using TBModels [6].
i s c r e t i z e ky d i r e c t i o n , 100 WCC c o m p u t a t i o n s w i l l be performed bands = [ 0 ] #C o n s i d e r i n g t h e o c c u p i e d bands ( can a l s o u s e bands=r a n g e ( 1 ) ) s y s t=c h e r n m o d e l #System we a r e c o n s i d e r i n g , eq .4 , a s PyBinding i n s t a n c e r v e c =2 * np .p i * np .d i a g ( np .o n e s ( 3 ) ) #R e c i p r o c a l l a t t i c e v e c t o r s WCCx=be .WSurf ( vec , s y s t , bands , ds , ds2 , r v e c )

FIG. 3 :
FIG. 3: (a) Density of states for 21 × 21 supercell as a function of vacancy size under periodic boundary conditions.(b)-(d) Supercell with vacancies at origin within a circle of radius r. (e)-(g) Wannier center charge spectra for supercells shown in (b)-(d) respectively.

Fig. ( 4
(a)) for the k z = 0 plane, disallowing topological classification as a strong topological insulator.At this point we compute the nested WCC spectra as a function of k z .In the BerryEasy package this performed as, import BerryEasy a s be NL= [ ] f o r kz i n np .l i n s p a c e ( 0 , 1 , 2 1 ) : vec=lambda t1 , t 2 : [ t1 , t2 , kz ] #i n t e g a t i o n a l o n g kx a s f u n c t .o f ky a t kz=0 ds =100 #d i s c r e t i z i n g i n t e g r a t i o n a l o n g kx i n t o 100 s t e p s ds2 =100 #d i s c r e t i z e i n t e g r a t i o n a l o n g ky i n t o 100 s t e p s bands = [ 0 , 1 ] #C o n s i d e r i n g t h e o c c u p i e d bands ( can a l s o u s e bands=r a n g e ( 2 ) ) s y s t=CHOTI model #System we a r e c o n s i d e r i n g , eq .9 , a s PyBinding i n s t a n c e r v e c =2 * np .p i * np .d i a g ( np .o n e s ( 3 ) ) # R e c i p r o c a l l a t t i c e v e c t o r s NL .append ( be .WNWL( vec , s y s t , bands , ds , ds2 , r v e c ) ) The above code example represents the final step in computing the nested Wilson loop, the code to define eq.(9) in the format of PyBinding, i.e. the syst parameter, is available on the corresponding webpage[27].The results seen in Fig. (4(b)), demonstrate that the nested Wilson loop shows a gapless spectra indicating that the surface resembles a non-trivial Chern insulator with gapless chiral edge states.These edge states are the famous chiral hinge modes.

FIG. 4 :
FIG. 4: (a) Gapped Wannier center charge spectra of eq.(9) at k z = 0 plane.(b) Spectra as a function of k z of Wannier Hamiltonian constructed via Wilson loop along the k x direction as a function of k y , also known as nested Wilson loop.(c) Spin-resolved Wannier center spectra detailing presence of non-trivial spin-Chern number in the k z = 0 plane of eq.(9).

FIG. 5 :
FIG. 5: (a) Average density of states at zero energy for eq.(10) upon introduction of chiral symmetry preserving disorder in a 100 × 100 system with periodic boundary conditions.(b) Spin-chern number of occupied ground-state as a function of disorder strength, averaging over 20 disorder configurations for a 21 × 21 supercell.(c)-(d) Momentum resolved Berry flux for (c) W = 1 and (d) W = 9 yielding C s = 1 and C s = 0 respectively.(e) Comparison of time for computation of Berry flux at a single k-point in momentum space utilizing a CPU or GPU as a function of the supercell dimensions.(f) Comparison of computation time in seconds to compute the Chern number of eq.(4) to identical accuracy in PythTB and BerryEasy in a supercell of linear system size L.
Fig. (5)(b), demonstrate that the spin-Chern number remains intact prior to closing of the average bulk gap.A sample output of this computation in the region before and after closing of the bulk gap is shown in Figs.(5)(c) and (5)(d) respectively.

FIG. 6 :
FIG. 6: (a) Band structure of bilayer-bismuth along high-symmetry path in the Brillouin zone.The Wannier center charge spectra is computed for the three lowest lying bands using both the (a) BerryEasy and (c) Z2Pack software package detailing identical results which classify the system as Z 2 non-trivial.

FIG. 7 :
FIG. 7: (a) Eigenvalues of the projected-spin operator when sampling the Brillouon zone along a 20 × 20 grid.The results detail the presence of a spectral gap allowing for topological classification of the negative(positive) eigenstates to determine C ↓ (C ↑ ).(b) Sum of the Wannier center charges for the negative eigenstates of the projected spin-operator at each value of k 2 .The spectral flow allows for determination of C ↓ = −1 utilizing eq.(6).