{'fraction_non_alphanumeric': 0.028423837093232183, 'fraction_numerical': 0.008509682571371531, 'mean_word_length': 4.351629273504273, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Pairing Correlations in Finite Systems: From the weak to the strong fluctuations regime

The Particle Number Projected Generator Coordinate Method is formulated for the pairing Hamiltonian in a detailed way in the projection after variation and the variation after projection methods. The dependence of the wave functions on the generator coordinate is analyzed performing numerical applications for the most relevant collective coordinates. The calculations reproduce the exact solution in the weak, crossover and strong pairing regimes. The physical insight of the Ansatz and its numerical simplicity make this theory an excellent tool to study pairing correlations in complex situations and/or involved Hamiltonians.


Introduction
The measurements of Black, Ralph and Tinkham [1,2] of discrete level spectra and spectroscopic gaps in nanometer Al isolated grains were interpreted as evidence of the superconductivity phenomenon. To understand the physics of such ultrasmall grains a great deal of theoretical effort was devoted to study such systems starting from grand canonical (BCS) and canonical ensembles [3] as well as very sophisticated theories [4] of the pairing Hamiltonian, see ref. [5] for a review. Later on the exact solution [6] of this naive model Hamiltonian was rediscovered. Some others studies treat the aspect of thermodynamic properties [7,8] while others, see for example [9], the question of persistence of pairing correlations above the BCS critical temperature is addressed. More recently some analytical results in special regimes have been obtained [10]. The main issue of all these studies is the proper description of the crossover between the few electron regime and the bulk one. To analyse this crossover several properties can be computed as a function of the mean electronic level spacing d (or the number of electrons N ) that characterizes the transition from one regime to the other. One of the findings of these studies was that the strong phase transition predicted in a grand canonical study was absent in more advanced theories as well as in the exact solution. In the BCS approach superconductivity is not possible for all d (N ) breaking down at a critical d value. This break down is number parity dependent and indicates that quantum fluctuations are not treated adequately by the BCS wave function. The knowledge of the exact solution for the simple-minded pairing Hamiltonian does not diminish importance to the theoretical approximations developed for the study of that Hamiltonian, see [11] for a review. These approximations are very general and allow the study of more sophisticated Hamiltonians for which no exact solution exists. The use of exactly solvable Hamiltonians [12], on the other hand, is very practical since it allows to check the accuracy of different approximations in the limiting cases represented by those Hamiltonians.
In a recent paper [13] we have proposed a new approach to study superconductivity in finite systems, namely the Generator Coordinate Method (GCM) [14], based on particle number projected BCS wave functions generated in a suitable way. In that paper the GCM approach was applied to superconducting grains described by the Pairing Hamiltonian and it was shown to provide an accurate description of these systems in perfect agreement with the exact Richardson solution. The purpose of this paper is two-fold, first, to present a detailed derivation of the relevant formula as well as the way to solve the Hill-Wheeler (HW) equations and, second, to analyse different generator coordinates in the context of pairing correlations. The derivation presented is comprehensive enough to allow for the application of the formalism to other pairing Hamiltonians. Furthermore since our theory is very general and not constrained by any requirement can be applied to more complex systems. As a matter of fact we have performed preliminary studies with the most general pairing Hamiltonians proposed in [12] and the results [15] are of the same quality as the ones presented in this investigation. Finite temperature effects are not considered in the present study.
In sect. 2 we derive the general formula of the GCM. In sect. 3 we discuss the different coordinates to be used in the calculations. The convergence and other issues concerning the numerical solution of the HW equations is analysed in sect. 4. Finally in sect. 5 the whole formal-ism is applied to study superconducting grains. The paper ends with the Conclusions and some numerical aspects discussed in the Appendices A and B.

Theory
The pairing Hamiltonian used in most calculations is given by where k+ (k−) labels the single particle level (time reversed) with energies ǫ k and c k , c † k destroys and creates electrons in their respective states. The interaction constant G is taken as λ d with d the level spacing and λ the BCS coupling constant whose value for Al is 0.224. The single particle energies ǫ k for simplicity take the values ǫ k = kd. The number N of electrons is equal to the number of levels and in the ground state they form N/2 Cooper pairs, so one works at half filling. This Hamiltonian allows the discussion of the crossover between the strong-coupling regime (d/∆ ≪ 1) that represents large grains and the weak-coupling regime (d/∆ ≫ 1) for small grains, in terms of the quantity d/∆ = 2sinh(1/λ)/N with ∆ the bulk gap, or equivalently in terms of the number of electrons N .
The simplest way to deal with pairing correlations is provided by the BCS theory [16]. Its Ansatz is given by the mean field wave function The variational parameters v k are related to the probability to find two electrons in the level k. The parameters u k are given by u 2 k + v 2 k = 1. The spontaneous particle number symmetry breaking mechanism implicit in Eq. 2 enlarges the available variational Hilbert space making the BCS approximation, in the case of large particle numbers, a very good one. The BCS state (2), on the other hand, undergoes strong particle number fluctuations and for finite systems like metal grains, the Ansatz (2) is unreliable and misses essential features. To correct this failure it is necessary to develop the BCS formalism in a canonical ensemble, where the particle number is fixed, rather than in a grand-canonical one. The restoration of the particle number in the BCS context was introduced by Dietrich and Mang [17] in a nuclear structure context and it was applied for the first time to superconducting grains by J. von Delft and F. Braun [18,19]. The projection method is based on the Anderson formulation of superconductivity [20] where projection onto good particle number is presented as an integration in the gauge variable ϕ, We assume the number of particles N to be even, the odd case is considered in appendix A.1. The formulation of the particle number projection (PNP) can be done in several ways [21]. Very compact formula are obtained in terms of the residuum integrals [17] defined by This definition holds for indices j 1 , · · · , j M such that j k = j p for all k and p. The integer M is simply a counter of the j's involved. In case that two or more indices are equal we define the corresponding residuum integral as zero. All expectation values can be easily calculated in terms of the residuum integrals. As an example we evaluate the matrix element N BCS|BCS N , direct substitution of Eq. 3 provides (5) In the same way the projected energy is given by The PNP energy, as the BCS one, depends only on the variational parameters u k , v k . Minimization with respect to these parameters leads to a set of N coupled non-linear equations The quantitiesǫ k , ∆ k and Λ k are defined bŷ The set of equations (7) resembles the ordinary BCS equations. In that equations, Λ k = 0 andǫ k = ǫ − Gv 2 k − µ. The Lagrange multiplier µ takes care, on the average, of the particle number conservation. Notice that in the projected equations the fields Λ k appear in addition. The solution of Eqs. (7) defines |BCS N . In the literature [18] this is usually called Projected BCS (PBCS) theory. Details of how the set of Eqs. (7) is numerically solved are given in appendix A.2.
To include additional correlations we consider a general superposition of different projected-BCS wave functions, The new wave function |Ψ N is based on the Generator Coordinate Method (GCM) developed by Hill and Wheeler in Nuclear Physics [14]. It has been also used by Peierls, Yoccoz and Thouless [22,23] among others to deal with the restoration of symmetries in mean field approaches as well as to deal with a variational approach to collective motion. It has also provided a variational derivation of the Random Phase Approximation [24]. The coordinate ξ refers to any parameter on which the BCS states may depend parametrically. In this way the superposition state |Ψ N takes care of the fluctuations associated to the parameter ξ. In principle, the variational quantities are the weights f (ξ) and the occupancies u k (ξ), v k (ξ) and should be determined invoking the variational principle. The final equations, however, result in an integro-differential set of equations very complicated to solve. In consequence some assumptions about occupancies are needed in order to facilitate the numerical implementation. If we assume that the quantities u k (ξ), v k (ξ) are known (see below) one deals only with the problem of calculating the weights f (ξ). This is accomplished by the Hill-Wheeler (HW) equation H is the Hamiltonian overlap, defined by (12) and N the norm overlap This equation is very similar to Eq. (5). The residuum integrals have now been generalized by In App. B we discuss the way to calculate these integrals and the Hamiltonian overlap of Eq. (12). It is important to notice that the solution of the HW equation provides not only the ground state but also the low-lying collective states. As we mentioned above, in the HW equation the quantities u k (ξ), v k (ξ) are supposed to be determined beforehand. They are usually fixed by the way the projected wave function |BCS(ξ) N is calculated, namely, whether |BCS(ξ) N is determined by projection after variation (PAV) or variation after projection (VAP). In the former (PAV), the occupancies are determined by the symmetryviolating wave function |BCS , i.e., by solving the ordinary BCS equations. In the VAP case the occupancies are given by the solution of the variational equations Eq.(7). In the BCS framework the VAP approach is known as PBCS. Obviously the VAP method is more involved but it is a fully self-consistent method that provides better results. We shall denote the first method GCMPAV and the second one GCMVAP.

SELECTION OF THE GENERATOR COORDINATE
The generator coordinate ξ is quite general and its selection is motivated by the physical problem. The BCS wave functions depend parametrically on the generator coordinate, its selection is therefore strongly related to the ways we have to characterize the wave function. Though there are many ways to choose the generator coordinate, we think that for the BCS case there are three relevant ones: The gap parameter ∆, the Lagrange parameter µ associated with the particle number of the BCS wave function and ∆N 2 = BCS|N 2 |BCS − BCS|N |BCS 2 , the fluctuations on the number of particles of the BCS wave function.
Instead of using directly the gap parameter as a coordinate it is numerically easier to generate BCS w.f.'s with different gap parameters by solving the corresponding BCS (PBCS) equations for different values G trial of the strength constant G. This method is easy to implement and very efficient. The second method is the simplest one. Now the generator coordinate is the chemical potential µ which in the ordinary BCS equations is used as Lagrange multiplier to fix the mean value of the particle number in the grand-canonical ensemble. In our case we solve the BCS equations for fixed µ and the use of different µ values allows to generate wave functions |BCS(µ) with different average particle number. The fact that |BCS(µ) does not have on the average the right particle number does not matter since later on we project on the right particle number. In this case one is looking for the fluctuations in the position of the Fermi level.
The last method, finally, considers fluctuations around the uncertainty in the particle number ∆N 2 . Now it is necessary to add a constraint to fix a given value of ∆N 2 . This is done by using the modified Hamiltonian H ′ = H −µN −µ 2 ∆N 2 , where the parameter µ 2 guaranties that the constraint is fulfilled. In principle we have set up six variational methods (PAV and VAP versions of each coordinate) but only five are feasible, because the VAP version of µ (by construction) is not possible. Although the ultimate test of the quality of the selection of the generator coordinate will be the eigenstates of the HW equation it is interesting to have a look on the Hilbert space generated by the different coordinates. The diagonal elements of the matrix H ξξ /N ξξ of the HW matrix, Eq. (11), are the projected total energies E N (ξ). This quantity is related to the condensation energy (CE) by E con (ξ) = E N (ξ) − E F with E F the uncorrelated energy of the Fermi sea, i.e., E F = 2 j ǫ j − GN/2. Though, as mentioned above, we have five ways to generate w.f. depending parametrically on ξ we shall concentrate in this section on the three PAV cases corresponding to the three different coordinates under study. In Fig.1 we display E con (ξ) as a function of the corresponding generator coordinate ξ and for different particle (level) numbers to cover the full range from weak to strong pairing regimes. For simplicity we plot only the curves for grains with an even number of particles. Let us first discuss the coordinate G trial . It is obvious that, for each number of particles N , a critical value G c (N ) of G trial exists such that no superconducting solution of the system is found below it. In panel a) we show the CE versus G trial − G c (N ). We find a parabolic behavior with the vertex moving to larger values of G trial − G c (N ) as the particle number decreases (as one would expect). The curves get softer with decreasing particle number with the curve N = 20 being specially soft. We also find that the value of the CE in the minimum is larger (in absolute value) as the particle number increases (as one also would expect).
In panel b) the quantity E con (µ) is plotted against µ − µ BCS (N ), µ BCS (N ) being the chemical potential of the BCS equation for the corresponding case. Because of the particle-hole symmetry of the model the subtraction of µ BCS (N ) provides symmetric curves around µ − µ BCS (N ) = 0 1 . For large N the solutions are parabola like curves which soften with decreasing particle number. For N ≥ 40 we find superconducting solutions for all µ values. This is not the case for N = 20 where for certain µ intervals we do not obtain any solution for the BCS equation, see below for more details. This is not surprising because the standard selfconsistent BCS equation does not provide a correlated solution in this case, see below. Lastly in panel c) E con (∆N 2 ) is plotted against ∆N 2 . Here we also obtain a parabolic behavior similar to the case a) with the difference that the minima shifted to large ∆N 2 correspond to the large particle numbers. The CE gets softer with larger particle number as one would expect. 1 It is important to notice that in this case the numbers 20, 40, etc correspond only to the number of levels and not to the number of particles of |BCS(µ) . Since we work without constraint on the particle number in general BCS(µ)|N |BCS(µ) = N . On the other hand since we are projecting on the particle number the wave functions |BCS(µ) N correspond to a system with N particles. It is clear that the energy minima of the different coordinates provide an approximation to an unconstrained VAP calculation. In Table 1 we have summarized the minima of the parabola as well as the VAP values and the exact ones. We find that all three coordinates do a good job for large particle numbers and that big differences appear for small particle numbers, i.e., in the weakly correlated regime. We find that in general and at this level the co- ordinate G trial is the most effective followed by ∆N 2 and µ. Of course this does not mean very much since the configuration mixing calculations will change these results.
Let's now analyse the wave functions generated with the different coordinates. To a given value of a coordinate, let say ξ 0 , corresponds a wave function |BCS(ξ 0 >. A simple way to characterize the physical content of this wave function is by the associated gap parameter ∆(ξ 0 ) = G k u k (ξ 0 )v k (ξ 0 ). In Fig. 2 we have represented the gap parameter ∆(ξ) associated to each wave function as a function of the coordinate ξ used to generate it. In panel a) we show the results for G trial . Of course the G entering into ∆ is the one of the original Hamiltonian, see Eq. (1), independently of the G trial used in the calculations. Taking into account the expression of ∆ we expect, in first order, a linear behavior with G trial and this is what we obtain. In general a very broad range of gap parameters is covered, which is the reason why the coordinate G trial can be considered equivalent to the gap parameter ∆. The case of the coordinate µ is considered in panel b), where we represent the corresponding gap parameter as a function of µ−µ N BCS . We find an oscillating behavior of ∆ with µ due to the symmetry of the model. Notice that the scale of the y-axis depends on the particle number considered, see the figure caption. For µ = kd, i.e., at the single particle energies ǫ k , we find maxima and for µ = k(d + 1/2) minima. The period and amplitude of the oscillations decrease with growing particle number because in this model d ∼ 1/N . For N ≥ 40 we obtain superconducting solutions for all µ values, in particular for µ = µ BCS , i.e., for the selfconsistent BCS equation. For N = 20, however, we observe that at and around µ = k(d + 1/2) we do not obtain correlated wave functions. As mentioned above this behavior is in agreement with the fact that the selfconsistent BCS solution does not have correlated solutions in this case. The situation is further illustrated in Fig. 3 for the N = 20 case. In the weak pairing regime we only find solutions for µ values corresponding to the hatched regions around a given level. In the region around the level k, the number of particles of the BCS w.f., i.e., the expectation value BCS(µ)|N |BCS(µ) , varies in a continuous way from 2(k − 1) to 2k. For example for the case of N = 20, i.e. k = 10, the BCS w.f. around the level 10 have average numbers of particles ranging from 18 to 20. In general from these w.f.'s it is always possible to project to 20 particles. In the regions between the hatched regions no BCS solution is found but only the Hartree-Fock (HF) one. The numbers of particles are obviously integer num- Finally, in panel c) the gap parameters corresponding to the ∆N 2 generator coordinate are plotted. The behavior is again linear, as for the coordinate G trial , but the range of the gap parameters involved in each wave function is the opposite one. In this case we obtain for small particle numbers a much larger range than for large particle numbers.

