Dual topological nodal line and nonsymmorphic Dirac semimetal in three dimensions

TL;DR

This paper introduces a new three-dimensional Dirac semimetal that combines topological and nonsymmorphic features, demonstrated through a minimal model and first-principles calculations, revealing rich topological phenomena.

Contribution

The study presents a novel 3D Dirac semimetal with combined topological and nonsymmorphic properties, supported by a minimal tight-binding model and material predictions.

Findings

Discovery of a 3D Dirac semimetal with both features

Identification of symmetry-enforced Dirac points and Weyl nodal lines

Prediction of real materials exhibiting this semimetal

Abstract

Previously known three-dimensional Dirac semimetals (DSs) occur in two types -- topological DSs and nonsymmorphic DSs. Here we present a novel three-dimensional DS that exhibits both features of the topological and nonsymmorphic DSs. We introduce a minimal tight-binding model for the space group 100 that describes a layered crystal made of two-dimensional planes in the $p4g$ wallpaper group. Using this model, we demonstrate that double glide-mirrors allow a noncentrosymmetric three-dimensional DS that hosts both symmetry-enforced Dirac points at time-reversal invariant momenta and twofold-degenerate Weyl nodal lines on a glide-mirror-invariant plane in momentum space. The proposed DS allows for rich topological physics manifested in both topological surface states and topological phase diagrams, which we discuss in detail. We also perform first-principles calculations to predict that the proposed DS is realized in a set of existing materials BaLa$X$B$Y_5$, where $X$ = Cu or Au, and $Y$ = O, S, or Se.

Tables

Table 1. Table 1: Lattice constants in a and c direction for unit cell of the BaLaXBY 5 family.

	BaLaCuBO5	BaLaCuBS5	BaLaCuBSe5	BaLaAuBO5	BaLaAuBS5	BaLaAuBSe5
a	5.4769 Å	6.5298 Å	6.8915 Å	5.7041 Å	6.6819 Å	7.0001 Å
c	7.4640 Å	8.6036 Å	8.9960 Å	7.7112 Å	8.6131 Å	9.0063 Å

Table 2. Table 2: Classification and representations of symmetry-lowering perturbations of SG 100 100 100 . The unit cell is illustrated by a box. Different color schemes are used for the atoms to represent different on-site energy.

Class	Strain	Perturbation	Band gap open
			$\Gamma$	$X$	$Y$	$M$
$A_{2u}$		${\tau }_z$	O	O	O	O
		$\sin \frac{k_x}{2}\cos \frac{k_y}{2}{\tau }_x{\sigma }_x+\cos \frac{k_x}{2}\sin \frac{k_y}{2}{\tau }_x{\sigma }_y$	$\times$	O	O	$\times$
$B_{1u}$		$\sin \frac{k_x}{2}\cos \frac{k_y}{2}{\tau }_x{\sigma }_y-\cos \frac{k_x}{2}\sin \frac{k_y}{2}{\tau }_x{\sigma }_x$	$\times$	O	O	$\times$
$B_{2g}$		$\sin \frac{k_x}{2}\sin \frac{k_y}{2}{\tau }_x$	$\times$	$\times$	$\times$	O
$B_{2u}$		$\sin \frac{k_x}{2}\cos \frac{k_y}{2}{\tau }_x{\sigma }_x-\cos \frac{k_x}{2}\sin \frac{k_y}{2}{\tau }_x{\sigma }_y$	$0$	${\tau }_x{\sigma }_x$	$-{\tau }_x{\sigma }_y$	${\tau }_y{\sigma }_z$
$E_g$		$\sin \frac{k_x}{2}\cos \frac{k_y}{2}{\tau }_y$	$\times$	O	$\times$	$\times$
		$\cos \frac{k_x}{2}\sin \frac{k_y}{2}{\tau }_y$	$\times$	$\times$	O	$\times$
		$\cos \frac{k_x}{2}\cos \frac{k_y}{2}{\tau }_y{\sigma }_x$	O	$\times$	$\times$	$\times$
		$\sin \frac{k_x}{2}\sin \frac{k_y}{2}{\tau }_y{\sigma }_y$	$\times$	$\times$	$\times$	O
$E_u$		$\sin \frac{k_x}{2}\cos \frac{k_y}{2}{\tau }_x{\sigma }_z$	$\times$	O	$\times$	$\times$
		$\cos \frac{k_x}{2}\cos \frac{k_y}{2}{\tau }_y{\sigma }_y$	O	$\times$	$\times$	$\times$
		$\sin \frac{k_x}{2}\sin \frac{k_y}{2}{\tau }_y{\sigma }_x$	$\times$	$\times$	$\times$	O
		$\cos \frac{k_x}{2}\sin \frac{k_y}{2}{\tau }_x{\sigma }_z$	$\times$	O	$\times$	$\times$

Equations

$${\mathcal{H}}^0(\mathbf{k})={\mathcal{H}}^t(\mathbf{k})+V(\mathbf{k}),$$

$${\mathcal{H}}^t(\mathbf{k})=$$

$$+$$

$$V(\mathbf{k})=$$

$$+$$

$${\mathcal{H}}^0(\mathbf{k}+\mathbf{G})=e^{-i\mathbf{d}({\tau }_z)\cdot \mathbf{G}}{\mathcal{H}}^0(\mathbf{k})e^{i\mathbf{G}\cdot \mathbf{d}({\tau }_z)}.$$

$${\mathcal{H}}^1(\mathbf{k})=$$

$$\displaystyle+m_{\scalebox{0.55}{$A_{\scalebox{1.0}{$2u$}}$}}\tau_{z}.$$

$${\mathcal{H}}_{+}^0=$$

$$+v_3\sin k_z{\tau }_z{\sigma }_z,$$

$${\mathcal{H}}_{-}^0=$$

$$+$$

$$\mathcal{W}(k_l,{\mathbf{k}}_0)=F({\mathbf{k}}_0)F({\mathbf{k}}_1)\cdots F({\mathbf{k}}_{N-2})F({\mathbf{k}}_{N-1}),$$

$${\left[\mathcal{W}(k_l,{\mathbf{k}}_0)\right]}_{nm}\equiv P{\left[\exp \left(i{\oint }_{k_l}d\mathbf{k}\cdot \mathbf{A}(\mathbf{k})\right)\right]}_{nm},$$

$$\mathcal{W}(k_y)=F(k_{x,0},k_y)F(k_{x,1},k_y)\cdots F(k_{x,N-2},k_y)F(k_{x,N-1},k_y),$$

