Relativistic k.p Hamiltonians for centrosymmetric topological insulators from ab initio wave functions

We present a method to microscopically derive a small-size k$\cdot$p Hamiltonian in a Hilbert space spanned by physically chosen ab initio spinor wave functions. Without imposing any complementary symmetry constraints, our formalism equally treats three- and two-dimensional systems and simultaneously yields the Hamiltonian parameters and the true $\mathbb{Z}_2$ topological invariant. We consider bulk crystals and thin films of Bi$_{2}$Se$_{3}$, Bi$_{2}$Te$_{3}$, and Sb$_{2}$Te$_{3}$. It turns out that the effective continuous k$\cdot$p models with open boundary conditions often incorrectly predict the topological character of thin films.

We present a method to microscopically derive a small-size k·p Hamiltonian in a Hilbert space spanned by physically chosen ab initio spinor wave functions. Without imposing any complementary symmetry constraints, our formalism equally treats three-and two-dimensional systems and simultaneously yields the Hamiltonian parameters and the true Z2 topological invariant. We consider bulk crystals and thin films of Bi2Se3, Bi2Te3, and Sb2Te3. It turns out that the effective continuous k·p models with open boundary conditions often incorrectly predict the topological character of thin films. Electronic structure of topological insulators (TIs) has been in focus of theoretical research regarding linear response, transport properties, Hall conductance, and motion of Dirac fermions in external fields [1,2]. These problems call for a physically justified model Hamiltonian of small dimension. As in semiconductors, it is thought sufficient that the model accurately reproduces the TI band structure near the inverted band gap [3]. The desired Hamiltonian is derived either from the theory of invariants [4] or within the k·p perturbation theory using the symmetry properties of the basis states [5].
In Ref. [3], along with the pioneering prediction of the topological nature of Bi 2 Se 3 , Bi 2 Te 3 , and Sb 2 Te 3 , a 4band Hamiltonian was first constructed from the theory of invariants, which is presently widely used to analyze the properties of bulk TIs as well as their surfaces and thin films [6][7][8][9][10][11][12][13][14]. The Hamiltonian parameters in Ref. [3] were obtained by fitting ab initio band dispersion curves. Later, an attempt was made [15] to recover the Hamiltonian of Ref. [3] by a k·p perturbation theory with symmetry arguments and to derive its parameters from the ab initio wave functions of the bulk crystals. Furthermore, in Ref. [15] the effective Landé g-factors for the Zeeman splitting [4,5] were introduced within the k·p theory of TIs.
To analyze how the properties of thin films are inherited from the bulk TI features, effective continuous models have been developed: they are based on the substitution k z → −i∂ z (originally introduced for slowly varying perturbations [16]) in the Hamiltonian of Ref. [3] and on the imposition of the open boundary conditions [15,[17][18][19]. These models predict a variety of intriguing phenomena at surfaces, interfaces, and thin films of TIs [20][21][22][23]. A fundamental issue here is the topological phase tran-sition between an ordinary 2D insulator and a quantum spin Hall insulator (QSHI). Apart from the theoretical prediction, the model parameters are fitted to the measured band dispersion to deduce the topological phase from the experiment [24,25]. By analyzing the signs and relative values of the parameters of the empirically obtained effective model a judgement is made on whether the edge states would exist in a given TI film, the logic being similar to that of Ref. [26]: The valence band should have a positive and conduction band a negative effective mass.
In order to avoid any ambiguity in deriving the model Hamiltonian and to treat 3D and 2D systems within the same formalism, one needs an ab initio and internally consistent scheme that realizes the k·p theory with the full inclusion of the microscopic structure of the system and generates a compact and physically transparent form of the Hamiltonian of a given size. A few attempts have been recently undertaken to predictably construct model Hamiltonians for classical bulk semiconductors [27] and graphene-based systems [28].
Here, we report a method to microscopically derive the relativistic Hamiltonian H kp accurate up to the second order in k from the spinor wave functions obtained with the all-electron full-potential extended linearized augmented plane wave method (ELAPW). The size of H kp is determined by the dimension of the subspace spanned by the physically chosen basis states. The form of the Hamiltonian is dictated by the symmetry of the wave functions unitary transformed to diagonalize the z-component of the total angular momentum and by a universal prescription to choose their phases. Here we apply this approach to centrosymmetric bulk crystals as well as to thin films of Bi 2 Se 3 , Bi 2 Te 3 , and Sb 2 Te 3 up to six quintuple layers (QLs). For each film, we calculate the topological invariant, indicating whether it is a QSHI. We conclude on the validity of the effective models by comparing the predictions by the k·p Hamiltonian with the actual properties of the film. Furthermore, within our approach, we derive for the first time the k·p Zeeman term for the films.
We construct the model Hamiltonian as a second-order k·p expansion around the point k = 0. To avoid any ambiguity, we obtain the expansion coefficients directly from ab initio eigenfunctions at Γ. For systems with inversion symmetry the energy bands E n are doubly degenerate with two orthogonal wave functions Ψ n1 and Ψ n2 that are parity eigenfunctions at the time reversal invariant momenta (TRIM). The k·p Hamiltonian is represented in the basis of these functions in terms of the matrix elements Ψ ni |π|Ψ mj of the velocity opera- [29]. Here σ is the vector of Pauli matrices and V (r) is the crystal potential. In the k·p expansion H kp = H (0) +H (1) +H (2) , the zeroorder term is just the band energy, H where α, β = x, y, z, and see Refs. [4,5]. Here m and n number the degenerate Kramers pairs, and i and j number the members of a pair. The index n ′ runs over all the bands excluding those forming the k·p basis (Löwdin's partitioning). Thus, when the dimension of H kp equals the dimension of the original full Hamiltonian the second order term H (2) vanishes [30] (hereafter, we refer to this case as the full-size k·p calculation). The ab initio band structure was obtained with the ELAPW method [31] using the full potential scheme of Ref. [32] within the local density approximation (LDA). The spin-orbit interaction is treated by a second variation method [33] including the scalar-relativistic bands up to at least 300 eV. This ensures a good convergence of the inverse effective mass, with a deviation from the second derivative of the E(k) curves within 3%. The experimental crystal lattice parameters were taken from Ref. [35] with the LDA relaxed atomic positions of Refs. [36][37][38]. Figure 1 compares the ab initio bands with those obtained by diagonalizing our k·p Hamiltonians of small size (4-and 8-band) and of full size. Note that the fullsize k·p calculation highly accurately reproduces the true bands: the error grows as k 2 [39], and at the Brillouin zone (BZ) boundary it is within 150 meV. For each Kramers-degenerate level n, the spinor wave functions Ψ ni form a two-dimensional basis. Numeri- cally obtained functions are arbitrarily ordered and have unphysical phases, which, however, affect the structure of H kp non-diagonal terms. In order to keep the same physically motivated ordering and to align the phases in different calculations we first transfer to the basis that diagonalizes the z-component of the total angular momentum J = L + S in the atomic sphere that has the largest weight in the n-th band (see Figs. S5-S7 in the Supplemental Material (SM) [34]). This establishes the numeration of the wave functions Ψ n1(2) → Ψ n↑(↓) . Next, we choose the phases of the new basis functions such that they become explicitly Kramers conjugate: Ψ n↓ =T Ψ n↑ , whereT = Kiσ y is the time reversal operator and K is the complex conjugation operator. Finally, for two pairs of different parity, n-th and m-th, we turn the phases such that iπ x(z) n↑m↓ be real. For the bulk TIs we choose the basis of four states Ψ v↑ , Ψ v↓ , Ψ c↑ , Ψ c↓ , where v and c stand for the valence and conduction band, respectively. This leads to the Hamiltonian [40] x + k 2 y . The 3D TI Hamiltonian (1) is the same (to within a unitary transformation) as in Refs. [3,15,[17][18][19] (the explicit matrix form is presented by Eq. (S1) in the SM [34]). The Pauli matrix τ operates in valence-conduction band space, whereas σ refers to the total angular momentum J. In Eq. (1) a direct product of these matrices is implied.
For the bulk crystals, the parameters in Eq. (1) are listed in Table I, and the eigenvalues E(k) of the result-  (1) is that they are very sensitive to details of the crystal geometry, as is the ab initio band structure [36][37][38]: even a small variation in atomic positions leads to considerable changes of the parameters (see Table S1 in the SM [34]). Furthermore, in all the 3D systems considered, see Table I, the parameters of H kp turned out to meet the conditions of the existence of the topological surface states [18]: the diagonal dispersion term M z( ) is positive and is larger than the electron-hole asymmetry: M z( ) > |C z( ) |, although C z and C are not negligible contrary to the assumption in Ref. [12]. Thus, our ab initio k·p Hamiltonian correctly predicts the topological character of these crystals in accord with the Z 2 topological invariant ν 3D obtained from the parities of the wave functions at the TRIM points [41]. We now use our second order perturbation theory to calculate the effective g-factors entering the Zeeman term (see Eq. (S2) in the SM [34]) that appears in the presence of static magnetic field [4,5]: The most important is that the calculated values are one or even two orders of magnitude larger than the free electron g-factor, g 0 ≈ 2. (The values obtained with the four-band k·p method are listed in Table S1 of the SM [34] for all the 3D TIs studied.) This result accords with the recent spin resonance measurements of the effective g-factor in Bi 2 Se 3 [42]: In contrast to the bulk TIs, for finite-thickness TI films an ambiguous behavior is observed. For a 2D system, in the basis Ψ slab v↑ , Ψ slab c↓ , Ψ slab c↑ , Ψ slab v↓ our k·p Hamiltonian reads (cf. Refs. [26,43]): where C = C 0 + C k 2 , M = M 0 + M k 2 , and the operator τ refers now to two decoupled sets of massive Dirac fermions. The last term in Eq. (2) ensures the characteristic spin-orbital texture of the TI surface states [44,45]. Figure 2 shows the film-thickness dependence of the parameters of the Hamiltonian (2) for the TIs considered (the plotted values are listed in Tables S2-S4 of the SM [34]). In contrast to the bulk TIs, for all the thicknesses the 4-band k·p spectrum does not have the absolute gap, see the red lines in Figs. 1(d)-1(f) for 2QL films and Figs. S2-S4 of the SM [34] for other thicknesses. Note that only for two QLs and only for Bi 2 Te 3 does the eightband Hamiltonian (blue lines) provide a quality close to that achieved in the 3D case.
Velocity V and electron-hole asymmetry C converge quite fast with the film thickness [Figs. 2(b) and 2(d)], demonstrating just slight changes starting from 3 QLs. The 6QL value of V can thus be compared with the Fermi velocity from the effective models (see, e.g., Ref. [17]) in the large-thickness limit. In this limit the  Fig. 2(d)]. The parameter C is positive everywhere, and it is notably larger than the absolute value of the diagonal dispersion term M .
For the same thickness, the parameter M may have different sign for different TIs, whereas for a given material M is found to "oscillate" with the number of QLs. It turns out that these oscillations do not correlate with the Z 2 topological invariant ν 2D obtained from the parities of the original wave functions at the TRIMs of the 2D BZ. This becomes evident from a comparison with the behavior of the gap parameter ∆, whose absolute value |∆| = −2M 0 yields the gap width atΓ and sign depends on ν 2D , with ∆ being negative for a topologically nontrivial film.
The Sb 2 Te 3 film becomes a 2D TI at 3 QLs and preserves this property up to 6 QLs, Fig. 2(c). For Bi 2 Te 3 , the 2 QL film is non-trivial, then 3 and 4 QLs are trivial, and 5 and 6 QLs are again non-trivial. The films of Bi 2 Se 3 are, on the contrary, trivial for all the thicknesses, while the k·p effective models predict them to be a QSHI at some of the thicknesses, see e.g. Ref. [17]. It should also be noted that for the same film the true invariant may depend on details of the crystal geometry and even on the band structure method, including the choice of exchange-correlation potential, see Refs. [46][47][48][49][50][51][52]. We emphasize that here the topological invariant and the k·p parameters are fully consistent because they are derived from the same band structure.
According to the effective continuous models [17,26], the relation between C , M 0 , and M one finds in Fig. 2 clearly predicts the absence of edge states. This means that a few-band k·p Hamiltonian does not provide a gen- eral and certain criterion of the topological character of 2D systems. Because the electron-hole asymmetry term is sometimes neglected in topological analysis [1], it is instructive to consider in more detail the 2QL films of Bi 2 Se 3 and Bi 2 Te 3 -the two thinnest films for which the sign of M correlates with the actual ν 2D . In Fig. 3, we see that with varying the van-der-Waals spacing (expansion for Bi 2 Te 3 and compression for Bi 2 Se 3 ) the parameters V , C and ∆ change steadily, and in Bi 2 Te 3 a transition from QSHI to the trivial state occurs (at 18% ∆ becomes positive), and at the same time M becomes negative, again following the true indicator ν 2D .
Finally, let us consider the effective-mass contribution to the g-factor for the films. In our approach, the static magnetic field B leads to the following Zeeman term: where g α = (g v α + g c α )/2 and ∆g α = (g v α − g c α )/2. A novel feature in Eq. (3) is the second term that contains the z-component of the cross product [σ ×B]. It resembles the Zeeman term for inversion-asymmetric quantum wells, where the "spin-momentum locking" term is also present [4]. As follows from Figs. 4(b) and 4(d), the parameters ∆g z and ∆g may "oscillate" with the thickness, and for Bi 2 Se 3 and Bi 2 Te 3 they tend to zero with increasing thickness. As a result, the leading contribution comes from the "conventional term" with g z and g , which at 6 QLs is already well converged, Figs. 4(a) and 4(c). For Sb 2 Te 3 , ∆g z and especially ∆g are rather big at 6 QLs, so the relevant term in H slab kp,Z should be taken into account at least up to this thickness. Note that in moving from 1 to 6 QLs ∆g z becomes negative for the first time when the film becomes a 2D TI, thus demonstrating a correlation with the topological invariant (see also the effect of expansion/compression for the 2 QL films in Fig. S8 of the SM [34]).
To summarize, we have developed a fully ab initio k·pperturbation approach to generate model Hamiltonians of a desired size. This ensures a physically meaningful behavior of the model Hamiltonian parameters with the continuously varying geometry (van-der-Waals spacing) and for different number of the building layers. By applying our approach to Bi 2 Se 3 , Bi 2 Te 3 , and Sb 2 Te 3 films, we have demonstrated that the widely used effective continuous models are not able to systematically predict the values and often even the relative sign of the model parameters. The failure to infer the general and certain criterion from H kp stems from its fundamental limitation: the topological properties of a crystal cannot be unambiguously determined from the behavior of a few bands in the vicinity of k = 0, even though the band inversion occurs just at that point.