Discovery of orbital ordering in Bi2Sr2CaCu2O8+x

The primordial ingredient of cuprate superconductivity is the CuO2 unit cell. Theories usually concentrate on the intra-atom Coulombic interactions dominating the 3d9 and 3d10 configurations of each copper ion. However, if Coulombic interactions also occur between electrons of the 2p6 orbitals of each planar oxygen atom, spontaneous orbital ordering may split their energy levels. This long-predicted intra-unit-cell symmetry breaking should generate an orbitally ordered phase, for which the charge transfer energy ε separating the 2p6 and 3d10 orbitals is distinct for the two oxygen atoms. Here we introduce sublattice-resolved ε(r) imaging to CuO2 studies and discover intra-unit-cell rotational symmetry breaking of ε(r). Spatially, this state is arranged in disordered Ising domains of orthogonally oriented orbital order bounded by dopant ions, and within whose domain walls low-energy electronic quadrupolar two-level systems occur. Overall, these data reveal a Q = 0 orbitally ordered state that splits the oxygen energy levels by ~50 meV, in underdoped CuO2.


1
Unforeseen and unexplained among the quantum matter states of hole-doped CuO2 is an electronic nematic phase 1 . Analyses of the three-orbital model [2][3][4][5][6][7][8][9] for the CuO2 chargetransfer insulator have long predicted a candidate mechanism for the nematic state involving an energy splitting between the two planar-oxygen 2p 6 orbitals within each unit cell.
Although never observed, such effects should be identifiable experimentally as a difference in charge-transfer energy ℰ( ) separating each oxygen 2p 6 orbital from the relevant copper 3d 10 orbital configuration. Moreover, if extant, such ordering of the oxygen 2p 6 orbitals may be compared to iron-based high temperature superconductors where intra-unit-cell rotational symmetry breaking between the iron and orbitals 10 -12 is pivotal 13,14 .

2
Soon after the discovery of Fe-based high temperature superconductivity the profound significance of ordering in the iron and orbitals was identified [10][11][12] , along with the subsequent realization that such orbital ordering (lifting energy degeneracy of and orbitals) was a fundamentally important phenomenon 13,14 . In these materials, such orbital ordering is pivotal for the tetragonal to orthorhombic phase transition into an electronic nematic phase 1,13 . This lifting of : energy degeneracy certainly has profound global effects 13,14 : sequentially with falling temperature, the orbital ordering renders the electronic structure strongly nematic along with a strongly anisotropic antiferromagnetic state; its spectrum of magnetic fluctuations exhibits equivalently reduced symmetry 13 ; finally a strongly anisotropic 14 and even orbital selective 15 form of high temperature superconductivity emerges. Hence, intra-unit-cell orbital ordering produces powerful, wide-ranging effects on the electronic structure, quantum magnetism, high temperature superconductivity and on the global phase diagram of the Fe-based high temperature superconductive materials 13,14, 16 , 17 . However, although long anticipated 2-9 , analogous intra-unit-cell orbital ordering for Cu-based high temperature superconductive materials, has never been observed.

3
In these materials, the planar Cu ++ ions are in the 3d 9 configuration with a singly occupied orbital, while the planar O -ions have filled 2p 6 orbitals. The Cu 3d 10 configuration, which is energetically disfavoured by the Coulomb energy U required to double occupy each orbital, is separated from the intervening oxygen 2p 6 Here , label any of the three orbitals, /ℏ are transition rates for electrons between orbitals , at sites i, j, are the orbital energies, and ↑ , ↓ are the d-orbital occupancies.
High temperature superconductivity is believed to emerge within this model driven by the charge-transfer superexchange pairing interactions between the Cu d-electrons (Methods Section A). In theory 1 , however, many other ordered phases could emerge upon hole doping CuO2, with a planar oxygen orbital-ordered phase 2-9 due to additional inter-oxygen repulsion terms ∑ in Eqn. 1, being a prime example.

4
One enigmatic phase that does emerge is the pseudogap (PG) regime. Its essential phenomenology 20,21 is that, for < * ( ) and < * ≈ 0.19 (yellow shading Fig. 1a), the Fermi surface becomes partially gapped thus diminishing the spectrum of electronic states, the magnetic susceptibility, the electronic specific heat, the c-axis electrical conductivity, and the spatially-averaged density of electronic states. Strikingly, however, in this same region of the phase diagram there is pervasive evidence for a nematic state in which the electronic structure breaks 90 o -rotational (C4) symmetry at wavevector Q=0, or equivalently within the CuO2 unit cell. Experimental signatures of this nematic state have been very widely reported based on multiple techniques as summarized in Fig. 1a, and this nematic phenomenology is observed for < * in the Bi2Sr2CaCu2O8+x, Bi2Sr2CuO6+x, Bi2-zPbzSr2-yLayCuO6+x, YBa2Cu3O7- x, La2-xBaxCuO4, and HgBa2CuO4+δ material families, thus strongly indicative of universality (Methods Section B). But no microscopic mechanism has yet been experimentally established for this CuO2 nematic phase.

