Topological Dirac Nodal-net Fermions in AlB$_2$-type TiB$_2$ and ZrB$_2$

TL;DR

This paper predicts a topological Dirac nodal-net semimetal state in TiB₂ and ZrB₂ using first-principles calculations, revealing unique topological features and surface states protected by symmetries.

Contribution

It identifies a novel Dirac nodal-net semimetal state in TiB₂ and ZrB₂, with detailed analysis of topological excitations and surface Dirac points, expanding understanding of topological semimetals.

Findings

Nodal-net state in TiB₂ and ZrB₂ identified without SOC

Presence of triple point, nexus, and nodal link topological excitations

Surface Dirac points analogous to graphene on TiB₂ surfaces

Abstract

Based on first-principles calculations and effective model analysis, a Dirac nodal-net semimetal state is recognized in AlB$_2$-type TiB$_2$ and ZrB$_2$ when spin-orbit coupling (SOC) is ignored. Taking TiB$_2$ as an example, there are several topological excitations in this nodal-net structure including triple point, nexus, and nodal link, which are protected by coexistence of spatial-inversion symmetry and time reversal symmetry. This nodal-net state is remarkably different from that of IrF$_4$, which requires sublattice chiral symmetry. In addition, linearly and quadratically dispersed two-dimensional surface Dirac points are identified as having emerged on the B-terminated and Ti-terminated (001) surfaces of TiB$_2$ respectively, which are analogous to those of monolayer and bilayer graphene.

Tables

Table 1. Table 1: Parity of the occupied bands of TiB 2 at TRIMs, Γ i subscript Γ 𝑖 \Gamma_{i} is in unit of the reciprocal lattice vectors. ’Total’ means the product of all parities of occupied bands. Γ 0 subscript Γ 0 \Gamma_{0} , Γ 1 ​ ( Γ 2 ) subscript Γ 1 subscript Γ 2 \Gamma_{1}(\Gamma_{2}) , Γ 3 subscript Γ 3 \Gamma_{3} and Γ 5 subscript Γ 5 \Gamma_{5} are the same as the notional Γ Γ \Gamma , M 𝑀 M , A 𝐴 A and L 𝐿 L respectively.

TRIM				parity			Total
${\Gamma }_0$ (0.0, 0.0, 0.0)	+	+	+	-	+		-
${\Gamma }_1$ (0.5, 0.0, 0.0)	+	-	-	+	+		+
${\Gamma }_2$ (0.0, 0.5, 0.0)	+	-	-	+	+		+
${\Gamma }_3$ (0.0, 0.0, 0.5)	-	+	-	-	+	+	-
${\Gamma }_4$ (0.5, 0.5, 0.0)	+	-	-	+	+		+
${\Gamma }_5$ (0.0, 0.5, 0.5)	+	-	+	+	-		+
${\Gamma }_6$ (0.5, 0.0, 0.5)	+	-	+	+	-		+
${\Gamma }_7$ (0.5, 0.5, 0.5)	+	-	+	+	-		+

Equations

$${ϵ}_1=\frac{1}{2}(\alpha +{\gamma }_1+{\gamma }_2+\sqrt{(\alpha +{\gamma }_1-{\gamma }_2{)}^2+8{\beta }^2})$$

$${ϵ}_2={\gamma }_1-\alpha$$

$${ϵ}_3=\frac{1}{2}(\alpha +{\gamma }_1+{\gamma }_2-\sqrt{(\alpha +{\gamma }_1-{\gamma }_2{)}^2+8{\beta }^2})$$

$$(C^2+A_2-A_1)k_x^2+(B_2-B_1-D)k_z^2=E_1-E_2$$

$$H(k)=(M+v{k}_x+{\alpha }_1{k}_x^2+{\alpha }_2{k}_y^2+{\alpha }_3{k}_z^2){\sigma }_z+\beta k_y{k}_z{\sigma }_x$$

$$E_{\pm }=\pm \sqrt{(v{k}_x+{\alpha }_1{k}_x^2+{\alpha }_2{k}_y^2+{\alpha }_3{k}_z^2{)}^2+{\beta }^2{k}_y^2{k}_z^2}$$

$$(-1{)}^{\nu }={\xi }_a{\xi }_b;{\xi }_a=\prod _n{\xi }_n({\Gamma }_a)$$

$$H(\mathbf{k})=V_F(k_x{\sigma }_x+k_y{\sigma }_y)$$

$$\displaystyle H(\bf{k})=\left(\begin{array}[]{cc}Bk^{2}&Ak_{+}+Ck_{-}^{2}\\ Ak_{-}+Ck_{+}^{2}&Bk^{2}\end{array}\right)$$

$$D(R)H(k)D^{†}(R)=H(R(k))$$

$$D(I)=\text{diag}\{1,1,1\}$$

$$D(c_{6z})=\text{diag}\{-e^{i2\pi /3},-e^{-i2\pi /3},1\}$$

$$\displaystyle D(c_{2y})=\left(\begin{array}[]{ccc}0&1&0\\ 1&0&0\\ 0&0&1\end{array}\right)$$

$$\displaystyle D(\text{TR})=\left(\begin{array}[]{ccc}0&1&0\\ 1&0&0\\ 0&0&1\end{array}\right)K$$

$$\displaystyle H(\bf{k})=\left(\begin{array}[]{ccc}\varepsilon_{1}({\bf k})&Ck_{-}^{2}&Dk_{-}k_{z}\\ Ck_{+}^{2}&\varepsilon_{1}({\bf k})&Dk_{+}k_{z}\\ Dk_{+}k_{z}&Dk_{-}k_{z}&\varepsilon_{2}({\bf k})\end{array}\right)$$

$$\displaystyle H(\bf{k})=\left(\begin{array}[]{ccc}\varepsilon_{1}({\bf k})&Ck_{-}^{2}+Fk_{+}^{4}&Dk_{-}k_{z}\\ Ck_{+}^{2}+Fk_{-}^{4}&\varepsilon_{1}({\bf k})&Dk_{+}k_{z}\\ Dk_{+}k_{z}&Dk_{-}k_{z}&\varepsilon_{2}({\bf k})\end{array}\right)$$

$$D({\sigma }_h)=\text{diag}\{-1,-1\}\oplus \text{diag}\{1,1\}$$

$$\displaystyle D(c_{2x})=\left(\begin{array}[]{cc}0&-1\\ -1&0\end{array}\right)\oplus\left(\begin{array}[]{cc}0&-1\\ -1&0\end{array}\right)$$

$$\displaystyle D(c_{3z})=\left(\begin{array}[]{cc}e^{-2i\pi/3}&0\\ 0&e^{2i\pi/3}\end{array}\right)\oplus\left(\begin{array}[]{cc}e^{2i\pi/3}&0\\ 0&e^{-2i\pi/3}\end{array}\right)$$

$$H(k)=$$

