Spectroscopic evidence for strong correlations between local superconducting gap and local Altshuler-Aronov density-of-states suppression in ultrathin NbN films

TL;DR

This study uses scanning tunneling spectroscopy to reveal that local superconducting gap variations in ultrathin NbN films are strongly linked to Altshuler-Aronov density-of-states suppression caused by electronic interactions, explained through a local Finkelstein mechanism.

Contribution

It demonstrates a purely fermionic explanation for superconducting inhomogeneities in disordered NbN films based on local LDOS suppression and resistivity variations.

Findings

Local superconducting gap correlates with LDOS depletion.

LDOS suppression can be modeled as local resistivity changes.

Resistivity variations induce local gap variations via Finkelstein mechanism.

Abstract

Disorder has different profound effects on superconducting thin films. For a large variety of materials, increasing disorder reduces electronic screening which enhances electron-electron repulsion. These fermionic effects lead to a mechanism described by Finkelstein: when disorder combined to electron-electron interactions increases, there is a global decrease of the superconducting energy gap $\Delta$ and of the critical temperature $T_c$, the ratio $\Delta$/$k_BT_c$ remaining roughly constant. In addition, in most films an emergent granularity develops with increasing disorder and results in the formation of inhomogeneous superconducting puddles. These gap inhomogeneities are usually accompanied by the development of bosonic features: a pseudogap develops above the critical temperature $T_c$ and the energy gap $\Delta$ starts decoupling from $T_c$. Thus the mechanism(s) driving the appearance of these gap inhomogeneities could result from a complicated interplay between fermionic and bosonic effects. By studying the local electronic properties of a NbN film with scanning tunneling spectroscopy (STS) we show that the inhomogeneous spatial distribution of $\Delta$ is locally strongly correlated to a large depletion in the local density of states (LDOS) around the Fermi level, associated to the Altshuler-Aronov effect induced by strong electronic interactions. By modelling quantitatively the measured LDOS suppression, we show that the latter can be interpreted as local variations of the film resistivity. This local change in resistivity leads to a local variation of $\Delta$ through a local Finkelstein mechanism. Our analysis furnishes a purely fermionic scenario explaining quantitatively the emergent superconducting inhomogeneities, while the precise origin of the latter remained unclear up to now.

Equations

$$\alpha (V)=\frac{R_{\square }}{2\pi R_Q}\ln \frac{\hslash {\omega }_0}{eV},$$

$$\delta \alpha (\mathbf{r},V)=\frac{\delta \rho (\mathbf{r})}{2\pi R_Q}\ln \frac{eV}{\hslash D/l_{\alpha }^2}.$$

$$\frac{\delta \Delta (\mathbf{r})}{⟨\Delta ⟩}\approx -\frac{\delta \rho (\mathbf{r})}{6\pi R_Q}{\ln }^3\frac{\hslash {\omega }_D}{\Delta (0)},$$

$${\rho }_{\Delta }(\mathbf{r})=\sum _{{\mathbf{r}}^{\prime }}\frac{(\Delta ({\mathbf{r}}^{\prime })-⟨\Delta ⟩)(\Delta ({\mathbf{r}}^{\prime }-\mathbf{r})-⟨\Delta ⟩)}{N{\sigma }_{\Delta }^2},$$

$${\rho }_{\text{cross}}(\mathbf{r})=\sum _{{\mathbf{r}}^{\prime }}\frac{(\Delta ({\mathbf{r}}^{\prime })-⟨\Delta ⟩)(\alpha ({\mathbf{r}}^{\prime }-\mathbf{r})-⟨\alpha ⟩)}{N{\sigma }_{\Delta }{\sigma }_{\alpha }}.$$

$$S(E)=\frac{2}{R_Q}{\int }_E^{1/\tau }\frac{d\omega }{\omega }R(\omega ),$$

$$R(\omega )=\frac{R_{\square }}{2\pi }\ln \frac{r_{\text{out}}(\omega )}{r_{\text{in}}(\omega )}=\frac{R_{\square }}{4\pi }\ln \frac{{\omega }_0}{\omega }.$$

$$S_{\text{homo}}(E)=\frac{R_{\square }}{4\pi R_Q}\left({\ln }^2\frac{{\omega }_0}{E}-{\ln }^2{\omega }_0\tau \right).$$

$$R(\mathbf{r},\omega )=\frac{\rho (\mathbf{r})}{2\pi }\ln \frac{a}{r_{\text{in}}(\omega )}+\frac{R_{\square }}{2\pi }\ln \frac{r_{\text{out}}(\omega )}{a}$$

$$S(\mathbf{r},E)=S_{\text{homo}}(E)+\frac{\delta \rho (\mathbf{r})}{4\pi R_Q}\left({\ln }^2\frac{1/\tau }{D/a^2}-{\ln }^2\frac{E}{D/a^2}\right),$$

$$\delta {\lambda }_{\epsilon ,{\epsilon }^{\prime }}=-\frac{1}{2\pi \nu }\int (dq)\frac{1}{D{q}^2+\epsilon +{\epsilon }^{\prime }},$$

$$\int (dq)f(q^2)\equiv \frac{1}{d}\int \frac{d{q}_{x}d{q}_y}{(2\pi {)}^2}\sum _{n=0}^{\infty }f[q_x^2+q_y^2+(\pi n/d{)}^2].$$

$$\frac{\delta \Delta }{{\Delta }_0}={\int }_{\Delta }^{{\omega }_D}\frac{d\epsilon }{\epsilon }\frac{d{\epsilon }^{\prime }}{{\epsilon }^{\prime }}\delta {\lambda }_{\epsilon ,{\epsilon }^{\prime }},$$

$$\delta {\lambda }_{\epsilon ,{\epsilon }^{\prime }}\approx \delta {\lambda }_{\epsilon ,{\epsilon }^{\prime }}^{\text{low}}+\delta {\lambda }^{\text{high}},$$

