text:
# Superconductors with Magnetic Impurities: Instantons and Sub-gap States

## I Introduction

While the spectral properties of a singlet $`s`$-wave superconductor are largely unaffected by weak non-magnetic impurities , the pair-breaking effect of magnetic impurities leads to the gradual destruction of superconductivity. Remarkably, the suppression of the quasi-particle energy gap is more rapid than that of the superconducting order parameter, admitting the existence of a narrow ‘gapless’ superconducting phase  in which the quasi-particle energy gap is destroyed while the superconducting order parameter remains non-zero. Now, according to the conventional (mean-field) description formulated in the seminal work of Abrikosov and Gor’kov (AG), an energy gap is maintained up to a critical concentration of magnetic impurities (at $`T=0`$, $`91`$% of the critical concentration at which superconductivity is destroyed). Yet, being unprotected by the Anderson theorem, it seems likely that the gap structure predicted by the mean-field theory is untenable and must be destroyed by ‘optimal’ fluctuations of the random impurity potential. Indeed, since the pioneering work of AG, several authors  have explored the nature of ‘sub-gap’ states in the superconducting system. The aim of this work is to present a detailed investigation of the spectrum and profile of sub-gap states in superconductors subject to a weak magnetic and non-magnetic impurity potential, thus systematically improving upon the mean-field theory of AG. Our preliminary findings have already been reported in a recent letter .

In the earliest works on the subject , attention was focussed on the the influence of strong magnetic impurities. In particular, in the unitarity limit, it was shown that a single magnetic impurity leads to the local suppression of the order parameter and creates a bound sub-gap quasi-particle state . For a finite impurity concentration, these intra-gap states broaden into a band  merging smoothly with the continuum bulk states.

By contrast, starting with a weak magnetic impurity distribution (i.e. one in which the magnetic scattering can be treated within the Born approximation), the mean-field theory of AG  predicts a gradual suppression of the quasi-particle energy gap. Defining the dimensionless parameter

$`\zeta {\displaystyle \frac{1}{\tau _s|\mathrm{\Delta }|}},`$ (1)
where $`|\mathrm{\Delta }|`$ represents the value of the homogeneous self-consistent order parameter, and $`\tau _s`$ denotes the Born scattering time due to magnetic impurities, the AG theory shows the gap to follow the relation

$`E_{\mathrm{gap}}(\tau _s)=|\mathrm{\Delta }|\left(1\zeta ^{2/3}\right)^{3/2}`$ (2)
showing an onset of the gapless region at $`\zeta =1`$ (note $`\mathrm{}=1`$ throughout). Within the same mean-field theory, for $`\zeta 1`$, the self-consistent order parameter varies as $`|\mathrm{\Delta }|=|\overline{\mathrm{\Delta }}|\mathrm{exp}[\pi \zeta /4]`$ where $`|\overline{\mathrm{\Delta }}|`$ represents the order parameter of the clean superconductor, confirming that the order parameter is finite at the onset of the gapless phase. The precise variation of the quasi-particle energy gap is compared to that of the self-consistent order parameter in Fig. 1. Staying within the framework of the mean-field theory, one can obtain a smooth interpolation between the strong and weak impurity scattering behaviors .

However, even for weak disorder, it is apparent that optimal fluctuations of the random potential must generate sub-gap states in the interval $`0<\zeta <1`$, and therefore provide non-perturbative corrections to the self-consistent Born approximation used by AG. Extending the arguments of Balatsky and Trugman , a fluctuation of the random potential which leads to an effective Born scattering rate $`1/\tau _s^{}`$ in excess of $`1/\tau _s`$ over a range set by the superconducting coherence length,

$`\xi =\left({\displaystyle \frac{D}{|\mathrm{\Delta }|}}\right)^{1/2},`$ (3)
induces quasi-particle states down to energies $`E_{\mathrm{gap}}(\tau _s^{})`$. These sub-gap states are localized, being bound to the region or ‘droplet’ where the scattering rate is large.

The situation bears comparison with band tail states in semi-conductors. Here rare or optimal configurations of the random impurity potential generate bound states, known as Lifshitz tail states , which extend below the band edge. However, the correspondence is, to some extent, superficial: band tail states in semi-conductors are typically associated with smoothly varying, nodeless wavefunctions. By contrast, the tail states below the superconducting gap involve the superposition of states around the Fermi level. As such, one expects these states to be rapidly oscillating on the scale of the Fermi wavelength $`\lambda _F`$, but modulated by an envelope which is localized on the scale of the coherence length $`\xi `$. This difference is not incidental. Firstly, unlike the semi-conductor, one expects the energy dependence of the density of states in the tail region below the mean-field gap edge to be ‘universal’, independent of the nature of the weak impurity distribution, but dependent only on the pair-breaking parameter $`\zeta `$. Secondly, as we will see, one can not expect a straightforward extension of existing theories  of the Lifshitz tails to describe the profile of tail states in the superconductor.

In the BCS approximation, the random system we consider is specified by the Gor’kov Hamiltonian

$$\widehat{H}=\left(\begin{array}{cc}\widehat{H}_0& |\mathrm{\Delta }|\sigma _2^{\mathrm{sp}}\\ |\mathrm{\Delta }|\sigma _2^{\mathrm{sp}}& \widehat{H}_0^T\end{array}\right)_{\mathrm{ph}}$$
(4)
where the matrix components index the particle/hole content, and

$`\widehat{H}_0={\displaystyle \frac{\widehat{𝐩}^2}{2m}}ϵ_F+V(𝐫)+J𝐒(𝐫)\sigma ^{\mathrm{sp}}`$
denotes the normal component of the random Hamiltonian. Here we take $`\mathrm{\Delta }=g_\mathrm{\Delta }c_{}c_{}`$ to be spatially homogeneous and determined self-consistently from the conventional AG mean-field theory. We will assess the validity of this assumption below. In addition to the weak potential impurity distribution $`V(𝐫)`$, the particles experience a quenched random magnetic impurity distribution $`J𝐒(𝐫)`$ where $`J`$ represents the exchange coupling and Pauli matrices $`\{\sigma _i^{\mathrm{sp}}\}`$ operate on the spin indices. The magnetic $`𝐒(𝐫)`$ and non-magnetic $`V(𝐫)`$ random impurity potentials are both taken to be Gaussian $`\delta `$-correlated with zero mean and variance

$`JS_\alpha (𝐫)JS_\beta (𝐫^{})_S`$ $`=`$ $`{\displaystyle \frac{1}{6\pi \nu \tau _s}}\delta ^d(𝐫𝐫^{})\delta _{\alpha \beta }`$
$`V(𝐫)V(𝐫^{})_V`$ $`=`$ $`{\displaystyle \frac{1}{2\pi \nu \tau }}\delta ^d(𝐫𝐫^{})`$
respectively, and $`\nu `$ represents the average density of states (DoS) per spin of the normal system.

To keep our discussion simple, and to make contact with the AG theory, we will take the quenched distribution of magnetic impurities to be ‘classical’ and non-interacting throughout. For practical purposes, this entails the consideration of structures where both the Kondo temperature  and, more significantly, the RKKY induced spin glass temperature  are smaller than the relevant energy scales of the superconductor. The remaining energy scales are arranged in the quasi-classical and dirty limits:

$`ϵ_F1/\tau (|\mathrm{\Delta }|,1/\tau _s)`$ (5)
where $`\tau `$ represents the transport time associated with non-magnetic impurities. We remark that these limits are not compatible with the situation in which the magnetic impurities provide the only source of scattering, i.e. $`\xi \mathrm{}=v_F\tau _s`$. We should not, therefore, expect a straightforward comparison with the analysis of Ref. .

Within the approximations above we will derive a quasi-classical field theory of the disordered superconducting system. Following the approach of Ref.  we will express spectral properties of the system in terms of an intermediate energy scale action which accommodates the quantum interference properties of the superconducting system. By investigating stationary inhomogeneous instanton field configurations of the action, we will expose the structure and profile of the sub-gap states, and thereby obtain an analytical expression for the DoS to exponential accuracy. Finally, we will comment on the universality of the present scheme by applying the same technique to investigate ‘gap fluctuations’ in a quantum dot coupled by open channels to a superconducting terminal.

Although the analysis is straightforward, the technology is somewhat involved. We have therefore decided to summarize the main conclusions of this investigation here in the introduction. In particular, we will show that, in the vicinity of the mean-field gap edge, the sub-gap DoS of the $`d`$-dimensional system is dominated by tail states which are confined to droplets of size

$`r_0(ϵ){\displaystyle \frac{\xi }{(1\zeta ^{2/3})^{1/4}}}\left({\displaystyle \frac{|\mathrm{\Delta }|}{E_{\mathrm{gap}}ϵ}}\right)^{1/4}.`$
The corresponding sub-gap DoS takes the form

$`{\displaystyle \frac{\nu (ϵ<E_{\mathrm{gap}})}{\nu }}\mathrm{exp}\left[4\pi g(\xi /L)^{d2}f_d(\zeta )\left({\displaystyle \frac{E_{\mathrm{gap}}ϵ}{|\mathrm{\Delta }|}}\right)^{(6d)/4}\right]`$
where $`g=\nu DL^{d2}`$ denotes the bare dimensionless conductance and $`f_d(\zeta )=a_d\zeta ^{2/3}(1\zeta ^{2/3})^{(2+d)/8}`$ represents a dimensionless function of the control parameter $`\zeta `$ ($`a_d`$ const.). When reparameterized in terms of the DoS just above the mean-field gap edge

$`\nu (ϵ>E_{\mathrm{gap}}){\displaystyle \frac{1}{\pi L^d}}\sqrt{{\displaystyle \frac{ϵE_{\mathrm{gap}}}{\mathrm{\Delta }_g^3}}}`$ (6)
where

$`\mathrm{\Delta }_g^{3/2}=4\pi \nu L^d\sqrt{{\displaystyle \frac{2}{3|\mathrm{\Delta }|}}}\zeta ^{2/3}(1\zeta ^{2/3})^{1/4}`$
the expression for the sub-gap DoS can be brought to the more compact form

$`{\displaystyle \frac{\nu (ϵ<E_{\mathrm{gap}})}{\nu }}\mathrm{exp}\left[\stackrel{~}{a}_d\left({\displaystyle \frac{r_0}{L}}\right)^d\left({\displaystyle \frac{E_{\mathrm{gap}}ϵ}{\mathrm{\Delta }_g}}\right)^{3/2}\right]`$
with $`\stackrel{~}{a}_d`$ some numerical constant.

In the zero-dimensional system, although the interpretation of the optimal fluctuation as a localized droplet is no longer appropriate, the expression above correctly interpolates to $`d=0`$ and coincides with the universal expression for gap fluctuations proposed in Ref. . The surprising dependence of the result on the dimensionless distance from the mean-field gap is one of the reasons why a Lifshitz argument appears difficult to construct for this problem.