5
Classically, a nematic metallic state may occur due to the Pomeranchuk instability of the Fermi surface 22 . But, for cuprates, theoretical analyses using a three-orbital band structure 2-9 predict the emergence of low energy electronic nematicity driven primarily by high energy orbital ordering. Specifically, when the Coulomb repulsions between electrons on nearest-neighbour O sites within the same unit cell are included in the strong-coupling limit 2 and in self-consistent mean-field theory 3 , an orbitally ordered nematic phase is predicted. On the same basis, such CuO2 orbital order is also predicted using diagrammatic perturbation theory 4 , perturbative expansion 5 , functional-renormalization-group techniques 7 , Hartree-Fock based models 8 , and in slave-boson mean-field theory 9 . Figure 1b represents conceptually such planar-oxygen orbital-order while Fig. 1c represents its schematic impact on the intra-unit-cell characteristics of charge-transfer energy ℰ. Here the 2 orbital of the oxygen atom along the CuO2 x-axis (Ox site) is separated from the upper Cu band by the charge transfer energy ℰ , while the notionally equivalent oxygen orbital along the y-axis (Oy site) exhibits a different charge transfer energy ℰ . Our objective is then a direct search for such rotational symmetry breaking at the charge-transfer energy scale, through visualization of ℰ within each CuO2 unit cell.

6
In practice, Bi2Sr2CaCu2O8+x samples with hole-density ≈ 0.17 and = 88 K are inserted into a spectroscopic imaging scanning tunneling microscope (SISTM) and cleaved in cryogenic ultra-high vacuum at T = 4.2 K. The technique of imaging charge transfer energies in cuprates has been established by several STM studies [23][24][25] of variations in the dI/dV spectra measured over the energy range -1.5 eV < E < 2 eV. However, intra-unit-cell resolution imaging of the charge transfer energies have proven challenging due to the high tunneling junction resistance (and thus large tip-sample separation) required for such high voltage imaging. Here we overcome such challenges and introduce sublattice resolved [26][27][28][29] imaging of ℰ( ) to CuO2 studies. A typical topographic image, T(r), of the terminal BiO surface from these studies, is shown in Fig. 2a. To visualize the intra-unit-cell structure of the charge-transfer energy scale, we use a recently developed technique 25 which analyzes high-voltage tunneling conductance spectra ( , ) ≡ / ( , ) so as yield the spatial variations in ℰ( ). Fig. 2b shows an exemplary high-voltage ( , ) spectrum, measured using junction-resistance ≈ 85 GΩ at a specific site within the CuO2 unit cell (yellow dot in Fig. 2a). By comparison to the spatially averaged ( ) in the same field-of-view (FOV) we see that they have different energy separation between the lower band (filled) and the upper band (empty) (compare red and blue arrows in Fig. 2b). Hence, by visualizing ( , ) in the (−1.6 V) ≤ ≤ (2 V) range at these junction-resistances, one can find departures in ℰ( ) from the average, throughout every unit cell in the FOV. To do so, we first use Fourier filtering to remove the structural supermodulation with wavevector from measured ( , ); this is necessary because the supermodulation is known to modulate the charge-transfer energy with an amplitude of approximately 150 meV 25 . We then integrate each measured ( , ) followed by evaluation of ( ) = ∫ ( , ) , the integrated DOS giving rise to the current for all positive energy states and similarly for the negative energy states, ( ) = ∫ ( , ) , where ± are the equivalent integrals for the spatially averaged spectrum ( ) . We normalize the current by the junction resistance, which is given , doing so leads to the energy variations of the charge transfer band. Both filled and empty states are used in this calculation, as explained in Methods Section C. Next to estimate the variations in the charge transfer energy away from its mean we use which averages over a wide range of tunnel conductances g between (0.01 nS) and (0.22 nS). This integration algorithm is demonstrated to provide an improved signalto-noise ratio (Methods Section A & C and Extended Data Fig. 1) in measuring the variations of charge transfer energy compared to the previous algorithm 25 . Finally, the Lawler-Fujita procedure 26 is used to morph the simultaneously measured topograph ( ) into a perfectly periodic image ( ) ( Fig. 2a and Methods Section D) in which the geometry of every unit cell is equivalent; the identical transformation carried out on ℰ ( ) yields an image of charge transfer energy variations ℰ( ) that is equally unit-cell periodic [26][27][28][29] . This ℰ( ) and its power spectral density Fourier transform ℰ( ) become the basis for studies of intraunit-cell symmetry breaking, with Fig. 2c showing the configuration of ℰ( ) measured in the FOV of Fig. 2a. The distribution of measured ℰ values of Fig. 2c is shown in Fig. 2d, and yields an RMS value ℰ ≈ 90 meV . As demonstrated below, this distribution is dominated by the spatial arrangements of an intra-unit-cell ordered state.

7
Next, we explore the symmetry of intra-unit-cell structure of ℰ( ) by studying the values of ℰ( ) measured at the two symmetry-inequivalent Bragg peaks = 2 / (1,0); = 2 / (0,1) of the CuO2 plane 26,27 . Figure  Bi2Sr2CaCu2O8+x. To double check, we also demonstrate directly by intra-unit-cell imaging in real space, that this rotational symmetry breaking is a valid empirical property of the charge transfer gap based in all unprocessed high voltage sublattice resolved dI/dV data (Methods Section I).

