Introduction

Recently, topological nodal line semimetals (TNLSMs) ZrSiX (X = S, Se, and Te) have attracted extensive attention. In TNLSMs, the conduction and valence bands cross at several points in k-space with forming a closed loop on the Fermi surface. These loops are protected by the nonsymmorphic symmetry1,2,3, and the electronic states along the nodal lines are topologically nontrivial. The crystalline and band structures of ZrSiX have been studied by means of X-ray diffraction4 and angle resolved photoemission spectroscopy (ARPES)5,6,7. The unique band structure of ZrSiX leads to various interesting physical properties, e.g., chiral anomaly, extremely large magnetoresistance, high-carrier mobility with high-carrier density, etc5,8. Due to the weak interlayer interaction9, ZrSiTe shows some two-dimensional (2D) characteristic in its band structure, revealed by quantum oscillation in magnetic susceptibility and magnetoresistance8. More interestingly, ZrSiTe was theoretically predicted to be a topological insulator when it is thinned down to the monolayer limit1,9,10. Although the physical properties of bulk ZrSiTe have been studied extensively, those of mono- and few-layers are still unclear.

For most 2D materials such as graphene, transition metal dichalchogenides (TMDs), and black phosphorus (BP), Raman spectroscopy provide a powerful tool to investigate the evolution of crystalline symmetry and interlayer interactions during a thinning down process11,12,13,14,15,16,17,18,19,20,21. Defects in 2D materials may induce a shift in frequency and an asymmetrical broadening of Raman peaks22,23. In addition, line defects, e.g., phonons at layer edges, can also be detected in black phosphorus24 and MoS225 by polarized Raman spectroscopy. However, the momentum conservation rule requires q ≈ 0 for a Raman process, thus phonon modes away from the Γ point in the q-space are hardly involved in a sole Raman scattering process26. For example, Hong et al.27 studied point-like defects in monolayer MoS2 and obtained nearly-perfect Raman spectra without showing any signature of defects, whereas the sample shows plenty of point defects in scanning transmission electron microscopic (STEM) images. Although it was recently demonstrated that photoluminescence (PL) might be another feasible method for detecting single-atom point defects28,29, it was designed for semiconductors and does not work for (semi-) metals. It therefore remains a huge challenge to quickly detect single-atom point defects using optical methods, especially in metallic, layered materials.

In this work, mono- to few-layered ZrSiTe films were exfoliated on a 285 nm-thick SiO2/Si substrate with the Scotch tape method. A systematic Raman study was performed in order to investigate the layer-dependent lattice vibration properties. Our results show significant thickness dependences of Raman peak position and intensity in ZrSiTe mono- and few-layers. In all thicknesses, three major Raman modes (around 95 (\(A_{1g}^2\)), 228 (\(A_{1g}^1\)), and 304 cm−1 (\(B_{1g}^1\))) and three minor modes are identified. Raman-active vibrational modes of 2D ZrSiTe were also modeled with density functional theory (DFT) calculations. All experimental observed modes except for the one located at 131.7 cm−1 are found to be the intrinsic vibration modes related to the crystal symmetry of ZrSiTe by our DFT calculations. The 131.7 cm−1 mode does not exist in the list of Raman-active modes predicted by the calculation of defect-free layers but can be understood as a combination of the intralayer shear modes at the Brillouin zone (BZ) boundary activated in thin layers. In long-stored samples, the defect-induced localized modes around 120 cm−1 also contribute to the 131.7 cm−1 mode and induce a red-shift of it. The total intensity of this mode can be further enhanced by 4.8 times upon increasing defect concentration.

Results and discussion

