Visualizing discrete Fermi surfaces and possible nodal-line to Weyl state evolution in ZrSiTe

Topological nodal line semimetals (TNLSMs) represent a quantum state of topological matter. When the crystal/time-reversal symmetry is broken, a nodal line state is expected to evolve into a Dirac semimetal, a Weyl semimetal, or other topological phases according to theoretical studies. Here, we report scanning tunneling microscopy (STM) based quasiparticle interference (QPI) measurements performed on the surface of TNLSM ZrSiTe single crystal. A discrete Fermi surface with multiple electron/hole pockets and the impurity-induced inter-/intra- pockets scatterings are directly visualized from QPI patterns. Moreover, the degenerated Dirac points at X point evolve into the pairs of Weyl nodes when Fe atoms are deposited, suggesting a possible phase transition from the nodal line to the Weyl state. The calculated band structures and the Weyl points by applying Zeeman splitting energies along x-direction, further confirm the existence of Weyl points in the Fe-doped ZrSiTe induced by the broken of time-reversal symmetry.


INTRODUCTION
Topological semimetals (TSMs) originate from the crossing of conduction and valence bands in the Brillouin zone (BZ), as shown in Fig. 1a-c. This series of materials have acquired tremendous interest recently because they extend the topological classification of matters from insulator to metal. Initially, it was found that TSMs are characterized by the Weyl/Dirac nodes, with the conduction and valence bands touching at discrete points [1][2][3][4][5][6][7][8] . The low energy excitation behavior of such materials can be described with the Dirac/Weyl equations. Recent studies found that the conduction and valence bands can also cross each other along a line or closed loop in BZ for certain materials, which is classified as the topological nodal line semimetals (TNLSMs) [9][10][11][12][13][14][15][16][17][18][19] .
In general, topological semimetals belong to symmetry-protected topological phases of matter. The crossing of the conduction and the valence bands is protected by the crystal and time-reversal symmetries (TRS). When either inversion symmetry or TRS is broken, a Dirac state can evolve into a Weyl state 20,21 . Likewise, the nodal lines can also break into discrete nodal points or even become fully gapped if certain symmetries are partially or fully broken. For example, without the consideration of spin-orbital coupling (SOC), theoretical calculations show that TaAs is a TNLSM with two nodal lines, protected by mirror reflection and spin rotation symmetries. However, if taking the SOC into consideration, each nodal line of TaAs is gapped into three pairs of Weyl nodes 22,23 . Another example is the breaking of time-reversal symmetry. It has been predicted theoretically that the broken TRS of the fourfold degenerate Dirac point leads to a transformation of the materials from the Dirac semimetals to Weyl semimetals 3,24,25 . Weyl fermions have been shown to exist in the inversion symmetry breaking compounds TaAs 26,27 , NbAs 28 , and TaP 29 . Recently, Weyl nodes were reported in the intrinsically TRS breaking compound, YbMnBi 2 30 . This motivates us to investigate the tunability of the topological nodal phase in the TNLSMs by the breaking of the symmetries. However, as far as we know, the evolution of a topological nodal line semimetal to a Weyl semimetal phase has not been achieved experimentally.
ZrSiS is found to be a new TNLSM material, with negligible SOC 2,12,31-35 . The SOC strength can be enhanced by replacing the S atom with heavier chalcogen ions, such as Se and Te. Such a chemical substitution can maintain a similar crystal structure. Theoretical calculations predicted an approximately 60 meV gap in ZrSiTe around the nodal line direction (i.e., M-Γ-M) in the BZ. Angleresolved photoemission spectroscopy (ARPES) has been utilized to investigate the nodal band structure in the ZrSiX (X = S, Se, Te) materials 2,18,36,37 . Scanning tunneling microscopy (STM) based quasiparticle interference (QPI) enables powerful insights into the k-space properties of superconducting, topological, Rashba, and other exotic electronic phases, and has been also performed to study these materials. Previous work reported the observation of the floating band state on the ZrSiSe (001) surface 38 and similar QPI measurements were also performed on the ZrSiS surface, which suggests the selective scattering from different subsets of the Zr 4dderived bands and pseudospin conservation at energies close to the line node 35,39 .
Here, we show a series of STM-based QPI measurements performed at the energies near Fermi level (E F ) on the vacuumcleaved surface of a high-quality ZrSiTe single crystal. By replacing S/Se with heavier chalcogen ion Te in the ZrSiX family, ZrSiTe possesses the heaviest chalcogen ions and hence the strongest SOC strength 2,18,36 . This SOC effect can effectively gap most of the band crossings inside the BZ, except for the ones along X-R directions 2,18,36,37 . Our DFT calculation in Fig. 2 indeed suggests that the Dirac crossings along ΓX and ΓM that form a line node around the Fermi level in the non-SOC case, are now gapped in the presence of the SOC, while the degeneracies at X, R, A, and M are unaffected due to the non-symmorphic protection. The impurityinduced QPI patterns in Fig. 3  surface with multiple electron/hole pockets and the impurityinduced inter-/intra-pockets scatterings. The bias-dependent QPI images provide clear evidence for the subtle band crossings near E F with astoundingly high resolution. Moreover, the breaking of TRS is achieved by introducing Fe atoms on the ZrSiTe surfaces. The degenerated Dirac crossings, as a consequence, develop into pairs of Weyl nodes near the X points. The Fermi arcs connecting the Weyl nodes were observed in the FT-QPI images. Our findings provide a practical scenario to tune the nodal line phase in the ZrSiX materials.
RESULTS AND DISCUSSION STM images, defects, and tunneling spectra on the ZrSiTe surface The crystalline structure of ZrSiX (X = S, Se and Te) has the tetragonal P4/nmm space group, adapting the PbFCl structure type. The Si square nets are separated by two layers of Zr-X, as depicted in the inset of Fig. 1d. The exposed (001)-oriented cleavage plane is then expected to be between the two adjacent X atomic layers 12,18 . Figure 1d gives the representative topographic STM image of the cleaved surface obtained at −20 mV. The square bright array represents the Te lattice on the surfaces, suggesting there is no obvious reconstruction during the cleaving and the subsequent cooling process. In addition to the periodical lattice, two types of defects are obviously visible on the surface, which we refer to as α and β types, respectively. Highly-resolved STM images of the two types of defects are given in Fig. 1e, f. Based on previous calculation 38 , defects α looks like dark holes, which can be attributed to the Te vacancies. Defects β, however, appear as bright protrusions. Since the exposed surface was prepared under ultrahigh vacuum (UHV), we assign defects β to the Te adatoms. Te adatom (f), respectively. The size is 4 nm × 4 nm for both images. g Hundreds of dI/dV spectra (green) measured at different locations on the bare surface areas (away from the defects) and the averaged spectrum (black) (V bias = 20 mV, I set = 300 pA, lock-in V mod = 3 mV). h STS curves were measured in different areas (marked as α, β, and γ in Fig. 1d). All the data are measured at 4.2 K.
Subsequent tunneling spectroscopic measurements are carried out in order to investigate the electronic structure of ZrSiTe. The typical dI/dV spectrum taken on the bare ZrSiTe surface (Fig. 1g), similar to the STS results on ZrSi(S/Se) 35,38,39 , except for a symmetrically-placed gap-like feature at E F , and the gap size is about 19 meV. A similar symmetric gap-like feature is also observed on the ZrSiSe surface, and the gap size is about 6 meV, as shown in Supplementary Fig. 1. Such tunneling spectra have been measured hundreds of times and the main feature maintains at different locations that are far from the defects (Fig. 1g), implying such density of states (DOS) is quite universal over the entire surface. By comparing such symmetrically-placed gap-like feature in dI/dV spectrum with the previous report in graphene 40 , we attribute the observed 19 meV gap at E F in Fig. 1g to a surface phonon-mediated inelastic tunneling process (±9.8 meV) 41 . For comparison, we measure the dI/dV spectra above the α and β types of defects, respectively. As shown in Fig. 1h, the overall shape of the spectra is similar to that obtained on the bare region (spectrum γ), except that defect α exhibits more unoccupied states to accept the tunneling electrons, while defects β has less unoccupied states. The evolution of the spectra across both defects can be seen in Supplementary Fig. 2. Previous ARPES measurements reveal that the non-symmorphic Dirac cones of ZrSiTe exhibit point-like   Supplementary Fig. 4c. d, f The QPI patterns are the same as that shown in b. The arrows represent the scattering q-vectors. e, g Fermi pockets with regions of high spectral density are marked as red dots and the arrows indicate the scattering pathways between these red dots. The scatterings between these dots are responsible for the q-vectors in the FT-STS images, and the corresponding relationship is d to e, f to g, respectively. h The dispersions of the q 1 , q 2 , q 3 , q 5, and q 9 vectors. The slopes of these dispersions are ranged from 1.5 to 4 eVÅ. The zoom-in high-resolution images (4 nm 2 × 4 nm 2 ) for a single Fe impurity were taken with +100 and +40 mV, respectively. d dI/dV image was taken simultaneously with the panel a with the bias voltage of +100 mV. e The corresponding Fourier transforms pattern. f dI/dV image taken simultaneously with the panel a with the bias voltage of +40 mV. g The corresponding Fourier transforms pattern. Note the Fourier transforms pattern in e, g are original without the symmetrized. The red square indicates the size of the first Brillouin zone. h Zoomed image of the yellow square region in FFT image of g. We use red and brown dashed curves to highlight the Fermi arcs. i Schematic drawing shows that each Dirac point located at X, splits into a pair of Weyl points (purple dots). Green Fermi arcs can be deducted directly from panel h. j The calculated band structure and Fermi contour (k) in the first BZ with two pairs of WPs, after 5 meV of the ZM splitting energy was added along the x-direction.
features at the X point, with the top cone located at about +40 meV and the lower one located at around −300 meV, respectively 36 . Among our dI/dV spectra, no prominent feature is observed at these characteristic energies, simply due to these Dirac points hidden inside of other bulk bands.

