Quasiparticle Interference of the van-Hove singularity in Sr$_2$RuO$_4$

The single-layered ruthenate Sr$_2$RuO$_4$ is one of the most enigmatic unconventional superconductors. While for many years it was thought to be the best candidate for a chiral $p$-wave superconducting ground state, desirable for topological quantum computations, recent experiments suggest a singlet state, ruling out the original $p$-wave scenario. The superconductivity as well as the properties of the multi-layered compounds of the ruthenate perovskites are strongly influenced by a van Hove singularity in proximity of the Fermi energy. Tiny structural distortions move the van Hove singularity across the Fermi energy with dramatic consequences for the physical properties. Here, we determine the electronic structure of the van Hove singularity in the surface layer of Sr$_2$RuO$_4$ by quasiparticle interference imaging. We trace its dispersion and demonstrate from a model calculation accounting for the full vacuum overlap of the wave functions that its detection is facilitated through the octahedral rotations in the surface layer.


I. INTRODUCTION
Strontium Ruthenate, Sr 2 RuO 4 , has played a leading role in discussions of unconventional superconductivity since its discovery almost three decades ago [1][2][3][4][5][6]. Much of the interest in the community centered on the possibility of chiral p-wave pairing, but the compound has also attracted attention simply because of its structural similarity to the cuprates, Fermi liquid behavior at low temperatures, and the availability of very clean samples with high quality surfaces. Recently, several new experimental results [7][8][9][10] have called into question the NMR results on which the traditional triplet pairing scenario was based [11], providing evidence for spin-singlet rather than triplet pairing and leading to a renaissance in the quest to identify the exact pairing state of this fascinating material. In principle, direct measurement of the superconducting gap by, e.g., Angular Resolved Photoemission Spectroscopy (ARPES), could provide important guidance, as it did in the cuprates. However the energy scales involved, such as the transition temperature of only 1.5K in Sr 2 RuO 4 , or the temperature at which the metamagnetic transitions in Sr 3 Ru 2 O 7 occur, are beyond the capabilities of current ARPES instruments. STM is a more appropriate tool, which due to its very high energy resolution that can be achieved at low temperatures and the ability to obtain information about the momentumand phase-resolved structure of the superconducting gap through quasiparticle interference (QPI) imaging [12,13] promises to resolve the most pressing questions about the superconducting properties of Sr 2 RuO 4 . * These authors contributed equally.
Significant progress has recently been made towards achieving this goal. The electronic structure in the bulk of Sr 2 RuO 4 near the Fermi energy is well-known to consist of weakly hybridized 1D sheets (α and β) of d xz/yz character, as well as a 2D d xy sheet (γ) that hybridizes with both. The γ-band has a van-Hove singularity which in the bulk is ∼ 14 meV above the Fermi energy [14], but whose energy depends sensitively on small structural changes [15] with significant consequences for the superconductivity [16]. This van Hove singularity plays not only an important role in the properties of Sr 2 RuO 4 , but also of the bilayer and trilayer ruthenates, where the van-Hove singularity has been suggested to be the origin of the metamagnetic behaviour [17].
In STM measurements, the situation has been less clear, and not all bands found in the bulk have been detected so far: Firmo et al. [19] argued that the gap observed in tunneling corresponds to that on the 1D d xz /d yz bands, with the d xy band not contributing to tunneling spectra but still exhibiting a sizeable gap due to proximity coupling. Tunneling to the STM tip is argued to happen primarily due to the d xz /d yz states in the sample, with the justification that d xy states associated with the γ-band have lobes that lie in the plane, while d xz /d yz states have lobes pointing out of the surface plane towards the tip (see also Fig. 1(a)). Similarly, in recent QPI experiments in the normal state of Sr 2 RuO 4 , the observed patterns were attributed to bands with d xz and d yz character [20]. The expected Bogoliubov-QPI in Sr 2 RuO 4 for a chiral order parameter has been previously investigated theoretically within a lattice Green's function framework that neglected the surface reconstruction and effect of the tunneling matrix elements [21]. A recent attempt to characterize the momentum-space struc- For an in-plane atomic orbital, such as dxy, the overlap with the tip orbitals is small, suggesting negligible contribution. (b) High-resolution differential conductance (dI/dV ) spectrum at T = 76mK, showing the characteristic gap-like structure at the surface of Sr2RuO4 (Vset = 8mV, Iset = 500.2pA, VL = 155µV). Four peaks can be clearly identified which are associated with the dxy-derived γ-band [18]. (c) Sketch of the structure of the reconstructed surface of Sr2RuO4 with RuO6 octahedra rotated by the angle ϑ = 6 • . The experimentally observed checkerboard charge order (right hand side) is equivalent to the breaking of C4-symmetry rendering the two Sr atoms (dark blue and light blue) and oxygen bonds along the horizontal and vertical direction inequivalent (indicated by orange and yellow circles) [18]. This symmetry breaking is described by a nematic order with d x 2 −y 2 symmetry (∆nem, red and blue indicate positive and negative sign, respectively, of the order parameter). The checkerboard charge order can be accounted for either through an additional staggered bond-centred order (∆ bond , colors as for the nematic order parameter), or through a staggered on-site order with different on-site energies for the dxy band at Ru (1) and Ru (2).
ture of the superconducting gap was made by Sharma et al. [22], who propose that their data was consistent with d-wave pairing and the signal from QPI due to the d xz /d yz bands. These results thus all raise the question what the role of the γ sheet is which has escaped detection in STM experiments so far. Detection of the γ band will allow one to to decide how large the gap is on this sheet, and understand whether it arises only from coupling to the α and β bands as argued in Ref. 19, or whether the STM tip is simply insensitive to states of d xy -character [20,22]. Recent STM experiments report signatures of the d xy band in tunneling spectra [18] (Fig. 1(b)) and an apparent absence of the superconducting gap, raising important questions about its role in superconductivity. Understanding and reconciling these seemingly contradictory interpretations is of primary importance in the effort to understand the QPI in this material and ultimately determine the momentum space structure of the superconducting gap in Sr 2 RuO 4 .
An important aspect of studies of clean surfaces of Sr 2 RuO 4 is the surface reconstruction, which arises as a spontaneous rotation of the RuO 6 octahedra similar to the crystal structure in the bulk of Sr 3 Ru 2 O 7 [23]. This reconstruction has recently been shown to influence QPI patterns on Sr 2 RuO 4 , and is possibly relevant for the low-energy electronic structure [18,20].
Here we resolve the mystery of the seeming absence of the γ-band from tunneling spectroscopy through combination of a phenomenological model of the low en-ergy electronic structure with ab-initio calculations of the tunneling matrix elements and ultra-low temperature scanning tunneling microscopy experiments. Our results demonstrate QPI of the γ-band and settle the normal state electronic structure and QPI in the surface layer of Sr 2 RuO 4 and thus provide a reference for modelling of the superconducting QPI. We show that for the unreconstructed surface tunneling into d xy states is suppressed primarily due to the alternating character of the Bloch d xy function at momenta near the van Hove point, and to a lesser extent due to the weaker extension of this Wannier function in the z direction. In the case of the reconstructed surface, we find that the amplitude of the d xy tunneling is enhanced due to an admixture of d z 2 orbital character mediated by the rotation of the oxygen octahedra.