NUMERICAL SOLUTION OF THE HILL-WHEELER EQUATION
For our purposes solving the HW equation, Eq. 11, is equivalent to the diagonalization of the Hamiltonian in the nonorthogonal basis of the generator states |BCS(ξ) N .
The usual procedure to deal with this equation [21] involves two diagonalizations. In a first step, the norm overlap N ξξ ′ is diagonalized with the functions u k (ξ) forming a complete orthonormal set in the space of the weights f (ξ). Its eigenvalues are never negative, n k ≥ 0, because the matrix N is definite positive. We shall keep the u k (ξ) with nonzero eigenvalues corresponding to the linearly independent states. In practice and due to numerical reasons we restrict the u k (ξ) to those with eigenvalues larger than a tolerance ε. For each of these functions there exist states |k , called the natural states, which span a collective subspace H C . In a second step the HamiltonianĤ is diagonalized in this space with The HW equations provide a set of wave functions, and energies E σ labeled by the index σ, the lowest one corresponding to the ground state and the others to excited states. In this work we are only interested in the ground state. Taking into account Eq. 10 and Eq. 19 one obtains Since the wave functions |BCS(ξ) N are not orthogonal, the weights f (ξ) cannot be interpreted as the probability amplitude to find the state |BCS(ξ) N in |Ψ N . It can be shown, however, that the functions are orthogonal and that they can be interpreted as probability amplitudes.
For numerical purposes all the expressions above involving integrals have to be replaced by sums discretizing the space ξ. In this form one deals with matrix equations easier to handle. The question that immediately arises is how to determine the optimal ξ-mesh to be used in the calculation. The border values of ξ are determined by energy arguments since the probability of mixing high-lying states is very small. The ξ-coordinate intervals used in the calculations are given in Table 2. The calculations depend furthermore on the mesh step used in the discretization. This parameter is chosen as to optimize the calculations, i.e., we take the largest mesh that includes all states with relevant information. This parameter is also related to the required accuracy. In the calculations performed we have not attempted to reproduce the exact results up to an unusual accuracy. In Fig. 4 the convergence of the condensation energy, as a function of the number of mesh points used in the calculations, is shown for different numbers of particles and for the coordinates G, µ and ∆N 2 . We observe that the number of mesh points needed for convergence depends on the ξ-coordinate and on the number of particles. The coordinate G trial , see top panel, provides the best convergence of the three calculations. For large particle numbers a very good convergence for relatively few mesh points is found. For small numbers of particle one has to go to larger mesh points to find the plateau. For the µ coordinate the situation is the reverse one, i.e., one finds earlier convergence for small numbers of particles. Finally, the situation for ∆N 2 is something in between the two former cases. We find that to reach convergence in energy 80 mesh points are sufficient for all coordinates. This is the number which we will use in all following numerical applications. A further check concerning the convergence is the number of natural states kept in the calculations. Since many of these states are linearly dependent some natural states |k will have a vanishing norm and must be excluded. In the calculation only those natural states with a norm larger than a given tolerance ε are kept. For a given tolerance we take as many states |k as needed to reach a good plateau. Now we analyze the energy convergence as a function of the number of natural states kept in the diagonalization of the HW equation or equivalently of the tolerance of the calculations. This is shown in Fig. 5 for the three coordinates and for grains with different numbers of particles. In the µ and ∆N 2 coordinates we find that for tolerances smaller than 10 −10 linear dependent states are introduced in the calculations providing unrealistic energy values. The interesting point is the nice plateau found for larger tolerances. The tolerance of 10 −10 corresponds, typically, to around 15 linearly independent states. The coordinate G is in this respect somewhat different. One observes that for tolerances of up to 10 −15 one still gets linearly independent states, which obviously correspond to highly excited states that do not affect the energy of the ground state. This tolerance typically amounts to 20 linearly independent states. From this respect we conclude that if one is interested in excited states the coordinate G is more effective than the other ones.
The diagonalization of the HW equation in the case of the µ coordinate requires some comments. As mentioned above in the weak pairing regime, for N = 20 for example, and for several µ intervals one does not find superconducting solutions. This circumstance, as explained above, shows up as "missing points" in the CE curves. These discontinuities do not affect however the solution of the HW equations, since these points have norm zero and do no mix with the other states.

