Magnetic-Field Tunable Intertwined Checkerboard Charge Order and Nematicity in the Surface Layer of Sr$_2$RuO$_4$

In strongly correlated electron materials, the electronic, spin, and charge degrees of freedom are closely intertwined. This often leads to the stabilization of emergent orders that are highly sensitive to external physical stimuli promising opportunities for technological applications. In perovskite ruthenates, this sensitivity manifests in dramatic changes of the physical properties with subtle structural details of the RuO$_6$ octahedra, stabilizing enigmatic correlated ground states, from a hotly debated superconducting state via electronic nematicity and metamagnetic quantum criticality to ferromagnetism. Here, it is demonstrated that the rotation of the RuO$_6$ octahedra in the surface layer of Sr$_2$RuO$_4$ generates new emergent orders not observed in the bulk material. Through atomic-scale spectroscopic characterization of the low-energy electronic states, four van Hove singularities are identified in the vicinity of the Fermi energy. The singularities can be directly linked to intertwined nematic and checkerboard charge order. Tuning of one of these van Hove singularities by magnetic field is demonstrated, suggesting that the surface layer undergoes a Lifshitz transition at a magnetic field of ~32T. The results establish the surface layer of Sr$_2$RuO$_4$ as an exciting 2D correlated electron system and highlight the opportunities for engineering the low-energy electronic states in these systems.


Introduction
The physical properties of strongly correlated electron materials often vary dramatically with comparatively modest external stimuli, [1] evidenced, for example, by magnetic-field driven metamagnetic transitions, [2,3] doping-induced metal-toinsulator transitions [4] and superconductivity [5] as well as a surprising sensitivity to uniaxial strain. [6] The changes in physical properties are usually accompanied by tiny structural distortions often reflecting, or even inducing, the lower symmetry of the new electronic states. This sensitivity of physical properties is exemplified in the perovskite ruthenates. The members of the Ruddlesden-Popper series of strontium ruthenate, Sr n+1 Ru n O 3n+1 , exhibit an exceptional variety of ground states ranging from superconductivity for n = 1 [7] via materials with a rich metamagnetic phase diagram for n = 2 [8] and 3 [9] to bulk ferromagnetism for n → ∞. [10] This wide variety of properties is intimately linked to small structural distortions of the RuO 6 cage. Already in the superconductor Sr 2 RuO 4 , with n = 1, a soft phonon mode associated with the octahedral rotation is observed. [11] Isoelectronic substitution of Sr by Ca results in a rich phase diagram with metallic and ferromagnetic phases until an antiferromagnetically ordered Mott insulator is reached, where the dominant change to the material is a small rotation and tilting of the RuO 6 cage. [1,4] Understanding the impact of that rotation on the electronic structure therefore provides a key to controlling electronic and physical properties of perovskite ruthenates. A prerequisite to link physical properties to details of the octahedral rotation is an understanding of the low-energy electronic structure. We here establish the surface layer of Sr 2 RuO 4 as a clean and well-controlled 2D model system where this can be achieved and explored.
The bulk material of Sr 2 RuO 4 has a tetragonal unit cell with undistorted RuO 6 octahedra [12,13] aligned with the highsymmetry directions of the crystal (left panel Figure 1a). It is a superconductor with a transition temperature T c ≈ 1.5K, [7] and is known to exhibit strong electron correlations, evidenced by effective masses of the bands between 6 and 17m e in the normal state. [14] Its bulk Fermi surface has been established by quantum oscillations [15] and confirmed by angle-resolved photoemission spectroscopy (ARPES). [14,16,17] In strongly correlated electron materials, the electronic, spin, and charge degrees of freedom are closely intertwined. This often leads to the stabilization of emergent orders that are highly sensitive to external physical stimuli promising opportunities for technological applications. In perovskite ruthenates, this sensitivity manifests in dramatic changes of the physical properties with subtle structural details of the RuO 6 octahedra, stabilizing enigmatic correlated ground states, from a hotly debated superconducting state via electronic nematicity and metamagnetic quantum criticality to ferromagnetism. Here, it is demonstrated that the rotation of the RuO 6 octahedra in the surface layer of Sr 2 RuO 4 generates new emergent orders not observed in the bulk material. Through atomic-scale spectroscopic characterization of the low-energy electronic states, four van Hove singularities are identified in the vicinity of the Fermi energy. The singularities can be directly linked to intertwined nematic and checkerboard charge order. Tuning of one of these van Hove singularities by magnetic field is demonstrated, suggesting that the surface layer undergoes a Lifshitz transition at a magnetic field of ≈32T. The results establish the surface layer of Sr 2 RuO 4 as an exciting 2D correlated electron system and highlight the opportunities for engineering the low-energy electronic states in these systems.
In Figure 1b, we show its Fermi surface: The d xy band (yellow line) has a van Hove singularity (vHs) at the M-point of the Brillouin zone (BZ) at about 14 meV above E F (Figure 1c). [19] Upon cleaving, the exposed surface layer exhibits a 6° rotation of the RuO 6 octahedra [20] (Figure 1d). In bulk Sr 2 RuO 4 , a similar octahedral rotation occurs on substitution of Sr by Ca quickly suppressing superconductivity. [4] Early STM studies [21,22] suggest that the rotation also suppresses superconductivity in the surface layer. The rotation further leads to orbital-dependent renormalizations close to the Fermi energy. [23,24] This rotation doubles the size of the unit cell and leads to a reconstruction of the Fermi surface (Figure 1d), and a shift of the vHs below E F (Figure 1e). [17] The vHs leads to a peak in the density of states (DOS). In Figure 1c,e) we show schematically the DOS of the unreconstructed (solid line) and reconstructed (dashed line) electronic structure in comparison. The energy of the vHs is seen to be highly sensitive to structural details of the RuO 6 octahedra. Due to the small interlayer coupling in Sr 2 RuO 4 , the surface reconstruction provides an opportunity to establish the influence of small structural distortions on the electronic structure in a very clean and chemically homogeneous 2D system, and hints at what the leading instabilities of the bulk material are. We here demonstrate that the reconstructed surface of Sr 2 RuO 4 by itself constitutes a strongly correlated system with its own unique electronic properties: checkerboard charge order, nematicity and four van Hove singularities (vHss) within a few millielectronvolts of the Fermi energy. We can link these phenomena through a phenomenological tight-binding model based on the bulk electronic structure and incorporating these emergent orders which yields a density of states in excellent qualitative agreement with tunnelling spectra. We demonstrate magnetic-field tuning of one of the vHss, and extrapolate from our measurements that the surface layer undergoes a Lifshitz transition at 32T.

