# Fixed-𝑁 Superconductivity: The Crossover from the Bulk to the Few-Electron Limit

## The model.—

We model the superconducting grain by a reduced BCS-Hamiltonian which has been used before to describe small superconducting grains  (it was phenomenologically successful for $`d\stackrel{~}{\mathrm{\Delta }}`$ , but probably is unrealistically simple for $`d\stackrel{~}{\mathrm{\Delta }}`$, for which it should rather be viewed as toy model):

$$H=\underset{j=0,\sigma }{\overset{N1}{}}\epsilon _jc_{j\sigma }^{}c_{j\sigma }\lambda d\underset{j,j^{}=0}{\overset{N1}{}}c_{j+}^{}c_j^{}c_j^{}c_{j^{}+}.$$
(1)   
The $`c_{j\pm }^{}`$ create electrons in free time-reversed single-particle-in-a-box states $`|j,\pm `$, with discrete, uniformly spaced, degenerate eigenenergies $`\epsilon _j=jd+\epsilon _0`$. The interaction scatters only time-reversed pairs of electrons within $`\omega _D`$ of $`\epsilon _F`$. Its dimensionless strength $`\lambda `$ is related to the two material parameters $`\stackrel{~}{\mathrm{\Delta }}`$ and $`\omega _D`$ via the bulk gap equation $`\mathrm{sinh}1/\lambda =\omega _D/\stackrel{~}{\mathrm{\Delta }}`$. We chose $`\lambda =0.22`$, close to that of Al . The level spacing $`d`$ determines the number $`N=2\omega _D/d`$ of levels, taken symmetrically around $`\epsilon _F`$, within the cutoff; electrons outside the cutoff remain unaffected by the interaction and are thus neglected throughout.

## Projected Variational Method.—

We construct variational ground states for $`H`$ by projecting BCS-type wave-functions onto a fixed electron number $`N=2n_0+b`$ , where $`n_0`$ and $`b`$ are the number of electron pairs and unpaired electrons within the cutoff, respectively. Considering $`b=0`$ first, we take

$`|0=C{\displaystyle _0^{2\pi }}𝑑\varphi e^{i\varphi n_0}{\displaystyle \underset{j=0}{\overset{N1}{}}}\left(u_j+e^{i\varphi }v_jc_{j+}^{}c_j^{}\right)|\text{Vac},`$ (2)   
where $`|\text{Vac}`$ is the vacuum state. Both $`v_j`$, the amplitude to find a pair of electrons in the levels $`|j,\pm `$, and $`u_j`$, the amplitude for the level’s being empty, can be chosen real  and obey $`u_j^2+v_j^2=1`$. The integral over $`\varphi `$ performs the projection onto the fixed electron pair number $`n_0`$, and $`C`$ is a normalization constant ensuring $`0|0=1`$.

Doing the integral analytically yields a sum over $`\left(\genfrac{}{}{0pt}{}{2n_0}{n_0}\right)`$ terms (all products in (2) that contain exactly $`n_0`$ factors of $`v_jc_{j+}^{}c_j^{}`$), which is forbiddingly unhandy for any reasonable $`n_0`$. Therefore we follow Ref.  and evaluate all integrals numerically instead. Introducing the following shorthand for a general projection integral,

$`R_n^{j_1\mathrm{}j_N}{\displaystyle _0^{2\pi }}{\displaystyle \frac{d\varphi }{2\pi }}e^{i(n_0n)\varphi }{\displaystyle \underset{jj_1\mathrm{}j_N}{}}(u_j^2+e^{i\varphi }v_j^2),`$   
the expectation value $`_0=0|H|0`$ can be expressed as

$`_0={\displaystyle \underset{j}{}}(2\epsilon _j\lambda d)|v_j|^2{\displaystyle \frac{R_1^j}{R_0}}\lambda d{\displaystyle \underset{j,k}{}}u_jv_ju_kv_k{\displaystyle \frac{R_1^{jk}}{R_0}}.`$ (3)   
Minimization with respect to the variational parameters $`v_j`$ leads to a set of $`2n_0`$ coupled equations,

