### CDW visualization in non-superconductive UTe_{2}

Table of Contents

#### Differential conductance imaging of CDW at *T* = 4.2 K

At *T* = 4.2 K and using superconducting tips to study the UTe_{2} (0–11) surface, we measure differential tunnelling conductance spectra *g*(** r**,

*V*) to visualize the CDW in the normal state of UTe

_{2}. Extended Data Fig. 1a–d shows

*g*(

**r**,

*V*) images

*V*= −7 mV, −15 mV, −23 mV and −29 mV with Fourier transform

*g*(

**q**,

*V*) shown as Extended Data Fig. 1e–h. Three CDW peaks at

**Q**

_{1,2,3}occur in all

*g*(

**q**,

*V*), representing incommensurate charge-density modulations with energy scale up to approximately 30 meV, consistent with ref.

^{24}.

#### CDW visualization at incommensurate wavevectors **Q**

_{1,2,3}

To calculate the amplitude \({g}_{{{\bf{Q}}}_{i}}({\bf{r}})\) of the CDW modulation represented by the peaks at **Q**_{i} (*i* = 1, 2, 3), we apply a two-dimensional computational lock-in technique. Here *g*(**r**) is multiplied by the term \({{\rm{e}}}^{{\rm{i}}{{\bf{Q}}}_{i}\cdot {\bf{r}}}\) and integrated over a Gaussian filter to obtain the lock-in signal

$${g}_{{{\bf{Q}}}_{i}}\left({\bf{r}}\right)=\frac{1}{\sqrt{2\pi }{\sigma }_{{\bf{r}}}}\int {\rm{d}}{\bf{R}}g\left({\bf{R}}\right){{\rm{e}}}^{{\rm{i}}{{\bf{Q}}}_{i}\cdot {\bf{R}}}{{\rm{e}}}^{\frac{{\left|{\bf{r}}-{\bf{R}}\right|}^{2}}{2{\sigma }_{{\bf{r}}}^{2}}}$$

(9)

in which *σ*_{r} is the cutoff length in the real space. In **q** space, this lock-in signal is

$${g}_{{{\bf{Q}}}_{i}}({\bf{r}})={{\mathscr{F}}}^{-1}{g}_{{{\bf{Q}}}_{i}}({\bf{q}})={{\mathscr{F}}}^{-1}\,\left[{\mathscr{F}}(\,g({\bf{r}}){{\rm{e}}}^{{\rm{i}}{{\bf{Q}}}_{i}\cdot {\bf{r}}})\cdot \frac{1}{\sqrt{2\pi }{\sigma }_{{\bf{Q}}}}{{\rm{e}}}^{-\frac{{q}^{2}}{2{{\sigma }_{{\bf{Q}}}}^{2}}}\right]$$

(10)

in which *σ*_{Q} is the cutoff length in **q** space. Here \({{\sigma }_{{\bf{r}}}=1/\sigma }_{{\bf{Q}}}\). Next, *g*_{Q}(**r**, *V*), the inverse Fourier transform of the combined **Q**_{i} (*i* = 1, 2, 3) CDWs, is presented in Extended Data Fig. 1i–l.

To specify the effect of filter size on the inverse Fourier transform, we show in Extended Data Fig. 1m–t the real-space density of states *g*(**r**, 10 mV), its Fourier transform *g*(**q**, 10 mV) and the evolution of inverse Fourier transform images as a function of the real-space cutoff length *σ*_{r}. The differential conductance map *g*_{Q}(**r**, 10 mV) is shown at a series of *σ*_{r}, including 10 Å, 12 Å, 14 Å, 18 Å, 24 Å and 35 Å. The distributions of the CDW domains in the filtered *g*_{Q}(**r**, 10 mV) images with cutoff lengths of 10 Å, 12 Å, 14 Å, 18 Å and 24 Å are highly similar. The cutoff length used in Fig. 2d is 14 Å, such that the domains of the CDW modulations are resolved and the irrelevant image distortions are excluded. The same filter size of 14 Å is chosen for all three **Q**_{i} vectors. Formally, the equivalent inverse Fourier transform analysis is carried out for Fig. 4c,d but with a filter size of 11.4 Å to filter both the CDW and PDW peaks.

