INTRODUCTION

Quasiparticle interference (QPI) has in recent years developed into a powerful tool for the characterization of electronic structure. Due to its unrivalled energy resolution, it has been highly successful in unravelling the detailed dispersion relation and structure of the superconducting order parameter in many correlated electron materials1,2,3. It is, however, inherently limited by imaging of the electronic states in two dimensions, and has consequently predominantly been applied to materials where the electronic structure can be approximated as two-dimensional, either due to a strong anisotropy, e.g., in layered materials1,2, or because the states are purely two-dimensional, as is the case for surface states4,5.

In a material with a crystal structure with high symmetry, for example, cubic, the interpretation of QPI will contain contributions from the full three-dimensional electronic structure, and scattering from defects that are below the surface layer will become non-negligible6. This makes the interpretation of the QPI patterns challenging as illustrated in Fig. 1a: a 2D measurement of the density of states at the surface contains information from a three-dimensional band structure and distribution of defects.

Fig. 1: Crystal structure and Brillouin zone of PbS.
figure 1

a Visualization of quasiparticle interference in an isotropic material. A typical STM topography (obtained from PbS), exhibiting an atomically flat surface with a wide range of defects, is shown (Vset = 0.8 V, Iset = 0.2 nA). The (schematic) differential conductance map, which is proportional to the density of states, shows clear signatures of quasiparticle interference (QPI), evidenced as oscillatory patterns around defects in the surface, as well as sub-surface layers. These patterns contain information from the full electronic structure. b Difference between the primitive cell of the bulk and surface. Left: the primitive cell of the bulk FCC crystal structure shown in a unit cell of PbS. Right: minimal supercell with lattice parameters a and \(c=\sqrt{2}a\) required to describe the surface. All crystallographic directions in the following refer to this supercell. c Schematic bulk Brillouin zone (BZ) of a cubic material with Fermi pockets centred at the zone boundaries, corresponding to the primitive cell shown on the left in b with a schematic electronic structure of an FCC material. The constant energy surface close to the top of the valence band in PbS is represented by pockets centred at the L-point (centre of the hexagonal faces). There are no states in the kz = 0 plane. d Brillouin zone corresponding to the minimal supercell shown on the right in (b), resulting in folding of the Brillouin zone and hence of the Fermi surface. e Cuts at different kz planes for the Fermi surface in (c). With decreasing kz the size of the contours increases. There is one dominant intra-pocket scattering vector q. f Decreasing kz planes for the Brillouin zone shown in (d). Due to the folding, there are two new intra-pocket scattering vectors, \({{{{{{{\bf{q}}}}}}}}^{\prime}\) and q″, that do not exist in the bulk Brillouin zone.

Full size image

To describe how the quasiparticle interference measured at the surface layer relates to the bulk electronic structure requires addressing what the correct mapping between surface and bulk reciprocal spaces for its description is. Figure 1b illustrates this point for a rocksalt crystal structure: the primitive cell in the bulk cannot be used to describe a (001) surface. The minimal unit cell required to describe the surface as well as its electronic structure is thus larger. As a consequence, the surface BZ in which the electronic structure is measured by QPI is not just a cut or projection of the bulk BZ, but a folded version of it.

In Fig. 1c, we show a Fermi surface in the bulk Brillouin zone and in the minimal Brillouin zone required to describe the surface (Fig. 1d). In many face-centred cubic systems, the Fermi surface has pockets at the zone boundary of the bulk BZ, as sketched in Fig. 1c. Considering a constant kz plane, the folding in the surface layer maps out-of-plane scattering vectors between hot spots on the Fermi surface into the same kz-plane (Fig. 1e, f).