Band structures, nodal lines, and Fermi surfaces of ZrSiTe
The calculated band structures of ZrSiTe without/with the SOC effect are shown in Fig. 2a, b, and the calculated nodal lines are shown in Fig. 2c, d correspondingly. By ignoring the SOC effect, the nodal lines, as shown in Fig. 2c, form a "nodal-cage" circle around the Г-Z line and two nodal lines along the X-R direction in the Brillouin zone (BZ). The 'nodal-cage' has an energy broadening at different k points, as indicated by the different colors of the energy bar from about −0.35 to 0.45 eV. However, the nodal lines at the BZ boundary along the X-R direction are very close to the Fermi level. When the SOC effect is included, the "nodal-cage" are all gaped, and only the nodal lines along the X-R direction that are protected by the non-symmorphic symmetry are preserved as shown in Fig. 2d, where all the nodes between ±0.5 eV around the Fermi level are counted. Along the X-R line k = (0, π, k z ), the little group is C 2v , and the generators of this little group are C ' 2z = { C 2z | (1/2, 1/2, 0) }, M ' x = { M x | (0, 1/2, 0) } and PT symmetry operation, and these operations satisfy (PT) 2 It is direct to prove that the irreducible representation of such a little group must be four-dimensional even in the double-space group. So that the Dirac nodal lines along the X-R survive when the SOC effect is included. The Dirac crossings along ΓX and ΓM that form a line node around the Fermi level in the non-SOC case, are now gapped in the presence of the SOC, while the degeneracies at X, R, A, and M are unaffected due to the non-symmorphic symmetries. Furthermore, the Fermi surfaces of ZrSiTe are calculated and shown in Fig. 2e, f, where the discontinuous Fermi surfaces in the first BZ and nodal lines at the BZ boundary along X-R direction are clearly seen. In order to compare with our QPI measurements, the Fermi contour at k z = 0 is also calculated and shown in Fig. 2g. Based on previous ARPES measurements 18,36 and our calculation, the Fermi contour of ZrSiTe in the first Brillouin zone is drawn schematically in Fig. 2h, which includes a diamond-shaped Fermi pocket around the Γ point and extra four semicircular-like pockets. Such Fermi contour includes multiple bands and small Fermi pockets crossing E F . As we see, the diamond-shaped Fermi surface consists of red pockets at the corners, yellow pockets in the middle edge of the diamond shape, and blue pockets contributing to extra four semicircles.

