Magnetic disorder in superconductors: Enhancement by mesoscopic fluctuations

TL;DR

This paper investigates how mesoscopic fluctuations in disordered superconducting films with magnetic impurities lead to a non-zero density of states within the superconducting gap, affecting the transition temperature and gap filling.

Contribution

It introduces a model accounting for mesoscopic fluctuations that enhances understanding of magnetic disorder effects beyond mean-field theory in superconductors.

Findings

Mesoscopic fluctuations induce a non-zero DOS within the superconducting gap.

The effect can surpass potential disorder fluctuations in filling the gap.

The suppression of $T_c$ varies depending on magnetic impurity strength.

Abstract

We study the density of states (DOS) and the transition temperature $T_c$ in a dirty superconducting film with rare classical magnetic impurities of an arbitrary strength described by the Poissonian statistics. We take into account that the potential disorder is a source for mesoscopic fluctuations of the local DOS, and, consequently, for the effective strength of magnetic impurities. We find that these mesoscopic fluctuations result in a non-zero DOS for all energies in the region of the phase diagram where without this effect the DOS is zero within the standard mean-field theory. This mechanism can be more efficient in filling the mean-field superconducting gap than rare fluctuations of the potential disorder (instantons). Depending on the magnetic impurity strength, the suppression of $T_c$ by spin-flip scattering can be faster or slower than in the standard mean-field theory.

Figures

Equations

$$H_{\mathrm{mag}}=J\sum _j{\psi }^{†}({\mathbf{r}}_j){\mathbf{S}}_j\mathbf{\sigma }\psi ({\mathbf{r}}_j).$$

$$\mathcal{S}={\mathcal{S}}_D+{\mathcal{S}}_{\Delta }+{\mathcal{S}}_{\text{mag}}.$$

$$\mathcal{S}_{D}=\frac{\pi\nu}{8}\int d^{2}\bm{r}\operatorname{tr}\bigl{[}D(\nabla Q)^{2}-4(\varepsilon\tau_{3}+\Delta\tau_{1})Q\bigr{]},$$

$$Q^2=1,Q=\overline{Q}\equiv {\tau }_1{\sigma }_2{Q}^{\text{T}}{\tau }_1{\sigma }_2.$$

$${\mathcal{S}}_{\Delta }=\frac{\nu }{\lambda T}\int d^2\mathbf{r}\sum _{a=1}^N|{\Delta }^a(\mathbf{r}){|}^2.$$

$${\mathcal{S}}_{\text{mag}}\approx \sum _j{s}_{\text{mag}}^{(j)}=-\frac{1}{2}\sum _j\mathrm{tr}\ln \left(1+i\sqrt{\alpha }Q({\mathbf{r}}_j){\tau }_3\mathbf{\sigma }{\mathbf{n}}_j\right).$$

$$\Bigl{\langle}\exp\sum_{j}f(\bm{r_{j}})\Bigr{\rangle}=\exp\left\{n_{s}\int d^{2}\bm{r}\Bigl{[}e^{f(\bm{r})}-1\Bigr{]}\right\},$$

$${\mathcal{S}}_{\text{mag}}\to -n_s\int d^2\mathbf{r}\left({⟨e^{\frac{1}{2}\mathrm{tr}\ln \left(1+i\sqrt{\alpha }Q(\mathbf{r}){\tau }_3\mathbf{\sigma }\mathbf{n}\right)}⟩}_{\mathbf{n}}-1\right).$$

$${\mathcal{S}}_{\text{mag}}=$$

$$+\frac{1}{2!}\sum _{m,n=1}^{\infty }\frac{(-1{)}^{m+n}}{4mn}\hat{C}{T}_{mn}$$

$$\displaystyle+\frac{1}{3!}\sum_{m,n,p=1}^{\infty}\frac{(-1)^{m+n+p-1}}{8mnp}\hat{C}T_{mnp}+\dots\Biggr{]}.$$

$$\hat{C}{T}_m=$$

$$\hat{C}{T}_{mn}=$$

$$\times \mathrm{tr}(Q{A}_{i_{m+1}}\dots Q{A}_{i_{m+n}}),$$

$${\mathcal{S}}_{\text{mag}}=$$

$$\displaystyle-\frac{1}{32}\hat{C}T_{22}+\frac{1}{32}\hat{C}T_{211}-\frac{1}{384}\hat{C}T_{1111}+\dots\Biggr{]}.$$

$$t(\epsilon )=\frac{1}{\pi g}\ln \frac{1}{|\epsilon |\tau },$$

$$\displaystyle\mathcal{S}_{\text{mag}}=n_{s}\int d^{2}\bm{r}\Biggl{[}\gamma_{2}\hat{C}T_{2}+\gamma_{11}\hat{C}T_{11}+\gamma_{4}\hat{C}T_{4}+\gamma_{31}\hat{C}T_{31}$$

$$\displaystyle+\gamma_{22}\hat{C}T_{22}+\gamma_{211}\hat{C}T_{211}+\gamma_{1111}\hat{C}T_{1111}+\dots\Biggr{]}.$$

$${\gamma }_{k_1{k}_2\dots k_q}(t)={\gamma }_{k_1{k}_2\dots k_q}(0)e^{n(n-1)t},$$

$${\mathcal{S}}_{\text{mag}}=-n_s\int d^2\mathbf{r}{⟨\exp \left\{\frac{1}{4}\mathrm{tr}\ln [1+\text{a}(Q{\tau }_3{)}^2]\right\}-1⟩}_{\text{a}},$$

$${\mathcal{P}}_{\alpha }(\text{a},t)=\frac{1}{4\text{a}\sqrt{\pi t}}\exp \left[-\frac{1}{4t}{\left(\frac{1}{2}\ln \frac{\text{a}}{\alpha }+t\right)}^2\right].$$

$$Q={\tau }_1\sin \theta +{\tau }_3\cos \theta .$$

$$\Delta =\pi \lambda T\sum _{\epsilon }\sin {\theta }_{\epsilon }.$$

$${\mathcal{S}}_{\text{mag}}^{\text{MF}}=-\frac{n_s}{4}\int d^2\mathbf{r}{⟨\mathrm{tr}\ln [1+\text{a}(Q{\tau }_3{)}^2]⟩}_{\text{a}}.$$

$$\mathcal{S}_{\text{mag}}^{\text{MF}}=-\frac{n_{s}}{2}\sum\limits_{\sigma=\pm}\int d^{2}\bm{r}\left\langle\ln\bigl{(}1+\textsf{a}\,e^{2i\sigma\theta}\bigl{)}\right\rangle_{\textsf{a}}.$$