$$\displaystyle\left(\begin{array}[]{cccc}\varepsilon(k)&Ck_{+}+Dk_{-}^{2}&iEk_{z}k_{+}&Fk_{z}\\ Ck_{-}+Dk_{+}^{2}&\varepsilon(k)&Fk_{z}&-iEk_{z}k_{-}\\ -iEk_{z}k_{-}&Fk_{z}&\varepsilon^{\prime}(k)&C^{\prime}k_{-}+D^{\prime}k_{+}^{2}\\ Fk_{z}&iEk_{z}k_{+}&C^{\prime}k_{+}+D^{\prime}k_{-}^{2}&\varepsilon^{\prime}(k)\end{array}\right)$$

$$\displaystyle H(k)=\left(\begin{array}[]{cccc}H_{11}&H_{12}\\ H_{12}^{{\dagger}}&H_{22}\end{array}\right)$$

$$H_{11}=$$

$$H_{22}=$$

$$H_{12}=$$

$$D(I)=\text{diag}\{1,1,1,1,1\}$$

$$\displaystyle D(c_{2y})=1\oplus\left(\begin{array}[]{cc}0&1\\ 1&0\end{array}\right)\oplus\left(\begin{array}[]{cc}0&-1\\ -1&0\end{array}\right)$$

$$D(c_{6z})=\text{diag}\{1,-e^{2i\pi /3},-e^{-2i\pi /3},e^{-2i\pi /3},e^{2i\pi /3}\}$$

$$\displaystyle D(\text{TR})=1\oplus\left(\begin{array}[]{cc}0&1\\ 1&0\end{array}\right)\oplus\left(\begin{array}[]{cc}0&-1\\ -1&0\end{array}\right)K$$

$$\displaystyle H(k)=\left(\begin{array}[]{ccccc}\varepsilon_{2}(k)&Dk_{z}k_{+}&Dk_{z}k_{-}&Ek_{+}^{2}&-Ek_{-}^{2}\\ Dk_{z}k_{-}&\varepsilon_{1}(k)&Ck_{-}^{2}&Fk_{z}k_{+}&0\\ Dk_{z}k_{+}&Ck_{+}^{2}&\varepsilon_{1}(k)&0&-Fk_{z}k_{-}\\ Ek_{-}^{2}&Fk_{z}k_{-}&0&\varepsilon_{3}(k)&Gk_{+}^{2}\\ -Ek_{+}^{2}&0&-Fk_{z}k_{+}&Gk_{-}^{2}&\varepsilon_{3}(k)\end{array}\right)$$

$$\displaystyle H(k)=\left(\begin{array}[]{cccc}\varepsilon_{1}(k)&Ck_{-}^{2}&0&0\\ Ck_{+}^{2}&\varepsilon_{1}(k)&0&0\\ 0&0&\varepsilon_{3}(k)&Gk_{+}^{2}\\ 0&0&Gk_{-}^{2}&\varepsilon_{3}(k)\end{array}\right)$$

$$\displaystyle D(c_{3z})=\left(\begin{array}[]{cc}e^{-2i\pi/3}&0\\ 0&e^{2i\pi/3}\end{array}\right),\sigma_{v}=\left(\begin{array}[]{cc}0&-1\\ -1&0\end{array}\right)$$

$$\displaystyle H(\bf{k})=\left(\begin{array}[]{cc}Bk^{2}&Ak_{+}+Ck_{-}^{2}\\ Ak_{-}+Ck_{+}^{2}&Bk^{2}\end{array}\right)$$

Full text

I Introduction

In the band theory of solids, a nodal point is a point of contact between conduction and valence bands. One way to classify nodal points is based on their dimensionality Chiu et al. (2016); Bzdušek et al. (2016); Kobayashi et al. (2017); Burkov et al. (2011), which can be divided into three categories. The first is zero-dimensional (0D) nodal points including Weyl points Wan et al. (2011); Weng et al. (2015a); Huang et al. (2015a); Xu et al. (2015); Lv et al. (2015); Soluyanov et al. (2015), Dirac points Novoselov et al. (2005); Wang et al. (2012); Liu et al. (2014); Young et al. (2012), triple points Zhu et al. (2016); Weng et al. (2016); Chang et al. (2016), and other higher degeneracy nodal points Bradlyn et al. (2016). The second is one-dimensional (1D) nodal line systems which include nodal ring Hu et al. (2016); Bian et al. (2016), nodal chain Bzdušek et al. (2016); Yu et al. (2017a) and nodal net Bzdušek et al. (2016). The third is two-dimensional (2D) nodal surfaces Liang et al. (2016). A semimetal with such nodal points is called a topological semimetal (TS), such as Weyl semimetal (WSM), Dirac semimetal (DSM), Dirac nodal line semimetal (DNLSM), nodal chain semimetal, et al.

The 0D nodal point system has been extensively studied during the past decade, due to its non-trivial topological properties and its exotic transport properties. For example, Dirac semimetal has been predicted to be a good candidate for quantum devices because of its massless Dirac fermion’s properties Novoselov et al. (2005); Zhang et al. (2005); Wang et al. (2012); Liu et al. (2014), Weyl semimetals have chiral anomaly which leads to a chiral magnetic effect Son and Yamamoto (2012); Song et al. (2016) i.e. an electric current parallel to an external magnetic field Vazifeh and Franz (2013); Fukushima et al. (2008); Ali et al. (2014); Zyuzin and Burkov (2012); Hosur and Qi (2013); Huang et al. (2015b), and a triple point metal could have topological Lifshitz transitions Zhu et al. (2016).

Investigations of 1D nodal line systems have been developing recent years, due to their great potential in diverse developments in materials science. A nodal line system has a non-trivial $\pi$ Berry phase around the nodal line which would shift the Landau level index by 1/2 Carmier and Ullmo (2008); Yan et al. (2017), and leads to drumhead surface states (SSs) Fang et al. (2016); Yu et al. (2017b). There are many types of nodal line systems, such as single and multiple Dirac nodal line (DNL) systems in the absence of spin-orbit coupling (SOC) Burkov et al. (2011); Bian et al. (2016); Hirayama et al. (2017); Cheng et al. (2017); Liang et al. (2016); Kim et al. (2015); Fang et al. (2015), nexus Heikkila and Volovik (2015); Hyart and Heikkila (2016); Zhu et al. (2016), nodal chain with SOC Bzdušek et al. (2016), and other crossing nodal lines Yu et al. (2017a); Weng et al. (2015b); Kim et al. (2015); Yu et al. (2015); Yan et al. (2017); Chang and Yee (2017); Chen et al. (2017); Kobayashi et al. (2017). All the proposed nodal line systems are protected by crystallographic symmetry. Some are protected by the coexistence of time-reversal symmetry $T$ and spatial-inversion symmetry $P$, namely $PT$ symmetry Cheng et al. (2017); Bian et al. (2016); Liang et al. (2016); Kim et al. (2015). Some are protected by mirror symmetry or glide symmetry Chiu and Schnyder (2014); Bzdušek et al. (2016).