Fermi surfaces and inter-/intra-pockets scatterings in QPI
The discrete Fermi surface and the scattering between multiple Fermi pockets can be revealed by QPI imaging. Figure 3a gives a representative dI/dV conductance mapping of the cleaved ZrSiTe surface with a bias voltage of −20 mV, acquired simultaneously with the topographic STM image shown in Fig. 1d. Fourier transforming of Fig. 3a displays a rich array of q-vectors in the momentum q-space, as shown in Fig. 3b, suggesting the existence of abundant scattering pathways. For clarity, we mark the scattering vectors with colored arrows, as illustrated in Fig. 3d, f. The scattering q-vectors connect the regions with high spectral density in the contours of constant energy (CCE). The constant energy used to acquire QPI patterns is chosen to be close to the E F . For brevity, we sort the scatterings into two groups, as shown respectively in Fig. 3e, g, in which the arrows represent these scattering vectors. After carefully examining the length and direction of these scattering q-vectors, all the q-vectors (colored arrows) can perfectly match with the corresponding high spectral density regions (red dots) as shown in Figs. 3e-g. This indicates that our QPI pattern matches well with the Fermi surface. The discrete high spectra density regions, implies the presence of electron and hole pockets after the 'nodal-cage' circle around the Г-Z line are gapped by SOC. Notably, these discrete high spectra density regions belong to different Fermi pockets, therefore the observed impurity-induced scattering processes involve the intraband and inter-band scattering.
To obtain the dispersion of these scattering q-vectors, a series of bias-dependent dI/dV images (see Supplementary Fig. 4), and the corresponding QPI patterns (see Supplementary Fig. 5) are acquired. The scattering spots are sharper and more discrete at lower bias voltages, suggesting that the nodal band crossing is more discontinuous at energies that are close to the E F . When the energies are away from E F , the DOS distributes more uniformly on the whole Fermi pocket in k-space. This phenomenon confirms the discrete characteristic of Fermi pockets in ZrSiTe. Besides, due to ultrahigh resolution in momentum space, some subtle band information can be visualized, e.g., comparing QPI patterns at −20 (Fig. 3b) and −60 meV (Fig. 3c), the q 4 vector moves outwards, indicating four semicircular pockets shrinks with the energy shifting away from E F . Furthermore, we extract the energy dispersion of q-vectors from the bias-dependent QPI images and systematically examine their energy dependence. As shown in Fig. 3h, the q 1 , q 2 , q 3 , q 5 , and q 9 vectors show the nearly linear dispersions within the energy windows from −60 to +60 meV. The slopes of these dispersions can be estimated in the range of 1.5 to 4 eVÅ, which is roughly in consistence with previous ARPES observations about the diamond-shaped Fermi surface formed by the bulk bands crossing in ZrSiX (X = S, Se and Te) and our calculations in such energy range. As reported in ref. 36 , the Fermi velocities of the surfacederived bands in ZrSiTe are ħv F = 3.9 and 5.2 eVÅ, respectively. Our slope values of the linear dispersion are roughly in consistence with these published data. It is worth noting that the q 1 , q 2 , q 3 , and q 9 vectors denote the inter-band scattering, while the q 5 is intra-band scattering.

