We see that for any $\Delta$, there is indeed a quantitative asymptotic
agreement between the bulk and surface-state spectra of the two models as the
node $(\epsilon,p_{\perp})=(0,0)$ is approached. For the Luttinger model, the
spectra are exactly quadratic in $p_{\perp}$, while for the Kane model, the
spectra have the same quadratic asymptotic behavior. Quantitative agreement
occurs within the validity range $|\epsilon|\ll|\Delta|$, $p_{\perp}\ll
p_{\Delta}$ [Eq. (72)] of the Luttinger model. We observe that the validity
range does expand upon increasing $|\Delta|$: for a fixed small enough
$p_{\perp}$, the difference between the surface-state bands of the Luttinger
(lighter solid red) and Kane (dashed red) models decreases upon increasing
$|\Delta|$.

Outside of the validity range of the Luttinger model, the spectra of the
Luttinger and Kane models deviate. The main qualitative difference is that for
the Luttinger model, the surface-state spectrum does not merge with the bulk-
band boundaries, since all scale quadratically, while the surface-state band
of the Kane model does merge with the bulk-band boundaries.

This establishes the first relation in the hierarchy of low-energy models
presented in Sec. V.

##  VIII Surfaces and symmetries

###  VIII.1 Surfaces

In the next Secs. IX, XI, and XIII, using the method of Sec. VI, we calculate
the surface states of the Dirac, line-node, and Luttinger semimetals, with
strain applied along the $z$ direction [Eq. (15)]. We do so for two
characteristic surface orientations: (i) A system occupying the half-space
$y<0$, with the $y=0$ surface parallel to the strain $z$ direction; with the
$(p_{x},p_{z})$ surface-momentum plane, for which the two Dirac points ${\bf
p}=(0,0,\pm_{u}p_{u})$ are projected onto different points
$(p_{x},p_{y})=(0,\pm_{u}p_{u})$, Fig. 1. (ii) A system occupying the half-
space $z>0$, with the $z=0$ surface perpendicular to the strain $z$ direction;
with the ${\bf p}_{\perp}=(p_{x},p_{y})$ surface-momentum plane, for which the
two Dirac points are projected onto the same point $(p_{x},p_{y})=(0,0)$, Fig.
1 [Fig. 14]. Specifying the half-space of the system for a given surface
(e.g., whether it is $z>0$ or $z<0$ for the $z=0$ surface; equivalently,
considering oriented surfaces) is important, since the chiral properties of
the surface states of the line-node and Weyl semimetals depend on that; for
the Dirac semimetal, the surface-state spectrum is the same though, due to the
symmetries, discussed below. Accordingly, the surface-state spectra will be
denoted as $\mathcal{E}_{y<0}(p_{x},p_{z})$ and
$\mathcal{E}_{z>0}(p_{x},p_{y})$, with the subscripts specifying the half-
space system and will be understood as sets of all the surface-state bands. As
with the Hamiltonians, the superscripts $L$, $D$, $LN$, and $W$ of the
surface-state spectra will be added to denote the Luttinger-, Dirac-, line-
node-, and Weyl-semimetal phases, respectively.

###  VIII.2 Symmetries

As mentioned in Secs. III.2, III.3, and III.4, the bulk Dirac semimetal has
${\bf D}_{4h}$ spatial symmetry, and the bulk line-node and Weyl semimetals
have ${\bf D}_{2d}$ spatial symmetry. For half-infinite samples, the
symmetries are deduced as follows.

Since the projected bulk-band boundaries $E_{y<0,\pm_{b}}(p_{x},p_{z})$ and
$E_{z>0,\pm_{b}}(p_{x},p_{y})$ [Eqs. (92) and (93)] are determined as extrema
of the bulk spectrum with respect to momentum component perpendicular to the
surface, their 2D symmetries are all of those of the 3D bulk system with the
perpendicular component simply dropped, regardless of whether it is changed or
not. For example, if for some $(p_{x},p_{y})$ the extremum of the bulk bands
$\varepsilon_{\pm_{b},\sigma}(p_{x},p_{y},p_{z})$ is reached at some
$p_{z}^{*}$, then, due to the $\pi$ rotation
$(p_{x},p_{y},p_{z})\rightarrow(p_{x},-p_{y},-p_{z})$ symmetry in ${\bf
D}_{2d}$ of the line-node or Weyl semimetal, the same value of the extremum is
reached at $(p_{x},-p_{y})$, which gives the effective 2D reflection symmetry
$E_{z>0,\pm_{b}}(p_{x},-p_{y})=E_{z>0,\pm_{b}}(p_{x},p_{y})$. Therefore, for
Dirac, line-node, and Weyl semimetals, the spatial symmetry of the bulk-band
boundaries is ${\bf C}_{2v}$ for the $y<0$ sample and ${\bf C}_{4v}$ for the
$z>0$ sample.

The surface-state spectrum, on the other hand, has only those symmetries of
the bulk that also preserve the sample geometry. Therefore, for the Dirac
semimetal, the spatial symmetry is ${\bf C}_{2v}$ for the $y<0$ sample and
${\bf C}_{4v}$ for the $z>0$ sample; for the line-node and Weyl semimetals,
the spatial symmetry is ${\bf C}_{2}$ for the $y<0$ sample and ${\bf C}_{2v}$
for the $z>0$ sample.

In addition to these spatial symmetries, time-reversal symmetry
$\mathcal{T}_{-}$ is still preserved for half-infinite samples, since the BCs
(86) are $\mathcal{T}_{-}$-symmetric.

##  IX Surface states of the Dirac semimetal

Figure 7:  The surface-state spectrum of the Dirac semimetal for the $y<0$
sample calculated using the “semi-analytical” method of Sec. VI within the
Luttinger model with compressive strain [$\hat{H}^{D}(\hat{{\bf p}})$, Eqs.
(22) and (86)] and absent cubic anisotropy $\alpha_{\square}=0$. The spectrum
can be presented in the universal dimensionless form
$\bar{\mathcal{E}}^{D}_{y<0}({\bar{p}}_{x},{\bar{p}}_{z};{\bar{\alpha}}_{0})$
[Eq. (109)], fully characterized by the dimensionless parameter
${\bar{\alpha}}_{0}=\alpha_{0}/\alpha_{z}$. The case ${\bar{\alpha}}_{0}=0$ of
particle-hole symmetry is presented. (a) The plane of the dimensionless
surface momentum $({\bar{p}}_{x},{\bar{p}}_{z})$ [Eq. (108)], showing the
Fermi contours (red) and the regions occupied by the surface-state band(s)
(shaded red with different intensity for the regions with one and two bands;
see, e.g., Fig. 8 for the difference), bordered by the merging contours
(purple). The projected $\pm_{u}$ Dirac points are shown in green. (b), (c),
(e), (f), (g) show the spectrum along straight-line paths, (d) and (h) show
the spectrum along circular paths, indicated in (a). In (d) and (h), the
vertical position of the horizontal axis is the Dirac-point energy
${\bar{\varepsilon}}_{u0}$. The surface-state bands
$\bar{\mathcal{E}}^{D}_{y<0}({\bar{p}}_{x},{\bar{p}}_{z};{\bar{\alpha}}_{0})$
are in red and the boundaries
${\bar{E}}^{D}_{y<0,\pm_{b}}({\bar{p}}_{x},{\bar{p}}_{z};{\bar{\alpha}}_{0})$
of the continuum (shaded light blue) of the bulk spectrum (25) are in blue.
See text for more explanations.  Figure 8:  The surface-state spectrum
$\bar{\mathcal{E}}^{D}_{y<0}({\bar{p}}_{x},{\bar{p}}_{z};{\bar{\alpha}}_{0})$
of the Dirac semimetal for ${\bar{\alpha}}_{0}=0.25$. Other conditions same as
in Fig. 7.

###  IX.1 Surfaces

In this section, we calculate the surface states of the Dirac semimetal, which
emerges when compressive strain is introduced to the Luttinger semimetal, for
the Luttinger model with the Hamiltonian $\hat{H}^{D}(\hat{{\bf p}})$ [Eq.
(22)] and BCs (86). We calculate the whole surface-state spectra
$\mathcal{E}_{y<0}^{D}(p_{x},p_{z})$ and $\mathcal{E}_{z>0}^{D}(p_{x},p_{y})$
of the Dirac semimetal with the semi-analytical method of Sec. VI for absent
cubic anisotropy, $\alpha_{\square}=0$, i.e., for the ${\bf O}(3)$ symmetry of
the Luttinger-semimetal part [Eq. (16)]. For $\alpha_{\square}=0$, due to the
axial rotation symmetry of the Dirac semimetal, the surface-state spectrum for
any surface with the normal perpendicular to the $z$ axis will be the same as
$\mathcal{E}^{D}_{y<0}(p_{x},p_{z})$ for the $y<0$ sample, and the surface-
state spectrum $\mathcal{E}^{D}_{z>0}(p_{\perp})$ for the $z>0$ sample will be
axially symmetric and depend on $p_{\perp}=\sqrt{p_{x}^{2}+p_{y}^{2}}$ only.
As in the Luttinger semimetal, including weak cubic anisotropy
$|{\bar{\alpha}}_{\square}|\ll 1$ [Eq. (106)] in the Dirac semimetal would
only lead to minor inconsequential warping (away from the borderline cases
${\bar{\alpha}}_{0}=\pm\tfrac{1}{2}$).

###  IX.2 Scaling

Figure 9:  The surface-state spectrum
$\bar{\mathcal{E}}^{D}_{y<0}({\bar{p}}_{x},{\bar{p}}_{z};{\bar{\alpha}}_{0})$
of the Dirac semimetal for ${\bar{\alpha}}_{0}=0.45$. Other conditions same as
in Fig. 7.  Figure 10:  The surface-state spectrum
$\bar{\mathcal{E}}^{D}_{y<0}({\bar{p}}_{x},{\bar{p}}_{z};{\bar{\alpha}}_{0})$
of the Dirac semimetal for the borderline case
${\bar{\alpha}}_{0}=\tfrac{1}{2}=0.5$ between the ranges
$-\tfrac{1}{2}<{\bar{\alpha}}_{0}<\tfrac{1}{2}$ and
$\tfrac{1}{2}<{\bar{\alpha}}_{0}<1$ of two and one surface-state bands,
respectively. Other conditions same as in Fig. 7.  Figure 11:  The surface-
state spectrum
$\bar{\mathcal{E}}^{D}_{y<0}({\bar{p}}_{x},{\bar{p}}_{z};{\bar{\alpha}}_{0})$
of the Dirac semimetal for ${\bar{\alpha}}_{0}=0.55$. Other conditions same as
in Fig. 7.  Figure 12:  The surface-state spectrum
$\bar{\mathcal{E}}^{D}_{y<0}({\bar{p}}_{x},{\bar{p}}_{z};{\bar{\alpha}}_{0})$
of the Dirac semimetal for ${\bar{\alpha}}_{0}=886$ of HgTe (Tab. 1). Other
conditions same as in Fig. 7.

For spherical symmetry ${\bf O}(3)$ of the Luttinger-semimetal part [Eq.
(16)], the bulk and surface-state spectra of the Dirac semimetal described by
the Luttinger model with strain [Eq. (22)] are still fully characterized by
one dimensionless parameter ${\bar{\alpha}}_{0}$. The presence of strain does
introduce energy $|u|$ and momentum

$p_{u0}=\sqrt{\frac{|u|}{\alpha_{z}}}$ (107)

scales, where the latter is the Dirac-point momentum $p_{u}$ [Eq. (23)] for
absent cubic anisotropy $\alpha_{\square}=0$; however, since there are no
other scales, its effect is that the spectra universally depend on the
relative dimensionless momentum

$\bar{\bf p}=({\bar{p}}_{x},{\bar{p}}_{y},{\bar{p}}_{z})=\frac{{\bf
p}}{p_{u0}}.$ (108)

The surface-state spectra can always be presented in the form

$\mathcal{E}^{D}_{y<0}(p_{x},p_{z};\alpha_{0},\alpha_{z};u)=|u|\bar{\mathcal{E}}^{D}_{y<0}({\bar{p}}_{x},{\bar{p}}_{z};{\bar{\alpha}}_{0}),$
(109)
$\mathcal{E}^{D}_{z>0}(p_{\perp};\alpha_{0},\alpha_{z};u)=|u|\bar{\mathcal{E}}^{D}_{z>0}({\bar{p}}_{\perp};{\bar{\alpha}}_{0})$
(110)

of the universal dimensionless functions
$\bar{\mathcal{E}}^{D}_{y<0}({\bar{p}}_{x},{\bar{p}}_{z};{\bar{\alpha}}_{0})$
and $\bar{\mathcal{E}}^{D}_{z>0}({\bar{p}}_{\perp};{\bar{\alpha}}_{0})$ of the
dimensionless surface momenta, fully characterized by just one dimensionless
curvature parameter ${\bar{\alpha}}_{0}$.

We stress that this exact scaling property of the surface-state spectrum of
the Dirac semimetal originating from the Luttinger semimetal under compressive
strain is specifically due to the scaleless property of the Luttinger model:
that its Hamiltonian (16) contains only the terms quadratic in momentum and
that the considered BCs (86) also have no scale (see discussion in Sec. VII).
For models that include higher-energy scales, like the Kane model (Sec. IV),
this is no longer exact, due to the presence of higher energy $\Delta$ and
momentum $p_{\Delta}$ scales. However, as long as the strain value $|u|$ is
within the validity range of the Luttinger model, $|u|\ll|\Delta|$, this
scaling property remains accurate as well.

###  IX.3 Asymptotic relation to the Luttinger semimetal

As the next general property, valid for any surface, at surface momentum much
larger than $p_{u}$, the surface-state spectra asymptotically approach those
[Eq. (100)] of the Luttinger semimetal:

$\bar{\mathcal{E}}^{D}_{y<0}({\bar{p}}_{x},{\bar{p}}_{z};{\bar{\alpha}}_{0})={\bar{\mathcal{A}}}({\bar{\alpha}}_{0})({\bar{p}}_{x}^{2}+{\bar{p}}_{z}^{2})+o({\bar{p}}_{x}^{2}+{\bar{p}}_{z}^{2}),\mbox{
}\sqrt{{\bar{p}}_{x}^{2}+{\bar{p}}_{z}^{2}}\gg 1,$ (111)
$\bar{\mathcal{E}}^{D}_{z>0}({\bar{p}}_{x},{\bar{p}}_{y};{\bar{\alpha}}_{0})={\bar{\mathcal{A}}}({\bar{\alpha}}_{0})({\bar{p}}_{x}^{2}+{\bar{p}}_{y}^{2})+o({\bar{p}}_{x}^{2}+{\bar{p}}_{y}^{2}),\mbox{
}\sqrt{{\bar{p}}_{x}^{2}+{\bar{p}}_{y}^{2}}\gg 1,$ (112)

where ${\bar{\mathcal{A}}}({\bar{\alpha}}_{0})$ denotes the set of curvatures
(101) of one or two surface-state bands, depending on ${\bar{\alpha}}_{0}$. We
exemplify this asymptotic relation for $y<0$ and $z>0$ samples, but the same,
of course, holds for any surface orientation. This is a manifestation of the
general annouced property that the low-energy perturbation, in this case,
strain, affects the spectrum of both the bulk and surface states only at lower
momentum, while at larger momentum, $p\gg p_{u}$, where the quadratic terms of
$\hat{H}^{L}(\hat{{\bf p}})$ in $\hat{H}^{D}(\hat{{\bf p}})$ [Eq. (22)]
dominate over the constant strain term, the spectrum asympotically approaches
that of the Luttinger semimetal.

###  IX.4 Surface states for $y<0$ sample

We describe the spectra
$\bar{\mathcal{E}}^{D}_{y<0}({\bar{p}}_{x},{\bar{p}}_{z};{\bar{\alpha}}_{0})$
by the 2D plots presented in Figs. 7-12 for the six parameter values
${\bar{\alpha}}_{0}=0,0.25,0.45,0.5,0.55,0.886$, respectively. Just like for
the Luttinger semimetal, the surface-state spectrum of the Dirac semimetal is
qualitatively different in the regimes $|{\bar{\alpha}}_{0}|<\tfrac{1}{2}$ and
$\tfrac{1}{2}<|{\bar{\alpha}}_{0}|<1$. We plot the surface-state spectra
$\bar{\mathcal{E}}^{D}_{y<0}({\bar{p}}_{x},{\bar{p}}_{z};{\bar{\alpha}}_{0})$
for the following values of ${\bar{\alpha}}_{0}$: the particle-hole symmetric
case ${\bar{\alpha}}_{0}=0$ (Fig. 7) and ${\bar{\alpha}}_{0}=0.25$ (Fig. 8)
and $0.45$ (Fig. 9), as three representative cases of the range
$|{\bar{\alpha}}_{0}|<\tfrac{1}{2}$; the borderline case
${\bar{\alpha}}_{0}=\tfrac{1}{2}=0.5$ (Fig. 10) between the two regimes; the
cases ${\bar{\alpha}}_{0}=0.55$ (Fig. 11) and $0.886$ (Fig. 12) (where the
latter is the value for HgTe), as two representative cases of the range
$\tfrac{1}{2}<|{\bar{\alpha}}_{0}|<1$. The cases ${\bar{\alpha}}_{0}=0.45$ and
${\bar{\alpha}}_{0}=0.55$, being close to the borderline case
${\bar{\alpha}}_{0}=\tfrac{1}{2}$, help better understand the transition
between the two regimes.

In each of Figs. 7-12, (a) shows the surface-momentum plane
$({\bar{p}}_{x},{\bar{p}}_{z})$, (b),(c),(e),(f), and (g) show the spectra
along the straight-line paths, and (d) and (h) show the spectra along the
circles
$({\bar{p}}_{x},{\bar{p}}_{z})=({\bar{p}}_{r}\sin\zeta,1+{\bar{p}}_{r}\cos\zeta)$,
$\zeta\in[0,2\pi)$, of radii ${\bar{p}}_{r}=0.2,0.5$, respectively, centered
at the projected Dirac point $({\bar{p}}_{x},{\bar{p}}_{z})=(0,1)$. The panel
(a) with the surface-momentum plane shows the following: regions (shaded red)
represent one or two surface-state bands according to the intensity, bounded
by the merging contours (purple), where the surface-state bands merge with the
bulk-band boundaries
${\bar{E}}^{D}_{y<0,\pm_{b}}({\bar{p}}_{x},{\bar{p}}_{z};{\bar{\alpha}}_{0})$,
labeled accordingly with $\pm_{b}$. Also, red curves show the fixed-energy
contours of the surface-state spectrum in the surface-momentum
$({\bar{p}}_{x},{\bar{p}}_{z})$ plane, obtained as crossing

