# Metallic nanograins: spatially nonuniform pairing induced by quantum
confinement

M. D. Croitoru1,3 A. A. Shanenko2 C. C. Kaun3 F. M. Peeters2 1Institut für
Theoretische Physik III, Universität Bayreuth, 95440 Bayreuth, Germany
2Departement Fysica, Universiteit Antwerpen, Groenenborgerlaan 171, B-2020
Antwerpen, Belgium  3Research Center for Applied Sciences, Academia Sinica,
11529 Taipei, Taiwan

###### Abstract

It is well-known that the formation of discrete electron levels strongly
influences the pairing in metallic nanograins. Here we focus on another effect
of quantum confinement in superconducting grains that was not studied
previously, i.e., spatially nonuniform pairing. This effect is very
significant when single-electron levels form bunches and/or a kind of shell
structure: in highly symmetric grains the order parameter can exhibit
variations with position by an order of magnitude. Nonuniform pairing is
closely related to a quantum-confinement induced modification of the pairing-
interaction matrix elements and size-dependent pinning of the chemical
potential to groups of degenerate or nearly degenerate levels. For
illustration we consider spherical metallic nanograins. We show that the
relevant matrix elements are as a rule enhanced in the presence of quantum
confinement, which favors spatial variations of the order parameter,
compensating the corresponding energy cost. The size-dependent pinning of the
chemical potential further increases the spatial variation of the pair
condensate. The role of nonuniform pairing is smaller in less symmetric
confining geometries and/or in the presence of disorder. However, it always
remains of importance when the energy spacing between discrete electron levels
$\delta$ is approaching the scale of the bulk gap $\Delta_{B}$, i.e.,
$\delta>0.1$-$0.2\,\Delta_{B}$.

metallic nanograins, nonuniform pairing, superconducting correlations, matrix
elements

###### pacs:

74.20.Fg, 74.78.Na

##  I Introduction

Quantum confinement plays a fundamental role in superconductors with nanoscale
dimensions. Interplay of quantum confinement and pairing correlations results
in important qualitative changes in the superconductor characteristics. parm ;
blat ; mush ; smit ; brau ; perali ; sier ; glad ; yuzb ; kres01 ; kres02 ;
sh01 ; sh02 ; sh03 ; sh04 ; sh05 ; gar  Because of technological reasons
quasi-0D superconducting structures (i.e., ensembles of small grains) were the
first where this interplay was investigated experimentally. Initial attempts
by Giaever and Zeller at the end of 60s used tunneling studies on large
ensembles of superconducting particles. gaiv01  Since that time most of the
studies on superconducting correlations in grains were performed with grain
powders li01 ; li02  or on films made of crystalline granules separated by
amorphous inter-granular space. bose01 ; bose02  In the pioneering work of
Ralph et al. ralph01 ; ralph02  the discrete electron spectrum was measured
for a single grain. Their technique (single-electron tunneling spectroscopy)
enabled them for the first time to probe superconducting correlations in an
individual ${\rm Al}$ grain. Very recently, STM was used to detect the
superconducting gap of an isolated ultra-small lead grain deposited onto a
silicon substrate (see e.g. Refs. bose03, ; bose04, ). These advances opened
new prospects to examine superconductivity in individual metallic nanograins
with unprecedented detail, e.g., to investigate the influence of the
confinement on the superconducting correlations.

The main feature of a superconducting nanograin that makes them different from
a bulk superconductor is the formation of discrete electron levels with
average energy spacing $\delta\approx 2\pi^{2}\hbar^{2}/(mk_{F}V)$, with
$k_{F}$ the bulk Fermi wave number and $V$ the system volume. It can be of the
same order as the bulk gap $\Delta_{B}$, or even larger in the case of ultra-
small nanograins. Therefore, size-quantization of the electron spectrum can
have a substantial impact on the basic superconducting characteristics of such
quasi-0D superconducting systems.

The understanding of the fundamental properties of superconducting
correlations in low-dimensional structures, in particular in isolated metallic
grains, has experienced a remarkable development in the last two decades.
Theoretical aspects, which have attracted the most attention are the
following. The problem of the breakdown of BCS superconductivity in ultra-
small metallic grains was addressed in several papers. delft ; smit ; mat  The
effect of the shell structure in the single-electron spectrum on the
superconducting correlations was pointed out for nanograins glad ; gar  and
ultrasmall metallic clusters. kres01 ; kres02  The ground state properties of
the BCS pairing Hamiltonian of ultra-small grains were considered beyond the
mean-field approximation using the Richardson exact solution. sier ; glad ;
yuzb

A spatially uniform pairing was assumed in these and other works and, as a
consequence, the matrix elements of the pairing interaction were taken
independent of the relevant single-electron quantum numbers, i.e., they were
set to $-g/V$, with $g>0$ the coupling constant and $V$ the volume. note1
This is, say, a bulk-like approximation recovered when the single-electron
wave functions are taken as plane waves. However, the translational invariance
is broken in nanograins, which leads to a position-dependent order parameter.
As a result, the pairing gap becomes strongly dependent on the relevant
quantum numbers, which is directly related to a confinement-induced
modification (as compared to $-g/V$) of the matrix elements controlling the
scattering of the time reversed states. Another important issue is that
single-electron levels can form bunches and even a kind of shell structure in
symmetric confining geometries. In this case the chemical potential $\mu$ can
be pinned to a group of nearly degenerate or degenerate levels. This is of
importance because the density of states in the vicinity of $\mu$ strongly
influences the superconducting correlations. In other words, such a pinning
plays the role of a filter that selects the contribution of a particular
single-electron shell (or of a group of close levels) to the superconducting
order parameter. Such a contribution is, as a rule, spatially nonuniform.

The aim of the present paper is to investigate effects related to a spatially
nonuniform pairing in metallic nanograins, which was not studied in previous
publications. For illustrative purposes we consider metallic spherical
nanograins, where the spatial dependence of the superconducting condensate is
pronounced (the order parameter can vary with position by an order of
magnitude). In less symmetric confining geometries and/or in the presence of
disorder spatial variations of the order parameter are reduced. However, our
study implies that nonuniform pairing remains of importance when the
interlevel spacing $\delta$ is approaching the scale of the order of the bulk
gap $\Delta_{B}$. Any remaining grouping of single-electron levels, that is
always present in real samples, even strengthens the effect of interest. We
work in the mean-field approximation and, thus, stay in the regime
$\delta\lesssim\Delta_{B}$. Below we consider ${\rm Sn}$ and ${\rm Al}$ with
$\Delta_{B}=0.616$ and $0.25\,{\rm meV}$, respectively (for the parameters
used below). Using the above values of $\Delta_{B}$, we find that the mean-
field approach is valid for $D>6$-$8\,{\rm nm}$, with $D$ the sphere diameter.

Our paper is organized as follows. In Sec. II, we outline the formalism how to
obtain a self-consistent solution to the problem. In Sec. III, we present our
numerical results. In particular, in Sec. III.1 we investigate the effects of
quantum confinement on pairing correlations through the modifications of the
matrix elements of the pairing interaction and the size-dependent pinning of
$\mu$ to single-electron shells. Sec. III.2 is focused on a spatial
distribution of the pair condensate and its relation to modifications of the
matrix elements and the size-dependent pinning of $\mu$. In Sec. III.3 we
discuss the interplay of Andreev reflection with quantum confinement,
resulting in the formation of Andreev-type states and significant dependence
of the pairing gaps on the relevant quantum numbers. A short summary and
discussion are given in Sec. IV.

##  II Formalism

The reduction of the system to the nanometer scale leads to the formation of a
discrete electron spectrum. Moreover, in the presence of quantum confinement,
the translational invariance of the system is broken, and the superconducting
order parameter is position dependent, i.e., $\Delta=\Delta({\bf r})$. For the
mean-field treatment of such a situation, it is appropriate to use the
Bogoliubov-de Gennes (BdG) equations, bogl ; degen  which can be written as