ZrSiTe has a crystalline structure in the tetragonal PbFCl type (P4/nmm space group No. 129)4,5,9,30. In the lattice of ZrSiTe, double layers of Zr–Te are separated by a layer of a Si square net. Two adjacent Te layers (Fig. 1b) are expected to be separated upon a cleavage along the (001)-orientation. Compared to its analogues, ZrSiS and ZrSiSe, the interlayer interactions in ZrSiTe are much weaker, evoking the possibility of isolating ZrSiTe monolayers from its bulk crystals by mechanical cleavage8,9. As a result, the thickness-dependent evolution of ZrSiTe lattice structure and physical properties could be studied with a layer-by-layer resolution. Figure 1a shows a typical atomic force microscopic (AFM) image and the corresponding optical image (inset of Fig. 1a) of the ZrSiTe mono- to few-layers, deposited on a 285 nm-thick SiO2/Si substrate. The ZrSiTe monolayer is light purple and the color quickly gets darker as the thickness of ZrSiTe gets larger from 2 to 6 layers. Then, the color of ZrSiTe film turns blue, green, beige and golden with respect to the increasing thicknesses. Similar color codes have also been found in TMDs such as MoS2, MoTe2 and WTe215,18,31,32, suggesting that the color of ZrSiTe on the SiO2/Si substrate is a good indicator of its thickness. The height profile (Fig. 1c, d) measured along the solid line marked in Fig. 1a indicates that the thickness of a single layer is 1.3 nm. This thickness value is slightly larger than the theoretical value 9.5 Å revealed by our DFT calculation for the bulk crystal, but is still reasonable16. In addition, the monolayer thickness of 1.5 nm is slightly larger than the interlayer distance of ZrSiTe. This phenomenon is usually observed in AFM measurements on exfoliated 2D materials, such as graphene33, MoS215,34,35, and other dichalcogenides18.

Fig. 1
figure 1

Topology of a typical ZrSiTe flake on a 285 nm SiO2/Si substrate. a The optical (the inset of a) and AFM images of a ZrSiTe flake on a 285 nm SiO2/Si substrate. The length of the scale bar is 1 μm. The variation of the thickness can be identified layer-by-layer both in optical and AFM topography images. b Crystalline structure of ZrSiTe. The cleavage plane between adjacent Te layers is indicated by an arrow. The line profile shown in c and d indicates the height profile along the blue and red solid line, respectively. The dashed lines indicate the average height of the corresponding plane of the layers. The thickness of a single layer can be estimated to be 1.3 nm from these line profiles

Full size image

Figure 2a compares the Raman spectra of freshly-cleaved ZrSiTe thin films of various thicknesses, in which the position and intensity of Raman signals consistently vary with the thickness. The sharpest peak is around 228 cm−1 (\(A_{1g}^1\)) in all samples except the monolayer one. Both this peak and the one around 304 cm−1 (\(B_{1g}^1\)), are symmetric and can be fitted with single Lorentzian peaks, indicating the absence of strong electron–phonon and spin–phonon couplings in these phonon modes. Our DFT calculations reveal the details of these three observed Raman-active modes. Figure 2e–g illustrates the lattice vibration displacements of these three Raman-active modes at 96.9 (\(A_{1g}^2\)), 226.4 (\(A_{1g}^1\)), and 299.5 cm−1 (\(B_{1g}^1\)) in a ZrSiTe bilayer, respectively. All of them are out-of-plane vibration modes.

Fig. 2
figure 2

Layer-dependent Raman spectra and phonon vibration modes in ZrSiTe. a The Raman spectra of fresh ZrSiTe with thickness ranging from 1L to bulk ZrSiTe and b the evolution of the representative Raman peak positions. The Raman modes in the monolayer ZrSiTe are distinct from the modes in the thicker ones. The sharpest mode in monolayer is at 131.7 cm−1 instead of around 228 cm−1 as in thicker layers. c The Raman spectra of ZrSiTe with various thickness after one-week storage in nitrogen gas-filled glove box, and d the representative Raman peak shifts. The Raman modes changed markedly in both monolayer and few layers. The modes at 131.7 and 142.7 cm−1 were enhanced in the monolayer, and the wave packet around 150 cm−1 arose in few layers. e–g The three representative eigen lattice vibration modes in 2L confirmed by the first-principle calculations. h, i The electronic and the phonon band structures of the monolayer ZrSiTe. The calculations were performed with considering spin–orbit coupling (SOC). The illustration in h shows the Brillouin zone and the path of high symmetry points. The three red horizontal lines in i represent the frequency of 66.7, 131.7, and 142.7 cm−1

Full size image