$$\delta {\lambda }_{\epsilon ,{\epsilon }^{\prime }}^{\text{low}}=-\frac{1}{2\pi \nu }\int (dq)\frac{\theta (q_D-q)}{D{q}^2+\epsilon +{\epsilon }^{\prime }},$$

$$\delta {\lambda }^{\text{high}}=-\frac{1}{2\pi \nu }\int (dq)\frac{\theta (q-q_D)}{D{q}^2},$$

$${\Delta }_1={\Delta }_0\exp [1/\lambda -1/(\lambda +{\lambda }^{\text{high}})].$$

$$\delta {\lambda }_{\epsilon ,{\epsilon }^{\prime }}^{\text{low}}=-\frac{\rho }{2\pi R_Q}\ln \frac{{\omega }_D}{\max (\epsilon ,{\epsilon }^{\prime })}$$

$$\frac{\delta {\Delta }^{\text{low}}}{{\Delta }_1}=-\frac{\rho }{6\pi R_Q}{\ln }^3\frac{{\omega }_D}{\Delta }.$$

$$\frac{\delta \Delta (\mathbf{r})}{⟨\Delta ⟩}=-\frac{\delta \rho (\mathbf{r})}{6\pi R_Q}{\ln }^3\frac{{\omega }_D}{\Delta }+\delta {\lambda }^{\text{high}}(\mathbf{r}){\ln }^2\frac{{\omega }_D}{\Delta },$$

$$\delta {\lambda }^{\text{high}}(\mathbf{r})=-\frac{2d\delta \rho (\mathbf{r})}{R_Q}\int (dq)\frac{\theta (q-q_D)}{q^2}.$$

Full text

I Introduction

Non-magnetic disorder was initially thought to have little effect on the superconducting properties of conventional thin films, as far as time-reversal symetry is preserved Anderson ; Abrikosov . Decades of intense experimental and theoretical works have ended to the opposite conclusion showing that beyond a critical disorder all systems transit either to a metallic or to an insulating state Gantmakher ; Feigelman . For thin films that are not single-crystals Brun2017 but consist of coupled nanocrystals two different classes of systems were historically considered. A qualitatively different superconducting behavior was observed between so-called granular and homogeneous disordered thin films. In granular ones a poor electrical coupling between the nanocrystals makes the film a disordered array of Josephson junctions Efetov ; Orr ; White ; Jaeger ; Gantmakher . Their superconducting properties are controlled by weak Josephson couplings between phases of local superconducting order parameters attributed to individual nanocrystals. In homogeneous ones the nanocrystals are much better electrically coupled to each other. Their superconducting properties are controlled by the subtle interplay between the non-magnetic disorder distribution and electron-electron interactions Graybeal ; Dynes ; Haviland ; Feigelman .

Nevertheless experimental and theoretical studies performed in the last ten years enabled to show that this historical categorization between granular and homogeneous thin films is lacking of an important link between them. At low disorder in homogeneously disordered films macroscopic measurements show that $T_{c}$ and $\Delta$ monotonously decrease with rising disorder corresponding to $k_{F}l$ decreasing from a value much larger than one toward unity Graybeal ; Dynes ; Haviland ; Feigelman ($k_{F}$ being the Fermi wavevector and $l$ the elastic electron mean-free path). This effect has been satisfactorily explained by a reduction of the electronic screening upon increasing disorder which reduces the diffusion constant $D$ and is called the “fermionic” mechanism Finkelstein ; Finkelstein1994 . Many films like, for example, Bi, Pb, MoGe, NbSi behave this way and follow the so-called Finkelstein scenario down to the almost complete destruction of superconductivity Dynes ; Valles ; Graybeal ; Crauste1 ; Crauste2 .

However a different behavior can also take place for sufficiently large disorder in nominally homogeneous films, i.e. typically for $1\leq k_{F}l\leq 3$: in amorphous InOx, TiN or NbN the single-particle gap probed by tunneling remains finite and decouples from the decreasing $T_{c}$, while a pseudogap emerges in a significant temperature range above $T_{c}$ Sacepe2008 ; Sacepe2010 ; Sacepe2011 ; Chand2012 . These latter features are believed to correspond to the so-called “bosonic” scenario, where disorder localizes single-electron wavefunctions and provoke an “emergent granularity” in the superconductor Ghosal ; Feigelman ; Bouadim . These emergent superconducting puddles were revealed by scanning tunneling microscopy/spectroscopy (STM/STS) in TiN, InOx and NbN Sacepe2008 ; Sacepe2011 ; Chand2012 ; Noat ; Kamlapure . From moderate to strong disorder, the long-range phase coherence between the emergent neighboring superconducting islands gradually weakens and eventually destroys, leaving the material with residual incoherent pairing properties Feigel2007 ; Feigel2010 .

From the above considerations one sees that in fact most films are expected to present an interplay between “fermionic” and “bosonic” effects for strong enough disorder before superconductivity disappears. Nevertheless there is a lack of experimental data combining macroscopic and local measurements on the same system to carefully explore this interplay. Such a detailed study was conducted in NbN Chand2012 and enabled to draw a phase diagram, where at low disorder the fermionic mechanism dominates, while bosonic effects develop at stronger disorder, as exposed above.

In this work we show that, in contrast to the common belief, the emergent gap inhomogeneities can be solely explained by the Finkelstein mechanism, where electron interactions play a crucial role. Our claim is supported by the fact that using scanning tunneling microscopy/spectroscopy in NbN thin films, we found a direct spatial cross-correlation between local superconducting gap maps and the mapping of the Altshuler-Aronov suppression of the local density-of-states (LDOS). Through appropriate modeling we demonstrate that both effects can be linked to local variations of the film’s resistance. Due to an interplay between disorder and electrons interaction, these resistivity variations drive corresponding spatial variations of the normal-state STS spectra and the local energy gap values.