$\displaystyle
E_{i}|u_{i}\rangle=\widehat{H}_{e}|u_{i}\rangle+{\widehat{\Delta}}|v_{i}\rangle,$
(1a) $\displaystyle
E_{i}|v_{i}\rangle={\widehat{\Delta}}^{\ast}|u_{i}\rangle-\widehat{H}_{e}^{\ast}|v_{i}\rangle,$
(1b)

where $E_{i}$ stands for the Bogoliubov-quasiparticle (bogolon) energy,
$\widehat{\Delta}=\Delta(\widehat{\bf r})$ (with $\widehat{\bf r}$ the
position operator) and the single-electron Hamiltonian is referred to the
chemical potential $\mu$, i.e.,

$\widehat{H}_{e}({\bf r})=\frac{\widehat{\bf p}^{2}}{2m_{e}}+V(\widehat{\bf
r})-\mu.$ (2)

We remark that any magnetic effects are beyond the scope of the present paper.
For simplicity, the confining interaction $V({\bf r})$ is taken as zero inside
the specimen and infinite outside: $V({\bf r})=V_{B}\;\vartheta(R-\rho)$ with
the barrier potential $V_{B}\to\infty$ ($R=D/2$ and $\rho$ is the radial
coordinate for the spherical confining geometry).

As a mean-field approach, the BdG equations should be solved in a self-
consistent manner

$\Delta({\bf r})=g\sum_{i}\langle{\bf r}|u_{i}\rangle\langle v_{i}|{\bf
r}\rangle\tanh(\frac{\beta E_{i}}{2}),$ (3)

where $g>0$ is the coupling constant for the effective electron-electron
interaction approximated by the delta-function potential, i.e., $\langle{\bf
r},{\bf r}^{\prime}|\Phi|{\bf r},{\bf r}^{\prime}\rangle=-g\delta({\bf r}-{\bf
r}^{\prime})$. The sum in Eq. (3) runs over the states with the single-
electron energy

$\xi_{i}=\bigl{[}\langle u_{i}|\widehat{H}_{e}|u_{i}\rangle+\langle
v_{i}|\widehat{H}_{e}|v_{i}\rangle\bigr{]}\in[-\hbar\omega_{D},\hbar\omega_{D}],$
(4)

with $\omega_{D}$ the Debye frequency. As is known, the solution of the BdG
equations has two branches: $(i,+)$ and $(i,-)$ (see Ref. swid, ) for which we
have $E_{i,+}>0$ and $E_{i,-}<0$. The sum in Eq. (3) should be taken over the
physical states [the $(i,+)$ branch], i.e., $E_{i}=E_{i,+}$.

For a given mean electron density $n_{e}$ the chemical potential $\mu$ is
determined from

$n_{e}=\frac{2}{V}\sum_{i}\bigl{[}f_{i}\langle
u_{i}|u_{i}\rangle+(1-f_{i})\langle v_{i}|v_{i}\rangle\bigr{]},$ (5)

with $V=\frac{4}{3}\pi R^{3}$ the volume of the spherical grain. For
conventional superconductors the energy gap is typically much smaller than the
chemical potential. As a result, $\mu$ stays nearly the same when passing from
the normal state to the superconducting one. degen  Therefore, one can solve
Eq. (5) in the absence of superconducting order ($\Delta({\bf r})=0$).

In a spherical nanograin, because of symmetry reasons, the order parameter
depends only on the radial coordinate, i.e., $\Delta=\Delta(\rho)$. Therefore
the pseudospinor in the particle-hole space can be characterized by the
quantum numbers of the angular momentum, i.e., ($l,m$). The angular part of
the pseudospinor $\Psi_{i}$ is given by the spherical harmonics
$Y_{lm}(\theta,\varphi)$ in polar coordinates $(\rho,\theta,\varphi)$, i.e.,

$\langle{\bf
r}|\Psi_{i}\rangle=Y_{lm}(\theta,\varphi)\left(\begin{array}[]{c}u_{jl}(\rho)\\\
v_{jl}(\rho)\end{array}\right),$ (6)

where $i=\\{j,l,m\\}$, with $j$ the radial quantum number associated with the
quantum-confinement boundary conditions

$u_{jl}(\rho)|_{\rho=R}=v_{jl}(\rho)|_{\rho=R}=0.$ (7)

To solve the BdG equations (1a) and (1b) numerically, $u_{jl}(\rho)$ and
$v_{jl}(\rho)$ are expanded in the eigenfunctions of the single-electron
Hamiltonian $\widehat{H}_{e}$ [see Eq. (2)]. In addition, iterations should be
invoked, to account for the self-consistency relation given by Eq. (3). This
program is significantly simplified by keeping only the pairing of the time-
reversed states, and  which is a standard approximation for the problem of
superconducting correlations in nanograins. In the framework of the BdG
equations this can be done through the so-called Anderson approximate solution
for which the particle- and hole-like wave functions are assumed to be
proportional to the single-electron wave function. It means that

$u_{jl}(\rho)={\cal U}_{jl}\,\chi_{jl}(\rho),\;u_{jl}(\rho)={\cal
V}_{jl}\,\chi_{jl}(\rho),$ (8)

with the radial part of the single-electron wave function given by

$\chi_{jl}(\rho)=\frac{\sqrt{2}}{R^{3/2}j_{l+1}(\alpha_{jl})}j_{l}(\alpha_{jl}\frac{\rho}{R}),$
(9)

with $j_{l}(x)$ the $l$-order spherical Bessel function of the first kind and
$\alpha_{jl}$ its $j$-node. The coefficients ${\cal U}_{jl}$ and ${\cal
V}_{jl}$ (taken as real) obey the standard constraint (see, e.g., Refs. ket, )

${\cal U}^{2}_{jl}+{\cal V}^{2}_{jl}=1.$ (10)

Then, inserting Eq. (8) into Eqs. (1a) and (1b) we find the following set of
coupled equations (here $E_{jlm}=E_{jl}$ and $\xi_{jlm}=\xi_{jl}$):

$\displaystyle[E_{jl}-\xi_{jl}]~{}{\cal U}_{jl}=\Delta_{jl}\,{\cal V}_{jl},$
(11a) $\displaystyle[E_{jl}+\xi_{jl}]~{}{\cal V}_{jl}=\Delta_{jl}\,{\cal
U}_{jl},$ (11b)

with

$\Delta_{jl}=\int\limits_{0}^{R}{\rm
d}\rho\,\rho^{2}\,\chi^{2}_{jl}(\rho)\Delta(\rho)$ (12)

and

$\xi_{jl}=\frac{\hbar^{2}}{2m_{e}}\frac{\alpha_{jl}^{2}}{R^{2}}-\mu.$ (13)

A nontrivial physical solution of Eqs. (11a) and (11b) exists only when

$E_{jl}=\sqrt{\xi_{jl}^{2}+\Delta^{2}_{jl}}.$ (14)

The Anderson prescription about the pairing of the time-reversed states allows
one to rephrase the self-consistency relation [see Eq. (3)] as follows:

$\Delta_{j^{\prime}l^{\prime}}=-\sum\limits_{jl}(2l+1)\;\frac{M_{j^{\prime}l^{\prime},jl}\;\Delta_{jl}}{2\sqrt{\xi^{2}_{jl}+\Delta^{2}_{jl}}}\tanh(\frac{\beta
E_{jl}}{2}),$ (15)

where

$\displaystyle
M_{j^{\prime}l^{\prime},jl}=-\frac{g}{4\pi}\int\limits_{0}^{R}\\!\\!{\rm
d}\rho\,\rho^{2}~{}\chi_{j^{\prime}l^{\prime}}^{2}(\rho)~{}\chi_{jl}^{2}(\rho).$

To derive Eq. (15), one should keep in mind the property of the spherical
harmonics $\sum_{m=-l}^{l}\;|Y_{lm}(\theta,\varphi)|^{2}=\frac{2l+1}{4\pi}$.
We remark that $M_{j^{\prime}l^{\prime},jl}$ is nothing else but the pairing-
interaction matrix element $\langle
i^{\prime},\bar{i^{\prime}}|\Phi|i,\bar{i}\rangle$ (with
$\bar{i}=\\{j,l,-m\\}$) averaged over the states with $m=-l,\ldots l$ and
$m^{\prime}=-l^{\prime},\ldots l^{\prime}$, i.e.,