Figure 2a
shows a topographic image of Sr 2 RuO 4 at a bias voltage V = 5 mV, recorded at a temperature of 76 mK, showing a SrO-terminated surface [25,26] with atomic resolution and demonstrating a low concentration of point defects (less than 0.1%). We find defects at the Ru site with two distinct orientations due to the octahedral rotation in the surface layer. The Fourier transformation of the topography (upper inset in Figure 2a), The d xy band exhibits a van Hove singularity (vHs) at the M point above E F , which results in a peak in the density of states (DOS, right). d) Structural model of the surface of Sr 2 RuO 4 with θ = 6°. The rotation is more apparent in the top view shown on the right and leads to a doubling of the unit cell. The unit cell with θ = 0° is shown as a blue square for comparison. e) Fermi surface with θ = 6° (yellow dashed line: d xy band; see Section S4, Supporting Information for details). The surface Brillouin zone is indicated by dotted red lines. f) Electronic structure corresponding to (e). The vHs at the M-point is pushed below E F with the octahedral rotation. The DOS is shown on the right, where the peak occurs below E F .
shows the presence of quasi-particle interference (QPI), consistent with previous reports, [27] as well as Bragg peaks at (0, 2π) and (2π, 0) due to the Sr lattice, and peaks with a low intensity at (π, π) corresponding to the periodicity of the surface reconstruction. When recording topographic STM images at bias voltages of −5 mV, Figure 2b, the appearance changes significantly. The image contrast is now dominated by the presence of a strong modulation of the charge density (lower inset in Figure 2b) which corresponds to a large increase in the intensity of the (π, π) peaks in the Fourier transformation. This modulation has been observed previously at the surface of Sr 2 RuO 4 , [25] as well as at that of Sr 3 Ru 2 O 7 . [28] However, the octahedral rotation due to the surface reconstruction cannot explain the charge modulation: adjacent Sr and Ru sites are in an equivalent structural environment that can be transformed into each other through a symmetry operation-a mirror operation for Ru and a 90° rotation for Sr. Detailed studies of the surface structure by low-energy electron diffraction (LEED) do not reveal any additional reconstruction of the surface layer apart from the octahedral rotation, [29] leaving the origin of the checkerboard charge order as an open question. We note that the checkerboard charge order is most prominent for small bias voltages, suggesting an electronic origin (see Section S2, Supporting Information). A typical differential conductance spectrum g(V) in the range ±95 mV is presented in Figure 2c. Here, kink-and gap-like features are observed at ±40 meV and ±5 meV, respectively. The latter is associated with a reduction of the differential conductance by almost 30% in relation to the value at 95 mV. The general shape of the spectrum, as well as the absence of a superconducting gap, is in agreement with g V R at π π = q ( , ) ckb obtained from maps taken at T = 2 K and 59 mK. Two peaks at −3.5 mV and +3.5 mV with opposite phase and a width of ≈1 mV can be seen (see Section S3, Supporting Information for details). Vertical dotted lines indicate the position of the four peaks observed in the average (r, ) g V spectrum at 59 mK. Insert: Fourier transformation (q, ) g V for the T = 2 K map, q ckb is indicated by a red circle.
previous reports. [21,22,27] However, from high-resolution spectra acquired at temperatures below 100 mK, we find that this gaplike feature actually exhibits four well defined peaks in the differential conductance within ±5 meV of the Fermi energy ( Figure 2d). The additional modulation of the charge density leads to pronounced signatures in spectroscopic maps: while topographic images obtained at positive bias voltage show predominantly the Sr lattice (top panel in Figure 2e), differential conductance maps exhibit a clear checkerboard charge modulation superimposed to this lattice at bias voltages close to E F (second and third panel in Figure 2e). The checkerboard exhibits a contrast inversion across the Fermi energy. We plot in Figure 2f the energy dependence of the intensity of the checkerboard analyzed through a phase-referenced Fourier transformation (PR-FT, see Section S3C, Supporting Information), which provides information about the amplitude and phase of the modulation. The amplitude of the checkerboard exhibits a pronounced maximum at −3.5 mV with a width of only about 1mV and a weaker maximum at +3.5 mV of similar width but with opposite phase. Both maxima occur at the same energy as the energy of the inner-most peaks in the tunneling spectrum (yellow dotted lines in Figure 2f). The phase change is consistent with what one would expect for a charge density wave. [30]

