Local density of states and its mesoscopic fluctuations near the transition to a superconducting state in disordered systems

Local density of states and its mesoscopic fluctuations near the transition to a superconducting state in disordered systems


We develop a theory of the local density of states (LDOS) of disordered superconductors, employing the non-linear sigma-model formalism and the renormalization-group framework. The theory takes into account the interplay of disorder and interaction couplings in all channels, treating the systems with short-range and Coulomb interactions on equal footing. We explore 2D systems that would be Anderson insulators in the absence of interaction and 2D or 3D systems that undergo Anderson transition in the absence of interaction. We evaluate both the average tunneling density of states and its mesoscopic fluctuations which are related to the LDOS multifractality in normal disordered systems. The obtained average LDOS shows a pronounced depletion around the Fermi energy, both in the metallic phase (i.e., above the superconducting critical temperature ) and in the insulating phase near the superconductor-insulator transition (SIT). The fluctuations of the LDOS are found to be particularly strong for the case of short-range interactions – especially, in the regime when is enhanced by Anderson localization. On the other hand, the long-range Coulomb repulsion reduces the mesoscopic LDOS fluctuations. However, also in a model with Coulomb interaction, the fluctuations become strong when the systems approaches the SIT.

72.15.Rn ,   71.30.+h ,   73.43.Nq  

I Introduction

Disordered superconductors show remarkable physics governed by interplay of superconductivity and Anderson localization. In particular, in two-dimensional (2D) systems, the competition between these two phenomena leads to a direct quantum phase transition between the insulating and superconducting states—the superconductor-insulator transition (SIT) [Goldman and Markovi, 1998,Gantmakher and Dolgopolov, 2010]. This is a zero-temperature transition that may be driven by varying the normal-state resistivity of a disordered film; experimentally, this is usually achieved by changing the film thickness. At a finite (but sufficiently low) temperature the insulating and superconducting phases of the film are separated by a metallic state. The physics of SIT and, more generally, of insulating, superconducting, and metallic states around SIT, has attracted a great deal of attention.

On the experimental side, two complementary approaches have been widely used to characterize the physics of disordered superconducting films under variation of temperature, film thickness, and magnetic field: (i) transport measurements and (ii) space-resolved tunneling spectroscopy. In the present paper, we focus on the second one and develop a theory of local density of states (LDOS) — including both its disorder-averaged value and fluctuations — as measured in space-resolved tunneling experiments.

Particularly intriguing experimental findings on tunneling spectroscopy of 2D disordered superconducting systems were provided by experiments on TiN and InO films [Sacépé et al., 2008, 2010, 2011; Sherman et al., 2014]. In short, it was found that (i) the pronounced soft gap in the tunneling spectrum survives across the superconductor-metal transition (i.e., with increasing temperature above ) and across SIT, and (ii) there are strong point-to-point fluctuations of the shape of the energy-dependence of LDOS on the superconducting side of the transition (i.e. below ). These results have been interpreted as evidence of (i) the existence of preformed Cooper pairs leading to a “pseudogap” in the non-superconducting states (metallic and insulating) [Sacépé et al., 2008, 2010, 2011; Sherman et al., 2014] and (ii) localization of some of Cooper pairs on the superconducting side of the transition, with the fraction of localized Cooper pairs increasing when the system approaches the SIT [Sacépé et al., 2011]. Qualitatively similar features, although considerably less pronounced, were observed in experiments on NbN films [Mondal et al., 2011,Noat et al., 2013]. Finally, a recent work on MoC films [Szab et al., 2016] did not discover any sizeable “pseudogap” or spatial fluctuations effects at all; the gap observed there was related to by the standard formula of the Bardeen-Cooper-Schrieffer (BCS) theory.

In order to understand the experimental findings — including features that are common for different materials as well as differences between the materials — one clearly needs the corresponding theory. In numerical works by Ghosal, Randeria, and Trivedi, Refs. [Ghosal et al., 1998,Ghosal et al., 2001], a solution of Bogoliubov-de Gennes equations for a 2D model with short-range interaction was carried out. It was found that the tunneling density of states shows a hard gap across the SIT and strong spatial fluctuations. More recently, these results were corroborated by Quantum Monte Carlo simulations [Bouadim et al., 2011]. While these results are very insightful, the numerical simulations for the inherently interacting problem are limited by relatively small system sizes. This makes it difficult to explore parametric dependences of observables in a sufficiently broad range, especially since the problem is characterized by a hierarchy of relevant length and energy scales. Such parametric dependences may be studied within analytical approaches, which are also expected to shed more light on underlying physical mechanisms. Feigel’man, Ioffe, Kravtsov, Yuzbashyan, and Cuevas, Refs. [Feigel’man et al., 2007,Feigel’man et al., 2010], studied the LDOS for a 3D system in the vicinity of Anderson-localization transition within a solution of the self-consistent BCS-type equation for the case of a short-range interaction. They found that a pseudogap develops when the superconductor is built out of localized single-particle states, and that this pseudogap increases when the system approaches the SIT.