$M_{j^{\prime}l^{\prime},jl}=\frac{1}{(2l^{\prime}+1)(2l+1)}\sum\limits_{m^{\prime}=-l^{\prime}}^{l^{\prime}}\sum\limits_{m=-l}^{l}\langle
i^{\prime},\bar{i^{\prime}}|\Phi|i,\bar{i}\rangle.$

As seen from Eq. (12), a spatially uniform order parameter means that the
pairing gaps $\Delta_{jl}$ do not depend on the quantum numbers $j$ and $l$.
This is compatible with Eq. (15) only when $M_{j^{\prime}l^{\prime},jl}$ does
not depend on $j^{\prime}$ and $l^{\prime}$. According to the definition given
by Eq. (LABEL:M_mat), we have
$M_{j^{\prime}l^{\prime},jl}=M_{jl,j^{\prime}l^{\prime}}$ and, so, if
$M_{j^{\prime}l^{\prime},jl}$ does not depend on $j^{\prime},l^{\prime}$, it
does not depend on $j,l$ either. So, we arrive at the standard simplified
approach of investigating the pairing correlations in metallic grains (see the
discussion in the Introduction). Below we show that the spatial dependence of
the order parameter can not be ignored in superconducting nanograins, which
implies significant variations of the matrix elements and pairing gaps with
the relevant quantum numbers. After a numerical solution of Eq. (15), the
position-dependent order parameter can be calculated from

$\Delta(\rho)=\sum\limits_{jl}\Delta^{(jl)}(\rho),$ (16)

with the shell-dependent contribution $\Delta^{(jl)}(\rho)$ given by

$\Delta^{(jl)}(\rho)=\frac{g}{8\pi}\,(2l+1)\,\frac{\chi_{jl}^{2}(\rho)\,\Delta_{jl}}{\sqrt{\xi^{2}_{jl}+\Delta^{2}_{jl}}}\tanh(\frac{\beta
E_{jl}}{2}).$ (17)

##  III Discussion of results

###  III.1 Enhanced intrashell matrix elements and quantum-size pinning of the
chemical potential

Figure 1: Critical temperature versus the grain diameter as calculated for:
(a) $M_{j^{\prime}l^{\prime},jl}\neq-g/V$ and $\mu\neq\mu_{B}$; (b)
$M_{j^{\prime}l^{\prime},jl}=-g/V$ and $\mu\neq\mu_{B}$; and (c)
$M_{j^{\prime}l^{\prime},jl}\neq-g/V$ and $\mu=\mu_{B}$. The dashed curves in
(a) show approximate lower and upper boundaries for the quantum-size
oscillations of $T_{c}$, both curves represent the same dependence
$T_{c}/T_{c,B}=1+a(D/D_{0})^{3/2}$, with $D_{0}=50\,{\rm nm}$ and $a=1$ (the
lower boundary) and $a=3.5$ (the upper one). The same curves are also given in
(b) and (c), for comparison.

Numerical calculations were performed with the set of parameters typical for
tin degen ; fett : $\hbar\omega_{D}/k_{B}=195~{}{\rm K}$, $gN(0)=0.25$, with
$N(0)$ the bulk density of states at the Fermi level (we use the bulk electron
density $n_{e}=148\,{\rm nm}^{-3}$, see, e.g., Ref. ash, ).

Figure 1(a) shows the critical temperature (in units of the bulk critical
temperature $T_{c,B}$) versus the nanograin diameter $D$ as calculated from
Eq. (15) when the matrix elements of the electron-electron interaction and the
size-dependent variation of the chemical potential have been fully taken into
account. Results in Fig. 1 are presented for a step $\Delta R=0.01~{}{\rm
nm}$. For each radius the critical temperature was defined as the temperature
above which the spatially-averaged order parameter
$\langle\Delta(\rho)\rangle$ becomes smaller than $0.01$ of its value at
$T=0$. Our numerical results exhibit two features typical of the size-
dependent pairing characteristics in high-quality superconducting nanograins
and nuclei. First, we observe an overall increase of $T_{c}$ with decreasing
$D$ (it is very pronounced due to the highly-symmetric confining geometry).
Second, $T_{c}$ oscillates wildly with $D$. This oscillatory behavior can be
understood in the following way. The pair correlations are nonzero only for
the states within a finite range (the Debye window) around the chemical
potential $\mu$. Moreover, the main contribution to the sum in Eq. (15) comes
from the states in the very vicinity of the Fermi level, because in this case
the expression $\Delta_{jl}/\sqrt{\xi^{2}_{jl}+\Delta^{2}_{jl}}\simeq 1$
$(\xi_{jl}\simeq 0)$. When varying the nanograin size, the number of states in
the Debye window changes. The smaller the diameter, the smaller the number of
relevant states contributing to the pairing characteristics and, as a result,
the more significant is such a change. This change is not monotonous but
rather oscillating due to a permanent competition between incoming and
outcoming states. As a consequence, all basic pairing characteristics, e.g.,
$T_{c}$ and pairing gaps $\Delta_{jl}$, exhibit quantum-size oscillations. It
is not only typical of nanograins with superconducting correlations (see,
e.g., the recent paper bose04 ) but it is also present in superconducting
nanowires sh01 ; sh02 ; sh03 ; sh04 ; sh05  and nanofilms. guo ; eom  Such
oscillations are pronounced for small diameters/thicknesses but decay with
increasing the characteristic size so that $T_{c}$ approaches the bulk
critical temperature $T_{c,B}$ (for our parameter $T_{c,B}=4.01\,{\rm K}$). It
is interesting to note that the overall increase of $T_{c}$ with decreasing
$D$ in Fig. 1(a) is similar to a size-dependent enhancement of the pairing gap
in nuclei, where it is proportional to $1/\sqrt{A}$ (see, e.g., Ref. sat, ),
with $A$ the number of nucleons. In particular, the two dashed curves in Fig.
1(a) show approximate upper and lower boundaries for $T_{c}$, highlighting the
magnitude of the quantum-size oscillations: both curves represent the same
dependence, i.e., $T_{c}/T_{c,B}=1+a\,(D_{0}/D)^{3/2}$, with $D_{0}=50\,{\rm
nm}$ and $a=1$ and $3.5$ for the lower and upper boundaries, respectively
[$(D_{0}/D)^{3/2}\propto N^{-1/2}_{e}$, with $N_{e}=n_{e}V$ the number of
electrons]. We remark that real samples exhibit inevitable shape and size
fluctuations that affect the high-degeneracy of single-electron levels. Hence,
measurements on an ensemble of nanograins will significantly smooth the
quantum-size oscillations in the critical temperature and reduce its overall
enhancement with decreasing nanograin size (see, also, Sec. IV). For instance,
in experimentally fabricated tin nanograins of a semi-spherical shape the
observed enhancement of the excitation gap over its bulk value is about bose04
$60\%$ for the particle heights $\approx 10$-$20\,{\rm nm}$. This is
significantly smaller than the enhancement of $T_{c}$ shown in Fig. 1(a).
However, detailed investigations of the enhancement of $T_{c}$ in
superconducting nanograins is beyond the scope of our present paper. Here we
are interested in a spatially nonuniform distribution of the pair condensate
which is of importance even in the presence of shape and size fluctuations and
disorder (see the discussion in Sec. IV).

In order to outline the role of the matrix elements
$M_{j^{\prime}l^{\prime},jl}$ [see Eq. (LABEL:M_mat)] of the electron-electron
interaction we also show what happens when the true matrix elements are simply
replaced by those of the bulk-like form: $M_{j^{\prime}l^{\prime},jl}=-g/V$,
which is what is usually done when investigating the superconducting
correlations in nanograins. The results are displayed in Fig. 1(b) and, as
seen, the difference with respect to Fig. 1(a) is significant. To simplify the
comparison, we show also in Fig. 1(b) two solid curves that represent the
radius-dependent upper and lower values of $T_{c}$ from Fig. 1(a).

