Localized Wannier function based tight-binding models for two-dimensional allotropes of bismuth

With its monoelemental composition, various crystalline forms and an inherently strong spin-orbit coupling, bismuth has been regarded as an ideal prototype material to expand our understanding of topological electronic structures. In particular, two-dimensional bismuth thin films have attracted a growing interest due to potential applications in topological transistors and spintronics. This calls for an effective physical model to give an accurate interpretation of the novel topological phenomena shown by two-dimensional bismuth. However, the conventional semi-empirical approach of adapting bulk bismuth hoppings fails to capture the topological features of two-dimensional bismuth allotropes because the electronic band topology is heavily influenced by crystalline symmetries as well as atom spacings. Here we provide a new parameterization using localized Wannier functions derived from the Bloch states in first-principles calculations. We construct new tight-binding models for three types of two-dimensional bismuth allotropes: a Bi (111) bilayer, bismuthene and a Bi(110) bilayer. We demonstrate that our tight-binding models can successfully reproduce the band structures, symmetries and topological features of these two-dimensional allotropes. We anticipate that these models can be extended to other similar two-dimensional topological structures such as antimonene and arsenene. Moreover, these models can serve as a starting point for investigating the electron/spin transport and electromagnetic response in low-dimensional topological devices.


I. INTRODUCTION
The discovery of spin-orbit coupling (SOC) induced topological phase transitions in electronic structure have led to a rapidly growing interest in the topological electronics and spintronics applications 1,2 . Due to an intrinsically strong SOC, bismuth and its compounds offer a rich playground for the development and application of topological band theories 3 . The electronic and topological structure of pure bismuth depends considerably on its crystal structure. Bulk bismuth is a topologically trivial semimetal despite a strong intrinsic SOC 4 . In contrast, when the bismuth lattice is confined to two dimensions, the resulting bismuth allotropes have been predicted to offer a rich spectrum of topologically non-trivial and distinct phases 5 . For instance, the Bi (111) bilayer-a monolayer bismuth arranged in a buckled configurationis predicted to have a quantum spin Hall (QSH) phase 6 . On the other hand, bismuthene-a graphene-like planar layer of bismuth arranged in a honeycomb latticeis reported to be a topological crystalline insulator 7,8 . Due to the diversity of non-trivial topological phases, two-dimensional bismuth allotropes have attracted extensive interest in recent years as a potential candidate for building novel topological electronic and spintronic devices [9][10][11] .
Recent experiments have demonstrated that manipulating the band topology of different 2D bismuth allotropes can lead to many exotic physical phenomena. For example, it has been shown that the ultrathin Bi (111) films can be used to tune topological edge states when interfaced with other 2D materials [12][13][14] . In a Bi (110) bilayer, elastic strains and external electric fields can significantly affect the stability of its topological phase due to its sensitivity to atomic buckling and charge doping 15,16 . For planar bismuthene, recent experiments and calculations have demonstrated a controllable orbital-filtering QSH effect due to a selective bonding with silicon carbide substrate 9 . Moreover, reports have also suggested the presence of exotic, higherorder topological hinge states in bismuth 10 . Despite the recent progress, a unified understanding of the correlation between the crystalline symmetries of 2D bismuth allotropes and their topological phases is still being developed.
Previous attempts to understand the influence of the crystal symmetries of 2D bismuth on electronic structure have largely been limited to semi-empirical tight-binding (TB) models of the bulk bismuth 4,17,18 . However, as we illustrate later, these models cannot be used for 2D allotropes of bismuth for several reasons. First, the symmetry of two-dimensional bismuth is different from that of the bulk-the semi-empirical TB model can not faithfully reflect the symmetry reduction from bulk to a surface-like structure. Moreover, the semi-empirical model does not consider the relaxation of atomic positions in the twodimensional layer relative to the bulk. A notable example is the planar honeycomb structure of bismuthene, which can not be directly related to the bulk bismuth symmetry. Consequently, these issues lead to a poor agreement between the band structure computed from the semiempirical model and first principles density functional theory (DFT) 19 .
To gain physical insights into correlation between crystal symmetry and topological phases of 2D allotropes of bismuth, here we develop effective TB models of Bi (111) and Bi (110) bilayers as well as bismuthene using localized Wannier functions constructed from first-principles calculations 20 . Wannier functions offer a natural choice of an orthonormal basis set due to the connection between the charge center of these functions and the Berry phase of Bloch states 21 . These functions can be exponentially localized to either atom centers or interstitial sites, and are therefore similar to the atomic orbitals. Furthermore, the Wannier basis set can be constructed without the need to fit free parameters to the band structures obtained from DFT or experiments. Wannier functions thereby allow us to construct a model Hamiltonian for each allotrope with relatively few parameters and yet still provide an accurate description of their band structure and its relationship with the crystal symmetry. These Wannier function based Hamiltonians can be further employed to investigate charge and spin carrier transport in electronic devices based on low-dimensional topological materials.
This article is organized as follows. In Sec. II we discuss the limitations in the applicability of the semiempirical TB parameters from bulk bismuth to a Bi (111) bilayer. We then introduce a method for parameterising new TB models from the first principles calculations in Sec. III. First we outline our density-functional methods in Sec. III A, then in Sec. III B we explain the construction of the TB models from localized Wannier functions, including constraints due to the symmetries of the crystallographic lattices. Finally, in Sec. IV we discuss the results of our parameterizations for a Bi (111) bilayer (Sec. IV A), bismuthene (Sec. IV B), and a Bi (110) bilayer (Sec. IV C). We conclude in Sec. V.