The paper is organized as follows: in section II a theory of the statistical properties of the Gor’kov Green function is developed within the framework of a supersymmetric field theory involving a non-linear $`\sigma `$-model functional. Here we follow closely the analysis of Ref.  and  (c.f. Ref. ). As a result we identify the conventional AG mean-field equation with the homogeneous saddle-point equation of the effective action. In section III we show that the non-vanishing of the DoS beneath the gap predicted by the standard AG theory is due to the appearance of inhomogeneous “instanton” saddle points of finite action. We identify the profile of these instantons with the envelope modulating the quasi-classical sub-gap states. The instanton configurations considered break the underlying supersymmetry of the action, and a detailed examination of fluctuations is required to understand how supersymmetry is restored by zero modes, as well as to appreciate fully how these configurations are able to contribute to the DoS. In section IV we examine the zero dimensional limit of the problem and compare our results to the recent literature . In doing so, we provide an explanation of the universal results reported in Ref. , and discuss the universality of the $`d>0`$ result. Finally, in section V we speculate on potential generalizations of the results presented here.

## II Field Theory of the Disordered Superconductor

The construction of the field theory of the disordered superconductor follows the quasi-classical method of Eilenberger  and Usadel  elevated to the level of an effective action. The starting point of the analysis is the generating functional for the single-particle Gor’kov Green function for the non-interacting quasi-particle Bogoluibov Hamiltonian. Here we borrow our notation from Ref. .

### A Generating Functional

Single-particle properties of the Gor’kov Hamiltonian (4) are obtained from the generating functional

$$𝒵[J]=D(\overline{\psi },\psi )e^{{\scriptscriptstyle 𝑑𝐫\left(i\overline{\psi }(\widehat{H}ϵ_{})\psi +\overline{\psi }J+\overline{J}\psi \right)}},$$
(7)
where $`ϵ_{}ϵi0`$ and the supervector fields have the internal structure $`\overline{\psi }=\left(\begin{array}{cccc}\overline{\psi }_{}& \overline{\psi }_{}& \psi _{}& \psi _{}\end{array}\right)`$, $`\psi ^T=\left(\begin{array}{cccc}\psi _{}& \psi _{}& \overline{\psi }_{}& \overline{\psi }_{}\end{array}\right)`$. As usual, by introducing both commuting and anticommuting elements the normalization $`𝒵[0]=1`$ is assured. (The generalization to the consideration of higher point functions, which involves the extension of the field space, follows straightforwardly.)

To condense the notation, it is convenient to perform the rotation $`\psi \psi ^{}=U\psi `$, $`\overline{\psi }\overline{\psi }^{}=\overline{\psi }U^{}`$ with

$`U=\left(\begin{array}{cc}1& 0\\ 0& i\sigma _2^{\mathrm{sp}}\end{array}\right)_{\mathrm{ph}},`$
after which the Gor’kov Hamiltonian takes the form

$`\widehat{H}=\left({\displaystyle \frac{\widehat{𝐩}^2}{2m}}+V(𝐫)ϵ_F\right)\sigma _3^{\mathrm{ph}}+J𝐒(𝐫)\sigma ^{\mathrm{sp}}+|\mathrm{\Delta }|\sigma _2^{\mathrm{ph}}.`$
The unusual phase coherence properties of the superconducting system rely on the particle/hole or charge conjugation symmetry

$$\widehat{H}=\sigma _2^{\mathrm{ph}}\sigma _2^{\mathrm{sp}}\widehat{H}^T\sigma _2^{\mathrm{sp}}\sigma _2^{\mathrm{ph}}.$$
(8)
To easily accommodate these effects, it is convenient to introduce the further space doubling,

$`2\overline{\psi }(\widehat{H}ϵ_{})\psi =\overline{\psi }(\widehat{H}ϵ_{})\psi +\psi ^T(\widehat{H}^Tϵ_{})\overline{\psi }^T`$ (9)
$`=\overline{\psi }(\widehat{H}ϵ_{})\psi +\psi ^T(\sigma _2^{\mathrm{ph}}\sigma _2^{\mathrm{sp}}\widehat{H}\sigma _2^{\mathrm{sp}}\sigma _2^{\mathrm{ph}}ϵ_{})\overline{\psi }^T`$ (10)
$`=2\overline{\mathrm{\Psi }}(\widehat{H}ϵ_{}\sigma _3^{\mathrm{cc}})\mathrm{\Psi }`$ (11)
where, defining the Pauli matrix $`\sigma _3^{\mathrm{cc}}`$ which operates in the charge-conjugation (cc) space,

$`\mathrm{\Psi }={\displaystyle \frac{1}{\sqrt{2}}}\left(\begin{array}{c}\psi \\ \sigma _2^{\mathrm{ph}}\sigma _2^{\mathrm{sp}}\overline{\psi }^T\end{array}\right)_{\mathrm{cc}},`$
$`\overline{\mathrm{\Psi }}={\displaystyle \frac{1}{\sqrt{2}}}\left(\begin{array}{cc}\overline{\psi }& \psi ^T\sigma _2^{\mathrm{ph}}\sigma _2^{\mathrm{sp}}\end{array}\right)_{\mathrm{cc}}.`$
This completes the formulation of the generating functional as a field integral involving $`16`$-component supervector fields $`\mathrm{\Psi }`$ and $`\overline{\mathrm{\Psi }}`$. The latter obey the symmetry relations

$`\mathrm{\Psi }=\sigma _2^{\mathrm{ph}}\sigma _2^{\mathrm{sp}}\gamma \overline{\mathrm{\Psi }}^T,\overline{\mathrm{\Psi }}=\mathrm{\Psi }^T\sigma _2^{\mathrm{ph}}\sigma _2^{\mathrm{sp}}\gamma ^1`$ (12)
with $`\gamma =E_{\mathrm{bb}}i\sigma _2^{\mathrm{cc}}E_{\mathrm{ff}}\sigma _1^{\mathrm{cc}}`$, where $`E_{\mathrm{bb}}=\mathrm{diag}(1,0)_{\mathrm{bf}}`$ and $`E_{\mathrm{ff}}=\mathrm{diag}(0,1)_{\mathrm{bf}}`$ project into the Boson-Boson and Fermion-Fermion sectors respectively.

### B Ensemble Averaging

Cast as a field integral, the impurity average of the generating functional over the Gaussian distributed random impurity potentials is straightforward. Separating the regular from the disordered components of the Hamiltonian $`\widehat{H}=\widehat{}_0+V(𝐫)\sigma _3^{\mathrm{ph}}+J𝐒(𝐫)\sigma ^{\mathrm{sp}}`$, an ensemble average over the random potentials obtains

$`𝒵[0]_{V,S}={\displaystyle D(\overline{\mathrm{\Psi }},\mathrm{\Psi })\mathrm{exp}\left[𝑑𝐫\left(i\overline{\mathrm{\Psi }}(\widehat{}_0ϵ_{}\sigma _3^{\mathrm{cc}})\mathrm{\Psi }\frac{1}{4\pi \nu \tau }(\overline{\mathrm{\Psi }}\sigma _3^{\mathrm{ph}}\mathrm{\Psi })^2\frac{1}{12\pi \nu \tau _s}(\overline{\mathrm{\Psi }}\sigma ^{\mathrm{sp}}\mathrm{\Psi })^2\right)\right]}.`$ (13)
The interactions generating by the impurity averaging can be decoupled by the introduction of a Hubbard-Stratonovich field. Beginning with the non-magnetic disorder, slow modes of the action are identified by rewriting the action in the approximate form

$`{\displaystyle \frac{1}{4\pi \nu \tau }}{\displaystyle 𝑑𝐫\left(\overline{\mathrm{\Psi }}\sigma _3^{\mathrm{ph}}\mathrm{\Psi }\right)^2}{\displaystyle \frac{1}{2\pi \nu \tau }}{\displaystyle \underset{|𝐪|<\mathrm{}^1}{}}\mathrm{str}\zeta (𝐪)\zeta (𝐪),`$
where $`\zeta (𝐪)=_𝐤\mathrm{\Psi }(𝐤)\overline{\mathrm{\Psi }}(𝐤+𝐪)\sigma _3^{\mathrm{ph}}`$. The latter can be decoupled by the slowly varying $`16\times 16`$ supermatrix field $`Q(𝐫)`$ according to the identity

$`\mathrm{exp}\left[{\displaystyle \frac{1}{2\pi \nu \tau }}{\displaystyle \underset{𝐪}{}}\mathrm{str}\zeta (𝐪)\zeta (𝐪)\right]`$
$`={\displaystyle DQ\mathrm{exp}\left[\underset{𝐪}{}\mathrm{str}\left(\frac{\pi \nu }{8\tau }Q(𝐪)Q(𝐪)\frac{1}{2\tau }Q(𝐪)\zeta (𝐪)\right)\right]}.`$
The symmetry properties of the fields $`Q(𝐫)`$ are inherited from the dyadic product $`\mathrm{\Psi }(𝐫)\overline{\mathrm{\Psi }}(𝐫)\sigma _3^{\mathrm{ph}}`$. Making use of Eq. (12) one finds the symmetry relations

$$Q=\sigma _1^{\mathrm{ph}}\sigma _2^{\mathrm{sp}}\gamma Q^T\gamma ^1\sigma _1^{\mathrm{ph}}\sigma _2^{\mathrm{sp}}.$$
(14)
The interaction generated by the magnetic impurity averaging can be treated  by performing all possible pairings and making use of the saddle-point approximation $`Q(𝐫)=2\mathrm{\Psi }(𝐫)\overline{\mathrm{\Psi }}(𝐫)\sigma _3^{\mathrm{ph}}_\mathrm{\Psi }/\pi \nu `$. This leads to the replacement

$`{\displaystyle \frac{1}{12\pi \nu \tau _s}}{\displaystyle 𝑑𝐫\left(\overline{\mathrm{\Psi }}\sigma ^{\mathrm{sp}}\mathrm{\Psi }\right)^2}`$
$`{\displaystyle \frac{\pi \nu }{24\tau _s}}{\displaystyle 𝑑𝐫\mathrm{str}\left(Q\sigma _3^{\mathrm{ph}}\sigma ^{\mathrm{sp}}\right)^2}.`$
Such an approximation, which neglects pairings at non-coincident points is allowed by the strong inequality $`(\mathrm{}/\xi )^d1`$. In addition we discard the contraction $`\overline{\mathrm{\Psi }}\sigma ^{\mathrm{sp}}\mathrm{\Psi }_\mathrm{\Psi }`$. The term generated by this procedure could in any case be decoupled by a slow Bosonic field $`𝐒(𝐫)`$ which would immediately be set to zero for the singlet saddle-points that will be the basis of this paper.

Gaussian in the fields $`\mathrm{\Psi }`$ and $`\overline{\mathrm{\Psi }}`$, the functional integration can be performed explicitly after which one obtains $`𝒵[0]_{V,S}=DQ\mathrm{exp}(S[Q])`$ where

