Topological hinge modes in Dirac semimetals

Abstract

Dirac semimetals (DSMs) are an important class of topological states of matter. Here, focusing on DSMs of band inversion type, we investigate their boundary modes from the effective model perspective. We show that in order to properly capture the boundary modes, k-cubic terms must be included in the effective model, which would drive an evolution of surface degeneracy manifold from a nodal line to a nodal point. Sizable k-cubic terms are also needed for better exposing the topological hinge modes in the spectrum. Using first-principles calculations, we demonstrate that this feature and the topological hinge modes can be clearly exhibited in β-CuI. We extend the discussion to magnetic DSMs and show that the time-reversal symmetry breaking can gap out the surface bands and hence is beneficial for the experimental detection of hinge modes. Furthermore, we show that magnetic DSMs serve as a parent state for realizing multiple other higher-order topological phases, including higher-order Weyl-point/nodal-line semimetals and higher-order topological insulators.

Graphical abstract

Keywords

topological / hinge / Dirac / semimetals

1 Introduction

The study of topological states and topological materials is an important research topic in the past two decades [1-11]. An important property of topological states is the bulk-boundary correspondence, i.e., the nontrivial topology in the bulk of a system would manifest as protected modes at the boundary. For example, a two-dimensional (2D) quantum anomalous Hall insulator features chiral zero-modes at its 1D edges [12]. As another example, 3D Weyl semimetals have protected surface Fermi arcs connecting the protections of bulk Weyl points with opposite chirality [10, 13]. The existence of surface Fermi arcs can be argued by considering a cylindrical surface in the Brillouin zone (BZ) that encloses one Weyl point [13]. By the Gauss Law, this 2D sub-system is essentially a 2D quantum anomalous Hall insulator, and the corresponding chiral zero-mode traces out a Fermi arc on the surface when we vary the radius of the cylinder.

Dirac semimetals (DSMs) are an important class of topological states that are closely related to Weyl semimetals [14-17]. In a DSM, the bands cross at isolated Dirac points at Fermi level. Each Dirac point is fourfold degenerate and can be regarded as formed by merging together a pair of Weyl points with opposite chirality. Because of this, a Dirac point does not have a net chiral charge (or a nontrivial Chern number). Previous works have shown that there are two types of Dirac points according to their formation mechanism [10, 14]. One type is the essential Dirac points, whose existence is enforced by certain nonsymmorphic space group symmetry [14, 18]. The other type is the accidental Dirac points, which is associated with band inversion in a region of the BZ [15, 16]. On the experimental side, the latter type attracted more interest, because it finds good material realizations, such as Na{}_3Bi and Cd{}_3As{}_2, and also because it hosts interesting boundary modes [19-26]. Initial first-principles calculations showed that Na{}_3Bi and Cd{}_3As{}_2 have surface Fermi arcs connecting the projections of bulk Dirac points [15, 16], similar to those in Weyl semimetals. However, subsequent studies pointed out that such surface arcs are not protected [27]. Recently, with the development of the concept of higher-order topology [28-39], Wieder et al. [40] found that these DSMs actually have a second-order topology with hinge Fermi arcs. A more recent work by Fang and Cano [41] presented a systematical classification of Dirac points with hinge modes.

In this work, we focus on this type of DSMs with band inversions and investigate the evolution of boundary modes from low-energy effective models. We show that in order to correctly capture the topology and boundary modes, the effective model must include terms beyond the second order in the momentum. Particularly, with the inclusion of k-cubic terms, there is an evolution of the surface degeneracy manifold from an open nodal line to a nodal point. This understanding offers guidance to search for materials with hinge modes that can be more readily probed in practice. We show that this is the case for \beta-CuI. Its hinge modes are directly exposed in first-principles calculations. We further extend the discussion to magnetic DSMs and show that the time reversal symmetry breaking can completely gap out the surface bands while maintaining the hinge modes, which could be beneficial for the experimental detection of hinge states. Furthermore, we show that magnetic DSMs are good platforms for realizing a variety of other higher-order topological states, such as magnetic higher-order Weyl semimetals, nodal-line semimetals, and topological insulators. Since effective models are widely used for understanding topological states, our findings have important implications on theoretical studies based on the such models. The results also point to concrete materials for which the topological hinge modes can be verified in experiment.

2 Effective model analysis

DSMs with band inversions such as Na{}_3Bi and Cd{}_3As{}_2 share similar low-energy band structures [15, 16]. They feature band inversion around a high-symmetry point (such as \mathrm{\Gamma }) in the BZ, and a pair of Dirac points are protected on a rotational axis that passes through the high-symmetry point. The commonly used low-energy effective model to study these DSMs is

(1)\begin{array}{r}{H}_0(\mathit{k})=\epsilon (\mathit{k})+M(\mathit{k}){\sigma }_z{s}_0+A{k}_x{\sigma }_x{s}_z-A{k}_y{\sigma }_y{s}_0,\end{array}

where the momentum \mathit{k} is measured from the band inversion high-symmetry point, {\sigma }_i and s_i are two sets of Pauli matrices, the functions \epsilon (\mathit{k})=C_0+C_1{k}_z^2 +C_2( k_x^2+k_y^2), M(\mathit{k})=M_0- M_1{k}_z^2- M_2(k_x^2+k_y^2), and C’s, M’s, and A are real model parameters. This model is expanded to the k-square order, which can describe the band inversion feature at the k=0 point if we require the M's share the same sign. Without loss of generality, we assume M_0,M_1 ,M_2> 0.

As we shall show in a while, the conventional model in Eq. (1) is not sufficient to capture the second-order topology and the correct boundary modes. To remedy this, expansion beyond the k square order is needed. Here, we shall include the k cubic terms, which are sufficient for the task.