$$\mathbf{k}=k_0\left(\cos \phi \sin \theta ,\sin \phi \sin \theta ,\cos \theta \right),$$

$$\mathcal{W}(\theta )=F(\theta ,{\phi }_0)F(\theta ,{\phi }_1)\cdots F(\theta ,{\phi }_{N-2})F(\theta ,{\phi }_{N-1}),$$

$${\mathcal{H}}^0(\mathbf{k})=$$

$$+v_0\cos \frac{k_x}{2}\cos \frac{k_y}{2}{\tau }_y{\sigma }_z+v_1\left(\sin k_x{\sigma }_x+\sin k_y{\sigma }_y\right){\tau }_z+v_2\left(\sin k_x{\sigma }_y-\sin k_y{\sigma }_x\right)+v_3\sin k_z{\tau }_z{\sigma }_z$$

$$+v_4\left(\sin \frac{k_x}{2}\cos \frac{k_y}{2}{\sigma }_y-\cos \frac{k_x}{2}\sin \frac{k_y}{2}{\sigma }_x\right){\tau }_x.$$

$$\mathcal{T}=i{\sigma }_y\mathcal{K},g_x={\tau }_y{\sigma }_x,C_{4z}={\tau }_z\exp \left[i\frac{\pi }{4}{\sigma }_z\right].$$

$${\mathcal{H}}^M(\mathbf{k})=k_x\left(c_1{\tau }_z{\sigma }_x+c_2{\tau }_x{\sigma }_x+c_3{\sigma }_y\right)+k_y\left(c_1{\tau }_z{\sigma }_y-c_2{\tau }_x{\sigma }_y-c_3{\sigma }_x\right)+k_z\left(c_4{\tau }_z{\sigma }_z+c_5{\tau }_x{\sigma }_z\right).$$

$${\mathcal{T}}^{-1}{\Gamma }_a\mathcal{T}=-{\Gamma }_a,$$

$${\mathcal{T}}^{-1}{\Gamma }_5\mathcal{T}={\Gamma }_5.$$

$${\mathcal{H}}^M(\mathbf{k})=k_x\left(c_1{\Gamma }_{23}+c_2{\Gamma }_2-c_3{\Gamma }_{24}\right)+k_y\left(c_1{\Gamma }_{13}-c_2{\Gamma }_1-c_3{\Gamma }_{14}\right)+k_z\left(-c_4{\Gamma }_{34}+c_5{\Gamma }_4\right).$$

$$H(k_x,0,0):\pm {\left(c_1^2+c_2^2+c_3^2\right)}^{1/2},$$

$$H(k_x,k_y,0):\pm {\left[{\left\{c_2{k}_x\pm {\left(c_1^2+c_3^2\right)}^{1/2}{k}_y\right\}}^2+{\left\{c_2{k}_y\pm {\left(c_1^2+c_3^2\right)}^{1/2}{k}_x\right\}}^2\right]}^{1/2},$$

$$H(k_x,0,k_z):\pm {\left\{c_3^2{k}_x^2+{\left(c_2{k}_x\pm c_4{k}_z\right)}^2+{\left(c_1{k}_x\mp c_5{k}_z\right)}^2\right\}}^{1/2},$$

$$H(k_x,k_y,k_z):$$

$$\pm \left[\left(c_1^2+c_2^2+c_3^2\right)\left(k_x^2+k_y^2\right)+\left(c_1^2+c_5^2\right)k_z^2\pm 2{\left\{4c_2^2\left(c_1^2+c_3^2\right)k_x^2{k}_y^2+\left[\left(c_2{c}_4\mp c_1{c}_5\right)k_x^2+\left(c_2{c}_4\pm c_1{c}_5\right)k_y^2\right]k_z^2\right\}}^{1/2}\right]$$

$$\displaystyle H^{1}({\bm{k}})=m_{\scalebox{0.55}{$A_{\scalebox{1.0}{$2u$}}$}}\tau_{x}+m_{\scalebox{0.55}{$E_{\scalebox{1.0}{$g$}}$}}^{(1)}\sin\frac{k_{x}}{2}\cos\frac{k_{y}}{2}\tau_{y}+m_{\scalebox{0.55}{$E_{\scalebox{1.0}{$g$}}$}}^{(2)}\cos\frac{k_{x}}{2}\sin\frac{k_{y}}{2}\tau_{y}+m_{\scalebox{0.55}{$B_{\scalebox{1.0}{$2g$}}$}}\sin\frac{k_{x}}{2}\sin\frac{k_{y}}{2}.$$

$${\mathcal{H}}_{+}^0(\mathbf{k})=t_1\cos \frac{k_x}{2}\cos \frac{k_y}{2}{\tau }_x+v_1\left(\sin k_x{\tau }_z{\sigma }_x+\sin k_y{\tau }_z{\sigma }_y\right)+v_3\sin k_z{\tau }_z{\sigma }_z,$$

$$\displaystyle\mathcal{H}^{0}_{-}({\bm{k}})=v_{-}\Bigg{[}v_{0}\cos\frac{k_{x}}{2}\cos\frac{k_{y}}{2}\tau_{y}\sigma_{z}+v_{2}\left(\sin k_{x}\,\sigma_{y}-\sin k_{y}\,\sigma_{x}\right)+v_{4}\left(\sin\frac{k_{x}}{2}\,\cos\frac{k_{y}}{2}\sigma_{y}-\cos\frac{k_{x}}{2}\sin\frac{k_{y}}{2}\sigma_{x}\right)\tau_{x}\Bigg{]}.$$

$${\nu }_n=\prod _i{\xi }_i,$$

$$\zeta=-{\rm Sign}\left(t_{1}\right)\times{\rm Sign}\left(m_{\scalebox{0.55}{$E_{\scalebox{1.0}{$g$}}$}}^{(1)}\right)\times{\rm Sign}\left(m_{\scalebox{0.55}{$E_{\scalebox{1.0}{$g$}}$}}^{(2)}\right)\times{\rm Sign}\left(m_{\scalebox{0.55}{$B_{\scalebox{1.0}{$2g$}}$}}\right).$$

Full text