$`S[Q]={\displaystyle 𝑑𝐫\left[\frac{\pi \nu }{8\tau }\mathrm{str}Q^2\frac{1}{2}\mathrm{str}\mathrm{ln}\left(\sigma _3^{\mathrm{ph}}(\widehat{}_0ϵ_{}\sigma _3^{\mathrm{cc}})+\frac{i}{2\tau }Q\right)\frac{\pi \nu }{24\tau _s}\mathrm{str}(Q\sigma _3^{\mathrm{ph}}\sigma ^{\mathrm{sp}})^2\right]}.`$
Further progress is possible only within a saddle-point approximation. Following Ref. , the saddle-point analysis will be conducted in a two-step process:

### C Saddle-point Approximation and the Non-linear $`\sigma `$-model

The first task is to make use of the quasi-classical parameter $`ϵ_F\tau 1`$ to construct an intermediate energy scale action. Dropping the symmetry breaking perturbations $`ϵ`$, $`|\mathrm{\Delta }|`$ and $`1/\tau _s`$, a variation of the action with respect to fluctuations of $`Q`$ obtains the saddle-point equation

$`Q(𝐫)={\displaystyle \frac{i}{\pi \nu }}𝐫|[i0\sigma _3^{\mathrm{cc}}\sigma _3^{\mathrm{ph}}+\widehat{h}_0+iQ/2\tau ]^1|𝐫,`$
where $`\widehat{h}_0=\widehat{𝐩}^2/2mϵ_\mathrm{f}`$. Taking into account the analytical properties of the Green function and the quasi-classical limit $`ϵ_\mathrm{f}\tau 1`$, one obtains the solution $`Q_0=\sigma _3^{\mathrm{ph}}\sigma _3^{\mathrm{cc}}`$. However, as usual, this saddle-point is not unique. In particular, $`Q=TQ_0T^1`$ is also a solution, for any constant matrix $`T`$ which is consistent with the fundamental symmetry (14) for $`Q`$. Transverse fluctuations of $`Q`$ away from the $`Q^2=𝟙`$ manifold may be integrated out within the saddle-point approximation due to the large parameter $`\nu L^d/\tau 1`$, where $`L`$ is the system size . Restricting attention to the manifold generated by the non-linear constraint $`Q^2=𝟙`$, an effective action is obtained by allowing $`T`$ to vary in space and expanding to second order in gradients of $`T`$, and first order in $`ϵ`$, $`|\mathrm{\Delta }|`$, and $`1/\tau _s`$

$$S[Q]=\frac{\pi \nu }{8}𝑑𝐫\mathrm{str}\left[D(Q)^24i\left(ϵ_{}\sigma _3^{\mathrm{cc}}+|\mathrm{\Delta }|\sigma _2^{\mathrm{ph}}\right)\sigma _3^{\mathrm{ph}}Q\frac{1}{3\tau _s}\left(Q\sigma _3^{\mathrm{ph}}\sigma ^{\mathrm{sp}}\right)^2\right].$$
(15)
Here $`D=v_F^2\tau /d`$ represents the classical diffusion constant associated with the non-magnetic impurities. In particular, the quasi-particle DoS is obtained from the functional integral

$`\nu (ϵ,𝐫)_{V,S}={\displaystyle \frac{\nu }{4}}\mathrm{Re}\mathrm{str}\left(\sigma _3^{\mathrm{bf}}\sigma _3^{\mathrm{ph}}\sigma _3^{\mathrm{cc}}Q(𝐫)\right)_Q.`$ (16)
The numerical factor leads to a DoS of $`4\nu `$ for the system as $`|ϵ|\mathrm{}`$. This is because both the particle-hole structure of the original Bogoliubov Hamiltonian and the spin each cause a doubling of the DoS.

This completes the derivation of the intermediate energy scale action. In the following section we will investigate the action (15) within mean-field theory and the influence of soft fluctuations on the low-energy properties of the system.

### D AG Mean-Field Theory and Fluctuations

The presence of symmetry-breaking terms in (15) originating from the order parameter and the magnetic impurity potential means that a mean-field analysis is already non-trivial. To assimilate the effect of these terms, and to establish contact with the AG theory, it is necessary to explore the saddle-point equation.

Varying the non-linear $`\sigma `$-model action with respect to fluctuations of $`Q`$, subject to the non-linear constraint, one obtains the saddle-point or mean-field equation,

$`D\left(QQ\right)+i[Q,ϵ_{}\sigma _3^{\mathrm{cc}}\sigma _3^{\mathrm{ph}}+i|\mathrm{\Delta }|\sigma _1^{\mathrm{ph}}]`$
$`+{\displaystyle \frac{1}{6\tau _s}}[Q,\sigma _3^{\mathrm{ph}}\sigma ^{\mathrm{sp}}Q\sigma _3^{\mathrm{ph}}\sigma ^{\mathrm{sp}}]=0.`$
With the Ansatz:

$`Q_{\mathrm{mf}}=\left[\sigma _3^{\mathrm{cc}}\sigma _3^{\mathrm{ph}}\mathrm{cosh}\widehat{\theta }+i\sigma _1^{\mathrm{ph}}\mathrm{sinh}\widehat{\theta }\right]𝟙^{\mathrm{sp}},`$
where the elements $`\widehat{\theta }=\mathrm{diag}(\theta _1,i\theta )_{\mathrm{bf}}`$ are diagonal in the superspace, the saddle-point equation decouples into Boson-Boson and Fermion-Fermion sectors, and takes the form

$`_{𝐫/\xi }^2\widehat{\theta }+2i\left(\mathrm{cosh}\widehat{\theta }{\displaystyle \frac{ϵ}{|\mathrm{\Delta }|}}\mathrm{sinh}\widehat{\theta }\right)\zeta \mathrm{sinh}(2\widehat{\theta })=0,`$ (17)
a result reminiscent of the Usadel equation of quasi-classical superconductivity . This is no coincidence: when subject to an inhomogeneous order parameter, the same effective action (15) describes the proximity effect in a hybrid normal/superconducting compound . In the present context, the spatially homogeneous form of Eq. (17) should be combined with the self-consistent equation for the order parameter

$`|\mathrm{\Delta }|={\displaystyle \frac{\pi \nu g_\mathrm{\Delta }}{\beta }}{\displaystyle \underset{n}{}}\mathrm{sin}\theta _n.`$
Here $`g_\mathrm{\Delta }`$ is the BCS coupling constant and $`\theta _n`$ indicates that the solution of the Usadel equation is taken at Matsubara frequency $`ϵiϵ_n`$. These two equations then coincide with the mean-field equations obtained by AG, and correctly reproduce the known form of $`E_{\mathrm{gap}}`$ as specified by Eq. (2).

The AG solution is not unique: for $`ϵ0`$, the saddle-point equation admits an entire manifold of homogeneous solutions parameterized by the transformations $`Q=TQ_{\mathrm{mf}}T^1`$ where now $`T=𝟙_{\mathrm{ph}}𝟙_{\mathrm{sp}}𝕥`$ and $`t=\gamma (t^1)^T\gamma ^1`$. Soft fluctuations of the fields, which are controlled by a non-linear $`\sigma `$-model defined on the manifold $`T\mathrm{OSp}(2|2)/\mathrm{GL}(1|1)`$ (symmetry class D in the classification of Ref. ), control the low-energy, long-range properties of the gapless system giving rise to unusual localization and spectral properties. (For a comprehensive discussion of the physics of the gapless phase, we refer to Refs. .) The phenomenology of Class D in the context of the gapless phase of the present system will be the subject of a forthcoming paper . Here we focus on the gapped phase, where Class D fluctuations are just one of a host of massive modes that are unimportant in the description of sub-gap states.

This completes the formal description of the bulk superconducting phase. The solution of the AG mean-field equation provides an adequate description of the bulk extended states. Soft fluctuations around the AG mean-field describe phase coherence effects due to quantum interference. However, within the present scheme it is not yet clear how to accommodate sub-gap states in the gapped phase of the AG theory. To identify such states, it is necessary to return to the saddle-point equation (17) and seek spatially inhomogeneous solutions.

## III Instantons and sub-gap states

Although the reduction and eventual destruction of the quasi-particle energy gap predicted by the AG mean-field theory can be reasonably justified on purely physical grounds, the integrity of the gap of the range $`0<\zeta <1`$ is less credible. Once time-reversal is symmetry is broken and the protection of Anderson’s theorem is lost, there remains no reason why a sharp gap should persist. Add to this the observation that the spin scattering rate must be subject to spatial fluctuations from the average value $`1/\tau _s`$, and one concludes that corrections to the DoS predicted by the AG theory must lead to the appearance of sub-gap states analogous to “band tails” in a disordered semiconductor .

This analogy is of course not new  nor, as far as practical calculation in the present formulation is concerned, is it particularly deep. This is because all averages have already been taken, so we can not look for an optimal fluctuation of some potential, as in the classic approaches to the study of band tail states in disordered semi-conductors . However, these studies hint at how one can proceed.

Band tail states in semi-conductors can be studied within the same functional integral formulation. In particular, the generating function of the single-particle Green function of a normal disordered conductor can be presented in the form of a supersymmetric field integral

$`𝒵[0]={\displaystyle D(\mathrm{\Psi },\overline{\mathrm{\Psi }})\mathrm{exp}\left[i𝑑𝐫\overline{\mathrm{\Psi }}\left(ϵ_+\frac{\widehat{𝐩}^2}{2m}V(𝐫)\right)\mathrm{\Psi }\right]},`$
where, once again, the random impurity distribution is drawn from a Gaussian $`\delta `$-correlated white-noise impurity potential. The optimal fluctuation method involves minimizing the action with respect to fluctuations in the fields $`\mathrm{\Psi }`$ and potential $`V`$. This involves seeking inhomogeneous solutions of the non-linear Schrödinger equation

$`\left(ϵ{\displaystyle \frac{\widehat{𝐩}^2}{2m}}V(𝐫)\right)\mathrm{\Psi }=0,`$
where the corresponding optimal potential is determined self-consistently by the relation $`V(𝐫)=|\mathrm{\Psi }(𝐫)|^2/2\pi \nu \tau `$. In the supersymmetric formulation, band tail states are identified with supersymmetry broken inhomogeneous solutions of the saddle-point equation (see Cardy  and Affleck ). Indeed, the anticipated exponential suppression of the DoS necessitates a breaking of supersymmetry to support a finite action. Here the phrase “supersymmetry breaking” is potentially misleading. We use it only to refer to *field configurations*, ubiquitous in the problems under discussion here, that do not respect the parity between Bose and Fermi degrees of freedom. However, any such configuration is just one member of a degenerate manifold differing by supersymmetric transformations. The latter maintain the invariance of the generating functional $`𝒵[0]`$ under global supersymmetric transformations.