#### Simulated UTe_{2} topography

To identify **q**-space peaks resulting from the (0–11) cleave-plane structure of UTe_{2}, we simulate the topography of the UTe_{2} cleave plane and Fourier transform. Subsequently, we can distinguish clearly the CDW signal from the structural periodicity of the surface. The simulation is calculated on the basis of the ideal lattice constant of the (0–11) plane of the UTe_{2}, **a*** = 4.16 Å and inter-Te-chain distance **b*** = 7.62 Å. Extended Data Fig. 2a is a simulated *T*(**r**) image in the FOV of 14.5 nm. The simulated topography *T*(**r**) is in good agreement with experimental *T*(**r**) images presented throughout. The Fourier transform, *T*(**q**), of the simulated *T*(**r**) in Extended Data Fig. 2b shows six sharp peaks, confirming that they are the primary peaks resulting from the cleave-plane structure. Most notably, the CDW peaks in Fig. 2c are not seen in the simulation. They are therefore not caused by the surface periodicity but instead originate from the electronic structure, as first demonstrated in ref. ^{24}.

### Normal-tip PDW detection at the NIS gap edge of UTe_{2}

Initial STM searches for a PDW on UTe_{2} were carried out using a normal tip at 280 mK. Extended Data Fig. 3a shows a typical line cut of the \({\rm{d}}I\,/\,{\rm{d}}V{| }_{{\rm{NIS}}}\) spectrum taken from the FOV in Extended Data Fig. 3b. There is a large residual density of states near the Fermi level. The gap depth *H* is defined as the difference between the gap bottom in the \({\rm{d}}I\,/\,{\rm{d}}V{| }_{{\rm{NIS}}}\) spectrum and the coherence peak height, that is, \(H\equiv {\rm{d}}I\,/\,{\rm{d}}V\,{{\rm{| }}}_{{\rm{NIS}}}(V\equiv {\Delta }_{{{\rm{UTe}}}_{{\rm{2}}}})-{\rm{d}}I\,/\,{\rm{d}}V\,{{\rm{| }}}_{{\rm{NIS}}}(V\equiv 0)\). Its modulation is extracted from the \({\rm{d}}I\,/\,{\rm{d}}V{| }_{{\rm{NIS}}}\) line cut and presented in Extended Data Fig. 3c; it modulates perpendicular to the Te atom chains.

Conventional, NIS tunnelling discloses superconducting energy-gap modulations as shown in Extended Data Fig. 3a. The superconducting energy gap is defined as half of the peak-to-peak distance in the \({\rm{d}}I\,/\,{\rm{d}}V\,{{\rm{| }}}_{{\rm{NIS}}}\) spectrum (Fig. 3a and Extended Data Fig. 3d). Its magnitude \(| {\Delta }_{{{\rm{UTe}}}_{{\rm{2}}}}| \) is found to lie approximately between 250 µeV and 300 µeV. We measure variations in the coherence peak position from the \({\rm{d}}I\,/\,{\rm{d}}V{| }_{{\rm{NIS}}}\) spectrum at each location **r**. The two energy maxima near \({\Delta }_{{{\rm{UTe}}}_{{\rm{2}}}}\) of each \({\rm{d}}I\,/\,{\rm{d}}V\,{{\rm{| }}}_{{\rm{NIS}}}\) spectrum are fitted with a second-order polynomial function (\({R}_{{\rm{RMS}}}^{2}=0.87\)). The energy gap is defined as the maxima of the fit, Δ_{+} for *V* > 0 and Δ_{−} for *V* < 0. The total gap map \({\Delta }_{{{\rm{UTe}}}_{{\rm{2}}}}({\bf{r}})\equiv [{\Delta }_{+}({\bf{r}})-{\Delta }_{-}({\bf{r}})]/2\) is derived from Δ_{+} and Δ_{−} (Extended Data Fig. 3e). The Fourier transform of \({\Delta }_{{{\rm{UTe}}}_{{\rm{2}}}}\left({\bf{r}}\right)\), \({\Delta }_{{{\rm{UTe}}}_{{\rm{2}}}}\left({\bf{q}}\right)\) (Extended Data Fig. 3f), shows three peaks at wavevectors **P**_{i=1,2,3}. They are the initial signatures of the energy-gap modulations of the three coexisting PDW states in UTe_{2}.