A material with these properties is the mineral galena (PbS). It naturally crystallizes in a rock-salt crystal structure and exhibits a perfect cleaving plane along the {100}-directions, thus providing a 3D band structure with cubic symmetry and atomically flat surfaces. It is a narrow direct bandgap semiconductor, whose properties differ from those of more conventional semiconductors (e.g., GaAs), where the bandgap occurs at the Γ-point in the centre of the Brillouin zone. In galena, the bandgap is at the L-point at the boundary of the first Brillouin zone (compare Fig. 2a) and its width decreases with decreasing temperature7,8 as well as under pressure9. The electronic structure is well known and features particle−hole symmetry across its semiconducting gap10. It undergoes topological changes as a function of energy when moving away from the valence band maximum (VBM) and from the conduction band minimum (CBm) with increasing absolute energy. At the VBM and CBm, there are no states in the kz = 0 plane but only at the edges of the folded BZ (analogous to Fig. 1d), making this material an ideal choice to establish the contribution from states with non-zero kz to the quasiparticle interference. It is only at about 0.6 eV below (0.8 eV above) the VBM (CBm) that states appear in the kz = 0 plane (compare Fig. 2a). From angular-resolved photoemission, there is no evidence for a surface state in the vicinity of the Fermi energy (Fig. 2c, also ref. 8), and so the detected QPI signal originates solely from the bulk electronic structure.

Fig. 2: Band structure and density of states.
figure 2

a Band structure of PbS from a tight-binding model for the bulk FCC BZ is shown on the left. The valence band is made up of sulfur p-orbitals, while the conduction band is mainly due to lead p-orbitals. The direct bandgap occurs at the L-point. Van-Hove singularities are indicated by M1 and M2 for the valence band and m1 and m2 for the conduction band. The panel on the right shows the total density of states. b Tunnelling spectroscopy, \({{{{{\rm{d}}}}}}\,{{{{{{\mathrm{ln}}}}}}}\,I/{{{{{\rm{d}}}}}}\,{{{{{{\mathrm{ln}}}}}}}\,V\), of galena, showing the gap Eg between conduction and valence band. A gap of 0.34 V is observed, with the bottom of the conduction band pinned to the Fermi level. Four peaks can be observed, which can be identified with the van Hove singularities in the band structure. Because of the strong increase in density of states at positive bias voltages, spectra have been recorded in three energy ranges indicated by dashed lines, which we probe with different setpoint conditions (Black: Vset = 1.2 V, Iset = 0.1 nA; dark grey: Vset = 0.6 V, Iset = 0.2 nA; light grey: Vset = 0.8 V, Iset = 0.2 nA). c Electronic band structure along the \(\overline{{{\Gamma }}}-\overline{{{{{{{{\rm{M}}}}}}}}}\) direction of the surface-projected Brillouin zone as measured by ARPES (left) and calculated by DFT (right), showing excellent agreement. The DFT calculation simulates the finite kz-resolution of the ARPES experiment using empirical parameters (see Supplementary Note 8). d Topographic images showing characteristic defects at lead (D1) and sulfur (D2) sites for positive and negative bias (top: Vset = 1 V, Iset = 150 pA; bottom: Vset = −1 V, Iset = 100 pA). e Spatial map of dlnI/dlnV at 0.85 V and f at −0.85 V. The circles identify the same types of defects as in (d).

Full size image

The interpretation of the quasiparticle interference detected at the surface of a material without a clear 2D anisotropy thus immediately raises several fundamental questions: (1) what is the effect of the three-dimensionality of the electronic structure on the quasiparticle interference detected at the surface? (2) How does the zone-folding from the bulk Brillouin zone (BZ) to the minimal cell required to describe the surface affect the scattering processes contributing to the QPI? and (3), how do sub-surface impurities contribute to the scattering signal?

Here, we address these questions through imaging of quasiparticle interference at the surface of the mineral galena (PbS), which has a bulk electronic structure with cubic symmetry. Comparison with theory reveals that in galena, the QPI signal observed at any given bias voltage is dominated by hotspots in certain kz planes of the bulk BZ, demonstrating 3D quasiparticle tomography of the electronic structure.

RESULTS

The galena sample was cleaved in-situ at cryogenic temperatures, producing a high-quality atomically-flat surface. A representative topographic STM image, measured at 20 K, shows atomic resolution, exhibiting the square atomic lattice of the lead atoms (Fig. 1a and Supplementary Fig. 1a). The surface exhibits a rich variety of defects found at both lead and sulfur lattice sites as a consequence of its natural origin and aiding in the observation of quasiparticle interference. In excellent agreement with the calculated total density of states (Fig. 2a) and with the previous reports8,11,12, typical tunnelling spectra (Fig. 2b) reveal the energy gap of galena on the order of Eg ≈ 0.34 eV. We find that the CBm is pinned to the Fermi level, indicating n-doping of the sample. Additionally, two sets of peaks are observed in tunnelling spectra: M1 and M2 at negative energies and m1 and m2 at positive energies, associated with van Hove singularities in the band structure related to band minima (m) and maxima (M) (compare Fig. 2a) in agreement with the known density of states13,14. Comparison of the calculated band structure with ARPES measurements of the electronic structure show excellent agreement between the two (Fig. 2c).