II. RESULTS
A. Open questions -coupling to dxy states and emergent orders As discussed in the introduction, we begin with the premise that understanding the QPI in the normal state [18,20,22] will be essential to identification of the symmetry and structure of the superconducting gap in Sr 2 RuO 4 by STM. There are several features of the measured patterns that are challenging to interpret. The first is the dramatic suppression of the features in the mea-sured spectrum that originate from bands with dominant d xy orbital content, compared to a calculation of N (q, ω) using the lattice Green's function with bulk electronic bands as done in Refs. [20,22], where QPI was modelled by ignoring any d xy contribution to the trace and hence density of states. This suppression of scattering features associated with the γ-band has been discussed in terms of d xy orbitals coupling weakly to the tip due to the location of their lobes in the xy plane [19], or to orbital-selective decoherence of these orbitals [24]. As a practical matter, QPI data has often been analyzed simply ignoring d xy contributions [20], which seems to work up to a point. Nevertheless, there are q-peaks observed in these patterns that correspond to scattering from Fermi surface points close to the van Hove singularities dominated by d xy states [18], so a complete theory needs to account for these features as well.
A second set of puzzles is associated with the checkerboard charge order and associated nematicity of the Sr 2 RuO 4 surface [18]. The phenomena associated with these include a chirality of impurity states emanating from defects on different sublattices, and the bias dependence of the intensity of the atomic peaks. While these have been definitively associated with the reconstructed surface, their origin is unclear. Here we show that all these phenomena can be explained in a natural way with a combination of two types of coexisting orders, a nematic order and a staggered bond order, and by accounting for the vacuum tail of the involved electronic states.