In the present paper, we develop a theory of the LDOS of disordered superconductors which employs the non-linear sigma-model (NLSM) formalism and the renormalization-group (RG) framework, and goes beyond the analysis of Refs. [Feigel’man et al., 2007,Feigel’man et al., 2010] in several important aspects. First, our theory takes into account mutual influence of disorder and interaction couplings in all channels. This influence leads to the renormalization that becomes strong for systems with sufficiently strong disorder (in particular those that are not too far from SIT). Second, we consider 2D systems with short-range and Coulomb interaction on equal footing. Third, we use the same formalism to explore (i) 2D systems that would be Anderson insulators in the absence of interaction and (ii) 2D or 3D systems that undergo Anderson transition in the absence of interaction. Fourth, we evaluate both the average LDOS and its mesoscopic fluctuations (which are related to the LDOS multifractality in normal disordered systems). Fifth, when calculating the average LDOS and its moments, we take into account renormalization effects originating from all interaction channels.

On the technical side, we exploit our recent works in two complementary directions: Refs. [Burmistrov et al., 2012,Burmistrov et al., 2015a], where the phase diagram and transport characteristics of 2D disordered systems around SIT were studied by means of the NLSM renormalization group, on the one hand, and Refs. [Burmistrov et al., 2013,Burmistrov et al., 2015b,Burmistrov et al., 2015c], where the LDOS multifractality was studied near Anderson metal-insulator transition (MIT) in a normal (i.e., not superconducting) system with Coulomb interaction, on the other hand. Application of a unified approach to the LDOS and its fluctuations near MIT (Refs. [Burmistrov et al., 2013,Burmistrov et al., 2015b,Burmistrov et al., 2015c]) and in disordered superconductors (this work) turns out to be very helpful for understanding similarities and differences between the two cases. We will return to this issue in the end of the paper.

The outline of the paper is as follows. In Sec. II we introduce the NLSM formalism and construct operators corresponding to the moments of LDOS. The anomalous dimensions of the moments of LDOS found within the two-loop approximation are presented in Sec. III. The obtained two-loop results are used in Sec. IV to analyze the scaling behavior of the disorder-averaged LDOS and of the LDOS moments for the following three cases: (i) superconducting transition in 2D system with weak short-ranged interactions; (ii) superconducting transition in 2D system with Coulomb interaction; (iii) superconducting transition in a system with weak short-ranged interactions which, in the absence of interactions, is close to the Anderson transition. Our results and conclusions are summarized in Sec. V. Several appendices contain technical details on the one- and two-loop RG equations and their analysis.

Ii Formalism

ii.1 NLSM action

We start with the description of the NLSM formalism to be used for the calculation of the local density of states and its fluctuations near the transition to the superconducting state. The action of the NLSM is given as a sum of the non-interacting part, [Wegner, 1979,Efetov et al., 1980], and terms arising from the interactions in the particle-hole singlet and triplet, and Cooper channels [Finkel’stein, 1983,Finkel’stein, 1984] (see Refs. [Finkelstein, 1990,Belitz and Kirkpatrick, 1994] for review):



Here we use notations from Ref. [Burmistrov et al., 2015a]. The Drude conductivity (including spin) in units is denoted as . The quantities and are interaction parameters in the singlet particle-hole, triplet particle-hole, and singlet Cooper channels, respectively. The parameter introduced by Finkelstein [Finkel’stein, 1983] describes the renormalization of the frequency term in the action (1).

The action (1) involves the following matrices


where stand for replica indices and integer numbers correspond to the Matsubara fermionic energies and . The sixteen matrices,


act in the spin (subsrcipt ) and particle-hole (subscript ) spaces. The corresponding Pauli matrices are defined in a standard form as follows


The vector combines three matrices, . The matrix field obeys the following constraints:


The charge conjugation matrix satisfies the following relation: . The matrix (as well as the trace operator ) acts in the replica, Matsubara, spin, and particle-hole spaces.

In the case of Coulomb interaction, the parameters and are related to each other, . This relation holds in the course of the renormalization [Finkel’stein, 1983]. This relation also reflects the symmetry of the NLSM action (1) under the spatially independent rotations of the matrix (so-called -invariance) [Baranov et al., 1999a,Burmistrov et al., 2015a].

