Surface termination dependent quasiparticle scattering interference and magneto-transport study on ZrSiS

The study on Dirac-like systems has progressed greatly in the past decade. While two-dimensional (2D) Dirac fermions were found in graphene [1] and on the surfaces of topological insulators [2], three-dimensional (3D) Dirac fermions that are protected by inversion symmetry and time reversal symmetry were demonstrated experimentally in Na₃Bi [3] and Cd₃As₂ [4]. Weyl fermions that stem from Dirac fermions due to the broken inversion symmetry or broken time reversal symmetry were also observed e.g. in TaAs [5]. Furthermore, it has been predicted that Dirac nodal line semimetal phase exists when the non-degenerate conduction and valence bands cross each other to form a closed loop [6]. Experimentally, several Dirac nodal line semimetals, within which their linear bands crossings connect to form continuous lines or arcs, were observed in PbTaSe₂ [7], PtSn₄ [8] and ZrSiS family [9–13]. However, it has been revealed that both PbTaSe₂ [7, 14] and PtSn₄ [8] exhibit the complex band structure, making them difficult to study the Dirac nodal bands separately. In contrast, the electronic band structure of ZrSiS is relatively simple. Indeed, recent ARPES studies on ZrSiS revealed the existence of non-symmorphic symmetry protected line nodes [9–13], which was first predicted by Kane and Young in 2015 [15]. In addition, there also exist a Dirac surface band at $\bar{{\rm{X}}}$ and a diamond-shaped Dirac band at $\bar{{\rm{\Gamma }}}$ from ARPES [9–13, 16] and scanning tunnelling microscope (STM) study [17, 18]. Both Dirac bands are also shown to disperse linearly up to few hundred meV even above E_F by quasiparticle scattering interference (QPI) imaging [17, 18], from which ZrSiS exhibits strongly band-selective scattering depending on impurity lattice site [17]. Moreover, transport measurements revealed extremely large, non-saturated and anisotropic magnetoresistance (MR) in ZrSiS [16, 19–21]. Several Shubnikov–de-Haas (SdH) oscillations measurements reported a major peak at ∼240 T, which is attributed to the quasi-2D Fermi surface [16, 19, 21–25] while a peak at ∼600 T observed by thermoelectric power [23] and high field SdH [24] measurement is assigned to a 3D Fermi surface. Recently, frequency independent flat optical conductivity [25], the rapid relaxation rate after photoexcitation [26] and Ag tip induced superconductivity [27] are also reported in ZrSiS. Thus, these fascinating phenomena demonstrated that the ZrSiS family is a great system for studying the Dirac nodal-lines [11, 12, 28, 29]. However, the experimental investigations on the detail electronic band structure above E_F and the orbital contribution are still inadequate. In addition, a proper model to extract the non-orbital MR contribution is needed to reveal the intrinsic MR contribution due to the Dirac fermions.

In this paper, we investigate the surface and bulk electronic structures in Dirac nodal line semimetal, ZrSiS, by density function theory (DFT) calculations, STM and magneto-transport measurements. Our bulk calculation reveals that the Fermi surface of ZrSiS primarily consists of a 3D diamond-shaped Dirac band at Γ and a quasi-2D tubular band between Γ and X. The slab calculation shows a surface band exists at $\bar{{\rm{X}}}$ with linear dispersion over several hundred meV. A fully spin-polarized surface state at E ∼ 0.6 eV is found to connect the diamond-shaped Dirac band and the bulk conduction band and is completely separated from the bulk bands. Our STM measurements reveal both S- and Si-surfaces and we deduce the observed atoms to be Zr, independent of the surface terminations. Our data from QPI imaging support the band structure from our DFT calculation. Our further analysis suggests that different surface terminations have negligible effect on the surface spectral weight from the Dirac bulk bands in ZrSiS but that of the surface bands can be greatly reduced due to the change in orbital composition. The enhanced surface spectral weight from the diamond-shaped band produces a stronger QPI pattern than that from the surface band. Our finding suggests an unusual scattering selection rule and highlights the important role of Zr-d orbitals. On the other hand, multiple SdH frequencies are clearly identified from the magneto-transport measurements. A model is proposed and applied on the angular dependent MR data to separate the dominant orbital MR contribution and the non-orbital MR that is likely associated with the Dirac bands. Our results thus demonstrate that ZrSiS and its related compounds are a valuable platform for exploring Dirac line nodes physics and their future applications.