In the rock salt crystal structure of PbS, defects at sulfur and lead sites are expected to lead to significantly different scattering patterns. Figure 2d shows characteristic defects on the lead and sulfur sites imaged at positive and negative bias voltage with a distinctly different appearance. From the chemical analysis of the sample (see Supplementary Note 1), the predominant defects are Bi (n-doping) and, with a smaller concentration, Ag (p-doping), leading to a net n-doping of the sample. Spatially-resolved maps of \({{{{{\rm{d}}}}}}\,{{{{{{\mathrm{ln}}}}}}}\,I/{{{{{\rm{d}}}}}}\,{{{{{{\mathrm{ln}}}}}}}\,V({{{{{{{\bf{r}}}}}}}},{V}_{b})\) at Vb = 0.85 V and −0.85 V are shown in Fig. 2e, f, respectively. The same type of defects as in Fig. 2d are highlighted on both maps. At both energies, clear QPI patterns are visible as spatial modulations around all defects, which appear in the Fourier transformation (Fig. 3a, e) as characteristic wave vectors along the [100] and [110] directions. The energies shown in Fig. 3a, e are close to the VBM and the CBm, respectively, where the constant energy surfaces have no states in the kz = 0 plane, yet we see strong signatures of quasiparticle scattering. Comparison with the bulk electronic structure reveals that the scattering we observe along the [100] direction can only originate from scattering between disconnected hole (electron)-pockets around the L-point at the boundary of the bulk Brillouin zone, as represented in Fig. 1c, d. The same pockets account for the scattering vectors visible along the [110] direction. Moving further away from the gap edge, these pockets become connected, as observed in Fig. 3g. The Fourier transformations show closed patterns whose characteristic wave vector decreases with increasing absolute energy.

Fig. 3: Comparison of experimental QPI patterns with T-matrix calculations.
figure 3

a−d QPI patterns at negative energy with decreasing energy. The panels on the left show the comparison between the experimental data and the T-matrix calculation of the QPI for the full bulk electronic structure. For each energy, the ratio pz/pxy of the orbital contributions required to describe the experimental data is indicated. The panels on the right show the QPI patterns originating from sulfur {px, py} (left) and pz states (right). e−h QPI patterns at positive energies with increasing energy with the orbital contributions from the Pb p-orbitals. Left and right panels as for (a−d). Near the band edge, disconnected QPI patterns are observed as expected for the top and bottom of the valence and conduction bands, respectively, which join as the absolute energy increases. For all panels, the q vector spans the reciprocal space from \((\frac{-2\pi }{a},\frac{-2\pi }{a})\) to \((\frac{2\pi }{a},\frac{2\pi }{a})\), where a is the lateral lattice parameter of the tetragonal supercell shown in Fig. 1b.

Full size image

In order to understand the origin of the features in QPI, we have developed a formalism based on T-matrix calculations to describe the QPI patterns from the full three-dimensional electronic structure. We use a tight-binding model for the bulk obtained from DFT calculations to obtain the bare Green’s function G0(k, ω) of the unperturbed system (see Methods section), and then evaluate the q-resolved density of states ρ(q, ω) from

$$\rho ({{{{{{{\bf{q}}}}}}}},\omega )=-\frac{1}{2\pi i}\mathop{\sum}\limits_{{{{{{{{\bf{k}}}}}}}}}\,{{\mbox{Tr}}}\,\{G({{{{{{{\bf{k}}}}}}}},{{{{{{{\bf{k}}}}}}}}-{{{{{{{\bf{q}}}}}}}},\omega )-{G}^{* }({{{{{{{\bf{k}}}}}}}},{{{{{{{\bf{k}}}}}}}}+{{{{{{{\bf{q}}}}}}}},\omega )\},$$
(1)