II Experimental results: local STS measurements of the gap and V-shape LDOS

Our samples consist of ultrathin NbN films grown ex-situ on sapphire, structured in nano-crystals of lateral size 2–5 nm Semenov . Special care was taken in order to minimize the thickness of the surface oxyde layer: the freshly grown samples were introduced into our STM set-up in less than 24 hours after their growth. We selected films having a $T_{c}\approx 0.25\,T_{c}^{\text{bulk}}=3.8$ K, so that seeming “bosonic” properties like gap inhomogeneities and pseudogap features have already developed. The nominal film thickness is about 2.1 nm. We studied their local superconducting properties by scanning tunneling microscopy/spectroscopy (STM/STS) at 300 mK, using PtIr tips, in an ultrahigh vacuum homemade setup. The presented $dI/dV$ measurements were obtained by numerical derivation of single $I(V)$ curves. $T_{c}$ was extracted from in situ 4-points electrical resistivity mesurements performed in the STM stage during the same run as the STM/STS measurements. The temperature dependence of the resitivity before the superconducting transition is typical of 2D disordered superconducting films and is shown in Fig. 1a.

STM topography measurements show that the film surface is very flat (see Fig. 2a). The size of the observed nanoscale structures correlates well with the size of the nanocrystals characterized independently by transmission electron microscopy (TEM) experiments. A characteristic $dI/dV$ spectrum measured locally is shown in Fig. 1b. It presents a fully gapped LDOS with well-defined superconducting coherence peaks. We also note that in-gap states are also present, most probably induced by paramagnetic defects at the NbN/sapphire interface. For energies larger than the ones of the coherence peaks, a strong background is seen instead of recovering a standard normal DOS, as in the Bardeen-Cooper-Schrieffer (BCS) case. Fig. 1c enables characterizing this background: it presents a $dI/dV$ spectrum measured on the same location but over a wider energy range about 20 times larger than $\Delta$. We see that the tunneling LDOS is strongly reduced over an energy range much larger than $\Delta$. This LDOS suppression is typical of systems where electronic correlations combined to disorder reduce the tunneling DOS at the Fermi level ($E_{F}$) Altshuler ; Altshuler1985 , an effect first theoretically explained by Altshuler and Aronov. In the following we refer to this effect as AA.

Fig. 1c shows that for $e|V|>\Delta$ the measured LDOS follows a power-law as a function of the applied bias voltage, $dI/dV\propto V^{\alpha}$. Please note that here the fit is truly performed using the functional dependence $V^{\alpha}$ without any offset: this curve goes to zero at zero voltage. Nevertheless since the fit is performed in the energy region outside the gap, data points exactly around zero voltages are not plotted in this figure. Such a behavior has also being evidenced for instance in highly disordered metallic thin films close to the insulating transition Richardella .

Theoretically such a dependence is expected in the framework of nonperturbative extension Ingold ; LevitovShytov ; Rollbuhler of the Altshuler-Aronov theory Altshuler ; Altshuler1985 . Let us note that this kind of effect is also referred to as dynamical Coulomb blockade effect in mesoscopic physics Ingold ; LevitovShytov ; Rollbuhler . Interestingly, it was shown that the effect of probing the tunneling DOS of an ultrasmall tunnel junction in a 2D system with electron-electron interaction is equivalent to dynamical Coulomb blockade formalism with an ohmic environmentRollbuhler . We have already used such an equivalence when modeling the superconducting proximity effect between a Pb single nano-crystal and a 2D metallic homogeneously disordered Pb film Brun2012 ; Serrier ; Roditchev .

The main finding of our work is a systematic spatial correlation between the local energy gap value $\Delta(\mathbf{r})$ and the strength of the AA LDOS suppression, quantified by the exponent $\alpha(\mathbf{r})$ in the power law $dI/dV(\mathbf{r})\propto V^{\alpha(\mathbf{r})}$. At each square nanometer of the $300\times 250$ nm2 area shown in Fig. 2a, the local gap value was extracted from the peak-to-peak distance in each $dI/dV(\mathbf{r})$ curve. The extracted map of the gap values is shown in Fig. 2b. The energy gap ranges between 1 and about 1.5 meV. The histogram of the gap values is shown in panel 2e and is close to a Gaussian distribution. The spatial distribution of $\Delta(\mathbf{r})$ shows emerging local inhomogeneities seen as lighter and darker puddles. Using the scale bar shown in the topography Fig. 2a, one sees that local puddles having a constant gap identified by a certain local color are smaller than 50 nm. Such a behavior was reported previously in similar systems for comparable disorder strength, but with comparatively less clean and less dense STS spectraSacepe2008 ; Noat ; Kamlapure ; Carbillet . $\Delta=1.2$ meV is used below as the average gap value of our film.

We call supergrain a local puddle with a nearly constant gap. We define $l_{\Delta}$ as the characteristic size of such supergrain. We found $l_{\Delta}\approx 27$ nm which is much larger than the low-$T$ coherence length $\xi(0)\approx 5$ nm. $l_{\Delta}$ corresponds to a correlation length characterizing the spatial extent over which local gap values are correlated. A precise way to extract it is to use the spatial dependence of the radial averaging $\rho_{\Delta}(r)$ of the 2D autocorrelation function of $\Delta(\mathbf{r})$ (details about the numerical computations of the correlation and cross-correlation functions can be found in the Appendices section). The spatial decay of $\rho_{\Delta}(r)$ and $l_{\Delta}$ are shown in Fig. 3b.

We provide now for the first time a local mapping of the disorder distribution combined to electron interactions which together are responsible for the spatial distribution of the gap inhomogeneities. At each location the gap value $\Delta(\mathbf{r})$ was measured, the local AA exponent $\alpha(\mathbf{r})$ was also extracted in the energy range [5;30] mV, separately for positive and negative energies. We then took for $\alpha(\mathbf{r})$ the average between the values fitted for positive and negative energies. The slight asymmetry in our spectra between positive and negative energies probably comes from a non-constant tip DOS.