Nematicity
The Sr-centred checkerboard charge order and nematicity are intimately linked in the surface layer through the octahedral rotation. The octahedral rotation itself preserves C 4 symmetry and does not give rise to an additional charge modulation or nematicity. However, the occurrence of a Sr-centred checkerboard, as we observe experimentally, and nematicity are equivalent: due to the checkerboard charge order centred at the Sr sites and the octahedral rotation, the oxygen atoms in the Ru plane become inequivalent between the horizontal, [10], and vertical, [01], directions. As shown in Figure 3a, the rotation means that the oxygen atoms connecting Ru atoms in the [01] (vertical) direction (colored in red) are closer to Sr atoms with a decreased charge density (shown dark), whereas oxygen atoms in the [10] (horizontal) direction (shown orange) are closer to Sr atoms with increased charge density. This leads to different hopping amplitudes across those oxygen atoms, indicated in Figure 3a by the yellow and red dashed lines, resulting in nematicity of the electronic states and a reduction from C 4 symmetry to C 2 symmetry ( Figure 3a). The converse holds true as well, nematicity on the oxygen sites (indicated by different coloured oxygen atoms in Figure 3a) results in their inequivalence and hence checkerboard charge order on the Sr sites -highlighting the equivalence of the two orders. Therefore, while the octahedral Adv. Mater. 2021, 33, 2100593 rotation itself gives rise to neither nematicity nor checkerboard charge order, the two become intimately related through the octahedral rotation-so if one occurs, the other will too. This emergent nematicity is confirmed experimentally through unidirectional modulations with atomic periodicity (Figure 3b,c) and anisotropy of the lowq q quasi-particle interference (Figure 3d-i). The atomic scale symmetry breaking reveals that for changes in the bias voltage V by about 1 mV the atomic-scale unidirectional periodicity changes direction between the highsymmetry directions [10] and [01] implying a small characteristic energy scale of the nematicity. This is also confirmed in the intensity of the atomic peaks ( Figure S5, Supporting Information). The symmetry breaking of long-wavelength quasi-particle interference is shown in a real-space map around four defects in Figure 3d-i (defects marked by black crosses) and reveals a change from predominant quasiparticle scattering along [10] to scattering along [01] and back as a function of energy.
We note that neither nematicity nor checkerboard charge order are captured by the surface structure determined by I(V)-LEED. [29] Nor does I(V)-LEED or our STM data provide any evidence for an orthorhombicity of the surface layer which could give rise to a lowered symmetry.