In this work, based on first-principles calculations, we theoretically study a new type of nodal structure, namely nodal net, in two AlB2-type diborides TiB2 and ZrB2, which present a unique combination of properties such as high bond strengths, high melting points, high thermal conductivities, low electrical resistance and low work functions and can be easily synthesized in the laboratory Waśkowska et al. (2011); Okamoto et al. (2010); Kumar et al. (2012); Wang et al. (2011). To better understand this complex nodal net structure, we find it including four classes nodal line structures: A, nodal ring in $k_{z}=0$ plane surrounding $K$ point; B, nodal line in three vertical mirror planes $\sigma_{v1}$, $\sigma_{v2}$ and $\sigma_{v3}$ (see definition on p. 273 of  Altmann and Herzig (1994)); C, nodal line along $\Gamma-A$ starting from a triple point; and D, a single isolated nodal ring at $k_{z}$=0.5 plane surrounding $A$ point. Furthermore, class-A and class-B nodal lines will cross at a k-point along the $\Gamma-K$ direction. All three class-B nodal lines in the vertical mirror planes terminate at $A$ point, which is also a termination of the class-C nodal line. So $A$ point is also called a neuxs point Cornwall (1999); Heikkila and Volovik (2015), which is the termination of several Dirac line nodes. This nodal net structure differs from previously reported 0D and 1D nodal structures which may lead to some new magnetic and electrical transport properties. In addition, the linearly and quadratically dispersed surface Dirac cones are found at the $\bar{K}$ point of the surface BZ for B-terminated and Ti-terminated surfaces of TiB2 respectively.

This paper is organized as follows. In Section.II, we elucidate the crystal and electronic structure of TiB2. In Section.III, we study the complex nodal net structure both using first-principles calculations and effective $k\cdot p$ model analysis. In Section.IV, we study the drumhead surface states both for B-terminated and Ti-terminated surfaces of TiB2, where a linearly and a quadratically dispersed surface Dirac points are found.

II Crystal and band structure of AlB2-type TiB2

The crystallographic data of TiB2 and ZrB2 are obtained from Ref.Kumar et al. (2012). TiB2 and ZrB2 have the same AlB2-type centro-symmetric crystal structures with the space group P6/mmm (191). As verified by calculation, we find TiB2 (Fig.2) and ZrB2 (Fig 6 in Appendix) have similar electronic structures, therefore, we take TiB2 as an example hereafter. As shown in Fig.1a, it is a layered hexagonal structure with alternating close-packed hexagonal layers of titanium and graphene-like boron layers. The optimized lattice constants are a=b=3.0335 $\r{A}$ and c=3.2263 $\r{A}$, which agree well with the experimental Post et al. (1954) and other theoretical Milman and Warren (2001) results.

To study the electronic properties of TiB2, the electronic band structure (BS) is calculated in absence of SOC (see details in Appendix A), as shown in Fig.2a. It shows that the valence and conduction bands near the Fermi level exhibit Dirac linear dispersion. There are six band crossing points (also called nodal points (NPs)) located along the H-$\Gamma$, $\Gamma$-A, A-H, K-$\Gamma$, M-K and L-A lines (marked as a to f points in Fig.2a). It is noticed that these six NPs deviate from the Fermi level about -0.16 eV, 0.5 eV, 0.35 eV, -0.01 eV, -0.28 eV and 0.32 eV, respectively. To investigate the formation mechanism of these crossing points, orbital-character analysis is performed. As shown in Fig.2a, the Ti 4d states in TiB2 are dominant feature for these six NPs. Specifically, a is dominated by Ti-dxz, Ti-dx2-y2, and Ti-dz2 orbitals; b is dominated by Ti-dyz and Ti-dz2 orbitals; c is dominated by Ti-dxy, Ti-dyz and Ti-dx2-y2 orbitals; d is dominated by Ti-dxz and Ti-dz2 orbitals; e is dominated by Ti-dxz, Ti-dx2-y2 and Ti-dz2 orbitals, and f is dominated by Ti-dyz and Ti-dx2-y2 orbitals, respectively. The Fermi surface (FS) of TiB2 is calculated and shown in Fig. 2b and 2c, which shows a lantern-like frame with compensated electron pockets and hole pockets, which is a feature of a topological semimetal and also will lead to non-saturated large positive magnetoresistance Ali et al. (2014). The FS calculated in this work agrees well with previous experimental Tanaka and Ishizawa (1980) and theoretical studies Kumar et al. (2012).

In presence of the SOC effect, the crossing points along the H-$\Gamma$, A-H, K-$\Gamma$, M-K and L-A lines are fully gapped, which is common in $PT$ protected systems Kim et al. (2015); Fang et al. (2016). SOC makes the crossing points open gaps about 26 meV, 18 meV, 25 meV, 23 meV, and 21 meV [Fig. 5 in Appendix], respectively. The millivolt level gaps indicated that the effect of SOC on the electronic band structure of TiB2 is quite weak and can be ignored in experimental work. One thing worth mentioning is that the SOC splitting would generate a Dirac point along the $\Gamma-A$ direction [Fig.5b in Appendix] although it happens between the (N+2)’th and (N+3)’th bands, where N is the number of occupied bands at $\Gamma$ point in the BZ .

III Nodal net structure

From previous studies, it is known that a nodal line system would have banana-shaped linked FSs Parker et al. (2009); Takane et al. (2016) as the NPs do not usually align at the same energy level. In other words, these banana-shape linked FSs are an indication of the existence of a nodal line structure. There is another clue to find NLSM for a $PT$ symmetry protected system, i.e., if there is a band touching point close to the Fermi level, there should be a nodal line including this point Fang et al. (2016). Based on these two clues, and the BS and FS shown in Fig.2, it is clear that there is a nodal line structure in the TiB2 system. By using the symmetrical Wannier tight binding model Zhu et al. (2016) and a software package WannierTools Wu et al. (2017), we found all k-points with zero local energy gap between the N’th and (N+1)’th energy bands, where N is the number of occupied bands at the $\Gamma$ point, i.e.: $\Delta({\bf k})=E_{N+1}({\bf k})-E_{N}({\bf k})=0$. The nodal points are plotted in Fig.3, from which, the energies of nodal points are not the same which leads to the lantern-like FSs as shown in Fig.2b.

TiB2 has $PT$ symmetry, which is enough to protect the existence of Dirac nodal lines in the absence of SOC. While, besides $PT$ symmetry, there are another four mirror symmetries $\sigma_{h}$, $\sigma_{v1}$, $\sigma_{v2}$ and $\sigma_{v3}$ in the $D_{6h}$ group. Such mirror symmetries will enforce the nodal lines embed on mirror planes. Thus, these nodal lines form a interconnected nodal net structure including four classes nodal lines: A, B, C, and D. In the Appendix, with the aid of DFT calculations, it is verified that class-A,B, and D nodal lines still exist but apart from the previous mirror plane by breaking the mirror symmetries, however the class-C nodal lines will disappear.