The 2D map of the exponent $\alpha(\mathbf{r})$ is presented in panel 2c and its histogram in panel 2f. Most of the values are situated between 0.26 and 0.43. It is seen that the exponent map also presents emerging granularity, with some resemblance with the gap map. Nevertheless the grainy structure of the exponent map looks sharper than the one of the gap map. This suggests that the correlation length $l_{\alpha}$ of $\alpha(\mathbf{r})$ should be smaller than $l_{\Delta}$. The correlation length $l_{\alpha}$ is extracted following a similar procedure as for $\Delta(r)$ above. The spatial dependence of the radial averaging $\rho_{\alpha}(r)$ of the 2D autocorrelation function of $\alpha(\mathbf{r})$ is shown in Fig. 3a. $l_{\alpha}$ also is indicated and yields $l_{\alpha}\approx 18$ nm. This value is indeed about two thirds of $l_{\Delta}$. In order to simplify the visual comparison between the gap map and the exponent map, we used the same color scale for both maps but with a reversed order, which makes it easier to see anti-correlations. This means that if statistically anti-correlations dominate between $\Delta(\mathbf{r})$ and $\alpha(\mathbf{r})$ one expects to see at the same place $\mathbf{r}$ bright puddles in both maps: where a large gap exists (bright puddle) a low $\alpha$ value should also be found (bright puddle). The situation is similar for darker puddles. Looking precisely at both maps Fig. 2b and Fig. 2c and trying to compare the shape and color of local puddles, one sees that anti-correlations are indeed present.

In order to mathematically quantify the cross-correlations existing between $\Delta(\mathbf{r})$ and $\alpha(\mathbf{r})$ we computed their 2D normalized cross-correlation function, which is shown in panel 2g. A very large anti-correlation of about $-0.55$ at coincident points is found. This confirms that the larger is the local energy gap value $\Delta(\mathbf{r})$, the smaller is the local AA exponent $\alpha(\mathbf{r})$ and vice-versa, as illustrated on panels 2h and i. The $dI/dV(\mathbf{r})$ spectrum with a large $\Delta(\mathbf{r})$ has a flatter AA background (red curve), and vice versa (blue curve). To extract the spatial extent over which local cross-correlations exist, we calculated the radial averaging $\rho_{\text{cross}}(r)$ of the 2D cross-correlation function presented in Fig. 2g. Its spatial dependence is shown on panel Fig. 3c. One sees that the typical correlation length of the cross-correlations is $l_{\text{cross}}\approx 25$ nm. This is larger than $l_{\alpha}$ and close to $l_{\Delta}$.

To show the robustness of our results we present in Fig. 4 a linecut of twelve $dI/dV(\mathbf{r})$ spectra measured along the line indicated by a white rectangle in Fig. 2c. This linecut starts at a bright puddle with larger local gap and smaller exponent. It corresponds to the bottom $dI/dV$ spectrum. 33 nm away this linecut ends at a dark puddle with smaller gap and larger exponent. The end corresponds to the top $dI/dV$ spectrum. Each spectrum is spaced by 3 nm from the previous one. It is seen that each individual $dI/dV(\mathbf{r})$ curve is well-fitted by a power law $dI/dV(\mathbf{r})=b\times V^{\alpha(\mathbf{r})}$, where $b$ is a constant.

One could naively think that a direct one-to-one correspondence could exist between the spatial distribution of the NbN nano-crystals encoded in the STM topography and the gap map or AA exponent map. In fact we found no cross-correlation between the local gap or exponent values and the topography of the probed area (see panel 2d). This is in agreement with our previous work Carbillet . It could be partly due to the inhomogeneous oxidation of the NbN surface if this process would change the one-to-one correspondence between the NbN crystallites distribution and the STM topography. This nevertheless seems unlikely regarding the size of the nanostructures seen in STM which agree well with the size of the nanocrystals measured from TEM experiments.

In contrast, the present results furnish a new way of characterizing locally and quantitatively the underlying disorder distribution combined to electron interaction effects which lead to the appearance of supergrains of size $l_{\Delta}$. Indeed our results provide a new type of information that was not reported before, to our knowledge. For instance in Refs. Sacepe2011, or Kamlapure, , the focus of the STS studies is set on the comparison between local gap values and local peak height values.

We finish this experimental section by giving the numerical values of relevant microscopic parameters that will be used hereafter in our theoretical analysis. The low-temperature coherence length can be estimated as $\xi(0)=\sqrt{\hbar D/\Delta}\approx 5.2$ nm. It corresponds to the diffusion coefficient $D\approx 0.5$ cm2/s obtained from $dB_{c2}(T)/dT$ data close to $T_{c}$, see Ref. Semenov, . It is also consistent with our own electrical resistivity measurements in perpendicular high-magnetic field Carbillet on similar films. Taking for the elastic mean free path $l\approx 0.5$ nm Semenov , we find for the elastic scattering time: $\tau=l^{2}/3D\approx 1.7\times 10^{-15}$ s.

III Theoretical modeling

Our theoretical explanation of the observed results is based on the standard understanding of:

• (i)

the Altshuler-Aronov anomaly Altshuler ; Altshuler1985 , properly generalized beyond perturbation theory for both weak and strong coupling LevitovShytov ; KamenevAndreev ;

• (ii)

the disorder-induced suppression of $T_{c}$ in superconducting films Finkelstein ; Finkelstein1994 .

Both effects have been described for a spatially homogeneous case. Both have thus to be properly generalized to an inhomogeneous situation, since this corresponds to the case of our experimental results.