B. Tight-binding model
We start from a model constructed from Wannier functions of the three Ru orbitals d xz , d yz and d xy which have been established as the relevant electronic states in the normal state Fermi surface of Sr 2 RuO 4 by ARPES experiments [25][26][27][28]. The surface Wannier states and the corresponding tight binding model are obtained from an ab-initio calculation (see Methods section for additional information). The lattice Hamiltonian is given by where t αβ R,R are hopping elements between the elementary cells described by the vectors R and R . See Supplementary Note 1 and Supplementary Fig. 1 for a plot of the band structure and Supplementary Note 2 and Supplementary Fig. 2 for a discussion of the Wannier functions. α and β are combined orbital and spin indices, the chemical potential enters as an on-site term and spinorbit coupling has been added as (complex-valued) onsite terms to represent H SOC = λL · S in the usual way, see Methods section for details.

C. Nematicity and checkerboard charge order
Unlike theoretical models for bulk Sr 2 RuO 4 , the surface reconstruction requires adoption of an elementary cell with two Ru atoms to describe the surface electronic structure as observed in STM and ARPES experiments. Therefore, we work in a basis that contains 6 Wannier states per spin, three of them centered at Ru(1) and three at Ru(2), see Fig. 1(c). As suggested by the crystal structure of the related material Sr 3 Ru 2 O 7 , in the surface layer the oxygen octahedra around Ru(1) atoms rotate anticlockwise and clockwise around Ru(2) (details on the resulting crystallographic structure are given in Appendix A 1), giving rise to the reconstructed surface and yielding qualitatively different tunneling properties as we discuss later. To fully describe the low energy electronic structure, two order parameters are required: a nematic term and a term describing the checkerboard charge order observed experimentally. Such terms could in principle be understood from microscopic models as instabilities of the electronic structure as worked out for the related material Sr 3 Ru 2 O 7 [29][30][31] and do modify hopping amplitudes of the d xy orbital. The nematic term affects the nearest neighbor (NN) hopping between Ru atoms. To describe the checkerboard charge order, we introduce a staggered bond order on the next nearest neighbor (NNN) hopping terms between Ru atoms, see Fig. 1(c). The symmetry properties of the two orders are identical once the octahedral rotation is present in the reconstructed surface layer [18]. The explicit form of the Bloch Hamiltonian as 6 × 6 matrix is not convenient to discuss at this point since the two sublattice basis already breaks the C 4 symmetry. We therefore discuss the symmetries in real space: The NN nematic term is of d x 2 −y 2 symmetry (see Appendix A 3) and induces a positive shift of the hopping amplitude along [0,1] (red ellipse connecting NN Ru atoms) and a negative shift of the hopping amplitude along [1,0] (blue ellipses) of amplitude ±∆ nem , preserving the translational symmetry between the Ru(1) and Ru(2) atoms.
The maximum contrast of the checkerboard pattern is on the Sr atoms, [18] suggesting that it is linked to the second-nearest neighbour (NNN) hopping parameter in the tightbinding model, motivating a description through a staggered bond order (see Fig. 1(c)). This NNN bond order breaks the translational symmetry on the Ru lattice, introducing a staggered next-nearest-neighbour interaction which is alternatingly strengthened and reduced by ∆ bond . The effect of the nematic and staggered bond order terms on the low energy electronic structure is comparable to the staggered on-site order proposed in Ref. 18   text, but show key results for the staggered on-site order in Supplementary Note 3 and Supplementary Fig. 3.

