Quantum-well states at the surface of a heavy-fermion superconductor – Nature

STM experiments

Single crystals of URu₂Si₂ were grown by the Czochralski technique in a Tetra Arc Furnace. We scanned samples for a low residual resistivity and a high superconducting critical temperature, close to 1.5 K. Such samples were then cut in a bar shape with dimensions 4 × 1 × 1 mm³, with the long distance parallel to the c axis. We mounted the samples on the sample holder of a scanning tunnelling microscope. The scanning tunnelling microscope was mounted in a dilution refrigerator. The resolution in energy of the setup was tested by measuring the superconducting tunnelling conductance with the tip and sample of s-wave superconductors Al and Pb down to 100 mK (ref. ⁴¹). Details of image-rendering software are provided in refs. ^42,43. The scanning tunnelling microscope head features a low-temperature movable sample holder, which is used to cleave the sample at cryogenic temperatures^41,44. At the same time, and importantly for this study, the sample holder allows modifying many times the scanning window. The terraces discussed in this work were found in three different samples, after studying hundreds of fields of view.

Surface termination in URu₂Si₂

We focus on U-terminated surfaces. In Extended Data Fig. 1a, we show the URu₂Si₂ crystal unit cell highlighting the U, Si and Ru planes; their inter-layer distances are indicated in units of the c-axis lattice parameter. In Extended Data Fig. 1b–d, we show STM images corresponding to different surface terminations. These surfaces are all obtained after cryogenic cleaving. On the surfaces full of square-shaped terraces, we find the results obtained in the main text. An example is shown in Extended Data Fig. 1b. All the observed terraces are separated by c/2 ≈ 4.84 Å. In Extended Data Fig. 1c,d, we show terraces with a triangular shape, in which we do not observe the phenomena discussed in the main text. Here the distance between consecutive terraces is about 0.11c, about 0.39c or about 0.61c, which correspond, respectively, to the three possible distances between U–Si planes (coloured arrows in Extended Data Fig. 1a). Therefore, we see that the surfaces with terraces having a triangular shape correspond to Si layers, sometimes with a U layer in between. By contrast, the surfaces with terraces having a square shape are U terminated. Atomically resolved images inside terraces (Extended Data Fig. 1e) provide the square atomic U lattice with an in-plane constant lattice of a = 4.12 Å . In Extended Data Fig. 1f, we show a typical atomic-sized image on Si-terminated surfaces. We do not observe atomic resolution and have sometimes seen circular defects. Defects in the U-terminated surfaces are very different, as shown in Extended Data Fig. 1g–o. We distinguish two distinct types of defect. The defects can be either point-like protrusions (Extended Data Fig. 1g,h) or troughs (Extended Data Fig. 1i). Sometimes, defects are arranged in small-sized square or rectangular structures (Extended Data Fig. 1j–o). Most of these defects are probably because of vacancies or interstitial atoms in layers below the U surface layer.

Tunnelling conductance in the HO state

The tunnelling conductance of the HO state has been discussed in refs. ^{45,46,47,48,49}. We have reproduced the results, as shown in Extended Data Fig. 2a,b. The tunnelling conductance results from simultaneous tunnelling into heavy and light bands, as in other heavy-fermion compounds^50,51,52. The red line in Extended Data Fig. 2a for T = 18 K follows a Fano function

in which A is a constant of proportionality, q is the ratio between two tunnelling paths and E_Fano is the Fano resonance energy with width ({Gamma }_{{rm{F}}}=2sqrt{{(pi {k}_{{rm{B}}}T)}^{2}+2{({k}_{{rm{B}}}{T}_{{rm{K}}})}^{2}}), T_K being the Kondo temperature^45,46. For the fit, we include an asymmetric linear background owing to the degree of particle–hole asymmetry in the light conduction band^45,53. To account for the thermal broadening, we convolute the result with the derivative of the Fermi–Dirac distribution. We find q = 0.8 ± 0.5, E_Fano = 3 ± 1 mV, Γ_F = 22 ± 1 mV and T_K = 90 ± 5 K, consistent with previous reports^45,46.

Inside the HO phase (red line in Extended Data Fig. 2b), we use the same Fano function, multiplied by an asymmetric BCS-like gap function with an offset δ_E

