Quasi-particle interference of the van Hove singularity in Sr2RuO4

The single-layered ruthenate Sr2RuO4 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 Sr2RuO4 by quasi-particle 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.


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 highquality 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 spinsinglet 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.5 K 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 momentum-and phaseresolved 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 .
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 (vHs) 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 vHs has 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 behavior 17 .
In scanning tunneling microscopy (STM) measurements, the situation has been less clear, and not all bands found in the bulk have been detected so far: Firmo et al. 18 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. 1a). 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 19 . 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. 20 A recent attempt to characterize the momentum-space structure of the superconducting gap was made by Sharma et al. 21 , who propose that their data were 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 of what the role of the γ sheet is which has escaped detection in STM experiments so far. Detection of the γ band will allow one 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. 18 , or whether the STM tip is simply insensitive to states of d xy -character 19,21 .
Recent STM experiments report signatures of the d xy band in tunneling spectra 22 (Fig. 1b) 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 19,22 .
Here, we resolve the mystery of the seeming absence of the γband from tunneling spectroscopy through a combination of a phenomenological model of the low-energy 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 modeling 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.

RESULTS
Open questions-coupling to d xy states and emergent orders As discussed in the introduction, we begin with the premise that understanding the QPI in the normal state 19,21,22 will be essential to the 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 measured spectrum that originate from bands with dominant d xy orbital content, compared with a calculation of N(q, ω) using the lattice Green's function with bulk electronic bands as done in refs. 19,21 , where QPI was modeled 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 18 , or to orbitalselective decoherence of these orbitals 24 . As a practical matter, QPI data has often been analyzed simply ignoring d xy contributions 19 , 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 22 , so a complete theory needs to account for these features as well.
The second set of puzzles is associated with the checkerboard charge order and associated nematicity of the Sr 2 RuO 4 surface 22 . 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. Although 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.

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). The lattice Hamiltonian is given by where t αβ R;R 0 are hopping elements between the elementary cells described by the vectors R and R 0 . (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) on-site terms to represent H SOC = λL ⋅ S in the usual way, see Methods section for details. The out-of-plane symmetry of a d yz orbital suggests a high probability of electron tunneling from a metallic tip and a significant contribution to the tunneling current. For an in-plane atomic orbital, such as d xy , the overlap with the tip orbitals is small, suggesting negligible contribution. b High-resolution differential conductance (dI/dV) spectrum at T = 76 mK, showing the characteristic gap-like structure at the surface of Sr 2 RuO 4 (V set = 8 mV, I set = 500.2 pA, V L = 155 μV). Four peaks can be clearly identified which are associated with the d xy -derived γ-band 22 . c Sketch of the structure of the reconstructed surface of Sr 2 RuO 4 with RuO 6 octahedra rotated by the angle ϑ = 6 ∘ . The experimentally observed checkerboard charge order (right-hand side) is equivalent to the breaking of C 4 -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) 22 . 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-centered order (Δ bond , colors as for the nematic order parameter), or through a staggered onsite order with different on-site energies for the d xy band at Ru(1) and Ru(2).
A. Kreisel et al.

Nematicity and checkerboard charge order
Unlike theoretical models for bulk Sr 2 RuO 4 , the surface reconstruction requires the 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 six Wannier states per spin, three of them centered at Ru (1) and three at Ru(2), see Fig. 1c. 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 the Methods section), 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-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. 1c. The symmetry properties of the two orders are identical once the octahedral rotation is present in the reconstructed surface layer 22 . 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 Methods section) 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, 22 suggesting that it is linked to the second nearest neighbor (NNN) hopping parameter in the tight-binding model, motivating a description through a staggered bond order (see Fig. 1c). This NNN bond order breaks the translational symmetry on the Ru lattice, introducing a staggered NNN 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. 22

