Origin of band inversion in topological Bi₂Se₃ Open Access

We appreciate the fact that topological materials form a broad class of materials ranging from insulators,¹ superconductors,² and semi-metals.³ To date, the topological idea in materials is still scanty and perhaps the best understood topological materials are topological insulators,⁴ which form the basis of our work herein. Topological insulators in bulk have insulating properties but often exhibit conducting edge states on the boundary protected by time reversal symmetry.⁵ An important aspect of a topological insulator is band inversion, which results from the crossing of the valence band and conduction band of different parities and a gap opening due to spin–orbit coupling.⁶ 

As we focus on novel topological insulators, it is important to shed light on their intrinsic character. These insulators are characterized by a Hamiltonian that is not, in any way, adiabatically connected to the atomic limit. This implies that when tuning external parameters to change the Hamiltonian, the process is extremely slow such that the material remains in the ground state throughout.⁷ A typical example of the atomic limit may be a case of diamond. Imagine that we pull the carbon atoms extremely far apart, yielding isolated carbon atoms.⁸ If we use diamond, it is possible to do this process without closing the gap implying that diamond is adiabatically connected to the atomic limit. This outcome is expected since diamond is a normal insulator.⁹ 

The same process can be applied to Bi₂Se₃, which in crystalline form is a topological insulator with a gap of 0.27 eV–0.30 eV.^10–12 If the previous procedure is repeated such that we pull the system apart to reach the atomic limit, it is impossible to do it without closing the bulk gap. This implies that Bi₂Se₃ is not adiabatically connected to the atomic limit and thus suits to be called a topological insulator.¹¹ 

The big question here is as follows: “What is the origin of such anomaly in topological insulators?” It is important to note that the electron wave function twists as you cross the Brillouin zone in a topological insulator.¹³ This is, however, not true in an atomic limit, and so the twist needs to be undone when getting to such a limit. This is what is achieved by gap closure. Equally important is the fact that these twists can be labeled by topological invariants¹⁴ whose mathematical form is highly dependent on the type of topology but are related to the Berry phase-like quantities¹⁵ that determine the evolution of the wave function as one crosses the Brillouin zone. In topological three-dimensional insulators, the topological invariant is usually a set of four numbers¹⁶ that take one of the two values (⁠$Z$ classification) and can be calculated by the Wannier charge center evolution across the Brillouin zone, as described in Refs. 5 and 17–19. Alternatively, if a system possesses inversion symmetry, the topological invariant is obtained by calculating the parity of eigenstates at special points in the Brillouin zone, as explained in Ref. 18.

Apart from the topological insulators, there is also another interesting class of insulators known as Chern insulators.²⁰ These are two-dimensional materials with topological and magnetic orders, as explained by Garrity and Vanderbilt.²⁰ In this class of materials, it is worth noting that the topological invariant is the Chern number attained by an integration of the Berry curvature over the Brillouin zone.

First-principles calculations have been employed in calculating topological invariants^21–23 by involving Berry phase-like quantities. Many studies have employed Wannier functions and packages that implement such calculations and have interfaces to some density functional theory (DFT) packages such as Wannier Tools.²⁴ There are also numerous databases such as the Topological Materials Database²⁵ that can aid in the topological classification of many materials based on semi-local DFT. If one accesses such a database, it forms a perfect starting point to figure out the topological order of a material. Caution has to be employed if higher levels of theory are applied such as hybrid functionals and GW. This application may easily lead to different results from those in the existing databases subsequently necessitating further analysis.

Having known that band inversion is a key ingredient of any topologically non-trivial material,²⁶ we will, therefore, focus on the fundamental definition of band inversion, how to recognize it in a band structure, and making inferences if band inversion is observed in a material. Our work will mainly focus on the topological insulator, Bi₂Se_3, and we show that indeed spin–orbit coupling is not a causative agent of band inversion.

This paper is organized as follows: In Sec. II, we account for the technicalities used in our calculations that may be critical for future reproducibility. In Sec. III, we present finer details of Bi₂Se₃ with and without spin–orbit coupling (SOC). We present the structural properties of Bi₂Se₃ in comparison to other works (Sec. III A), the various orbital contributions in Bi₂Se₃ and their specific locations in the band-structure (Sec. III C), the band-structure of Bi₂Se₃ (Sec. III B), the demonstration of band inversion in Bi₂Se₃ with and without spin–orbit (Sec. III D), and the calculation of topological invariants (Sec. III E). Conclusion and future perspectives of this article are highlighted in Sec. IV.

