Superconductivity in a disordered metal with Coulomb interactions

We study electronic densities of states (DOS) of strongly disordered superconducting thin ﬁlms of TiN. We ﬁnd through scanning tunneling microscopy (STM) measurements that the DOS decreases towards the Fermi level in the normal phase obtained by applying magnetic ﬁelds. The DOS shows spatial ﬂuctuations whose length scale is related to the energy-dependent decrease in the DOS close to the Fermi level and is similar in normal and superconducting phases. This shows that Coulomb interaction in disordered metals strongly inﬂuences the normal DOS and suggests that reduced Coulomb screening leads to spatial variations in the superconducting DOS. DOI: 10.1103/PhysRevResearch


I. INTRODUCTION
Superconductor-insulator transition (SIT) is an exemplary quantum phase transition occurring in disordered superconducting films [1][2][3].The SIT manifests a dramatic change in the ground state of the system from superconducting (zero resistance) to insulator (infinite resistance) often occurring in a small parameter range.The resulting phase diagram bears similarities to cuprates and interface superconductors, which also show zero temperature phase transitions between insulating and superconducting phases [4][5][6][7][8][9].In those systems, the tuning parameter is chemical or electrical doping [7,10,11].Disorder remains very difficult to identify and study as it is probably intertwined with modifications of an anisotropic electronic band structure.In the SIT, the transition occurs purely as a function of the degree of disorder, and in the metallic side disorder has washed out all anisotropies [1,12].The disorder is characterized through the temperature dependence of the normal state resistivity, which increases with decreasing temperature in the insulating phase.
Studies of the SIT show some similar aspects as for instance studies in the cuprates, such as a pseudogap above the superconducting critical temperature T c [13][14][15][16][17]. Disorder modifies the ratio between Coulomb and superconducting coupling energies [18].Suppressed Coulomb screening results in the depletion of the normal state density of states (DOS) near the Fermi level, which is often termed the zero bias anomaly (ZBA) [19,20].In disordered Be films, the DOS evolves upon cooling from a Coulomb-like gap [21] to a well opened, semiconducting like gap which is often termed a hard gap [22].This points out that there can be a self-induced electronic granular structure [23,24] comprising superconducting droplets hosting a superconducting gap [18,25,26].In presence of strong disorder, the superconducting gap does not show coherence peaks, suggesting that Cooper pairs are localized [27].These experiments mostly address the properties of systems with very strong disorder.However, it has been suggested that the role of disorder might crucially modify the temperature dependence of the DOS in a metallic system showing an incipient decrease of Coulomb screening.Bartosch and Kopietz showed that, as a consequence of reduced Coulomb screening, a disordered metal loses the Fermi liquid regime at sufficiently low temperatures.Instead of the energy and temperature independent DOS of a usual metal, the disordered metal presents at very low temperature a DOS vanishing at the Fermi level with N (E ) ∝ |E | that is washed out into an energy independent DOS when increasing temperature [28,29].This suggests that the onset of the reduction in Coulomb screening occurs relatively far from the transition and leads to identifiable zero bias anomalies in disordered metals.Here we make a comparative study of the low energy normal and superconducting DOS in disordered metallic TiN.We find that the DOS in the normal phase has a strong V-shape anomaly which follows closely analytical expressions provided in Refs.[28,29].We also show that the spatially fluctuating DOS of the normal state leads to spatial fluctuations in the superconducting state, suggesting that the energy and spatial DOS dependencies are interlinked in presence of disorder.transport properties were chosen to be as close as possible to film TiN2 and TiN3 in Ref. [25].A detailed analysis of the morphology of the film using x-ray, transmittion electron microscopy (TEM), and atomic scale scanning tunneling microscopy (STM) is provided in Appendix A. The film consists of a network of strongly interconnected single crystal like structures oriented randomly on the substrate and separated by only small angle grain boundaries.The crystal lattice parameters are the same as bulk TiN.The STM results were carried out using a home made dilution refrigerator with a voltage resolution of about 15 μV [30][31][32].The junction resistance was always well above a M , much larger than the film's sheet resistance.We comment on details about STM measurements in high resistance samples in Appendix D and note that we normalize the results at a fixed bias voltage and that its value does not influence the obtained bias voltage dependence.We also note that generally we do not maintain the same field of view when modifying the magnetic field.We term our sample D15.We also provide in Appendix B results on the same sample, but with an additional heat treatment (2 minutes heating to 200 • C on air).The heat treatment lead to more insulating behavior.In D15, we performed atomic scale STM and spectroscopic studies at a magnetic field of 4 T and at zero field.We studied a grid similar to the one in Ref. [25].After the STM measurements we performed transport experiments obtaining the results expected from Ref. [33].For the heated sample, discussed in Appendix B, we were not able to scan or obtain atomic resolution, probably due to surface degradation from heat treatment.However, we could obtain tunneling conductance curves at a few locations and the results are consistent with the transport experiments.In the following, we focus mainly on results in D15.