The resulting function is convoluted with the derivative of the Fermi–Dirac distribution function. We find δ_E(4.1 K) = 1.5 ± 0.5 meV and Δ_HO(4.1 K) = 4.0 ± 0.5 meV, consistent with previous reports^45,46.

Note that we also observe further features at lower temperatures and smaller bias voltages (Extended Data Fig. 2c). The red line in Extended Data Fig. 2c is a fit described below. The features above the superconducting gap can also be roughly obtained by using two Lorentzians at ε₋ and ε₊ and an asymmetric background. Probably, the peaks at ε₋ and ε₊ are because of avoided crossings in the band structure of the 2DHF at very low energies. We notice that the small feature at ε₊ occurs at a very similar energy range as a kink in the band structure found at the surface of Th-doped URu₂Si₂ (refs. ^48,49). In the calculations we show below, we can identify features in surface f-derived bands that can be associated to such peaks in the tunnelling conductance. However, such features can form as a result of correlations elsewhere in the Brillouin zone as well.

Quantum-well states at terraces between steps

In Extended Data Fig. 3a, we show the tunnelling conductance background subtracted from Fig. 2a to obtain Fig. 2b. To carry out the background subtraction, we first identify the features at ε₋ and ε₊ in the conductance map. These are the light-blue regions centred at ε₊ and the red–yellow region centred at ε₋ in Extended Data Fig. 3a. We then identify the edge states occurring at the steps, given by the red areas at the sides of Extended Data Fig. 3a. Similar peaks are obtained on steps separating different terraces. The nature and shape of the 1DESs is discussed below. We then model these features by a set of Lorentzians and obtain the pattern shown in Extended Data Fig. 3a. We subtract this pattern from the experiment (Fig. 2a) to obtain the pattern shown in Extended Data Fig. 3b.

To model the quantum-well states, we use the Fabry–Pérot interferometer expression for the density of states g_FP(x, E) given by

with (k=sqrt{2{m}^{* }(E-{E}_{0})/{hbar }^{2}}), m* the electronic effective mass, r the reflection amplitude, ϕ the phase and L the width of the terrace²³. The Fabry–Pérot interferometer is an optical resonator made of semireflecting mirrors and provides a simple and insightful way to model electronic wavefunctions confined between two wells. More information on surface band structure and on quantum-well states by confinement is provided in refs. ^{3,9,54,55,56,57,58,59,60}. We assume a symmetric potential well with L = 20 nm, r = 0.5 and ϕ = −π. The pattern generated by equation (3) is shown in Extended Data Fig. 3c. White points provide the positions of quantized levels as in Fig. 2a,b. It is not difficult to see that the structure of quantized levels is renormalized together with the electronic band structure. Smaller electronic effective masses imply larger quantized level width and separation, and vice versa. However, the simple Fabry–Pérot model does not take into account relaxation by electron–electron interactions, which lead to the extra level broadening discussed in the main text and in refs. ^27,28,61.

The black lines in Fig. 2d are fits to the equation (3). To account for the behaviour at the edges, we add the equation (5) for the 1DES. We use the parameters extracted for the terrace L₃, discussed in Extended Data Table 1. We show further examples in Extended Data Fig. 4. Note that, in Fig. 2d, we use equation (3) along with the contribution from the 1DES and contributions for ε₋ and ε₊ at the bias voltages at which these features are observed in the tunnelling conductance.

The 2DHF quantization was observed on the surfaces of different URu₂Si₂ samples. In Extended Data Fig. 5a–c, we show the result on another sample. Notice here that terraces have different sizes. We show in Extended Data Fig. 5a the STM topography image. In Extended Data Fig. 5b, we show a height profile through the white line in Extended Data Fig. 5a. In Extended Data Fig. 5c, we represent the tunnelling conductance along the central terrace (L ≈ 57 nm) of this profile. We observe similar tunnelling conductance curves as those presented in the main text. Notice the features at ε₋ and ε₊. The quantized levels are also readily observed. These occur, however, at different energy values, as the size of the terrace L is different to that of the terrace in the main text. In Extended Data Fig. 5d, we represent the values of the quantized levels found in terraces of different sizes L by different colours; we show the dispersion relation of the 2DHF as a magenta line.