Dirac semimetals (DS) refer to a class of topological semimetals, characterized by hosting massless Dirac fermions in momentum space Armitage et al. (2018). First identified in graphene with the vanishingly weak spin-orbit coupling (SOC), the massless Dirac fermion system has attracted a surge of interest, exhibiting exotic properties and potential applications for future electronic devices Geim (2009); Allen et al. (2010). Notably, with the advent of topological insulators Hasan and Kane (2010); Qi and Zhang (2011), the three-dimensional (3D) DS with strong spin-orbit coupling has reinforced their status as an important class of topological semimetals. It was first noted that a 3D DS can occur at the phase boundary between the topological and the normal insulators in the presence of inversion symmetry Fu et al. (2007); Murakami (2008). Later, Young et al. found that the 3D DS can be stabilized by crystalline symmetries and time-reversal symmetry Young et al. (2012), and Wang et al. theoretically proposed the material realizations in Na3Bi and Cd3As2 Wang et al. (2012, 2013), which were confirmed experimentally Liu et al. (2014a); Xu et al. (2013); Borisenko et al. (2014); Neupane et al. (2014); Liu et al. (2014b). Currently, the DSs are expected to exist in a variety of forms, such as two-dimensional (2D) DSs Wieder and Kane (2016), double DSs Wieder et al. (2016), type-@slowromancapii@ DSs Chang et al. (2017), and Dirac-Weyl semimetals Gao et al. (2018).

In spite of this variety, it is surprising to notice that all the previously known DSs fall into two disjoint classes, dubbed topological and nonsymmorphic DSs, respectively Yang and Nagaosa (2014). The nonsymmorphic class of the DSs is characterized by hosting Dirac points (DPs) that are pinned at the time-reversal invariant momenta (TRIMs) of the Brillouin zone (BZ). On the other hand, the topological class of the DSs distinguish themselves from the nonsymmorphic DSs by having a pair of DPs off TRIMs. Another distinguishing feature of the topological DSs is the coexistence of nontrivial band topology in the bulk, manifested as gapless excitations on the surface Kargarian et al. (2016); Bednik (2018). In contrast, the bulk bands of the nonsymmorphic DSs are expected to be topologically trivial. Instead, a topological nature of the nonsymmorphic class is reflected in topological phase transitions, driven by symmetry-lowering perturbations from the nonsymmorphic DS into either a topological insulator or a normal insulator Young et al. (2012); Steinberg et al. (2014); Yang and Nagaosa (2014); Schoop et al. (2016); Yang et al. (2017).