D. Tunneling into dxy states
In Fig. 2(a), we plot cuts through the Wannier functions obtained via downfolding a DFT calculation of the unreconstructed surface of Sr 2 RuO 4 onto a low-energy band structure consisting of three d-orbitals. The full isosurfaces of these rather complicated functions are shown in Supplementary Fig. 2, but at the location of the STM tip, some 4−5Å above the SrO surface plane, they resemble atomic d-orbitals. Note that the d xz Wannier function has one maximum roughly half way to the NN Ru atomic positions and the d xy Wannier function is much smaller in magnitude and vanishes by symmetry above the NN Ru atom. The octahedral rotation in the surface layer leads to a number of important changes in the electronic states associated with the d xy band: (1) the van Hove singularity in the d xy band shifts below the Fermi energy, (2) the d xy band acquires d z 2 character with increased octahedral rotation and (3) the Wannier functions in the vacuum become chiral, with opposite chirality on Ru(1) and Ru(2) atoms. In Fig. 2(b), we show how the Wan-nier functions appear at the tip position in the presence of a 6 • octahedral rotation. While the d xz,yz states are not qualitatively altered, the Wannier functions associated with the γ band acquire a chiral character such that they no longer vanish above the NN Ru positions (light red dot). The Wannier functions shown correspond to a Ru(2) position, with the function associated with Ru(1) having the opposite chirality. Fig. 2(c) shows the electronic structure for an octahedral rotation of ϑ = 6 • and the projected density of states for different octahedral rotations. The van-Hove singularity in the d xy -derived γ band has moved across the Fermi level compared to the unrotated case and has acquired a significant d z 2 character, especially close to the M point, as a consequence of the octahedral rotation. The van-Hove singularity at the M point does not have any d z 2 character in its projected density of states (PDOS) without rotation. These findings do not change in a fully relativistic ab-initio calculation. It is through this admixture of d z 2 character with the octahedral rotation that tunneling into the γ band is facilitated and gives rise to peaks in the cLDOS from the vHss. These are absent when employing Wannier functions obtained without oxygen rotation (dashed lines in Fig. 2(e)).  tions, we plot in Fig. 3(a) the square of the Wannier function integrated over the x − y plane, |W | 2 , as a function of z, corresponding to the tip height in a tunneling experiment. Once sufficiently away from the surface, for z > 2Å, the Wannier functions show the expected exponential decay, i.e. |W | 2 ∝ exp(−z/α). For the case without octahedral rotation, the d xy weight in vacuum at values of z relevant for tunneling (assumed here to be typically at 5Å above the surface, but the exact height is irrelevant for our analysis) is an order of magnitude smaller than the weight of d xz,yz , as anticipated in Firmo et al. [19] However, once the octahedral rotation is considered, the weight of the d xy orbital is only about 3 times smaller and exhibits a decay length which is 10% larger compared to the d xz/yz orbitals (decay length α xy = 2.2Å vs. α xz,yz = 2.0Å).  The decay length for the d xy orbital thus changes from a value smaller than that of the d xz,yz orbitals to a larger value due to the rotation. The suppression of the vacuum overlap in the unreconstructed state alone is therefore not sufficient to explain the lack of most d xy features in QPI. We will show below that the d xy states contribute most strongly near the van Hove point, close to k = ( 1 2 , 0), thus tunneling should be proportional to the value of the Wannier functions at the centre of the NN Ru atom. Without rotation of the oxygen octahedra, this is zero by symmetry, see Fig. 2(a), i.e. no tunneling is expected from states close to k = ( 1 2 , 0). This is also seen by the absence of any features due to the vHs in the cLDOS when calculated without rotation of the oxygen octahedra, see Fig. 2(e).