Obviously, the form of the k cubic terms depends on the crystal symmetry of the material to be considered. To be specific, let’s consider the constraint of D_{6h} point group symmetry, which applies to the material Na{}_3Bi. In Appendix B, we also present the analysis for the D_{4h} point group (applying to Cd{}_3As{}_2), which leads to slightly different terms, but the qualitative results regarding their influence on the topology are not affected. Considering the constraint from time-reversal symmetry \mathcal{T}=-\mathrm{i}{\sigma }_0{s}_y\mathcal{K} (\mathcal{K} the complex conjugation) and the generators of the D_{6h} group: \mathcal{P}={\sigma }_z{s}_0, {\mathcal{M}}_x =\mathrm{i}{\sigma }_0{s}_x and C_{6z}={\mathrm{e}}^{\mathrm{i}(\pi /3)\widehat{J_z}/\hslash }={\mathrm{e}}^{\mathrm{i}(\pi /3)(2{\sigma }_0-{\sigma }_z)s_z}, the symmetry-allowed k cubic terms include

(2)\begin{array}{r}{H}_1(\mathit{k})=B{k}_z[(k_x^2-k_y^2){\sigma }_x{s}_x+2k_x{k}_y{\sigma }_x{s}_y],\end{array}

with B a real parameter. Note that besides H_1, there are additional k-cubic terms proportional to the last two terms in Eq. (1) timed by k_i^2. However, these terms are not important for our discussion, so they are neglected here.

The spectrum of the effective model H=H_0+H_1 can be readily solved, which is given by

(3)\begin{array}{r}{E}_{\pm }(\mathit{k})=\epsilon (\mathit{k})\pm \sqrt{M(\mathit{k}{)}^2+A^2{k}_{+}{k}_{-}+|B{k}_z{k}_{+}^2{|}^2},\end{array}

where k_{\pm }=k_x\pm \mathrm{i}{k}_y, and each band is doubly degenerate due to the \mathcal{P}\mathcal{T} symmetry. The bands cross at two Dirac points located at (0,0,\pm k_D) on the high-symmetry axis, with k_D=\sqrt{M_0/M_1}. Around each Dirac point, the band dispersion is linear in k at the leading order. For example, expanding the dispersion at (0,0,+ k_D), we have E(\mathit{q})=\pm \sqrt{A(q_x^2+q_y^2)+4M_0{M}_1{q}_z^2}\sim q, where \mathit{q} and the energy are measured from the Dirac point.

Now, we analyze the boundary modes of this effective model. First, it is noted that the k-cubic terms in H_1 do not affect the bulk Dirac point features. For instance, if we put B=0 in Eq. (3), one finds that the location of the Dirac points and the leading order dispersion are not affected at all, which seems to imply that H_1 is inessential. Hence, let us first consider the surface spectrum by neglecting the H_1 term. The calculation results are presented in Fig.1. Here, to study a surface, we discretize the model on a hexagonal lattice as in Fig.1(a). The bulk band structure in Fig.1(c) captures the low-energy features, particularly the Dirac points on the \mathrm{\Gamma }-A path. In Fig.1(d, e), we plot the calculated spectrum for the side surface normal to \hat{y}, where the projections of the two bulk Dirac points can be well distinguished. In Fig.1(f), one clearly observes a pair of surface Fermi arcs connecting the two projected Dirac points, which are similar to the previous first-principles results on Na{}_3Bi and Cd{}_3As{}_2 [15, 16]. These Fermi arcs are formed by the cutting of Fermi energy with the surface bands indicated in Fig.1(e). One can see that the surface bands linearly cross on the \overline{\mathrm{\Gamma }}-\overline{Z} path in the surface BZ between the surface protections of Dirac points at (0,\pm k_D), which form a nodal line connecting the projected Dirac points in the surface band structure. This picture can be better visualized in Fig.1(d), which maps out the surface band dispersion.

The surface spectrum for H_0 can be understood in the following. Consider a slice in the BZ with constant k_z=\lambda for H_0, which constitutes a 2D sub-system {\tilde{H}}_0^{\lambda }(k_x,k_y) labeled by \lambda. We have

(4)\begin{array}{rl}{\tilde{H}}_0^{\lambda }=& \epsilon (k_x,k_y,\lambda )+M(k_x,k_y,\lambda ){\sigma }_z{s}_0+A{k}_x{\sigma }_x{s}_z\\ & -A{k}_y{\sigma }_y{s}_0.\end{array}

One finds that this 3D model has exactly the same form as the famous Bernevig-Hughes-Zhang model [42] for 2D topological insulators. Particularly, the model is topologically nontrivial for , i.e., for a 2D slice in the region between the two Dirac points, which is consistent with the assumed band inversion feature around \mathit{k}=0. Thus, each constant k_z slice between the two Dirac points is effectively a 2D topological insulator, which has a pair of 1D edge bands forming a Dirac type crossing. The crossing traces out the surface nodal line connecting the two Dirac points on a side surface. This clarifies the origin of the surface spectrum of H_0 in Fig.1(d)−(f).

It must be noted that a conventional topological insulator requires the protection of the time reversal symmetry. In the 2D model {\tilde{H}}_0^{\lambda }, we have an anti-unitary symmetry {\mathcal{T}}^{\ast }=- \mathrm{i}{\sigma }_0{s}_y \mathcal{K}, which resembles but is not the true time reversal symmetry for k_z\ne 0, because in the 3D system, time reversal operation should also reverse the sign of k_z. It follows that the surface spectrum in Fig.1(d) is enabled by an emergent symmetry ({\mathcal{T}}^{\ast }) limited to H_0, but not protected by any true symmetry of the original system. As a result, the surface bands with a nodal line represents a critical state susceptible to perturbations from higher-order terms.

Next, we show that restoring the k-cubic terms in H_1 helps to capture the correct topology. Note that by putting k_z=\lambda in H_1, we obtain its contribution to the 2D sub-system of a constant k_z slice:

(5)\begin{array}{r}{\tilde{H}}_1^{\lambda }=B\lambda [(k_x^2-k_y^2){\sigma }_x{s}_x+2k_x{k}_y{\sigma }_x{s}_y].\end{array}

Clearly, {\tilde{H}}_1^{\lambda } breaks the emergent symmetry {\mathcal{T}}^{\ast } of {\tilde{H}}_0^{\lambda }. In other words, if we treat {\tilde{H}}_0^{\lambda } as describing a 2D {\mathcal{T}}^{\ast }-invariant topological insulator, {\tilde{H}}_1^{\lambda } can be regarded as perturbations that break the effective time reversal symmetry. Consequently, the 1D Dirac type crossing in the edge bands for {\tilde{H}}_0^{\lambda } would open a gap. This is confirmed by the calculated surface spectrum in Fig.2(b) by including the H_1 term, which destroys the surface nodal line. It should be noted that the k_z=0 slice is special as it preserves the true time reversal symmetry, so it remains a 2D topological insulator with gapless edge bands. For the 3D system, this means that although the surface nodal line is destroyed, there is still a robust nodal point of the surface bands at \overline{\mathrm{\Gamma }}.

This feature can also be understood from another perspective. Note that the bulk Dirac points are protected by the rotational symmetry on the k_z axis. They can be gapped out by breaking the rotational symmetry while preserving \mathcal{T}. Then the system would transform to a 3D strong topological insulator because of the assumed band inversion at \mathrm{\Gamma }. It is well known that a 3D topological insulator features Dirac-cone type surface bands. This explains the Dirac type surface dispersion in Fig.1(b−d), and the nodal point is just the neck point of the surface Dirac cone. This discussion clarifies the important role played by H_1, under which the surface bands evolve from Fig.1(d) with a nodal line to Fig.2(b) with a Dirac cone. Inspecting the Fermi contour at the surface, the Fermi arcs in Fig.1(f) would generally transform into a closed loop in Fig.2(d), similar to that in a 3D strong topological insulator.

We have shown that by including the k-cubic terms, the 2D sub-system described by {\tilde{H}}^{\lambda }(k_x,k_y)= {\tilde{H}}_0^{\lambda }+{\tilde{H}}_1^{\lambda } with and \lambda \ne 0 is no longer a 2D conventional topological insulator. Both its bulk and its edge are gapped. Nevertheless, the band inversion feature is still maintained in the model, and we will show that {\tilde{H}}^{\lambda } corresponds to a 2D second-order topological insulator. The second-order topology can be inferred from the nested Wilson loop calculation [29]. In Fig.3(b), we plot the obtained nested Berry phase as a function of \lambda. One observes that the result is nontrivial (trivial) for (>k_D). Thus, each constant k_z( \ne 0) slice of the BZ between the two Dirac points is effectively a 2D second-order topological insulator.

A 2D second-order topological insulator should have protected corner modes. We implement {\tilde{H}}^{\lambda } on a hexagonal lattice and plot the calculated spectrum for a nanodisk geometry in Fig.3(c) and (d). Here, we take \lambda =\pi /3. When we put B=0, i.e., drop the k-cubic terms, the zero-modes are distributed throughout the edge of the disk [Fig.3(c)]. This is the critical state, for which the system resembles the conventional topological insulator with gapless edge modes. As soon as we turn on the k-cubic terms, the edge becomes gapped and the zero-modes are localized at the corners of the disk [Fig.3(d)], confirming the second-order topology.

Since {\tilde{H}}^{\lambda } is a constant k_z slice of the DSM, its corner modes would constitute the hinge modes at hinges between the side surfaces of a 3D DSM. To explicitly demonstrate this, we consider a tube geometry as shown in Fig.2(f). The obtained spectrum in plotted in Fig.2(e), in which the hinge modes are marked with red color. In Fig.2(f), we verify that these modes are indeed distributed at the hinges between the side surfaces of the system.

From the model study, we have seen that the k-cubic terms are indispensable for describing the correct boundary modes of the DSM. On the 2D surface, the generic Fermi contour is a Fermi loop from the Dirac-cone surface bands. The bulk band inversion leads to second-order topology with hinge modes bounded by the projected Dirac points on the 1D hinges between side surfaces.

3 Material example

The analysis in the last section shows that to better visualize the hinge modes, the system needs to have sizable k-cubic terms. In the effective model we discussed, the k-cubic terms correspond to a kind of spin-orbit coupling in the next neighbor hopping process. This suggests that we should look for materials with strong spin-orbit coupling in the low-energy bands. In materials Na{}_3Bi and Cd{}_3As{}_2, the cubic terms are relatively small, which makes the surface Fermi contour still close to Fermi arcs. And the hinge modes there coexist in energy with the surface modes for a fixed k_z, making it difficult to resolve the hinge modes in the spectrum.

Here, we show that \beta-CuI is a good candidate to probe the hinge modes. The previous work by Le et al. [43] has revealed \beta-CuI as a DSM formed by band inversion. Here, we find that this material has sizable k-cubic terms, and we shall directly investigate its hinge modes.

As illustrated in Fig.4(a), the structure of \beta-CuI belongs to the space group R\overline{3}m (No. 166), same as the famous topological insulator Bi{}_2Se{}_3 family. From the crystal field environment, one observes that the iodine atoms can be classified as two types denoted as I{}_1 and I{}_2, where I{}_1 is octahedrally coordinated by six Cu atoms forming a sandwich ABC tri-layer stacking, while I{}_2 connects two Cu atoms parallel to the c axis separating the Cu−I{}_1−Cu sandwich layer. The relaxed lattice constants are a=4.3710 Å and c=20.8611 Å in the hexagonal lattice description (see Appendix A for the computation approach), which are in good agreement with the experimental results (a=4.2986 Å, c=21.4712 Å) [44]. The Wyckoff positions of Cu, I{}_1 and I{}_2 are 6c (0, 0, 0.1246), 3a (0, 0, 0) and 3b (0, 0, 0.5), respectively.