$\bar{\mathcal{E}}^{D}_{y<0}({\bar{p}}_{x},{\bar{p}}_{z};{\bar{\alpha}}_{0})={\bar{\varepsilon}}_{u0}$

with the level [Eq. (24)]

${\bar{\varepsilon}}_{u0}=\frac{\varepsilon_{u}|_{\alpha_{\square}=0}}{|u|}={\bar{\alpha}}_{0}$

of the Dirac points (for ${\bar{\alpha}}_{\square}=0$) in dimensionless units.
Assuming the Fermi level is at that Dirac-point energy, these fixed-energy
contours are often called Fermi contours.

In the range $|{\bar{\alpha}}_{0}|<\tfrac{1}{2}$, there are two surface-state
bands, which, however, exist not at all momenta. The two surface-state bands
are present everywhere except for finite-size regions in-between the Dirac
points, enclosed by the merging contours with the $\pm_{b}$ bulk-band
boundaries that pass through both Dirac points ${\bar{p}}_{z}=\pm 1$. At
${\bar{\alpha}}_{0}=0$, due to particle-hole symmetry, the $\pm_{b}$ merging
contours are the same [Fig. 7(a)]; for any other
$0<|{\bar{\alpha}}_{0}|<\tfrac{1}{2}$ [Fig. 8(a) and Fig. 9(a)], they differ.
At large momentum, the two surface-state bands asymptotically do recover those
of the Luttinger semimetal, as expected according to Eq. (111).

In the range $\tfrac{1}{2}<{\bar{\alpha}}_{0}<1$, on the other hand, there is
only one surface-state band, which exists everywhere, except, interestingly,
for one point $({\bar{p}}_{x},{\bar{p}}_{z})=(0,0)$, where it touches the
$-_{b}$ bulk-band boundary. The band also connects to both Dirac points. In
this range $\tfrac{1}{2}<{\bar{\alpha}}_{0}<1$, at large momenta, this one
surface-state band also asymptotically recovers that of the Luttinger
semimetal [Eq. (111)].

The transition between the two regimes occurs continuously as follows: upon
increasing ${\bar{\alpha}}_{0}$ in the range
$0<{\bar{\alpha}}_{0}<\tfrac{1}{2}$, the region of absent upper surface-state
band expands, and eventually disappears at infinity at
${\bar{\alpha}}_{0}=\tfrac{1}{2}$. The region of absent lower surface-state
band shrinks along the ${\bar{p}}_{x}$ direction, so that at
${\bar{\alpha}}_{0}=\tfrac{1}{2}$ the two parts of the merging contour at
${\bar{p}}_{x}\gtrless 0$ collapse onto one straight-line segment at
${\bar{p}}_{x}=0$, occupying $-1<{\bar{p}}_{z}<1$, connecting the projected
Dirac points. Upon increasing ${\bar{\alpha}}_{0}=\tfrac{1}{2}+0$
infinitesimally, the now only surface-state band detaches from this line
everywhere except for the point $({\bar{p}}_{x},{\bar{p}}_{z})=(0,0)$.

Figure 13:  The surface-state spectrum
$\bar{\mathcal{E}}^{D}_{y<0}({\bar{p}}_{r}\sin\zeta,{\bar{p}}_{r}\cos\zeta;{\bar{\alpha}}_{0}=\tfrac{1}{2})$
(lighter red) of the Dirac semimetal along the circular paths of larger radii
${\bar{p}}_{r}=3$ [(a) and (b)] and $5$ [(c) and (d)], centered at
$({\bar{p}}_{x},{\bar{p}}_{z})=(0,0)$, for the borderline case
${\bar{\alpha}}_{0}=\tfrac{1}{2}$. Other conditions same as in Fig. 7. The
spectrum fits very well to the trial function (113) (dashed red), suggesting a
possible analytical solution.

The borderline case ${\bar{\alpha}}_{0}=\tfrac{1}{2}$ that continuously
connects these two regimes is itself particularly interesting. For the
Luttinger semimetal obtained from the Luttinger model, one surface-state band
merges with the bulk continuum, while the curvature
${\bar{\mathcal{A}}}_{-}({\bar{\alpha}}_{0}=\tfrac{1}{2})=0$ of the other
quadratic band vanishes and it becomes completely flat [Fig. 5(e)]. For the
Dirac semimetal obtained from the Luttinger model, we believe that, based on
the general asymptotic argument, one can only conclude that the surface-state
energy
$\bar{\mathcal{E}}_{y<0}^{D}({\bar{p}}_{x},{\bar{p}}_{z};{\bar{\alpha}}_{0}=\tfrac{1}{2})=o({\bar{p}}_{x}^{2}+{\bar{p}}_{z}^{2})$
grows slower than a quadratic power in this limit, as written in Eqs. (111)
and (112). The actual calculation shows that the limit
$\bar{\mathcal{E}}_{y<0}^{D}({\bar{p}}_{x},{\bar{p}}_{z};{\bar{\alpha}}_{0}=\tfrac{1}{2})$
at large $\sqrt{{\bar{p}}_{x}^{2}+{\bar{p}}_{z}^{2}}$ is finite, but depends
on the direction of $({\bar{p}}_{x},{\bar{p}}_{z})$. In Fig. 13, we plot the
surface-state spectrum
$\bar{\mathcal{E}}_{y<0}^{D}({\bar{p}}_{r}\sin\zeta,{\bar{p}}_{r}\cos\zeta;{\bar{\alpha}}_{0}=\tfrac{1}{2})$
along the circles of relatively large radii ${\bar{p}}_{r}=3$ and $5$ centered
at $({\bar{p}}_{x},{\bar{p}}_{z})=(0,0)$. It appears to fit very accurately
the trial function

$\tfrac{1}{2}+\tfrac{3}{4}(\cos 2\zeta-1).$ (113)

This suggests that this could be the exact next-order term of the large-
momentum expansion [the leading term of which is the Luttinger-model result
Eqs. (111) and (112)], which could perhaps be demonstrated analytically, e.g.,
via a variant of perturbation theory. We also observe that the spectrum along
the directions
$\bar{\mathcal{E}}^{D}_{y<0}(0,{\bar{p}}_{z};{\bar{\alpha}}_{0}=\tfrac{1}{2})\equiv\tfrac{1}{2}={\bar{\varepsilon}}_{u0}$
(for $|{\bar{p}}_{z}|>1$) [Fig. 10(b)] and
$\bar{\mathcal{E}}^{D}_{y<0}({\bar{p}}_{x},0;{\bar{\alpha}}_{0}=\tfrac{1}{2})\equiv-1$
[Fig. 10(e)] is exactly flat.

Regarding the Fermi contours [panel (a) of Figs. 7-12], in the whole range
$-\tfrac{1}{2}<{\bar{\alpha}}_{0}<\tfrac{1}{2}$, the two surface-state bands
never cross the Dirac-point energy level, and hence, there are no such Fermi
contours. In the range $\tfrac{1}{2}<{\bar{\alpha}}_{0}<1$, the one surface-
state band does cross the Dirac-point energy level and therefore, there is a
Fermi contour. There are two symmetric smooth parts of the contour in the
${\bar{p}}_{x}\gtrless 0$ regions, each connected to both Dirac points. We
discuss these properties below in Sec. IX.6, in relation to previous works.

The qualitatively different situations in the regimes
$|{\bar{\alpha}}_{0}|<\tfrac{1}{2}$ and $\tfrac{1}{2}<|{\bar{\alpha}}_{0}|<1$
is clearly manifested in the spectra along circular paths around the projected
Dirac point [panels (d) and (h) of Figs. 7-12 and Fig. 15], which further
confirm the general properties of the surface-state bands (in particular,
their merging behavior). For $|{\bar{\alpha}}_{0}|<\tfrac{1}{2}$, there are
two surface-state bands, both occupy finite segments of the circle and merge
with respective bulk-band boundaries. For
$\tfrac{1}{2}<|{\bar{\alpha}}_{0}|<1$, there is one surface-state band, which
exists on the whole circle and has no merging points. The transition between
these regimes occurs at ${\bar{\alpha}}_{0}=\tfrac{1}{2}$ as follows: one band
shrinks and disappears as its two merging points meet; the other band extends
and detaches from the bulk-band boundary as its two merging points meet.

Due to the asymptotic linear scaling of the surface-state spectrum in the
vicinity of the Dirac points, for small enough ${\bar{p}}_{r}$, the surface-
state spectrum

$\bar{\mathcal{E}}^{D}_{y<0}({\bar{p}}_{r}\sin\zeta,1+{\bar{p}}_{r}\cos\zeta;{\bar{\alpha}}_{0})={\bar{\varepsilon}}_{u0}+{\bar{\mathcal{V}}}^{D}_{y<0}(\zeta){\bar{p}}_{r}+o({\bar{p}}_{r}),\mbox{
}{\bar{p}}_{r}\rightarrow 0,$

is fully characterized by its velocity ${\bar{\mathcal{V}}}^{D}_{y<0}(\zeta)$
as a function of the angle $\zeta$. Therefore, the surface-state spectrum
along one circular path of small enough radii, being proportional to these
velocities, is already sufficient to fully characterize the whole local
surface-state spectrum as a function of 2D surface momentum
$({\bar{p}}_{x},{\bar{p}}_{z})$. In particular, it captures the merging
behavior with the bulk bands and the topological properties in the vicinity of
the node, as will be demonstrated in detail in Sec. XIII.

We see that the behaviors of the surface-state spectrum (controlled by a
single dimensionless parameter ${\bar{\alpha}}_{0}$) at large momenta [Eqs.
(111) and (112)], where it asymptotically approaches that of the Luttinger
semimetal, and in the vicinity of the Dirac points are related in a consistent
way. There is a qualitatively different behavior in the regimes
$|{\bar{\alpha}}_{0}|<\tfrac{1}{2}$ and $\tfrac{1}{2}<|{\bar{\alpha}}_{0}|<1$
in both limits of momentum.

###  IX.5 Surface states for $z>0$ sample

Figure 14:  The surface-state spectrum
$\bar{\mathcal{E}}^{D}_{z>0}({\bar{p}}_{\perp};{\bar{\alpha}}_{0})$ [Eq.
(110)] of the Dirac semimetal for the $z>0$ sample for the six cases
${\bar{\alpha}}_{0}=0,0.25,0.45,0.5,0.55,0.886$. Other conditions same as in
Fig. 7 (including absent cubic anisotropy $\alpha_{\square}=0$). The two Dirac
points are projected onto one point ${\bar{p}}_{\perp}=0$. In the vicinity of
${\bar{p}}_{\perp}=0$, the surface-state bands
$\bar{\mathcal{E}}^{D}_{z>0}({\bar{p}}_{\perp};{\bar{\alpha}}_{0})$ (always
two at smaller $p_{\perp}$) have a quadratic asymptotic behavior, while the
bulk-band boundaries
${\bar{E}}^{D}_{z>0,\pm_{b}}({\bar{p}}_{\perp};{\bar{\alpha}}_{0})$ have a
linear one.

In Fig. 14, we present the surface states
$\bar{\mathcal{E}}^{D}_{z>0}({\bar{p}}_{\perp};{\bar{\alpha}}_{0})$ of the
Dirac semimetal for the $z>0$ sample for the six cases
${\bar{\alpha}}_{0}=0,0.25,0.45,0.5,0.55,0.886$. For the ${\bf O}(3)$ symmetry
of the Luttinger-semimetal part ($\alpha_{\square}=0$), the surface-state
spectrum has axial symmetry and depends only on the absolute value
${\bar{p}}_{\perp}$ (including moderate cubic anisotropy $\alpha_{\square}$
would only lead to some warping of the spectrum). The bulk-band boundaries
${\bar{E}}^{D}_{z>0,\pm_{b}}({\bar{p}}_{\perp};{\bar{\alpha}}_{0})$ are
asymptotically linear in the vicinity of the projected double Dirac point at
${\bar{p}}_{\perp}=0$. We find that in the whole Luttinger-semimetal regime
$-1<{\bar{\alpha}}_{0}<1$, the structure of the surface-state spectrum near
${\bar{p}}_{\perp}=0$ is the same: there are two surface-state bands in the
vicinity of the projected Dirac point with quadratic asymptotic behavior.
Since at large momentum the behavior has to agree with that of the Luttinger
semimetal [Eq. (112)], the behavior is different at intermediate momenta
${\bar{p}}_{\perp}\sim 1$ in the two regimes. For
$|{\bar{\alpha}}_{0}|<\tfrac{1}{2}$, the two surface-state bands originating
from the Dirac point evolve into the two surface-state bands of the Luttinger
semimetal at large momenta. For $\tfrac{1}{2}<|{\bar{\alpha}}_{0}|<1$, one of
the bands originating from the Dirac point merges with the bulk-band boundary
at intermediate momentum ${\bar{p}}_{\perp}\sim 1$ [as seen in Fig. 14(k)],
while the other band evolves into the one surface-state band of the Luttinger
semimetal at large momenta. In the borderline case
${\bar{\alpha}}_{0}=\frac{1}{2}=0.5$, one of the bands is exactly flat [Fig.
14(g),(h)].

###  IX.6 Discussion

We now discuss the above-established properties of the surface states of the
Dirac semimetal that originates from the Luttinger semimetal under compressive
strain, calculated within the Luttinger model, in relation to the previous
works.

Surface states of Dirac semimetals are often compared to those of Weyl
semimetals. The topological properties of Weyl semimetals are directly tied to
those of quantum anomalous Hall systems and their surface states obey quantum
Hall topology via bulk-boundary correspondence. We discuss the topological
properties of the Dirac, line-node, and Weyl semimetal in more detail in Sec.
XIII, after we introduce the necessary concepts. In the vicinity of a single
projected Weyl point, the surface-state band is guaranteed to exist, since the
band is chiral along the path enclosing the Weyl point and its chirality is
equal to the Chern number $\pm 1$ of the Weyl point, as has been recently
shown for the most general linear-in-momentum model [54]. Hence, the Fermi
contour of the surface-state band (its fixed-energy cut, called so assuming
the Fermi energy is at the Weyl-point level) is also always guaranteed to
exist in the vicinity of a single projected Weyl point.

A Dirac point has the structure of two related Weyl points with opposite Chern
numbers $\pm 1$ and has the net zero $+1-1=0$ Chern number. Based on this, one
can conclude that surface states in Dirac semimetals are not enforced to be
present, at least not by quantum Hall topology, since a Dirac semimetal is a
topologically trivial system in terms of quantum Hall topology. Nonetheless,
in various theoretical models [9], surface-state bands are typically present
in the vicinity of Dirac points. Here, we also see that surface-state band(s)
are present generically [at least for the used BCs (86)] in the whole
Luttinger-semimetal regime $-1<{\bar{\alpha}}_{0}<1$.

By analogy with the surface-state Fermi contours in Weyl semimetals, which are
guaranteed to exist, in Refs. 56, 57, a particular attention was paid to the
presence or absence of the surface-state Fermi contours in Dirac semimetals,
especially in the vicinity of Dirac points. In one [56] presented four-band
lattice model of a Dirac semimetal, Fermi contours existed, but were not
connected to the Dirac points; in the other [57] eight-band model, Fermi
contours were completely absent.

Here, we present a model (Luttinger model with strain) with the minimal number
of bands (four, double-degenerate valence and conduction bands) that does not
have surface states at the Dirac-point energy at all (for the $y<0$ sample,
for the strain along the $z$ direction), neither connected to nor disconnected
from the Dirac points, in the substantial range
$-\tfrac{1}{2}<{\bar{\alpha}}_{0}<\tfrac{1}{2}$ (half) of the parameter space
$-1<{\bar{\alpha}}_{0}<1$. We stress that, nonetheless, two surface-state
bands do exist in this regime and they are also connected to both Dirac
points; they just never cross the Dirac-point energy level. (Therefore, one
should always make it clear whether one is talking about the presence or
absence of the surface-state bands altogether, or just whether they cross the
Fermi level or not.) We also fully characterize the evolution with
${\bar{\alpha}}_{0}$ between the regimes of two surface-state bands with
absent Fermi contours ($-\tfrac{1}{2}<{\bar{\alpha}}_{0}<\tfrac{1}{2}$) and
one surface-state band with a present Fermi contour
($\tfrac{1}{2}<|{\bar{\alpha}}_{0}|<1$). The circular-path plots (Fig. 15)
make it particularly clear how the Fermi contours disappear: the band simply
moves from crossing to not crossing the Dirac-point energy level; in the
borderline case ${\bar{\alpha}}_{0}=\tfrac{1}{2}$ between the two regimes, the
extremum of the band is exactly at that level. We see that, at least in this
model, this evolution also coincides with the other surface-state band
disappearing and the band in question detaching from the bulk-band boundary.

Therefore, the absence of the surface-state Fermi contours connected to the
Dirac points (which, again, does not at all rule out the presence of the
surface-state bands) does not seem to be rare at all and could be expected in
more systems and models (especially the behavior in the vicinity of the Dirac
points seems to be a more generic low-energy feature than the model from which
it was been derived).

##  X Surface states of the Dirac semimetal within the linear-in-momentum
model

In this section, we calculate the surface states of the Dirac semimetal in the
vicinity of the Dirac points within the linear-in-momentum model. The
Hamiltonian (31) for the latter has already been derived in Sec. III.2.1.
Here, we first derive the BCs for the linear-in-momentum model and then
analytically calculate the surface-state spectrum. The main goal is to
demonstrate the expected asymptotic agreement in the vicinity of the Dirac
points with the surface states calculated from the Luttinger model in Sec. IX.

###  X.1 Boundary conditions for the linear-in-momentum model of the Dirac
semimetal

####  X.1.1 Systematic procedure

We want to derive the BCs for the wave function $\hat{\Psi}({\bf r})$ [Eq.
(28)] of the linear-in-momentum low-energy model [Eq. (31)] of the Dirac
semimetal from the BCs (86) of the Luttinger model. This is another instance
of the general systematic low-energy-expansion procedure [8, Samokhin, 53, 46,
41, 55] of deriving BCs from the underlying “microscopic” models, which for
the system in question is as follows.

For any surface, consider the stationary Schrödinger equation

$\hat{H}^{D}(\hat{{\bf p}})\hat{\psi}({\bf r})=\varepsilon_{u}\hat{\psi}({\bf
r})$ (114)

for the Luttinger Hamiltonian of the Dirac semimetal at the Dirac-point energy
$\epsilon=\varepsilon_{u}$ [Eq. (24)]. There are eight particular solutions
with real momenta ${\bf p}=(0,0,\pm_{u}p_{u})$, the general linear
combinations of which reads