Fermi arcs and the broken of TRS induced by Fe deposition
In order to further tune the topological nodal phase in ZrSiX materials, we deposited a small amount of Fe atoms onto the cleaved ZrSiTe surface, by which the TRS could be effectively broken by the magnetic moments of Fe impurities. Figure 4a-c shows a topographic STM image on which several isolated Fe atoms can be identified. No other kinds of point defects or step edges are visible in this region so that Fe atoms act as scattering centers solely. Figure 4d, f give the dI/dV images acquired at +100 mV and +40 meV, respectively. By Fourier transforming these dI/dV images, we find that the QPI image taken at 100 mV exhibits a featureless pattern in its FFT image (Fig. 4e), except for Bragg spots of a surface lattice. However, the situation is dramatically different for the QPI pattern taken at +40 mV (Fig. 4g), which is known to be the energy of the Dirac crossing of ZrSiTe 36 . As shown in Fig. 4g, h, four sections of arcs are clearly observed in the FFT image, and we highlight two of them with red color and the other two with brown color. We present several cuts of the arcs along the yellow dashed lines in Supplementary Fig.  10a, as shown in Supplementary Fig. 10b. The black arrows denote the intensity peaks clearly associated with the appearance of arcs. Note these arcs can only be observed at the energies close to the energy of the Dirac cone (see Supplementary Fig. 6 for details). Subsequent converting the q-vectors into k-vectors in the momentum space leads to four Fermi arcs (green arcs) appearing around the X/Y point of the BZ boundary, as schematically shown in Fig. 4i. The detail of conversion from q-vector to k-vector is shown in Supplementary Fig. 10c-f. Contributions to the QPI features from both arcs were superimposed in Supplementary  Fig. 10f 42 . Such features represent the shape of a ninja star, which is consistent with the experimental imaging in Fig. 4h. Such QPI measurements were repeated on different Fe-doped samples, and the consistent QPI pattern with better quality and similar Fermi arcs are observed on the >0.33 at.% Fe-doping sample (see Supplementary Note 9 and Supplementary Fig. 11). In our case, the introduction of magnetic iron atoms breaks TRS and thus lifts the fourfold degeneracies of Dirac points. Each Dirac point located at X point of the BZ boundary, therefore, evolves into the pair of Weyl points, given by the endpoints of Fermi arcs.
In order to further confirm the emergence of these Weyl points, we perform a DFT calculation (Fig. 4j, k and Supplementary Fig. 7). We first construct the tight-binding (TB) Hamiltonian by the symmetry-adapted Wannier functions for the Zr d, Si p, and Te p orbitals without performing the procedure of maximizing localization in coded WANNIER90. Then the Zeeman (ZM) splitting energy along the x-direction is added to the TB Hamiltonian to search the Weyl points (WPs) by using WannierTools. It is worth noting that, several pairs of Weyl points are really produced near the X point in BZ with a moderate Zeeman split energy, which exhibits a good agreement between experiment and theoretical prediction. Figure 4j, k shows the calculated band structures and the appearance of two pairs of Weyl points (WPs) after applying the 5 meV of Zeeman (ZM) splitting energy along the x-direction. When 7 meV of the ZM splitting energy is added, as shown in Supplementary Fig. 7a, six pairs of WPs have been found within ±0.1 eV relative to the Fermi level. We note that the splitting of WPs in k z = 0 plane in Supplementary Fig. 7a is very small. When increasing the ZM splitting energy, more WPs would appear according to our calculations. For example, 12 pairs of WPs within ±0.1 eV of the Fermi level are found, as shown in Supplementary  Fig. 7b.
Our method provides a practical attempt to tune the degenerated Dirac crossings in ZrSiTe into the pairs of Weyl nodes by doping with magnetic impurity (broken of TRS). One may question whether a rather low coverage of magnetic impurities would be sufficient to provide some kinds of magnetic order which would cause a global breaking of TR-symmetry. A small amount of adsorbed magnetic adatoms may only break locally the TR-symmetry of electrons through the scattering near these magnetic defects. In fact, scientists have previously found the ability to break the TRS in topological insulators with magnetic doping concentrations well below 1% [43][44][45] . Both Cr-doped and V-doped topological insulator (Bi, Sb) 2 Te 3 thin films have been used to realize The quantum anomalous Hall (QAH) effect 46,47 , which proves that a low coverage magnetic doping could induce the global magnetic ordering. Indeed, intrinsically TR-symmetry breaking (i.e., intrinsically magnetically ordered) topological semimetallic materials, will be more promising candidates for realizing the phase transition from the nodal line to the Weyl state.
In summary, we have investigated the surface topography and electronic structures of topological nodal line semimetals ZrSiTe, utilizing low-temperature scanning tunneling microscopy and Fourier transform scanning tunneling spectroscopy. Our QPI patterns match the Fermi contour perfectly and reflect the detail of electron and hole pockets with high resolution in k-space. Furthermore, the breaking of TRS by introducing Fe atoms on the ZrSiTe surface leads to the degenerated Dirac crossings developing into the pairs of Weyl nodes near the X points, which are evidenced by the observation of Fermi arcs at the Dirac energy. Our findings provide a practical route to operate potentially the nodal phases between nodal lines, Dirac points, and Weyl points in the ZrSiX material system, by turning on the SOC and further the breaking the time-reversal symmetry.