$$\epsilon \sin {\theta }_{\epsilon }-\Delta \cos {\theta }_{\epsilon }+\frac{n_s}{\pi \nu }{⟨\frac{\text{a}\sin 2{\theta }_{\epsilon }}{1+{\text{a}}^2+2\text{a}\cos 2{\theta }_{\epsilon }}⟩}_{\text{a}}=0,$$

$${\theta }_{\epsilon }\approx \Delta /(\epsilon +1/{\tau }_s),$$

$$\frac{1}{{\tau }_s}=\frac{2n_s}{\pi \nu }{⟨\frac{\text{a}}{(1+\text{a}{)}^2}⟩}_{\text{a}}.$$

$$\frac{{\tau }_{s0}}{{\tau }_s}\approx 1+\frac{2-16\alpha +6{\alpha }^2}{(1+\alpha {)}^2}t(T_c)+O(t^2).$$

$$\ln \frac{T_{c0}}{T_c}=\psi \left(\frac{1}{2}+\frac{1}{2\pi T_c{\tau }_s}\right)-\psi \left(\frac{1}{2}\right),$$

$$F_E(\psi )=\sinh \psi -\frac{E}{\Delta }\cosh \psi -\frac{[\alpha n_s/(\pi \nu \Delta )]\sinh 2\psi }{1+{\alpha }^2-2\alpha \cosh 2\psi }.$$

$$\rho (E)/2\nu =\sqrt{2/3}{\eta }^{-2/3}\text{Re}\sqrt{ϵ+i{ϵ}_{\ast }},$$

Full text

I Introduction

The properties of superconductors in the presence of impurities have remained at the focus of intense theoretical and experimental research during the past half-century. It is generally accepted that the potential scattering in $s$-wave superconductors affects neither the transition temperature, $T_{c}$, nor the density of states (DOS), $\rho(E)$. This statement usually referred to as Anderson’s theorem Abrikosov and Gor’kov (1959a, b); Anderson (1959) is valid for sufficiently good metals. As the potential disorder increases, the emergent inhomogeneity due to the interplay of quantum interference (Anderson localization) and interaction leads to modification of $T_{c}$ Maekawa and Fukuyama (1982); Maekawa et al. (1984); Anderson et al. (1983); Finkel’stein (1987); Feigel’man et al. (2007, 2010); Feigel’man and Skvortsov (2012); Burmistrov et al. (2012, 2015) and $\rho(E)$ Hurault and Maki (1970); Abrahams et al. (1970); Di Castro et al. (1990); Burmistrov et al. (2016), with the effect being controlled by the parameter $1/(k_{F}l)\ll 1$ (where $k_{F}$ is the Fermi momentum and $l$ is the mean free path).

Magnetic impurities violating the time-reversal symmetry affect superconductivity much stronger, already at $k_{F}l\to\infty$. Classical magnetic impurities lead both to suppression of $T_{c}$ and to reduction of the superconducting gap in $\rho(E)$ with the increase of their concentration $n_{s}$ Abrikosov and Gor’kov (1960). Beyond the Born limit, magnetic impurities produce degenerate subgap bound states (see Fig. 1a). Their hybridization results in the formation of an energy band giving rise to a nontrivial DOS structure Yu (1965); Soda et al. (1967); Shiba (1968); Rusinov (1969). The account for the Kondo effect Müller-Hartmann and Zittartz (1971); Matsuura (1977); Bickers and Zwicknagl (1987), the indirect exchange interaction between magnetic impurities Ruvalds and Liu (1981), or the spin-flip scattering assisted by the electron-phonon interaction Jarrell (1990) can lead to the reentrant behavior of $T_{c}$ vs. $n_{s}$ (see Ref. Balatsky et al. (2006) for a review).

A hard gap in $\rho(E)$ obtained for superconductors with magnetic impurities in the mean-field approximation is smeared by inhomogeneity. This can be due to rare fluctuations of a potential disorder Lamacraft and Simons (2000, 2001); Meyer and Simons (2001); Marchetti and Simons (2002), $n_{s}$ Silva and Ioffe (2005), or superconducting order parameter Larkin and Ovchinnikov (1971). A combined theory of these mechanisms has been developed in Refs. Skvortsov and Feigel’man (2013); Fominov and Skvortsov (2016).

In this Letter we describe a novel mechanism for smearing of the superconducting gap. We reconsider the problem of rare classical magnetic impurities with the Poissonian statistics in a dirty superconductor. The key point that distinguishes our work from the previous ones is that we take into account mesoscopic fluctuations of the local DOS in a potential disorder. Physically, this implies that the energies of subgap bound states become dependent on the spatial positions of magnetic impurities (see Fig. 1b). Averaging over these bound states results in a non-zero homogenous DOS at all energies in the region of the phase diagram where in the absence of this effect $\rho(E)$ is zero within the mean-field approximation. Motivated by the recent experiment on magnetic Gd impurities in superconducting MoGe films Kim et al. (2012), in this Letter we develop the theory of the enhancement of magnetic disorder by mesoscopic fluctuations in the case of a dirty superconducting film.

The outline of the paper is as follows. In Sec. II we present description of dirty superconductors with rare magnetic impurities in terms of the nonlinear sigma model and its renormalization. Our results for the renormalized spin-flip rate, superconducting transition temperature, and the density of states are given in Sec. III. We end the paper with discussions (Sec. IV) and conclusions (Sec. V) Some details of calculations are presented in Appendices.

II Nonlinear sigma model for paramagnetic impurities

We consider a two-dimensional (2D) dirty $s$-wave superconductor in the presence of both potential (spin-preserving) and magnetic disorder. Scattering off the former is responsible for the dominant contribution to the momentum relaxation rate $1/\tau$. Much weaker spin-flip scattering rate is related with the exchange interaction between magnetic impurities and electrons described by the Hamiltonian

[TABLE]

We shall treat rare magnetic disorder under standard assumptions Yu (1965); Soda et al. (1967); Shiba (1968); Rusinov (1969); Marchetti and Simons (2002): (i) impurity positions $\bm{r}_{j}$ have the Poisson distribution; (ii) impurity spins $\bm{S}_{j}$ are classical statistically independent vectors with the flat distribution over their orientations, $\prod_{j}\delta(\bm{S}_{j}^{2}-S^{2})$.

The low energy description of two-dimensional disordered superconductors with rare paramagnetic impurities can be conveniently formulated in terms of the replicated verion of a nonlinear sigma-model Finkelstein (1990); Belitz and Kirkpatrick (1994). Its action can be written as

[TABLE]

Here $\mathcal{S}_{D}$ is the standard diffusive action

[TABLE]

where $\nu$ and $D$ denote the density of states at the Fermi energy (per one spin projection) and the diffusive coefficient in the normal state, respectively. The matrix $Q$ operates in the spin, Nambu, replica, and Matsubara energy spaces. It is subject to the following constraints Houzet and Skvortsov (2008):