$\hat{\psi}^{\text{bulk}}(z)=\hat{\Psi}_{+_{u}}e^{+_{u}ip_{u}z}+\hat{\Psi}_{-_{u}}e^{-_{u}ip_{u}z},$
(115)

where $\hat{\Psi}_{\pm_{u}}$ are arbitrary constant four-component vectors.
These vectors represent the low-energy wave function $\hat{\Psi}({\bf r})$
[Eq. (28)].

For a given surface, we look for a general solution to Eq. (114) that does not
grow exponentially into the bulk. Such general solution is a linear
combination of particular solutions with real momenta for the directions along
the surface, but with generally complex momenta for the directions
perpendicular to the surface. Among these particular solutions, there will
necessarily be the above plane-wave bulk solutions (115) with real momenta,
representing the low-energy wave function. However, for some surfaces, there
will be additional particular solutions, exponentially decaying into the bulk
over spatial scales much shorter that the large spatial scales of interest
considered in the low-energy model. Applying the BCs (86) of the original
model to such general solution and excluding the coefficients of the latter
decaying solutions, we obtain the relations for the constants
$\hat{\Psi}_{\pm_{u}}$, which represent the sought BCs for the low-energy wave
function.

Below we apply this procedure specifically to the $z>0$ and $y<0$ samples.

####  X.1.2 $z>0$ sample

For the $z>0$ sample, the Dirac points $(0,0,\pm_{u}p_{u})$ are projected onto
the same momentum $(p_{x},p_{y})=(0,0)$ in the surface-momentum plane. We look
for the general solution $\hat{\psi}(z)$ to

$\hat{H}^{D}(0,0,\hat{p}_{z})\hat{\psi}(z)=\varepsilon_{u}\hat{\psi}(z)$

that does no grow exponentially into the bulk, as $z\rightarrow+\infty$. The
Hamiltonian $\hat{H}^{D}(0,0,\hat{p}_{z})$ is diagonal and we immediately find
that the bulk plane-wave solutions [Eq. (115)] at the nodal points
$p_{z}=\pm_{u}p_{u}$ are the only particular solutions. So, the general
solution reads

$\hat{\psi}(z)=\hat{\psi}^{\text{bulk}}(z)=\hat{\Psi}_{+_{u}}e^{+ip_{u}z}+\hat{\Psi}_{-_{u}}e^{-ip_{u}z}.$

Applying the BCs (86), we obtain

$\hat{\psi}(z=0)=\hat{\Psi}_{+_{u}}+\hat{\Psi}_{-_{u}}=\hat{0}.$

To leading order, these relations hold for present coordinate dependence
$\hat{\Psi}({\bf r})$ of the low-energy wave function. Therefore, upon the
substitution $\hat{\Psi}_{\pm_{u}}\rightarrow\hat{\Psi}_{\pm_{u}}(x,y,z=0)$,
these represent the BCs

$\hat{\Psi}_{+_{u}}(x,y,z=0)+\hat{\Psi}_{-_{u}}(x,y,z=0)=\hat{0}$ (116)

for the low-energy wave function of the linear model of the Dirac semimetal
for $z>0$ sample. We see that different $j_{z}$ components are not coupled by
such BCs, while the two Dirac points $\pm_{u}$ for each $j_{z}$ component are
coupled.

####  X.1.3 $y<0$ sample

For the $y<0$ sample, the Dirac points $(0,0,\pm_{u}p_{u})$ are projected onto
the different momenta $(p_{x},p_{z})=(0,\pm_{u}p_{u})$ in the surface-momentum
plane. We look for the general solution of the form

$\hat{\psi}(y,z)=\hat{\psi}_{+_{u}}(y)e^{+_{u}ip_{u}z}+\hat{\psi}_{-_{u}}(y)e^{-_{u}ip_{u}z}$
(117)

to the Schrödinger equation at $\epsilon=\varepsilon_{u}$ that does not grow
exponentially into the bulk, as $y\rightarrow-\infty$. Due to the linear
independence of the functions $e^{\pm_{u}ip_{u}z}$, we have the decoupled
equations

$\hat{H}^{D}(0,\hat{p}_{y},\pm_{u}p_{u})\hat{\psi}_{\pm_{u}}(y)=\varepsilon_{u}\hat{\psi}_{\pm_{u}}(y)$
(118)

for the functions $\hat{\psi}_{\pm_{u}}(y)$.

In this case, the situation is technically more complicated. In the
Hamiltonians $\hat{H}^{D}(0,\hat{p}_{y},\pm_{u}p_{u})$ with nonzero
$\hat{p}_{y}$, i.e., for coordinate-dependent $\hat{\psi}_{\pm_{u}}(y)$, all
four wave-function components with definite $j_{z}$ are coupled. For each
Dirac point $\pm_{u}$, we find that there are six momentum $p_{y}$ solutions
to the characteristic equation

$\det[\hat{H}^{D}(0,p_{y},\pm_{u}p_{u})-\varepsilon_{u}\hat{1}_{4}]=0,$

for which the particular solutions are not growing exponentially into the
bulk: there is a four-fold-degenerate solution $p_{y}=0$, representing the
bulk low-energy wave function, and a double-degenerate solution

$p_{y}=-i\kappa,\mbox{
}\kappa=2p_{u}\frac{\frac{\sqrt{3}}{2}(\alpha_{z}+\tfrac{2}{5}\alpha_{\square})}{\sqrt{(\alpha_{z}-\frac{3}{5}\alpha_{\square})^{2}-\alpha_{0}^{2}}}=2p_{u}\frac{v_{\perp}}{\sqrt{v_{z}^{2}-v_{0}^{2}}}$

for which the particular solutions decay into the bulk. The corresponding
general solution to Eq. (118) not growing into the bulk reads

$\hat{\psi}_{\pm_{u}}(y)=\left(\begin{array}[]{c}\Psi_{\pm_{u},+\frac{3}{2}}\\\
\Psi_{\pm_{u},+\frac{1}{2}}\\\ \Psi_{\pm_{u},-\frac{1}{2}}\\\
\Psi_{\pm_{u},-\frac{3}{2}}\end{array}\right)+\left[b_{\pm_{u},1}\left(\begin{array}[]{c}\pm_{u}\sqrt{v_{z}^{2}-v_{0}^{2}}\\\
v_{0}+\tfrac{1}{2}v_{z}\\\ 0\\\
-\frac{\sqrt{3}}{2}v_{z}\end{array}\right)+b_{\pm_{u},2}\left(\begin{array}[]{c}v_{0}-\tfrac{1}{2}v_{z}\\\
\mp_{u}\sqrt{v_{z}^{2}-v_{0}^{2}}\\\ -\frac{\sqrt{3}}{2}v_{z}\\\
0\end{array}\right)\right]e^{\kappa y}.$ (119)

Here, $\Psi_{\pm_{u},j_{z}}$,
$j_{z}=+\frac{3}{2},+\frac{1}{2},-\frac{1}{2},-\frac{3}{2}$, and
$b_{\pm_{u},1}$, $b_{\pm_{u},2}$ are the six constant free coefficients.
Applying the BCs (86) to the general solution (117), we have

$\hat{\psi}(y=0,z)=\hat{0}.$

These BCs have to be satisfied at every point $(x,z)$ on the $y=0$ surface.
Due to the linear independence of $e^{\pm_{u}ip_{u}z}$, the BCs have to be
satisfied independently for the functions (119):

$\hat{\psi}_{\pm_{u}}(y=0)=\hat{0}.$

For each Dirac point $\pm_{u}$, this gives four linear homogeneous relations
between the six coefficients. Excluding $b_{\pm_{u},1}$ and $b_{\pm_{u},2}$
from these relations, we arrive at the two relations

$\displaystyle\tfrac{\sqrt{3}}{2}v_{z}\Psi_{\pm_{u},+\frac{3}{2}}+(v_{0}-\tfrac{1}{2}v_{z})\Psi_{\pm_{u},-\frac{1}{2}}\pm_{u}\sqrt{v_{z}^{2}-v_{0}^{2}}\Psi_{\pm_{u},-\frac{3}{2}}=0,$
(120)
$\displaystyle\tfrac{\sqrt{3}}{2}v_{z}\Psi_{\pm_{u},+\frac{1}{2}}+(v_{0}+\tfrac{1}{2}v_{z})\Psi_{\pm_{u},-\frac{3}{2}}\mp_{u}\sqrt{v_{z}^{2}-v_{0}^{2}}\Psi_{\pm_{u},-\frac{1}{2}}=0$
(121)

for the four coefficients. These become the BCs for the low-energy wave
function upon the identification

$\hat{\Psi}_{\pm_{u}}\rightarrow\hat{\Psi}_{\pm_{u}}(x,y=0,z)$

of the coefficients with the wave-function components at the surface. It is
also instructive for further calculations of the surface states to express
these BCs in terms of the velocities $v_{0,z}$ [Eq. (33)], since these will
combine naturally with those in the linear Hamiltonian [Eq. (31)].

####  X.1.4 Summary

We summarize how the BCs for the wave function $\hat{\Psi}({\bf r})$ of the
linear model of the Dirac semimetal originate for the $z>0$ and $y<0$ samples.
This wave function has eight components and there are four BCs for any
surface, since its Hamiltonian $\hat{\mathcal{H}}^{D}(\hat{{\bf k}})$ [Eq.
(31)] is linear in momentum $\hat{{\bf k}}$.

The eight plane-wave solutions (four for each Dirac point) are present for any
surface. For the $z>0$ sample, these were the only relevant particular
solutions to the stationary Schrödinger equation for the Luttinger Hamiltonian
at the Dirac-point energy $\epsilon=\varepsilon_{u}$. All of these solutions
have the same (absent) coordinate dependence on $(x,y)$ along the $z=0$
surface; as a result, the four BCs for the wave function $\hat{\psi}({\bf r})$
of the Luttinger model translate directly to four BCs (116) for the eight
components of the low-energy wave function $\hat{\Psi}({\bf r})$.

One the other hand, for the $y<0$ sample, in addition to the eight bulk plane-
wave solutions, there are extra four (two for each Dirac point) decaying
solutions; so, there are twelve free coefficients (six for each Dirac point).
Here, the solutions for the two nodes have different dependencies
$e^{\pm_{u}ip_{u}z}$ along the $y=0$ surface; as a result, the four functional
BCs for the Luttinger-model wave function $\hat{\psi}({\bf r})$, due to their
translation symmetry, reduce to eight relations for the twelve free
coefficients. Excluding the four coefficients of the decaying solutions, we
again obtain four BCs (121) for the eight components of the low-energy wave
function $\hat{\Psi}({\bf r})$.

###  X.2 Analytical calculation of the surface states

We now calculate the surface states analytically in the vicinity of the Dirac
points from the linear model of the Dirac semimetal derived above.

####  X.2.1 $y<0$ sample

For energies

$\epsilon\in(E_{y<0,\pm_{u},-_{b}}^{D,\text{lin}}(k_{x},k_{z}),E_{y<0,\pm_{u},+_{b}}^{D,\text{lin}}(k_{x},k_{z}))$

within the gap between the bulk-band boundaries

$E_{y<0,\pm_{u},\pm_{b}}^{D,\text{lin}}(k_{x},k_{z})=\varepsilon_{u}\pm_{u}v_{0}k_{z}\pm_{b}\sqrt{v_{\perp}^{2}k_{x}^{2}+v_{z}^{2}k_{z}^{2}},$
(122)

the characteristic equation

$\det[\hat{\mathcal{H}}^{D}_{\pm_{u},\pm_{j_{z}}}({\bf
k})-\epsilon\hat{1}_{2}]=0$ (123)

with respect to $k_{y}$ for each $\pm_{u}$, $\pm_{j_{z}}$ block has one
solution (the same for $\pm_{j_{z}}$)

$k_{y,\pm_{u}}(\epsilon,k_{x},k_{z})=-i\sqrt{k_{x}^{2}+\frac{v_{z}^{2}k_{z}^{2}-[\pm_{u}v_{0}k_{z}-(\epsilon-\varepsilon_{u})]^{2}}{v_{\perp}^{2}}}$
(124)

that corresponds to a particular solution to the stationary Schrödinger
equation decaying into the bulk, as $y\rightarrow-\infty$. For each $\pm_{u}$,

$\hat{\Psi}^{+_{u}}(y;\epsilon,k_{x},k_{z})=\left(\begin{array}[]{c}\hat{\Psi}_{+_{u}}^{+_{u}}(y;\epsilon,k_{x},k_{z})\\\
\hat{0}_{4}\end{array}\right),\mbox{
}\hat{\Psi}^{-_{u}}(y;\epsilon,k_{x},k_{z})=\left(\begin{array}[]{c}\hat{0}_{4}\\\
\hat{\Psi}_{-_{u}}^{-_{u}}(y;\epsilon,k_{x},k_{z})\end{array}\right)$

the general solution is their linear combination

$\hat{\Psi}_{\pm_{u}}^{\pm_{u}}(y;\epsilon,k_{x},k_{z})=\left[c_{\pm_{u},+_{j_{z}}}\left(\begin{array}[]{c}\hat{\chi}^{\pm_{u},+_{j_{z}}}(\epsilon,k_{x},k_{z})\\\
\hat{0}_{2}\end{array}\right)+c_{\pm_{u},-_{j_{z}}}\left(\begin{array}[]{c}\hat{0}_{2}\\\
\hat{\chi}^{\pm_{u},-_{j_{z}}}(\epsilon,k_{x},k_{z})\end{array}\right)\right]e^{ik_{y,\pm_{u}}(\epsilon,k_{x},k_{z})y},$

where

$\hat{\chi}^{\pm_{u},+_{j_{z}}}(\epsilon,k_{x},k_{z})=\left(\begin{array}[]{c}v_{\perp}k_{x}-\sqrt{v_{\perp}^{2}k_{x}^{2}+v_{z}^{2}k_{z}^{2}-[v_{0}k_{z}-(\epsilon-\varepsilon_{u})]^{2}}\\\
(v_{0}-v_{z})k_{z}\mp_{u}(\epsilon-\varepsilon_{u})\end{array}\right),$ (125)
$\hat{\chi}^{\pm_{u},-_{j_{z}}}(\epsilon,k_{x},k_{z})=\left(\begin{array}[]{c}-v_{\perp}k_{x}+\sqrt{v_{\perp}^{2}k_{x}^{2}+v_{z}^{2}k_{z}^{2}-[v_{0}k_{z}-(\epsilon-\varepsilon_{u})]^{2}}\\\
(v_{0}+v_{z})k_{z}\mp_{u}(\epsilon-\varepsilon_{u})\end{array}\right)$ (126)

are the nontrivial solutions to

$[\hat{\mathcal{H}}^{D}_{\pm_{u},\pm_{j_{z}}}(k_{x},k_{y,\pm_{u}}(\epsilon,k_{x},k_{z}),k_{z})-\epsilon\hat{1}_{2}]\hat{\chi}=\hat{0}_{2}$

and $c_{\pm_{u},\pm_{j_{z}}}$ are arbitrary coefficients. Regarding notation,
here and below, the superscript labels of the wave functions denote the
quantum numbers, while the subscripts are reserved for their components. In
$\hat{0}_{4}$ and $\hat{0}_{2}$, the sizes of the null vectors are indicated
for clarity.

Inserting these wave functions into the BCs (121) yields the equation for the
energy $\epsilon$ of the surface states. Solving it, we obtain the surface-
state spectrum for the linear model of the Dirac semimetal

$\mathcal{E}^{D,\text{lin}}_{y<0,\pm_{u},\pm}(k_{x},k_{z})=\varepsilon_{u}+\tfrac{\sqrt{3}}{2}[\pm_{u}\sqrt{3}v_{0}k_{z}\pm\sqrt{(v_{z}^{2}-v_{0}^{2})(\tfrac{4v_{\perp}^{2}}{3v_{z}^{2}}k_{x}^{2}+k_{z}^{2})}].$
(127)

Due to the linear scaling, the surface-state spectrum (127), as well as the
bulk-band boundaries (122), are linear in the absolute value
$\sqrt{k_{x}^{2}+k_{z}^{2}}$ of the surface momentum relative to the projected
Dirac points and are fully characterized by their dependence on the polar
angle $\zeta$,

$(k_{x},k_{z})=\sqrt{k_{x}^{2}+k_{z}^{2}}(\sin\zeta,\cos\zeta).$

The surface-state spectrum (127) is shown in Fig. 15, along with that
calculated within the Luttinger model, for comparison; their quantitative
asymptotic agreement is discussed below in Sec. X.4. The properties do fully
agree with those obtained from Luttinger model. The solution is fully
characterized by the ratio

$\frac{v_{0}}{v_{z}}=\frac{\alpha_{0}}{\alpha_{z}-\frac{3}{5}\alpha_{\square}}.$

In the Luttinger-semimetal regime from which this Dirac semimetal originates,
$-1<v_{0}/v_{z}<1$. Because of the linear scaling, the merging contours in
$(k_{x},k_{z})$ are straight half-lines with specific merging angles
$\zeta_{m}$. The latter can be obtained by equating the surface-state energy
(127) to the bulk-band boundaries (122), which gives

$\cos\zeta_{m}=\pm\frac{2\frac{v_{0}}{v_{z}}\frac{v_{\perp}}{v_{z}}}{\sqrt{1-4\frac{v_{0}^{2}}{v_{z}^{2}}(1-\frac{v_{\perp}^{2}}{v_{z}^{2}})}}.$

Depending on $v_{0}/v_{z}$, there are two regimes, with one and two surface-
state bands (127). The borderline case can be found be demanding that the
merging angle(s) $\zeta_{m}$ equal $0$ or $\pi$, which is when one of the two
bands as functions of $\zeta$ shrinks and disappears and the other one
detaches from the bulk-band boundaries. This happens at

$\frac{v_{0}}{v_{z}}=\pm\tfrac{1}{2}.$

For $-\tfrac{1}{2}<v_{0}/v_{z}<\tfrac{1}{2}$, both surface-state bands (127)
are present. For $-1<v_{0}/v_{z}<-\tfrac{1}{2}$, only the $+$ band is present;
for $\tfrac{1}{2}<v_{0}/v_{z}<1$, only the $-$ band is present.

We also observe an interesting relation: the surface-state spectrum (127) of
the linear-in-momentum model of the Dirac semimetal can be expressed (at least
for $k_{x}=0$) in terms of the curvatures (105) of the surface-state spectrum
of the Luttinger model of the Luttinger semimetal (also for finite cubic
anisotropy):