Class-A nodal lines. Those nodal rings surrounding $K$ point, embed in the $k_{z}=0$ plane, which is a mirror plane $\sigma_{h}$ of $D_{6h}$, and is shown as six arcs around $K$ point in Fig.3. The effective $k\cdot p$ model at $K$ point was constructed within the little group $D_{3h}$, and shown in Eq. (45). When $k_{z}=0$, $H_{12}=0$, which leads to two uncoupled blocks. Those two eigenvalues of each block lead to one upward and one downward parabola. So, if the energies of the two blocks at $K$ point are different, there will be a nodal ring surrounding $K$ point.

Class-B nodal lines. Those Weyl nodes sitting on the vertical mirror planes $\sigma_{v1}$, $\sigma_{v2}$ and $\sigma_{v3}$, are shown as edges of the lantern in Fig.3a. We could use the effective $k\cdot p$ model at $\Gamma$ point shown in Eq.(22) to prove the existence of a nodal line on such k-planes. For simplicity, we choose $\sigma_{v1}(\sigma_{xz})$, a mirror perpendicular to the y-axis as an example. On this plane, $k_{y}=0$, then the eigenvalues of Eq.(22) are

[TABLE]

where $\alpha=Ck_{x}^{2}$, $\beta=Dk_{x}k_{z}$ and $\gamma_{i}=E_{i}+A_{i}k_{x}^{2}+B_{i}k_{z}^{2}$ with i=1, 2. At $\Gamma$ point, $k_{x}=0,k_{z}=0$, $\epsilon_{1}=E_{1}$, $\epsilon_{2}=E_{1}$, $\epsilon_{3}=E_{2}$. It is clear that, along $\Gamma-A$, where $k_{x}=0$, $\epsilon_{1}=\epsilon_{2}$, which is evidence of the existence of class-C nodal lines. The nodal line in classes other than class-C could exist only if $\epsilon_{1}=\epsilon_{3}$ or $\epsilon_{2}=\epsilon_{3}$. $\epsilon_{1}=\epsilon_{3}$ would leads to constraint $E_{2}-E_{1}=Ck_{x}^{2}$ and $k_{x}k_{z}=0$, which results in a touching point $k_{z}=0,k_{y}=0,k_{x}=\sqrt{(E_{2}-E_{1})/C}$. However, the DFT fitting results shows that $E_{2}-E_{1}<0$. So there is no touching point between $\epsilon_{1}$ and $\epsilon_{3}$. Another possibility, $\epsilon_{2}=\epsilon_{3}$, would lead to:

[TABLE]

There are three possibilities arising from the relationship between the parameters in Eq.1

If $(C^{2}+A_{2}-A_{1})(E_{1}-E_{2})>0$ and $(B_{2}-B_{1}-D)(E_{1}-E_{2})>0$, then there will be an elliptic nodal ring centered at $\Gamma$ point. 2. 2.

If $(C^{2}+A_{2}-A_{1})(B_{2}-B_{1}-D)<0$, then there will be an hyperbolic nodal ring. 3. 3.

If $(C^{2}+A_{2}-A_{1})(E_{1}-E_{2})<0$ and $(B_{2}-B_{1}-D)(E_{1}-E_{2})<0$, there will be no nodal line.

According to the DFT fitting parameters, TiB2 belongs to the first class, which will have a elliptical nodal ring on three vertical mirror planes $\sigma_{v1}$, $\sigma_{v2}$ and $\sigma_{v3}$.

Actually, there are another three mirror symmetries $\sigma_{d1}$, $\sigma_{d2}$ and $\sigma_{d3}$ in $D_{6h}$. Take $\sigma_{d1}(\sigma_{yz})$ for instance, we could perform the same analysis above by setting $k_{x}=0$; however, we have to admit that the results are the same as in the $k_{y}=0$ plane, i.e. there will be nodal line in the $k_{x}=0$ plane. The reason for this is that the $k\cdot p$ model in Eq.(22) is up to second order, which will lead to isotopic effects on $k_{x}$ and $k_{y}$. Eventually, there are not only nodal lines in the mirror planes, but also a nodal surface encompassing the $\Gamma$ point. So, in order to degenerate the nodal surface and distinguish between the $k_{x}=0$ and $k_{y}=0$ planes, we have to include higher order terms in the $k\cdot p$ model, as discussed in the Appendix.D.2.

Nodal link of class-A and class-B NLs. Fig.3a shows that the class-A NL is linked with class-B NL at O point along $\Gamma-K$ direction, which is also shown in Fig.3c. $\Gamma-K$ is the overlap of $\sigma_{v}$ and $\sigma_{h}$ mirror planes. Eventually, the little group of O is $c_{2v}$ with an additional $PT$ symmetry. The 2-band effective $k\cdot p$ model of O up to 2nd order of ${\bf k}$ is given as

[TABLE]

where $\sigma_{x}$ and $\sigma_{z}$ are Pauli-matrices. Since O is a nodal link point, so the mass term $M$ should be zero. Eigenvalues of Eq.2 are

[TABLE]

The nodal point exists only if $k_{y}k_{z}=0$ and $vk_{x}+\alpha_{1}k_{x}^{2}+\alpha_{2}k_{y}^{2}+\alpha_{3}k_{z}^{2}=0$, which leads to two class of NLs, one embeds in $k_{z}=0$ plane, which belongs to class-A NL, the other one embeds in $k_{y}=0$ plane, which belongs to class-B NL. It is clearly that both NLs link together at O. Taking class-A as an example, $k_{z}=0$, and the other condition becomes $\alpha_{1}(k_{x}+\frac{v}{2\alpha_{1}})^{2}+\alpha_{2}k_{y}^{2}=\frac{v^{2}}{4\alpha_{1}}$. It is learned that $v\neq 0$ is the only condition that keeps the NL. $\alpha_{1}\alpha_{2}>0$ would lead the NL to be elliptic (class-A NL shown as red ring in Fig.3c), while $\alpha_{1}\alpha_{2}<0$ would lead the NL to be hyperbolic (class-B NL shown as blue line in Fig.3c). It is easily to prove that the breaking down of the mirror symmetry would only move the nodal link point O and distort the NLs. However, the breaking down of $PT$ symmetry would gap out all nodal points, because $\sigma_{y}$ matrix would be introduced to Eq.2 then. So the nodal link of class-A and class-B NLs are robust under the protection of $PT$ symmetry.