APPLICATION TO SUPERCONDUCTING GRAINS
In this section we present a systematic study of properties of superconducting grains. The results of the GCM calculation of different quantities are compared with the PBCS approximation and the exact Richardson solution.

Ground State Condensation energies
Condensation energies characterize the presence of pairing correlations. The crossover between superconducting and fluctuation dominated regimes can be described through this quantity.
As in the former cases the condensation energy E con is defined as the difference between the total energy in the corresponding approximation and the energy of the uncorrelated Fermi sea. For example in the GCM approaches it is given by E con = E σ=0 − E F , see Eq. 19 and below. This quantity is displayed in Fig. 6 for even grains (up to 600 electrons) and odd grains (up to 601 electrons) as a function of the particle number N . In both plots we give numerical results for the approximations discussed above, BCS, PBCS and the GCMPAV and GCMVAP approaches. The GCMPAV results are presented for the coordinates G trial , µ and ∆N 2 and the GCMVAP for G trial and ∆N 2 . The grand-canonical (BCS) calculation of E con predicts vanishing correlations in the few-electron regime in the even and odd systems. The PBCS condensation energies, on the other hand, though always negative predict an unrealistic sharp crossover between the fluctuation dominated regime and the bulk which is more pronounced in odd grains. This artifact is not present neither in the GCM approaches nor in the exact calculations. The simpler GCMPAV approaches already predict a smooth crossover for odd and even grains. The more involved GCMVAP approaches not only predict a smooth crossover but their predictions coincide with the exact results. Concerning the GCMPAV calculations we find that the µ coordinate is the most effective of all of them followed by the ∆N 2 one. Paradoxically the calculation with the µ coordinate is the simplest one from the numerical point of view.
The reason why the µ coordinate is the most successful one is probably due to the fact that using this coordinate one has the right inertia parameter for the rotations in the gauge space associated with the operatorN . As a matter of fact this was demonstrated by Peierls and Thouless [23] in the context of the translational invariance and a PAV approach by the double projection technique (see Eq.10 above). In this work they show that the right inertia parameter of the collective motion associated with the linear momentum operatorP (N in our case) is obtained when the GCM coordinates are the position (ϕ in our case ) and the velocity (µ in our case). That means the dynamics associated with Eq. 11 has the right inertial parameter. In a VAP approach one always obtain the right mass parameter [21].