In Extended Data Fig. 5e, we show as coloured points the bias voltage dependence of the energy spacing ΔE between consecutive quantized levels for terraces L₃, L₄, the terrace with length L = 57 nm (shown in Extended Data Fig. 5c) and a terrace with length L = 27 nm (not shown). We can write that (Delta E={E}_{n+1}-{E}_{n}=left(frac{{hbar }^{2}{pi }^{2}}{2{m}^{* }{L}^{2}}right)left({(n+1)}^{2}-{n}^{2}right)), with n = 1, 2, 3,…. This gives a square-root dependence of ΔE on the energy, (Delta E={E}_{n+1}-{E}_{n}=left(frac{{hbar }^{2}{pi }^{2}}{2{m}^{* }{L}^{2}}right)left(2n+1right)propto sqrt{E}), shown in Extended Data Fig. 5e. In Extended Data Fig. 5f, we plot the average value of (frac{Delta E}{2n+1}) for each terrace as a function of L. We find the expected (frac{1}{{L}^{2}}) dependence.

Note that in Extended Data Fig. 5 and Fig. 1d, we do not perform any symmetrization. We can see a tendency of the quantized states to shift towards the sides of the terrace, giving intensity patterns that are slightly asymmetric. We have calculated the expected patterns for different reflection coefficients at each side of the terrace. This produces asymmetric patterns similar to those observed experimentally. However, it is difficult to separate such an asymmetry from the signal coming from the edge states in the tunnelling conductance. Lateral symmetrization thus remains the best way to analyse and understand quantum states in the terraces observed here in URu₂Si₂.

We use (hbar /tau ={Gamma }_{0}+left(| {E}_{0}| /4pi right){left[(E+| {E}_{0}| )/{E}_{0}right]}^{2}) to fit the energy dependence of the lifetime⁶¹. The parameter |E₀|/4π is a prefactor that fits our experiment well. The prefactor is sometimes provided as a number²⁸ and has also been estimated as (frac{{e}^{2}{k}_{{rm{F}}}^{2}pi }{4pi {{epsilon }}_{0},32,hbar {q}_{{rm{TF}}}}) or, equivalently, (frac{e{m}^{3/2}}{32times {3}^{5/6}{pi }^{2/3}{{epsilon }}_{0}^{1/2}{hbar }^{4}{n}^{5/6}}) (with ϵ₀ the dielectric constant, n the electron density and q_TF the Thomas–Fermi screening length)²⁷. These estimations provide similar values to |E₀|/4π.

To obtain the energy dependence of the reflection coefficient, r(E), we used the model described in ref. ²⁴. To this end, we consider a one-dimensional periodic array of scattering objects, each modelled by a square potential well of width b < L (L is the width of the terrace). A constant complex potential W provides confinement and coupling to the bulk states. We can then write

with (k=sqrt{frac{2{m}^{* }E}{{hbar }^{2}}}) and (q=sqrt{frac{2{m}^{* }(E-W)}{{hbar }^{2}}}). The reflection coefficient r(E) is given by r(E) = |R(E)|². For the dashed line in Fig. 2e, we use typical parameters for Cu, with m* = 0.46m₀ and W = (−2 − 1i) eV. We shift the obtained curve in energy to obtain a result within the energy range of our data. For the continuous line in Fig. 2e, we use m* = 17m₀ and W = (−18 − 5i) meV.

Band-structure calculations at U-terminated surfaces of URu₂Si₂

The band structure of bulk URu₂Si₂ has been analysed previously in detail using DFT calculations^62,63,64. Relevant results coincide with angle-resolved photoemission, STM and quantum oscillation studies^{22,65,66,67,68,69,70,71}.