As shown in Fig. 2b, two Raman modes, i.e., \(A_{1g}^2\) and \(B_{1g}^1\), are significantly stiffened in thicker samples while the \(A_{1g}^1\) mode slightly softens. Layer-dependent frequency shifts of Raman modes were frequently reported in other 2D materials11,12,13,14,15,17,18,19,20,21,31,32. According to our DFT calculations, all these three modes are assigned to out-of-plane vibrations (Fig. 2e–g). While it is straightforward that the stiffening of modes \(A_{1g}^2\) and \(B_{1g}^1\) with respect to increasing sample thickness is ascribed to stronger interlayer interactions in thicker layers, the softening of the \(A_{1g}^1\) mode is rather complicated. This mode mainly represents the out-of-plane vibrations of Zr atoms, which is influenced by the neighboring surface Te atoms. Surface effects, which weaken with the increase of thickness, are significant in varying vibrational frequencies36,37 of 2D materials. These surface Te atoms lead to appreciable softening of the \(A_{1g}^1\) mode in ZrSiTe with thickness from 1L to 4L, and the softening becomes not significant above 4L, as shown in Fig. 2a. Table 1 summarizes the experimental and theoretical frequencies of the Raman-activated modes in 1L to 4L ZrSiTe. All experimental frequencies were well reproduced by the theoretical calculations within an error of 3%, indicating the reliability of the DFT calculations. An exception lies in the \(A_{1g}^2\) mode in the monolayer, where the experimental frequency is 4.5% higher than the theoretical value. It is worth noting that in thicker samples, e.g. 4L, more Raman modes could be observed. And the numbering of modes based on their phonon energies may vary. In this manuscript, however, we use the labling of modes in the bilayer ZrSiTe for all thicknesses for convenience. In addition to the peak position shifts, the overall intensity of Raman signals gets lower if the layer thickness increases, which is most likely ascribed to stronger light absorption in thicker ZrSiTe layers. In fact, it is difficult to acquire effective Raman signals in ZrSiTe samples thicker than 100 nm. The spectrum shown in Fig. 2a is measured from a roughly 50 nm-thick film, denoted as a bulk result.

Table 1 Calculated and experimentally energies and symmetries of typical Raman modes in 1 to 4 layer ZrSiTe
Full size table

It is worth to note that the Raman spectrum of the ZrSiTe monolayer is distinct from those of other thicknesses. In the monolayer, the strongest mode is centered at 131.7 cm−1 instead of 304 cm−1. This mode fades out with respect to the increasing layer thickness, and cannot be clearly observed in the samples thicker than 3L. According to the DFT calculations, this mode is not supposed to be Raman-active in ideal defect-free ZrSiTe monolayers. It was reported that tellurium thin films show two representative Raman modes at 120 and 140 cm−1, respectively38,39,40. However, the peak position of the mode at 131.7 cm−1 we found in ZrSiTe is appreciably away from the frequencies of both tellurium modes. In addition, the two representative tellurium modes can be observed in all thicknesses, inconsistent with the fact that the 131.7 cm−1 mode fades out with increasing thickness of the ZrSiTe layers. Also, the degraded tellurium is more likely to be amorphous instead of forming the crystal. As a result, the 131.7 cm−1 peak observed only in 1L and 2L ZrSiTe samples is unlikely to originate from the modes of α-phase tellurium few-layers. In 2D materials, defects often have significant influences on Raman spectra11,12,41. It is reasonable to infer that the 131.7 cm−1 mode in thin layer ZrSiTe samples is related to defects.

We performed scanning tunneling microscopic (STM) measurements on an in situ cleaved ZrSiTe surface. Figure 3a shows a topographic STM image of the cleaved surface. The surface Te lattice was visualized as a bright dot array. The nearest Te–Te distance was measured to be 3.7 Å, which is consistent with the theoretical value of 3.68 Å4, suggesting no obvious surface reconstruction occurred during the cleavage. In addition, three types of defects, labeled as α, β, and γ respectively, were clearly observed on the surface. High-resolution STM images of these three types of defects are given in the right panels of Fig. 3a. The α defect appears as a dark hole. The β defect appears as a bright protrusion with a star-like shape. And the γ defect shows a similar bright protrusion with β, but in a square shape. Our statistics show that the β defect, with a concentration of 2.6‰ estimated by counting the number of defects in all scanned areas, is the most common defect on the ZrSiTe surface.