Pairing correlations.
In a canonical ensemble the BCS order parameter is identically zero. For this reason it is necessary to define another quantity to characterize pair correlations in a state of fixed numbers of electrons. We choose the pairing parameter used in ref. [18] where the subindex b indicates the number parity of the grain. The quantities C k 's are defined by and form a set of correlators which measure the fluctuations in the occupation numbers. The expectation values are to be calculated with the wave functions of the corresponding approach using the formula developed in Appendices A and B. In an uncorrelated or in a blocked state one has C k = 0. In the grand-canonical case the C k 's reduce to C k = u k v k and ∆ b coincides with the usual superconducting order parameter.
In Fig. 7 we show our results for the pairing parameter in units of∆ for even (upper panel) and odd (lower panel) systems respectively. As we can see in both plots the sharp transition occurring in the BCS and PBCS methods is absent in the GCM approaches as well as in the exact solution. The peculiar behaviors of ∆ b in the exact and GCM approximations before and after the BCS breakdown are related to the change of a pairing delocalized in energy (weak pairing regime) to a localized one (strong pairing regime). The rough decrease of ∆ b with N is connected to the special feature of the model, for which the constant G of Eq. (22) is inverse proportional to the number of electrons. The fact that ∆ b converges monotonically to the final value∆ , in the even case from above and in the odd one from below, is due to the blocking effect. In these plots we observe again that the GCMVAP approaches provide solutions closer to the exact one than the GCMPAV approaches.