In \beta-CuI, the p orbitals of I{}_1 atoms and the p_z orbitals of I{}_2 atoms are strongly affected by the crystal fields from the surrounding Cu atoms and are repelled away from the Fermi level. Meanwhile, due to the positive valence of Cu, the d orbitals of Cu are completely filled and are located at around -2.5 eV. Therefore, near the Fermi level, the valence and conduction bands are mainly contributed by the I{}_2-5p_{x,y} and Cu-4s orbitals. Our first-principles result confirms this analysis. Fig.4(c) shows the band structure and projected density of states (PDOS) of \beta-CuI without spin-orbit coupling (SOC). Around the Fermi energy, there is an energy band inversion at the \mathrm{\Gamma } point, caused by the Cu-4s and the I{}_2-5p_{x,y} orbitals. The Cu-4s bands are about 0.47 eV lower than the I{}_2-5p_{x,y} bands, and there is a band crossing point along the \mathrm{\Gamma }-\mathrm{Z} line. After turning on SOC, the band inversion at \mathrm{\Gamma } is enhanced to 0.77 eV, and the band crossing along \mathrm{\Gamma }-\mathrm{Z} still exists [Fig.4(d)]. Each band here is doubly degenerate due to \mathcal{P}\mathcal{T}. The irreducible representations of the two crossing bands belong to {\mathrm{\Gamma }}_4 and {\mathrm{\Gamma }}_5\oplus {\mathrm{\Gamma }}_6 of C_{3v} group along \mathrm{\Gamma }-\mathrm{Z}, respectively. Therefore, the crossing point is a fourfold Dirac point, consistent with the previous result [43].

Now, we turn to the surface spectrum of \beta-CuI. Fig.5(a) shows the calculated surface spectrum for the (100) surface. One observes features similar to those in Fig.2(c). Particularly, one can see the large splitting of the nodal line on the \overline{\mathrm{\Gamma }}-\overline{\mathrm{Z}} path between the projected Dirac points, and the Fermi contour takes the form of a loop rather than arcs [Fig.5(b)]. These evidences indicate sizable k-cubic terms which break the effective {\mathcal{T}}^{\ast } symmetry.

The surface spectrum in Fig.5(a) suggests that there is a good chance to resolve the hinge modes in the system. To calculate the hinge spectrum, we consider the tube geometry shown in Fig.5(d). The result is plotted in Fig.5(c). Indeed, we find two hinge bands within the surface band gap bounded by the projected Dirac points. (The appearance of a small gap of the bulk spectrum in Fig.5(c) is due to the finite size effect.) By checking the wave function distribution, we verify that these modes are located at the hinges of the sample, as shown in Fig.5(d). The two bands split in energy because they are not connected by symmetry. These hinge modes manifest the second-order topological character of \beta-CuI.

Finally, let’s construct the \mathit{k}\cdot \mathit{p} effective model for \beta-CuI. \beta-CuI has the D_{3d} point group symmetry. The symmetry-constrained model is slightly more complicated than that discussed in the last section, but the qualitative features are the same. Choosing the basis at \mathrm{\Gamma } as | S_{1/2}^{+},\pm 1/ 2⟩ and | P_{3/2}^{-},\pm 3/ 2⟩, the symmetry generators are represented as \mathcal{P}={\sigma }_z{s}_0, {\mathcal{M}}_x =\mathrm{i}{\sigma }_0{s}_x, C_{3z}={\mathrm{e}}^{\mathrm{i}(2\pi /3)\widehat{J_z}/\hslash }= {\mathrm{e}}^{\mathrm{i}(2\pi /3)(2{\sigma }_0-{\sigma }_z)s_z} and \mathcal{T}=-\mathrm{i}{\sigma }_0{s}_y\mathcal{K}.

Then, the symmetry allowed effective model can be obtained as

(6)\begin{array}{r}H(\mathit{k})=H_0+H_1,\end{array}

where H_0 contains terms up to k-square order

(7)\begin{array}{rl}{H}_0& =\epsilon (\mathit{k})+M(\mathit{k}){\sigma }_z{s}_0+A_0\left(k_x{\sigma }_x{s}_z-k_y{\sigma }_y{s}_0\right)\\ & +D_0\left(k_x{\sigma }_x{s}_x-k_y{\sigma }_x{s}_y\right),\end{array}

and H_1 contains k-cubic terms

(8)\begin{array}{rl}{H}_1& =B_1{k}_z[(k_x^2-k_y^2){\sigma }_x{s}_x+2k_x{k}_y{\sigma }_x{s}_y]\\ & +B_2{k}_z[(k_x^2-k_y^2){\sigma }_x{s}_z+2k_x{k}_y{\sigma }_y{s}_0]\\ & +[A_1{k}_z^2+A_2(k_x^2+k_y^2)](k_x{\sigma }_x{s}_z-k_y{\sigma }_y{s}_0)\\ & +[D_1{k}_z^2+D_2(k_x^2+k_y^2)](k_x{\sigma }_x{s}_x-k_y{\sigma }_x{s}_y).\end{array}

Here, \epsilon (\mathit{k}) and M(\mathit{k}) have the same expression as in Eq. (1). The model parameters can be obtained from fitting the first-principles band structure in Fig.4(d). We obtain that C_0=-0.0518\text{eV}, C_1=0.6661\text{eV}\cdot {\text{Å}}^2, C_2=3.1243\text{eV}\cdot {\text{Å}}^2, M_0=0.1930\text{eV}, M_1=4.9640 \text{eV}\cdot {\text{Å}}^2, M_2=0.8866 \text{eV}\cdot {\text{Å}}^2, A_0=-1.5556\text{eV}\cdot \text{Å}, A_1= -0.0937\text{eV}\cdot {\text{Å}}^3, A_2=-0.8030\text{eV}\cdot {\text{Å}}^3, D_0=0.2264 \text{eV}\cdot \text{Å}, D_1=-0.0570\text{eV}\cdot {\text{Å}}^3, D_2= 8.9368\text{eV}\cdot {\text{Å}}^3, B_1=6.7806\text{eV}\cdot {\text{Å}}^3, and B_2= 1.4844 \text{eV}\cdot {\text{Å}}^3. The result shows that the k-cubic terms are sizable for \beta-CuI.

