{'fraction_non_alphanumeric': 0.035635399793844794, 'fraction_numerical': 0.012928876454130466, 'mean_word_length': 4.192842942345925, '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': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Number-conserving theory of nuclear pairing gaps: a global assessment

We study odd-even mass staggering of nuclei, also called pairing gaps, using a Skyrme self-consistent mean-field theory and a numerically exact treatment of the pairing Hamiltonian. We find that the configuration-space Monte Carlo method proposed by Cerf and Martin offers a practical computational procedure to carry out the numerical solutions in large-dimensional model spaces. Refitting the global strength of the pairing interaction for 443 neutron pairing gaps in our number-conserving treatment, we find the correction to the pairing correlation energies and pairing gaps to have rms values of 0.6 MeV and 0.12 MeV, respectively. The exact treatment provides a significant improvement in the fit to experimental gaps, although it is partially masked by a larger rms error due to deficiencies in other aspects of the theory such as the mean-field energy functional.


I. INTRODUCTION
Computer resources now make it possible to test theories of nuclear structure using the entire body of nuclear data. One particular aspect of nuclear structure is pairing, which is important for determining stability and dynamical properties of nuclei. The Bardeen-Cooper-Schrieffer (BCS) theory [1] has been a paradigm for treating nuclear pairing, but it is not well justified in finite nuclei. Besides its violation of particle-number conservation, the condensate may collapse in finite systems. A recent global study of nuclear pairing gaps [2] found that ∼ 25% of nuclei lacked a BCS pairing condensate because of the weakness of the interaction or a low singleparticle density of states. The observed smoothness of nuclear binding energies calls for a theory that does not force a discontinuous jump between ground states with and without pairing condensates. In Ref. 2 it was found that a small but significant overall improvement in theory could be achieved by using the Lipkin-Nogami (LN) extension of BCS to correct for particle-number violation [3]. However, the LN treatment has its own limitations. For example, it becomes inaccurate near closed shells when implemented in the usual way [4,5]. On a practical level, iterative BCS-LN solvers often have convergence problems near closed shells. We note that there are many methods other than the LN extension of BCS to treat the pairing interaction more accurately, including direct diagonalization in truncated spaces [6]. For methods that emphasize particle number conservation, see Ref. [7] and references therein.
Here we address the question of the importance of a better treatment of pairing correlations by carrying out a global survey using a numerically exact technique to calculate pairing correlation energies at fixed particle number. In particular, we employ the configuration-space Monte Carlo (CSMC) algorithm of Cerf and Martin [8]. Numerically exact solutions can also be obtained by direct diagonalization of the pairing Hamiltonian in config-uration spaces of fixed seniority [9,10], but the CSMC method is more efficient and can be implemented in much larger spaces.
Since our aim is to assess the relative performance of the theory with and without an exact treatment of pairing, we will avoid introducing extraneous elements and closely follow the methodology of Ref. [2]. In that work the performance of various self-consistent meanfield (SCMF) methods was tested on the neutron and proton pairing gaps for odd-A nuclei; here we use the same 443 odd-neutron gaps to assess the importance of an exact treatment of the pairing interaction.
The neutron pairing gap for an (odd) neutron number N is defined by the second-order energy difference in neutron number ∆ (3) o (N ) = − where E(N ) is the ground-state energy of the nucleus with N neutrons and Z protons. The proton number Z is the same for all three nuclei in Eq. (1) and is not indicated explicitly in the formula. As a prototype of an SCMF theory we use the energy density functional constructed from the SLy4 Skyrme functional [11] for the normal density part and a densitydependent contact interaction for the pairing part. The Hartree-Fock+BCS (HF+BCS) equations are solved using the EV8 code [12]. We construct a pairing Hamiltonian whose single-particle orbital energies and pairing matrix elements are extracted from the SCMF calculation of Ref. [2]. Next, we solve this Hamiltonian exactly using the CSMC method, which is free of a sign problem when all pairing matrix elements are attractive. The SCMF interaction energies are also taken from the calculations of Ref. [2]. The main difference here is in the treatment of pairing correlation energies and the strength of the pairing interaction. The performance of the theory is measured by the root-mean-square (rms) of the residuals with respect to the experimental data set after making a least-squared fit of the overall pairing interaction strength. The differences between the Monte Carlo treatment and the BCS approximation are used to estimate the importance of a particle-number-conserving exact treatment of pairing.
The outline of this paper is as follows. In Sec. II we discuss our methodology of constructing a pairing Hamiltonian from the SCMF results and how we use its exact solution to obtain an improved estimate of the pairing gap. In Sec. III we describe the CSMC method used to solve the pairing Hamiltonian. In Sec. IV we present our results for the pairing gaps. Our conclusions are given in Sec. V.

II. METHODOLOGY
The leading approach in the search for a computationally tractable theory of nuclear structure starts with an SCMF theory to construct a set of configurations and then mixes the configurations through a residual interaction to restore broken symmetries and add a correlation contribution to the total energy. By itself, the mean-field theory is straightforward. However, there are different ways to introduce correlations, even if we limit ourselves to pairing correlations. The BCS treatment is the simplest way to introduce pairing correlations and can be easily implemented once the single-particle wave functions and energies have been obtained from the meanfield theory. The Hartree-Fock-Bogoliubov (HFB) approximation is an extension that is required when the orbital properties depend on pairing, but like BCS it violates particle-number conservation. To gain the benefit of the HFB approximation, the pairing Hamiltonian must be defined in very large model spaces and with a more general interaction than can be treated with the present CSMC method. Here we construct pairing Hamiltonians in which the orbitals are fixed from the mean-field calculation, as in the BCS approximation, and have interaction matrix elements that are all attractive. The HFB might be required in the dripline region. However, since very few of the known experimental gaps are in this region, our conclusions should apply to the vast majority of nuclei for which data exist.
For construction of the pairing Hamiltonian, we follow closely the treatment of Ref. [13] as implemented in EV8. The single-particle energies are taken directly from the eigenvalues of the single-particle SCMF Hamiltonian and the interaction is chosen as a density-dependent contact interaction together with an energy cutoff factor described below. In Eq. (2) ρ 0 = 0.16 fm −3 is the conventional saturation density of nuclear matter and the parameter η controls the specific density dependence. We will use η = 0.5, called the "mixed" density-dependent pairing interaction. The strength V 0 is determined by minimizing the rms of the residuals of the calculated pairing gaps from their experimental counterparts [2]. As implemented in EV8, the mean field is invariant under time reversal and the self-consistent single-particle orbitals appear in degenerate time-reversed pairs i and i with energy ε i . The total number of orbital pairs is Ω. The antisymmetrized pairing matrix elements V ij are taken to be where V is given by Eq.
(2) and f i are energy cutoff factors [13] with a = 5 MeV and b = 0.5 MeV. Denoting the singleparticle orbitals by φ i (r, σ), we have Next we construct the pairing Hamiltonian (6) Note that the Hamiltonian (6) does not include diagonal matrix elements V ii . We assume that they have already been incorporated into the mean-field part of the energy density functional.
An Hamiltonian of the form Eq. (6) was used in a recent study comparing the BCS approximation with exact matrix diagonalization results in model spaces of size Ω = 16 [14]. However the sizes of the model space required for a global survey are prohibitively large for direct matrix diagonalization methods to be practical. We therefore use the CSMC method which scales much more gently as a function of Ω. This method can be used to find the exact ground-state energy E CSMC of the Hamiltonian in Eq. (6) to within a statistical error. We note that the pairing Hamiltonian can be solved algebraically for special forms of the interaction following Richardson's method [15], but these are not applicable to more general interactions such as in Eq. (2).
Our improved estimate for the total ground-state energy is given by where E SCMF is the SCMF energy calculated with SLy4 plus the density-dependent contact pairing interaction (2), and E BCS is the BCS ground-state energy of the Hamiltonian in Eq. (6). In Eq. (7) we are essentially replacing the BCS energy of the Hamiltonian in Eq. (6) with its exact CSMC ground-state energy. Both E SCMF and E BCS are calculated with an interaction strength V 0 determined by minimizing the rms deviation of the SCMF gaps from the experimental gaps. However, E CSMC is calculated with a renormalized strength, determined by minimizing the rms residuals (with respect to experiment) of the theoretical gaps calculated from Eq. (7) and Eq. (1). One problem of using a pairing Hamiltonian from a theory such as the one discussed in Ref. [2] is that there are diagonal interaction matrix elements both in the mean-field part as well as in the pairing part of the energy functional. Since the pairing interaction is added to describe correlations beyond those obtained in the SCMF with a single Slater determinant, it should not add diagonal interactions beyond the mean field. At an extreme, if the BCS condensate collapses, the BCS correlation energy should be zero. Because of these considerations we do not include diagonal interaction matrix elements V ii in our Hamiltonian (6) for either the CSMC or the BCS calculations. Such matrix elements remain, however, in the pairing part of the SCMF theory.

III. CONFIGURATION SPACE MONTE CARLO SOLVER
The CSMC Hamiltonian solver has been applied to individual isotope chains [16] but our work here is its first use in a global survey. To introduce the various parameters of the method and make our presentation self-contained, we review the algorithm in Sec. III A. In Sec. III B, we discuss the statistical Monte Carlo error and demonstrate the computational scaling properties of the method.

A. The Monte Carlo method
In the following, we assume the particle number N to be even. The algorithm applies to Hamiltonians of the form in Eq. (6) for which all pairing interaction matrix elements satisfy V ij ≤ 0. The overall operation of the algorithm is similar to many other Monte Carlo methods where a trial state |Φ is evolved in imaginary time whereĤ is the system's Hamiltonian and E t is an energy parameter adjusted to keep the normalization of |Ψ(τ ) approximately fixed. The initial evolution filters out the ground-state component of the initial trial state, and subsequent evolution is used to obtain better statistics for the ground-state energy. The remaining details of the Monte Carlo method relies on the representation of Ψ(τ ) as a superposition of pure paired configurations [8]. Let us label the fully paired eigenstates ofĤ 1 , in Eq. (6), by |n = |n 1 , n 2 , . . . n Ω where n i = 0 or 1 is the pair occupation number of the i-th two-fold degenerate level, and The wave function |Ψ(τ ) can be written as a linear combination of these paired configurations |n The coefficients α n (τ ) can be chosen to be all positive and normalized as n α n (τ ) = 1 .
Thus, the wave function |Ψ(τ ) can be represented by an ensemble of paired configurations |n that are distributed with probability α n (τ ). The evolution in imaginary time is carried out as a series of time evolutions, each of which is over a small time step ∆τ Using the Suzuki-Trotter symmetric decomposition [17][18][19], we write the short-time propagator as We are interested in calculating the matrix elements of Eq. (14) between two paired configurations |n and |n ′ . SinceĤ 1 is diagonal in the |n basis, the only non-trivial part is the matrix elements of e −Ĥ2∆τ . Expanding this propagator in a Taylor series, we have [8] where ω = N 2 (Ω − N 2 + 1) and ν = ωV ∆τ defines a dimensionless time step withV = ij V ij /Ω 2 . Here W κ (n → n ′ ) represents the weight of a path of L pair hops that takes the configuration n to n ′ . Each pair hop describes the transition of a pair of particles from an occupied two-fold level to an empty two-fold level. The probability to have L pair hops in the time interval ∆τ is a Poisson distribution P (L) = e −ν ν L /L!. The parameter ω represents the total number of possible pair hops and ν is the average number of pair hops in the time interval ∆τ . The weight W κ (n → n ′ ) is given by where (i m ,ī m ) and (j m ,j m ) are the orbital pairs whose occupations are swapped at the m-th step of the L-step pair hop process.
In practice, we carry out the Monte Carlo evolution as follows. We take the initial state |Φ to be the groundstate configuration ofĤ 1 and replicate it N e times to generate the initial ensemble. Subsequently, the members of the ensemble are evolved independently. For each time step ∆τ , the time evolution is done stochastically using Eqs. (15) and (16). The number of pair hops L is drawn from the Poisson distribution P (L), and an L-step pair hop process is carried out. At each step we choose an occupied pair orbital (i m ,ī m ) from a uniform distribution, and an orbital (j m ,j m ) which is either unoccupied or equal to (i m ,ī m ), again from a uniform distribution. The occupation numbers of (i m ,ī m ) and (j m ,j m ) are swapped. The resulting new configuration is replicated stochastically with a weight of We adjust the normalization energy E t during the time evolution to keep the ensemble size stable. At the k-th time step we define where N e (k) is the size of the ensemble at the k-th time step. We repeat the above process N τ times for the initial evolution. The resulting ensemble of fully paired configurations n m (m = 1, . . . , N e ) can be expressed as the wave function and is our first representative of the ground-state ensemble.
The ground-state energy E can be calculated from the ground state Ψ using E = n ′ n ′ |H|Ψ (where we have used n ′ n ′ |Ψ = 1). Approximating Ψ by Ψ 1 in Eq. (18), we estimate the ground-state energy to be where E sp (n) is given by Eq. (10) and The prime on the summation in Eq. (20) denotes that the sum is restricted to those combinations ij where the orbital pair (i,ī) is occupied in |n and the orbital pair (j,j) is either unoccupied in |n or the same as (i,ī). Additional representatives |Ψ i of the ground-state wave function are generated by evolving the ensemble an additional number of time steps N T and taking a representative every N c steps to ensure uncorrelated ensembles. Using relations similar to Eq. (19), we obtain . , E NE for the groundstate energy. The final estimate for the CSMC groundstate energy is and its corresponding statistical error is So far, we have discussed a system with an even number of particles N . The generalization to an odd number N o is straightforward. We put a single particle in one of the orbitals of a degenerate pair. This pair of orbitals becomes effectively blocked, i.e., it cannot participate in the pair transitions between orbitals. The energy of the remaining N o − 1 particles is found by applying CSMC to the reduced space in which the blocked orbital pair is excluded. The total energy of the N o -particle system is then given by where ε b is the single-particle energy of the blocked orbital and E b (N o − 1) is the energy of N o − 1 particle system in the reduced space. The ground-state energy of the odd-N system is found by minimizing Eq. (23) over different choices of the blocked orbital b.
For the calculations in this work, we have taken N τ = 5000, N T = 50, 000, N c = 500 and N e (0) = 25, 000. This gives N E = 100 estimators E i of the the ground-state energy and its error in Eqs. (21) and (22), respectively.

B. Statistical error
The statistical error σ in the CSMC energy estimate can be written as where σ 2 in is the intrinsic variance of the energy, i.e., the variance of the quantity E sp (n) + E v (n) for paired configurations n that are distributed according to α n . Here N E = N T /N c is the number of uncorrelated ensembles of size N e used in the CSMC calculation. The replication process described in Sec. III A introduces correlations between configurations in the ensemble at a given time step and N e /χ e (with χ e > 1) represents the effective number of uncorrelated configurations.
In the following, we provide an estimate for the intrinsic standard deviation σ in assuming a constant pairing interaction. In this case E v (n) in Eq. (20) is a constant, and σ 2 in is determined solely by the variance of the E sp (n). With E sp (n) given by Eq. (10), its average valueĒ sp over the various configurations in the ensemble is wheren i = n n i α n . 1 Eq. (25) holds for any constant λ but we choose λ to be the chemical potential to minimize the particle number fluctuations.
In the Appendix we use the BCS wave function to estimate the fluctuations in E sp [see Eq. (A.4)]. For a uniform single-particle spectrum with a bandwidth E c ≫ ∆ (∆ is the BCS pairing gap), we find We can use this expression to estimate the scaling of σ in with the size Ω of the single-particle space. For weak to moderate pairing, ∆ ∝ Ω, and σ in ∝ Ω 3/2 . Our simple estimate (26) is accurate to within a factor of ∼ 2 (see Fig. 8 in the Appendix). We have checked that even in cases when the single-particle spectrum is non-uniform and the pairing interaction is orbital-dependent (e.g., the nuclear pairing Hamiltonians used in this work), expression (26) (where ∆ is taken to be an average pairing gap) provides a reasonable estimate for the intrinsic error in E sp . If all N e configurations of the ensemble at a given time step were to be uncorrelated, the CSMC error would have been σ in / √ N E N e . Since these configurations are correlated in the CSMC calculation, the actual statistical error is larger by a factor of √ χ e [see Eq. (24)]. In Fig. 1 we show (solid circles) this enhancement factor √ χ e (as determined empirically from the CSMC statistical error for ν = 0.1) versus Ω for a uniform single-particle spectrum with level spacing of 1 MeV at half filling (N = Ω) and a constant pairing strength of V ij = 0.3 MeV. The dashed line is a fit to χ e ∼ 1 + (Ω/Ω 0 ) 3 . In general we find that the scaling of χ e with Ω depends on the strength of the pairing interaction.
To illustrate the scaling of the CSMC computational time with the size Ω of the single-particle space, we consider the same example as in Fig. 1. The results, shown 1 Note thatn i differs from the quantum mechanical expectation value of the pair occupation operatorn i . The lower panel of Fig. 2 shows the statistical error calculated from Eq. (22). It appears to scale as Ω 3 (dashed line) 2 . Spaces as large as Ω = 30 are easily computed and we shall argue below that the accuracy achieved is adequate for our purposes. In contrast, if the calculations were done by conventional matrix diagonalization, one would have to deal with a matrix of dimension 1.6 × 10 8 and Ω = 50 (corresponding to a matrix of dimension 1.3 × 10 14 ) would be completely out of reach.
The CSMC calculations for the global survey were carried out using two values of ν (0.05 and 0.10) and averaging their respective energy estimates. We checked that the corresponding time steps were sufficiently small to avoid a significant systematic error from the Suzuki-Trotter decomposition in Eq. (14). To test for other biases in the CSMC algorithm, we compared with the matrix diagonalization of the Hamiltonian (6) for two cases presented in Ref. [14], namely 118 Sn and 206 Pb. The pairing Hamiltonians were obtained from Ref. [2]; H 1 is derived from SCMF with the Skyrme SLy4 energy functional andĤ 2 is of the contact form as in Eq. (2) with η = 0 (and no cutoff factors f i ). The single-particle space in these examples has a size Ω = 16, which requires matrices of dimension ∼ 13, 000 for the direct diagonalization. The calculated correlation energies (measured relative to the HF ground-state energy) are shown in Ta   For the CSMC calculation in our global survey, it is important that the Monte Carlo statistical error does not degrade the accuracy of the calculated pairing gaps to a point where the performance measure would be affected. The maximal permissible statistical error is estimated as follows. We take a typical rms of the residuals (deviations between theory and experiment) in the range 0.25-0.30 MeV, and demand that the Monte Carlo statistical contribution, calculated in quadratures, be less than 0.01 MeV. This requires that the average statistical error for the pairing gap be smaller than √ 0.25 2 − 0.24 2 ≈ 0.07 MeV. In fact, with our choices of the numerical parameters in the global calculation, the maximal statistical error (∼ 0.05 MeV) satisfies this upper bound for all cases.

IV. RESULTS
For our global survey of odd neutron gaps, we take the same nuclei as in Ref. [2], where 443 odd neutron pairing gaps were calculated and compared with experiment. Our procedure for obtaining a new set of theoretical gaps involves the following steps: 1. We start with the full SLy4+pairing energies as calculated in Ref. [2] and construct the pairing Hamil-tonianĤ in Eq. (6) using the converged SCMF single-particle energies and wave functions. The diagonal interaction matrix elements V ii are not included in Eq. (6).
2. We calculate the BCS ground-state energy ofĤ, taking the same interaction strength V 0 as in the original calculations (i.e., V 0 = 700 MeV-fm 3 ), and subtract it off the total SCMF energy.
3. We calculate the exact ground-state energy ofĤ by CSMC and add it back to to obtain our new estimate of the ground-state energy. The overall interaction strength is renormalized in the CSMC calculation, and its value is determined by minimizing the rms residuals of the newly calculated pairing gaps.
4. To make a fair comparison with the BCS, we repeat step 3 with BCS energy (excluding diagonal interaction matrix elements as in the CSMC), refitting the overall strength of the interaction to minimize the rms residuals. 3 Some remarks are in order regarding our refit. The CSMC calculations were performed at two different values of V 0 (560 and 700 MeV fm 3 ). We used a linear interpolation to obtain the ground-state energy for interaction strengths between these two values. The new value of V 0 is determined by minimizing the rms of the residuals using a linear least-squared fit. The model space for our CSMC calculation consists of all orbitals for which f 2 i > 0.01 [see Eq. (4)]. We verified the convergence of our calculations by repeating them for a model space with f 2 i > 0.001. The largest model space used in these calculations is Ω = 64 for the nucleus N = 156 and Z = 100. For this nucleus each CSMC calculation takes about 30 minutes on a single processor.
The results for our refits are shown in Table II along with the fit reported in Ref. [2]. First we note that the fitted value of the interaction strength V 0 is smaller for the CSMC gaps than its fitted value for the BCS gaps by about 6%. It is not surprising that the required strength is higher in a theory (e.g., BCS) that is subject to pairing collapse and gives a zero correlation energy in some of the nuclei. In fact, the differences between correlation energies comparing the CSMC and BCS can be quite large; their rms difference is ∼ 0.6 MeV for the more than 900 nuclei in our data set. However, the observable quantities are not the correlation energies but the pairing gaps, for which the differences are much smaller. The rms of the differences between the CSMC and the BCS pairing gaps is 0.12 MeV. Given that the total rms residuals of the theory with respect to experiment is of the order 0.25-0.30 MeV, this difference between CSMC and BCS appears to be quite significant. However, one must realize that when there are independent sources of error, the larger ones can effectively mask the others. This can be seen in the third column of table II, reporting the rms residuals of the CSMC and BCS pairing gap with respect to experiment. The corresponding values, 0.28 and 0.24 MeV, only differ by 0.04 MeV. However, this is close to what one would expect when adding in quadratures the error of the BCS approximation (0.12 MeV) and the other sources of error (0.24 MeV). We also note that the values for the fitted strength and the rms residuals reported in Ref. [2] are somewhat higher than their corresponding values in our BCS fit, to which it should be compared.  To illustrate the performance of the theory locally, we compare in Fig. 3 the theoretical and experimental pairing gaps for the chain of Sn isotopes (Z = 50). All three theories (SCMF, BCS and CSMC) overestimate the gap but follow correctly its overall dependence on neutron number, including the the dip at N = 65 and the sharp drop near the N = 82 shell closure. When compared with the experimental gaps, the CSMC shows a modest but systematic improvement over the SCMF and BCS theories, except for the N = 77 − 81 nuclei in the vicinity of the shell closure. It would be useful to know whether there are any systematic criteria for identifying nuclei for which the improved treatment of pairing has the most benefit. One criterion could be the magnitude of the error (i.e., residual) comparing the SCMF or BCS pairing gaps with their experimental values. To examine the dependence on the SCMF error we take subsets of gaps whose SCMF residual is larger (in absolute value) than some given value and calculate the rms of the subset as a function of the lower cutoff. The results are shown in Fig. 4. The rms error of the subset increases with cutoff but in the CSMC approach it does so at a lower rate than in the BCS treatment. For example, when we keep only those nuclei whose SCMF rms residual is greater than 0.5 MeV, we find that the BCS rms increases to 0.68 MeV while the CSMC rms is only 0.52 MeV, an improvement of 0.14 MeV. The inset of Fig. 4 shows the rms of the CSMC correction to the BCS residual versus the lower cutoff of the SCMF residual. This rms exhibits a gradual increase from about 0.12 MeV when all the nuclei are included to about 0.19 MeV when we include only those nuclei whose SCMF residual is greater than 0.5 MeV. Thus, there is a mild increase in the benefit derived from the exact treatment when the residual error is large.
To narrow further the conditions under which an exact treatment is beneficial, we go back to the symmetries that are broken in mean field and BCS theory, namely particle-number conservation and rotational symmetry. A measure of the BCS violation of particle-number conservation is given by where v 2 i are the BCS occupation numbers. We divide the nuclei with odd number of neutrons into bins of width 1 according to their particle-number fluctuation ∆N (nuclei with ∆N = 0 have their own). Fig. 5 shows the rms of the residuals for the nuclei in each bin versus the midpoint of the bin. The bin with ∆N = 0 consists of all nuclei for which the BCS pairing has collapsed. Clearly the CSMC treatment is needed in that situation. The CSMC also gives an improvement for the bin centered at ∆N = 3.5 having the strongest pairing condensate. This is likely due to the too-large pairing strength V 0 required for the global BCS fit. Thus, when compared to the exact CSMC results, the BCS approximation seems to be adequate when 1 ≤ ∆N ≤ 3. If we use in the BCS treatment the lower value of the CSMC interaction strength, we find the BCS performance to improve gradually with increasing ∆N and to become comparable to the CSMC performance for ∆N ≥ 2.
The violation of rotational symmetry in a nucleus is often characterized by the mass quadrupole deformation parameter β 2 defined by where A, r 2 and Q 0 are, respectively, the mass number, rms radius and intrinsic quadrupole moment of the nucleus. We divide the odd-N nuclei into bins of width 0.1 according to their deformation β 2 in the SCMF treatment. The rms of the residuals are calculated for the nuclei in each bin and their values (in both BCS and CSMC) are plotted versus the bin centers in Fig. 6. We observe that there is almost no difference between the two treatments for oblate nuclei. For spherical and strongly deformed prolate nuclei, the CSMC gives a moderate improvement over BCS while for moderately deformed prolate nuclei there is a significant improvement in the CSMC method as compared with the BCS approximation.
In Fig. 7 we show the average value of the ratio between the CSMC correction to the BCS correlation energy and the CSMC correlation energy, (δE CSMC − δE BCS )/δE CSMC , versus ∆N (here δE BCS ≡ E HF −E BCS and δE CSMC ≡ E HF − E CSMC ). At ∆N = 0 the BCS solution collapses to the HF solution and this ratio is just 1. We observe the above ratio to decrease monotonically versus ∆N as the BCS approximation becomes better.

V. CONCLUSION
Starting from an SCMF theory of the pairing gaps and treating pairing correlations exactly beyond the BCS approximation (with a renormalized pairing interaction strength), we found a significant improvement in the theory as measured by the rms residuals of the pairing gaps. The exact calculations were carried out by constructing a pairing Hamiltonian from the SCMF output and using the configuration space Monte Carlo (CSMC) method.
We find the improvement in the rms residuals of the pairing gaps to be most significant in nuclei for which the BCS condensate was weak, as measured by the smallness of particle-number fluctuation ∆N . Based on our results, the BCS seems to be adequate if one limits the theory to nuclei for which 1 ≤ ∆N ≤ 3. The artificially high value of the BCS interaction strength in a global fit leads to pairing gaps that are on average too large in nuclei with ∆N ≥ 3. We also found the improvement to be larger in moderately deformed (β 2 ∼ 0.1 − 0.2) prolate nuclei.
The total residual in the SCMF treatment of pairing gaps can be thought of as coming from two parts, the inadequacy in the mean field and the approximate treatment of pairing. Since in the CSMC method the pairing part is treated exactly, the CSMC residuals are wholly due to the inadequacy of the mean field. The rms of the residuals of the pairing gaps in the BCS treatment over their values in the CSMC treatment is then an estimate of the error involved in the BCS approximation. We find this rms value to be ∼ 0.12 MeV (see inset of Fig. 4), and propose it as the bound on the accuracy that can be achieved in an SCMF theory that treats pairing correlations approximately.
From the computational perspective, the most notable aspect of this work is the use of the CSMC algorithm, which has not been used previously in a global survey. where ∆ i are the level-dependent BCS pairing gaps. The BCS estimate (A.4) for σ in is shown in Fig. 8 by solid squares. The BCS expression (A.4) for the intrinsic variance can be further simplified for a uniform single-particle spectrum (in which case the gap ∆ is level-independent) with a bandwidth E c ≫ ∆. In this limit most single-particle levels satisfy |ǫ i − λ| ≫ ∆ and the BCS amplitudes can be replaced by their non-interacting values. This leads to the simple expression in Eq. (26).
In Fig. 8 we also compare the BCS estimate (A.4) (solid squares) with its simplified version (26) (open squares). The latter slightly overestimates the BCS expression; both underestimate the exact intrinsic variance but provide a reasonable estimate within a factor of ∼ 2. All cases in Fig. 8 scale as Ω 3/2 (the respective fits are shown by lines).
In the remaining part of this Appendix we show that Eq. (A.1) (with λ being the chemical potential) is a good approximation, i.e., that the covariance contribution to σ 2 in in (A.1) is small. Eq. (25) holds for an arbitrary parameter ζ replacing λ. Since the number of particles N is conserved, we have in general and cov(n i , n j ) = n i n j −n inj is the covariance of n i and n j . The covariance contribution to (A.5) vanishes for values of ζ for which A(ζ) = 0, a quadratic equation in ζ. Using particle number conservation 2 j n i = N = 2 jn i , we have is the midpoint between the two solutions.
As an example we show in Fig. 9 the quantity A/σ 2 in as a function of ζ for a uniform single-particle spectrum.
In general the zeros ζ ± of A are not known without performing a full CSMC calculation. However, we observe that |A/σ 2 in | is rather small in the region between ζ − and ζ + as compared with its typical value outside this region. Thus taking ζ ≈ ζ 0 in Eq. (A.5) and ignoring A leads to a good approximation to σ 2 in . The sums on the r.h.s. of Eq. (A.10) are dominated by those levels i for whichn i is close to 1/2, i.e., by levels in the vicinity of the chemical potential λ. Thus we expect ζ 0 to be in proximity to the chemical potential. We can also estimate ζ 0 directly from Eq. (A.10) using the BCS expressions forn i in Eq (A.3).