Our ZrSiS single crystals are grown by two-step chemical vapour transport process as described in detail elsewhere [30] and have the PbFCl-type crystal structure in tetragonal space group P4/nmm, as showed in figure 1(a). The lattice constants, determined by x-ray diffraction, are a = 3.544 Å and c = 8.055 Å. In ZrSiS, square Zr₂S₂ layers are alternatively stacked along the c-axis with double layers of Si atoms. We calculate the electronic structure based on DFT theory. The bulk band structure without and with SOC is shown in figures 1(b) and (c), respectively. Several Dirac bands can be found in the bulk band structure without SOC (circled in red and blue). The Dirac band crossings also form two one-dimensional (1D) line nodes between R and X at E ∼ −0.45 to −0.67 eV and also between A and M at E ∼ −2.2 eV. When the SOC is taken into account in the calculation, the Dirac bands near the Fermi level (marked as red circles in figures 1(b) and (c)) open a small gap (∼37 meV) but their dispersion remains linear as illustrated in figure 1(c). In contrast, the 1D line nodes (marked as blue circles in figures 1(b) and (c)) in the occupied states remain degenerate due to the protection from non-symmorphic symmetry [9–13, 17, 31]. These Dirac points are connected to form a 3D diamond-shaped Fermi surface at Γ. Moreover, quasi-2D tubular Fermi surface is also found along ΓX as shown in figure 1(f). These results are consistent with recent ARPES [9–13, 16, 31] and SdH [16, 19, 21, 22, 24, 25] measurements on ZrSiS. To better compare the theoretical calculation with STM and ARPES measurements, we also carry out a slab calculation and the resultant band structure with the surface component is shown in figure 1(d). By comparing the electronic structure of ZrSiS from the bulk and slab calculations, we find the surface states with linear dispersion around $\bar{X}$ at E_F, and this is consistent with previous ARPES [9–13] and STM [17] experiments. Interestingly, we also find that the SOC opens a gap of ∼100 meV (orange circles in figures 1(b) and (c)) on the Dirac crossings around E ∼ 0.6 eV and a surface state emerges to connect the top of inner Dirac band and the bottom of conduction bands (orange arrow along $\overline{{\rm{\Gamma }}M}$ in figure 1(d)). Our DFT calculation also shows that this surface state is fully spin-polarized and the spins are all aligned in-plane clockwise as showed in figure 1(e). The surface state is still protected by time reversal symmetry but not by topology since our calculation shows ZrSiS is topologically trivial. From the calculated band structure, the linear dispersion of all Dirac bands in ZrSiS extends over a wide energy range of several hundred meV, which has a great advantage over other Dirac semimetals with linear dispersion of much smaller energy range.

Zoom In Zoom Out Reset image size

