{'fraction_non_alphanumeric': 0.048410912670510474, 'fraction_numerical': 0.017947546055407115, 'mean_word_length': 4.196291560102302, 'pattern_counts': {'":': 0, '<': 28, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 19, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Multifractally-enhanced superconductivity in thin films

The multifractal superconducting state originates from the interplay of Anderson localization and interaction effects. In this article we overview the recent theory of the superconductivity enhancement by multifractality and extend it to describe the spectral properties of superconductors on the scales of the order of the superconducting gap. Specifically, using the approach based on renormalization group within the nonlinear sigma model, we develop the theory of a multifractal superconducting state in thin films. We derive a modified Usadel equation that incorporates the interplay of disorder and interactions at energy scales larger than the spectral gap and study the effect of such an interplay on the low-energy physics. We determine the spectral gap at zero temperature which occurs to be proportional to the multifracally enhanced superconducting transition temperature. The modified Usadel equation results in the disorder-averaged density of states that, near the spectral gap, resembles the one obtained in the model of a spatially random superconducting order parameter. We reveal strong mesoscopic fluctuations of the local density of states in the superconducting state. Such strong mesoscopic fluctuations imply that the interval of energies in which the superconducting gap establishes is parametrically large in systems with multifractally-enhanced superconductivity.


Introduction
Anderson localization [1] and superconductivity are two fundamental quantum phenomena. Semiclassicaly, i.e., without taking into account the quantum interference effects ( and also discarding the influence of disorder on interaction), superconductivity is not affected by the electron scattering on non-magnetic disorder that is compatible with the symmetry of the superconducting order parameter (so-called "Anderson theorem" [2][3][4]). However, at the next level of sophistication -with quantum effects included, the Anderson theorem meets Anderson localization. At this level, superconductivity and Anderson localization were considered as antagonists: strong localization [5][6][7][8] or interplay of weak disorder and Coulomb interaction [9][10][11][12][13][14][15][16][17][18] were predicted to destroy superconductivity 1 . This view was supported by a discovery [19] and the further study of the superconductor-insulator transition (SIT) in thin films (see Refs. [20][21][22] for a review).
The paradigm of suppression of superconductivity by Anderson localization has been challenged in Refs. [23,24], 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 threedimensional disordered systems) has been demonstrated for the situation when the long-ranged Coulomb repulsion between electrons is not effective. Later, the multifractal enhancement of T c in thin superconducting films has been predicted by the present authors [25,26]. Recently, these theoretical predictions have been supported by numerical solution of the disordered attractive Hubbard model on a two-dimensional [27] and three-dimensional [28] lattices. Recently, an enhancement of the superconducting transition temperature by disorder observed in monolayer niobium dichalcogenides has been explained by the multifractality [29,30].
The hallmark of the multifractally-enhanced superconductivity in disordered systems is strong mesoscopic fluctuations of the local order parameter and the local density of states [24]. The point-to-point fluctuations of the local density of states have been observed in many experiments on tunneling spectroscopy on thin superconducting films [31][32][33][34][35][36]. Recently, these fluctuations have been retrieved from numerical solution of disordered attractive two-dimensional Hubbard model [37,38].
For temperatures above the critical temperature, T > T c , the point-to-point fluctuations of the local density of states observed in the tunneling experiments are in qualitative agreement with the theory [39] developed by the present authors for the mesoscopic fluctuations of the local density of states in the normal phase of disordered superconducting films. Our theory employs the renormalization-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 disorder and interactions in all channels (particle-hole as well as particle-particle) and to describe systems 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 as recent numerical data on the disordered Hubbard model represent a new challenge for an analytical theory. In the present paper, we develop a theory of the mesoscopic fluctuations of the local 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 dielectric constant, and consider the regime of an intermediate disorder strength [25,26], in which the superconducting transition temperature is parametrically enhanced by multifractality compared to the conventional Bardeen-Cooper-Schrieffer (BCS) result. We calculate the average density of states and its mesoscopic fluctuations in the low-temperature limit in the presence of the interplay between disorder and interactions.
We demonstrate that the gap in the disorder-averaged density of states at zero temperature is proportional to the superconducting transition temperature, i.e., it is also enhanced by multifractality. In spite of the enhancement of the spectral gap, we find that the combined effect of disorder and interaction results in the suppression of the coherence peaks. Their height is finite and controlled by the ratio of a bare disorder and interaction. We also show that the mesoscopic fluctuations of the local density of states in the superconducting state are enhanced in the same way as in the normal phase for temperatures above T c .
On the technical side, in order to describe the superconducting state, we derive the Usadel equation modified by renormalization effects due to interplay of disorder and interactions (in both particle-hole and particle-particle channels) at short scales. This modified Usadel equation can 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 derived by means of the renormalization group in the normal phase. We discuss the relation between the approaches based on the renormalization of the sigma model and the self-consistent solution of the 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 the Finkel'stein NLSM as applied to superconducting systems [15,18,25,26]. The mean-field description of the superconducting state is briefly outlined in Sec. 3. The description of the superconducting state beyond the mean-field approximation is developed in Sec. 4. The developed theory is used to estimate the superconducting transition temperature (Sec. 5) and the energy dependence of the gap function (Sec. 6). The results for the local density of states and its mesoscopic fluctuations in the superconducting state are presented in Sec. 7 and Sec. 8, respectively. Finally, our results and conclusions are summarized in Sec. 9.

Finkel'stein NLSM action in the normal state
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 singlet channel), S (σ) int (the particle-hole triplet channel), and S (c) int (the particle-particle channel) (see Refs. [40][41][42] for a review): where dr Tr I α n t r0 Q Tr I α −n t r0 Q, (2b) n r=0,3 j=1,2,3 dr Tr I α n t r j Q Tr I α −n t r j Q, Here, Q(r) is the matrix field in the replica, Matsubara, spin, and particle-hole spaces. The trace Tr acts in the same spaces. The matrix field obeys the nonlinear constraint, as well as chargeconjugation symmetry relation: The nonlinear constraint on the matrix field Q suggests the following parametrization: where α, β = 1, . . . , N r stand for replica indices and integers n, m correspond to the Matsubara fermionic frequencies ε n = πT (2n + 1). The representation (4) implies that the diffusive fluctuations in the normal state are associated with smooth variations around the metallic saddle point Q = Λ.
The action (1) contains three constant matrices: The sixteen matrices, t r j = τ r ⊗ s j , r, j = 0, 1, 2, 3, operate in the particle-hole (subscript r) and spin (subscript j) spaces. The matrices τ 0 , τ 1 , τ 2 , τ 3 and s 0 , s 1 , s 2 , s 3 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 2 /h and including spin). The parameter Z ω describes the frequency renormalization upon the renormalization group flow [40]. The bare value of Z ω is equal to πν/4, where ν denotes the bare density of states at the Fermi level including a spin-degeneracy factor. The interaction amplitude Γ s (Γ t ) encodes interaction in the singlet (triplet) particle-hole channel. The interaction in the Cooper channel is denoted as Γ c . Its negative magnitude, Γ c < 0, corresponds to an attraction in the particle-particle channel. In what follows, it is convenient to introduce the dimensionless interaction parameters, γ s,t,c = Γ s,t,c /Z ω . If Coulomb interaction is present, the following relation holds: γ s = −1. This condition remains intact under the action of the renormalization group flow [26,40]. In the case of the BCS model and, under an assumption of strong disorder, ω D τ 1 (ω D is the Debye frequency), the bare values of interaction amplitudes are related as −γ s0 = γ t0 = γ c0 . In what follows, we shall refer to such a relation between interaction amplitudes as the BCS line.