Several surface states have been observed by angle-resolved photoemission spectroscopy^19,20,21,22. The surface state discussed in refs. ^19,20 is formed by a hole-like band with its maximum at −35 meV and is thus far from what we observe here. At the X point of the Brillouin zone, there are no bulk states. Angle-resolved photoemission spectroscopy measurements show hints of surface-like bands with two-dimensional character at these points²¹. We have taken a closer look at the X point through DFT calculations. To this end, we built a U-terminated supercell consisting of 37 atomic layers, giving a total of ten U layers (Extended Data Fig. 6a). We performed DFT calculations using the full-potential linearized augmented plane-wave method with local orbitals as implemented in the WIEN2k package⁷². Atomic spheres radii were set to 2.5, 2.5 and 1.9 Bohr radii for U, Ru and Si, respectively. We used a 19 × 19 × 1 mesh of k-points in the first Brillouin zone, reduced by symmetry to 55 distinct k-points. The RK_max parameter was set to 6.5, resulting in a basis size of approximately 5,400 (more than 100 basis functions per atom). Spin–orbital coupling was included in the second variational step⁷³ and relativistic local orbitals were included for U 6p_1/2 and Ru 4p_1/2 states. The basis for calculations of the spin–orbital eigenvalue problem consisted of scalar-relativistic valence states of energies up to about 5 Ry, resulting in a basis size of about 3,800. The local density approximation was used for the treatment of exchange and correlation effects^63,74.

In Extended Data Fig. 6b, we highlight in particular the U spin-up character of the obtained surface-projected band structure. The spin-down character is much less pronounced within the shown energy range. There are several bands inside gaps of the bulk band structure, but only those around the X point of the simple tetragonal Brillouin zone, X_st (see Extended Data Fig. 6c), are sufficiently shallow to provide large effective masses.

We find a surface state (upper inset of Extended Data Fig. 6b) that consists of two hybridized hole bands, forming an M-shaped feature close to the Fermi level. The dispersion relation found in our experiment (magenta line in the upper inset of Extended Data Fig. 6b) is compatible with the central part of the M-shaped feature.

1DES and HO within U layers

To analyse the 1DES at the step between two terraces, we use a one-dimensional Dirac-function-like potential at the step, V(x) = U₀δ(x − x_1DES), in which x_1DES is the position of the 1DES. We take U₀ = b₀V₀, with b₀ the width of the potential well and V₀ the energy depth (V₀ < 0). We add a complex potential, V(x) → (U₀ − iU₁)δ(x − x_1DES) to simulate the coupling of the 1DES to the bulk of the crystal. A schematic representation of this model is shown in Extended Data Fig. 7a. Solving the Schrödinger equation for E < 0, we obtain the Green’s function of the states in the potential well

in which ({lambda }_{x}=frac{{hbar }^{2}}{{m}^{* }}| {U}_{0}{| }^{-1}) is the decay length with m* the effective mass. E_1DES and η_1DES are the energy position and the energy broadening of the 1DES given by

in which δV is the height of the potential barrier of the well relative to the Fermi level.

We can now fit the tunnelling conductance at the 1DES using

convoluted with the derivative of the Fermi–Dirac distribution function. Extended Data Table 1 shows the extracted fitting parameters E_1DES, η_1DES, λ_x and x_1DES for the four different terraces L₁ to L₄ from Fig. 1.

From Extended Data Table 1, we see that the energy position and the energy broadening of the 1DES are independent of the terrace size, with average values of E_1DES = −0.52 ± 0.14 meV and η_1DES = 0.45 ± 0.06 meV. We also see that all the spatial features are always at the same position with respect to the step, x_1DES ≈ 4.0a₀, a₀ being the in-plane lattice constant, with a decay length λ_x ≈ 0.9 nm ≈ 2a₀. The latter indicates that 1DESs and 2DHFs couple when the decay length reaches a few interatomic distances. With the extracted average values from Extended Data Table 1 for λ_x, η_1DES and E_1DES, we obtain U₀ = 5.4 meVÅ, U₁ = 0.38 meVÅ and δV = 3.1 meV.

We can analyse the 1DES through the tunnelling conductance at a step (Extended Data Fig. 7b,c). At low bias voltages, we find a dip in the tunnelling conductance of a few nanometres at the upper side of the step (blue lines in Extended Data Fig. 7b; for example, at −1.2 mV). The dip fills with the 1DES at about E_1DES (red lines in Extended Data Fig. 7b at −0.4 mV) and empties again at higher bias voltages. This shows that charge depletion close to the step opens a gap in the band structure. The gap is filled at the resonant energy of the 1DES, as observed previously in metals^30,75,76. By normalizing the tunnelling conductance to its shape far from the step (Extended Data Fig. 7c), we can follow the decay of the 1DES into the quantum-well states of the 2DHF with the model described above (Extended Data Fig. 7a). The decay length is on the order of the inverse of the wavevector of the 2DHF.

