This note is heavily based on Yoichi Ando’s review paper: https://arxiv.org/pdf/1304.5693.pdf # Abstract In this work, we aim to give a pedagogical introduction to the $Z_2$ [[Z2 Topological Invariant|topological invariants]] of [[Time-Reversal Symmetric Topological Insulator|time-reversal symmetric topological insulators]] (TRS-TIs). We present a very detailed derivation of the theoretical framework, as well as discuss some of the possible methods to calculate the [[Topological Invariant|invariants]] starting from *ab initio* or [[Tight-Binding Model|tight-binding]] calculations. This work is not intended as a literature review of [[Topological Insulator|topological insulators]], rather, as a solid introduction to the underlying theory and the available computational methods. # Introduction Since the seminal theory papers on [[Topological Insulator|topological insulators]] (TIs) in 2007, both the experimental and theoretical sides of the field have rapidly grown. The [[topology]] in the name comes from the [[Wavefunction|wavefunctions]] of TIs: they span a topologically nontrivial [[Hilbert Space]]. Since these systems possess a bulk [[Band Gap|gap]], and [[topology]] cannot change without closing a gap, then TIs must have [[gapless]] [[Surface State|surface states]] at the interface with [[vacuum]] (or any other [[trivial insulator]]). In 3D [[Strong Topological Insulator|strong TIs]], this manifests as [[gapless]] [[Surface State|surface states]], and in 2D, [[gapless]] [[Edge State|edge states]]. This is known as the [[bulk-boundary correspondence]]. In [[Strong Topological Insulator|strong TIs]], the [[gapless]] [[surface state]] is protected by [[Time-Reversal Operator|time-reversal symmetry]]. While most research is focused on systems with both [[Time-Reversal Operator|time-reversal]] and [[Inversion Operator|inversion symmetries]], [[Inversion Operator|inversion symmetry]] is not a necessary condition to have a [[Topological Insulator|TI]]. It is also possible to have topological insulators that are protected by [[Spatial Symmetry|spatial symmetries]] (e.g. [[Mirror Symmetry]]). These are known as [[Topological Crystalline Insulator|topological crystalline insulators]] (TCIs) and are characterized by higher order invariants and are thus beyond the scope of this work. It is also worth noting that most current research on 3D TIs is concerned with [[Strong Topological Insulator|strong TIs]]. However, there are a few papers on [[Weak Topological Insulator|weak TIs]], but they are much harder to study experimentally. In general, TIs are characterized by [[band inversion]] and they commonly have strong [[spin-orbit coupling]] (SOC), a necessary condition for [[band inversion]]. The strong SOC leads to what's known as helical spin polarization or [[spin-momentum locking]], where the spin is locked to the surface momentum. This paper is organized as follows: in section II we discuss some preliminary topics required for a full understanding of the theory, including [[Berry phase]], [[Time-Reversal Operator|time reversal symmetry]] and [[Inversion Operator|inversion symmetry]]. In section III, we cover the 1D time-periodic model introduced by Fu and Kane. We consider the simple case with only two bands and then extend it to $2N$ bands. We additionally discuss the extension to 2D systems. In section IV, we extend the discussion to 3D systems and note the different types of [[Time-Reversal Symmetric Topological Insulator|TIs]] in 3 spatial dimensions. In Section V we treat systems with [[Inversion Operator|inversion symmetry]] in full detail, and derive the simplified [[Fu-Kane formula]]. In Section VI, we discuss some of the available computational methods, including the computation of [[Inversion Operator|parity]] eigenvalues for [[Centrosymmetric]] systems, and the [[Wannier Charge Center]] (WCC) method for general systems. We derive some of the equations for the WCC method and discuss its strengths and weaknesses. For both methods, we note some of the available open-source software packages. We note that this work aims to provide an educational introduction to the $Z_2$ [[Z2 Topological Invariant|topological invariants]] that characterize [[Time-Reversal Symmetric Topological Insulator|TIs]], with an emphasis on how the equations are derived and how to calculate the invariants computationally. Thus, we will only very briefly discuss the properties of these systems, and direct readers to the many excellent papers available on the experimental side of the field and other properties of TIs. # Preliminaries ## Berry Phase The [[Berry phase]] is especially important for studying [[topological systems]], and we show some important relations here. Consider a system that depends on some time-dependent [[vector]] $\pmb{x}(t)$. We assume the [[Hamiltonian]] is $H[\pmb{x}(t)]$, and the Hamiltonian's $n^{\text{th}}$ [[eigenstate]] is $\ket{n, \pmb{x}(t)}$. The [[Schrödinger equation]] for this system is: $ \begin{align} H[\pmb{x}(t)] \ket{n, \pmb{x}(t)} = E_n[\pmb{x}(t)] \ket{n, \pmb{x}(t)} .\end{align} $ We assume $\pmb{x}(t=0) = \pmb{x}_0$ and that $\pmb{x}(t)$ changes [[Adiabatic|adiabatically]] in time. From the time-dependent [[Schrödinger equation]]: $ \begin{align} H[\pmb{x}(t)] \ket{n, t} = i\hbar \frac{\partial}{\partial t} \ket{n, t}, \end{align} $ we can obtain the [[eigenstate]] at time $t$: $ \begin{align} \ket{n,t} = &\exp\left(\frac{i}{\hbar}\int_0^t \text{d} t' E_n[\pmb{x}(t')]\right) \exp(-\gamma[\pmb{x}(t)]) \ket{n, \pmb{R}(t)} , \end{align} $ where $ \begin{align} \gamma[\pmb{x}(t)] \equiv \int_0^t \text{d} t'\left(i \pmb{\dot{x}}(t') \cdot \bra{n, \pmb{x}(t')} \nabla_{\pmb{x}} \ket{n, \pmb{x}(t')}\right). \end{align} $ The phase of the first term in the equation for $\ket{n, t}$ is the trivial dynamical phase, and $\gamma$ is known as the nontrivial geometric or [[Berry phase]]. Now assume that the [[vector]] $\pmb{x}(t)$ traces a closed [[contour]] $C$ from $t = 0$ to $t=T$, i.e. $\pmb{x}(t=0) = \pmb{x}(t=T) = \pmb{x}_0$. We label the area enclosed by the contour $S$. Then we can simplify $\gamma$ as follows: $ \begin{align} \gamma[C] &\equiv \int_0^T \text{d} t\left(i \pmb{\dot{x}}(t) \cdot \bra{n, \pmb{x}(t)} \nabla_{\pmb{x}} \ket{n, \pmb{x}(t)}\right), \\ &= \oint_C \text{d} \pmb{x} \cdot i \bra{n, \pmb{x}(t)} \nabla_{\pmb{x}} \ket{n, \pmb{x}(t)},\\ &\equiv -\oint_C \text{d} \pmb{x} \cdot \pmb{A}_n(\pmb{x}),\\ &\equiv -\int_S \text{d} \pmb{S} \cdot \pmb{F}_n(\pmb{x}), \end{align} $ where in the last equation we used [[Stokes' theorem]]. Here we have defined the [[Berry connection]] $\pmb{A}_n(\pmb{x})$ and the [[Berry curvature]] $\pmb{F}_n(\pmb{x})$ of the $n^{\text{th}}$ [[Band Structure|band]]: $ \begin{align} \pmb{A}_n(\pmb{x}) &\equiv -i \bra{n,\pmb{x}} \nabla_{\pmb{x}} \ket{n,\pmb{x}},\\ \pmb{F}_n(\pmb{x}) &\equiv \nabla_{\pmb{x}} \times \pmb{A}_n(\pmb{x}) .\end{align} $ ## Time-Reversal Operator Consider a system preserving [[Time-Reversal Operator|time-reversal symmetry]] (TRS). The [[time-reversal operator]] for a spin-1/2 system is $\Theta = -i \sigma_y K$ where $K$ denotes complex conjugation and $\sigma_\mu$ ($\mu = 1,2,3$) are the [[Pauli matrices]]. Note that $\Theta$ is an [[anti-unitary operator]], i.e. $\Theta^{-1} = \Theta^{\dagger}$. For half-integer spin, we also have $\Theta^2 = -1$. ## The Bloch Hamiltonian Consider a periodic system with Hamiltonian $\mathcal{H}$: $ \begin{align} \mathcal{H} \ket{\psi_{n\pmb{k}}} = E_{n\pmb{k}} \ket{\psi_{n\pmb{k}}} .\end{align} $ Using [[Bloch's theorem]] we can write: $ \begin{align} \ket{\psi_{n\pmb{k}}} = e^{i \pmb{k}\cdot\pmb{r}} \ket{u_{n\pmb{k}}} .\end{align} $ which leads to the reduced [[Schrödinger equation]] for $\ket{u_{n\pmb{k}}}$: $ \begin{align} H(\pmb{k}) \ket{u_{n\pmb{k}}} = E_{n\pmb{k}} \ket{u_{n\pmb{k}}} ,\end{align} $ where $H(\pmb{k}) = e^{-i \pmb{k}\cdot\pmb{r}} \mathcal{H} e^{+i \pmb{k}\cdot\pmb{r}}$ is the [[Bloch's Theorem|Bloch Hamiltonian]]. If the system preserves TRS, i.e. $[\Theta, \mathcal{H}] = 0$, then: $ \begin{align} \Theta \mathcal{H} &= \mathcal{H} \Theta, \\ \Theta \mathcal{H} \Theta^{-1} &= \mathcal{H}, \\ e^{+i \pmb{k}\cdot\pmb{r}} \Theta \mathcal{H} \Theta^{-1} e^{-i \pmb{k}\cdot\pmb{r}} &= e^{+i \pmb{k}\cdot\pmb{r}} \mathcal{H} e^{-i \pmb{k}\cdot\pmb{r}}, \\ \Theta e^{-i \pmb{k}\cdot\pmb{r}} \mathcal{H} e^{+i \pmb{k}\cdot\pmb{r}} \Theta^{-1} &= H(-\pmb{k}), \\ \Theta H(\pmb{k}) \Theta^{-1} &= H(-\pmb{k}) .\end{align} $ Thus, the energy states of a system with [[Time-Reversal Operator|time-reversal symmetry]] come in pairs ([[Kramers Theorem|Kramers' pairs]]), where the $+\pmb{k}$ and the $-\pmb{k}$ have the same energy. At the [[time-reversal invariant momenta]] (TRIM), i.e. the points where $+\pmb{k}$ and $-\pmb{k}$ are equivalent due to the periodicity of the [[Brillouin Zone|BZ]], the [[Kramers Theorem|Kramers' pairs]] is [[Degeneracy|degenerate]]. There are eight [[Time-Reversal Invariant Momenta|TRIM]] in 3D and four in 2D, as shown in Fig. \ref{fig:TRIM} ![[Pasted image 20210617162026.png]] ![[Pasted image 20210617162842.png]] # TR Operator in Bloch Basis We can write down a matrix in the [[Bloch's Theorem|Bloch basis]] that relates states at $-\pmb{k}$ and $\pmb{k}$: $ \begin{align} w_{mn}(\pmb{k}) = \bra{u_{m,-\pmb{k}}}\Theta\ket{u_{n,\pmb{k}}} .\end{align} $ This is known as the [[sewing matrix]]. We can compute $\ket{u_{n,-\pmb{k}}}$ as: $ \begin{align} \sum_n w^*_{mn}(\pmb{k})\Theta \ket{u_{n,\pmb{k}}} = \ket{u_{m,-\pmb{k}}} ,\end{align} $ $w(\pmb{k})$ also has the following property: $ \begin{align} w_{mn}(\pmb{k}) &= - w_{nm}(-\pmb{k}) .\end{align} $ It is also easy to show that $w(\pmb{k})$ is a [[Unitary Operator|unitary matrix]]: $ \begin{align} \left(w^{\dagger}(\pmb{k}) w(\pmb{k})\right)_{mn} &= \sum_{l} w^{\dagger}_{ml}(\pmb{k}) w_{ln}(\pmb{k}) \\ &=\sum_{l} \bra{u_{l,\pmb{k}}}\Theta^{\dagger}\ket{u_{m,-\pmb{k}}}\times\bra{u_{l,-\pmb{k}}}\Theta\ket{u_{n,\pmb{k}}},\\ &= \delta_{mn} .\end{align} $ Thus $w^\dagger(\pmb{k}) w(\pmb{k}) = I$ and $w(\pmb{k})$ is unitary. ## U(2N) Berry Connection Matrix For a system with $2N$ bands, the $U(2N)$ [[Berry connection]] matrices are defined as: $ \begin{align} \pmb{a}_{mn}(\pmb{k}) = - i \bra{u_{m,\pmb{k}}}\nabla_{\pmb{k}}\ket{u_{n,\pmb{k}}} .\end{align} $ We can relate $\pmb{a}(-\pmb{k})$ and $\pmb{a}(\pmb{k})$ as follows: $ \begin{align} \pmb{a}_{mn}(-\pmb{k}) &= - i \bra{u_{m,-\pmb{k}}}\nabla_{-\pmb{k}}\ket{u_{n,-\pmb{k}}} ,\\ &= i \bra{u_{m,-\pmb{k}}}\nabla_{\pmb{k}}\ket{u_{n,-\pmb{k}}} , \\ &= i \sum_{l,l'} \bra{u_{l',\pmb{k}}}\Theta^\dagger (w_{ml'}^*(\pmb{k}))^\dagger \nabla_{\pmb{k}} \times(w_{nl}^*(\pmb{k}))\Theta\ket{u_{l,-\pmb{k}}} , \\ &= -i \sum_{l,l'} \bra{u_{l',\pmb{k}}}\Theta w_{l'm}(\pmb{k}) \times ([\nabla_{\pmb{k}} w_{nl}^*(\pmb{k})]\Theta\ket{u_{l,-\pmb{k}}} +w_{nl}^*(\pmb{k})\nabla_{\pmb{k}}\Theta\ket{u_{l,-\pmb{k}}}) ,\\ %&= -i \sum_{l,l'} \bra{u_{l',\pmb{k}}}\Theta w_{l'm}(\pmb{k}) [\nabla_{\pmb{k}} w_{nl}^*(\pmb{k})]\Theta\ket{u_{l,-\pmb{k}}} + \bra{u_{l',\pmb{k}}}\Theta w_{l'm}(\pmb{k}) w_{nl}^*(\pmb{k})\nabla_{\pmb{k}}\Theta\ket{u_{l,-\pmb{k}}} , \\ \end{align} $ $ \begin{align} \pmb{a}_{mn}(-\pmb{k})&= -i \sum_{l,l'} w_{l'm}(\pmb{k}) \bra{u_{l',\pmb{k}}}\Theta^2\ket{u_{l,-\pmb{k}}} \times \nabla_{\pmb{k}} (w^\dagger(\pmb{k}))_{ln} \\ &+ \bra{u_{l',\pmb{k}}}\Theta w_{l'm}(\pmb{k}) w_{nl}^*(\pmb{k})\nabla_{\pmb{k}}\Theta\ket{u_{l,-\pmb{k}}} , %&= i \sum_l w_{lm}(\pmb{k})\nabla_{\pmb{k}} (w^\dagger(\pmb{k}))_{ln} \end{align} $ Thus, we get: $ \begin{align} \pmb{a}(-\pmb{k}) = w(\pmb{k}) \pmb{a}^*(\pmb{k}) w^\dagger(\pmb{k}) + i w(\pmb{k}) \nabla_{\pmb{k}} w^\dagger(\pmb{k}) \end{align} $ Observe that: $ \begin{align} \pmb{a}^*_{mn}(\pmb{k}) &= (- i \bra{u_{m,\pmb{k}}}\nabla_{\pmb{k}}\ket{u_{n,\pmb{k}}} )^* \\ %&= i \bra{u_{n,\pmb{k}}}\nabla_{\pmb{k}}\ket{u_{m,\pmb{k}}}\\ &= \pmb{a}_{nm} .\end{align} $ And thus: $ \begin{align} \text{tr}[\pmb{a^*}(\pmb{k})] = \sum_n \pmb{a^*}_{nn}(\pmb{k}) = \sum_n \pmb{a}_{nn}(\pmb{k}) = \text{tr}[\pmb{a}(\pmb{k})] .\end{align} $ Moreover: $ \begin{align} \nabla_{\pmb{k}} (w w^\dagger) &= \nabla_{\pmb{k}}(1) = 0 ,\\ &= (\nabla_{\pmb{k}} w) w^\dagger + w (\nabla_{\pmb{k}} w^\dagger) = 0 \end{align} $ Thus: $ \begin{align} w \nabla_{\pmb{k}} w^\dagger = -(\nabla_{\pmb{k}} w) w^\dagger .\end{align} $ Taking the [[trace]] of the equation for $\pmb{a}(-\pmb{k})$ above and using the invariance of the trace under [[Similarity Transformation|unitary transformations]]: $ \begin{align} \text{tr}[\pmb{a}(-\pmb{k})] &=\text{tr}[w(\pmb{k}) \pmb{a}^*(\pmb{k}) w^\dagger(\pmb{k})]\\ &+ i \text{tr}[w(\pmb{k}) \nabla_{\pmb{k}} w^\dagger(\pmb{k})]\nonumber ,\\ &=\text{tr}[\pmb{a}^*(\pmb{k})] + i \text{tr}[w(\pmb{k}) \nabla_{\pmb{k}} w^\dagger(\pmb{k})] .\end{align} $ Now replacing $\pmb{k} \to -\pmb{k}$ and using what we have derived so far, the cyclic property of the [[trace]], as well as $\text{tr}[AB] = \text{tr}[(AB)^T] = \text{tr}[B^T A^T] = \text{tr}[A^T B^T]$, we get: $ \begin{align} \text{tr}[\pmb{a}(\pmb{k})] &=\text{tr}[\pmb{a}^*(-\pmb{k})] + i \text{tr}[w(-\pmb{k}) \nabla_{-\pmb{k}} w^\dagger(-\pmb{k})] ,\\ &=\text{tr}[\pmb{a}(-\pmb{k})] - i \text{tr}[w^T(\pmb{k}) \nabla_{\pmb{k}} (w^\dagger(\pmb{k}))^T] ,\\ &=\text{tr}[\pmb{a}(-\pmb{k})] - i \text{tr}[w(\pmb{k}) \nabla_{\pmb{k}} (w^\dagger(\pmb{k}))] ,\\ &=\text{tr}[\pmb{a}(-\pmb{k})] + i \text{tr}[(\nabla_{\pmb{k}} w(\pmb{k})) w^\dagger(\pmb{k})] ,\\ &=\text{tr}[\pmb{a}(-\pmb{k})] + i \text{tr}[w^\dagger(\pmb{k}) \nabla_{\pmb{k}} w(\pmb{k})] .\end{align} $ ## Parity Operator Consider a state $\ket{\pmb{x}, \sigma}$ where $\pmb{x}$ is the spatial coordinate and $\sigma$ is the spin. The parity or spatial [[inversion operator]] $\Pi$ acts on this state as: $ \begin{align} \Pi \ket{\pmb{x}, \sigma} = \ket{-\pmb{x}, \sigma} .\end{align} $ In $\pmb{k}$ space, $\Pi$ acts on a state $\ket{\pmb{k}, \sigma}$ as: $ \begin{align} \Pi \ket{\pmb{k}, \sigma} = \ket{-\pmb{k}, \sigma} .\end{align} $ One can easily show that $\Pi$ is an [[idempotent operator]], i.e. $\Pi^2 = 1$ and the [[Eigenvalue|eigenvalues]] of $\Pi$ are $\pm 1$. Note that $\Pi$ is also [[unitary operator]], i.e. $\Pi^\dagger = \Pi^{-1}$ If we write down a [[Hamiltonian]] as: $ \begin{align} H = \sum_{\pmb{k},\sigma,\sigma'} \ket{\pmb{k},\sigma} H_{\sigma\sigma'} \bra{\pmb{k},\sigma'} .\end{align} $ Then parity acts as: $ \begin{align} \Pi H \Pi^{-1} &= \sum_{\pmb{k},\sigma,\sigma'} \Pi \ket{\pmb{k},\sigma} H(\pmb{k})_{\sigma\sigma'} \bra{\pmb{k},\sigma'}\Pi^\dagger,\\ &= \sum_{\pmb{k},\sigma,\sigma'} \ket{-\pmb{k},\sigma} H(\pmb{k})_{\sigma\sigma'} \bra{-\pmb{k},\sigma'},\\ &= \sum_{\pmb{k},\sigma,\sigma'} \ket{\pmb{k},\sigma} H(-\pmb{k})_{\sigma\sigma'} \bra{\pmb{k},\sigma'} .\end{align} $ Thus, if the system satisfies $H(\pmb{k}) = H(-\pmb{k})$ then $\Pi H \Pi^{-1} = H$ and the system preserves [[Inversion Operator|inversion symmetry]]. # 1D, time-periodic model In this section, we consider a one-dimensional system with [[Time-Reversal Operator|time-reversal symmetry]] and [[lattice constant]] $a = 1$. We additionally assume that the system's parameters are periodically time-dependent, with a period $T$, i.e.: $ \begin{align} H(t+T) = H(t) .\end{align} $ We additionally assume that: $ \begin{align} H(-t) = \Theta H(t) \Theta^{-1} .\end{align} $ Note that for this system, the [[time-reversal invariant momenta]] are at $k = 0, \pi/2$ and $t = 0, T/2$ (since $H(-T/2) = H(-T/2+T) = H(T/2$, and thus $H(T/2) = \Theta H(-T/2) \Theta^{-1} = \Theta H(T/2) \Theta^{-1}$. Refer to Fig. \ref{fig:kramers} for an example of Kramers pairs in a 1D system. ![[Kramers.png]] ## 2-band model We will consider a system with only two [[Band Structure|bands]], a Kramers pair, with [[Bloch's Theorem|Bloch wavefunctions]] $\ket{u_{1k}(t)}$ and $\ket{u_{2k}(t)}$. In the modern theory of polarization, the [[charge polarization]] $P_\rho$ is defined in terms of the [[Berry connection]] as follows: $ \begin{align} P_\rho = \int_{\text{BZ}} \frac{\text{d}\pmb{k}}{(2\pi)^3} \sum_{n=1}^{n_b} A_{n}(\pmb{k}) .\end{align} $ where $A_{n}(\pmb{k}) = -i \bra{u_{n\pmb{k}}} \nabla_{\pmb{k}} \ket{u_{n\pmb{k}}}$ is the [[Berry connection]] of the $n^\text{th}$ band and the sum is over the $n_b$ occupied bands. The integral is done over the [[Brillouin zone]]. For our 1D, 2-band system, this reduces to: $ \begin{align} P_\rho(t) &= \int_{-\pi}^{\pi} \frac{\text{d} k}{2\pi} (-i \bra{u_{1k}(t)} \nabla_k \ket{u_{1k}(t)} \\ &\phantom{=\int_{-\pi}^{\pi} \frac{\text{d} k}{2\pi}}-i \bra{u_{2k}(t)} \nabla_k \ket{u_{2k}(t)}), \nonumber\\ &= \int_{-\pi}^{\pi} \frac{\text{d} k}{2\pi} (a_{11}(k,t) + a_{22}(k,t)), \\ &= \int_{-\pi}^{\pi} \frac{\text{d} k}{2\pi} \text{tr}[a(k,t)], .\end{align} $ We can write $P_\rho(t)$ as $P_\rho(t) = P_1(t) + P_2(t)$ where the partial polarization $P_i(t)$ for $i = 1,2$ is defined as: $ \begin{align} P_i(t) &= \int_{-\pi}^{\pi} \frac{\text{d} k}{2\pi} a_{ii}(k,t) .\end{align} $ The [[time-reversal polarization]] $P_\theta$ is then defined as: $ \begin{align} P_\theta(t) \equiv P_1(t) - P_2(t) = 2P_1(t) - P_\rho(t) .\end{align} $ This is a measure of the charge polarization between the two bands. Since the two bands are a [[Kramers Theorem|Kramers pair]], $P_\theta(t)$ then measures the difference in charge polarization between the spin-up and spin-down bands. Since the system is [[Time-Reversal Operator|time-reversal symmetric]] at $t = 0$ and $t = T/2$ as shown above, then it must be [[Kramers Theorem|Kramers’ degenerate]] at all $k$ points at these two times. Thus, at $t=0$ and $t=T/2$, we must have that $\ket{u_{1,-k}}$ is equal to $\Theta \ket{u_{2k}}$ up to a phase, i.e., $ \begin{align} \Theta \ket{u_{2k}(t)} = e^{-i \chi(k)} \ket{u_{1,-k}(t)}, \quad t = 0, T/2 .\end{align} $ We can also relate $\Theta \ket{u_{1,k}}$ to $\ket{u_{2,-k}}$ as follows: $ \begin{align} \Theta \ket{u_{2k}} &= e^{-i \chi(k)} \ket{u_{1,-k}}, \\ \Theta^2 \ket{u_{2k}} &= \Theta e^{-i \chi(k)} \ket{u_{1,-k}}, \\ -\ket{u_{2k}} &= e^{+i \chi(-k)} \Theta \ket{u_{1,-k}}, \\ \Theta \ket{u_{1k}} &= -e^{-i \chi(-k)} \ket{u_{2,-k}}, \quad t = 0, T/2 .\end{align} $ Thus, we can compute the $w$ matrix for this system at $t = 0$ and $t = T/2$: $ \begin{align} w_{11}(k,t) &= \bra{u_{1,-k}} \Theta \ket{u_{1,k}} \nonumber\\ &= \bra{u_{1,-k}} -e^{-i \chi(-k)} \ket{u_{2,-k}}\nonumber\\ &= 0 = w_{22}(k,t), \\ w_{12}(k,t) &= \bra{u_{1,-k}} \Theta \ket{u_{2,k}} \nonumber \\ &= \bra{u_{1,-k}} e^{-i \chi(k)} \ket{u_{1,-k}} = e^{-i \chi(k)}, \\ w_{21}(k,t) &= \bra{u_{2,-k}} \Theta \ket{u_{1,k}} \nonumber \\ & = \bra{u_{1,-k}} -e^{-i \chi(-k)} \ket{u_{2,-k}} = -e^{-i \chi(-k)}, \end{align} $ And finally: $ \begin{align} w(k) = \begin{pmatrix} 0 & e^{-i \chi(k)} \\ -e^{-i \chi(-k)} & 0\end{pmatrix} , \quad t = 0, T/2 .\end{align} $ Noting that $ \begin{align} a_{11}(-k) &= -i \bra{u_{1,-k}} \frac{\partial}{\partial (-k)} \ket{u_{1,-k}} ,\\ &= i \bra{u_{2,k}}\Theta^\dagger e^{-i\chi(k)} \frac{\partial}{\partial k} e^{+i\chi(k)} \Theta \ket{u_{2,k}},\\ &= i \bra{u_{2,k}}\Theta^\dagger e^{-i\chi(k)} ( e^{+i\chi(k)} \times \left[i \frac{\partial \chi(k)}{\partial k} + \frac{\partial}{\partial k}\right]\Theta \ket{u_{2,k}} ),\\ &= \bra{u_{2,k}}\Theta^\dagger \frac{\partial \chi(k)}{\partial k} \Theta \ket{u_{2,k}} -i \bra{u_{2,k}}\Theta^\dagger\frac{\partial}{\partial k}\Theta \ket{u_{2,k}} ,\\ %&= \frac{\partial \chi(k)}{\partial k} \bra{u_{2,k}}\Theta^\dagger \Theta \ket{u_{2,k}} \nonumber \\ %&\phantom{=}-i \bra{u_{2,k}}\Theta^\dagger \frac{\partial}{\partial k}\Theta \ket{u_{2,k}} ,\\ &= a_{22}(k) - \frac{\partial}{\partial k} \chi(k) .\end{align} $ Similarly: $ \begin{align} a_{22}(-k) &= a_{11}(k) - \frac{\partial}{\partial k} \chi(-k) .\end{align} $ Thus, we can write $P_1(t)$ at $t=0, T/2$ as: $ \begin{align} P_1(t) &= \frac{1}{2\pi} \int_{-\pi}^{\pi} a_{11}(k,t) \\ &= \frac{1}{2\pi} \left( \int_{0}^{\pi} dk a_{11}(k,t) + \int_{-\pi}^{0} dk a_{11}(k,t) \right) ,\\ &= \frac{1}{2\pi} \int_{0}^{\pi} dk \left(a_{11}(k,t) + a_{11}(-k,t)\right) ,\\ &= \frac{1}{2\pi} \int_{0}^{\pi} dk \bigg(a_{11}(k,t) + a_{22}(k,t) \nonumber \\ &\phantom{\frac{1}{2\pi} \int_{0}^{\pi} dk} - \frac{\partial \chi(k)}{\partial k}\bigg),\\ &= \frac{1}{2\pi} \int_{0}^{\pi} dk \, \text{tr}[a(k,t)] \nonumber \\ &- \frac{1}{2\pi}\left(\chi(\pi) - \chi(0)\right), \quad t = 0, T/2 .\end{align} $ %% %And similarly for $P_2(t)$: %\begin{align} % P_2(t) &= \frac{1}{2\pi} \int_{-\pi}^{\pi} a_{22}(k,t) = \frac{1}{2\pi} \left( \int_{0}^{\pi} dk a_{22}(k,t) + \int_{-\pi}^{0} dk a_{22}(k,t) \right) ,\\ % &= \frac{1}{2\pi} \int_{0}^{\pi} dk \left(a_{22}(k,t) + a_{22}(-k,t)\right) ,\\ % &= \frac{1}{2\pi} \int_{0}^{\pi} dk \left(a_{22}(k,t) + a_{11}(k,t) - \frac{\partial \chi(-k)}{\partial k}\right),\\ % &= \frac{1}{2\pi} \int_{0}^{\pi} dk \, \text{tr}[a(k,t)] - \frac{1}{2\pi}\left(\chi(\pi) - \chi(0)\right), \quad t = 0, T/2 %.\end{align} %\hl{Wrong, this should not be the same as $P_1(t)$} %% Rewriting $\chi(k) = i \log(w_{12}(k))$, we get: $ \begin{align} P_1(t) = \frac{1}{2\pi} \int_{0}^{\pi} dk \, \text{tr}[a(k,t)] - \frac{i}{2\pi} \log\frac{w_{12}(\pi,t)}{w_{12}(0,t)}, \quad t = 0, T/2 .\end{align} $ Now, we can compute the [[time-reversal polarization]] for the system, and use the equation for the [[trace]] of $\pmb{a}(\pmb{k})$ above , and the Jacobi identity $\nabla_{\pmb{k}} \log(\text{det}[w(\pmb{k}]) = \text{tr}[w^\dagger(\pmb{k})\nabla_{\pmb{k}}w(\pmb{k})]$ to simplify the result: $ \begin{align} P_\theta(t) &= 2 P_1(t) - P_\rho (t) ,\\ &= \int_0^\pi \frac{dk}{2\pi} \left(\text{tr}[a(k,t)] - \text{tr}[a(-k,t)]\right) \nonumber \\ & \quad \quad - \frac{i}{\pi}\log\frac{w_{12}(\pi,t)}{w_{12}(0,t)},\\ &= \int_0^\pi \frac{dk}{2\pi} \left(i\text{tr}[w^\dagger(k,t) \nabla_{k} w(k,t)]\right) \nonumber \\ & \quad \quad - \frac{i}{\pi}\log\frac{w_{12}(\pi,t)}{w_{12}(0,t)},\\ &= \frac{i}{2\pi} \int_0^\pi dk \frac{\partial}{\partial k} \log(\text{det}[w(k,t)]) \nonumber \\ & \quad \quad - \frac{i}{\pi}\log\frac{w_{12}(\pi,t)}{w_{12}(0,t)},\\ &= \frac{i}{\pi} \log \sqrt{\frac{\text{det}[w(\pi,t)]}{\text{det}[w(0,t)]}} - \frac{i}{\pi}\log\frac{w_{12}(\pi,t)}{w_{12}(0,t)} .\end{align} $ For this model, the [[determinant]] is simple, $\text{det}[w(k,t)] = w_{12}(k,t)^2$, so we get: $ \begin{align} P_\theta(t) = \frac{1}{i\pi}\log\bigg(&\frac{\sqrt{w_{12}(0,t)^2}}{w_{12}(0,t)}\cdot \frac{w_{12}(\pi,t)}{\sqrt{w_{12}(\pi,t)^2}}\bigg), \quad t=0,T/2 .\end{align} $ The argument of the log is $\pm 1$. Thus, $P_\theta = 0$ or $1 \mod 2$ if we restrict the angle of the complex log to $0, 2\pi$ in the complex plane. Note that the [[Hilbert Space]] forms a [[torus]] in this case due to the [[periodic boundary conditions]], and the [[Bloch's Theorem|Bloch wavefunction]] $\ket{u_{nk}(t)}$ is a [[map]] from the phase space $(k,t)$ to the [[Hilbert Space]]. Based on the values of $P_\theta$ between times $t=0$ and $t=T/2$, we can classify the [[Hilbert Space]] into two classes, "twisted" and trivial. The twisted case occurs when $P_\theta(t=T/2) \neq P_\theta(t=0)$. We can define a quantity $\nu$: $ \begin{equation} \nu = P_\theta(t=T/2) - P_\theta(t=0), \end{equation} $ $\nu$ is only defined mod 2. Using the equation for $P_\theta$ above, we can write: $ \begin{align} \nu = \frac{1}{i\pi}\log\bigg(&\frac{\sqrt{w_{12}(0,T/2)^2}}{w_{12}(0,T/2)}\cdot\frac{w_{12}(\pi,T/2)}{\sqrt{w_{12}(\pi,T/2)^2}}\nonumber \\ &\cdot \frac{w_{12}(0,0)}{\sqrt{w_{12}(0,0)^2}}\cdot\frac{\sqrt{w_{12}(\pi,0)^2}}{w_{12}(\pi,0)}\bigg) \end{align} $ Or, if we identify $\Lambda_i = (0,0), (\pi,0), (0, T/2), (\pi,T/2)$ for $i=1\ldots4$, we get: $ \begin{equation} (-1)^\nu = \prod_{i=1}^4 \frac{w_{12}(\Lambda_i)}{\sqrt{w_{12}(\Lambda_i)^2}} \end{equation} $ ## General 1D Case and extension to 2D Assume we have 2N occupied bands, where bands $n$ and $n+1$ ($n = 1,3,\ldots,2N-1$) are [[Kramers Theorem|Kramers’ pairs]]. At the time-reversal symmetric times $t=0$ and $t=T/2$, we have: $ \begin{align} \Theta \ket{u_{n+1,k}(t)} &= e^{-i\chi_n(k)} \ket{u_{n,-k}(t)} , \\ \Theta \ket{u_{nk}(t)} &= e^{-i\chi_n(-k)} \ket{u_{n+1,-k}(t)} \end{align} $ Thus, the $w$ matrix is: $ \begin{align} &w\left(k, t\in\{0,T/2\}\right)= \nonumber \\ &\begin{pmatrix} 0 & e^{-i \chi_1(k)} & 0 & 0 & \ldots\\ -e^{-i \chi_1(-k)} & 0 & 0 & 0 & \ldots \\ 0 & 0 & 0 & e^{-i \chi_2(k)} & \ldots\\ 0 & 0 & -e^{-i \chi_2(-k)} & 0 & \ldots \\ \vdots & \vdots & \vdots & \vdots & \ddots \end{pmatrix} \end{align} $ Thus, at $t=0$ and $t=T/2$ we get: $ \begin{align} w_{12}(\Lambda_i)w_{34}(\Lambda_i) \ldots w_{2N-1,2N}(\Lambda_i) \nonumber \\ = \prod_{n=1}^{N} e^{-i\chi_n(\Lambda_i)} = e^{-i\sum_{n=1}^{N}\chi_n(\Lambda_i)} \nonumber \\ = \text{Pf}[w(\Lambda_i)] = \sqrt{\text{det}[w(\Lambda_i)]}. \end{align} $ Where in the last step we used the property of the [[Pfaffian]] of [[skew-symmetric matrix|antisymmetric matrices]] $\text{Pf}[A]^2 = \text{det}[A]$. Going through a calculation identical to that in the previous subsection, we get: $ \begin{align} P_1(t) &= \int_0^\pi A(k) - \frac{i}{2\pi}\log\left(\frac{\text{Pf}[w(\pi,t)]}{\text{Pf}[w(0,t)]}\right), \\ & \quad \quad t=0, T/2 \nonumber .\end{align} $ Which leads to the time-reversal polarization: $ \begin{align} P_\theta(t) &= \frac{i}{\pi}\log\left(\frac{\sqrt{\text{det}[w(0,t)]}}{\text{Pf}[w(0,t)]}\frac{\text{Pf}[w(\pi,t)]}{\sqrt{\text{det}[w(\pi,t)]}}\right), \quad \quad t=0, T/2 .\end{align} $ And finally $\nu$: $ \begin{align} (-1)^\nu = \prod_{i=1}^{4} \frac{\text{Pf}[w(\Lambda_i)]}{\sqrt{\text{det}[w(\Lambda_i)]}} .\end{align} $ Observe that if we set $T=2\pi$, then the [[Time-Reversal Invariant Momenta|TRIM]] points are $\Lambda_i = (0,0), (\pi,0), (0,\pi), (\pi,\pi)$. In this case, the periodic [[Phase Space]] $(k,t)$ of this 1D system is [[Isomorphism|isomorphic]] to the [[Brillouin Zone|BZ]] of a 2D system $(k_x, k_y)$. # Extension to 3D In three-dimensions, there are 8 [[Time-Reversal Invariant Momenta|TRIM]] points, forming 6 "time-reversal invariant planes." One would thus expect 6 $Z_2$ [[Z2 Topological Invariant|topological invariants]] analogously to the 2D case. However, as shown in by Moore and Balents, there are only four independent $Z_2$ invariants. We relabel the TRIM points as: $ \begin{align*} \Lambda_{n_1 n_2 n_3} = \frac{n_1}{2} \pmb{b}_1 + \frac{n_2}{2} \pmb{b}_2 + \frac{n_3}{2} \pmb{b}_3 .\end{align*} $ where $n_i = 0, 1, i = 1,\ldots 3$ and $\pmb{b}_i$ is the $i^{\text{th}}$ [[primitive reciprocal lattice vector]]. For each TRIM we define: $ \begin{align*} \delta(\Lambda_{n_1 n_2 n_3}) \equiv \frac{\text{Pf}[w(\Lambda_{n_1 n_2 n_3})]}{\sqrt{\text{det}[w(\Lambda_{n_1 n_2 n_3})]}} .\end{align*} $ And the [[Z2 Topological Invariant|topological invariants]] $\nu_i$, $i=0,1,2,3$ are given by: $ \begin{align*} (-1)^{\nu_0} = \prod_{n_1,n_2,n_3=0}^{1} \delta(\Lambda_{n_1 n_2 n_3}), \\ (-1)^{\nu_i} = \prod_{n_i=1,n_{j\neq i}=0}^{1} \delta(\Lambda_{n_1 n_2 n_3}) .\end{align*} $ The final $Z_2$ invariant is often written as $Z_2=(\nu_0;\nu_1 \nu_2 \nu_3)$. The first index is known as the strong or global index, and the other three are known as weak indices. There are three classes of topological insulators in 3D: * Strong TI: $\nu_0 = 1$, $\nu_{i\neq 0}$ can be zero or one, e.g., $Z_2 = (1;000)$ * Weak TI: $\nu_0 = 0$, At least one $\nu_{i\neq 0} = 1$ e.g., $Z_2 = (0;001)$ * Trivial Insulator: $\nu_i = 0$, $i = 0, \ldots 3$, i.e., $Z_2 = (0;000)$ # Systems with Inversion Symmetry In this section, we work on a general system with [[Time-Reversal Operator|time-reversal]] as in the previous system. However, we also assume the system is [[Inversion Operator|inversion symmetric]], i.e. $\Pi H \Pi^{-1} = H$. The [[Bloch's Theorem|Bloch functions]] are $\ket{u_\alpha(\pmb{k}}$, where $\alpha$ is a band index. At the [[Time-Reversal Invariant Momenta|TRIM]] $\Lambda_i$, the [[Inversion Operator|parity]] [[eigenvalue]] are $\xi_\alpha(\Lambda_i) = \pm 1$, i.e. $\Pi \ket{u_\alpha(\Lambda_i)} = \xi_\alpha(\Lambda_i) \ket{u_\alpha(\Lambda_i)}$. As before, we assume the system has $2N$ bands, where bands $n$ and $n+1$ are Kramers pairs ($n = 1, 3, \ldots, N-1$). The [[Berry connection]], $\pmb{A}(\pmb{k}) = \text{tr}[\pmb{a}(\pmb{k})]$ is related to the [[Berry curvature]] $\pmb{F}(\pmb{k})$ via the relation: $ \begin{align} \pmb{F}(\pmb{k}) = \nabla_{\pmb{k}} \times \pmb{A}(\pmb{k}) .\end{align} $ In systems with [[Time-Reversal Operator|time-reversal symmetry]], it is well known that $\pmb{F}(-\pmb{k}) = -\pmb{F}(\pmb{k})$. In systems with [[Inversion Operator|inversion symmetry]], we additionally have $\pmb{F}(-\pmb{k}) = +\pmb{F}(\pmb{k})$. Thus, in systems with both TRS and IS, the [[Berry curvature]] is identically zero $\pmb{F}(\pmb{k}) \equiv 0$. Therefore, we can always pick a [[gauge]] in which $\pmb{A}(\pmb{k} = 0$ for all $\pmb{k}$. We can define a new [[Matrix]] $v_{\alpha \beta}(\pmb{k})$: $ \begin{align} v_{\alpha \beta}(\pmb{k}) \equiv \bra{u_{\alpha}(\pmb{k})} \Pi \Theta \ket{u_{\beta}(\pmb{k})} .\end{align} $ where $\alpha$ and $\beta$ are [[Band Structure|band]] indices. It is straightforward to show that $v(\pmb{k})$ is [[skew-symmetric matrix|antisymmetric]] in the band indices (i.e. $v^T = - V$ or $v_{\alpha \beta}(\pmb{k}) = - v_{\beta \alpha}(\pmb{k})$) and [[Unitary Operator|unitary]] (i.e. $v(\pmb{k})^\dagger v(\pmb{k}) = I$ or $\sum_{\beta} v_{\alpha \beta}(\pmb{k}) v_{\gamma \beta}^*(\pmb{k}) = \delta_{\alpha \gamma}$. Moreover, $v(\pmb{k})$ is related to the [[Berry connection]] via: $ \begin{align} \pmb{A}(\pmb{k}) &= \frac{i}{2} \text{tr}[v^\dagger \nabla_{\pmb{k}} v]\\ &= \frac{i\nabla_{\pmb{k}}}{2} \log\left(\text{det}[v(\pmb{k})]\right)\\ &= i \nabla_{\pmb{k}} \text{Pf}[v(\pmb{k})] .\end{align} $ Thus, to set the [[Berry connection]] $\pmb{A}(\pmb{k}) = 0$, we need to pick a gauge in which $\text{Pf}[v(\pmb{k})] = 1$, i.e. pick the phase of $\ket{u_{\alpha}(\pmb{k})}$ to guarantee that condition. Back to the [[sewing matrix]], we have: $ \begin{align} w_{\alpha \beta}(\Lambda_i) &= \bra{u_\alpha(\Lambda_i)}\Theta\ket{u_\beta(\Lambda_i)}, \\ &= \bra{u_\alpha(\Lambda_i)}\Pi \Pi \Theta\ket{u_\beta(\Lambda_i)}, \\ &= \xi_\alpha(\Lambda_i) \bra{u_\alpha(\Lambda_i)} \Pi \Theta\ket{u_\beta(\Lambda_i)}, \\ &= \xi_\alpha(\Lambda_i) v_{\alpha \beta}(\Lambda_i) .\end{align} $ Note that [[index notation|repeated indices]] *do not* imply [[Einstein Summation Convention|summation]] in the last equation. Thus: $ \begin{align} \text{Pf}[w_{\alpha \beta}(\Lambda_i)] = \prod_{n=1}^{2N}\xi_{2n}(\Lambda_i) \text{Pf}[v_{\alpha \beta}(\Lambda_i)] .\end{align} $ However, since we have fixed the gauge to ensure that $\text{Pf}[v_{\alpha \beta}(\pmb{k})] = 1$, then $\text{Pf}[w_{\alpha \beta}(\Lambda_i)] = \prod_{n=1}^{2N}\xi_{2n}(\Lambda_i)$. Thus, we finally end up with the well-known [[Fu-Kane formula]] for systems with [[Inversion Operator|inversion symmetry]]: $ \begin{align} \delta(\Lambda_i) &= \frac{\text{Pf}[w(\Lambda_i)]}{\sqrt{\text{det}[w(\Lambda_i)]}} = \text{sign}({\text{Pf}[w(\Lambda_i)]}) \nonumber \\ &= \prod_{n=1}^{2N}\xi_{2n}(\Lambda_i) .\end{align} $ # Computational Methods ## Systems with IS: parity eigenvalues In systems with [[Inversion Operator|inversion symmetry]], the easiest method to calculate the $Z_2$([[Cyclic Group of Order 2]]) invariants is to compute the parity [[eigenvalue]]. Starting from a well-converged *ab initio* [[Self-Consistent Field Calculation|self-consistent field]] (SCF) calculation, one can compute the [[Irreducible Representation|irreducible representation]] of the [[double group|double group]] at the $k$ points of interest (i.e. the [[Time-Reversal Invariant Momenta|TRIM]]). There are multiple methods to compute the [[Irreducible Representation]], but they are well beyond the scope of this paper. An example of a software that implements this function is `bands.x` distributed with `Quantum ESPRESSO`. ## Wannier Charge Center While there are a few other methods to compute the $Z_2$ invariant in non-[[centrosymmetric]] systems, the most robust is the [[Wannier charge center]] (WCC) method, which is also known as the [[Wilson loop]] method. This method relies on computing the [[maximally localized Wannier functions]] (MLWFs) of the system before starting the calculation of the $Z_2$ invariants. For a review of MLWFs, we refer the readers to [[Soluyanov2011]]. In this section, we will focus on Fu and Kane's 2 band 1D model, and direct readers to [[Pizzi2020]] for details on the computational implementation and the straightforward extension to higher dimensions. We assume the system has $2N$ bands (i.e, $N$ [[Kramers Theorem|Kramers' pairs]]). Recall that [[maximally localized Wannier functions|Wannier functions]] for 1D systems are defined as: $ \begin{align} \ket{W_n(x-x_j)} = \frac{1}{2\pi} \int_{-\pi}^{\pi} \text{d} k e^{i k (x - x_j} \ket{u_{nk}} ,\end{align} $ where $x_j$ is the position of the [[unit cell]] of interest, and $n = 1,\ldots2N$. We then define the [[Wannier charge center]] (WCC): $ \begin{align} \overline{x}_n &= \bra{W_n(x)}\hat{x}\ket{W_n(x)} \\ &= \frac{i}{2\pi} \int_{-\pi}^{\pi} \text{d} k \bra{u_{nk}} \partial_k \ket{u_{nk}} .\end{align} $ That is, the [[Wannier Charge Center|WCC]] $\overline{x}_n$ is the [[Expectation Value]] of the [[position operator]] in the home cell ($x_j = 0$) for band $n$. We relabel the bands such that $\ket{u_{nk}} \to \ket{u_{\alpha k}^S}$ where $S = 1, 2$ and $\alpha = 1,\ldots N$. This means that $\ket{u_{\alpha k}^1}$ and $\ket{u_{\alpha k}^2}$ are a [[Kramers Theorem|Kramers] pair]]. Note that: $ \begin{align} \sum_{\alpha} \overline{x}_{\alpha}^S = \frac{i}{2\pi} \oint_{\text{BZ}} \text{d} k \sum_\alpha \bra{u_{\alpha k}^S} \partial_k \ket{u_{\alpha k}^S} = \frac{1}{2\pi} \oint_{\text{BZ}} A^S .\end{align} $ Thus, we obtain the [[charge polarization]] and [[time-reversal polarization]]: $ \begin{align} P_\rho = P^1 + P^2 &= \frac{1}{2\pi} \oint \dd k (A^{1}(k) + A^{2}(k)) \\ &= \sum_\alpha \overline{x}_{\alpha}^1 + \overline{x}_{\alpha}^{2}\\ P_\theta = P^1 = P^2 &= \sum_\alpha \overline{x}_{\alpha}^1 - \overline{x}_{\alpha}^{2} .\end{align} $ Therefore, we can finally get the 1D topological invariant: $ \begin{align} \delta = P_\theta(T/2) - P_\theta(0) &= \sum_\alpha \overline{x}_{\alpha}^1(T/2) - \overline{x}_{\alpha}^{2}(T/2) \nonumber \\&- \sum_\alpha \overline{x}_{\alpha}^1(0) - \overline{x}_{\alpha}^{2}(0) .\end{align} $ Note that there are issues with [[gauge]] and [[Continuous Map|continuity]] discussed in more detail in the work by Soluyanov. The implementation of this method is straightforward in theory: one needs to compute the [[Wannier Charge Center|WCCs]] and plot them at every value of $t$, and visually track their evolution. However, this is tedious and the authors of [[Soluyanov2011]] came up with an automated method that relies on tracking the largest gap in the WCC spectrum. This method is generally quite robust, extremely quick, and works on all systems irrespective of [[Inversion Operator|inversion symmetry]]. However, the drawbacks is the requirement to "wannierize" (i.e. find the [[maximally localized Wannier functions|MLWFs]]) the system first. Wannierization can be a tricky process for some materials and, while it is not computationally expensive, requires significant time for trial and error. There are several available software packages for this method. For example, `WannierTools` is a very robust package that can calculate the $Z_2$ [[Z2 Topological Invariant|invariants]] as well as quite a few other quantities. It is parallelized and quite efficient. However, it is not easily extensible and difficult to interface with. Another alternative is `Z2Pack`, which functions in a very similar way but is built using Python. Thus, it is slower but much easier to interface with and quite well documented. Both packages rely on Wannier90 for construction of the MLWFs.