8
This motivates development of an intra-unit-cell orbital order parameter for CuO2.
We focus specifically on the two planar oxygen atom sites within each CuO2 unit cell. First, we sample the values of ℰ( ) surrounding oxygen atom sites along the x-axis from Cu atom and label them ℰ − /2 : ℰ + /2 ; equivalently we sample the two nearest planar oxygen atom sites along the y-axis from and label them ℰ − /2 : ℰ + /2 , as shown inset to Fig. 3c. Using this approach, a nematic order parameter ℰ ( ) may be defined 26,30 as the difference between the average ℰ ( ) and ℰ ( ) inside each unit cell: Strong breaking of intra-unit-cell C4 symmetry is now observed in ℰ( ) , while the predominant disorder is revealed as Ising domains of opposite sign (a phenomenology that cannot be the result of anisotropy of the scan tip). The relative preponderance of the two orbital order domains is about 2:1 in these ≈ 0.17 hole-doped samples (Methods Section E). This intra-unit-cell structure in charge-transfer energy implies a redistribution of electric charge that breaks rotational symmetry in the form of a charge quadrupole. But each such charge quadrupole has two possible orientations relative to the local environment, with two different energies ( Fig. 3b): it is thus an electronic quadrupolar two-level system. Figure   3c shows the histogram of measured ℰ ( ) in Fig. 3a, whose RMS value is approximately 25 meV. The histogram of measured | ℰ ( )| < 5 meV within proposed Ising domain walls is presented in Fig. 3d, representing a dense reservoir of thermally active quadrupolar twolevel systems which should undergo rapid thermal fluctuations at finite temperatures. Such two-level systems occupy an area fraction of approximate 25% in this sample (Methods Section F). Another issue is the size of the Ising domains, the largest of which (when averaged over all FOV (Methods Section E)) typically has a size of 30 nm 2 , equivalent to approximately 200 CuO2 unit cells. In the same FOV of these studies we can identify the locations of the interstitial oxygen dopant ions (Methods Section H). When superimposed on the simultaneously measured ℰ ( ) images (Fig. 3e), we find a propensity for these ions to lie in the domain walls (Fig. 3f) between the orbitally ordered regions (Methods Section H).

9
To epitomize the internal electronic structure of orbitally ordered CuO2, from Fig. 3a we select all regions with ℰ > +5 meV (domain ) and all regions with ℰ < −5 meV