Table 1: Matrix elements $M_{j^{\prime}l^{\prime},jl}=M_{jl,j^{\prime}l^{\prime}}$ in units of $-g/V$ calculated at $D=7.1\,{\rm nm}$ for quantum numbers such that $\xi_{j^{\prime}l^{\prime}},\xi_{jl}<\hbar\omega_{D}$: $M_{j^{\prime}l^{\prime},jl}$ | $j^{\prime}$ | $l^{\prime}$ | $j$ | $l$
---|---|---|---|---
10.62 | 31 | 11 | 31 | 11
1.9 | 31 | 11 | 23 | 29
1.33 | 31 | 11 | 19 | 39
0.64 | 31 | 11 | 8 | 71
0.41 | 31 | 11 | 1 | 101
4.71 | 23 | 29 | 23 | 29
1.7 | 23 | 29 | 19 | 39
0.7 | 23 | 29 | 8 | 71
0.43 | 23 | 29 | 1 | 101
3.72 | 19 | 39 | 19 | 39
0.77 | 19 | 39 | 8 | 71
0.46 | 19 | 39 | 1 | 101
2.69 | 8 | 71 | 8 | 71
0.69 | 8 | 71 | 1 | 101
3.61 | 1 | 101 | 1 | 101

To clarify the physical reason why using the true matrix elements leads to
significant deviations from the results found for
$M_{j^{\prime}l^{\prime},jl}=-g/V$, we show in Table I the numerical values of
$M_{j^{\prime}l^{\prime},jl}$ (calculated in units of $-g/V$) for
$D=14.2\,{\rm nm}$ (only the states within the Debye window are given). As
seen, the diagonal (intrashell) matrix elements $M_{jl,jl}$ are strongly
enhanced as compared to $-g/V$. However, the matrix elements controlling the
scattering of the time reversed states between different shells (intershell)
are often decreased in absolute value with respect to $-g/V$. So, the question
arises why the superconducting correlations are enhanced for the true matrix
elements? The point is that the intershell interactions are of less importance
due to a size-dependent pinning of the chemical potential to the groups of
degenerate or nearly degenerate levels (shells can be often close to each
other in energy), see the next paragraph. When $\mu$ is pinned to a particular
shell, then the single-electron energy measured from $\mu$ is zero for the
states from this shell. These states make a major contribution to
superconducting correlations unless diameters are not large enough
$D<20$-$30\,{\rm nm}$, in other words, the number of contributing shells is
less than $10$-$15$. In this case the superconducting correlations are nearly
determined by the pairing gap $\Delta_{jl}$ associated with the shell pinned
to $\mu$. From Eq. (15) it is seen that $\Delta_{jl}$ for the states with
$\xi_{jl}=0$ is mainly governed by the intrashell matrix element $M_{jl,jl}$.
For instance, when ignoring the contribution of all other states one simply
obtains (at $T=0$)

$\Delta_{jl}\approx-(l+\frac{1}{2})\,M_{jl,jl}.$

When the diameter increases beyond $20$-$30\,{\rm nm}$, then the intershell
matrix elements approach $-g/V$ while the intrashell matrix elements are still
significantly different from the bulk-like behavior. However, the role of the
states with $\xi_{jl}=0$ is becoming less and less important for larger
diameters due to the presence of larger and larger number of shells making a
contribution to the pairing correlations. As a consequence, the difference
between the data in Figs. 1(a) and (b) decreases when approaching
$D=35$-$40\,{\rm nm}$, together with the amplitude of the quantum-size
oscillations of $T_{c}$.

Figure 2: (a) Size variations of the chemical potential, accompanied by an
overall shift of $\mu$ to upper values with decreasing $D$. (b) Details of the
quantum-size pinning of $\mu$ (filled squares) to the single-electron levels
(solid curves), small diameters are shown for simplicity.

Figure 3: Spatial distribution of the pair condensate in spherical nanograins:
$\Delta(\rho)$ (calculated at $T=0$) versus $\rho$ for diameters $D=12\,{\rm
nm}$ (a), $13.52\,{\rm nm}$ (b), $14.2\,{\rm nm}$ (c), $16\,{\rm nm}$ (d),
$16.4\,{\rm nm}$ (e) and $17.54\,{\rm nm}$ (f).

In the fully self-consistent scheme the chemical potential is determined in
such a way that the mean electron density $n_{e}$ is constant [see Eq. (5)].
However, size-dependent variations of $\mu$ are of importance not only because
they simply prevent the mean electron density from deviations. In fact, such
deviations are almost insignificant: our calculations for $\mu=\mu_{B}$ show
that $n_{e}$ decreases by a few percent when $D$ reduces to $10$-$20\,{\rm
nm}$. A more interesting thing is that the size-dependent variations of $\mu$
have a pronounced effect on the superconducting correlations. In particular,
this can be seen from Fig. 2(c), where $T_{c}$ is calculated for the true
matrix elements and $\mu=\mu_{B}$. What is the reason for this suppression of
$T_{c}$? In the presence of the formation of strongly degenerate electron
levels or bunches of electron levels with almost negligible spacing between
them, the chemical potential lies mostly at the highest partly-filled
degenerate level (see, e.g., Refs. kres01, ; kres02, ). Pairing correlations
are significant only within the Debye window around the chemical potential
$\mu$ and are strongest delft  exactly at $\mu$. Hence, when $\mu$ is pinned
to a shell level, this favors the pairing correlations at this level and, in
turn, through the self-consistency relation, favors the pairing correlations
at neighboring shells. In other words, if the level to which the chemical
potential is pinned is highly degenerate than the phase space for the
strongest pair scattering is enlarged and, consequently, the system gains in
interaction energy and, as a result, superconducting correlations are strongly
enhanced. It is different when $\mu$ is not pinned to a shell, which is mostly
the case for a constant chemical potential, e.g., for $\mu=\mu_{B}$. Here the
relevant shells entering the Debye window are as a rule specified by
$\xi_{jl}\not=0$ and, so, their contributions are diminished.

The above discussion is further illustrated by our numerical results for $\mu$
in Fig. 2. As seen from panel (a), when keeping the electron density of the
system constant, $\mu$ slightly shifts systematically up with decreasing $D$
and exhibits size-dependent oscillations, as seen from Fig. 2(a). These
oscillations are a signature of the size-dependent pinning of $\mu$ to groups
of degenerate or nearly degenerate single-electron levels. This is clearly
seen from Fig. 2(b), where variations of $\mu$ (filled squares) are plotted
versus $D$ together with the single-electron energies measure from the band
bottom, i.e., $\frac{\hbar^{2}}{2m_{e}}\frac{\alpha^{2}_{jl}}{R^{2}}$ (solid
curves). For the sake of simple illustration, panel (b) shows the data for
extremely small diameters, where the energy spacing between the shell levels
is pronounced and, as a result, the size-dependent oscillations of $\mu$ are
not so wild as it happens for higher diameters. As follows from Fig. 2(b)
$\mu$ is pinned to a shell level in most cases, which, as mentioned above,
represents incomplete shells. Sometimes $\mu$ can be found between two
neighboring shell levels, which corresponds to the case of a fully occupied
lower shell.

###  III.2 Spatially nonuniform pair condensate

In the previous paragraph we considered the effect of quantum confinement on
pairing correlations through the matrix elements and quantum-size pinning of
$\mu$. As discussed at the end of Sec. II, a framework which incorporates both
issues appears to be only consistent when the position-dependent
superconducting order parameter is taken into consideration. Thus, our results
discussed in the previous section suggest that the spatial variations of
$\Delta(\rho)$ will be pronounced even in nanograins with diameters up to
$D=20$-$30\,{\rm nm}$. However, it is usually argued that spatial variations
of $\Delta(\rho)$ cost significant extra energy and, so, they are strongly
suppressed when $D\ll\xi$, with $\xi$ the bulk coherence length (see, for
instance, Ref. mat, ). In addition, $D$ should be larger than $\lambda_{F}$:
in practice, $k_{F}D\sim 10$ is assumed to be sufficient to ignore any spatial
dependence of the order parameter. kres01 ; kres02  For typical metallic
parameters $k_{F}D\sim 200$-$400$ for $D=10$-$20\,{\rm nm}$ and this is the
reason why the spatial dependence of the order parameter was ignored in most
papers on superconducting correlations in nanograins. To go in a more detail
on this point, we below discuss our numerical results on $\Delta(\rho)$.