4 Magnetic Dirac semimetal

Since time-reversal symmetry is not a necessary condition for the existence of Dirac points, in this section we discuss hinge modes in DSMs with broken \mathcal{T}, i.e., in magnetic DSMs [45, 46]. Compared to the nonmagnetic DSMs discussed so far, magnetic DSMs exhibit an important difference in the surface spectra. As discussed in Section 2, a nonmagnetic DSM has Dirac-cone type surface bands protected by the \mathcal{T} symmetry. In a magnetic DSM, the \mathcal{T} symmetry is broken, so the surface Dirac cone is generally gapped.

To explicitly demonstrate this point, we construct a four-band lattice model which follows the P4/m^{\prime }mm magnetic space group symmetry (No. 123.341 in Belov-Neronova−Smirnova notation). As shown in Fig.6(a), we take a simple tetragonal lattice with two sites in a unit cell, labeled as A and B sites. At A site, we put two basis orbitals |s\uparrow ⟩ and |s\downarrow ⟩; and at B site, we put |p_{-}\uparrow ⟩ and |p_{+}\downarrow ⟩ as basis (p_{\pm }=p_x\pm \mathrm{i}{p}_y). In these four bases, the generators of the space group are represented by \mathcal{P}\mathcal{T}=-\mathrm{i}{\sigma }_z{s}_y\mathcal{K}, {\mathcal{M}}_x=\mathrm{i}{\sigma }_0{s}_x, C_{4z}={\mathrm{e}}^{\mathrm{i}(\pi /4)(2{\sigma }_0-{\sigma }_z)s_z}. Then, we construct the following minimal model that respects these symmetries:

(9)\begin{array}{rl}\mathcal{H}=& \epsilon (\mathit{k}){\sigma }_0{s}_0+m(\mathit{k}){\sigma }_z{s}_0+v(\sin k_x{\sigma }_x{s}_z+\sin k_y{\sigma }_y{s}_0)\\ & +w(\cos k_x-\cos k_y){\sigma }_x{s}_x,\end{array}

where \epsilon (\mathit{k})=2 {\epsilon }_1\cos k_z+4{\epsilon }_2\cos k_x\cos k_y+{\epsilon }_3\sin k_z, and m(\mathit{k})=m_0 +2m_1\cos k_z+ 4m_2\cos k_x\cos k_y+m_3 \sin k_z. Here, {\epsilon }_0 and m_0 represent on-site energy, with v, w, \epsilon’s and m’s being real parameters. With properly chosen parameters, this model describes a DSM state as shown in Fig.6(b), which has a pair of Dirac points along the k_z axis. Here, the \mathcal{T} symmetry is broken by the w term. If we drop the w term, the model would reduce to a nonmagnetic DSM similar to the ones discussed in Section 2. To see this, we expand the lattice model (9) at the \mathrm{\Gamma } point for small k (the diagonal term \sim {\sigma }_0{s}_0 is dropped since it does not affect the topology). Then, we obtain the following k\cdot p model up to k-quadratic terms:

(10)\begin{array}{rl}{H}_{\text{eff}}=& M(\mathit{k}){\sigma }_z{s}_0+A_z{k}_z{\sigma }_z{s}_0+A{k}_x{\sigma }_x{s}_z+A{k}_y{\sigma }_y{s}_0\\ & +B(k_x^2-k_y^2){\sigma }_x{s}_x,\end{array}

where M(\mathit{k})=M_0-M_1{k}_z^2 -M_2( k_x^2+k_y^2), M_0=-m_0+2m_1+4m_2, M_1=m_1, M_2=2m_2, A_z=m_3, A=v, and B=w/2. This model is very similar to model (1) except for the last term. Importantly, unlike the k-cubic term in Eq. (2), the B( k_x^2-k_y^2){\sigma }_x{s}_x term opens a gap in the 2D subsystem H(k_x, k_y) for any fixed k_z, including the k_z=0 slice, because this term derives from the \mathcal{T}-symmetry breaking w term. It follows that the surface Dirac cone [as in Fig.2(b)−(d)] will be gapped out.

This feature is confirmed by our numerical results shown in Fig.6(c) and (d). One observes that as expected, the Dirac cone at \overline{\mathrm{\Gamma }} is removed and the surface bands are gapped. Meanwhile, the existence of the hinge modes, as corresponding to the second-order topology, is not affected. As shown in [Fig.6(e, f)], due to the absence of the surface Dirac cone, the hinge modes can be more clearly observed in the spectrum. This could be an advantage for the detection of hinge modes. To look for suitable material candidates of magnetic HODSMs, the system should possess proper symmetries (such as \mathcal{P}\mathcal{T} and rotations) to protect the Dirac points. Meanwhile, symmetries like the effective time-reversal symmetry that protect gapless surface states should be broken.

DSMs represents a kind of symmetry-protected topological state. By suitably breaking the protecting symmetry for the Dirac points, we find that DSMs can be transformed into other types of higher-order topological states. We will discuss three examples below.

Starting from the effective model in Eq. (10), let’s first consider adding a perturbation term \mathrm{\Delta }H =D{\sigma }_z{s}_z, which breaks the \mathcal{P}\mathcal{T} symmetry. One finds that each Dirac point would split into two Weyl points located on the k_z axis, as illustrated in Fig.7(a) and (b). Correspondingly, there appear Fermi arcs on the side surfaces, connecting the projections of these Weyl points. Since the original band inversion at \mathrm{\Gamma } is preserved, the hinge modes are maintained in this Weyl semimetal state, as shown in Fig.7(e) and (f). Therefore, the resulting state represents a magnetic higher-order Weyl semimetal. This is similar to the recently reported higher-order Weyl semimetals but is realized with broken time reversal symmetry.