10
These data may be compared to pioneering X-ray scattering studies of orbital ordering across the La2-xBaxCuO4 materials family. In these crystals, signatures of intra-unitcell symmetry breaking should occur as anisotropy in the scattering tensor at planar Cu sites.
Here, it is the intensity of the (001) scattering peak at the Cu L-edge, which directly measures the electronic nematicity of the Cu 3d states within a single CuO2 plane, and provides the indications of intra-unit-cell nematicity 31  inside the unit-cell due to the density wave, average to 0 over one period 26,27 . However, higher harmonics of a CDW can, in principle, produce composite order parameters which break symmetry at Q = 0, as can "vestigial" CDW states in the presence of disorder 32 . Another important subject is the effect of electron-lattice interactions: while the predominant theoretical motivation for such a state has been based on the repulsive potential between the oxygen atoms 2-9 (Fig. 1b), it is inevitable that such orbital order will also couple to the B1g phonon mode. If this mode softens in Bi2Sr2CaCu2O8+x, the O atom elevated above the CuO2 plane experiences a different potential than its neighbor which is depressed below it 33 .
One surprising new revelation is the distribution of spatially localized, low activation energy two-level systems (TLS) within the boundaries between orbital ordered domains. These are the electronic quadrupoles within each unit cell (Fig. 4c, d) whose energy barrier between the two orientations of ℰ tends to zero (

11
Nevertheless, significant consequences do emerge from these studies. First, when physically realistic models of CuO2 electronic structure are used (e.g. Eqn. 1 with parameters relevant to real materials), Coulomb interactions between electrons at the two crystalequivalent oxygen sites should generate intra-unit-cell rotational symmetry breaking i.e.
Spatially, the state is arranged in Ising domains, which appear to be bounded by randomly sited oxygen dopant-ions. Moreover, within domain walls as ℰ → 0 (Fig. 3d) we find an ensemble of low energy-barrier electronic quadrupolar two-level systems. Most fundamentally, the microscopic mechanism proposed theoretically 2-9 for the cuprate nematic phase i.e. orbital order between oxygen orbitals at the two separate oxygen sites of CuO2, is highly consistent with the observed intra-unit cell rotational symmetry breaking of ℰ( ) which splits the energy between the two oxygen atoms by ∼ 50 meV in Bi2Sr2CaCu2O8+x.
Cu 3d Averaged spectrum Typical spectrum  was -750 mV at 25 pA. The bulk incommensurate crystal supermodulation is seen clearly in Fig. 2a. It is at 45˚ to, and therefore is the mirror plane between, the x-and y-axes as always 26 .
The orbital ordering domains are not influenced by supermodulation for this symmetry reason. Furthermore, we compare three independent orbital ordering domains analyzed with and without the supermodulation in Methods Section K. Despite crystal structure strongly modulated by the supermodulation, the orbital ordering domains remain consistent. Therefore the supermodulation has no discernable influence on the intra-unit-cell symmetry breaking of ℰ( ) discussed in this paper. (b) Typical high-voltage differential conductance spectra ( , ) from (a) show as a solid black curve while the spatially averaged spectrum ( ) is shown as a dashed curve. The exemplary spectrum is measured at location a denoted as a yellow dot in (a). Such high junction resistances of 85 GΩ or large tip-sample distances preclude effects on g(V) of the tip-sample electric field. Separation between the lower band and the upper band is clearly very distinct for the exemplary spectrum (blue arrows) and for the average spectrum (red double-headed arrow). We estimate that energy difference ℰ using Eqn. 2. Setpoint for the differential conductance map ( , ) was 600 mV at 7 pA. The peaks are removed in ℰ( ). (h) The linecuts from = (0, 0) to (1, 0)2 / and from = (0, 0) to (0, 1)2 / in the ℰ( ) from (g). The power spectral density ℰ( ) breaks C4 symmetry at its Bragg peaks because the plots of ℰ( ) show anisotropy at = (1, 0)2 / and = (0, 1)2 / . This is the first direct evidence of intra-unit-cell rotational symmetry breaking at the charge transfer energy scale in cuprates.     Fig. 3a where ℰ > +5 meV. (b) Unit-cell averaged structure of ( ) averaged over all regions in Fig. 3a where ℰ < −5 meV.
In both (a) and (b) the C4 symmetry is preserved. (c) Unit-cell averaged structure of ℰ( ) averaged over all regions in Fig. 3a where ℰ > +5 meV. The charge transfer energy strongly breaks C4 rotation symmetry about every Cu site (yellow dots) and consequently there is an energy splitting of approximately 50 meV between the charge transfer energies at the two crystal-equivalent oxygen sites (indicated by the crosses). Unprocessed high voltage dI/dV images reveal directly the subtending data (Methods Section I). Correspondence to: jcseamusdavis@gmail.com

Methods
The spectroscopic imaging scanning tunneling microscope (SISTM) measurements of Bi2Sr2CaCu2O8+x crystals are carried out in a custom-built SISTM. The sample is grown using the floating zone method with doping controlled by oxygen depletion. The hole doping level studied in this paper is p~17%. The samples are cleaved in a cryogenic ultrahigh vacuum at a temperature of ~4.2 K so that the atomically flat BiO termination is revealed. The samples are then immediately inserted into the STM head. All measurements are carried out using tungsten tips at a base temperature of ~4.2 K.

A. Visualization of Spatially Modulating Charge Transfer Energy
The spatially resolved charge transfer energy ℰ( ) may be visualized in cuprates using high voltage 1-3 differential tunnel conductance ( , ) imaging of Bi2Sr2CaCu2O8+x 4 . Differential conductance ( , ) was visualized in the range −1.6 V ≤ ≤ 2 V, and at the extremely high junction resistance ≈ 85 GΩ ( = 600 mV: = 7 pA) necessary to suppress all tip-induced electric field effects. Fig. 2b in the main text shows a typical example of a ( , ) spectrum measured at T = 4.2 K 4 . The edges of the filled lower band and the empty upper band can be identified from the appearance of extremely rapid increase in the density of states. Here the value of the charge transfer energy ℰ( ) is estimated at every location by subtracting these band edges at a constant differential conductance ≈ 20 pS as follows: Although charge transfer measurements have been successful in STM studies, the real space distribution nor atomic resolution has yet been achieved. The difficulty here is maintaining a stable tunneling junction at high bias requires an unusually high junction resistance. In this paper, we achieve intra unit cell resolution in the charge transfer gap at such high junction resistance and thus resolve variations of the charge transfer energy on the Ox and Oy sites.
The crystal supermodulation is a quasi-periodic lattice modulation along the (1,1) direction in the crystal structure at the wavevector . The resulting ℰ( ) map for Bi2Sr2CaCu2O8+x shows strong modulations at the same wavevector, with an amplitude of ≈ 0.3 eV and has a spatial average value of ⟨ℰ( )⟩ ≈ 1.2 eV which agrees well with the charge-transfer energy of Bi2Sr2CaCu2O8+x as measured independently using a variety of different experimental techniques 5,7,8 . Analysis of the supermodulation influence on ℰ( ) and the corresponding supermodulation influence on the electron-pair density ( ) measured by scanned Josephson tunnelling microscopy have revealed that the charge-transfer energy ℰ( ) controls the density of condensed electron-pairs ( ) as / ℰ ≈ −0.81 ± 0.17 eV in Ref. 4. This is in excellent quantitative agreement with cluster Dynamical mean field theory (cDMFT) solutions of the Emery three-band model which predict / ℰ = 0.93 ± 0.1 eV for Bi2Sr2CaCu2O8+x 9 .
The quantitative agreement between theory and experiment here serves a dual purpose. Firstly, it provides strong evidence that charge-transfer superexchange is the electronpairing mechanism of the cuprate high temperature superconductors and secondly, it provides further justification of the experimental ℰ( ) visualization technique described above and used throughout this study. Hence, Extended Data Fig. 1 compares the images of charge transfer energy ℰ( ) derived from the simple algorithm in Eqn. M1 and the more general algorithm for charge transfer energy variation ℰ( ) used in the main text. The images of ℰ and ℰ bear high similarities in real space and reciprocal space. Their linecuts show identical and strong symmetry breaking at the Bragg peaks where the ratio of the intensity of the Bragg peak in the y-direction to the Bragg peak in the x-direction reveals, = 1.9. Thus, the general algorithm calculating the charge transfer energy variation ℰ( ) throughout this paper is consistent with the earlier simple algorithm 4 for calculating ℰ( ). Figs. 1 d, e show the histograms of the distribution of the charge transfer energy ℰ( ) and variation of the charge transfer energy ℰ( ) respectively. Clearly, the ℰ( ) histogram features a narrower σ and therefore features a higher signal-to-noise ratio. Thus, we implement the algorithm to measure charge transfer energy variation ℰ( ) throughout our paper due to its enhanced signal-to-noise ratio.

B. Experimental evidence of Q=0 rotational symmetry breaking in cuprates.
The existence of a nematic state in cuprates, which breaks rotational symmetry at Q = 0 has been seen experimentally. This state has been widely reported based on multiple techniques, (Main Text Fig. 1a). Neutron scattering (NS) is the most widely used [10][11][12][13][14][15] , here polarized neutrons are used to detect the anisotropy at the Bragg peaks, indicating IUC symmetry breaking. Since translational symmetry is found to be preserved, this can be explained by symmetry breaking inside each unit cell. Other experimental methods used are angleresolved photoemission spectroscopy (ARPES) [16][17][18] , where above T* left-circularly polarized (LCP) and right-circularly polarized (RCP) light have the same intensity at the mirror plane, but when cooled below T*, a difference in LCP and RCP intensity at the mirror plane emerges, which is a signature of a = 0 rotational symmetry breaking state. Torque-magnetometry (MT) studies 19,20 have reported a kink of the in-plane anisotropy of the susceptibility at T* which relates to rotational symmetry breaking. Optical anisotropy (ORA) 21 measurements show a significant change in the second-harmonic of the optical response near T*, whereas the linear response remains unchanged in the same region, this could be explained if bulk inversion symmetry was broken, alongside this nematicity is also seen. In polarization resolved Raman scattering 22 suppression of the susceptibility near T* is observed. Elastoresistance (Transport) measurements 23 of in-plane anisotropy which onsets near T* also indicates a nematic state. Resonant ultrasound spectroscopy (RUS) 24 finds a phase transition at T* by noting discontinuities in the frequencies and widths of the vibrational normal modes of a crystal. STM experiments 25 have detected the intra-unit-cell rotational symmetry breaking representing nematicity in the low energy density of states, and this nematicity exhibits Ising domains which diminish in size and intensity approaching the pseudogap endpoint.
The motivation to search for a charge transfer energy splitting between the Ox and Oy site is well predicted by leading researchers who report, for example in Kivelson et al at the strong coupling limit, three-band Emery model produces a nematic phase 31 or Fischer et al study the charge transfer energy using the three-orbital model for the cuprates 32 ; or Hardy et al report a nematic metal phase which arises from the spontaneous partial orbital polarization of the multiorbital non-Fermi liquid 33 .
This well-established theory is the motivation for us to search for real space, cuprate orbital ordering in Bi2Sr2CaCu2O8+x using STM. We have atomically resolved the orbital ordering in Bi2Sr2CaCu2O8+x and confirmed that the orbital asymmetry is located on the oxygen sites.

C. Sublattice-resolved Charge Transfer Energy Visualization
The technical details of measuring the variation of the charge transfer energy ℰ( ) are covered in this section. The averaged differential conductance ( ) is shown as a dashed curve in Extended Data V. The variation in each integral ( ); ( ) from the average values, ; ± , occurs due to the variation in energy separation ℰ( ) between the lower band and the upper band from its average value. To evaluate the energy splitting between the point spectrum and the averaged spectrum, the integral difference is normalized by the difference between the maximum differential conductance and the minimum difference conductance . = 0.22 nS is given by the maximum value of the field of view averaged spectrum and = 0.009 nS is the minimum differential conductance of the averaged spectrum.  Fig. 2). In Case 1 the separation in energy between lower and upper bands in the point spectrum is below the total field of view averaged. The area integrated below the averaged spectrum is larger than the area below the point spectrum thereby the energy variation is positive ℰ( )>0. In Case 2 the separation in energy between lower and upper bands in the point spectrum is above that of the averaged spectrum. The area integrated below the averaged spectrum is smaller than the area below the point spectrum and this case gives rise to a negative energy variation ℰ( )<0.
The choice of the integration range of conductances in Eq. M2 does not alter the conclusion that C4 symmetry is broken in ℰ. Extended Data Fig. 3