Class-C nodal lines. Actually, such nodal lines are formed by two degenerated bands along the $\Gamma-A$ direction, are protected by $c_{3}$ symmetry, and are shown as a segment between $\Gamma$ and $A$ points in Fig.3a. Along $\Gamma-A$, there is another topological Fermion, namely the triple point Zhu et al. (2016), which is a type-A triple point according to the Ref.Zhu et al. (2016)’s definition, because the Berry phase around $\Gamma-A$ is zero. Such a triple point would evolve into a Dirac point in the presence of SOC [Fig.5b]. However, the topologically induced SSs of triple points would only happen in a system with SOC, and would not happen in the absence of SOC.

Class-D nodal lines. This is an isolated nodal ring at $k_{z}$=0, and is protected by the $\sigma_{h}$ mirror symmetry. It is easily deduced from the $k\cdot p$ model (Eq.77). The eigenvalues of Eq.77 are $\epsilon_{1}^{\pm}=E_{1}+(A\pm C)(k_{x}^{2}+k_{y}^{2})$ and $\epsilon_{3}^{\pm}=E_{3}+(B_{3}\pm G)(k_{x}^{2}+k_{y}^{2})$. since $E_{1}<0$ and $E_{3}>0$, but $(A\pm C)>0$ and $(B_{3}+G)*(B_{3}-G)<0$, $\epsilon_{1}^{\pm}$ would have a cross-point with $\epsilon_{3}^{\pm}$ due to its typical band inversion. The resulting nodal line is shown as an orange circle at the top plane of Fig.3a.

Topological number for the nodal net. The topological number of DNL is characterized by a quantized $\mathbb{Z}_{2}$ topological charge $\nu$ Kim et al. (2015), which is given by the parity of the Berry phase along a loop $S$ that interlinks with the Dirac ring [red loop in Fig.3a]. It is verified that $\nu$ is 1 for red loops $S^{1}$ and $S^{2}$. The topological charge $\nu$ of red loops $S^{1}$ is identical to the topological charge of the green circle in Fig.3a which is composed of lines H2-H1, H1-L1, L1-L2, and L2-H2. Due to the mirror symmetry, the summation of Berry phases along H1-L1 and L2-H2 is zero. We could define a topological charge $\nu_{H}$ and $\nu_{L}$ for H2-H1 and L1-L2 respectively since they form a closed loop in k-space. From previous studies Kim et al. (2015), The topological number $\nu$ for a time reversal invariant loop which links two parity-invariant momenta $\Gamma_{a}$ and $\Gamma_{b}$ is related to the parity of $\Gamma_{a}$ and $\Gamma_{b}$ as

[TABLE]

where $\xi_{n}(\Gamma_{a})$ is the parity for the occupied bands.