where \(G({{{{{{{\bf{k}}}}}}}},{{{{{{{{\bf{k}}}}}}}}}^{\prime},\omega )\) is the retarded Green’s function of the system including the impurities (see Methods section). Note that here, unlike in the QPI calculations done usually, the wave vectors k and \({{{{{{{{\bf{k}}}}}}}}}^{\prime}\) are three-dimensional wave vectors. We account for scattering processes through the T-matrix which encodes the details of the scattering potential \({V}_{{{{{{{{\bf{k}}}}}}}},{{{{{{{{\bf{k}}}}}}}}}^{\prime}}\) (see Methods section, eq. (4)).

To obtain the surface-projected QPI, \(\tilde{\rho }({q}_{x},{q}_{y},\omega )\), we first divide the folded bulk BZ into a series of (qx, qy) planes parallel to the surface, and for each plane calculate the bulk QPI \(\rho ({q}_{x},{q}_{y},\omega ){| }_{{k}_{z}}\), using the Green’s function defined in Eq. (1). While we calculate the QPI in planes with qz = 0, scattering vectors with non-zero qz are accounted for implicitly through the folding. We then integrate over all kz components to obtain the QPI, \(\tilde{\rho }({q}_{x},{q}_{y},\omega )\). The resulting QPI pattern obtained from a calculation describing the kz dependence and the orbital character of the wave functions is shown in Fig. 3 next to the experimental data. To fully describe the experimental data, we manually adjust the ratio of the contributions from px,y, and pz orbitals. The extracted ratio pz/px,y is shown in Supplementary Fig. 8. This change in orbital contribution with bias voltage reflects the diversity of scattering processes from different defects in the bulk and their directional anisotropy and energy dependence. The comparison reveals excellent and nearly quantitative agreement.

Our analysis reveals several important points: (1) For both positive and negative bias voltages, pz- and {px, py}-orbitals create distinct QPI patterns, despite the orbital degeneracy of bulk electronic bands. (2) As a function of bias voltage, we observe a change in the relative intensity with which bands with pz- and {px, py}-orbital character appear in the observed QPI: the ratio is found to continuously increase with increasing energy both above and below the Fermi energy. (3) The experimental QPI pattern can only be reproduced when accounting for the full kz-dependence of the electronic structure. A quasi-2D model would fail spectacularly for the present system.

DISCUSSION

Our framework to describe the QPI allows us to provide an understanding of these observations and disentangle the contributions from different kz planes to make the connection from the QPI patterns to the bulk electronic structure.

To better understand the energy dependence of the relative contributions of bands with different orbital character to the observed QPI patterns, a clear understanding of symmetry and dimensionality of the electronic bands and their interaction with the existing impurities is necessary. Bulk PbS has an isometric crystal structure with the highest possible symmetry. As such, the p-orbital manifold is subject to an isotropic cubic crystal field, protecting the orbital degeneracy of the energy bands. At the surface, this is no longer the case, due to the surface discontinuity, suppressing the electron hopping along the surface normal. The (001) surface termination creates a tetragonal crystal field near the surface, causing an energy splitting between the pz and {px, py} orbitals. As the latter are spatially distributed along [100] and [010] directions, they contribute predominantly to the QPI features at the zone centre or along the corresponding in-plane surface reciprocal wave vectors, [10] and [01] (see Fig. 3). On the other hand, the pz orbital, due to its apical distribution, can isotropically contribute to scattering in any in-plane direction. As a result, it becomes the main channel for impurity scattering within and near the [110] plane, appearing as diagonal features in the observed and calculated QPI. Furthermore, both Pb-pz and S-pz can strongly interact with the impurities beneath the surface, as they are the only orbitals that can form strong interlayer σ-type bonding; {px, py} can only form π-type bonding between the PbS layers, which is much weaker than σ bonding, leading to weaker scattering patterns.

Such a characteristic difference between these two orbital sets and their contribution to the QPI reveals a fundamental distinction in the orbital character of the underlying scattering processes. To further explain how they respond to the bias variation, we also need to consider the energetic alignment of the impurity states and their relative coupling to the host energy bands.