Collective wave functions.
We now look at the structure of the GCM states in the space of the collective parameter ξ. The collective weights f (ξ) can not be interpreted as probability amplitudes because the generating states |BCS(ξ) N are not, in general, orthogonal to each other. The amplitudes G(ξ) of Eq. (21), on the other hand, play the role of "collective wave functions", they are orthogonal and their modules squared have the meaning of a probability.
The quantities |G(ξ)| 2 are plotted in Figs. 8, 9 and 10 as a function of the parameters G, µ and ∆N 2 , respectively. The behavior of |G(ξ)| 2 as a function of ξ indicates which are the most relevant components of the states |Ψ N (ξ) in terms of the parameter ξ. To guide the eye we have also plotted in these figures the projected energy E N (ξ) of Eq. (6). For simplicity we restrict our discussion to even systems. In Fig. 8 we represent these quantities for the coordinate G trial in the GCMPAV and GCMVAP approaches. Since the wave functions of both approaches do not differ qualitatively we shall discuss both cases together. The fact that the projected energies are lower in the PAV than in the VAP approach for a given G trial − G c is due to the fact that G c is almost zero for all particle numbers in the VAP approach while it varies considerably with the particle number for the PAV case, see Table 2. We find broad potential energy curves for small particle numbers and narrower ones with increasing N . Interestingly the potential energy curves for the GCMPAV and GCMVAP approaches are rather different for small particle numbers and they become similar for the large ones. Concerning the wave functions for N = 20 we obtain very broad distributions corresponding to a situation of weak  pairing dominated by fluctuations in the order parameter ∆, see also panel a) of Fig. 2. For N ∼ 40 we find that the wave functions are not that extended anymore but they still present a two peak distribution, wit the first peak around the non-superconducting solution and the other around a superconducting one. For larger particle numbers (N ∼ 86, 172, 400) a one peak distribution emerges with the width of the peak getting smaller for increasing particle number. Looking at Fig. 2a we see that at large N the distribution peaks around the wave function with a gap very close to the bulk one.
The results for the µ coordinate are presented in Fig. 9. For N = 20, in the weak pairing region, we obtain a flat potential which shape corresponds to the physics already discussed in relation with Fig. 3. The collective wave function, also according to the discussion of Fig. 2b and Fig. 3, displays an oscillating behavior with maxima around µ = kd and minima (zeroes) around µ = k(d+1/2). The height of the maxima decreases considerably for k values different from 10 and 11, i.e., the collective wave function is mainly formed by the HF solution around the Fermi level and the N=20 components of the BCS solution of the levels above and below the Fermi level. For N = 40 the weak pairing regime persists and the wave function displays a structure similar to the N = 20 case but with the strength much more concentrated owing to the fact that the level spacing decreases with increasing number of particles. For N = 86 the two peak structure is just a reminiscence of the weak pairing situation and for N = 172 and 400 a one peak structure emerges indicating the strong pairing situation. The potential energy curves get steeper with increasing N and the localization of the peak around µ N BCS sharpens in the same way. As one can see in Fig. 2b the range of the ∆ parameter covered by the wave functions diminishes with increasing N .
Lastly we discuss the ∆N 2 coordinate in Fig. 10. The potential energy curves are easy to understand. In the small particle number limit the BCS solution does not provide a superconducting solution and therefore ∆N 2 = 0. On the other hand the BCS approximation provides the exact solution in the bulk limit, i.e. there ∆N 2 >> 1. Accordingly we expect minima in the projected energies at small ∆N 2 for low N and at large ∆N 2 for large N. The potential energy curves get softer with growing N because for increasing ∆N 2 it is energetically easier to change this value. As it should be the potential energy curve for the GCMVAP approach lies below the GCMPAV one. Concerning the collective wave functions, their behaviors correspond to the shape of the potentials. For N ≤ 86 there is a finite probability of having an uncorrelated HF solution as a component of the collective wave functions and only for N ≥ 172 we obtain Wigner-like functions with zero amplitude for the HF component. The range of ∆'s covered by the wave functions can be read from Fig. 2c. Finally, we would like to mention that a very detailed comparison of our wave functions and the exact ones has been made in ref. [13]. We find that the physical content of the GCMPAV wave functions and the exact ones is identical.