Tight-Binding Model for Surface Electronic Structure
To understand the microscopic consequences nematicity and the observed charge modulation have on the electronic structure of the surface layer, we have developed a minimal tightbinding model for the band structure at the surface of Sr 2 RuO 4 . This model was generated by taking a tight-binding model of the bulk of Sr 2 RuO 4 , [18] and then altering the position of the d xy vHs such that it is shifted toward the Fermi level (see Section S4, Supporting Information for details). We then include the emergent orders at a phenomenological level. The surface reconstruction and the checkerboard charge order ( Figure 2) are accounted for by including a weak intraband hybridization (Δ hyb ). The nematicity and inequivalence of the [10] and [01] lattice directions ( Figure 3) are included through a nematic order parameter Δ nem = δ nem (cos(k x ) − cos(k y )), which we apply only to the d xy orbital. A similar order parameter has previously been introduced to describe nematicity in Sr 3 Ru 2 O 7 . [31] The full Hamiltonian is then where H ( ) Ru k k is the Hamiltonian for the ruthenium 4d t 2g bands in the unreconstructed Brillouin zone, [18] and π π = ( , ) Q Q accounts for the doubling of the unit cell (for details see Section S4, Supporting Information). The inclusion of the emergent orders leads to the formation of four vHss around the M point, shown in Figure 4a which shows the band structure obtained from Equation (1). We follow the notation by van Hove [32] to label the four vHss: nematicity results in two saddle points, 1 S and 2 S , on the high symmetry axis at the M point, whereas hybridization leads to a partial gap with a band maximum M and an additional saddle point around E F . The resulting complex low-energy electronic structure around the M point is shown in a 3D representation in Figure 4b. The partial drop in the g(V) spectrum around E F , observed in Figure 2d, is naturally explained by the hybridization of the Ru d xy -bands caused by the doubling of the unit cell. The four singularities lead to maxima in the density of states (see Figure 4c) in excellent agreement with the structure of the low-energy g(V) spectrum. Differences, such as the magnitude of the drop, are likely a consequence of tunneling matrix element effects, which are neglected in the calculation.