ii.2 Moments of the local density of states

As usual, the local density of states is expressed via the exact single-particle Green function. Within the NLSM formalism, the disorder-averaged LDOS is determined by the operator which is linear in :


Here symbol denotes the trace in spin and particle-hole spaces only. The index denotes a fixed replica and denotes the averaging with the NLSM action (1). The density of states at energy of the order of inverse elastic scattering time is denoted by . We remind that plays a role of the high-energy (ultraviolet) cutoff of the NLSM theory. The disorder-averaged LDOS can be obtained after the analytic continuation of to the real energies, .

Next, let us introduce the irreducible two-point correlation function


which allows us to find the second moment of the LDOS. In the NLSM approach, the correlator at coinciding spatial points is related to the following bilinear in operator:


The correlation function defined for real energies can be obtain from the following Matsubara counterpart


after analytic continuation: , , and . We use the convention that and .

The following comments are in order here.

  • The condition that replica indices and are nonequal in Eq. (9) stems from the fact that the two-point correlation function measures mesoscopic fluctuations of the LDOS. This forbids interaction lines between two fermionic loops corresponding to the LDOS in the diagrammatic approach.

  • The bilinear-in- operator (9) is the eigenoperator under the action of the RG (below we shall explicitly prove this statement by means of the two-loop calculations).

  • Disorder-averaged higher moments of the LDOS can be expressed in terms of higher-order irreducible correlation functions of the -field similarly to Eq. (9). Explicit examples of operators corresponding to the third and forth moments of the LDOS can be found in Ref. [Burmistrov et al., 2015c]. The corresponding operators are also eigenoperators of the renormalization group.

Iii Renormalization group for LDOS

In this section we outline the RG formalism in the context of the calculation of the moments of the distribution of the local DOS as derived from the NLSM.

iii.1 Perturbative expansion

To resolve the nonlinear constraint (5) we adopt the square-root parametrization:


We use the following notations: and with and . As a consequence of the charge-conjugation constraint (5), the blocks and obey the following relations:


These relations imply that elements in the expansion are real or purely imaginary. The perturbative (in ) analysis of the NLSM action (1) is performed by expanding the action in powers of .

From the expansion the NLSM action (1) to the second order in , we find the bare propagators (see Ref. [Belitz and Kirkpatrick, 1994]). The propagators in the particle-hole channel (diffusons) read ( and )


where and . The standard diffusive propagator is given by


with . The diffusons renormalized by interaction in the particle-hole channels read


The propagators in the particle-particle channel (cooperons) can be written as ( and )


where and . The propagator stands for the standard superconducting-fluctuation propagator:


Here we have introduced . The quantity is the diffusion coefficient, is the ultraviolet energy scale, and denotes the digamma function. We note that the fluctuation propagator (16) is written under the assumption that the infrared energy scale of the theory is determined by temperature.

For the purpose of regularization in the infrared, we add the extra term


to the NLSM action (1). The presence of this term in the action results, in particular, in the substitution of for in the propagators (13), (14), and (16).

iii.2 The disorder-averaged LDOS

For the sake of completeness, and in order to set notations, we remind the reader the result of one-loop renormalization of the disorder-averaged LDOS. Expanding the matrix to the second order in , we derive from Eq. (6) the following expression:


Next, using Eqs. (12) and (15), we find


Finally, performing the analytic continuation to real frequencies, , we obtain that the disorder-averaged LDOS can be written as


where the renormalization factor is given by (see Refs. [Al’tshuler and Aronov, 1985,Kamenev and Levchenko, 2009] for a review):


Here , , , and are retarded propagators obtained from the corresponding Matsubara propagators. The Keldysh part of the fluctuation propagator is related to the retarded one via the bosonic distribution function as follows: . The fermionic distribution function is denoted as . We adopt the following short-hand notation:


We note that the definition (20) of coincides with the definition of the field renormalization constant in Ref. [Belitz and Kirkpatrick, 1994] and the wavefunction renormalization constant in Ref. [Castellani et al., 1984a]. We stress that differs from the Finkel’stein’s frequency renormalization factor and should not be confused with the latter.

To derive the RG equation for , we set temperature and energy to zero and study the dependence of on the infrared regulator . Then in dimensions we obtain


Here are dimensionless interaction amplitudes and denotes dimensionless resistivity, with and being the area of a -dimensional sphere. As usual, Eq. (23) determines the anomalous dimension of the disorder-averaged LDOS. Using the minimal subtraction scheme (see e.g., Ref. [Amit, 1993]), we obtain in the one-loop approximation [Al’tshuler and Aronov, 1979a,Finkel’stein, 1984,Castellani et al., 1984a]:


where is the running RG length scale. We note that a more accurate treatment of the term with the fluctuation propagator in Eq. (21) (see Appendix A of Ref. [Burmistrov et al., 2015a]) results in exactly the same RG equation as Eq. (24). Thus, this one-loop RG equation is formally exact in all three interaction couplings , , and . In the case of fully broken spin-rotational symmetry, Eq. (24) holds with the contribution of the triplet particle-hole channel in the right hand side being omitted.

iii.3 The second moment of the LDOS

Now we consider the renormalization of the second moment of the LDOS. We restrict our consideration by one- and two-loop orders in .

One-loop results

In the one-loop approximation one obtains




Hence, we find


where . Setting and using as the infrared regulator, we get


Two-loop results

Details of the calculation of the two-loop contribution to the correlation function are presented in Appendix A. Using Eqs. (121) and (128) of Appendix A, we find the following two-loop contribution to the irreducible two-point correlation function:


Here the function denotes the polylogarithm,




Anomalous dimension of the second moment of LDOS

Above we have derived the dependence of on the momentum scale within the two-loop approximation. However, itself acquires renormalization [Baranov et al., 1999b]. The renormalized momentum scale is defined as follows


where denotes the renormalized conductivity at the momentum scale . As a consequence of Eq. (32), satisfies the following relation: . Using Eq. (23) and the one-loop result for the conductivity [Al’tshuler and Aronov, 1979b,Al’tshuler et al., 1980,Finkel’stein, 1983,Castellani et al., 1984a]


we find


We note that within the one-loop approximation, there is no contribution to Eq. (34) due to interaction in the Cooper channel.

By using Eqs. (28), (29), and (34), we can rewrite the second moment of the LDOS in terms of the renormalized momentum scale and the renormalization factor as follows:




Here we omit the terms that are finite in the limit . The bare value of is unity, , and


Next, we introduce a dimensionless quantity . With the help of Eqs (33) and (36), we express , and as follows:


The interaction parameters are renormalized at the one-loop level [Finkel’stein, 1983]. However this does not affect the two-loop result for the anomalous dimension of since is independent of . To the lowest orders in , the renormalization parameters become




The anomalous dimension of is derived from standard conditions that (as well as and ) does not depend on the momentum scale . In this way, we obtain the following two-loop result for the anomalous dimension of :


Here we omit ‘prime’ and ‘bar’ signs for brevity. We remind the reader that the function is defined in Eq. (30). We note that as it is known for free electrons [Wegner, 1979]. In the case of Coulomb interaction, one has . The interaction affects the anomalous dimension at the two-loop order only.

We emphasize that coefficients , and satisfy the relation


This guaranties the absence in Eq. (41) of terms divergent in the limit , i.e. the renormalizability of . In addition, this relation proves that the operator corresponding to is the RG eigenoperator. Indeed, if the operator corresponding to is the linear combination of several eigenoperators, the relation (42) would imply non-linear system of equations which has no non-trivial solutions in general.

Combining the above results, the second moment of the LDOS can be written as


The scaling behavior of and are governed by the anomalous dimensions  (24) and (41), respectively.

iii.4 The -th moment of the LDOS

In this section, we generalize the results obtained in the previous section for the -th moment of the LDOS. The important observation is that the irreducible -th moment of the LDOS, , involves connected contributions from averages of the number of matrices . Therefore, has no one- and two-loop contributions for . Consequently, as in the case of noninteracting electrons [Hf and Wegner, 1986; Wegner, 1987a, b], the anomalous dimension for the -th moment of the LDOS becomes proportional to the factor within one- and two-loop approximation (see details in Ref. [Burmistrov et al., 2015c]). Thus, we find


where the behavior of is determined by the following RG equation:


Here the function is defined in Eq. (30). We note that for the special cases and , one can demonstrate that Eq. (45) holds as well (see Appendix C of Ref. [Burmistrov et al., 2015c]).

As we have already mentioned above, the two-loop contributions to can be interpreted as the renormalization of diffusions and cooperons involved in the one-loop term (see Appendix B). Therefore, within the two-loop approximation, the corrections due to fluctuating Cooper pairs to comes from the term only. In the one-loop approximation the fluctuation corrections to the -th moment of the LDOS are fully determined by those in the average LDOS via the factor .