Conclusions
In this paper we have presented a detailed formulation of the particle number projected Generator Coordinate Method. We have discussed two different coordinates to generate wave functions for the variation after projection method and three for the projection after variation one. The theory has been applied to study superconducting grains with a pairing Hamiltonian. We have shown that the GCMVAP calculations with both proposed coordinates reproduce the exact results in the weak, crossover and bulk regimes. Concerning the GCMPAV calculations we find that all three proposed coordinates, in spite of not being able to reproduce the exact results, describe qualitatively the correct physics washing out the phase transition found in the BCS and the PBCS approaches. Concerning the degree of accuracy we find that the µ coordinate is the most effective of all the three followed by the ∆N 2 one.
We think that these results are rather general and apply to many more complex Hamiltonians than the naive pairing one considered here. Since the GCM Ansatz includes explicitly fluctuations in the wave function it is very appropriate to deal with finite systems where phase transitions may take place. The method, contrary to other approximations, applies equally well to systems with very few or very large particle number. On the other hand the GCM Ansatz is very versatile to be adapted to other physical situations by considering additional coordinates to the ones discussed in this work.
This work has been supported in part by DGI,

A Peculiarities of the PBCS approximation
In this appendix we discuss some numerical aspects of the solution of the equations used in this article.

A.1 Odd particle number case
The PBCS and HW methods can be extended to systems with an odd number of particles by blocking one of the available states. A system with an odd particle number is described by the state (24) with N an even number and l the blocked state. The residuum integral of Eq. (e −iϕ/2 u 2 k + e iϕ/2 v 2 k ), (25) in an obvious notation. As before all expectation values can easily be calculated in terms of the residuum integrals, for example the norm matrix element is given by and the projected energy by The superindex is somewhat superfluous and it can be suppressed. Blocking different states l one obtains different excited states. The lowest of these energies corresponds to the ground state of the system with an odd particle number. The PBCS and HW variational equations obtained with theses states can be guessed without further calculations: In all sums and products in the equations presented for even systems the term corresponding to the blocked state l has to be excluded. In the same way one can write down states with 3,5,.. etc. blocked states.

