Local density of states and its mesoscopic fluctuations near the transition to a superconducting state in disordered systems
Abstract
We develop a theory of the local density of states (LDOS) of disordered superconductors, employing the nonlinear sigmamodel formalism and the renormalizationgroup framework. The theory takes into account the interplay of disorder and interaction couplings in all channels, treating the systems with shortrange 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 superconductorinsulator transition (SIT). The fluctuations of the LDOS are found to be particularly strong for the case of shortrange interactions – especially, in the regime when is enhanced by Anderson localization. On the other hand, the longrange 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.
pacs:
72.15.Rn , 71.30.+h , 73.43.NqI Introduction
Disordered superconductors show remarkable physics governed by interplay of superconductivity and Anderson localization. In particular, in twodimensional (2D) systems, the competition between these two phenomena leads to a direct quantum phase transition between the insulating and superconducting states—the superconductorinsulator transition (SIT) [Goldman and Markovi, 1998,Gantmakher and Dolgopolov, 2010]. This is a zerotemperature transition that may be driven by varying the normalstate 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) spaceresolved 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 disorderaveraged value and fluctuations — as measured in spaceresolved 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 superconductormetal transition (i.e., with increasing temperature above ) and across SIT, and (ii) there are strong pointtopoint fluctuations of the shape of the energydependence 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 nonsuperconducting 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 BardeenCooperSchrieffer (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 Bogoliubovde Gennes equations for a 2D model with shortrange 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 Andersonlocalization transition within a solution of the selfconsistent BCStype equation for the case of a shortrange interaction. They found that a pseudogap develops when the superconductor is built out of localized singleparticle 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 nonlinear sigmamodel (NLSM) formalism and the renormalizationgroup (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 shortrange 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 metalinsulator 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 twoloop approximation are presented in Sec. III. The obtained twoloop results are used in Sec. IV to analyze the scaling behavior of the disorderaveraged LDOS and of the LDOS moments for the following three cases: (i) superconducting transition in 2D system with weak shortranged interactions; (ii) superconducting transition in 2D system with Coulomb interaction; (iii) superconducting transition in a system with weak shortranged 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 twoloop 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 noninteracting part, [Wegner, 1979,Efetov et al., 1980], and terms arising from the interactions in the particlehole singlet and triplet, and Cooper channels [Finkel’stein, 1983,Finkel’stein, 1984] (see Refs. [Finkelstein, 1990,Belitz and Kirkpatrick, 1994] for review):
(1) 
where
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 particlehole, triplet particlehole, 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
(2) 
where stand for replica indices and integer numbers correspond to the Matsubara fermionic energies and . The sixteen matrices,
(3) 
act in the spin (subsrcipt ) and particlehole (subscript ) spaces. The corresponding Pauli matrices are defined in a standard form as follows
(4) 
The vector combines three matrices, . The matrix field obeys the following constraints:
(5) 
The charge conjugation matrix satisfies the following relation: . The matrix (as well as the trace operator ) acts in the replica, Matsubara, spin, and particlehole 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 (socalled 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 singleparticle Green function. Within the NLSM formalism, the disorderaveraged LDOS is determined by the operator which is linear in :
(6) 
Here symbol denotes the trace in spin and particlehole 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 highenergy (ultraviolet) cutoff of the NLSM theory. The disorderaveraged LDOS can be obtained after the analytic continuation of to the real energies, .
Next, let us introduce the irreducible twopoint correlation function
(7) 
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:
(8) 
The correlation function defined for real energies can be obtain from the following Matsubara counterpart
(9) 
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 twopoint 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 bilinearin operator (9) is the eigenoperator under the action of the RG (below we shall explicitly prove this statement by means of the twoloop calculations).

Disorderaveraged higher moments of the LDOS can be expressed in terms of higherorder 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 squareroot parametrization:
(10) 
We use the following notations: and with and . As a consequence of the chargeconjugation constraint (5), the blocks and obey the following relations:
(11) 
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 particlehole channel (diffusons) read ( and )
(12) 
where and . The standard diffusive propagator is given by
(13) 
with . The diffusons renormalized by interaction in the particlehole channels read
(14) 
The propagators in the particleparticle channel (cooperons) can be written as ( and )
(15) 
where and . The propagator stands for the standard superconductingfluctuation propagator:
(16) 
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.
iii.2 The disorderaveraged LDOS
For the sake of completeness, and in order to set notations, we remind the reader the result of oneloop renormalization of the disorderaveraged LDOS. Expanding the matrix to the second order in , we derive from Eq. (6) the following expression:
(18) 
Next, using Eqs. (12) and (15), we find
(19) 
Finally, performing the analytic continuation to real frequencies, , we obtain that the disorderaveraged LDOS can be written as
(20) 
where the renormalization factor is given by (see Refs. [Alâtshuler and Aronov, 1985,Kamenev and Levchenko, 2009] for a review):
(21) 
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 shorthand notation:
(22) 
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
(23) 
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 disorderaveraged LDOS. Using the minimal subtraction scheme (see e.g., Ref. [Amit, 1993]), we obtain in the oneloop approximation [Al’tshuler and Aronov, 1979a,Finkel’stein, 1984,Castellani et al., 1984a]:
(24) 
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 oneloop RG equation is formally exact in all three interaction couplings , , and . In the case of fully broken spinrotational symmetry, Eq. (24) holds with the contribution of the triplet particlehole 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 twoloop orders in .
Oneloop results
In the oneloop approximation one obtains
(25) 
and
(26) 
Hence, we find
(27) 
where . Setting and using as the infrared regulator, we get
(28) 
Twoloop results
Details of the calculation of the twoloop contribution to the correlation function are presented in Appendix A. Using Eqs. (121) and (128) of Appendix A, we find the following twoloop contribution to the irreducible twopoint correlation function:
(29) 
Here the function denotes the polylogarithm,
(30) 
and
(31) 
Anomalous dimension of the second moment of LDOS
Above we have derived the dependence of on the momentum scale within the twoloop approximation. However, itself acquires renormalization [Baranov et al., 1999b]. The renormalized momentum scale is defined as follows
(32) 
where denotes the renormalized conductivity at the momentum scale . As a consequence of Eq. (32), satisfies the following relation: . Using Eq. (23) and the oneloop result for the conductivity [Al’tshuler and Aronov, 1979b,Al’tshuler et al., 1980,Finkel’stein, 1983,Castellani et al., 1984a]
(33) 
we find
(34) 
We note that within the oneloop 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:
(35) 
where
(36) 
Here we omit the terms that are finite in the limit . The bare value of is unity, , and
(37) 
Next, we introduce a dimensionless quantity . With the help of Eqs (33) and (36), we express , and as follows:
(38) 
The interaction parameters are renormalized at the oneloop level [Finkel’stein, 1983]. However this does not affect the twoloop result for the anomalous dimension of since is independent of . To the lowest orders in , the renormalization parameters become
(39) 
and
(40) 
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 twoloop result for the anomalous dimension of :
(41) 
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 twoloop order only.
We emphasize that coefficients , and satisfy the relation
(42) 
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 nonlinear system of equations which has no nontrivial solutions in general.
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 twoloop 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 twoloop approximation (see details in Ref. [Burmistrov et al., 2015c]). Thus, we find
(44) 
where the behavior of is determined by the following RG equation:
(45) 
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 twoloop contributions to can be interpreted as the renormalization of diffusions and cooperons involved in the oneloop term (see Appendix B). Therefore, within the twoloop approximation, the corrections due to fluctuating Cooper pairs to comes from the term only. In the oneloop 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 spinrotational symmetry, the anomalous dimension can be obtained from Eq. (45) as follows: (i) one omits the contribution of the triplet particlehole 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:
(46) 
Iv Scaling analysis
iv.1 Weak coupling RG equations in 2D
Recently [Burmistrov et al., 2015a], the full set of oneloop RG equations describing the renormalization of resistivity and interactions has been derived by means of the background field renormalization of the NLSM (1):
(47)  
(48)  
(49)  
(50)  
(51) 
Here stands for the number of triplet diffusive modes. In the case of preserved spinrotational symmetry, all three triplet diffusons contribute to the RG equations, . If spinrotational 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 disorderindependent and disorderinduced terms in the right hand side of Eq. (50) demonstrates that oneloop 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 twoloop contributions in Eqs. (45) and (46).
iv.2 Weak shortranged interactions in 2D
We start our analysis of the fluctuations of LDOS with the case of weak shortranged 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 spinrotational symmetry. For , the set of RG equations (47)  (50) can be simplified to
(52) 
(53) 
Equation (52) yields usual weak localization behavior:
(54) 
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:
(55) 
where is doubly degenerate. We emphasize that the eigenvalue coincides with the oneloop 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 socalled BCS line
(56) 
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
(57) 
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
(58) 
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 BCStype mechanism). The temperature of this superconducting transition can be estimated as
(59) 
In the considered regime of sufficiently strong disorder, , the temperature given by Eq. (59) is much higher than the clean BCS value
(60) 
When becomes smaller than , Eq. (59) crosses over into Eq. (60). Therefore, shows a nonmonotonous 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 superconductorinsulator transition. The latter takes place at , i.e. . It is worth noting that at the superconductorinsulator 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