Second, we consider the perturbation \mathrm{\Delta }H =D{k}_z{\sigma }_y{s}_y, which breaks the C_{4z} symmetry. This again destroys the Dirac points. In this case, one finds that each Dirac point evolves into a nodal ring, as illustrated in Fig.8(b). The hinge modes manifesting the higher-order topology are also explicitly shown in Fig.8(c) and (d).

In recent works, the higher-order nodal-line semimetal (HONLSM) state was predicted for \mathcal{P}\mathcal{T} symmetric systems and found its realization in 3D graphdiyne [38, 47]. Note that different from those cases which are spinless, the example here is in a spinful system. It follows that the two nodal rings here are actually not protected by any exact symmetry of the system (so that they would generally be weakly gapped). Nevertheless, the protection may come from an approximate chiral symmetry that could emerge at the low energy.

Finally, the DSM could be fully gapped in the bulk and turned into a magnetic higher-order topological insulator. This could be achieved, e.g., by the perturbation term \mathrm{\Delta }H =D{k}_z{\sigma }_x{s}_y. The result is shown in Fig.9. One observes that the system becomes an insulator, with chiral propagating hinge modes in the bulk gap.

From the above discussion, we can see that DSMs, especially the magnetic DSMs, can serve as a good parent platform to realize various other higher-order topological states and to study the topological phase transitions among them.

5 Conclusion

In this work, we have discussed how to capture the topological boundary modes in the effective model approach to DSMs and how to better expose the hinge modes in the spectrum. We show that the k-cubic terms, which are often neglected in such models, are essential for capturing the correct boundary-mode topology. Using the effective model, we can understand the evolution of surface spectrum driven by the k-cubic terms. Based on such understanding, we show that the surface Dirac cone and the topological hinge modes can be clearly exhibited in \beta-CuI. We show that in magnetic DSMs, the breaking of \mathcal{T} symmetry can gap out the surface Dirac cone while preserving the hinge modes. This could be an advantage for the detection of hinge modes. Furthermore, we show the transitions from magnetic DSMs to other higher-order topological phases, which makes magnetic DSMs a nice platform for investigating a range of states with higher-order topology. Our finding clarifies the key features of the topological boundary modes of DSMs. It has important implications on theoretical studies on DSMs using the effective model approach and on the search of suitable concrete materials with higher-order topology.

6 Appendix A Computation method

The first-principles calculations have been carried out based on the density-functional theory (DFT) as implemented in the Vienna ab initio simulation package (VASP) [49, 50], using the projector augmented wave method [51] and Perdew−Burke−Ernzerhof (PBE) [52] exchange-correlation functional approach. The plane-wave cutoff energy was set to 500 eV. The Monkhorst-Pack k-point mesh [53] of size 8\times 8\times 8 was used for the BZ sampling in bulk calculations. The surface spectrum of \beta-CuI was calculated by constructing the maximally localized Wannier functions (MLWF) [54, 55] and surface Green’s function methods [56, 57] implemented in WannierTools [58].

7 Effective model with D_4h symmetry

Here, we consider the effective model constrained by the D_{4h} symmetry: \mathcal{P}={\sigma }_z{s}_0, {\mathcal{M}}_x =\mathrm{i}{\sigma }_0{s}_x and C_{4z}= {\mathrm{e}}^{\mathrm{i}(\pi /2)\widehat{J_z}/\hslash }={\mathrm{e}}^{\mathrm{i}(\pi /4)(2{\sigma }_0-{\sigma }_z)s_z}. Using the approach discussed in the main text, we find that the model expanded up to k-cubic order reads

\begin{array}{}\text{({rm{B}}1)}& \begin{array}{rl}H& =H_0+H_1,\\ H_0& =\epsilon (\mathit{k})+M(\mathit{k}){\sigma }_z{s}_0+A(\mathit{k})k_x{\sigma }_x{s}_z-A(\mathit{k})k_y{\sigma }_y{s}_0,\\ H_1& =B_1{k}_z(k_x^2-k_y^2){\sigma }_x{s}_x+2B_2{k}_x{k}_y{\sigma }_y{s}_x.\end{array}\end{array}

The functions \epsilon, M, and A have the same form as in model [Eq. (1)]. One can see that the main difference from model [Eq. (1)] is that there is one more independent parameter in H_1. The qualitative features of the surface and hinge spectra are the same as discussed in the main text (see ).

References

[1] M. Z. Hasan, C. L. Kane. Topological insulators. Rev. Mod. Phys., 2010, 82(4): 3045

[2] X. L. Qi, S. C. Zhang. Topological insulators and superconductors. Rev. Mod. Phys., 2011, 83(4): 1057

[3] S.Q. Shen, Topological Insulators, Vol. 174, Springer Berlin Heidelberg, Berlin, Heidelberg, 2012

[4] B.A. BernevigT.L. Hughes, Topological Insulators and Topological Superconductors, Princeton University Press, 2013

[5] A. Bansil, H. Lin, T. Das. Topological band theory. Rev. Mod. Phys., 2016, 88(2): 021004

[6] C. K. Chiu, J. C. Y. Teo, A. P. Schnyder, S. Ryu. Classification of topological quantum matter with symmetries. Rev. Mod. Phys., 2016, 88(3): 035005

[7] S. A. Yang. Dirac and Weyl materials: Fundamental aspects and some spintronics applications. Spin, 2016, 6(2): 1640003

[8] X. Dai. Weyl fermions go into orbit. Nat. Phys., 2016, 12(8): 727

[9] A. A. Burkov. Topological semimetals. Nat. Mater., 2016, 15(11): 1145