$\mathcal{E}^{D,\text{lin}}_{y<0,+_{u},\pm}(k_{x}=0,k_{z})=\varepsilon_{u}+{\bar{\mathcal{A}}}_{\pm}^{{\bf
O}_{h}}({\bar{\alpha}}_{0},{\bar{\alpha}}_{\square})\alpha_{z}\cdot
2p_{u}k_{z},\mbox{ }k_{z}>0,$

and similar for the $-_{u}$ Dirac point or $k_{z}<0$.

Figure 15:  The surface-state spectrum
$\mathcal{E}^{D,\text{lin}}_{y<0}(k_{x},k_{z})$ [Eq. (127)] of the Dirac
semimetal in the vicinity of the $+_{u}$ projected Dirac point for the $y<0$
sample, calculated analytically within the linear-in-momentum model (dashed
red) with the Hamiltonian $\hat{\mathcal{H}}^{D}(\hat{{\bf k}})$ [Eq. (31)]
and BCs (121). The spectrum is presented in the dimensionless momentum and
energy units [Eqs. (108) and (109)] and for absent cubic anisotropy
$\alpha_{\square}=0$. and the comparison with spectrum calculated within the
Luttinger model (solid lighter red). The bulk-band boundaries of the linear-
in-momentum and Luttinger models are in dashed blue and solid lighter blue,
respectively. (a),(b),(c),(e),(f),(g) are spectra along the circular paths of
radius ${\bar{p}}_{r}=0.05$; (d) and (h) are spectra along straight-line paths
passing through the projected Dirac point. Quantitative asymptotic agreement
between the two models as the projected Dirac point is approached is observed.

####  X.2.2 $z>0$ sample

Figure 16:  (a) The flat surface-state spectrum
$\mathcal{E}^{D,\text{lin}}_{z>0}(k_{\perp})\equiv\varepsilon_{u}$ [Eq. (132)]
of the Dirac semimetal for the $z>0$ sample, calculated within the linear-in-
momentum model with the Hamiltonian $\hat{\mathcal{H}}^{D}(\hat{{\bf k}})$
[Eq. (31)] and BCs (116). This flat spectrum and the quadratic asymptotic
behavior of the surface-state bands (Fig. 14) calculated within the Luttinger
model as the projected double Dirac point is approached constitute
quantitative asymptotic agreement between the two models, see Sec. X.4. (b)
The flat surface-state spectrum
$\mathcal{E}^{LN,\text{lin}}_{z>0}(k_{\perp})\equiv\varepsilon_{u}$ [Eq.
(143)] of the line-node semimetal for the $z>0$ sample, calculated within the
linear-in-momentum model with the Hamiltonian
$\hat{\mathcal{H}}^{LN}(\hat{{\bf k}})$ [Eq. (39)] and BCs (116).

For the $z>0$ sample, the two Dirac points are projected onto the same surface
momentum $(k_{x},k_{y})=(0,0)$, as a result, the $\pm_{u}$ parts
$\hat{\Psi}_{\pm_{u},\pm_{j_{z}}}({\bf r})$ are coupled by the BCs (116).
However, the $\pm_{j_{z}}$ parts are decoupled in both the Hamiltonian and
BCs.

For the energies

$\epsilon\in(E_{z>0,-_{b}}^{D,\text{lin}}(k_{\perp}),E_{z>0,+_{b}}^{D,\text{lin}}(k_{\perp}))$

within the gap between the bulk-band boundaries

$E^{D,\text{lin}}_{z>0,\pm_{b}}(k_{\perp})=\varepsilon_{u}\pm_{b}v_{\perp}k_{\perp}\sqrt{1-\tfrac{v_{0}^{2}}{v_{z}^{2}}},$

the characteristic equation (123), with respect to $k_{z}$ for each
$\pm_{u},\pm_{j_{z}}$ block has one solution (the same for $\pm_{j_{z}}$)

$k_{z,\pm_{u}}(\epsilon,k_{\perp})=\frac{\mp_{u}v_{0}(\epsilon-\varepsilon_{u})+i\sqrt{k_{\perp}^{2}v_{\perp}^{2}(v_{z}^{2}-v_{0}^{2})-v_{z}^{2}(\epsilon-\varepsilon_{u})^{2}}}{v_{z}^{2}-v_{0}^{2}}$
(128)

that corresponds to a particular solution to the stationary Schrödinger
equation decaying into the bulk, as $z\rightarrow+\infty$. The general
solution for each $\pm_{j_{z}}$ is their linear combination

$\hat{\Psi}^{+_{j_{z}}}(z;\epsilon,k_{\perp},\phi)=c_{+_{u},+_{j_{z}}}\left(\begin{array}[]{c}\hat{\chi}^{+_{u},+_{j_{z}}}(\epsilon,k_{\perp},\phi)\\\
\hat{0}_{2}\\\
\hat{0}_{4}\end{array}\right)e^{ik_{z,+_{u}}(\epsilon,k_{\perp})z}+c_{-_{u},+_{j_{z}}}\left(\begin{array}[]{c}\hat{0}_{4}\\\
\hat{\chi}^{-_{u},+_{j_{z}}}(\epsilon,k_{\perp},\phi)\\\
\hat{0}_{2}\end{array}\right)e^{ik_{z,-_{u}}(\epsilon,k_{\perp})z},$ (129)
$\hat{\Psi}^{-_{j_{z}}}(z;\epsilon,k_{\perp},\phi)=c_{+_{u},-_{j_{z}}}\left(\begin{array}[]{c}\hat{0}_{2}\\\
\hat{\chi}^{+_{u},-_{j_{z}}}(\epsilon,k_{\perp},\phi)\\\
\hat{0}_{4}\end{array}\right)e^{ik_{z,+_{u}}(\epsilon,k_{\perp})z}+c_{-_{u},-_{j_{z}}}\left(\begin{array}[]{c}\hat{0}_{4}\\\
\hat{0}_{2}\\\
\hat{\chi}^{-_{u},-_{j_{z}}}(\epsilon,k_{\perp},\phi)\end{array}\right)e^{ik_{z,-_{u}}(\epsilon,k_{\perp})z},$
(130)

where

$\hat{\chi}^{\pm_{u},+_{j_{z}}}(\epsilon,k_{\perp},\phi)=\left(\begin{array}[]{c}\mp_{u}v_{\perp}k_{\perp}e^{-i\phi}\\\
\epsilon-\varepsilon_{u}\mp_{u}(v_{0}-v_{z})k_{z,\pm_{u}}(\epsilon,k_{\perp})\end{array}\right),$
$\hat{\chi}^{\pm_{u},-_{j_{z}}}(\epsilon,k_{\perp},\phi)=\left(\begin{array}[]{c}\pm_{u}v_{\perp}k_{\perp}e^{-i\phi}\\\
\epsilon-\varepsilon_{u}\mp_{u}(v_{0}+v_{z})k_{z,\pm_{u}}(\epsilon,k_{\perp})\end{array}\right)$

are the nontrivial solutions to

$[\hat{\mathcal{H}}^{D}_{\pm_{u},\pm_{j_{z}}}({\bf
k}_{\perp},k_{z,\pm_{u}}(\epsilon,k_{\perp}))-\epsilon\hat{1}_{2}]\hat{\chi}=\hat{0}_{2}$

and $c_{\pm_{u},\pm_{j_{z}}}$ are arbitrary coefficients.

Applying the BCs (116), we obtain the homogeneous system of equations [akin to
Eq. (97)] for these coefficients, which has nontrivial solutions when its
determinant is zero [akin to Eq. (98)], which yields the equation

$\chi^{+_{u},+_{j_{z}}}_{+\frac{3}{2}}(\epsilon,k_{\perp},\phi)\chi^{-_{u},+_{j_{z}}}_{+\frac{1}{2}}(\epsilon,k_{\perp},\phi)-\chi^{-_{u},+_{j_{z}}}_{+\frac{3}{2}}(\epsilon,k_{\perp},\phi)\chi^{+_{u},+_{j_{z}}}_{+\frac{1}{2}}(\epsilon,k_{\perp},\phi)=0$
(131)

for the energy $\epsilon$ of the surface states for the $+_{j_{z}}$ block. For
every surface momentum $(k_{x},k_{y})$, the only solution to this equation is
a nondegenerate $\epsilon=\varepsilon_{u}$. Similar situation holds for the
$-_{j_{z}}$ block. Hence, the surface-state spectrum (Fig. 16) for $z>0$
sample obtained from the linear-in-momentum model of the Dirac semimetal is
flat and double-degenerate:

$\mathcal{E}^{D,\text{lin}}_{z>0,\pm_{j_{z}}}(k_{\perp})\equiv\varepsilon_{u}.$
(132)

The momentum solutions of the surface states equal

$k_{z,\pm_{u}}(\epsilon=\varepsilon_{u},k_{\perp})=i\varkappa,\mbox{
}\varkappa=\frac{v_{\perp}k_{\perp}}{\sqrt{v_{z}^{2}-v_{0}^{2}}},$

the coefficients are $c_{+_{u},\pm_{j_{z}}}=c_{-_{u},\pm_{j_{z}}}=1$, and the
wave functions of the surface states reads

$\hat{\Psi}^{+_{j_{z}}}(z;\epsilon=\varepsilon_{u},k_{\perp},\phi)=\left(\begin{array}[]{c}\hat{\chi}^{+_{u},+_{j_{z}}}(\epsilon=\varepsilon_{u},k_{\perp},\phi)\\\
\hat{0}_{2}\\\
\hat{\chi}^{-_{u},+_{j_{z}}}(\epsilon=\varepsilon_{u},k_{\perp},\phi)\\\
\hat{0}_{2}\end{array}\right)e^{-\varkappa z},\mbox{
}\hat{\chi}^{\pm_{u},+_{j_{z}}}(\epsilon=\varepsilon_{u},k_{\perp},\phi)=\left(\begin{array}[]{c}\mp_{u}v_{\perp}k_{\perp}e^{-i\phi}\\\
\mp_{u}(v_{0}-v_{z})i\varkappa\end{array}\right),$ (133)
$\hat{\Psi}^{-_{j_{z}}}(z;\epsilon=\varepsilon_{u},k_{\perp},\phi)=\left(\begin{array}[]{c}\hat{0}_{2}\\\
\hat{\chi}^{+_{u},-_{j_{z}}}(\epsilon=\varepsilon_{u},k_{\perp},\phi)\\\
\hat{0}_{2}\\\
\hat{\chi}^{-_{u},-_{j_{z}}}(\epsilon=\varepsilon_{u},k_{\perp},\phi)\end{array}\right)e^{-\varkappa
z},\mbox{
}\hat{\chi}^{\pm_{u},-_{j_{z}}}(\epsilon=\varepsilon_{u},k_{\perp},\phi)=\left(\begin{array}[]{c}\pm_{u}v_{\perp}k_{\perp}e^{-i\phi}\\\
\mp_{u}(v_{0}+v_{z})i\varkappa\end{array}\right).$ (134)

###  X.3 Emergent chiral symmetry of the linear-in-momentum model

We find a symmetry explanation for the obtained flat behavior of the surface-
state spectrum (132) of the linear model of the Dirac semimetal for the $z>0$
sample. We notice that both the Hamiltonian (31) and BCs (116) possess chiral
symmetry. Namely, we notice that the Hamiltonian (31) changes its sign
relative to the Dirac-point energy $\varepsilon_{u}$ under the unitary
transformation

$\hat{S}[\hat{\mathcal{H}}^{D}(\hat{{\bf
k}})-\varepsilon_{u}\hat{1}_{8}]\hat{S}^{\dagger}=-[\hat{\mathcal{H}}^{D}(\hat{{\bf
k}})-\varepsilon_{u}\hat{1}_{8}],$ (135)
$\hat{S}=\hat{\tau}_{x}\otimes\hat{1}_{4},$

since the blocks $\hat{\mathcal{H}}_{\pm_{u}}^{D}(\hat{{\bf
k}})-\varepsilon_{u}\hat{1}_{8}$ differ only by the sign. This means that the
Hamiltonian satisfies chiral symmetry and $\hat{S}$ is its operator [58]. The
operation $\hat{S}$ interchanges $\pm_{u}$ components of the wave function:

$\hat{S}\hat{\Psi}({\bf r})=\left(\begin{array}[]{c}\hat{\Psi}_{-_{u}}({\bf
r})\\\ \Psi_{+_{u}}({\bf r})\end{array}\right).$

This transformed wave function also satisfies the BCs (116) for the $z>0$
sample, which means that they also satisfy chiral symmetry [8, 46, 54]. Hence,
the whole $z>0$ system with the boundary has chiral symmetry. As a result, its
spectrum, including the surface states, satisfies chiral symmetry. Regarding
the double degeneracy of the flat surface-state band, for these Hamiltonian
and BCs it is guaranteed by the decoupling of the $\pm_{j_{z}}$ parts in both
of them. Moreover, we notice that both $\pm_{j_{z}}$ surface-state solutions
have the same eigenvalue $-1$ of the chiral symmetry operator:

$\hat{S}\hat{\Psi}^{\pm_{j_{z}}}(z;\epsilon=\varepsilon_{u},k_{\perp},\phi)=-\hat{\Psi}^{\pm_{j_{z}}}(z;\epsilon=\varepsilon_{u},k_{\perp},\phi).$

This means that the double degeneracy is also protected by chiral symmetry,
meaning that double degeneracy will persist if the Hamiltonian or BCs are
modified (for example, the decoupling of $\pm_{j_{z}}$ blocks could be broken)
in a way that still preserves their chiral symmetry.

For the $y<0$ sample, the BCs (121) of the linear-in-momentum model do not
obey chiral symmetry and therefore, neither does the surface-state spectrum.

We note that this is an emergent symmetry of specifically the linear-in-
momentum model for the $z>0$ sample that arises in this low-energy limit. In
the Luttinger model with strain, there is no such chiral symmetry. And indeed
the surface-state spectrum is not flat; its quadratic dispersion at small
$p_{\perp}$ is, however, in accord with the flat behavior in the linear model.

###  X.4 Asymptotic agreement between the linear-in-momentum and Luttinger
models of the Dirac semimetal

Comparing the surface-state spectra of the Dirac semimetal in the vicinity of
the Dirac points, $|{\bf p}\mp_{u}{\bf p}_{u}|\ll p_{u}$, calculated within
the linear-in-momentum and Luttinger models, Fig. 15 for the $y<0$ sample and
Figs. 14 and 16 for the $z>0$ sample (flat versus quadratic behavior), we
observe a clear quantitative asymptotic agrement, as the projected Dirac
points are approached. This establishes the next relation in the hierarchy of
low-energy models presented in Sec. V. Just as the Luttinger model is
sufficient to fully capture the asymptotic behavior of the surface states of
the Kane model in the vicinity of the quadratic node in the Luttinger-
semimetal phase (Sec. VII.2), the linear-in-momentum model is sufficient to
fully capture the asymptotic behavior of the surface states of the Luttinger
model in the vicinity of the Dirac nodes in the Dirac-semimetal phase.

##  XI Surface states of the line-node semimetal

In this section, we calculate and explore within the Luttinger model the
surface states of the line-node semimetal that arises in the Dirac semimetal
upon introducing the linear BIA term; the Dirac semimetal itself, studied in
the previous Sec. IX, arises from the Luttinger semimetal upon introducing
compressive strain. The bulk properties of this line-node semimetal have been
studied in Sec. III.3.

The Luttinger model (35) of the line-node semimetal is fully characterized by
three dimensionless parameters: two of the Luttinger-semimetal part,
${\bar{\alpha}}_{0}$, ${\bar{\alpha}}_{\square}$, and the ratio
$p_{\beta_{1}}/p_{u}$ of the characteristic momenta of the linear BIA [Eq.
(43)] and strain [Eq. (23)] terms. We present the line-node (this section) and
Weyl (next Sec. XIII) semimetal cases for the parameters
$\alpha_{0,z,\square}$ and $\beta_{1}$ of HgTe and the value $u=-3\text{meV}$
of the compressive strain. Although we include the nonzero cubic anisotropy
$\alpha_{\square}$, it does not seem to play any notable role for the regime
of other parameters we explore: we have compared to the case of absent cubic
anisotropy ($\alpha_{\square}=0$) and found only minor inessential
quantitative differences.

For the chosen value of strain, the hierarchy presented in Sec. V is satisfied
well. As demonstrated in Sec. III.3, upon introducing the linear BIA term into
the Dirac semimetal, the Dirac points of the bulk spectrum transform into line
nodes. As we show below, the surface-state spectra of both the $y<0$ and $z>0$
samples also undergo rather nontrivial transformations at the $p_{\beta_{1}}$
scale around the former projected Dirac points.

###  XI.1 $y<0$ sample

Figure 17:  The surface-state spectrum $\mathcal{E}^{LN}_{y<0}(p_{x},p_{z})$
of the line-node semimetal for the $y<0$ sample, calculated within the
Luttinger model with compressive strain and linear BIA term [Eqs. (35) and
(86)] for the parameters of HgTe (Tab. 1) and the strain value
$u=-3~{}\text{meV}$. The structure of the figure is the same as in Fig. 12 for
the Dirac semimetal. Compared to the surface-state spectrum of the Dirac
semimetal, the key changes occur in the vicinity of the (former) projected
Dirac points, which transform into projected line nodes (green). The spectrum
is presented in the dimensionless momentum and energy units [Eqs. (108) and
(109)] of the Dirac-semimetal case for better comparison.  Figure 18:  The
surface-state spectrum $\mathcal{E}^{LN}_{y<0}(p_{x},p_{z})$ (red) of the
line-node semimetal along $p_{x}=\text{const}$ straight-line paths in the
vicinity $p_{z}\approx p_{z}^{LN}(0)=p_{z0}$ of the projected $+_{u}$ line
node (additional characterization of Fig. 17). For $-p_{\perp
0}<p_{x}<p_{\perp 0}$, $p_{\perp 0}=p_{\perp}^{LN}(0)$, the path crosses the
projected line node; the bulk-band boundaries
$E^{LN}_{y<0,\pm_{b}}(p_{x},p_{z})$ (blue) along this path have a linear node
and the surface-state band $\mathcal{E}^{LN}_{y<0}(p_{x},p_{z})$ also has a
linear asymptotic behavior on both sides.  Figure 19:  Additional
characterization of the surface-spectrum $\mathcal{E}_{y<0}^{LN}(p_{x},p_{z})$
of the line-node semimetal (Fig. 17) in the vicinity of the projected $+_{u}$
line node: (a),(b),(c) along $p_{z}=\text{const}$ straight-line paths and
(d),(e),(f) along circular paths. The latter are useful for comparison with
those around the projected Weyl points of the Weyl semimetal, presented in
Figs. 23, 24, and 25. The vertical position of the horizontal axis is the
energy $\varepsilon^{LN}(0)$ of the line node.