Finkel'stein NLSM action in the superconducting state
The symmetry breaking that leads to the appearance of the superconducting state changes the saddle point of the NLSM [43][44][45][46]. We shall follow the approach of Ref. [46] which allows us to take into account the effects of disorder and residual quasiparticle interactions beyond the mean-field level.
Let us single out the term with n = 0 from the expression for S (c) int , Eq. (2d). Then, upon the Hubbard-Stratonovich decoupling by spatially dependent fields ∆ α r (r), r = 1, 2, we find whereS (c) Since, by construction, the action depends quadratically on ∆ α r (r), we can solve the equation of motion instead of performing the functional integral. The equation of motion reads We note that the charge-conjugation symmetry, Eq. (3), guarantees that ∆ α 1,2 (r) are real functions. 4 The presence of nonzero ∆ α r (r) changes the saddle-point equation for Q. Performing variation of the NLSM action (1) with S (c) int replaced byS (c) int , we derive the following saddle-point equation: Here, γ 0 ≡ γ s and γ 1 = γ 2 = γ 3 ≡ γ t , and we have introduced the diffusion coefficient D = g/(16Z ω ). We mention that in the absence of interaction Eq. (10) is nothing but the Usadel equation, 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 spatially dependent solutions for Q(r) which mimic fluctuations of a disorder potential in the original microscopic formulation of the problem. To enumerate all these solutions and to perform summation over them seem like a daunting task. In order to circumvent this difficulty we shall use the smallness of the bare resistance, 1/g 1. Then, we can sum over spatially dependent solutions for Q(r) by treating the fluctuations around the spatially independent solution of the Usadel equation by means of the renormalization group. We start from splitting ∆ α r into constant and spatially dependent part, The spatially varying part, δ∆ α r (r), contains information about the mesoscopic fluctuations of the superconducting order parameter. These mesoscopic fluctuations turn out to be large (see Appendix A) [24,37,38,52]. In order to take them into account, we perform a formally exact integration over δ∆ α r (r). As a result, the spatial fluctuations of the Hubbard-Stratonovich field corresponding to the order parameter get fully encoded in additional correlations of the Q field. This procedure results in a modification of S (c) int , where V denotes the volume of a superconductor, After integrating over δ∆ α r (r), the saddle-point equation (10) gets modified, It is worth noting that the information about the disorder-induced spatial fluctuations of the Hubbard-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 coupling of∆ and Q and the fluctuations of field Q. Before considering the effect of such fluctuations, in the following section we analyse the mean-field solution of the modified Usadel equation.
The self-consistency equation (16) transforms into the following relation: Since the superconductor order parameter ∆ is spatially independent by construction, we consider a spatially independent solution for θ ε n . Then, we find: where ∆ satisfies the self-consistency relation This is nothing but the standard self-consistency condition for the gap ∆ in the BCS theory [47]. The saddle-point solution (17) can be conveniently written as rotation around the matrix Λ, (R −1 ) αβ nm = δ ε n ,ε m cos(θ ε n /2) − t φ δ ε n ,−ε m sgn ε n sin(θ ε n /2) δ αβ .
We note that the matrix R satisfies the relation, CR T = R −1 C. We emphasize that although the spatially independent saddle point ansatz (17) solves Usadel equation (15), the solution (17) is completely insensitive to disorder and residual interaction between quasiparticles in accordance with the "Anderson theorem" [2][3][4]. 6

Beyond the mean-field approximation
In order to see the effect of disorder and residual interactions between quasiparticles, one needs to go beyond the mean-field approximation of the previous section. The fluctuations of the matrix Q around the saddle-point ansatz (17) modify the effective potential for the spectral angle θ ε . In what follows, we establish a perturbative renormalization group approach for accounting these fluctuations.