III. RESULTS
Figure 1 shows the temperature evolution of tunneling conductance G(V, T ) together with the temperature dependence of resistance R(T ) [solid line in Fig. 1(b)] under magnetic fields in the normal phase.Here, T c = 1.1 K and H c2 = 2.6 T (see Appendix B, Fig. 5).The tunneling conductance G(V ) [Fig.1(a)] shows a very strong suppression near zero bias which evolves when decreasing temperature.The suppression goes in parallel to the increase of the resistance R when decreasing temperature.The noticeable suppression of G(V ) near E F persist up to 10 K, as well as the negative slope in R(T ).
In Fig. 2 (a), we show the tunneling conductance G(V ) for different temperatures.Note that there is a strong decrease of the tunneling conductance towards zero bias.To analyze this behavior in more detail we remind that the sample can be considered as a 2D metal in this temperature range [20,33].There is indeed a crossover from 3D to 2D behavior when the diffusion length, 2π hD/E , approaches the film thickness [20].Taking the film thickness d and diffusion coefficient D (see Table I in Appendix B), we obtain E 5 meV, which is the voltage range we work with.We use the model of Bartosch and Kopietz (BK) [28,29] which, as we will show, provides an accurate description of our data.The model treats a lowdimensional metal with long-range Coulomb interactions in the limit of small eV .The BK model was used for various 1D systems, in which a dip around zero bias was observed in the tunneling conductance [34][35][36].The model starts describing a usual metal and takes into account the absence of Coulomb screening when correlations are enhanced by disorder.The model provides the DOS as a function of temperature, FIG. 2. In (a), we show the tunneling conductance G(V ) as grey filled symbols for sample D15 (a) normalized to about 3 mV for different temperatures, marked in the figure and a magnetic field of 7 T (magnetic field dependence is shown in Appendix F, Fig. 12).Dashed lines provide the fits of the tunneling conductance discussed in the text.Lines show the corresponding DOS, N (eV ), obtained after eliminating the temperature induced smearing in G(V ).In (b), we show the derivative of the DOS N N = dN (E )/dE as a function of temperature.showing a dip which develops when decreasing temperature.
The energy scale for the dip is related to the scattering parameters, particularly τ 0 , the defect scattering time and an additional time related to electron-electron interactions, τ 1 .We can write, following Refs.[28,29], that the tunneling conductance is given by where we have taken into account temperature induced smearing of the Fermi function in the DOS.r 0 is a dimensionless measure for system's resistance [29] and τ 1 = 4τ 0 /(κl ) 4 , where κ is the Thomas-Fermi screening wave vector.For a good metal, the Thomas-Fermi screening length is short compared with the mean free path (κl 1) so that τ 1 τ 0 and we find flat G(V ) and an energy independent DOS [29].In our case, τ 1 and τ 0 are of the same order, leading to the observed ZBA [29].
In Fig. 2(a), we show as dashed lines the temperature evolution of G(V ) together with curves calculated with Eq. (1).Notice the similarity to the actual DOS, shown by the lines in Fig. 2(a), indicating that temperature smearing is weak as compared to the effect of Coulomb correlations.We use just two parameters, r 0 and τ 0 (taking τ 1 = τ 0 ), which already provides an excellent fit.We find τ 0 = 7.3 ± 0.3 × 10 −15 s.With this, we find R = 3 ± 0.05 k [from r 0 = R/(h/e 2 )], which is very close to sample's sheet resistance at room temperature R = 2.94 k .Furthermore, this value also describes the suppression of the superconducting critical temperature T c , following an expression originally proposed by Finkel'stein [37,38] (see also Appendix B).The suppression of T c is due to scattering and electron density fluctuations produced by the reduction of Coulomb screening.As we will see below, the length scale associated to spatial DOS fluctuations we measure here is related to τ 0 .BK model exemplarily shows how the Fermi liquid regime is lost by Coulomb interactions in low dimensionality.In essence, the Coulomb gap N (E ) ∝ |E | appears in the DOS in a system with finite conductance at a finite temperature at energies sufficiently close to the Fermi level [28].
Notice that we are in the peculiar situation where the normal state DOS N (E ) varies in the same energy scale as the superconducting gap.The shift in the chemical potential occuring in BCS theory below T c (μ 0 − 2 μ 0 where μ 0 is the chemical potential of the normal phase) is then modified by a factor N N (E F ), giving μ 0 − c N N (E F ) 2 (with c being a constant) [39].E F is taken in this simple approach to be the Fermi level.In Fig. 2(b), we show N N (E ) for different temperatures.We see that N N (E ) is a nonmonotonic function with a well developed maximum that rapidly increases with cooling.The energy where the maximum occurs is decreasing with temperature.For temperatures below T c , the maximum is very close to the Fermi level.We can take the point where N N (E ≈ E F ) is not zero as determining the shift in the chemical potential.This provides in our case a value that is considerably larger than the one found for noninteracting electrons.Furthermore, it increases most when crossing T c towards lower temperatures.We should note that the value of the gap in the DOS below T c (or the position of the quasiparticle peaks) is also considerably larger than the BCS superconducting gap value ≈1.76 k B T c [25].This suggests that there could be a link between the increased shift in the chemical potential due to the shape of the normal DOS and the gap in the DOS in the superconducting phase.On the other hand, at temperatures well above T c , the maximum in N N (E ) is shallow and separated from the Fermi level.Thus, the reduction of Coulomb screening and the decrease in the Fermi level DOS prevents the formation of Cooper pairs, explaining the decrease in T c when increasing disorder.The same tendency is obtained for sample D15e, as shown in Appendix B.
The normal phase DOS, according to BK, might also show small spatial fluctuations, if the scattering is spatially dependent, which implies a spatially dependent r 0 [28,29].We indeed find fluctuations in the normal phase tunneling conductance which connect to the zero field DOS fluctuations.
In Figs.3(a)-3(d), we show the spatially fluctuating tunneling conductance acquired in the normal and superconducting phases.The fluctuations in the tunneling conductance amount to small variations of a few percent as a function of the position.To obtain a length scale associated to these fluctuations we calculate the autocorrelation function of the tunneling conductance images.The autocorrelation function ACF( r) provides the statistical correlation between any two points in an image separated by a vector r (see Refs. [11,40] and Appendix C).As there are no symmetric patterns in our images, we plot the radially averaged ACF(r) in Figs.3(e) and 3(f).We observe that ACF(r) decreases exponentially with r.When we average over all bias voltages, we find Thus we identify a length scale d ≈ 20 nm, which is shared by normal and superconducting phases.Taking a value of v F = 2 × 10 6 m/s consistent with bandstructure calculations for bulk TiN [41] and the above value of τ 0 = 7.3 ± 0.3 × 10 −15 s, we find that d ≈ v F τ 0 .Thus the depression of the DOS induced by Coulomb interactions also leads to spatial fluctuations in the DOS.
The detailed spatial distribution of the conductance suffers however changes when entering the superconducting phase.The absolute value of ACF(r) is slightly larger at 4 T than at 0 T, indicating that the variations of the conductance are stronger in the normal phase.To characterize the spatial distribution of these variations, we have calculated the distribution function of multifractal exponents F (α).
In a usual metallic system, the amplitude of the wave function scales as the inverse of the size of the system L d where L is the lateral size and d the dimension, e.g., in 2D as ∝ L −2 .In a fully localized electronic system, wavefunctions decay however exponentially.It was shown that at the Anderson transition wavefunctions have a divergent localization length following power laws that are scale independent [42,43].Such a situation can be characterized by calculating the areas where the wavefunctions decay as ∝ L −α and obtaining from it F (α) (see also Appendix E).When the DOS shows just random changes, F (α) is only defined around α = 2.The function F (α) is very useful to identify situations when the DOS (which is proportional to the square of the electronic wavefunction) forms areas that have either a small or a large DOS, instead of totally random distributions.When the spatial variations show some structure and deviate from random behavior, F (α) broadens and becomes an inverted parabola whose maximum value is no longer located at α = 2. Often, the deviations from 2 are in the 0.1% range and have been associated to critical behavior at the localization transition [42][43][44][45], multifractal superconductivity in 2D layers of transition metal dichalcogenides [46][47][48][49] and to fractal superconductivity close to the SIT [50,51].
We observe a deviation from 2 in a similar, but smaller range in the normal phase under magnetic fields [Fig.3(h)].We are far from the SIT in the metallic site and the correlations we observe are of very small magnitude.Inside the superconducting phase we can however identify spatial structure.This occurs particularly for bias voltages of about half the quasiparticle peak position or at the same value of the BCS superconducting gap [Fig.3(h)] and suggests that there is an emergent granularity inside the superconducting phase.We can speculate that the tendency continues when approaching and crossing the SIT, probably with a shift in the maximum granularity towards the gap edge [50,51].