D. Lawler-Fujita Symmetrization of Charge Transfer Energy Images
Here we describe how to process the image in removing a pico-meter distortion due to the piezo-electronic drift of the STM tip scanner. First, we apply the affine transformations to the topograph, T'( ) and the simultaneously taken electronic structure images, such that the two sets of Bragg wavevectors of Bi atoms in BiO plane, and in T'( ) satisfy | | = | | and • = 0 after the transformation. Then, we apply the Lawler-Fujita (LF) algorithm 34 to remove the STM tip scanner drift in the data, in which the following steps are taken. The T'( ), in which the STM tip drift is embedded, is given by where ( ) is an amplitude of the complex field ( ) for the wavevector and ( ) is the displacement field associated with the STM tip scanner drift that distorts the T'( ) 34 . To get the displacement field ( ) for a periodic modulation at , the following equations of the two-dimensional lock-in technique are evaluated, where  is a coarse-graining length and Φ ( ) is a spatial phase shift of the modulation at . In this study, we use  = 5. The Lawler-Fujita algorithm is key for visualizing intra-unit-cell electronic structure. To obtain this goal, the dataset with atomic resolution and the dataset with high voltage (-1.6 V to 2 V) must be registered with atomic accuracy. The two datasets are measured in the same field of view. The topographs of the two datasets, ′( ) in the atomic resolution dataset and ′( ) in the high voltage dataset, are used to register the identical atoms in one dataset with those in the other. The Lawler-Fujita algorithm is applied to both datasets to remove the lattice distortions of the image introduced by the distortions of the piezo-tube, and to produce the corrected images of ( ) and ( ) . Subsequently ( ) and ( ) are registered by spatial translations, whose accuracy is evaluated by the cross-correlation between ( ) and ( ). The spatial translations correct any small difference in field-ofview between the measurements. ( ) and ( ) are now registered. All transformation parameters applied to ′( ) are subsequently applied to the high-voltage differential conductance map ′( , ) that is measured simultaneously with the topography. The differential conductance map is ( , ) registered.
Extended Data Fig. 5 shows two topographs of a Bi2Sr2CaCu2O8+x surface taken with atomic resolution and high voltage. They are corrected using the Lawler-Fujita algorithm and registered with atomic resolution. The precision of the image registration is shown in Extended Data Fig. 5e. The maximum of the cross-correlation between ( ) and ( ) coincides with the (0,0) cross-correlation vector. The offset of the two registered images are within three pixels. Taking an uncertainty of half the size of the three pixels provides the precision of the registration method is better than 80 pm everywhere in the whole field-ofview.
To demonstrate that evidence of nematicity is present in the unprocessed charge transfer energy variation ℰ′( ), Extended Data Fig. 6 compares the unprocessed data of ℰ′( ) and the processed ℰ( ). ℰ′( ) is calculated from the unprocessed ′( , ) dataset and ℰ( ) is calculated from the LF-corrected and registered ( , ) dataset. The two images show many features in common. The Fourier transform ℰ′( ) of the unprocessed data shows anisotropy at the Bragg peaks. The processed data ℰ( ) shows the same anisotropic Bragg peaks and the background noise is much lower in the drift-corrected data than in the unprocessed data. Thus, the LF algorithm removes the pico-meter distortion due to the STM tip scanner and increases the signal-to-noise ratio. The LF algorithm does not alter the conclusion that C4 symmetry is broken in ℰ.