What does this tell us about the identification of optimal fluctuations and sub-gap states in the superconductor? Following the analysis above, one might guess that sub-gap states are associated with inhomogeneous configurations of the $`\mathrm{\Psi }`$ field action. However, we anticipate that optimal solutions corresponding to sub-gap states are localized on a length scale in excess of the superconducting coherence length. In the dirty limit, $`\xi \mathrm{}\lambda _F`$, this implies that the localized sub-gap states are quasi-classical in nature. Their existence on the level of the $`\mathrm{\Psi }`$ field action will be buried in the fast $`\lambda _F`$ oscillations of the wavefunction. To reveal the sub-gap states, we must first remove the fast short length scale fluctuations of the quasi-classical Green function and look for an equation of motion for the slowly varying envelope of the wavefunction. But this is just the program of the usual quasi-classical method.

The term “sub-gap states” is a little misleading in this context. Band tails are bound states of some rare potential that sit by themselves below the bulk of the spectrum. Each rare configuration that make the gap soft in the present case will give rise to many states beneath the AG gap. Thus the term “gap fluctuation”, used in Ref.  to describe the zero dimensional SN system, may be more appropriate.

As well as being quasi-classical in nature, the existence of sub-gap states is not affected by working in the dirty limit. As such, their existence must be accommodated in the non-linear $`\sigma `$-model functional (15) since the validity of this description relied only on the quasi-classical parameter $`ϵ_F\tau 1`$ and the dirty limit assumption. To identify sub-gap states in the present formalism, we should therefore investigate inhomogeneous solutions of the low-energy saddle-point equation in $`Q`$ — the Usadel equation . Such a solution should be thought of as defining an envelope for the quasi-classical sub-gap states.

Therefore, let us revisit the mean-field equation and look for inhomogeneous solutions at energies $`ϵ<E_{\mathrm{gap}}`$. To focus our discussion, let us begin by restricting attention to the quasi one-dimensional geometry. To stay firmly within the diffusive regime, we therefore impose the requirement that the system size $`L`$ be much smaller than the localization length of the normal system $`\xi _{\mathrm{loc}.}\nu L_\mathrm{w}D`$, where $`L_\mathrm{w}`$ denotes the cross-section. Later, in section III C, we will generalize our discussion to encompass systems of higher dimension. Furthermore, since, over the interval $`0<\zeta <1`$, the quasi-particle energy gap varies more rapidly than the superconducting order parameter, we will neglect self-consistency of the order parameter. Taking self-consistency into account will not alter our qualitative findings, and will only weakly affect the quantitative results.

### A Instantons in the Quasi One-dimensional Geometry

To investigate inhomogeneous solutions of the mean-field equation (17) it is convenient to recast the equation in terms of its first integral

$$(_{x/\xi }\widehat{\theta })^2+V(\widehat{\theta })=\mathrm{const},$$
(18)
where

$`V(\widehat{\theta })=4i\left(\mathrm{sinh}\widehat{\theta }{\displaystyle \frac{ϵ}{\mathrm{\Delta }}}\mathrm{cosh}\widehat{\theta }\right)\zeta \mathrm{cosh}2\widehat{\theta }`$
denotes the complex potential. Let us denote by $`\theta _{\mathrm{ag}}`$ the values of $`\theta _1`$ and $`i\theta `$ at the conventional saddle point, and focus on an energy $`ϵ`$ below the gap predicted by the AG theory. Here $`\mathrm{Im}\theta _{\mathrm{ag}}=\pi /2`$ such that the mean-field DoS $`\nu _{\mathrm{ag}}(ϵ)=4\nu \mathrm{Re}\mathrm{cosh}\theta _{\mathrm{ag}}`$ vanishes. The corresponding value of $`\mathrm{Re}\theta _{\mathrm{ag}}`$ depends sensitively on the energy, with $`\mathrm{Re}\theta _{\mathrm{ag}}=0`$ for $`ϵ=0`$.

Considering the Boson-Boson sector only, if we require that $`\theta _1(x\pm \mathrm{})=\theta _{\mathrm{AG}}`$, what kind of inhomogeneous solution is possible? The values of $`\theta _1`$ at which $`_x\theta _1=0`$ can be identified by considering the complex (dimensionless) potential function $`V(\theta _1)`$ from which we can determine the endpoints of the ‘motion’ in the complex plane, just as one would use a real potential normally. By inspection one may see that, on the line $`\mathrm{Im}\theta _1=\pi /2`$, the potential is purely real. This is not the only contour where $`\mathrm{Im}V=0`$, but, by considering forces, it is not hard to see that either $`\mathrm{Im}\theta _1=\pi /2`$ always during the motion, or $`\theta _1`$ follows a trajectory with an endpoint at $`\mathrm{Im}\theta _1<0`$. For reasons outlined below, we will discount this latter possibility. The former case amounts to considering “bounce” trajectories in the *real* potential $`V(i\pi /2+\varphi )=V_\mathrm{r}(\varphi )`$ where

$$V_\mathrm{r}(\varphi )4\left(\mathrm{cosh}\varphi \frac{ϵ}{|\mathrm{\Delta }|}\mathrm{sinh}\varphi \right)+\zeta \mathrm{cosh}2\varphi .$$
(19)
A typical potential is shown in Fig. 2.

Now integration over the angles $`\widehat{\theta }`$ is constrained to certain contours . Is the bounce solution accessible to both? As usual, the contour of integration over the Boson-Boson field $`\theta _1`$ includes the entire real axis, while for the Fermion-Fermion field, $`i\theta `$ runs along the imaginary axis from $`0`$ to $`i\pi `$. With a smooth deformation of the integration contours, the AG saddle-point is accessible to both the angles $`\widehat{\theta }`$ . By contrast, the bounce solution and the AG solution can be reached simultaneously by a smooth deformation of the integration contour only for the Boson-Boson field $`\theta _1`$ (see Fig. 3). The bounce solution is therefore associated with a breaking of supersymmetry at the level of the saddle point.

Thus we have identified an inhomogeneous saddle-point configuration for which the supersymmetry is broken: $`\theta _1`$ executes a bounce whilst $`i\theta `$ remains at the mean-field value $`\theta _{\mathrm{ag}}`$. The symmetry broken solution then incurs the (finite) real action

$`S=4\pi \nu L_W(D|\mathrm{\Delta }|)^{1/2}S_\varphi (ϵ/|\mathrm{\Delta }|,\zeta )`$
where, defining $`\varphi ^{}`$ as the endpoint of the motion,

$$S_\varphi _{\varphi _{\mathrm{AG}}}^\varphi ^{}𝑑\varphi \sqrt{V_\mathrm{r}(\varphi _{\mathrm{ag}})V_\mathrm{r}(\varphi )}.$$
(20)
Now, as mentioned above, there exists a second possibility for a bounce solution in which one moves away from $`\theta _{\mathrm{ag}}`$ parallel to the imaginary axes. Indeed, such a solution would seem to be a natural candidate for the Fermion-Fermion field $`i\theta `$. However, since the endpoint for this trajectory lies at $`\mathrm{Re}\theta <0`$ outside the integration domain which runs from $`0`$ to $`\pi `$, this would seem to be excluded.

As $`ϵ`$ approaches $`E_{\mathrm{gap}}`$ from below, the potential (19) becomes more shallow, with the maximum merging with one of the minima when we reach the gap. Near the edge, up to an irrelevant constant, an expansion of the potential in powers of $`(\varphi \varphi _{\mathrm{ag}})`$ leads to the cubic form

$`V_\mathrm{r}[\varphi ]\alpha \left({\displaystyle \frac{E_{\mathrm{gap}}ϵ}{|\mathrm{\Delta }|}}\right)^{1/2}(\varphi \varphi _{\mathrm{ag}})^2`$ (21)
$`+\beta (\varphi \varphi _{\mathrm{ag}})^3`$ (22)
where the dimensionless coefficients are specified by

$`\alpha =6\sqrt{{\displaystyle \frac{2}{3}}}\left({\displaystyle \frac{E_{\mathrm{gap}}}{|\mathrm{\Delta }|}}\right)^{1/6},\beta =2\left({\displaystyle \frac{\zeta E_{\mathrm{gap}}}{|\mathrm{\Delta }|}}\right)^{1/3}.`$ (23)
Note that, making use of Eq. (2), both of these coefficients depend solely on the dimensionless parameter $`\zeta `$. From this expansion, one can obtain an analytic solution for $`S_\varphi `$. To leading order in $`(E_{\mathrm{gap}}ϵ)/|\mathrm{\Delta }|`$ one finds

$$S_\varphi =\frac{4}{15}\frac{\alpha ^{5/2}}{\beta ^2}\left(\frac{E_{\mathrm{gap}}ϵ}{|\mathrm{\Delta }|}\right)^{5/4}.$$
(24)
This approximation is shown in Fig. 4 along with the exact result obtained by numerical integration. Note that the action vanishes exactly at the gap. For completeness we give the explicit form of the bounce solution

$`\varphi (x)\varphi _{\mathrm{ag}}={\displaystyle \frac{\alpha }{\beta }}{\displaystyle \frac{1}{\mathrm{cosh}^2(x/2r_0)}},`$
where the extent of the instanton is set by

$$r_0(ϵ)=\frac{\xi }{\alpha ^{1/2}}\left(\frac{|\mathrm{\Delta }|}{E_{\mathrm{gap}}ϵ}\right)^{1/4}.$$
(25)
Indeed the size of the instanton is easily understood from the quadratic “stiffness” term in Eq. (22). Thus one finds that, while the overall scale is set by the superconducting coherence length $`\xi `$, the size of the droplet diverges both as $`ϵ`$ approaches $`E_{\mathrm{gap}}`$ and, noting that $`\alpha (1\zeta ^{2/3})^{1/4}`$, as one approaches the gapless phase $`\zeta 1`$.

This completes the analysis of the saddle-point solution together with the corresponding action. However, as this level we are presented with two problems:

* the contribution of a second saddle point would seem to spoil the normalization condition $`𝒵[0]_{V,S}=1`$, which should be preserved within the saddle point approximation;
* confined to the line $`\mathrm{Im}\theta =\pi /2`$, when substituted into the DoS source (16), the bounce configuration does not appear to generate states!

The resolution of both problems lies in the nature of the fluctuations around the symmetry broken mean-field solution. These field fluctuations can be separated into “radial” and “angular” contributions. The former involve fluctuations of the diagonal elements $`\widehat{\theta }`$, while the latter describe rotations including those Grassmann transformations which mix the bf sector . Both classes of fluctuations play a crucial role.

### B Fluctuations

Before turning to the technical analysis, let us outline qualitatively the influence of the fluctuations around the mean-field. As usual, associated with radial fluctuations around the bounce, there exists a zero mode (due to translational invariance of the solution), and a negative energy mode. The latter, which necessitates a $`\pi /2`$ rotation of the corresponding integration contour to follow the line of steepest descent (c.f. Ref.  and see below), has two effects: firstly it ensures that the non-perturbative contributions to the local DoS are non-vanishing, and secondly, that they are positive. Turning to the angular fluctuations, the breaking of supersymmetry is accompanied by the appearance of a Grassmann zero mode separated by a gap from higher excitations which restores the global supersymmetry (c.f. spin symmetry breaking in a ferromagnet of finite extent). The zero mode ensures that the symmetry broken inhomogeneous saddle-point configurations respect the normalization condition $`𝒵[0]_{V,S}=1`$.