Magnetic-Field Tuning of the Electronic Structure
The presence of multiple vHss in the immediate vicinity of the Fermi energy offers the opportunity to drive the material through a Lifshitz transition using magnetic field and Zeeman splitting of the bands. We investigate the influence of an applied magnetic field, focusing on the behaviour of the most prominent vHs, M, at −3.4 mV. Figure 4d shows tunneling spectra recorded in magnetic fields from 0 T to 13.4 T, applied perpendicular to the surface, at T = 76 mK. The spectra reveal a clear splitting of the vHs with increasing magnetic field. In Figure 4e, the fielddependence of the energy of the vHs reveals a linear behavior as expected for a Zeeman-like splitting. We find the slope to be +0.10 mV T −1 for the peak moving toward the Fermi energy and −0.07 mV T −1 for the one moving away. The difference in slope can be attributed to an overall chemical potential shift of the d xy band of 17 μV T −1 . We find a g-factor of g* ≈ 3 (with the splitting ΔE = g*μ B B z ). This value is consistent with the known Wilson ratio W = g*/g = 1.5 of bulk Sr 2 RuO 4 , [33,34] a surprisingly good agreement given that the Wilson ratio contains contributions from the whole Fermi surface, whereas we determine g* only for one band. From our experimental data, we can extrapolate that the vHs will cross the Fermi energy at a magnetic field of about 32 T. At this field, the surface layer is expected to undergo a Lifshitz transition. As the energy of the vHs is associated with the checkerboard charge order, once it is split in a magnetic field, the charge order becomes spin-polarized. We show in Figure 4f that indeed also the energy at which the checkerboard charge order appears most prominently shifts with the change in the energy of the vHs, demonstrating that the two are intimately linked. Our measurements demonstrate that the electronic structure of the reconstructed surface layer of Sr 2 RuO 4 has all the ingredients for a field-tuned Lifshitz transition: 1) a vHs close to the Fermi energy, 2) the energy of the vHs can be tuned by magnetic field toward the Fermi energy, and 3) a Lifshitz transition of the electronic structure within reach of available magnetic fields.

Discussion
Our measurements show how the low-energy electronic structure of a 2D layer of Sr 2 RuO 4 is impacted by octahedral rotations and stabilizes new emergent phases close to a magnetic-field induced Lifshitz transition. We provide a comprehensive description of the surface electronic structure through a phenomenological tight-binding model that incorporates the two orders we detect, nematicity and checkerboard charge order and reproduces all features we observe in our data. It captures the pronounced gap-like feature seen in tunneling spectra. The quick changes in the appearance of the nematicity in spectroscopic maps within a narrow energy interval are naturally explained by the vHss in the k x -and k y -directions becoming inequivalent. While checkerboard charge order and nematicity are intimately intertwined, it is not clear which of the two dominates, if either does. Possible mechanisms include 1) nematicity driven by electronic correlations, 2) charge order driven by nesting or a lattice distortion, 3) a cooperative effect of both, and 4) an antiferromagnetic order or fluctuations which stabilize the checkerboard order and nematicity. For the first three scenarios, one would expect a significant structural distortion accompanying these, as is seen for charge density waves or nematicity in other systems. A structural distortion in the surface layer of Sr 2 RuO 4 beyond the octahedral rotation has not been reported. [29] In a scenario where magnetic order or fluctuations drive nematicity and charge order, one would still expect a structural distortion, though much smaller. Recently, evidence for nematicity in thin films of Sr 2 RuO 4 has been reported, suggesting that it may not be limited to the surface layer. [35] A comparison with Sr 3 Ru 2 O 7 reveals intriguing parallels: Sr 3 Ru 2 O 7 exhibits a similar (albeit larger) octahedral rotation, and the energy of the vHs as detected by ARPES is found at ≈4 meV below E F , [36] close to the energy at which we find the dominant vHs in the surface layer of Sr 2 RuO 4 . Due to stronger correlations and a g-factor of g* ≈ 14.6 in Sr 3 Ru 2 O 7 , [37] the vHs can be tuned to the Fermi energy at 8T, whereas in the surface layer of Sr 2 RuO 4 , we find a g-factor that is almost 4 times smaller. Consequently, the vHs is expected to reach the Fermi energy only at 32T. It remains to be seen whether the parallels go any further when the surface layer in Sr 2 RuO 4 undergoes the Lifshitz transition. There are also important differences, though, including that the system we report here is strictly 2D which puts the criticality of the Lifshitz transition into a different universality class than what is expected for the bulk of  reveal a clear shift of a peak as a function of magnetic field, [28] possibly because quantum fluctuations play a much larger role and already have influence on the line shape of the vHs in the density of states. Previous studies of the reconstructed surface have not been able to detect signatures of superconductivity, [21,22] except a recent study [38] where a gap-like structure in the vicinity of the Fermi energy has been attributed to superconductivity. Tunneling spectra showing clear signatures of superconductivity with a temperature and/or magnetic field dependence consistent with bulk Sr 2 RuO 4 [39][40][41] were either obtained from surfaces that were exposed to air leading effectively to a dirty surface much like in ARPES experiments probing the bulk electronic structure or on a different surface reconstruction than the one studied here. We do not observe any evidence for superconductivity even in tunneling spectra acquired at temperatures well below 100 mK. We have ensured that the spectroscopic resolution of our instrument is sufficient to detect superconducting gaps of materials with a similar transition temperature and gap sizes on the order of 200 μeV. [42,43] This indicates that the emergent electronic order suppresses superconductivity in the surface layer. This suppression may imply that the density of states of the d xy band around E F , which becomes gapped out at the surface, and the fluctuations associated with the reconstruction play a crucial factor in superconducting pairing. This scenario naturally results in a competition of nematicity and the charge density modulations with superconductivity, reminiscent to what is found in other strongly correlated electron systems. [5,44] Understanding the leading instability in this highly debated material therefore undoubtedly will have to account for this susceptibility toward density wave formation as observed in other unconventional superconductors.