[TABLE]

Here the transposition ${}^{\textsf{T}}$ acts in both the Matsubara energy space and the replica space. The Pauli matrices $\tau_{j}$ ($\sigma_{j}$) act in the Nambu (spin) spaces. The matrix $\varepsilon$ is the diagonal matrix with the elements $\varepsilon_{n}=\pi T(2n+1)$.

The superconducting correlations are described by the order-parameter matrix $\Delta$ which is diagonal in the Nambu space with matrix elements $\Delta^{a}(\bm{r})$. In the absence of a supercurrent, $\Delta$ is chosen to be real. The action $\mathcal{S}_{\Delta}$ reads

[TABLE]

Here $N$ stands for the number of replica and $\lambda>0$ denotes the attraction amplitude in the Cooper-channel.

We consider the case of rare classical magnetic impurities with the concentration $n_{s}$ [the precise condition on $n_{s}$ see below], when the magnetic part of the action, $\mathcal{S}_{\text{mag}}$, becomes separable in the individual magnetic impurities Meyer and Simons (2001):

[TABLE]

Here $\bm{n}_{j}$ stands for the three-dimensional unit vector and the dimensionless parameter $\alpha=(\pi\nu JS)^{2}$ is expressed in terms of the impurity spin $S$ and exchange constant $J$. We note that approximation (6) of the full action $\mathcal{S}_{\text{mag}}$ is equivalent to the self-consistent $T$-matrix approximation for magnetic scattering which treats all orders in scattering off a single magnetic impurity but neglects diagrams with intersecting impurity lines.

Performing the Poisson averaging over positions of the magnetic impurities with the help of the following relation Friedberg and Luttinger (1975)

[TABLE]

we find that the contribution to the nonlinear sigma model action due to magnetic impurities becomes

[TABLE]

Here $\langle\dots\rangle_{\bm{n}}$ stands for the averaging over direction of the unit vector $\bm{n}$. Expanding $\mathcal{S}_{\text{mag}}$ in powers of $\sqrt{\alpha}$, we find

[TABLE]

Here we introduced the operators

[TABLE]

and so on. The operator $\hat{C}$ acts as the symmetric tensor: $C_{i_{1}\dots i_{m}}=\langle n_{i_{1}}\dots n_{i_{m}}\rangle$ . For convenience we defined the self-dual matrix $\bm{A}=i\sqrt{\alpha}\tau_{3}\bm{\sigma}=\overline{\bm{A}}$ . Since operators $T_{n\dots}$ are symmetric with respect to its indices, the expansion can be written in the following form:

[TABLE]

The nonlinear sigma model action (2) with the magnetic part given by Eq. (11) provides full description of quantum effects for a dirty superconductor in the diffusive regime. These effects (weak localization and Aronov-Altshuler-type corrections) are responsible for the renormalization of system’s parameters, e.g. the diffusion coefficient and the attraction amplitude. In the 2D case, the magnitude of quantum corrections at the energy scale $\varepsilon$ is governed by the parameter

[TABLE]

where $g=h/(e^{2}R_{\square})\gg 1$ is the bare dimensionless conductance of the film. In a superconductor, renormalization stops at $\varepsilon\sim\max\{T_{c},|\Delta|\}\sim T_{c}$. Assuming that the transition temperature is not too low, $t(T_{c})\ll 1$, one can neglect the renormalization of the conductance and interaction parameters between the energy scales $1/\tau$ and $T_{c}$ (see Refs. Finkelstein (1990); Belitz and Kirkpatrick (1994) for a review). In contrast, renormalization of the magnetic-impurity part $S_{\rm mag}$ of the nonlinear sigma model is essential.

Treating this renormalization in the one-loop approximation, we find that after the renormalization this part of the action can be written as (see Appendix A)

[TABLE]

Here the coefficients $\gamma_{k_{1}k_{2}\dots k_{q}}$, where $k_{1}+k_{2}+\dots+k_{q}=n$, are given as follows:

[TABLE]

where initial values of the coefficients $\gamma_{k_{1}k_{2}\dots k_{q}}(0)$ follow from Eq. (11).

In what follows we are interested in the singlet sector of the theory. Therefore, one can operate with $Q$ matrix which is the unit matrix in the spin space, $Q=Q_{0}\sigma_{0}$. Then we can average over directions of the impurity magnetization $\bm{n}$ in operators $\hat{C}T_{k_{1}k_{2}\dots k_{q}}$. Then the renormalized magnetic-impurity part of the nonlinear sigma-model action (13) can be written in the following convenient short-hand notation (see Appendix B):

[TABLE]

where the averaging $\langle\dots\rangle_{\textsf{a}}$ is defined with respect to the following log-normal distribution function:

[TABLE]

Comparing Eq. (15) with Eq. (8) we may interpret the effect of renormalization as follows: Now instead of a single value of $\alpha$ there is a log-normal distribution of the effective strength of impurity a, schematically shown in Fig. 1.

III Results

In the mean-field approximation, a dirty superconductor in the diffusive regime is described by two coupled equations: the self-consistency equation for the superconducting order parameter $\Delta$ and the Usadel equation for the quasiclassical Green’s function Usadel (1970); Belzig et al. (1999); Kopnin (2001). These equations can be derived as the saddle-point equations of the nonlinear sigma model described in the previous section. The mean-field solution for the $Q$ matrix can be parametrized as

[TABLE]

Performing variation of the action $\mathcal{S}$ on the configuration (17) with respect to $\Delta$, we find the following mean-field self-consistency equation:

[TABLE]

For the study of space-averaged configurations at the mean-field level, it is sufficient to retain the term of the first order in trace only in the renormalized action (15):

[TABLE]

Since the eigenvalues of $(Q\tau_{3})^{2}$ are $e^{\pm 2i\theta}$, we find explicitly

[TABLE]

Performing variation of the action $\mathcal{S}$ [with $\mathcal{S}_{\text{mag}}$ given by Eq. (20)] on the configuration (17) with respect to $\theta_{\varepsilon}$, we find the following modified Usadel equation:

[TABLE]

where $\theta_{\varepsilon}$ is the energy-dependent spectral angle, $\nu$ is the normal DOS at the Fermi energy per one spin projection, and $\varepsilon=\pi T(2n+1)$ denotes the fermionic Matsubara frequencies. We remind that the averaging $\langle\dots\rangle_{\textsf{a}}$ in Eq. (21) is defined with respect to the log-normal distribution function (16). Since $\mathcal{P}_{\alpha}(\textsf{a},t\to 0)\to\delta(\textsf{a}-\alpha)$, Eq. (21) at $t=0$ coincides with the standard Usadel equation in the case of magnetic impurities Abrikosov and Gor’kov (1960); Shiba (1968); Rusinov (1969); Marchetti and Simons (2002). The linearity of Eq. (21) in $n_{s}$ is justified for small concentration of magnetic impurities: $n_{s}\xi^{2}/g\ll 1$, where $\xi=l/\sqrt{T_{c}\tau}$ is the dirty superconducting coherence length.