The main physical mechanism leading to the suppression of $T_{c}$ in dirty superconductors is through the enhancement of Coulomb repulsion, which in turn reduces the effective attraction in the Cooper channel. An important distinction should here be made between the case of 2D and 3D superconductors. In the 3D case, the whole effect can be simply described in terms of the Coulomb pseudopotential AMR . However, our NbN ultrathin film is a 2D superconductor. Here, in contrast to the 3D situation, one cannot reduce the effect of Coulomb interaction to the modification of just a single Cooper coupling constant. In 2D one should instead track the whole energy dependence of the effective Cooper coupling, which changes considerably within the relevant energy scale of $T_{c}$. This requires a more sophisticated analysis, which can be done via a renormalization group technique developed by Finkelstein Finkelstein ; Finkelstein1994 .

In 2D, the effects (i) and (ii) listed above are essentially determined by the film resistivity. Since we want to describe an inhomogeneous situation for these two effects, we assume a model, where the local resistivity $\rho(\mathbf{r})$ fluctuates with the correlation length $l_{\alpha}$. With particular details relegated to Appendices, we summarize below the main findings.

III.1 Link between the tunneling LDOS suppression and the local resistance

According to Refs. Altshuler, ; Altshuler1985, ; Ingold, ; LevitovShytov, ; KamenevAndreev, the tunneling DOS suppression in an interacting diffusive normal metal, $\nu(E)$, is determined by the spreading resistance between the diffusive scale and the field propagation scale (see Appendix B.1). Relevant $eV$ values (dozens of mV) are much smaller than $\hbar/\tau\approx 0.33$ eV, justifying well the diffusive regime. In the 2D geometry, the log-normal dependence of $\nu(E)$ LevitovShytov ; KamenevAndreev can be approximated by a power law with a weakly energy-dependent exponent

[TABLE]

where $R_{Q}=h/e^{2}=25.8$ k$\Omega$ is the resistance quantum and $\omega_{0}$ is the plasma frequency. Taking $\hbar\omega_{0}\approx 10$ eV from Ref. Semenov, , $V=10$ mV (typical voltage scale for experimental determination of $\alpha$, see Fig. 1c) and $R_{\Box}=6.8$ k$\Omega$ (maximal sheet resistance $R_{\Box}^{\text{max}}$ before the transition, see Fig. 1a), we obtain $\alpha=0.29$. This value is quite close to the mean value 0.35 measured experimentally. A small discrepancy could be related to the use of an underestimated $R_{\Box}^{\text{max}}$ that is significantly affected by thermal amplitude superconducting fluctuations, whereas the Altshuler-Aronov effect involves larger frequencies, where this suppression is less efficient. We conclude that the mean value of the AA exponent $\langle\alpha(\mathbf{r})\rangle$ is well described by the theory developed in Refs. LevitovShytov, ; KamenevAndreev, .

The key point in our explanation of inhomogeneity seen in $\alpha(\mathbf{r})$ and $\Delta(\mathbf{r})$ is an assumption that they both originate from fluctuations of the local 2D resistivity $\rho(\mathbf{r})=R_{\Box}+\delta\rho(\mathbf{r})$. Using the concept of the spreading resistance to calculate $\delta\alpha(\mathbf{r},V)$, we obtain (see Appendix B.2):

[TABLE]

Using the experimentally measured dispersion of the AA exponent $\sigma_{\alpha}=0.028$ extracted from panel 2f, its correlation length $l_{\alpha}\approx 18$ nm and taking $V=10$ mV, we find $eV/(\hbar D/l_{\alpha}^{2})\approx 70$ and hence obtain the sheet resistance dispersion $\sigma_{\rho}=1.1$ k$\Omega$. This value corresponds to 16 % of $R_{\Box}^{\text{max}}$. This dispersion value among local resistivities seems quite reasonable regarding the moderate but not too strong disorder existing in our NbN thin films.

We would like here to comment on the justification and meaning associated with the definition of a local resistivity. The correlation length $l_{\alpha}\approx 18$ nm of the fluctuations of $\alpha$ exceeds well the electronic mean free path $l\approx 0.5$ nm. This implies that our sample is deep in the diffusive regime for the considered energy scale corresponding to $\alpha$ and its variations, thus the local resistivity can be defined. Since the spatial fluctuations of the AA exponent are mainly controlled by the diffusive length scale (the diffusive length $L$ is related to the energy $E$ through the relation $L=\sqrt{\hbar D/E}$), these fluctuations enable us to probe the local film resistivity.

III.2 Link between the local gap inhomogeneities and the local resistances

The additional suppression of the superconducting order parameter in regions with larger local resistivity is explained by the Finkelstein mechanism. In general, one should distinguish between spatial fluctuations coming from small ($q<q_{D}$) and large ($q>q_{D}$) momenta, where $q_{D}=\sqrt{\omega_{D}/D}$ with $\omega_{D}$ the Debye frequency. For our NbN film we estimate $\hbar\omega_{D}\approx 300$ K, $q_{D}\approx 1$ nm*-1*. The low-momentum contribution is given by the usual log-cube 2D expression

[TABLE]

cut off at the Debye frequency (see Appendix C.1). The high-momentum contribution can be expressed in terms of the renormalization of the BCS coupling constant $\lambda(\textbf{r})$. Substituting the resistivity dispersion $\sigma_{\rho}=1.1$ k$\Omega$ into Eq. (3), we obtain for the relative gap dispersion $\sigma_{\Delta}/\left<\Delta\right>=0.074$, whereas the experimental value is 0.057 (see panel 2e). Note that Eq. (3) contains a large degree of uncertainty: the precise value of $\hbar\omega_{D}$ being unknown, $\sigma_{\Delta}$ can only be determined with roughly 50 % accuracy.

Nevertheless, the fact that Eq. (3) nearly describes the measured gap dispersion indicates that the high-momentum contribution to spatial fluctuations of $\Delta(\mathbf{r})$ is absent. This case could physically correspond to a situation, where the 2D resistivity would spatially fluctuate due to the variations of interface resistance between neighboring crystallites (that would lead to small $q<q_{D}$ momenta contribution), while the 3D resistivity within each crystallite would remain identical (no change in the contribution of the $q>q_{D}$ momenta).