As shown in Fig. 3, at negative bias voltages close to the Fermi energy, both pz and {px, py} contribute almost equally to the observed QPI, but the latter becomes increasingly dominant at more negative Vb. This means that scattering in the {px, py} channel at surface impurities becomes dominant at high negative bias voltages. Given that the existing impurities are mainly charge-donors (cationic), this finding is not unexpected. Cations tend to have impurity states at higher energies than anions. At negative bias voltages, we are probing anionic S-derived bands, further away from the impurity states, lowering the probability of scattering in total and from defects in the sub-surface layers in particular. Thus the resulting QPI tends to become more dominated by the surface {px, py} orbitals as Vb becomes more negative. For a positive Vb, the situation is the opposite. Since now both the impurity states and Pb-derived bands are cationic, near the band edge, we see a more significant contribution from {px, py} orbitals. This can be interpreted as a higher probability of scattering from the top surface layer as the Pb-site impurities in this layer couple more strongly to the neighbouring Pb ions through {px, py} orbitals. As Vb increases, the scattering from defects in subsurface layers through pz orbitals appears to gain an increasing strength, leading to a dramatic enhancement of the ratio pz/px,y. Such a dimensional anisotropy in quasiparticle scattering and its dependence on the chemical character of impurities underlines the importance of our bulk-based analysis of the QPI patterns observed at the surface.

The relation between the kz planes dominating the QPI signal and the electronic structure is illustrated in Fig. 4. Figure 4a shows a cut in the (010) plane through constant energy surfaces of the electronic structure at increasing absolute energies E1…E4. It is obvious that the quasiparticle interference patterns near the band edge cannot originate from the kz = 0 plane as there are no states—demonstrating the need for a description of the QPI that accounts for the full 3D band structure. For surface defects, which will make the dominant contribution, the scattering vectors that contribute most connect pockets of the band structure with group velocities in the surface plane, indicated by the red dashed line in Fig. 4a.15 To identify the scattering vectors dominating the quasiparticle interference patterns and efficiently calculate the QPI signal, we evaluate the flatness of the combined pockets in the electronic structure in the minimal cells required for the surface (folded) and the bulk (primitive) Brillouin zone, Fig. 4b. In the folded zone, the points with in-plane group velocity map into a single kz-plane, and hence the dominant scattering vectors have a z-component which is zero, making the calculation of the QPI signal more efficient. To determine the dominant kz-plane, we define the flatness as the inverse of the derivative of the area with respect to kz, i.e., \({({{{{{\rm{d}}}}}}A/{{{{{\rm{d}}}}}}{k}_{z})}^{-1}\). While using the bulk Brillouin zone shows a dominant kz vector at the zone boundary, in the folded Brillouin zone, corresponding to the minimal surface unit cell, the dominant kz plane is away from the high-symmetry planes—consistent with the experiment. The dominant kz-plane extracted in this way is shown in Fig. 4c, with the resulting scattering vectors. We note that this dominant kz plane is energy-dependent, as the band pockets grow, the dominant kz plane moves towards the zone centre (see Fig. 4a). The flatness curve in Fig. 4b shows additional maxima at its edges which are of no significance in QPI calculations, as the Fermi surfaces at these kz cuts have a negligible cross-sectional area.

Fig. 4: kz selectivity of quasiparticle interference.
figure 4

a Representation of (010)-planes at ky = 0 with increasing absolute energy. The first frame depicts an energy close to the VBM. Lines indicate constant energy contours of the bulk electronic structure, orange lines in the bulk Brillouin zone, and black lines of folded bands in the surface Brillouin zone. The red dashed lines indicate the kz plane where the bands exhibit an in-plane group velocity. The folding maps the scattering vectors connecting hot spots in different kz planes into the same kz plane. With increasing absolute energy (indicated by the arrow), the dominant kz-plane for scattering moves towards kz = 0. b To determine the dominant kz plane, we use the surface area A (left panel, showing the merged area enclosed by the union of the two pockets for the folded zone) to determine the flatness ∣dA/dkz∣−1 (right panel) as a function of kz for the pockets in (a) (shown here for E = −0.85 eV). Once folding is accounted for, the points with the largest flatness occur at kz = 0.88π/c, indicated by the dashed red line. c Cut along the red line in (a) for the first layer. The vectors q1 and q2 connect unfolded Fermi surfaces. d Maps of the q-resolved density of states \(\tilde{\rho }({{{{{{{\bf{q}}}}}}}},E)\) due to QPI, taking into account the energy dependence of the orbital contributions, as a function of energy E for different kz planes. e Corresponding experimental data along the same path as in (d). The circles indicate the position of the peaks (error bars represent q-space resolution of the dI/dV-maps). f QPI layer at Vb = −0.85 V, with the arrows as in (c).

Full size image