### Superconductive-tip PDW visualization at the SIS gap edge of UTe_{2}

#### Tip preparation

Atomic-resolution Nb superconducting tips are prepared by field emission. To determine the tip gap value during our experiments, we measure conductance spectrum on UTe_{2} at 1.5 K (*T*_{c} = 1.65 K), in which the UTe_{2} superconducting gap is closed. The tip gap \(| {\Delta }_{{\rm{tip}}}| \cong 1.37\,{\rm{meV}}\) is the energy of the coherence peak (Extended Data Fig. 4a). The measured spectrum is fitted using a Dynes model^{40}. The typical \({\rm{d}}I\,/\,{\rm{d}}V{| }_{{\rm{SIS}}}\) measured at 280 mK on UTe_{2} (Fig. 3b) shows the total gap value \(E=| {\Delta }_{{\rm{tip}}}| +| {\Delta }_{{{\rm{UTe}}}_{{\rm{2}}}}| \approx 1.62\,{\rm{meV}}\). Therefore, we estimate \(| {\Delta }_{{{\rm{UTe}}}_{{\rm{2}}}}| \approx 0.25\,{\rm{meV}}\), consistent with a previous report^{9} and the \({\rm{d}}I\,/\,{\rm{d}}V{| }_{{\rm{NIS}}}\) shown in Fig. 3a and Extended Data Fig. 3.

#### SIS tunnelling amplification of energy-gap measurements

To determine the energy of *E*_{+}(**r**) and *E*_{−}(**r**) at which the maximum conductance in \({\rm{d}}I\,/\,{\rm{d}}V{| }_{{\rm{SIS}}}(V)\) occurs, we fit the peak of the measured \({\rm{d}}I\,/\,{\rm{d}}V{| }_{{\rm{SIS}}}(V)\) spectra using a second-order polynomial fit:

$$g(V)=a{V}^{2}+bV+c$$

(11)

This polynomial closely fits the experimental data. Extended Data Fig. 4b,c shows two typical \({\rm{d}}I\,/\,{\rm{d}}V{| }_{{\rm{SIS}}}(V)\) spectra measured at +*V* and −*V* along the trajectory indicated in Fig. 3c. The evolution of fits *g*(*V*) in Extended Data Fig. 4d,e shows a very clear energy-gap modulation.

#### Shear correction and Lawler–Fujita algorithm

To register several images to precisely the same FOV, the Lawler–Fujita algorithm is applied to the experimental data. Then, to recover the arbitrary hexagon of the Te lattice, shear correction is applied to correct any image distortions caused by long-range scanning drift during days of continuous measurement.

To correct against picometre-scale distortions of the lattice, we apply the Lawler–Fujita algorithm. Let \(\widetilde{T}(\widetilde{{\bf{r}}})\) represent a topograph of a perfect UTe_{2} lattice without any distortion. Three pairs of Bragg peaks **Q**_{1}, **Q**_{2} and **Q**_{3} can be obtained from Fourier transform of \(\widetilde{T}(\widetilde{{\bf{r}}})\). Hence \(\widetilde{T}(\widetilde{{\bf{r}}})\) is expected to take the form

$$\widetilde{T}(\widetilde{{\bf{r}}})=\mathop{\sum }\limits_{i=1}^{3}{T}_{i}\cos \left({{\bf{Q}}}_{i}\cdot \widetilde{{\bf{r}}}+{\theta }_{i}\right)$$

(12)

The experimentally obtained topography *T*(**r**) may suffer from a slowly varying position-dependent spatial phase shift *θ*_{i}(**r**), which can be given by

$$T({\bf{r}})=\mathop{\sum }\limits_{i=1}^{3}{T}_{i}\cos \left({{\bf{Q}}}_{i}\cdot {\bf{r}}+{\theta }_{i}({\bf{r}})\right)$$