IV Discussion

We have shown that local probe spectroscopic studies are able to demonstrate a direct relation between the properties of the electronic excitation spectrum pertinent to the lowest-energy scale: the superconducting gap ($\sim 1$ meV) on one hand, and to much higher energy scales (few dozens of meV) relevant to the Coulomb anomaly induced by the Altshuler-Aronov effect on the other hand. This supports the relevance of the “fermionic” mechanism and Finkelstein theory for superconductivity suppression by moderate disorder in NbN thin films not only at the global scale but also locally.

We characterized the nanoscale low-temperature granularity that develops in the superconducting NbN film and have shown that electron-electron interactions play a crucial role in determining this precise granularity. A natural generalization then would be to study and understand the temperature evolution of the inhomogeneity of $\Delta(\textbf{r})$ and $\alpha(\textbf{r})$, in particular around and above $T_{c}$, to connect to the pseudogap phenomenon and its current understanding in various materials. Indeed, we point out that for the present ultrathin film the ratio $\Delta(0)/k_{B}T_{c}\approx 3.7$ is twice larger than the BCS theory prediction 1.76, and significantly larger than the value 2.1 known for thicker and less disordered NbN films Semenov , fabricated using the same set-up, procedures and recipees than in the present work.

The increase of the $\Delta(0)/k_{B}T_{c}$ ratio with the sheet resistance can be understood within the usual mechanism of strong thermal phase fluctuations EmeryKivelson , which are known to be important in 2D superconducting systems with small superfluid density (suppressed due to strong disorder). Indeed, the width of the Gaussian thermal fluctuational region in 2D, as seen from conductivity measurements, can be estimated by comparing the Aslamazov-Larkin paraconductivity $\sigma_{\text{AL}}^{\text{2D}}=(e^{2}/16\hbar)\ln(T/T_{c0})$ and normal-state conductivity as $\mathrm{Gi}=(e^{2}/16\hbar)R_{\Box}$, where $T_{c0}$ is the mean-field transition temperature. Taking into account that phase fluctuations are about three times stronger than GiLV-book , we obtain an estimate for the renormalized $T_{c}$, which coincides with the Berezinskii-Kosterlitz-Thouless transition: $(T_{c0}-T_{c})/T_{c0}\approx 3\,\text{Gi}$. With the normal-state resistance $R_{\Box}\approx 7$ KOhm (see Fig. 1a), this gives $T_{c0}\approx 1.4T_{c}$. Thus we conclude that the major part of the increase of the $\Delta(0)/k_{B}T_{c}$ ratio can be understood as being due to fluctuation suppression of the measured $T_{c}$ with respect to the mean-field value $T_{c0}$, which should be compared to the gap width $\Delta(0)$. The same conclusion was reached in Ref. Chand2012, , where much thicker (actually 3D) NbN films were studied.

The above explanation of the enhanced $\Delta(0)/k_{B}T_{c}$ and pseudogap phenomenon (discussed in more details in the previous paper Carbillet ) refers to the “bosonic” mechanism of phase fluctuations, while local correlation between Altshuler-Aronov exponent $\alpha(\textbf{r})$ and the low-temperature gap value $\Delta(\textbf{r})$ are nicely explained within the “fermionic” mechanism. In that sense NbN films seem to demonstrate simultaneous action of different superconductivity suppression mechanisms. It would be interesting to see if strongly disordered NbN films exhibit a double-gap behavior found recently in InOx filmsDubouchet2019 . Other methods such as the ones addressing local coherence through scanning Josephson microscopy Bergeal2008 by mapping the local critical current might also shed important new light Cho2019 , when compared to other local spectroscopic properties such as the local gap map or peak height map. Finally, addressing the LDOS measurements and its mesoscopic fluctuations near $T_{c}$ in disordered thin films similar to ours would also help testing and comparing with recent theoretical predictionsSF2005 ; we2012 ; Burmistrov .

V Conclusion

In conclusion, we have demonstrated sharp spatial anti-correlations in highly resistive NbN thin films between the local magnitude of the superconducting energy gap $\Delta(\textbf{r})$ and the value of the exponent $\alpha(\textbf{r})$ characterizing the local density-of-states suppression at the Fermi level due to the Altshuler-Aronov effect induced by Coulomb repulsion. The correlation lengths characteristic of the emergent inhomogeneities at $T\ll T_{c}$ are $l_{\Delta}\approx 27$ nm for $\Delta(\textbf{r})$ and $l_{\alpha}\approx 18$ nm for $\alpha(\textbf{r})$, while the cross-correlation length $l_{\text{cross}}$ between $\Delta(\textbf{r})$ and $\alpha(\textbf{r})$ is $l_{\text{cross}}\approx 25$ nm. Our results are in agreement with the predictions of a local version of the Finkelstein theory for the “fermionic” mechanism of superconductivity suppression by disorder. This behavior is generally expected for any 2D superconducting materials satisfaying: (a) a strong Coulomb interaction between conduction electrons, (b) a spatial scale $l_{\alpha}$ characteristic of the inhomogeneities of the Altshuler-Aronov density-of-states suppression much larger than the superconducting coherence length $\xi(0)$. Our results suggest that in NbN films the emerging granularity can be fully described by a “fermionic” mechanism. Additional theoretical and experimental work is needed to see if the same claim is valid for other recently discovered 2D disordered materials, including monolayer or bilayer dichalcogenides Costanzo ; Xi , bilayer graphene Cao or interface oxydes Caprara .

Acknowledgements.

We thank D. S. Antonenko for useful discussions. This work was supported by the French ANR project SUPERSTRIPES under the reference number ANR-15-CE30-0026, by the Russian Scientific Foundation under Grant No. 20-12-00361 and by the Russian Academy of Sciences Program “Low temperature physics”.