IV. DISCUSSION
Recently, authors of Ref. [52] analyzed the DOS for energies above the superconducting gap and found a connection between the Coulomb-like variation of the normal DOS and the size of the superconducting gap in thin NbN films.Since the upper critical field was too high, the low energy DOS could not be established.Nevertheless, they found a neat anticorrelation between the power law of the ZBA and the size of the superconducting gap, which is not related to the topography of the sample.As we state above, we do not correlate images at different magnetic fields, because we take the images at different fields of view.Here we focus instead on the size of the relevant length scales.In Ref. [52], they showed that the Coulomb interactions modify the value of the superconducting gap and are responsible for the granularity in the superconducting phase.Here we find instead that there are random variations of the ZBA which show a tendency to form spatial structures in the superconducting phase [Fig.3(h)].Other important measurements in NbN showed the presence of a pseudogap at high magnetic fields, with a nearly magnetic field independent temperature onset when cooling [17].The magnetic field induced a gradual suppression of superconductivity which became more granular with field, suggesting the presence at high magnetic fields of fluctuating Cooper pair islands connected among each other and leading to a pseudogap feature in the DOS [16,17].We observe a similar DOS at high magnetic fields, although the structure is established at zero field.Furthermore, we identify a connection between energy-and spatial-dependent DOS, which is well accounted for by taking exclusively into account the reduction of Coulomb screening.Finally, the magnetic field dependence, as shown in Appendix F, is very weak.Thus, while we cannot completely rule out temperature induced superconducting fluctuations, we can explain all of our findings using Coulomb screening.