The quantity a in Eq. (21) plays a role of the renormalized impurity strength. Since the bare impurity strength $\alpha$ is proportional to the local DOS which is subjected to mesoscopic fluctuations, $\mathcal{P}_{\alpha}(\textsf{a},t)$ reflects the log-normal distribution of the local DOS in 2D weakly disordered systems Altshuler et al. (1986); Lerner (1988). Contrary to naive expectations, one should average over a the Usadel equation rather than physical observables, e.g. the DOS. This is a consequence of the Poisson distribution of impurity positions $\bm{r}_{j}$.

III.1 Effective spin-flip rate

In the vicinity of the thermal transition, $\Delta\to 0$ and we can linearize Eq. (21) with respect to $\theta_{\varepsilon}$. This procedure yields

[TABLE]

where the effective spin-flip rate is given by

[TABLE]

At $t=0$, one recovers the standard expression for the bare spin-flip rate due to magnetic impurities, $1/\tau_{s0}=2\alpha n_{s}/[\pi\nu(1+\alpha)^{2}]$ Rusinov (1969). In the limiting cases $\alpha\to 0$ and $\alpha\to\infty$, the spin-flip rate (23) becomes enhanced in comparison with the bare one: $1/\tau_{s}=\exp(2t)/\tau_{s0}$ and $1/\tau_{s}=\exp(6t)/\tau_{s0}$, respectively. For an arbitrary value of $\alpha$, the asymptotic expansion at $t\ll 1$ has the form (see Appendix C.1):

[TABLE]

At small $t$, the spin-flip rate is suppressed (enhanced) for $\alpha_{0}<\alpha<1/(3\alpha_{0})$ (otherwise), where $\alpha_{0}=1/(4+\sqrt{13})\approx 0.13$. The overall behavior of the ratio $\tau_{s0}/\tau_{s}$ as a function of $t$ and $\alpha$ is illustrated in Fig. 2.

III.2 Transition temperature

Since the spin-flip rate (23) is the only characteristic of magnetic disorder that enters the linearized solution for $\theta_{\varepsilon}$, we obtain the standard equation for the superconducting transition temperature (see Appendix C.2):

[TABLE]

where $T_{c0}$ denotes the transition temperature in the absence of magnetic impurities, and $\psi(z)$ stands for the di-gamma function. Equation (25) was derived by Abrikosov and Gor’kov (AG) in the Born limit ($\alpha\to 0$) Abrikosov and Gor’kov (1960), and later was shown to describe the suppression of $T_{c}$ for arbitrary values of $\alpha$ Rusinov (1969). For a scale-independent spin-flip time, $\tau_{s}=\tau_{s0}$, Eq. (25) defines a universal function $T_{c}^{\text{AG}}(1/\tau_{s0})$ shown by the black dashed line in Fig. 3. Superconductivity is eventually destroyed at the critical spin-flip rate $1/\tau_{sc}=2\pi e^{\psi(1/2)}T_{c0}\approx 0.88\,T_{c0}$ Abrikosov and Gor’kov (1960). This standard approach corresponds to the limit $t=0$, when mesoscopic fluctuations can be neglected.

An essential modification introduced by the log-normal distribution of the impurity strength (16) is that now the spin-flip rate $1/\tau_{s}$ depends on the parameter $t(T_{c})$, i.e. on the conductance $g$ and the transition temperature $T_{c}$ itself. This leads to an unusual behavior illustrated in Fig. 3, where we present the numerical solutions of Eq. (25) for fixed values of $g$ and $T_{c0}\tau$ and for various values of $\alpha$. At finite $t$, dependence of $1/\tau_{s}$ on $T_{c}$ renders the curves $T_{c}(1/\tau_{s0})$ sensitive to a particular value of $\alpha$. In the range $\alpha_{0}<\alpha<1/(3\alpha_{0})$, the spin-flip rate decreases monotonously down to zero with increasing $t$. Therefore the reduction of $T_{c}$ with the increase of $1/\tau_{s0}$ is slower than for $t=0$. This agrees qualitatively with the slowdown of $T_{c}$ suppression with increasing the film resistance measured in Ref. Kim et al. (2012). In the opposite case, for $\alpha<\alpha_{0}$ and $\alpha>1/(3\alpha_{0})$, the dependence of $T_{c}$ on $1/\tau_{s0}$ is qualitatively different since the ratio $\tau_{s0}/\tau_{s}$ can be larger than unity and is a non-monotonous function of $t$. Since the spin-flip rate is enhanced, the reduction of $T_{c}$ with the increase of $1/\tau_{s0}$ is faster than in the case $t=0$. The non-monotonicity of $\tau_{s0}/\tau_{s}$ results in the existence of two solutions of Eq. (25) for $T_{c}$. Formally, it admits the solution with nonzero $T_{c}$ for any value of the parameter $\tau_{sc}/\tau_{s0}$. However, we remind that our approach is valid provided the inequality $T_{c}\gg\exp(-\pi g)/\tau$ holds.

The dependence of the spin-flip rate on $g$ transforms into the dependence of $T_{c}$ on the film conductance. To illustrate this effect, we fix the value of the parameter $\tau_{sc}/\tau_{s0}$ and plot the ratio $T_{c}/T_{c}^{\text{AG}}(1/\tau_{s0})$ on the film resistance $R_{\square}$ for some values of $\alpha$ in the inset to Fig. 3. Since for $\alpha_{0}<\alpha<1/(3\alpha_{0})$ the spin-flip rate decreases monotonously with the increase of $t$, $T_{c}$ is enhanced with respect to $T_{c}^{\text{AG}}$. The non-monotonous dependence of $1/\tau_{s}$ on $t$ obtained for $\alpha<\alpha_{0}$ and $\alpha>1/(3\alpha_{0})$ leads to the reentrant behavior of $T_{c}$ on $R_{\square}$.

It is worthwhile to mention that not only the suppression of $T_{c}$ by magnetic impurities but also the reduction of $\Delta$ is modified at finite $g$ due to the log-normal distribution of the effective impurity strength els .

III.3 Density of states