E. Comparison Intra-unit-cell Charge Transfer Energy Images from Multiple FOVs
Here we compare the result of multiple experiments to evaluate the repeatability of the phenomena reported. Firstly, the Fourier analysis was repeated on three different field-ofviews of the same sample. The results of this analysis are shown in Extended Data Figs. 7-12. The topograph of each FOV is shown in Extended Data Figs. 7, 9, 11 a, the supermodulation can be clearly identified as expected for Bi2Sr2CaCu2O8+x with ≈ 0.17 . The ( ) is an atomically-resolved topograph that is registered with the differential conductance map ( , ). Extended Data Figs. 7, 9, 11 b then show the variations of the charge transfer energy, ℰ( ). All FOVs show high levels of disorder in these images, and also all span the same energy ranges of -300 meV to +200 meV. The power spectral density (PSD) Fourier transform of the topograph is then shown in Extended Data Figs. 7, 9, 11 c. Highlighted in the blue dashed circle is the Bragg peaks and in the red dashed circle is . The Bragg peaks are of similar magnitude in the topograph, therefore C4 symmetry is preserved indicating that neither the lattice nor the scanning tip produces IUC symmetry breaking. Long-range disorder is seen at = 0. ℰ( ), the PSD Fourier transform of charge transfer energy variation, is then shown in Extended Data Figs. 7, 9, 11 d, again with the Bragg peaks highlighted with dashed circles. In all FOVs it is clear the Bragg peaks at (0, 1)2 / are more intense than the peaks at (1, 0)2 / , meaning IUC C4 rotational symmetry is broken. Linecuts of ( ) are then taken in Extended Data Figs. 7, 9, 11 e. The red linecut is taken from = (0, 0) to (1, 0)2 / and the blue linecut is taken from = (0, 0) to (0, 1)2 / . These linecuts allow quantitative analysis of Bragg peaks. In all FOVs it is clear ≈ in ( ). Confirming that the tip and lattice preserve C4 symmetry. Extended Data Figs. 7, 9, 11 f are linecuts of ℰ( ), again taken from = (0, 0) to (1, 0)2 / and the blue linecut is taken from = (0, 0) to (0, 1)2 / . In all FOVs, there is a clear anisotropy between the intensities of the Bragg peaks. Therefore, IUC rotational symmetry is broken in ℰ. The inset here is a histogram of the distribution of ℰ( ). The histograms are asymmetric and the peak of the histogram is shifted from 0 in all FOVs.
To sum up Extended Data Figs. 7,9,11: the anisotropy in the Bragg peaks of the charge transfer energy variations ℰ are repeatable in multiple experiments at the same hole density. The ℰ from independent FOVs show similar statistics. The nematicity is not generated by the crystallography nor the scanning tip. Therefore, the Fourier analysis of ℰ appears robust and reliable.
Subsequently the IUC charge transfer order parameter and the oxygen site-specific imaging of the charge transfer are repeated for the three field-of-views in Extended Data Figs. 8, 10, 12, respectively. Extended Data Figs. 8, 10, 12 a report the nematic domains for the data in Extended Data Figs. 7, 9, 11, respectively. ℰ ( ) is the intra-unit-cell d-symmetry breaking of the charge transfer energy variations and its calculation details of ℰ ( ) are discussed in the following Section F. The dashed circle shows the r-space radius of the gaussian smoothing used in calculating ℰ ( ). The three FOVs show similar disordered Ising nematic domains. The relative strength between the two nematic domains is defined as the area ratio between the two Ising domains / . The area ratio / is approximately 2.1 ± 0.2, which is consistent with the relative strength between the two Bragg peaks in the power spectral density Fourier transform of ℰ, ℰ ℰ( )~1 .9 ± 0.3. This result demonstrates that the relative preponderance of the two orbital order domains is 2:1 in the 17% hole-doped sample.
Extended Data Figs. 8, 10, 12 f display the oxygen site-specific nematic order parameter ℰ ( ) images (Equation 3 in the Main Text). The spatial distribution of ℰ ( ) is highly similar to that of ℰ ( ) with a cross-correlation coefficient of 0.93. This focuses attention on the charge transfer symmetry breaking mainly occurring on the planar oxygen sites.
The microscopic structure inside the Ising domains is visualized from the unit-cell-averaged images of each domain, as shown in Extended Data Fig. 8, 10, 12  To determine the fraction of the total area occupied by the two-level systems (TLS) in the sample, and to validate the estimated domain wall cut-off at | ℰ | = 5 meV, we calculate the domain boundary trajectory and broaden it by the correlation length. The correlation length of ℰ ( ) is used to estimate the distance from the domain walls where these TLS are likely to be thermally active. First, the autocorrelation of ℰ ( ) is given by where R is the location being integrated at, while r is varied. This is shown in Extended Data Fig. 14a. The decay length of this autocorrelation function is defined as the correlation length ξ. This can be determined by plotting the autocorrelation as a function of radial angle and then fitting, where C is a fitting constant  Fig. 14d), and the area fraction of approximately 25% TLS is found inside these domain walls.
In the main text the domain walls are defined as | ℰ ( )| < 5 meV (the white area defined in Fig. 3d). This leads to an area fraction of 25%, consistent with the method used above. Thus, the cut-off of | ℰ ( )| < 5 meV is mathematically validated using the above algorithm.
The unit cells inside one Ising domain are averaged to increase the signal-to-noise ratio of the microstructure of the IUC symmetry breaking order 35 . Firstly, each CuO2 unit cell is identified in the topography ( ) (gray grid in Extended Data Fig. 15a). Each unit cell has a size of 10×10 pixels and is labelled as ( , , ). The image series ( , , ) is categorized into two datasets ( , , ) and ( , , ). Where I1 is inside a positive domain with ( ) > 5 and I2 is inside a negative domain with ℰ ( ) < −5 ( Fig. Extended Data 15b). The images , ( , , ) have lateral dimensions x and y and a counting index n. ( ) has already been atomically registered and the lattice displacements have been corrected using the Lawler-Fujita algorithm, making each lattice square and periodic. The image set of , ( , , ) ) is made up with n images of , ( , ). The images are cropped independently from ( ), where each image has a size of 11 × 11 pixels (5.9 × 5.9 Å) and corresponds to one CuO2 unit cell rotated by 45°. The single CuO2 unit cell is defined by the location of Bi/Cu atoms (Extended Data Fig. 15). Here the sublattice is a product of two cosine window functions because the product has a more localized focus on the sublattice than the sum of two cosine window functions. The negative values in the window functions ( * ± ) * and ( * ) * ± are set to 0. Like the real-space construction, the oxygen sublattices are displaced by half a unit cell and one oxygen atom is set as the origin of this sublattice. The sublattice decomposition from the reciprocal space rules out the error of identifying the position of the atomic sites in the real space.