V. CONCLUSIONS
To sum up, we connect the energy dependence of the suppression of DOS near Fermi energy in a two-dimensional metal with long-range Coulomb interactions [28,29].Our results show how the Fermi liquid regime is lost in presence of disorder due to Coulomb correlations and that this leads to a fluctuating DOS with a length scale shared by normal and superconducting phases.Future work might address the detailed temperature and field dependence.More detailed measurements at higher magnetic fields might provide further information fluctuations and their role.Furthermore, it is of interest to study the vortex lattice in such a disordered superconductor.
It is also interesting to ask about the possible consequence of our findings to crystalline systems with doping or interface superconductors, in particular those where a pseudogap is observed, as for instance cuprates.Quite generally, our results suggest that a reduction of Coulomb screening can be accounted for within analytical expressions [28,29] and that it is accompanied by spatial fluctuations.It would be good to examine deviations from a Fermi liquid picture in the light of our measurements.For instance, the formation of spatial structures such as stripes or fluctuating charge order in crystalline systems [8,9].The role of disorder induced by dopants, particularly in materials with a small electronic density, might also influence Coulomb interactions in systems that mostly show metallic properties.

APPENDIX A: STRUCTURAL CHARACTERIZATION OF TIN THIN FILMS
We made a characterization of the morphology of the film studied in the main text (D15) by combining x-ray, TEM and STM results.
If there is crystalline order at some length scale, x-ray scattering leads to rings located at Bragg peaks of the lattice constant.Indeed, x-ray scattering of our thin films leads to ring like patterns, which are superposed to the fourfold crystalline Bragg peaks of the Si substrate [Figs.4(a) and 4(b)].We can index all peaks with the planes corresponding to the bulk crystal structure [Fig.4(a)].This shows that growth favors the nucleation of the bulk crystalline structure.
When analyzing the films with TEM, we observe a distribution of many small angle boundaries between small crystals of slightly differing orientations [Fig.4(c)].The corresponding Fourier transform [Fig.4(d)] shows again the Bragg peaks of the substrate and those of the crystalline structure of TiN.Because the field of view is limited, we observe many Bragg peaks, all located at the same distance, instead of the ring that is observed for the whole film in x-ray scattering.Notice the small difference between orientations of crystals.Analyzing the TEM image, we can obtain the distribution of the distance between boundaries.As we show in the inset of [Fig.4(c)], we find a broad distribution, starting at about 2 nm and ending at about 7 nm.The average distance between small angle boundaries is of 5 nm.
In a STM experiment, we can zoom into a small single crystal, as shown in Fig. 4(e).We can observe the atomic lattice.The square lattice is somewhat distorted due to disorder or tip related effects during scanning, but in the Fourier transform, Fig. 4(f), we identify again peaks at the interatomic distances.
Thus, in all, we have a highly disordered system.Scattering centers are unlikely coinciding always with small angle boundaries, so that the relevant length scale can be much larger, as we actually show from our experiments and discuss in the main text.