Perturbative expansion
Taking into account the fluctuations of Q, we renormalize the NLSM action. For this purpose, we shall develop a peturbation expansion around the saddle point Q, using the square-root parametrization of the matrix field Q: We adopt the following notations: W n 1 n 2 = w n 1 n 2 and W n 2 n 1 = w n 2 n 1 , where n 1 0 and n 2 < 0. The blocks w and w satisfy the charge-conjugation constraints: These constraints imply that some elements (w αβ n 1 n 2 ) r j in the expansion w αβ n 1 n 2 = r j (w αβ n 1 n 2 ) r j t r j are purely real and the others are purely imaginary.
We note that we have not yet fixed θ ε n at this stage.

Renormalization of the NLSM action
The interaction of the diffusive modes encoded in W renormalizes the NLSM action. Since we consider the spatially independent saddle point, Q, and S (ρ) int , S (σ) int , andŜ (c) int are zero at this saddle point, we are interested in modifications of these terms in the NLSM action by the quantum corrections that can be expressed through the diffusive propagators (25). To the lowest order in disorder it is enough to approximate Q as This results in Then, using Eq. (25), we retrieve Here, in the last line, we omitted the term proportional to Γ s +3Γ t that involves θ ε n in the diffusion propagator only, as it yields a negligible correction with respect to the term we retain. In a similar way, one can renormalize the interaction in the Cooper channel, n δ ε n 1 −ε n 2 ,ω n + δ ε n 1 +ε n 2 ,ω n − δ ε n 1 −ε n 2 ,ω n × cos 2 θ ε n 1 /2 cos 2 θ ε n 2 /2 + sin 2 θ ε n 1 /2 As above, we have neglected here the terms that involve θ ε n in the diffusion propagator only. As one can see, the renormalization ofŜ (c) int results in the non-extensive term, proportional to V 0 . In what 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 in the normal state but with the specific slow field Q (see Supplemental Materials of Ref. [25] and Appendix A of Ref. [26] 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 quantum fluctuations (treated to the lowest order in 1/g). Including this correction, we obtain the perturbatively modified action in the following form: This action will serve as a basis for the modified equations of motion with the quantum interference and interaction effects included.

Modified Usadel equation
Minimizing Eq. (30) with respect to θ ε and ∆, we find (ε > 0) It is convenient to introduce the energy-dependent quantity that we shall term the gap function. Then, we can write the modified Usadel equation in the following concise form (ε > 0), Solving Eq. (33), we find the spectral angle, Supplemented by Eq. (34), Eq. (32) becomes a self-consistent equation for the gap function. We note that, contrary to the gap function, the Hubbard-Stratonovich field (the order parameter) ∆ remains independent of the energy, A few remarks are in order here. First, we note that a similar form of the modified Usadel equation has been derived in Refs. [52,53] in the case of a residual Coulomb interaction and without the exchange interaction, i.e., for γ s = −1 and γ t = 0, by means of the diagrammatic technique 2 . Second, the accuracy of our approximation does not allow us to see that ∆ in the definition of E ε , Eq. (26), is changed to ∆ ε . However, there is an argument in favor of such transformation from ∆ to ∆ ε . Let us make a shift in the spectral angle θ ε → θ ε + δθ ε . Then, on one hand, δθ ε can be reabsorbed into some W matrix. Therefore, the propagator for δθ ε is related with the propagator for diffusive modes, Eq. (25). On the other hand, the very same propagator can be obtained from the linear variation of the modified Usadel equation. The latter fact implies the presence of ∆ ε in the expression for the propagator of diffusive modes. Therefore, we make the following substitution in the diffusive propagator, 2 The Cooperon screening factor of Ref. [52] is given as w(ε) = ∆ ε /∆.

9
where L ε = √ D/ε. We note that for energies larger than the spectral gap edge the transformation from ∆ 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 originating from the interplay of disorder and interactions in the particle-hole channel. Since this correction appeared in a procedure similar to the background field renormalization, we are allowed to rewrite Eq. (32) as the following self-consistency equation for ∆ ε , where the dependence of γ c on L is governed by the following equation Here y = ln L/ and t = 2/(π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 at energy of the order of 1/τ, i.e. ∆ = ∆ ε∼1/τ .

Critical temperature for the transition to the superconducting state
The modified self-consistency equation (37) allows us to determine the superconducting transition beyond the mean-field approximation. The latter yields the BCS result, see Eq. (21).

Estimate for T c from the renormalization group equations in the normal phase
We start by reminding the reader on how the superconducting transition temperature can be estimated from the renormalization group equations in the normal phase. In the case of weak short-range interactions, these renormalization group equations, to the lowest order in t, read [25,26,48] dt It is worthwhile to mention that, within the lowest order in t, the dependence of the renormalization-group functions on the interaction parameters can be computed exactly [40,41]. There are subtleties (extensively discussed in the past, see e.g. Ref. [16][17][18]) related with the renormalization of the Cooper-channel interaction. Recently, the difficulties with the Cooper-channel renormalization have been resolved within the background-field renormalization scheme, and the renormalization-group functions for interaction parameters γ s,t,c have been computed exactly (including dependence on γ c ) within the lowest order in t [26] (see also Ref. [46]). For the purpose of this paper, it will be sufficient to use the equations (40a) -(40e) which are of the lowest order in γ s,t,c (and are not affected by the above subtleties).
The renormalization group equations are supplemented by the initial conditions t = t 0 , γ s = γ s0 , γ t = γ t0 , and γ c = γ c0 at y = 0 (L = ). We assume that the initial values of resistance and interaction parameters are small, t 0 1 and |γ s0 |, |γ t0 |, |γ c0 | 1. In addition, we limit our discussion to the case of t 0 |γ c0 |. We mention that for t 0 |γ c0 | disorder modifies the BCS transition temperature only slightly. Under the above assumptions, we can neglect the γ 2 c term in the r.h.s. of Eq. (40d) at the initial stage of the renormalization group flow. Then, the interaction parameters flow towards the BCS line, The initial value of γ is equal to γ 0 = (3γ t0 + 2γ c0 − γ s0 )/6. We assume that γ 0 < 0. Projecting renormalization group equations (40a) -(40d) onto the BCS line, we find Let us emphasize that the disorder-induced correction to the renormalization group equation (the term 2tγ in Eq. (41) is intimately related to the multifractality of wave functions in the noninteracting theory. The coefficient 2t coincides with the lowest-order term in the expansion of the absolute value of the multifractal exponent |∆ 2 | in powers of the resistance t [49]. As known [50], the exponent ∆ 2 < 0 governs the scaling behavior of the disorder-averaged fourth moment of the wave functions. In the NLSM formalism, the scaling of the fourth moment of the wave functions is governed by an operator bilinear in Q. Therefore, to the linear order in interaction, the renormalization-group equations of the interacting theory are controlled by scaling dimensions at the noninteracting fixed point [25,51]. According to the renormalization group equations (41), γ flows towards negative values and diverges eventually at a finite length scale L c . For |γ 0 | t 2 0 one can estimate this lengthscale as [25] L where The divergence of γ at the finite length scale signals instability of the normal phase towards superconducting order at temperature We emphasize that T RG c is enhanced parametrically in comparison with the clean BCS transition temperature T BCS c given by Eq. (39). This enhancement occurs due to the positive sign of the disorder-induced term 2tγ in Eq. (41). As explained above, this is a consequence of the multifractality in the noninteracting theory. Therefore, the origin of the enhancement of the superconducting transition temperature is deeply rooted into the multifractality of two-dimensional disordered electrons. Also we note that at t 0 ∼ |γ 0 | transition into the insulating phase occurs [25].

Relation between renormalization group flow and the modified self-consistency equation
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. As usual, in order to estimate the superconducting transition temperature, we can linearize Eq. (37). Then we find In order to make connection of Eq. (45) with the renormalization group equations in the normal state, following Refs. [52,53], we substitute L ε+ε by min{L ε , L ε }. We also substitute the sum over Matsubara energies by an integral over continuous energies (see below for arguments as to why it is possible). Then Eq. (45) can be reduced to the following differential equation, where ∆ ε = −(1/2)dΥ ε /dy ε and y ε = ln L ε / . Equation (46) is worthwhile to compare with the equation that 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 Υ in Eq. (47), the condition Υ = 0 appears automatically, since γ c diverges at the length scale L c . Therefore, we can identify Υ ε with Υ(L ε ): Then the condition Υ ε=πT = 0 determines a certain temperature T c for which L πT c = L c . Temperature T c defined in this way is given by Eq. (44). Thus, the superconducting transition temperature determined from the self-consistency equation coincides with the one found from the renormalization group equations in the normal phase. It is worthwhile to mention that there is a subtle difference between Eq. (46) and Eq. (47). The variables t, γ s , and γ t in Eq. (47) are governed by renormalization group equations (40a) -(40d), whereas in Eq. (46) they obey Eqs. (40a) -(40c) and (38). In the arguments given above we neglected this subtlety.

Estimate for T c from the modified self-consistency equation
Now we shall compare the prediction for the transition temperature (44) that stems from the renormalization group flow in the normal state with the one obtained from the accurate analysis of the modified self-consistency equation (45). In order to analyse this equation, we project the renormalization group equations onto the BCS line −γ s = γ t = γ c = γ and neglect interaction 12 corrections 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 group equations: Solving the above equation for t(L) and γ(L), we find Equation (45) can be viewed as an eigenvalue problem, such that the transition temperature is determined by the maximum eigenvalue of a corresponding matrix. We parametrize the transition temperature as T c = (2πτ −1 ) exp(−2/t 0 + 2/t c ) where t c satisfies inequality 1 t c t 0 . Then, rewriting Eq. (45) with the help of Eq. (50), we obtain The above approximate analysis of the self-consistency equation suggests that in the case ζ 1, the maximum eigenvalue of M(ζ) behaves as λ (M) max ∝ 1/ζ. The numerical data for λ (M) max can be viewed in Fig. 1a. The numerical results support the above expectation. For ζ 1, one finds λ (M) max = u M /(2ζ) with u M ≈ 1. 26. We note that this result satisfies the Perron-Frobenius bound at ζ 1, Therefore, the self-consistency equation (51) results in the following expression for the superconducting transition temperature, cf. Eq. (44): We note that the assumption ζ = 2/t c 1 is indeed satisfied.

13
The fact that ζ is large allows us to find the eigenvector ∆ n that corresponds to the maximum eigenvalue λ (M) max . Let us use the Euler-Maclaurin formula for the summation over n in Eq. (51). This is justified by the condition ζ 1. In addition, we approximate ln(n + n + 1) by ln(max{n, n } + 1) (as we shall see below, such simplification provides qualitatively correct results). With these approximations, we arrive at the following equation: where The quantity u ∞ = |γ 0 |/t 0 1 corresponds to u n with the maximally allowed index n 1/(2πT c τ). Let us parametrize ∆(u n ) as ∆(u n ) = ∆(u 0 ) f (u n ), where a function f (u) satisfies the normalization condition f (u 0 ) = 1. Then the function f (u) satisfies the following differential equation, where C = 2t 2 0 /|γ 0 | 1. Using the smallness of u ∞ , the solution of the above equation can be written as Here, J 2 (x) denotes the Bessel function of the first kind. The yet unknown parameter u 0 can be found from the boundary condition at u = u 0 , f (u 0 ) = Cu 0 . To the lowest order in C 1, the parameter u 0 can be found as where u c ≈ 0.92 is the solution of the equation F (u c ) = 0. Hence, neglecting the difference between u 0 and u c (proportional to t 2 0 /|γ 0 | 1), we find that the superconducting transition temperature is given by [cf. Eqs. (44) and (53)] We note that the 25% discrepancy in the numerical values of u c and u M indicates that the approximation ln(n + n + 1) by ln(max{n, n } + 1) is not justified by any small parameter. Nevertheless, our approximate solution for the eigenvector with the maximum eigenvalue, Eq. (56), is in good agreement with the eigenvector determined numerically, as shown in Fig. 1b. The estimate (44) for the superconducting temperature from the analysis of the renormalization group equations in the normal phase and the estimates (53) and (58) on the basis of the self-consistency equation are essentially the same except for different numerical values of the constant 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. 14 6. The energy dependence of the gap function

The gap function near T c
The non-trivial energy dependence of the gap function at the transition temperature suggests that ∆(z 0 ) vanishes at T = T c . The dependence of the gap function on temperatures at T c −T T c can be found in a way similar to the BCS theory. We expand the modified self-consistency equation (37) to the third order in the gap function: We note that the expansion of γ L √ in powers of ∆ ε and ∆ ε leads to the terms which are proportional to a small factor t 4 0 /γ 2 0 1. Such terms can be safely neglected. Let us parameterize the gap function as ∆ ε = ∆ 0 (T ) f (u ε ) where u ε = |γ 0 |t(L ε )/t 2 0 . Then taking into account that the sum over ε in the last term of the right-hand side of Eq. (59) is dominated by ε ∼ T c , we find where u T = |γ 0 |t(L πT )/t 2 0 . The above integral equation is reduced to Eq. (55) but with u T instead of u 0 . Using Eq. (57) and expressing u T and u 0 in terms of T and T c , respectively, we obtain the following well-known result of the BCS theory: We note that the corrections to the BCS-type temperature dependence are controlled by the small parameter t 2 0 /|γ 0 | 1.

The gap function at low temperatures T
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 introduce the energy ε 0 which is given by the solution of the equation ε 0 = ∆ ε 0 . Then, we assume that for ε < ε 0 the gap function ∆ ε is close to its value ∆ 0 at ε = 0. For ε > ε 0 the gap function ∆ ε is a monotonously decreasing function that reaches the value ∆ at ε ∼ 1/τ. At ε < ε 0 the self-consistency equation (37) can be approximately written as 15 Since the dependence of γ on L, governed by Eq. (50), is only logarithmical, the integral over ε in Eq. (62) is dominated by ε ∼ ε 0 . We then find: This result implies that ∆ ε is constant at ε < ε 0 up to corrections of the order of t 4 0 /γ 2 0 . With the same accuracy the energy scale ε 0 coincides with ∆ 0 . Further, we can substitute u ∆ 0 for u ∆ 0 + √ 2ε 0 in Eq. (63). Therefore, at ε < ∆ 0 the gap function behaves as At ε ε 0 ∆ 0 we approximate the self-consistency equation (37) as follows: Here, the dependence of γ on L is governed by Eq. (50). We note that Eq. (65) is justified for ε ε 0 . For ε ε 0 , a more accurate treatment of the self-consistency equation (37) results in corrections of the order of t 2 0 /|γ 0 |. Parametrizing the energy dependence of the gap function at ε ε 0 as ∆ ε = ε 0 f (u ε ), we find that the integral equation (65) results in Eq. (55) with u ε 0 entering in place of u 0 and with the value of the constant C = 2c 1 t 2 0 /|γ 0 |, where c 1 = arcsinh(1). We note that, in order to find the precise value of c 1 , one needs to take into account the t 2 0 /|γ 0 | 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 0 . With the same accuracy, we have ε 0 = ∆ 0 T c , cf. Eq. (58).
Expressing the right-hand side of Eq. (65) in terms of the function f (z) and setting ε ∼ 1/τ, we retrieve We note that ∆ ∆ 0 . All in all, we can summarize the energy dependence of the spectral gap at T = 0 as The overall dependence of the gap function ∆ ε on ε is shown in Fig. 2. We emphasize that within our approximate treatment, the values of ∆ 0 and T c coincide with the exponential accuracy, i.e. the exponential factor in Eq. (67) coincides with the exponential factor in Eq. (58). However, we cannot exclude a possibility that the magnitude of the ratio ∆ 0 /T c differs from the one in the BCS theory. This can result from the corrections of the order of t 2 0 /|γ 0 | that determine the coefficient C in Eq. In order to consider the effect of nonzero but low temperatures, T T c , we need to perform summation over Matsubara frequencies in the interval 0 < ε n ∆ 0 more accurately. We assume that the parameter ∆ 0 becomes temperature dependent, ∆ 0 → ∆ 0 (T ). Then Eq. (65) should be modified as follows: Using the parametrization, ∆ ε = ∆ 0 (T ) f (u ε ), one can check that the function f (u) satisfies Eq. (55) with u 0 replaced by u 0 (T ) = |γ 0 |t(L ∆ 0 (T ) )/t 2 0 . The constant C becomes a T -dependent function given by Using Eq. (57), we can express u 0 (T ) in terms of u 0 , yielding Therefore, the dependence of the magnitude of the gap function on temperature T T c at ε < ∆ 0 is the same as in the BCS theory. We note that there are corrections to this result which are of the order of t 2 0 /|γ 0 |.

Disorder-averaged density of states
In this section, we analyze the disorder-averaged density of states of single-particle excitations in the superconducting state. Within the NLSM formalism, the average density of states can be obtained as where the analytic continuation to real energies, ε n → E + i0, is performed. Here symbol 'sp' denotes the trace over particle-hole and spin spaces. Plugging in the parametrization (23) into Eq. (71), one finds The factor Z ε can be written as Projecting Eq. (73) to the BCS line, one finds that This factor encodeds the interaction-induced corrections to the density of states of the type that leads to the zero-bias anomaly in the normal state [39].
As one can see, in order to compute the average density of states, we need to make an analytic continuation for the gap function ∆ ε from Matsubara energies to real energies. We emphasize that 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.

The real and imaginary parts of ∆(E)
We restrict our consideration to the case of zero temperature, T = 0. Performing analytic continuation in Eq. (64), we find that the gap function at E ∆ 0 is purely real, At energies E ∆ 0 , the imaginary part of ∆(E) appears. To show that the nonzero imaginary part does exist, one can consider the perturbative expression (32). Performing analytic continuation to the real energy, we find that Im ∆(E ∼ 1/τ) ∼ ∆t 0 . 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 ε ∆ 0 , we energy dependence of the gap function is described by Eq. (67). After analytic continuation, iε → E + i0, in Eq. (67), we obtain Here, we have used the fact that t(L E ) 1. The real and imaginary parts of ∆(E) read: Expression (78a) has the same range of validity as Eq. (67). It provides us with an estimate for the real part of the gap function at E ∆ 0 withcorrections of the order of t 2 0 /|γ 0 | neglected. The situation with the imaginary part of the gap function is more delicate. Equation (78b) means that Im ∆(E ∼ ∆ 0 ) ∼ ∆ 0 t 4 0 /γ 2 0 . We note that the contribution of the same order will be given by the continuation of Eq. (75) to energies larger than ∆ 0 . This fact implies that, in order to compute Im ∆(E) at energies E ∼ ∆ 0 , one needs to retain terms of the order of t 4 0 /γ 2 0 in the solution of Eq. (37). However, at large energies, E ∆ 0 , Eq. (78b) provides the leading result to the imaginary part of the gap function. For example, at E ∼ 1/τ, Eq. (78b) yields Im ∆(E ∼ 1/τ) ∼ ∆t 0 that matches with the perturbative result. Additional argument that the more accurate (to the order t 4 0 /γ 2 0 ) expression for the real part of ∆(E) is needed in order to determine Im ∆(E) is given by the 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) = 0 for E E g and Im ∆(E) > 0 for E > E g . Making analytic continuation to the real frequencies in the selfconsistency equation (37), we obtain For energies close to the energy E g , we can expand the right hand side of Eq. (79) as (80b) We shall demonstrate below that the parameter α ∼ t 4 0 /γ 2 0 1. Solving Eq. (80a) for α 1, we find the dependence of the gap function for real energies close to E g , In agreement with our assumptions, Im ∆(E) = 0 for E E g . We note that ∆(E) E for E E g .
For E > E g the imaginary part of the gap function is non-zero, Im ∆(E) = −iα E 2 − E 2 g . We note that negative sign of Im ∆(E) at energies E > E g is needed for positivity of the density of states (see below). However, away from E g the imaginary part of Im ∆(E) has to change sign in order to be consistent with the asymptotic at large energies, Eq. (78b). The change of the sign of Im ∆(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) implies that E g is of the order of ∆ 0 . The sum over ε in Eq. (80b) is dominated by ε ∆ ε ∆ 0 . Then we find Here, we have taken into account that ∆ 2 (E g ) − E 2 g = αE g ∆ 0 .

The average density of states
With the help of the above results for ∆(E), we can derive the expression for disorderaveraged density of states as a function of energy. For large energies, E ∆ 0 , the imaginary part of ∆(E) can be neglected in comparison with its real part, see Eqs. (78a) and (78b). Since for E ∆ 0 the energy is always larger than the real part of ∆(E), we find Here and hereinafter we neglect the factor Z ε , see Eq. (72), since it provides negligible corrections of the order of t 2 0 /|γ 0 | in the parametric regime that we consider. In order to compute ρ(E) at energies close to E g ∆ 0 , we use the expression (81) for ∆(E).
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: The square-root growth of ρ(E) turns into the maximum at E = E max = E g (1 + α 2 /2). The magnitude at the maximum is given as ρ(E max ) = ν/(2α) ∼ νγ 2 0 /t 4 0 ν. For energies larger than E max , the density of states decays, 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. (See Sec. 8.2 for discussion of the mesoscopic fluctuations of the local density of states.) We note that the profile of ρ(E) is very similar to the density of states in the presence of the depairing term in the Usadel equation, introduced, e.g., by the scattering off magnetic impurities [54]. The comparison of Eq. (85), as well as the position of the maximum and its magnitude, with the Abrikosov-Gor'kov theory suggests that the effective depairing parameter η ∼ α 3 . We note that the depairing term in the Usadel equation appears also in the model with a spatially varying, random interaction in the Cooper channel [55] and with a spatially varying random order parameter [56]. 20 Figure 3: Sketch of the dependence of the local density of states ρ(E, r)/ν on the energy E for |γ 0 | t 0 |γ 0 |. The red solid curve corresponds to the disorder-averaged density of states (see text in Sec. 7.2). The dashed red curves and the shaded region illustrates the mesoscopic fluctuations of the local density of states. There are no quasiparticle states below E g . The energy E g * separates the region of strong and weak mesoscopic fluctuations. The energy E * demarcates the energy regions where the mesoscopic fluctuations are cut off by the length scale L * (E g < E < E * ) and the dephasing length 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).
By means of Eqs. (34) and (36), we find (for arbitrary signs of ε and ε ) 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: At energies below the gap edge, E E g , the gap function is real and satisfy ∆(E) E. Therefore, the fluctuations of the local density of states are zero identically, Above the superconducting gap, E > E g the mesoscopic fluctuations are non-zero. At energies close to the gap, 0 < E − E g E g , we find [δρ(E, r)] 2 2t 0 ρ 2 (E) ln min{L g , L} , 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 Finally, at large energies, E ∆ 0 , the general expression (92) can be simplified as We have neglected Im ∆(E) everywhere except the denominator of the first diffusive propagator on the right-hand side of Eq. (92). Here, the length L * is defined through Since Im ∆(E) is small, the length scale L * is large, L * L 2∆ 0 . In particular, at E ∼ ∆ 0 one can find the following estimate: L * ∼ L 2∆ 0 |γ 0 |/t 2 0 ∼ L g . 22 We note that the diffusive propagators in Eq. (92) are affected by dephasing due to electronelectron interactions [39]. This results in appearance of the dephasing length L φ in Eq. (96). For energies E ∆ 0 the dephasing length can be estimated as 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 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 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 0 ln(min{L * , L g , L φ }/ ) ∼ 1 for energies of the order of ∆ 0 . This indicates that the perturbative treatment of mesoscopic fluctuations has to be extended to a more elaborated renormalization group approach.

Renormalization group approach
The perturbative result (99) has exactly the same form as the perturbative correction to the dispersion of the local density of states above the superconducting transition. Essentially, the superconducting state leads to the specific infrared cut off length scale, min{L * , L g , L φ }, only. This suggests that one can use the one-loop renormalization group equation for m 2 ≡ [ρ(E, r)] 2 /ρ 2 (E) derived in Ref. [39], upto the length scale min{L * , L g , L φ }. Solving Eqs. (100) and (49), we find We note that the length scale L * starts from the value of the order of L g and, then, it grows with increasing E. The dephasing length L φ has the opposite behavior: it decreases with increasing energy. Therefore, there is an energy E * at which the lengths L * and L φ become of the same order. Using where u E = |γ 0 |t(L E )/t 2 0 , we can find the following estimate, 23 We note that E * ∆ 0 . 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 of the local density of states (in the case L = ∞), Interestingly, in accordance with Eq. (104), the disorder-averaged normalized second moment of the local density of states has the maximum at E = E * . The magnitude of this maximum can be estimated as Since the right hand side of the relation (105) is much smaller than one, we can approximate the disorder-averaged second moment of the local density of states at E g < E < E * as 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 local density of states decreases. Using Eq. (104), we find Thus, the mesoscopic fluctuations of ρ(E, r) are suppressed at the ultraviolet energy scale, E ∼ 1/τ. The dependence of [ρ(E, r)] 2 on the energy is illustrated in Fig. 3, where the functions ρ(E)± [ρ(E, r)] 2 are shown by dashed red curves. The fluctuations (not normalized) decreases with increase of the energy. They are strong in the energy interval, E g < E < E * . At E * the dependence of [ρ(E, r)] 2 on the energy is changed, as it was explained above. We note that there is a certain energy E g * which is the solution of the following equation: [ρ(E, r)] 2 = ρ(E). Using Eq. (104), we find that 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 local density of states below its average value.

Distribution function for the local density of states
In a similar way as it was done for the temperatures above T c [39], one can generalize the result (104) to the higher moments of the local density of states, We note that according to Eq. (104), the function A(E) > 1 in the energy interval E g < E < E g * , i.e. for region of strong mesoscopic fluctuations. The expression (109) implies the following 24

Summary and conclusions
To summarize, we have developed the theory of the multifractal superconducting state in thin films. Treating the fluctuations around the mean-field spatially homogeneous solution, we derived the modified Usadel equation that incorporates the interplay of disorder and interactions at 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, yields parametrically the same estimate for the multifractally enhanced superconducting transition temperature as the one derived by considering the instability in renormalization group equations in the normal phase.
(ii) The mutual effects of disorder and interactions result in strong dependence of the superconducting gap function on energy (see Fig. 2): at energies of the order of the spectral gap the gap function ∆ ε is parametrically enhanced in comparison with its magnitude at ultraviolet energies ∼ 1/τ.
(iii) The spectral gap at zero temperature is multifractally enhanced in the same way as the superconducting transition temperature.
(iv) The energy dependence of the gap function in the Usadel equation results in the profile of the disorder-averaged density of states that, near the spectral gap, resembles the one derived 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 effects. On the one hand, it results in the enhancement of the spectral gap, but on the other hand, it induces the effective depairing parameter that cuts off the coherence peaks in the average density of states. The corresponding depairing parameter is estimated as ∼ (t 2 0 /γ 0 ) 6 1.
(v) The mesoscopic fluctuations of the local density of states in the superconducting state are strong at 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 0 /|γ 0 | 1 that is similar to the estimate in the normal phase at temperatures close to the superconducting transition temperature [39]. It is worth emphasizing that strong spatial fluctuations of the local density of states (including the local value of the gap and the amplitude of coherence peaks) emerge in our theory despite the fact the model does not involve any macroscopic inhomogeneities. Indeed, our starting point is a model with short-range disorder and with all parameters being spatially uniform. Emergent strong fluctuations are a mesoscopic effect resulting from quantum interference in a disordered system.
Based on these findings for the statistics of fluctuations of the local density of states in thin films with multifractally-enhanced superconductivity, we conclude that disorder-induced interference effects dramatically affect spectral properties of these superconducting films, by fully governing the physics in a wide energy interval where the spectral gap for single-particle excitations 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 and non-zero in the other. The distribution function of the local density of states has a log-normal form, such that the typical value of the density of states is lower that the average value. In other words, 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 qualitative agreement with tunneling spectroscopy data on thin superconducting films [31][32][33][34][35][36] and with numerical solution of disordered attractive two-dimensional Hubbard model [37,38]. The strong mesoscopic fluctuations of the local density of states are accompanied by strong mesoscopic fluctuations of the superconducting order parameter.
On the technical side, within our approach, we have integrated over the spatial fluctuations of the superconducting order parameter from the very beginning. Therefore, they are hidden in the termŜ (c) int that describes the interaction in the Cooper channel. The renormalization group analysis of the NLSM in the normal phase suggests that such a procedure is more convenient, since it allows one to take into account the interplay between the spatial fluctuations of the superconducting 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 the NLSM action at length scales of the order of L T c . In our approach, the mesoscopic fluctuations of the local superconducting order parameter manifests themselves in the course of the renormalization group flow as the energy dependence of the gap function.
Since we considered a macroscopically homogeneous sample, it was natural to address fluctuations around the spatially homogeneous mean-field solution. The emergent fluctuations of the order parameter encoded in the mesoscopic fluctuations of the density of states, as observed in experiments of the "gap tomography", are due to mesoscopic quantum interference effects above the gap-edge energy. At the same time, it is interesting to see whether nonperturbative tails of the density of states inside the gap would appear in spatially homogeneous disordered films studied within the NLSM formalism. For this purpose, one should consider a possibility of existence of non-trivial saddle-point solutions of the NLSM action. A more straightforward way to obtain such tails would be through the introduction of macrocroscopic inhomogeneities (say, in the local concentration of impurities) that could induce "background" fluctuations of the superconducting gap.
We restricted our consideration to the case of a weak short-ranged interaction (suppressed by, e.g., high dielectric constant of the substrate), for which a strong effect of the enhancement of superconductivity by multifractality was predicted. The approach developed in this paper can be generalized in several directions. Our theory can be extended to include a strong shortranged interaction, as well as the Coulomb interaction. In particular, it would be interesting to study how the multifractality affects the BCS -BEC crossover [61,62]. One can also study the multifractal superconducting state that occurs in the system that without superconducting instability is near the interacting metal-insulator transition. Our theory can also be extended to the case of an applied magnetic field that destroys the superconducting state and may produce an intermediate phase with giant magnetoresistance [63] (for studying the giant magnetoresistance near the transition within the NLSM formalism, see Ref. [26]). Finally, our approach does not take into account phase fluctuations of the order parameter, giving rise to the Berezinskii-Kosterlitz-Thouless phenomena in superconducting 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. [46]. We note, however, that for sufficiently high conductance of the normal state, such effects are expected to influence the results described in this article only slightly.

Acknowledgements
The authors are grateful to M. Skvortsov  In this Appendix we estimate the mesoscopic fluctuations of the superconducting order parameter. While the quantity that is actually measured in experiments on the tomography of the superconductors is the local density of states (Sec. 8), the results of the scanning-tunnelingmicroscopy measurements are frequently translated into the maps of the fluctuating local order parameter by fitting the density of states with the BCS "ansatz". Here, instead of translating our results for the mesoscopic fluctuation of the density of states into the fluctuations of the order parameter, we calculate the latter directly within the NLSM formalism. The relation (9) suggests that the variance of the order parameter can be written as follows: Tr t r0 L α 1 0 Q(r) · Tr t r0 L α 2 0 Q(r) − 2 Tr t r0 L α 1 0 Q(r)t r0 L α 2 0 Q(r) , (A.1) where α 1 α 2 are some fixed replica indices. We note that the operator in the brackets in Eq. that follow from Eq. (25). Here, Y = (t/2) ln(L/ ) and we have neglected the energy dependence of the diffusive propagators at scales L L T c . Using parametrization (23), expanding the right hand side of Eq. (A.1) to the second order in W, and applying Eq. (25), we obtain the following perturbative result: Here we have used the self-consistency equation (31). We note that the fluctuations δ∆ do not lead to a finite single-particle density of states below the gap edge E g but rather correspond to the spatial fluctuations of the superconducting condensate. The hard spectral gap for the quasiparticles results from an effective averaging over such fluctuations with the self-consistency or renormalization group procedure.With the logarithmic accuracy, we can estimate the righthand side of Eq.  4) in a similar manner, to obtain yet stronger fluctuations of the superconducting order parameter in thin films with multifractally-enhanced superconductivity [60].

Appendix B. Kramers-Kronig relations for ∆(E)
In this appendix, we demonstrate how the Kramers-Kronig relations for ∆(E), We note that here we have neglected the contribution to Im ∆(E) from energies close to the spectral gap E g . As one can see, when the approximation (B.3), which is valid for E ∆ 0 , is used for all energies, an additional contribution to the real part of ∆(E) appears, which is of the order of (t 4 0 /γ 2 0 )z 2 (z 2 f (z)) . The latter is of the order of t 4 0 /γ 2 0 at E ∼ ∆ 0 . This implies that, in order to be able to compute the imaginary part of ∆(E) from the Kramers-Kronig relation, one needs to derive expression for Re ∆(E) with the accuracy of the order of t 4 0 /γ 2 0 .