(13)

To get *θ*_{i}(**r**), we use a computational two-dimensional lock-in technique to the topography

$${A}_{{\bf{Q}}}\,({\bf{r}})=\int {\rm{d}}{\bf{R}}T({\bf{R}}){{\rm{e}}}^{{\rm{i}}{\bf{Q}}\cdot {\bf{R}}}{{\rm{e}}}^{-\frac{{\left({\bf{r}}-{\bf{R}}\right)}^{2}}{2{\sigma }^{2}}}$$

(14)

$${A}_{{{\bf{Q}}}_{i}}\left(\,{\bf{r}}\right)={{\mathscr{F}}}^{-1}{A}_{{{\bf{Q}}}_{i}}\left(\,{\bf{q}}\right)={{\mathscr{F}}}^{-1}\,\left[{\mathscr{F}}(T\left({\bf{r}}\right){{\rm{e}}}^{{\rm{i}}{{\bf{Q}}}_{i}\cdot {\bf{r}}})\cdot \frac{1}{\sqrt{2\pi }{\sigma }_{{\rm{Q}}}}{{\rm{e}}}^{-\frac{{q}^{2}}{2{{\sigma }_{{\rm{Q}}}}^{2}}}\right]$$

(15)

$$\left|{A}_{{\bf{Q}}}\left(\,{\bf{r}}\right)\right|=\sqrt{{\left({\rm{Re}}{A}_{{\bf{Q}}}\left({\bf{r}}\right)\right)}^{2}+{\left({\rm{Im}}{A}_{{\bf{Q}}}\left({\bf{r}}\right)\right)}^{2}}$$

(16)

$${\theta }_{i}({\bf{r}})={\tan }^{-1}\frac{{\rm{Im}}{A}_{{\bf{Q}}}\left(\,{\bf{r}}\right)}{{\rm{Re}}{A}_{{\bf{Q}}}\left(\,{\bf{r}}\right)}$$

(17)

for which *σ* is chosen to capture the lattice distortions. In the Lawler–Fujita analysis, we use *σ*_{q} = 3.8 nm^{−1}. Mathematically, the relationship between the distorted and the perfect lattice for each **Q**_{i} is \({{\bf{Q}}}_{i}\cdot {\bf{r}}{\boldsymbol{+}}{\theta }_{i}({\bf{r}})={{\bf{Q}}}_{i}\cdot \widetilde{{\bf{r}}}+{\theta }_{i}\). We define another global-position-dependent quantity, the displacement field \({\bf{u}}\left({\bf{r}}\right)={\bf{r}}-\widetilde{{\bf{r}}}\), which can be obtained by solving equations

$${\bf{u}}\left({\bf{r}}\right)={\left(\begin{array}{c}\begin{array}{c}{{\bf{Q}}}_{1}\\ {{\bf{Q}}}_{2}\end{array}\\ {{\bf{Q}}}_{3}\end{array}\right)}^{-1}\left(\begin{array}{c}\begin{array}{c}{\theta }_{1}-{\theta }_{1}\left({\bf{r}}\right)\\ {\theta }_{2}-{\theta }_{2}\left({\bf{r}}\right)\end{array}\\ {\theta }_{3}-{\theta }_{3}\left({\bf{r}}\right)\end{array}\right)$$

(18)

Finally, a drift-corrected topography, \(\widetilde{T}\left(\widetilde{{\bf{r}}}\right)\) is obtained by

$$\widetilde{T}\left(\widetilde{{\bf{r}}}\right)=T({\bf{r}}-{\bf{u}}\left({\bf{r}}\right))$$

(19)

By applying the same correction of **u**(**r**) to the simultaneously taken differential conductance map *g*(**r**), we can get

$$\widetilde{g}\left(\widetilde{{\bf{r}}}\right)=g({\bf{r}}-{\bf{u}}({\bf{r}}))$$

(20)

in which \(\widetilde{g}\left(\widetilde{{\bf{r}}}\right)\) is the drift-corrected differential conductance map.