$$2(\widehat{\epsilon }_j+\mathrm{\Lambda }_j)u_jv_j=\mathrm{\Delta }_j(u_j^2v_j^2),$$
(4)   
where the quantities $`\widehat{\epsilon }_j`$, $`\mathrm{\Lambda }_j`$ and $`\mathrm{\Delta }_j`$ are defined by

$`\widehat{\epsilon }_j`$ $``$ $`(\epsilon _j\lambda d/2){\displaystyle \frac{R_1^j}{R_0}},\mathrm{\Delta }_j\lambda d{\displaystyle \underset{k}{}}u_kv_k{\displaystyle \frac{R_1^{jk}}{R_0}},`$   
$`\mathrm{\Lambda }_j`$ $``$ $`{\displaystyle \underset{k}{}}\left(\epsilon _j{\displaystyle \frac{\lambda d}{2}}\right)v_k^2\left[{\displaystyle \frac{R_2^{jk}R_1^{jk}}{R_0}}{\displaystyle \frac{R_1^k}{R_0}}{\displaystyle \frac{R_1^jR_0^j}{R_0}}\right]`$   
$`{\displaystyle \frac{\lambda d}{2}}{\displaystyle \underset{k,\mathrm{}}{}}u_kv_ku_{\mathrm{}}v_{\mathrm{}}\left[{\displaystyle \frac{R_2^{jk\mathrm{}}R_1^{jk\mathrm{}}}{R_0}}{\displaystyle \frac{R_1^k\mathrm{}}{R_0}}{\displaystyle \frac{R_1^jR_0^j}{R_0}}\right].`$   
We obtain an upper bound on the ground state energy and a set of $`v_j`$’s, i.e. an approximate wave function, by solving these equations numerically. To this end we use a formula of Ma and Rasmussen  to express any $`R_n^{j_1\mathrm{}j_N}`$ in terms of $`R_0`$ and all $`R_0^j`$’s, and evaluate the latter integrals using fast Fourier transform routines.

Next consider states with $`b`$ unpaired electrons, e.g. states with odd number parity or excited states: Unpaired electrons are “inert” because the particular form of the interaction involves only electron pairs. Thus the Hilbert subspace with $`b`$ specific levels occupied by unpaired electrons, i.e. levels “blocked” to pair scattering , is closed under the action of $`H`$, allowing us to calculate the energy, say $`_b`$, of its ground state $`|b`$ by the variational method, too. To minimize the kinetic energy of the unpaired electrons in $`|b`$ we choose the $`b`$ singly occupied levels, $`jB`$, to be those closest to the Fermi surface . Our variational ansatz for $`|b`$ then differs from $`|0`$ only in that $`_j`$ is replaced by $`\left(_{jB}c_{j+}^{}\right)_{jB}`$. Thus in all products and sums over $`j`$ above, the blocked levels are excluded (the $`u_j`$ and $`v_j`$ are not defined for $`jB`$) and the total energy $`_b`$ has an extra kinetic term $`_{jB}\epsilon _j`$.

In the limit $`d0`$ at fixed $`n_0d`$, the PBCS theory reduces to the g.c. BCS theory of Ref.  (proving that the latter’s $`N`$-fluctuations become negligible in this limit): The projection integrals can then be approximated by their saddle point values ; since $`\varphi =0`$ at the saddle, the $`R`$’s used here are all equal, thus $`\mathrm{\Lambda }_j`$ vanishes, the variational equations decouple and reduce to the BCS gap equation, and the saddle point condition fixes the mean number of electrons to be $`2n_0+b`$. To check the opposite limit of $`d\stackrel{~}{\mathrm{\Delta }}`$ where $`n_0`$ becomes small, i.e. the f.d. regime, we compared our PBCS results for $`_0`$ and $`_1`$ with MFF’s exact results , finding agreement to within $`1\%`$ for $`n_012`$. This shows that “superconducting fluctuations” (as pairing correlations are traditionally called when, as in this regime, the g.c. pairing parameter vanishes ) are treated adequately in the PBCS approach. Because it works so well for $`d\stackrel{~}{\mathrm{\Delta }}`$ and $`d\stackrel{~}{\mathrm{\Delta }}`$, it seems reasonable to trust it in the crossover regime $`d\stackrel{~}{\mathrm{\Delta }}`$ too, though here, lacking any exact results for comparison, we cannot quantify its errors.