We carried out scalar relativistic and fully relativistic first-principles calculations employing the density functional theory^27,28 formalism as implemented in the siesta method.²⁹ We used a single-ζ basis set for the semi-core states and double-ζ plus polarization for the valence states of all the two atoms. Exchange and correlation functionals were treated within the generalized gradient approximation.³⁰ Core electrons were replaced by ab initio norm-conserving fully separable,³¹ Troullier–Martin pseudopotentials.³² In siesta, the one-electron eigenstates are expanded in a set of numerical atomic orbitals. A Fermi–Dirac distribution with a temperature of 0.075 eV was used to smear the occupancy of the one-particle electronic eigenstates. To get a converged system, a two-step procedure was performed: First, to relax the atomic structure and the one-particle density matrix with a sensible number of k-points (9 × 9 × 5 Monkhorst–Pack³³ k-point mesh); second, freezing the relaxed structure and density matrix, a non-self consistent band structure calculation was performed with a much denser sampling of 60 × 60 × 60 in the real space integrations, and a uniform grid with an equivalent plane-wave cutoff of 600 Ry was used. In this system, we equally computed the fat bands (⁠$Fi,n,σ,k→$⁠), which are defined as the periodic equivalent of the Mulliken population,

where $Ci,n,σ,k→$ and $Si,j,k→$ are the orbital coefficients and the overlap matrix elements in that order, respectively. The indices i and j denote basis functions, n is the band index, σ is the spin index, and $k→$ is a reciprocal vector in the Brillouin zone.

Since Bi₂Se₃ has inversion symmetry,⁶ we calculated the topological invariants from the symmetry of Bloch functions at the Brillouin zone points by employing the newly developed interface between siesta and wannier 90.³⁴ 

In this calculation, all the atomic coordinates were relaxed until the forces were smaller than 0.01 eV/Å, and the stress tensor components were below 0.0001 eV/ų. For spin–orbit coupling (SOC) in siesta, we employed the recipe described in Refs. 35 and 36.

In Fig. 1, we illustrate Bi₂Se₃ with a rhombohedral crystal structure. In addition, it has a space group $D3d5(R3¯m)$⁠. Usually, it is displayed in quintuple layers with an Se₁–Bi₁–Se₂–Bi₁′–Se₁′ arrangement. Bi₂Se₃ portrays layered structures with a triangular lattice within its single layer. It is important to note that this material consists of five-atom layers arranged along the z-direction, known as quintuple layers (QLs). A QL consists of five atoms with two equivalent Se atoms, two equivalent Bi atoms, and a third Se atom. Some scholars¹¹ denote the two Se atoms as Se₁ and Se_1′, the two Bi atoms as Bi₁ and Bi_1′, while the third Se as Se₂.

It is usually important to perform geometrical relaxations in any first-principles calculation to avoid errors and for an accurate prediction of other dependent quantities such as the electronic properties that are extremely sensitive. Table I presents the fully relaxed lattice parameters of Bi₂Se₃ with a rhombohedral crystal. Our calculations, with and without taking into account the effects of spin–orbit coupling, are in good agreement with the available theoretical and experimental results. It is clear that our lattices are slightly overestimated. This is a well-known occurrence when employing GGA, as described in Ref. 37.

When using bare GGA, we found that the lattices a and c were overestimated by 1.47% and 0.73%, respectively, in comparison to the experimental values. On employing spin–orbit coupling, we realized that it increased the error in lattice parameters. This kind of discrepancy can easily be corrected by using van der Waals (vdW) corrections, as employed in Ref. 38. The inclusion of vdW in such calculations will definitely give matching experimental results as far as the lattices are concerned.

