Multifractally-enhanced superconductivity in thin films
MMultifractally-enhanced superconductivity in thin films
I. S. Burmistrov a,b, ∗ , I. V. Gornyi c,d,e , A. D. Mirlin c,d,f,a a L. D. Landau Institute for Theoretical Physics, Semenova 1a, 142432, Chernogolovka, Russia b Laboratory for Condensed Matter Physics, National Research University Higher School of Economics, 101000 Moscow, Russia c Institute for Quantum Materials and Technologies, Karlsruhe Institute of Technology, 76021 Karlsruhe, Germany d Institut f¨ur Theorie der Kondensierten Materie, 76128 Karlsruhe, Germany e Io ff e Institute, 194021 St. Petersburg, Russia f Petersburg Nuclear Physics Institute, 188300 St. Petersburg, Russia
Abstract
The multifractal superconducting state originates from the interplay of Anderson localizationand interaction e ff ects. In this article we overview the recent theory of the superconductivity en-hancement by multifractality and extend it to describe the spectral properties of superconductorson the scales of the order of the superconducting gap. Specifically, using the approach based onrenormalization group within the nonlinear sigma model, we develop the theory of a multifractalsuperconducting state in thin films. We derive a modified Usadel equation that incorporates theinterplay of disorder and interactions at energy scales larger than the spectral gap and study thee ff ect of such an interplay on the low-energy physics. We determine the spectral gap at zerotemperature which occurs to be proportional to the multifracally enhanced superconducting tran-sition temperature. The modified Usadel equation results in the disorder-averaged density ofstates that, near the spectral gap, resembles the one obtained in the model of a spatially randomsuperconducting order parameter. We reveal strong mesoscopic fluctuations of the local densityof states in the superconducting state. Such strong mesoscopic fluctuations imply that the inter-val of energies in which the superconducting gap establishes is parametrically large in systemswith multifractally-enhanced superconductivity. Keywords:
Anderson localization, multifractality, superconductivity
1. Introduction
Anderson localization [1] and superconductivity are two fundamental quantum phenomena.Semiclassicaly, i.e., without taking into account quantum interference e ff ects, superconductivityis not a ff ected by the electron scattering on non-magnetic disorder that is compatible with thesymmetry of the superconducting order parameter (so-called “Anderson theorem” [2–4]). How-ever, at the next level of sophistication – with quantum e ff ects included, the Anderson theoremmeets Anderson localization. At this level, superconductivity and Anderson localization wereconsidered as antagonists: strong localization [5–8] or interplay of weak disorder and Coulombinteraction [9–15] were predicted to destroy superconductivity. This view was supported by a ∗ Corresponding author. Fax: +7-495-702-9317 e-mail: [email protected]
Preprint submitted to Annals of Physics February 1, 2021 a r X i v : . [ c ond - m a t . m e s - h a ll ] J a n iscovery [16] and the further study of the superconductor-insulator transition (SIT) in thin films(see Refs. [17–19] for a review).The paradigm of suppression of superconductivity by Anderson localization has been chal-lenged in Refs. [20, 21], where the enhancement of the superconducting transition temperature, T c , due to the multifractal behavior of wave functions near the Anderson transition (e.g., in three-dimensional disordered systems) has been demonstrated for the situation when the long-rangedCoulomb repulsion between electrons is not e ff ective. Later, the multifractal enhancement of T c in thin superconducting films has been predicted by the present authors [22, 23]. Thesetheoretical predictions have been supported by numerical solution of the disordered attractiveHubbard model on a two-dimensional lattice [24]. Recently, an enhancement of the supercon-ducting transition temperature by disorder observed in monolayer niobium dichalcogenides hasbeen explained by the multifractality [25, 26].The hallmark of the multifractally-enhanced superconductivity in disordered systems isstrong mesoscopic fluctuations of the local order parameter and the local density of states [21].The point-to-point fluctuations of the local density of states have been observed in many exper-iments on tunneling spectroscopy on thin superconducting films [27–32]. Recently, these fluc-tuations have been retrieved from numerical solution of disordered attractive two-dimensionalHubbard model [33, 34].For temperatures above the critical temperature, T > T c , the point-to-point fluctuations of thelocal density of states observed in the tunneling experiments are in qualitative agreement withthe theory [35] developed by the present authors for the mesoscopic fluctuations of the local den-sity of states in the normal phase of disordered superconducting films. Our theory employs therenormalization-group (RG) framework within the non-linear sigma-model (NLSM) formalism.The advantage of such approach is the possibility to take into account mutual influence of dis-order and interactions in all channels (particle-hole as well as particle-particle) and to describesystems with short-range and Coulomb interaction on equal footing.The experimental data on the fluctuations of the local density of states below T c , as well asrecent numerical data on the disordered Hubbard model represent a new challenge for an ana-lytical theory. In the present paper, we develop a theory of the mesoscopic fluctuations of thelocal density of states in the multifractal superconducting state. For the sake of concreteness,we shall focus on the case of thin superconducting films. We assume the absence of long-range(Coulomb) interaction, which can be achieved in structures put on a substrate with a high dielec-tric constant, and consider the regime of an intermediate disorder strength [22, 23], in which thesuperconducting transition temperature is parametrically enhanced by multifractality comparedto the conventional Bardeen-Cooper-Schrie ff er (BCS) result. We calculate the average density ofstates and its mesoscopic fluctuations in the low-temperature limit in the presence of the interplaybetween disorder and interactions.We demonstrate that the gap in the disorder-averaged density of states at zero temperatureis proportional to the superconducting transition temperature, i.e., it is also enhanced by multi-fractality. In spite of the enhancement of the spectral gap, we find that the combined e ff ect ofdisorder and interaction results in the suppression of the coherence peaks. Their height is finiteand controlled by the ratio of a bare disorder and interaction. We also show that the mesoscopicfluctuations of the local density of states in the superconducting state are enhanced in the sameway as in the normal phase for temperatures above T c .On the technical side, in order to describe the superconducting state, we derive the Usadelequation modified by renormalization e ff ects due to interplay of disorder and interactions (inboth particle-hole and particle-particle channels) at short scales. This modified Usadel equation2an be reformulated as a standard Usadel equation, but with an energy dependent gap function,and yields the same estimate for the superconducting transition temperature as the one derivedby means of the renormalization group in the normal phase. We discuss the relation between theapproaches based on the renormalization of the sigma model and the self-consistent solution ofthe gap equation in the presence of multifractal correlations of matrix elements.The structure of the paper is as follows. In Sec. 2 we overview the formalism of theFinkel’stein NLSM as applied to superconducting systems [14, 15, 22, 23]. The mean-fielddescription of the superconducting state is briefly outlined in Sec. 3. The description of the su-perconducting state beyond the mean-field approximation is developed in Sec. 4. The developedtheory is used to estimate the superconducting transition temperature (Sec. 5) and the energydependence of the gap function (Sec. 6). The results for the local density of states and its meso-scopic fluctuations in the superconducting state are presented in Sec. 7 and Sec. 8, respectively.Finally, our results and conclusions are summarized in Sec. 9.
2. Formalism of Finkel’stein NLSM
The action of the Finkel’stein NLSM is given as the sum of the non-interacting NLSM, S σ , and contributions resulting from electron-electron interactions, S ( ρ )int (the particle-hole singletchannel), S ( σ )int (the particle-hole triplet channel), and S ( c )int (the particle-particle channel) (see Refs.[36–38] for a review): S = S σ + S ( ρ )int + S ( σ )int + S ( c )int , (1)where S σ = − g (cid:90) d r Tr( ∇ Q ) + Z ω (cid:90) d r Tr ˆ ε Q , (2a) S ( ρ )int = − π T Γ s (cid:88) α, n (cid:88) r = , (cid:90) d r Tr I α n t r Q Tr I α − n t r Q , (2b) S ( σ )int = − π T Γ t (cid:88) α, n (cid:88) r = , (cid:88) j = , , (cid:90) d r Tr I α n t r j Q Tr I α − n t r j Q , (2c) S ( c )int = − π T Γ c (cid:88) α, n (cid:88) r = , (cid:90) d r Tr t r L α n Q Tr t r L α n Q . (2d)Here, Q ( r ) is the matrix field in the replica, Matsubara, spin, and particle-hole spaces. The traceTr acts in the same spaces. The matrix field obeys the nonlinear constraint, as well as charge-conjugation symmetry relation: Q ( r ) = , Tr Q = , Q = Q † = − CQ T C , C = it . (3)The nonlinear constraint on the matrix field Q suggests the following parametrization: Q = T − Λ T , T † = T − , CT T = T − C , Λ αβ nm = sgn ε n δ ε n ,ε m δ αβ t , (4)where α, β = , . . . , N r stand for replica indices and integers n , m correspond to the Matsubarafermionic frequencies ε n = π T (2 n + ff usive fluctua-tions in the normal state are associated with smooth variations around the metallic saddle point Q = Λ . 3he action (1) contains three constant matrices:ˆ ε αβ nm = ε n δ ε n ,ε m δ αβ t , ( I γ k ) αβ nm = δ ε n − ε m ,ω k δ αβ δ αγ t , ( L γ k ) αβ nm = δ ε n + ε m ,ω k δ αβ δ αγ t . (5)The sixteen matrices, t r j = τ r ⊗ s j , r , j = , , , , (6)operate in the particle-hole (subscript r ) and spin (subscript j ) spaces. The matrices τ , τ , τ , τ and s , s , s , s are the standard sets of the Pauli matrices.The bare value of the coupling constant g coincides with the Drude value of the conductance(in units e / h and including spin). The parameter Z ω describes the frequency renormalizationupon the renormalization group flow [36]. The bare value of Z ω is equal to πν/
4, where ν denotesthe bare density of states at the Fermi level including a spin-degeneracy factor. The interactionamplitude Γ s ( Γ t ) encodes interaction in the singlet (triplet) particle-hole channel. The interac-tion in the Cooper channel is denoted as Γ c . Its negative magnitude, Γ c <
0, corresponds toan attraction in the particle-particle channel. In what follows, it is convenient to introduce thedimensionless interaction parameters, γ s , t , c = Γ s , t , c / Z ω . If Coulomb interaction is present, thefollowing relation holds: γ s = −
1. This condition remains intact under the action of the renor-malization group flow [23, 36]. In the case of the BCS model and, under an assumption of strongdisorder, ω D τ (cid:28) ω D is the Debye frequency), the bare values of interaction amplitudes arerelated as − γ s = γ t = γ c . In what follows, we shall refer to such a relation between interactionamplitudes as the BCS line . The symmetry breaking that leads to the appearance of the superconducting state changesthe saddle point of the NLSM [39–42]. We shall follow the approach of Ref. [42] which allowsus to take into account the e ff ects of disorder and residual quasiparticle interactions beyond themean-field level.Let us single out the term with n = S ( c )int , Eq. (2d). Then, upon theHubbard–Stratonovich decoupling by spatially dependent fields ∆ α r ( r ), r = ,
2, we find S ( c )int = Z ω π T γ c (cid:90) d r (cid:88) α (cid:88) r = , (cid:2) ∆ α r ( r ) (cid:3) + Z ω (cid:90) d r (cid:88) α (cid:88) r = , ∆ α r ( r ) Tr t r L α Q + ˜ S ( c )int , (7)where ˜ S ( c )int = − π T Γ c (cid:88) α, n (cid:44) (cid:88) r = , (cid:90) d r (cid:16) Tr t r L α n Q (cid:17) . (8)Since, by construction, the action depends quadratically on ∆ α r ( r ), we can solve the equation ofmotion instead of performing the functional integral. The equation of motion reads ∆ α r ( r ) = π T | γ c | Tr t r L α Q ( r ) , r = , . (9)We note that the charge-conjugation symmetry, Eq. (3), guarantees that ∆ α , ( r ) are real functions.4he presence of nonzero ∆ α r ( r ) changes the saddle-point equation for Q . Performing variationof the NLSM action (1) with S ( c )int replaced by ˜ S ( c )int , we derive the following saddle-point equation: − D ∇ ( Q ∇ Q ) + [ ˆ ε, Q ] + (cid:88) α (cid:88) r = , ∆ α r ( r )[ t r L α , Q ] − π T (cid:88) α, n (cid:88) r = , (cid:88) j = γ j [ I α − n t r j , Q ] Tr I α n t r j Q − π T γ c (cid:88) α, n (cid:44) (cid:88) r = , [ t r L α n , Q ] Tr t r L α n Q = , (10)Here, γ ≡ γ s and γ = γ = γ ≡ γ t , and we have introduced the di ff usion coe ffi cient D = g / (16 Z ω ). We mention that in the absence of interaction Eq. (10) is nothing but the Usadelequation, well known in the theory of dirty superconductors.Substitution of ∆ α r ( r ) from the self-consistency equation (9) into the Usadel equation (10)results in the nonlinear matrix equation for Q ( r ). In general, such equation has many spatiallydependent solutions for Q ( r ) which mimic fluctuations of a disorder potential in the originalmicroscopic formulation of the problem. To enumerate all these solutions and to perform sum-mation over them seem like a daunting task. In order to circumvent this di ffi culty we shall usethe smallness of the bare resistance, 1 / g (cid:28)
1. Then, we can sum over spatially dependent solu-tions for Q ( r ) by treating the fluctuations around the spatially independent solution of the Usadelequation by means of the renormalization group. We start from splitting ∆ α r into constant andspatially dependent part, ∆ α r ( r ) = ∆ α r + δ ∆ α r ( r ) , (cid:90) d r δ ∆ α r ( r ) = . (11)The spatially varying part, δ ∆ α r ( r ), contains information about the mesoscopic fluctuations ofthe superconducting order parameter. These mesoscopic fluctuations turn out to be large (seeAppendix A) [21, 33, 34, 45]. In order to take them into account, we perform a formally exactintegration over δ ∆ α r ( r ). As a result, the spatial fluctuations of the Hubbard-Stratonovich fieldcorresponding to the order parameter get fully encoded in additional correlations of the Q field.This procedure results in a modification of S ( c )int , S ( c )int = Z ω V π T γ c (cid:88) α (cid:88) r = , (cid:2) ∆ α r (cid:3) + Z ω V (cid:88) α (cid:88) r = , ∆ α r Tr t r L α Q + ˆ S ( c )int , (12)where V denotes the volume of a superconductor, Q = V (cid:90) d r Q ( r ) , (13)and ˆ S ( c )int = − π T Γ c (cid:88) α, n (cid:88) r = , (cid:90) d r (cid:104) Tr t r L α n (cid:16) Q − Q δ n , (cid:17)(cid:105) . (14)After integrating over δ ∆ α r ( r ), the saddle-point equation (10) gets modified, − D ∇ ( Q ∇ Q ) + [ ˆ ε, Q ] + [ ˆ ∆ , Q ] − π T (cid:88) α, n (cid:88) r = , (cid:88) j = γ j [ I α − n t r j , Q ] Tr I α n t r j Q − π T γ c (cid:88) α, n (cid:88) r = , [ t r L α n , Q ] Tr t r L α n ( Q − Q δ n , ) = , (15)5here ˆ ∆ = (cid:88) α (cid:88) r = , ∆ α r t r L α , ∆ α r = π T | γ c | Tr t r L α Q . (16)It is worth noting that the information about the disorder-induced spatial fluctuations of theHubbard-Stratonovich field ∆ α r , and, hence, about the spatial fluctuations of the order parameter,is not lost after the integration over δ ∆ α r ( r ) and remains encoded in Eq. (15) through the couplingof ˆ ∆ and Q and the fluctuations of field Q . Before considering the e ff ect of such fluctuations, inthe following section we analyse the mean-field solution of the modified Usadel equation.
3. Mean-field description of the superconducting state
Let us seek the solution of the Usadel equation (15) and the self-consistency equation (16) inthe following form Q αβ nm = (cid:16) cos θ ε n sgn ε n δ ε n ,ε m + t φ sin θ ε n δ ε n , − ε m (cid:17) δ αβ , t φ = cos φ t + sin φ t , ∆ α = ∆ cos φ, ∆ α = ∆ sin φ. (17)Here, we assume that the spectral angle θ ε n is an even function of ε n . We note that the aboveansatz, Q , satisfies the charge-conjugation condition. Substituting Eq. (17) into the Usadelequation (15), we obtain D ∇ θ ε n − | ε n | sin θ ε n + ∆ cos θ ε n = . (18)The self-consistency equation (16) transforms into the following relation: ∆ = π T | γ c | (cid:88) ε n sin θ ε n . (19)Assuming an infinitely large superconductor, we consider a spatially independent solution for θ ε n . Then, we find: cos θ ε n = | ε n | (cid:112) ε n + ∆ , sin θ ε n = ∆ (cid:112) ε n + ∆ , (20)where ∆ satisfies the self-consistency relation ∆ = π T | γ c | (cid:88) ε n ∆ (cid:112) ε n + ∆ . (21)This is nothing but the standard self-consistency condition for the gap ∆ in the BCS theory [43].The saddle-point solution (17) can be conveniently written as rotation around the matrix Λ , Q = R − Λ R , R αβ nm = (cid:104) δ ε n ,ε m cos( θ ε n / − t φ δ ε n , − ε m sgn ε m sin( θ ε n / (cid:105) δ αβ . (22)( R − ) αβ nm = (cid:104) δ ε n ,ε m cos( θ ε n / − t φ δ ε n , − ε m sgn ε n sin( θ ε n / (cid:105) δ αβ . We note that the matrix R satisfies the relation, CR T = R − C . We emphasize that although thespatially independent saddle point ansatz (17) solves Usadel equation (15), the solution (17) iscompletely insensitive to disorder and residual interaction between quasiparticles in accordancewith the “Anderson theorem” [2–4]. 6 . Beyond the mean-field approximation In order to see the e ff ect of disorder and residual interactions between quasiparticles, oneneeds to go beyond the mean-field approximation of the previous section. The fluctuations of thematrix Q around the saddle-point ansatz (17) modify the e ff ective potential for the spectral angle θ ε . In what follows, we establish a perturbative renormalization group approach for accountingthese fluctuations. Taking into account the fluctuations of Q , we renormalize the NLSM action. For this pur-pose, we shall develop a peturbation expansion around the saddle point Q , using the square-rootparametrization of the matrix field Q : Q = R − (cid:16) W + Λ √ − W (cid:17) R , W = (cid:32) ww (cid:33) . (23)We adopt the following notations: W n n = w n n and W n n = w n n , where n (cid:62) n < w and w satisfy the charge-conjugation constraints: w = − Cw T C , w = − Cw ∗ C . (24)These constraints imply that some elements ( w αβ n n ) r j in the expansion w αβ n n = (cid:80) r j ( w αβ n n ) r j t r j arepurely real and the others are purely imaginary.In this paper, we restrict our considerations to the expansion of the renormalization groupequations to the lowest order in residual electron-electron interactions. This is justified in the caseof weak short-ranged interaction, which corresponds to small magnitudes of the bare interactionparameters, | γ s | , | γ t | , | γ c | (cid:28)
1. Then, in order to study the e ff ect of fluctuations, we shall needthe propagators for di ff usive modes in the noninteracting theory. In other words, the propagatorsare determined by the NLSM action in which terms S ( ρ )int , S ( σ )int , and ˆ S ( c )int are omitted. Within thisperturbative-in-interaction scheme, one finds (cid:68) [ w r j ( p )] α β n n [ w r j ( − p )] β α n n (cid:69) = g δ α α δ β β δ ε n ,ε n δ ε n ,ε n D p ( i ε n , i ε n ) , r , j = , . . . , , D p ( i ε n , i ε n ) = Dp + E ε n + E ε n , (25) E ε n = | ε n | cos θ ε n + ∆ sin θ ε n . (26)We note that we have not yet fixed θ ε n at this stage. The interaction of the di ff usive modes encoded in W renormalizes the NLSM action. Sincewe consider the spatially independent saddle point, Q , and S ( ρ )int , S ( σ )int , and ˆ S ( c )int are zero at thissaddle point, we are interested in modifications of these terms in the NLSM action by the quan-tum corrections that can be expressed through the di ff usive propagators (25). To the lowest orderin disorder it is enough to approximate Q as Q → Q + R − WR . S ( ρ )int + S ( σ )int → − π T (cid:90) d r (cid:42)(cid:88) α, n (cid:88) r = , (cid:88) j = Γ j Tr (cid:2) RI α n t r j R − W (cid:3) Tr (cid:2) RI α − n t r j R − W (cid:3)(cid:43) . (27)Then, using Eq. (25), we retrieve S ( ρ )int + S ( σ )int = − π T N r Vg (cid:88) n n (cid:90) d d q (2 π ) d D q ( i ε n , i ε n ) × (cid:88) n (cid:40) ( Γ s + Γ t ) δ ε n + ε n ,ω n + ( Γ s + Γ t ) (cid:0) δ ε n − ε n ,ω n − δ ε n + ε n ,ω n (cid:1)(cid:34) cos (cid:16) θ ε n / (cid:17) × cos (cid:16) θ ε n / (cid:17) + sin (cid:16) θ ε n / (cid:17) sin (cid:16) θ ε n / (cid:17)(cid:35) − ( Γ s − Γ t ) δ ε n + ε n ,ω n sin θ ε n sin θ ε n (cid:41) → π T N r Vg ( Γ s − Γ t ) (cid:88) ε,ε (cid:48) > sin θ ε sin θ ε (cid:48) (cid:90) d d q (2 π ) d D q ( i ε, − i ε (cid:48) ) . (28)Here, in the last line, we omitted the term proportional to Γ s + Γ t that involves θ ε n in the di ff usionpropagator only, as it yields a negligible correction with respect to the term we retain. In a similarway, one can renormalize the interaction in the Cooper channel,ˆ S ( c )int → − π T Γ c (cid:88) α, n (cid:88) r = , (cid:90) d r (cid:28)(cid:16) Tr (cid:2) Rt r L α n R − W (cid:3)(cid:17) (cid:29) + π T Γ c V (cid:88) α (cid:88) r = , (cid:42)(cid:32)(cid:90) d r Tr (cid:2) Rt r L α n R − W (cid:3)(cid:33) (cid:43) = − π T Γ c V N r g (cid:88) n , n (cid:90) d d q (2 π ) d D q ( i ε n , i ε n ) (cid:88) n (cid:40) δ ε n − ε n ,ω n + (cid:0) δ ε n + ε n ,ω n − δ ε n − ε n ,ω n (cid:1) × (cid:34) cos (cid:16) θ ε n / (cid:17) cos (cid:16) θ ε n / (cid:17) + sin (cid:16) θ ε n / (cid:17) sin (cid:16) θ ε n / (cid:17)(cid:35)(cid:104) − (2 π ) d δ ( q ) δ n , / V (cid:105)(cid:41) → − π T Γ c N r g (cid:88) ε> D q = ( i ε, − i ε ) sin θ ε . (29)As above, we have neglected here the terms that involve θ ε n in the di ff usion propagator only. Asone can see, the renormalization of ˆ S ( c )int results in the non-extensive term, proportional to V . Inwhat follows, we shall safely neglect this term in the thermodynamic limit V → ∞ .We mention that the above computation is similar to the background field renormalization inthe normal state but with the specific slow field Q (see Supplemental Materials of Ref. [22] andAppendix A of Ref. [23] for details). We also note that to the lowest order in interactions γ s , t , c the renormalization of the parameter Z ω coincides with the renormalization of the Q matrix, i.e.,with the Z factor that determines the renormalization of the average density of states.Equations (28) – (29) determine the correction δ S fl [ Q ] to the classical action due to quan-tum fluctuations (treated to the lowest order in 1 / g ). Including this correction, we obtain theperturbatively modified action in the following form: S [ Q ] + δ S fl [ Q ] = π T Z ω N r V (cid:40) ∆ π T γ c + (cid:88) ε> (cid:104) ε cos θ ε + ∆ sin θ ε (cid:105) + π T ( γ s − γ t ) g (cid:88) ε,ε (cid:48) > sin θ ε sin θ ε (cid:48) (cid:90) d d q (2 π ) d D q ( i ε, − i ε (cid:48) ) (cid:41) . (30)8his action will serve as a basis for the modified equations of motion with the quantum interfer-ence and interaction e ff ects included. Minimizing Eq. (30) with respect to θ ε and ∆ , we find ( ε > − ε sin θ ε + (cid:40) ∆ + π T ( γ s − γ t ) g (cid:88) ε (cid:48) > (cid:90) d d q (2 π ) d (cid:104) D q ( i ε, − i ε (cid:48) ) + D q ( i ε (cid:48) , − i ε ) (cid:105) sin θ ε (cid:48) (cid:41) cos θ ε = , ∆ = π T | γ c | (cid:88) ε (cid:48) > sin θ ε (cid:48) . (31)It is convenient to introduce the energy-dependent quantity ∆ ε = − π T (cid:88) ε (cid:48) > (cid:40) γ c − ( γ s − γ t ) g (cid:90) d d q (2 π ) d (cid:104) D q ( i ε, − i ε (cid:48) ) + D q ( i ε (cid:48) , − i ε ) (cid:105)(cid:41) sin θ ε (cid:48) , (32)that we shall term the gap function . Then, we can write the modified Usadel equation in thefollowing concise form ( ε > − ε sin θ ε + ∆ ε cos θ ε = . (33)Solving Eq. (33), we find the spectral angle,sin θ ε = ∆ ε (cid:112) ε + ∆ ε , cos θ ε = ε (cid:112) ε + ∆ ε . (34)Supplemented by Eq. (34), Eq. (32) becomes a self-consistent equation for the gap function. Wenote that, contrary to the gap function, the Hubbard-Stratonovich field (the order parameter) ∆ remains independent of the energy, ∆ = π T | γ c | (cid:88) ε> ∆ ε (cid:112) ε + ∆ ε . (35)A few remarks are in order here. First, we note that a similar form of the modified Usadelequation has been derived in Refs. [45, 46] in the case of a residual Coulomb interaction andwithout the exchange interaction, i.e., for γ s = − γ t =
0, by means of the diagrammatictechnique . Second, the accuracy of our approximation does not allow us to see that ∆ in thedefinition of E ε , Eq. (26), is changed to ∆ ε . However, there is an argument in favor of suchtransformation from ∆ to ∆ ε . Let us make a shift in the spectral angle θ ε → θ ε + δθ ε . Then, onone hand, δθ ε can be reabsorbed into some W matrix. Therefore, the propagator for δθ ε is relatedwith the propagator for di ff usive modes, Eq. (25). On the other hand, the very same propagatorcan be obtained from the linear variation of the modified Usadel equation. The latter fact impliesthe presence of ∆ ε in the expression for the propagator of di ff usive modes. Therefore, we makethe following substitution in the di ff usive propagator, D − q ( i ε, − i ε (cid:48) ) → q + L − √ ε +∆ ε + √ ε (cid:48) +∆ ε (cid:48) , (36) The Cooperon screening factor of Ref. [45] is given as w ( ε ) = ∆ ε / ∆ . L ε = √ D /ε . We note that for energies larger than the spectral gap edge the transformationfrom ∆ to ∆ ε is not essential.Finally, the form of the second term in the curly brackets in the right hand side of Eq. (32)coincides exactly with the perturbative correction to the interaction in the Cooper channel orig-inating from the interplay of disorder and interactions in the particle-hole channel. Since thiscorrection appeared in a procedure similar to the background field renormalization, we are al-lowed to rewrite Eq. (32) as the following self-consistency equation for ∆ ε , ∆ ε = − π T (cid:88) ε (cid:48) > γ c (cid:18) L √ ε +∆ ε + √ ε (cid:48) +∆ ε (cid:48) (cid:19) ∆ ε (cid:48) (cid:113) ε (cid:48) + ∆ ε (cid:48) , (37)where the dependence of γ c on L is governed by the following equation d γ c dy = − t γ s − γ t ) . (38)Here y = ln L /(cid:96) and t = / ( π g ) denotes the dimensionless resistance. Comparison of Eq. (37)and Eq. (35) suggests that the order parameter ∆ coincides with the value of the gap function atenergy of the order of 1 /τ , i.e. ∆ = ∆ ε ∼ /τ .
5. Critical temperature for the transition to the superconducting state
The modified self-consistency equation (37) allows us to determine the superconducting tran-sition beyond the mean-field approximation. The latter yields the BCS result, T BCS c ∼ τ − exp( − / | γ c | ) , (39)see Eq. (21). c from the renormalization group equations in the normal phase We start by reminding the reader on how the superconducting transition temperature can beestimated from the renormalization group equations in the normal phase. In the case of weakshort-range interactions, these renormalization group equations, to the lowest order in t , read[22, 23, 44] dtdy = t (cid:0) − γ s / − γ t / − γ c (cid:1) , (40a) d γ s dy = − t (cid:0) γ s + γ t + γ c (cid:1) , (40b) d γ t dy = − t (cid:0) γ s − γ t − γ c (cid:1) , (40c) d γ c dy = − t (cid:0) γ s − γ t (cid:1) − γ c , (40d) d ln Z ω dy = t (cid:0) γ s + γ t + γ c (cid:1) . (40e)The renormalization group equations are supplemented by the initial conditions t = t , γ s = γ s , γ t = γ t , and γ c = γ c at y = L = (cid:96) ). We assume that the initial values of resistance and10nteraction parameters are small, t (cid:28) | γ s | , | γ t | , | γ c | (cid:28)
1. In addition, we limit ourdiscussion to the case of t (cid:29) | γ c | . We mention that for t (cid:28) | γ c | disorder modifies the BCStransition temperature only slightly. Under the above assumptions, we can neglect the γ c term inthe r.h.s. of Eq. (40d) at the initial stage of the renormalization group flow. Then, the interactionparameters flow towards the BCS line, − γ s = γ t = γ c = γ. The initial value of γ is equal to γ = (3 γ t + γ c − γ s ) /
6. We assume that γ <
0. Projectingrenormalization group equations (40a) – (40d) onto the BCS line, we find dtdy = t , d γ dy = t γ − γ / . (41)According to these renormalization group equations γ flows towards negative values and divergeseventually at a finite length scale L c . For | γ | (cid:29) t one can estimate this lengthscale as [22] L c = (cid:96) exp (cid:0) / t − / t RG c (cid:1) , (42)where t RG c = t / (2 | γ | ) . (43)The divergence of γ at the finite length scale signals instability of the normal phase towardssuperconducting order at temperature T RG c = D / L c . Hence, we obtain T RG c ∼ τ exp (cid:32) − t + t RG c (cid:33) . (44)We emphasize that T RG c is enhanced parametrically in comparison with the clean BCS transitiontemperature T BCS c given by Eq. (39). We note that at t ∼ (cid:112) | γ | transition into the insulatingphase occurs [22]. We note that Eq. (38) has a striking resemblance with the renormalization group equation for γ c in the normal state, Eq. (40d). As we shall demonstrate below, this fact is not occasional. Asusual, in order to estimate the superconducting transition temperature, we can linearize Eq. (37).Then we find ∆ ε = − π T (cid:88) ε (cid:48) > γ c ( L ε + ε (cid:48) ) ∆ ε (cid:48) ε (cid:48) . (45)In order to make connection of Eq. (45) with the renormalization group equations in thenormal state, following Refs. [45, 46], we substitute L ε + ε (cid:48) by min { L ε , L ε (cid:48) } . We also substitute thesum over Matsubara energies by an integral over continuous energies (see below for argumentsas to why it is possible). Then Eq. (45) can be reduced to the following di ff erential equation, d Υ ε dy ε = − t ( γ s − γ t ) (cid:2) Υ ε − Υ ε = π T (cid:3) , d Υ ε dy ε (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ε = π T = , d ln Υ ε dy ε (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ε = /τ = γ c , (46)where ∆ ε = − (1 / d Υ ε / dy ε and y ε = ln L ε /(cid:96) . Equation (46) is worthwhile to compare with theequation d Υ dy = − t ( γ s − γ t ) Υ , d ln Υ dy (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) y = = γ c , (47)11hat can be obtained from Eq. (40d) upon the substitution γ c = ( d Υ / dy ) / (2 Υ ). Equations (46) –(47) become identical if we impose the additional condition Υ ε = π T =
0. Such condition for Υ ε in Eq. (46) is needed in order to have one-to-one correspondence between Υ ε and ∆ ε . For Υ inEq. (47), the condition Υ = γ c diverges at the length scale L c .Therefore, we can identify Υ ε with Υ ( L ε ): Υ ε ≡ Υ ( L ε ) . (48)Then the condition Υ ε = π T = T c for which L π T c = L c . Tempera-ture T c defined in this way is given by Eq. (44). Thus, the superconducting transition temperaturedetermined from the self-consistency equation coincides with the one found from the renormal-ization group equations in the normal phase. It is worthwhile to mention that there is a subtledi ff erence between Eq. (46) and Eq. (47). The variables t , γ s , and γ t in Eq. (47) are governed byrenormalization group equations (40a) – (40d), whereas in Eq. (46) they obey Eqs. (40a) – (40c)and (38). In the arguments given above we neglected this subtlety. c from the modified self-consistency equation Now we shall compare the prediction for the transition temperature (44) that stems from therenormalization group flow in the normal state with the one obtained from the accurate analysisof the modified self-consistency equation (45). In order to analyse this equation, we project therenormalization group equations onto the BCS line − γ s = γ t = γ c = γ and neglect interactioncorrections to the renormalization of t (see Sec. 5.1). This implies substitution of γ c by γ in Eq.(45), where the length-scale dependence of γ is governed by the following renormalization groupequations: dtdy = t , d γ dy = γ t . (49)Solving the above equation for t ( L ) and γ ( L ), we find t ( L ) = t − t ln L /(cid:96) , γ ( L ) = γ t ( L ) t . (50)Equation (45) can be viewed as an eigenvalue problem, such that the transition temperature isdetermined by the maximum eigenvalue of a corresponding matrix. We parametrize the transitiontemperature as T c = (2 πτ − ) exp( − / t + / t c ) where t c satisfies inequality 1 (cid:29) t c (cid:29) t . Then,rewriting Eq. (45) with the help of Eq. (50), we obtain ∆ n = | γ | t (cid:88) n (cid:48) (cid:62) M n , n (cid:48) (2 / t c ) ∆ n (cid:48) , M n , n (cid:48) ( ζ ) = ζ + ln( n + n (cid:48) + n (cid:48) + / . (51)The above approximate analysis of the self-consistency equation suggests that in the case ζ (cid:29) M ( ζ ) behaves as λ ( M )max ∝ /ζ . The numerical data for λ ( M )max can beviewed in Fig. 1a. The numerical results support the above expectation. For ζ (cid:29)
1, one finds λ ( M )max = u M / (2 ζ ) with u M ≈ .
26. We note that this result satisfies the Perron-Frobenius bound at ζ (cid:29) λ ( M )max < min n (cid:88) n (cid:48) (cid:62) M n , n (cid:48) ( ζ ) = /ζ . (52)12a) � � � � � + � � � � � ζ � � � � � �� ������������ ζ � / λ � � � ( � ) (b) � ��� � ∞ � � ��� ��� ��� ��������������� � / � � � ( � ) Figure 1: (a) Dependence of the inverse maximum eigenvalue, 1 /λ ( M )max on ζ . The dots are numerical values, the red solidline is a fit. Extrapolation to the infinite size of the matrix M is performed. (b) Comparison between the function f ( z ),Eq. (56), (red curve), the function f ( z ) (blue dot-dashed curve) and the eigenvector (blue dots) corresponding to themaximal eigenvalue of the matrix M for t = .
14 and γ = − . Therefore, the self-consistency equation (51) results in the following expression for the super-conducting transition temperature, cf. Eq. (44): T SC c ∼ τ exp (cid:32) − t + t SC c (cid:33) , t SC c = u M t / | γ | . (53)We note that the assumption ζ = / t c (cid:29) ζ is large allows us to find the eigenvector ∆ n that corresponds to the max-imum eigenvalue λ ( M )max . Let us use the Euler-Maclaurin formula for the summation over n (cid:48) inEq. (51). This is justified by the condition ζ (cid:29)
1. In addition, we approximate ln( n + n (cid:48) + { n , n (cid:48) } +
1) (as we shall see below, such simplification provides qualitatively correctresults). With these approximations, we arrive at the following equation: ∆ ( u n ) = u n u (cid:90) u n duu ∆ ( u ) + u n (cid:90) u ∞ du ∆ ( u ) + t | γ | u n ∆ ( u ) , (54)where u n = | γ | t [2 / t c + ln( n + . The quantity u ∞ = | γ | / t (cid:28) u n with the maximally allowed index n (cid:39) / (2 π T c τ ). Let us parametrize ∆ ( u n ) as ∆ ( u n ) = ∆ ( u ) f ( u n ), where a function f ( u ) satisfiesthe normalization condition f ( u ) =
1. Then the function f ( u ) satisfies the following di ff erentialequation, u f (cid:48)(cid:48) ( u ) − f (cid:48) ( u ) + f ( u ) = , f (cid:48) ( u ) = Cu , f (cid:48) ( u ∞ ) = f ( u ∞ ) / u ∞ , (55)where C = t / | γ | (cid:28)
1. Using the smallness of u ∞ , the solution of the above equation can bewritten as f ( u ) = F ( u ) F ( u ) , F ( u ) = uJ (4 √ u ) . (56)Here, J ( x ) denotes the Bessel function of the first kind. The yet unknown parameter u can befound from the boundary condition at u = u , f (cid:48) ( u ) = Cu . To the lowest order in C (cid:28)
1, theparameter u can be found as u = u c − Cu c , (57)13here u c ≈ .
92 is the solution of the equation F (cid:48) ( u c ) =
0. Hence, neglecting the di ff erencebetween u and u c (proportional to t / | γ | (cid:28) T c ∼ τ exp (cid:32) − t + t c (cid:33) , t c = u c t / | γ | . (58)We note that the 25% discrepancy in the numerical values of u c and u M indicates that the approx-imation ln( n + n (cid:48) +
1) by ln(max { n , n (cid:48) } +
1) is not justified by any small parameter. Nevertheless,our approximate solution for the eigenvector with the maximum eigenvalue, Eq. (56), is in goodagreement with the eigenvector determined numerically, as shown in Fig. 1b.The estimate (44) for the superconducting temperature from the analysis of the renormal-ization group equations in the normal phase and the estimates (53) and (58) on the basis of theself-consistency equation are essentially the same except for di ff erent numerical values of theconstant u . In all three cases its value is of the order of unity (1 .
5, 1 .
26, and 0 .
92, respectively).We emphasize also that significant dependence of the gap function ∆ ε on the Matsubara energy ε appears just at the superconducting transition temperature.
6. The energy dependence of the gap function c The non-trivial energy dependence of the gap function at the transition temperature suggeststhat ∆ ( z ) vanishes at T = T c . The dependence of the gap function on temperatures at T c − T (cid:28) T c can be found in a way similar to the BCS theory. We expand the modified self-consistencyequation (37) to the third order in the gap function: ∆ ε = π T (cid:88) ε (cid:48) > | γ ( L ε + ε (cid:48) ) | ∆ ε (cid:48) ε (cid:48) − π T (cid:88) ε (cid:48) > | γ ( L ε + ε (cid:48) ) | ∆ ε (cid:48) ε (cid:48) . (59)We note that the expansion of γ (cid:18) L √ ε +∆ ε + √ ε (cid:48) +∆ ε (cid:48) (cid:19) in powers of ∆ ε and ∆ ε (cid:48) leads to the termswhich are proportional to a small factor t /γ (cid:28)
1. Such terms can be safely neglected. Letus parameterize the gap function as ∆ ε = ∆ ( T ) f ( u ε ) where u ε = | γ | t ( L ε ) / t . Then taking intoaccount that the sum over ε (cid:48) in the last term of the right-hand side of Eq. (59) is dominated by ε (cid:48) ∼ T c , we find f ( u ε ) = u ε u T (cid:90) u ε duu f ( u ) + u ε (cid:90) u ∞ du f ( u ) + C u ε , C = t | γ | − ζ (3) ∆ ( T )8 π T c f ( z T ) , (60)where u T = | γ | t ( L π T ) / t . The above integral equation is reduced to Eq. (55) but with u T insteadof u . Using Eq. (57) and expressing u T and u in terms of T and T c , respectively, we obtain thefollowing well-known result of the BCS theory: ∆ ( T ) = (cid:34) π ζ (3) T c ( T c − T ) (cid:35) / , T c − T (cid:28) T c . (61)We note that the corrections to the BCS-type temperature dependence are controlled by the smallparameter t / | γ | (cid:28)
1. 14 .2. The gap function at low temperatures T (cid:28) T c Let us first analyze the energy dependence of the gap function ∆ ε at the zero temperature.We expect that the function ∆ ε has a form similar to the one at T = T c , see Fig. 1b. We introducethe energy ε which is given by the solution of the equation ε = ∆ ε . Then, we assume that for ε < ε the gap function ∆ ε is close to its value ∆ at ε =
0. For ε > ε the gap function ∆ ε is amonotonously decreasing function that reaches the value ∆ at ε ∼ /τ .At ε < ε the self-consistency equation (37) can be approximately written as ∆ ε = /τ (cid:90) d ε (cid:48) ∆ ε (cid:48) (cid:113) ε (cid:48) + ∆ ε (cid:48) (cid:12)(cid:12)(cid:12) γ ( L ∆ + √ ε (cid:48) +∆ ε (cid:48) ) (cid:12)(cid:12)(cid:12) − (cid:112) ε + ∆ ε − ∆ ∆ + (cid:113) ε (cid:48) + ∆ ε (cid:48) t ( L ∆ + √ ε (cid:48) +∆ ε (cid:48) ) =∆ − (cid:2) (cid:113) ε + ∆ ε − ∆ (cid:3) /τ (cid:90) d ε (cid:48) ∆ ε (cid:48) (cid:113) ε (cid:48) + ∆ ε (cid:48) (cid:12)(cid:12)(cid:12) γ ( L ∆ + √ ε (cid:48) +∆ ε (cid:48) ) (cid:12)(cid:12)(cid:12) t ( L ∆ + √ ε (cid:48) +∆ ε (cid:48) ) ∆ + (cid:113) ε (cid:48) + ∆ ε (cid:48) . (62)Since the dependence of γ on L , governed by Eq. (50), is only logarithmical, the integral over ε (cid:48) in Eq. (62) is dominated by ε (cid:48) ∼ ε . We then find: ∆ ε (cid:39) ∆ − t γ u ∆ + √ ε (cid:20) (cid:113) ε + ∆ ε − ∆ (cid:21) ∞ (cid:90) dx √ + x ( ∆ /ε + √ + x ) . (63)This result implies that ∆ ε is constant at ε < ε up to corrections of the order of t /γ . With thesame accuracy the energy scale ε coincides with ∆ . Further, we can substitute u ∆ for u ∆ + √ ε in Eq. (63). Therefore, at ε < ∆ the gap function behaves as ∆ ε (cid:39) ∆ − t γ u ∆ (cid:20) (cid:113) ε + ∆ − ∆ (cid:21) . (64)At ε (cid:62) ε (cid:39) ∆ we approximate the self-consistency equation (37) as follows: ∆ ε = (cid:12)(cid:12)(cid:12) γ ( L ε ) (cid:12)(cid:12)(cid:12) ε (cid:90) d ε (cid:48) ∆ (cid:113) ε (cid:48) + ∆ + (cid:12)(cid:12)(cid:12) γ ( L ε ) (cid:12)(cid:12)(cid:12) ε (cid:90) ε d ε (cid:48) ∆ ε (cid:48) ε (cid:48) + /τ (cid:90) ε d ε (cid:48) ∆ ε (cid:48) ε (cid:48) (cid:12)(cid:12)(cid:12) γ ( L ε (cid:48) ) (cid:12)(cid:12)(cid:12) . (65)Here, the dependence of γ on L is governed by Eq. (50). We note that Eq. (65) is justifiedfor ε (cid:29) ε . For ε (cid:39) ε , a more accurate treatment of the self-consistency equation (37) resultsin corrections of the order of t / | γ | . Parametrizing the energy dependence of the gap functionat ε (cid:62) ε as ∆ ε = ε f ( u ε ), we find that the integral equation (65) results in Eq. (55) with u ε entering in place of u and with the value of the constant C = c t / | γ | , where c = arcsinh(1).We note that, in order to find the precise value of c , one needs to take into account the t / | γ | corrections to Eq. (65). However, we are not interested in such accuracy and set C to zero. Then,we find that the function f ( u ) is given by Eq. (56) with u c replacing u . With the same accuracy,we have ε = ∆ (cid:39) T c , cf. Eq. (58).Expressing the right-hand side of Eq. (65) in terms of the function f ( z ) and setting ε ∼ /τ ,we retrieve ∆ = c ∆ ( γ / t ) ∆ , c ∆ = u c (cid:90) duu f ( u ) ≈ . . (66)15 �� � �� �� �� �� ���������������� �� ( ε / Δ � ) Δ ε / Δ � Figure 2: Dependence of ratio ∆ ε / ∆ on ln( ε/ ∆ ) at T = γ = − .
005 and t = .
03. The monotonously decreasingpart of the curve corresponds to the function f ( z ε ). We note that ∆ (cid:28) ∆ . All in all, we can summarize the energy dependence of the spectral gap at T = ∆ ε = ∆ , ε < ∆ , F (cid:0) | γ | t ( L ε ) / t (cid:1) F ( u c ) , ε (cid:62) ∆ , ∆ ∼ τ exp − t + | γ | u c t . (67)The overall dependence of the gap function ∆ ε on ε is shown in Fig. 2.We emphasize that within our approximate treatment, the values of ∆ and T c coincide withthe exponential accuracy, i.e. the exponential factor in Eq. (67) coincides with the exponentialfactor in Eq. (58). However, we cannot exclude a possibility that the magnitude of the ratio ∆ / T c di ff ers from the one in the BCS theory. This can result from the corrections of the orderof t / | γ | that determine the coe ffi cient C in Eq. (55), see Eq. (57).In order to consider the e ff ect of nonzero but low temperatures, T (cid:28) T c , we need to performsummation over Matsubara frequencies in the interval 0 < ε n (cid:46) ∆ more accurately. We assumethat the parameter ∆ becomes temperature dependent, ∆ → ∆ ( T ). Then Eq. (65) should bemodified as follows: ∆ ε = | γ ( L ε ) | π T (cid:88) ε (cid:48) > ∆ ( T ) ε (cid:48) + ∆ ( T ) − /τ (cid:90) ∆ ( T ) d ε (cid:48) ∆ ( T ) (cid:113) ε (cid:48) + ∆ ( T ) + ε (cid:90) ε d ε (cid:48) ε (cid:48) | γ ( L ε ) | ∆ ε (cid:48) + /τ (cid:90) ε d ε (cid:48) ε (cid:48) | γ ( L ε (cid:48) ) | ∆ ε (cid:48) . (68)Using the parametrization, ∆ ε = ∆ ( T ) f ( u ε ), one can check that the function f ( u ) satisfies Eq.(55) with u replaced by u ( T ) = | γ | t ( L ∆ ( T ) ) / t . The constant C becomes a T -dependent func-tion given by C ( T ) (cid:39) C ( T = − t | γ | ∞ (cid:90) dE (cid:113) E + ∆ e − √ E +∆ / T = C ( T = − t | γ | (cid:114) π T ∆ e − ∆ / T . (69)16sing Eq. (57), we can express u ( T ) in terms of u , yielding ∆ ( T ) = ∆ − (2 π T ∆ ) / e − ∆ / T , T (cid:28) T c . (70)Therefore, the dependence of the magnitude of the gap function on temperature T (cid:28) T c at ε < ∆ is the same as in the BCS theory. We note that there are corrections to this result which are ofthe order of t / | γ | .
7. Disorder-averaged density of states
In this section, we analyze the disorder-averaged density of states of single-particle excita-tions in the superconducting state. Within the NLSM formalism, the average density of statescan be obtained as ρ ( E ) = ν (cid:104) Q αα nn (cid:105) , (71)where the analytic continuation to real energies, ε n → E + i
0, is performed. Here symbol ‘sp’denotes the trace over particle-hole and spin spaces. Plugging in the parametrization (23) intoEq. (71), one finds ρ ( E ) = ν Re Z / ε cos θ ε (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ε →− iE + = ν Re Z / ε ε (cid:112) ε + ∆ ε (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ε →− iE + . (72)The factor Z ε can be written as Z ε = Z ( L √ ε +∆ ε ), where [35] d ln Zd ln y = t ( γ s + γ t + γ c ) . (73)Projecting Eq. (73) to the BCS line, one finds that Z ε (cid:39) + γ ( L √ ε +∆ ε ) . (74)This factor encodeds the interaction-induced corrections to the density of states of the type thatleads to the zero-bias anomaly in the normal state [35].As one can see, in order to compute the average density of states, we need to make an analyticcontinuation for the gap function ∆ ε from Matsubara energies to real energies. We emphasizethat the nontrivial dependence of ∆ ε on ε implies that the spectral gap at real energies, ∆ ( E ),has both real and imaginary parts. The symmetry implies that Re ∆ ( E ) is an even function of E whereas Im ∆ ( E ) is an odd function of E . ∆ ( E )We restrict our consideration to the case of zero temperature, T =
0. Performing analyticcontinuation in Eq. (64), we find that the gap function at E (cid:54) ∆ is purely real, ∆ ( E ) (cid:39) ∆ − t u γ (cid:20) (cid:113) ∆ − E − ∆ (cid:21) . (75)At energies E (cid:38) ∆ , the imaginary part of ∆ ( E ) appears. To show that the nonzero imaginary partdoes exist, one can consider the perturbative expression (32). Performing analytic continuation17o the real energy, we find that Im ∆ ( E ∼ /τ ) ∼ ∆ t . We note that the imaginary part of ∆ ( E )is nonzero due to the imaginary part of the renormalized interaction in the Cooper channel, γ c ,which corresponds to the superconducting fluctuation propagator.For ε (cid:38) ∆ , we energy dependence of the gap function is described by Eq. (67). Afteranalytic continuation, i ε → E + i
0, in Eq. (67), we obtain ∆ ( E ) = ∆ f (˜ u E ) , (76)where ˜ u E = | γ | ˜ t ( L E ) t = | γ | / t + ( t /
2) ln( E τ ) − i π t / (cid:39) u E + i π t | γ | u E . (77)Here, we have used the fact that t ( L E ) (cid:28)
1. The real and imaginary parts of ∆ ( E ) read:Re ∆ ( E ) =∆ f | γ | t ( L E ) t , (78a)Im ∆ ( E ) =∆ π | γ | t ( L E )4 t f (cid:48) | γ | t ( L E ) t . (78b)Expression (78a) has the same range of validity as Eq. (67). It provides us with an estimatefor the real part of the gap function at E (cid:62) ∆ withcorrections of the order of t / | γ | neglected.The situation with the imaginary part of the gap function is more delicate. Equation (78b) meansthat Im ∆ ( E ∼ ∆ ) ∼ ∆ t /γ . We note that the contribution of the same order will be given bythe continuation of Eq. (75) to energies larger than ∆ . This fact implies that, in order to computeIm ∆ ( E ) at energies E ∼ ∆ , one needs to retain terms of the order of t /γ in the solution of Eq.(37). However, at large energies, E (cid:29) ∆ , Eq. (78b) provides the leading result to the imaginarypart of the gap function. For example, at E ∼ /τ , Eq. (78b) yields Im ∆ ( E ∼ /τ ) ∼ ∆ t thatmatches with the perturbative result. Additional argument that the more accurate (to the order t /γ ) expression for the real part of ∆ ( E ) is needed in order to determine Im ∆ ( E ) is given bythe Kramers-Kronig relations for imaginary and real parts of ∆ ( E ) (see Appendix B).The above results suggest that there is an energy E g such that Im ∆ ( E ) = E (cid:54) E g and Im ∆ ( E ) > E > E g . Making analytic continuation to the real frequencies in the self-consistency equation (37), we obtain ∆ ( E ) = π T (cid:88) ε (cid:48) > (cid:12)(cid:12)(cid:12) γ ( L √ ∆ ( E ) − E + √ ε (cid:48) +∆ ε (cid:48) ) (cid:12)(cid:12)(cid:12) ∆ ε (cid:48) (cid:113) ε (cid:48) + ∆ ε (cid:48) . (79)For energies close to the energy E g , we can expand the right hand side of Eq. (79) as ∆ ( E ) = ∆ ( E g ) − α (cid:34) (cid:112) ∆ ( E ) − E − (cid:113) ∆ ( E g ) − E (cid:35) , (80a) α = π T (cid:88) ε (cid:48) > (cid:12)(cid:12)(cid:12) γ ( L √ ∆ ( E g ) − E + √ ε (cid:48) +∆ ε (cid:48) ) (cid:12)(cid:12)(cid:12) t ( L √ ∆ ( E g ) − E + √ ε (cid:48) +∆ ε (cid:48) ) ∆ ε (cid:48) ε (cid:48) + ∆ ε (cid:48) . (80b)We shall demonstrate below that the parameter α ∼ t /γ (cid:28)
1. Solving Eq. (80a) for α (cid:28)
1, wefind the dependence of the gap function for real energies close to E g , ∆ ( E ) = ∆ ( E g ) − α (cid:113) E − E , ∆ ( E g ) = E g (1 + α / , | E − E g | (cid:28) E g . (81)18n agreement with our assumptions, Im ∆ ( E ) = E (cid:54) E g . We note that ∆ ( E ) (cid:62) E for E (cid:54) E g .For E > E g the imaginary part of the gap function is non-zero, Im ∆ ( E ) = − i α (cid:113) E − E . Wenote that negative sign of Im ∆ ( E ) at energies E > E g is needed for positivity of the density ofstates (see below). However, away from E g the imaginary part of Im ∆ ( E ) has to change sign inorder to be consistent with the asymptotic at large energies, Eq. (78b). The change of the sign ofIm ∆ ( E ) occurs at the energy of the order of a few E g .Now we can estimate the parameter α . The comparison of Eq. (81) with Eq. (75) impliesthat E g is of the order of ∆ . The sum over ε (cid:48) in Eq. (80b) is dominated by ε (cid:48) (cid:39) ∆ ε (cid:48) (cid:39) ∆ . Thenwe find α (cid:39) (cid:12)(cid:12)(cid:12) γ ( L ∆ ) (cid:12)(cid:12)(cid:12) t ( L ∆ ) ∞ (cid:90) d ε (cid:48) ∆ ε (cid:48) + ∆ = π u t γ . (82)Here, we have taken into account that (cid:113) ∆ ( E g ) − E = α E g (cid:28) ∆ . With the help of the above results for ∆ ( E ), we can derive the expression for disorder–averaged density of states as a function of energy. For large energies, E (cid:29) ∆ , the imaginarypart of ∆ ( E ) can be neglected in comparison with its real part, see Eqs. (78a) and (78b). Sincefor E (cid:29) ∆ the energy is always larger than the real part of ∆ ( E ), we find ρ ( E ) = ν E (cid:113) E − ∆ f (cid:16) | γ | t ( L E ) / t (cid:17) , ∆ (cid:28) E (cid:46) /τ . (83)Here and hereinafter we neglect the factor Z ε , see Eq. (72), since it provides negligible correc-tions of the order of t / | γ | in the parametric regime that we consider.In order to compute ρ ( E ) at energies close to E g (cid:39) ∆ , we use the expression (81) for ∆ ( E ).Introducing ∆ E = E − E g , we find for | ∆ E | (cid:28) E g , ρ ( E ) = να , ∆ E < , Im (cid:20) − ∆ E α E g − i (cid:113) ∆ E α E g (cid:21) − , ∆ E (cid:62) . (84)As expected, the energy E g determines the gap edge in the disorder-averaged density of states.Using Eq. (84), we derive a square-root dependence of the density of state near the gap edge: ρ ( E ) = √ α (cid:115) E − E g E g , < E − E g (cid:28) α E g . (85)The square-root growth of ρ ( E ) turns into the maximum at E = E max = E g (1 + α / ρ ( E max ) = ν/ (2 α ) ∼ νγ / t (cid:29) ν . For energies larger than E max , the density of states decays, ρ ( E ) = ν √ (cid:115) E g E − E g , α E g (cid:28) E − E g (cid:28) E ∗ . (86)19 � � � � �� � ( � - � � ) α � � � ���� ρ ( � ) ρ ( � ��� ) � � � * � � * ��� ~ γ �� / � �� ρ ( � )/ ν Figure 3: Sketch of the dependence of the local density of states ρ ( E , r ) /ν on the energy E for | γ | (cid:28) t (cid:28) (cid:112) | γ | . Thered solid curve corresponds to the disorder-averaged density of states (see text in Sec. 7.2). The dashed red curves andthe shaded region illustrates the mesoscopic fluctuations of the local density of states. There are no quasiparticle statesbelow E g . The energy E g ∗ separates the region of strong and weak mesoscopic fluctuations. The energy E ∗ demarcatesthe energy regions where the mesoscopic fluctuations are cut o ff by the length scale L ∗ ( E g < E < E ∗ ) and the dephasinglength L φ ( E ∗ < E ) (see text in Sec. 8.2). Inset: the region close to the gap edge E g . The curve corresponds to Eq. (84). We note that Eq. (86) matches the high energy asymptotics, Eq. (83), at energy of the order of E g . The dependence of the disorder-averaged density of states on energy is shown in Fig. 3. (SeeSec. 8.2 for discussion of the mesoscopic fluctuations of the local density of states.) We notethat the profile of ρ ( E ) is very similar to the density of states in the presence of the depairingterm in the Usadel equation, introduced, e.g., by the scattering o ff magnetic impurities [47]. Thecomparison of Eq. (85), as well as the position of the maximum and its magnitude, with theAbrikosov-Gor’kov theory suggests that the e ff ective depairing parameter η ∼ α . We note thatthe depairing term in the Usadel equation appears also in the model with a spatially varying, ran-dom interaction in the Cooper channel [48] and with a spatially varying random order parameter[49].
8. Mesoscopic fluctuations of the local density of states in the superconducting state
We start analysis of the mesoscopic fluctuations of the local density of states from the calcu-lation of its dispersion. It can be expressed in terms of bilinear in Q operators [50], (cid:104) [ δρ ( E , r )] (cid:105) = ν
32 Re (cid:104) P ( i ε n , i ε n ) − P ( i ε n , i ε n ) (cid:105)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) i ε n , n → E + i i ε n → E − i , (87)where δρ ( E , r ) = ρ ( E , r ) − (cid:104) ρ ( E , r ) (cid:105) and P ( i ε n , i ε m ) = (cid:104)(cid:104) sp Q α α nn ( r ) · sp Q α α mm ( r ) (cid:105)(cid:105) − (cid:104) sp (cid:2) Q α α nm ( r ) Q α α mn ( r ) (cid:3) (cid:105) . (88)20ere α (cid:44) α are some fixed replica indices and (cid:104)(cid:104) A · B (cid:105)(cid:105) = (cid:104) AB (cid:105)−(cid:104) A (cid:105)(cid:104) B (cid:105) . Using parametrization(23) for Q in Eq. (88) and expanding to the second order in W , we obtain P ( i ε n , i ε n ) = − g (cid:90) d d q (2 π ) d (cid:40) sin (cid:32) θ ε n (cid:33) cos (cid:32) θ ε n (cid:33) D q ( i ε n , − i ε n ) + cos (cid:32) θ ε n (cid:33) sin (cid:32) θ ε n (cid:33) D q ( i ε n , − i ε n ) (cid:41) (89)and P ( i ε n , i ε n ) = − g (cid:90) d d q (2 π ) d (cid:40) cos (cid:32) θ ε n (cid:33) cos (cid:32) θ ε n (cid:33) D q ( i ε n , i ε n ) + sin (cid:32) θ ε n (cid:33) sin (cid:32) θ ε n (cid:33) D q ( − i ε n , − i ε n ) (cid:41) . (90)By means of Eqs. (34) and (36), we find (for arbitrary signs of ε and ε (cid:48) ) P ( i ε, i ε (cid:48) ) = − g (cid:90) d d q (2 π ) d − ε (cid:112) ε + ∆ ε ε (cid:48) (cid:113) ε (cid:48) + ∆ ε (cid:48) DDq + (cid:112) ε + ∆ ε + (cid:113) ε (cid:48) + ∆ ε (cid:48) . (91)Making an analytic continuation to real frequencies in accordance with the prescription in Eq.(87), we obtain the following result for the dispersion of fluctuations of the local density of states: (cid:104) [ δρ ( E , r )] (cid:105) = ν g Re (cid:90) d d q (2 π ) d (cid:34)(cid:32) + E | ∆ ( E ) − E | (cid:33) DDq + (cid:112) ∆ ( E ) − E − (cid:32) + E ∆ ( E ) − E (cid:33) DDq + (cid:112) ∆ ( E ) − E (cid:35) . (92)At energies below the gap edge, E (cid:54) E g , the gap function is real and satisfy ∆ ( E ) (cid:62) E . There-fore, the fluctuations of the local density of states are zero identically, (cid:104) [ δρ ( E , r )] (cid:105) = , E (cid:54) E g . (93)Above the superconducting gap, E > E g the mesoscopic fluctuations are non-zero. At ener-gies close to the gap, 0 < E − E g (cid:28) E g , we find (cid:104) [ δρ ( E , r )] (cid:105) (cid:39) t ρ ( E ) ln min { L g , L } (cid:96) , (94)where L stands for the system size and the average density of states, ρ ( E ) is given by Eq. (84).The length L g is defined as L g = (cid:113) D / ( α E g ) ∼ ( | γ | / t ) L ∆ . (95)Finally, at large energies, E (cid:29) ∆ , the general expression (92) can be simplified as (cid:104) [ δρ ( E , r )] (cid:105) (cid:39) t ρ ( E ) ln min { L ∗ , L φ , L } (cid:96) − (Re ∆ ( E )) E ln min { L ∗ , L φ , L } L √ E − (Re ∆ ( E )) . (96)21e have neglected Im ∆ ( E ) everywhere except the denominator of the first di ff usive propagatoron the right-hand side of Eq. (92). Here, the length L ∗ is defined through D / L ∗ = | Im ∆ ( E ) | / (cid:112) E − (Re ∆ ( E )) . (97)Since Im ∆ ( E ) is small, the length scale L ∗ is large, L ∗ (cid:29) L ∆ . In particular, at E ∼ ∆ one canfind the following estimate: L ∗ ∼ L ∆ | γ | / t ∼ L g .We note that the di ff usive propagators in Eq. (92) are a ff ected by dephasing due to electron-electron interactions [35]. This results in appearance of the dephasing length L φ in Eq. (96). Forenergies E (cid:29) ∆ the dephasing length can be estimated as D / L φ ∼ t ( L E ) γ ( L E ) E . (98)At the spectral edge, E = E g , L φ diverges due to restriction of the phase volume for quasiparticles(there are no quasiparticles below E g ). Therefore, for energies close to the spectral edge, E − E g (cid:28) E g , the dephasing length exceeds the length L g . This is the reason why L φ is absent in Eq.(94).The above results (93)–(96) can be summarized as (cid:104) [ δρ ( E , r )] (cid:105) (cid:39) t ρ ( E ) ln min { L , L ∗ , L g , L φ } (cid:96) . (99)We remind the reader that the lengths L g , L ∗ , and L φ are defined in Eqs. (95), (97), and (98),respectively. We mention that, with a logarithmic accuracy, t ln(min { L ∗ , L g , L φ } /(cid:96) ) ∼ ∆ . This indicates that the perturbative treatment of mesoscopic fluctuationshas to be extended to a more elaborated renormalization group approach. The perturbative result (99) has exactly the same form as the perturbative correction tothe dispersion of the local density of states above the superconducting transition. Essentially,the superconducting state leads to the specific infrared cut o ff length scale, min { L ∗ , L g , L φ } ,only. This suggests that one can use the one-loop renormalization group equation for m ≡(cid:104) [ ρ ( E , r )] (cid:105) /ρ ( E ) derived in Ref. [35], d ln m dy = t + O ( t ) (100)upto the length scale min { L ∗ , L g , L φ } . Solving Eqs. (100) and (49), we find (cid:104) [ ρ ( E , r )] (cid:105) ρ ( E ) (cid:39) (cid:34) t (min { L , L ∗ , L g , L φ } ) t (cid:35) . (101)We note that the length scale L ∗ starts from the value of the order of L g and, then, it grows withincreasing E . The dephasing length L φ has the opposite behavior: it decreases with increasingenergy. Therefore, there is an energy E ∗ at which the lengths L ∗ and L φ become of the sameorder. Using DL ∗ = ∆ E t | γ | u E f ( u E ) f (cid:48) ( u E ) , DL φ = t | γ | u E E , (102)22here u E = | γ | t ( L E ) / t , we can find the following estimate, E ∗ ∼ ∆ t | γ | ln | γ | t / . (103)We note that E ∗ (cid:29) ∆ . Substituting expressions (102) for the length scales L ∗ , and L φ , into Eq.(101), we obtain the following energy dependence of the disorder-averaged second moment ofthe local density of states (in the case L = ∞ ), (cid:104) [ ρ ( E , r )] (cid:105) = u c t γ ρ ( E ) + (2 u c t / | γ | ) ln( | γ | / t ) , E ∼ ∆ , − ( u c t / | γ | ) ln[ ∆ u E f ( u E ) f (cid:48) ( u E ) t / ( E | γ | )] , ∆ (cid:28) E (cid:28) E ∗ , (cid:16) + [ u c t / (2 | γ | )] ln[ t u E E / ( | γ | ∆ )] (cid:17) − , E ∗ (cid:28) E (cid:28) /τ . (104)Interestingly, in accordance with Eq. (104), the disorder-averaged normalized second moment ofthe local density of states has the maximum at E = E ∗ . The magnitude of this maximum can beestimated as (cid:104) [ ρ ( E ∗ , r )] (cid:105) ρ ( E ∗ ) ρ ( E ∼ ∆ ) (cid:104) [ ρ ( E ∼ ∆ , r )] (cid:105) − ∼ ( u c t / | γ | ) ln( | γ | / t ) . (105)Since the right hand side of the relation (105) is much smaller than one, we can approximate thedisorder-averaged second moment of the local density of states at E g < E < E ∗ as (cid:104) [ ρ ( E , r )] (cid:105) = ( u c t /γ ) ρ ( E ) . (106)This indicates the strong mesoscopic fluctuations of the local density of states at energies E g < E < E ∗ . With further increase of energy, the disorder-averaged second moment of the localdensity of states decreases. Using Eq. (104), we find (cid:104) [ ρ ( E ∼ /τ, r )] (cid:105) (cid:39) ρ ( E ∼ /τ ) . (107)Thus, the mesoscopic fluctuations of ρ ( E , r ) are suppressed at the ultraviolet energy scale, E ∼ /τ .The dependence of (cid:104) [ ρ ( E , r )] (cid:105) on the energy is illustrated in Fig. 3, where the functions ρ ( E ) ± (cid:112) (cid:104) [ ρ ( E , r )] (cid:105) are shown by dashed red curves. The fluctuations (not normalized) decreaseswith increase of the energy. They are strong in the energy interval, E g < E < E ∗ . At E ∗ thedependence of (cid:104) [ ρ ( E , r )] (cid:105) on the energy is changed, as it was explained above. We note that thereis a certain energy E g ∗ which is the solution of the following equation: (cid:112) (cid:104) [ ρ ( E , r )] (cid:105) = ρ ( E ).Using Eq. (104), we find that E g ∗ ∼ (cid:32) ∆ τ (cid:33) / ∼ τ exp − t + | γ | u c t . (108)As shown in Fig. 3, for energies E < E g ∗ the fluctuations of the local density of states are strong.For energies E > E g ∗ , the mesoscopic fluctations are not enough to reduce significantly the localdensity of states below its average value. 23 .3. Distribution function for the local density of states In a similar way as it was done for the temperatures above T c [35], one can generalize theresult (104) to the higher moments of the local density of states, (cid:104) [ ρ ( E , r )] q (cid:105) = ρ q ( E ) (cid:2) A ( E ) (cid:3) q ( q − / , A ( E ) = (cid:104) [ ρ ( E , r )] (cid:105) /ρ ( E ) . (109)We note that according to Eq. (104), the function A ( E ) > E g < E < E g ∗ ,i.e. for region of strong mesoscopic fluctuations. The expression (109) implies the followinglog-normal distribution for the normalized local density of states, ˜ ρ = ρ ( E , r ) /ρ ( E ), (see Ref.[51] for the case of a normal metal) P ( ˜ ρ ) = A ( E ) √ π ln A ( E ) exp −
12 ln A ( E ) (cid:32) ln ˜ ρ +
32 ln A ( E ) (cid:33) , E g < E < E g ∗ . (110)The log-normal distribution (110) predicts that the most probable value for the local density ofstates is given as ρ mode ( E ) = ρ ( E ) (cid:104) [ ρ ( E , r )] (cid:105) / . (111)This result implies that the most probable height of the coherence peak can be estimated as ρ mode ( E ∼ ∆ ) /ν ∼ | γ | / t . We note that for | γ | (cid:28) t (cid:28) | γ | / the magnitude of the coherencepeak in ρ mode is much larger than ν . In the region | γ | / (cid:28) t (cid:28) | γ | / the coherence peak isabsent.We also introduce the typical value of the normalized local density of states, ˜ ρ typ = exp (cid:104) ln ˜ ρ (cid:105) where (cid:104) . . . (cid:105) denotes the average with respect to the distribution (110). Then we obtain, ρ typ ( E ) = ρ ( E ) (cid:104) [ ρ ( E , r )] (cid:105) / . (112)This result leads to the following estimate: ρ typ ( E ∼ ∆ ) /ν ∼ | γ | / t . For | γ | (cid:28) t (cid:28) | γ | / ,the magnitude of the coherence peak in ρ type is much larger than the bare value of the density ofstates. For | γ | / (cid:28) t (cid:28) | γ | / , the coherence peak is absent. We note that, strictly speaking,the typical value of ˜ ρ should be determined not from the distribution (110) but from the modifieddistribution in which rare events of extremely small or extremely large values of ˜ ρ are suppressed(see Ref. [52]). However, as one can check, such a modification does not significantly modifythe result (112).We mention that the following relation holds: ρ mode ( E ) (cid:28) ρ typ ( E ) (cid:28) ρ ( E ) . (113)This means that for energies E g < E < E g ∗ , an experimentally measured local density of stateswill be typically much smaller than its average value given by Eqs. (83) and (84). Nevertheless,for | γ | (cid:28) t (cid:28) | γ | / a measured local density of states will typically have a high coherencepeak and the spectral gap of the order of E g . For | γ | / (cid:28) t (cid:28) | γ | / , the coherence peak ina measured local density of states will be completely suppressed and the density of states willhave a soft spectral gap of the order of E g ∗ . In both ranges of t , a measured local density ofstate at some spatial points will have much higher coherence peak than in the typical spatialregions. For energies E > E g ∗ , an experimentally measured local density of states will be closeto the disorder-averaged density of states (83). In this energy interval, E > E g ∗ , the mesoscopic24uctuations of the local density of states are small. This physical picture for the local density ofstates is illustrated in Fig. 3.Following Ref. [35], one can also generalize the result (104) to correlation functions of thelocal density of states at di ff erent spatial points and di ff erent energies [53]. Finally, we note thatthe fact of strong mesoscopic fluctuations of the local density of states is consistent with thestrong mesoscopic fluctuations of the superconducting order parameter, see Appendix A.
9. Summary and conclusions
To summarize, we have developed the theory of the multifractal superconducting state inthin films. Treating the fluctuations around the mean-field spatially homogeneous solution, wederived the modified Usadel equation that incorporates the interplay of disorder and interactionsat energy scales larger than the spectral gap.Our key findings are as follows:(i) The modified Usadel equation, in combination with the self-consistency equation, yieldsparametrically the same estimate for the multifractally enhanced superconducting transi-tion temperature as the one derived by considering the instability in renormalization groupequations in the normal phase.(ii) The mutual e ff ects of disorder and interactions result in strong dependence of the super-conducting gap function on energy (see Fig. 2): at energies of the order of the spectralgap the gap function ∆ ε is parametrically enhanced in comparison with its magnitude atultraviolet energies ∼ /τ .(iii) The spectral gap at zero temperature is multifractally enhanced in the same way as thesuperconducting transition temperature.(iv) The energy dependence of the gap function in the Usadel equation results in the profileof the disorder-averaged density of states that, near the spectral gap, resembles the onederived in the model of a spatially random superconducting order parameter (see Fig. 3).We stress that the interplay of disorder and interactions leads to two opposite e ff ects. Onthe one hand, it results in the enhancement of the spectral gap, but on the other hand, itinduces the e ff ective depairing parameter that cuts o ff the coherence peaks in the averagedensity of states. The corresponding depairing parameter is estimated as ∼ ( t /γ ) (cid:28) E < E g ∗ (see Eq. (108)). In the energy interval E g < E < E ∗ (see Fig. 3 and Eq.(103)) their relative amplitude is of the order of t / | γ | (cid:29) ff ect resultingfrom quantum interference in a disordered system.Based on these findings for the statistics of fluctuations of the local density of states inthin films with multifractally-enhanced superconductivity, we conclude that disorder-induced25nterference e ff ects dramatically a ff ect spectral properties of these superconducting films, byfully governing the physics in a wide energy interval where the spectral gap for single-particleexcitations establishes. Specifically, we have identified a parametrically large energy range, E g < E < E g ∗ , where the quasiparticle spectral gap can be zero in some spatial regions andnon-zero in the other. The distribution function of the local density of states has a log-normalform, such that the typical value of the density of states is lower that the average value. In otherwords, the system may locally look as superconducting at energies much higher than E g .Our results for strong mesoscopic fluctuations of the local density of states are in qualitativeagreement with tunneling spectroscopy data on thin superconducting films [27–32] and withnumerical solution of disordered attractive two-dimensional Hubbard model [33, 34]. The strongmesoscopic fluctuations of the local density of states are accompanied by strong mesoscopicfluctuations of the superconducting order parameter.On the technical side, within our approach, we have integrated over the spatial fluctuations ofthe superconducting order parameter from the very beginning. Therefore, they are hidden in theterm ˆ S ( c )int that describes the interaction in the Cooper channel. The renormalization group anal-ysis of the NLSM in the normal phase suggests that such a procedure is more convenient, sinceit allows one to take into account the interplay between the spatial fluctuations of the supercon-ducting order parameter (the interaction in the Cooper channel) and charge and spin fluctuations(interaction in the particle-hole channel). This interplay leads to strong renormalization of theNLSM action at length scales of the order of L T c . In our approach, the mesoscopic fluctuationsof the local superconducting order parameter manifests themselves in the course of the renormal-ization group flow as the energy dependence of the gap function.Since we considered a macroscopically homogeneous sample, it was natural to address fluc-tuations around the spatially homogeneous mean-field solution. The emergent fluctuations of theorder parameter encoded in the mesoscopic fluctuations of the density of states, as observed inexperiments of the “gap tomography”, are due to mesoscopic quantum interference e ff ects abovethe gap-edge energy. At the same time, it is interesting to see whether nonperturbative tails of thedensity of states inside the gap would appear in spatially homogeneous disordered films studiedwithin the NLSM formalism. For this purpose, one should consider a possibility of existenceof non-trivial saddle-point solutions of the NLSM action. A more straightforward way to obtainsuch tails would be through the introduction of macrocroscopic inhomogeneities (say, in the localconcentration of impurities) that could induce “background” fluctuations of the superconductinggap.We restricted our consideration to the case of a weak short-ranged interaction (suppressedby, e.g., high dielectric constant of the substrate), for which a strong e ff ect of the enhancementof superconductivity by multifractality was predicted. The approach developed in this paper canbe generalized in several directions. Our theory can be extended to include a strong short-rangedinteraction, as well as the Coulomb interaction. One can also study the multifractal superconduct-ing state that occurs in the system that without superconducting instability is near the interactingmetal-insulator transition. Our theory can also be extended to the case of an applied magneticfield that destroys the superconducting state and may produce an intermediate phase with giantmagnetoresistance [54] (for studying the giant magnetoresistance near the transition within theNLSM formalism, see Ref. [23]). Finally, our approach does not take into account phase fluc-tuations of the order parameter, giving rise to the Berezinskii-Kosterlitz-Thouless phenomena insuperconducting films, as well as the existence of the vortices with a normal state in the core.Such fluctuations can be incorporated into our theory in the way similar to the one in Ref. [42].We note, however, that for su ffi ciently high conductance of the normal state, such e ff ects are26xpected to influence the results described in this article only slightly.
10. Acknowledgements
The authors are grateful to M. Skvortsov and K. Tikhonov for useful discussions. I.S.B.is grateful to M. Stosiek and F. Evers for collaboration on a related project. The research waspartially supported by the Russian Foundation for Basic Research (grant No. 20-52-12013) –Deutsche Forschungsgemeinschaft (grant No. EV 30 / Appendix A. Mesoscopic fluctuations of the superconducting order parameter
In this Appendix we estimate the mesoscopic fluctuations of the superconducting order pa-rameter. While the quantity that is actually measured in experiments on the tomography of thesuperconductors is the local density of states (Sec. 8), the results of the scanning-tunneling-microscopy measurements are frequently translated into the maps of the fluctuating local orderparameter by fitting the density of states with the BCS “ansatz”. Here, instead of translating ourresults for the mesoscopic fluctuation of the density of states into the fluctuations of the orderparameter, we calculate the latter directly within the NLSM formalism. The relation (9) suggeststhat the variance of the order parameter can be written as follows: (cid:104) δ ∆ r ( r ) (cid:105) = (cid:18) π T γ c (cid:19) (cid:34)(cid:68)(cid:68) Tr (cid:2) t r L α Q ( r ) (cid:3) · Tr (cid:2) t r L α Q ( r ) (cid:3)(cid:69)(cid:69) − (cid:68) Tr (cid:2) t r L α Q ( r ) t r L α Q ( r ) (cid:3)(cid:69)(cid:35) , (A.1)where α (cid:44) α are some fixed replica indices. We note that the operator in the brackets in Eq.(A.1) is the eigenoperator under renormalization group flow. This can be easily checked with thehelp of the following identities: (cid:104) Tr AW Tr BW (cid:105) = Y Tr (cid:104) AB − Λ A Λ B − ACB T C + A Λ CB T C Λ (cid:105) , (A.2a) (cid:104) Tr AW BW (cid:105) = Y (cid:104) Tr A Tr B − Tr Λ A Tr Λ B + Tr ACB T C − Tr A Λ CB T C Λ (cid:105) , (A.2b)that follow from Eq. (25). Here, Y = ( t /
2) ln( L /(cid:96) ) and we have neglected the energy dependenceof the di ff usive propagators at scales L (cid:29) L T c .Using parametrization (23), expanding the right hand side of Eq. (A.1) to the second orderin W , and applying Eq. (25), we obtain the following perturbative result: (cid:104) ( δ ∆ ) (cid:105) ∆ = g (cid:88) ε> sin θ ε − (cid:88) ε,ε (cid:48) > sin θ ε sin θ ε (cid:48) (cid:90) d d q (2 π ) d D q ( i ε, − i ε (cid:48) ) (A.3)Here we have used the self-consistency equation (31). We note that the fluctuations δ ∆ do notlead to a finite single-particle density of states below the gap edge E g but rather correspondto the spatial fluctuations of the superconducting condensate. The hard spectral gap for thequasiparticles results from an e ff ective averaging over such fluctuations with the self-consistencyor renormalization group procedure.With the logarithmic accuracy, we can estimate the right-hand side of Eq. (A.3) as (cid:104) ( δ ∆ ) (cid:105) ∆ (cid:39) t ln min { L T , L ∆ } (cid:96) . (A.4)27sing Eqs. (53) and (67), the perturbative result (A.4) implies that for temperatures T (cid:46) T c themesoscopic fluctuations of the superconducting order parameter are large, (cid:104) ( δ ∆ ) (cid:105) ∼ ∆ . Thisis consistent with our results from Sec. 8.1, where the mesoscopic fluctuations of the density ofstates were calculated, cf. Eq. (99). Including renormalization e ff ects as in the calculation ofthe mesoscopic fluctuations of the density of states in Sec. 8.2, one can generalize Eq. (A.4) ina similar manner, to obtain yet stronger fluctuations of the superconducting order parameter inthin films with multifractally-enhanced superconductivity [53]. Appendix B. Kramers–Kronig relations for ∆ ( E ) In this appendix, we demonstrate how the Kramers-Kronig relations for ∆ ( E ),Im ∆ ( E ) = − p . v . ∞ (cid:90) −∞ d ωπ Re ∆ ( ω ) ω − E , Re ∆ ( E ) = p . v . ∞ (cid:90) −∞ d ωπ Im ∆ ( ω ) ω − E . (B.1)are satisfied by the approximate expressions for the real and imaginary parts of ∆ ( E ). Let usrewrite the Kramers-Kronig relation for Re ∆ ( E ) asRe ∆ ( E ) = p . v . ∞ (cid:90) d ωπω (cid:2) Im ∆ ( ω + E ) + Im ∆ ( ω − E ) (cid:3) . (B.2a)Expressing the imaginary part of ∆ ( E ) in terms of its real part in accordance with Eq. (78b),Im ∆ ( ω ) = − ( π/ ω∂ ω Re ∆ ( ω ) , (B.3)we findRe ∆ ( E ) = − p . v . ∞ (cid:90) d ω (cid:104) ∂∂ω (cid:2) Re ∆ ( ω + E ) + Re ∆ ( ω − E ) (cid:3) + E ω ∂∂ E (cid:2) Re ∆ ( ω + E ) − Re ∆ ( ω − E ) (cid:3)(cid:105) = Re ∆ ( E ) + π E ∂∂ E Im ∆ ( E ) . (B.4)We note that here we have neglected the contribution to Im ∆ ( E ) from energies close to thespectral gap E g . As one can see, when the approximation (B.3), which is valid for E (cid:29) ∆ , isused for all energies, an additional contribution to the real part of ∆ ( E ) appears, which is of theorder of ( t /γ ) z ( z f (cid:48) ( z )) (cid:48) . The latter is of the order of t /γ at E ∼ ∆ . This implies that, inorder to be able to compute the imaginary part of ∆ ( E ) from the Kramers-Kronig relation, oneneeds to derive expression for Re ∆ ( E ) with the accuracy of the order of t /γ .28 eferences [1] P.W. Anderson, Absence of di ff usion in certain random lattices , Phys. Rev. , 1492 (1958).[2] A. A. Abrikosov and L. P. Gor’kov, On the theory of superconducting alloys, I. The electrodynamics of alloys atabsolute zero , Sov. Phys. JETP , 1090 (1959).[3] A. A. Abrikosov and L. P. Gor’kov, Superconducting alloys at finite temperatures , Sov. Phys. JETP , 220 (1959).[4] P. W. Anderson, Theory of dirty superconductors , J. Phys. Chem. Solids , 26 (1959).[5] L. N. Bulaevskii and M. V. Sadovskii, Localization and superconductivity , JETP Lett. , 640 (1984).[6] M. Ma and P. A. Lee, Localized superconductors , Phys. Rev. B , 5658 (1985).[7] A. Kapitulnik and G. Kotliar, Anderson localization and the theory of dirty superconductors , Phys. Rev. Lett. ,473 (1985).[8] G. Kotliar and A. Kapitulnik, Anderson localization and the theory of dirty superconductors. II , Phys. Rev. B ,3146 (1986).[9] S. Maekawa and H. Fukuyama, Localization e ff ects in two-dimensional superconductors , J. Phys. Soc. Jpn. ,1380 (1982).[10] H. Takagi and Y. Kuroda, Anderson localization and superconducting transition temperature in two-dimensionalsystems , Solid State Comm. , 643 (1982).[11] S. Maekawa, H. Ebisawa, and H. Fukuyama, Theory of dirty superconductors in weakly localized regime , J. Phys.Soc. Jpn. , 2681 (1984).[12] P. W. Anderson, K. A. Muttalib, and T. V. Ramakrishnan, Theory of the “universal” degradation of T c in high-temperature superconductors , Phys. Rev. B , 117 (1983).[13] L. N. Bulaevskii and M. V. Sadovskii, Anderson localization and superconductivity , J. Low Temp. Phys. , 89(1985).[14] A. M. Finkel’stein, Superconducting transition temperature in amorphous films , JETP Lett. , 46 (1987).[15] A. M. Finkel’stein, Suppression of superconductivity in homogeneously disordered systems , Physica B 197, 636(1994).[16] D. B. Haviland, Y. Liu, and A. M. Goldman,
Onset of superconductivity in the two-dimensional limit , Phys. Rev.Lett. , 2180 (1989).[17] A. M. Goldman and N. Markovi´c, Superconductor- insulator transitions in the two-dimensional limit , Phys. Today51, 39 (1998)[18] V. F. Gantmakher and V. T. Dolgopolov,
Superconductor-insulator quantum phase transition , Physics-Uspekhi , 1 (2010).[19] B. Sac´ep´e, M. Feigel’man, T. M. Klapwijk, Quantum breakdown of superconductivity in low-dimensional mate-rials , Nat. Phys. , 734 (2020).[20] M. V. Feigel’man, L. B. Io ff e, V. E. Kravtsov, and E. A. Yuzbashyan, Eigenfunction fractality and pseudogap statenear the superconductor-insulator transition , Phys. Rev. Lett. , 027001 (2007).[21] M. V. Feigel’man, L. B. Io ff e, V. E. Kravtsov, and E. Cuevas, Fractal superconductivity near localization thresh-old , Ann. Phys. , 1390 (2010).[22] I.S. Burmistrov, I.V. Gornyi, A.D. Mirlin,
Enhancement of the critical temperature of superconductors by Ander-son localization , Phys. Rev. Lett. , 017002 (2012).[23] I. S. Burmistrov, I. V. Gornyi, A. D. Mirlin,
Superconductor-insulator transitions: Phase diagram and magne-toresistance , Phys. Rev. B , 014506 (2015).[24] M. N. Gastiasoro and B. M. Andersen, Enhancing superconductivity by disorder , Phys. Rev. B , 184510 (2018).[25] K. Zhao , H. Lin, X. Xiao, W. Huang, W. Yao, M. Yan, Y. Xing, Q. Zhang, Z.-X. Li, S. Hoshino, J. Wang, S.Zhou, L. Gu, M. S. Bahramy, H. Yao, N. Nagaosa, Q.-K. Xue, K. T. Law, X. Chen, and S.-H. Ji, Disorder-inducedmultifractal superconductivity in monolayer niobium dichalcogenides , Nat. Phys. , 904 (2019).[26] C. Rubio-Verd´u, A. M. Garc´ıa-Garc´ıa, H. Ryu, D.-J. Choi, J. Zald´ıvar, S. Tang, B. Fan, Z.-X. Shen, S.-K. Mo, J. I.Pascual, and M. M. Ugeda, Visualization of multifractal superconductivity in a two-dimensional transition metaldichalcogenide in the weak-disorder regime , Nano Lett. , 5111 (2020).[27] B. Sac´ep´e, C. Chapelier, T. I. Baturina, V. M. Vinokur, M. R. Baklanov, and M. Sanquer, Disorder-inducedinhomogeneities of the superconducting state close to the superconductor-insulator transition , Phys. Rev. Lett. , 157006 (2008).[28] B. Sac´ep´e, C. Chapelier, T. I. Baturina, V. M. Vinokur, M. R. Baklanov, and M. Sanquer,
Pseudogap in a thin filmof a conventional superconductor , Nat. Commun. , 140 (2010).[29] B. Sac´ep´e, Th. Dubouchet, C. Chapelier, M. Sanquer, M. Ovadia, D. Shahar, M. Feigel’man, and L. Io ff e, Local-ization of preformed Cooper pairs in disordered superconductors, Nat. Phys. , 239 (2011).[30] D. Sherman, B. Gorshunov, S. Poran, N. Trivedi, E. Farber, M. Dressel, and A. Frydman, E ff ect of Coulomb in-teractions on the disorder-driven superconductor-insulator transition , Phys. Rev. B , 035149 (2014).[31] M. Mondal, A. Kamlapure, M. Chand, G. Saraswat, S. Ku- mar, J. Jesudasan, L. Benfatto, V. Tripathi, and . Raychaudhuri, Phase fluctuations in a strongly disordered s-wave NbN superconductor close to the metal-insulator transition , Phys. Rev. Lett. , 047001 (2011).[32] Y. Noat, V. Cherkez, C. Brun, T. Cren, C. Carbillet, F. Debontridder, K. Ilin, M. Siegel, A. Semenov, H.-W.H¨ubers, and D. Roditchev,
Unconventional superconductivity in ultrathin superconducting NbN films studied byscanning tunneling spectroscopy , Phys. Rev. B , 014503 (2013).[33] B. Fan and A. M. Garc´ıa-Garc´ıa, Enhanced phase-coherent multifractal two-dimensional superconductivity , Phys.Rev. B , 104509 (2020).[34] M. Stosiek, B. Lang, and F. Evers,
Self-consistent-field ensembles of disordered Hamiltonians: E ffi cient solverand application to superconducting films , Phys. Rev. B , 144503 (2020).[35] I. S. Burmistrov, I. V. Gornyi, and A. D. Mirlin, Local density of states and its mesoscopic fluctuations near thetransition to a superconducting state in disordered systems , Phys. Rev. B , 205432 (2016).[36] A. M. Finkel’stein, Electron Liquid in Disordered Conductors , vol. 14 of Soviet Scientific Reviews, ed. by I.M.Khalatnikov, Harwood Academic Publishers, London, (1990).[37] D. Belitz and T. R. Kirkpatrick,
The Anderson-Mott transition , Rev. Mod. Phys. , 261 (1994).[38] I. S. Burmistrov, Finkel’stein nonlinear sigma model: interplay of disorder and interaction in 2D electron systems ,JETP, , 669 (2019).[39] M. V. Feigel’man, A. I. Larkin, and M. A. Skvortsov, Phys. Rev. B , 12361 (2000).[40] I. V. Yurkevich and I. V. Lerner, Phys. Rev. B , 064522 (2001).[41] A. Levchenko and A. Kamenev, Phys. Rev. B , 094518 (2007).[42] E. J. Koenig, A. Levchenko, I. V. Protopopov, I. V. Gornyi, I. S. Burmistrov, and A. D. Mirlin, B erezinskii–Kosterlitz–Thouless transition in homogeneously disordered superconducting films, Phys. Rev. B , 214503(2015).[43] J. Bardeen, L.N. Cooper, and J.R. Schrie ff er, Phys. Rev. , 1175 (1957).[44] L. Dell’Anna, Enhancement of critical temperatures in disordered bipartite lattices , Phys. Rev. B , 195139(2013).[45] M. A. Skvortsov and M. V. Feigel’man, Superconductivity in Disordered Thin Films: Giant Mesoscopic Fluctua-tions , Phys. Rev. Lett. , 057002 (2005).[46] M. V. Feigel’man and M. A. Skvortsov, Universal Broadening of the Bardeen-Cooper-Schrie ff er Coherence Peakof Disordered Superconducting Films , Phys. Rev. Lett. , 147002 (2012).[47] A. A. Abrikosov and L. P. Gor’kov, Contribution to the theory of superconducting alloys with paramagneticimpurities , Zh. Eksp. Teor. Fiz. , 1781 (1960).[48] A. I. Larkin and Yu. N. Ovchinnikov, Density of states in inhomogeneous superconductors , Sov. Phys. JETP ,1144 (1972).[49] M. A. Skvortsov and M. V. Feigel’man, Subgap states in disordered superconductors , JETP , 487 (2013).[50] I. S. Burmistrov, I. V. Gornyi, and A. D. Mirlin,
Multifractality at Anderson transitions with Coulomb interaction ,Phys. Rev. Lett. , 066601 (2013).[51] I. V. Lerner,
Distribution functions of current density and local density of states in disordered quantum conductors ,Phys. Lett. A , 253 (1988).[52] M. S. Foster, S. Ryu, and A.W. W. Ludwig,
Termination of typical wavefunction multifractal spectra at the Ander-son metal-insulator transition: Field theory description using the functional renormalization group , Phys. Rev. B , 075101 (2009).[53] I. S. Burmistrov, I. V. Gornyi, and A. D. Mirlin, to be published elsewhere.[54] A. Datta, A. Banerjee, N. Trivedi, and A. Ghosal, New paradigm for a disordered superconductor in a magneticfield , arXiv:2101.00220 ., arXiv:2101.00220 .