Consider now the superconducting phase with a finite $\Delta$. The DOS can be obtained from the solution of Eq. (21) after analytic continuation to real energies $E$: $\rho(E)=2\nu\operatorname{Re}\cos\theta_{-iE+0}$. It is convenient to parametrize the spectral angle as $\theta=\pi/2+i\psi$. Without renormalization ($t=0$), the angle $\psi(E)$ should be determined from equation $F_{E}(\psi)=0$, where Abrikosov and Gor’kov (1960); Shiba (1968); Rusinov (1969); Marchetti and Simons (2002)

[TABLE]

This leads to a complicated structure of the DOS at energies $|E|<\Delta$, which depends on the values of $\alpha$ and $\eta=1/\tau_{s0}\Delta$ (see Ref. Fominov and Skvortsov (2016) for a review). In the case $\eta>(\frac{1-\alpha}{1+\alpha})^{2}$, the impurity band touches the Fermi energy, leading to a finite DOS at $E=0$. Below we shall consider the opposite regime, $\eta<(\frac{1-\alpha}{1+\alpha})^{2}$, in which $\rho(E)$ has a finite gap $E_{g0}$ for $t=0$. The gap opens since $F_{E}(\psi)=0$ possesses only real solutions at energies $|E|<E_{g0}$.

Typical modification of the DOS at finite $t$ is illustrated in Fig. 4, where we plot $\rho(E)$ obtained by numerical solution of Eq. (21) at $\varepsilon\to-iE+0$. For $t\ll 1$, mesoscopic fluctuations of magnetic disorder affect the DOS in two ways. (i) On the perturbative level, they shift the position of the gap: $E_{g0}\to E_{g}$, with $(E_{g0}-E_{g})/E_{g0}\propto t$, but the gap remains hard. (ii) A finite DOS below the renormalized gap is then generated nonperturbatively in $t$, due to the tail of the distribution $\mathcal{P}_{\alpha}(\textsf{a},t)$. In Fig. 4, the smearing of $E_{g}$ can be clearly seen for $t=0.03$, whereas for larger $t$ the smearing and gap shift cannot be separated. A profound feature of the DOS is its finite value right at the Fermi energy.

In the limit of weak renormalization, $t\ll 1$, the DOS can be obtained analytically. The general expression is quite cumbersome (see Appendix C.3), so we present here only the results in the regime of weak magnetic impurities ($\alpha\ll\eta^{2/3}\ll 1$). The gap smearing at $E\to E_{g}$ is described by

[TABLE]

where $\epsilon=(E-E_{g})/\Delta$ and

[TABLE]

The subgap DOS (27) decays as power law. The residual DOS at the Fermi energy it is determined by the probability $\mathcal{P}_{\alpha}(1,t)$ to find $\textsf{a}=1$ and in the limit $t\ll 1$ reads

[TABLE]

This result is non-perturbative in both $t$ and $\alpha$. The dependence of $\rho(0)$ on $t$ for some values of $\ln\alpha$ is shown in the inset to Fig. 4. Its non-monotonicity is related to that of $\mathcal{P}_{\alpha}(1,t)$ as a function of $t$. At a fixed value of $t$, $\rho(0)$ behaves non-monotonically with the impurity strength $\alpha$ at a given value of $\eta$.

IV Discussions

Our main Eq. (21) could be derived for a toy model of Poissonian magnetic impurities with the strength independently distributed according to $\mathcal{P}_{\alpha}(\textsf{a},t)$. We emphasize however that in a disordered film the log-normal distribution is generated intrinsically due to mesoscopic fluctuations of the local DOS.

The log-normal distribution $\mathcal{P}_{\alpha}(\textsf{a},t)$ predicts an exponentially small probability for realization of very small and very large values of the effective impurity strength a. As well-known from the theory of mesoscopic fluctuations of the local DOS and wave function multifractality, this implies that typically the impurity strength $\textsf{a}_{-}<\textsf{a}<\textsf{a}_{+}$ is realized Mirlin (2000). Using results of Ref. Foster et al. (2009), we obtain the following estimate for the termination points: $\textsf{a}_{\pm}=\alpha\exp[\pm(4/\sqrt{\pi g})\ln 1/(T_{c}\tau)]$ (see Appendix D). In order our result (29) for $\rho(0)$ were applicable to a typical sample, vicinity of $\textsf{a}=1$ should be inside the interval $(\textsf{a}_{-},\textsf{a}_{+})$. It is fulfilled provided $(4/\sqrt{\pi g})\ln 1/(T_{c}\tau)\gg\ln(1/\alpha)$.

We emphasize the difference with the instanton analysis Marchetti and Simons (2002); Fominov and Skvortsov (2016), where the effect of mesoscopic fluctuations on magnetic disorder was not taken into account: (i) in our approach, the DOS is modified already at the mean-field level and (ii) our results (27)–(29) involve a spreading resistance $\ln(\xi/l)/2\pi g$ which is parametrically larger than the sheet resistance $1/g$ emerging in the instanton analysis. As a result, our mechanism predicts a larger DOS at the Fermi level, whereas near $E_{g}$ it prevails provided $\ln^{2}[\eta^{2/3}/(4\alpha)]\ll 16t\ln[g\eta^{1/2}/(\alpha^{3/4}\sqrt{t})]$ (see Appendix C.3).

Although our results were derived for a weak disorder, $t\ll 1$, they can be extended to the case of a moderate disorder, $t\sim 1$ (provided $g\gg 1$) els . In this situation the mean-field equation for $\theta_{\varepsilon}$ remains the same as Eq. (21), but the distribution function $\mathcal{P}_{\alpha}(\textsf{a},t)$ must be found taking Fermi-liquid renormalizations into account.

The enhancement of magnetic disorder due to mesoscopic fluctuations is not restricted to classical magnetic impurities. It is known Kettemann and Mucciolo (2006); Micklitz et al. (2006); Kettemann and Mucciolo (2007); Micklitz et al. (2007) that the Kondo effect in the disordered electron systems is also modified by mesoscopic fluctuations of the local DOS. Therefore, the theory for the interplay of the Kondo effect and superconductivity developed in Refs. Müller-Hartmann and Zittartz (1971); Matsuura (1977); Bickers and Zwicknagl (1987) needs to be modified for disordered films els .

The dependence of $T_{c}$ on the film conductance can be caused by a variety of reasons, among which are the dependence of the DOS at the Fermi energy on disorder, renormalization of the Cooper channel attraction in ballistic and diffusive regimes, Berezinskii-Kosterlitz-Thouless transition etc. Finkelstein (1990); Belitz and Kirkpatrick (1994) The sensitivity of the spin-flip rate on the conductance is a new mechanism providing a nontrivial dependence of $T_{c}$ on $g$.

V Conclusions