In Fig. 3 the radial dependence of the superconducting order parameter is
shown as calculated from Eq. (16) for $D=12\,{\rm nm}$ (a), $13.52\,{\rm nm}$
(b), $14.2\,{\rm nm}$ (c), $16\,{\rm nm}$ (d), $16.4\,{\rm nm}$ (e) and
$17.54\,{\rm nm}$ (f). The shells making a contribution to the superconducting
correlations are also displayed in each panel, and the quantum numbers of the
shell level pinned to $\mu$ are underlined. As seen, we in general have a
nonuniform distribution of the pair condensate for diameters $D=10$-$20\,{\rm
nm}$, which is in agreement with our expectations. For example, let us
consider the results plotted in panel (c). Here $\mu$ is pinned to the shell
level $(l,j)=(101,1)$ and, so, single-electron states with $j=1$ and $l=101$
make a major contribution to $\Delta(\rho)$, which results in a significant
enhancement of the order parameter next to the edge, i.e., for
$\rho/R=0.9$-$1.0$. The profile of this enhancement is determined by the
radial wave function $\chi^{2}_{1,101}(\rho)$ with two pronounced local maxima
($\Delta/\Delta_{B}=14.3$ and $7.2$ at $\rho/R=0.9$ and $0.97$, respectively)
and one node (recall that $j$ is the number of the nodes of the radial wave
function). All the other shells displayed in Fig. 3(c) are specified by
$\xi_{jl}\not=0$ and, as a result, their contributions is much less
significant. The local maximum $\Delta(\rho)/\Delta_{B}=2.3$ at $\rho/R=0.1$
is due to states $(j,l)=(31,11)$. The shells with $(j,l)=(23,29)$ and
$(19,29)$ are responsible for local enhancements of the order parameter up to
$1.5$-$2.0\Delta_{B}$ at $\rho/R=0.27$ and $0.36$, respectively. At last, the
shell $(8,71)$ produces the local maximum at $\rho/R=0.64$. In general, the
larger the angular momentum, the larger the values of $\rho/R$ at which the
corresponding single-electron states have an effect on the profile of
$\Delta(\rho)$.

It is worth noting that typically, the order parameter is strongly suppressed
in the center ($\rho=0$) except of rare cases when states with zero angular
momentum contribute to the pairing correlations. One such example is given in
Fig. 3(d), where a narrow pick can be seen at $\rho=0$ due to the contribution
of the shell with $(j,l)=(41,0)$.

Figure 4: The order parameter $\Delta(\rho)$ for sufficiently large diameters
$D=33.6\,{\rm nm}$ (a), $D=35\,{\rm nm}$ (b) and $D=36\,{\rm nm}$ (c).

From Fig. 3 it follows that the radial distribution of the pair condensate
remains strongly nonuniform even for $D\approx 20\,{\rm nm}$. We would like to
note that when selecting concrete values of $D$ for Fig. 3, we did not even
take diameters for which $T_{c}$ is close to the upper dashed curve in Fig.
1(a). In the case of a strong enhancement of $T_{c}$ the radial distribution
of the pair condensate is as a rule strongly nonuniform. The points selected
for Fig. 3 are mainly in a vicinity of the lower dashed curve in panel (a) of
Fig. 1: for $D=13.52,\,14.2,\,16$ and $16.4\,{\rm nm}$ we have
$T_{c}/T_{c,B}=6.41,\,6.48,\,4.082$ and $3.78$, respectively. However, even in
this case the order parameter can vary with position by an order of magnitude.
Spatial variations of $\Delta(\rho)$ are significantly relaxed only when $D$
approaches $30$-$40\,{\rm nm}$, as seen from Fig. 4.