Fig. 3
figure 3

Topographic STM images for the various surface defects in ZrSiTe. a A typical 30 × 30 nm2 atom-resolved STM image of in situ exfoliated ZrSiTe surface with defects. The zoom-in of the high-resolution STM image for the defects of types of α, β, and γ are shown in the right. The corresponding defects are marked with blue squares in the large image (a). The length of the scale bar is 6 nm in the main panel and 1 nm in the right zoom-in panels. b–d Top views of the ball-stick models of a 3 × 3 monolayer ZrSiTe with a Te vacancy, a Te adatom, and a Zr vacancy, respectively. e–g Simulated STM images of the defects shown in b–d, respectively. Compare a with e–g, the defect types of α, β, and γ can be evidently ascribed to Te vacancy, Te adatom, and Zr vacancy, respectively

Full size image

In order to identify these atomic defects, three most probable types of defects, i.e. Te vacancy, Te adatom and Zr vacancy, were modeled by our DFT calculations. Figure 3b–d shows the top view of a 3 × 3 lattice ball-stick model of the ZrSiTe monolayer with these three types of defects. Meanwhile, Fig. 3e–g exhibits the corresponding simulated STM images, which are highly consistent with the experimental images of the α, β, and γ defects, respectively. Moreover, in the tellurium-rich limit, the formation energy is 1.3 eV for the Te adatom, 2.8 eV for the Te vacancy and 3.8 to 4.6 eV for the Zr vacancy. In contrast, in the zirconium-rich (tellurium-deficient) limit, the formation energy is 2.1 to 1.3 eV, 2.0 to 2.8 eV and 4.6 eV for these three types of defects, respectively. Therefore, the most common β defect is assigned to the Te adatom, because of the high consistency between the simulated and experimental results, and the lowest formation energy. Based on similar derivations, we assess defects α and γ to the Te vacancy and the Zr vacancy, respectively.

In order to further understand the relationship between the Raman signals and the defect related lattice vibration modes of the ZrSiTe layers, we performed density functional perturbation theory (DFPT) calculations for ZrSiTe mono- and few-layers. According to our DFPT results, the 131.7 cm−1 mode found in thinner layers is most likely resulted from atomic defects assisted or induced Raman scattering, rather than an eigen-vibrational mode at the Γ point (see Supplementary Information Fig. S1a). It was reported that atomic defects could provide additional inelastic scattering processes, leading to the zone boundary phonons modes being observable in Raman spectra, like the D peak in graphene11. Here, we show the phonon dispersion relation of the ZrSiTe monolayer in Fig. 2i, which indicates a mode with a rather small dispersion along X-M exhibiting a frequency of 133.0 cm−1 at the X(0.5,0,0) point on the boundary of the BZ. Figure S2b illustrates the vibrational displacements of this mode, which suggest an intralayer shear vibration inducing rumpling of the layer. This rumpling is activated in thinner layers where the perpendicular interlayer interactions are not stiff enough. It might be squeezed to a nearly in-plane mode when the thickness is significantly large, resulting in severely depressed Raman signals around 130 cm−1 in thicker layers (Fig. 2a). Similarly, some ‘forbidden’ modes can be also observed in the 2D materials with the CdI2 crystal structure (e.g., HfS2, SnS2, PtSe2)25,42,43,44,45,46. However, unlike defect-induced Raman modes in our case, they are infrared-active modes or secondary-order Raman modes presented in the spectra owing to resonant Raman scattering. We also found several localized modes caused by defects in the adjacent frequency range (Fig. S2 and S3). Thus, the 131.7 cm−1 peak could be a combined signal of those adjacent localized modes since these defects also exist in as-prepared samples, which calls for a control experiment.