Appendix A Numerical computations of the auto-correlation and cross-correlation functions

The auto-correlation functions calculated from the gap map presented in Fig. 2b and the exponent map in Fig. 2c are computed using the 2D expression (4) with periodic boundary conditions to treat the finite size spectroscopic data. We chose such periodic boundary conditions for simplicity. Nevertheless, these latter conditions do not affect the global behavior of the auto-correlation function since the correlation lengths that are extracted from the dependence of the auto-correlation functions are about $1/10-1/20$ of the image size. Thus the spatial dependence of the auto-correlation function is not affected by boundary conditions in our case.

In order to ensure proper normalization, the auto-correlation functions are calculated by subtracting the mean value of the variable of interest and by dividing by the standard deviation $\sigma$ of this variable. For instance for the local gap values $\Delta(\mathbf{r})$ the 2D auto-correlation function $\rho_{\Delta}(\mathbf{r})$ is defined by

[TABLE]

where $N$ is the total number of pixels in the $\Delta(\mathbf{r})$ gap map. The 2D auto-correlation functions were calculated in this way for both $\Delta(\mathbf{r})$ and $\alpha(\mathbf{r})$ data presented in Fig. 2. Since $\rho_{\Delta}(\mathbf{r})$ and $\rho_{\alpha}(\mathbf{r})$ were rather isotropic we decided to circularly average them to extract the correlation lengths $l_{\Delta}$ and $l_{\alpha}$. The result of the circular averaging of $\rho_{\Delta}(\mathbf{r})$ and $\rho_{\alpha}(\mathbf{r})$ is presented in Fig. 3, where the spatial dependence of $\rho_{\Delta}(r)$ and $\rho_{\alpha}(r)$ is plotted as a function of the distance $r$. The correlation lengths are determined in such a way that the correlation functions remain above the noise level close to zero.

The 2D cross-correlation function $\rho_{\text{cross}}(\mathbf{r})$ between $\Delta(\mathbf{r})$ and $\alpha(\mathbf{r})$ is calculated in the same way as above using also periodic boundary conditions with the following formula:

[TABLE]

$\rho_{\text{cross}}(\mathbf{r})$ is plotted in panel g of Fig. 2. The 2D cross-correlation function between $\Delta(\mathbf{r})$ and the topography $Z(\mathbf{r})$ presented in panel a of Fig. 2 is calculated in the same way, $\alpha(\mathbf{r})$ in (5) being replaced by $Z(\mathbf{r})$. Finally since $\rho_{\text{cross}}(\mathbf{r})$ has a rather isotropic behavior, as seen on Fig. 2g, it was circularly averaged in order to extract the cross-correlation length $l_{\text{cross}}$ from the spatial dependence of $\rho_{\text{cross}}(r)$.

Appendix B Altshuler-Aronov anomaly in inhomogeneous films

We consider the tunneling density of states $\nu(E)$ in a diffusive normal film, where the local resistivity fluctuates with the correlation length $a$. The mean free path $l$, the film thickness $d$ and other parameters satisfy the experimentally relevant set of inequalities: $l\ll d\ll\sqrt{D/E}\ll a$. Here and below we use the system of units with $k_{B}=\hbar=1$.

B.1 Homogeneous film

The theory of the Altshuler-Aronov anomaly, also called zero-bias anomaly, Altshuler ; Altshuler1985 ; Ingold ; Fin1983 ; LevitovShytov ; KamenevAndreev predicts that the Coulomb-induced suppression of the tunneling density of states can be written in the form $\nu(E)=\nu_{0}e^{-S(E)}$, where $\nu_{0}$ is the bare density of states at the Fermi level in the absence of interaction. The function $S(E)$ can be considered as the gauge phase correlation function Fin1983 ; KamenevAndreev or the instanton action LevitovShytov . In the case of diffusive electron motion it can be written in the form

[TABLE]

where $R_{Q}=h/e^{2}$ is the resistance quantum and $R(\omega)$ is the spreading resistance between the diffusive length $r_{\text{in}}(\omega)=\sqrt{D/\omega}$ and electromagnetic field propagation length $r_{\text{out}}(\omega)=\sqrt{D\omega_{0}/\omega^{2}}$, where $\omega_{0}$ is the plasma frequency. In the 2D geometry with a constant sheet resistance $R_{\Box}$,

[TABLE]

Taking the integral in Eq. (6) we arrive at the standard result for homogeneous films:

[TABLE]

For practical purposes, the slow log-normal energy dependence of $\nu(E)$ predicted by Eq. (8) can be replaced by the power law $\nu(E)\propto E^{\alpha}$, with an almost energy-independent index $\alpha(E)=-\partial S_{\text{homo}}(E)/\partial\ln E$, leading to Eq. (1).

B.2 Inhomogeneous film

Consider now the Altshuler-Aronov anomaly in an inhomogeneous film. We will treat inhomogeneity in a simple model, where the local sheet resistance can be written as $\rho(\mathbf{r})=R_{\Box}+\delta\rho(\mathbf{r})$, with $\delta\rho(\mathbf{r})$ being a Gaussian noise with zero average and the correlation length $a$. The spreading resistance $R(\mathbf{r}_{0},\omega)$ now depends on the tunneling position $\mathbf{r}_{0}$ and can be estimated using logarithmic dependence of resistance in 2D. Therefore, with logarithmic accuracy one may calculate $R(\mathbf{r}_{0},\omega)$ replacing $\rho(\mathbf{r})$ by $\rho(\mathbf{r}_{0})$ for $|\mathbf{r}-\mathbf{r}_{0}|<a$ and by $R_{\Box}$ for $|\mathbf{r}-\mathbf{r}_{0}|>a$.

In the experimentally relevant case, $E\gg D/a^{2}$, and therefore $r_{\text{in}}(\omega)\ll a$. Then the spreading resistance $R(\mathbf{r},\omega)$ can be estimated as