In this paper we provide an exception to this a priori classification of 3D DSs. Developing a minimal tight-binding model for space groups (SGs) $P4bm$ (# 100), we establish the existence of a novel type of 3D DSs, characterized by featuring both the topological nodal lines and nonsymmorphic DSs. It is shown that the DS hosts the DPs that reside at TRIMs, which is a characteristic feature of the nonsymmorphic DSs. Simultaneously, the bulk bands carry nontrivial band topology, giving rise to topological surface states, which is unexpected from the previously known nonsymmorphic Dirac semimetals. A striking consequence of this dual nonsymmorphic and topological nature of the DS is the rich topological physics manifested not only in the surface energy spectrum but also in topological phase transitions driven by symmetry-breaking perturbations. Drumhead-like topological surface states arise due to the nontrivial band topology in the bulk, characterized by hosting Weyl nodal lines (WNLs). Moreover, symmetry-lowering perturbations derive a topological phase transition from the proposed DS to distinct topological phases, including a weak topological insulator (WTI) and Weyl and double Weyl semimetal (WS) phases. Using first-principles calculations, we also discuss its material realization in an existing compound, BaLaCuBO5.

Let us begin with elucidating the role of symmetries in SG 100 to protect degeneracies of the Bloch states. SG 100 has the distinguishing feature that it is generated by a glide-mirror $g_{x}$ and a fourfold rotation $C_{4z}$ without inversion symmetry. As emphasized in Zaheer (2014); Wieder et al. (2018), the double glide-mirrors, $g_{x}$ and $g_{y}=C_{4z}^{-1}g_{x}C_{4z}$, together with time-reversal symmetry $\mathcal{T}$, span four-dimensional irreducible representations (FDIRs) at $M=(\pi,\pi,0)$ and $A=(\pi,\pi,\pi)$, where $g_{x}$ and $g_{y}$ satisfy the minimal algebras for a FDIR, $g_{x}^{2}=g_{y}^{2}=1$ and $[\mathcal{T},g_{x(y)}]=\{g_{x},g_{y}\}=0$. Moreover, the linear dispersion of the bands is generic at $M$ and $A$ since a $\mathcal{T}$-odd vector representation of the point group at the $M$ and $A$ points is present in the tensor product of the FDIRs Young et al. (2012); Zaheer (2014); Wieder et al. (2016). Therefore, the presence of DPs are enforced in SG 100 when the filling is an odd multiple of four. In addition, $g_{x}$ ($g_{y}$) and $\mathcal{T}$ further give rise to a constraint to the connectivity of the bands, such that the Kramers pairs at $\Gamma$ and $Y(X)$ should exchange their partners from $\Gamma$ to $Y(X)$ without opening a band gap, leading to hourglass-like connectivity Young and Kane (2015); Wang et al. (2016); Wieder et al. (2018). As a consequence, it is guaranteed that additional twofold-degenerate WNLs are present on the $k_{x}=0$ ($k_{y}=0$) plane, protected by glide-mirror $g_{x}$ ($g_{y}$).

The above symmetry-analysis provides a guiding principle to design the DS hosted in SG 100. Since SG 100 and the $p4g$ wallpaper group are equivalent, generated by $C_{4z}$ and $g_{x}$, a minimal four-band tight-binding model can be constructed from an infinite stack of the identical layer in the $p4g$ wallpaper group as illustrated in Fig. 1. The constructed lattice model is presented in Fig. 2(a). A unit cell comprises two sublattices $A$ and $B$ (labeled by $\tau_{z}=\pm 1$), which are coordinated at ${\bm{d}}(\tau_{z})=\frac{1}{4}\left[\left(2+\tau_{z}\right){{\bm{a}}}_{x}+\left(2-\tau_{z}\right){{\bm{a}}}_{y}\right]$, respectively. The corresponding tight-binding Hamiltonian is given as

[TABLE]

where

[TABLE]

describes the nearest hopping of electrons, and

[TABLE]

describes the potential terms that lower the transnational symmetry of $\mathcal{H}^{t}({{\bm{k}}})$ into SG 100. $V({\bm{k}})$ is constructed, such that it preserves the generators of SG 100, $g_{x}=\tau_{x}\exp\left(-i\frac{\pi}{2}\sigma_{x}\right)$ and $C_{4z}=\exp\left(-i\frac{\pi}{4}\sigma_{z}\right)$, and time-reversal symmetry $\mathcal{T}=i\sigma_{y}K$, where $\{\sigma_{i}\}_{i=x,y,z}$ are the Pauli matrices for spins. We adopted a gauge, in which the Hamiltonian $\mathcal{H}^{0}({\bm{k}})$ transforms under the translation of a reciprocal lattice vector ${\bm{G}}$ according to

[TABLE]

Figure 2(c) shows the electronic energy bands calculated from the tight-binding model Eq. 1. Without the inversion symmetry, each band along the high-symmetry $M$-$\Gamma$-$X$-$R$ line is non-degenerate, thus forming a fourfold degeneracy at $M$. Note that the bands are linearly dispersing in the vicinity of $M$ point. Therefore, the bands feature a nonsymmorphic DP at $M$. We also have confirmed that an additional DP is present at $A$, as we expected from the symmetry-analysis. Based on the Wilson bands calculations Yu et al. (2011), we find that the DPs carry the zero Chern number 111See the Supplemental Materials for the details of the Wilson bands calculations.. This indicates that the fourfold degeneracy is a genuine DP, in the sense that it is a composite of two WPs with $\pm 1$ Chern numbers, respectively.

Besides the fourfold degeneracy, we find that the bands also feature twofold-degeneracy WNLs, the hourglass-like band connectivity on the high-symmetry $\Gamma$-$X$ ($Z$-$R$) line. This guarantees the presence of a twofold degeneracy on the $k_{y}=0$ plane, as shown in Fig. 2(c). A close inspection in the entire BZ reveals that one-dimensional nodal lines are present in the vicinity of $X$ ($R$) lying on the $k_{y}=0$ plane. We have confirmed that a Weyl line node carries the $\pi$ Berry phase, calculated along a $C_{2}\mathcal{T}$-invariant path that threads the nodal line [See the left panel of Fig. 2(f).], where the Berry phase is $\mathbb{Z}_{2}$-quantized. As a consequence of the $\pi$ Berry phase, drumhead-like states emerge on the surface where the projected interior region of a nodal line has non-zero area, such as the (100) surface. As shown in Fig. 2(d), the slab band calculation results in the topological surface states at $E=0$ near the $\overline{\Gamma}$ ($\overline{R}$), which constitute a part of the drumhead-like surface states on the (010) surface.

The WNL hosted in SG 100 is of a hourglass-type Bzdušek et al. (2016); Wang et al. (2017a, b), which is robust against the band inversion at $M$ ($R$). As illustrated in Fig. 2(e), the band inversion at $X$ shrinks the size of the nodal line into a fourfold-degenerate DP. However, instead of annihilating it, the band inversion reverts the DP to a nodal line due to the hourglass-like band connectivity. We assert that this type of WNLs can be characterized by a non-trivial $\mathcal{Z}_{2}$ topological invariant, calculated on the time-reversal invariant sphere that encloses a nodal line [See the right panel of Fig. 2(f).]. The Wilson bands calculation results in the same connectivity of the Wilson bands as those of a 3D Dirac point 222See the Supplemental Materials for the detailed results of the topological invariants.. The nontrivial $\mathcal{Z}_{2}$ invariant, again, reveals that the WNL can be shrunk to form a 3D DP.

Having demonstrated the topological aspect of the DS, we now move to its nonsymmorphic aspect, captured in a topological phase diagram shown in Fig. 3. From the DS phase, we consider symmetry-lowering perturbations 333See the Supplemental Materials for the classification of the perturbations by the point group D4h.. Among diverse possibilities, as a representative example, here we consider a combination of the inversion symmetric ${E_{g}}$- and $B_{2g}$-mode strains and an m_{\scalebox{0.55}{A_{\scalebox{1.0}{$2u$}}}}-mode staggered potential. These perturbations are described by a perturbed Hamiltonian $\mathcal{H}^{1}({\bm{k}})$, where

[TABLE]

For simplicity, we assume the mass parameters are equivalent between the inversion-symmetric perturbations (m_{s}\equiv m_{\scalebox{0.55}{E_{\scalebox{1.0}{$g$}}}}=m_{\scalebox{0.55}{B_{\scalebox{1.0}{$2g$}}}}). Furthermore, we decompose the pristine Hamiltonian $\mathcal{H}^{0}({\bm{k}})$ (Eq. 1) into the inversion-symmetric part $\mathcal{H}^{0}_{+}({\bm{k}})$ and inversion-asymmetric part $\mathcal{H}^{0}_{-}({\bm{k}})$, where

[TABLE]

and

[TABLE]

Here, $v_{-}$ is introduced to parametrize the overall strength of the inversion-asymmetric part.

Figure 3(a) shows a topological phase diagram that is obtained from $\mathcal{H}^{0}+\mathcal{H}^{1}$ in the ($v_{-}$,$m_{s}$,m_{\scalebox{0.55}{A_{\scalebox{1.0}{$2u$}}}}) space. We first note that the DS phase in SG 100 resides along the $|v_{-}|$ (red) axis. From this DS phase, a centrosymmetric strain, described by $m_{s}$, drives a topological phase transition; positive (negative) $m_{s}$ induces a WTI (normal insulator), characterized by $\mathbb{Z}_{2}$ topological indices $(\mu_{0};\mu_{1},\mu_{2},\mu_{3})=(0;001)$[(0;000)]. Therefore, the DS phase defines a phase boundary between the normal and topological insulator phases Murakami (2008), thus exhibiting the nonsymmorphic nature of the DS. In addition to the WTI phase, we find that a Weyl semimetal (WS) can also be induced from the DS phase by applying the staggered potential (|m_{\scalebox{0.55}{A_{\scalebox{1.0}{$2u$}}}}|>0). Interestingly, we find that the three distinctive WS phases are allowed: (1) one having regular (single) WPs with the Chern number $|\mathcal{C}|=1$, (2) another having double WPs with $|\mathcal{C}|=2$, and (3) the other having both single and double WPs. We also note that an archetypal centrosymmetric DS phase is restored from the DS phase by turning off the noncentrosymmetric interactions $v_{-}=0$, from which a WNL semimetal phase is induced by the $m_{s}$ strains, represented by a vertical green line in the figure.

Figure 3(b) illustrates the detailed process of topological phase transition via the creation and annihilation of WPs along the vertical (yellow) path indicated in Fig 3(a). When varying m_{\scalebox{0.55}{A_{\scalebox{1.0}{$2u$}}}} from 0.1 to 0.4 in the unit of $t_{1}$, the in-plane $\mathcal{C}=1$ WP near $X$ ($Y$) and the in-plane $\mathcal{C}=-1$ WP near the DP of $M$ fuse and annihilate eventually, while the WPs residing on the $k_{z}$-axis find their anti-chiral partners by moving along the $k_{z}$ axis. This inter-TRIM WP annihilation results in the trivial insulator phase. On the other hands, Fig. 3(c) illustrates the evolution of the WPs during the topological phase transition from the WS with $\mathcal{C}=2$ to the WTI phase that occurs along the horizontal (yellow) path indicated in Fig. 3(a). Apart from zero, $m_{s}>0$ splits a WP with $\mathcal{C}=2$ into two $\mathcal{C}=1$ WPs off the $k_{z}$-axis. One of the two $\mathcal{C}=1$ WPs encounters with other two $\mathcal{C}=-1$ WPs from the $k_{z}=0$ plane. This event results in a single $\mathcal{C}=-1$ WP, indicated by a solid green circle. The resultant $\mathcal{C}=-1$ WP is eventually annihilated on the $k_{z}=0$ plane by meeting with another WP with $\mathcal{C}$=1, which originates from the double WP on the $k_{z}$-axis. This annihilation results in a WTI 444See Supplemental Material at http:// for the detailed calculations of the associated topological invariants..

Finally, searching for materials that realize the DS in SG 100, we have found an existing material BaLaCuBO5 Norrestam et al. (1994). BaLaCuBO5 is a layered system in SG 100 as shown in Fig. 4(a). It comprises $p4g$ multilayers with an each layer preserving $C_{4z}$ rotation and double glide-mirrors $g_{x}$ and $g_{y}$ symmetries. Our first-principles calculations, performed using Quantum Espresso package Giannozzi et al. (2009), support that BaLaCuBO5 realizes the proposed DS phase 555See Supplemental Material at http://xxxxx for the details of computational methods and the first-principles results for other candidates.. Figs. 4(b) and 4(c) show the first-principles electronic energy bands of BaLaCuBO5. The sticking of four bands is clear from the band structure, featuring filling-enforced gaplessness Watanabe et al. (2016). A fourfold-degenerate DP is present at $M$, and the hourglass-like band connectivity appear on the $\Gamma$-$X$ line. The hourglass-like band connectivity leads to a band crossing on the $\Gamma$-$X$ line, as shown in the magnified view in Fig. 4(d). The presence of band crossing signals the presence of a WNL that encircles the $M$ point lying on the $k_{x}=0$ plane, which we have confirmed throughout the band calculations performed in the entire BZ. Our results is in good agreement with the a time-reversal invariant topological encyclopedia online, which indicates BaLaCuBO5 as a high-symmetry point topological semimetal Zhang et al. (2019).

In conclusion, we have established the existence of a novel type DSs in three dimensions, characterized by hosting topological surface states and mediating topological phase transitions. Hosting topological surface states in the nonsymmorphic DSs, the proposed DS features a unique topological character unlike archetypal 3D DSs. The surface energy spectrum should give rise to drumhead-like topological surface states, which should be feasible to observe in the BaLaCuBO5 compound using a known experimental technique, such as angle-resolved photoemission spectroscopy (ARPES). Moreover, defining a symmetry-tuned topological critical point between a normal insulator and a WTI, the proposed DS can transform to diverse topological phases by symmetry-lowering perturbations.

I Supplementary Material for “Dual topological nodal line and nonsymmorphic Dirac semimetal in three dimensions”

I.1 First-principles calculations

Our first-principles calculations were performed based on density functional theory (DFT) as implemented in the Quantum Espresso package Giannozzi et al. (2009). We used the Perdew–Burke–Ernzerhof exchange-correlation functional Perdew et al. (1996) and norm–conserving, optimized, designed nonlocal pseudopotentials Rappe et al. (1990). The spin-orbit coupling was fully considered for the electronic structure calculations using a noncollinear scheme. The electronic wave functions were expanded in terms of a discrete set of plane-waves basis within the energy cutoff of 680 eV. The 8$\times$8$\times$4 Monkhorst-Pack ${\bm{k}}$-points were sampled from the first Brillouin zone (BZ) Monkhorst and Pack (1976). The atomic structures were fully relaxed within a force tolerance of 0.005 eV/Å. The lattice constants for relaxed unit cells of BaLaXBY5 family are given in Table 1. To highlight the elementary band representation (EBR) of our interest in the DFT calculation of BaLaCuBO5, the Wannier90 package was exploited to construct the tight-binding Hamiltonian by using maximally-localized Wannier function for dxy orbitals of Cu Mostofi et al. (2008).

I.2 Chern number calculations

In this section, we introduce two computational methods to calculate the Chern number. First one is to efficiently find a Weyl point (WP) during the phase transition, which carries a non-zero Chern number. Then, we divide the BZ into the cubic-grids and calculate the Berry phase on each surface to track the path of the WP. The other one is a standard Wilson loop method that we used to determine the Chern number of the time-reversal-invariant plane of the weak-topological insulator (WTI) or the Dirac point (DP).

Figure S1 illustrates the methods that we employed to calculate the Chern number of WPs during the phase transition. We track the position of WP during the phase transition between trivial and topological insulators, by calculating the Berry phase on the surface of the cubic-grid of BZ as in Fig. S1(a). Since WP plays the role of the monopole of the Berry connection, a cubic-grid with non-zero Berry phase manifests the WPs of net charge equals to its non-zero Berry phase. We divide BZs into $400\times 400\times 400$ cubic-grids to track down the path of WPs as the parameters are changed to complete the phase transition via creation and annihilation of WPs.

On the other hands, the non-abelian Wilson loop calculation provides the technical venue to determine the vanishing Chern number of the DP, as well as $\mathcal{Z}_{2}$ topological invariant Yu et al. (2011). The Wilson loop has a mathematical structure given by

[TABLE]

where the overlap matrix $\left[F({\bm{k}}_{i})\right]_{nm}\equiv\langle\psi_{n}(vvk_{i})|\psi_{m}({\bm{k}}_{i+1})\rangle$ is defined by the inner product of the occupied states at two adjacent momenta ${\bm{k}}_{i}$ and ${\bm{k}}_{j}$ on the closed path $k_{l}$. For equal spacing slides of the closed loop $k_{l}$ with infinitesimal spacing $\Delta{\bm{k}}$, the Wilson matrix is mathematically equivalent to the non-abelian Berry phase:

[TABLE]

where P represents the integral on the closed loop $k_{l}$ is path-ordered, ${\bm{k}}_{0}$ is the starting point of the closed path integral, and $\left[{\bm{A}}({\bm{k}})\right]_{nm}\equiv i\langle\psi_{n}({\bm{k}})|\partial_{{\bm{k}}}\psi_{n}({\bm{k}})\rangle$ is the non-abelian Berry connection on the momentum space.

By sweeping the specific momentum plane with the Wilson matrix, one can determine the $Z$ Chern number or $Z_{2}$ invariant of the plane. Particularly, in the case of the TR symmetrical plane, the even-number-crossing and odd-number-crossing of $\phi(k_{y})$, the phases of eigenvalues of Wilson matrix, indicate the trivial and non-trivial $Z_{2}$ invariant of the plane Yu et al. (2011). In Fig. S1(b-1), the Wilson loops on the $k_{z}=0$ plane is illustrated. The Wilson matrix $\mathcal{W}(k_{y})$ is calculated in the closed loops aligned in $k_{x}$ direction;

[TABLE]

where $k_{x,j}=-\pi+2\pi j/N$. Fig. S1(b-2) shows a clear odd-number-crossing of the eigenphase $\phi(k_{y})$, which manifests that the phase is a topological insulator (TI) phase. We confirm that there exists the WTI phase on the phase diagram generated by the $m_{S}$ perturbation by implementing the Wilson band calculation on equally separated Wilson loop of $N=200$.

Additionally, we calculate the Wilson band on the spherical surface enclosing the Weyl nodal line (WNL) to investigate the relation to the DP which is achieved by shrinking the WNL into a point via restoring the inversion symmetry. The Wilson band sphere enclosing the WNL is parameterized by

[TABLE]

as illustrated in Fig. S1(c-1). The Wilson matrix as a function of $\theta$ is calculated via

[TABLE]

where $\left[(\theta,\varphi_{i})\right]_{n,m}=\langle\theta,\varphi_{i};n|\theta,\varphi_{i+1};m\rangle$ is the overlap matrix of occupied states between two neighboring azimuth angles $\varphi_{i}=2\pi\,i/N_{\varphi}$ and $\varphi_{i+1}=2\pi\,(i+1)/N_{\varphi}$. In Fig. S1(c-2), the phases of eigenvalues of the Wilson matrix, $\phi(\theta)$, are illustrated with adjusting small parameter $k_{0}$. The flowing pair of Wilson bands from [math] to [math] ($=2$) via $-\pi$ and $\pi$ shows an identical winding structure with the Wilson bands of DP as illustrated in Fig. S4(c).

I.3 Band BaLaXBY5 families

The DFT band structures of BaLa$X$B$Y_{5}$ family with $X$ = Cu and Au and $Y$ = O, S, and Se, of which the space symmetry belongs to the space group (SG) 100, are shown in Fig. S2. Every member of this family features the nodal structures with four-band sticking forming and elementary band representations Bradlyn et al. (2017). We found that the four bands near the Fermi level mainly comprise the $d$ orbitals of the $X$ atoms. The Fermi level is well separated from the bands other than the four bands near the Fermi level. The four-band sticking enforces the nodal structure, resulting in a filling-enforced semimetal Watanabe et al. (2016) when the filling is $4\mathbb{Z}+2$ as in the cases of these compounds. Near the $X$ point on the $\Gamma$-$X$ line, the hourglass-like crossing appears forming two-fold degenerate node, which constitutes twofold-degenerate WNL forming along the $k_{z}$ direction. At the same time, the fourfold degenerate DP exists at the $M$ point in all the members of the material family.

Figure S3 shows the comparison between the tight-binding (TB) and first-principles band structure for BaLaCuBO5, where readers can find that the tight-binding (TB) model well reproduces the four DFT bands near the Fermi level. We provided the TB model in the main text as

[TABLE]

Here, the simple first and second nearest-hoppings designated by $t_{1}$, $t_{2}$ are considered in the plane, and along the $k_{z}$ direction, the nearest-hopping designated by $t_{3}$ is considered. The parameter set of $\mathcal{H}^{0}({\bm{k}})$ that we find to reproduce the DFT results best is found as $t_{0}=0.124$, $t_{1}=1.626$, $t_{2}=0.119$, $t_{3}=0.0$, $v_{0}=0.0$, $v_{1}=0.0$, $v_{2}=0.04$, $v_{3}=0.01$, and $v_{4}=0.013$ in the order of eV. We first matched the energy eigenvalues at the time-reversal-symmetrical momenta (TRIMs), then adjusted the rest of parameters to maximally reduce the mismatches via introducing the interactions beyond the nearest- and next-nearest-neighbor interactions.

I.4 Double-glide Dirac node decomposition

The $M$ point hosts a fourfold degenerate DP, protected by double-glide-mirrors. Here we develop the low-energy effective theory for the Dirac point that respects the symmetries of SG $100$ and the time-reversal (TR) symmetry, represented by

[TABLE]

Note that the phase factors of each representation are arranged, such that they satisfy the commutation relations of SG 100 and possess the corresponding eigenvalues. The $k\cdot p$ model respecting these symmetries together with the TR symmetry is given by

[TABLE]

We introduce five anti-commuting $\Gamma$-matrices and ten combinations of them to express the $SU(4)$ invariant Hilbert space. Our choice is $\Gamma_{(1,2,3,4,5)}=(\tau_{x}\sigma_{y},\tau_{x}\sigma_{x},\tau_{y},\tau_{x}\sigma_{z},\tau_{z})$, which is followed by their ten combinations of $\Gamma_{ab}\equiv-i\Gamma_{a}\Gamma_{b}~{}(a\neq b)$. Note that $\Gamma_{1\sim 4}$ are $\mathcal{T}$-odd

[TABLE]

while $\Gamma_{5}$ is $\mathcal{T}$-even

[TABLE]

With the $\Gamma$-matrix set, one can rewrite the Hamiltonian of Eq. (S8) as following.

[TABLE]

The eigenenergies for the Hamiltonian in particular momentum space are given as follows.

[TABLE]

We find that the symmetry-preserving $k\cdot p$ model is identical with the Taylor-expanded TB result in the vicinity of $M$ point, which is clearly shown by setting the coefficients to $c_{1}=-v_{1}$, $c_{2}=v_{4}$, $c_{3}=-v_{2}$, $c_{4}=v_{3}$ , and $c_{5}=0$. The term $c_{5}\neq 0$ is responsible for a higher-order hopping beyond the next-nearest-neighbor interaction, which we excluded in our TB model. Figure S4(a) shows the energy-momentum relationship on $k_{x}-k_{y}$ plane, obtained from the $k\cdot p$ Hamiltonian. Unlike a conventional centrosymmetric Dirac cone, in which all the bands doubly degenerate away from the DP, our double-glide Dirac semimetal (DGDS) DP exhibits bifurcation of the conduction and valence bands away from the DP except on the glide-invariant $k_{x}=\pi$ and $k_{y}=\pi$ lines, due to the absence of inversion symmetry.

By calculating the Chern number via the Wilson band calculation, we convince ourselves that the DP protected by double-glide-mirrors without inversion symmetry carries zero Chern number, as it is supposed to be since a Dirac cone is a composite of two WPs with opposite Chern numbers $\pm 1$. The detailed description for the Wilson matrix calculations is provided in Section I.2. To evaluate the Chern number of the DP, which is turned out to be zero, the Wilson band calculation is implemented on the sphere enclosing the double-glide DP parameterized as in Eq. (S4), where the center of the sphere is set to the DP. The Wilson matrix as a function of $\theta$ is calculated via Eq. (S5). In Fig. S4(c), we plot the phases of eigenvalues of the Wilson matrix, $\phi(\theta)$, with adjusting small parameter $k_{0}$. The pair of Wilson bands flows from [math] to [math] ($=2$) via $-\pi$ and $\pi$, respectively, which is a typical Chern number calculation for a DP that comprises two WP with opposite Chern number $\pm 1$. Our calculation proves that the genuine three-dimensional DP occurs from the $p4g$ multiplayer in SG 100, in which inversion is absence.

I.5 Symmetry lowering perturbation

In this Section, we provide detailed information about symmetry lowering perturbations of SG 100. In Table 2, we classify the possible symmetry lowering perturbations in the minimal model of SG 100. Note that we consider only the perturbations that open the gap of at least one TRIM. In the first column, we classify the perturbations by the symmetry representations of the point group D4h. Note that the subscript $u$ ($g$) is introduced to designate the inversion symmetric (asymmetric) perturbation. We provide the atomistic illustrations in the second column for the perturbations that cause the corresponding symmetry-lowering perturbation, which can be considered as a combination of uni-axial strain, buckling, and staggered potential. We also provide the representation of the perturbation in the third column. In the last four columns, we inform whether the perturbation opens an energy band gap at the corresponding TRIM, $\Gamma$, $X$, $Y$, and $M$. In the case where multiple perturbations are allowed in a class, they are separate with horizontal lines.

The perturbation of the first class $A_{2u}$ breaks double-glide-mirrors $g_{x}$ and $g_{y}$ but preserves the $C_{4z}$ symmetry. The different on-site energies of the two sublattices can be achieved by substituting one of the two sublattices. On the other hands, $B_{1u}$ class breaks $C_{4z}$, while preserving double-glide-mirrors $g_{x}$ and $g_{y}$, which can be achieved by applying a uniaxial strain. Both $B_{2g}$ and $B_{2u}$ classes break $C_{4z}$ and double-glide-mirrors $g_{x}$ and $g_{y}$. The only difference between $B_{2g}$ and $B_{2u}$ is inversion symmetry; $B_{2g}$ additionally breaks the inversion symmetry, while $B_{2u}$ preserves it. The combination of uniaxial strain and inversion-preserving buckling or inversion-breaking substitution can generate $B_{2g}$ and $B_{2u}$, respectively. $E_{g}$ and $E_{u}$ classes preserve one of the double-glide-mirrors, $g_{x}$ or $g_{y}$ and break the other glide as well as $C_{4z}$ symmetry. $E_{g}$ can be achieved applying shear stress in the $x$-$z$ or $y$-$z$ plane, while additional glide-preserving buckling generates $E_{u}$ by breaking inversion symmetry.

For the phase diagram, $A_{2u}$, $B_{2g}$ and $E_{g}$ classes of the perturbations are considered. To study the effect of breaking glide-mirrors of the DGDS, we consider the class $A_{2u}$, which breaks double-glide-mirrors selectively, and m_{\scalebox{0.55}{A_{\scalebox{1.0}{$2u$}}}}\tau_{x} is chosen since it affects the entire BZs. Also, $B_{2g}$ and $E_{g}$ classes are considered with the expectation that they drive DGDS to the (weak) TI phase. The perturbed Hamiltonian for these perturbations can be written as

[TABLE]

By adding the perturbation $H^{1}({\bm{k}})$ to the pristine Hamiltonian $H^{0}({\bm{k}})=H^{0}_{+}({\bm{k}})+H^{0}_{-}({\bm{k}})$, we study the phase diagram achievable via applying the corresponding symmetry-lowering perturbations. Note that we divide the $H^{0}_{+}({\bm{k}})$ and $H^{0}_{-}({\bm{k}})$ are the centrosymmetric and non-centrosymmetric parts of the pristine Hamiltonian are given respectively by

[TABLE]

and

[TABLE]

Similar to a regular nonsymmorphic Dirac semimetal (DS), which defines a critical point between normal insulator (NI) and topological insulator (TI) phases Murakami (2008); Young et al. (2012); Young and Kane (2015), the DGDS occurs at a phase boundary between NI and TI phases tuned by symmetry-lowering perturbations. We demonstrated this by explicitly calculating $\mathcal{Z}_{2}$ topological invariant when a band gap is opened by symmetry-lowering perturbations. To efficiently calculate the $\mathcal{Z}_{2}$ topological invariant, we start from the centrosymmetric limit and push it into the noncentrosymmetric limits until a band gap closes and reopens, which signals a topological phase transition. In the centrosymmetric limit, one can easily determine the TI phase by calculating the Fu-Kane $\mathbb{Z}_{2}$ index $(\nu_{1},\nu_{2},\nu_{3};\nu)$, which can be calculated as

[TABLE]

where product $i$ runs over TRIMs at high-symmetrical plane $j$ of BZ, and $\xi_{i}$ corresponds to the inversion-symmetry eigenvalues Fu et al. (2007). The absence of band gap closing guarantees the same topological insulator state that we find from the centrosymmetric limit in the noncentrosymmetric region. The inversion symmetric limit can be obtained by setting $v_{-}=0$ and m_{\scalebox{0.55}{A_{\scalebox{1.0}{$2u$}}}}=0, where every band is doubly degenerate due to the Kramers theorem. One can open the energy gap at TRIMs by setting $t_{1}$, m_{\scalebox{0.55}{E_{\scalebox{1.0}{$g$}}}}^{(1)}, m_{\scalebox{0.55}{E_{\scalebox{1.0}{$g$}}}}^{(2)} and m_{\scalebox{0.55}{B_{\scalebox{1.0}{$2g$}}}} nonzero. Without considering any accidental band crossing between the occupied and conduction bands off high-symmetry momenta, the Kane-Mele $\mathbb{Z}_{2}$ invariant becomes $(0,0,\zeta;0)$, where $\zeta$ is given by

[TABLE]

Here, the representations for inversion symmetry at TRIMs $\Gamma$, $X$, $Y$ and $M$ are given by

[TABLE]

We confirm that Eq. (S16) is consistent with the $\mathbb{Z}_{2}$ invariant calculated by the non-abelian Wilson loop as illustrated in Fig. S1(b). Since only the signs of the parameters m_{\scalebox{0.55}{E_{\scalebox{1.0}{$g$}}}}^{(1)}, m_{\scalebox{0.55}{E_{\scalebox{1.0}{$g$}}}}^{(2)} and m_{\scalebox{0.55}{B_{\scalebox{1.0}{$2g$}}}} matters to the $\mathbb{Z}_{2}$ invariant, we set m_{\scalebox{0.55}{E_{\scalebox{1.0}{$g$}}}}^{(1)}=m_{\scalebox{0.55}{E_{\scalebox{1.0}{$g$}}}}^{(2)}=m_{\scalebox{0.55}{E_{\scalebox{1.0}{$g$}}}} as a representative system for simplicity. Then, the perturbation in the Hamiltonian becomes

[TABLE]

Also, we further simplify the Hamiltonian by setting \mathcal{H}^{1}({\bm{k}})=m_{\scalebox{0.55}{E_{\scalebox{1.0}{$g$}}}}=m_{\scalebox{0.55}{B_{\scalebox{1.0}{$2g$}}}}=m_{S}, of which the results are presented in the main text. In the following Section, the detail description of the phase diagrams in terms of parameters $v_{-}$, $m_{s}$, m_{\scalebox{0.55}{A_{\scalebox{1.0}{$2u$}}}} is provided.

I.6 $v_{-}$, $m_{s}$, m_{\scalebox{0.55}{A_{\scalebox{1.0}{$2u$}}}} phase diagram

Here, we provide detailed results of the perturbed Hamiltonian $\mathcal{H}^{\prime}({\bm{k}})=H^{0}({\bm{k}})+H^{1}({\bm{k}})$ in the DGDS and establish phase diagrams generically accessible from the DGDS, where $H^{0}({\bm{k}})$ and $H^{1}({\bm{k}})$ are defined in Eqs. (S14) and (S18). In Fig. S5, 2D version of the phase diagram for parameters $v_{-}$, $m_{s}$ and m_{\scalebox{0.55}{A_{\scalebox{1.0}{$2u$}}}} are given. As illustrated in Fig. S5(a) and S5(b), the m_{s}-m_{\scalebox{0.55}{A_{\scalebox{1.0}{$2u$}}}} phase diagrams for $v_{-}=0$ and $v_{-}=1$ share a similar sketch, with few differences in detail. In common, the white, gray, and blue colored areas correspond to normal insulator (NI), $\mathcal{C}=1$ WS, and WTI phases, respectively. In contrast, nonsymmorphic DS (orange dot) and WNL (green line) phase are found in $v_{-}=0$ limit, while DGDS phase replaces the nonsymmorphic DS and the WNL phase vanishes when $v_{-}=1$. The remaining uncolored white area represents the NI phase. All the phases, NI, WS, WNL, nonsymmorphic DS, and DGDS phases, are recovered in the v_{-}-m_{\scalebox{0.55}{A_{\scalebox{1.0}{$2u$}}}} phase diagram in Fig. S5(c), only WS phase is distinguished by a different Chern number.

We analyze the symmetry group of which the perturbed system satisfied. We start by the unperturbed system in the inversion symmetric limit, where $\tau_{x}$ is chosen for the representation of (restored) inversion. The corresponding symmetry group belongs to SG 125, where the nonsymmorphic DS phase resides. The perturbation m_{\scalebox{0.55}{A_{\scalebox{1.0}{$2u$}}}} breaks the inversion and glide-mirrors in the way that the two-fold rotations $C_{2x,2y}$ survive, and brings the system to SG 89 where WNL phase can exist. On the other hands, increasing $m_{S}$ from SG 129 brings the system to SG 13, of which the generators are the inversion $I$ and the glide $g_{z}$. In SG 13, the WTI phase can be obtained by band inversion.

In Fig. S6, the sketches of the nodal structures of each phase in Fig. S5 deforming each other are given. Fig. S6(a) shows the entire nodal structures of the whole BZ, while detailed local illustrations in the vicinity of $X$ ($R$) and $M$ ($R$) points are given in Figs. S6(c) and S6(d). The DPs can be achieved by compressing the WNLs, which corresponds to turning off m_{\scalebox{0.55}{A_{\scalebox{1.0}{$2u$}}}} from the WNL phase as illustrated in (b) and (c). The DP at $X$ ($R$) point is the compression of single WNL while DP at $M$ ($A$) is the compression of two WNL crossing at $k_{z}$-axis. Note that WNL at $X$ ($R$) point the DGDS phase deforms to the DP of DS phase in the same way of (b).

On the other hands, increasing $v_{-}$ from the WNL phase shrinks the WNL into WPs of which the positions are located on the trace of the WNL. WNL at $X$ ($R$) is deformed into two $\mathcal{C}=-1$ ($\mathcal{C}=1$) on the crossing points of WNL and $k_{z}$-axis and two $\mathcal{C}=1$ ($\mathcal{C}=-1$) WPs on the $k_{z}=0$ ($k_{z}=\pi$) plane, while WNL at $M$ ($A$) point is deformed into two $\mathcal{C}=2$ ($\mathcal{C}=-2$) WPs on the crossing points of the WNL and $k_{z}$-axis and four $\mathcal{C}=-1$ ($\mathcal{C}=1$) WPs on the crossing points of the WNL and $k_{z}=0$ ($k_{z}=\pi$) plane. Increasing $m_{s}$ form the WS phase separate each $\mathcal{C}=2$ WP into two $\mathcal{C}=1$ WPs, respectively, then yields the complex movements of WPs as illustrated in the main text and ends up with the WTI phase. Indeed, in general, the phase transition occurs along a path other than those of in Fig. S6, but the basic formula of the deformation of the nodal structures follows the sketches in Fig. S6.