The quasiparticle dispersion obtained from the calculation for different kz planes is shown in Fig. 4d together with the experimentally observed one (Fig. 4e). The dispersion at energies close to the VBM and CBm is visible only at kz planes close to the zone boundary, near kz = π/c, whereas planes close to kz = 0 become more relevant for higher absolute energies. This becomes apparent from comparing the scattering vectors obtained from the experiment with the calculation (white dotted lines in Fig. 4d), showing agreement for some scattering vectors for \({k}_{z}=\frac{\pi }{c}\) (q1 and q2), whereas other vectors clearly do not fit (\({q}_{2}^{\prime}\)) but are better described by scattering from a different kz-plane (see Supplementary Fig. 5 for detailed comparison). The dominant scattering vectors in the kz plane with the largest contribution at −0.85 V are shown in Fig. 4f superimposed to the experimental quasi-particle interference.

Our results highlight how quasiparticle interference imaging can be used to probe the three-dimensional electronic structure. As a basis, we have chosen a system that has a cubic symmetry and a three-dimensional electronic structure, and does not exhibit electron correlations, which enables for a complete and accurate theoretical description of the electronic structure and the wave functions, making this an ideal test system. The high precision of the band structure calculation is seen in the comparison with the ARPES data.

Different to quasi-two-dimensional materials, where it is often sufficient to only consider the kz = 0 plane, in this system, the dispersion in the direction normal to the surface cannot be ignored: the QPI patterns observed in STM cannot be simply explained by the quasiparticle signal coming from a particular kz-plane, or even just the kz = 0 plane. What we rather find is that the QPI signal is dominated by bias-dependent kz planes, determined by the plane of dominant in-plane group velocity. Taking into account the surface Brillouin zone, we can understand the QPI signal once the correct kz-plane is determined, resulting in excellent agreement between theory and experiment of both the dominant QPI features as well as their dispersion.

The picture that emerges has far-reaching consequences for the interpretation of quasiparticle interference in materials where the electronic structure is not quasi-2D. In these cases, the dominant scattering plane, which can be multiple planes in multi-band systems, needs to be determined to be able to extract information about the electronic structure from the quasiparticle interference patterns. We note that in many anisotropic systems, the quasiparticle interference will be dominated by scattering vectors in high-symmetry planes, where the group velocity has to be in-plane by symmetry, for example k-vectors in the kz = 0 plane and at the zone face16. The present example, however, shows that this is not generally true, and the dominant kz plane can differ from high symmetry planes and even become energy-dependent. Correct identification of the dominant kz-planes and their energy dependence enables tomographic characterization of the bulk electronic structure from quasiparticle interference imaging.

In conclusion, we have demonstrated quasiparticle tomography of the bulk electronic structure from a surface-sensitive measurement. Modelling the scattering vectors due to the three-dimensional electronic structure using the zone folding that occurs in the surface region, we find that the dominant contribution to the quasiparticle interference imaging comes from an energy-dependent kz-plane. We introduce a theoretical framework to efficiently compute the resulting QPI patterns. By accounting for the full 3D electronic structure and the orbital composition of the surface electronic bands, we are able to fully reconstruct the experimentally observed quasiparticle interference patterns from theory. Our results introduce a method to extract information about the 3D bulk electronic structure from a quasi-particle interference measurement in 2D.

METHODS

Sample details

The sample used in this study is a mineral sample of galena collected in Malawi. Galena is the natural ore mined for lead and silver. Our sample is an argentiferous galena sample, i.e., has significant silver content. We have analyzed the sample composition by Inductively Coupled Plasma Mass Spectrometry (ICPMS), for details of the analysis see Supplementary Note 1. Galena has a cubic crystal symmetry with well-defined cleavage planes in the [001] directions, exposing mirror-like surfaces.

STM measurements

We performed STM measurements in a home-built ultra-low temperature STM, mounted in a dilution refrigerator17, working in a cryogenic vacuum. The STM tip was made of Pt−Ir wire, prepared prior to measurements on a single crystal of Au(111) by field emission. The PbS sample was cleaved in situ at low temperatures. All measurements were performed at 20 K to thermally excite a sufficiently large number of charge carriers so that STM imaging becomes possible. The temperature was stabilized by a heater on the STM head. The bias voltage was applied to the sample, with the tip at virtual ground. Tunnelling spectroscopy was performed using a lock-in amplifier to measure differential conductance, dI/dV, where the bias voltage V was modulated at a frequency of ν = 437Hz. Due to the increase in differential conductance at large bias voltages, different regions of the measured bias voltages had to be probed with different bias and current setpoints, Vset and Iset, so that both the semiconducting gap and the features at negative bias could be resolved. In this way, the measured bias was split into three intervals: [−1.8, −0.6], [−0.6, 0.6], and [0.6, 2.0], in volts, as shown in Fig. 2b. The spectroscopic and QPI data shown in this work correspond to18

