Interplay of hidden orbital order and superconductivity in CeCoIn5

Visualizing atomic-orbital degrees of freedom is a frontier challenge in scanned microscopy. Some types of orbital order are virtually imperceptible to normal scattering techniques because they do not reduce the overall crystal lattice symmetry. A good example is dxz/dyz (π,π) orbital order in tetragonal lattices. For enhanced detectability, here we consider the quasiparticle scattering interference (QPI) signature of such (π,π) orbital order in both normal and superconducting phases. The theory reveals that sublattice-specific QPI signatures generated by the orbital order should emerge strongly in the superconducting phase. Sublattice-resolved QPI visualization in superconducting CeCoIn5 then reveals two orthogonal QPI patterns at lattice-substitutional impurity atoms. We analyze the energy dependence of these two orthogonal QPI patterns and find the intensity peaked near E = 0, as predicted when such (π,π) orbital order is intertwined with d-wave superconductivity. Sublattice-resolved superconductive QPI techniques thus represent a new approach for study of hidden orbital order.

In a crystalline metal, strong electronic correlations occurring between electrons derived from different orbitals in the same atom can yield an orbital-selective Hund's metal state 1,2 , or even orbital-selective superconductivity 3,4,5 . Similarly, symmetry breaking orbital order may occur, with one of the most famous cases being in the Fe-based high temperature superconductors 6,7 . However, some types of orbital order are almost indiscernible because they do not occur with any lattice distortion which reduces the overall crystal lattice symmetry. For example, () orbital order in a tetragonal array of transition-metal atoms occurs when the degeneracy of dxz and dyz orbitals is lifted and each predominates energetically over the other at alternating lattice sites (Fig. 1a). This subtle state does not alter the crystal lattice symmetry meaning that it is virtually invisible to normal photon and neutron scattering techniques, since these techniques are mainly sensitive to the core electron scattering and the nuclear scattering, respectively 8 , 9 . By contrast, conventional scanning tunneling microscopy (STM) has reported evidence for () orbital order on the surface of CeCoIn5 10 , revealing an opportunity for quasiparticle scattering interference (QPI) imaging , a powerful technique for detecting subtle orbital selective effects 3,11 .
( ) is the so-called T-matrix, representing the possible scattering processes between states | ⟩ and | + ⟩ for a point like s-wave scatterer. Atomic scale imaging of these interference patterns ( , ) is achieved using spatial mapping of the differential conductance, ( , ) 14 .

Modeling the () orbital order
As a concrete model, we consider orbital order of dxz/dyz-orbitals on a 2D square-lattice ( Fig. 1a). To accommodate the () orbital order, the unit cell is enlarged to a two-sublattice basis allowing for the incorporation of a staggered, nematic orbital order preserving the translational and global 4 -symmetry. Including superconductivity the model Hamiltonian takes the form, where the Nambu spinor is defined as with , , ( ) annihilating an electron with momentum and spin at sublattice in orbital .Here ℋ 0 ( ) contains intra-and interorbital nearest-and next-nearest-neighbour hoppings allowed by the d-wave symmetry of the orbitals, ℋ ( ) introduces the on-site anti-ferro-orbital order and Δ d ( ) contains nearest-and next-nearest-neighbour intraorbital d-wave pairings as introduced in Ref. 15. Here for generality we consider the simplest model Hamiltonian (ℋ 0 ( ), ℋ ( )) rather than specific Hamiltonian of CeCoIn5. As the model is spin-independent we suppress the spin index below. To separate the energy scales of the orbital and superconducting orders, the orbital order is introduced at an energy scale well above the superconducting gap, i.e. Δ oo ≫ Δ d . The Hamiltonian in (2) is chosen as a minimal model approach where ℋ 0 ( ) describes the simplest band dispersion allowing for the implementation of local 4 -symmetry breaking. A detailed description of the model and parameters can be found in SI Section 1.
To simulate QPI, a non-magnetic impurity is introduced as a point-like potential. We choose an on-site impurity as the scattering center, because this kind of impurity widely exists in the crystals and is located at a high symmetry point required to detect the local symmetry breaking caused by the orbital order. The impurity, either a different element or lattice vacancy, is assumed to exhibit a trivial spatial structure leading to identical potential strengths in the orbital degree of freedom. The local density of states (LDOS) is computed using a T-matrix approach as where is the real-space position of the two-ion unit cell, ∈ { = , ; = , }, the T-matrix is given by is the bare, retarded Greens function obtained from (2). We always insert the impurity at one of the two sites in the unit cell positioned at = for simplicity. Note that ( , , ) contains four components for the unit cell at , corresponding to the orbital and sublattice degrees of freedom. The position of a single lattice site, , is uniquely mapped from the set { , } enabling a straightforward transition to the site-resolved LDOS. To allow for reliable comparison to experimental data, we calculate the local density of states above the surface of the material following a simplified method that takes the Wannier orbitals into account 16,17 and basically weigh the computed ( , ) by atomic-like dxz/dyz orbitals. To account for the experimental resolution of 100 eV, we perform an additional Gaussian energy convolution, details on these calculations can be found in SI Section 2.