We first investigate the surface electronic structure with STM measurements, which are carried out in a homemade low temperature ultra-high vacuum STM system at T = 4.5 K with electrochemically etched tungsten tips [33]. The STM tips are cleaned and reshaped first by field-emission on the gold surface. The ZrSiS single crystal is in situ cleaved at the sample stage of the STM head at around T = 5 K right before STM measurements. Due to the nature of a layered crystal structure with a glide plane between S layers, it is generally assumed that the cleavage should occur naturally here and yield large flat and charge neutral S-surface. Interestingly, we find that ZrSiS can also be cleaved to reveal Si-surface. Figure 2(a) shows a 80 × 80 nm² topographic image of ZrSiS surface with three large terraces. The linecut shows the step height of 8.46 and 3.61 Å, from which we can identify two S-surfaces and a Si-surface as illustrated in figure 2(a) and appendix A3. The high resolution topographic images acquired on S- and Si-surface are shown in figures 2(b) and (c), respectively. Surprisingly, both terminations show the 1 × 1 atomic structure with the same lattice constant of 3.54 ± 0.11 Å and the same orientation, which matches that of S-layer or Zr-layer. When the STM tip is moved closer to the surface by taking the topographic image with higher current, each surface shows different features. Additional atoms become visible at the centre of the square lattice on the S-surface as shown in figure 2(d) while the atoms on the Si-surface become less clear due to hybridization between atoms along a-axis in figure 2(e). Thus, we deduce that the Zr atoms are observed and the surface atoms are invisible from STM in figures 2(b) and (c). The typical differential tunnelling spectra measured on the clean regions of both surfaces show almost identical characteristic semimetallic spectra as illustrated in figure 2(c). The local density of states (LDOS) near Fermi level is the lowest due to Dirac crossing in the band structure but not zero. The measured tunnelling spectra appear to be in excellent agreement with our orbital-projected density of states (PDOS) calculation, as displayed in figure 2(f), in which Zr-d orbital dominates the LDOS near E_F over S-p and Si-p orbitals on the surface. This explains why it is easier to observe Zr atoms than S or Si atoms in STM topographic images and confirms the previous impurities study of ZrSiS by STM and DFT study [17]. The simulated STM images also support this view [17]. A dip exists at E ∼ 0.6 eV in PDOS calculation and is also observed in tunnelling spectrum at E ∼ 0.74 eV, which could be related to the change in energy dispersion as indicated by black and orange arrows in figure 1(d). This dip feature will also allow us to better compare the experimental data with DFT results at the high energy in the following section.

Zoom In Zoom Out Reset image size

To further investigate the dependence of the electronic structure on differently terminated surfaces, we utilize Fourier transform-scanning tunnelling spectroscopy (FT-STS), which can relate the real space characteristic wavevectors resulting from QPI to the band structure in momentum space. First, energy resolved differential tunnelling conductance maps, dI/dV(r, E = eV) are taken in a field of view (FOV) of 60 × 60 nm² on S- and Si-surfaces as shown in figures 3(a)–(e) and (k)–(o), respectively. The standing waves due to QPI are visible around the defects and impurities and show very similar patterns at the same energy. The main difference is that many impurities on the Si-surface do not contribute to any QPI (possibly S clusters, as shown in appendix A4) and the QPI signals diminish at E < ∼125 meV. We then take FT of the dI/dV(r, E) maps and obtain the corresponding dI/dV(q, E) maps on S- and Si-surfaces, as shown in figures 3(f)–(j) and (p)–(t), respectively (the movies, available in the supplementary material, are available online at stacks.iop.org/NJP/20/103025/mmedia). Three pronounced dispersive C₄-symmetric q-vectors (q₁, q₂ and q₃) can be observed in the dI/dV(q, E) maps (figure 3). After analysing the direction, wavelength and energy dispersion of q-vectors and comparing with the calculated band structure, we can identify q₁ as the intra-pocket scattering between the diamond-shaped band at $\bar{{\rm{\Gamma }}},$ q₂ as the intra-pocket scattering of the surface band at $\bar{{\rm{X}}}$ and q₃ as the inter-pocket scattering between the surface bands at $\bar{{\rm{X}}}$ as shown in figure 4(a). Our q-vectors assignments are consistent with previous QPI studies on ZrSiS [17, 18]. However, on the Si-surface, all q-vectors become weaker with decreasing energy and the QPI from the surface bands (q₂ and q₃) diminished at E < ∼125 meV. This is mainly because the disorder from the large amount of impurities on the Si-surface suppresses the QPI. In addition, the surface band near $\bar{{\rm{X}}}$ has significant S-p orbital contribution and the dispersion towards $\bar{{\rm{\Gamma }}}$ and $\bar{M}$ is linked to a change in its Zr-d orbitals [31], which reduces the spectral weight on the surface bands and in term, weakens the QPI near $\bar{{\rm{X}}}$ on the Si-surface. To quantify the QPI data, we then extract q₁ and q₂ from the dI/dV(q, E) maps taken on both surfaces (see figures 4(b) and (c)), in which both q-vectors exhibit the identical energy dispersion on both surfaces within the overlapped energy range and agree well with our calculated band structure. The small difference between our QPI and DFT data at the high energy can be accounted if we take into account an energy-shift by comparing the dip at E ∼ 0.74 eV in tunnelling spectrum with PDOS. Figures 4(b) and (c) show that the diamond-shaped Dirac bulk bands at $\bar{{\rm{\Gamma }}}$ exhibit the linear dispersion from E ∼ −0.1 eV to E ∼ +0.5 eV and the surface bands at $\bar{{\rm{X}}}$ from E ∼ −0.1 eV to E ∼ 0.8 eV. The derived Fermi velocity of the bulk Dirac band is ∼1.6 × 10⁶ m s⁻¹ and of the surface band is ∼8.2 × 10⁵ m s⁻¹, both of which are comparable with that in graphene. Our results demonstrate that the surface termination does not affect the bulk Dirac bands in ZrSiS; however, the spectral weight of the surface bands is greatly reduced due to the change in orbital composition.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To verify the existence of the predicted spin-polarized surface at higher energy, we take dI/dV(r, E) maps between 500 and 1000 meV with higher q-resolution on the same FOV of the S-surface (figure 3(a)) because it will be difficult to resolve smaller q₁ as shown in figures 5(a)–(c). The corresponding dI/dV(q, E) maps are shown in figures 5(d)–(f), in which q₁ vanishes at E > ∼660 meV and appears again at E > ∼860 meV. This energy range coincides with the dip in tunnelling spectrum and PDOS in figure 2(f). Thus, we can conclude that the disappearance of q₁ is consistent with the time reversal symmetry protection from the predicted fully spin-polarized surface band in figures 1(d) and (e). We are unable to make the comparison with data from Si-surface because of the large background at q = 0 due to the large amount of surface impurities. Although such a surface state is not protected by the non-symmorphic symmetry nor the band topology, its full in-plane spin polarization with time reversal protection may be useful in a related compound that has the Fermi level within its proximity.