$$\frac{{{{{{{{\rm{d}}}}}}}}\,{{{{{{\mathrm{ln}}}}}}}\,I}{{{{{{{{\rm{d}}}}}}}}\,{{{{{{\mathrm{ln}}}}}}}\,V}(V)=\frac{({{{{{{{\rm{d}}}}}}}}I/{{{{{{{\rm{d}}}}}}}}V)(V)}{I(V)/V},$$
(2)

to allow consistent comparison between measurements taken with different bias and current setpoints.

Angle-resolved photoemission spectroscopy

We performed angle-resolved photoemission spectroscopy (ARPES) on the same sample, at the APE beamline of the Elettra synchrotron in Trieste, Italy19. The samples were cleaved in situ and measured at 78 K, using p-polarized synchrotron light with the photon energy of 70 eV, and the data were collected using a Scienta-Omicron DA30 electron energy analyzer.

Calculations

The QPI in momentum space, i.e., Eq. (1), was computed using the retarded Green’s function G defined as

$$G({{{{{{{\bf{k}}}}}}}},{{{{{{{\bf{{k}}}}}}}^{\prime}}},\omega )={G}_{0}({{{{{{{\bf{k}}}}}}}},\omega ){\delta }_{{{{{{{{\bf{k}}}}}}}},{{{{{{{\bf{{k}}}}}}}^{\prime}}}}+{G}_{0}({{{{{{{\bf{k}}}}}}}},\omega ){T}_{{{{{{{{\bf{k}}}}}}}},{{{{{{{\bf{{k}}}}}}}^{\prime}}}}(\omega ){G}_{0}({{{{{{{\bf{{k}}}}}}}^{\prime}}},\omega ),$$
(3)

with T being a matrix encompassing all scattering processes between the wave-vectors k and \({{{{{{{\bf{k}}}}}}}}^{\prime} ={{{{{{{\bf{k}}}}}}}}-{{{{{{{\bf{q}}}}}}}}\) by the impurity potential \({V}_{{{{{{{{\bf{k}}}}}}}},{{{{{{{\bf{{k}}}}}}}^{\prime}}}}\),

$${T}_{{{{{{{{\bf{k}}}}}}}},{{{{{{{\bf{{k}}}}}}}^{\prime}}}}(\omega )={V}_{{{{{{{{\bf{k}}}}}}}},{{{{{{{\bf{{k}}}}}}}^{\prime}}}}+\mathop{\sum}\limits_{{{{{{{{\bf{k}}}}}}}}^{\prime\prime} }{V}_{{{{{{{{\bf{k}}}}}}}},{{{{{{{\bf{k}}}}}}}}^{\prime\prime} }{G}_{0}({{{{{{{\bf{k}}}}}}}}^{\prime\prime} ,\omega ){T}_{{{{{{{{\bf{k}}}}}}}}^{\prime\prime} ,{{{{{{{\bf{{k}}}}}}}^{\prime}}}}(\omega ),$$
(4)

and G0(k, ω) being the bare Green’s function

$${G}_{0}({{{{{{{\bf{k}}}}}}}},\omega )=\mathop{\sum}\limits_{n}\frac{\left|n,{{{{{{{\bf{k}}}}}}}}\right\rangle \left\langle n,{{{{{{{\bf{k}}}}}}}}\right|}{\omega +{{{{{\rm{i}}}}}}\eta -{\epsilon }_{n}({{{{{{{\bf{k}}}}}}}})},$$
(5)

resulting from the Bloch eigenvalues ϵn(k) and eigenvectors \(\left|n,{{{{{{{\bf{k}}}}}}}}\right\rangle\) with the broadening factor η.