####  XI.1.1 Bulk-band boundaries and projected line node

First, let us discuss the bulk-band boundaries

$E_{y<0,+_{b}}^{LN}(p_{x},p_{z})=\min_{p_{y},\sigma}\varepsilon_{+_{b},\sigma}^{LN}(p_{x},p_{y},p_{z}),$
(136)
$E_{y<0,-_{b}}^{LN}(p_{x},p_{z})=\max_{p_{y},\sigma}\varepsilon_{-_{b},\sigma}^{LN}(p_{x},p_{y},p_{z}),$
(137)

and the projection of the line node onto the surface-momentum plane
$(p_{x},p_{z})$, which are somewhat nontrivial.

Because of the symmetric orientation of the $y=0$ surface relative to the bulk
symmetry axes, the boundaries $E^{LN}_{y<0,\pm_{b}}(p_{x},p_{z})$ of the
conduction and valence bulk bands meet at the curve

$(p_{x},p_{z})=(p_{\perp}^{LN}(\phi)\cos\phi,p_{z}^{LN}(\phi)),\mbox{
}\phi\in[0,\pi],$ (138)

in the surface momentum plane ($\phi$ becomes merely a parameter here), which
is the projection of the 3D line (45), and at energy
$\epsilon=\varepsilon^{LN}(\phi)$:

$E_{y<0,+_{b}}^{LN}(p_{\perp}^{LN}(\phi)\cos\phi,\pm_{u}p_{z}^{LN}(\phi))=E_{y<0,-_{b}}^{LN}(p_{\perp}^{LN}(\phi)\cos\phi,\pm_{u}p_{z}^{LN}(\phi))=\varepsilon^{LN}(\phi),\mbox{
}\phi\in[0,\pi]$ (139)

The two line nodes $\pm_{u}$ have nonoverlapping projections, but for each of
them, the two points at $\pm\phi$, $\phi\in(0,\pi)$, of the 3D line node are
projected onto one point.

As explained in Sec. III.3, the line node has an approximately circular shape,
but does have a weak dependence of the momenta $p_{\perp}^{LN}(\phi)$ and
$p_{z}^{LN}(\phi)$ and energy $\varepsilon^{LN}(\phi)$ on the polar angle
$\phi$. Because of the dependence $\varepsilon^{LN}(\phi)$, there is no single
energy at which the region in the surface-momentum plane $(p_{x},p_{z})$
occupied by the bulk states would be equal to the projection (138) of the line
node. Rather, at any energy, the region occupied by the bulk states is solid.
In Fig. 17(a) and Fig. 18(a), we plot the bulk and surface-state spectrum at
the fixed energy $\epsilon=\varepsilon^{LN}(0)$. Since the dependencies
$\varepsilon^{LN}(\phi)$ and $p_{x}^{LN}(\phi)$ are weak, the region of the
bulk states has only a narrow width along the $p_{z}$ direction. Since
$\varepsilon^{LN}(\phi)\geq\varepsilon^{LN}(0)$, the crossing with this level
occurs for $-_{b}$ bulk states, the equation

$E_{y<0,-_{b}}^{LN}(p_{x},p_{z})=\varepsilon^{LN}(0).$

defines the boundary of that region. Because the $p_{z}^{LN}(\phi)$ and
$\varepsilon^{LN}(\phi)$ dependencies are numerically really weak, these
differences cannot be resolved in the scale of Fig. 17(a) and Fig. 18(a) and
both the projected line node and the fixed-energy region appear as a straight-
line segment.

####  XI.1.2 Surface-state spectrum

We describe the surface-state spectrum $\mathcal{E}^{LN}_{y<0}(p_{x},p_{z})$
of the line-node semimetal for the $y<0$ sample by the 2D plots presented in
Figs. 17-19. Fig. 17 has the same structure as Figs. 7-12 for the Dirac
semimetal. Figs. 18 and 19 further characterize the surface-state spectrum in
the vicinity of the projected line node. Fig. 18 shows the spectrum along the
$p_{x}=\text{const}$ paths and Fig. 19 – along some additional paths.

In accord with Sec. VIII.2, as the Dirac semimetal is transformed into the
line-node semimetal upon introducing the linear BIA term, the ${\bf C}_{2v}$
spatial symmetry of the bulk-band boundaries remains, while the spatial
symmetry of the surface-state spectrum is lowered from ${\bf C}_{2v}$ to ${\bf
C}_{2}$.

For $p_{\beta_{1}}\ll p_{u}$, the key nontrivial changes occur in the vicinity
of the former Dirac points, which transform into the line nodes. Away, the
surface-state spectrum is qualitatively similar to that of the Dirac semimetal
for the most part, compare Fig. 17 to Fig. 12, except for the region around
the touching point at $(p_{x},p_{z})=(0,0)$ in the Dirac semimetal for
$\tfrac{1}{2}<|{\bar{\alpha}}_{0}|<1$, due to the proximity of the surface-
state band to the bulk-band boundary of the Dirac semimetal.

As in the Dirac semimetal in the regime $\tfrac{1}{2}<|{\bar{\alpha}}_{0}|<1$,
there is still only one surface-state band in the line-node semimetal at
${\bar{\alpha}}_{0}=0.886$ of HgTe. However, unlike in the Dirac semimetal,
there appear patches of absent surface states. The touching point
$(p_{x},p_{z})=(0,0)$ in the Dirac semimetal grows into a finite-size region
of absent surface states of size $\sim p_{\beta_{1}}$ in the line-node
semimetal. In the vicinity of the $+_{u}$ projected line node, two patches of
absent surface states emerge, both scaling as $\sim p_{\beta_{1}}$, with
merging contours connected to the two ends of the projected line node: a
smaller patch around $(p_{x},p_{z})=(-p_{\perp 0},p_{z0})$ and a larger one
around $(p_{x},p_{z})=(p_{\perp 0},p_{z0})$; we denote $p_{\perp
0}=p_{\perp}^{LN}(0)$ and $p_{z0}=p_{z}^{LN}(0)$ for brevity. The smaller
patch is visible in Fig. 18(a), while in Fig. 17(a) it is not well-resolved.
Other than these regions, the surface-state band merges with the bulk-band
boundaries at the projected line node from both sides $p_{z}\gtrless
p_{z}^{LN}(\phi(p_{x}))$, as clearly seen in Fig. 18 [$\phi(p_{x})$ is the
inverse of $p_{x}=p_{\perp}^{LN}(\phi)\cos\phi$]. The spectrum around the
$-_{u}$ projected line node is the same by symmetry.

The weak dependence of the momentum $p^{LN}_{z}(\phi)$ of the line node is
manifested in the spectrum along the $p_{z}=p_{z0}$ path, shown in Fig. 19(b):
the “gaps” in the bulk spectrum are only due to $p_{z}^{LN}(\phi(p_{x}))$
depending on $p_{x}$.

In Fig. 18, the spectrum along a series of straight-line paths with
$p_{x}=\text{const}$ is shown. Except for the close proximity to the ends
$p_{x}=\pm p_{\perp 0}$ of the projected line node, the surface-state band
$\mathcal{E}_{y<0}^{LN}(p_{x},p_{z})$ at fixed $p_{x}\in(-p_{\perp
0},+p_{\perp 0})$ is asymptotically linear in $p_{z}$ as it approaches the
projected line node from both sides. Due to the variations of
$p_{z}^{LN}(\phi(p_{x}))$ being much weaker than the $0.1p_{\perp 0}$ scale,
the velocities

$\partial_{p_{z}}\mathcal{E}_{y<0}^{LN}(p_{x},p_{z}\rightarrow
p_{z}^{LN}(\phi(p_{x}))\pm 0)$

can be effectively deduced from the straight-line paths at fixed
$p_{z}=p_{z0}-0.1p_{\perp 0}$ and $p_{z0}+0.1p_{\perp 0}$, shown in Fig. 19(a)
and (c) [the level of the horizontal axis in Fig. 19 is
$\varepsilon^{LN}(0)$]. On the $p_{z}^{LN}(\phi(p_{x}))<p_{z}$ side of larger
$p_{z}$, around the right end of the projected line node, the velocity changes
from negative to positive as $p_{x}$ decreases from $p_{\perp 0}$; the
velocity depends quite linearly on $p_{x}$. The termination point of the Fermi
contour, defined by $\mathcal{E}^{LN}_{y<0}(p_{x},p_{z})=\varepsilon^{LN}(0)$,
on the $p_{z}^{LN}(\phi(p_{x}))<p_{z}$ side [Fig. 18(a)] is also in accord
with the switch of the sign of the velocity
$\partial_{p_{z}}\mathcal{E}_{y<0}^{LN}(p_{x},p_{z})$ with $p_{x}$. On the
$p_{z}<p_{z}^{LN}(\phi(p_{x}))$ side of smaller $p_{z}$, on the other hand,
the velocity stays negative and close to that of the lower bulk-band boundary
in the whole range of $p_{x}$.

In Fig. 19(d),(e),(f), we also plot the spectrum along the circular paths of
radius $p_{r}=0.1p_{\perp 0}$ centered at $(p_{x},p_{z})=(-p_{\perp
0},p_{z0})$, $(0,p_{z0})$, and $(+p_{\perp 0},p_{z0})$. Such circles cross the
projected line node. These circular paths will become particularly insightful
for the Weyl semimetal, where their centers turn into projected Weyl points;
we will discuss the comparison of these spectra between the line-node and Weyl
semimetals in Sec. XIII.1.

###  XI.2 $z>0$ sample

Figure 20:  The surface-state spectrum $\mathcal{E}^{LN}_{z>0}(p_{x},p_{y})$
of the line-node semimetal for the $z>0$ sample. Other conditions same as in
Fig. 17. (a) Surface-momentum plane. (b),(c),(d) The spectrum along the
$\phi=0,+\frac{\pi}{4},-\frac{\pi}{4}$ paths, respectively.

####  XI.2.1 Bulk-band boundaries and projected line node

For the $z>0$ sample, the two line nodes $\pm_{u}$ [Eq. (45)] project onto the
same curve

$(p_{x},p_{y})=(p_{\perp}(\phi)\cos\phi,p_{\perp}(\phi)\sin\phi),\mbox{
}\phi\in[0,2\pi),$

in the surface-momentum plane. The bulk-band boundaries

$E_{z>0,+_{b}}^{LN}(p_{x},p_{y})=\min_{p_{z},\sigma}\varepsilon_{+_{b},\sigma}^{LN}(p_{x},p_{y},p_{z}),$
(140)
$E_{z>0,-_{b}}^{LN}(p_{x},p_{y})=\max_{p_{z},\sigma}\varepsilon_{-_{b},\sigma}^{LN}(p_{x},p_{y},p_{z}).$
(141)

meet at this projected line node at the line-node energy:

$E_{z>0,+_{b}}^{LN}(p_{\perp}^{LN}(\phi)\cos\phi,p_{\perp}^{LN}(\phi)\sin\phi)=E_{z>0,-_{b}}^{LN}(p_{\perp}^{LN}(\phi)\cos\phi,p_{\perp}^{LN}(\phi)\sin\phi)=\varepsilon^{LN}(\phi),\mbox{
}\phi\in[0,2\pi).$ (142)

Because of the energy dependence $\varepsilon^{LN}(\phi)$ of the line node, at
a fixed energy $\varepsilon^{LN}(0)$, there is also a finite-width region of
bulk states, whose contours are defined by

$E_{z>0,-_{b}}^{LN}(p_{x},p_{y})=\varepsilon^{LN}(0).$

Due to the weak numerical dependence of the line-node energy, these deviations
from the circle are not visible in the scale of Fig. 20.

####  XI.2.2 Surface-state spectrum

The surface-state spectrum $\mathcal{E}^{LN}_{z>0}(p_{x},p_{y})$ of the line-
node semimetal for the $z>0$ sample is presented in Fig. 20. Again, compared
to the Dirac semimetal, the key changes occur at the scale $p_{\perp}\sim
p_{\beta_{1}}$ set by the linear BIA term. In the Dirac semimetal, there are
two surface-state bands in the vicinity of the projected double Dirac point,
both connected to the latter. In the line-node semimetal, the two surface-
state bands are degenerate at $p_{\perp}=0$, due to time reversal symmetry
$\mathcal{T}_{-}$. For a given $\phi$, as a function of $p_{\perp}$, the upper
surface-state band passes through the projected line node. The lower surface-
state band merges with the lower bulk-band boundary and then reappears in the
outer region, and we notice that it does so virtually unhybridized with the
bulk bands, as if passing them through. In fact, it visually appears that the
main qualitative effect of linear BIA term is to shift the two surface-state
bands of the Dirac semimetal in opposite directions along $p_{\perp}$; we
demonstrate below that in the linear-in-momentum model this is a precise
mathematical property.

While for every $\phi$ the dependence of the spectrum on $p_{\perp}$ is
qualitatively the same, there is some quantitative dependence on $\phi$ as
well, i.e., axial rotation symmetry is broken, which is predominantly due to
the linear BIA term rather than cubic anisotropy. We plot the spectrum for
fixed $\phi=0,\pm\frac{\pi}{4}$ as a function of $p_{\perp}$ (the spectrum for
$\phi=\frac{\pi}{2}$ is the same as for $\phi=0$). The spectrum is most
stretched along the $\phi=-\frac{\pi}{4}$ direction and most squeezed along
the $\phi=\frac{\pi}{4}$ direction, as also clearly evident in the in the
Fermi contour (red) of the lower surface-state band in the outer region.

####  XI.2.3 Analytical calculation for the linear-in-momentum model

We also find the surface-state spectrum analytically for the linear-in-
momentum model (39), which includes the linear BIA term to leading order. This
provides some additional insights. We have seen that in the basis (36) the
Hamiltonian of the model is decoupled [Eq. (40)] into two $2\times 2$ blocks,
labeled $\pm_{m}$. The BCs (116) for the $z=0$ surface have been derived in
Sec. X.1. Clearly, this form of BCs holds upon any same change of basis within
the $\pm_{u}$ parts of the wave function. Repeating the steps of Sec. X.2.2,
we obtain that the double-degenerate flat-band solution

$\mathcal{E}^{D,\text{lin}}_{z>0,\pm_{m}}(k_{\perp})\equiv\varepsilon_{u},$
(143)

remains also in the presence of the leading BIA term. Indeed, its effect
solely amounts to shifting the spectra of the two decoupled blocks in momentum
$k_{\perp}$ by $\pm m/v_{\perp}$, compare Eqs. (41) and (42) to Eq. (34).
Chiral symmetry still holds in the presence of the BIA term and still provides
an explanation for this behavior. The surface-state spectrum is presented in
Fig. 15(b). The bulk-band boundaries [Eqs. (41) and (42)] are

$E^{LN,\text{lin}}_{z>0,\pm_{b}}(k_{\perp})=\varepsilon_{u}\pm_{b}\sqrt{1-\tfrac{v_{0}^{2}}{v_{z}^{2}}}|v_{\perp}k_{\perp}-m|.$

Another conclusion that can be drawn is that all the deviations (lifting of
the degeneracy and nonzero dispersion) of the surface-state bands (Fig. 20)
from this double-degenerate flat-band behavior are due to the next-order terms
of the Luttinger model that have been neglected in the linear-in-momentum
model.

##  XII Bulk topology and bulk-boundary correspondence of a Weyl semimetal

In Sec. XIII, we calculate and analyze the surface states of the Weyl
semimetal. One of our main goals in that section is exploring the
manifestation of bulk topology in the behavior of the surface states; such
relation is known as bulk-boundary correspondence [58]. In this section, we
remind the basic topological properties of Weyl semimetals and formulate a
type of bulk-boundary correspondence that directly relates the characteristics
of the surface-state spectrum to the Weyl points, which we will actively use
in Sec. XIII.

###  XII.1 Bulk topology of a Weyl semimetal

The topological properties of 3D Weyl semimetals are directly tied to those of
a 2D quantum anomalous Hall (QAH) systems (Chern insulators) [9]. In the
general topological classification scheme [59, 58], both systems belong to
class A with no assumed symmetries in their respective dimensions. Any surface
$\Sigma$ in the 3D momentum space that is closed or extends to infinity may be
regarded as an effective 2D QAH bulk system. The Chern number

$C(\Sigma)=\frac{1}{4\pi}\int_{\Sigma}{\text{d}}^{2}{\mbox{\boldmath{$\sigma$}}}_{\bf
p}\,{\bf B}({\bf p})$ (144)

of this system is given by the flux of the wave-function-basis-independent
Berry curvature vector field ${\bf B}({\bf p})$ through that surface
(${\mbox{\boldmath{$\sigma$}}}_{\bf p}$ is the normal area vector at point
${\bf p}\in\Sigma$ on the surface).

Weyl points are singularities of the Berry curvature, at which the orientation
of the field depends on the direction from which the Weyl point is approached.
The Chern number

$C^{W}({\bf p}_{0})=C(\Sigma_{{\bf p}_{0}})$ (145)

of a Weyl point at ${\bf p}_{0}$ is defined as the Chern number over a surface
$\Sigma_{{\bf p}_{0}}$ enclosing only that point, e.g., a sphere of a
sufficiently small radius centered at the Weyl point. Weyl points can be seen
as sources and drains of the Berry curvature field, depending on the signs of
their Chern numbers. Note that, as a flux integral, the sign of the Chern
number depends on the assigned orientation of the surface. A unified
convention must be chosen for the orientations of the surfaces.

Since the Berry curvature ${\bf B}({\bf p})$ is a rotorless vector field,
$[\nabla_{\bf p}\times{\bf B}({\bf p})]={\bf 0}$, the Chern numbers of
different surfaces obey conservations laws, which follow from the Gauss
theorem. The difference of the Chern numbers of the two surfaces $\Sigma_{1}$
and $\Sigma_{2}$ not crossing any Weyl poins is given by the sum of the Chern
numbers of all Weyl points enclosed in the region between the two surfaces,

$C(\Sigma_{1})-C(\Sigma_{2})=\sum_{{\bf p}_{\alpha}\text{between
{\hbox{\Sigma_1}} and {\hbox{\Sigma_2}}}}C^{W}({\bf p}_{\alpha}).$ (146)

In particular, the Chern number of all surfaces that can be continuously
deformed to one another without crossing any of the Weyl points is the same.
Sets of such surfaces can be regarded as topological phases of effective 2D
QAH systems. The formula (146) can be interpreted as the change as the Chern
numbers across a topological phase transition realized as the surface $\Sigma$
is continuously deformed from $\Sigma_{1}$ and $\Sigma_{2}$. A surface
$\Sigma$ passing through at least one Weyl point can be seen as a topological
phase transition, where the effective 2D QAH system becomes a gapless nodal
semimetal.