## Ground state energies.—

Figure 1(a) shows the ground state condensation energies $`E_b=_bF_b|H|F_b`$ for even and odd grains ($`b=0`$ and $`1`$, respectively), which is measured relative to the energy of the respective uncorrelated Fermi sea ($`|F_0=_{j<n_0}c_{j+}^{}c_j^{}|\text{Vac}`$ or $`|F_1=c_{n_0+}^{}|F_0`$), calculated for $`N600`$ using both the canonical ($`E_b^C`$) and g.c. ($`E_b^{GC}`$)  approaches. The g.c. curves suggest a “breakdown of superconductivity”  for large $`d`$, in that $`E_b^{GC}=0`$ above some critical $`b`$-dependent level spacing $`d_b^{GC}`$. In contrast, the $`E_b^C`$’s are (i) significantly lower than the $`E_b^{GC}`$’s, thus the projection much improves the variational ansatz; and (ii) negative for all $`d`$, which shows that the system can always gain energy by allowing pairing correlations, even for arbitrarily large $`d`$. As anticipated in , the “breakdown of superconductivity” is evidently not as complete in the canonical as in the g.c. case. Nevertheless, some remnant of it does survive in $`E_b^C`$, since its behaviour, too, changes markedly at a $`b`$ (and $`\lambda `$) dependent characteristic level spacing $`d_b^C`$ ($`<d_b^{GC}`$): it marks the end of bulk BCS-like behavior for $`d<d_b^C`$, where $`E_b^C`$ is extensive ($`1/d`$), and the start of a f.d. plateau for $`d>d_b^C`$, where $`E_b^C`$ is practically intensive (almost $`d`$ independent). The standard heuristic interpretation  of the bulk BCS limit $`\stackrel{~}{\mathrm{\Delta }}^2/(2d)`$ (which is indeed reached by $`E_b^C`$ for $`d0`$) hinges on the scale $`\stackrel{~}{\mathrm{\Delta }}`$: the number of levels strongly affected by pairing is roughly $`\stackrel{~}{\mathrm{\Delta }}/d`$ (those within $`\stackrel{~}{\mathrm{\Delta }}`$ of $`\epsilon _F`$), with an average energy gain per level of $`\stackrel{~}{\mathrm{\Delta }}/2`$. To analogously interpret the $`d`$ independence of $`E_b^C`$ in the f.d. regime, we argue that the scale $`\stackrel{~}{\mathrm{\Delta }}`$ loses its significance – fluctuations affect all $`n_0=\omega _D/d`$ unblocked levels within $`\omega _D`$ of $`\epsilon _F`$ (this is made more precise below), and the energy gain per level is proportional to a renormalized coupling $`\stackrel{~}{\lambda }d`$ (corresponding to the $`1/N`$ correction of  to the g.c. BCS result). The inset of Fig. 1(a) shows the crossover to be quite non-trivial, being surprisingly abrupt for $`E_1^C`$.

## Parity Effect.—