Samples and STM measurements
The ZrSiSe/ZrSiTe crystals were synthesized using a chemical vapor transport method as described elsewhere 12 . The studied crystals have the typical sizes of 2 mm × 2 mm × 0.5 mm. Samples were cleaved in situ at room temperature (RT) under a vacuum with pressure better than 1 × 10 −10 torr. The cleaved sample was quickly transferred into an Unisoko-1300 commercial STM for measurement at a temperature of 4.2 K. Fe atoms were deposited onto the surface with the evaporation temperature at 1075°C with a water-cooling K-cell evaporator. STM measurements were conducted using an electrochemically etched tungsten tip. dI/dV spectroscopy curves were acquired by employing a lock-in technique with a bias modulation.

QPI imaging
To obtain the QPI imaging, differential conductance dI/dV (r, eV) ≡ g (r, E = eV) was firstly acquired at locations r and bias voltage V. After scanning a certain area, we get a differential conductance mapping g (r, E = eV) proportional to the local density of electronic states LDOS (r, E = eV). The mapping displays spatial modulations arising from an elastic scattering of quasiparticles. Subsequently, this mapping can be Fourier transformed to extract scattering vectors (q-vectors), which connect electronic states in kspace via a set of selection rules. The symmetrization of Fourier transform patterns is done by using four-quadrants averaging.

DFT calculation
The first-principles calculation was performed by the density functional theory as implemented in the Vienna ab initio simulation package (VASP) 48 . Perdew-Burke-Ernzerhof (PBE) 49 type of the generalized gradient approximation (GGA) was used as the exchange-correlation potential. The cutoff energy for the wave function expansion was set as 400 eV, and 11 × 11 × 9 k-meshes were used. Spin-orbit coupling (SOC) was considered consistently in the calculations. A real-space tightbinding Hamiltonian was obtained by constructing the symmetryadapted Wannier functions for the Zr d, Si p, and Te p orbitals without performing the procedure for maximizing localization in coded WANNIER90 50 , which was used to calculate the topological properties of this system by using WannierTools 51 .

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.

CODE AVAILABILITY
The codes that were used here are available upon request to the corresponding author.