Since the impurities in our PbS samples are expected to be non-magnetic and randomly distributed in the whole system, we neglect the spin- and momentum-dependence of the scattering potential and regarded it as \({V}_{{{{{{{{\bf{k}}}}}}}},{{{{{{{\bf{{k}}}}}}}^{\prime}}}}={V}_{0}{\mathbb{1}}\), where \({\mathbb{1}}\) is the identity matrix and V0 is a constant potential, here fixed at V0 = 0.1 eV added to the on-site energy at the defect site. Similarly, the resulting T matrix was assumed to be only a function of ω and independent of k. We then used the convolution technique developed by Kohsaka et al.20 to compute G(k, k − q, ω) over a fine 600 × 600 q-mesh, spanning the entire surface BZ. To further accelerate the calculations and, more importantly, to account for the vacuum overlap, the Bloch eigenvalues ϵn(k) and eigenvectors \(\left|n,{{{{{{{\bf{k}}}}}}}}\right\rangle\) were extracted from a low-energy tight-binding (TB) model, via a set of atomic-orbital (AO)-like Wannier functions21. We used the Fourier representation of these Wannier functions to transform all the Green’s functions from their periodic lattice model to a continuum model22,23,24. The resulting QPI was subsequently projected onto a plane with a fixed height z = 5 Å, corresponding to the distance between the STM tip and the surface (in our model this corresponds to a plane 5 Å above the topmost PbS layer in the minimal bulk supercell).

To construct the TB model, we first performed a fully-relativistic DFT calculation using Perdew−Burke−Ernzerhof exchange-correlation functional as implemented in WIEN2K25, for the tetragonal bulk supercell containing two formula units of PbS shown in Fig. 1b. The corresponding BZ was sampled by 20 × 20 × 16 k-points. To accurately account for the correlation effects, an additional non-local potential term was added to the DFT Kohn-Sham Hamiltonian using the modifed Becke−Johnson scheme26. From this, we then downfolded a 12 × 12 TB model using twelve AO-like Wannier functions21, representing the top six valence bands, dominated by S-2p orbitals, and bottom six conduction bands, mainly of Pb-6p character. The resulting TB Hamiltonian,

$$H=\mathop{\sum}\limits_{{{{{{{{\bf{R}}}}}}}},{{{{{{{\bf{R}}}}}}}}^{\prime} }{t}_{{{{{{{{\bf{R}}}}}}}},{{{{{{{\bf{R}}}}}}}}^{\prime} }^{\mu \nu ,\sigma \sigma ^{\prime} }{c}_{{{{{{{{\bf{R}}}}}}}}\mu \sigma }^{{{{\dagger}}} }{c}_{{{{{{{{\bf{R}}}}}}}}^{\prime} \nu \sigma ^{\prime} }$$
(6)

describes the creation of an electron by the operator \({c}_{{{{{{{{\bf{R}}}}}}}}\mu \sigma }^{{{{\dagger}}} }\) at the lattice site R with the respective AO and spin characters μ and σ, after the annihilation of another electron with AO and spin characters ν and \(\sigma ^{\prime}\) by the operator \({c}_{{{{{{{{\bf{R}}}}}}}}^{\prime} \nu \sigma ^{\prime} }\) at site \({{{{{{{\bf{R}}}}}}}}^{\prime}\). The corresponding hopping parameter is denoted by \({t}_{{{{{{{{\bf{R}}}}}}}},{{{{{{{\bf{R}}}}}}}}^{\prime} }^{\mu \nu ,\sigma \sigma ^{\prime} }\). Here, μ and ν can be any of the outermost {px, py, pz} orbitals of Pb and S.

By aligning the c-axis of the TB supercell along the crystalline [001] direction, we can ensure that all surface projections occur within a square superlattice exactly as found in the STM experiment. This results in a double folding of the wave functions in k-space along the kz direction, as shown in Fig. 1c, d. Correspondingly, any kz cut of the Fermi surface contains the energy contours from the primitive cell as well as their folded counterparts. An important outcome of such folding is that it enables us to selectively determine all major inter- and intra-pocket scattering vectors q without knowledge of their qz component (see Fig. 1e, f). One can identify which kz component of bulk wave functions contribute the most to the QPI pattern observed at the surface and how it evolves by varying the bias voltage. Also, since in our methodology we use AO-like Wannier functions, we can easily decompose the calculated QPI in terms of different orbital components. As such, it is possible to specify how strongly each orbital contributes to the observed QPI.