An environment-controlled experiment was carried out to study the effects of degradation of ZrSiTe thin films. We kept the samples in a nitrogen gas-filled glove box, isolated from oxygen and moisture. After one-week storage in the glove box, the Raman spectra of ZrSiTe mono- and few-layers show marked variations. Figure 2c shows the thickness dependence of Raman spectra of the ZrSiTe thin films after one-week storage. For the monolayer, the intensity of the 131.7 cm−1 mode is enhanced by 4.8 times and the full width at half maximum (FWHM) is 2.04 times that of a fresh-cleaved ZrSiTe monolayer, which indicates a higher density of defects in the long-stored monolayer. An appreciable red-shift of this mode was also observed after the one-week storage. Our calculation reveals that several defect-induced localized modes exist around 120 cm−1. Here, we show the vibration displacements of a mode at 122.2 cm−1 in Fig. 4a–c as an example, which is induced by a Te adatom in the monolayer. More defect-induced localized modes are illustrated in Fig. S2c–j. In addition, the phonons between the X and M points can be involved in Raman scattering processes via defect mediation, and activate the rumpling modes. The energy of these phonons is slightly lower than the ones at the X point. Given that the Raman signal around 131.7 cm−1 is a combination of the intralayer shear mode and the defect-induced localized modes, the red-shifts indicate an increasing weight of the localized modes. On the other hand, the \(A_{1g}^1\) and \(B_{1g}^1\) modes are strongly weakened, which most likely result from inversion symmetry breaking in the perpendicular direction because of the defects generated on the topmost Te layer.

Fig. 4
figure 4

Vibrational displacements of two defect-induced Raman modes. a–c Top, front and bottom view of the vibration modes induced by a Te adatom at 122.2 cm−1 and b–f Top view, front view and right view of the vibration modes induced a Te vacancy at 142.9 cm−1, respectively, as shown in a 3 × 3 1L ZrSiTe lattice

Full size image

In the ZrSiTe few-layers, modes\(A_{1g}^2\), \(A_{1g}^1\), and \(B_{1g}^1\) are still dominant. However, the 131.7 cm−1 mode disappears and a rather wide peak consisting of several different modes fades in at around 150 cm−1. All these emerging modes are also results of various kinds of defects. Figure 4d–f illustrates the vibrational displacements of a Te vacancy induced local mode at 142.7 cm−1 and Fig. S3 summarizes other defect-induced Raman-active vibrational modes near 150 cm−1. The absence of the 131.7 cm−1 mode in ZrSiTe few-layers could be, most likely, understood by the scenario that the increasing of interlayer interactions with increasing layer thickness stabilize the ZrSiTe crystalline structure and suppress the layer rumpling motions. And the arising of 143.6 cm−1 modes could be ascribed to the Te-vacancies.

Degradation of the samples was further investigated with X-ray photoelectron spectroscopy (XPS). Detailed experiment results are discussed in Fig. S6. It is clear that roughly one third of the surface Te atoms degraded during a two-week storage in the glove box. It is, thus, reasonable to presume that after one-week storage in N2 atmosphere, a significant number of Te defects appear in the ZrSiTe thin films and the defect-induced mode at 143.6 cm−1 arises afterwards.

In order to confirm these defect-induced vibrational modes, we acquired the secondary differential conductance (d2I/dV2) spectra at defect α and on the pristine surface, as shown in Fig. 5. Such inelastic tunneling spectroscopy (IETS) has been used to detect the surface phonon modes, for instance, on Au(111) surface47. The vibrational spectrum taken on Au(111) shows the symmetric maxima centered at an energy of around ± 9 meV, well within the energy range of phonons on the Au(111) surface47. A giant phonon-mediated inelastic tunneling process (63 ± 2 meV) was also found in graphene by d2I/dV2 spectra48. In our case, a symmetric-placed pair of peak/dip at ± 7.8 meV is also observed on the pristine surface, as shown in Fig. 5b. At the location of defect α, the vibrational energy increases to ±12 meV. According to our DFT calculation, 7.8 (62.9 cm−1) and 12 meV (96.8 cm−1) are corresponding to the energy of the in-plane E2g mode (Fig. S2a) and the Te vacancy-induced intralayer modes (Fig. S4), respectively, which is decided by the selection rules of inelastic tunneling processes. The surface phonon and point-defect-induced phonon modes measured by our IETS agree quite well with the Raman spectra and the DFT calculations, which directly evidence the existence of defect-induced phonon modes. In order to further investigate the structural and symmetry evolution of ZrSiTe few-layers as a function of temperature, we did temperature-dependent Raman measurement and it is clear that no significant crystal phase transition happens during the cooling down despite of consistent stiffening and strengthening of the Raman modes (see Supplementary Information Fig. S5).

