⇦Chinese Physics Letters, 2021, Vol. 38, No. 7, Article code 077105 A Programmable k$\cdot$p Hamiltonian Method and Application to Magnetic Topological Insulator MnBi$_2$Te$_4$ Guohui Zhan (占国慧)1†, Minji Shi (施敏吉)1†, Zhilong Yang (杨志龙)1, and Haijun Zhang (张海军)1,2* Affiliations 1National Laboratory of Solid State Microstructures, School of Physics, Nanjing University, Nanjing 210093, China 2Collaborative Innovation Center of Advanced Microstructures, Nanjing University, Nanjing 210093, China Received 8 May 2021; accepted 2 June 2021; published online 18 June 2021 Supported by the Fundamental Research Funds for the Central Universities (Grant No. 020414380185), the Natural Science Foundation of Jiangsu Province (Grant No. BK20200007), the National Natural Science Foundation of China (Grant Nos. 12074181 and 11834006), and the Fok Ying-Tong Education Foundation of China (Grant No. 161006).
†G. Zhan and M. Shi contributed equally to this work.
*Corresponding author. Email: zhanghj@nju.edu.cn Citation Text: Zhan G H, Shi M J, Yang Z L, and Zhang H J 2021 Chin. Phys. Lett. 38 077105    Abstract In the band theory, first-principles calculations, the tight-binding method and the effective $\boldsymbol{k}$$\cdot$$\boldsymbol{p}$ model are usually employed to investigate electronic structures of condensed matters. The effective $\boldsymbol{k}$$\cdot$$\boldsymbol{p}$ model has a compact form with a clear physical picture, and first-principles calculations can give more accurate results. Nowadays, it has been widely recognized to combine the $\boldsymbol{k}$$\cdot$$\boldsymbol{p}$ model and first-principles calculations to explore topological materials. However, the traditional method to derive the $\boldsymbol{k}$$\cdot$$\boldsymbol{p}$ Hamiltonian is complicated and time-consuming by hand. We independently developed a programmable algorithm to construct effective $\boldsymbol{k}$$\cdot$$\boldsymbol{p}$ Hamiltonians for condensed matters. Symmetries and orbitals are used as the input information to produce the one-/two-/three-dimensional $\boldsymbol{k}$$\cdot$$\boldsymbol{p}$ Hamiltonian in our method, and the open-source code can be directly downloaded online. At last, we also demonstrated the application to MnBi$_2$Te$_4$-family magnetic topological materials. DOI:10.1088/0256-307X/38/7/077105 © 2021 Chinese Physics Society Article Text Recently, research of topological states and topological materials, such as topological insulators, topological semimetals and topological superconductors, has become an important topic in condensed matter physics, and great progresses have been achieved in both experiment and theory.[1–3] The developed topological band theory with first-principles calculations played a key role in the past years. For example, topological invariants, such as $\mathcal{Z}_2$, were defined through the band theory.[4] The concept of band inversion of topological states stemmed from the band theory.[5,6] Also, the topological boundary states were successfully predicted through first-principles calculations.[7] However, the first-principles band structures seem like a black box with much-hidden information, which results in a barrier to deeply understand the essential physical pictures of topological states. Differently, the $\boldsymbol{k}$$\cdot$$\boldsymbol{p}$ Hamiltonian has a simple form and a clear physical picture, which is a necessary and useful supplement for the first-principles calculations.[8–11] With the $\boldsymbol{k}$$\cdot$$\boldsymbol{p}$ Hamiltonian, Fu predicted an unconventional hexagonal warping term in surface states of topological insulator Bi$_2$Te$_3$.[12] Therefore, combination of first-principles calculations and the $\boldsymbol{k}$$\cdot$$\boldsymbol{p}$ Hamiltonian has become a standard paradigm to theoretically study topological states and topological materials.[6,7,13] Based on the group theory and the $\boldsymbol{k}$$\cdot$$\boldsymbol{p}$ perturbation theory, Liu et al. provided an instructive demonstration showing how to construct the $\boldsymbol{k}$$\cdot$$\boldsymbol{p}$ Hamiltonian of Bi$_2$Se$_3$-family three-dimensional topological insulators.[11] However, though the derivation is a standard method, the detailed process is quite troublesome and time-consuming by hand, especially for dealing with high-order terms of $\boldsymbol{k}$$\cdot$$\boldsymbol{p}$ Hamiltonians. Therefore, it is practicable and highly necessary to develop automated programmable methods to efficiently construct $\boldsymbol{k}$$\cdot$$\boldsymbol{p}$ Hamiltonians. In this context, there have been some innovative proposes, such as kdotp-symmetry code developed by Gresch et al.,[14,15] Qsymm Python package developed by Varjas et al.[16,17] In this Letter, based on group theory, we independently develop a programmable algorithm to construct $\boldsymbol{k}$$\cdot$$\boldsymbol{p}$ Hamiltonians for all one-/two-/three-dimensional materials (the open-source code can be download from https://github.com/shimj/Model-Hamiltonian). With this automated algorithm, only the crystal symmetry and atomic orbitals involved are needed to produce the $\boldsymbol{k}$$\cdot$$\boldsymbol{p}$ Hamiltonian, which effectively avoids time-consuming calculations and latent mistakes. The schematic of our programmable method is shown in Fig. 1. We have successfully applied our method to MnBi$_2$Te$_4$-family magnetic topological materials.[18–20] Basic $\boldsymbol{k}$$\cdot$$\boldsymbol{p}$ Hamiltonian Theory. In the band theory, the wavefunction in a periodic lattice is written as the Bloch wavefunction, $\psi_{\boldsymbol{k}}(\boldsymbol{r})=u_{\boldsymbol{k}}(\boldsymbol{r})e^{i\boldsymbol{k}\cdot\boldsymbol{r}}$. Here, $u_{\boldsymbol{k}}(\boldsymbol{r})$ is a periodic function with $u_{\boldsymbol{k}}(\boldsymbol{r})=u_{\boldsymbol{k}}(\boldsymbol{r}+\boldsymbol{R})$, where $\boldsymbol{R}$ indicates the lattice vector. By substituting the Bloch wavefunction $\psi_{\boldsymbol{k}}$ into the Schrödinger equation, we have $$ \mathcal{H}_{\boldsymbol k}u_{\boldsymbol k} = E_{\boldsymbol k}u_{\boldsymbol k},~~ \tag {1} $$ where the Hamiltonian has the $\boldsymbol{k}$-dependent form $\mathcal{H}_{\boldsymbol k}=p^2/2m+V+\hbar\boldsymbol{k}\cdot \boldsymbol{p}/m_0+\hbar^2k^2/2m$. Generally, the representation matrix of Hamiltonian $\mathcal{H}_{\boldsymbol k}$ near a $\boldsymbol{k_0}$, with $\boldsymbol{k}=\boldsymbol{k_0}+\delta\boldsymbol{k}$, can be approached by Taylor series as $$ \mathcal{H}_{\boldsymbol{k_0}+\delta\boldsymbol{k}}\doteq\sum_{\alpha+\beta+\gamma\leq n} (\delta k_x)^{\alpha}(\delta k_y)^{\beta}(\delta k_z)^{\gamma}\varGamma_{\alpha\beta\gamma},~~ \tag {2} $$ where $\varGamma_{\alpha\beta\gamma}$ are constant matrices, $n$ indicates the order of Taylor series expansion. In the following, we will derive the representation matrix of $\mathcal{H}_{\boldsymbol k}$ from the commutation relation between Hamiltonian and symmetry operations.
cpl-38-7-077105-fig1.png
Fig. 1. Workflow for constructing the $\boldsymbol k \cdot \boldsymbol p$ effective Hamiltonian. The selected generators of the group and atomic orbitals are the input information to run the program. In the process, the key point is to obtain the representation matrix $\mathcal{F}$ and $\mathcal{M}$. By solving $(\mathcal{F}\otimes\mathcal{M}-I)\alpha=0$, the Hermitian matrices can be obtained from the solution $A$.
When a periodic condensed matter system preserves an $S$ symmetry, its Hamiltonian $\mathcal{H}$ satisfies the commutation relation $[\mathcal{H},S]=0$. In $\boldsymbol{k}$-space, we have $$ S\mathcal{H}_{\boldsymbol k}S^{-1} = \mathcal{H}_{\gamma\boldsymbol k},~~ \tag {3} $$ where $$ \gamma = \begin{cases} \tilde{g}, & S = Q,\\ -1, & S = \mathcal{T},\\ -\tilde{g}, & S = Q\mathcal{T}. \end{cases}~~ \tag {4} $$ Here, $Q$ represents an operator corresponding to the spatial operation $g=\{\tilde{g} | \tau\}$, $\mathcal{T} = (-i\sigma_y\otimes I) \mathcal{K}$ is the time reversal operator, $\tilde{g}$ is a rotation or reflection operator, and $\tau$ is a translation vector. It is noteworthy that $2\pi$-rotation for a spinful system is not included in our discussion, since it leads to an identity $(-1)\mathcal{H}_{\boldsymbol k}(-1)=\mathcal{H}_{\boldsymbol k}$ without giving any constraint to the Hamiltonian $\mathcal{H}_{\boldsymbol k}$. In principle, any complete set of lattice-periodic functions, for example, $\{e^{i\boldsymbol G\cdot \boldsymbol r}\}$, can be used as basis to accurately solve the Schrödinger equation $\mathcal{H}_{\boldsymbol k}|u_{\boldsymbol k}\rangle = E_{\boldsymbol k}|u_{\boldsymbol k}\rangle$. Since we are interested in the bands around a certain point $\boldsymbol k_0$ in the Brillouin zone (BZ), it allows us to use $\{|u_i^{\boldsymbol k_0}\rangle\}$ as a basis to approximately solve the Schrödinger equation in the framework of the perturbation theory, where $i$ indicates the band index. Once we know the symmetry transformations of these states, the matrix representation of $\mathcal{H}_{\boldsymbol k}$ under this basis can then be determined in form. Around a high-symmetry point $\boldsymbol k_0$, for a symmetry operation $S$ meeting $\gamma \boldsymbol k_0 = \boldsymbol k_0$, we can rewrite Eq. (3) as $$ S \mathcal{H}_{\boldsymbol k_0+\delta \boldsymbol k}S^{-1} = \mathcal{H}_{\boldsymbol k_0+\gamma\delta\boldsymbol k}.~~ \tag {5} $$ First of all, the symmetry operations under consideration have to keep the linear space spanned by $\{|u_i^{\boldsymbol k_0}\rangle\}$ invariant. Otherwise, we must either drop some symmetry operations or add more states into the basis. With the representation matrix $D(S)$ for symmetry operations $S$ defined as $$\begin{alignat}{1} &S^u (|u^{\boldsymbol k_0}_1\rangle,|u^{\boldsymbol k_0}_2\rangle,\dots)=(|u^{\boldsymbol k_0}_1\rangle,|u^{\boldsymbol k_0}_2\rangle,\dots)D(S^u),~~ \tag {6a} \end{alignat} $$ $$\begin{alignat}{1} &S^a (|u^{\boldsymbol k_0}_1\rangle,|u^{\boldsymbol k_0}_2\rangle,\dots)=(|u^{\boldsymbol k_0}_1\rangle,|u^{\boldsymbol k_0}_2\rangle,\dots)D(S^a)\mathcal{K},~~ \tag {6b} \end{alignat} $$ where we have explicitly marked the symmetry operation $S$ by $u$ for “unitary” and $a$ for “antiunitary”, Eq. (5) immediately gives $$\begin{alignat}{1} &D(S^u)\mathcal{H}_{\boldsymbol k_0}(\gamma^{-1}\delta\boldsymbol k)D(S^u)^{-1} = \mathcal{H}_{\boldsymbol k_0}(\delta\boldsymbol k),~~ \tag {7a} \end{alignat} $$ $$\begin{alignat}{1} &D(S^a)\mathcal{H}_{\boldsymbol k_0}(\gamma^{-1}\delta\boldsymbol k)^*D(S^a)^{-1} = \mathcal{H}_{\boldsymbol k_0}(\delta\boldsymbol k),~~ \tag {7b} \end{alignat} $$ where $\mathcal{H}_{\boldsymbol k_0}(\delta\boldsymbol k)$ denotes $\mathcal{H}_{\boldsymbol k_0+\delta\boldsymbol k}$. LCAO Basis. In group theory, for a linear space invariant under a group, one can always find a set of basis in which all the vectors can form irreducible representations of that group and are orthogonal with each other. For a general vector, by analyzing its components, one can easily find all other vectors related by the group operations and all the irreducible representations. Since the first-principles calculations can give the atomic-orbital projection of the wavefunction, it would be useful to set up a database containing the representations of all low-level linear combinations of atomic orbitals (LCAOs). First of all, we ignore the spin. The periodic part of linear combination of atomic orbitals has the form $$\begin{align} u^{\boldsymbol k}_{m\alpha} (\boldsymbol r)= \frac{1}{\sqrt{N}}\sum_{\boldsymbol R_n}e^{-i\boldsymbol k\cdot(\boldsymbol r-\boldsymbol r_\alpha-\boldsymbol R_n)}\phi_m(\boldsymbol r-\boldsymbol r_\alpha - \boldsymbol R_n),~~ \tag {8} \end{align} $$ where $\alpha$ and $m$, respectively, refer to the $\alpha$-th atom in a unit cell and its $m$-th atomic orbital. Here we also define $\phi'_m(\boldsymbol r)=e^{-i\boldsymbol k\cdot \boldsymbol r}\phi_m(\boldsymbol r)/\sqrt{N}$. It should be noted that if we collect the periodic wavefunctions for all the same atoms in a unit cell and all the atomic orbitals with the same quantum number $\ell$ to span a linear space, it must be invariant under symmetry operations. To simplify the calculation, we analyze the indices $m$ and $\alpha$ separately. So, we need two linear spaces: a complex linear space $V_1$ spanned by $\{\phi'_m(\boldsymbol r)\}$ and $V_2=\mathbb{R}^{\alpha_{\max}}$, where $\alpha_{\max}$ means the number of atoms in our consideration of a unit cell and the vectors in it indicate the spatial distribution of atoms. Then with the correspondence as $$\begin{align} &\phi'_m(\boldsymbol r)\otimes (\overbrace{0,\dots ,0}^{\alpha-1},1,0,\dots ,0)^{\rm T} \leftrightarrow \sum_{\boldsymbol R_n}\phi'_m(\boldsymbol r-\boldsymbol r_\alpha-\boldsymbol R_n),~~ \tag {9} \end{align} $$ we know that the tensor product of these two spaces $V_1\otimes V_2$ is isomorphic to the space $V_0$ spanned by $\{u^{\boldsymbol k}_{m\alpha}(\boldsymbol r)\}$ for a certain $\boldsymbol k$. Finally, with the straightforward correspondence, $$\begin{alignat}{1} &S\sum_{\boldsymbol R_n}\phi'_m(\boldsymbol r-\boldsymbol r_\alpha-\boldsymbol R_n) \\ \leftrightarrow& S^{\rm P}\phi'_m(\boldsymbol r)\otimes S(\overbrace{0,\dots ,0}^{\alpha-1},1,0,\dots ,0)^{\rm T},~~ \tag {10} \end{alignat} $$ where $S^{\rm P}$ is the translation-deleted part of $S$. One can easily find the resulting state of any symmetry operations on $u^{\boldsymbol k}_{m\alpha} (\boldsymbol r)$. Above all, the direct product between the representations of $V_1$ and $V_2$ is that of $V_0$. When taking the spin into account, we could regard $V$ as $V_1\otimes V_2\otimes V_{\rm spin}$. In this way, the above discussion can be immediately generalized to the situation with spin, where we must consider the effect of $S^{\rm P}$ on the spin part. For different systems, the representations formed by $V_2$ differ from each other, but they are mostly direct sum of pure one-dimensional representations. Thus in the database, we only collect matches between representations of the atomic orbital with or without spin and all one-dimensional representations. Expansion with Hermitian Matrices. Considering the hermiticity of $\mathcal{H}_{\boldsymbol k_0}(\delta\boldsymbol k)$, it can be expanded by $n^2$ linearly independent Hermitian matrices with real coefficients $h(\delta\boldsymbol k)$, and then we can obtain $n^2$ independent equations. We denote these Hermitian matrices by $\{B_i,i=1,2,\dots ,n^2\}$. They span a linear space on real number field. Then $\mathcal{H}_{\boldsymbol k_0}(\delta\boldsymbol k)$ can be expressed as $$ \mathcal{H}_{\boldsymbol k_0}(\delta\boldsymbol k) = \sum_i h_i(\delta\boldsymbol k)B_i.~~ \tag {11} $$ Therefore, Eq. (7) can be expressed as $$\begin{alignat}{1} &\sum_j h_j(\gamma^{-1}\delta\boldsymbol k) D(S^u)B_jD(S^u)^{-1} = \sum_i h_i(\delta\boldsymbol k)B_i,~~~~~~~~ \tag {12a} \end{alignat} $$ $$\begin{alignat}{1} &\sum_j h_j(\gamma^{-1}\delta\boldsymbol k) D(S^a)B^*_jD(S^a)^{-1} = \sum_i h_i(\delta\boldsymbol k)B_i.~~~~~~~~ \tag {12b} \end{alignat} $$ Now we define two new matrices $\mathcal{M}$ and $\mathcal{N}$ as follows: $$\begin{alignat}{1} &D(S^u)B_{j}D(S^u)^{-1} = \sum_m \mathcal{M}_{mj}B_m,~~ \tag {13a} \end{alignat} $$ $$\begin{alignat}{1} &D(S^a)B^*_{j}D(S^a)^{-1} = \sum_m \mathcal{M}_{mj}B_m,~~ \tag {13b} \end{alignat} $$ $$\begin{alignat}{1} &h_j(\gamma^{-1}\delta\boldsymbol k) = \sum_n \mathcal{N}_{nj}h_n(\delta\boldsymbol k).~~ \tag {14} \end{alignat} $$ Comparing the coefficients of $B_i$ on both sides of Eq. (12) yields $\sum_n (\mathcal{MN}^{\rm T})_{\rm in}h_n(\delta\boldsymbol k) = h_i(\delta\boldsymbol k)$ for $i=1,2,\dots ,n^2$, that is, $$ \mathcal{N}^{\rm T}h(\delta\boldsymbol k) = \mathcal{M}^{-1}h(\delta\boldsymbol k).~~ \tag {15} $$ On the one hand, if the basis $\{|u^{\boldsymbol k_0}_i\rangle, i =1,2,\dots ,n\}$ is orthonormal, it is obvious that $D$ represents unitary matrices ($D^{-1}=D^†$), which implies that $D(S^u)B_{j}D(S^u)^{-1}$ are vectors in Hermitian matrix space on real number field. Thus, $\mathcal{M}$ must be real. On the other hand, if we define the inner product of Hermitian matrices as ${\rm tr}(B_iB_j^†)$, we can prove $$\begin{align} &{\rm tr}\left[D(S^u)B_{i}D(S^u)^{-1}[D(S^u)B_{j}D(S^u)^{-1}]^†\right]\\ ={}&{\rm tr}(B_{i}B_{j}^†),~~ \tag {16} \end{align} $$ which means that the inner product remains unchanged (the same for $S^a$). This implies that if $\{B_i\}$ is an orthonormal basis, $\mathcal{M}$ must be orthogonal ($\mathcal{M}^{-1}=\mathcal{M}^{\rm T}$). Especially, if the combination of inversion and time reversal symmetry $\mathcal{PT}$ is in consideration, then we have $$ \mathcal{M}(\mathcal{PT})h(\delta\boldsymbol k) = h(\delta\boldsymbol k),~~ \tag {17} $$ which allows us to simplify the calculation. If we assume that the solution of $[\mathcal{M}(\mathcal{PT})-I] X = 0$ is $\{b_i, i=1,2,\dots ,m\}$, where $m < n^2$ always holds, then $h(\delta\boldsymbol k)=\sum_ih'_i(\delta\boldsymbol k)b_i$ and $\mathcal{H}_{\boldsymbol k_0}(\delta\boldsymbol k) = \sum_ih'_i(\delta\boldsymbol k) B'_i$, where we denote $\sum_j(b_i)_jB_j$ by $B'_i$. Now we can recalculate $\mathcal{M}'$ with $\{B'_i, i=1,2,\dots ,m\}$, but this time only the $m\times m$ submatrix on the top left corner of $\mathcal{M}'$ is meaningful and actually we can only get this submatrix with the incomplete set of Hermitian matrices $\{B'_i\}$. Expansion of $h(\delta\boldsymbol k)$. Since $h(\delta\boldsymbol k)$ is a polynomial of $\delta\boldsymbol k$, we could expand $h(\delta\boldsymbol k)$ by a monomial basis of $\delta\boldsymbol k$. Suppose that $\{f_i(\delta\boldsymbol k), i=1,2,\dots ,d\}$ is the basis, then $h(\delta\boldsymbol k)$ can be expanded as $$ h_i(\delta\boldsymbol k) = \sum_j f_j(\delta\boldsymbol k)A_{ji}.~~ \tag {18} $$ Then with the definition of matrix $\mathcal{F}$ below, $$ f_i(\gamma^{-1}\delta\boldsymbol k) = \sum_j \mathcal{F}_{ji}f_j(\delta\boldsymbol k).~~ \tag {19} $$ Equation (14) gives $\mathcal{F}A=A\mathcal{N}$, and then Eq. (15) gives $A=\mathcal{F}A\mathcal{M}^{\rm T}$. If defining $\alpha$ as a vector from stacking up columns of $A^{\rm T}$, we can obtain $$ (\mathcal{F}\otimes \mathcal{M}-I)\alpha=0.~~ \tag {20} $$ Assuming that the solution is $\alpha = \sum_l c_l\tilde{\alpha}(l)$, where $\{\tilde{\alpha}(l)\}$ is a basis of the solution space and $c_l$ are real coefficients, then $A = \sum_l c_l\tilde{A}(l)$, and finally $$\begin{alignat}{1} {\cal H}_{\boldsymbol k_0}(\delta\boldsymbol k) = \sum_ih_i(\delta\boldsymbol k)B_i = \sum_{ijl}c_lf_j(\delta\boldsymbol k){\tilde{A}(l)}_{ji}B_i. ~~~~~~~ \tag {21} \end{alignat} $$ Note that $\mathcal{F}$ is a block diagonal matrix, since polynomials with different orders cannot transform to each other. We could separately calculate $\mathcal{F}^{(n)}$ and $A^{(n)}$ for each Taylor series expansion order $f^{(n)}(\delta\boldsymbol k)$. Finally, we stack $f^{(n)}(\delta\boldsymbol k)$ and $A^{(n)}$, respectively, to get $f(\delta\boldsymbol k)$ and $A$, that is, $$ f(\delta\boldsymbol k) = \begin{pmatrix} f^{(0)}(\delta\boldsymbol k)\\ f^{(1)}(\delta\boldsymbol k)\\ f^{(2)}(\delta\boldsymbol k)\\ \vdots \end{pmatrix} ,~~ \tag {22} $$ $$\begin{align} &\tilde{A}= \sum_{l=1}^{l^{(0)}} c_{l}\begin{pmatrix} \tilde{A}^{(0)}(l)\\ 0\\ 0\\ \vdots \end{pmatrix} +\sum_{l=l^{(0)}+1}^{l^{(0)}+l^{(1)}} c_{l}\begin{pmatrix} 0\\ \tilde{A}^{(1)}(l)\\ 0\\ \vdots \end{pmatrix}\\ &+\sum_{l=l^{(0)}+l^{(1)}+1}^{l^{(0)}+l^{(1)}+l^{(2)}}c_{l}\begin{pmatrix} 0\\ 0\\ \tilde{A}^{(2)}(l)\\ \vdots \end{pmatrix}+\cdots~~ \tag {23} \end{align} $$ The Case of $\boldsymbol k_0$ on the BZ Boundary under a Symmetry Operation $\gamma $ with $\gamma \boldsymbol k_0-\boldsymbol k_0\neq0$. When the $\boldsymbol k_0$ is on the BZ boundary, we should take into account the following equation: $$ \mathcal{H}_{\boldsymbol k+\boldsymbol G} = e^{-i\boldsymbol G\cdot \boldsymbol r}\mathcal{H}_{\boldsymbol k}e^{i\boldsymbol G\cdot \boldsymbol r},~~ \tag {24} $$ where $\boldsymbol G$ is an arbitrary reciprocal lattice vector. With this equation, the symmetry operation meeting $\gamma \boldsymbol k_0-\boldsymbol k_0=\boldsymbol G\neq0$ can also limit the form of the Hamiltonian around $\boldsymbol k_0$, since $$\begin{alignat}{1} S\mathcal{H}_{\boldsymbol k_0+\delta \boldsymbol k}S^{-1} = \mathcal{H}_{\gamma\boldsymbol k_0+\gamma \boldsymbol \delta k} =e^{-i\boldsymbol G\cdot \boldsymbol r}\mathcal{H}_{\boldsymbol k_0+\gamma \delta \boldsymbol k}e^{i\boldsymbol G\cdot \boldsymbol r}.~~ \tag {25} \end{alignat} $$ We can define $\tilde{S}$ as $\tilde{S} = e^{i\boldsymbol (\gamma\boldsymbol k_0 - \boldsymbol k_0)\cdot \boldsymbol r}S$, which leads to an equation with the same form as Eq. (5), $$ \tilde{S} \mathcal{H}_{\boldsymbol k_0+\delta\boldsymbol k}\tilde{S}^{-1} = \mathcal{H}_{\boldsymbol k_0+\gamma\delta \boldsymbol k}.~~ \tag {26} $$ Although the set of symmetry operations $\{\tilde{S}\}$ actually describes the full symmetry of $\mathcal{H}_{\boldsymbol k_0}$, it is not always a group. The product of two arbitrary elements $\tilde{S}_1$ and $\tilde{S}_2$ is $$ \tilde{S}_1\tilde{S}_2 = e^{i\gamma_1(\boldsymbol k_0-\gamma_2\boldsymbol k_0)\cdot \boldsymbol \tau_1}e^{i(\gamma_1\gamma_2\boldsymbol k_0-\boldsymbol k_0)\cdot \boldsymbol r}S_1S_2,~~ \tag {27} $$ which is still in $\{\tilde{S}\}$ and corresponding to $S_1S_2$ only if $e^{i\gamma_1(\boldsymbol k_0-\gamma_2\boldsymbol k_0)\cdot \boldsymbol \tau_1}=1$. Due to the existence of the symmetry $\gamma$ with $\gamma \boldsymbol k_0\neq\boldsymbol k_0$, only if all symmetries in consideration are symmorphic ($\tau=0$), $\{\tilde{S}\}$ becomes a group. Otherwise, we either drop the symmetry operation $\gamma$ with $\gamma \boldsymbol k_0\neq\boldsymbol k_0$ or the nonsymmorphic symmetry operations in our code. Application to Magnetic Topological Material MnBi$_2$Te$_4$. MnBi$_2$Te$_4$ is a versatile magnetic topological materials to realize to quantum anomalous Hall (QAH) state, antiferromagnetic (AFM) topological insulator, magnetic axion insulator, tunable dynamical axion field and ferromagnetic (FM) Weyl semimetal.[18,21–28] Here, with the above method, we take MnBi$_2$Te$_4$ as an example to demonstrate how to construct the $\boldsymbol{k}$$\cdot$$\boldsymbol{p}$ Hamiltonian.
cpl-38-7-077105-fig2.png
Fig. 2. Crystal structure and band structure of magnetic topological insulator MnBi$_2$Te$_4$. (a) The unit cell of AFM MnBi$_2$Te$_4$ consists of two SLs. The red arrows represent the spin moment of Mn atom. The green arrow denotes for the half translation operator $\tau_{1/2}$. [(b), (c)] The fat band structures of FM and AFM states. The blue and red dots denote the characters of Te and Bi $p_z$-orbitals, respectively. (d) The unit cell of FM MnBi$_2$Te$_4$ has one SL. (e) The process of construct FM MnBi$_2$Te$_4$ $\boldsymbol{k}$$\cdot$$\boldsymbol{p}$ Hamiltonian.
As shown in Fig. 2(a), MnBi$_2$Te$_4$ has a layered crystal structure with a triangle lattice. The trigonal axis (threefold rotation symmetry $C_{3z}$) is defined as the $z$ axis, a binary axis (twofold rotation symmetry $C_{2x}$) is defined as the $x$ axis and a bisectrix axis (in the reflection plane) is defined as the $y$ axis for the coordinate system. The material consists of septuple layers (SL) (e.g., Te–Bi–Te–Mn–Te–Bi–Te) arranged along the $z$ direction. The structural of nonmagnetic MnBi$_2$Te$_4$ is described by the space group $D^5_{3d}$ ($R\bar{3}m$) whose generators are $C_{3z}$, $P$, $C_{2x}$ and pure translations which can be ignored since their representations formed by LCAOs are identities. The first-principles calculations of FM and AFM states [Figs. 2(b) and 2(c)] show that four bands at the $\varGamma$ point near the Fermi-level all contains $p_z$ orbitals of Bi and Te atom, the basis for both FM and AFM states can be expressed as $\{|{\rm Bi}, p_z,\uparrow\rangle^+, i|{\rm Te}, p_z,\uparrow\rangle^-, |{\rm Bi}, p_z,\downarrow\rangle^+, -i|{\rm Te}, p_z,\downarrow\rangle^-\}$, where the superscripts $+/-$ indicate parities. For the FM state, the symmetries are three-fold rotation $C_{3z}$, inversion symmetry $P$, and the combination of two-fold rotation and time reversal symmetry $C_{2x}\varTheta$, the representation matrices can be obtained by our code: $$\begin{align} D(C_{3z})={}&{\exp}\Big(-i\frac{\pi}3\sigma_z\otimes 1_{2\times2}\Big),\\ D(P)={}&1_{2\times2}\otimes\tau_z,\\ D(C_{2x}\varTheta) ={}& {\exp}\Big(i\frac\pi2\sigma_z\otimes 1_{2\times2}\Big)\mathcal{K},~~ \tag {28} \end{align} $$ where $\mathcal{K}$ is the complex conjugation operator. Then the code calculates the Matrices $\mathcal{F}$ and $\mathcal{M}$ for each symmetry automatically, and finally gives the $\boldsymbol k\cdot\boldsymbol p$ Hamiltonian up to the 2nd order as follows, $$\begin{alignat}{1} \mathcal{H}_{\rm FM}(\boldsymbol k)= \begin{pmatrix} M_1(\boldsymbol k) & A_1k_z & 0 & A_2 k_-\\ A_1k_z & M_2(\boldsymbol k) & A_4k_- & 0\\ 0 & A_4k_+ & M_3(\boldsymbol k) & A_3k_z\\ A_2 k_+ & 0 & A_3k_z & M_4(\boldsymbol k) \end{pmatrix},~~~ \tag {29} \end{alignat} $$ where $k_\pm = k_x\pm ik_y$ and $M_n(\boldsymbol k) = M_0^n+B_1^nk_z^2+B_2^n(k_x^2+k_y^2)$ with $n=1,2,3,4$. Note that the off-diagonal two-order terms are omitted, since they contribute higher order of $\boldsymbol k$ to the energy. To clearly show the effect induced by the FM exchange field, we consider the nonmagnetic (NM) state and assume that both $C_{2x}$ and $\varTheta$ are preserved. In addition to $D(C_{3z}),D(P)$, the representation matrices also include $D(C_{2x}) = {\exp}[-i(\pi/2)\sigma_x\otimes 1_{2\times2}]$, $D(\varTheta)=i\sigma_y\otimes 1_{2\times2}\mathcal{K}$. Then we can obtain the NM $\boldsymbol k\cdot\boldsymbol p$ Hamiltonian $$\begin{alignat}{1} &\mathcal{H}_{\rm NM}(\boldsymbol k) \\ ={}& \epsilon_0(\boldsymbol k)+\begin{pmatrix} M(\boldsymbol k) & \tilde{A}_1k_z &0 & \tilde{A}_2k_-\\ \tilde{A}_1k_z & -M(\boldsymbol k) & \tilde{A}_2k_- & 0\\ 0 & \tilde{A}_2k_+ & M(\boldsymbol k) & -\tilde{A}_1k_z\\ \tilde{A}_2 k_+ & 0& -\tilde{A}_1k_z & -M(\boldsymbol k) \end{pmatrix}, ~~~~~~~ \tag {30} \end{alignat} $$ where $\epsilon_0(\boldsymbol k) = C + D_1k_z^2+ D_2(k_x^2+k_y^2)$ and $M(\boldsymbol k) = M_0+B_1 k_z^2+B_2(k_x^2+k_y^2)$. Due to the same symmetry and the orbital basis, the $\mathcal{H}_{\rm NM}(\boldsymbol k)$ of NM MnBi$_2$Te$_4$ is the same as that of topological insulator Bi$_2$Se$_3$.[7] Comparing $\mathcal{H}_{\rm NM}(\boldsymbol k)$ and $\mathcal{H}_{\rm FM}(\boldsymbol k)$, we can obtain the FM term $\delta \mathcal{H}_{\rm FM}(\boldsymbol k)$ induced by the FM order, $$\begin{alignat}{1} &\delta \mathcal{H}_{\rm FM}(\boldsymbol k)\\ ={}&\begin{pmatrix} M_{13}(\boldsymbol k) & A_1^{\prime}k_z &0 &A_2^{\prime} k_-\\ A_1^{\prime}k_z & M_{24}(\boldsymbol k) & -A_2^{\prime} k_-&0\\ 0 & -A_2^{\prime} k_+ & -M_{13}(\boldsymbol k) & A_1^{\prime}k_z\\ A_2^{\prime} k_+ & 0 & A_1^{\prime}k_z & -M_{24}(\boldsymbol k) \end{pmatrix},~~~~~~ \tag {31} \end{alignat} $$ where $M_{13}(\boldsymbol k) = [M_1(\boldsymbol k)- M_3(\boldsymbol k)]/2$, $M_{24}(\boldsymbol k) = [M_2(\boldsymbol k)- M_4(\boldsymbol k)]/2$, $A_1^{\prime}=(A_1+A_3)/2$ and $A_2^{\prime}=(A_2-A_4)/2$. The combination of NM $\mathcal{H}_{\rm NM}(\boldsymbol k)$ and FM $\delta\mathcal{H}_{\rm FM}(\boldsymbol k)$ can give a clear physical picture of the FM Weyl semimetal MnBi$_2$Te$_4$, schematically shown in Fig. 3. The mass term $M_{13,24}(\boldsymbol k)$, due to the FM exchange field, lifts the band degeneracy to produce a pair of Weyl points along the $k_z$ direction.
cpl-38-7-077105-fig3.png
Fig. 3. Schematic of FM Weyl semimetal from NM insulator. Due to the time reversal symmetry and inversion symmetry, the bands are double degeneracy for the NM state, described by $\mathcal{H}_{\rm NM}(\boldsymbol k)$, schematically shown on the left. Through considering the FM exchange with $\delta\mathcal{H}_{\rm FM}(\boldsymbol k)$, the band degeneracy is lifted by the mass term $M_{13,24}$ to produce a pair of Weyl points.
For the A-type AFM state, the time reversal symmetry $\varTheta$ is broken, but the combination operation $\varTheta\tau_{1/2}$ of $\varTheta$ and a translation $\tau_{1/2}$ [shown in Fig. 2(a)] is preserved. Compared with nonmagnetic state, the representation of the combination symmetry $\varTheta\tau_{1/2}$ in AFM state is the same as that of the time reversal symmetry $\varTheta$. Thus their Hamiltonians are the same in form. In summary, we have developed a programmable algorithm to efficiently construct the one-/two-/three-dimensional effective $\boldsymbol{k}$$\cdot$$\boldsymbol{p}$ Hamiltonians for condensed matters, which avoids time-consuming calculations and possible latent mistakes. The corresponding code is open source and it can be downloaded through the github platform. Our method has been successfully applied to the recent works on antiferromagnetic MnBi$_2$Te$_4$-family axion insulators.[18–20] Therefore, our algorithm and open-source code provide a useful tool for theoretical studies of condensed matters. References Colloquium : Topological insulatorsTopological insulators and superconductorsWeyl and Dirac semimetals in three-dimensional solidsTopological insulators with inversion symmetryQuantum Spin Hall Effect and Topological Phase Transition in HgTe Quantum WellsTopological insulators from the perspective of first-principles calculationsTopological insulators in Bi2Se3, Bi2Te3 and Sb2Te3 with a single Dirac cone on the surfaceMotion of Electrons and Holes in Perturbed Periodic FieldsSemiconductors and SemimetalsModel Hamiltonian for topological insulatorsHexagonal Warping Effects in the Surface States of the Topological Insulator Bi 2 Te 3 Chern Semimetal and the Quantized Anomalous Hall Effect in HgCr 2 Se 4 Qsymm: algorithmic symmetry finding and symmetric Hamiltonian generationTopological Axion States in the Magnetic Insulator MnBi 2 Te 4 with the Quantized Magnetoelectric EffectLarge Dynamical Axion Field in Topological Antiferromagnetic Insulator Mn 2 Bi 2 Te 5Dynamical axion state with hidden pseudospin Chern numbers in MnBi 2 Te 4 -based heterostructuresIntrinsic magnetic topological insulators in van der Waals layered MnBi 2 Te 4 -family materialsExperimental Realization of an Intrinsic Magnetic Topological Insulator *Prediction and observation of an antiferromagnetic topological insulatorUnique Thickness-Dependent Properties of the van der Waals Interlayer Antiferromagnet MnBi 2 Te 4 FilmsQuantum anomalous Hall effect in intrinsic magnetic topological insulator MnBi 2 Te 4Robust axion insulator and Chern insulator phases in a two-dimensional antiferromagnetic topological insulatorIntrinsic magnetic topological insulator phases in the Sb doped MnBi2Te4 bulks and thin flakesTunable 3D/2D magnetism in the (MnBi2Te4)(Bi2Te3)m topological insulators family
[1] Hasan M Z and Kane C L 2010 Rev. Mod. Phys. 82 3045
[2] Qi X L and Zhang S C 2011 Rev. Mod. Phys. 83 1057
[3] Armitage N P, Mele E J, and Vishwanath A 2018 Rev. Mod. Phys. 90 015001
[4] Fu L and Kane C L 2007 Phys. Rev. B 76 045302
[5] Bernevig B A, Hughes T L, and Zhang S C 2006 Science 314 1757
[6] Zhang H and Zhang S C 2013 Phys. Status Solidi RRL 7 72
[7] Zhang H, Liu C X, Qi X L, Dai X, Fang Z, and Zhang S C 2009 Nat. Phys. 5 438
[8] Luttinger J M and Kohn W 1955 Phys. Rev. 97 869
[9] Kane E O 1966 The ${\boldsymbol k}\cdot{\boldsymbol p}$ Method (Amsterdam: Elsevier) chap 3 p 75
[10]Voon L C L Y and Willatzen M 2009 The ${\boldsymbol k}$$\cdot$${\boldsymbol p}$ Method: Electronic Properties of Semiconductors (Berlin: Springer Science & Business Media)
[11] Liu C X, Qi X L, Zhang H, Dai X, Fang Z, and Zhang S C 2010 Phys. Rev. B 82 045122
[12] Fu L 2009 Phys. Rev. Lett. 103 266801
[13] Xu G, Weng H, Wang Z, Dai X, and Fang Z 2011 Phys. Rev. Lett. 107 186806
[14]Gresch D 2018 Identifying Topological Semimetals (PhD Dissertation) (Troyer, Matthias, ETH Zürich)
[15]Gresch D kdotp-Symmetry Code
[16] Varjas D et al. 2018 New J. Phys. 20 093026
[17]Akhmerov A Qsymm Code
[18] Zhang D, Shi M, Zhu T, Xing D, Zhang H, and Wang J 2019 Phys. Rev. Lett. 122 206401
[19] Zhang J, Wang D, Shi M, Zhu T, Zhang H, and Wang J 2020 Chin. Phys. Lett. 37 077304
[20] Wang H, Wang D, Yang Z, Shi M, Ruan J, Xing D, Wang J, and Zhang H 2020 Phys. Rev. B 101 081109
[21] Li J, Li Y, Du S, Wang Z, Gu B L, Zhang S C, He K, Duan W, and Xu Y 2019 Sci. Adv. 5 eaaw5685
[22] Gong Y, Guo J, Li J, Zhu K, Liao M, Liu X, Zhang Q, Gu L, Tang L, Feng X, Zhang D, Li W, Song C, Wang L, Yu P, Chen X, Wang Y, Yao H, Duan W, Xu Y, Zhang S C, Ma X, Xue Q K, and He K 2019 Chin. Phys. Lett. 36 076801
[23] Otrokov M M, Klimovskikh I I, Bentmann H, Estyunin D, and Chulkov E V 2019 Nature 576 416
[24] Otrokov M M, Rusinov I P, Blanco-Rey M, Hoffmann M, Vyazovskaya A Y, Eremeev S V, Ernst A, Echenique P M, Arnau A, and Chulkov E V 2019 Phys. Rev. Lett. 122 107202
[25] Deng Y, Yu Y, Shi M Z, Guo Z, Xu Z, Wang J, Chen X H, and Zhang Y 2020 Science 367 895
[26] Liu C, Wang Y, Li H, Wu Y, Li Y, Li J, He K, Xu Y, Zhang J, and Wang Y 2020 Nat. Mater. 19 522
[27] Chen B, Fei F, Zhang D, Zhang B, Liu W, Zhang S, Wang P, Wei B, Zhang Y, and Zuo Z 2019 Nat. Commun. 10 4469
[28] Klimovskikh I I, Otrokov M M, Estyunin D, Eremeev S V, and Chulkov E V 2020 npj Quantum Mater. 5 54