#### Lattice registration of UTe_{2} energy gap \({{\boldsymbol{\Delta }}}_{{{\bf{UTe}}}_{{\bf{2}}}}{\boldsymbol{(}}{\bf{r}}{\boldsymbol{)}}\)

We measure two separate \({\rm{d}}I\,/\,{\rm{d}}V{| }_{{\rm{SIS}}}({\bf{r}},V)\) maps separated by several days and in two overlapping FOVs, with energy ranges −1.68 meV < *E* < −1.48 meV and 1.5 meV < *E* < 1.7 meV. Therefore, we obtain two datasets, *T*_{+}(**r**) with the simultaneous \({{{\rm{d}}I/{\rm{d}}V| }_{{\rm{SIS}}}}_{+}({\bf{r}},V)\) at positive bias and *T*_{−}(**r**) with the simultaneous \({{{\rm{d}}I/{\rm{d}}V| }_{{\rm{SIS}}}}_{-}({\bf{r}},V)\) at negative bias.

After the shear and Lawler–Fujita corrections are applied, the lattice in the corrected topographs of *T*_{+}(**r**) and *T*_{−}(**r**) become nearly perfectly periodic. Next, we perform rigid spatial translations to register the two topographs to the exact same FOV with a lateral precision better than 27 pm. Extended Data Fig. 5a,b shows two topographs of registered *T*_{+}(**r**) and *T*_{−}(**r**). Cross-correlation (XCORR) of two images *I*_{1} and *I*_{2}, *X*(**r**, *I*_{1}, *I*_{2}) at **r** is obtained by sliding two images **r** apart and calculating the convolution,

$$X({\bf{r}},{I}_{1},{I}_{2})=\frac{\int {I}_{1}^{\ast }({{\bf{r}}}_{1}){I}_{2}({\bf{r}}+{{\bf{r}}}_{1}){\rm{d}}{{\bf{r}}}_{1}}{\sqrt{\int {|{I}_{1}({{\bf{r}}}_{1})|}^{2}{\rm{d}}{{\bf{r}}}_{1}\int {|{I}_{2}({{\bf{r}}}_{2})|}^{2}{\rm{d}}{{\bf{r}}}_{2}}}$$

(21)

in which the denominator is a normalization factor such that, when *I*_{1} and *I*_{2} are exactly the same image, we can get *X*(**r** = **0**, *I*_{1}, *I*_{2}) = 1 with the maximum centred at (0, 0) cross-correlation vector. Extended Data Fig. 5c shows that the maximum of XCORR between *T*_{+}(**r**) and *T*_{−}(**r**) coincides with the (0, 0) cross-correlation vector. The offset of the two registered images are within one pixel. The multiple-image registration method is better than 0.5 pixels = 27 pm in the whole FOV and the maxima of the cross-correlation coefficient between the topographs is 0.93. All transformation parameters applied to *T*_{+}(**r**) and *T*_{−}(**r**) to yield the corrected topographs are subsequently applied to the corresponding \({\rm{d}}I\,/\,{\rm{d}}V{| }_{{\rm{SIS}}}({\bf{r}},V)\) maps obtained at positive and negative voltages.

#### Particle-hole symmetry of the superconducting energy gap \({{\boldsymbol{\Delta }}}_{{{\bf{UTe}}}_{{\bf{2}}}}{\boldsymbol{(}}{\bf{r}}{\boldsymbol{)}}\)

The cross-correlation map in Extended Data Fig. 5f provides a two-dimensional measure of agreement between the positive and negative \({\rm{d}}I\,/\,{\rm{d}}V{| }_{{\rm{SIS}}}\left(V\right)\) energy-maxima maps in Extended Data Fig. 5d,e. The inset of Extended Data Fig. 5f shows a line cut along the trajectory indicated in Extended Data Fig. 5f. It shows a maximum of 0.92 and coincides with the (0, 0) cross-correlation vector. Thus, it shows that gap values between positive bias and negative bias are highly correlated.

#### PDW visualization at incommensurate wavevectors **P**

_{1,2,3}