Zoom In Zoom Out Reset image size

To further study the bulk electronic structure and transport properties, we carry out magneto-transport measurements by standard four-probe geometry. The resistivity equals about $9\,\mu {\rm{\Omega }}\,{\rm{cm}}$ at room temperature and gradually decreases with decreasing temperature, showing a metallic behaviour (figure 6(a)). Below 20 K, the resistivity is nearly independent of the temperature with a residual resistivity of ∼0.23 $\mu {\rm{\Omega }}\,{\rm{cm}},$ giving a residual resistivity ratio (RRR $\equiv {\rho }_{300{\rm{K}}}/{\rho }_{10{\rm{K}}}$) of ∼38.8. Figure 6(b) shows the MR $({\rm{MR}}\equiv [\rho ({\rm{H}})/\rho (0)]-1={\rm{\Delta }}\rho /\rho )$ data at T = 2 K and at several different $\theta$ values with a magnetic field up to 15 Tesla, where $\theta$ is defined as the angle between the magnetic field and the normal direction of the sample surface as illustrated in the upper inset of figure 6(b). Basically, it exhibits a positive MR, but its magnitude does not vary monotonically with $\theta .$ The largest MR turns out to appear at about $\theta$ = 30°, giving a ${\rm{\Delta }}\rho /\rho$ ≈ 173 at 15 Tesla. As $\theta$ further increases, the MR gradually drops to about ${\rm{\Delta }}\rho /\rho$ = 2.7 at $\theta$ = 90° and 15 Tesla. The MR data at $\theta$ = 90° in figure 6(b) are multiplied by 10 for visibility. We note that, in lower fields, MR curves show crossover from concave curvature to convex curvature at $\theta$ = 90°, where practical field linear MR is observed for ${\mu }_{0}H\,\leqslant 1\,{\rm{Tesla}}$ at $\theta$ = 90°. The MR data are replotted in polar figure shown in figure 6(c) at several different fields, which shows a butterfly-like pattern with two-fold symmetry.