Whereas the ground state energies are not observable by themselves, the parity-dependent spectral gaps $`\mathrm{\Omega }_0=_2_0`$ and $`\mathrm{\Omega }_1=_3_1`$ are measurable in RBT’s experiments by applying a magnetic field . Figure 1(b) shows the canonical $`(\mathrm{\Omega }_b^C)`$ and g.c. $`(\mathrm{\Omega }_b^{GC})`$ results for the spectral gaps. The main features of the g.c. predictions are : (i) the spectral gaps have a minimum, which (ii) is at a smaller $`d`$ in the odd than the even case, and (iii) $`\mathrm{\Omega }_1<\mathrm{\Omega }_0`$ for small $`d`$, which was argued to constitute an observable parity effect. Remarkably, the canonical calculation reproduces all of these qualitative features, including the parity effect, differing from the g.c. case only in quantitative details: the minima are found at smaller $`d`$, and $`\mathrm{\Omega }_0^{GC}<\mathrm{\Omega }_0^C`$ for large $`d`$. The latter is due to fluctuations, neglected in $`_b^{GC}`$, which are less effective in lowering $`_b^C`$ the more levels are blocked, so that $`|_b^C_b^{GC}|`$ decreases with $`b`$.

## Wave functions.—

Next we analyze the variationally-determined wave-functions. Each $`|b`$ can be characterized by a set of correlators

$$C_j^2(d)=c_{j+}^{}c_{j+}c_j^{}c_jc_{j+}^{}c_{j+}c_j^{}c_j,$$
(5)   
which measure the amplitude enhancement for finding a pair instead of two uncorrelated electrons in $`|j,\pm `$. For any blocked single-particle level and for all $`j`$ of an uncorrelated state one has $`C_j=0`$. For the g.c. BCS case $`C_j=u_jv_j`$ and the $`C_j`$’s have a characteristic peak of width $`\stackrel{~}{\mathrm{\Delta }}`$ around $`\epsilon _F`$, see Fig. 2(a), implying that pairing correlations are “localized in energy space”. For the BCS regime $`d<\stackrel{~}{\mathrm{\Delta }}`$, the canonical method produces $`C_j`$’s virtually identical to the g.c. case, vividly illustrating why the g.c. BCS approximation is so successful: not performing the canonical projection hardly affects the parameters $`v_j`$ if $`d\stackrel{~}{\mathrm{\Delta }}`$, but tremendously simplifies their calculation (since the $`2n_0`$ equations in (2) then decouple). However, in the f.d. regime $`d>d_b^C`$, the character of the wave-function changes: weight is shifted into the tails far from $`\epsilon _F`$ at the expense of the vicinity of the Fermi energy. Thus pairing correlations become delocalized in energy space (as also found in ), so that referring to them as mere “fluctuations” is quite appropriate. Fig. 2(b) quantifies this delocalization: $`C_j`$ decreases as $`(A|\epsilon _j\epsilon _F|+B)^1`$ far from the Fermi surface, with $`d`$-dependent coefficients $`A`$ and $`B`$; for the g.c. $`d=0`$ case, $`A=2`$ and $`B=0`$; with increasing $`d`$, $`A`$ decreases and $`B`$ increases, implying smaller $`C_j`$’s close to $`\epsilon _F`$ but a slower fall-off far from $`\epsilon _F`$. In the extreme case $`dd_b^C`$ pair-mixing is roughly equal for all interacting levels.