Consequences for QPI of () orbital order
The consequences of this () orbital order for QPI experiments are intriguing.
Surprisingly, theoretical modeling for the r-space QPI patterns, ( , E), around the impurity at sublattice a (Fig. 2a) and sublattice b (Fig. 2d) at energy | | > Δ well outside the superconducting gap, show almost identical ( , ) . At energies | | < Δ , however, the situation is radically different. Here ( , ) around chemically identical impurity atoms at sublattice a (Fig. 2b) and sublattice b (Fig 2e) are vividly different. The key consequence is that the amplitude of scattering interference is far more intense along one axis than along the other axis, depending on which sublattice the impurity atom resides. The interference pattern breaks C4-rotational symmetry, indicating the existence of the hidden () orbital order which breaks C4-symmetry locally. We stress that the impurity potential itself is pointlike and of identical strength on both orbitals. The C4-symmetry breaking takes place because the impurity chooses a specific sublattice, in conjunction with the underlying orbital order.
To quantify this local symmetry breaking effect, we define a dimensionless value Meanwhile, at the energy | | > Δ, still within the energy scale of the orbital order (Δ ), ( , ) is less than 2% (Fig.S2). Thus, the orbital order can be clearly unraveled below the energy scale of the superconducting order parameter because of the opening of the superconducting gap selectively enhances its visibility. For comparison, the QPI simulation of the normal-state model at the energies | | < Δ can be seen Fig.S3, which is equivalent to the | | > Δ case of the superconducting model.