Zoom In Zoom Out Reset image size

As shown in figure 6(b), apparent SdH oscillations are observed as the field becomes larger than 3 Tesla. The corresponding fast Fourier transformation (FFT) spectra at different $\theta$ values ranging from 0° to 90° are shown in figure 7(a), where the curves are shifted for clarity and the curves for $\theta$ = 75° and 90° are further amplified by 8-fold and 10-fold, respectively. At $\theta$ = 0°, three peaks of F₁, F₂ and F₃ can be clearly identified from the FFT spectrum as denoted in figure 7(a). The major SdH frequency F₁ equals 240.3 Tesla. F₂ is nearly twice of F₁ and likely to be the second harmonic contribution. F₃, on the other hand, gives a small SdH frequency of about 18.4 Tesla. The corresponding dρ/dH data versus $1/{\mu }_{0}{H}_{\perp }\equiv 1/{\mu }_{0}H\,\cos \,\theta$ at $\theta$ = 0°, 15°, 30°, 45° and 60° are shown in figure 7(b), where beating patterns are observed and justified the multiple frequencies in the FFT spectra. The SdH oscillation period for F₁ is found to scale as $1/{\mu }_{0}{H}_{\perp },$ indicating a 2D-like band structure that is likely deriving from the hole-type quasi-2D tubular Fermi surface shown in figure 1(f). This 2D-like feature is further demonstrated by the dashed line in figure 7(a) and also the excellent linear fitting to the F₁-1/cosθ curve in the lower-left inset in figure 7(b), which is in big contrast to the weak θ dependence of F₃. For magnetic fields along c-axis ($\theta$ = 0°), the hole-type quasi-2D tubular Fermi surface gives an extremal Fermi surface area of A_F = 0.023 084 Ų, which corresponds to a SdH frequency of about 241.8 Tesla according to the Onsager relation (F = (h/4π²e)A_F). This turns out to agree well with the observed major SdH frequency of F₁ = 240.3 Tesla. We note that, as $\theta \geqslant 30^\circ ,$ the major frequency F₁ appears to get broader and then split into multiple peaks, where the peak separation grows larger as $\theta$ increases. This behaviour may be related to the fine structures in the 2D tubular Fermi surface shown as the blue area in figure 1(f). The observed major SdH frequency F1 agrees with several previous reports [24, 29, 34], but we do note that the 3D diamond-shaped band (the green area in figure 1(f)) was not seen in our SdH measurements, which is likely due to its reported unusual mass enhancement [29, 35] and our limited field strength.

Zoom In Zoom Out Reset image size

The presence of multiple bands conduction in ZrSiS is evident from the nonlinear field dependence of Hall resistivity (not shown here), and multiple SdH frequencies. In addition, the calculated band structure also suggests the existence of both electron and hole bands at the Fermi level. To find out the intrinsic MR from the bulk bands, we carry out the following analysis. According to the two-band model, the corresponding orbital MR contribution can be expressed as

$$\begin{eqnarray}&&{({\rm{\Delta }}\rho /\rho )}_{{\rm{o}}{\rm{r}}{\rm{b}}}=\,\displaystyle \frac{{n}_{1}{\mu }_{1}{n}_{2}{\mu }_{2}({\mu }_{2}-{\mu }_{1}){({\mu }_{0}H)}^{2}}{{({n}_{1}{\mu }_{1}+{n}_{2}{\mu }_{2})}^{2}+{({n}_{1}+{n}_{2})}^{2}{\mu }_{1}^{2}{\mu }_{2}^{2}{({\mu }_{0}H)}^{2}},\end{eqnarray} \tag{ 1 }$$

where ${\mu }_{1(2)}$ and ${n}_{1(2)}$ are the mobility and density, respectively, of the first (second) band. ${({\rm{\Delta }}\rho /\rho )}_{{\rm{o}}{\rm{r}}{\rm{b}}}$ is only sensitive to the perpendicular component of the field, and it scales as ${({\mu }_{0}{H}_{\perp })}^{2}$ in lower fields. Based on the Matthiessen's rule, the total MR can thus be expressed as