APPENDIX B: TRANSPORT IN THE SAMPLE OF THE MAIN TEXT AND IN ANOTHER, MORE INSULATING, SAMPLE
It is useful to show and discuss the whole temperature and magnetic field dependence of the resistance.In Fig. 5, we show the result for the sample discussed in the main text.Notice that the upper critical field is of approximately 2.6 T and that the variation of the resistance with the magnetic FIG. 4. In a we show the intensity vs reciprocal distance obtained in x-ray scattering of a thin film of TiN.We index each peaks with their corresponding lattice planes.The space group is given in the legend.The structure is depicted in the upper right part of the figure.In (b), we show the actual x-ray scattering experiment.The blue line provides the scan made to obtain (a).The green circles show the single crystalline Si substrate.We mark with a red dashed line circle the reciprocal size of the lattice constant a, as defined in the upper right inset of (a).In (c), we show a transmission electron microscopy (TEM) image of the TiN film.We observe that the film consists of a TiN crystals separated by small angle grain boundaries.We measured the distance between boundaries and show the result as a histogram in the inset.In (d), we show the Fourier transform of (c).We remark the peaks corresponding to the Si substrate by green circles.The spots coming from the crystallites are distributed on a circle.We mark the reciprocal lattice constant a by red circles.In (e), we show a STM image of a single grain of the TiN film.We observe the atomic lattice.The corresponding Fourier transform is shown in (f).We highlight the reciprocal distance a again with red circles and a line.FIG. 5.In (a), we show the temperature and magnetic field dependence of the sheet resistance of the TiN film discussed in the main text.In (b), we show the same data with the resistance plotted in a color scale.The bar providing the scale for the resistance is given in the right of (b).In (c), we show the resistance vs temperature for different magnetic fields.The color corresponds to the value of the magnetic field, following the bar on the upper right.In the inset, we show the conductance as a function of temperature (normalized to the quantum of conductance G 00 ), in a logarithmic scale.field strongly flattens out for fields above the critical field.The resistance [inset of Fig. 5(c)] shows an increase with decreasing temperature in the normal phase above 4 T.
In Fig. 6, we show the results other film mentioned in the text, D15e, which was obtained from the same film after an additional heat treatment.We observe that the critical temper- FIG. 6.In (a), we show the sheet resistance vs temperature and magnetic field of a thin film with a larger normal phase resistance than the one discussed in the main text.In (b), we show the corresponding colormap, with the resistance given as a color.The scale follows the bar on the right.In (c), we show the resistance as a function of temperature for different magnetic fields, with the magnetic field given by the bar on the bottom right.The top right inset is the normalized conductance G/G 00 vs. temperature at zero magnetic field and at B = 7 T. ) FIG. 7. In (a), we show the temperature evolution of the tunneling conductance G(V ) for a sample which is more resistive than the one discussed in the main text.The sample is termed D15e.We show the data as black circles.The fits to the tunneling conductance G(V ) are given as dashed lines.The DOS according to Bartosh-Kopietz is given by the lines.In (b), we show the derivative of the DOS for different temperatures.
ature and the upper critical field are considerably smaller.At the same time, the resistance shows a stronger increase with decreasing temperature in the normal phase under magnetic fields.
In the sample D15e, we took several tunneling conductance curves.We show the result in Fig. 7.We only show data in the normal phase obtained by applying a magnetic field of 7 T.The superconducting gap opens as a tiny feature in the center of the curve and we could not analyze this in detail.We see that in this sample the tunneling conductance drops to zero at low temperatures much more strongly than in D15.The behavior is also captured by BK model, which provides a temperature variation of N N [Fig.7(b)] leading to a peak at energies that are slightly smaller, but of the same order as in D15.
We summarize some parameters obtained previously in both samples in Table I.
In the Fig. 8, we provide the critical temperature versus sheet resistance in TiN films, from [38].There, we calculated the dependence of T c using ln (

R (kΩ)
FIG. 8. We reproduce data and model of Ref. [38].We show the critical temperature vs the sheet resistance in different TiN samples as black points.The magenta line is the same quantity calculated using Ref. [53]., where r = G 00 R, R is the normal resistance per square, γ = ln[ h/(kT c0 τ )].This expression was derived by Finkel'stein to explain the decrease in T c in disordered thin films due to scattering and electron density fluctuations (Ref.[24]).We see that this expression explains the variation of T c with R. By taking γ = 5.73 we found τ = 7.3 × 10 −15 s, the value we obtain in the main text from our fits to the conductance and from the measured spatial dependence of the DOS fluctuations in the normal phase.

APPENDIX C: CALCULATION OF THE AUTOCORRELATION FUNCTION
In absence of visible patterns, the most efficient way to define a length scale is to calculate the autocorrelation function ACF.This was made in Refs.[11,40] to obtain spatial dependencies related to the opening of the pseudogap in cuprates and to analyze patterns in an oxyde sample.
We should distinguish between the autocorrelation function often used to analyze pictures or images and the statistical ACF we use here.
The usual method starts by defining two arbitrary matrices A and B [54].The cross-correlation matrix C is a measure of similarity between A and B as a function of the displacement of one relative to the other.The elements of C are calculated by displacing A over B a given vectorial lag, then calculating the element to element product of the overlapping elements and finally taking the sum of them.If we take a 2 × 2 matrix , Thus the self correlation of a matrix is just the crosscorrelation of the matrix with itself.However, the statistical ACF measures correlations between two points separated by a distance r.The spatial autocorrelation function ACF(r) of an image is given by the correlation of any two pixels i and j of the image, separated by a given vector r = (r, θ ) = r i − r j , where r i and r j the position of those pixels [11,40].This leads to the following definitions: where δ r,(r i −r j ) , (C2) ) The radial averaged spatial autocorrelation function ACF(r) is the result of averaging the autocorrelation function value for all the vectors with the same magnitude r = |r|.It gives an idea of the spatial extension of the correlated regions in the image.A white noise image has a flat and close to zero ACF(r), with however a sharp peak at r = 0, because every pixel is correlated with itself.In real images, peaks ACF(r) mean that there are regions where pixel intensities are spatially correlated, within the same region or with neighboring regions.
The often used correlation method for images is computationally very fast, since it just needs to "slide, multiply, and sum" the matrix over itself a total of (2d − 1) 2 times, being d the dimension of the matrix.However, it is extremely sensitive to offsets and does not produce the statistical correlation among points exactly.For the calculation of one element, it blindly multiplies and sums neighboring pixels regardless of their relative intensities, orientation or distance.This can cause that two actually correlated elements or regions vanish or are obscured because of negative-plus-positive sums of their surrounding region.Also, it has no renormalization or consideration for the bright center.When images are only slightly displaced they have much more overlapping pixels contributing to the sum than when they are barely touching, making the edge of the image self correlation matrix fainter.
The calculation of the statistical ACF(r) [11,40] by contrast is an offset-independent method, i.e. the result depends only on the relative difference of intensity between two pixels, not on their absolute value.It does not have the drawbacks mentioned in the previous paragraph and provides the connections present between points of a certain distance.The main drawback is of course that it is computationally much more demanding, because it individually checks, counts, multiplies and classifies every possible pair of pixels on the image several times, in order to calculate every parameter needed for the sums.
We can see the results of both methods in Fig. 9.We take a simple vortex lattice image taken in β − Bi 2 Pd from Ref. [55].There is experimental noise in the image and vortices appear as blue circles.The statistical correlation is in essence the size of a vortex.We expect that ACF(r) decreases with r, with a length scale that is the vortex size, and becomes negative at the intervortex distance.
We produce three columns where we have arbitrarily moved the colorscale zero vale.At one side of the histogram containing the values of the pixels [Figs.9(a), 9(d), and 9(g)], at the center of the pixel values [Figs.9(b), 9(e), and 9(h)] and at the average pixel of the histogram [Figs.9(c), 9(f), and 9(i)].
In that way, we see that the result of the calculation is independent of the choice [Figs.9(d), 9(e), and 9(f)] only when we use the statistical ACF(r).As expected, we can determine the vortex core size from the decay at small r, and find a value which coincides with those provided in literature by analyzing the shape of isolated vortices [55].The ACF(r) becomes negative at the intervortex distance and oscillates, with a spatially decaying amplitude.The latter is due to the ratio of regions with zero DOS with respect to those with a finite DOS.The former is particularly large in this system [55].
By contrast, the computationally less demanding method of multiplying matrices only leads to a similar result when the zero of the histogram of the values in image is centered at the average of the histogram [Fig.9(i)].Furthermore, the obtained length scale is strongly distorted.
Of course, a multiplication of 2D images also leads to a 2D image Figs.9(g), 9(h), and 9(i).When using a statistical ACF, the result just depends on the polar coordinates (r, φ).In images that show no in-plane symmetry, as those discussed in the text, there is no angular dependence either and everything is in the radial dependence.When treating a vortex image with a sixfold symmetry, the radial dependence is as shown in the insets of Figs.9(d), 9(e), and 9(f).The angular dependence is given however as a function of a single coordinate, φ.For the purpose of comparison, we have stretched this φ dependence into a 2D matrix in Figs.9(d), 9(e), and 9(f).The relevant parameter is, though, the radial dependence of the correlations [insets of Figs.9(d)-9(i)].
We note again that we normalize ACF in such a way as to obtain a percentage variation of the conductance over its average value.In particular, the variations we see in the main text are between 5% and 10% of the conductance, providing accordingly the y-axis in Fig. 3(e)-3(g).
We provide in Fig. 10 the whole ACF as a function of the position.This is radially averaged in Fig. 3.There are a few remarks in place.First, the central peak is removed, as it is divergent due to the conversion from cartesian to radial coordinates.Second, we should note that the ACF always provides a variation relative to the average conductance value, whatever the latter is.That means that each panel shows the conductance variations with a maximized contrast.Therefore the minimal conductance variations found close to zero bias at zero field Fig. 10(a) are highlighted in such a plot.Third, as we plot the ACF with respect to the average value of the conductance, it can become negative and even oscillate.This is quite apparent at some bias voltage some images in Fig. 10.There is clearly no ordered pattern in the ACF (the images consist of bumps and stripes which vary as a function of the bias voltage, as we see in Fig. 10), although the orientation and level of correlation certainly varies in different positions.
It will be interesting to analyze in future the statistics of such a behavior over very large areas and many images, as it might provide further clues about the deviation from random behavior induced by electronic correlations.