Fig. 5
figure 5

The secondary differential conductance spectroscopy of defect α. a STM image for defect α and b the secondary differential conductance (d2I/dV2) spectroscopy taken on the defect α (blue) and the pristine surface (black), respectively. The length of the scale bar in a is 10 nm

Full size image

To confirm the reproducibility of the thickness dependence of ZrSiTe Raman spectra, we did a Raman mapping over a ZrSiTe flake with its thickness varying from monolayer to four layers. Figure 6 shows intensity maps of the representative Raman peaks over the ZrSiTe flake. The monolayer region can be clearly identified where the 131.7 and 143.6 cm−1 signals are significant (bottom-right corner of the sample), whereas the other four modes appear in thicker regions but show nearly no intensity in the monolayer region. The 131.7 cm−1 mode could also be seen in the bilayer region in the bottom-left corner with much lower intensity. These results support our assessment that the 131.7 cm−1 mode is related to the Te adatoms which exhibit a much higher concentration in the ZrSiTe monolayer. We could also identify signals of the 143.6 cm−1 mode at some spots along the edge of the ZrSiTe bilayer, which is consistent with the trend that the edges of ZrSiTe thin films are more suspect to be defective. Moreover, the highest intensity of the two characteristic Raman modes (\(A_{1g}^1\) and \(B_{1g}^1\)) appears on the bilayer region located at the bottom-left corner of the map. In addition, in the center of the bottom of the map, three-layered ZrSiTe shows higher Raman intensity than the thicker surrounding areas, which is consistent with the results that thinner ZrSiTe shows stronger Raman signals.

Fig. 6
figure 6

Intensity mapping of Raman scattering over a ZrSiTe flake. a The optical image of a ZrSiTe flake with various thicknesses spreading over its area. The length of the scale bar in the optical image is 5 μm. b–g Intensity maps of six Raman modes which can be observed in the sample. The monolayer can be clearly identified in the intensity map of the defect related modes at 131.7 and 143.6 cm−1. Generally, the intensity of Raman scattering decreases with the increasing thickness of ZrSiTe from 2L to higher, as shown in the map of the \(A_{1g}^1\) and \(B_{1g}^1\) modes

Full size image

In summary, we have performed a systematic investigation on the Raman-activated vibrational modes in ZrSiTe mono- and few-layers by combining Raman spectroscopy, STM and DFT calculations. We observed three representative Raman modes located at 95 (\(A_{1g}^2\)), 228 (\(A_{1g}^1\)), and 304 cm−1 (\(B_{1g}^1\)) in all thicknesses. Exceptionally, single-atom point defects activate BZ edge modes around 130 cm−1 and induce localized defect modes near 120 or 150 cm−1. The edge modes were assigned to intralayer shear modes near the BZ edge, corresponding to the surface rumpling of the layers, which is most significant in the monolayer but severely suppressed in thicker layers. The localized modes near 120 cm−1 are believed to be a fingerprint of Te adatoms. Their signals are negligible in freshly prepared samples, however, the longer the sample is stored, the stronger their Raman intensity shows and the poorer the vertical symmetry of the layer holds; which results in an appreciable red-shift of the edge modes. All these results evidence that Raman spectroscopy is able to observe some ‘forbidden’ modes in 2D materials, owing to defect-mediated scattering in the q-space or defect-induced local vibrational modes, which allows the detection of various single-atom point defects in metallic layers.

Methods

Sample preparation and Raman measurements