QPI signature of () orbital order in CeCoIn5
To explore these predictions we studied CeCoIn5, a prototypical heavy-fermion superconductor, whose crystal unit cell has dimensions a=b=4.6Å, c=7.51Å and with A standard Co terminated surface topography ( ) is shown in Fig. 3a with sublattices marked by red dots and blue dots, respectively. The Co terminated surface is identified from both the tunneling conductance spectrum and the domain boundaries (SI section4). In this field of view (FOV), we find two single atom defects allocated at sublattice a and b, respectively. These two defects are nearly identical in topography image (Fig. 3a). gives the comparison of ( , ) surrounding the same defect in the superconducting state ( < ) (Fig. S5a,b,c) and in the normal state( > ) (Fig. S5d,e,f). The local anisotropy ( , ) is only enhanced at E=0 in the superconducting state while has no apparent change at E=0 in the normal state, in agreement with the theoretical prediction in Next, we study the local anisotropy ( , ) around the defects at the two sublattices.  (Fig. 4d). Obviously, the conductance anisotropy is rotated by 90-degrees for scattering centres at the different sublattice sites. To analyse the energy-dependence of ( , )we plot in Fig. 4b,e, the line profiles of ( , ) along the high symmetry directions (0,1) and (1,0) versus bias. The anisotropy is very weak (light blue and light red) at the energies outside the superconducting gap, while, inside the superconducting gap, the anisotropy rapidly increases (dark blue and dark red) and its maxima are indistinguishable from E~0. Moreover, the curves of ( , ) at the second atom site away from the defect center (region marked by black squares in Fig. 4a,d) also exhibits this property (Fig. 4c, f). For comparison, we plot the theoretical curve of ( , ) along the same high symmetry directions at each energy in Fig. S4. The theory curve features the same tendency as the experimental curve that ( , ) is significantly enhanced inside the superconducting gap and the maximum of ( , ) is indistinguishable from E~0.
Finally, we use a multi-atom (MA) averaging technique resolved by sublattice to establish the repeatability of these phenomena for all equivalent impurity atoms. Figure 5a,b and 5e,f indicate the scattering centers at sublattice a ( Fig. 5a,b) and sublattice b (Fig. 5e,f) marked by red circles that are involved in the MA analysis. The MA technique averages the mapping data over several defects located at the same sublattice 26 .
Since multiple sites are averaged, the random local distortion and noise are suppressed, and the common feature surrounding the defects is enhanced (SI Section 3). Figures 5c,g present the MA-averaged topography and differential conductance map ( , ) at E=0 for impurity atoms on sublattice a and b, respectively. The similar features seen in MA-averaged differential conductance map ( , ) (Fig. 5c,g) and single-defect differential conductance map ( , ) (Fig. 3c) reveals that the scattering interferences from the two sublattices are highly distinct and repeatably rotated by 90-degrees relative to each other. One advantage of the MA process is that, since the random features in the mapping image are suppressed, the averaged image can be regarded as a single defect that resides in a defect-free large FOV, even though the actual sample is defective. This advantage allows us to Fourier transform the interference signal surrounding the single defect with high resolution in q-space. We set the real-space origin ( = ) at the defect site and focus on the real part of Fourier  (Fig. S7). Note that our theoretical model is based on a simple band dispersion, as described in SI Section 1.

Discussion
In this work we have explored the QPI signatures of () orbital order in CeCoIn5. The subtlety of such orders is in their preservation of crystal lattice symmetries which makes them undetectable by traditional scattering techniques 8,9 . On the other hand, pioneering STM visualization studies of anisotropic electron density due to orbital order has been reported 27, 10 . Such experiments must be carried out under extreme tunneling conditions, for example at currents >10nA that, according to the Tersoff-Hamann theory 28 , require a miniscule tip-sample distance. Such tip-surface distances usually challenge the stability of the STM junction and, moreover, the tip-sample interaction may then become so intense as to alter the sample properties. By contrast, taking CeCoIn5 ( , ) orbital order as an example, we have explored the possibility of using conventional junction QPI to detect the local symmetry breaking orbital order. From theory, it was predicted that, even with an isotropic impurity, the underlying orbital order should reveal itself as a sublattice-selective anisotropy in the surrounding QPI pattern, due to the different effective coupling of the impurity to the two orbitals. This is because although the impurity is described as a simple point like potential with no spatial or orbital structure, the scattering T-matrix reflects the orbital order. This suggests strongly that the specific type of impurity Information for additional details on data treatment and extraction.

Theory:
The two-dimensional square lattice including staggered orbital order and superconductivity has been modelled by the Bogoliubov-de Gennes Hamiltonian in Equation Data Availability All data are available in the main text on Zenodo 31 . Additional information is available from the corresponding author upon reasonable request.      c,g, Simultaneous MA-averaged differential conductance map ( , ) at ~0 around the defects at sublattice a (c) and sublattice b (g).

Interplay of Hidden Orbital Order and Superconductivity in CeCoIn5
Weijiong where and σ ( i = 1,2,3 ) are the Pauli matrices in sublattice and orbital space, respectively, and Δ = 0.25| 1 | is the orbital order parameter. Note that because of the In this work, we only consider the simplest model Hamiltonian including staggered orbital order and it is not identical to the real Fermi surface of CeCoIn5. We do not discuss a more complete model including both Ce and In atoms and the superconductivity originating from Ce atoms, since such issues are both beyond the scope of our current work and not relevant to its conclusions. Nevertheless, as shown in Fig. S7, the overall pattern of the real part of BQPI is still present in a good agreement between the calculation and the experiment except some inconsistencies in the exact period of the Friedel oscillations. This implies that our model indeed captures the key ingredients of symmetry-breaking QPI induced by the orbital order.