APPENDIX D: ASPECTS OF STM IN HIGH RESISTANCE FILMS
Let us remind a few basic aspects of tunneling STM experiments.In a usual STM experiment, there is often a small resistance in series with the tunnel junction.This helps controlling the noise level by producing low-pass filters together with capacitors.The tunneling current is given by where R is the resistance in series, V 0 the voltage drop in that resistance, V 1 the voltage drop at the tip-sample junction and σ Junction the actual conductance of the junction.Usually, 1 σ Junction R, so that V 0 ≈ 0 and V ≈ V 1 and we can write the widely used relation between the tunneling current in a junction and the DOS: Its derivative is σ Junction (V ) = dI dV is zero for V < , with being the superconducting gap, finite V > and eventually diverges exactly at the quasiparticle peaks for a conventional s-wave BCS superconductor.In our case, R ≈ 20 k , and 1 σ Junction is well above a M , so that the condition 1 σ Junction R is maintained for every purpose during the voltage ramp.However, the current flows through the sample before reaching the tip-sample junction.Thus, there is an additional resistance R S which adds to R due to the sample.It is important to see that the value of R S is the one found in macroscopic transport experiments and does not considerably modify R. For example, a sample resistance of 10 k with a current flow of 1 nA leads to a voltage drop of 10 μV.Thus, the effect of R + R S is of at most a voltage shift of a few tens of μV .If we assume a non-Ohmic R S , the voltage drop in the sample might become larger.However, by varying the tip-sample distance we can modify σ Junction , and thus the relative role of R S in the tunneling experiment.With a voltage drop in the sample, we expect to see shifts of features in the DOS of the sample.We do not observe such shifts, showing that R S is ohmic.
On the other hand, a DOS strongly varying with energy has further consequences that we need to discuss.To see this, we remind that the tip is positioned through a feedback mechanism that maintains a constant tunneling current at a bias of several mV, which is usually of order of a nA.The feedback mechanism thus imposes a relation between I t (V ) which can be written as where V B is the bias voltage at which the control system of the STM is working [57].Thus, when modifying the bias voltage is accordingly modified too.However, normalizing the current or its derivative at, say V 0 < V B , eliminates this factor and provides curves that are comparable to each other with different bias voltages.
In the case we consider here, I t (V 1 ) is linked to the DOS, but we measure I t (V ).The feedback mechanism acts in the same way, with modified normalization constants, which are eliminated by normalizing the current to I at another bias voltage, I t (V 0 ).The same applies for the tunneling conductance G(V ).Although the influence of the feedback loop becomes more important in more resistive samples (as the heated D15), normalization at a fixed voltage can be used to obtain results that can be compared among samples and measurement conditions.