Tunneling into d xy states
In Fig. 2a, 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 dorbitals. 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 Wannier functions appear at the tip position in the presence of a 6 ∘ octahedral rotation. Although 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.   (2) atoms. c DFT band structure for the reconstructed surface with ϑ = 6 ∘ employed in this work, with d z 2 orbital weight highlighted in blue. The right panel shows the projected d z 2 density of states for ϑ = 6 ∘ as in the panel on the left, as well as for ϑ = 3 ∘ and 0 ∘ (red, yellow). d Lattice density of states from the model described in the text (black: total; blue: d xy on the Ru(1) and Ru (2) atoms exhibiting peaks at the positions marked with arrows, green: d xz/yz (approximately fourfold degenerate) with only tiny structure within the energy range shown). e Continuum density of states (cLDOS) evaluated 5 Å above the surface above the two inequivalent Ru atoms and above the two inequivalent Sr atoms, see inset. Solid lines show results using Wannier functions with octahedral rotation and exhibit features at the peaks of the d xy lattice DOS, whereas a calculation using Wannier functions without octahedral rotation yields a completely flat cLDOS above the Ru positions within this energy range.
shows the electronic structure for an octahedral rotation of ϑ = 6 ∘ and the projected density of states for different octahedral rotations. The vHs in the d xy -derived γ band has moved across the Fermi level compared with 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 vHs at the M point does not have any d z 2 character in its projected density of states 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 continuum local density of states (cLDOS) from the vHss. These are absent when employing Wannier functions obtained without oxygen rotation (dashed lines in Fig. 2e).
As a rough estimate of possible tunneling contributions, we plot in Fig. 3a 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., jWj 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. 18 However, once the octahedral rotation is considered, the weight of the d xy orbital is only about three times smaller and exhibits a decay length which is 10% larger compared with 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 center of the NN Ru atom. Without rotation of the oxygen octahedra, this is zero by symmetry, see Fig. 2a, 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. 2e.

Checkerboard charge order
In experiments, one of the most prominent features of the surface electronic structure is a pronounced Sr-centered checkerboard charge order, which in our model is accounted for through the staggered bond order. In Fig. 4a, b, we show measured differential conductance maps at positive and negative bias voltages, respectively, in comparison to calculated maps (Fig. 4c, 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 centered on the Sr atoms as found experimentally, and reproduces the contrast inversion between positive and negative bias voltages (compare Fig. 4c, d). A notable difference between the experimental and calculated maps is that the experimental data are 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. 4c, d the cLDOS is taken at a 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 R 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. 4e, 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. 4a, c), whereas at negative bias voltage the agreement is not quite as good: the LDOS map in Fig. 4d reproduces the checkerboard charge order as seen experimentally, but the calculated differential conductance map (Fig. 4f) shows the dominant contrast on top of the ruthenium atoms. This difference is likely an artifact because the simulated tip-sample 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 QPI, 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 QPI significantly. For comparison, we show in the supplementary material (Suppl. Fig. 3 and Suppl. influence on the norm of the Wannier orbitals associated with the d xz/yz bands, it enhances the value of the Wannier orbital associated with the γ band, which has predominantly d xy character, significantly, and also increases the decay length in jWj 2 / expðÀz=αÞ from α xy = 1.9 Å to α xy = 2.2 Å. b Surface geometry of Sr 2 RuO 4 with the Ru and O atoms on the surface at z ≈ 0 (green and red dotted lines for these planes) and the STM tip~5 Å above (black dashed line).

A. Kreisel et al.
Note 3) the calculations for the staggered on-site order, showing very similar results.