To formally investigate the fluctuation determinant, let us implement the rational parameterization

$$Q=R\sigma _3^{\mathrm{ph}}\sigma _3^{\mathrm{cc}}\frac{1+iP}{1iP}R^1,$$
(26)
where the condition $`P=\sigma _2^{\mathrm{ph}}\mathrm{\Gamma }P^T\mathrm{\Gamma }^1\sigma _2^{\mathrm{ph}}`$ with

$`\mathrm{\Gamma }\sigma _3^{\mathrm{cc}}\gamma =E_{\mathrm{bb}}\sigma _1^{\mathrm{cc}}iE_{\mathrm{ff}}\sigma _2^{\mathrm{cc}},`$
is imposed by (14). With this choice, a variant of that used in Ref. , the measure is trivial . $`R`$ rotates $`Q`$ from the metallic saddle point $`\sigma _3^{\mathrm{ph}}\sigma _3^{\mathrm{cc}}`$ to the bounce configuration:

$`R(x)=\mathrm{exp}\left[{\displaystyle \frac{1}{2}}\sigma _2^{\mathrm{ph}}\sigma _3^{\mathrm{cc}}\widehat{\theta }(x)\right].`$
Defining

$`P=\left(\begin{array}{cc}C& A\\ B& \mathrm{\Gamma }C^T\mathrm{\Gamma }^1\end{array}\right)_{\mathrm{ph}},`$
where $`A=\mathrm{\Gamma }A^T\mathrm{\Gamma }^1`$ and $`B=\mathrm{\Gamma }B^T\mathrm{\Gamma }^1`$, the condition $`Q^2=𝟙`$ requires $`[\sigma _3^{\mathrm{ph}}\sigma _3^{\mathrm{cc}},P]_+=0`$, or

$`[\sigma _3^{\mathrm{cc}},C]_+=[\sigma _3^{\mathrm{cc}},A]_{}=[\sigma _3^{\mathrm{cc}},B]_{}=0.`$
Thus the matrix field $`C`$ is off-diagonal in the cc space, while $`A`$ and $`B`$ are diagonal. In fact, the field fluctuations $`C`$ describe the low-energy quantum interference effects — the soft Class D modes. However, the fluctuations contained within the fields $`C`$ are oblivious to the supersymmetry breaking. We will therefore deal only with the fluctuations that are parameterized by $`A`$ and $`B`$. (Moreover, we will neglect the massive spin triplet fluctuations.)

With the explicit parameterization

$`A`$ $`=`$ $`\left(\begin{array}{cccc}s_\mathrm{b}& 0& \overline{\eta }_A& 0\\ 0& s_\mathrm{b}& 0& \eta _A\\ \eta _A& 0& s_\mathrm{f}& 0\\ 0& \overline{\eta }_A& 0& s_\mathrm{f}\end{array}\right)`$
$`B`$ $`=`$ $`\left(\begin{array}{cccc}s_\mathrm{b}^{}& 0& \overline{\eta }_B& 0\\ 0& s_\mathrm{b}^{}& 0& \eta _B\\ \eta _B& 0& s_\mathrm{f}^{}& 0\\ 0& \overline{\eta }_B& 0& s_\mathrm{f}^{}\end{array}\right),`$
an expansion of the action (15) to quadratic order around the bounce configuration

$`\widehat{\theta }(x)=\left(\begin{array}{cc}\theta _1(x)& 0\\ 0& \theta _{\mathrm{ag}}\end{array}\right)_{\mathrm{bf}}.`$
obtains

$`S=4\pi \nu L_W(D|\mathrm{\Delta }|)^{1/2}(S_\varphi +S_\mathrm{q})`$
where $`S_\varphi `$ represents the contribution from the saddle-point alone (20), and

$`S_\mathrm{q}`$ $`=`$ $`{\displaystyle 𝑑u\left[(|s_\mathrm{b}|^2+|s_\mathrm{f}|^2)+V^{}(\varphi )s_{\mathrm{b}}^{}{}_{}{}^{2}+V^{\prime \prime }(\varphi )s_{\mathrm{b}}^{\prime \prime }{}_{}{}^{2}+V^{}(\varphi _{\mathrm{ag}})s_{\mathrm{f}}^{}{}_{}{}^{2}+V^{\prime \prime }(\varphi _{\mathrm{ag}})s_{\mathrm{f}}^{\prime \prime }{}_{}{}^{2}\right]}`$ (28)
$`+{\displaystyle 𝑑u\left[(\overline{\xi }_+\xi _++\overline{\xi }_{}\xi _{})+V_+\overline{\xi }_+\xi _++V_{}\overline{\xi }_{}\xi _{}\right]},`$
where $`s_{\mathrm{b}/\mathrm{f}}s_{\mathrm{b}/\mathrm{f}}^{}+is_{\mathrm{b}/\mathrm{f}}^{\prime \prime }`$, $`\xi _\pm (\eta _B\pm \eta _A)/\sqrt{2}`$, $`\overline{\xi }_\pm (\overline{\eta }_B\pm \overline{\eta }_A)/\sqrt{2}`$, $`ux/\xi `$ and $`\theta _{\mathrm{ag}}=i\pi /2+\varphi _{\mathrm{ag}}`$. Here the various potentials are given by

$`V^{}(\varphi )`$ $`=`$ $`2\left(\mathrm{cosh}\varphi {\displaystyle \frac{ϵ}{|\mathrm{\Delta }|}}\mathrm{sinh}\varphi \zeta \mathrm{cosh}2\varphi \right)`$
$`V^{\prime \prime }(\varphi )`$ $`=`$ $`(\varphi )^2+2\left(\mathrm{cosh}\varphi {\displaystyle \frac{ϵ}{|\mathrm{\Delta }|}}\mathrm{sinh}\varphi \zeta \mathrm{sinh}^2\varphi \right),`$
and

$`V_\pm (\varphi )={\displaystyle \frac{(\varphi )^2}{4}}+\mathrm{cosh}\varphi +\mathrm{cosh}\varphi _{\mathrm{ag}}`$
$`{\displaystyle \frac{ϵ}{|\mathrm{\Delta }|}}(\mathrm{sinh}\varphi +\mathrm{sinh}\varphi _{\mathrm{ag}}){\displaystyle \frac{\zeta }{2}}\left(\mathrm{cosh}2\varphi +\mathrm{cosh}2\varphi _{\mathrm{ag}}\right)`$
$`\zeta \left(\mathrm{sinh}\varphi \mathrm{sinh}\varphi _{\mathrm{AG}}\mathrm{cosh}\varphi \mathrm{cosh}\varphi _{\mathrm{AG}}\right).`$
As a check, let us consider the conventional saddle point. At $`ϵ=0`$ (for simplicity), $`\varphi =\varphi _{\mathrm{ag}}=0`$, and the quadratic action assumes the form

$`S_\mathrm{q}={\displaystyle 𝑑u\left(|s_\mathrm{b}|^2+|s_\mathrm{f}|^2+\overline{\xi }_+\xi _++\overline{\xi }_{}\xi _{}\right)}`$
$`+2{\displaystyle }du[|s_\mathrm{b}|^2+|s_\mathrm{f}|^2+\overline{\xi }_+\xi _++\overline{\xi }_{}\xi _{}`$
$`\zeta (s_{\mathrm{b}}^{}{}_{}{}^{2}+s_{\mathrm{f}}^{}{}_{}{}^{2}+\overline{\xi }_{}\xi _{})].`$
This action is manifestly supersymmetric and, therefore, performing the integrations over all fields gives unity. (Moreover, since $`\zeta <1`$ in the gapped phase, the integral is manifestly convergent.) This in why the usual AG DoS is just given by evaluating the source (16) at the supersymmetric AG saddle point.

Since $`V^{}(\varphi _{\mathrm{ag}})`$ and $`V^{\prime \prime }(\varphi _{\mathrm{ag}})`$ are both positive definite, integration over the Fermion-Fermion degrees of freedom merely generates some (weakly energy dependent) positive prefactor. The Boson-Boson sector is more interesting. In particular, it is straightforward to verify that the action for the components $`s_\mathrm{b}^{}`$ simply reflects the longitudinal variation

$`{\displaystyle \frac{1}{2}}{\displaystyle 𝑑x𝑑x^{}\mathrm{\Delta }\theta _1(x)\frac{\delta ^2S}{\delta \theta _1(x)\delta \theta _1(x^{})}\mathrm{\Delta }\theta _1(x^{})}`$
with $`2s_\mathrm{b}^{}=\mathrm{\Delta }\theta _1`$. Thus the action for $`s_\mathrm{b}^{}`$ is the one that could have been written down from the outset: it is that of the fluctuations of $`\theta _1`$ around the bounce instanton discussed in section III A. To address the influence of this class of fluctuations we can draw on the standard literature .

To perform the Gaussian integration over the fields $`s_\mathrm{b}^{}`$, we form the expansion

$`s_\mathrm{b}^{}(x)={\displaystyle \underset{\mathrm{n}}{}}a_\mathrm{n}\phi _\mathrm{n}(x)`$
in terms of the eigenfunctions $`\phi _\mathrm{n}`$ of the quadratic operator in the action for $`s_\mathrm{b}^{}`$. Now the action for $`s_\mathrm{b}^{}`$ exhibits a zero mode $`\phi _1\theta _1`$ due to translational invariance of the action (the bounce can be positioned anywhere in space). This zero mode must be accommodated by the introduction of a collective coordinate which in turn introduces a Jacobian factor associated with the change of variables from $`a_1`$ to the collective coordinate. Furthermore, since the instanton is a bounce, the zero mode $`\phi _1`$ has a node, and hence there exists a ground state of $`V^{}`$ with negative energy. This requires the contour for the integration variable $`a_0`$ to be deformed away from the real axis — it takes a right turn at zero and heads in the negative imaginary direction. In the context of the contour drawn in Fig. 3, this is because we have to get back to the real axis for $`\theta _1`$.

This deformation of the integration contour has a profound consequence. While the contribution to the sub-gap states from the bounce solution alone vanishes (recall that the instanton was confined to the line $`\mathrm{Im}\theta _1=\pi /2`$, so that $`\mathrm{cosh}\theta _1`$ remains imaginary), the rotation of the integration contour introduces a factor of $`i`$ resulting in an imaginary contribution to the Green function below the gap. This is turn would signify a non-zero sub-gap DoS. The mechanism operates in the related context of Landau levels broadened by disorder and discussed by Efetov and Marikhin  (see also an earlier paper by Affleck ).

It is reassuring to note that the deformation of the integration contour, which is constrained to return to the undeformed Bosonic contour, is unambiguous, and gives rise to only positive definite contributions to the DoS. Finally, as usual, since the contour only runs over half a Gaussian for $`a_0`$, its contribution yields a factor of $`1/2`$.