Quasiparticle Interference Simulation of Two Sublattice Scatterings
The local density of states (LDOS) is computed using a T-matrix approximation as ( , , ) = − 1 Im( ( , ) + ( , ) ( , ) (− , )) (S10) where is the real-space position of the two-ion unit cell, ∈ { = , ; = , }, the T-matrix is given by To improve the direct comparison to experiment we implement two modifications to the calculated LDOS. First, we take into account that the tunnelling process to the STM tip is in the exponential limit, i.e. the STM tip is several Å above the surface and the tunnelling conductance is proportional to the local density of states at that position. . σ = Δ/12 corresponds to the experimental energy resolution ~100meV. Fig. S1 gives the direct comparison between calculated LDOS ( ) and measured density of state ( ) far from the impurities and at the impurity/defect site to reveal that our choice of model parameters allows to describe the spectral properties of the impurities in the experiment. In Fig. S2 we present the anisotropy in real space (see main text) as calculated in the superconducting state showing that the anisotropy at zero energy (E=0) is much larger than at a finite energy above the superconducting gap. In contrast, Fig. S3 shows the corresponding result as obtained without superconducting order (normal state) where the anisotropy is very small for both E = 0 and at finite energy. As a consequence, the orbital order would be difficult to detect in this local probe. The same conclusions can be read off from Fig. S4 where the anisotropy is plotted as function of energy with and without superconducting order.
We point out that the authors of Ref. [4] report a similar symmetry-breaking QPI caused by the nematicity in the FeSeTe system. Their observation is distinct from our result in two aspects. First, the global QPI analysis they perform discovers the order that breaks overall crystal lattice symmetry but should not yield anti-ferro-orbital order which perserves the global C4 symmetry as in CeCoIn5. Second, the nematicity they discover is only observed in the non-superconductive state at high energy beyond the coherence peak.
This is also distinct from our present picture where the anisotropy from orbital order is significantly enhanced within the superconductive quasiparticles below the superconducting gap.

Multi-atom Technique
The multi-atom technique 5 is overlapping and averaging the same type of defects to suppress the random noise and highlight the common features of the defects. we first identify the coordinates of the centers of the defects from the topography. The defects are chosen in the topography by eye selecting only those of the same type as in Fig. 3 and 4, since they are well allocated at sublattice site a/b. Then, the selected defects at sublattice a/b are distinguished by the surrounding scattering pattern in ( , = 0) map.
Here we choose one defect as example. The chosen defects are marked by a red circle in Fig. 5a,e (topography) and 5b,f ( ( , = 0)). The exact coordinate = ( , ) for the center of the defect is figured out by calculating the center of mass of with intensity as the weights in the subsidence region around the impurity as seen in the topography.
The shift operation is done in the q-space on Fourier transformed map ( ) and

Identification of Termination Surface
As reported by Ref. [6], three different cleaved surfaces can be found in CeCoIn5: Co surface, CeIn surface and In2 surface. We first rule out the In2 surface because this surface is reconstructed. Both Co surface and CeIn surface show the same lattice constant ~4.6Å. Here, we identify our measured surface as Co surface by two observations.
First, Fig. S8 shows the measured scanning tunneling spectrum on our sample surface. The spectrum presents a dip at ~5meV corresponding to the heavy-fermion hybridization gap. This spectrum is exactly the same as the spectrum measured on Co surface in Ref. [7], except that our energy resolution is better. However, the spectrum of the CeIn surface in Ref. [7] displays that the density of states at -30meV is larger than that at 30meV, different from what we observed.
Second, since the orbital order breaks the equivalence of Co sites in sublattice a and b, two degenerate states should appear on the surface. At the interface of these two degenerate states, domain boundaries should form. Ref. [7] reports that the domain boundaries only appear on Co surface, implying that the orbital order occurs only on Co surface. We also observe many domain boundaries on our measured surface (Fig. S9), indicating our cleaved surface is Co-terminated.
Furthermore, Fig. S10 shows two nearby domains with several defects close to the domain boundary. In the ( , = 0) map (Fig. S10c)