APPENDIX E: MULTIFRACTAL ANALYSIS OF IMAGES
The distribution function of multifractal exponents F (α) provides the dimensions α that are found in a given field of view and is often used in characterization of 2D images.We describe our method in detail in Appendix F of Ref. [58] and the code is available at Ref. [59].Here we provide a few details as a guide.We first remind the properties of F (α). F (α) consists of just a point at α = 2 for a fully random image.For an image containing a monofractal distribution of points, F (α) is just a point at α equal to the fractal dimension.In a multifractal distribution F (α) is an inverted parabola with a maximum F (α max ) at α max .The parabola is defined around a range of α, as shown in Fig. 11, where we show F (α) in our images.
We the counting method to obtain F (α), described in Ref. [60].We first calculate the sum of all conductances in a box of length , normalized to the sum over all possible boxes, P i ( ).We vary the box size between 64 × 64 and 1 × 1.Then we introduce a set of exponents q and calculate the sum over all possible P i ( ) q , where i is a counting index (we use −10 < q < 10).This has an exponential dependence as a function of the box size with a set of exponents τ (q), which can be found from the logarithmic derivative of the sum over the box size over q.The set of dimensions α present in the patterns forming the image are distributed with 12.We plot as points the normalized tunneling conductance vs bias voltage at zero field (red points) and under magnetic fields of 4 T (light blue points) and 7 T (blue points).Data are taken at 0.3 K.We provide as lines the fits to BK model under magnetic fields.At zero field, we use the same fit, numerically modified as described in the text.a function F (α) which is obtained from τ (Q), α(Q) and the said logarithmic derivative as explained in Refs.[58,60].
Thus the x axis in Fig. 11 provides the set of dimensions that are found in the image and the most probable dimension occurs at the maximum of F (α). Thus the position of the maximum F (α) at α max deviates from 2 as the image deviates from a random distribution.In Fig. 3(h), we plot α max − 2.