Conclusions
We show that the surface layer of Sr 2 RuO 4 provides a 2D model system to study the intricate structure-property relationships of a strongly correlated electron system. We demonstrate the equivalence of checkerboard charge order and nematicity in this system, and find that the reconstructed electronic structure leads to four vHss within 5 meV of the Fermi energy. Magnetic-field tuning of one of these vHs implies that the surface layer can be used as a well-controlled test system to study magnetic-field induced Lifshitz transitions, enabling detailed comparison with microscopic theories. Because the emergent surface phase is strictly 2D, limited to the surface layer, all relevant information about the electronic states is accessible spectroscopically. Magnetic-field tuned Lifshitz transitions have been proposed to be at the heart of the quantum critical behavior in a range of heavy fermion and strongly correlated electron materials, [45][46][47] yet the ability to spectroscopically trace the electronic structure across a field-tuned quantum phase transition has remained elusive. Given the sensitivity of the energy of the vHs in Sr 2 RuO 4 to uniaxial strain, [6] we expect that the magnetic field at which the Lifshitz transition is extrapolated to occur here can be reduced substantially by combining uniaxial strain with magnetic field. This creates the opportunity to spectroscopically verify the role of quantum fluctuations across a magnetic field-tuned Lifshitz transition through a detailed study of the line shape of the vHs as it is tuned across the Fermi energy. The sensitivity of the electronic structure of ruthenates to tiny structural modifications found here shows opportunities for tailoring correlated electronic phases in 2D and exploring their physics for novel electronic devices.

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.

S1. MATERIALS AND METHODS
A. Single crystal growth.
The Sr 2 RuO 4 crystals used in this work were grown by the floating-zone technique with Ru self-flux, using a commercial image furnace with double-elliptical mirrors and two 2.0kW halogen lamps (S1). Morphological and elemental characterization of the crystals was carried out using a Zeiss Leo EVO 50 scanning electron microscope (SEM) equipped with an Oxford INCA Energy 300 energy dispersive X-ray spectroscopy (EDS) system. The structure and crystalline quality of the samples were assessed by a high-resolution X-ray diffractometer (Panalytical, X Pert MRD), with a Cu K-α source.

Resistance measurements
Transport measurements were performed using a four probe technique with a 3 He refrigerator. The resistance as a function of temperature (Fig. S2) shows a superconducting transition at T c = 1.5K, in agreement with data reported in literature for good quality crystals(S2).
We determine a residual-resistance ratio (RRR) lim T →0K R(300K) R(T ) ≈ 666, extrapolated to 0K from fitting the resistance above T c to R(T ) = R res + AT 2 ( Figure S2(a)). This value is comparable to that of high-purity crystals reported in the literature(S3).