G. Oxygen-site Specific Charge Transfer Order: Real Space and Reciprocal Space
The real space nematic order parameter can now be constructed Similarly, the nematic order parameter is calculated from the reciprocal space analysis The images of the segregated Ox and Oy sublattices and , the discrete images of the nematic order parameter ℰ ( ) image and the continuous ℰ ( ) are presented in Extended Data Fig. 16. The sublattice decomposition calculated from real space and from reciprocal space are almost identical, which validates the sublattice decomposition algorithm.
While the order parameters ℰ ( ) and ℰ ( ) originate from two slightly different physical definitions, they are clearly similar visually, as expected. The cross-correlation coefficient of these two images is 0.93 (where 1 means the two images are identical) this confirms that both methods are able to visualize these orbital ordering Ising domains.

H. Dopant Oxygen Ion Pinning of Orbital Ordered Domains
To explore causes for the obvious heterogeneity on the orbital order domains, we search for the dopant oxygen ions and their relationship with the domains. Each dopant oxygen ion is identified as a maximum located at -900 meV in the differential conductance spectrum (inset of Extended Data Fig. 17a). A total of ~70 oxygen dopants are found in the differential conductance map ( , −900 meV) with an image size of 17 × 17 nm 2 (Extended Data Fig.  17a). The locations of the oxygen dopant are overlaid onto the simultaneously measured IUC oxygen-specific nematic order parameter ℰ ( ) in Extended Data Fig. 17d. Visually we see that the oxygen dopants are near the ℰ domain walls (yellow contours).
The distance between each oxygen dopant and their nearest location on the domain walls is calculated and a distribution is shown in Extended Data Fig. 17g. To validate that the oxygen dopants are located near the ℰ domain walls, the distribution of is compared to the distance between randomly generated points and the domain walls. The random points are generated from a two-dimensional Poisson disc sampling function and the random points are separated from each other. The random points have no spatial correlation with the orbital ordering domains. The expected averaged distance to the domain walls is calculated and compared to in Extended Data Fig. 17g.
The distribution of is different from from two aspects. While distribution has a sharp peak at 1.6 Å, distribution has a blunt plateau. The deviation in the distance distribution indicates clearly that the oxygen dopants are located near the domain walls of the orbital order domains, providing statistical evidence that they are pinning the nematic domains.
We repeated the above statistical analysis for three independent FOVs (Extended Data Figs. 17a-c). Extended Data Fig. 17f is chosen to be presented in the main text Fig. 3e. The sum of the three histograms (Extended Data Fig. 17g-i) are presented in the main text Fig. 3f. A total of 237 oxygen dopants are studied in the total histogram.
The field of view used in the main text contains 70 dopant atoms. The dopant radius is quantified as the spatial decay length of the dopant resonance 1 away from the dopant site. We take spatially averaged linecuts over typical dopant maxima and fit with an exponential decay function. The average decay length gives us a dopant size of approximately 2 Å. A ratio of dopant area to the total field-of-view area can be calculated to be on the order of 10%. Therefore, dopant atoms do not make a significant difference to overall ℰ( ) imaging procedures.