[TABLE]

Substituting that into Eq. (6), we get

[TABLE]

where $S_{\text{homo}}(E)$ is given by Eq. (8). Taking the derivative, we see that the AA exponent fluctuates according to Eq. (2), where we identify $a$ with $l_{\alpha}$.

Appendix C Gap fluctuations in inhomogeneous films

If the correlation length of the local resistivity fluctuations exceeds the superconducting coherence length, $a\gg\xi(0)$, then the spectral gap at a given point $\mathbf{r}$ is determined by the order parameter $\Delta$ in a homogeneous system with the resistivity $\rho=\rho(\mathbf{r})$. The dependence of $\Delta$ on $\rho$ results from the Coulomb repulsion enhanced in a diffusive media Altshuler1985 . Assuming the effect can be treated perturbatively, we will rederive the known expression Ovchina ; Maekawa ; Takagi ; Fukuyama focusing on proper cutoffs of energy and momentum integrals.

C.1 Order parameter suppression in homogeneous films

Coulomb suppression of superconductivity Fin1983 can be described in terms of the energy-dependent correction to the Cooper channel interaction constant we2012 :

[TABLE]

where $\nu$ is the 3D density of states at the Fermi energy per one spin projection, both Matsubara energies, $\varepsilon$ and $\varepsilon^{\prime}$, are assumed to be larger than $T_{c}$, and

[TABLE]

Perturbative account of $\delta\lambda_{\varepsilon,\varepsilon^{\prime}}$ in the normal state reproduces the log-cube shift of $T_{c}$ Ovchina ; Maekawa ; Takagi ; Fukuyama . In the superconducting state, it leads to the shift of the order parameter compared to the clean case. In the low-temperature limit ($T\ll T_{c}$), this shift is given by we2012

[TABLE]

where $\Delta_{0}$ is the gap value in the clean system. This formula generalizes the BCS relation $\delta\Delta/\Delta=\delta\lambda/\lambda^{2}$ to the case of an energy-dependent coupling constant.

Since the energies $\varepsilon$ and $\varepsilon^{\prime}$ in Eq. (13) are limited by $\omega_{D}$, one can represent the general expression (11) as

[TABLE]

where $\delta\lambda_{\varepsilon,\varepsilon^{\prime}}^{\text{low}}$ is the energy-dependent contribution from low ($q<q_{D}$) momenta and $\delta\lambda^{\text{high}}$ is the energy-independent contribution from high ($q>q_{D}$) momenta:

[TABLE]

and $q_{D}=\sqrt{\omega_{D}/D}$ is determined by the Debye frequency.

The high-momentum contribution $\delta\lambda^{\text{high}}$ is just a correction to the Cooper channel constant. According to Eq. (13) it leads to the perturbative suppression of the order parameter: $\delta\Delta^{\text{high}}/\Delta=\delta\lambda^{\text{high}}/\lambda^{2}$. This expression can be extended beyond the perturbation theory using the BCS relation that results in the renormalization of the bare gap: $\Delta_{0}\to\Delta_{1}$, where

[TABLE]

Once the high-momentum contribution is incorporated into $\Delta_{1}$, the perturbative shift of the gap due to the low-momentum part $\delta\lambda^{\text{low}}_{\varepsilon,\varepsilon^{\prime}}$ should be determined by Eq. (13) with $\Delta_{0}$ replaced by $\Delta_{1}$.

In order to estimate the low-momentum contribution $\delta\lambda^{\text{low}}_{\varepsilon,\varepsilon^{\prime}}$ for our NbN films, we have to compare $\omega_{D}\approx 300$ K and the transverse Thouless energy $\hbar D(\pi/d)^{2}\approx 950$ K. Since the latter is larger, one can leave only the $n=0$ mode in Eq. (12), leading to the 2D contribution (the 2D resistivity is $\rho=R_{Q}/4\pi\nu Dd$)

[TABLE]

and hence to the following shift of $\Delta$:

[TABLE]

For our NbN films with $\rho=R_{\Box}^{\text{max}}=6.8$ k$\Omega$, $\hbar\omega_{D}\approx 300$ K and $\Delta=12$ K, the low-momentum contribution to the gap suppression according to Eq. (19) is given by $\delta\Delta^{\text{low}}/\Delta_{1}\approx 0.55$, which is at the border of applicability of the perturbation theory. At the same time, $\lambda^{\text{high}}$ in Eq. (15) is determined by the 3D diffusion which should be cut off at $q\sim 1/l$: $\delta\lambda^{\text{high}}\sim-(\rho/\pi^{2}R_{Q})(d/l)\sim 0.1$, yet with an unknown prefactor. Depending on the latter, $\Delta_{1}/\Delta_{0}$ is in the range of several tens of percent, that is a bit too small to describe experimental data.

C.2 Inhomogeneous films

In Sec. C.1, we calculated the gap suppression $\delta\Delta$ with respect to the clean case ($\rho\to 0$). To get gap fluctuations $\delta\Delta(\mathbf{r})$ in an inhomogeneous system with long-range fluctuations of $\rho(\mathbf{r})$, we simply replace $\rho$ by $\delta\rho(\mathbf{r})$:

[TABLE]

where $\delta\lambda^{\text{high}}(\mathbf{r})$ contains the large-momenta contribution:

[TABLE]

Equation (20) is derived in the model, where the resistivity [both the 3D $\rho_{3}(x,y)$ and the 2D $\rho(x,y)$ related via $\rho_{3}=\rho/d$] exhibits small spatial fluctuations. In a more realistic model presumably applicable to our NbN films, these are the crystalline interfaces that are the sources of resistance fluctuations at scales larger that a few nm. For such a model of structural disorder, high-momentum fluctuations of $\delta^{\text{high}}(\mathbf{r})$ are absent, and one arrives at Eq. (3).