Taking a closer look at the steps, we surprisingly find a notable in-plane anisotropy of the 1DES. As we see in Extended Data Fig. 7b, the 1DES is observed when crossing steps along the dashed red line in Extended Data Fig. 7b but not along the dashed blue line, as discussed in the main text. It is useful here to take a closer look at the HO as well. As discussed in refs. ^{14,34,64,77,78,79}, there is no dipolar (magnetic) or structural order related to the HO phase. Instead, the U lattice can present some sort of long-range electronic ordering, whose actual symmetry and shape is considered as a relevant and open mystery¹⁴.

Previous NMR measurements indicated the absence of fourfold in-plane symmetry breaking in bulk URu₂Si₂ (ref. ³⁵), whereas our STM data clearly show an in-plane symmetry breaking in the 1DES. At the surface, there can be marked changes in the electronic structure owing to a modification of the valence of uranium atoms^21,80. Rather, the breaking of the in-plane symmetry observed here in the 1DES suggests that a fundamental breaking of the near-surface electronic properties of the U lattice is at play in the HO phase.

Interplay between superconductivity and the 2DHF

We consider several parallel conduction channels between the tip and the surface. For simplicity, we take into account tunnelling into the 2DHF and into the feature of largest size at ε₋ (Extended Data Fig. 8a). The first channel, t₁, connects the tip with the 2DHF. The 2DHF is superconducting by proximity from the bulk superconductor, which we model using a coupling t_s. With the second channel, t₂, we connect the tip to other surface features. We can write the retarded Green function ({hat{G}}^{{rm{r}}}) as

with

in which E_2DHF is the energy associated to the 2DHF and includes the shift of energy owing to HO and Fano resonance, W is an energy scale related to the normal density of states at the Fermi level by ρ(E_F) = 1/(πW), Δ is the superconducting gap and η is a small energy relaxation rate. We have added the self energies iΓ_j/2, (j = 2DHF, −, +), with Γ_j the width of the resonance j.

The differential conductance is calculated as

with

Here f(E, T) is the Fermi–Dirac distribution at the energy E and temperature T and ({sigma }_{0}=frac{2{e}^{2}}{h}) is the quantum of conductance (with h being Planck’s constant and e the elementary charge). Notice that T is the transmission, not the density of states often used to discuss STM measurements in superconductors. Notice also that we take into account tunnelling into the 2DHF (({G}_{11}^{{rm{r}}})) and into ε₋ (({G}_{33}^{{rm{r}}})), with mixed contributions (({G}_{31}^{{rm{r}}}) and ({G}_{13}^{{rm{r}}})).

To fit the tunnelling conductance curves shown in Fig. 4a,b and Extended Data Fig. 8b, we subtracted an asymmetric background (see Extended Data Figs. 2 and 8b). In Extended Data Table 2, we show the parameters used to obtain the tunnelling conductance with temperature (shown as black lines in Fig. 4a,b) from equation (11). We do not vary the parameters η = 0.018 meV, ε_2DHF = 12 meV, Γ_2DHF = 1 meV, Γ₋ = 0.55 meV and Γ₊ = 0.14 meV with temperature. We see that the superconducting lifetime itself is practically negligible, η = 0.018 meV ≪  Δ = 0.2 meV. The large zero-bias conductance is not because of the incomplete coupling to the bulk, as t_s is close to one. There is further smearing coming from features at ε₊ and ε₋. Γ₊ and Γ₋ provide the smeared superconducting density of states and a finite tunnelling conductance at zero bias. Notice that the superconducting gap vanishes at the critical temperature T_c but that the strongly bias-voltage-dependent tunnelling conductance remains up to higher temperatures (Extended Data Fig. 8).

There are numerous evidences for d-wave or more complex superconducting properties in URu₂Si₂. The differences in the superconducting density of states between these superconducting states and s-wave superconductivity are relatively small, particularly because the tunnelling conductance obtained with the model described here provides smeared conductance curves. The same occurs for the temperature dependence of the superconducting gap. Instead, the shape and asymmetry of the tunnelling conductance at defects might provide the connection between the properties of the superconducting 2DHF and the unconventional superconductivity of the bulk.