E. Checkerboard charge order
In experiments, one of the most prominent features of the surface electronic structure is a pronounced Srcentred checkerboard charge order, which in our model is accounted for through the staggered bond order. In Fig. 4(a) and (b) we show measured differential conduc-tance maps at positive and negative bias voltages, respectively, in comparison to calculated maps (Fig. 4(c, d)) of the continuum local density of states at a constant height above the surface, fully accounting for the vacuum tail of the wave functions. The LDOS maps demonstrate that the staggered bond order indeed leads to a checkerboard charge order centred on the Sr atoms as found experimentally, and reproduces the contrast inversion between positive and negative bias voltages (compare Fig. 4(c,  d)). A notable difference between the experimental and calculated maps is that the experimental data is dominated by the checkerboard charge order, whereas the calculated maps show the checkerboard charge order as a subdominant contribution superimposed to the atomic contrast. We attribute this difference between the experimental and calculated maps to the different treatment of the tip-sample distance: in our measurements, the tip-sample distance is set at each point independently to yield a constant current, whereas in the calculations shown in Fig. 4(c, d) the local density of states is taken at constant height. To faithfully reproduce the experimental data requires calculating differential conductance maps where the tip height is locally adjusted to maintain a constant integral Eset 0 ρ(E)dE, emulating the effect of the feedback loop of an STM regulating on a constant tunneling current before the spectrum is recorded. Such maps are shown in Fig. 4(e, f) for the same energies as in (c, d), showing a complete suppression of the atomic contrast. At positive bias voltages, we find excellent agreement between the calculated and measured differential conductance maps ( Fig. 4(a, c)), whereas at negative bias voltage the agreement is not quite as good: The LDOS map in Fig. 4(d) reproduces the checkerboard charge order as seen experimentally, but the calculated differential conductance map (Fig. 4(f)) shows the dominant contrast on top of the ruthenium atoms. This difference is likely an artifact because the simulated tipsample distance is significantly smaller than the one expected for the experiment. The theoretically tractable tip-sample distances are limited due to technical reasons related to the accuracy of the Wannier functions at large distances, indications for quantitative changes are given by the (slightly) larger decay length of the d xy Wannier function. Our results show that the setpoint effect, which is normally considered detrimental to the interpretation of spectroscopic maps, suppresses the atomic corrugation in the differential conductance maps. For the following comparison of the quasi-particle interference, we have verified that the main impact of the setpoint effect is a suppression of the atomic contrast in the differential conductance maps, otherwise not affecting the signal due to quasi-particle interference significantly. For comparison, we show in the supplementary material (Suppl.

F. Quasiparticle Interference
In order to provide a full picture of the low energy electronic structure around the Fermi energy, we use QPI imaging and compare our experimental QPI maps to the simulated continuum LDOS maps. Continuum LDOS maps in real space and simulated topographies exhibit chiral QPI patterns around Ru-site defects, as also found experimentally. The rotational sense of these patterns depends on the position of the defect in a Ru(1) or Ru(2) site (see Supplementary Fig. 4 and Supplementary Note 4). For comparison with the experiment, we average over both types of defects, as also the experimental data is acquired over fields of view with defects in both sites. Comparison of the simulated QPI, fully taking into account the tunneling matrix elements through the Wannier functions (for details on the method see Appendix A 4), reveals excellent overall agreement between experiment and theory (see Fig. 5). We note that the positions of features in q-space deviate slightly between theory and experiment. For example the outer dominant ring-like structure is larger in the calculation. This is because the bare electronic structure from the first principles calcu-lations does not to match exactly the true Fermi surface, a deviation which is not relevant for the following discussion. The key features are qualitatively consistent between theory and experiment: the outer square-shaped scattering from the quasi-1D bands, and the inner square coming from the scattering processes crossing the zoneboundary. At low q-vectors, an intensity distribution with C 2 symmetry is seen which is reproduced in the calculation. As pointed out previously [20], the appearance of spectral weight in some parts of the BZ can only be understood accounting for the surface reconstruction; these are essentially the structures parallel to q x and q y in Fig. 5, which can be traced back to bands with predominantly d xz /d yz character (see also orbital decomposition in fig. 5). By contrast, as can be seen in Fig. 5 there are scattering processes close to the atomic and reconstruction peaks associated with the d xy -derived γ-band, the band which exhibits the van Hove singularity close to the Fermi energy. In the following, we will use the model calculations to establish the signatures of the γ-band and the van-Hove singularity in QPI.

QPI of the van Hove singularity
The QPI signal of the γ-band and the vHs is dominated by scattering vectors connecting the tips of the constant energy contours close to the vHs. Due to the background near q → 0 from the impurity distribution, the small q vectors connecting the tips near the vHs are difficult to detect reliably. We have identified two ways to still accurately detect the dispersion close to the vHs and determine the energy of the vHs: (1) the scattering vectors connecting the points of highest density of states include ones which cross the Brillouin zone, leading to QPI features around the atomic Bragg peaks, where the noise background is much lower. (2) At the van Hove singularity, the scattering vector becomes commensurate with the atomic contrast, resulting in a resonant enhancement of the atomic contrast when the energy becomes equal to that of the vHs. Fig. 6(a) shows the Fermi surface extracted from our tight-binding model, with one of the scattering vectors with high joint density of states connecting points near the van-Hove singularity leading to QPI around the atomic peaks. This scattering vector is already apparent from the QPI map as a feature in close proximity to the atomic peaks, see Fig. 6(b). From line cuts through the three-dimensional energy-momentum data along the q x and q y direction, Fig. 6(c, d), the dispersion of these peaks can be tracked. A clear hole-like dispersion is observed with a band maximum a few millivolts above the Fermi energy. Fig. 6(f, g) show for comparison the same cuts obtained from the calculations. While the calculations exhibit more fine structure than seen in the experiment, the main feature of a hole-like dispersion around the atomic peak is faithfully reproduced. As expected from the nematicity of the electronic structure, the van-Hove singularities occur at slightly different energies along the q x -and q y -directions, providing an estimate of the magnitude of the nematic term ∆ nem in the Hamiltonian. The nematicity also leads to a pronounced anisotropy of the low-q QPI. As a function of energy, the contributions from small q scattering vectors crossing the zone boundary are expected to evolve according to the touching of the (reconstructed) bands. In real space, these correspond to interference patterns with long wavelength and rotation of the dominant wave vector as a function of energy, as noted in Ref. 18. Calculated maps of the continuum LDOS confirm this feature and allow us to unequivocally assign it to tunneling into the d xy -derived γ-band since the van Hove singularity responsible for this low-q scattering occurs only in this orbital channel. Notably, these results confirm our theoretical conjecture that the γ band becomes detectable in tunneling, facilitated by the octahedral rotations.
The crossing of the QPI signal of the dispersion of the γ-band through the atomic peak seen in Fig. 6(c, d) provides an alternative measure of the van Hove singularities in the electronic structure: at the energy of a vHs at the zone boundary, the quasi-particle scattering becomes commensurate with the atomic periodicity leading to a significant increase in the intensity of the atomic peaks in spectroscopic maps. In Fig. 6(c, d), this becomes apparent as saturation of the contrast at distinct points along the q at = (1, 0) and (0, 1) line. For clarity, we plot the energy dependence of the QPI signalg(q at , V ) as a function of energy eV in Fig. 6(e). Traditional (lattice-only) T-matrix calculations are unable to capture this feature because they exhibit the same periodicity as the Brillouin zone, and the QPI signal at the atomic peak is thus identical to the one at the zone centre. The continuum LDOS is able to describe this intensity modulation: in our calculations, the intensities of the atomic peaks at q at in maps of the LDOSρ(q, E) show sharp peaks when the γ band crosses the zone boundary (compare Fig. 6(h)), providing an alternative way to determine the energy of the van Hove singularities in the electronic structure without the necessity of undertaking a full QPI mapping.

III. DISCUSSION
Our theoretical modelling and measurements provide a comprehensive picture of the low energy electronic structure of the surface layer of Sr 2 RuO 4 , and identify clear signature of the γ-band in QPI with potential implications for its superconducting state. Quasi-particle interference of the γ-band which is predominantly of d xycharacter had hitherto been assumed to only contribute negligibly to the tunneling signal. Our measurements show a clear QPI signal from this band through comparison with theory. From the calculations for a hypothetical unreconstructed surface, without octahedral rotation, we indeed find that tunneling to the γ band would be about an order of magnitude smaller than for the d xz and d yz bands, leading to a negligible contribution to the tunneling conductance. The small tunneling probability is due to the real space properties of the d xy Wannier function and the oscillatory nature of the Bloch wave function near k = ( 1 2 , 0). In the reconstructed surface, the octahedral rotation leads to additional d z 2 weight for the γ-band, making it accessible by QPI.
In summary we have identified the physical ingredients necessary to describe the electronic structure at the surface of Sr 2 RuO 4 which include spin-orbit coupling, the nematic and staggered charge orders as well as the rotation of the oxygen octathedra. While the terms in the Hamiltonian outlined in detail in the Methods section can in principle be derived from a microscopic description, the effective parameters are subject to renormalizations due to the strongly correlated nature of the material. Despite the excellent agreement we observe between the experimental data and theoretical modelling, one notable difference remains: while in tunneling spectroscopy the gap-like structure around the Fermi energy leads to a significant suppression of differential conductance by about 40%, this is not accurately captured in the calculation, where the suppression remains significantly smaller. This can have a number of origins, including that the calculations are carried out for smaller tip-sample distances than used in the experiments, potential additional relaxation of the surface layer in the z-direction and that a larger part of the Fermi surface becomes gapped out than is captured in the model. Nevertheless, our measurements clearly demonstrate that all three bands of the t 2g manifold (d xz , d yz and d xy ) which are expected to be present at the Fermi energy can be detected in QPI. We also note that our measurements, taken at 59 mK, do not show evidence of a superconducting gap. This is particularly surprising given that all three bands which contribute to the Fermi surface are clearly detected. The absence of spectral features from such a gap in many high-resolution STM experiments on the SrO-terminated surface of Sr 2 RuO 4 [18,32,33] remains an important open puzzle. One possibility highlighted by our analysis is that superconductivity is suppressed if the surface is reconstructed by octahedral rotation, as assumed here. It is conceivable that disorder, or other subtle surface effects, may lead to other reconstructions that do not suppress superconductivity -calling for new ways to suppress the surface reconstruction, possibly through adsorbate layers to facilitate a detection of the superconducting gap in tunneling experiments and thus determination of the superconducting order parameterproviding a resolution to the long-standing mystery of the symmetry of the superconducting order parameter in Sr 2 RuO 4 .
Appendix A: Methods

First principles Wannier calculations
Density functional theory calculations [34] were performed with the projected augmented wave (PAW) method as implemented in the Vienna ab initio simulation package (VASP) [35,36]. The generalized gradient approximation of Perdew, Burke and Ernzerhof was used for the exchange correlation functional [37]. To be able to capture the Wannier functions high above the surface we performed a monolayer calculation of perovskite Sr 2 RuO 4 . The lattice constant was taken a = 3.87Å, and the vacuum length was chosen to be about 21Å. To incorporate the rotations we construct a √ 2 × √ 2 supercell of the Sr 2 RuO 4 monolayer unit cell. We perform two calculations, one without a rotation and one with a ϑ = 6 • rotation of the Ru-O in-plane bonds. In the rotation of the O atoms we keep the lattice constants, the Ru positions and the Ru-O bond angles fixed suth that the O position is then given by (x, y, 0) = 0.25(tan ϑ + 1, tan ϑ − 1, 0). The energy cutoff of the plane waves was chosen as 650 eV. The total energy was converged to 10 −7 eV. The Brillouin zone integration was sampled by using a 7 × 7 × 1 Γ-centered Monkhorst-Pack grid. To construct the Wannier functions and the tight-binding models, the d xz , d yz and d xy orbitals were projected on the low energy bands, employing the Wannier90 code package [38] with input parameters num_iter=0 and dis_num_iter=10000. The outer energy window was taken as [-3,1] eV and the inner frozen energy window as [-1.7,0.2] eV, both relative to the Fermi level.
2. Surface tight-binding Hamiltonian with spin-orbit coupling As described in the main text, the Hamiltonian matrix in momentum space is given by where H tb (k) = 1 spin ⊗t(k) ab is the Bloch Hamiltonian as obtained from the ab initio Wannier calculation (including the ϑ = 6 • rotation of the Ru-O in-plane bonds) with an overall band renormalization of Z = 1/4 to match experimentally observed Fermi velocities [25]. The 6 × 6 matrix t(k) ab is the Fourier representation of the in-plane hoppings of our two dimensional model for the orbitals a and b with a, b = (d xy , d (2) xz , d (2) yz , d xy ) and the superscript denotes the sublattice index of the Ru atoms.
We introduce spin-orbit coupling via an onsite spinorbit coupling term proportional to the unit matrix in the sublattice space, where the 3 × 3 matrices stem from the product of the spin operator S = 1 2 (σ x , σ y , σ z ) and the representation of the angular momentum matrices L µ in the basis of the real-valued d orbitals, H We use λ = 20 meV to yield agreement of the splittings from the avoided crossings of the quasi 1 dimensional d xz and d yz bands along the path from Γ to the M point (see Fig. 6) with the measured spectral function of the surface bands [25].

Nematicity and Staggered orders
We introduce the nematic and staggered orders through an additional term H nem (k) in the Hamiltonian, The nematic and bond order contributions are diagonal in spin space, H nem (k) = 1 spin ⊗ t nem (k). The matrix t nem (k) has only nonzero components in the d xy orbital components such that the sub-matrix in this subspace, spanned by the elements 3 and 6, reads t nem (k)| [3,6], [3,6] = where the momenta are in the Brillouin zone obtained from an elementary cell with two Ru atoms, see Fig. 1(c) and we chose ∆ nem = 2.5 meV, ∆ bond = 5 meV. To describe the staggered onsite order, we use the same expression with f bond (k) = 1.

Wannier method for simulations of tunneling
Here we adopt a Wannier-function based approach which allows us to relate the tunneling rate to the local density of states in vacuum a fewÅ above the surface of the SrO layer, where the STM tip is positioned. The current at voltage V (or differential conductance) as measured in an STM experiment can be calculated by where A 0 is a constant containing the tip density of states and the tunneling matrix element, and ρ(r, ω) is the continuum local density of states (cLDOS) at the tip position r = (x, y, z) which is assumed to be at severalÅ above the surface atoms of Sr 2 RuO 4 . The cLDOS can be calculated conveniently from the continuum Green's function G(r, r ; ω) = ψ(r) † ψ(r ) ω (where ψ(r) = R,µ c Rµ w Rµ are the continuum electron operators) via [39,40] ρ(r, ) = − 1 π Im G(r, r; ).
Usually, the electronic structure is discussed using the lattice (tight-binding) model with a lattice Green's func-tionĜ µ,ν R,R (ω), a matrix in the combined orbital and spin space µ, lattice position R. The continuum Green function can be calculated using a basis transformation as where the matrix elements w Rµ (r) are the Wannier functions which are obtained in the tight-binding downfolding for H tb (k) as well. Finally, let us mention that the differential conductance (at constant tip height) is obtained by taking the derivative of Eq. (A8) with respect to the bias voltage, yielding the proportionality dI/dV ∝ ρ(r, eV ). A topographic map z(x, y) as obtained experimentally by keeping the current I 0 constant for a given bias voltage V 0 , can be calculated by solving the equation for z(x, y) which requires the evaluation of the continuum LDOS within a height range and for all energies to carry out the integral. Finally, for realistic calculations of conductance maps as obtained in topographic mode, one evaluates continuum LDOS at the height profile [41], i.e, ρ t (x, y, eV ) = ρ((x, y, z(x, y)), eV ) .
In a homogeneous system, the lattice Green function is translation invariant,Ĝ R,R =Ĝ 0 R−R (ω), where the r.h.s. can be calculated by Fourier transform from the Green function in momentum spacê For calculations including impurities, we consider a simple potential scatterer at a Ru atom with lattice position i * , diagonal in the combined orbital (sublattice) and spin space with the Hamiltonian Within the T-matrix approach, the lattice Green function is given by [40,41]  . The local Green function is just given byĜ(ω) =Ĝ 0 R=0 (ω). Note that all frequency arguments in the Green functions are shorthand notations for ω + iη with an energy broadening η which we choose to be sub meV to achieve satisfactory energy resolution by use of k grids of size 2500 × 2500 (or 3500 × 3500 for Fig. 2(d,e)) . Finally, lattice density of states can be calculated as trace, N (q, ω) = − 1 π Im TrĜ 0 q (ω).

Ultra-low temperature STM
Quasi-particle interference imaging has been performed using a dilution-refrigerator based low temperature STM operating at temperatures down to 10mK [42]. Spectroscopic maps g(r, V ) shown here, where a tunneling spectrum is recorded at each point of a topographic image z(r), were recorded by stabilizing the tip sample distance at a tunneling setpoint V set , I set before switching the feedback loop off to record the spectrum. The differential conductance was recorded through a lock-in technique, adding a voltage modulation V L to the bias voltage. While some of the maps shown here where recorded in magnetic field (as indicated in figure captions), the qualitative behaviour of the features reported here is not affected by the field. The magnetic field shifts some features slightly in energy (as reported in [18]); however the shifts are small compared to the overall energy scale of the comparison with theory.

Sample growth and characterization
Single crystals of Sr 2 RuO 4 were grown by a flux feeding floating zone (FFFZ) with Ru self-flux using a commercial image furnace equipped with double elliptical mirrors and two 2.0 kW halogen lamps (NEC Machinery, model SC1-MDH11020). Details of the FFFZ crystal growth are described in detail elsewhere [43][44][45]. Several techniques, including x-ray diffraction, energy dispersive spectroscopy (EDS) and polarized light optical microscopy (PLOM) analysis, have been used to fully characterize the structure, quality, and purity of the crystals.
Data availability: The data underpinning the findings of this study are available online [46].
Code availability: The computational data and source code are available upon reasonable request to A. Kreisel (kreisel@itp.uni-leipzig.de).  Figure 1(a)) and a fat-band plot showing the orbital character of the bands crossing the Fermi energy (Supp. Fig. 1(b)).  (e, f) Calculated differential conductance map ρt(r, E) at E = ±5meV, for Eset = 10meV for the same model. The intensity of all images has been normalized by the spatial average, color bars indicate relative intensity.

Supplementary Note 3. STAGGERED ON-SITE ORDER
In this supplementary note, we discuss the effect of different staggered order parameters on the calculated tunneling. Employing the staggered on-site order (see Methods section of the main text), there are few differences in the resulting spectra and maps as we demonstrate in Supplementary Figure 3 which compares the experimental data as presented in the main text ( Fig. 4) with the theoretical results for the continuum LDOS ρ(r, E) and the simulated differential conductance maps ρ t (r, E) yielding qualitatively similar results, with minor quantitative differences.