To quantify how the total amount of pairing correlations, summed over all states $`j`$, depends on $`d`$, Fig. 1(c) shows the pairing parameter $`\mathrm{\Delta }_b(d)=\lambda d_jC_j`$ proposed by Ralph , calculated with the canonical $`(\mathrm{\Delta }_b^C)`$ and grand-canonical $`(\mathrm{\Delta }_b^{GC})`$ approaches. By construction, both $`\mathrm{\Delta }_b^{GC}`$ and $`\mathrm{\Delta }_b^C`$ reduce to the bulk BCS order parameter $`\stackrel{~}{\mathrm{\Delta }}`$ as $`d0`$, when $`C_ju_jv_j`$. $`\mathrm{\Delta }_b^{GC}`$ decreases with increasing $`d`$ and drops to zero at the same critical value $`d_b^{GC}`$ at which the energy $`E_b^{GC}`$ vanishes , reflecting again the g.c. “breakdown of superconductivity”. In contrast, $`\mathrm{\Delta }_b^C`$ is non-monotonic and never reaches zero; even the slopes of $`\mathrm{\Delta }_e^C`$ and $`\mathrm{\Delta }_e^{GC}`$ differ as $`d0`$, , illustrating that the $`1/N`$ corrections neglected in the g.c. approach can significantly change the asymptotic $`d0`$ behavior (this evidently also occurs in Fig. 1.b). Nevertheless, $`\mathrm{\Delta }_b^{GC}`$ does show a clear remnant of the g.c. breakdown, by decreasing quite abruptly at the same $`d_b^C`$ at which the plateau in $`E_b^C`$ sets in. For the odd case this decrease is surprisingly abrupt, but is found to be smeared out for larger $`\lambda `$. We speculate that the abruptness is inversely related to the amount of fluctuations, which are reduced in the odd case by the blocking of the level at $`\epsilon _F`$, but increased by larger $`\lambda `$. $`\mathrm{\Delta }_b^C`$ increases for large $`d`$, because of the factor $`\lambda d`$ in its definition, combined with the fact that (unlike in the g.c. case) the $`C_j`$ remain non-zero due to fluctuations.

Our quantitative analysis of the delocalization of pairing correlations is complimentary to but consistent with that of MFF . Despite being limited to $`n_012`$, MFF also managed to partially probe the crossover regime from the f.d. side via an ingenious rescaling of parameters, increasing $`\lambda `$ at fixed $`\omega _D`$ and $`d`$, thus decreasing $`d/\stackrel{~}{\mathrm{\Delta }}`$; however, the total number of levels $`2\omega _d/d`$ stays fixed in the process, thus this way of reducing the effective level spacing, apart from being (purposefully) unphysical, can only yield indirect and incomplete information about the crossover, since it captures only the influence of the levels closest to $`\epsilon _F`$. Our method captures the crossover fully without any such rescalings.

## Matveev-Larkin’s parity parameter.—

ML  have introduced a parity parameter, defined to be the difference between the ground state energy of an odd state and the mean energy of the neighboring even states with one electron added and one removed: $`\mathrm{\Delta }_{\text{ML}}=_1\frac{1}{2}(_0^{\text{add}}+_0^{\text{rem}})`$. Figure 3 shows the canonical and g.c. results for $`\mathrm{\Delta }_{\text{ML}}`$, and also the large-$`d`$ approximation given by ML, $`\mathrm{\Delta }_{\text{ML}}=d/(2\mathrm{log}(\alpha d/\stackrel{~}{\mathrm{\Delta }}))`$, where the constant $`\alpha `$ (needed because ML’s analysis holds only with logarithmic accuracy) was used as fitting parameter (with $`\alpha =1.35`$). As for the spectral gaps, the canonical and g.c. results are qualitatively similar, though the latter of course misses the fluctuation-induced logarithmic corrections for $`d>d^C`$.

In summary, the crossover from the bulk to the f.d. regime can be captured in full using a fixed-$`N`$ projected BCS ansatz. With increasing $`d`$, the pairing correlations change from being strong and localized within $`\stackrel{~}{\mathrm{\Delta }}`$ of $`\epsilon _F`$, to being mere weak, energetically delocalized “fluctuations”; this causes the condensation energy to change quite abruptly, at a characteristic spacing $`d^C\stackrel{~}{\mathrm{\Delta }}`$, from being extensive to intensive (modulo small corrections). Thus, the qualitative difference between “superconductivity” for $`d<d^C`$, and “fluctuations” for $`d>d^C`$, is that for the former but not the latter, adding more particles gives a different condensation energy; for superconductivity, as Anderson put it, “more is different”.

We would like to thank R. Fazio, G. Falci, A. Mastellone for sending us their numerical data, and K. Likharev, T. Pohjola, D. Ralph, A. Rosch, G. Schön, and A. Zaikin for discussions. This research was supported by “SFB 195” of the Deutsche Forschungsgemeinschaft and the German National Merit Foundation.