Scanning tunneling microscopy (STM)
Experiments were performed in a home-built ultra-low temperature STM operating in a dilution refrigerator.(S4) Samples were prepared by in-situ cleaving at low temperatures (∼ 20K) in cryogenic vacuum. We used STM tips cut from PtIr wire, and prepared them in-situ by field emission on a Au(111) single crystal. Bias voltages were applied to the sample, with the tip at virtual ground.
Spectroscopic measurements were performed using a lock-in amplifier to measure differential conductance g(V ). The bias voltage V was modulated at a frequency of ν = 397Hz, the lock-in modulation V L , specified as amplitude, is given in figure captions where applicable.
The bias and current setpoints, V set and I set , are indicated for both topographies and spectroscopic maps.
After cleaving and inserting the samples into the STM head, we have always observed the reconstructed SrO termination as observed in previous reports(S5-S8) and as can be inferred from the appearance of the defects and the additional Fourier peaks due to the checkerboard charge order. sponding to the checkerboard charge order at q ckb = (π, π) for bias voltages between +/ − 70mV.
The inset shows a typical Fourier transformation taken with −5mV, where the peak at q ckb is indicated by a red circle. This data was taken at 2K and the tunneling resistance was kept at 100MΩ for all measurements.

POGRAPHIC STM IMAGES
The appearance of the checkerboard modulation in topographies is bias dependent, as seen in Figures 2(a) and (b) of the main manuscript. Here we show the bias dependence of the appearance of the checkerboard charge order in the voltage range between −70mV and +70mV. We plot the intensityz(q ckb ) of the peak in the Fourier transformation associated with the checkerboard order (q ckb = (π, π)) as a function of applied bias voltage V in Figure S3. The intensity of the checkerboard charge order has a sharp maximum around −5mV, with its intensity decreasing by a factor of 25 at −70mV and of 8 at +70mV. Using this information, for spectroscopic maps shown here, we have chosen a tunneling setpoint at small positive bias voltages where the checkerboard charge order is weak to minimize the setpoint effect.

S3. DIFFERENTIAL CONDUCTANCE MAPS
In this section, we show the real space conductance maps underpinning the data shown in Fig. 3(d) and how it has been processed.