For the considered generalized Luttinger model, for the Dirac semimetal, the
Chern number of any surface not crossing the Dirac points is zero by
cancellation, because a Dirac node consists of two Weyl nodes with opposite
Chern numbers. Upon introducing linear and cubic BIA terms, by continuity, the
Chern number remains zero for any surface $\Sigma$ that is not crossed by the
created line node or Weyl points during this process. Knowing this and the
Chern numbers of the Weyl points, one can then easily determine the Chern
numbers of any other surface using Eq. (146). The Chern numbers $C^{W}({\bf
p}_{\alpha})$ of the Weyl points can be calculated using the expansion of the
Hamiltonian about the Weyl points. We choose the sign convention so that the
Chern numbers of the four Weyl points of the two groups ${\bf p}=(\pm
p_{\perp}^{W},0,\pm p_{z}^{W})$ and ${\bf p}=(0,\pm p_{\perp}^{W},\pm
p_{z}^{W})$ equal $-1$ and $+1$, respectively. This corresponds to the
convention as in Ref. 54, with the outward sphere orientation and the flux is
of the local pseudospin of the lower (filled) bulk band.

Among possible surfaces $\Sigma$ in 3D momentum space, it is common to
consider families of parallel planes with the same common orientation and
different positions. For the considered Weyl semimetal, for the
$p_{x}=\text{const}$ family of planes $\Sigma_{p_{x}}$, there are four
“topological phases” separated by the planes
$p_{x}=-p^{W}_{\perp},0,+p^{W}_{\perp}$ of three “phase transitions”. The sum
of the Chern numbers of the two Weyl points in the $p_{x}=-p^{W}_{\perp}$
plane equals $-2$, that of the four Weyl points in the plane $p_{x}=0$ equals
$+4$, and that of the two WPs in the $p_{x}=+p^{W}_{\perp}$ plane equals $-2$.
According to Eq. (146), these determine the changes of the Chern numbers of
the planes as the phase-transition plane is crossed. The Chern numbers of the
$p_{x}=\text{const}$ planes in the phases $p_{x}<-p^{W}_{\perp}$,
$-p^{W}_{\perp}<p_{x}<0$, $0<p_{x}<p^{W}_{\perp}$, and $p^{W}_{\perp}<p_{x}$
therefore equal $C(\Sigma_{p_{x}})=0,-2,+2,0$, respectively, for the normal of
the plane pointing in the positive $p_{x}$ direction. For the
$p_{z}=\text{const}$ planes, since the sum of the Chern numbers of the four
Weyl points in the $p_{z}=\pm p_{z}^{W}$ planes is zero, the Chern number of
any $p_{z}=\text{const}$ plane is zero.

These bulk topological relations are also manifested in the momentum planes of
specific sample surfaces, as the surfaces in 3D momentum space project as
paths onto the 2D surface-momentum plane, as shown in Fig. 21(a) for $y<0$
sample and Fig. 26(a) for $z>0$ sample and will be discussed in more detail in
the respective sections.

###  XII.2 Generalized bulk-boundary correspondence of a Weyl semimetal

####  XII.2.1 Bulk-boundary correspondence for paths in the surface-momentum
plane

In Ref. 54, the formulation of the bulk-boundary correspondence for Weyl
semimetals was presented that is more general than commonly considered. Also,
as its special case, a version of the bulk-boundary correspondence follows
that fully characterizes the vicinity of the Weyl point: the chirality of the
surface-state spectrum along a path enclosing the projected Weyl point is
equal to the Chern number of the Weyl point. In Ref. 54, it was presented for
one Weyl point. The generalization to any number of Weyl points is
straightforward, as we now present.

To be specific, we consider the $z>0$ sample; the formulation is, of course,
applicable to any surface orientation. Consider an arbitrary path ${\bf
p}_{\perp}^{\gamma}(t)$, parameterized by $t$, in the surface-momentum plane
${\bf p}_{\perp}=(p_{x},p_{y})$ that is either closed or terminates with both
ends at infinity $|{\bf p}_{\perp}|=+\infty$, and does not cross any projected
Weyl points. Each such path defines an effective 2D QAH system with the
Hamiltonian $\hat{H}^{W}({\bf p}_{\perp}^{\gamma}(t),\hat{p}_{z})$ on a
generalized cylinder $\Sigma^{\gamma}$ in the 3D momentum space, obtained by
passing straight $p_{z}$ lines (perpendicular to the surface) through this
path $\gamma$. The Chern number $C(\Sigma^{\gamma})$ [Eq. (144)] of such
system is given by the flux through this generalized cylinder.

The surface-state spectrum $\mathcal{E}^{W}_{z>0}({\bf
p}_{\perp}^{\gamma}(t))$ along this path becomes the “edge-state” spectrum as
a function of $t$ of this effective 2D QAH system defined on the space
$(p_{z},t)$. The bulk-boundary correspondence for this system states that the
signed number

$N(\Sigma^{\gamma})=N[\mathcal{E}^{W}_{z>0}({\bf p}_{\perp}^{\gamma}(t))]$
(147)

of the effective chiral edge-state bands is given by the Chern number
$C(\Sigma^{\gamma})$ of the cylinder $\Sigma^{\gamma}$:

$N(\Sigma^{\gamma})=C(\Sigma^{\gamma}),$ (148)

An edge-state band is called chiral if it connects lower $-_{b}$ (“valence”)
and upper $+_{b}$ (“conduction”) bulk bands. The signed number $\pm 1$ of one
chiral band, also called “chirality”, is determined by the order of the
merging points $t_{\pm_{b}}$ with these bulk bands on the $t$ axis. If there
are several chiral bands, then the total signed number (chirality)
$N(\Sigma^{\gamma})$ is given by the sum of their individual chiralities.

The bulk-boundary correspondence of Weyl semimetals is commonly expressed in
terms of “conventional” 2D QAH systems defined on infinite 2D planes in 3D
momentum space. This now becomes a special case of the above more general
formulation, when the paths ${\bf p}_{\perp}^{\gamma}(t)$ are straight lines,
and the generalized cylinders on which the QAH systems are defined are such
planes. For example, for a plane with the fixed $p_{x}=P_{x}=\text{const}$, to
be denoted as $\Sigma_{P_{x}}$, the parameter $t=p_{y}$ is the momentum
component itself. The surface-state spectrum
$\mathcal{E}^{W}_{z>0}(P_{x},p_{y})$ is interpreted in this picture as
effective edge-state spectra (functions of $p_{y}$) of a family a 2D QAH
systems defined on the planes $\Sigma_{P_{x}}$. According to the bulk-boundary
correspondence [Eqs. (152) and (148)], the signed number

$N(\Sigma_{P_{x}})=N[\mathcal{E}^{W}_{z>0}(P_{x},p_{y})]=C(\Sigma_{P_{x}})$
(149)

of the chiral edge-state bands is determined by the Chern number
$C(\Sigma_{P_{x}})$ of the plane $\Sigma_{P_{x}}$. The Chern numbers
$C^{W}({\bf p}_{\alpha})$ of the Weyl points determine via (146) the
corresponding changes in the chiralities of the effective edge-state spectra
upon crossing the Weyl points by the plane:

$N(\Sigma_{P_{x1}})-N(\Sigma_{P_{x2}})=\sum_{{\bf p}_{\alpha}\text{ between
}\Sigma_{P_{x1}}\text{ and }\Sigma_{P_{x2}}}C^{W}({\bf p}_{\alpha}).$ (150)

However, the above more general formulation (148) allows us to formulate, as
its special case, a different type of bulk-boundary correspondence of Weyl
semimetals that is much better suited for the characterization of the surface-
state spectrum in the vicinity of projected Weyl points. Consider a closed
path $\gamma$ that encloses some number of projected Weyl points ${\bf
p}_{\alpha}$. The Chern number

$C(\Sigma^{\gamma})=\sum_{{\bf p}_{\alpha}\text{ inside
}\Sigma^{\gamma}}C^{W}({\bf p}_{\alpha})$ (151)

of the 2D QAH system defined on the cylinder $\Sigma^{\gamma}$ generated by
such path is equal to the net sum of the Chern numbers [Eq. (145)] of the Weyl
points enclosed by the cylinder $\Sigma^{\gamma}$, since such cylinder and the
set of spheres around the Weyl points can be continuously deformed into each
other. According to the bulk-boundary correspondence (148), the signed number
(chirality)

$N(\Sigma^{\gamma})=N[\mathcal{E}^{W}_{z>0}({\bf
p}_{\perp}^{\gamma}(t))]=\sum_{{\bf p}_{\alpha}\text{ inside
}\Sigma^{\gamma}}C^{W}({\bf p}_{\alpha})$ (152)

of the effective chiral edge-state bands $\mathcal{E}^{W}_{z>0}({\bf
p}_{\perp}^{\gamma}(t))$ along the path $\gamma$ is given directly by this net
sum of the Chern numbers of the enclosed Weyl points.

We will explicitly demonstrate the manifestation of both of these types of
bulk-boundary correspondence for the studied Weyl-semimetal model.

####  XII.2.2 Bulk-boundary correspondence for the oriented Fermi contours

In fact, for a given surface, the essential information about the topological
properties of the surface-state spectrum $\mathcal{E}^{W}_{z>0}({\bf
p}_{\perp})$ and its bulk-boundary correspondence for any path ${\bf
p}_{\perp}^{\gamma}(t)$ can be captured in one 2D plot of the oriented fixed-
energy contours in the surface-momentum plane ${\bf p}_{\perp}$. Such fixed-
energy contours (also called Fermi contours, assuming the Fermi level is at
the Weyl-point energy) where the surface-state bands cross the energy
$\varepsilon^{W}$ level of the Weyl points, are defined by the equation

$\mathcal{E}^{W}_{z>0}({\bf p}_{\perp})=\varepsilon^{W}.$

The chirality $N[\mathcal{E}^{W}_{z>0}({\bf p}^{\gamma}_{\perp}(t))]$ of the
effective edge-state spectrum $\mathcal{E}_{z>0}^{W}({\bf
p}^{\gamma}_{\perp}(t))$ along an arbitrary path ${\bf p}_{\perp}^{\gamma}(t)$
can alternatively be determined as the net signed number of its crossings of a
certain energy level, such as $\varepsilon^{W}$, where the sign of the
crossing is determined by the sign of the “velocity”

$\partial_{t}\mathcal{E}^{W}_{z>0}({\bf p}^{\gamma}_{\perp}(t))={\bf
v}^{W}_{z>0}({\bf p}^{\gamma}_{\perp}(t))\cdot\partial_{t}{\bf
p}_{\perp}^{\gamma}(t)$

with respect to $t$. This sign is therefore determined by the scalar product
of the path’s tangent vector $\partial_{t}{\bf p}_{\perp}^{\gamma}(t)$ and the
velocity vector

${\bf v}^{W}_{z>0}({\bf p}_{\perp})=\partial_{{\bf
p}_{\perp}}\mathcal{E}^{W}_{z>0}({\bf p}_{\perp})$

at the intersection of the path and the Fermi contour. The velocity ${\bf
v}^{W}_{z>0}({\bf p}_{\perp})$ is orthogonal to the Fermi contour and which of
the two options of its orientation is realized can be specified without
plotting the vector by assigning an orientation to the Fermi contour. It makes
sense to coordinate the contour orientation with the Chern numbers of the
projected Weyl points so that the points with positive and negative Chern
numbers acts as effective sources and drains, respectively.

Summarizing, for any path $\gamma$, one can determine the chirality
$N(\Sigma^{\gamma})$ of the effective edge state-spectrum along it by
considering its intersection points with the Fermi contour and counting the
numbers of these crossings with the signs of the velocities deduced from the
local orientation of the contour.

Generally, there may be several branches of the Fermi contours connected that
meet at one projected Weyl point, in which case there is a conservation law
according to the “source and drain” picture, that the net number of the
oriented branches has to equal the net Chern number of the projected Weyl
point. In particular, each single projected Weyl point with the Chern number
$\pm 1$ emanates or terminates one Fermi contour. A single Fermi contour
connecting two projected single Weyl points with the Chern numbers $\pm 1$ is
typically referred to as a Fermi arc and is the simplest situation. For $0$ or
$\pm 2$ net Chern numbers of the projected Weyl points, the Fermi-contour
structures are more complicated, as we will see below.

##  XIII Surface states of the Weyl semimetal

In this section, we calculate and explore the surface states of the Weyl
semimetal that arises in the Luttinger semimetal upon introducing compressive
strain and linear and cubic BIA terms. Without the cubic BIA terms, the system
is a line-node semimetal, whose surface states were studied in the previous
Sec. XI. We calculate the surface states within the Luttinger model, described
by the Hamiltonian $\hat{H}^{W}({\bf p})$ [Eq. (49)] and BCs (86), for the
half-infinite $y<0$ and $z>0$ samples, using the semi-analytical method
presented in Sec. VI. We present the Weyl-semimetal case for the parameters
$\alpha_{0,z,\square}$ and $\beta_{1}$ of HgTe provided in Tab. 1 and
compressive strain $u=-3\text{meV}$, the same values as for the line-node
semimetal case. The choice of the cubic BIA parameters has been explained in
Sec. IV.3. Such values obey well the hierarchy $p_{\beta_{3}}\ll
p_{\beta_{1}}\ll p_{u}$ of scales described in Sec. V.

###  XIII.1 $y<0$ sample

Figure 21:  The surface-state spectrum $\mathcal{E}^{W}_{y<0}(p_{x},p_{z})$ of
the Weyl semimetal for the $y<0$ sample, calculated within the Luttinger model
with compressive strain and linear and cubic BIA terms [Eqs. (49) and (86)]
for the parameters of HgTe (Tab. 1) and the strain value $u=-3~{}\text{meV}$.
The structure of the figure is the same as in Fig. 12 for the Dirac semimetal
and Fig. 17 for the line-node semimetal. The projected Weyl points are shown
in green. Away from the (former) projected line nodes, the bulk and surface-
state spectra are very similar to that of the line-node semimetal (Fig. 17).
Figure 22:  The surface-state spectrum $\mathcal{E}^{W}_{y<0}(p_{x},p_{z})$ of
the Weyl semimetal along the $p_{x}=\text{const}$ straight-line paths in the
vicinity of the former $+_{u}$ projected line node $p_{z}\approx
p_{z}^{W}\approx p_{z}^{LN}(0)$ (additional characterization of Fig. 21). The
main difference from the analogous Fig. 18 for the line-node semimetal is that
the line node gaps out everywhere expect for the Weyl points at $p_{x}=0,\pm
p_{\perp}^{W}$. and the surface-state spectrum adjusts accordingly. The edge-
state spectra $\mathcal{E}^{W}_{y<0}(p_{x},p_{z})$ of the effective 2D QAH
systems defined on $p_{x}=\text{const}$ planes obey the bulk-boundary
correspondence (149) in all four “gapped” regions of $p_{x}$. Note that the
aspect ratio in (a) differs from 1 to resolve the features along the $p_{z}$
direction.  Figure 23:  The surface-state spectrum
$\mathcal{E}^{W}_{y<0}(p_{x},p_{z})$ of the Weyl semimetal along the circular
paths around the projected Weyl point
$(p_{x},p_{z})=(-p_{\perp}^{W},p_{z}^{W})$ with the Chern number $C^{W}=-1$
(additional characterization of Fig. 21). The bulk-boundary correspondence
[54] (152) characterizing the vicinity of the projected Weyl point is
satisfied and, for small enough radius, the emergent linear regime (153) is
observed, as identified by the asymptotically straight Fermi and merging
contours as the projected Weyl point is approached; the same concerns Figs. 24
and 25.  Figure 24:  The surface-state spectrum
$\mathcal{E}^{W}_{y<0}(p_{x},p_{z})$ of the Weyl semimetal along the circular
paths around the projected double Weyl point $(p_{x},p_{z})=(0,p_{z}^{W})$
with the net Chern number $C^{W}=+2$ (additional characterization of Fig. 21).
Note that the aspect ratio in (e) differs from 1 in order to resolve the
features along the $p_{z}$ direction.  Figure 25:  The surface-state spectrum
$\mathcal{E}^{W}_{y<0}(p_{x},p_{z})$ of the Weyl semimetal along the circular
paths around the projected Weyl point
$(p_{x},p_{z})=(p_{\perp}^{W},p_{z}^{W})$ with the Chern number $C^{W}=-1$
(additional characterization of Fig. 21).

We describe the surface-state spectrum $\mathcal{E}^{W}_{y<0}(p_{x},p_{z})$ of
the Weyl semimetal for the $y<0$ sample by the 2D plots presented in Figs.
21-25. Fig. 21 shows the overall behavior and has the same structure as Fig.
12 for the Dirac semimetal and Fig. 17 for the line-node semimetal. Figs.
22-25 fully characterize the vicinity of the former projected line node: Fig.
22 shows the spectrum along the series of $p_{x}=\text{const}$ paths,
analogous to Fig. 18 for the line-node semimetal, and Figs. 23, 24, and 25 –
along the circular paths around the projected Weyl points
$(p_{x},p_{z})=(-p_{\perp}^{W},p_{z}^{W})$, $(0,p_{z}^{W})$, and
$(+p_{\perp}^{W},p_{z}^{W})$ of the $+_{u}$ group, respectively. The surface-
state spectrum along some other characteristic paths, not presented here, can
be deduced using the symmetry
$\mathcal{E}_{y<0}^{W}(-p_{x},-p_{z})=\mathcal{E}_{y<0}^{W}(p_{x},p_{z})$.

####  XIII.1.1 Projected bulk spectrum

As explained in Sec. III.4, the cubic BIA terms transform the line-node
semimetal into a Weyl semimetal by “gapping out” each line node in the bulk
spectrum everywhere except for the four Weyl points. For the $y=0$ surface,
there is a total of six projected Weyl points in the surface-momentum plane
$(p_{x},p_{z})$ [Fig. 21(a) and Fig. 22(a)]: two Weyl points with the same
Chern numbers are projected onto each of the points $(p_{x},p_{z})=(0,\pm
p^{W}_{z})$, resulting in their net Chern number $C^{W}=+2$; one Weyl point is
projected onto each of the points $(p_{x},p_{z})=(\pm p^{W}_{\perp},\pm
p^{W}_{z})$, with the Chern number $C^{W}=-1$. The upper and lower bulk-band
boundaries $E_{y<0,\pm_{b}}^{W}(p_{x},p_{z})$ of the Weyl semimetal touch only
at these six projection Weyl points, at the common energy $\varepsilon^{W}$.
For the used parameters, the differences between the values
$(p_{\perp}^{LN}(0),p_{z}^{LN}(0),\varepsilon^{LN}(0))$ of the line node and
$(p_{\perp}^{W},p_{z}^{W},\varepsilon^{W})$ of the Weyl points are negligible.