We calculated the band structure of Bi₂Se₃ along high-symmetry points Γ-Z-F-Γ-L in the first Brillouin zone. In Fig. 2(a), we set the Fermi level at 0 eV, and it is represented by the red dotted line. For a plain GGA calculation, an energy separation at Γ of 0.05 eV is observed between the top of the valence band and the bottom of the conduction band. This is also in agreement with 0.08 eV obtained by Aguilera and co-workers⁴⁰ when using LDA.⁴¹ This hints to the fact that Bi₂Se₃ is a direct bandgap material. This energy gap was severely underestimated due to the DFT limitation that arises from the approximations used in the exchange and correlation functional.

When we considered the effect between the orbital motion and electron spin as shown in Fig. 2(b), we found the energy gap to be of the order of 0.28 eV but no longer direct as also reported in the literature.⁴⁰ The gap was found between Z and F. This is consistent with the previous DFT and experimental findings, as shown in Table II. From Fig. 2(b), we show that spin–orbit coupling has a significant impact on the band structure and alters the bandgap at the Γ point. Since we no longer have a direct gap, it is noted that SOC induces gap opening between the conduction and the valence band, as shown in Fig. 4 and elaborated in Fig. 5. The inversion seen illustrates non-trivial topology in Bi₂Se₃ due to the band inversion in momentum space.

Another exciting feature in this material is that the topological surface states for this crystal are seen to be simple and are described by a single gap-less Dirac cone at the $k→$ = 0, Γ point in the surface Brilloiun zone.⁴⁰ 

One important feature in Fig. 2(b) is the lifting of the degeneracy of one-electron energy levels in the band structure. This has been seen to be significantly dominant at the Γ point. Without the employment of spin–orbit coupling as depicted in Fig. 2(a), we found doubly degenerated bands. This type of degeneracy is usually broken by a spin-dependent term in the Hamiltonian, as described in Refs. 35 and 36, of the siesta method. Apart from the increased bandgap, the bandwidth also increased.

An investigation of the total and partial densities of states always helps to clarify the nature of the bandgaps. In principle, we get information about the origin of various orbitals in the band structure. From Fig. 3, it can clearly be seen that the s-orbitals of Bi and Se contribute majorly in the core states, while the p-orbitals of Se are seen to dominate the valence states. A close look at the conduction states shows domination by the p-orbitals of Bi. At the Fermi level, the p-orbitals of Se and Bi are fully responsible for the properties seen herein. A similar orbital arrangement in Bi₂Se₃ was reported in Ref. 38. There is also a very strong hybridization in the Se p and Bi p states.

A second look at Fig. 3 shows that the spin–orbit coupling suppresses the density of states. The total density of states from ≈30 states/eV in the scalar relativistic scenario decreases to ≈3 states/eV in a fully relativistic case. This should not raise an alarm since similar observations were seen in Ref. 42. One common effect of spin–orbit coupling is that it always tends to lower the energy of a system.⁴³ 

It is well known that in Bi₂Se₃, the topological index distinguishing ordinary insulating from topological insulating behavior is controlled by band inversion at the Γ point,⁴⁴ as shown in Fig. 5. The goal of this work was to demonstrate band inversion in Bi₂Se₃ and, thereafter, seek to answer the question of the origin of band inversion. The usual strategy to get a topological insulator, and already used by Kane and Mele,¹⁶ is to induce the wave function twist using spin–orbit coupling (SOC). This is known for Bi₂Se₃ but may be present in other materials.

Without spin–orbit coupling (SOC), as shown in Fig. 5(a), the “conduction band” in green is majorly made of Bi p_z orbitals with a little admixture of Se p_z and the “valence band”, shown in red, is made of Se p_z orbitals. However, due to DFT underestimation, the bands tend to overlap, giving the system a gap of 0.05 eV, as explained in Sec. III B. When spin–orbit coupling (SOC) is included in Fig. 5(b), a gap opens at Γ point, and now, we have a proper conduction band that has contributions from the band that made the valence band originally (red) and vice versa. This is called a band inversion.

A band inversion like this may suggest that the material has a topological order, but the only way to confirm it is by calculating the topological invariant and later to calculate the Weyl chirality to prove the topological nature. The actual cause of band inversion in materials is still a subject of discussion as far as scholars are concerned.