So $(-1)^{\nu_{L}}=\xi_{M_{1}}\xi_{L_{1}}$. Since closed loop L1-L2 links $M_{1}$ and $L_{1}$ which are the parity-invariant momenta, it is verified by DFT calculations that $\xi_{M_{1}}=1$ and $\xi_{L_{1}}=1$ (see details in Appendix.E, i.e. $\nu_{L}=0$. So the topological number for H2-H1 is $\nu_{H}=(\nu_{S^{1}}-\nu_{L})\ \text{mod}\ 2=1$. Eventually, the topological numbers are $\nu=0$ and $1$ for regions 1 and 2 shown in Fig.3b with different colors respectively. There will be odd number of nodal lines between regions 1 and 2 due to the topological number change from 0 to 1.

IV Drumhead Surface States

Based on previous studies Kim et al. (2015); Chan et al. (2016); Kobayashi et al. (2017), the 1D $\mathbb{Z}_{2}$ invariant $\nu$ partially guarantees the presence of drumhead surface states. In this section the drumhead surface states on B-terminated and Ti-terminated (001) cleavage surfaces of TiB2 are studied.

B-terminated surface structure is shown in Fig.4d, which is a honeycomb lattice like graphene. By using WannierTools Wu et al. (2017) and the method of iterative Green’s function Sancho et al. (1985) solution based on a tight-binding model (TBM), the surface state spectrum was calculated, as shown in Fig.4a-c. In Fig.4a, it is shown that between $\bar{M}-\bar{K}$ there is a drumhead surface state coming from a ”Dirac” point which is the projection of the nodal line. In region 1, those drumhead surface states form a linear Dirac cone at $\bar{K}$ point, which is analogous to the Dirac cone in Graphene. SSs obtained from the TB model usually are used to explain the topological properties. To compare the SS with ARPES experimental data for further use, we use a first-principle calculations for a slab system due to that fact that a real surface system would have charge reconstructions which cannot be described by TBM. So we simulated a 20-layer slab of TiB2 with VASP. The BS is shown in Fig.4f, in which the color denotes the contributions from the $p_{z}$ orbital of the surface’s boron atoms. Basically, the DFT results are close to the TBM results. The difference is that the zero-energy point of the surface Dirac cone of DFT results is about 0.45 eV higher than that of the TBM results. Since the lattice and the orbital of the boron surface are the same as that of graphene, the effective $k\cdot p$ models at $\bar{K}$ point are the same, and are given as

[TABLE]

where $V_{F}$ is the Fermi velocity. By fitted to the DFT calculations, $V_{F}$ is estimated to $1.5$\mathrm{eV}\text{,}\mathrm{\SIUnitSymbolAngstrom}$$, which is about $2.28\times 10^{5}m/s$ in SI units.

Ti-terminated surface structure is shown in Fig.4e, which is a hexagonal lattice of the titanium atom. The DFT calculated BS of a 20-layer slab of TiB2 with titanium in the outer surface is shown in Fig.4g, in which the color denotes contributions from $d_{zx},d_{yz}$ orbitals of the surface’s titanium atoms. There are drumhead SSs coming from the nodal line in this case as same as B-terminated surface; however, the dispersion of SS at $\bar{K}$ is different from the B-terminated surface. In Fig.4m, it seems that there is a 2D quadratic dispersed Dirac cone at $\bar{K}$ point. However, by zooming into the dispersion at $\bar{K}$ point, it turns out that it is a faked quadratic Dirac cone. Not only there is a linearly dispersed Dirac cone at $\bar{K}$ point, but there is also a linearly dispersed Dirac cone along $\bar{K}-\bar{\Gamma}$. Regarding the point group $C_{3v}$ describing the surface structure, the effective $k\cdot p$ model is obtained as

[TABLE]

where $k_{\pm}=k_{x}\pm ik_{y}$ and $k^{2}=k_{x}^{2}+k_{y}^{2}$. When the quadratic terms $B$ and $C$ are missing, then Eq.(8) reduces to Eq.(5) which leads to massless linear Dirac dispersion $E(k)=\pm|k|$. In the absence of linear term $A$, Eq.(8) describes a Dirac point with parabolic energy dispersion. Such a quadratic Dirac point is unstable due to its topological number $\nu$ being trivial. It will split into four linear Dirac points with nontrivial topological number $\nu=1$ under the $C_{3}$ symmetry  Heikkila and Volovik (2015) which leads to the linear term in Eq.(8) . Eventually, one is centered at $\bar{K}$ point, the other three are connected by the $C_{3}$ symmetry (see Fig.4j). Fitting to the DFT calculated BS, we obtain the following parameters for a Ti-terminated surface: $A=0.08eV\cdot\AA,B=-0.46eV\cdot\AA^{2}$, and $C=2.3eV\cdot\AA^{2}$, in which the linear term is much weaker than the quadratic term.

The surface Dirac cone at $\bar{K}$ of TiB2 with B-terminated and Ti-terminated surfaces are very similar to monolayer and bilayer graphene, which have linearly dispersed and quadratically dispersed Dirac cones respectively. The Berry phase around the linear Dirac cone is $\pi$, while the Berry phase around the quadratically dispersed Dirac cone is $2\pi$. Such differences in the Berry phase would lead to different quantum oscillations. In Ref. Zhang et al. (2014), it was proposed that a Dirac point could be observed in monolayer TiB2; however, the Dirac point proposed in this paper could be observed on the surface of a thick slab system rather than a monolayer system.

V Conclusion

In this paper, based on first-principles calculations and model analysis, a novel $PT$ symmetry protected Dirac nodal net state is recognized in AlB2-type TiB2 and ZrB2 in the absence of SOC. This complex nodal net structure is composed with four classes of NLs: A, B, C, and D, in which, class-A and class-B NLs link together at O along the $\Gamma-K$ direction, three class-B NLs in the vertical mirror planes terminate at A point, which is also a termination of the class-C NL. Several $k\cdot p$ models for these four classes of NLs are constructed under the constraint of their symmetry which confirmed the formation of this nodal net. The topological numbers $\nu$ for different regions in BZ are calculated. It is noted that there are two different dispersed drumhead-like Dirac cones emerging on B-terminated and Ti-terminated surfaces, which are analogous to those of monolayer and bilayer graphene, indicating some novel surface transport properties in TiB2 and ZrB2. We believe that this work will guide further progress in understanding the novel properties of TiB2 and ZrB2, and two different terminated surfaces are good platforms to study 2D Dirac fermions. In addition, AlB2-type TiB2 and ZrB2 can be easily synthesized and provides two prototype materials to study the topological nodal net structure.

Note added.- Recently, Ref. Zhang et al. (2017) appeared, discussing some of the topological properties of metal-diboride where the nodal line around $K$ is the Class-A nodal line in the nodal net of this work.

VI Acknowledge

We acknowledged helpful discussions with X.Dai. X.F and B.W were supported by the National Natural Science Foundation of China (NSFC-51372215), C.Y and Z.S were supported by National Natural Science Foundation of China, the National 973 program of China (Grant No. 2013CB921700), Q.W was supported by Microsoft Research, and the Swiss National Science Foundation through the National Competence Centers in Research MARVEL and QSIT.

Appendix A Computational methods

In this work, the electronic properties for AlB2-type TiB2 are studied by using density functional theory (DFT) Hohenberg and Kohn (1964); Kohn and Sham (1965) as implemented in the Vienna Ab initio Simulation Package (VASP) Kresse and Furthmüller (1996, 1996); Kresse and Joubert (1999). The exchange correlation functional of Perdew-Burke-Emzerhof generalized gradient approximation (GGA-PBE) Blöchl (1994); Perdew et al. (1996) are performed. The standard version of PBE pseudo-potential is adopted in this work explicitly treating four valence electrons for the Ti atom (4d35s1) and three valence electrons for the B atom (3s23p1). A cutoff energy of 500 eV and an 11$\times$11$\times$9 k-mesh are used to perform the bulk calculation. The conjugate-gradient algorithm is used to relax the ions, and the convergence thresholds for total energy and ionic force component are chosen as 1$\times$10*-7* eV and 0.001 eV/$\r{A}$.

For the slab calculations (Fig.4f-g), the thickness of the B-terminated slab is 20 layers of titanium and 21 layers of boron, while the thickness of the Ti-terminated slab is 20 layers of titanium and 19 layers of toron. A $26\times 26\times 1$ $\Gamma$ centered k mesh and a 14$\AA$-thick vacuum are used in the DFT simulations. The surface are fully relaxed with energy convergence up to 1$\times$10*-7* eV and force up to 0.001 eV/$\r{A}$.

The nodal-net searching and surface states spectrum calculations shown in Fig.4a-c are done using the open-source software WannierTools Wu et al. (2017) which is based on Wannier tight binding model (WTBM) constructed with Wannier90 Mostofi et al. (2008). Ti $s,p,d$, and B $s,p$ orbitals are used as initial projectors for WTBM construction. WTBMs constructed with Wannier90 do not exactly fulfil all crystal symmetries which is very important for nodal points searching because usually nodal points are protected by crystal symmetries except Weyl points. The WTBM is symmetrized to be compatible with the crystal symmetry using the method described in Red. Zhu et al. (2016).

Appendix B band structure and Fermi surface of AlB2-type ZrB2

In general, ZrB2 is often compared with TiB2. The structure of AlB2-type ZrB2 also has hexagonal structure with a space group of P6/mmm (No.191). Its optimized lattice constants are a=b=3.1748(3) $\r{A}$ and c=3.5579(7) $\r{A}$, which is slightly larger than those of TiB2. As shown in Fig.6a, there are also six band crossing points along the high-symmetry path similar with TiB2. Among them, Zr 4d states are the main contribution orbitals for these six band crossing points. The Zr-dxz and Zr-dxy orbitals are much higher than the Fermi level. From Figs. 6b and c, the Fermi surface of ZrB2 is sightly different from that of TiB2. The circling surface around A point is disappeared in the Fermi surface of ZrB2. The main part of the lantern-like Fermi surface of ZrB2 is basically consistent with the Fermi surface of TiB2, which also shows a nodal-net feature.

Appendix C Nodal net structure of TiB2 without vertical mirror planes

Since the nodal net of TiB2 is $PT$ symmetry protected, any perturbations that preserve the inversion symmetry can only distort the nodal-net but not destroy them. As an example,we apply a uniaxial strain (compression 1$\%$) along the [100] crystal direction. After deformation, the structure of AlB2-type TiB2 belongs to space group of C2/m (No. 12), which just preserves the Mz mirror reflection symmetry. As a result, as shown in Fig.7, it is found that the class-A nodal line still embeds in the $k_{z}=0$ plane, one of the class-B nodal line become isolated with the other two class-B nodal lines which link with the class-D nodal line that embeds in the $k_{z}=0.5$ plane, while the class-D nodal line along $\Gamma-A$ disappears because it is protected by $C_{3}$ symmetry which is destroyed under such strain. Thus, the nodal net in TiB2 is robustly stable, which does not requires the protection of mirror symmetry.

Appendix D Effective $k\cdot p$ model

We derive several $k\cdot p$ models describing the bulk bands in the vicinity of $\Gamma$, A, and K in the 3D BZ, and a $k\cdot p$ for surface states at $\bar{K}$ point in 2D BZ. The $k\cdot p$ models are used to get a better understanding of the surface states and the nodal net structures, and would be useful for further investigations of Landau level and quantum transport properties.

The $k\cdot p$ models were calculated using the kdotp_symmetry code, which implements the method described in Ref. Gresch et al. (2017), basically, there are two things that should be prepared before applying kdotp_symmetry. Firstly, we identify the little group $G$ of a high symmetry point $K_{0}$ of which we want to construct a low energy effective model, and get the generators $R$ of little group $G$. Secondly, we identify the representations of $R$ on the basis of the eigenvectors at the selected high-symmetry points. Then kdotp_symmetry will produce the $k\cdot p$ model under the constraint

[TABLE]

where D(R) is the representation matrix of symmetry operator $R$.

D.1 $k\cdot p$ model at $\Gamma$ point

The little point group at Gamma point of bulk TiB2 is $D_{6h}$ plus TR symmetry (see Table 76 of Ref.Koster (1962)). There are three generators of $D_{6h}$ including spatial-inversion $I(-x,-y,-z)$, 2-fold rotation $c_{2y}(-x,y,-z)$ and 6-fold rotation $c_{6z}(\frac{x}{2}+\frac{\sqrt{3}}{2}y,-\frac{\sqrt{3}}{2}x+\frac{y}{2},z)$. From the fat-band analysis of Fig.2, the relevant bands come from $d_{z^{2}}$ which belongs to the $A_{1g}$ representation, $d_{zx}$ and $d_{yz}$ orbitals which form the basis of its $E_{1g}$ representation. The representations of group generators according to the symmetrical basis {($d_{zx}$+i$d_{yz}$)/$\sqrt{2}$, ($d_{zx}$-i$d_{yz}$)/$\sqrt{2}$ , $d_{z^{2}}$} are given by

[TABLE]

where, $K$ is the complex conjugation operator. Considering these symmetries, and the constraint of Eq.(9), a 3-band model up to second order of k around $\Gamma$ point for bulk TiB2 is given by

[TABLE]

where $\varepsilon_{1}({\bf k})=E_{1}+A_{1}(k_{x}^{2}+k_{y}^{2})+B_{1}k_{z}^{2}$, $\varepsilon_{2}({\bf k})=E_{2}+A_{2}(k_{x}^{2}+k_{y}^{2})+B_{2}k_{z}^{2}$.

As mentioned in the main text, Eq.(22) with only second order momentum $k$ would lead to a nodal surface other than the nodal line structure. To distinguish the difference between $\Gamma-K$ and $\Gamma-M$ directions, we have to introduce the sixth order of $k_{x}$, $k_{y}$ in the $k_{z}=0$ plane, and introduce a fourth order of $k_{x}$, $k_{y}$ in the off-diagonal part, so eventually, the new $k\cdot p$ model is given by

[TABLE]

where $\varepsilon_{1}({\bf k})=E_{1}+A_{1}(k_{x}^{2}+k_{y}^{2})+B_{1}k_{z}^{2}$, $\varepsilon_{2}({\bf k})=E_{2}+A_{2}(k_{x}^{2}+k_{y}^{2})+B_{2}k_{z}^{2}+L(k_{x}^{2}+k_{y}^{2})^{2}+M(k_{+}^{6}+k_{-}^{6})$.

By fitting to the DFT band structure of TiB2, the parameters in Eq.(26) are obtained: $E_{1}=1.787$\mathrm{eV}, $B_{1}=-3.8$\mathrm{eV}\text{\,}\mathrm{\SIUnitSymbolAngstrom}$^{2}$, $A_{1}=2.6$\mathrm{eV}\text{\,}\mathrm{\SIUnitSymbolAngstrom}$^{2}$, $E_{2}=-2.12$\mathrm{eV}, $A_{2}=1.63$\mathrm{eV}\text{,}\mathrm{\SIUnitSymbolAngstrom}$^{2}$, $B_{2}=5.1$\mathrm{eV}\text{,}\mathrm{\SIUnitSymbolAngstrom}$^{2}$, $L=1.3$\mathrm{eV}\text{,}\mathrm{\SIUnitSymbolAngstrom}$^{4}$, $C=3.55$\mathrm{eV}\text{,}\mathrm{\SIUnitSymbolAngstrom}$^{2}$, $M=0.65$\mathrm{eV}\text{,}\mathrm{\SIUnitSymbolAngstrom}$^{6}$, $F=1.83$\mathrm{eV}\text{,}\mathrm{\SIUnitSymbolAngstrom}$^{4}$, and $D=5.1$\mathrm{eV}\text{,}\mathrm{\SIUnitSymbolAngstrom}$^{2}$. The comparison between DFT bands and $k\cdot p$ bands is shown in Fig.8

The nodal line structure around $\Gamma$ point calculated from Eq.26 with the DFT fitted parameters is shown in Fig.9a. It is shown that the nodal net close to $\Gamma$ point is very similar to Fig.3, the class-A, class-B, and class-C nodal lines are captured successfully; however the nodal line beyond the nexus point $A$ is not captured. That is because the nexus point $A$ is at the boundary of the BZ, which is related to an infinity in the $k\cdot p$ model. The position of the nexus point $A$ could be tuned by changing $F$ in Eq.26. In Fig.9c, it is shown that the nexus point would disappear if $F$ was very large. The nodal line in the $\sigma_{d}$ mirror plane could be shown if the fourth and sixth order terms in Eq.26 become smaller (Fig.9b).

D.2 $k\cdot p$ model at K point

The little group at K point in the TiB2 is $D_{3h}$ (see Table 65 of Ref.Koster (1962)), of which there are three generators including horizontal mirror $\sigma_{h}(x,y,-z)$, 2-fold rotation $c_{2x}(x,-y,-z)$ and 3-fold rotation $c_{3z}(-\frac{x}{2}-\frac{\sqrt{3}}{2}y,\frac{\sqrt{3}}{2}x-\frac{y}{2},z)$. Around K point, the relevant representations are $\Gamma_{5},and\Gamma_{6}$, of which the basis are $\{d_{zx},d_{yz}\}$ and $\{d_{xy},d_{x^{2}-y^{2}}\}$. Taking the symmetrical orbitals { $|Y_{2}^{1}\rangle$, -$|Y_{2}^{-1}\rangle$} and { $-|Y_{2}^{2}\rangle,|Y_{2}^{-2}\rangle$} as a basis, where $|Y_{l}^{m}\rangle$ is the complex Spherical harmonic function, the representations of the generators are given by

[TABLE]

Considering these symmetries above and the constraint of Eq.(9), a 4-band model up to second order of k around ${K}$ point in bulk TiB2 is given by

[TABLE]

where $\varepsilon(k)=E_{0}+A(k_{x}^{2}+k_{y}^{2})+Bk_{z}^{2}$ and $\varepsilon^{\prime}(k)=E_{0}^{\prime}+A^{\prime}(k_{x}^{2}+k_{y}^{2})+B^{\prime}k_{z}^{2}$. The fitted parameters of TiB2 are $E_{0}=-1.8651$\mathrm{eV}, $E_{0}^{\prime}=1.3387$\mathrm{eV}, $A=2.55$\mathrm{eV}\text{,}\mathrm{\SIUnitSymbolAngstrom}$^{2}$, $A^{\prime}=-16.5$\mathrm{eV}\text{,}\mathrm{\SIUnitSymbolAngstrom}$^{2}$, $C=3.7$\mathrm{eV}\text{,}\mathrm{\SIUnitSymbolAngstrom}, $C^{\prime}=-1.66$\mathrm{eV}\text{\,}\mathrm{\SIUnitSymbolAngstrom}, $D=-0.58$\mathrm{eV}\text{,}\mathrm{\SIUnitSymbolAngstrom}$^{2}$, $D^{\prime}=18.8$\mathrm{eV}\text{,}\mathrm{\SIUnitSymbolAngstrom}$^{2}$, $B=8.73$\mathrm{eV}\text{,}\mathrm{\SIUnitSymbolAngstrom}$^{2}$, $B^{\prime}=-6.37$\mathrm{eV}\text{,}\mathrm{\SIUnitSymbolAngstrom}$^{2}$, $E=24.9$\mathrm{eV}\text{,}\mathrm{\SIUnitSymbolAngstrom}$^{2}$, $F=6.84$\mathrm{eV}\text{,}\mathrm{\SIUnitSymbolAngstrom}$$. The comparison between DFT bands and $k\cdot p$ bands is shown in Fig. 13.

In order to analyse the nodal line structure, Eq.(42) can be written into two blocks

[TABLE]

where

[TABLE]

With the fitted parameters, we find that our $k\cdot p$ model not only describes the nodal line surrounding $K$ point, but can also predict part of the nodal line in the vertical mirror plane.

D.3 $k\cdot p$ model at A point

For the little group at A $D_{6h}$ (see Table 76 of Ref.Koster (1962)), relevant representations are $A_{1g}$ $E_{1g}$ and $E_{2g}$, which constitutes {$d_{z^{2}}$}, {$d_{zx}$, $d_{yz}$} and{$d_{xy}$, $d_{x^{2}-y^{2}}$}. By symmetrization of those orbitals according to the $D_{6h}$ group, the symmetrical basis are chosen as $|Y_{2}^{0}\rangle$, $|Y_{2}^{1}\rangle$, -$|Y_{2}^{-1}\rangle$, -$|Y_{2}^{2}\rangle$ and $|Y_{2}^{-2}\rangle$, and the related representations of its generators are given by

[TABLE]

The constructed $k\cdot p$ model is given by

[TABLE]

where $\varepsilon_{i}(k)=E_{i}+A_{i}(k_{x}^{2}+k_{y}^{2})+B_{i}k_{z}^{2}$ with $i=1,2,3$. The fitted parameters are $E_{1}=-0.2426$\mathrm{eV}, $A_{1}=170$\mathrm{eV}\text{\,}\mathrm{\SIUnitSymbolAngstrom}$^{2}$, $B_{1}=3.53$\mathrm{eV}\text{\,}\mathrm{\SIUnitSymbolAngstrom}$^{2}$, $E_{2}=1.575$\mathrm{eV}, $A_{2}=1.27$\mathrm{eV}\text{,}\mathrm{\SIUnitSymbolAngstrom}$^{2}$, $B_{2}=-6.16$\mathrm{eV}\text{,}\mathrm{\SIUnitSymbolAngstrom}$^{2}$, $E_{3}=0.3822$\mathrm{eV}$$, $A_{3}=-2.39$\mathrm{eV}\text{,}\mathrm{\SIUnitSymbolAngstrom}$^{2}$ $B_{3}=26.55$\mathrm{eV}\text{,}\mathrm{\SIUnitSymbolAngstrom}$^{2}$, $C=0.47$\mathrm{eV}\text{,}\mathrm{\SIUnitSymbolAngstrom}$^{2}$, $D=0.0$, $E=0.0$, $F=0.1$\mathrm{eV}\text{,}\mathrm{\SIUnitSymbolAngstrom}$^{2}$, $G=5.62$\mathrm{eV}\text{,}\mathrm{\SIUnitSymbolAngstrom}$^{2}$. The fitted band is shown in Fig.12.

Particularly, when the effective $k_{z}=0$, this means that bulk $k_{z}=0.5$, Eq.(72) decouples into three block diagonal matrices. From Fig.2, it is shown that along $A-H$, band crossing happens between $d_{z}^{2}$ and the combination of $d_{xz}$ and $d_{yz}$. i.e. we only have to consider the following block of Eq.(72)

[TABLE]

D.4 $k\cdot p$ model for surface states at $\bar{K}$ points

The little group at $\bar{K}$ point of slab system TiB2 is $C_{3v}$, which has two generators $c_{3z}(-\frac{x}{2}-\frac{\sqrt{3}}{2}y,\frac{\sqrt{3}}{2}x-\frac{y}{2},z)$ and $\sigma_{v}(x,-y,z)$. According to the DFT calculations, it was determined that the surface state at $\bar{K}$ belongs to $E$ (see Table 49 of Ref.Koster (1962)) representation of $C_{3v}$. On the basis of complex orbitals, the related representations of its generators are given by

[TABLE]

Considering these symmetries above and the constraint of Eq.(9), a 2-band model up to the second order of k around $\bar{K}$ point for surface states is given by

[TABLE]

where $k_{\pm}=k_{x}\pm ik_{y}$ and $k^{2}=k_{x}^{2}+k_{y}^{2}$. The linear part of Eq.(85) leads to massless Dirac dispersion $E(k)=\pm|k|$. The combination of the linear term and the quadratic term leads to 3-fold rotation symmetry of the energy dispersion. Fitting to DFT calculation band structure, We obtain the following parameters for Ti-terminated surface: $A=0.08$\mathrm{eV}\text{,}\mathrm{\SIUnitSymbolAngstrom}$,B=-0.46$\mathrm{eV}\text{,}\mathrm{\SIUnitSymbolAngstrom}$^{2},C=2.3$\mathrm{eV}\text{,}\mathrm{\SIUnitSymbolAngstrom}$^{2}$, and for the B-terminated surface A=1.5$\mathrm{eV}\text{\,}\mathrm{\SIUnitSymbolAngstrom}$, B=0, C=0.

Appendix E Parities at TRIMs

The parities of the occupied bands of TiB2 at TRIMs are listed in Table.1. It is noted that there are six occupied bands at $A$, and five occupied bands at other TRIMs. The product of the parities of the occupied bands at $M$ and $L$ are 1, which leads to $\xi_{M}=1$ and $\xi_{L}=1$.