[10] N. P. Armitage, E. J. Mele, A. Vishwanath. Weyl and Dirac semimetals in three-dimensional solids. Rev. Mod. Phys., 2018, 90(1): 015001

[11] J. Qi, H. Liu, H. Jiang, X. C. Xie. Dephasing effects in topological insulators. Front. Phys., 2019, 14(4): 43403

[12] F. D. M. Haldane. Model for a quantum Hall effect without Landau levels: Condensed-matter realization of the “parity anomaly”. Phys. Rev. Lett., 1988, 61(18): 2015

[13] X. Wan, A. M. Turner, A. Vishwanath, S. Y. Savrasov. Topological semimetal and Fermi-arc surface states in the electronic structure of pyrochlore iridates. Phys. Rev. B, 2011, 83(20): 205101

[14] S. M. Young, S. Zaheer, J. C. Y. Teo, C. L. Kane, E. J. Mele, A. M. Rappe. Dirac semimetal in three dimensions. Phys. Rev. Lett., 2012, 108(14): 140405

[15] Z. Wang, Y. Sun, X. Q. Chen, C. Franchini, G. Xu, H. Weng, X. Dai, Z. Fang. Dirac semimetal and topological phase transitions in A₃Bi (A = Na, K, Rb). Phys. Rev. B, 2012, 85(19): 195320

[16] Z. Wang, H. Weng, Q. Wu, X. Dai, Z. Fang. Three-dimensional Dirac semimetal and quantum transport in Cd₃As₂. Phys. Rev. B, 2013, 88(12): 125427

[17] S. Li, Z. M. Yu, Y. Yao, S. A. Yang. Type-II topological metals. Front. Phys., 2020, 15(4): 43201

[18] J. A. Steinberg, S. M. Young, S. Zaheer, C. L. Kane, E. J. Mele, A. M. Rappe. Bulk Dirac points in distorted spinels. Phys. Rev. Lett., 2014, 112(3): 036403

[19] Z. K. Liu, B. Zhou, Y. Zhang, Z. J. Wang, H. M. Weng, D. Prabhakaran, S. K. Mo, Z. X. Shen, Z. Fang, X. Dai, Z. Hussain, Y. L. Chen. Discovery of a three-dimensional topological Dirac semimetal, Na₃Bi. Science, 2014, 343(6173): 864

[20] Z. K. Liu, J. Jiang, B. Zhou, Z. J. Wang, Y. Zhang, H. M. Weng, D. Prabhakaran, S. K. Mo, H. Peng, P. Dudin, T. Kim, M. Hoesch, Z. Fang, X. Dai, Z. X. Shen, D. L. Feng, Z. Hussain, Y. L. Chen. A stable three-dimensional topological Dirac semimetal Cd₃As₂. Nat. Mater., 2014, 13(7): 677

[21] M. Neupane, S. Y. Xu, R. Sankar, N. Alidoust, G. Bian, C. Liu, I. Belopolski, T. R. Chang, H. T. Jeng, H. Lin, A. Bansil, F. Chou, M. Z. Hasan. Observation of a three-dimensional topological Dirac semimetal phase in high-mobility Cd₃As₂. Nat. Commun., 2014, 5(1): 3786

[22] S. Jeon, B. B. Zhou, A. Gyenis, B. E. Feldman, I. Kimchi, A. C. Potter, Q. D. Gibson, R. J. Cava, A. Vishwanath, A. Yazdani. Landau quantization and quasiparticle interference in the three-dimensional Dirac semimetal Cd₃As₂. Nat. Mater., 2014, 13(9): 851

[23] S. Borisenko, Q. Gibson, D. Evtushinsky, V. Zabolotnyy, B. Büchner, R. J. Cava. Experimental realization of a three-dimensional Dirac semimetal. Phys. Rev. Lett., 2014, 113(2): 027603

[24] T. Liang, Q. Gibson, M. N. Ali, M. Liu, R. J. Cava, N. P. Ong. Ultrahigh mobility and giant magnetoresistance in the Dirac semimetal Cd₃As₂. Nat. Mater., 2015, 14(3): 280

[25] S. Y. Xu, C. Liu, S. K. Kushwaha, R. Sankar, J. W. Krizan, I. Belopolski, M. Neupane, G. Bian, N. Alidoust, T. R. Chang, H. T. Jeng, C. Y. Huang, W. F. Tsai, H. Lin, P. P. Shibayev, F. C. Chou, R. J. Cava, M. Z. Hasan. Observation of Fermi arc surface states in a topological metal. Science, 2015, 347(6219): 294

[26] J. Xiong, S. K. Kushwaha, T. Liang, J. W. Krizan, M. Hirschberger, W. Wang, R. J. Cava, N. P. Ong. Evidence for the chiral anomaly in the Dirac semimetal Na₃Bi. Science, 2015, 350(6259): 413

[27] M. Kargarian, M. Randeria, Y. M. Lu. Are the surface Fermi arcs in Dirac semimetals topologically protected. Proc. Natl. Acad. Sci. USA, 2016, 113(31): 8648

[28] F. Zhang, C. L. Kane, E. J. Mele. Surface state magnetization and chiral edge states on topological insulators. Phys. Rev. Lett., 2013, 110(4): 046404

[29] W. A. Benalcazar, B. A. Bernevig, T. L. Hughes. Quantized electric multipole insulators. Science, 2017, 357(6346): 61

[30] J. Langbehn, Y. Peng, L. Trifunovic, F. von Oppen, P. W. Brouwer. Reflection-symmetric second-order topological insulators and superconductors. Phys. Rev. Lett., 2017, 119(24): 246401

[31] Z. Song, Z. Fang, C. Fang. (d−2)-dimensional edge states of rotation symmetry protected topological states. Phys. Rev. Lett., 2017, 119(24): 246402