II. THE APPLICABILITY OF TIGHT BINDING PARAMETERS FROM BULK TO A BILAYER
The electronic properties of Bi (111) bilayers 6,19 are typically modeled using semi-empirical TB parameters derived for bulk bismuth 4 . For example, this strategy was adopted in Ref. 6 and 19, by using the hopping parameters of bulk bismuth to bismuth bilayer, with the SOC strength (λ) increased from its bulk value to better fit the energy splittings between the bands of a bilayer when compared to DFT calculations. In the following we will investigate the validity of this approach, as outlined in Ref. 6 and 19, and determine if using bulk tight-binding parameters to two dimensional structure can faithfully reproduce all critical electronic features of the topologically nontrivial Bi (111) bilayer.
We use TB parameters for bulk bismuth to calculate the band structure of a Bi (111) bilayer but trun- cate these parameters by retaining only the hoppings inside the bilayer 6 . All calculations were performed using NanoNet, an extendable Python framework for electronic structure computations based on the tight-binding method 22 . The resulting band structure is shown in the colored lines of Fig. 1. The SOC strength is taken from the bulk value (1.5 eV). In this figure we have also plotted a band structure for the bilayer computed from DFT (black lines) for reference (see Sec. III A for method). There is good agreement between the valence bands of TB and DFT for a SOC strength of 1.5 eV but the agreement between the conduction bands is poor. The lowest conduction band of the TB model has energy minima at points Γ and M whereas DFT predicts two additional minima, one at K and another along the high symmetry segment M → Γ.
Another property of interest is band inversion, i.e. the exchange of a crystal's electronic properties (e.g. parity) between its conduction and valence bands. The parity of each eigenstate at the Γ point in Fig. 1 is labeled as either even (+) or odd (−). The parities of the bands shown in the figure are even for the valence bands and odd for the conduction bands with no SOC (λ = 0.0 eV). When we use the bulk SOC value of 1.5 eV, the parity of the second highest valence band and lowest conduction band are exchanged. This is in good agreement with a previous study that assumes a Bi (111) film retains the bulk lattice parameter at about 20-50 nm of film thickness 23 . However, several reports have shown that when bismuth films are reduced to only one bilayer thickness, the change of parities occurs between the highest valence band and lowest conduction band, contrary to the data presented in Fig. 1. 24,25 This contrasting behavior has been associated with the change of the lattice parameter from its bulk value. 24,25 SOC also significantly affects the band gap in a Bi (111) bilayer. Fig. 2 shows the variation of the band gap at Γ point between the highest valence band and the lowest conduction band as a function of the SOC strength. The maximum energy gap that can be achieved at the Γ point by increasing the strength of the SOC is approximately 0.4 eV. This gap is smaller than the gap computed from DFT (see inset of Fig. 2) which typically underestimates band gaps compared to their experimental values 26 . We have found that the band splitting agrees with DFT results if we increase the SOC strength in the TB model to 1.8 eV, as suggested by Ref. 19. However, the band shape and the exchange of parities still do not agree with ab initio calculations and experimental observations 24,27 . Based on the results above, we can conclude that an effective tight-binding model of two-dimensional bismuth layered materials can not be achieved by using bulk hopping parameters, even by tuning the strength of SOC. Consequently, a new model is needed for 2D allotropes of bismuth that adapts to the structural changes of the material when reduced from the bulk to a thin film.

III. TIGHT BINDING MODELS FROM WANNIER FUNCTIONS
A. Density-functional method We first obtained the electronic structure of the 2D allotropes of bismuth through DFT as implemented in the Vienna Ab initio Simulation Package (VASP) 28,29 . Exchange and correlation effects were captured by the generalized gradient approximation and the Perdew-Burke-Ernzerhoff functional 30 . The atomic structures of the allotropes were first optimized prior to a self-consistent convergence of the electronic structure using a Γ-centered 21 × 21 × 1 k-point grid. The energy cutoff for the planewave basis used was 600 eV. Wavefunctions obtained using DFT were utilized to construct TB models based on Wannier functions using the code Wannier90 31 , as discussed in the next section.