To summarize, we reconsidered the problem of rare classical magnetic impurities with the Poissonian statistics in a dirty superconducting film. We took into account renormalization of the multiple spin-flip scattering due to mesoscopic fluctuations of the local DOS in a potential disorder. This effect results in the log-normal distribution of the effective magnetic impurity strength rendering the energies of quasiparticle bound states position dependent (see Fig. 1). In the superconducting state, this results in the smearing of the hard gap (obtained in the absence of spin-flip renormalization) and emergence of a non-zero DOS for all energies already at the mean-field level. Depending on the bare magnetic impurity strength, the superconducting transition temperature is suppressed by the spin-flip scattering slower or faster than in the absence of renormalization. Finally, we mention that our results can be extended to the model with an arbitrary distribution of magnetic impurities, the vicinity of a superconductor-insulator transition, the case with Coulomb repulsion in addition to attraction, the presence of Zeeman splitting, etc. els .

Acknowledgements.

We thank M. V. Feigel’man, Ya. V. Fominov, and A. D. Mirlin for useful discussions. The research was supported by Skoltech NGP Program (Skoltech–MIT joint project).

Appendix A Renormalization of the action $\mathcal{S}_{\text{mag}}$

In this Appendix we present details of the one-loop renormalization of the magnetic-impurity part of the action.

To renormalize $\mathcal{S}_{\text{mag}}$ we write $Q=\Lambda(1+W+\dots)$, where $W$ obeys two linear constraints: $\Lambda W+W\Lambda=0$ and $W=-\overline{W}$, and the convergency condition $W=-W^{\dagger}$. The matrix $\Lambda$ is assumed to be self-dual: $\Lambda=\overline{\Lambda}$. Then the quadratic part of the action reads:

[TABLE]

The quadratic part of the action determines the following contraction rules:

[TABLE]

where $t=[2/(\pi g)]\ln(L/l)$ with $L$ denoting the infrared length scale.

Next we write the matrix $Q$ as $Q=U^{-1}\Lambda(1+W+\dots)U$ where the slow field $U$ obeys the condition $\overline{U}=U^{-1}$. Using contraction rules (31), we find

[TABLE]

The action for the slow modes after integration over fast modes $W$ can be found as

[TABLE]

Here $\langle\dots\rangle_{W}$ denotes the averaging over fast modes $W$. In what follows, in the expansion in the right-hand side of Eq. (33) we neglect all terms except the lowest order one in the impurity concentration, $\langle\mathcal{S}_{\text{mag}}\rangle_{W}$. The smallness of the omitted terms is controlled by the condition $n_{s}\xi^{2}/g\ll 1$, where $\xi\sim\sqrt{D/T_{c}}$ is the superconducting coherence length in the dirty limit. As we shall see below, it is the first term in the right-hand side of Eq. (33) that is responsible for the logarithmic renormalization of $\mathcal{S}_{\text{mag}}$.

A.1 Operators of the second order in $Q$

In the Born approximation (first order in $\alpha$) we need to consider the operators $T_{2}$ and $T_{11}$ with two $Q$ matrices involved. Their contribution to $\mathcal{S}_{\text{mag}}$ is controlled by the coefficients $\gamma_{2}$ and $\gamma_{11}$ with the initial conditions following from (11): $\gamma_{2}(0)=1/4$ and $\gamma_{11}(0)=-1/8$.

Using the contraction rules (32), we find

[TABLE]

The operators $T_{2}$ and $T_{11}$ transform into each other under the renormalization. We note that under renormalization the operators with the same or fewer number of $Q$ matrices are generated only.

The eigenvalues of $M_{2}$ are equal to $2$ and $-1$. The eigenvalue $2$ corresponds to the operator $\hat{C}{T}_{2}-\hat{C}T_{11}/2$:

[TABLE]

The operator $\hat{C}{T}_{2}-\hat{C}T_{11}/2$ is known to be a pure scaling operator beyond the lowest order perturbation theory Wegner (1980); H$\mathrm{\ddot{o}}$f and Wegner (1986); Wegner (1987a, b).

We emphasize that the operators of the second order in $Q$ enter the magnetic part of the action, Eq. (11), precisely in combination $\hat{C}{T}_{2}-\hat{C}T_{11}/2$. This implies that

[TABLE]

A.2 Operators of the fourth order in $Q$

The next nontrivial order in $\alpha$ involves operators which are of the fourth order in $Q$. Their flow is described by the system

[TABLE]

The operators of the forth order in $Q$ are mixed under the renormalization. In addition, the operators of the second order in $Q$ are generated.

The system of equations (37) can be cast in the matrix form

[TABLE]

Here we used Eq. (34). We emphasize that the matrix $M_{4}$ is the upper triangular block matrix. This reflects the fact that under renormalization the operators with the same or fewer number of $Q$ matrices are generated only.

The matrix $M_{4}$ has the following eigenvalues: 12, 5, 2, 2, $-1$, $-1$, and $-6$. The largest eigenvalue $12$ corresponds to the operator $\hat{C}T_{4}-(2/3)\hat{C}T_{31}-(1/4)\hat{C}T_{22}+(1/4)\hat{C}T_{211}-(1/48)\hat{C}T_{1111}$. It is known that this operator is the pure scaling operator from arguments based on the group representation theory Wegner (1980); H$\mathrm{\ddot{o}}$f and Wegner (1986); Wegner (1987a, b). It is worth emphasizing that the operators of the forth order in $Q$ enter the magnetic part of the action, Eq. (11), precisely in the combination $\hat{C}T_{4}-(2/3)\hat{C}T_{31}-(1/4)\hat{C}T_{22}+(1/4)\hat{C}T_{211}-(1/48)\hat{C}T_{1111}$. This implies that the coefficients in the action (13) are simply

[TABLE]

A.3 Renormalization of operators of arbitrary order in $Q$

In general, one can derive the following set of renormalization group equations:

[TABLE]

[TABLE]

[TABLE]

and so on. Using these equations we find for the renormalization of the action:

[TABLE]

Note that all terms in Eq. (43) which contain $\alpha$ cancel each other.

Appendix B The renormalized action $\mathcal{S}_{\text{mag}}$

In this Appendix we present the details of derivation of Eq. (8).

All in all, we find from Eq. (43) that the coefficients $\gamma_{k_{1}k_{2}\dots k_{q}}$, where $k_{1}+k_{2}+\dots+k_{q}=n$, behaves in the same way:

[TABLE]

In what follows we are interested in the mean-field analysis of the renormalized action (13) for which the singlet sector of the theory is important only. Therefore, one can operate with $Q$ matrix which is the unit matrix in the spin space, $Q=Q_{0}\sigma_{0}$. Then averaging over directions of the impurity magnetization $\bm{n}$ becomes trivial. We find (all indices, $m$, $n$, $\dots$ are even)

[TABLE]