[32] F. Schindler, A. M. Cook, M. G. Vergniory, Z. Wang, S. S. P. Parkin, B. A. Bernevig, T. Neupert, B. Andrei Bernevig, T. Neupert. Higher-order topological insulators. Sci. Adv., 2018, 4(6): eaat0346

[33] F. Schindler, Z. Wang, M. G. Vergniory, A. M. Cook, A. Murani, S. Sengupta, A. Y. Kasumov, R. Deblock, S. Jeon, I. Drozdov, H. Bouchiat, S. Guéron, A. Yazdani, B. A. Bernevig, T. Neupert. Higher-order topology in bismuth. Nat. Phys., 2018, 14(9): 918

[34] X. L. Sheng, C. Chen, H. Liu, Z. Chen, Z. M. Yu, Y. X. Zhao, S. A. Yang. Two-dimensional second-order topological insulator in graphdiyne. Phys. Rev. Lett., 2019, 123(25): 256402

[35] H. X. Wang, Z. K. Lin, B. Jiang, G. Y. Guo, J. H. Jiang. Higher-order Weyl semimetals. Phys. Rev. Lett., 2020, 125(14): 146401

[36] S. A. A. Ghorashi, T. Li, T. L. Hughes. Higher-order Weyl semimetals. Phys. Rev. Lett., 2020, 125(26): 266804

[37] H. Qiu, M. Xiao, F. Zhang, C. Qiu. Higher-order Dirac sonic crystals. Phys. Rev. Lett., 2021, 127(14): 146601

[38] C. Chen, X. T. Zeng, Z. Chen, Y. X. Zhao, X. L. Sheng, S. A. Yang. Second-order real nodal-line semimetal in three-dimensional graphdiyne. Phys. Rev. Lett., 2022, 128(2): 026405

[39] H. D. Scammell, J. Ingham, M. Geier, T. Li. Intrinsic first- and higher-order topological superconductivity in a doped topological insulator. Phys. Rev. B, 2022, 105(19): 195149

[40] B. J. Wieder, Z. Wang, J. Cano, X. Dai, L. M. Schoop, B. Bradlyn, B. A. Bernevig. Strong and fragile topological Dirac semimetals with higher-order Fermi arcs. Nat. Commun., 2020, 11(1): 627

[41] Y. Fang, J. Cano. Classification of Dirac points with higher-order Fermi arcs. Phys. Rev. B, 2021, 104(24): 245101

[42] B. A. Bernevig, T. L. Hughes, S. C. Zhang. Quantum spin Hall effect and topological phase transition in HgTe quantum wells. Science, 2006, 314(5806): 1757

[43] C. Le, X. Wu, S. Qin, Y. Li, R. Thomale, F. C. Zhang, J. Hu. Dirac semimetal in β-CuI without surface Fermi arcs. Proc. Natl. Acad. Sci. USA, 2018, 115(33): 8311

[44] Y. Shan, G. Li, G. Tian, J. Han, C. Wang, S. Liu, H. Du, Y. Yang. Description of the phase transitions of cuprous iodide. J. Alloys Compd., 2009, 477(1−2): 403

[45] P. Tang, Q. Zhou, G. Xu, S. C. Zhang. Dirac fermions in an antiferromagnetic semimetal. Nat. Phys., 2016, 12(12): 1100

[46] G. Hua, S. Nie, Z. Song, R. Yu, G. Xu, K. Yao. Dirac semimetal in type-IV magnetic space groups. Phys. Rev. B, 2018, 98: 201116(R)

[47] K. Wang, J. X. Dai, L. B. Shao, S. A. Yang, Y. X. Zhao. Boundary Criticality of PT-invariant topology and second-order nodal-line semimetals. Phys. Rev. Lett., 2020, 125(12): 126403

[48] S.NieJ. ChenC.YueC.LeD.Yuan W.ZhangH. Weng, Tunable Dirac semimetals with higher-order Fermi arcs in Kagome lattices Pd₃Pb₂X₂ (X = S, Se), arXiv: 2203.03162 (2022)

[49] G. Kresse, J. Hafner. Ab initio molecular-dynamics simulation of the liquid-metal–amorphous-semiconductor transition in germanium. Phys. Rev. B, 1994, 49(20): 14251

[50] G. Kresse, J. Furthmüller. Efficient iterative schemes for ab initio total-energy calculations using a plane-wave basis set. Phys. Rev. B, 1996, 54(16): 11169

[51] P. E. Blöchl. Projector augmented-wave method. Phys. Rev. B, 1994, 50(24): 17953

[52] J. P. Perdew, K. Burke, M. Ernzerhof. Generalized gradient approximation made simple. Phys. Rev. Lett., 1996, 77(18): 3865

[53] H. J. Monkhorst, J. D. Pack. Special points for Brillouin-zone integrations. Phys. Rev. B, 1976, 13(12): 5188

[54] N. Marzari, D. Vanderbilt. Maximally localized generalized Wannier functions for composite energy bands. Phys. Rev. B, 1997, 56(20): 12847

[55] I. Souza, N. Marzari, D. Vanderbilt. Maximally localized Wannier functions for entangled energy bands. Phys. Rev. B, 2001, 65(3): 035109

[56] M. P. L. Sancho, J. M. L. Sancho, J. Rubio. Quick iterative scheme for the calculation of transfer matrices: application to Mo(100). J. Phys. F Met. Phys., 1984, 14(5): 1205

[57] M. P. L. Sancho, J. M. L. Sancho, J. M. L. Sancho, J. Rubio. Highly convergent schemes for the calculation of bulk and surface Green functions. J. Phys. F Met. Phys., 1985, 15(4): 851

[58] Q. Wu, S. Zhang, H. F. Song, M. Troyer, A. A. Soluyanov. WannierTools: An open-source software package for novel topological materials. Comput. Phys. Commun., 2018, 224: 405

RIGHTS & PERMISSIONS

Higher Education Press