B. Construction of Wannier tight binding models and symmetry constraints
For periodic systems with translational symmetry, the one-particle states can be expressed by the Bloch state |ψ nk with band index n and crystal momentum k. When deriving parameters for the hopping of an electron from one orbital in a crystal to another it is more convenient to consider localized orbitals rather than Bloch states because the latter are delocalised. Wannier functions are one such choice of localised orbitals and they can be constructed by the inverse Fourier transform of a Bloch state. In this way we can obtain a real-space Wannier TB Hamiltonian by a discrete Fourier transformation 32 : where R is a real-space lattice vector, 0 is the realspace lattice vector that defines the home cell [i.e. R = (0, 0, 0)], N is the number of points in the k-point grid and E q are eigenvalues from DFT. U q is a unitary transformation that takes the Bloch state at point q in k-space to a rotated Bloch state in the Wannier gauge 32 , i.e.
where n and m are band indices. The corresponding reciprocal space TB Hamiltonian H αβ (k) can then be expressed as a Fourier transformation from the real space Hamiltonian H αβ (R) to k space: where the subscript αβ denotes a Hamiltonian matrix element that corresponds to hopping from orbital α at τ α in the home cell to orbital β at τ β within a cell at R. H αβ (R) can be further expanded as: where t αβ (R − 0) represents the hopping parameter between neighboring atomic orbitals α in the home cell and β in the cell R extracted from the Wannier TB Hamiltonian. For simplicity, we have used simple numerical subscript t n in later presentation of our TB models. A detailed discussion of the extraction process and the role of each hopping parameter can be found in the Appendix. Next we discuss the symmetry constraints on this Hamiltonian. The representations generated from the method implemented in Wannier90 33 often contain small numerical errors that can break the symmetries of the crystal's energy bands. We correct these errors and restore the corresponding constraints on the crystal's symmetry. Our Hamiltonians satisfy the following symmetry constraints in k-space and real-space respectively 33 : D(g) is the matrix representation of the symmetry operation g, which is an element of symmetry group G. The matrix representations used for the 2D allotropes of bismuth (D 3d for the Bi (111) bilayer, D 6h for bismuthene and D 2h for the Bi (110) bilayer) are given in the Appendix. The subscripts α and β are the orbital indices before the symmetry operation and i and j are the corresponding indices after this operation. R is a lattice vector that defines the cell's position after the symmetry operation.
To ensure our Hamiltonian matrix satisfies each symmetry constraint for a particular allotrope, we take a group average over all Hamiltonians transformed by the symmetry operations in G through Eq. 6. This yields a symmetrized Hamiltonian 33 : where l is the number of elements in the symmetry group.

IV. RESULTS
To constrain our basis functions to crystal symmetries and fix Wannier orbitals to atomic positions, we construct our TB Hamiltonian from a spinless case without performing maximal localization on Wannier Functions. We have also treated SOC independently, evaluating it in the basis of atomic orbitals and then fitting the SOC parameters by comparing to the bands of DFT. In the following subsections we will discuss the symmetry of the crystal structures, properties of the basis functions, and construction of the TB models for each allotrope.
It should also be noted that we will emphasize the necessity of determining the symmetry properties of electronic bands in constructing the TB model for each allotrope. The SOC effects are also treated differently for each bismuth allotrope considering the changes in crystalline symmetries. We calculate the matrix representations of all the symmetry operations for eigenstates at high symmetry points. Finally we confirm the accuracy of our TB models by calculating the basis functions for these irreducible representations using the projection operator method 34 and check against orbital characters obtained using DFT calculations for consistency.