Following the recommendations in Ref. 26, it is argued that band inversion is not necessarily induced by spin–orbit coupling (SOC) but may also occur when the strength of some other external parameter such as structural distortion increases. This is demonstrated in Fig. 6(b). We found out that after slight distortions of the lattice (a by 0.47% and c by 2.36%), band inversion occurred at F and Γ points. Unfortunately, the bands did not open as portrayed in Fig. 5(b), confirming that spin–orbit is necessary for the opening of the gap in our case. With the lattice distortion in the system, we also noted an increase in the bandwidth. The lifting of degeneracy around the Γ point was also seen just as it was predicted in Fig. 5(b), giving a notion that degeneracy can be lifted using lattice distortions.

The authors herein argue that while what has been illustrated in Fig. 6(b) is not enough to have a topological insulator, it is still vital to open the bandgap to get to a situation similar to the one depicted in Fig. 4 or Fig. 5(b), which usually still needs spin–orbit coupling (SOC).

Often, it is argued that spin–orbit is responsible for band inversion in topological Bi₂Se₃. This is not true at all based on our lattice strained band inversion demonstrated in Fig. 6(b). Since we performed a fully relativistic optimization, an elongation of the c lattice was observed, while the a lattice shrank. We had to repeat a scalar relativistic calculation with a distorted lattice and similar effects were seen. The point herein is that it is possible to witness band inversion (i) in materials whose spin–orbit coupling energy is too low, as reported in Ref. 26, and (ii) in a strained lattice similar to the case of Fig. 6(b).

It is therefore prudent to note that band inversion is caused by a couple of factors: (i) scalar relativistic effects such as mass velocity and Darwin terms and (ii) lattice distortion in the crystal. Spin–orbit coupling in Bi₂Se₃ is only vital in opening the gap and making it topological. The mass velocity and Darwin terms in this case can be written as −P⁴/8m³c² and $−ℏ2/4m2c2⋅dV/dr⋅∂/∂r→$⁠, respectively, as elaborated in Ref. 45. We note herein that even when no spin–orbit is considered, the scalar relativistic terms are included in our scheme since we have a heavy element such as Bi. Such contributions will not split bands, but they are rather dependent on the angular momentum behavior of the band, particularly for p states.⁴⁶ 

We note that the response of Bi₂Se₃ toward a given perturbation is characterized by topological invariants.^18,19 This means that a given topological invariant describes the fundamental properties of the band structure at preferred k-points. Usually, it is used to describe strong topological insulators and if present, the system will have topologically protected bands at the surface. Since we have inversion symmetry in Bi₂Se₃, we obtained the topological invariants from the symmetry of the Bloch function at six special Brillouin zone points,

where the reciprocal lattice vectors are b_i, while n_i = 0, 1. If we consider ψ_i,n to be the nth occupied Bloch function at γ_i, we defined the symmetry function as

where Θ is the inversion operator. Since we know the symmetry functions in Bi₂Se₃, it follows that the strong topological invariant ν is given by

There also exist three weak invariants ν₁, ν₂, and ν₃, given by the products of four δ_i’s for which Γ_i reside in the same plane,

Usually, the weak invariants are not so much robust toward atomic disorder. We found out that the $Z2$ topological index in Bi₂Se₃ was [1, 0, 0, 0], a signature of a strong topological material.

In this work, we have demonstrated the fundamental band inversion in topological Bi₂Se₃, which is a result of both scalar relativistic and fully relativistic effects. Our work indicates a critical role played by spin–orbit coupling in such materials on the opening of the band structure. Our bandgap energy of 0.28 eV was also in a perfect agreement with the existing scholarly work. One strong finding in this work is the fact that spin–orbit coupling is not the driving force for band inversion in Bi₂Se₃, but scalar relativistic terms and lattice distortions do. It is also important to note that the basic features of topological materials can equally be probed and understood by the use of simple tight-binding models using the newly developed interface between siesta and wannier 90, as illustrated in Ref. 47. This will give us room to deal with complex dynamics including electron–phonon and electron–electron interactions that will permit us to simulate polaron charge transport or excitonic phenomena in Bi₂Se₃. It would be interesting to see if band inversion can still be observed at the Γ-point if higher theories such as hybrid functionals and GW are employed.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The authors gratefully acknowledge the computer resources, technical expertise, and assistance provided by the Centre for High Performance Computing (CHPC - MATS862 and MATS0712), Cape Town, South Africa.