A. Real space maps
The topographic image acquired simultaneously with the map taken at 2K is shown in Figure S4(a) and the corresponding real space differential conductance map, g(r, V ) = dI/dV (r, V ), for V = −3.4mV in Figure S4(b). The checkerboard charge order is clearly visible. In addition quasiparticle interference (QPI) effects are observed around the defect seen in the topography.
B. Processing of differential conductance maps In Figure S4 we show images following each step of the data processing. The raw Fourier transformation of the g(r, V ) map is displayed in Fig. S4( breaking shown in Fig. 4(b) and (c) in the main text and is linked to the nematicity.

C. Phase-referenced Fourier transformation
In figs. 3(d) and 6(c) of the main text, we show phase-referenced Fourier transformations to deduce the characteristic energy scale of the checkerboard charge order and the relative phase at positive and negative bias voltages. The Fourier transformationg(q, V ) of a differential conductance mapg(r, V ) which shows the atomic lattice exhibits two pairs of Bragg peaks at q at = (±2π, 0) and (0, ±2π), indicated by a black circle in Figure S6(a).
The checkerboard order shown in Figure 3 for the data at 2K. The intensity of the peak at q CKB = (π, π) due to the checkerboard charge order is strongly energy dependent. The layer corresponding to the energy of the vHs is highlighted by a red box, where the peaks at q CKB exhibit the highest intensity. Additionally, the intensity of the Bragg peaks q at = (0, ±2π) and (±2π, 0) switches with energy due to the nematicity (cf. fig. 4(b), (c) of the main text). For each layer, the higher intensity peak is highlighted by a red circle, and the lower intensity peak with a dashed white circle. The range of the color bar is the same for all images. √ 2 times larger, with the unit vector 45 • rotated. The periodicity due to the checkerboard charge order shows up as four peaks at q CKB = (±π, ±π) and q CKB = (±π, ∓π), one of which is indicated by the red circle in Fig. S6(a). The amplitude of this checkerboard order is reflected in the intensity of the Fourier peak, |g(q CKB , V )|, as shown in Fig. S5. We can determine the characteristic energy scale of the checkerboard modulation as a function of energy by plotting, |g(q CKB , V )| as a function of applied bias V , however losing the phase information φ(q, V ) contained in the Fourier transformation, where ∆x and ∆y are the size of the map in along the x and y directions and N x and N y are the number of pixels in each direction. While analyzing the phase φ(q, V ) itself is possible, it suffers from an arbitrary global phase factor. To remove this global phase factor, we use a phase-referenced Fourier transformation (PR-FT)(S9) In this PR-FT, the phase at each q-vector is referenced to the phase at a specific energy V 0 , removing the global phase factor. This allows tracking of the change in phase as a function of energy relative to that layer by simply plotting the real part of the PR-FT, Here, we reference the phase to the map layer at V 0 = −3.4mV. The real part of the PR-FT images, Figure S6(c), Re[g R ((π, π), V )] allows us to determine if there is a phase shift between the checkerboard order at different energies as well as the relative amplitude. A phase reversal means a change in sign. At the reference energy, V 0 , Re[g R ((π, π), V 0 ] will be positive by definition. Fig. S6(c) shows Re[g R ((π, π), V = +3mV], where the peaks at q CKB appear with a negative sign, evidencing a phase shift with respect to the charge modulation compared to panel A (map parameters as in fig. S4).

S4. TIGHT-BINDING MODEL
The tight-binding model for the surface of Sr 2 RuO 4 may be defined, in an analogous manner to that of the bulk, via hopping associated with the three bands derived from the ruthenium t 2g orbitals (S10), This places the vHs just above the Fermi level and slightly decreases the Fermi wave vector, k F , of the d xz and d yz bands as suggested by ARPES measurements (S11). We note that the shape of the d xy related pockets observed in both ARPES and DFT are slightly larger than in our model, however this difference does not affect the conclusions drawn here.
We then introduce a hybridisation between the two Ru sites, ∆ hyb = 3meV, as off-diagonal elements in Eq. S7, and include a phenomenological C 4 symmetry breaking term ∆ nem (k) = δ nem (cos(k x ) − cos(k y )) specifically to the d xy orbital in H Ru , with δ nem = 2meV, to produce the full Hamiltonian defined in Eq. 1 of the main text. A similar nematic term has been discussed for Sr 3 Ru 2 O 7 previously.(S12) The density of states presented in Fig. 5(c) has been calculated via N 0 (ω) = − 1 π Tr Im k G(k, ω) .
Here, G(k, ω) is the Green's function defined as G(k, ω) = 1 ω−H(k)+iΓ , where ω is the energy and Γ a broadening parameter. We use a k-grid of 4096x4096 lattice points and Γ = 0.1meV. Fig. S7 shows the Fermi surface, band structure and DOS given by the tight-binding model described above. Fig. S7(a)-(c) show the case where only the doubling of the unit cell is taken into account (Eq. S7). The vHs was put above, but very close to, E F and it appears as a sharp peak with logarithmic divergence, cut off by the broadening parameter Γ. Fig. S7(d)-(f) show the case were only the C 4 -symmetry breaking term is included. The d xy band becomes C 2 -symmetric, and the vHs splits into two peaks, one above E F and another below, although no gap opens around the Fermi energy. Fig. S7(g)-(i) show the case of eq. S9, where both the nematic term and the hybridization potential are included. A gap is opened between the d xy bands, creating four vHs. The calculated DOS reproduces the measured differential conductance spectrum, as shown in Fig. 5 of the main text.
In Fig. S8, we present the DOS calculated using the surface tight-binding model in the presence of a magnetic field. To simulate this, we introduce a Zeeman splitting term to the Hamiltonian Here σ = +1 for up spin states and −1 for down spin states. g * has been set to 3, as determined experimentally and discussed in the main text, µ B = 5.788 · 10 −5 eVT −1 and B z is the magnetic field strength in the z-direction. The introduction of a magnetic field reduces the peak height of the vHs's. line. For the data at other fields, we fit the values of a and c to the background using the same arctangent function, but keep V 0 and Γ fixed at the values of the fit at 13.4T, Figure   S9(a). To determine the energy of the van Hove singularity, we fit a Lorentzian to the peak in the spectrum, Figure S9(b).