The “gapping out” of the line node is best seen in the spectrum along the
straight-line path $p_{z}=p_{z}^{W}$ passing through the $+_{u}$ group of the
three projected Weyl points, shown in Fig. 21(j), where the “gap”
$E^{W}_{y<0,+_{b}}(p_{x},p_{z}^{W})-E^{W}_{y<0,-_{b}}(p_{x},p_{z}^{W})\sim\varepsilon_{\beta_{3}}$
[Eq. (53)] in the bulk spectrum is now due to the cubic BIA terms. Compare
this to a much smaller scale of the “gap”
$E^{LN}_{y<0,+_{b}}(p_{x},p_{z0})-E^{LN}_{y<0,-_{b}}(p_{x},p_{z0})$ in the
bulk spectrum of the line-node semimetal along the
$p_{z}=p^{LN}_{z}(0)=p_{z0}$ path, shown in Fig. 19(b), which is only due to
the very weak $p_{z}^{LN}(\phi)$ dependence.

####  XIII.1.2 Surface-state spectrum, general properties

In the surface-momentum plane $(p_{x},p_{z})$, changes occur in the regions of
width $p_{\beta_{3}}$ [Eq. (54)] around the former projected line nodes, which
for the considered values is much smaller than the typical radius $\sim
p_{\beta_{1}}$ [Eq. (43)] of the line-node. Away, the bulk and surface-state
spectrum is very similar to that of the line-node semimetal, compare Fig. 21
to Fig. 17 and Fig. 19. Already for the $p_{z}=p_{z}^{W}\pm 0.1p_{\perp}^{W}$
paths, positioned close to the $p_{z}=p_{z}^{W}$ line, on which the $+_{u}$
group of projected Weyl points resides, the bulk-state and surface-state
spectra of the Weyl semimetal are essentially identical quantitatively to
those of the line-node semimetal: compare Fig. 21(i) and (k) to Fig. 19(a) and
(c).

The key changes upon introducing the cubic BIA terms are as follows. As
demonstrated in Sec. XI, the surface-state band
$\mathcal{E}^{LN}_{y<0}(p_{x},p_{z})$ of the line-node semimetal merges with
the bulk states at both sides of the projected line node. Upon introducing the
cubic BIA terms and opening the small gap $\sim\varepsilon_{\beta_{3}}$ [Eq.
(84)] in the bulk spectrum at the $+_{u}$ projected line node, the two sides
$p_{z}\gtrless p_{z}^{LN}(\phi(p_{x}))$ of the surface-state band
$\mathcal{E}^{LN}_{y<0}(p_{x},p_{z})$ adjust, so that in the Weyl semimetal
the surface-state band $\mathcal{E}^{W}_{y<0}(p_{x},p_{z})$ remains connected
to either the upper or lower bulk-band boundary
$E^{W}_{y<0,\pm_{b}}(p_{x},p_{z})$. This happens differently in the
$p_{x}\in(-p_{\perp}^{W},0)$ and $p_{x}\in(0,+p_{\perp}^{W})$ regions. Namely,
in the region $p_{x}\in(-p_{\perp}^{W},0)$, the part of the surface-state band
approaching the former projected line node from smaller $p_{z}$ and smaller
energies connects to the lower $-_{b}$ bulk-band boundary and the part thereof
approaching from larger $p_{z}$ and higher energies connects to the upper
$+_{b}$ bulk-band boundary. As a result, the surface-state band does not cross
the opened bulk gap in the region $p_{x}\in(-p_{\perp}^{W},0)$ around
$p_{z}=p_{z}^{W}$. In the region $p_{x}\in(0,+p_{\perp}^{W},0)$, the behavior
is reversed. As a result, the surface-state band crosses the opened bulk gap
twice in the region $p_{x}\in(0,+p_{\perp}^{W})$ around $p_{z}=p_{z}^{W}$.

These differences in the bulk and surface-state spectra of the Weyl and line-
node semimetals are best seen in the spectrum along the series of
$p_{x}=\text{const}$ paths, presented in Fig. 22 for the vicinity of the
$+_{u}$ group of projected Weyl points, as compared to Fig. 18 for the line-
node semimetal. The spectrum in the vicinity of the $-_{u}$ group is related
to $+_{u}$ via the symmetry
$\mathcal{E}_{z>0}^{W}(-p_{x},-p_{z})=\mathcal{E}_{z>0}^{W}(p_{x},p_{z})$.

These differences are also seen in the spectrum along the circular paths of
radius $p_{r}=0.1p_{\perp}^{W}$ around the projected Weyl points
$(p_{x},p_{z})=(-p_{\perp}^{W},p_{z}^{W})$, $(0,p_{z}^{W})$, and
$(+p_{\perp}^{W},p_{z}^{W})$, presented in Fig. 23(a), Fig. 24(b), Fig. 25(a)
for the Weyl semimetal, respectively, as compared to Fig. 19(d),(e),(f) for
the line-node semimetal, where the spectrum along the circular paths of radius
$p_{r}=0.1p_{\perp 0}$ centered at the points $(p_{x},p_{z})=(-p_{\perp
0},p_{z0})$, $(0,p_{z0})$, and $(+p_{\perp 0},p_{z0})$ of the projected line
node are presented.

These differences are also manifested in the Fermi
[$\mathcal{E}^{W}_{y<0}(p_{x},p_{z})=\varepsilon^{W}$, red] and merging
[$\mathcal{E}^{W}_{y<0}(p_{x},p_{z})=E^{W}_{y<0,\pm_{b}}(p_{x},p_{z})$,
purple] contours of the Weyl semimetal, presented in Fig. 22(a) and Figs. 23,
24, and 25. We see that the regions of absent surface states around the
$(p_{x},p_{z})=(\pm p_{\perp}^{W},p_{z}^{W})$ projected Weyl points are still
present, inherited from the line-node semimetal, controlled mainly by the
linear BIA term. However, additional Fermi and merging contours appear in the
region of the former projected line node, according to the restructuring of
the bulk and surface-state bands explained above.

Away from the regions of the former projected line nodes, the Fermi contours
of the Weyl semimetal have the same shape as those of the line-node semimetal
[Fig. 21(a)]. As the Fermi contour connected to the
$(p_{x},p_{z})=(-p_{\perp}^{W},p_{z}^{W})$ projected Weyl point deviates from
the $+_{u}$ region, it becomes the same as that of the line-node semimetal
until it eventually reaches the $-_{u}$ region, where it connects to the
$(0,-p_{z}^{W})$ projected double Weyl point. The Fermi contour connected to
the $(p_{\perp}^{W},p_{z}^{W})$ projected Weyl point connects to the
$(0,p_{z}^{W})$ projected double Weyl point along an almost straight line. The
other Fermi contour connected to $(0,p_{z}^{W})$ also moves towards the
$(p_{\perp}^{W},p_{z}^{W})$ point; however, it does not connect to it, but
only passes closely; after which it deviates from the $+_{u}$ region and
follows closely the Fermi contour of the line-node semimetal until it reaches
to $-_{u}$ group, which it connect to the $(p_{\perp}^{W},-p_{z}^{W})$ point.
Altogether, there are two disconnected contours, in accord with absent surface
states crossing the $\varepsilon^{W}$ level in the
$p_{x}\in(-p_{\perp}^{W},0)$ region around $p_{z}=+p_{z}^{W}$ and in the
$p_{x}\in(0,+p_{\perp}^{W})$ region around $p_{z}=-p_{z}^{W}$. Note that,
similar to the discussion for the Dirac semimetal in Sec. IX.6, this only
means that there is no connectivity at this energy: all six projected Weyl
points are connected by the surface-state band.

####  XIII.1.3 Surface-state spectrum, topological properties

We now explore the topological aspect of the Weyl semimetal, as well as of the
Dirac and line-node semimetals. As per the general topological framework
presented in Sec. XII.2.1, the surface-state spectrum
$\mathcal{E}^{W}_{y<0}(p_{x}^{\gamma}(t),p_{z}^{\gamma}(t))$ along various
paths $(p_{x}^{\gamma}(t),p_{z}^{\gamma}(t))$ presented in Figs. 21-25
(straight lines and circles), can be regarded as edge-state spectra of
effective generalized 2D QAH systems, defined on the respective cylinders
$\Sigma^{\gamma}$. All these spectra obey the bulk-boundary correspondence
(148), as we now demonstrate.

According to the changes in the spectrum in the vicinity of the former
projected line nodes, the contribution from the region around
$p_{z}=p_{z}^{W}$ to the chirality of the effective edge-state bands oriented
with growing $p_{z}$ is $0$ in the $p_{x}\in(-p_{\perp}^{W},0)$ region and
$+2$ in the $p_{x}\in(0,+p_{\perp}^{W})$ region. By the symmetry
$\mathcal{E}_{y<0}^{W}(-p_{x},-p_{z})=\mathcal{E}_{y<0}^{W}(p_{x},p_{z})$, the
contribution to the chirality from the region around $p_{z}=-p_{z}^{W}$, is
$-2$ in the $p_{x}\in(-p_{\perp}^{W},0)$ region and $0$ in the
$p_{x}\in(0,+p_{\perp}^{W})$ region. Adding these together, the chiralities
$N(\Sigma_{p_{x}})$ of the effective edge-state spectra of the
$p_{x}=\text{const}$ QAH systems (Fig. 22) are indeed in accord with the Chern
numbers $C(\Sigma_{p_{x}})=0,-2,2,0$, of the four topological phases at
$p_{x}<-p_{\perp}^{W}$, $-p_{\perp}^{W}<p_{x}<0$, $0<p_{x}<p_{\perp}^{W}$, and
$p_{\perp}^{W}<p_{x}$, respectively, manifesting bulk-boundary correspondence
(Sec. XII.2.1). The planes with $p_{x}=-p_{\perp}^{W},0,p_{\perp}^{W}$ that
cross the projected Weyl points represent the phase transitions between these
phases.

The chirality of the effective edge-state spectra of the $p_{z}=p_{z}^{W}\pm
0.1p_{\perp}^{W}$ QAH systems [Fig. 21(i),(k),(f)] is zero, in accord with the
zero Chern number $C(\Sigma_{p_{z}})=0$ of the topological phases
$p_{z}<-p_{z}^{W}$, $-p_{z}^{W}<p_{z}<p_{z}^{W}$, $p_{z}^{W}<p_{z}$. We remind
that for $p_{z}=\text{const}$ QAH systems, although there are phase
transitions at the planes $p_{z}=\pm p^{W}_{z}$ containing the $\pm_{u}$
groups of Weyl points, there is no change of the Chern number, since the sum
of the Chern numbers of the Weyl points at $p_{z}=\pm p^{W}_{z}$ is zero:
$1-1-1+1=0$.

Also, we have assigned orientation according to the source-and-drain picture
to all Fermi contours in the surface-momentum planes in all Figs. 21-25 for
the Weyl semimetal. As explained in Sec. XII.2.2, this allows one to deduce
the chirality of effective edge-state spectrum along any path. We observe that
the so-determined chiralities are indeed in full accord with the bulk-boundary
correspondence (148).

We can now also comment on the topological aspect of the surface-state
spectrum of the Dirac and line-node semimetals. These systems can be regarded
as special cases of Weyl semimetals and their surface states can also be
analyzed in terms of quantum Hall topology and its bulk-boundary
correspondence. They are trivial Weyl semimetals with composite, net-neutral
nodes, since for any surface not crossing the nodes the effective 2D QAH
system is trivial with zero Chern number.

Of particular interest are the circular paths enclosing one projected Dirac
point [(d) and (h) of Figs. 7-12 and Fig. 15]. The Chern number of the Dirac
point is zero, and in all six presented cases of ${\bar{\alpha}}_{0}$, the
surface-state spectrum along such circular paths is in accord with bulk-
boundary correspondence. Although the surface-state spectrum is qualitatively
different in the two regimes $|{\bar{\alpha}}_{0}|<\tfrac{1}{2}$ and
$\tfrac{1}{2}<|{\bar{\alpha}}_{0}|<1$, it is topologically the same and has
the same zero chirality: for $|{\bar{\alpha}}_{0}|<\tfrac{1}{2}$, there are
two edge-state bands, one merging with the lower bulk-band boundary
$E^{D}_{y<0,-_{b}}(p_{x},p_{z})$, and the other with the upper bulk-band
boundary $E^{D}_{y<0,+_{b}}(p_{x},p_{z})$; for
$\tfrac{1}{2}<|{\bar{\alpha}}_{0}|<1$, there is one edge-state band not
connected to the bulk-band boundaries. Of course, the net chirality of the
effective edge-state spectrum is also zero for any other path (like straight
lines) not crossing the projected Dirac points.

This topological property of the surface states is preserved by continuity
when the Dirac semimetal is transformed into the line-node semimetal and when
the latter is in turn transformed into the Weyl semimetal. For the circular
paths enclosing the $+_{u}$ Dirac point of the Dirac semimetal [Fig. 12(d) and
(h)], for the circular paths enclosing the $+_{u}$ line node of the line-node
semimetal [Fig. 17(d) and (h)], and for the circular paths enclosing the whole
$+_{u}$ group of projected Weyl points of the Weyl semimetal [Fig. 21(d) and
(h)], the effective edge-state spectra are topologically equivalent (moreover,
those of the line-node and Weyl semimetals are very close quantitatively),
their chirality is zero in accord with the zero total net Chern number of the
enclosed singularities (nodal points or lines), manifesting bulk-boundary
correspondence. This is because all three contours and systems can be
continuously deformed to each other without crossing the nodes (and therefore
while the effective QAH systems remain gapped).

####  XIII.1.4 Topological properties and emergent linear scaling near
projected Weyl points

Next, we consider circular paths enclosing individual projected Weyl points.
Circular paths of sufficiently small radii centered around the projected Weyl
points are particularly useful, as they provide a direct and complete
characterization of the whole surface-state spectrum in their vicinity,
serving two purposes: capturing the local topological properties and the
emergent linear scaling. For each projected Weyl point
$(p_{x},p_{z})=(-p_{\perp}^{W},p_{z}^{W})$, $(0,p_{z}^{W})$, and
$(+p_{\perp}^{W},p_{z}^{W})$ of the $+_{u}$ group, we present the surface-
state spectra along several such circles with decreasing radii $p_{r}$ in
Figs. 23, 24, and 25, respectively.

First, we explicitly see the bulk-boundary correspondence (152) manifest: the
signed number (chirality) of the effective edge-state bands connecting the
lower and upper bulk-band boundaries is indeed in full accord with the net
Chern number of the projected Weyl points enclosed within the circle. Although
the spectra around the same point differ quantitatively for different radii,
their topological properties remain the same.

Second, as the radius $p_{r}$ is decreased, the emergence of the momentum
scale $p_{\beta_{3}}$ [Eq. (54)] of the cubic BIA terms becomes evident. Below
this scale around the projected Weyl points, a universal linear scaling
behavior of the surface-state spectrum emerges, e.g.,

$\mathcal{E}_{y<0}^{W}(p_{\perp}^{W}+p_{r}\sin\zeta,p_{z}^{W}+p_{r}\cos\zeta)=\varepsilon^{W}+\mathcal{V}_{y<0}^{W}(\zeta)p_{r}+o(p_{r}),$
(153)

for the projected Weyl point $(p_{x},p_{z})=(p_{\perp}^{W},p_{z}^{W})$ and
similar for the other points. In this regime, the spectrum is fully
characterized by its velocity $\mathcal{V}^{W}_{y<0}(\zeta)$ as a function of
the angle $\zeta$. Therefore, for small enough radius $p_{r}$, the spectrum
along the circle essentially maps out the velocity
$\mathcal{V}^{W}_{y<0}(\zeta)$.

In Ref. 54, the general surface-state structure in the vicinity of one
projected Weyl point was calculated for the most general form of the linear-
in-momentum Hamiltonian and BCs. The main general properties are: there is one
surface-state band, which exist in a half-plane of the surface momentum; the
merging half-lines with the upper and lower $\pm_{b}$ bulk-band boundaries
form one straight line. Hence, for any model, the surface-state properties
must asymptotically reduce to this behavior in the vicinity of a single
projected Weyl point. The linear regime can be best identified by the Fermi
and merging contours, which become asymptotically straight as the projected
Weyl point is approached (and also by the linear spectrum along the straight-
line paths passing through the Weyl points). We do observe this asymptotic
behavior for the considered generalized Luttinger model of the Weyl semimetal
in the vicinity of the four single projected Weyl points $(\pm
p_{\perp}^{W},\pm p_{z}^{W})$ [Figs. 23 and 25]. For the
$(-p_{\perp}^{W},p_{z}^{W})$ point, the linear regime is reached at slightly
less than $p_{r}=0.001p_{\perp}^{W}$ [Fig. 23(g)]. For the
$(p_{\perp}^{W},p_{z}^{W})$ point, the linear regime is reached at around
$p_{r}=0.01p_{\perp}^{W}$ [Fig. 25(f)]. For the $(0,p_{z}^{W})$ double
projected Weyl point, the linear regime is reached earlier, already at
$p_{r}=0.1p_{\perp}^{W}$. [Fig. 24(d) and (e)]. For the double projected Weyl
point, the general linear-in-momentum model has not yet been derived and
studied, but, of course, linear scaling is also anticipated in this case.

The emergence of the linear scaling proves that the vicinities of projected
Weyl points could be described by a linear-in-momentum model, which would be
the last low-energy model in the hierarchy of Sec. V. Although we do not
derive these models from the generalized Luttinger model here, the clearly
observed linear scaling proves that they exist and could be derived.

We also point out the following. Although the surface-state spectrum along the
$p_{x}=0$ straight-line path is qualitatively the same for the Dirac [Fig.
12(b)] and Weyl [Fig. 21(b) and Fig. 22(g)] semimetals, the full spectra in
the vicinity of the projected Dirac point $(p_{x},p_{z})=(0,p_{u0})$ and the
projected double Weyl point $(p_{x},p_{z})=(0,p_{z}^{W})$ are topologically
and qualitatively completely different, as the spectra along the circular
paths show: for the Dirac point with the $C^{W}=0$ Chern number, the two
branches at $p_{x}=0$ are cuts of the same topologically trivial band at
$\zeta=0,\pi$ [Fig. 12(d) and Fig. 15]; for the projected Weyl point with the
$C^{W}=+2$ Chern number, the two branches are cuts of two different chiral
subbands [Fig. 24(a),(b),(c)]. Therefore, it is important to realize that
single straight-line paths alone do not provide a complete picture of the
surface-state spectrum and such similarities can be misleading. The plots
along circular paths of sufficiently small radii, on the other hand, fully
characterize the surface-state spectrum in the vicinity of the nodes.