This completes the discussion of the contribution to the functional integral from the Bosonic degrees of freedom. Finally we turn to the role of the Grassmann fluctuations. Considering the supersymmetric structure of the action, it is immediately clear that supersymmetry breaking must be accompanied by the existence of zero mode in the Grassmann sector which restores the supersymmetry. Indeed, an inspection of the potential $`V_{}`$ identifies a zero energy eigenfunction gapped from the others (see Fig. 5.

This situation may be compared with the analysis of Andreev and Altshuler , who identified a stationary phase saddle-point that determines the oscillatory part of energy level correlation functions for normal diffusive conductors. Again this saddle-point is more correctly a manifold of points related by supersymmetry, but the fact that it is spatially uniform guarantees that the spectrum of Grassmann fluctuations is truly gapless.

This zero mode indicates the existence of a degenerate manifold parameterized by a Grassmann coordinate. The integration of unity over this Grassmann coordinate yields zero, which ensures the necessary normalization condition $`𝒵[0]_{V,S}=1`$ is met. A non-zero result can be obtained only if there is a prefactor (or source term) that breaks supersymmetry. Making use of the DoS source (16), we expand in $`P`$ to find the lowest order part quadratic in the Grassmanns (at order $`P^2`$)

$`\mathrm{str}\left(\sigma _3^{\mathrm{bf}}\sigma _3^{\mathrm{ph}}\sigma _3^{\mathrm{cc}}Q(x)\right)`$
$`8i(\mathrm{sinh}\varphi (x)\mathrm{sinh}\varphi _{\mathrm{ag}})\times \left(\overline{\xi }_+(x)\xi _+(x)+\overline{\xi }_{}(x)\xi _{}(x)\right).`$
From this result, we see that the local sub-gap DoS remains non-zero only in the vicinity of the instanton where the supersymmetry is broken.

Thus, taking into account Gaussian fluctuations and zero modes, one obtains the non-perturbative, one instanton contribution to the sub-gap DoS:

$`{\displaystyle \frac{\nu (ϵ)_{V,S}}{4\nu }}(i|K|){\displaystyle 𝑑xi(\mathrm{sinh}\varphi (x)\mathrm{sinh}\varphi _{\mathrm{ag}})|\phi _0^{}(x)|^2\sqrt{\frac{LS_\varphi }{\xi }}\mathrm{exp}\left[4\pi \nu L_\mathrm{w}\sqrt{D|\mathrm{\Delta }|}S_\varphi \right]},`$ (29)
where the factor $`\sqrt{LS_\varphi /\xi }`$ represents the Jacobian associated with the introduction of the collective coordinate , $`i|K|`$ is the overall factor arising from the non-zero modes, and the Grassmann zero mode wavefunction $`\phi _0^{}`$ is normalized such that

$`{\displaystyle 𝑑x|\phi _0^{}|^2}=1.`$
Here we have assumed that the $`s_\mathrm{b}^{\prime \prime }`$ integration only contributes to the positive prefactor. We have checked this numerically for a few cases and expect a general statement could be made by moving to a different parameterization.

Eq. (29) is the main result of this section. Note the non-perturbative nature of the result, both in the coupling constant $`g^1`$ of the $`\sigma `$-model, and the (dimensionless) spin scattering rate $`\zeta `$.

### C Sub-gap States in Dimensions $`d>1`$

The calculation above was tailored to the consideration of the quasi one-dimensional geometry. The generalization to higher dimensions follows straightforwardly. In particular, it is necessary to seek inhomogeneous solutions of the saddle-point equation (18) where the gradient operator must be interpreted as the higher dimensional generalization. Generally, this equation must be solved numerically. However, for energies $`ϵ`$ in the vicinity of the gap $`E_{\mathrm{gap}}`$, an analytic expression for the energy scaling can be obtained.

Using the approximation to $`V_\mathrm{r}[\varphi ]`$ (22) valid when $`(E_{\mathrm{gap}}ϵ)/|\mathrm{\Delta }|1`$, the exponential dependence of the sub-gap DoS can be deduced in higher dimension. In this limit, dimensional analysis of the cubic equation of motion yields the scaling form

$`\varphi (𝐫)\varphi _{\mathrm{ag}}(ϵ)={\displaystyle \frac{\alpha }{\beta }}f(𝐫/r_0),`$
where $`r_0`$ is the characteristic length defined by Eq. (25). When substituted back into the action, one finds that the DoS depends exponentially on the parameter $`4\pi g(\xi /L)^{d2}S_\varphi `$ where

$`S_\varphi =a_d\zeta ^{2/3}(1\zeta ^{2/3})^{(2+d)/8}\left({\displaystyle \frac{E_{\mathrm{gap}}ϵ}{|\mathrm{\Delta }|}}\right)^{(6d)/4}.`$ (30)
Here $`g=\nu DL^{d2}`$ denotes the bare dimensionless conductance of the normal system, and $`a_d`$ is a numerical constant ($`a_1=8{}_{}{}^{4}\sqrt{24}/5`$) In particular, the exponent depends linearly on the energy separation from the gap in two dimensions.

### D Numerics

To assess the validity of the approximations used in obtaining the results above we have investigated numerically the DoS in the vicinity of the mean-field energy gap using a numerical diagonalisation of a non-interacting tight-binding Bogoluibov Hamiltonian with disorder in both the on-site potential matrix elements and in the spin impurity scattering potential. Taking a two-dimensional lattice of size $`22\times 22`$ site with an on-site disorder taken from the range of $`ϵ[3,3]`$ measured in units of the hopping matrix element, $`|\mathrm{\Delta }|=1`$ measured in the same units, and various strengths of the spin impurity potential, the DoS is shown in Fig. 6. Notice that, as the strength of the spin impurity potential is increased, the quasi-particle energy gap is quenched.

For the weakest magnetic impurity potential, we have expanded the region in the vicinity of the mean-field energy gap. Fig. 7 shows a exponential scaling of the ‘sub-gap’ DoS with an exponent which depends linearly on energy as predicted by the theory above.

## IV Zero dimensional problems and universality

In the previous section we considered the instanton contribution to the sub-gap DoS in the infinite system. For completeness let us now consider the zero dimensional case that obtains when $`r_0`$, the size of the instanton, exceeds the system size $`L`$, which will happen when $`ϵ`$ approaches close enough to $`E_{\mathrm{gap}}`$ from below in any finite system. In this limit one can clearly not fit an instanton inside the system. Leaving aside the practical relevance of this situation, theoretical motivation is provided by a recent paper  that explored this regime using a random matrix analysis.

However, before turning to the consideration of the zero-dimensional limit of the present problem let us first try to draw some intuition from a closely related investigation of a different system comprised of a normal quantum dot contacted to a superconducting terminal (as shown in Fig. 8). In such a geometry it is well established (see e.g. Ref. ) that the proximity effect induces a gap in the DoS of the normal dot. Indeed, in Ref. , the integrity of the gap is proposed as a signature of irregular or chaotic dynamics inside the dot. Now near the gap edge the DoS of the dot takes the singular form

$$\nu (ϵ>E_{\mathrm{gap}})\frac{1}{\pi L^d}\sqrt{\frac{ϵE_{\mathrm{gap}}}{\mathrm{\Delta }_\mathrm{g}^3}},$$
(31)
with $`E_{\mathrm{gap}}=cN\delta `$ and $`\mathrm{\Delta }_\mathrm{g}=c^{}N^{1/3}\delta `$, where $`c0.048`$ and $`c^{}0.068`$, $`\delta =1/L^d\nu `$ denotes the single particle level spacing, and $`N`$ is the number of fully transmitting channels in the contact.

However, the location of the gap edge relies on a mean-field analysis of the coupled system. In Ref. Vavilov et al. have argued that optimal fluctuations of the impurity potential give rise to gap fluctuations. The hypothesis introduced in Ref.  is that the spectral statistics near a gap edge are universal. This allows a random matrix theory analysis of gap fluctuations and leads to the following expression for the sub-gap DoS,

$`{\displaystyle \frac{\nu (ϵ<E_{\mathrm{gap}})}{\nu }}\mathrm{exp}\left[{\displaystyle \frac{2}{3}}\left({\displaystyle \frac{E_{\mathrm{gap}}ϵ}{\mathrm{\Delta }_\mathrm{g}}}\right)^{3/2}\right].`$ (32)
Now the AG mean-field solution for a superconducting quantum dot with magnetic impurities also predicts the existence of a square root edge (see Eq. (6) and inset of Fig. 4). Then, when recast in the form of Eq. (31), it is pertinent to ask whether the expression for the sub-gap DoS coincides with Eq. (32) in the zero dimensional limit. This is the situation addressed in Ref..

In the following we will show that the universal result (32) is explicitly recovered by the present theory. Moreover, in doing so, we will expose the origin of the universal structure reported in Ref.  and describe its implications for universality of the $`d>0`$ problem.

### A Superconducting dot with magnetic impurities

Let us therefore consider explicitly the action of a superconducting grain in the presence of a weak magnetic impurity potential. When $`r_0L`$ only the zero spatial mode contributes significantly to the action (15). In this limit, the action assumes the zero dimensional form

$`S[Q]={\displaystyle \frac{i\pi }{2\delta }}\mathrm{str}\left[\left(ϵ_{}\sigma _3^{\mathrm{cc}}+|\mathrm{\Delta }|\sigma _2^{\mathrm{ph}}\right)\sigma _3^{\mathrm{ph}}Q\right]`$ (33)
$`+{\displaystyle \frac{\pi }{24\tau _s\delta }}\mathrm{str}\left(Q\sigma _3^{\mathrm{ph}}\sigma ^{\mathrm{sp}}\right)^2,`$ (34)
where, as usual, $`\delta `$ denotes the single-particle level-spacing. As in the higher dimensional problem, a variation of the action with respect to $`Q`$ obtains a mean-field equation, now without spatial variation. Parameterizing the saddle-point equation as in section II D, we obtain the zero-dimensional Usadel or AG mean-field equation (17)

$`2i\left(\mathrm{cosh}\widehat{\theta }{\displaystyle \frac{ϵ}{|\mathrm{\Delta }|}}\mathrm{sinh}\widehat{\theta }\right)\zeta \mathrm{sinh}(2\widehat{\theta })=0.`$
From this equation, we can identify the usual AG solution which in turn recovers the AG phenomenology.

The inclusion of bounce configurations in the previous calculations was based upon the observation that, although the contribution they make is exponentially small, they are the least action configurations on the part of the contour that gives a finite sub-gap DoS. In the zero-dimensional case we are spared having to think about the problem in function space. The action is proportional to the potential of Fig. 2. The correct contour thus passes through the maximum of the potential (minimum of the action) corresponding to the usual AG saddle point, and turns away from the real $`\varphi `$ axes (i.e. the line $`\mathrm{Im}\theta _1=\pi /2`$) at the minimum of the potential

$`\varphi _{0\mathrm{d}}(ϵ)=\varphi _{\mathrm{ag}}(ϵ)+{\displaystyle \frac{2\alpha }{3\beta }}\sqrt{{\displaystyle \frac{E_{\mathrm{gap}}ϵ}{|\mathrm{\Delta }|}}},`$
(marked in Fig. 2). This part of the contour, parallel to the imaginary axes, gives a contribution to the DoS, and the second saddle point is in fact a *maximum* on this portion by analyticity. Following the same arguments as in section III A, this solution is inaccessible to the Fermionic contour. We find that the sub-gap DoS near $`E_{\mathrm{gap}}`$ in the zero-dimensional case scales as

$`{\displaystyle \frac{\nu (ϵ<E_{\mathrm{gap}})}{\nu }}\mathrm{exp}\left[{\displaystyle \frac{4\pi }{27}}{\displaystyle \frac{|\mathrm{\Delta }|}{\delta }}{\displaystyle \frac{\alpha ^3}{\beta ^2}}\left({\displaystyle \frac{E_{\mathrm{gap}}ϵ}{|\mathrm{\Delta }|}}\right)^{3/2}\right],`$
where $`\alpha `$ and $`\beta `$ are the coefficients defined in Eq. (23). We note that the general result for the energy dependence of the exponent written down in section III C for dimensions $`d1`$ applies also for $`d=0`$.

To establish contact with the universal result given in Eq. (31) it is helpful to recast the result in a modified form. To do this we note that, in the vicinity of the mean-field gap edge, the DoS can be expanded as

$`{\displaystyle \frac{\nu _{\mathrm{ag}}(ϵ>E_{\mathrm{gap}})}{4\nu }}`$
$`\sqrt{{\displaystyle \frac{2}{3}}}\zeta ^{2/3}(1\zeta ^{2/3})^{1/4}\left({\displaystyle \frac{ϵE_{\mathrm{gap}}}{|\mathrm{\Delta }|}}\right)^{1/2}.`$
Then, if we define

$`\mathrm{\Delta }_\mathrm{g}^{3/2}{\displaystyle \frac{4\pi }{\delta }}\sqrt{{\displaystyle \frac{2}{3|\mathrm{\Delta }|}}}\zeta ^{2/3}(1\zeta ^{2/3})^{1/4},`$
the mean-field DoS can be brought to the form of Eq. (31), and the sub-gap DoS takes the universal form

$`{\displaystyle \frac{\nu (ϵ<E_{\mathrm{gap}})}{\nu }}\mathrm{exp}\left[{\displaystyle \frac{4}{3}}\left({\displaystyle \frac{E_{\mathrm{gap}}ϵ}{\mathrm{\Delta }_\mathrm{g}}}\right)^{3/2}\right].`$ (35)
Leaving aside a spurious (yet systematic — see below) factor of $`2`$ , the sub-gap DoS obtained above coincides with the universal expression shown in Eq. (32).

Once again, to assess the validity of the approximation scheme that leads to the universal result above, we have looked for the scaling in a numerical investigation of a random matrix Hamiltonian with the same symmetry. In the sub-gap region, Fig. 9 shows a good fit of the DoS to the predicted exponential scaling.

The rescaling of the DoS above and the appearance of the universal form suggests that we should revisit the $`d`$-dimensional result and look for a similar rescaling. From Eq. (30) it is straightforward to verify that in this case

$`{\displaystyle \frac{\nu (ϵ<E_{\mathrm{gap}})}{\nu }}\mathrm{exp}\left[\stackrel{~}{a}_d\left({\displaystyle \frac{r_0(ϵ)}{L}}\right)^d\left({\displaystyle \frac{E_{\mathrm{gap}}ϵ}{\mathrm{\Delta }_\mathrm{g}}}\right)^{3/2}\right],`$
where $`\stackrel{~}{a}_d`$ represents some numerical constant, and $`r_0`$ is the characteristic length defined by Eq. (25). Finally, by defining

$`\stackrel{~}{\mathrm{\Delta }}_\mathrm{g}^{3/2}(ϵ){\displaystyle \frac{\delta }{\stackrel{~}{\delta }(ϵ)}}\mathrm{\Delta }_g^{3/2},`$
where $`\stackrel{~}{\delta }(ϵ)=1/(\nu r_0^d(ϵ))`$ is the level spacing inside a region of size $`r_0`$, the volume dependent prefactor can be absorbed into the expression and DoS can be brought to the form

$`{\displaystyle \frac{\nu (ϵ<E_{\mathrm{gap}})}{\nu }}\mathrm{exp}\left[\stackrel{~}{a}_d\left({\displaystyle \frac{E_{\mathrm{gap}}ϵ}{\stackrel{~}{\mathrm{\Delta }}_\mathrm{g}}}\right)^{3/2}\right],`$
revealing a simple relation between the $`d=0`$ and $`d>0`$ problems.

The coincidence of Eqs. (32) and (35) indeed suggests that gap fluctuations are universal. To expose the origin of the universal scaling within the present formalism, let us consider the quantum dot geometry of Fig. 8 within the $`\sigma `$-model scheme.

### B Quantum dot contacted to a superconductor

Following Ref. , let us consider a normal quantum dot contacted to a bulk superconducting terminal with order parameter $`\mathrm{\Delta }`$. Taking $`N`$ fully open channels between the dot and the lead, it is straightforward to show that, in the zero-dimensional limit, the effective action of the hybrid SN system takes the form

$`S[Q]=i{\displaystyle \frac{\pi ϵ_+}{2\delta }}\mathrm{str}\left[\sigma _3^{\mathrm{cc}}\sigma _3^{\mathrm{ph}}Q\right]`$ (36)
$`{\displaystyle \frac{N}{2}}\mathrm{str}\left[\mathrm{ln}(1+Q_\mathrm{L}Q)\right],`$ (37)
where $`Q_\mathrm{L}`$ represents the $`Q`$-matrix of the superconducting terminal (see below), and $`Q`$ represents the zero-dimensional supermatrix field within the normal region. The term containing the logarithm describes the coupling to the lead — see Ref.  for a derivation. Even for a diffusive dot, we note that retaining the full logarithm is the essential element in the correct treatment of the zero-dimensional limit that obtains when $`D/L^2N\delta `$, a point to which we will return later. For a clean superconducting lead, at energies $`ϵ`$ much less than the bulk superconducting gap $`\mathrm{\Delta }`$, one can set $`Q_\mathrm{L}=\sigma _1^{\mathrm{ph}}`$.

As usual, to obtain the mean-field expression for the DoS it is necessary to minimize the action with respect to variations in $`Q`$. Doing so, one obtains the saddle-point equation

$`i{\displaystyle \frac{\pi ϵ_+}{2\delta }}[Q,\sigma _3^{\mathrm{cc}}\sigma _3^{\mathrm{ph}}]+{\displaystyle \frac{N}{2}}[Q,(1+Q_LQ)^1Q_L]=0`$
As usual, applying the Ansatz that the saddle-point solution is contained within the diagonal parameterization,

$`Q=\sigma _3^{\mathrm{cc}}\sigma _3^{\mathrm{ph}}\mathrm{cosh}\widehat{\theta }+i\sigma _1^{\mathrm{ph}}\mathrm{sinh}\widehat{\theta },`$
the saddle-point equation takes the form

$$\frac{\pi ϵ}{\delta }\mathrm{sinh}\widehat{\theta }+\frac{N}{2}\frac{\mathrm{cosh}\widehat{\theta }}{1+i\mathrm{sinh}\widehat{\theta }}=0.$$
(38)
A straightforward analysis of the symmetric saddle-point solution leads to the mean-field result for the DoS shown in Eq. (31) (remembering that now $`\nu (ϵ)=2\nu \mathrm{Re}\mathrm{cosh}\theta _{\mathrm{mf}}(ϵ)`$ as there are no spin degrees of freedom in this problem). In particular, one can straightforwardly determine $`E_{\mathrm{gap}}`$ by setting $`\mathrm{cosh}\theta _{\mathrm{mf}}`$ to be imaginary. Thus $`\mathrm{sinh}\theta _{\mathrm{mf}}ib`$ for real $`b`$ and (38) gives

$`ϵ(b)={\displaystyle \frac{N\delta }{2\pi }}{\displaystyle \frac{1}{b}}\sqrt{{\displaystyle \frac{b1}{b+1}}}.`$
The extremum of this function gives the largest energy corresponding to a real value of $`b`$. This occurs at $`b=(1+\sqrt{5})/2=1+\gamma `$, where $`\gamma `$ is the golden mean, and yields $`E_{\mathrm{gap}}=(N\delta /2\pi )\gamma ^{5/2}0.048N\delta `$ as required.

As with the case of magnetic impurities, to explore the influence of gap fluctuations, it is necessary to seek the symmetry broken saddle-point configuration. However, for reasons outlined below, it is possible to identify a pattern which implies the universality of the resulting analysis:

### C Universalities

At first glance the situations considered in the two proceeding sections are rather different. The actions (34) and (37) would seem not to have much in common. However, a simple and general argument may be established to reveal the universal character. As before, defining $`\theta _{\mathrm{mf}}(ϵ)=i\pi /2+\varphi _{\mathrm{mf}}(ϵ)`$, the mean-field DoS for the SN device is given by $`\nu _{\mathrm{mf}}(ϵ)=2\nu \mathrm{Im}\mathrm{sinh}\varphi _{\mathrm{mf}}(ϵ)`$, where $`\varphi _{\mathrm{mf}}`$ is determined by the condition $`\delta S/\delta \varphi [\varphi =\varphi _{\mathrm{mf}}]=0`$. Since the DoS displays a square root singularity described by Eq. (31), the (saddle-point) action near the edge is constrained to be of the form

$`S[\widehat{\varphi }]=k\mathrm{str}\left[{\displaystyle \frac{1}{3}}\widehat{s}^3+\left({\displaystyle \frac{\delta }{2\pi }}\right)^2\left({\displaystyle \frac{ϵ_+E_{\mathrm{gap}}}{\mathrm{\Delta }_\mathrm{g}^3}}\right)\widehat{s}\right],`$
where $`\widehat{s}(ϵ)=\mathrm{sinh}\widehat{\varphi }(ϵ)\mathrm{sinh}\widehat{\varphi }_{\mathrm{mf}}(E_{\mathrm{gap}})`$. Here, the elements $`\widehat{\varphi }=\mathrm{diag}(\varphi _{\mathrm{bb}},\varphi _{\mathrm{ff}})`$ and $`\widehat{s}=\mathrm{diag}(s_{\mathrm{bb}},s_{\mathrm{ff}})`$ are diagonal in the superspace. (As one may check, a variation of the action for $`ϵ>E_{\mathrm{gap}}`$ obtains the symmetric mean-field solution

$`s_{\mathrm{bb}}=s_{\mathrm{ff}}=i{\displaystyle \frac{\delta }{2\pi }}\sqrt{{\displaystyle \frac{ϵE_{\mathrm{gap}}}{\mathrm{\Delta }_g^3}}}`$
which in turn recovers the expression (31) for $`\nu (ϵ)`$.) Moreover, since the term containing $`\widehat{s}`$ is linear in the energy, we can determine the value of $`k`$ from the knowledge that $`ϵ`$ appears in the action as $`(2\pi ϵ_+/\delta )\mathrm{sinh}\widehat{\varphi }`$. (It is this term that can more generally contain the Dyson index ‘$`\beta `$’, which therefore appears in the general expression for gap fluctuations described in Ref. .) In the present case, we thus have $`k=(2\pi \mathrm{\Delta }_\mathrm{g}/\delta )^3`$.

Now, as discussed in the previous section, when $`ϵ<E_{\mathrm{gap}}`$ there exists two saddle-point solutions at

$`s_\pm =\pm {\displaystyle \frac{\delta }{2\pi }}\sqrt{{\displaystyle \frac{E_{\mathrm{gap}}ϵ}{\mathrm{\Delta }_g^3}}}.`$
As before, one of these solutions ($`s_{}(ϵ)\varphi _{\mathrm{mf}}(ϵ)`$) is associated with the conventional symmetric mean-field solution while the other represents a second saddle-point accessible only to the Bosonic contour. Taking this second, symmetry broken saddle-point into account (i.e. setting $`s_{\mathrm{bb}}=s_+`$ and $`s_{\mathrm{ff}}=s_{}`$), one obtains the saddle-point action

$`S[\widehat{\varphi }]={\displaystyle \frac{4}{3}}\left({\displaystyle \frac{E_{\mathrm{gap}}ϵ}{\mathrm{\Delta }_\mathrm{g}}}\right)^{3/2}.`$
It is this symmetry broken saddle-point which controls the sub-gap DoS and leads to the universal scaling form proposed in Ref. . This generalizes the arguments applied to the superconducting dot with magnetic impurities.

### D Discussion

Following on from this discussion, to conclude this section, let us make two remarks which bare on the universality of the general scheme. The first of these remarks concerns the integrity of the scaling of the sub-gap DoS when different impurity distributions are taken into account. The second remark concerns the extension of the ideas above to the consideration of the hybrid SN system beyond the zero-dimensional regime.

Firstly, for the superconductor with magnetic impurities, one can generalize the arguments above to show that the energy scaling of the sub-gap DoS even in the $`d`$-dimensional case is insensitive to the nature of the random impurity distribution. This is in contrast to Lifshitz band tail states in semi-conductors where the energy scaling depends sensitively on this distribution. To understand this, let us suppose that the distribution of magnetic impurities $`J𝐒(𝐫)`$ is not Gaussian $`\delta `$-correlated, as we assumed throughout, but obeys some arbitrary statistics defined by a probability functional $`P[J𝐒(𝐫)]`$. When the ensemble average over $`J𝐒(𝐫)`$ is performed one would obtain in the $`\mathrm{\Psi }`$-field action a contribution of the form

$`\mathrm{ln}\mathrm{exp}[i{\displaystyle }d𝐫\overline{\mathrm{\Psi }}J𝐒\sigma ^{\mathrm{sp}}\mathrm{\Psi }]_PC[\overline{\mathrm{\Psi }}\sigma ^{\mathrm{sp}}\mathrm{\Psi }(𝐫)],`$
which defines $`C[\mathrm{}]`$, the generating functional of connected correlators of $`J𝐒(𝐫)`$. Though this is in general a very complicated and indeed non-local functional of $`\overline{\mathrm{\Psi }}\sigma \mathrm{\Psi }(𝐫)`$, one can in principle find a *local* $`Q`$-field action by including pairings only at coincident points, justified by the assumption $`(\mathrm{}/\xi )^d1`$ about the non-magnetic disorder. The mean-field description of this system then follows from the homogeneous solution of the saddle-point equation, an Usadel equation like (17) with some potential. Generally this potential will have the same characteristics as the real potential of (19) plotted in Fig. 2 on the line $`\mathrm{Im}\theta _1=\pi /2`$. The central maximum is due to the $`|\mathrm{\Delta }|`$ term; the upturn at large $`\varphi `$ arises from the small pair-breaking part, and the asymmetry comes from the $`ϵ`$ term. Now, if mean-field theory leads to a square-root singularity in the DoS (a circumstance which can be avoided only by a special tuning of parameters), one can expect that increasing the energy leads to the maximum merging with one of the minima according to

$`V_\mathrm{r}[\varphi ]\alpha \left({\displaystyle \frac{E_{\mathrm{gap}}ϵ}{|\mathrm{\Delta }|}}\right)^{1/2}(\varphi \varphi _{\mathrm{mf}})^2+\beta (\varphi \varphi _{\mathrm{mf}})^3`$
with $`\alpha `$ and $`\beta `$ chosen appropriately. Then the analysis of section III applies. In particular the scaling of the exponent with $`((E_{\mathrm{gap}}ϵ)/|\mathrm{\Delta }|)^{(6d)/4}`$ is expected to be *universal and independent of the details of the magnetic impurity potential*.

Now let us turn to the generality of the present scheme in describing ‘gap fluctuations’ in extended hybrid superconductor/normal systems. The latter has been discussed in a very recent paper by Ostrovsky et al. . In this work, the authors developed an instanton approach analogous to that employed in the magnetic impurity system here to estimate the profile of gap fluctuations in the $`d`$-dimensional SNS system. Now from the discussion above, it is possible to expose the relation between these two works: in the SNS system, the energy gap induced in the normal region due to the proximity effect is determined by the Thouless energy defined as $`E_\mathrm{T}1/\tau _{\mathrm{dwell}}`$, where $`\tau _{\mathrm{dwell}}`$ is the time required for electrons in the normal region to feel the presence of the superconductor . The Thouless energy is determined by $`E_\mathrm{T}\mathrm{min}\{D/L^2,\mathrm{\Gamma }N\delta \}`$ where $`\mathrm{\Gamma }`$ is the transparency of the contact to the superconductor ($`\mathrm{\Gamma }=1`$ in the analysis of the zero-dimensional system above).

In the diffusive limit $`D/L^2\mathrm{\Gamma }N\delta `$, at the mean-field level, the position of the quasi-particle energy gap is found by solving the Usadel equation with the appropriate boundary conditions . As a result one obtains a square root singularity in the DoS. In this case the mean-field solution is itself inhomogeneous. The sub-gap correction is found by identifying a second inhomogeneous instanton configuration that breaks supersymmetry at the level of the action . Both solutions merge at the mean-field gap. Now, following the arguments above, it is simple to see how the phenomenology of Ref.  fits into the same general scheme: in this case the relevant coordinate of $`Q`$ interpolates between the inhomogeneous mean-field solution and the instanton. The result is a sub-gap DoS which assumes the familiar form of Eq. (35), with appropriately defined geometry dependent parameters $`E_{\mathrm{gap}}`$ and $`\mathrm{\Delta }_\mathrm{g}`$. Naturally the introduction of $`d_{}`$ transverse dimensions gives the expected energy dependence of $`(E_{\mathrm{gap}}ϵ)^{(6d_{})/4}`$ in the exponent.

In the opposite limit $`D/L^2\mathrm{\Gamma }N\delta `$ (not considered in Ref. ), gradients of $`Q`$ are heavily penalized and the coupling to the lead must be retained in its ‘logarithmic form’, with $`Q`$ being taken as constant in the dot. (Indeed, the logarithm is crucial to reproduce even the mean-field expression for the DoS (31) with the correct coefficients.) This is the true zero-dimensional limit treated above. As we have seen, with this action, one recovers the known universal expression for the spectrum of gap fluctuations below the mean-field edge.

## V Conclusions

To conclude, we have developed a quasi-classical field theory describing a superconductor in the dirty limit with weak magnetic impurity scattering. The Abrikosov-Gor’kov mean-field treatment of this system, showing a diminished but hard gap in the DoS, can be straightforwardly recovered as a homogeneous saddle-point of the effective action. Where zero DoS is predicted by the mean-field theory, there exist spatially inhomogeneous saddle-point configurations that break supersymmetry at the level of the action. A careful analysis of fluctuations around these instanton configurations demonstrate how supersymmetry is restored by a manifold of equivalent configurations parametrized by a Grassmann coordinate, but more importantly how the configurations give rise to a finite, though exponentially small, DoS. In contrast to band tail states in semi-conductors, the quasi-particle nature of the sub-gap states leads to universality of their properties.

Finally, let us remark on the connection of the results presented above to related problems in the literature. The resulting expression for the DoS (30) was found to be non-perturbative in the $`\sigma `$-model coupling $`1/g`$, which measures the strength of non-magnetic disorder. We note that other non-perturbative results in disordered systems have been obtained by related instanton calculations. As well as the investigation of tail states in semi-conductors , a supersymmetric field theory was developed by Affleck  (see also Refs. ) to investigate tail states in the lowest Landau level. There it was shown that tail states correspond to instanton configurations of the $`\mathrm{\Psi }`$-field action (c.f. Ref. ). It is also interesting to compare the present scheme with the study of ‘anomalously localized states’  (see also, Ref. ). There one finds that long-time current relaxation in a disordered wire is also associated with instanton configurations of the $`\sigma `$-model action. Finally, a Lifshitz argument has been applied on the level of the Usadel equation in the study of gap fluctuations due to inhomogeneities of the BCS interaction .

Although we have focussed largely on the question of tail states below the mean-field quasi-particle gap of a superconductor with magnetic impurities, we expect the instanton approach developed here to be more widely applicable. Indeed, in this work we have shown the intimate connection between the study of sub-gap states in the magnetic impurity problem and gap fluctuations in hybrid superconductor/normal structures. Furthermore, the same instanton approach describes gap fluctuations in superconductors with a quenched inhomogeneous distribution of the BCS coupling constant , as well as quasi-two dimensional superconducting films subject to strong in-plane magnetic fields . In both cases the latter are described by a mean-field theory which assumes the Abrikosov-Gor’kov form.

More speculatively, it seems likely that the same general scheme can be employed to study the influence of optimal fluctuations on the nature of bulk transitions. For example, in the present system, the transition to bulk superconductivity with reducing magnetic impurity concentration will be preempted by the nucleation of superconducting islands or droplets within the metallic/insulating phase (c.f. Ref. ). Similarly, the Stoner transition to a bulk itinerant ferromagnet in a disordered system will be mediated by the formation of a droplet phase in which islands become ferromagnetic . In both cases, we expect these ‘droplet phases’ to be associated with inhomogeneous instanton configurations of the corresponding low-energy action.

Acknowledgements: We are grateful to Alexander Altland, Alexander Balatsky, Alex Kamenev, Julia Meyer, and Martin Zirnbauer for valuable discussions. We are also deeply indebted to Dima Khmel’nitskii for bringing this general subject to our attention, and to John Chalker for crucial discussions, particularly those made at an early stage of this work. One of us (AL) would like to acknowledge the financial support of Trinity College.