APPENDIX F: MAGNETIC FIELD DEPENDENCE
In Fig. 12, we compare tunneling spectroscopy data obtained at zero field with data obtained above the critical field.We can see that the magnetic field induced modifications in the normal phase above H c2 are small.This is not surprising.As we can see in Fig. 5, the resistance varies only slightly with magnetic field at low temperatures in the normal phase above H c2 .We also include the fit discussed in the text.To account for the zero field data, we use the fit obtained under magnetic fields multiplied by an arbitrary function σ (V ) which is chosen to maintain the same area below the conductance vs voltage curves at zero field and under magnetic fields.

3 FIG. 1 .
FIG. 1. (a)The tunneling conductance vs. voltage normalized to its value at 5 mV is shown as a color scale (scale bar on the right) and is plotted as a function of temperature at a magnetic field of 7 T. (b) Sheet resistance vs temperature in the same sample at the same magnetic field (orange line).The result at zero magnetic field is shown by the black dashed line.Vertical dashed line provides T c at zero magnetic field.The magnetic field is applied perpendicular to the film, well above the critical field.Inset in (b) shows the resistance vs 1/T in a log scale.All resistance vs temperature measurements are provided in Appendix B, the sample is named D15.Well above H c2 (of about 2 T at zero temperature), i.e., above about 4 T, the resistance vs temperature and the tunneling conductance shows a weak magnetic field dependence.

FIG. 3 .
FIG.3.In (a) and (b), we show the tunneling conductance normalized to its value at 1 mV as a function of the position along a line (black arrow) at zero field (a) and with a magnetic field of 4 T (b).Data are taken at 0.3 K.The color scale provides the normalized conductance following the scale bar in the right of each panel.In (c) and (d), we show tunneling conductance maps over an area of 400 nm at zero magnetic field (c) and with a magnetic field of 4 T (d).The field of view is different, due to a modification of the position when applying the magnetic field.We show maps at the bias voltages 0, 0.15, and 0.3 mV, from left to right.The normalized conductance is shown by the bars in the right side of each figure.The arrows provide the scan over which we trace the tunneling conductance in (a) and (b).Note that we make one tunneling conductance curve approximately each 5 nm.In (e) and (f), we show the ACF(r) at zero field and at 4 T. The amplitude of the ACF is given by the color bar on the right of each figure and has been normalized to show directly the spatial variations in the conductance.In (g), we show ACF(r) averaged for all bias voltages at zero magnetic field (green squares) and at 4 T (red circles).In (h), we show the position of the maximum of the distribution function of multifractal exponents F (α) as a function of the bias voltage.We mark by a black vertical dashed line the position of the quasiparticle peaks and by a blue dashed line the position of the maximum in F (α).

a 21 a
22 and a 3 × 3 matrix B =

FIG. 9 .
FIG. 9.In (a)-(c), we show a test image with a vortex lattice, from Ref. [56].The color scale is given by the scale bar at the right in arbitrary units.The inset shows the color scale (from blue to yellow) and in white the color histogram from the image.The red line shows the position of zero in the image.In (a), the bottom of the scale is given the zero value, in (b) the center of the scale and in (c) the position of the average pixel value.In (d)-(f), we show the result of the statistical ACF(r, φ) following the model described in the text for each of the images in (a)-(c).The bar in the right of each image provides the amplitude of the ACF(r).In the insets we show ACF(r).We see that all figures (d)-(f) are identical.In (g)-(h), we show the correlation function obtained applying usual image analysis software to the images (a)-(c).Insets provide the radial average and the color bar at the right the amplitude.The result is strongly dependent on the average value of the color scale of the image.Black scale bars are 100 nm long.

FIG. 10 .
FIG. 10.We show the ACF as a function of the position in (a)-(f) for different bias voltages.The color scale is give by the scale bars on the right.Bias voltages and magnetic fields are given in each panel.We obtain from these images the radial average discussed in the Figs.3(e) and 3(f).We show the radial average of each image in the lower right insets.Dashed vertical line in the insets and circle in the main panels provide the distance discussed in the text, 20 nm, as a guide.
FIG.11.In (a), we show as lines the distribution fractal at zero magnetic field and in (b), at a magnetic field of 4 T. The bias voltage dependence is given by the color scale on the right.Random behavior produces a point at α = 2 and F (α = 2) = 2. Deviations from random behavior increase the range of F (α) in α and displace the maximum of F (α) from 2. Notice that the distribution at zero field is close to random at low and high bias voltage.The deviations from random behavior slightly increase when decreasing bias voltage at 4 T.In Fig.3, we plot the position of the maximum in F (α), α max .
FIG.12.We plot as points the normalized tunneling conductance vs bias voltage at zero field (red points) and under magnetic fields of 4 T (light blue points) and 7 T (blue points).Data are taken at 0.3 K.We provide as lines the fits to BK model under magnetic fields.At zero field, we use the same fit, numerically modified as described in the text.

TABLE I .
The sample D15e was made from a sample identical to D15 by heating.R is the resistance of the film per square at room temperature; B c2 (0) is the upper critical field at T = 0; and