and so on. Then the renormalized action for magnetic impurities becomes

[TABLE]

Decoupling the Gaussian part with an auxiliary integral over $\lambda$ we obtain

[TABLE]

Now all summations become trivial:

[TABLE]

where

[TABLE]

Finally, we find

[TABLE]

where $\operatorname{tr}$ still includes summation over the spin space. This equation is equivalent to Eq. (8).

Appendix C The spin-flip rate, the transition temperature, and the density of states

In this Appendix we present the details of derivation of results for the spin-flip rate, the transition temperature, and the density of states.

C.1 The spin-flip rate near the transition temperature

According to Eq. (23), contrary to the usual case of magnetic disorder Abrikosov and Gor’kov (1960), the spin-flip rate in the presence of mesoscopic fluctuations acquires a weak logarithmic (in 2D) dependence on energy through the function $t=t(\varepsilon)$:

[TABLE]

Here we introduced $\mu=(2t-\ln\alpha)/(4\sqrt{t})$ and $\beta=4\sqrt{t}$. Expanding the denominator in the last integral in the right hand side of Eq. (51) in powers of $\exp(-2\pi|u|/\beta)$ and, then, performing integration over $u$, we find

[TABLE]

Here we introduce the function $f(z)=\exp(z^{2})[1-\operatorname{erf}(z)]$. At $\beta/\pi\ll 1$ we can use the expansion of the function $f(z)$ in series in $1/z$:

[TABLE]

Performing summation over $k$ in (52), we obtain

[TABLE]

In the limiting cases Eq. (54) reduces to

[TABLE]

Expanding the right hand side of Eq. (54) to the first order in $t$, we find

[TABLE]

where $1/\tau_{s0}=2\alpha n_{s}/[\pi\nu(1+\alpha)^{2}]$.

For $\beta/\pi\gg 1$, the sum $k$ in Eq. (52) reduces to the integral. Then, we find

[TABLE]

Next, for $\mu\beta\gg 1$, which holds for $t\gg 1$, we obtain

[TABLE]

C.2 The transition temperature

Knowledge of the effective spin-flip rate allows one to compute the dependence of the superconducting transition temperature on the spin-flip rate, $1/\tau_{s0}$, and potential disorder. Using Eqs. (18) and (22), we find the following equation for the transition temperature:

[TABLE]

where $\varepsilon=2\pi T_{c}(n+1/2)$. Performing formal expansion in the right hand side of Eq. (59) with respect of the difference $\tau_{s}^{-1}(\varepsilon)-\tau_{s}^{-1}(T_{c})$, we obtain

[TABLE]

where

[TABLE]

Since the effective spin-flip rate depends on the Matsubara energy $\varepsilon$ via $t(\varepsilon)=t(T_{c})-\frac{1}{\pi g}\ln[\pi(2n+1)]$, we can represent $\tau_{s}^{-1}(\varepsilon)$ as follows

[TABLE]

In the case $1/(2\pi T_{c}\tau_{s})\ll 1$, the sum in Eq. (66) is dominated by the term with $k=1$. The condition $g\gg 1$ allows us to consider in Eq. (66) the term with $l=1$ only. Therefore, we find

[TABLE]

where numerical constant

[TABLE]

Using Eq. (59), we find that the suppression of $T_{c}$ for the case $1/(2\pi T_{c}\tau_{s})\ll 1$ is given as

[TABLE]

As one can see, the correction to the expression for $T_{c}$ due to the dependence of the effective spin-flip rate on the Matsubara energy is negligible provided the conductance is large enough, $g\gg 1$.

In the opposite case, $1/(2\pi T_{c}\tau_{s})\gg 1$, we can integrate over $n$ in Eq. (66) instead of summation and find

[TABLE]

Provided the following condition

[TABLE]

holds, we can neglect the term $X_{\rm rest}$ in the right hand side of Eq. (59) in comparison with the $\ln 1/(2\pi T_{c}\tau_{s})$ which appears due to the di-gamma function.

All in all, the correction $X_{\rm rest}$ to the mean-field equation (59) for $T_{c}$ which stems from the energy dependence of the effective spin-flip rate can be neglected if the following inequality holds:

[TABLE]

Since our theory is valid for $t(T_{c})=[1/(\pi g)]\ln[{1}{(2\pi T_{c}\tau)}]\ll 1$ and $1/\tau_{s}\ll 1/\tau$, the condition (68) is always satisfied.

C.3 The density of states

The average DOS can be extracted from $\langle Q_{\varepsilon\varepsilon}\rangle$ analytically continued to the real energies $E$: $i\varepsilon\to E+i0^{+}$. The mean-field equation (21) can be written as

[TABLE]

where

[TABLE]

It is convenient to parametrize the spectral angle as $\theta_{\varepsilon}=\pi/2+i\psi$ such that the density of states becomes:

[TABLE]

For arbitrary values of $t$ and $\alpha$, Eq. (69) is a complicated integral equation which can be solved numerically. Below, we demonstrate how its solution and, consequently, the density of states, can be found analytically at $t\ll 1$.

At first, we rewrite the function $\mathcal{F}(\theta,\mu,\beta)$ as follows

[TABLE]

Here we used the relation 3.514.2 from the book Gradsteyn and Ryzhik (2000). Expanding the denominator in the last integral in the right hand side of Eq. (72) in powers of $\exp(-2\pi|u|/\beta)$ and, then, performing integration over $u$, we find

[TABLE]

where $f(z)=\exp(z^{2})[1-\operatorname{erf}(z)]$. Since we are interested in $\beta/\pi\lesssim 1$ we can use expansion of the function $f(z)$ in powers of $1/z$ (see Eq. (53)). Then performing summation over $k$ in the right hand side of Eq. (73), we find

[TABLE]

where

[TABLE]

While deriving Eq. (74) we used the following relation for the Euler di-gamma functions:

[TABLE]

We note that both the real and imaginary parts of the function $H(z)$ are exponentially small at $z\gg 1$.

Using the result (74) and making transformation $\varepsilon\to-iE$, we obtain the following form of the mean-field equation (21):

[TABLE]

where the function (cf. Eq. (26))

[TABLE]

is defined in such a way that the mean-field equation at $t=0$ is given as

[TABLE]

In what follows, we focus at the case

[TABLE]

in which the density of states has a finite gap $E_{g0}$ at $t=0$ Fominov and Skvortsov (2016). In this case, Eq.  (79) has a real solution $\psi$ for $|E|<E_{g0}$. The energy $E_{g0}$ and the corresponding value $\psi_{g0}$ are determined from the following equations:

[TABLE]

Since for $t\ll 1$ the arguments of the functions $H$ in Eq. (77) are large we can use the asymptotic expression for $H(z)$ at $z\gg 1$. In this way, we find