ZrSiTe single crystals were prepared using a chemical vapor transport method. The stoichiometric mixture of Zr, Si and Te powder was sealed in a quartz tube with iodine used as transport agent (2 mg cm−3). Plate-like single crystals with metallic luster can be obtained via the vapor transport growth with a temperature gradient from 950 to 850 °C. The mono- to few-layer ZrSiTe films were achieved by mechanical exfoliation with the Scotch® magic tape in a glove box in N2 atmosphere. During the whole experiment process, the water and oxygen rate was kept below 1 ppm in the glove box. The ZrSiTe thin films were transferred onto the 285 nm SiO2/Si substrate. And a microscope was used to locate the flakes. Then the samples were kept in covered chip holders in the glove box before further measurements. The samples were transferred into aluminum sample boxes with a glass window that are sealed by rubber O-rings before they were taken out of the glove box for Raman spectroscopic measurement. The Raman spectroscopy was performed with a Horiba Jobin Yvon® LabRAM HR800 confocal Raman microscopy. The system is equipped with a Nd:YAG laser (wavelength 532 nm), a 1800 g mm−1 grating, and a Nikon ×50 (NA = 0.5) long-working-distance objective lens. The laser spot size on the sample is close to 1 μm and the spectral resolution is ≈0.55 cm−1. Each spectrum was acquired with an integration time of 8 seconds. Such experimental conditions are identical for all spectral measurements. For polarized Raman measurements, the cross-polarization configuration was used. A half-wave plate was placed in the incident laser path to rotate the polarization of the laser. During the room temperature Raman measurement, the samples were kept in the sealed aluminum sample boxes. And the measurement was performed soon after the location of the ZrSiTe films was identified, normally within 1 h after the films were exfoliated from the bulk crystal. The light went through the glass window and focused on the sample surface to get clear signals from the ZrSiTe thin films. A 532 nm laser was used as the activation light and the laser power was 5 mW during the experiment. For the temperature-dependent Raman measurement, the sample was loaded into the special sample holder used for the measurement as soon as it was taken out from the sealed sample box to eliminate the exposure of the samples to the air. All the samples were stored in a covered sample box in the glove box when they were not used for measurement. A Bruker® Dimension Edge AFM was used to determine the thickness of ZrSiTe thin film. The STM image was captured with a Unisoku STM system. A bulk ZrSiTe crystal with typical 0.5 × 2 mm2 size was attached to the sample holder and the cleavage was done in-situ after the system was pumped to a vacuum better than 1 × 10−10 Torr and cooled down to 4.3 K.

Theoretical methods

DFT calculations were performed using the generalized gradient approximation for the exchange-correlation potential, the projector augmented wave method49,50 and a plane-wave basis set as implemented in the Vienna ab-initio simulation package (VASP)51,52,53, and Quantum Espresso (QE)54. DFPT was employed to calculate phonon-related properties, including Raman intensity (QE), activity (QE) and shifts (VASP), vibrational frequencies at the Gamma point (VASP) and other vibration-related properties (VASP). The kinetic energy cut-off for the plane-wave basis set was set to 700 eV for geometric and vibrational properties for primitive cells from 1 to 4L and 550 eV for the calculations of 3 × 3 1L supercells with defects. A k-mesh of 20 × 20 × 1 was adopted to sample the first Brillouin zone of the conventional unit cell of the few-layers and a 6 × 6 × 1 k-mesh was used for the 3 × 3 × 1 supercell. A Gaussian smearing of 0.05 eV was used in our calculations. In optimizing the system geometry and vibration calculations, van der Waals (vdW) interactions were considered at the vdW-DF55,56 level with the optB86b exchange functional (optB86b-vdW)57,58,59, which was proved to be accurate in describing the structural properties of layered materials27,37,60,61.The shape and volume of each supercell were fully optimized and all atoms in the supercell were allowed to relax until the residual force per atom was less than 1 × 10−3 eV Å−1. The electronic band structure is calculated using Perdew–Burke–Ernzerhof (PBE) functional with spin–orbit coupling (SOC)62.

A 6 × 6 1L supercell was used for the STM simulation with the 400 eV energy cut-off and 2 × 2 × 1 k-point mesh. The defect-deformation energy ΔEi was induced to explain the concentration difference of the three types of defects. The defect-deformation energy can be defined as63