$$\begin{eqnarray}&&({\rm{\Delta }}\rho /\rho )(\theta ,H)=C(\theta ){({\mu }_{0}H\cos \theta )}^{2}+D(\theta ,H),\end{eqnarray} \tag{ 2 }$$

where ${({\rm{\Delta }}\rho /\rho )}_{{\rm{o}}{\rm{r}}{\rm{b}}}$ = $C(\theta )\,{({\mu }_{0}H\cos \theta )}^{2},$ $C(\theta )$ is a θ-dependent prefactor and $D(\theta ,H)$ term represents the non-orbital contribution to the Mr. Assuming that $D(\theta ,H)$ term grows with field at an exponent less than 2, the value of $({\rm{\Delta }}\rho /\rho )(\theta ,H)/\,{({\mu }_{0}H\cos \theta )}^{2}$ should gradually approach ${\rm{C}}(\theta )$ in higher fields, where $D(\theta ,H)/\,{({\mu }_{0}H\cos \theta )}^{2}$ term becomes vanishingly small. Figure 8(a) shows the $({\rm{\Delta }}\rho /\rho )(\theta ,H)/\,{({\mu }_{0}H\cos \theta )}^{2}$ versus $\theta$ curves at different fields ranging from 0.3 to 15 Tesla, where the curve indeed gradually approaches a field independent function that is shown as a thick red line for ${\mu }_{0}H$ = 15 Tesla. By taking the thick red line as ${\rm{C}}(\theta )$ and using equation (2), we can extract the non-orbital part MR contribution $D(\theta ,H).$ The corresponding $D(\theta ,H)$ versus $H$ curves at different $\theta$ values are shown in figure 8(b). Surprisingly, the $D(\theta ,H)$ is practically linear in H field for ${\mu }_{0}H\,\leqslant 4$ Tesla with a drop in the slope by about 4-fold as $\theta$ value goes from 0° to 75°. Figure 8(c) replots the $D(\theta ,H)$ versus $H$ in a log-scale, where the orange dashed line is the linear fit to the data. We remark that H-linear dependency in $D(\theta ,H)$ persists down to a very weak field of about 0.3 Tesla, and the extracted exponent ${\rm{\alpha }},$ defined as $D(\theta ,H)\propto {H}^{\alpha },$ as a function of $\theta$ is shown in figure 8(d). The exponent $\alpha$ is found to be close to 1.0 and nearly independent of the $\theta$ values, which also justifies the assumption we made earlier ($\alpha$ < 2) for the analysis using equation (2). We do note the failure of H-linearity in $D(\theta ,H)$ for ${\mu }_{0}H$ > 4 Tesla as can be seen in figure 8(c), where we attribute such deviation from H-linearity to the additional contribution from quantum oscillations in the quantum regime. Therefore, the intrinsic MR from the bulk band effect is likely to be best revealed in the low field regime as shown in figure 8(b).

Zoom In Zoom Out Reset image size

The H-linear MR phenomenon has been observed in topological Dirac and Weyl semimetals, such as Dirac semimetal Cd₃As₂ [36] and Weyl semimetal NbP [37, 38]. Combining the observation of linear dispersion band from QPI analysis and also special Dirac line node in the calculated band in ZrSiS, the observed H-linear $D(\theta ,H)$ may provide further support for the existence of the Dirac-like bulk band from the magneto-transport measurements. Nevertheless, the correct theoretical description for the H-linear MR due to Dirac-like band remains a debatable issue [38–41]. On the other hand, it is also suggested that similar H-linear MR effect can be observed in a narrow band gap semiconductor with inhomogeneous conductivity [40, 41]. As shown in figure 8(b), $D(\theta ,H)$ appears to be at maximum at $\theta \,$ = 0°, and the slope d $D(\theta ,H)$/dH seems to roughly scales as cos $\theta ,$ which may rule out the possible Zeeman-like spin spitting effects on Dirac electrons [38]. We also note that the H-linear in $D(\theta ,H)$ remains valid down to weak field regime, suggesting the quantum MR due to the filling of single Landau level may not apply either [39]. Further investigation is needed to elucidate the intrinsic mechanism for the H-linear MR effect in ZrSiS.