[TABLE]

We note that for $t\ll 1$, we can write

[TABLE]

C.3.1 The density of states near the band gap $E_{g0}$

The solution of Eq. (77) for $t\ll 1$ depends on the energy interval we are interested in. We start from the energies close to the bare gap edge $E_{g0}$. The function $\widetilde{F}_{E}(\psi,\alpha,t)$ has similar behaviour as the function $F_{E}(\psi,\alpha)$. Although at nonzero $t$ the density of states is finite at some energy, it is convenient to define the characteristic energy $E_{g}$ and corresponding angle $\psi_{g}$ which are the solutions of the following set of equations:

[TABLE]

For $t\ll 1$ we find that the difference between the characteristic energy $E_{g}$ and the bare gap $E_{g0}$ is given as

[TABLE]

In the Abrikosov-Gor’kov regime, $\alpha\ll\eta^{2/3}\ll 1$, where $\eta\equiv 1/(\tau_{s0\Delta})$, the above expression for the shift of the bare gap acquires the following simple form:

[TABLE]

Here we took into account that $\cosh\psi_{g0}=1/\eta^{1/3}$ and $E_{g0}=\Delta(1-\eta^{2/3})^{3/2}$.

Now we can find the dependence of the density of states on energy near $E_{g}$. Expanding the left hand side of Eq. (82) in $\epsilon=(E-E_{g})/\Delta$ and $\psi-\psi_{g}$, we find the following result for the density of states:

[TABLE]

Here we introduced the energy scale

[TABLE]

In the regime $\alpha\ll\eta^{2/3}\ll 1$ the result for the density of states for $|\epsilon|\ll\eta^{2/3}$ becomes

[TABLE]

Now it is instructive to compare our results for the density of states with the results of the instanton analysis Fominov and Skvortsov (2016); Marchetti and Simons (2002). The density of states due to instantons near the band gap $E_{g}$ is given as

[TABLE]

As one can see there is the characteristic energy scale $\Gamma={|\partial_{\psi_{g}}^{2}\widetilde{F}_{E_{g}}(\psi_{g},\alpha,t)|}/({g\cosh\psi_{g}})$ in Eq. (90). Using Eq. (89), we find

[TABLE]

Therefore, our contribution to the density of states dominates the instanton one near the band gap $E_{g}$ provided $\epsilon_{*}\gg\Gamma$. In the Abrikosov-Gor’kov regime, $\alpha\ll\eta^{2/3}\ll 1$, this condition becomes

[TABLE]

At the Fermi level our contribution to the density of states dominates the result due to instanton analysis since the latter involves the sheet resistance $1/g$ which is parametrically smaller than spreading resistance $t=\ln(\xi/l)/(2\pi g)$.

C.3.2 The density of states at low energies

At energies which are much smaller than the characteristic energy, $|E|\ll E_{g}$, the equation (82) without the right hand side has the real solutions only. We substitute $\psi=\psi^{\prime}+i\psi^{\prime\prime}$ with $\psi^{\prime\prime}\ll 1$ into Eq. (82) and splitting into the real and imaginary parts. Then we find

[TABLE]

The density of states can be found as

[TABLE]

We present the comparison between the density of states found from numerical solution of Eq. (69) and analytical result (93) in Fig. 5. To plot the curves in this figure we neglect the difference between $\psi_{g}$ and $\psi_{g0}$ as well as between $\widetilde{F}_{E}(\psi^{\prime},\alpha,t)$ and $F_{E}(\psi^{\prime},\alpha)$.

At $E=0$, $\psi^{\prime}=0$ is the solution of Eq. (C.3.2). Then from Eq. (93) we find the density of states at zero (Fermi) energy

[TABLE]

Appendix D The effect of termination of the multifractal spectrum

In this Appendix we discuss how the termination of the multifractal spectrum affects our results.

The result (44) for the coefficients $\gamma_{k_{1}k_{2}\dots k_{q}}$ is derived by consideration of the contributions related with $\langle\mathcal{S}_{\textrm{mag}}\rangle_{W}$. In this approximation operators $T_{k_{1}k_{2}\dots k_{q}}$ with given $n=k_{1}+\dots+k_{q}$ always transform under the renormalization group into linear combinations of operators $T_{l_{1}l_{2}\dots l_{q}}$ with $m=l_{1}+\dots+l_{q}\leqslant n$. Therefore, the renormalization group equations remain linear in coefficients $\gamma_{k_{1}k_{2}\dots k_{q}}$. In general, one needs to take into account terms which are nonlinear in $\mathcal{S}_{\textrm{mag}}$, e.g. $\langle[\mathcal{S}_{\textrm{mag}}]^{2}\rangle_{W}$. Then the fusion of two operators $T_{k_{1}k_{2}\dots k_{q}}$ and $T_{l_{1}l_{2}\dots l_{q}}$ into a single operator $T_{s_{1}s_{2}\dots s_{q}}$ with $s_{1}+\dots+s_{q}=n+m$ is possible. This renders the renormalization group equations for $\gamma_{k_{1}k_{2}\dots k_{q}}$ nonlinear Foster et al. (2009). This nonlinearity results in termination of the multifractal spectrum Mirlin (2000) which implies the following modification of Eq. (44):

[TABLE]

Here $n_{c}=\sqrt{2/t_{0}}\gg 1$ and $t_{0}=2/(\pi g)\ll 1$ denotes the bare resistance. The function $\textsf{y}_{n}$ obeys the following symmetry property: $\textsf{y}_{1-n}=\textsf{y}_{n}$ Gruzberg et al. (2013). Let us now define the function $\mathcal{G}(\lambda)$ as

[TABLE]

Then we find

[TABLE]

where $X$ is given by Eq. (49).

At $t\ll 1$ and $t/t_{0}\gg 1$, the function $\mathcal{G}(\lambda)$ can be written as

[TABLE]

where $\theta(z)$ denotes the Heaviside step function and $\lambda_{c}=t(2\sqrt{2/t_{0}}-1)\approx 2t\sqrt{2/t_{0}}$. This form of the function $\mathcal{G}(\lambda)$ implies that the integration over a in Eq. (21) is restricted to the range $\textsf{a}_{-}<\textsf{a}<\textsf{a}_{+}$, where $\textsf{a}_{\pm}=\alpha\exp(\pm 2\lambda_{c})$. Since for the existence of a finite density of states near the Fermi energy, vicinity of $\textsf{a}=1$ is important, this point should be within the range of integration over a in Eq. (21), i.e. $\textsf{a}_{-}<1<\textsf{a}_{+}$. The latter condition is fulfilled provided $4t\sqrt{2/t_{0}}\gg\ln(1/\alpha)$.