The inverse Fourier transform analysis for PDW state in Fig. 4c is implemented using the same technique described here in Methods. The filter size chosen to visualize the PDW is 11.4 Å. The inverse Fourier transform of the CDW in Fig. 4d is calculated using an identical filter size of 11.4 Å.

#### Independent PDW visualization experiments

To confirm that the PDW discovered is present in several FOVs, we show a typical example of the gap modulation Δ_{+}(**r**) from one different FOV in Extended Data Fig. 6. The \({\rm{d}}I\,/\,{\rm{d}}V{| }_{{\rm{SIS}}}({\bf{r}},V)\) map is measured in the voltage region surrounding the positive Nb-UTe_{2} energy maxima near 1.6 meV. The spectra in this FOV are fitted with a second-order polynomial and shear corrected as described here in Methods. The resulting gap map, *δ*Δ_{+}(**r**), is presented in Extended Data Fig. 6b. The Fourier transform of this map, *δ*Δ_{+}(**q**), is presented in Extended Data Fig. 6c. *δ*Δ_{+}(**q**) features the same PDW wavevectors (**P**_{1}, **P**_{2}, **P**_{3}) reported in the main text.

### Energy modulations of subgap Andreev resonances

Surface Andreev bound states must occur in *p*-wave topological superconductors^{41}. Moreover, based on the phase-changing quasiparticle reflections at the *p*-wave surface, finite-energy Andreev resonances should also occur in the junction between a *p*-wave and an *s*-wave superconductor^{42} and are observed in UTe_{2}. Inside the SIS gap, we measure the \({\rm{d}}I\,/\,{\rm{d}}V{| }_{{\rm{SIS}}}\left({\bf{r}},V\right)\) map in the energy range from −500 µeV to 500 µeV. The map is measured in the FOV in Extended Data Fig. 7a, the same FOV as in Figs. 3 and 4. Three conductance peaks are resolved at approximately −300 µeV, 0 and 300 µeV, annotated with green arrows in the typical subgap spectrum in Extended Data Fig. 7b. The energy maximum of the positive subgap states between 200 µeV to 440 µeV is assigned as *A*_{+}. The energy maximum of the negative subgap states between −440 µeV to −200 µeV is assigned as *A*_{−}. The averaged energy of the Andreev subgap states is defined as \({\Delta }_{{\rm{A}}}({\bf{r}})\equiv \left[{A}_{+}({\bf{r}})-{A}_{-}({\bf{r}})\right]/2\), which ranges from 300 µeV to 335 µeV (Extended Data Fig. 7c). Fourier transform of the subgap energies Δ_{A}(**q**) exhibit two sharp peaks at the PDW wavevectors **P**_{1} and **P**_{2} (Extended Data Fig. 7d).

In the case of two superconductors with very different gap magnitude, when the sample bias voltage shifts the smaller gap edge to the chemical potential, the Andreev process of electron (hole) transmission and hole (electron) reflectional plus electron-pair propagation may produce an energy maximum in d*I*/d*V*|_{SIS} at the voltage of smaller gap energy. Hence, the observations in Extended Data Fig. 7d may be expected if the UTe_{2} superconducting energy gap is modulating at the wavevectors **P**_{1} and **P**_{2}. Extended Data Fig. 7e shows that the energy of the Andreev states modulate in space with a peak-to-peak amplitude near 10 µeV (see histogram in Extended Data Fig. 7f).

### Enhancement of signal-to-noise ratio using superconductive tips

Superconducting STM tips provide an effective energy resolution beyond the Fermi–Dirac limit. They have therefore been widely used as a method of enhancing the energy resolution of STM spectra^{26,27,28,29,30,31}.

To better quantify the signal-to-noise ratio improvement of the measured energy-gap modulations, we compare the fitting quality of the superconducting gap maps obtained using a normal tip (Extended Data Fig. 3) and a superconducting tip (Fig. 4). The fitting quality is defined using the coefficient