$${{\Delta E}}_{} = {{E}}_{\mathrm{i}}\left( {{\mathrm{defect}}} \right) - {{E}}\left( {{\mathrm{supercell}}} \right) + \mathop {\sum }\limits_{{{s}} = 1}^{N{\mathrm{species}}} {{n}}_{{s}}^{\mathrm{i}}{\mathrm{\mu }}_{{s}}$$

Ei(defect) is the total energy of the system with defect while the E(supercell) is the total energy of the supercell without defect. The \({{n}}_{{s}}^{\mathrm{i}}\) mains the number of the atoms of the s defect: this parameter is negative to adatoms and is positive to vacancies. The μs is the chemical potential of the s defect. Considering it is difficult to obtain the exact value of μs, the system was assumed to vary between two limiting cases: tellurium-rich [\({\mathrm{\mu }}_{{\mathrm{Te}}} = {\mathrm{\mu }}_{{\mathrm{Te}}}^ \ast\) and \({\mathrm{\mu }}_{{\mathrm{Zr}}}^ \ast + {{\Delta H}} \,<\, {\mathrm{\mu }}_{{\mathrm{Zr}}} \,<\, {\mathrm{\mu }}_{{\mathrm{Zr}}}^ \ast\), where \({\mathrm{\mu }}_{{\mathrm{Te}}}^ \ast\) is the chemical potential of Te in the α-phase crystal form, \({\mathrm{\mu }}_{{\mathrm{Zr}}}^ \ast\) is the chemical potential of Zr in the bulk form and ΔH is the formation enthalpy of ZrSiTe] and zirconium-rich (tellurium-deficient) \(\left[ {{\mathrm{\mu }}_{{\mathrm{Zr}}} = {\mathrm{\mu }}_{{\mathrm{Zr}}}^ \ast \,{\mathrm{and}}\,{\mathrm{\mu }}_{{\mathrm{Te}}}^ \ast + {{\Delta H}} \,<\, {\mathrm{\mu }}_{{\mathrm{Te}}} \,<\, {\mathrm{\mu }}_{{\mathrm{Te}}}^ \ast } \right]\). The 3 × 3 bilayer supercell was used to calculate the defect-deformation energy.

Vibrational frequencies and polarization vectors at the Brillouin zone center were calculated using the DFPT. In a rigid-layer vibrational mode, the whole layer can be treated as one rigid body64,65. The projected interlayer force constant (p-ILFC) ki (i stands for the projected direction, e.g., x, y, or z) was constructed by summing inter-atomic force constants over all atoms from each of the two adjacent layers, as \({{k}}_{{\mathrm{AB}}}^{{i}} = \mathop {\sum }\limits_{{{a}},{{b}}} {{D}}_{{{a}},{{b}}}^{{i}}\left( {{{a}} \in \left[ {{\mathrm{atoms}}\,{\mathrm{in}}\,{\mathrm{layer}}\,{\mathrm{A}}} \right],{{b}} \in \left[ {{\mathrm{atoms}}\,{\mathrm{in}}\,{\mathrm{layer}}\,{\mathrm{B}}} \right]} \right)\). The matrix of inter-atomic force constants, essentially the same as the Hessian matrix of Born–Oppenheimer energy surface, is defined as the energetic response to a distortion of atomic geometry in DFPT66. It reads, \({{D}}_{{{a}},{{b}}}^{{i}} = \frac{{\partial ^2{{E}}({{R}})}}{{\partial {{R}}_{{a}}^{{i}}\partial {{R}}_{{b}}^{{i}}}}\), where R is the coordinate of each ion, E(R) is the ground-state energy. The force on an individual ion can thus be expressed to \(F_i = - \frac{{\partial E({\mathbf{R}})}}{{\partial {\mathbf{R}}_i}}\). Vibrational frequency ω are related to the eigenvalues of the Hessian matrix and atomic masses. It reads \({\mathrm{det}}\left| {\frac{1}{{4{\mathrm{\pi }}^2{{c}}^2}}\frac{{{{D}}_{{{a}},{{b}}}}}{{{{m}}_{{{a}},{{b}}}}} - {\mathrm{\omega }}^2} \right| = 0\), where ma,b is the effective mass, defined as \(\sqrt {{{m}}_{{a}}{{m}}_{{b}}}\).