A.2 Numerical solution of the PBCS equations.
The set of coupled non-linear equations (7) is usually solved by an iterative procedure. In order to speed up such procedure it is convenient to introduce the variable χ k through the relation where the normalization condition u 2 k + v 2 k = 1 has been taken into account by construction. In terms of χ k the variational equations look like The additional transformation allows to isolate the new variable in terms of the fieldŝ ǫ k , Λ k and ∆ k The right hand side of this equation does not depend explicitly on θ k because the fieldsǫ k , Λ k and ∆ k are independent of θ k . This fact is very useful to solve Eq. (30) by numerical iteration. We start with a guess of θ k (for instance the grand-canonical solution) and solve (30) until the convergence of the energy, to a given tolerance, has been reached.
The variational equations (30) involve the computation of the residuum integrals R ν1...νM n . These integrals can be calculated analytically using the existing closed analytical expression [25]. Their evaluation, however, requires the addition of many terms making the whole computation a very time consuming approach. As the integrals must be computed many times for different sets of u k , v k , the numerical solution of (30) in an efficient way requires the computation of the residuum integrals in a fast and accurate way. For this purpose we have implemented Fast Fourier Transform routines to evaluate the integrals reducing the number of different integrals as much as possible to minimize the computational effort. For this purpose the following two identities can be used. The first one was found by Dietrich et al [17]. It can be shown that the residuum integrals satisfy the following recursion relation This formula allows to calculate all residuum integrals if the integrals R 0 0 and R j 0 are known. This relation reduces to N + 1 the overall number of numerical integrations for a given set of v j 's and u j 's.
In the PBCS one only needs Ma's relation to calculate three terms: R jk 1 ,R jkl 2 − R jkl 1 and R jk 2 − R jk 1 . If R 0 0 and R j 0 are known, R j 1 can be obtained by Dietrich's recursion relation. First we consider R jk m with m = 1 or m = 2. Ma and Rasmussen's formula reduces to where we have used the identity ζ j = 1 v 2 j −1. The difference R jk 2 − R jk 1 can be written in a simplified way as follows The calculation of R jkl 2 − R jkl 1 is a bit more complicated, the result is Since the indices of the residuum integrals can be permuted, R jkl 2 − R jkl 1 can be expressed for all possible combinations of v j , v k , v l by the equations above.

B The generalized residuum integrals
The Hamiltonian overlap H ξξ ′ , Eq. (12), can be calculated using the generalized Dietrich's recursion relation, The residuum integrals R i 1 , R ij 1 needed to calculate H ξξ ′ can be written in terms of R j 0 and R jk 0 using Dietrich's relation twice .