###  XIII.2 $z>0$ sample

Figure 26:  The surface-state spectrum $\mathcal{E}^{W}_{z>0}(p_{x},p_{y})$ of
the Weyl semimetal for the $z>0$ sample. Other conditions same as in Fig. 17.
(a) Surface-momentum plane. (b),(e),(h) The spectrum along the
$\phi=0,+\frac{\pi}{4},-\frac{\pi}{4}$ paths, respectively, with the insets
around the former projected line node. (k),(l),(m),(p),(q),(r) The spectrum
along the circular paths around the projected Weyl point
$(p_{x},p_{y})=(p_{\perp}^{W},0)$ and (n),(o),(s),(t) the surface momentum
plane at the same scale. In (o) and (t), the aspect ratio differs from 1 to
resolve the features along the $p_{z}$ direction.

The surface-state spectrum $\mathcal{E}^{W}_{z>0}(p_{x},p_{y})$ of the Weyl
semimetal for the $z>0$ sample is presented in Fig. 26. As a consequence of
the bulk band structure explained in Sec. XII.1, the Weyl points with equal
Chern numbers are projected onto each other in pairs in the surface-momentum
plane $(p_{x},p_{y})$. This leads to four projected Weyl points
$(p_{x},p_{y})=(\pm p_{\perp}^{W},0)$ and $(0,\pm p_{\perp}^{W})$ [Fig. 26(a)]
with the net Chern numbers $C^{W}=\mp 2$, respectively. Away from the region
of width $p_{\beta_{3}}$ around the former projected line node, the spectrum
is very similar to that of the line-node semimetal, compare to Fig. 20.

The same main changes upon introducing the cubic BIA terms and transforming
the line-node semimetal to the Weyl semimetal occur in the surface-state
spectrum for the $z>0$ sample as for the $y<0$. As demonstrated in Sec. XI,
the upper band of the surface-state spectrum
$\mathcal{E}^{LN}_{z>0}(p_{x},p_{y})$ of the line-node semimetal merges with
the projected line node from both the inner $p_{\perp}<p_{\perp}^{LN}(\phi)$
and outer $p_{\perp}^{LN}(\phi)<p_{\perp}$ regions. Upon introducing the cubic
BIA terms and opening the gap in the bulk spectrum in the four angular sectors
of $\phi$, the two sides sides of the upper surface-state band adjusts, so
that in the Weyl semimetal the upper band of the surface-state spectrum
$\mathcal{E}^{W}_{z>0}(p_{x},p_{y})$ remains connected to either the upper or
the lower bulk-band boundary $E^{W}_{z>0,\pm_{b}}(p_{x},p_{y})$. This happens
differently in different angular sectors between the projected Weyl points. In
the sectors $0<\phi<\frac{\pi}{2}$ and $\pi<\phi<\frac{3\pi}{2}$ [Fig.
26(e),(f), and (g)], the inner band growing with $p_{\perp}$ merges with the
$+_{b}$ upper bulk-band boundary and the growing outer band merges with the
$-_{b}$ lower bulk-band boundary. As a result, in these sectors, these parts
of the surface-state band cross the gaps, each contributing $+1$ chirality for
any path oriented with growing $p_{\perp}$. In the remaining sectors
$\frac{\pi}{2}<\phi<\pi$ and $\frac{3\pi}{2}<\phi<2\pi$ [Fig. 26(h),(i), and
(j)], the situation is reversed: the inner band growing with $p_{\perp}$
merges with the $-_{b}$ lower bulk-band boundary and the growing outer band
merges with the $+_{b}$ lower bulk-band boundary. As a result, in these
sectors, these parts of the surface-state band do not cross the gap and
contribution nothing to chirality. These properties are also manifested in the
Fermi and merging contours [Fig. 26(a),(n),(o),(s),(t)]. From each projected
Weyl point $(p_{x},p_{y})=(0,\pm p_{\perp}^{W})$ with $C^{W}=+2$, two inner
Fermi contours originate and terminate at the projected Weyl point
$(p_{x},p_{y})=(\pm p_{\perp}^{W},0)$ with $C^{W}=-2$, respectively.

The net chirality of the effective edge-state bands can be deduced from this
for any path and will be equal to the Chern number of the effective QAH
system, manifesting the bulk-boundary correspondence. In particular, for
straight-line paths passing through ${\bf p}_{\perp}={\bf 0}$, with
$\phi=\text{const}\neq 0,\pm\frac{\pi}{2},\pi$ the chirality is always zero,
in accord with the zero Chern number, but there are two different
realizations. For $0<\phi<\frac{\pi}{2}$, there is a cancellation of $+2$ and
$-2$ contributions. For $\frac{3\pi}{2}<\phi<2\pi$, there are no edge-state
bands crossing the gap. The zero chirality and Chern number for any
$\phi=\text{const}\neq 0,\pm\frac{\pi}{2},\pi$ path can also be explained in
terms of $\mathcal{T}_{-}$ time-reversal symmetry, see Sec. XIV.1.

In Fig. 26(k),(l),(m),(p),(q),(r), we also present the spectrum along the
circular paths enclosing individual projected Weyl points for the radii
$p_{r}=0.1p_{\perp}^{W},0.01p_{\perp}^{W}$. The bulk-boundary correspondence
is manifested. From the plots Fig. 26 (n),(o),(s),(t) of the surface-momentum
plane, linear scaling, analogous to Eq. (153), is seen emerging at about
$p_{r}=0.001p_{\perp}^{W}$ as the Fermi and merging contours become
asymptotically straight as the projected Weyl point is approached.

##  XIV Others aspects

###  XIV.1 Effective 2D $\mathcal{T}_{-}$ insulators

Just like any surface (in particular, any plane) in the 3D momentum space of a
Weyl semimetal can be regarded as an effective 2D QAH system (class A), for
present $\mathcal{T}_{-}$ symmetry (which is the case for the considered
model), any $\mathcal{T}_{-}$-invariant surface can be regarded as an
effective 2D $\mathcal{T}_{-}$-symmetric topological insulator (if no Weyl
points are crossed) (class AII). Since the BC (86) is
$\mathcal{T}_{-}$-symmetric, the corresponding effective edge-state spectrum
satisfies the class-AII bulk-boundary correspondence.

Any plane passing through ${\bf p}={\bf 0}$ is $\mathcal{T}_{-}$-symmetric.
For the $y<0$ sample, relevant $\mathcal{T}_{-}$-symmetric planes are those
that project onto the lines $(p_{x},p_{z})=p_{r}(\sin\zeta,\cos\zeta)$ with
$\zeta=\text{const}$. There are three phases as a function of $\zeta$. The
effective edge-state spectrum, which for $|\zeta|\ll 1$ can be deduced from
the $p_{x}=\text{const}$ paths around $p_{z}=p_{z}^{W}$ in Fig. 22 and
symmetry, is trivial in terms of class-AII topology in all three phases,
although differs qualitatively. For the $z>0$ sample, relevant
$\mathcal{T}_{-}$-symmetric planes are those that project onto the lines
$(p_{x},p_{y})=p_{\perp}(\cos\phi,\sin\phi)$ with $\phi=\text{const}$. There
are two phases: $\phi\in(-\frac{\pi}{2},0)$ and $\phi\in(0,\frac{\pi}{2})$.
The effective edge-state spectrum is nontrivial in terms of class-AII topology
in both phases [Fig. 26(e) and (h)], although differs qualitatively.

###  XIV.2 Relation between different surfaces

Figure 27:  Relation between the surface-state spectra of the Weyl semimetal
for different surface orientations. For an effective 2D QAH system, specified
by a plane in 3D momentum space of the Weyl semimetal, exemplified here with
$p_{x}=p_{\perp}^{W}/2$ plane of the $0<p_{x}<p_{\perp}^{W}$ topological
phase, the effective edge-state spectra $\mathcal{E}^{W}_{y<0}(p_{x},p_{z})$
for the $y<0$ sample and $\mathcal{E}^{W}_{z>0}(p_{x},p_{y})$ for the $z>0$
sample have the same chirality $N=2$, equal to the Chern number
$C(\Sigma_{p_{x})})=2$ of the plane, in accord with the bulk-boundary
correspondence (149). This can be understood from the visualization of the
chiral quantum Hall edge states following the edge of one sample of any shape,
such as a sample occupying the quadrant region specified by both $y<0$ and
$z>0$ conditions. Other conditions are the same as in Fig. 21.  Figure 28:
Comparison of the surface-state spectra of the Weyl semimetal calculated
within the generalized Luttinger (first three rows) and Kane (last three rows)
models, for the parameters of HgTe (Tabs. 1 and 2) and for different values of
compressive strain $u$. For each model, the first row presents the Fermi
contours of the surface-state band at the Weyl-point energy $\varepsilon^{W}$,
the second and third rows present the spectrum along the
$\phi=\pm\frac{\pi}{4}$ paths, respectively. For small enough strain $u$ (like
$u=-3\text{meV}$), there is qualitative and good quantitative agreement at
momenta on order $p_{\perp}\sim p_{\beta_{1}}$. For larger strain, outside of
the validity range of the Luttinger model, the spectrum of the Kane model
undergoes qualitative transformations, while the spectrum of the Luttinger
model remains qualitatively the same. The units of momentum and energy in
these plots are $\mu\text{m}^{-1}$ and meV.

Although considering effective 2D QAH systems defined on cylinders in the 3D
momentum space of a Weyl semimetal is particularly useful for characterizing
the vicinity of the Weyl points, this construction is tailored to a specific
surface orientation. Considering QAH systems defined on planes, on the other
hand, has the advantage that different surfaces become different effective
edges of the same QAH system, the topological properties of the edge-state
spectra of which are related.

For instance, for $p_{x}=\text{const}$ planes, the surface-state spectra
$\mathcal{E}_{y<0}^{W}(p_{x},p_{z})$ and $\mathcal{E}_{z>0}^{W}(p_{x},p_{y})$
can be regarded as effective edge-state spectra at fixed $p_{x}=\text{const}$
as functions of $p_{z}$ and $p_{y}$, respectively. The bulk-boundary
correspondence [Eq. (149)] has to hold for any edge of a given QAH system;
hence, the chiralities of these edge-state spectra have to be the same and
equal to the bulk Chern number $C(\Sigma_{p_{x}})$. This can be understood
from the chiral properties of the edge states of a QAH system (Fig. 27) with
an edge consisting of straight-line sections with different orientations: due
to their chiral nature, the edge states cannot terminate and persist with
preserved chirality into the next section.

We illustrate this in Fig. 27 for the effective QAH system defined on the
$p_{x}=p_{\perp}^{W}/2$ plane, which exemplifies the $0<p_{x}<p_{\perp}^{W}$
phase with the Chern number $C(\Sigma_{p_{x}})=+2$. The surface states for
$y<0$ and $z>0$ samples in this case have already been presented in Figs. 21
and 26. According to the considerations in Secs. XIII.1 and XIII.2, the net
chirality of the effective edge-state spectra is $N=+2$ for both surfaces, as
seen in Fig. 27. Note that it is important here that both $(-y,x,z)$ and
$(z,x,y)$ form right-handed systems, which is why the chiralities of the
spectra for $y<0$ and $z>0$ samples are the same (and not opposite).

###  XIV.3 Comparison between the Luttinger and Kane models

In this section, we perform the comparison of the surface-state spectra
calculated within the generalized Luttinger and Kane models in the Weyl
semimetal phase, i.e., with all perturbations present: strain and linear and
cubic BIA terms. As explained in Secs. IV and V, the Luttinger model is the
rigorous low-energy limit of the Kane model, when all the relevant energy
scales are much smaller than the level spacing $\Delta$ between the
$j=\frac{3}{2}$ and $j=\frac{1}{2}$ states. For the hierarchy of scales we
consider, the largest scale in the Luttinger model is the strain $|u|$;
therefore, agreement is expected when the condition $|u|\ll|\Delta|$ is
satisfied. We now explore what this means quantitatively by comparing the
surface-state spectra for various values of $|u|$ while $\Delta$ and other
parameters remains fixed and equal to those of HgTe (Tabs. 1 and 2).

We calculate the spectra for the $z>0$ sample along the characteristic
$\phi=\pm\frac{\pi}{4}$ paths and their Fermi contours, presented in Fig. 28.

We notice that for the bulk spectrum at the momentum scale $p_{\beta_{1}}$
[Eq. (43)] of the linear BIA term, there is no appreciable quantitative
difference between the Luttinger and Kane models for the presented values of
strain and there are no qualitative differences as the strain is changed: the
values $p_{\perp}^{W}$ of the Weyl points remain almost the same (they are
determined mainly by the $p_{\beta_{1}}$ scale, which does not depend on
strain), and only the energy position and splitting of the bulk bands have
some dependence on $|u|$.

Further, there are also no qualitative changes in the surface-state spectrum
of the Luttinger model upon changing strain. This is understandable, since
there is no larger scale in the Luttinger model that could compete with
strain.

At the smallest presented value $u=-3\text{meV}$ (which is the value used
throughout for the line-node and Weyl semimetals), the surface-state spectra
of the Kane and Luttinger models agree qualitatively and reasonably well
quantitatively, confirming that at such value the Luttinger model (with the
values of other parameters those of HgTe) is applicable. (Upon decreasing the
strain value, quantitative agreement would only improve.) The agreement is
best seen in the Fermi contours: there is a closed Fermi contour in the outer
region and two “double” Fermi contours connected to the pairs of the projected
double Weyl points.

However, we observe that, upon increasing the strain magnitude $|u|$ and thus
making it more comparable to the level spacing $|\Delta|$, the surface-state
spectrum of the Kane model undergoes not just quantitative, but qualitative
transformations. The main “driver” of these changes is the surface-state
energy $\mathcal{E}^{W,K}_{z>0,\pm}({\bf 0})$ at the double-degeneracy point
${\bf p}_{\perp}={\bf 0}$ moving upwards relative to the bulk-band boundaries
upon increasing $|u|$; this effect does not happen in the Luttinger model.
Upon this evolution, first, the connectivity to the Weyl points switches from
the lower to the upper surface-state band between $u=-3\text{meV}$ and
$u=-10\text{meV}$. The energy of the ${\bf p}_{\perp}={\bf 0}$ point still
remains below the Weyl-point (Fermi) level $\varepsilon^{W,K}$. The closed
Fermi contour shrinks and passes through the (former) line-node region. For
smaller strain magnitude $|u|$, this Fermi contour in the outer region is due
to the lower surface-state band $\mathcal{E}^{W,K}_{z>0,-}({\bf p}_{\perp})$
crossing the Fermi level and for larger $|u|$, this Fermi contour in the inner
region is due to the upper surface-state band $\mathcal{E}^{W,K}_{z>0,+}({\bf
p}_{\perp})$. Further, at the strain value very close to $u=-14\text{meV}$,
the surface-state energy $\mathcal{E}^{W,K}_{z>0,\pm}({\bf 0})$ crosses the
Fermi level $\varepsilon^{W,K}$, where the Fermi contour shrinks into a point.
At larger strain magnitude, $u=-16,-17,-20\text{meV}$, the closed Fermi
contour in the inner region is again due to the crossing of the lower surface-
state band. The qualitative difference in the Fermi contours at
$u=-16\text{meV}$ and $u=-17,-20\text{meV}$ is whether along the
$\phi=\frac{\pi}{4}$ path (and in a range of $\phi$ around it) the lower
surface-state band crosses the Fermi level twice or not at all, respectively,
whether it has a minimum in $p_{\perp}$ below or above the Fermi level.

##  XV Conclusion

We have presented the theoretical analysis of the surface states of the
Luttinger semimetal in the presence of several perturbations: compressive
strain and linear- and cubic-in-momentum bulk-inversion-asymmetry (BIA) terms.
There are several physical, technical, and methodological points presented in
this work, which we summarize and emphasize in this conclusion.

It has been demonstrated that the Luttinger semimetal without perturbations
can indeed be regarded as the parent, highest-symmetry system for multiple
nontrivial phases. For the considered perturbations, the system has a peculiar
structure, where each successive perturbation creates a new semimetal phase
via the modification of the nodal structure of the previous phase, resulting
in a sequence of four semimetal phases: the parent Luttinger semimetal and
Dirac, line-node, and Weyl semimetals.

Calculation of the surface states of these semimetal phases has been the main
practical goal of this work. Since already the unperturbed Luttinger semimetal
exhibits surface states [8], the surface states of the new phases can be
regarded as their evolution.

Among the theoretical technical points, we have shown that the semi-analytical
method of calculating surface states within the continuum models with BCs
allows for calculations for a true half-infinite system with any desired
resolution. This allowed us to capture very fine features of the surface-state
spectrum, e.g., in the vicinity of the nodes. We have also used an efficient
and insightful method of characterizing surface-state spectra through a set of
2D plots, which consists of the plot of the surface-momentum plane with Fermi
and merging contours and shaded regions occupied by the bands and 2D plots of
the surface-state spectrum along a set of paths: key straight-line paths and
circular paths around the projected point nodes. Circular paths provide a full
characterization of the vicinity of the projected point nodes, including the
demonstration of bulk-boundary correspondence. This set of 2D plots provides a
comprehensive description of the whole 3D plot of the surface-state spectrum
as a function of the 2D surface momentum.

Among the interesting physical features of the system, for the considered Weyl
semimetal, we observed multiple situations where the effective edge-states
bands of different QAH systems with the same Chern number (either within one
gapped phase or in different phases with the same Chern number), while having
the same chirality and being topologically equivalent, nonetheless differ
quite a bit qualitatively (such as the situations of absent edge states
crossing the bulk gap and present edge states with cancelled chiralities for
the zero Chern number).

The main theoretical methodological point has been the application of
continuum models with BCs for the study of surface states. In accord with the
hierarchy of the perturbations and semimetal phases they create, there is a
hierarchy of low-energy continuum models describing the vicinities of the
nodal structures. We derive most of these models and demonstrate quantitative
asymptotic agreement between some of them, which explicitly proves that
continuum models with proper BCs are perfectly applicable for the study of
surface states.

###### Acknowledgements.

M. K. acknowledges the financial support by the Grant No. KH 461/1-1 of the
Deutsche Forschungsgemeinschaft (DFG, German Research Foundation). J.-B. M.
and E. M. H. acknowledge funding by the DFG through SFB 1170, Project-ID
258499086, through the Würzburg-Dresden Cluster of Excellence on Complexity
and Topology in Quantum Matter ct.qmat (EXC2147, Project-ID 390858490), as