Results at point defects

We analyse in more detail here the tunnelling conductance at defects. We plot the tunnelling conductance obtained on two different defects in Extended Data Fig. 9a as blue and red lines. Notice the pronounced electron–hole asymmetry, which provides curves that widely differ from the curves far from defects (black line in Extended Data Fig. 9a). As discussed above, we can identify two kinds of defects: protrusions with height increases of around 15 pm, probably because of interstitial atoms located beneath the surface (Extended Data Fig. 9b,d,f,h), and troughs of around 15 pm in size, probably because of vacancies beneath the surface (Extended Data Fig. 9c,e,g,i). The defects visibly affect the tunnelling conductance. We plot the tunnelling conductance at ε₊, 0 mV and ε₋ for both types of defect in Extended Data Fig. 9f,g. In Extended Data Fig. 9h,i, we show the spatial dependence of the tunnelling conductance along a crystalline axis for both types of defect. At the site of the defect, there is a pronounced electron–hole asymmetry, which is opposite for each kind of defect. Protrusions provide a substantially enhanced tunnelling conductance for empty states above the Fermi level, whereas troughs provide the opposite (Extended Data Fig. 9a). At zero bias, the troughs show a pronounced in-plane anisotropy and a reduction of the superconducting gap, whereas the protrusions are in-plane isotropic and show an opened superconducting gap. When leaving the defect, the usual behaviour is recovered after several nanometres (Extended Data Fig. 9h,i).

To explain the pronounced modification of the electron–hole asymmetry of the tunnelling conductance, let us consider electron correlations creating a large Fermi surface scenario at some portion of the band structure. At low temperatures, correlations provide an avoided band crossing owing to hybridized heavy and light bands. We can assume that it happens somewhere in the band structure of URu₂Si₂, for example, close to the Γ point, as suggested in refs. ^14,18,47. Then, we obtain at low temperatures a large Fermi surface with heavy electrons and a close-lying light band with smaller wavevectors above the Fermi level. We can then assume that both kinds of defect couple to different parts of such a correlated band structure. Following our experiments, protrusions create 2DHF coupling to the small light band and a large density of states for empty states above the Fermi level (red curve in Extended Data Fig. 9a). Troughs can instead favour coupling to the heavy band, giving the opposite behaviour (blue curve in Extended Data Fig. 9a). The superconducting gap is disturbed at troughs. Generally, we expect all sorts of defects to be pair breaking, as URu₂Si₂ is not a s-wave superconductor. However, two-dimensional quantized states screen pair-breaking interactions from the underlying superconducting bulk⁸¹. This suggests that the absence of in-gap states in protrusions is because of screening by the 2DHF. Troughs, however, produce a strong coupling to the heavy bands that carry the heavy superconducting state and thus also contribute to pair breaking.

In-gap states at troughs have a certain in-plane anisotropy (Extended Data Fig. 9g, middle panel). The anisotropy of the superconducting gap has been analysed with macroscopic measurements^82,83,84,85. There are indications for gap nodes, for instance, from specific heat measurements^82,83,85. Local vortex dynamics shows an in-plane fourfold or twofold pinning potential⁸⁶. When measured as a function of the angle, there is a pronounced out-of-plane anisotropy, which has been associated to nodes along the c axis, but there is little in-plane variation⁸⁵. Following such a nodal structure, several experiments propose a chiral k_z(k_x + ik_y) superconducting wavefunction^84,85,87. The surface states of a (k_x ± ik_y) superconductor are predicted to be chiral arc states connecting the Weyl nodes⁸⁸. It is so far unclear how these surface states of the superconducting phase adapt to the surface states of the normal phase, which appear as a consequence of the surface-induced perturbation of the crystalline potential. As mentioned, we observe an asymmetry at zero bias (Extended Data Fig. 9g, middle panel). To increase the signal-to-noise ratio, we were forced to apply a fourfold symmetrization. Thus, the anisotropy might be twofold or fourfold. The chiral k_z(k_x + ik_y) state is in-plane isotropic. However, the shape of in-gap states is determined by the anisotropy of the normal-state wavefunctions, together with the symmetry of the superconducting wavefunction, so that the observed in-plane anisotropy of the in-gap states probably reflects the normal-state in-plane anisotropy.