A. Bi (111) bilayer
The Bi (111) bilayer has a quasi-2D honeycomb structure characterized by an out-of-plane buckling with an intra-layer spacing of approximately 0.87Å. There are two inequivalent atoms in the primitive cell as shown in Fig. 3. The crystal structure of the Bi (111) bilayer has a symmorphic space group symmetry of P 3m1 (D 3 3d , SG164). The symmetry generators for its point subgroup include the identity (E), inversion (P ), threefold rotation along the z direction (C 3z ), and three twofold rotations (C 2 ). From the orbital characters obtained from our DFT calculations, we find the s orbitals form an isolated set of bands well below the Fermi level. Closer to the Fermi level, p orbitals form two sets of isolated bands: three valence bands (VBs) and three conduction bands (CBs). We therefore choose to include only p orbitals in our TB model.
The number of nearest neighbor (NN) hoppings to include is a critical parameter for a TB model. To investigate the relationship between neighbor interactions and the agreement of our model with DFT, we calculate the average energy difference between the two models along a path of high symmetry in the first Brillouin zone (FBZ). We find that a Hamiltonian with second NN interactions can well reproduce the topological properties of the Bi (111) bilayer, with good agreement to the occupied states. However, the inclusion of up to fourth NNs is needed to recover the features of unoccupied bands (see Fig. 4 (b)), but at the cost of requiring more hopping parameters. In the following we will focus more on the details of building a TB model based on second NN interactions. The model considering fourth NNs can be simply constructed in a similar way by incorporating more neighbor interactions.
Next we discuss the steps involved in constructing our TB model in more detail. We start from a spinless Hamiltonian written in a basis of p-like orbitals. The basis set is {p (1) and (2) denote each of the two atoms in the primitive cell, as labelled in Fig. 3 (c). We use a convention where the real space Hamiltonian is dependent on the lattice vectors R and sublattice vectors τ α and τ β that correspond to particular atomic sites. We construct hopping parameters within a group of NN interactions by applying symmetry operations to the Hamiltonian (which in turn depends on R, τ α , and τ β ). It is thereby possible to generate all hopping parameters within a given set of NN interactions [i.e. first, second, third, or fourth, as shown in Fig. 3 (a)] from the interactions between one pair of neigbor atoms in the set. We denote r 1r 3 as nearest neighbor hopping vectors, r 4r 6 as second-nearest neighbor hopping vectors, r 7r 9 as third-nearest neighbor hopping vectors, and r 10r 15 as fourth-nearest neighbor hopping vectors The symmetry operations not only generate hopping parameters but also enforce the symmetry constraints that are defined in Eqs. 5 and 6.
The k-space TB Hamiltonian can then be constructed from the real-space Hamiltonian using Equation 3. The k-space tight-binding Hamiltonian consists of sub-blocks related to hoppings between two sub-lattice atoms as shown in Equation 8: Here we use the superscript to denote interactions be-tween atoms and subscript for interactions between orbitals in Hamiltonian H. The diagonal blocks H 11 (k) and H 22 (k) are on-site energies of the two bismuth atoms (Bi1 and Bi2) in the home cell, while H 12 (k) and H 21 (k) are hopping matrix between the two bismuth atoms. The matrix elements in every H ij (k) are therefore denoted as H ij (k), which represents the interactions between orbitals. The indices 1, 2, ..., 6 correspond to the obtials {p zBi1 , p xBi1 , p yBi1 , p zBi2 , p xBi2 , p yBi2 }. We have dropped k dependence in the sub-blocks for brevity: H 22 (k) and H 21 (k) can be obtained by applying an inversion operation to H 11 (k) and Hermit adjoint to H 12 (k) respectively: Constrained by the inversion symmetry and Hermitian condition, H 12 (k) is symmetric. Hence, we only need to solve for H 11 (k) and H 12 (k). H ii (k) represents the selfinteraction energy of an orbital with itself, which can be expressed as: . Here t i denote real space hopping parameters as defined in Eq.

4.
For the remaining off-diagonal elements H 11 (k), they are given by: where the constants h (m) ij are related to hopping For off-diagonal hopping sub-block H 12 (k), the matrix elements can be similarly expressed as: Again, h . Following this procedure, we obtain the spinless TB Hamiltonian with 12 independent parameters (t 1 -t 15 , except t 10 , t 13 and t 14 ) for Bi (111) bilayer with up to the next-nearest neighbor interactions as summarized in Table I. If one desires a better description of the conduction bands, third and fourth NN interactions can be included. As shown in Table I, a TB model involving up to four NN interactions included needs 10 additional independent hopping parameters (t 16 -t 27 , except t 18 and t 20 ). In the Appendix, we discuss the symmetry constraints imposed on the real space Hamiltonian and give the explicit forms of the real space Hamiltonian up to fourth NN intaractions.
Next, we addSOC interactions to the spinless Hamiltonian. The SOC contribution to the spinor Hamiltonian is expressed as: where λ is the on-site SOC strength. This can be expanded further as: where L is the angular momentum operator on orbitals and S is the angular momentum operator on spin states, L ± and S ± are the raising operators defined as L x ± iL y and S x ± iS y , respectively. The on-site SOC strength λ has a large effect on inducing band inversion at Γ point. When λ = 1.1 eV, the system experiences a gap-closing phase transition, with exchange of parity character between occupied and unoccupied states. The VB and CB near the Fermi level will be separated further apart if λ is further increased. We notice a particular feature in DFT bands, which shows a small segment of nonzero curvature in the highest VB near Γ just below the Fermi level ( Fig. 4 (b)). We can obtain this feature with large λ at the expense of band mismatches in other parts of the occupied bands. This issue can be resolved when we introduce a next-nearest neighbor SOC term. This additional next-nearest neigbor SOC term reflects the strong SOC strength of bismuth and the threefold rotation symmetry of the honeycomb lattice. The form of this term is expressed as 35 : where H I is the intrinsic SOC originated from the buckling of of the (111) bilayer, λ mn is effective SOC strength associated with hopping from orbital m in home cell to orbital n in next-nearest neigbor cell, and [σ z ] αα is the element (α, α) of the Pauli matrix. As discussed by Kochan et al. 35 , only diagonal elements in the spin basis are nonzero for intrinsic SOC. Combining symmetry considerations and the requirement of Hermiticity, we obtain the intrinsic next-nearest SOC term λ mn in the form of a 3 × 3 matrix with 6 independent parameters (more detailed derivation of this term is given in the Appendix): We obtain the parameters from λ 1 to λ 6 from first principles again and the on-site SOC term λ is tuned to fit the energy bands at Γ near the Fermi level. Therefore the SOC part adds another 7 parameters to our TB model as shown in Table I.

IRs
Basis This concludes the construction process of our spinless and spinor TB models. An advantage of these models is that they use minimal parameters to accurately reproduce the electronic structure of the bismuth allotropes studied here. We now examine the symmetry properties and basis functions of our model to verify this claim. Fig. 4 shows the irreducible representations (IRs) at Γ obtained from our TB models and corresponding orbital characters from DFT calculations shown in circles. Table II summarizes the basis functions of IRs. The notation we use for IRs is consistent with Koster's 36 . The upper part of the table is the basis functions for the spinless case while the bottom part corresponds to the case with spin. For degenerate states, mixing between basis functions within the degenerate subspace is allowed as long as the two states are related by time reversal (TR) symmetry.
In Fig. 4 the highest VB (with IR Γ + 3 and basis functions of {|p + x , p + y }) has a double degeneracy protected by the threefold rotation symmetry of the crystal structure and the lowest CB belongs to the Γ − 2 IR, which consists of an equal mixing of p z orbitals from the two atoms in the primitive cell. These orbital characters are consistent with those from DFT calculations. Therefore, the eigenstates in our TB models have the correct symmetry properties. After taking electron spins into account, the new double group IRs near the Fermi level are: Γ + 3 ⊗ Γ + 6 = Γ + 4 ⊕ Γ + 5 ⊕ Γ + 6 and Γ − 2 ⊗ Γ + 6 = Γ − 6 . In other words, the product of the basis functions for the highest VB (corresponding to Γ + 3 ) with spinor basis functions can be decomposed into direct sum of states Γ + 4 , Γ + 5 and Γ + 6 , and a product of the lowest CB Γ − 2 state with a spinor yields an anti-symmetric Γ − 6 state. As a result of both TR and inversion symmetry, all the eigenstates are doubly degenerate. In addition, one dimensional IRs Γ ± 4 and Γ ± 5 are related by TR symmetry and form a Kramer pair.
The band inversion is characterized by the exchange of states between Γ + 4 ⊕ Γ + 5 and Γ − 6 . The accompanying exchanging of parity between occupied and unoccupied states gives rise to the topological nontriviality of the system. The Z 2 topological invariant can be calculated simply by taking parity values at four TR invariant momentum points in the FBZ 37 . Inversion parity eigenvalues of occupied states at Time Reversal Invariant Momenta (TRIM) are listed in Table III. The Z 2 invariant of 1 calculated from our TB for Bi (111) bilayer is consistent with the values reported in literature using DFT methods 24 .
Bismuthene is a 2D allotrope of bismuth arranged in a planar honeycomb crystal structure. The crystal structure has symmetry of space group P 6/mmm (D 1 6h , SG191). The symmetry generators for its point group include identity (E), inversion (P ), twofold rotations along the z direction (C 2z ), and two threefold rotations (C 3 ). By analyzing the orbital character obtained from our DFT calculations (see Fig. 5) we find that there is negligible hybridization between s and p orbitals, and therefore we choose to include only p-orbitals in our TB model of bismuthene. We have found that the third and fourth NN hoppings have negligible contribution to the band structure and the TB model with next-nearest neighbors already shows satisfactory agreement to DFT results. Therefore, we only include up to next-nearest hoppings in the TB model. First we start with Hamiltonian without spin. The basis set we take is {p zBi1 , p xBi1 , p yBi1 , p zBi2 , p xBi2 , p yBi2 } where the 1 and 2 subscripts distinguish the sublattice atoms.
Similar to the Bi (111) bilayer, the k-space Hamiltonian of bismuthene also consists of four sub-blocks: The form of the sub-blocks are: Again, we only need to solve for H 11 (k) and H 12 (k) since the other two sub-blocks can be obtained by inversion and Hermitan adjoint. The H ij (k) of bismuthene has the same components H ij (k) as the Bi (111) bilayer except those zero elements in H 11 (k) and H 22 (k), because bismuthene has all the symmetry element generators of Bi (111) bilayer with an extra horizontal mirror symmetry σ h element. The general form of the on-site hopping term H ii (k) is: respectively. The only off-diagonal term H 23 (k) in H 11 (k) has the form: where the constants h 23 = t 8 , and h Again, h . Finally, we obtain 9 independent hopping parameters as shown in Table IV.
We then include the effect of SOC in the TB Hamiltonian by duplicating the basis set to introduce a spin degree of freedom and adding an on-site SOC term H soc . The SOC strength for bismuthene (λ) is taken directly from our ab initio calculation as defined in Eq. 23. This approach is different from that we used for a Bi (111) bilayer where we have added the next-nearest SOC terms to capture the non-negligible buckling effect on the SOC. In comparison, we found that for bismuthene the onsite atomic SOC term alone can well reproduce the band structure. The λ value for bismuthene is presented in Table IV.
We check the symmetry properties of bismuthene in our TB model by identifying the IRs of the eigenstates at the Γ point and projecting them onto their corresponding basis functions, listed in Table V. We compare these IRs with projected orbital characters predicted by DFT as shown in Fig. 5. Without SOC, the eigenstates formed by spinless p-orbitals can be categorized into representations Γ − 2 , Γ + 3 , Γ + 6 and Γ − 5 . Γ + 3 and Γ + 6 states correspond to symmetric states formed by out-of-plane p-orbitals and in-plane p-orbitals respectively. Γ − 2 and Γ − 5 correspond to anti-symmetric states. In the Fig. 5, it can be seen that the orbital characters from TB model are consistent with those from DFT calculations.
To analyze the symmetry of the Bloch states with SOC, we introduce spinor representation Γ + 7 and generate IRs of spinful eigenstates from previous IRs of spinless states by again taking the direct product with spinor and decomposing onto IRs of the double group for occupied states: SOC has several effects on the electronic states of bismuthene. First we notice the degeneracy corresponding to the Γ + 6 state and Γ − 5 state is lifted. Band exchange happens between the upper two bands in the VB and another two high-energy bands in the CB. Even though SOC exchanges the parity of these two VB states, it does not affect the topological property of bismuthene because there is no exchange between unoccupied and occupied states. We also see that with SOC the electronic structure of bismuthene undergoes a transition from semimetallic phase to an insulator phase after a band gap opens near the K point.
We confirm this further by calculating the Z 2 invariant using parity values at TRIM listed in Table VI. We have numbered bands from low energy to high energy and neglected spin degeneracy for brevity. The result that bismuthene has Z 2 = 0, is consistent with the previous studies. 38 TABLE V. Basis functions of IRs at Γ point (point group D 6h ), obtained from projection operator method. 34 . There are two sets of p-orbital combinations allowed for Γ 8+ and Γ 7− IRs. The bands with these two IRs consist of linear combination of the in-plane and out-of-plane p-orbitals, in a similar way to Γ 6+ and Γ 6− IRs from D 3d group as listed in

C. Bi (110) bilayer
A Bi (110) bilayer consists of two layers of planar bismuth. Upon relaxation, we find this 2D allotrope has a non-symmorphic crystal symmetry unlike the (110) surface of bulk bismuth, which has been reported experimentally 39 . The absence of buckling in each planar layer leads to a crystal symmetry of space group P mna (D 7 2h , SG53). The symmetry generators for this group include two twofold rotations C 2y , C 2z and two corresponding mirror operations σ 2y , σ 2z with a glide translation vector of ( 1 2 , 1 2 , 0). The presence of glide plane symmetry usually results in additional band degeneracies in the electronic structure. The projected orbital characters from our DFT calculations show that electronic bands of Bi (110) bilayer are mainly from p-orbitals (see Fig. 6). The contribution of s-orbitals is to energy bands far below the Fermi level and in this case isolated from the rest of the band structure. Therefore we only consider p-orbitals as basis functions in our TB model for Bi (110) bilayer.
It is more convenient to use neighboring cells instead of neighboring atoms in the construction of TB models for Bi (110) bilayer since it has a low-symmetry crystal structure and an accurate description of its electronic structure requires inclusion of more neighboring atoms compared to the other two allotropes. In this way, we have used interactions between orbitals in the home cell and orbitals in all eight neighbor cells around it as shown in Fig. 3 (f).
As before, we start by constructing the spinless Hamiltonian. The real-space Hamiltonian matrix describes hopping between the home cell and neighbor cells at R. The entries of the Hamiltonian can be divided into 16 sub-blocks to represent hopping between four inequivalent atoms in the primitive cells. Therefore, we analyze the form of the real-space Hamiltonian in the unit of a sub-block. For    We now examine the symmetry properties of electronic bands near the Fermi level. For spinless eigenstates, the electronic bands at Γ consist of mixing between the sets of p y and p z orbitals, while mixing of p x orbitals is disallowed, as shown in Table VIII. After SOC is included, the orbital characters at Γ are a mixture of all p-orbitals. The effect of SOC also induces a band inversion at Γ between the two bands at either side of the Fermi level, implying a gap-closing topological phase transition.
The most notable feature in the electronic bands of the Bi (110) bilayer is the additional degenerate states at X and Y due to presence of glide mirror symmetries. For a spinless Hamiltonian, the bands at X and Y are doubly existence of symmetry-protected band degeneracies. Taking spins into consideration, the orbital characters at X and Y are from mixing of the p-orbitals across all atoms. The 2-fold band degeneracy imposed by the glide mirror symmetries is still retained. Then combining with TR and inversion symmetry, the bands will become fourfold degenerate with SOC included. The representations of these fourfold degenerate bands are related to each other by TR and σ z symmetry. Finally, we calculate inversion eigenvalues of the occupied states and list them in Table X. Z 2 =1 indicates Bi (110) bilayer is a topological insulator with non-symmorphic symmetry 16 . In conclusion, we have constructed tight binding models for three different allotropes of 2D bismuth, namely, the Bi (111) bilayer, bismuthene and the Bi (110) bilayer, by projecting Bloch states calculated from DFT onto localized atomic orbitals as Wannier functions. We have utilized a minimal set of independent hopping parameters in the TB model to represent the energy bands near the Fermi level. These TB models can accurately reproduce the band topology of these three allotropes compared to previous semi-empirical TB models.
In all cases considered here, we find that the crystalline symmetry plays an important role in constructing Wannier-based TB models. We have tailored the form of the TB model to meet the symmetry constraints imposed on each type of bismuth allotrope. These symmetry constraints can eliminate the numerical errors in the Wannier Hamiltonian, leading to a significant reduction in the number of independent TB parameters. Symmetrized TB models can faithfully recover the band representations near the Fermi level for each bismuth allotrope and also help us to understand the physical origin of specific features in the band structures.
In summary, we have shown that TB models of twodimensional bismuth allotropes can be effectively derived from the Wannier basis. This symmetry-based approach could be conveniently applied to other topologically nontrivial two-dimensional materials such as antimonene and arsenene 40,41 . Another advantage of these simplified but accurate models is that they can be easily modified to incorporate external factors such as strain, electric and magnetic fields 42,43 . Moreover, these models can be a cornerstone for designing new topological structures exhibiting exotic quantum transport properties. For example, the tight-binding models proposed by Su, Schrieffer and Heeger 44 have inspired experimental realization of stable topological quantum states in graphene nanoribbons 45 . These models could also be effective in investigating structure complexities in topological systems, e.g. interactions with defects and disorders 46 . Therefore we expect our simplified and symmetrized tight-binding models could play a critical role in multiscale modeling of low-dimensional topological structures, potentially offering an accurate picture of device physics at lowered computational cost 22 . Here we present how to generate a real space Hamiltonian from the interaction effect of symmetries operations on atomic orbitals and lattice sites. The crystal structure of a Bi (111) bilayer has the D 3d point group symmetry. The symmetry constraints on the real space basis function include Identity (E), inversion (I), threefold rotation (C 3 ) and mirror symmetry (σ). We can write down the matrix representation of these symmetry operations on a set of p orbital basis {p z , p x , p y } from the same atom: Other symmetry elements can be generated using: We also introduce matrices to describe the symmetry operations on lattice sites. In the primitive cell of Bi (111) bilayer, there are two atoms. We can therefore write these operations as 2×2 matrices in the basis of lattice sites 1,2 as: We can readily obtain the form of real space spinless Hamiltonian by Equation 5. The D(g) matrices in the Equation 5 are 6×6 matrices as a result of the tensor product between orbital representation matrices and lattice site representation matrices. In addition, the real space Hamiltonian should be Hermitian: where R is a lattice vector. This can be derived from a Fourier transform of the Bloch Hamiltonian.
In the following, we give the form of different real space Hamiltonians of a Bi (111) bilayer by considering up to 4 th atomic neighbors in detail. Each real space Hamiltonian consists of 4 3×3 sub-blocks corresponding to four possible ways of hoppings between two lattice sites. Therefore we can analyze different atomic neighbor hoppings by focusing on specific sub-blocks of a real space Hamiltonian. The hopping terms in all real space Hamiltonians are significantly constrained by symmetry operations. These terms are given in Table I.
First, for investigating the on-site and the nearest atomic NNs, we derive the real space Hamiltonian for the home cell [000]: The diagonal sub-blocks are on-site hoppings while the off-diagonal sub-blocks are hopping between lattice sites in neighboring cells and the home cell. In the diagonal sub-blocks, C 3 , σ, and Hermiticity lead to non-zero hopping terms in the diagonal entries of H [000] . As the two lattice sites are related by inversion symmetry, the two diagonal sub-blocks should have the same form, resulting in the only independent parameters in the diagonal sub-blocks to being t 1 and t 2 . Similarly, for the offdiagonal sub-blocks, the symmetry constraints also limits the number of independent parameters to four. Consequently, H [000] can be written with only six independent hopping terms.
We now analyze the form of next-nearest atomic neighbor hopping terms by examining the interactions between the home cell R = [000] and the cell at R = [100]. There are also four types of hopping terms here: hopping from lattice site 1 (2) in home cell to lattice site 1 (2) in the cell at R = [100]. The real space Hamiltonian of cell R = [100] is expressed in the form: By applying mirror symmetry σ (100) and Hermiticity, further constraints can be found on H [100] : Hence there are only six independent parameters in the next-nearest neighbor part of H [100] . The off-diagonal blocks are zero since we only include the terms for nextnearest atomic neighbors here.
As for 3 rd atomic neighbor interactions, we use the off-diagonal part of the real space Hamiltonian of cell R = [110], which has the form: There also exist restrictions on entries in the sub-blocks due to C 010 2 and Hermiticity. These are: Consequently there are only four independent parameters for 3 rd atomic neighbors.
Finally , for 4 th atomic neighbor interactions, we look at the off-diagnoal terms of H 4 100 : There are six independent entries in the sub-block of the real space Hamiltonian H 4 100 due to inversion symmetry and Hermiticity. The rest of the hopping terms can be generated easily using the other C 2 and C ± 3 rotations and taking the Hermitian adjoint. Now we have finished construction of spinless Hamiltonian with up to 4 th neighbor interaction included using the form of real space Hamiltonian.
2 nd Nearest neighbors:  In the Bi (111) bilayer we consider two types of SOC terms, namely on-site SOC and a next-nearest atomic neighbor SOC referred to as intrinsic. The on-site SOC can be directly calculated using Eq. 23 and has the form: Next we turn to the intrinsic SOC. The form of this SOC can be obtained by using the approach in the previous section for deriving hopping terms. The crystal symmetries, TR symmetry and Hermicity play important roles in limiting the number of independent SOC terms. We where σ z is a Pauli matrix and H I, 11 [110] represents the SOC between atom 1 in the home cell and atom 1 in cell R = [110]. H I, 11 [110] has the form: Using the findings in the previous section, when the symmetry constraints act on the next-nearest hopping terms, there will only be six independent terms. We can therefore include the contribution of SOC in a Bi (111) bilayer through an on-site SOC strength λ and next-nearest SOC strengths λ 1 to λ 6 .

c. Bismuthene
In the case of bismuthene, we construct a real space Hamiltonian by including hopping terms to the nextnearest atomic neighbor. Similarly to a Bi (111) bilayer, we give the form of Hamiltonian for cell R = [000] as: The next-nearest neighbor hopping terms can be obtained by looking at diagonal sub-blocks of H [100] : Applying the same symmetry operations of a (111) bilayer listed in the previous section on the real space Hamiltonian of Bismuthene, the real space Hamiltonian including up to NN interactions can be obtained. We describe the electronic structure of bismuthene with nine independent parameters as shown in Table IV. The number of hopping terms are reduced by the presence of crystal symmetries and Hermicity. When considering SOC, the form of on-site SOC Hamiltonian of bismuthene is identical to that of Bi (111) bilayer.

d. Bi (110) bilayer
For a Bi (110) bilayaer, we use a neighboring cell instead of the neigboring atoms, but the explicit form of the real space Hamiltonian of cell R = [uvw] still follows the definition as described for the Bi (111) bilayer and bismuthene. For cell R = [000], we write the sub-blocks representing the interactions between one lattice site in the home cell R = [000] and another lattice site in the home cell as, with superscript indicating the sub-lattice sites and subscript the cell that interacts with the home cell:

Character table for symmetry groups of three allotropes
The character tables for a Bi (111) bilayer and bismuthene can be derived from their corresponding point groups. Those character tables can be readily found in Koster's book 36 . Here we show the character table for a Bi (111) bilayer (D 3d ) and bismuthene (D 6h ) in Table  XI and Table XIII, respectively. The character table for a Bi (110) bilayer involves nonsymmporhic symmetries and the characters for Γ, M , X and Y point are listed in Table XIV to Table XVII.        1