Quasi-particle 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 see Methods section), 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 calculations does not match exactly the true Fermi surface, a deviation that 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 zone boundary. At low-q vectors, an intensity distribution with C 2 symmetry is seen which is reproduced in the calculation. As pointed out previously 19 , 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 vHs close to the Fermi energy. In the following, we will use the model calculations to establish the signatures of the γ-band and the vHs in QPI.
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 the highest density of states include ones that 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. Figure 6a shows the Fermi surface extracted from our tightbinding 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. 6b. From line cuts through the threedimensional energy-momentum data along the q x and q y direction, Fig. 6c, 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. Figure 6f, g shows for comparison the same cuts obtained from the calculations. Although the calculations exhibit more fine structure than seen in the experiment, the main feature of a hole-like dispersion A. Kreisel et al. around the atomic peak is faithfully reproduced. As expected from the nematicity of the electronic structure, the vHss 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. 22 . Calculated maps of the cLDOS confirm this feature and allow us to unequivocally assign it to tunneling into the d xy -derived γ-band since the vHs 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. 6c, 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. 6c, 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. 6e. 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 center. The cLDOS is able to describe this intensity modulation: in our calculations, the intensities of the atomic peaks at q at in maps of the cLDOSρðq; EÞ show sharp peaks when the γ band crosses the zone boundary (compare Fig. 6h, providing an alternative way to determine the energy of the vHss in the electronic structure without the necessity of undertaking a full QPI mapping.

DISCUSSION
Our theoretical modeling and measurements provide a comprehensive picture of the low-energy electronic structure of the surface layer of Sr 2 RuO 4 , and identify a clear signature of the γband in QPI with potential implications for its superconducting state. QPI of the γ-band which is predominantly of d xy -character 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 realspace 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. Although 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 modeling, one notable difference remains: while in tunneling spectroscopy the gap-like structure around the Fermi energy leads to significant suppression of differential conductance by~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 SrOterminated surface of Sr 2 RuO 4 22,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 detection of the superconducting gap in tunneling experiments and thus determination of the superconducting order parameter-providing a resolution to the long-standing mystery of the symmetry of the superconducting order parameter in Sr 2 RuO 4 .

First-principles Wannier calculations
Density functional theory calculations 34 were performed with the projected augmented wave method as implemented in the Vienna ab initio simulation package 35,36 . The generalized gradient approximation of 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 one-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 . Fig. 6 QPI close to the atomic peaks. a Fermi surface from the model with surface reconstruction. The color represents the orbital character of the bands, with d xz /d yz character shown in green/red, respectively, and d xy character in blue. The gray arrow indicates the q-vector connecting the tips of the pockets of the γ-band close to the van Hove singularity. The QPI dispersion from this q-vector is expected near the atomic peaks. b Fourier transformation of a differential conductance map at V = 2.52 mV. QPI features due to the q-vector shown in (a) (marked by an arrow) can be observed close to the Bragg peaks at (1, 0) (in units of 2π/a). c, d Energy-momentum cuts through the QPI map along q x (c) and q y (d) close to (1, 0) and (0, 1). A clear dispersing feature is seen which collapses onto the atomic peak (white arrows). The dispersion differs in the q x and q y directions as a consequence of nematicity (V set = 5.6 mV, I set = 225 pA, V L = 300 μV, T = 76 mK, B = 6.5 T). e Differential conductancegðq at ; VÞ for the atomic peaks q at = (1, 0) and (0, 1), showing prominent features at V = 2 mV and 3 mV where the vHs crosses the zone boundary. f, g Corresponding energy-momentum cuts from the model along q x (f) and q y (g), showing the signatures of the dispersion of the scattering vectors associated with the γ band around the atomic peaks. h As in (e), local density of statesρðq at ; EÞ at q at = (1, 0) and (0, 1) as a function of E. As in the experimental data, two maxima are seen due to the vHs at the zone boundary.

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 with f bond ðkÞ ¼ 1 2 ðcos k x þ cos k y Þ and f nem ðkÞ ¼ 1 4 ½1 þ exp½iðk x À k y Þ À expðÀik y Þ À expðik x Þ where the momenta are in the Brillouin zone obtained from an elementary cell with two Ru atoms, see Fig. 1c and we chose Δ nem = 2.5 meV, Δ bond = 5 meV. To describe the staggered on-site 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 that allows us to relate the tunneling rate to the local density of states in a 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 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 0 ; ωÞ ¼ hψðrÞ y ψðr 0 Þi ω (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 (tightbinding) model with a lattice Green's functionĜ μ;ν R;R 0 ðωÞ, a matrix in the combined orbital and spin space μ, lattice position R. The continuum Green function can be calculated using a basis transformation as Gðr; r 0 ; ωÞ ¼ X R;R 0 ;μνĜ μ;ν R;R 0 ðωÞw Rμ ðrÞw R 0 ν ðr 0 Þ; (11) 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. (9) 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 0 dω ρðx; y; zðx; yÞ; ωÞ ; for z(x, y) that requires the evaluation of the continuum LDOS (cLDOS) 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 cLDOS 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 ¼Ĝ 0 RÀR 0 ðωÞ, where the r.h.s. can be calculated by Fourier transform from the Green function in momentum spacê G 0 k ðωÞ ¼ ½HðkÞ À ω À1 ,Ĝ 0 R ðωÞ ¼ P kĜ 0 k ðωÞe iRÁk . 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 wherê is the T-matrix andV imp ¼ V imp 1 spin Ŝ is the matrix representation of Eq. (14). For an impurity on a Ru(1) atomic position, we useŜ ¼ 1 3 0 0 0 , whereas for the impurity on a Ru(2) atomic positionŜ ¼ 0 0 0 1 3 . 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. 2d, e). Finally, the lattice density of states can be calculated as trace, Nðq; ωÞ ¼ À 1 π ImTrĜ 0 q ðωÞ.
Ultra-low temperature STM QPI imaging has been performed using a dilution-refrigerator-based lowtemperature STM operating at temperatures down to 10 mK 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. Although some of the maps shown here were recorded in a magnetic field (as indicated in figure captions), the qualitative behavior of the features reported here is not affected by the field. The magnetic field shifts some features slightly in energy (as reported in ref. 22 ); however, the shifts are small compared with 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, and polarized light optical microscopy 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.