In summary, we have reported a combined study on the bulk and surface electronic structures of ZrSiS by using DFT, STM and magneto-transport measurements. Our band calculations show that ZrSiS exhibits 3D Dirac line nodes protected by non-symmorphic symmetry, Dirac surface bands at $\bar{X},$ a fully spin-polarized surface band at $\bar{{\rm{\Gamma }}}$ at E ∼ 0.6 eV. ZrSiS has a 3D electron-like diamond-shaped Fermi surface and a quasi-2D hole-like tubular Fermi surface. Our results from QPI imaging and SdH measurements show overall excellent agreement with the calculated band structures. Moreover, both FSs have strong Zr-d orbital contribution and high carrier mobility, which can contribute the butterfly-shaped MR We analyse our angular dependent MR data with a two-band model and deduce the H-linear non-orbital part of MR at H < 4 Tesla is most likely intrinsic. Our surface termination dependent topographs and QPI imaging demonstrate that the Zr-d orbitals contribution to the surface spectral weight of the Dirac diamond-shaped bulk band is robust against the surface modification whereas the orbital composition of the Dirac surface bands depends on the type of surface atoms. Our results support the importance of Zr-d orbitals on the electronic structures and the unusual Dirac line nodes phase in ZrSiS, which can also be generalized to the related compounds. Thus, the ZrSiS family provides a rich platform to explore the effect of non-symmorphic symmetry and orbital degree of freedom on their electronic structures.

We thank Andreas W Rost, Hsin Lin, Peng-Jen Chen, Tay-Rong Chang, Minn-Tsong Lin, Christopher J Butler, Horng-Tay Jeng and Sungkit Yip for helpful discussions. This work is supported by Academia Sinica (AS-iMATE-107-12), Ministry of Science and Technology, Taiwan and National Taiwan University. T-MC and G-YG are grateful for the generous support from Golden Jade Fellowship of Kenda Foundation.

A1. Computational details

A2. Slab calculation of orbital decomposed band structure from surface atoms on S-surface

A3. The determination of surface terminations based on step height analysis

Topographic images taken in the same FOV and the lines profile along the diagonal direction with (a) bias = −200 mV and (b) bias = +200 mV reveal the identical step heights. Topographic images taken on the yellow and the red terrace exhibit identical bias dependence as shown in figures 2(b) and (d). The measured step height between them is 8.46 Å, which is approximately to 1 unit-cell (uc) of the bulk lattice. The measured step height between the yellow and the blue terrace is 3.61 Å (or 0.43 uc). (c) Three possible scenarios for the cleaved surfaces observed in figure 2(a). The grey regions represent the atoms being removed after the cleavage. From the comparison between the measured step height and these scenarios, we deduce the observed terraces in figure 2 agree with the Case I, which is also consistent with the STM simulation of the S-surface [17]. The Zr-termination may also exist on the cleaved surface since it is on the opposite side of our cleave plane; however, it is absent within our scan range. We note that the bond lengths among Si-Zr tetrahedron, Zr–S within quasi-2D square, and Zr–S across the gap of quintuple layers are 2.812 Å, 2.650 Å, and 2.827 Å, respectively. Thus, the observation of cleaved Si-surface is not surprising from the same level of weaker bond strength reflected on the longer bond length. In fact, the observation of cleaved Si-surface also reflects the semimetal nature of ZrSiS that the Zr–S gap between quintuple layers are not of simple van der Waals type. In addition, the sigma bonding among SiZr₄-tetrahedron requires sp³ orbital hybridization for Si, which is the key controlling factor for the non-symmorphic symmetry of ZrSiS.

A4. Topographic images of S- and Si-surfaces

High resolution topographic images taken in the same 60 × 60 nm² FOV as the differential conductance maps in figure 2. (a) S-surface and (b) Si-surface.

A5. Magnetic field dependent tunnelling spectra on the S-surface