$${R}^{2}\left({\bf{r}}\right)=1-\frac{{\sum }_{i=1}^{N}{\left[g\left({\bf{r}},{V}_{i}\right)-{\rm{d}}I/{\rm{d}}V\left({\bf{r}},{V}_{i}\right)\right]}^{2}}{{\sum }_{i=1}^{N}{\left[g\left({\bf{r}},{V}_{i}\right)-\bar{g}\left({\bf{r}}\right)\right]}^{2}}$$

(22)

in which \({\rm{d}}I\,/\,{\rm{d}}V\left(V\right)\) is the measured spectrum, *g*(**r**, *V*) is the fitted spectrum and \(\bar{g}\left({\bf{r}}\right)\) is the averaged fitted spectrum. Extended Data Fig. 8a shows a typical spectrum measured using a superconductive tip, \({\rm{d}}I\,/\,{\rm{d}}V{| }_{{\rm{SIS}}}\) from the FOV in Fig. 3c. Extended Data Fig. 8d is a typical \({\rm{d}}I\,/\,{\rm{d}}V{| }_{{\rm{NIS}}}\) spectrum measured using a normal tip from the FOV in Extended Data Fig. 3. The energy-maximum noise level is decisively lower in \({\rm{d}}I\,/\,{\rm{d}}V{| }_{{\rm{SIS}}}\) spectra than in \({\rm{d}}I\,/\,{\rm{d}}V{| }_{{\rm{NIS}}}\) spectra and the fitting quality \({R}_{{\rm{SIS}}}^{2}\) is substantially higher than \({R}_{{\rm{NIS}}}^{2}\).

Extended Data Fig. 8b,c shows maps of the fitting parameter *R*^{2} calculated from fitting the d*I*/d*V*|_{SIS} energy-maxima map obtained using a superconductive tip, that is, the \({\Delta }_{{{\rm{UTe}}}_{{\rm{2}}}}({\bf{r}})\) images presented in Fig. 3e,f. Extended Data Fig. 8e,f shows maps of *R*^{2} calculated from the coherence peak fitting of d*I*/d*V|*_{NIS} obtained using a normal tip, that is, the \({\Delta }_{{{\rm{UTe}}}_{{\rm{2}}}}({\bf{r}})\) images presented in Extended Data Fig. 3e. Comparing these *R*^{2} quality-of-fit parameter maps, we find that a much larger fraction of normal-tip coherence peak maps have poor correspondence with the fitting procedures used. For superconducting tips, the root-mean-square values of the fitting parameter, \({R}_{{\rm{RMS}}}^{2}\), are 0.98 and 0.99 for the positive and negative coherence peak fitting, respectively. The normal-tip \({R}_{{\rm{RMS}}}^{2}\) values are 0.87 and 0.86 for the positive and negative coherence peak fitting, respectively. The superconducting tip therefore demonstrably achieves a marked signal-to-noise ratio enhancement for evaluation of \({\Delta }_{{{\rm{UTe}}}_{{\rm{2}}}}\left({\bf{r}}\right)\) images.

As the signal-to-noise ratio is increased in the SIS-convoluted coherence peaks measured using a superconducting tip, it has been possible to resolve the UTe_{2} energy-gap modulations of order approximately 10 μV. Fundamentally, the energy resolution is associated with the ability of the superconductive tip to resolve the energy at which the d*I*/d*V*|_{SIS} coherence peak reaches its maximum amplitude. Consequently, we determine our energy resolution to be 10 μV.

Thus, the same superconductor energy-gap modulations in \({\Delta }_{{{\rm{UTe}}}_{{\rm{2}}}}\left({\bf{r}}\right)\) of UTe_{2} can be observed using either a superconducting tip or a normal tip. However, the former substantially increases the SIS conductance at \(\left|E\right|={\Delta }_{{{\rm{UTe}}}_{{\rm{2}}}}+{\Delta }_{{\rm{tip}}}\) and allows for considerably better imaging of these energy maxima and thus \({\Delta }_{{{\rm{UTe}}}_{{\rm{2}}}}\left({\bf{r}}\right)\).

### Interplay of subgap quasiparticles and PDW