I. Direct visualization of orbital ordering from original data
Here we demonstrate directly that the rotational symmetry breaking we report is a true property of the charge transfer gap from unprocessed data (Extended Data Figs. [18][19][20]. The repeatability of charge transfer energy splitting in three additional experiments using a different setup condition (-750 mV, 25 pA) in the unprocessed data demonstrates that the evidence of orbital ordering is independent of setup conditions and data analysis, and orbital ordering is a robust experimental result.
The two datasets are complementary, but it's crucial to emphasize the indispensable role of the original dataset in measuring charge transfer energy variations. Measured at 85 GΩ and spanning a wide energy range of -1.6 V to 2 V, this dataset provides a precise and rigorous mapping of charge transfer energy variations. This is a direct indicator of orbital ordering, a critical aspect that underpins the entire study.
This independent experiment robustly and repeatedly shows the -50 meV ~ -30 meV energy splitting between Ox and Oy sites within the unit cell in the -750 mV and 25 pA setpoint experiments, all in unprocessed data. Thus, in Bi2Sr2CaCu2O8+x intra-unit cell rotational symmetry breaking at the charge transfer energy scale definitely does occur on the intraunit-cell oxygen sites, and this effect occurs oppositely in two Ising domains.

J. Sample Quality Characterization
It is important to consider sample quality, to ensure a high-quality measurement was taken. The single crystal of Bi2Sr2CaCu2O8+x sample studied throughout this study is shown in Extended Data Fig. 22a. The sample is 17% hole-doped. This is an isotope substituted Bi2Sr2CaCu2O8+x such that the oxygens ( 16 O) in the sample are almost fully substituted by 18 O. The sample quality is characterized by a magnetic susceptibility measurement M|H versus temperature across the superconducting transition, shown in Extended Data Fig. 22b. M is the magnetization and H is magnetic field strength. Here, the samples magnetic susceptibility was measured from 50 to 110 K, with a sharp transition seen at 88 K, whose behavior indicates a superconducting transition. Therefore, from this data it is clear a high-quality sample of Bi2Sr2CaCu2O8+x was studied.

K. Orbital ordering domains are unaffected by supermodulation.
To further epitomize the lack of effect the supermodulation has on Q = 0 orbital ordering, we have compared the topographs and the orbital ordering domains obtained with and without the supermodulation from three independent FOVs (Extended Data Fig. 23). The three FOVs are randomly chosen from Extended Data Fig. 21. The first column (Extended Data Fig. 23a, e, i) shows the topographs where the supermodulation is removed via Fourier filtering at the QSM peak. Their orbital ordering domains are presented in the third column (Extended Data Fig. 23c, g, k). The second column presents the topographs where the supermodulation is kept in the analysis (Extended Data Fig. 23b, f, j). Although the atomic structure is strongly modulated by the supermodulation, their orbital ordering domains in the fourth column (Extended Data Fig. 23d, h, l) remain virtually the same as the orbital ordering domains where the QSM peak is filtered (the third column of Extended Data Fig. 23).
Our experimental results align well with our expectations as illustrated in Main Text Fig. 2a, where the bulk incommensurate crystal supermodulation is clearly visible. This supermodulation is oriented at 45˚, serving as the mirror plane between the x and y axes, consistent with previous findings. Importantly, this supermodulation does not impact the intra-unit-cell symmetry breaking of ℰ( ) discussed in this paper, owing to the specific symmetry involved 34 .
In conclusion the orbital ordering domains are independent from the supermodulation. The orbital ordering domains are not a vestige of the crystal supermodulation but a true property of Bi2Sr2CaCu2O8+x. The narrower distribution in ℰ( ) demonstrates an improvement in signal-to-noise ratio over the ℰ( ) algorithm. From ED. Fig. 1, the histogram of ℰ( ) has a narrower σ than the histogram of ℰ( ), demonstrating the signal-to-noise ratio is higher in the algorithm used here to measure the variations in charge transfer energy ℰ( ). Thus this ℰ( ) algorithm with higher SNR was used throughout this paper.  . Clearly these point spectra show the Ox spectra are shifted by -50 meV~-30 meV with respect to their intra-unit-cell Oy spectra. The locations of the Ox and Oy sites are shown as circles in (a-b). (d) 10 point spectra sampled on the oxygen sites from Domain 2 reveal that the Oy spectra are shifted by -50 meV~-30 meV with respect to the Ox spectra from the same unit cell. The locations of the unit cell and oxygen sites are shown in circles in (a-b). The dI/dV splitting to opposite directions in the two different domains, proves charge transfer energy splitting occurs in unprocessed differential conductance spectra. Figure 19. Intra-unit-cell splitting of unprocessed differential conductance spectra (FOV2). (a) Topograph of the FOV where the spectra are measured. (b) Image of orbital-order intra-unit-cell order parameter ( ) sampled on oxygen sites. Two Ising domains are labelled. (c) 10 dI/dV point spectra sampled on the oxygen sites from Domain 1 in (b) reveal that the Ox spectra are shifted by -50 meV~-30 meV with respect to the Oy spectra from the same unit cell. (d) 10 point spectra sampled on the oxygen sites from Domain 2 reveal that the Oy spectra are shifted by -50 meV~-30 meV with respect to the Ox spectra from the same unit cell.    . (b, f, j) present the topographs where the supermodulation is kept. (c, g, k) show the orbital ordering domains calculated from the topographs in the first column where the supermodulation is removed. (d, h, l) show the orbital ordering domains calculated from the topographs in the second columns where the supermodulation remains. The orbital ordering domains from the third (without supermodulation) and fourth columns (with supermodulation) are virtually identical. The supermodulation has virtually no effect on the orbitalorder parameter used in this paper.