For our parameters $k_{F}=16.4\,{\rm nm}$ and, so, we obtain $k_{F}D\approx
300$ for $D\approx 20\,{\rm nm}$. Hence, the criterion $k_{F}D\gg 1$ is not
very useful in order to estimate the effect of spatial variations of the pair
condensate. Based on our numerical study, we would like to suggest another
criterion related to a more sensitive energy scale, which in the
superconducting state is governed by the bulk pairing gap $\Delta_{B}$. The
spatial distribution of the order parameter is always strongly inhomogeneous
when $\delta\sim\Delta_{B}$ (here it is even better to replace $\Delta_{B}$ by
the size-dependent pairing gap). The spatial variations decay with a decrease
in the ratio of the mean interlevel spacing to the bulk order parameter, i.e.,
$\delta/\Delta_{B}$, and our numerical results suggest that such variations
are significantly reduced only when $\delta/\Delta_{B}<0.05$-$0.1$ (recall
that effects of a magnetic field are beyond the scope of our paper). For ${\rm
Sn}$ spherical superconducting grains this regime is achieved when
$D>40$-$50\,{\rm nm}$ (note that $\delta\approx 2\pi^{2}\hbar^{2}/(mk_{F}V)$
underestimates the intershell spacing for spherical confining potential).
Despite that our results are for a highly symmetric confining geometry, we can
expect that the order parameter will be always spatially nonuniform for
$\delta/\Delta_{B}>0.1$, even when shape imperfections and disorder dissolve a
shell structure. The reason is that the number of contributing states (i.e.,
the states in the energy interval $\approx[\mu-\Delta_{B},\mu+\Delta_{B}$) is
not very large for $\delta/\Delta_{B}>0.1$. In this case the states pinned to
$\mu$ always make a major contribution to the order parameter and, so, the
profile of the squared absolute value of the corresponding wave function will
mainly determine the spatial distribution of the condensate. Thus, the domain
$\delta/\Delta_{B}=0.1$-$1.0$ is in general characterized by strong effects
due to the spatially nonuniform pairing.

We remark that our conclusions do not contradict the usual argument that
spatial variations of the order parameter cost extra energy. Let us compare a
bulk superconductor with a superconducting nanograin. In bulk the relevant
matrix elements controlling the scattering of the time-reversed states are
$-g/V$ and the order parameter is spatially uniform (in the absence of a
magnetic field). As opposed to bulk, the pair condensate significantly varies
with position in nanograins, which results, of course, in an increase of the
kinetic energy. However, the intrashell matrix elements are now enhanced in
absolute value as compared to $-g/V$ due to quantum confinement. This
compensates energy costs of spatial variations of the order parameter.

The discussion in the previous paragraph is also related to arguments that
invoke the conventional Ginzburg-Landau theory. According to this arguments
the order parameter is uniform in samples with size smaller than the bulk
coherence length. When applying this to nanograins, one can conclude that the
pair condensate should not vary with its position. However, this is not true.
It is well-known that one should be careful when applying the conventional
Ginzburg-Landau theory to superconductors with characteristic size smaller
than the zero-temperature (BCS) coherence length $\xi_{0}$. Strictly speaking,
Gor’kov’s derivation of the conventional Ginzburg-Landau formalism from the
BCS approach is not applicable on a scale smaller than $\xi_{0}$ (see,.e.g.,
Ref. fett, ). For ${\rm Sn}$ we have $\xi_{0}\approx 230\,{\rm nm}$ (see,
e.g., Ref. degen, ). Thus, in the case of interest $D\ll\xi_{0}$, and one can
hardly invoke the conventional Ginzburg-Landau formalism to check whether or
not $\Delta(\rho)$ varies with $\rho$.

###  III.3 Confinement-induced Andreev-type states

Here we would like to discuss one more issue related to a spatially nonuniform
pairing in nanograins. This is the formation of Andreev-type states induced by
quantum confinement sh03 ; sh06  (see also a similar paper torm  discussing
Andreev-type states in an ultracold trapped superfluid Fermi gas). Since the
60s (see Refs. saint, ; car, ; andr, ) it is known that quasiparticles can
“feel” a spatial variation of the superconducting order parameter as a kind of
potential barrier. This physical mechanism (referred to as Andreev mechanism
below) is the basis for Andreev quantization investigated previously for the
core of a single vortex for the mixed state of a type-II superconductor car
and for an isolated normal region of the intermediate state of a type-I
superconductor andr  (or for a similar case of SNS contacts saint ). Based on
our consideration of Sec. III.2, one can expect that Andreev-type states can
play a remarkable role in superconducting nanograins due to significant
spatial variations of the superconducting order parameter. This is very
similar to recently investigated Andreev-type states in superconducting
nanowires/nanofilms sh03 ; sh06 , where the pair condensate is position
dependent in the direction perpendicular to the nanowire/nanofilm due to the
quantization of the perpendicular electron motion. In Ref. sh06,  it was shown
that

$\Delta_{i}=\int\\!{\rm d}^{3}r\,\Delta({\bf r})\,\Bigl{[}|u_{i}({\bf
r})|^{2}+|v_{i}({\bf r})|^{2}\Bigr{]},$ (18)

which means that the pairing energy gap $\Delta_{i}$ is the averaged value of
the order parameter ”watched” by the quasiparticles with quantum numbers $i$.
Note that $|u_{i}({\bf r})|^{2}+|v_{i}({\bf r})|^{2}$ can be interpreted as
the spatial distribution of quasiparticles according to the well-known
constraint $\int{\rm d}^{3}r(|u_{i}({\bf r})|^{2}+|v_{i}({\bf r})|^{2})=1$
[see, e.g., Ref. ket,  and Eq. (10)]. When inserting Eqs. (6) into Eq. (18),
one can easily obtain Eq. (12) with $\Delta_{i}=\Delta_{jl}$. If
quasiparticles avoid the domains of enhanced pair condensate, the
corresponding integral in the right-hand-side of Eq. (18) becomes smaller and,
hence, such quasiparticles have smaller pairing gaps $\Delta_{jl}$. They can
be referred to as Andreev-type states.

Our numerical study of quantum-number dependent pairing gaps $\Delta_{jl}$ for
metallic nanograins reveals a significant role of Andreev mechanism. Let us
consider $D=13.52\,{\rm nm}$, the corresponding spatial distribution of the
pair condensate is given in Fig. 3(b). To show how different species of
quasiparticles are distributed in the radial direction in this case, the
radial-dependent shell contributions (at $T=0$) $\Delta^{(jl)}(\rho)$ [see
Eqs. (16) and (17)] are plotted in Fig. 5(a). We remark that such a
representation is more informative than simply a plot of
$|u_{jl}(\rho)|^{2}+|v_{jl}(\rho)|^{2}$. First, the radial dependence of
$\Delta^{(jl)}\propto\chi^{2}_{jl}(\rho)$ is the same as that of
$|u_{jl}(\rho)|^{2}+|v_{jl}(\rho)|^{2}\propto\chi^{2}_{jl}(\rho)$ [see Eq.
(8)]. Second, a plot of $\Delta^{(jl)}(\rho)$ gives also information how the
corresponding states contribute to $\Delta(\rho)$. From Fig. 3(c) we can see
that a significant enhancement of the order parameter occurs at
$\rho/R=0.45$-$0.7$. From Fig. 5(a) it is clear that this enhancement is due
to the states with $(j,l)=(14,48)$ and $(9,63)$. Other shells, i.e., $(27,16)$
and $(23,25)$ contribute less, and the corresponding quasiparticles,
representing Andreev-type states, are mainly located beyond the domain
$\rho=0.45-1.0$. As a result, they have smaller pairing gaps, i.e.,
$\Delta_{27,16}=2.65\,\Delta_{B}$ and $\Delta_{23,25}=2.81\,\Delta_{B}$, as
compared to $\Delta_{14,48}=4.098\,\Delta_{B}$ and
$\Delta_{9,63}=5.77\,\Delta_{B}$. As seen, the quasiparticles with
$(j,l)=(27,16)$ are most successful in avoiding the local enhancement of
$\Delta(\rho)$ at $\rho/R=0.45$-$0.7$ and, so, $\Delta_{27,16}$ is the
smallest pairing gap. Such a manifestation of Andreev mechanism is not a
particular feature of $D=13.52\,{\rm nm}$. In general, $\Delta_{jl}$ strongly
varies with $j$ and $l$ for diameters $<30-40\,{\rm nm}$, i.e., where spatial
variations of the order parameter are still pronounced. Quite often such
variations can be an order of magnitude, as, e.g., for $D=14.2\,{\rm nm}$ (see
$\Delta(\rho)$ given in Fig. 2(c)). At this diameter a great enhancement of
$\Delta(\rho)$ takes place at $\rho/R=0.9$. This is due to the contribution of
the shell with $(j,l)=(1,101)$ [see Fig. 5(b)]. Other shells make much less
important inputs and the corresponding quasiparticles are mainly distributed
beyond the domain $\rho/R=0.9$-$1.0$. So, as compared to
$\Delta_{1,101}=9.32\,\Delta_{B}$, they have significantly smaller pairing
gaps, i.e., $\Delta_{31,11}=1.6\Delta_{B}$, $\Delta_{23,29}=1.62\,\Delta_{B}$,
$\Delta_{19,39}=1.72\,\Delta_{B}$ and $\Delta_{8,71}=2.35\,\Delta_{B}$. Thus,
the interplay of Andreev mechanism and quantum confinement is responsible for
variations of $\Delta_{jl}$ with the relevant quantum numbers.

Figure 5: Shell-dependent contributions to the order parameter
$\Delta_{jl}(\rho)$ for relevant shells: (a), $D=13.52\,{\rm nm}$,
$(j,l)=(27,16),\,(23,25),\,(14,48)$ and $(9,63)$; (b) $D=14.2\,{\rm nm}$,
$(j,l)=(31,11),\,(23,29),\,(19,39),\,(8,71)$ and $(1,101)$.

One could expect that such a serious difference in pairing gaps of different
quasiparticle species can result in a pronounced drop of the ratio of
$\Delta_{E}$ (the minimal energy gap) to the critical temperature
$k_{B}T_{c}$, similar to the case for quantum superconducting nanowires. sh03
The main idea here is that $\Delta_{E}$ is governed by Andreev-type states
and, hence, is decreased. Unlike $\Delta_{E}$, $T_{c}$ is controlled by the
quasiparticles making a major contribution to $\Delta(\rho)$ and, so, $T_{c}$
is coupled to their higher pairing gaps. As a result, $\Delta_{E}/k_{B}T_{c}$
can be significantly smaller than in bulk. For instance, one can expect that
$\Delta_{E}=\Delta_{31,11}=1.6\,\Delta_{B}$ at $D=14.2\,{\rm nm}$ while
$T_{c}$ is governed by $\Delta_{1,101}=9.32\,\Delta_{B}$. However, this is not
correct for nanograins. The point is that $\Delta_{E}$ is a spectroscopical
gap which is probed by STM. It is defined as $\Delta_{E}=\min_{jl}E_{jl}$. For
nanowires the subband-dependent pairing gap is always the minimal
quasiparticle energy due to a quasi-free spectrum in the direction parallel to
the nanowire. For nanograins this is different. In particular, for
$D=14.2\,{\rm nm}$ we have the following single-electron energies (absorbing
$\mu$) of the relevant shells: $\xi_{31,11}=-18.6\,\Delta_{B}$,
$\xi_{23,29}=-26.03\,\Delta_{B}$, $\xi_{19,39}=20.9\,\Delta_{B}$,
$\xi_{8,71}=22.6\,\Delta_{B}$ and $\xi_{1,101}=0$. Hence, one can calculate
that $\Delta_{E}=E_{1,101}=\Delta_{1,101}$ in spite of the fact that
$\Delta_{1,101}$ is the largest pairing gap. Thus, although Andreev mechanism
plays a significant role in superconducting nanograins, it can hardly be
probed by STM-measurements due to the nonzero interlevel spacing, unlike
quantum superconducting nanowires.

##  IV Conclusions and discussion

In conclusion, we have shown that the spatial distribution of the pair
condensate is essentially nonuniform in metallic nanograins. In particular,
the spatially nonuniform pairing can proliferate in nanograins even when
$k_{F}D\sim 300$ and, so, the usual criterion to neglect variations of the
superconducting condensate with position, i.e., $k_{F}D\gg 1$, is not very
useful and can result in wrong conclusions. This is the reason why effects due
to spatially nonuniform pairing in superconducting grains were previously
overlooked. Our study suggests that a new criterion should be based on a more
delicate energy scale (as compared to the Fermi energy), which, in the
superconducting state, is given by the bulk order parameter $\Delta_{B}$. It
turns out that the pairing becomes spatially nonuniform when the interlevel
spacing $\delta$ exceeds $0.1$-$0.2\,\Delta_{B}$. Variations of the order
parameter with position exhibit a pronounced enhancement with an increase of
$\delta/\Delta_{B}$. When $\delta\sim\Delta_{B}$, such variations can be
almost an order of magnitude in highly symmetric grains. At first sight, this
seems impossible because it costs extra energy for such spatial variations.
However, a nonuniform distribution of the pair condensate is accompanied by
enhanced pairing interaction matrix elements, which compensates the energy
cost for an inhomogeneous distribution of the condensate. Another point is the
size-dependent pinning of the chemical potential to groups of degenerate or
nearly degenerate energy levels. Such a pinning plays the role of a filter
that increases the contribution of the single-electron levels in the vicinity
of the chemical potential and suppresses contributions of other states. This
results in an additional mechanism favoring spatially nonuniform pairing in
metallic nanograins.

Figure 6: (a) Single-electron energies $\xi_{i}$ (given in units of the Debye
energy $\hbar\omega_{D}$) ordered in ascending manner versus the ordering
number $N$ for the rectangular-shaped aluminum nanograin with
$L_{x}=7.06\,{\rm nm},\,L_{y}=L_{x}/1.1,\,L_{z}=1.1\,L_{x}$ (squares) and for
a cubic nanograin with $L_{x}=7.06\,{\rm nm}$ (triangles). (b) The size-
dependent excitation gap $\Delta_{E}/\Delta_{B}$ versus $L_{x}$ (in steps of
$\delta L_{x}=0.01\,{\rm nm}$) for an aluminum nanograin of the rectangular
shape with the dimensions $L_{x},L_{y}=L_{x}/1.1,L_{z}=1.1\,L_{x}$: squares
represent the results calculated with the modified matrix elements and with
proper variations of $\mu$; stars are the data obtained for the bulk-like
matrix elements $-g/V$ and $\mu=\mu_{B}$. (c) The spatial distribution of the
pair condensate in the rectangular grain with $L_{x}=7.06\,{\rm nm}$ and
$L_{y}=L_{x}/1.1,\,L_{z}=1.1\,L_{x}$,
$\Delta(x)=\Delta(x,y,z)|_{y=x/1.1,z=1.1x}$.

In this paper we investigated a highly symmetric confining geometry. Due to
this feature the problem becomes effectively one-dimensional (the order
parameter depends only on the radial coordinate) and, so, sufficiently large
diameters up to $D\approx 40\,{\rm nm}$ can be investigated. This size is
almost impossible to reach theoretically for grains with the order parameter
depending on three relevant coordinates due to time consuming numerical
calculations. Such an effectively one-dimensional problem has large
degeneration factors for the corresponding shell structure, resulting in a
significant enhancement of the pairing correlations. In reality there can be
several issues that may lead to a splitting of the shell levels. It will
decrease the degeneration factors and, so, reduce the pairing correlations,
since the main contribution to the sum in the gap equation comes from the
transitions within the same shell pinned to the chemical potential. Among such
issues is the Jahn-Teller deformation, i.e., the transformation of a spherical
nanograin with incompletely filled shells to an ellipsoidal shape. In
addition, the surface imperfections and impurities can significantly change
the distribution of single-electron levels. However, our qualitative results
are quite generic and do not depend on a particular shape of a nanograin and
the presence of possible imperfections. For instance, when
$\delta\sim\Delta_{B}$ the pair condensate will always be spatially nonuniform
because only a few single-electron levels enters the energy interval
$\approx[\mu-\Delta_{B},\mu+\Delta_{B}]$. Due to the dominant contribution of
such levels to $\Delta({\bf r})$, one can expect that the pair condensate
acquires a profile governed by the squared absolute value of the wave function
for the single-electron state closest in energy to $\mu$. This is
significantly strengthened by an increase (in absolute value) of the diagonal
matrix elements $\langle i,\bar{i}|\Phi|i,\bar{i}\rangle$ and, in addition, by
the pinning of the chemical potential to the single-particle levels.

We remark that the diagonal matrix elements, i.e., $\langle
i,\bar{i}|\Phi|i,\bar{i}\rangle$ [see the definition for $\Phi$ below Eq. (3)]
are always enhanced as compared to $-g/V$ in the presence of quantum
confinement, whatever disorder and shape imperfections. This can be seen from
the following simple arguments. Introducing $\varphi_{i}({\bf r})$, the wave
function associated with state $i$, one can write

$\langle i,\bar{i}|\Phi|i,\bar{i}\rangle=-g\int\\!\\!{\rm
d}^{3}r\;|\varphi_{i}({\bf r})|^{4}.$

Due to the normalization condition we have $|\varphi_{i}({\bf
r})|^{2}=\frac{1}{V}+d_{i}({\bf r})$, where $\int{\rm d}^{3}rd_{i}({\bf
r})=0$. Then, the above matrix element can be rearranged as

$\langle
i,\bar{i}|\Phi|i,\bar{i}\rangle=-\frac{g}{V}\Bigl{[}1+V\\!\int\\!\\!{\rm
d}^{3}r\,d^{\,2}_{i}({\bf r})\Bigr{]}.$

The second term in the brackets is always positive in the presence of quantum
confinement, i.e., when $d_{i}({\bf r})\not=0$. It is zero only when
$\varphi_{i}({\bf r})$’s are chosen in the form of plane waves, which results
in $\langle i,\bar{i}|\Phi|i,\bar{i}\rangle=-g/V$.

The above discussion can be supplemented by our numerical results calculated
from the BCS-like equation similar to Eq. (15) but now for aluminum nanograins
of rectangular shape with dimensions
$L_{x},\,L_{y}=L_{x}/1.1,\,L_{z}=1.1L_{x}$. For aluminum we have degen ; fett
; ash : $\hbar\omega_{D}/k_{B}=375~{}{\rm K}$, $gN(0)=0.18$, and
$\mu_{B}=11.67\,{\rm eV}$, which corresponds to the electron density
$n_{e}=181\,{\rm nm}^{-3}$. In Fig. 6(a) single-electron levels arranged in
the ascending order are shown within the Debye window for the rectangular
nanograin with $L_{x}=7.06\,{\rm nm}$ (squares). The same is also given here
for a cubic aluminum nanograin with $L_{x}=L_{y}=L_{z}=7.06\,{\rm nm}$
(triangles). As seen, single-electron levels for the rectangular shape are
distributed in a nearly equidistant manner (with $\delta\approx
0.2$-$0.3\,{\rm meV}\sim\Delta_{B}$) contrary to the states corresponding to
the cubic geometry. It is well-known that an almost equidistant distribution
smit  of single-electron levels near $\mu$ is also expected in the presence of
significant imperfections such as the surface roughness and/or impurities. So,
our results in Fig. 6 give a feeling about the role of the spatially
nonuniform pairing in disordered metallic grains. The excitation energy gap
$\Delta_{E}$ for the rectangular nanograin is shown in units of $\Delta_{B}$
in Fig. 6 as a function of $L_{x}$ in the interval $L_{x}=7$-$8\,{\rm nm}$.
Here squares represent our results calculated with the modification of the
matrix elements and with $\mu$ varying with $L_{x}$; stars are the results
found for the bulk-like matrix elements $-g/V$ and $\mu=\mu_{B}$. As seen,
$\Delta_{E}$ ($\propto T_{c}$) is now two-times enhanced as compared to
$\Delta_{B}$ (on average), which is much less significant than for highly
symmetric grains (compare with Fig. 1) due to a splitting of the shell levels.
However, the effect of interest is still pronounced: $\Delta_{E}$ calculated
for the modified matrix elements and with account of size variations of $\mu$
is generally larger by a factor of $1.5$-$2.0$. The spatial profile of the
order parameter is nonuniform with local enhancements over its average value
by about $100\%$, see, e.g., Fig. 6(c). For rectangular grains with
$L_{x}=7$-$8\,{\rm nm}$ we have $\delta\sim\Delta_{B}$. However, as we
checked, the spatially nonuniform pairing and the related effects of the
modification of the relevant matrix elements and the size-dependent pinning of
$\mu$ are of significance even for smaller $\delta$’s, i.e., when
$\delta>0.1$-$0.2\,\Delta_{B}$ ($L_{x}<14$-$15\,{\rm nm}$). For instance, at
$L_{x}=11\,{\rm nm}$ the order parameter exhibits variations of about
$30$-$40\%$ of its averaged value. These results are in agreement with our
expectations based on the investigation of the highly symmetric spherical
grains.

###### Acknowledgements.

This work was supported by the Alexander von Humboldt Foundation, the Flemish
Science Foundation (FWO-Vl) and the Belgian Science Policy (IAP).

## References

  * (1) R. Parmenter, Phys. Rev. 166, 392 (1968).
  * (2) J. M. Blatt and C. J. Thompson, Phys. Rev. Lett. 10, 332 (1963).
  * (3) B. Mühlschlegel, D. J. Scalapino, and R. Denton, Phys. Rev. B 6, 1767 (1972).
  * (4) R. A. Smith, and V. Ambegaokar, Phys. Rev. Lett. 77, 4962 (1996).
  * (5) A. Perali, A. Bianconi, A. Lanzara and M. L. Saini, Solid State Commun. 100, 181 (1996).
  * (6) F. Braun and J. von Delft, Phys. Rev. B 59, 9527 (1999).
  * (7) G. Sierra, J. Dukelsky, G. G. Dussel, J. von Delft, and F. Braun, Phys. Rev. B 61, R11890 (2000).
  * (8) V. N. Gladilin, V. M. Fomin, and J. T. Devreese, Solid State Commun. 121, 519 (2002).
  * (9) E. A. Yuzbashyan, A. A. Baytin, and B. L. Altshuler, Phys. Rev. B 68, 214509 (2003).
  * (10) Y. N. Ovchinnikov and V. Z. Kresin, Eur. Phys. J. B 45, 5 (2005); ibid, Eur. Phys. J. B 47, 333 (2005).
  * (11) V. Z. Kresin and Y. N. Ovchinnikov, Phys. Rev. B 74, 024514 (2006).
  * (12) A. A. Shanenko and M. D. Croitoru, Phys. Rev. B 73, 012510 (2006).
  * (13) M. D. Croitoru, A. A. Shanenko, and F. M. Peeters, Phys. Rev. B 76, 024511 (2007).
  * (14) A. A. Shanenko, M. D. Croitoru, R. G. Mints, and F. M. Peeters, Phys. Rev. Lett. 99, 067007 (2007).
  * (15) A. A. Shanenko, M. D. Croitoru, and F. M. Peeters, Phys. Rev. B 78, 024505 (2008).
  * (16) M. D. Croitoru, A. A. Shanenko, C. C. Kaun, and F. M. Peeters, Phys. Rev. B 80, 024513 (2009).
  * (17) A. M. García-García, J. D. Urbina, E. A. Yuzbashyan, K. Richter, and B. L. Altshuler, Phys. Rev. Lett. 100, 187001 (2008).
  * (18) I. Giaever and H. R. Zeller, Phys. Rev. Lett. 20, 1504 (1968); H. R. Zeller and I. Giaever, Phys. Rev. 181, 789 (1969).
  * (19) W.-H. Li, C. C. Yang, F. C. Tsao, S. Y. Wu, P. J. Huang, M. K. Chung, and Y. D. Yao, Phys. Rev. B 72, 214516 (2005)
  * (20) W.-H. Li, C.-W. Wang, C.-Y. Li, C. K. Hsu, C. C. Yang, and C.-M. Wu, Phys. Rev. B 77, 094508 (2008);
  * (21) S. Bose, P. Raychaudhuri, R. Banerjee, P. Vasa, and P. Ayyub, Phys. Rev. Lett. 95, 147003 (2005)
  * (22) S. Bose, C. Galande, S. Chockalingam, R. Banerjee, P. Raychaudhuri, and P. Ayyub, J. Phys.: Condens. Matter 21,205702 (2009)
  * (23) D. C. Ralph, C. T. Black, and M. Tinkham, Phys. Rev. Lett. 74, 3241 (1995).
  * (24) C. T. Black, D. C. Ralph, and M. Tinkham, Phys. Rev. Lett. 76, 688 (1996).
  * (25) I. Brihuega, S. Bose, M. M. Ugeda, C.H. Michaelis, and K. Kern, arXiv: 0904.0354, April 2009.
  * (26) S. Bose, A. M. García-García, M. M. Ugeda, J. D. Urbina, C. H. Michaelis, and K. Kern, Nature Mat. 9, 550 (2010).
  * (27) J. von Delft, A. D. Zaikin, D. S. Golubev, and W. Tichy, Phys. Rev. Lett. 77, 3189 (1996).
  * (28) K. A. Matveev and A. I. Larkin, Phys. Rev. Lett. 78, 3749 (1997).
  * (29) Small corrections to $-\frac{g}{V}$ of about $1\%-3\%$ were considered in Ref. gar,  for chaotic grains.
  * (30) Impact of a shell structure on superconducting correlations were previously considered within the grand canonical formalism in Refs. kres01, ; kres02,  and gar, ; bose03, . However, the authors of Refs. kres01 ; kres02  did not investigate the dependence of the relevant matrix elements on quantum numbers (and, so, ignored the spatial dependence of the pair condensate entirely) but took account of the pinning of the chemical potential to incomplete shells. In Refs. gar, ; bose03,  niether modifications of the matrix elements due to quantum confinement nor the quantum-size pinning of the chemical potential were considered.
  * (31) N. N. Bogoliubov, Sov. Phys.-Usp. 2, 236 (1959) [see, also, N. N. Bogoliubov, _Selected Works, Part II, Quantum and Classical Statistical Mechanics_(Gordon and Breach, Amsterdam, 1991)].
  * (32) P. G. de Gennes, _Superconductivity of Metals and Alloys_ (W. A. Benjamin, New York, 1966).
  * (33) A. V. Swidzinsky, _Spatially Inhomogeneous Problems in the Theory of Superconductivity_ (Nauka, Moscow, 1982).
  * (34) P. W. Anderson, J. Phys. Chem. Solids 11, 26 (1959).
  * (35) J. B. Ketterson and S. N. Song, Superconductivity (Cambridge Univ. Press., Cambridge, 1999).
  * (36) A. L. Fetter and J. D. Walecka, _Quantum Theory of Many-Particle Systems_ (Dover, New York, 2003).
  * (37) N. W. Ashcroft and N. D. Mermin, Solid State Physics (Saunders, 1976).
  * (38) Y. Guo, Y. F. Zhang, X. Y. Bao, T. Z. Tang, L. X. Zhang, W. G. Zhu, E. G. Wang, Q. Niu, Z. Q. Qiu, J. F. Jia, Z. X. Zhao, and Q. K. Xue, Science 306, 1915 (2004).
  * (39) D. Eom, S. Qin, M. Y. Chou, and C. K. Shih, Phys. Rev. Lett. 96, 027005 (2006); S. Y. Qin, J. Kim, Q. Niu, and C. K. Shih, Science 324, 1314 (2009).
  * (40) W. Satula, J. Dobaczewski, and W. Nazarewicz, Phys. Rev. Lett. 81, 3599 (1998).
  * (41) J.-P. Martikainen and P. Törma, Phys. Rev. Lett. 95, 170407 (2005).
  * (42) P. G. de Gennes, and D. Saint-James, Phys. Lett. 4, 151 (1963).
  * (43) C. Caroli, P. G. de Gennes, and J. Matricon, Phys. Lett. 9, 307 (1964).
  * (44) A. F. Andreev, Sov. Phys. JETP 22, 455 (1966).
  * (45) A. A. Shanenko, M. D. Croitoru, and F. M. Peeters, Phys. Rev. B 78, 054505 (2008).