Here we show simultaneous normal-tip-measured modulations of the UTe_{2} subgap states and \({\Delta }_{{{\rm{UTe}}}_{{\rm{2}}}}\left({\bf{r}}\right)\) at *T* = 280 mK, to study their interplay. Extended Data Fig. 9a shows the integrated differential conductance from −250 μV to 250 μV, \({\sum }_{-250\,{\rm{\mu }}{\rm{V}}}^{250\,{\rm{\mu }}{\rm{V}}}g\left({\bf{r}},E\right)\). Inverse Fourier transform of the three wavevectors **Q**_{1,2,3} from \({\sum }_{-250\,{\rm{\mu }}{\rm{V}}}^{250\,{\rm{\mu }}{\rm{V}}}g\left({\bf{r}},E\right)\) and **P**_{1,2,3} from the simultaneous \({\Delta }_{{{\rm{UTe}}}_{{\rm{2}}}}\left({\bf{r}}\right)\) in Extended Data Fig. 3e are compared in Extended Data Fig. 9c,d. Clearly, from the highly distinct spatial structure of these images, there is no one-to-one correspondence between the subgap density-of-states modulations and the simultaneously measured PDW energy-gap modulations in UTe_{2}. Overall, there is a very weak anticorrelation, with a cross-correlation value of −0.23 ± 0.05 that is not inconsistent with coincidence. Hence we demonstrate that there is no deterministic influence of the subgap density-of-states modulations on the PDW energy-gap modulations in superconducting UTe_{2}.

### Visualizing the interplay of PDW and CDW in UTe_{2}

The analysis of phase difference between PDW and CDW at three different wavevectors is shown in Extended Data Fig. 10. The inverse Fourier transforms of each CDW and PDW wavevector demonstrate a clear half-period shift between the two density waves (Extended Data Fig. 10a–f). This shift motivates the statistical analysis of the phase difference. The phase map of \({g}_{{{\rm{Q}}}_{1}}({\bf{r}},-9\,{\rm{mV}})\), \({\phi }_{1}^{{\rm{C}}}\left({\bf{r}}\right)\), and the phase map of \({\Delta }_{{{\rm{P}}}_{1}}({\bf{r}})\), \({\phi }_{1}^{{\rm{P}}}\left({\bf{r}}\right)\), are calculated. The phase difference between two corresponding maps is defined as \(| \delta {\phi }_{1}| ={\phi }_{1}^{{\rm{C}}}\left({\bf{r}}\right){\boldsymbol{-}}{\phi }_{1}^{{\rm{P}}}\left({\bf{r}}\right)\) for the **P**_{1}:**Q**_{1} wavevectors. Identical procedures are carried out for **P**_{2}:**Q**_{2} and **P**_{3}:**Q**_{3}. The histograms resulting from this procedure show that the statistical distributions of the phase shift \(|\delta {\phi }_{i}{\rm{|}}\) are centred around π (Extended Data Fig. 10j–l). Although the distribution varies, this π phase shift reinforces the observation of the spatial anticorrelation between CDW and PDW.

As shown in the inset of Extended Data Fig. 10g, the three PDW wavevectors are related by reciprocal lattice vectors: **P**_{2} = **P**_{1} − **G**_{3} and **P**_{3} = **G**_{1} − **P**_{1}. Nevertheless, the three UTe_{2} PDWs seem to be independent states when analysed in terms of the spatial modulations of the amplitude of the **P**_{1,2,3} peaks from Fig. 4 using equation (16). The amplitude of **P**_{1,2} has a domain width beyond 10 nm in the real space (Extended Data Fig. 10g,h). The amplitude of **P**_{3} is short-ranged, of which the averaged domain width is approximately 5 nm (Extended Data Fig. 10i). The one-pixel shift of **P**_{3} from the central axis is within the error bar of experimental measurements. The spatial distributions of the three PDWs are negligibly correlated with cross-correlation values of their amplitude of *X*(**P**_{1}, **P**_{2}) = −0.3, *X*(**P**_{1}, **P**_{3}) = 0.9 and *X*(**P**_{2}, **P**_{3}) = 0.28. The weak cross-correlation relationships indicate that the three PDWs are independent orders.