In the absence of spin-rotational symmetry, the anomalous dimension can be obtained from Eq. (45) as follows: (i) one omits the contribution of the triplet particle-hole channel and (ii) one multiples the right hand side of Eq. (45) by the factor (see Ref. [Burmistrov et al., 2015c]). Thus, in the case of broken spin rotational symmetry, Eq. (45) takes the following form:


Iv Scaling analysis

iv.1 Weak coupling RG equations in 2D

Recently [Burmistrov et al., 2015a], the full set of one-loop RG equations describing the renormalization of resistivity and interactions has been derived by means of the background field renormalization of the NLSM (1):


Here stands for the number of triplet diffusive modes. In the case of preserved spin-rotational symmetry, all three triplet diffusons contribute to the RG equations, . If spin-rotational symmetry is broken, all triplet modes are suppressed at long lengthscales, . In addition, in this case RG equation for should be ignored.

The above RG equations are derived in the lowest order in . Extending the previous result [Finkel’stein, 1984], in this order they contains all contributions due to the interaction in the Cooper channel, . Comparison of the disorder-independent and disorder-induced terms in the right hand side of Eq. (50) demonstrates that one-loop RG equations can be used towards the superconducting instability upto the length scale at which [Burmistrov et al., 2015a]. We note that similar conclusion follows from comparison of one and two-loop contributions in Eqs. (45) and (46).

iv.2 Weak short-ranged interactions in 2D

We start our analysis of the fluctuations of LDOS with the case of weak short-ranged interactions, , in two dimensions. The phase diagram and transport properties for this case have been discussed in details in Refs. [Burmistrov et al., 2012,Burmistrov et al., 2015a]. The existence of a large region of the superconducting phase with transition temperature higher than standard expression of BCS theory has been predicted.

For the sake of convenience, we briefly remind the reader the main steps of analysis of Ref. [Burmistrov et al., 2012]. Let us focus on the case of preserved spin-rotational symmetry. For , the set of RG equations (47) - (50) can be simplified to


Equation (52) yields usual weak localization behavior:


Here denotes the bare value of resistance. In the absence of interactions, Eq. (54) suggests that a strong Anderson insulator emerges at the scale . In the case of not too weak interactions (or sufficiently weak disorder), , Eqs. (53) reduce to the standard RG equation for the BCS instability in the clean case. Then the superconducting transition occurs at the scale of the order of .

For the case of disorder which is strong compared to the interaction, , the renormalization proceeds in two steps. At the first step, we can neglect the term. This brings us to a linear system of equations. The corresponding matrix has two eigenvalues:


where is doubly degenerate. We emphasize that the eigenvalue coincides with the one-loop result for . This occurs since the interactions in the NLSM action (1) are described by the operators bilinear in . Thus, three couplings are tend to the eigenvector corresponding to the positive eigenvalue . This eigenvector parametrizes the so-called BCS line


It is this relation between the couplings that one obtains starting from a standard BCS Hamiltonian with the attraction only. Projecting Eqs. (53) onto the BCS line, we get


We note that the RG flow in directions perpendicular to the eigenvector does not affect the results in any essential way (see Ref. [Burmistrov et al., 2012] for details).

We are interested in the case when the dominant bare interaction is the Cooper attraction such that the initial value . Equation (57) describes two distinct scenarios. For , the resistance becomes of order unity when is still small. This means that, with further increase of the length scale, the system flows toward a strong Anderson insulator.

In the opposite case, , increases under the RG transformation due to the first term in the r.h.s. of Eq. (57). The attractive interaction overtakes and reaches unity at the scale . We note that


i.e. strong attraction is arisen in the region of good metal. In this situation, we expect that with further increase of the length scale, the RG flow develops a superconducting instability (due to standard BCS-type mechanism). The temperature of this superconducting transition can be estimated as


In the considered regime of sufficiently strong disorder, , the temperature given by Eq. (59) is much higher than the clean BCS value


When becomes smaller than , Eq. (59) crosses over into Eq. (60). Therefore, shows a non-monotonous dependence on the disorder strength and gets strongly enhanced in the intermediate range of resistivity, . For a given interaction strength , the superconducting critical temperature is the largest when the system approaches the superconductor-insulator transition. The latter takes place at , i.e. . It is worth noting that at the superconductor-insulator transition the critical temperature on the superconducting side is given by , where is the localization length defined such as .

The above analysis [Burmistrov et al., 2012] for was extended in Ref. [Burmistrov et al., 2015a] to the case of strong interaction couplings. The numerical integration of the full RG equations (47) - (50) provided the phase diagram of the SIT. In particular, it was found that, when the system is initially close to the BCS line, the enhancement of occurs for