text:
# Nuclear isotope shifts within the local energy-density functional approach

## 1 Introduction

For a long time, the fine structure of the isotopic dependence, i.e. dependence on the neutron number, of nuclear charge radii could not be explained satisfactorily. Because of the apparent lack of relevant experimental information, the effective particle–particle (pp) force which leads to the observed pairing properties of nuclei, has mostly been assumed in a very simple form, just sufficient to produce reasonable gap parameters $`\overline{\mathrm{\Delta }}`$ extracted from the observed odd-even staggering in the nucleon separation energies. Even more sophisticated effective interactions, e.g. the Gogny force , do not include a density dependence in the pp channel.

It then has been demonstrated in Hartree-Fock-Bogolyubov (HFB) type calculations that three- or four-body forces are essential to reproduce the experimental data on isotope shifts of nuclear charge radii , and indeed that these isotope shifts give indirect experimental information on the effective in-medium many-body force, or what is equivalent, on the density dependence of the effective interaction, particularly in the pp channel.

It should be emphasized that a better knowledge of the density dependence of the nuclear pairing force would give the possibility to predict the pairing gap as a function of nuclear matter density which is, in particular, of great importance for understanding the pairing phenomena in neutron stars. It is well known that at present the pairing gap can not be obtained with sufficient accuracy from nuclear matter calculations based on bare NN interaction. The empirical information gained from the studies of real nuclei, specifically from the combined analyses of the nucleon separation energies and isotope shifts in charge radii, seems to be indispensable in this respect. A more general remark is that, as one may notice, a “good” microscopic theory, which could supply an effective interaction for describing the nuclear ground states and low-energy nuclear structure on a satisfactory level, is still lacking. The presently most successful simultaneous description of the bulk nuclear properties, such as binding energies and radii, throughout the periodic chart is achieved with phenomenological density-dependent interaction such as the Skyrme force . This can be traced to the fact that some major effects produced by the density dependence of the effective interaction in the particle–hole (ph) channel, which is the second variational derivative of the corresponding energy-density functional with respect to normal density, are now established fairly well111One may still notice that a universal and unique parametrization of the ph force has not yet been found. In particular, considerable effort is still continuing to optimize the HF part of the Skyrme-type functional .. On the same ground, one expects that simultaneous description of the differential observables, such as odd-even mass differences and the odd-even effects in radii, would shed light on the density dependence of the effective interaction in the pp channel.

In  an effective interaction has been chosen which consisted of two- and three- (or four-) body parts, with parameters adjusted independently of those of the single-particle potential used for the shell model description of the reference nucleus. Only differences with respect of this reference nucleus have been calculated in a self-consistent way.

This procedure corresponds to the philosophy of the Landau-Migdal theory of finite Fermi systems , where there is no simple connection between the single particle well parameters and the effective interaction of quasiparticles close to the Fermi surface. However, there are consistency requirements  and since the parameters have been fitted to a large number of data, a quite reliable set of force parameters is available . The interaction of Refs.  has been restricted by the requirement that it should reduce to the Migdal force in the ph channel, the three-body part leading to the density dependence.

In the meantime, since Migdal had formulated his theory of finite Fermi systems (FFS), there has been much progress in Hartree-Fock calculations , as well as in the self-consistent version of the FFS theory , demonstrating that initially Migdal may have been too pessimistic concerning the possibility to use the same interaction for ground state properties and for low-lying excited states.

In deriving the HF or HFB equations, the starting point is the minimization of the expectation value of the energy, expressed by an effective interaction. The latter is considered as a substitute for the $`G`$-matrix derived from the true NN interaction, appearing in some kind of “effective Hamiltonian”. Therefore, ideally, the same interaction should be used in the pp (hh) and the particle-hole (ph) channel. The energy is then obtained from the interaction and the state vector which is assumed to be a Slater determinant.

The self-consistent FFS theory or energy density functional (EDF) method starts from relations between the total energy, expressed as a functional of the density, the single particle potential, the effective interaction, and the quasiparticle density: Instead of parameterizing the interaction and deriving the energy from it in a certain approximation, an ansatz for the energy functional is assumed from which the other quantities can be derived. The interaction in the pp (hh) channel is obtained from the energy functional by a different procedure than that in the ph channel, and therefore these two interactions are different, which is of great practical importance. This feature is shared with Migdal’s theory of finite Fermi systems.

Otherwise, from a technical point of view, there is little difference between the EDF method and the HF or, with pairing, HFB method. The latter correspond to special choices of the energy functional. On the other hand in the EDF method, to compute the density, a Slater determinant or HFB-type state vector is employed. Thus the differences which exist in principle, due to the necessary approximations are of little practical importance.

In a recent paper  the EDF method has been applied to magic nuclei 40,48Ca and 208Pb, investigating especially the influence of the spin–orbit interaction on collective states. Besides demonstrating the importance of including the two-body spin orbit force for the consistency of the method in describing the excited states, with the chosen set of parameters good results on the ground state properties have been obtained.

Here an analogous EDF approach is applied to open-shell nuclei with special attention to the pairing part of the functional. First results obtained mainly in diagonal-pairing approximation (i.e. with a HF–BCS formalism) have already been given in a few short publications  (the paper  gives the first application of the EDF approach to deformed nuclei and contains the results of a HFB-type calculation for dysprosium isotopes). In the present paper we give the results of more elaborate calculations based on the coordinate-space technique developed in  which allows us, for a local pairing field $`\mathrm{\Delta }(\stackrel{}{r})`$, to solve the Gor’kov equations in spherical finite systems without any approximations . The latter approach, even in the case of contact pairing force, corresponds, to a good approximation as will be shown in the present paper, to a full treatment of the HFB problem by applying the general variational principle to the EDF with a fixed energy cutoff $`ϵ_\mathrm{c}>ϵ_\mathrm{F}`$ ($`ϵ_\mathrm{F}`$ is the Fermi energy). It is known from nuclear matter calculations that $`s`$-wave pairing vanishes at or slightly above the equilibrium density $`0.16`$ fm-3; and so, in our local EDF approach, the cutoff is chosen to be larger than the corresponding Fermi energy $`ϵ_{0\mathrm{F}}37`$ MeV. In actual calculations we shall use $`ϵ_\mathrm{c}=40`$ MeV.

In Section 2 we give some theoretical foundation of the local EDF method for superfluid nuclei. In Section 3 the normal part of the energy-density functional with the parameters used in the calculations is described. In Section 4 the pairing part of the functional is presented. There, a possible extrapolation to uniform nuclear matter and the dilute limit with weak and strong pairing is also considered. Section 5 is devoted to the analysis of specific numerical results for a few isotopic chains of spherical nuclei and to the comparison of these results with experimental observations and other theoretical approaches. The influence of the density-dependent pairing on the evolution of the mean energies of multipole excitations is also studied using the self-consistent sum rule approach. In Section 6, the contribution of the ground state (“phonon”) correlations to the charge radii is analyzed. In Section 7 the conclusions are summarized. Some important theoretical and technical aspects are presented in the Appendices A, B, and C. Appendix A contains a thorough formulation of the generalized variational principle for the systems with pairing correlations which could be described by a local cutoff EDF. In Appendix B, the expression for the pairing energy is derived by using the Green’s function formalism. In Appendix C, we give a detailed description of the coordinate-space technique which we apply here to spherical nuclei.

## 2 Method of the local energy functional with pairing

In this section we give some foundation of the local-energy-density functional method which we apply here. Let us first briefly mention a few approaches used to calculate the nuclear ground state energy which is the expectation value of the Hamiltonian $`\widehat{H}`$ (with a “bare” NN-interaction) over the exact ground state vector $`|\mathrm{\Phi }`$:

$$E=\mathrm{\Phi }|\widehat{H}|\mathrm{\Phi }.$$
(2.1)
In general, the vector $`|\mathrm{\Phi }`$ is a functional of the particle density $`\varrho ^\mathrm{p}`$  given by

$$\varrho ^\mathrm{p}(\stackrel{}{r}_1)=\mathrm{Tr}_{st}\mathrm{\Phi }|\mathrm{\Psi }^{}(1)\mathrm{\Psi }(1)|\mathrm{\Phi },$$
(2.2)
where $`s`$ and $`t`$ are the spin and isospin indices and $`\mathrm{\Psi }^{}(1)`$ and $`\mathrm{\Psi }(1)`$ are the particle creation and annihilation operators, respectively, $`(1)=(\{\stackrel{}{r},s_z,\tau _3\}_1)`$. Therefore, the energy is a functional of $`\varrho ^\mathrm{p}`$:

$$E=E[\varrho ^\mathrm{p}].$$

It is a difficult problem to construct such a functional and calculate the energy for many body systems. Formally, one can introduce single-particle degrees of freedom and, assuming that there is no pairing condensate, extract the single-particle motion on the HF level by writing222Of course, in the nuclear case, this expression with bare $`\widehat{H}`$ is meaningless since $`\mathrm{HF}|\widehat{H}|\mathrm{HF}`$ becomes infinite due to the repulsive core. It is understood then that, for the first term in (2.3), the $`\widehat{H}`$ is taken with some kind of microscopically derived effective interaction, for example with the Brueckner $`G`$-matrix (see, e.g., the discussion in ). Because of inevitable approximations, such a replacement would never lead to the vanishing of the second term, $`E_{\mathrm{corr}}`$.

$$E[\widehat{\rho }]=\mathrm{HF}|\widehat{H}|\mathrm{HF}+E_{\mathrm{corr}}[\widehat{\rho }],$$
(2.3)
where $`|\mathrm{HF}`$ is the HF vacuum — a Slater determinant of single-particle wave functions $`\varphi _i`$. The latter are defined in a self-consistent manner by the equation $`\widehat{h}\varphi _i=ϵ_i\varphi _i`$ in which the single-particle Hamiltonian $`\widehat{h}`$ is given by the variational derivative of the energy functional with respect to the single-particle density matrix $`\widehat{\rho }`$: $`\widehat{h}=\delta E/\delta \widehat{\rho }`$. The second term, $`E_{\mathrm{corr}}`$, is the dynamical correlation energy (short range correlations, RPA correlations, many-particle correlations, parquet diagrams, etc.). The density matrix $`\widehat{\rho }`$ may be associated with the Landau-Migdal quasiparticle density matrix defined as

$$\widehat{\rho }(1,2)=\mathrm{HF}|\psi ^{}(2)\psi (1)|\mathrm{HF},$$

with $`\psi ^{}`$ ($`\psi `$) the quasiparticle creation (annihilation) operators. It is in turn a functional of the particle density $`\varrho ^\mathrm{p}`$:

$$\widehat{\rho }=\widehat{\rho }[\varrho ^\mathrm{p}].$$

In homogeneous infinite nuclear matter, because of equality of the particle and quasiparticle numbers, the following relation between the two densities holds:

$$\varrho ^\mathrm{p}=\rho =\mathrm{Tr}_{st}\widehat{\rho }(1,1).$$
(2.4)
In an inhomogeneous system, the local difference between the densities could be attributed to the quasiparticle form factor . Using the effective radius approximation one may write

$$\varrho ^\mathrm{p}(\stackrel{}{r})=(1+\frac{1}{6}R_\mathrm{q}^2\stackrel{}{}^2)\rho (\stackrel{}{r}),$$
(2.5)
where $`R_\mathrm{q}`$ is the effective quasiparticle radius.

Hence the mean square nuclear radius for the particles may differ from that for the quasiparticles:

$$r^2^\mathrm{p}=r^2+R_\mathrm{q}^2.$$
(2.6)
Some additional effects may be related with a possible change of the internal structure of nucleons in nuclear medium. Such effects are not included explicitly in the ordinary approaches, based either on a “bare” Hamiltonian $`H`$ or on some kind of effective interaction, in which the nucleons or quasiparticles are considered as point-like objects333To get the charge density distribution, and the charge radius of a nucleus, the point nucleon densities are usually folded with the free nucleon charge form factors. These may also change in the nuclear medium.. Within the EDF approach, however, tracing back to the Hohenberg–Kohn density-functional theory , if the functional were known, its minimization would determine both the ground-state energy and the correct particle density $`\varrho ^\mathrm{p}(\stackrel{}{r})`$. Anyway, for the differences between mean square radii of nuclei, which is of our main interest here, one has

$$\delta r^2^\mathrm{p}=\delta r^2.$$
(2.7)
The common problem of any many-body theory is how to calculate the dynamical correlation energy $`E_{\mathrm{corr}}`$ or how to take it into account in some effective way. In nuclear physics a few self-consistent methods are very popular:

* The HF method with effective density-dependent interaction (e.g., with zero-range Skyrme  or finite-range Gogny force ). In this method the energy is calculated as a ground state expectation value of an effective Hamiltonian $`H_{\mathrm{eff}}`$ taken over a Slater determinant,

$$E[\widehat{\rho }]=\mathrm{HF}|\widehat{H}_{\mathrm{eff}}|\mathrm{HF}.$$
(2.8)

* The relativistic mean field model  based on a phenomenological relativistic Lagrangian which includes mesonic and nucleonic degrees of freedom. In most applications to the ground states of finite nuclei, the relativistic Hartree formalism with static meson fields, and with no-sea approximation for nucleons, has been used (see  and references therein).

* The quasiparticle Lagrangian method (QLM) which is an extension  of the Landau-Migdal quasiparticle concept. In this method an effective local quasiparticle Lagrangian is constructed taking into account the first-order energy-dependence of the nucleon mass operator and the requirements imposed by self-consistency relations . The QLM can be reformulated in terms of an effective Hamiltonian  in which the (linear) energy dependence of the nucleon mass operator is completely hidden so that this approach becomes equivalent to the EDF method (see also the discussion of this point in ).

* The method of the effective quasiparticle local energy density functional (which is referred to as EDF here). The possibility to use this method is based on the existence theorem of Hohenberg and Kohn . With the Kohn-Sham quasiparticle formalism , this theorem allows one to write the nuclear EDF in the form

$$E[\widehat{\rho }]=\mathrm{Tr}(t\widehat{\rho })+E_{\mathrm{int}}[\rho ],$$
(2.9)
where the first term is taken with the free kinetic energy operator $`t=p^2/2m`$ ($`m`$ is the bare nucleon mass). The ground state energy is then determined by making this EDF stationary with respect to infinitesimal variations of the single-particle wave functions belonging to the class of the Slater determinants from which the density matrix $`\widehat{\rho }`$ (and the density $`\rho `$) can be calculated self-consistently. This method is flexible in the sense that the EDF does not (and should not) expose the same symmetry properties for the effective force, which is the second variational derivative of the EDF with respect to $`\rho `$, as the underlying bare NN interaction or effective Skyrme-type forces. Particularly, the force in the ph channel obtained from the EDF does not come out to be antisymmetrized, and the effective force in the pp channel needed for describing the nuclear pairing properties has to be derived in a different way. As a result, the matrix elements of the pp force, although they should be antisymmetrized, have no direct connection, say by means of simple angular momentum recoupling like the Pandya transformation, with those of the ph force444This was explicitly demonstrated, already on the level of the “ladder” approximation, in . We note the term “force”, used to name the second variational derivative of the EDF with respect to the density matrix, may be somewhat confusing in this context. In fact, this derivative, for infinite systems, has strict correspondence to the quasiparticle interaction amplitude on the Fermi surface introduced in the Landau theory of Fermi liquids . The total quasiparticle scattering amplitude is, of course, antisymmetric due to the Pauli principle resulting in sum rules for the Landau parameters. These parameters should satisfy the Pomeranchuk stability conditions  in order to prevent the collapse of the system with respect to the low-energy ph excitations and provide the ground state energy to be a minimum rather than simply stationary. For finite systems, this means that all the RPA solutions obtained with the self-consistent ph interaction, taken as the second variational derivative of EDF at the stationary point, should have real frequencies $`\omega `$ such that $`\omega ^2>0`$.. This is, as discussed in the Introduction, in accord with the philosophy of Migdal theory of finite Fermi systems . (For an example of the phenomenological nuclear EDF see ).

In most of the applications of the conventional functionals to finite nuclei the pairing correlations are introduced either at the BCS or the HFB level (see, e.g., Ref.  and references therein). When the contact pairing interaction is used, the pairing part of the functional is usually evaluated within a truncated space of the quasiparticle levels imposing an energy cutoff which is often chosen with some freedom but not too far from the Fermi surface. A more rigorous approach should be based on the general variational principle for the cutoff functional which uses the same quasiparticle basis both in the ph and pp channel.

We proceed as follows. In nuclei with pairing correlations one works with approximate state vectors which are not eigenstates of the particle number. One assumes nonzero anomalous expectation values, $`\widehat{\nu }(1,2)=\mathrm{\Phi }|\mathrm{\Psi }(1)\mathrm{\Psi }(2)|\mathrm{\Phi }0`$, and the energy of a superfluid nucleus may be given by a functional of the generalized density matrix $`\widehat{R}`$ containing both normal, $`\widehat{\rho }`$, and anomalous, $`\widehat{\nu }`$, components (see Appendix A). In analogy to eq. (2.3), the energy of a superfluid nucleus can be written as

$$E[\widehat{R}]=\mathrm{HFB}|\widehat{H}|\mathrm{HFB}+E_{\mathrm{corr}}[\widehat{R}]$$
(2.10)
with $`|\mathrm{HFB}`$ the HFB quasiparticle vacuum and $`E_{\mathrm{corr}}`$ the dynamical correlation energy.

The generalized density matrix $`\widehat{R}`$, in turn, is a functional of $`\varrho ^\mathrm{p}`$ and therefore

$$E=E[\widehat{R}[\varrho ^\mathrm{p}]]=E[\varrho ^\mathrm{p}].$$
(2.11)
The method we follow here is based on an effective quasiparticle energy functional of the generalized density matrix:

$$E[\widehat{R}]=E_{\mathrm{kin}}[\widehat{\rho }]+E_{\mathrm{int}}[\widehat{\rho },\widehat{\nu }],$$
(2.12)
where $`E_{\mathrm{kin}}[\widehat{\rho }]=\mathrm{Tr}(t\widehat{\rho }),`$ and $`E_{\mathrm{int}}[\widehat{\rho },\widehat{\nu }]=E_{\mathrm{int}(\mathrm{normal})}[\widehat{\rho }]+E_{\mathrm{anomal}}[\widehat{\rho },\widehat{\nu }].`$ The anomalous energy $`E_{\mathrm{anomal}}`$ is chosen such that it vanishes in the limit $`\nu 0`$.

It should be emphasized that the weak-pairing approximation $`|\overline{\mathrm{\Delta }}|ϵ_\mathrm{F}`$ is assumed throughout this paper. That means we need to retain only the first-order term $`\nu ^2`$ in the $`\nu `$-dependent, anomalous part of the energy functional. One may then write

$$E_{\mathrm{anomal}}[\widehat{\rho },\widehat{\nu }]=\frac{1}{4}\left(\widehat{\nu }^{}\widehat{}_\mathrm{a}^{\mathrm{pp}}[\widehat{\rho }]\widehat{\nu }\right),$$
(2.13)
where $`\widehat{}_\mathrm{a}^{\mathrm{pp}}`$ is an antisymmetrized effective interaction in the pp channel and the round brackets (…) imply integration and summation over all variables.

To calculate the ground state properties, one can now use the general variational principle with two constraints,

$$\mathrm{HFB}|\widehat{N}(\mu )|\mathrm{HFB}N(\mu )=N,$$
(2.14)
$$\widehat{R}^2=\widehat{R},$$
(2.15)
leading to the variational functional of the form

$$I[\widehat{R}]=E[\widehat{R}]\mu N(\mu )\mathrm{Tr}\widehat{\mathrm{\Lambda }}(\widehat{R}\widehat{R}^2),$$
(2.16)
where $`N`$ is the particle number, $`\mu `$ the chemical potential, and $`\widehat{\mathrm{\Lambda }}`$ the matrix of Lagrange parameters (see, e.g. Ref. , and Appendix A for more details).

The major difficulty in implementing this general variational principle in practice is connected with the anomalous part of the energy functional. The anomalous energy (2.13) can be calculated if one knows a solution of the gap equation for the pairing field $`\widehat{\mathrm{\Delta }}`$ and the anomalous density matrix $`\widehat{\nu }`$. In general, the gap equation,

$$\widehat{\mathrm{\Delta }}=\frac{1}{2}\widehat{}_\mathrm{a}^{\mathrm{pp}}\widehat{\nu },$$
(2.17)
is nonlocal and its solution (starting, for example, from a realistic bare NN interaction ), even for uniform nuclear matter, poses serious problems. First attempts were made very recently to construct an effective pairing interaction for semi-infinite nuclear matter  and for finite nuclei  with an approximate version of the Brueckner approach. Those studies are in an initial stage, and so far they do not give any guidance how to choose an effective pairing interaction, in particular its density dependence, that could be used in nuclear structure calculations. It is assumed that a simple universal effective interaction in the pp-channel can be invented to correctly describe the nuclear pairing properties. Our intention is to show that, to a good approximation in the case of weak pairing $`|\overline{\mathrm{\Delta }}|ϵ_\mathrm{F}`$, the EDF method and the general variational principle can be used with an effective density-dependent contact pp-interaction.

We proceed with the formal development by using the Green’s function formalism as described in detail in Appendices A and B. Here we give only a brief account of the main issues.

We introduce an arbitrary cutoff $`ϵ_\mathrm{c}`$ in the energy space, but such that $`ϵ_\mathrm{c}>ϵ_\mathrm{F}`$, and split the generalized density matrix into two parts,

$$\widehat{R}=\widehat{R}_\mathrm{c}+\delta _\mathrm{c}\widehat{R},$$
(2.18)
where $`\delta _\mathrm{c}\widehat{R}`$ is related to the integration over energies $`|ϵ|>ϵ_\mathrm{c}`$.

The gap equation is renormalized to yield

$$\widehat{\mathrm{\Delta }}=\frac{1}{2}\widehat{}_\mathrm{a}^\xi \widehat{\nu }_\mathrm{c},$$
(2.19)
where $`\widehat{\nu }_\mathrm{c}`$ is the cutoff anomalous density matrix and $`\widehat{}_\mathrm{a}^\xi `$ is the effective antisymmetrized pp-interaction in which the contribution coming from the energy region $`|ϵ|>ϵ_\mathrm{c}`$ far from the Fermi surface is included by renormalization (see Appendix A).

For homogeneous infinite matter with weak pairing ($`|\overline{\mathrm{\Delta }}|ϵ_\mathrm{F}`$), it is shown that the variational functional $`E\mu N`$ does not change in first order in $`|\overline{\mathrm{\Delta }}|^2/ϵ_\mathrm{F}`$ upon variation with respect to $`\delta _\mathrm{c}\widehat{R}`$. The total energy of the system and the chemical potential also remain the same in first order in $`|\overline{\mathrm{\Delta }}|^2/ϵ_\mathrm{F}`$ if one imposes the particle number constraint (2.14) for the cutoff functional. To a good approximation, as discussed in Appendix A, this should also be valid for finite (heavy) nuclei.

Such an outcome may be understood by noting that the major pairing effects, if $`|\overline{\mathrm{\Delta }}|ϵ_\mathrm{F}`$ which is generally the case in nuclear matter and in finite nuclei, are developed near the Fermi surface. In this connection we mention the well-known fact that the pairing energy, in the BCS approximation, is defined by a sum concentrated near the Fermi surface (see, for example, Ref. ):

$$E_{\mathrm{pair}}\frac{1}{2}\overline{n}\overline{\mathrm{\Delta }}^2,$$
(2.20)
with $`\overline{n}`$ the average level density and $`\overline{\mathrm{\Delta }}`$ the average energy gap in the vicinity of the Fermi surface. In infinite matter, this corresponds to the pairing energy per particle (see Appendix B):

$$\frac{E_{\mathrm{pair}}}{N}=\frac{3}{8}\frac{\mathrm{\Delta }^2(\stackrel{}{p}_\mathrm{F})}{ϵ_\mathrm{F}}.$$
(2.21)
It follows that, with the cutoff functional, this leading pairing contribution to the energy of the system is exactly accounted for.

Thus we find that nuclear ground state properties can be described by applying the general variational principle to minimize the cutoff functional which has exactly the same form as in eq. (2.16) with the constraint (2.14) but with $`\widehat{R}`$ replaced by $`\widehat{R}_\mathrm{c}`$:

$$I_\mathrm{c}[\widehat{R}]=E_{\mathrm{kin}}[\widehat{\rho }_\mathrm{c}]+E_{\mathrm{int}}^\mathrm{c}[\rho _\mathrm{c},\nu _\mathrm{c}]\mu _\mathrm{c}N_\mathrm{c}(\mu _\mathrm{c})\widehat{\mathrm{\Lambda }}(\widehat{R}\widehat{R}^2),$$
(2.22)
$$\mathrm{HFB}|\widehat{N}_\mathrm{c}(\mu _\mathrm{c})|\mathrm{HFB}N_\mathrm{c}(\mu _\mathrm{c})=N.$$
(2.23)
Here $`E_{\mathrm{kin}}[\widehat{\rho }_\mathrm{c}]=\mathrm{Tr}(t\widehat{\rho }_\mathrm{c})`$ and $`E_{\mathrm{int}}^\mathrm{c}[\rho _\mathrm{c},\nu _\mathrm{c}]=E_{\mathrm{int}(\mathrm{normal})}^\mathrm{c}[\widehat{\rho }_\mathrm{c}]+E_{\mathrm{anomal}}^\mathrm{c}[\widehat{\rho }_\mathrm{c},\widehat{\nu }_\mathrm{c}]`$ with

$$E_{\mathrm{anomal}}^\mathrm{c}[\widehat{\rho }_\mathrm{c},\widehat{\nu }_\mathrm{c}]=\frac{1}{4}(\widehat{\nu }_\mathrm{c}^{}\widehat{}_\mathrm{a}^\xi [\widehat{\rho }_\mathrm{c}]\widehat{\nu }_\mathrm{c}).$$
(2.24)
The above consideration has been carried out without any pre-assumptions concerning the density functional $`E_{\mathrm{int}}[\widehat{R}]E_{\mathrm{int}}^\mathrm{c}[\rho _\mathrm{c},\nu _\mathrm{c}]`$. Now, recalling the Hohenberg-Kohn theorem , we specify that the density functional can be chosen to be of a local form, i.e. dependent on the normal and anomalous local real densities $`\rho (\stackrel{}{r},\tau )`$ and $`\nu (\stackrel{}{r},\tau )`$ defined in Appendix A by (A.61) and (A.64), respectively. Then $`E_{\mathrm{normal}}^\mathrm{c}`$ is a functional of the normal densities (isoscalar, isovector, spin–orbital, etc.), and $`E_{\mathrm{anomal}}^\mathrm{c}`$ is a functional of the normal and anomalous densities, the latter acquiring the simple form

$$E_{\mathrm{anomal}}^\mathrm{c}[\rho _\mathrm{c},\nu _\mathrm{c}]=\underset{\tau =\mathrm{n},\mathrm{p}}{}\stackrel{}{r}\nu _\mathrm{c}^{}(\stackrel{}{r},\tau )^\xi (\stackrel{}{r},\tau ;[\rho _\mathrm{c}])\nu _\mathrm{c}(\stackrel{}{r},\tau ).$$
(2.25)
The pairing potential is defined by

$$\mathrm{\Delta }(\stackrel{}{r},\tau )=\frac{1}{2}\frac{\delta E_{\mathrm{anomal}}^\mathrm{c}[\rho _\mathrm{c},\nu _\mathrm{c}]}{\delta \nu _\mathrm{c}(\stackrel{}{r},\tau )}.$$
(2.26)
That corresponds to the gap equation which takes on now a very simple multiplicative form:

$$\mathrm{\Delta }(\stackrel{}{r},\tau )=^\xi (\stackrel{}{r},\tau ;[\rho _\mathrm{c}])\nu _\mathrm{c}(\stackrel{}{r},\tau ).$$
(2.27)
The normal mean-field potential is given by

$$U(\stackrel{}{r},\tau )=\frac{\delta E_{\mathrm{int}}^\mathrm{c}[\rho _\mathrm{c},\nu _\mathrm{c}]}{\delta \rho _\mathrm{c}(\stackrel{}{r},\tau )}.$$
(2.28)
Due to the density dependence of $`^\xi `$ it includes also a contribution arising from the variation of the pairing interaction energy:

$$U_{\mathrm{pair}}(\stackrel{}{r},\tau )=\frac{\delta E_{\mathrm{anomal}}^\mathrm{c}[\rho _\mathrm{c},\nu _\mathrm{c}]}{\delta \rho _\mathrm{c}(\stackrel{}{r},\tau )}.$$
(2.29)
We emphasize that having found the gap function $`\mathrm{\Delta }(\stackrel{}{r})`$, and the mean-field potential $`U(\stackrel{}{r})`$, the Gor’kov equations  may be solved, for spherical nuclei, exactly by using the coordinate-space technique (see Refs.  and Appendix C). The generalized Green’s function obtained this way can be integrated over energy up to $`ϵ_\mathrm{c}`$ to yield both the normal and anomalous densities $`\rho _\mathrm{c}`$ and $`\nu _c`$ which are used then to compute the energy of the system. This corresponds to a full HFB treatment of the nuclear ground state properties. It remains, of course, an art to find a “good” local functional, in particular its anomalous part.

## 3 The normal part of the functional

In this and the next section we describe in some detail the energy-density functional which we use in the present calculations. We shall omit in the following the cutoff index c for simplicity. The functional is local in the sense that it contains only normal and anomalous densities, not the density matrix. The interaction part of the energy density is

$$\epsilon _{\mathrm{int}}=\epsilon _{\mathrm{main}}+\epsilon _{\mathrm{Coul}}+\epsilon _{sl}+\epsilon _{ss}+\epsilon _{\mathrm{anomal}}.$$
(3.1)
The terms $`\epsilon _{\mathrm{main}}`$ and $`\epsilon _{\mathrm{Coul}}`$ have been described in Refs. , where also the connection of the parameters involved with the characteristics of nuclear matter has been given. Two different versions of the spin–orbit part $`\epsilon _{sl}`$ have been used: the one of Ref. , and also, after , the variant555As a rule, the convention $`\mathrm{}=1`$ is used throughout the paper but $`\mathrm{}`$ will appear in some expressions. Hopefully this should not lead to any confusion.

$`\epsilon _{sl}(\stackrel{}{r})`$ $`=`$ $`{\displaystyle \frac{1}{2}}C_0r_0^2{\displaystyle \underset{i,k=\mathrm{p},\mathrm{n}}{}}\{\kappa ^{ik}[\rho ^i{\displaystyle \underset{ss^{}}{}}\psi ^{}(\stackrel{}{r},s)[\stackrel{}{p}\times \stackrel{}{\sigma }]\psi (\stackrel{}{r},s^{})^k`$ (3.2)
$`+`$ $`\rho ^k{\displaystyle \underset{ss^{}}{}}\psi ^{}(\stackrel{}{r},s)[\stackrel{}{p}\times \stackrel{}{\sigma }]\psi (\stackrel{}{r},s^{})^i]`$
$`+`$ $`g_1^{ik}{\displaystyle \underset{\alpha ,\beta }{}}{\displaystyle \underset{ss^{}}{}}\psi ^{}(\stackrel{}{r},s)\sigma _\alpha p_\beta \psi (\stackrel{}{r},s^{})^i{\displaystyle \underset{ss^{}}{}}\psi ^{}(\stackrel{}{r},s)\sigma _\alpha p_\beta \psi (\stackrel{}{r},s^{})^k\}`$
where $`\psi ^{}`$ and $`\psi `$ are field operators, the brackets denote the ground state expectation value, and a superscript $`k`$ indicates that a projection on particles of type $`k`$ is performed.

As given by (3.2), $`\epsilon _{sl}`$ corresponds to the two-body spin–orbit interaction

$$_{sl}^\omega =C_0r_0^2(\kappa +\kappa ^{}\stackrel{}{\tau }_1\stackrel{}{\tau }_2)[_1\delta (\stackrel{}{r}_1\stackrel{}{r}_2)\times (\stackrel{}{p}_1\stackrel{}{p}_2)](\stackrel{}{\sigma }_1+\stackrel{}{\sigma }_2),$$
(3.3)
and to the first–order velocity harmonic of the spin–dependent part of the Landau interaction amplitude

$$_1^\omega =C_0r_0^2(g_1+g_1^{}\stackrel{}{\tau }_1\stackrel{}{\tau }_2)\delta (\stackrel{}{r}_1\stackrel{}{r}_2)(\stackrel{}{\sigma }_1\stackrel{}{\sigma }_2)(\stackrel{}{p}_1\stackrel{}{p}_2),$$
(3.4)
which has to be symmetrized in such a way that the momentum operators $`\stackrel{}{p}_{1,2}`$ never act on the delta-function $`\delta (\stackrel{}{r}_1\stackrel{}{r}_2)`$.

In these expressions, the multiplier $`r_0^2`$ has been introduced to make the parameters $`\kappa `$ and $`g_1`$ dimensionless. It is given by $`r_0^2=(3/8\pi \rho _0)^{2/3}`$ where $`\rho _0`$ is the equilibrium density of one kind of particles in symmetric nuclear matter. The factor $`C_0`$ is the inverse density of states at the Fermi energy, $`ϵ_{0\mathrm{F}}=k_{0\mathrm{F}}^2/2m^{}`$, in saturated nuclear matter. It is given by $`C_0=2ϵ_{0\mathrm{F}}/3\rho _0=\pi ^2/k_{0\mathrm{F}}m^{}`$.

Defining spin–orbit densities by

$`\rho _{sl}^k(\stackrel{}{r})={\displaystyle \underset{ss^{}}{}}\psi ^{}(\stackrel{}{r},s)(\stackrel{}{\sigma }\stackrel{}{l})\psi (\stackrel{}{r},s^{})^k,`$
then in spherically symmetric nuclei, the spin–orbit energy density (3.2) can be simplified to

$$\epsilon _{sl}=C_0r_0^2\underset{i,k=\mathrm{n},\mathrm{p}}{}\left(\frac{1}{r}\rho _{sl}^i\kappa ^{ik}\frac{\rho ^k}{r}+\frac{1}{4r^2}\rho _{sl}^ig_1^{ik}\rho _{sl}^k\right).$$
(3.5)
In deformed nuclei, the expression for $`\epsilon _{sl}`$ cannot be written in such a simple form. In the axially symmetric case the derivation can be performed in cylindrical coordinates; the resulting formulae are given in .

We just mention some features of those parts of $`\epsilon _{\mathrm{int}}`$ which are not given here in detail: The main contribution $`\epsilon _{\mathrm{main}}`$ consists of a volume part and of a surface term both containing simple fractional-linear functions of the normal isoscalar density $`\rho _+=\rho _\mathrm{n}+\rho _\mathrm{p}`$. This leads to the effective density-dependent forces in the scalar-isoscalar and scalar-isovector channels (see e.g. Refs.  for more details). There is no momentum dependence in these parts, therefore we have the effective mass $`m^{}`$ equal to the bare mass $`m`$. The surface term involves a finite, but small Yukawa range parameter $`R=0.35`$ fm; in the limit $`R0`$ it could be written with the help of the Laplace operator acting on functions of the density. It vanishes for constant density. The Coulomb part $`\epsilon _{\mathrm{Coul}}`$ is taken in the usual form  with an exchange contribution in the Slater approximation.

Note that in the more complicated parts of the force like those of eqs. (3.3) and (3.4) we did not introduce any density dependence. In general, all of the constants $`\kappa `$, $`g_1`$, and others to follow below could be density dependent as well. Until now there are no experimental data which would make it necessary to introduce such a dependence for the laboratory nuclei within or not too far from the stability valley. However, going to the nucleon drip lines and beyond, a more complicated density-dependent spin–orbit force might be of relevance. This is indicated by the relativistic mean field models for very neutron-rich nuclei  and by the exact Monte Carlo methods for small pure neutron drops . Both approaches predict a much smaller (by a factor of 2–3) spin–orbit potential than the usual Skyrme-type ansatz. The role played by the density dependence of the effective spin–orbit interaction is left to be studied in future work.

The spin–spin term $`\epsilon _{ss}`$ has hitherto been omitted since its expectation value in spin-zero nuclei vanishes. We include it here for completeness; it might become important for magnetic excitations or polarization. A possible simple form of it is

$$\epsilon _{ss}(\stackrel{}{r})=\frac{1}{2}C_0\underset{i,k=\mathrm{n},\mathrm{p}}{}g_0^{ik}\stackrel{}{\sigma }^i\stackrel{}{\sigma }^k.$$
(3.6)
Note that the strength $`g_0^{ik}`$ can be considered as the zeroth term in a decomposition of the spin–spin interaction amplitude at the Fermi surface in Legendre polynomials of the scattering angle. The next (first) term, proportional to $`g_1^{ik}`$, has been included in the spin–orbit part of the functional; however, the normalization chosen in eq. (3.2) is different from that of Ref. .

As long as there is no static spin polarization ($`\stackrel{}{\sigma }=0`$) in the ground states of even nuclei and our interest is only in electric and mass quantities we can continue to neglect $`\epsilon _{ss}`$. One delicate point should be addressed in this connection. Usually, when deriving the Hartree-Fock equations, one starts with a Hamiltonian containing a well behaved and properly antisymmetrized interaction. In the case of contact forces, antisymmetry prevents self-interaction of the particles, so there is no need for precautions against that. Going over to an effective force which, as in our case, is not antisymmetrized in the ph channel, the self-interaction is not automatically excluded. Dealing with a not-too-small number of particles, by fitting to the experimental data the self-interaction is compensated (or accounted for) in the average by the choice of the parameters.

This does not hold for the spin–spin part of the force. For an odd nucleus, consisting of a core with $`J=0`$ plus one particle, there is a nonvanishing spin density, and from eq. (3.6) there is also a contribution to the total energy which to lowest order is just the self-interaction of the odd particle and should therefore be excluded. On the other hand, considering the contribution of this part of the functional to the single particle potential when determining the wave function of the core, will give the polarization of the core which is a physical quantity. — Usually, on the level of Hartree-Fock type calculations the spin–spin interaction is left out; it may be considered in a second step by perturbation methods.

The basic set of functional parameters in the notations of Refs.  used in the present calculations is the set DF3:

$$\begin{array}{cc}a_+^\mathrm{v}=6.422,\hfill & a_{}^\mathrm{v}=5.417,\hfill \\ h_{1+}^\mathrm{v}=0.163,h_{2+}^\mathrm{v}=0.724,\hfill & h_1^\mathrm{v}=0,h_2^\mathrm{v}=3.0,\hfill \\ a_+^\mathrm{s}=11.1,h_+^\mathrm{s}=0.31,\hfill & a_{}^\mathrm{s}=4.00,h_{}^\mathrm{s}=0,\hfill \\ \kappa ^{\mathrm{pp}}=0.285,\kappa ^{\mathrm{pn}}=0.135,\hfill & g_1^{\mathrm{pp}}=g_1^{\mathrm{pn}}=0.12,\hfill \\ r_0=1.147\mathrm{fm},\hfill & R=0.35\mathrm{fm}.\hfill \end{array}\}$$
(3.7)
This set set has been specifically deduced in  to reproduce not only the already known ground-state properties of magic nuclei but also the very recent experimental data  on the single-particle energies near the “magic cross” at 132Sn. The set (3.7) corresponds to the following nuclear matter characteristics at saturation: compression modulus $`K_0=200`$ MeV, chemical potential (binding energy per nucleon) $`\mu _0=16.05`$ MeV, asymmetry energy parameter $`\beta _0=28.7`$ MeV, saturation density $`2\rho _0=0.1582`$ fm-3 ($`k_{0\mathrm{F}}=1.328`$ fm-1, $`ϵ_{0\mathrm{F}}=36.57`$ MeV, $`C_0=\pi ^2/k_{0\mathrm{F}}m=308.2`$ MeV$``$fm3). As in Ref. , the functional contains zero-range isoscalar spin–orbit interaction $`\kappa `$ and velocity-dependent spin-isospin interaction $`g_1^{}`$ (i.e. the first Landau-Migdal harmonic in the $`\sigma \tau `$ channel). This is in contrast to  where a finite range spin–orbit force had been used.

With the above parameter set of the normal part of the density functional the ground states of magic nuclei are described fairly well. For example, the calculated rms charge radii $`r_{\mathrm{ch}}^2^{1/2}`$ of 40Ca, 48Ca and 208Pb are 3.480 fm, 3.478 fm, and 5.500 fm, respectively, which agree nicely with experimental values deduced very recently  from the combined analysis of optical, muonic and elastic electron scattering data: $`r_{\mathrm{ch}}^2^{1/2}=3.4767(8)`$ fm for 40Ca, 3.4736(8) fm for 48Ca, and 5.5013(7) fm for 208Pb. The predictions for some other double magic nuclei are: $`r_{\mathrm{ch}}^2^{1/2}=3.721`$ fm (56Ni), 3.951 fm (78Ni), 4.453 fm (100Sn), and 4.705 fm (132Sn).

## 4 The pairing part of the functional

The pairing energy density $`\epsilon _{\mathrm{anomal}}`$ in eq. (3.1) is chosen in the form prescribed by (2.25):

$$\epsilon _{\mathrm{anomal}}(\stackrel{}{r})=\underset{\tau =n,p}{}^{\xi ,\tau \tau }(\stackrel{}{r};[\rho ])|\nu _\tau (\stackrel{}{r})|^2.$$
(4.1)
As discussed in Refs.  and mentioned in the introduction, effective three- or more-body forces are indispensable if one wants to reproduce isotope shifts in charge radii which means that $`\epsilon _{\mathrm{anomal}}`$ must depend on the normal density , and therefore we introduce

$$^{\xi ,\mathrm{nn}}=^{\xi ,\mathrm{pp}}=C_0f^\xi (x),$$
(4.2)
The dimensionless strength $`f^\xi `$ is supposed to be, within the EDF approach, a local functional of the isoscalar dimensionless density $`x=(\rho _\mathrm{n}+\rho _\mathrm{p})/2\rho _0`$, with $`\rho _{\mathrm{n}(\mathrm{p})}`$ the neutron (proton) density; $`C_0`$ and $`\rho _0`$ are defined in the previous section. We shall use the parametrization of $`f^\xi `$ suggested in Ref. :

$$f^\xi (x)=f^\xi (x(\stackrel{}{r}))=f_{\mathrm{ex}}^\xi +h^\xi x^q(\stackrel{}{r})+f_{}^\xi r_0^2\left(\stackrel{}{}x(\stackrel{}{r})\right)^2.$$
(4.3)
The parameter $`f_{\mathrm{ex}}^\xi `$ is negative and simulates an attraction in the pp channel in the far nuclear exterior, $`h^\xi `$ and $`f_{}^\xi `$ are taken to be positive . The exponent $`q`$ in the second term is introduced to have a more flexible parametrization. A repulsive short-range part of the $`G`$-matrix effective interaction with a $`x^{2/3}`$ dependence has been discussed, e.g., by Bethe many years ago  while developing a Thomas-Fermi theory for large finite nuclei. The choice $`q=\frac{2}{3}`$ seems thus to be reasonable, and indeed our calculations showed that some improvements in reproducing the isotopic shifts may be achieved with this choice compared to the linear case $`q=1`$. We shall use $`q=\frac{2}{3}`$ in the present paper. The three-body force of Refs.  corresponds to a linear dependence of $`^\xi `$ on $`\rho `$ ($`q=1`$, $`f_{}^\xi =0`$). As shown in Ref. , the self-consistent EDF (HF+BCS) calculations with the density-gradient term $`f_{}^\xi `$ in pairing force provide desirable size of isotopic shifts and right order of odd-even staggering observed in lead isotopes, the coupling of the proton mean field with neutron pairing being the major effect in this case. The variation of the anomalous energy (4.1) incorporating the above force with respect to normal densities gives the following contribution to the central mean-field potential:

$$U_{\mathrm{pair}}^{\mathrm{n},\mathrm{p}}=\frac{C_0}{4\rho _0}\left[qh^\xi x^{q1}2r_0^2f_{}^\xi \left(\mathrm{\Delta }x+\stackrel{}{}x\stackrel{}{}\right)\right]\left(|\nu _\mathrm{n}|^2+|\nu _\mathrm{p}|^2\right).$$
(4.4)
As will be shown in the next section, different choices of the parameters of the particle–particle force (4.3) are possible to reasonably describe the neutron separation energies and isotopic shifts in charge radii. In particular, the following sets are deduced for the lead isotopes:

$$\begin{array}{cccc}f_{\mathrm{ex}}^\xi =0.56,\hfill & h^\xi =0,\hfill & f_{}^\xi =0\hfill & \hfill (\mathrm{a})\\ f_{\mathrm{ex}}^\xi =1.20,\hfill & h^\xi =0.56,\hfill & f_{}^\xi =2.4\hfill & \hfill (\mathrm{b})\\ f_{\mathrm{ex}}^\xi =1.60,\hfill & h^\xi =1.10,\hfill & f_{}^\xi =2.0\hfill & \hfill (\mathrm{c})\\ f_{\mathrm{ex}}^\xi =1.79,\hfill & h^\xi =1.36,\hfill & f_{}^\xi =2.0\hfill & \hfill (\mathrm{d})\\ f_{\mathrm{ex}}^\xi =2.00,\hfill & h^\xi =1.62,\hfill & f_{}^\xi =2.0\hfill & \hfill (\mathrm{e})\\ f_{\mathrm{ex}}^\xi =2.40,\hfill & h^\xi =2.16,\hfill & f_{}^\xi =2.0\hfill & \hfill (\mathrm{f})\end{array}\}.$$
(4.5)
Let us, in the rest of this section, study the behavior of $`\mathrm{\Delta }`$ as a function of density $`\rho =2\rho _0x`$ (or the Fermi momentum $`k_\mathrm{F}=(3\pi ^2\rho /2)^{1/3}k_{0\mathrm{F}}x^{1/3}`$) in symmetric ($`\rho _\mathrm{n}=\rho _\mathrm{p}`$) uniform infinite nuclear matter with pairing interaction of eqs. (4.2), (4.3), and also discuss how the pairing affects the equation of state (EOS) and the position of the equilibrium point. The density-gradient term $`f_{}^\xi `$ vanishes in this case, thus only the terms with parameters $`f_{\mathrm{ex}}^\xi `$ and $`h^\xi `$ are left. For uniform matter the gap equation (2.27) reduces to

$$\mathrm{\Delta }(x)=^\xi (x)_{kk_\mathrm{c}}\frac{\stackrel{}{k}}{(2\pi )^3}\frac{\mathrm{\Delta }(x)}{2\sqrt{(ϵ_kϵ_\mathrm{F}(x))^2+\mathrm{\Delta }^2(x)}},$$
(4.6)
where666Note that, in our approach, the upper limit in the truncated momentum or energy space depends on density since $`ϵ_\mathrm{c}`$ is chosen to be measured from the Fermi level position determined by the chemical potential $`\mu `$ (see Appendix C). At any given density, $`\mu `$ is introduced as a Lagrange multiplier and therefore, in varying the cutoff functional, the phase space is kept fixed. $`k_\mathrm{c}=\sqrt{2m(ϵ_\mathrm{F}+ϵ_\mathrm{c})}/\mathrm{}`$ and $`ϵ_k=\mathrm{}^2k^2/2m`$. Canceling $`\mathrm{\Delta }(x)`$ on both sides of this equation, to find a nontrivial solution $`\mathrm{\Delta }0`$, one gets

$$\frac{1}{4}x^{1/3}f^\xi (x)_0^{t_\mathrm{c}(x)}t\frac{\sqrt{t}}{\sqrt{(t1)^2+\delta ^2(x)}}=1,$$
(4.7)
where $`\delta (x)=\mathrm{\Delta }(x)/ϵ_\mathrm{F}(x)`$ and $`t_\mathrm{c}(x)=1+ϵ_\mathrm{c}/ϵ_\mathrm{F}(x)`$ with $`ϵ_\mathrm{c}`$ the energy cutoff measured from the Fermi energy $`ϵ_\mathrm{F}=ϵ_{0\mathrm{F}}x^{2/3}`$. The solution of the last equation in the weak pairing approximation, for $`f^\xi (x)<0`$, is given by (see, e.g., ):

$$\mathrm{\Delta }(x)=8ϵ_{0\mathrm{F}}x^{2/3}\sqrt{\frac{s(x)1}{s(x)+1}}\mathrm{exp}\left(s(x)2+\frac{2}{f^\xi (x)x^{1/3}}\right),$$
(4.8)
with $`s(x)=\sqrt{1+ϵ_\mathrm{c}/ϵ_\mathrm{F}(x)}k_\mathrm{c}/k_\mathrm{F}`$. In the far nuclear exterior when $`x0`$, $`^\xi `$ should be determined by the free NN scattering. However, approaching the nuclear surface, finite range and nonlocal in-medium effects may already considerably modify the force even at low density, and therefore the best parametrization of an effective contact force at $`x1`$ might not necessarily agree with the vacuum value which could be extracted from the free NN interaction, as discussed by Migdal many years ago .

Shown in Fig. 1 are the values of dimensionless strength $`f^\xi `$ (upper panel) and the results for $`\mathrm{\Delta }`$ (lower panel) in infinite matter with the parameter sets (a)–(f) of eq. (4.5) deduced for the lead chain. The calculations are performed with cutoff $`ϵ_\mathrm{c}=40`$ MeV and Fermi energy $`ϵ_{0\mathrm{F}}=36.57`$ MeV of saturated nuclear matter from the parametrization (3.7) of the normal part of the density functional.

One should notice that $`ϵ_\mathrm{F}`$ entering the integrand of the gap equation (4.6) can be expressed directly through density by $`ϵ_\mathrm{F}=\mathrm{}^2k_\mathrm{F}^2/2m`$ with $`k_\mathrm{F}=(3\pi ^2\rho /2)^{1/3}k_{0\mathrm{F}}x^{1/3}`$ only if the pairing is weak indeed so that the dependence of the Fermi energy (and the chemical potential $`\mu `$) on $`\mathrm{\Delta }`$ can be disregarded. Otherwise one should introduce the particle number condition

$$x=\frac{2}{\rho _0}_{kk_\mathrm{c}}\frac{\stackrel{}{k}}{(2\pi )^3}n_k(x),$$
(4.9)
where

$$n_k(x)=\frac{1}{2}\left(1\frac{ϵ_kϵ_\mathrm{F}(x)}{\sqrt{(ϵ_kϵ_\mathrm{F}(x))^2+\mathrm{\Delta }^2(x)}}\right),$$
(4.10)
and solve the system of the two equations (4.6) and (4.9) with respect to $`\mathrm{\Delta }`$ and $`ϵ_\mathrm{F}`$. The results shown in Fig. 1 by full lines correspond to such a direct solution while those shown by dashed lines to the weak pairing approximation, eq. (4.8).

As seen in Fig. 1, the approximation (4.8) works well in the entire range of $`k_\mathrm{F}`$ for the set (a), but for the other sets this is true only at $`k_\mathrm{F}`$ greater than $``$1.2 fm-1 and also, for the sets (b), (c) and (d), at $`k_\mathrm{F}`$ less than 0.42, 0.14 and 0.042 fm-1, respectively (in these regions the ratio $`\mathrm{\Delta }/ϵ_\mathrm{F}`$ does not exceed 0.1). For the sets (e) and (f), at lower densities, eq. (4.8) can not be used any more to estimate the gap even by order of magnitude since, as seen from the behavior of the corresponding dashed lines, it becomes divergent.

All parameter sets (4.5) except (a) reproduce the neutron separation energies and the isotope shifts of charge radii of lead isotopes reasonably well (see next section). For the sake of comparison also shown in Fig. 1 are the values of the $`{}_{}{}^{1}S_{0}^{}`$ pairing gap in nuclear matter obtained for the CD-Bonn potential without medium effects (using free single-particle spectrum $`ϵ_k=k^2/2m`$ and for the Gogny D1 force in the HFB framework . The agreement between the two latter calculations is relatively good while both deviate noticeably from our predictions. The curve for contact density-independent pairing force, set (a), stands by itself with a positive derivative $`\mathrm{\Delta }(x)/x`$ everywhere; no acceptable description of $`r_{\mathrm{ch}}^2`$ could be obtained in this case (see Fig. 5). Interestingly, for the sets (b)–(e) which reproduce satisfactorily both the neutron separation energies $`S_\mathrm{n}`$ and mean squared charge radii $`r_{\mathrm{ch}}^2`$ there exists a pivoting point at $`k_\mathrm{F}1.15`$ fm-1 (at $``$0.65 of the equilibrium density) with the same value of $`\mathrm{\Delta }_{\mathrm{piv}}3.3`$ MeV ($`f_{\mathrm{piv}}^\xi 0.78`$).

One can see in Fig. 1 that for parameter sets with bigger absolute values of $`f^\xi `$ and $`h^\xi `$ the region where $`\mathrm{\Delta }(x)`$ varies strongly moves towards lower densities, the slope becomes steeper and around $`x1`$ the pairing gap tends to be very small indeed. One may think that in finite nuclei the effect of the strong dependence of $`\mathrm{\Delta }`$ on $`x`$, especially at small densities, would not influence the nuclear properties noticeably. But even if $`\mathrm{\Delta }`$ would be concentrated only in the surface region and outside of the nucleus, the anomalous density, due to quantum effects, would not vanish inside the nucleus where $`x1`$. In this case the effect of density dependence should persist in the nuclear volume and should be more pronounced at larger values of $`h^\xi `$. This point was confirmed by our calculations (cf. Ref.  and below).

Consider in more detail the behavior of $`\mathrm{\Delta }`$ at very low densities when $`k_\mathrm{F}0`$ assuming the weak pairing regime, i.e. $`f_{\mathrm{ex}}^\xi >f_{\mathrm{cr}}^\xi `$ where $`f_{\mathrm{cr}}^\xi `$ is a critical strength constant given by

$$f_{\mathrm{cr}}^\xi =2\frac{k_{0\mathrm{F}}}{k_{0\mathrm{c}}},$$
(4.11)
with $`k_{0\mathrm{c}}=\sqrt{2mϵ_\mathrm{c}}/\mathrm{}`$. This constant is determined by the condition

$$1+\frac{f_{\mathrm{cr}}^\xi }{4\sqrt{ϵ_{0\mathrm{F}}}}_0^{ϵ_\mathrm{c}}\frac{ϵ}{\sqrt{ϵ}}=0,$$
(4.12)
obtained from the gap equation (4.6) at $`\mathrm{\Delta }=ϵ_\mathrm{F}=0`$. For the chosen energy cutoff and parametrization (3.7) of our density functional we have $`f_{\mathrm{cr}}^\xi =1.912`$. To leading order, on the other hand, at $`k_\mathrm{F}0`$ from eq. (4.8) we obtain

$$\mathrm{\Delta }=cϵ_\mathrm{F}\mathrm{exp}\left(\frac{\pi }{2k_\mathrm{F}a}\right),a<0,$$
(4.13)
where $`c=8\mathrm{e}^21.083`$ and where we have introduced the quantity

$$a=\frac{\pi }{2k_{0\mathrm{F}}}\left(\frac{\sqrt{2mϵ_\mathrm{c}}}{\mathrm{}k_{0\mathrm{F}}}+\frac{2}{f_{\mathrm{ex}}^\xi }\right)^1\frac{\pi }{4k_{0\mathrm{F}}}\left(\frac{1}{f_{\mathrm{ex}}^\xi }\frac{1}{f_{\mathrm{cr}}^\xi }\right)^1.$$
(4.14)
It can be shown that this quantity is nothing but the singlet nucleon–nucleon scattering length. To be more specific, consider first the two-neutron problem in the vicinity of the critical point, $`f_{\mathrm{ex}}^\xi f_{\mathrm{cr}}^\xi `$, when the attraction is strong enough to produce a bound pair state (the scattering problem at any $`f^\xi `$ will be considered below). The bound-state wave function for the relative motion of two neutrons in case of contact interaction $`C_0f_{\mathrm{ex}}^\xi \delta (\stackrel{}{r})`$ is determined by (cf. Ref. ):

$$\psi (\stackrel{}{r})=C_0f_{\mathrm{ex}}^\xi G_0(\stackrel{}{r},ϵ_\mathrm{b})\psi (0),$$
(4.15)
where $`G_0`$ is the free Green’s function in the truncated space,

$$G_0(\stackrel{}{r},ϵ_\mathrm{b})=_{ϵ_k^{\mathrm{nn}}ϵ_\mathrm{c}^{\mathrm{nn}}}\frac{\stackrel{}{k}}{(2\pi )^3}\frac{\mathrm{e}^{\mathrm{i}\stackrel{}{k}\stackrel{}{r}}}{ϵ_k^{\mathrm{nn}}ϵ_\mathrm{b}},$$
(4.16)
with $`ϵ_\mathrm{b}`$ the binding energy and $`ϵ_k^{\mathrm{nn}}=\mathrm{}^2k^2/m`$ (the reduced mass is $`m/2`$). In the rest system of a nucleus, the center-of-mass energy $`ϵ^{\mathrm{nn}}`$ in the scattering problem would correspond to the energy $`ϵ^{\mathrm{nn}}/2`$ of each of the two nucleons in the $`s`$-wave pairing problem. This implies that the cutoff in the $`k`$-space in eq. (4.16), $`k_\mathrm{c}^{\mathrm{nn}}=\sqrt{mϵ_\mathrm{c}^{\mathrm{nn}}}/\mathrm{}`$, should be the same as in the gap equation. Thus, in the energy space, one has $`ϵ_\mathrm{c}^{\mathrm{nn}}=2ϵ_\mathrm{c}`$ and from (4.15) at $`\stackrel{}{r}=0`$ one obtains the equation to determine $`ϵ_\mathrm{b}`$:

$$1=\frac{1}{4\sqrt{ϵ_{0\mathrm{F}}}}f_{\mathrm{ex}}^\xi _0^{ϵ_\mathrm{c}}ϵ\frac{\sqrt{ϵ}}{ϵϵ_\mathrm{b}/2},$$
(4.17)
which reduces to

$$\frac{f_{\mathrm{ex}}^\xi }{f_{\mathrm{cr}}^{\xi \mathrm{nn}}}=\left(1\sqrt{\frac{ϵ_\mathrm{b}}{2ϵ_\mathrm{c}}}\mathrm{arctan}\sqrt{\frac{2ϵ_\mathrm{c}}{ϵ_\mathrm{b}}}\right)^1,$$
(4.18)
where $`f_{\mathrm{cr}}^{\xi \mathrm{nn}}`$ is a critical value of $`f_{\mathrm{ex}}^\xi `$ at which eq. (4.15) has a bound state solution $`ϵ_\mathrm{b}=0`$. Comparing eq. (4.17) at $`ϵ_\mathrm{b}=0`$ with eq. (4.12) one finds that $`f_{\mathrm{cr}}^{\xi \mathrm{nn}}`$ coincides with the value defined above by eq. (4.11). Now, in the vicinity of $`f_{\mathrm{cr}}^\xi `$ we can set $`\mathrm{arctan}\sqrt{2ϵ_\mathrm{c}/|ϵ_\mathrm{b}|}\pi /2`$ and use $`ϵ_\mathrm{b}=\mathrm{}^2/ma_{\mathrm{nn}}^2`$ in eq. (4.18). Then it is easy to see that the solution for the scattering length $`a_{\mathrm{nn}}`$ is exactly the same as given by (4.14). The expression (4.13) for $`\mathrm{\Delta }`$ agrees with the results of Ref.  based on a general analysis of the gap equation at low densities when $`k_\mathrm{F}|a|1`$. But we should stress that (4.13) is valid only in the weak-coupling regime corresponding to negative $`a`$. In the opposite case the gap in the dilute limit has to be found in a different way.

At $`f_{\mathrm{ex}}^\xi >f_{\mathrm{cr}}^\xi `$ the scattering length is negative, and from eq. (4.13) it follows that at low densities the pairing gap is exponentially small and eventually $`\mathrm{\Delta }(k_\mathrm{F}0)=0`$. Such a weak pairing regime with Cooper pairs forming in a spin singlet $`l=0`$ state exists up to the critical point at which the attraction becomes strong enough to change the sign of the scattering length. Then the strong pairing regime sets in, eqs. (4.8) and (4.13) are not valid any more, and $`\mathrm{\Delta }`$ should be determined directly from the combined solution of the gap equation (4.6) and the particle number condition (4.9). In the dilute systems, $`ϵ_\mathrm{F}`$ plays the role of the chemical potential $`\mu `$. The latter is defined by $`\mu =ϵ_\mathrm{F}(k_\mathrm{F})+U(k_\mathrm{F})`$ with $`U(k_\mathrm{F})`$ the HF mean field at the Fermi surface which is negligible for the fermion gas. At the critical point $`\mu `$ becomes negative and a bound state of a single pair of nucleons with the binding energy $`ϵ_\mathrm{b}=2\mu `$ becomes possible . This can be easily seen from the gap equation (4.6) written in the form

$$\left(\frac{k^2}{m}2\mu \right)\varphi _k=\mathrm{sgn}(ϵ_k\mu )\sqrt{1\varphi _k^2}_{k^{}k_\mathrm{c}}\frac{\stackrel{}{k}^{}}{(2\pi )^3}^\xi \varphi _k^{},$$
(4.19)
where we have introduced the functions $`\varphi _k=\mathrm{\Delta }/\sqrt{(ϵ_k\mu )^2+\mathrm{\Delta }^2}`$ and replaced $`ϵ_\mathrm{F}`$ by $`\mu `$. In the strong coupling regime, $`\mu <0`$, and in the dilute limit, $`|\varphi _k|1`$, this equation reduces to the Schrödinger equation for a single bound pair where $`2\mu `$ plays the role of the eigenvalue. It is equivalent to the coordinate-space equation (4.15). Thus $`2\mu =ϵ_\mathrm{b}`$ where the binding energy $`ϵ_\mathrm{b}`$ is determined from (4.18). The pairing gap $`\mathrm{\Delta }`$ at low densities in this regime can be found from the particle number condition (4.9) which in the leading order now reads

$$x=\frac{3\mathrm{\Delta }^2}{8ϵ_{0\mathrm{F}}^{3/2}}_0^{ϵ_\mathrm{c}}\frac{\sqrt{ϵ}ϵ}{(ϵ+|\mu |)^2}.$$
(4.20)
In the vicinity of the critical point from this equation we find

$$\mathrm{\Delta }^2=\frac{16}{3\pi }xϵ_{0\mathrm{F}}^2\frac{1}{k_{0\mathrm{F}}a},$$

which gives

$$\mathrm{\Delta }=\frac{\mathrm{}^2}{m}\left(\frac{2\pi \rho }{a}\right)^{1/2},a>0.$$
(4.21)
From this consideration it follows that, in the dilute case, the energy needed to break a condensed pair goes smoothly from $`2\mathrm{\Delta }`$ to $`2\mu =ϵ_\mathrm{b}`$ as a function of the coupling strength when the regime changes from weak to strong pairing. But as seen from (4.13) and (4.21), the behavior of $`\mathrm{\Delta }`$ at low densities is such that the derivative $`\mathrm{\Delta }/\rho `$ at $`\rho 0`$ as a function of $`f_{\mathrm{ex}}^\xi `$ exhibits a discontinuity from 0 to $`\mathrm{}`$. This is illustrated in Fig. 2 where we have plotted the gap $`\mathrm{\Delta }(\rho )`$ at very low densities for the sets (c)–(e) of eq. (4.5) embracing both regimes. We note also that the analytical expressions (4.13) and (4.21) give a purely imaginary gap at the critical point when the scattering length changes sign ($`k_\mathrm{F}a\sqrt{2m\mu }/\mathrm{}\mathrm{i}`$ if $`\mu `$ becomes negative). In this connection it is instructive to write the weak coupling expression (4.8) for $`\mathrm{\Delta }`$ in the following form :

$$\mathrm{\Delta }(k_\mathrm{F})=cϵ_\mathrm{F}\mathrm{exp}\left[\frac{\pi }{2}\mathrm{cot}\delta (k_\mathrm{F})\right],$$
(4.22)
where $`c=8\mathrm{e}^2`$ and where we have introduced the Fermi level phase shift $`\delta (k_\mathrm{F})`$ defined by

$$k_\mathrm{F}\mathrm{cot}\delta (k_\mathrm{F})=\frac{4k_{0\mathrm{F}}}{\pi }\left(\frac{1}{f^\xi (k_\mathrm{F})}+\frac{k_\mathrm{c}(k_\mathrm{F})}{2k_{0\mathrm{F}}}\right)\frac{k_\mathrm{F}}{\pi }\mathrm{ln}\left(\frac{k_\mathrm{c}(k_\mathrm{F})k_\mathrm{F}}{k_\mathrm{c}(k_\mathrm{F})+k_\mathrm{F}}\right),$$
(4.23)
with $`k_\mathrm{c}(k_\mathrm{F})=\sqrt{k_{0\mathrm{c}}^2+k_\mathrm{F}^2}`$. Eq. (4.23) corresponds to an exact solution of the nn scattering problem at the relative momentum $`k=k_\mathrm{F}`$ with the states truncated by a momentum cutoff $`k_\mathrm{c}=k_\mathrm{c}(k_\mathrm{F})`$ for contact interaction $`C_0f^\xi (k_\mathrm{F})\delta (\stackrel{}{r})`$ (see, e.g., Ref. ). In the dilute limit, from eq. (4.22) one notices again that the pairing gap becomes pure imaginary at the critical point since $`\mathrm{cot}\delta \mathrm{i}`$ when $`|a_{\mathrm{nn}}|\mathrm{}`$. At very low densities, with the parametrization (4.3) of the pairing force and with the chosen density-dependent cutoff, eq. (4.23) reduces to

$$k_\mathrm{F}\mathrm{cot}\delta (k_\mathrm{F})\frac{1}{a_{\mathrm{nn}}}+\frac{1}{2}r_{\mathrm{nn}}k_\mathrm{F}^2\frac{2k_\mathrm{F}}{\pi }\left[\frac{k_\mathrm{F}}{2k_{0\mathrm{c}}}\frac{2h^\xi }{(f_{\mathrm{ex}}^\xi )^2}\left(\frac{k_\mathrm{F}}{k_{0\mathrm{F}}}\right)^{3q1}\right],$$
(4.24)
where $`a_{\mathrm{nn}}`$ is the scattering length defined by eq. (4.14) and $`r_{\mathrm{nn}}`$ is the effective range, $`r_{\mathrm{nn}}=4/\pi k_{0\mathrm{c}}`$. The first two terms in this equation would describe low-energy behavior of the nn $`s`$-wave phase shift through an expansion of $`k\mathrm{cot}\delta `$ in powers of the relative momentum $`k=k_\mathrm{F}`$ if the interaction were density-independent – in our case, if the coupling strength and momentum cutoff were fixed by $`f^\xi =f_{\mathrm{ex}}^\xi `$ and $`k_\mathrm{c}=k_{0\mathrm{c}}`$, respectively. It follows that with a density-dependent effective force, such an expansion contains additional terms which for the parametrization used here are of the same order as the effective range term777Our effective pairing interaction with the choice $`q=1/3`$ would lead in the dilute limit to the expression for $`\mathrm{\Delta }`$ of the form of eq. (4.13) but with a different prefactor $`c`$ depending on the value of $`h^\xi `$. If, furthermore, we define the latter parameter by $`h^\xi =(1+2\mathrm{ln}2)(f_{\mathrm{ex}}^\xi )^2/6`$ we get in the leading order $`\mathrm{\Delta }(k_\mathrm{F})=(2/\mathrm{e})^{7/3}ϵ_\mathrm{F}\mathrm{exp}(\pi /2k_\mathrm{F}a)`$, i.e. the result obtained in Ref.  for a non-ideal Fermi gas by studying the singularities of the interaction amplitude (vertex function $`\mathrm{\Gamma }`$) and taking into account the terms up to the second order in $`k_\mathrm{F}a`$.. This simply demonstrates that, for reproducing the pairing gap, the effective interaction even at very low densities need not necessarily coincide with the bare NN interaction.

When the Fermi momentum $`k_\mathrm{F}`$ approaches from below the upper critical point at which the pairing gap closes, $`k_\mathrm{F}^{\mathrm{cr}}=k_{0\mathrm{F}}(f_{\mathrm{ex}}^\xi /h^\xi )^{1/3q}`$ defined by $`f^\xi (k_\mathrm{F})=0`$, we get from eq. (4.22)

$$\mathrm{\Delta }(k_\mathrm{F})cϵ_\mathrm{F}\mathrm{exp}\left[\frac{2k_{0\mathrm{F}}}{3qf_{\mathrm{ex}}^\xi (k_\mathrm{F}k_\mathrm{F}^{\mathrm{cr}})}\right].$$
(4.25)
Thus, at higher densities, the pairing gap becomes exponentially small when $`k_\mathrm{F}`$ approaches $`k_\mathrm{F}^{\mathrm{cr}}`$, in agreement with general analysis of the gap solutions . In weak coupling, as we have already shown, $`\mathrm{\Delta }(k_\mathrm{F})`$ is also exponentially small at low densities, in the vicinity of $`k_\mathrm{F}=0`$. It is noteworthy that near both critical points the pairing potential $`\mathrm{\Delta }`$ can be found from the behavior of the phase shift by using eq. (4.22) with a smooth density-dependent prefactor $`c(k_\mathrm{F})`$ (the details will be given elsewhere ). In the present paper, as an illustration, we show in Fig. 1 by the dotted line the values of $`\mathrm{\Delta }(k_\mathrm{F})`$ obtained from eq. (4.22), with a prefactor $`c=8\mathrm{e}^2`$, by using “experimental” nn phase shifts, without electromagnetic effects888We thank Rupert Machleidt for providing us with these nn phase shifts.. It is seen that $`\mathrm{\Delta }`$ obtained this way at low densities closely follows the solution of the gap equation with the CD-Bonn potential . The nn phase shift passes zero at the relative momentum $`k1.71`$ fm-1, and the gap should vanish at the corresponding Fermi momentum. Unfortunately, the solutions for $`\mathrm{\Delta }`$ are given in Ref.  only in the region up to $`k_\mathrm{F}=1.4`$ fm-1 which is rather far from the upper critical point.

For symmetric nuclear matter, with the functional DF3 used in the present paper, the energy per particle is

$$\frac{E}{A}(x)=\frac{2}{\rho _0x}_{kk_\mathrm{c}}\frac{\stackrel{}{k}}{(2\pi )^3}\frac{\mathrm{}^2k^2}{2m}n_k(x)+\frac{1}{3}ϵ_{0\mathrm{F}}a_+^\mathrm{v}f_+^\mathrm{v}(x)x+\frac{3\mathrm{\Delta }^2(x)}{2f^\xi (x)xϵ_{0\mathrm{F}}},$$
(4.26)
where $`f_+^\mathrm{v}(x)=(1h_{1+}^\mathrm{v}x)/(1+h_{2+}^\mathrm{v}x)`$ with the parameters given by (3.7); here the “particle–hole” term $`f_+^\mathrm{v}`$ vanishes in the dilute limit linearly in density. The chemical potential is

$$\mu (x)=ϵ_\mathrm{F}(x)+\frac{1}{3}ϵ_{0\mathrm{F}}a_+^\mathrm{v}[f_+^\mathrm{v}(x)x^2+2f_+^\mathrm{v}(x)x]+\frac{3f^\xi (x)}{2f^{\xi 2}(x)}\frac{\mathrm{\Delta }^2(x)}{ϵ_{0\mathrm{F}}},$$
(4.27)
where the prime denotes the derivative with respect to the dimensionless density $`x`$. The Fermi energy $`ϵ_\mathrm{F}(x)`$ and the pairing gap $`\mathrm{\Delta }(x)`$ entering these equations are determined from (4.6) and (4.9). The last two terms in (4.27), even in strong pairing regime, vanish in the dilute limit at least as $`x^q`$ if $`0<q<1`$ or linearly in $`x`$ if $`q1`$ (in our case, the exponent in the density-dependent pairing force is $`q=\frac{2}{3}`$). Thus, we see again that, in strong coupling, in the leading order $`\mu =ϵ_\mathrm{F}=ϵ_\mathrm{b}/2<0`$.

The calculated energy per nucleon as a function of the isoscalar density $`\rho `$ is shown in the upper panel in Fig. 3 together with the results of the nuclear matter calculations  for the UV14 plus TNI model. It is seen that DF3 gives qualitatively reasonable description of the nuclear matter EOS and that pairing could contribute noticeably to the binding energy especially at lower densities. In the lower panel in Fig. 3 we have plotted the pairing energy per nucleon, $`(E/A)_{\mathrm{pair}}`$, obtained by subtracting from eq. (4.26) the corresponding value of $`E/A`$ at $`\mathrm{\Delta }=0`$. This pairing energy is a sum of the positive contribution coming from the kinetic energy, i.e. from the first term in (4.26), and the negative anomalous energy – the last term in (4.26) (an expression for $`(E/A)_{\mathrm{pair}}`$ in the case of weak pairing is derived in Appendix B). The pairing contribution is small for the density-independent force with $`f^\xi =0.56`$, set (a) of eq. (4.5), and increases, as expected, for the sets (b)–(f) with a shift to lower densities as $`f_{\mathrm{ex}}^\xi `$ becomes gradually more attractive. For the sets (e) and (f) the attraction is strong, $`f_{\mathrm{ex}}^\xi <f_{\mathrm{cr}}^\xi `$. In these cases a nonvanishing binding energy in the dilute limit is solely due to Bose-Einstein condensation of the bound pairs, the spin-zero bosons, when all the three quantities, $`\mu `$, $`E/A`$ and $`(E/A)_{\mathrm{pair}}`$, reach the same value $`ϵ_\mathrm{b}/2`$ ($`ϵ_b=0.0646`$ and $`1.616`$ MeV for the set (e) and (f), respectively). This is illustrated in Fig. 4 where we have plotted $`E/A`$ and $`(E/A)_{\mathrm{pair}}`$ as functions of $`\rho `$ at very low densities. Analytically, for the kinetic energy term in (4.26) in the leading order we find

$$\frac{E_{\mathrm{kin}}}{A}=\frac{3\mathrm{\Delta }^2}{8xϵ_{0\mathrm{F}}^{3/2}}_0^{ϵ_\mathrm{c}}\frac{ϵ\sqrt{ϵ}ϵ}{(ϵ+|\mu |)^2}\frac{3\mathrm{\Delta }^2}{2xϵ_{0\mathrm{F}}}\left(\frac{1}{f_{\mathrm{cr}}^\xi }\frac{3\pi }{8}\frac{1}{k_{0\mathrm{F}}a}\right),$$
(4.28)
where we have used the definition (4.11) and the relation $`|\mu |=\mathrm{}^2/2ma^2`$. Combining this with the last term in (4.26), one gets

$$\frac{E}{A}=\frac{3\pi \mathrm{\Delta }^2}{16xϵ_{0\mathrm{F}}}\frac{1}{k_{0\mathrm{F}}a}.$$
(4.29)
By using eq. (4.21), this expression reduces exactly to $`\mathrm{}^2/2ma^2=ϵ_\mathrm{b}/2`$. It follows that the model just described is in fact parameter-free: in the strong coupling regime near the critical point, the ground state properties of nuclear matter in the dilute limit are completely determined by the scattering length.

We have considered subsaturated nuclear matter with $`{}_{}{}^{1}S_{0}^{}`$ pairing within our local EDF framework and demonstrated some results, including extrapolation to the low density limit, with a few possible parameter sets of the pairing force deduced from experimental data for lead isotopes. At low densities, in the $`T=0`$ case of symmetric $`N=Z`$ matter, the $`{}_{}{}^{3}S_{1}^{}^3D_1`$ pairing correlations leading to a Bose deuteron gas formation might be more important since the n–p force is more attractive than that in the p–p or n–n pairing channels (see Ref.  and references therein). Thus, our approach, with $`{}_{}{}^{1}S_{0}^{}`$ pairing only, would be more relevant for an asymmetric $`NZ`$ case and for pure neutron systems. From this point of view the best choice of the effective contact pairing force for the EDF calculations seems to be the set (d) of eq. (4.5) since it gives the singlet scattering length $`a_{\mathrm{nn}}17.2`$ fm which corresponds to a virtual state at $`140`$ keV known experimentally999The contact interaction in a truncated space of states, with two parameters $`f_{\mathrm{ex}}^\xi `$ and $`k_{0\mathrm{c}}`$, can be calibrated to produce not only the empirical scattering length $`a_{\mathrm{nn}}`$ but also a realistic effective range $`r_{\mathrm{nn}}2.8`$ fm which is directly related to the momentum space cutoff by $`r_{\mathrm{nn}}=4/\pi k_{0\mathrm{c}}`$, see Ref. . It would yield a rather low value for the energy cutoff $`ϵ_\mathrm{c}9`$ MeV. However, as we have already mentioned, even at low densities the force can be considerably modified by nonlocality and polarization effects so that the energy dependence of the Fermi-level $`s`$-wave nn phase shift, eq. (4.24), at low relative momenta $`k`$ governed by the effective range might be different from the free two-nucleon case.. As seen in Fig. 1, with this choice the behavior of $`\mathrm{\Delta }`$ at low densities agrees well with the calculations based on realistic NN forces. At higher densities, however, our predictions for $`\mathrm{\Delta }`$ with the set (d) go much higher reaching a maximum of $`4.84`$ MeV at $`k_\mathrm{F}0.92`$ fm-1 while the calculations of Ref.  give a maximum of about 3 MeV at $`k_\mathrm{F}0.82`$ fm-1. With a bare NN interaction, assuming charge independence and $`m_\mathrm{n}=m_\mathrm{p}`$ in the free single-particle energies, the pairing gap would be, at given $`k_\mathrm{F}`$, exactly the same both in symmetric nuclear matter and in neutron matter. As shown in Refs. , if one includes medium effects in the effective pairing interaction, most importantly the polarization RPA diagrams in the cross channel, the pairing gap in neutron matter would be substantially reduced to values of the order of 1 MeV at the most. Whether such a mechanism works in the same direction for symmetric nuclear matter is still an open question. If the gap were smaller than that obtained with the Gogny D1 force, which is also shown in Fig. 1, it would be difficult to explain the observed nuclear pairing properties. The effective contact density-dependent force (4.3) with our preferable parameter set (d) of eq. (4.5) yields larger pairing energy in nuclear matter than the Gogny force, but in finite nuclei this is compensated by the repulsive gradient term. The phenomenological pairing force used here contains dependence on the isoscalar density only since we have analyzed the existing data on separation energies and charge radii for finite nuclei with a relatively small asymmetry characterized by $`(NZ)/A0.25`$. An extrapolation to neutron matter with such a simple force would give a larger pairing gap than for nuclear matter. This suggests that some additional dependence on the isovector density $`\rho _\mathrm{n}\rho _\mathrm{p}`$ might be present in the effective pairing force. This possibility is planned to be tested in our future work, with a more careful analysis of experimental data though the relevant data base is not rich enough.

Now consider how the density changes when the pairing gap appears in nuclear matter. To leading order, one finds the following expression for the energy per nucleon near the saturation point:

$$\frac{E}{A}=\frac{E_0}{A}+\frac{K_0}{18}\frac{(xx_0)^2}{x_0^2}+\beta (x)I^2\frac{3}{8}\frac{\mathrm{\Delta }^2(x)}{ϵ_\mathrm{F}(x)}$$
(4.30)
where $`K_0`$ is the compression modulus at saturation density, $`\beta (x)I^2`$ is the asymmetry energy, with $`I=(\rho _\mathrm{n}\rho _\mathrm{p})/2\rho _0x(NZ)/A`$. This expression is valid at $`|xx_0|1,|I|1`$ and $`\mathrm{\Delta }ϵ_\mathrm{F}`$. Higher order effects connected with $`I`$-dependence of $`K,\beta ,\mathrm{\Delta }`$ and $`ϵ_\mathrm{F}`$ are neglected. The derivation of the last term in (4.30) is discussed in detail in Appendix B. Due to this pairing term, the position of the equilibrium point may be shifted to lower or higher densities depending on the behavior of $`\mathrm{\Delta }(x)`$ near $`x=1`$. If $`\mathrm{\Delta }`$ does not depend on $`x`$ in the vicinity of $`x=1`$ then the equilibrium density decreases due to presence of $`ϵ_\mathrm{F}(x)x^{2/3}`$ in the denominator of the pairing term (the system gains more binding energy). This means an expansion of the system. The effect is enhanced if $`\mathrm{\Delta }`$ becomes larger during such an expansion, i.e. when the derivative $`\mathrm{\Delta }(x)/x`$ at $`x=1`$ is negative. This point may be illustrated by the following simple consideration. At equilibrium the pressure $`P=0`$ which means

$$\frac{}{x}(E/A)=0.$$
(4.31)
For saturated symmetric nuclear matter without pairing the dimensionless density is, by definition, $`x_0=1`$. When $`I0`$ and $`\mathrm{\Delta }0`$, from eq. (4.30) with condition (4.31) the new equilibrium density can be found (in units of $`2\rho _0`$), which is the solution of the equation

$$xx_0=\frac{9x_0^2}{K_0}\left[\frac{\mathrm{\Delta }(x)}{4ϵ_\mathrm{F}(x)}\left(3\frac{\mathrm{\Delta }(x)}{x}\frac{\mathrm{\Delta }(x)}{x}\right)\frac{\beta (x)}{x}I^2\right].$$
(4.32)
From this equation one can see that the density should be sensitive indeed to the derivatives of $`\mathrm{\Delta }`$, and a negative slope in $`\mathrm{\Delta }`$ should cause a decrease of the density. The obtained relations can be used to estimate the influence of pairing on the charge radii for heavy nuclei as was done in Ref. . It was shown that pairing interaction with strong $`\rho `$-dependence at $`x1`$ might significantly change the equilibrium density. For the parametrization used here the size of this effect is controlled by the parameter $`h^\xi `$. The shift of the saturation point is relatively small: $`|\delta \rho |/2\rho _00.8`$ % for all the six sets of eq. (4.5) ($`\delta \rho /2\rho _0+0.4`$ % and $`0.4`$ % for the set (a) and (d), respectively). In finite nuclei, the surface term $`f_{}^\xi `$ is equally important to produce a kink in the radius evolution along isotope chain at magic neutron number and, especially, to explain the observed odd-even staggering in $`r_{\mathrm{ch}}^2`$ . Although, as seen in Fig. 1, the calculations with bare NN interaction or with Gogny force give a negative sign for $`\mathrm{\Delta }(x)/x`$ near $`x1`$ ($`k_\mathrm{F}1.33`$ fm-1), but the slope might be not steep enough. This is the probable reason why the HFB calculations with the Gogny force, which give a good description of the global pairing properties of nuclei, could not reproduce the kink in lead isotopes . The density dependence of the pairing force leads to the direct coupling between the neutron anomalous density and the proton mean field as given by (4.4). The suppression of $`|\nu _\mathrm{n}|^2`$ in the odd neutron subsystem because of the blocking effect influences the potential $`U_{\mathrm{pair}}^\mathrm{p}`$ of the proton subsystem through the volume, $`h^\xi `$, and surface, $`f_{}^\xi `$, couplings, and this moves the behavior of the proton radii towards the desired regime .

## 5 Numerical results for some isotope chains

The calculations for finite spherical nuclei are performed using the coordinate-space technique which is described in detail in Appendix C. In the construction of the Gor’kov Green’s functions, the four linear-independent solutions $`(u_i,v_i),i=1`$–4, were found by the Numerov method with a radial step of 0.1 fm and physical boundary conditions imposed for each $`lj`$ channel at the origin and at $`r=25`$ fm. The contour in the complex energy plane with Im$`ϵ=\pm 8`$ MeV along the horizontal sections and with energy cutoff $`ϵ_\mathrm{c}=40`$ MeV measured from the Fermi level was used (see Figs. 21, 22 in Appendix C). The integration along the contour was performed by the Simpson method with automatic step selection. The convergence of the iteration procedure is controlled by a few criteria: for two successive steps $`i`$ and $`i+1`$, the conditions $`|\rho _{i+1}(r)\rho _i(r)|<10^6`$ fm-3 and $`|\mathrm{\Delta }_{i+1}(r)\mathrm{\Delta }_i(r)|<50`$ keV should be achieved for all $`r`$, the chemical potential $`\mu `$ should finally satisfy the condition $`|N(\mu )N|<0.01`$ and the contribution $`N_{lj}`$ to the total particle number $`N(\mu )`$ from the states with higher angular momenta $`lj`$ should not exceed 0.01. Such criteria guarantee an accuracy not worse than 0.1% for all the calculated quantities of our interest. The mean square charge radii were computed from ground state charge densities obtained by folding the point nucleon distributions with nucleon charge form factors, including the relativistic electromagnetic spin–orbit correction, in the same way as in Ref. .

Let us discuss first the results obtained for the lead chain. As shown in the upper panel of Fig. 5, the experimental neutron separation energies are reproduced equally well for all the chosen parameter sets of the pairing force. The curves (a)–(f) correspond to the sets marked with the same letters in eq. (4.5). The curve (a) is obtained for the “constant” pairing force, without density dependence. The other sets (b)–(f) differ from that by non-zero values of $`h^\xi `$ and $`f_{}^\xi `$ parameters. As already mentioned in Section 4, for each $`f_{\mathrm{ex}}^\xi `$ value it is possible to find such values of $`h^\xi `$ and $`f_{}^\xi `$ that the pairing energy for a given nucleus is practically the same. These parameters are kept fixed for all isotopes of the lead chain.

In the lower panel of Fig. 5 the calculated isotope shifts of mean squared charge radii, $`\delta r_{\mathrm{ch}}^2`$, for the lead chain with respect to 208Pb as a reference nucleus are presented. These curves except (a) do not differ from each other significantly and all of them reproduce qualitatively the kink and the average size of odd-even staggering. The curve (a) corresponds to a simple pairing delta-force; in this case the behavior of $`\delta r_{\mathrm{ch}}^2`$ is rather smooth and neither kink nor staggering are reproduced. The variant (a) of eq. (4.5) was included just to give an example that without any density dependence in the contact pairing effective interaction it is also possible to describe the neutron separation energies but not the evolution of charge radii. With the parametrization of eq. (4.3) used here for the pairing force, the size of the staggering and the kink in charge radii of Pb isotopes can be reproduced only if the “gradient” parameter $`f_{}^\xi `$ is $`2`$. Without gradient term the results for radii would be in between the curves (a) and (b) shown in the lower panel of Fig. 5 (an example is given in the lower panel of Fig. 6). On the whole, as seen in Fig. 5, the set (d) yields a somewhat better description than the other sets, particularly for the lighter Pb isotopes. For this reason, and due to the fact that the set (d) corresponds to the correct value of the singlet scattering length in the dilute limit (see previous Section), this set will be taken as our preferable choice, and the corresponding functional will be referred to as DF3(d) in the following.

In Fig. 6 we compare the DF3(d) predictions for Pb isotopes with the HFB results obtained for the Gogny D1S force  and for two state-of-the-art models based on Skyrme density-dependent forces SkP and SLy4 . The calculations with Gogny force for the odd isotopes are done in the blocking approximation, and thus we can show the results both for even and odd Pb isotopes101010We thank Jean-François Berger, Jacques Dechargé and Sophie Peru for providing us with these results.. The Skyrme-HFB calculations111111We thank Jacek Dobaczewski for providing us with a numerical file of these Skyrme–HFB calculations. were done, unfortunately, without blocking, and therefore we show the radii only for even isotopes; the $`S_n`$ values for odd isotopes were obtained by correcting the no-blocking HFB energy by adding an average pairing gap, this was considered to be a very good approximation. The details of the HFB+SkP and HFB+SLy4 calculations can be found in  (see also the references therein), and here we only briefly discuss how the pairing correlations have been implemented in these models. The HFB+SkP calculations have been done with the same SkP force in the ph and pp channels, and therefore in these calculations the pairing force has density dependence $`x^{1/6}`$. The energy cutoff has been chosen by including the quasi-particle states up to the energy equal to the depth of the effective single-particle potential . This implies that the cutoff depends not only on the particle numbers $`N`$ and $`Z`$ but also on the quantum numbers $`j`$ and $`l`$; for $`l`$=0 it varies in the range of about 40–50 MeV. The HFB+SLy4 calculations have been done with the density-dependent zero-range pairing force in the pp channel, and the SLy4 force in the ph channel. The form of this pairing force is identical to the velocity-independent piece of the Skyrme interaction (see Ref. ) with parameters $`V_0=t_0=2488.913`$ MeV fm3 and $`V_3=19990`$ MeV fm3+1/6, and it also has the $`x^{1/6}`$ density dependence. The prescription for the energy cutoff has been identical as in the HFB+SkP calculations. We remark that such a prescription is much different from that which we use here; with such a recipe it is not easy to construct the corresponding effective contact pairing force with a fixed phase-space cutoff, hence no simple extrapolation to the pairing properties of uniform nuclear matter.

One can see in the upper panel of Fig. 6 that the HFB+D1S, HFB+SkP, SkLy4 and our DF3(d) calculations describe the neutron separation energies with more or less the same quality, some deviations are observed in the region above 207Pb. For the lighter Pb isotopes, the Gogny force slightly overestimates the odd-even effect in $`S_n`$ while with DF3(d) it is slightly underestimated.

In the lower panel of Fig. 6 it is seen that the Skyrme-HFB calculations with SkP and SLy4 forces give too small kink in radii. Calculations with the Gogny D1S force for $`\delta r_{\mathrm{ch}}^2`$ could not reproduce the kink either. The latter model yields sizeable staggering in radii, but the effect is too small, at least by a factor of two smaller in amplitude than experimentally observed; moreover, its sign in isotopes below $`A`$=194 is reversed with respect to the observations. As for the Skyrme functionals, one suspects that they could not give the correct size of staggering as well. Such a conclusion may be supported by our DF3 calculations with parameters $`f_{\mathrm{ex}}^\xi =1.79,h^\xi =1.66,f_{}^\xi =0`$, i.e. with the pairing force of (4.5) which contains only a $`x^{2/3}`$ dependence, without gradient of density. The external strength constant $`f_{\mathrm{ex}}^\xi =1.79`$ is taken from our preferable set (d), and the parameter $`h^\xi =1.66`$ is found to get practically the same $`S_n`$ values as with the variant (d) itself. The results for the radii are shown in the lower panel of Fig. 6 by the dotted line. It is seen that the staggering appears but it is too small indeed. With a $`x^{1/6}`$ dependence, which is used with the SkP and SLy4 parametrization, the effect would be even smaller. Thus, to reproduce the scale of the odd-even effects in radii (and the scale of the kinks) within the local energy-density functional approach, the blocking mechanism should be enhanced by some peculiar density dependence of the effective pairing force, the linear dependence on $`\rho ^\sigma `$ with $`\sigma `$=1, 1/3, 1/6… could not give the desirable size.

This fact clearly demonstrates the importance of including odd-mass nuclei in the fitting procedure of force constants. Now, in general, a complete description of the ground state of an odd nucleus is quite difficult which reflects itself e.g. in the difficulty to reproduce experimental magnetic moments. These depend sensitively on the time-reversal-odd components of the polarization induced by the odd particle. However, the operator $`r^2`$ is time-reversal-even, and the small, not-well-known, odd components of the polarization enter the expectation value only quadratically. If, with given interaction, $`r^2`$ is not well described by a HFB-type ground state, there is no hope to reproduce it with a more sophisticated state vector. As can be seen from Figs. 5 and 6, it is possible to reproduce the order of magnitude of the staggering of $`r^2`$ with density-dependent pairing force (4.3) containing a gradient term.

Similar arguments apply to the nuclei off magic numbers. Looking only at even neutron isotopes, the slope of $`r^2`$ as a function of $`N`$ changes at the magic numbers, producing a “kink”. With simple interactions in the pp-channel, this is not reproduced in HFB-type calculations121212A kink has, however, been obtained in relativistic mean-field calculations, with simple pairing in the Pb chain , at the expense of significantly too weak binding of the neutron single particle states above the Fermi level of 208Pb which increases the polarization effects through neutron–proton interaction (see Ref.  for a detailed discussion of this point, and also our commentary to Figs. 9 and 10 of the present paper).. This relative increase of $`r^2`$ off magic numbers is thought to be connected with “dynamical deformation” due to ground state fluctuations of the surface. These ground state fluctuations of the surface, being time-independent, actually can not be separated from a static diffuseness of the surface, and are present in an independent particle model too. In Section 6 below we will give some estimates of the contribution of RPA-type ground state correlations to this effect. It emerges that the main part of the effect should—and can—be obtained on the HFB level.

In Fig. 7 the $`r`$-dependence of $`\mathrm{\Delta }`$ for two isotopes 190Pb and 214Pb (the ends of the measured chain) is shown. One can see that with increasing external attraction $`f_{\mathrm{ex}}^\xi `$, $`\mathrm{\Delta }`$ becomes smaller in the volume and gradually concentrates in the outer part of the nuclear surface. This is demonstrated more clearly in Fig. 8 were $`\mathrm{\Delta }`$ is plotted as a function of dimensionless isoscalar density $`x`$. The curves (c)–(f) are crossing very nearly at the same point which has approximately the same coordinates in both cases: $`x_{\mathrm{piv}}0.7`$ and $`\mathrm{\Delta }_{\mathrm{piv}}0.7`$$`0.8`$ MeV. This “pivoting” point corresponds to the one seen in Fig. 1 for $`\mathrm{\Delta }`$ in nuclear matter but in finite nuclei the value of $`\mathrm{\Delta }_{\mathrm{piv}}`$ turns out to be lower by a factor of 4 due to the repulsive gradient term $`f_{}^\xi `$.

It is of interest to compare the predictions obtained with different mean-field approaches not only for the proton radii but also for the neutron ones. An example is given in Fig. 9 where the results for the lead isotopes obtained with the functional DF3(d) are shown in comparison with the HFB+SkP and HFB+SLy4 calculations, and with recently published relativistic mean-field Hartree calculations  based on the effective force NL3 and the constant pairing gap prescription of Ref.  within the Bardeen-Cooper-Schriffer formalism (RMF–HBCS). It is seen that the EDF model with gradient pairing yields staggering in the evolution of the proton and neutron radii, and produces a kink in both of them at A=208. The Skyrme-type functionals produce very small kinks in proton radii. The calculations with the RMF–HBCS model are given in Ref.  for even-even nuclei only, but one would guess that this model is not able to produce noticeable odd-even staggering too. At the same time, as seen in Fig. 9, both Skyrme–HFB and the RMF–HBCS functionals predict bigger neutron radii than our EDF calculations, and all of them yield a distinctive kink in the evolution of $`r^2_\mathrm{n}`$ at $`A=208`$. One observes rather close agreement between the proton radii obtained with all considered models but a wide spread in the neutron radii. The RMF–HBCS neutron radii are particularly large. We remark at this point that the “experimental” rms neutron radius in 208Pb of 5.593 fm deduced from the analysis of the high-energy polarized proton scattering in Ref.  is significantly lower than the prediction of the RMF(NL3) model; this fact has been already mentioned in Refs. . But one should bear in mind that the result of this analysis is model-dependent, and the errors in the extracted neutron radii could be quite large. The difference between neutron and proton rms radii for 208Pb, +0.14$`\pm `$0.04 fm, has been also deduced in Ref. , and we just put the corresponding error bars for the $`r^2_n^{1/2}`$ value in this nucleus as shown in Fig. 9.

Considering general trends in the behavior of a given sequence of nuclei, we notice that, because of the strong pn attraction, the evolution of the proton radius should be driven by that of the neutron one through self-consistency, and vice versa. Moreover, in neutron-rich nuclei, an anticorrelation should exists between neutron radii and neutron separation energies: the larger $`r^2_\mathrm{n}`$, the smaller $`S_{2\mathrm{n}}`$. One may suspect that too big neutron radii with a kink at $`N=126`$, which is correlated with a kink in proton radii, would mean too small $`S_{2\mathrm{n}}`$ values beyond the shell closure. This is indeed the case as clearly demonstrated in Fig. 10. This figure displays two-neutron separation energies $`S_{2\mathrm{n}}`$ (upper panel) and differences of mean square charge radii $`\delta r_{\mathrm{ch}}^2`$ with respect to 208Pb (lower panel) for the even lead isotopes. The results obtained with the RMF(NL3), Skyrme-SLy4 and EDF-DF3(d) models are compared there with experiment. It is seen that all these models reproduce $`S_{2\mathrm{n}}`$ values for lighter Pb isotopes below $`A=208`$ fairly well, with more or less the same quality, the deviation from experiment being within 1 MeV or so (one should bear in mind that $`S_{2\mathrm{n}}`$ for the $`A=190`$$`198`$ lead isotopes are derived in Ref.  from systematic trends, and the models discussed here describe these trends reasonable well). But for the heavier lead isotopes the predictions are different. Our EDF calculations with gradient pairing agree with experimental $`S_{2\mathrm{n}}`$ values within 0.5 MeV while Skyrme-SLy4 and RMF(NL3) models give deviations up to 1.5 and 2.0 MeV, respectively. It is also seen in the lower panel of Fig. 10 that the isotope shifts in charge radii for lighter Pb isotopes are nicely reproduced with our EDF approach, the $`\delta r_{\mathrm{ch}}^2`$ values for the $`A>208`$ nuclei and the kink being slightly underestimated. The Skyrme-SLy4 and RMF(NL3) models work worse and yield results of equal quality for the lighter isotopes, but produce much different kinks and $`\delta r_{\mathrm{ch}}^2`$ values beyond $`A=208`$. The kink obtained with the RMF model incorporating simple BCS pairing, with constant pairing gap, is impressive indeed, and good description of this anomaly in the isotope shifts of Pb nuclei has been reputed as a considerable success of the relativistic mean-field approach. The above consideration, which has much in common with that of Ref. , points out, however, the importance of simultaneous description of both physical quantities, i.e. the energetic ($`S_{2\mathrm{n}}`$) and geometrical ( $`\delta r_{\mathrm{ch}}^2`$) differential observables, in order to get a deeper insight into the nature of this anomaly in the isotope shifts.

The calculated neutron separation energies and differences of mean squared charge radii with respect to 116Sn for tin isotopes with mass number from A=99 to A=136 are shown in Fig. 11 in comparison with experimental data and with the Skyrme–HFB calculations based on the SkP and SLy4 forces. In the lower panel of this figure, we also show the RMF–HBCS $`\delta r_{\mathrm{ch}}^2`$ values . The results for radii from these three models are available for even isotopes only. Of the six parameter sets (4.5) of the pairing force, only the results for our preferable set (d), with a scaling factor of 1.05, are shown. Similar to the lead chain, the $`S_\mathrm{n}`$ values for tin isotopes are reproduced reasonably well. However, on both ends of the considered chain, near the magic nuclei 100Sn and 132Sn, the common trend for all mean-field calculations presented in the upper panel of Fig. 11 is that the size of the odd-even effect in neutron separation energies becomes noticeably smaller than the experimental one. As seen in the lower panel in Fig. 11, the evolution of the charge radii in the Skyrme-type and RMF–HBCS calculations is rather smooth, the SkP, SLy4 and RMF functionals slightly overestimate $`\delta r_{\mathrm{ch}}^2`$ in the heavier tin isotopes above $`A=120`$. The desirable size of staggering in the behavior of $`\delta r_{\mathrm{ch}}^2`$ for tin isotopes could be obtained with gradient pairing as shown by open circles in this figure. At the same time, the EDF calculations with functional DF3 significantly underestimate $`\delta r_{\mathrm{ch}}^2`$ in lighter isotopes below A=114. This fact seems to be correlated with the weakening of the pairing when approaching A=100 and points out the need of both more careful adjustment of the normal isovector part of the energy-density functional and, probably, invoking the dependence of the pairing force on the isovector density $`\rho _{}=\rho _\mathrm{n}\rho _\mathrm{p}`$.

The calculated two–neutron separation energies for even tin isotopes up to the neutron drip line are shown in Fig. 12. Again, the set (d) of the pairing force with a scaling factor of 1.05 was used in the EDF calculations with the parametrization DF3. One can observe noticeable differences between the Skyrme–HFB calculations with the SkP and SLy4 forces, the RMF–HBCS predictions, and the EDF results. The RMF-HBCS calculations  are performed with deformed code, and available for the tin isotopes only up to $`A=156`$; these calculations produce some irregularities in the $`S_{2\mathrm{n}}`$ values near $`A=100`$ and $`A=132`$ (the nuclei beyond 140Sn are predicted to be well deformed). At the same time, the three other spherical microscopic calculations give surprisingly similar predictions when approaching the neutron drip line, the heaviest bound tin nucleus being 172Sn, 174Sn and 176Sn for the EDF, SkP and SLy4, respectively. One should notice that the Skyrme–HFB calculations were performed with a discretized continuum in a spherical box of finite size  which allows to artificially obtain a stable solution for the nuclear ground state even for positive value of the chemical potential. Such a solution is not possible with the coordinate-space technique used in the present paper because of the physical boundary conditions imposed for the scattering states.

The evolution of the isoscalar and isovector densities along the tin isotope chain is illustrated in Fig. 13. It is seen that isovector density $`\rho _{}=\rho _\mathrm{n}\rho _\mathrm{p}`$ increases both in the volume and at the surface as the neutron number $`N`$ becomes larger. The isoscalar (matter) density tends to be slightly decreased with $`N`$ in the volume and became more diffused at the surface. The matter distributions in the two magic nuclei, 100Sn and 132Sn, have a steeper slope at the surface compared to the non-magic isotopes. The radial dependence of the neutron pairing field $`\mathrm{\Delta }(r)`$, of the anomalous density $`\nu (r)`$ and of the normal neutron density $`\rho _\mathrm{n}(r)`$ for the stable nucleus 116Sn, for the unstable isotope 150Sn and for the drip-line nuclide 172Sn are shown in Fig. 14. One can see that in the heavier isotopes, the anomalous density becomes more diffuse at the surface with a longer tail compared to the normal densities, in agreement with the discussion given in Appendix C. It is also seen that the pairing field $`\mathrm{\Delta }`$ in the outer part of the nuclear surface has a prominent maximum which moves outwards and gets a larger tail near the drip line. The concentration of $`\mathrm{\Delta }`$ outside of the half density point is related to the specific cancellation between the external attraction and the repulsive gradient term in the effective pairing force .

The staggering phenomenon observed in the behavior of charge radii plotted as a function of neutron number is the one of the prominent odd-even effects which has been systematically measured and widely discussed. Similar effects are expected to exist for other quantities such as neutron and matter radii, centroid energies of multipole excitations (position of the giant resonances), etc.

Shown in the upper panel of Fig. 15 are the root-mean-square proton, neutron and matter radii for the tin isotope chain calculated with the functional DF3 and with the set (d) of the pairing force. In the middle panel of this figure we show the mean energies $`\omega _{\mathrm{dip}}`$ of dipole isovector excitations. These are calculated as the square root of the ratio $`m_3/m_1`$ with $`m_1`$ and $`m_3`$ the linear and cubic energy-weighted sum rule, i.e. the first and the third moment of the corresponding RPA strength distribution, respectively :

$$\omega _{\mathrm{dip}}=\left[\frac{\mathrm{}^2}{3m}\frac{A}{NZ}\stackrel{}{r}\stackrel{}{r}^{}\frac{\rho _\mathrm{n}}{\stackrel{}{r}}^{\mathrm{np}}(\stackrel{}{r},\stackrel{}{r}^{})\frac{\rho _\mathrm{p}}{\stackrel{}{r}^{}}\right]^{1/2}.$$
(5.1)
Here $`^{\mathrm{np}}`$ is the effective neutron–proton interaction obtained from the energy-density functional as the second variational derivative $`\delta ^2E_{\mathrm{int}}/\delta \rho _\mathrm{n}\delta \rho _\mathrm{p}`$. In the lower panel of Fig. 15 we show the mean energies $`\omega _{\mathrm{Fermi}}`$ of the charge-exchange $`0^+`$ excitations, i.e. the Fermi transitions, in tin nuclei with $`N>Z`$. These energies are calculated within the sum rule approach as the ratio $`(m_1^++m_1^{})/(m_0^+m_0^{})`$ with $`m_1^+`$ and $`m_1^{}`$ the first moment of the strength distribution of the Fermi transitions in the $`\beta ^+`$ and $`\beta ^{}`$ channel, respectively (energy-weighted sum rules) and with $`m_0^+`$ and $`m_0^{}`$ the corresponding non-energy weighted sum rules. The resulting expression for $`\omega _{\mathrm{Fermi}}`$ is given by (see, for example, ):

$$\omega _{\mathrm{Fermi}}=\frac{1}{NZ}\stackrel{}{r}U_{\mathrm{Coul}}(\stackrel{}{r})\left(\rho _\mathrm{n}(\stackrel{}{r})\rho _\mathrm{p}(\stackrel{}{r})\right),$$
(5.2)
where $`U_{\mathrm{Coul}}`$ is the Coulomb mean field potential.

It is seen in Fig. 15 that the rms neutron, proton and matter radii as functions of the mass number $`A`$ reveal a kink at magic 132Sn, and at larger $`A`$ the difference between rms neutron and proton radii starts to increase more rapidly. The staggering in radii is present mostly in the region between the two magic nuclei, from 100Sn to 132Sn, and this effect is practically washed out beyond $`A=140`$. The mean energy of dipole transitions occurs to be in anticorrelation with such a behavior in radii, and this seem to be in qualitative agreement with experimental data (one should mention that eq. (5.1) overestimates the position of the giant dipole resonance by $`0.5`$ MeV because of the cubic energy-weighted sum rule $`m_3`$ used in the derivation of $`\omega _{\mathrm{dip}}`$). One also observes a distinct kink in the behavior of $`\omega _{\mathrm{dip}}`$ at $`A=132`$. Beyond this magic number the mean dipole energy decreases rather fast and then nearly saturates when approaching $`A=172`$. This might be connected with an enhancement of the low-energy dipole transitions  and also with possible appearance of the so-called soft dipole mode  in nuclei near the neutron drip line. The anticorrelations in the staggering behavior of $`\omega _{\mathrm{dip}}`$ and the rms radii $`r^2_\mathrm{n}^{1/2}`$, $`r^2_\mathrm{p}^{1/2}`$ can easily be understood by considering the influence of pairing on the gradients of neutron and proton densities entering eq. (5.1). The odd-even effect in $`\omega _{\mathrm{dip}}`$ is constructive and more pronounced in the $`A132`$ region where both gradients strongly overlap. Beyond 132Sn, the larger differences between neutron and proton rms radii imply lower overlap between neutron and proton density gradients at the nuclear surface, hence a smaller mean dipole energy. The situation with mean energy of the Fermi charge-exchange transitions is different. From eq. (5.2) one expects that the correlation between $`\omega _{\mathrm{Fermi}}`$ and the neutron and proton rms radii should be destructive. As seen in the lower panel in Fig. 15, the staggering in the evolution of $`\omega _{\mathrm{Fermi}}`$ with $`A`$ is very weak indeed and almost invisible. The kink at the magic mass number $`A=132`$ is also much less pronounced compared to the dipole case. Remarkable enough, the theoretical self-consistent sum-rule predictions for $`\omega _{\mathrm{Fermi}}`$ in tin isotopes are in excellent agreement with the available experimental data on the position of the isobaric analog states.

Fig. 16 displays the neutron separation energies and isotope shifts in mean square charge radii for the calcium isotope chain. In this case, the EDF calculations with the preferable set (d), eq. (4.5), deduced from the Pb isotopes yield noticeable but rather small staggering both in $`S_\mathrm{n}`$ and $`\delta r_{\mathrm{ch}}^2`$. A possible way to get the desirable size of odd-even effects is to enhance the pairing correlations by introducing a scaling factor 1.35 for the pairing force (all parameters being scaled in the same way). It is seen then in Fig. 16 that the neutron separation energies are reproduced fairly well, with more or less the same quality as with the Skyrme SkP and SLy4 functionals, though the effect of the scaling is rather mild. The effect on the behavior of charge radii is more drastic: the $`\delta r_{\mathrm{ch}}^2`$ values are increased by a factor of 3 to 4 yielding a very good agreement with experiment. This example demonstrates that the anomalous $`A`$-dependence of charge radii in Ca isotopes could be correctly described within the EDF approach with renormalized pairing force. However, the possibility to determine a single parametrization of the effective interaction in the pp channel, which would be valid both for heavy and light nuclei, remains an open problem. In any case, we got an indication that the staggering in $`\delta r_{\mathrm{ch}}^2`$ of Ca isotopes, which is a challenge for nuclear models, could be explained with a peculiar density dependence of the pairing force. One can see in Fig. 16 that the Skyrme SkP and SLy4 functionals produce a nearly smooth increase of radii with $`A`$, with weak irregularities below and above $`A=48`$, but still with rather sizeable kinks at double magic 48Ca. The RMF–HBCS predictions for $`\delta r_{\mathrm{ch}}^2`$, Ref. , are much different: the charge radius is nearly a constant up to 48Ca, with a sudden increase for 50Ca. Interestingly, the functional DF3 yields, in agreement with experiment, a small negative isotope shift in charge radii between 40Ca and 48Ca, which does not depend on pairing.

The EDF calculations, as demonstrated in Fig. 17, show that the trends and the staggering effects in the mean energy of isovector dipole excitations and charge-exchange Fermi transitions for calcium isotopes are in strong anticorrelation with evolution of the ground state mean squared radii. Similar to the case of the tin isotopes discussed above, the anticorrelation in $`\omega _{\mathrm{Fermi}}`$ is weaker than in $`\omega _{\mathrm{dip}}`$. Unfortunately, there are no systematic experimental data on dipole excitations in Ca isotopes at our disposal, from which the $`\omega _{\mathrm{dip}}`$ values could be deduced. As for the mean energies of the Fermi transitions, which are calculated with the self-consistent EDF sum rule approach and plotted as a function of $`A`$ in the lower panel in Fig. 17, they may be compared with the experimental positions of the isobaric analog states. It is seen that the sum rule expression, eq. (5.2), underestimates the energy of the analog states in calcium isotopes by some 200 keV but reproduces qualitatively the observed trend.

As, in the past, RPA ground state correlations have been invoked to explain details of the isotopic dependence of nuclear charge radii, we estimate the influence of RPA ground state correlations.

## 6 Ground state fluctuations

An important point is the issue of RPA-type ground state correlations. In several papers  the influence of the average mean square deformations on the mean square charge radii has been discussed in a phenomenological way which is based on the following sum rule: Take the identity for the multipole operator $`Q_{\lambda ,\mu }=\stackrel{}{r}r^\lambda Y_{\lambda ,\mu }(\mathrm{\Omega })\psi ^{}(\stackrel{}{r})\psi (\stackrel{}{r})`$

$$\underset{\mu }{}0||Q_{\lambda ,\mu }|^2|0=\underset{\mu ,s}{}|s|Q_{\lambda ,\mu }|0|^2.$$
(6.1)
The right hand side is equal to the summed transition strength $`_sB(\mathrm{E}\lambda ;0s)`$, the left hand side is the expectation value of the squared multipole operator $`Q_\lambda `$ in the ground state. In semiclassical nuclear models, this expectation value is proportional to the square of the deformation of multipolarity $`\lambda `$. Thus, even for spherical nuclei the rms deformation does not vanish. In the liquid drop model this is interpreted as due to zero point surface vibrations, but it is to be emphasized that there is no observable time dependence, and no breaking of the spherical symmetry due to this “dynamical” deformation. Experimentally, the effect can not be distinguished from diffuseness of the surface arising from any other mechanism.

Our method (as well as any HF or HFB calculation) incorporates only the dynamical deformation effects corresponding to a gas of noninteracting particles (magic nuclei) or of a gas with only pairing (nonmagic nuclei). To include the RPA ground state correlations, the calculated charge radii should be increased by

$$\delta r_{\mathrm{ch}}^2=\overline{\delta }r_{\mathrm{ch}}^2_{\mathrm{RPA}}\overline{\delta }r_{\mathrm{ch}}^2_{\mathrm{HFB}}.$$
(6.2)
However, the parameters have been adjusted to reproduce the experimental density distribution and therefore the average RPA correlations are included due to the choice of parameters. In a high precision fit, this should be corrected for, by subtracting the value given by (6.2) for the reference nucleus.

A microscopic method to calculate the corrections to $`r_{\mathrm{ch}}^2`$ due to low-energy surface vibrations has been proposed in . We use the liquid drop model and eq. (6.1) to estimate the influence of isoscalar collective excitations on the charge radii, as proposed by Esbensen and Bertsch .

As already noted in , the inclusion of noncollective states in eq. (6.1), which can not be considered to be surface oscillations will not do much harm, because the contribution of these states will be essentially the same to $`\overline{\delta }r_{\mathrm{ch}}^2_{\mathrm{RPA}}`$ and $`\overline{\delta }r_{\mathrm{ch}}^2_{\mathrm{HFB}}`$ and will drop out. According to , for multipolarity $`\lambda >0`$

$$\overline{\delta }r_{\mathrm{ch}}^2r^2_\mathrm{d}r^2_\mathrm{s}=\frac{5}{4\pi }r^2_\mathrm{s}\underset{\lambda }{}\beta _\lambda ^2,$$
(6.3)
where the suffixes d and s mean (dynamically) deformed and spherical, respectively.

Within the same model

$$S^0(\mathrm{E}\lambda )\underset{s}{}B(E\lambda ;0s)=\underset{\mu }{}0||Q_{\lambda ,\mu }|^2|0=\left(\frac{3ZeR_0^\lambda }{4\pi }\right)^2\beta _\lambda ^2.$$
(6.4)
With the help of (6.4) and (6.1) the sum over $`\beta _\lambda ^2`$ in eq. (6.3) can be replaced by the sum over the total $`B(\mathrm{E}\lambda )`$-values, and we reach the following expression:

$$\delta r_{\mathrm{ch}}^2=\frac{4\pi }{3Z^2e^2}\underset{\lambda }{}\frac{1}{R_0^{\lambda 2}}[S_{\mathrm{RPA}}^0(\mathrm{E}\lambda )S_{\mathrm{HFB}}^0(\mathrm{E}\lambda )].$$
(6.5)
Of course, for $`\lambda =1`$ one has to exclude the spurious translation. For $`\lambda =0`$ we have to use a compressible drop model, and obtain similarly

$$\delta r_{\mathrm{ch}}^2=\frac{140\pi }{108Z^2e^2R_0^2}[S_{\mathrm{RPA}}^0(\mathrm{E0})S_{\mathrm{HFB}}^0(\mathrm{E0})].$$
(6.6)
where $`S(\mathrm{E0})`$ is calculated according to eq. (6.4) with $`Q_{00}=r^2Y_{00}`$.

The excited states have been calculated in quasiparticle-RPA (QRPA) which has been used in the form (we adopt the same symbolic notations as in Ref. ):

$$\widehat{V}=e_q\widehat{V}_0+\widehat{F}\widehat{A}\widehat{V},S_{\mathrm{RPA}}^0=\frac{1}{\pi }\mathrm{Im}(e_q\widehat{V}_0\widehat{A}\widehat{V}),$$
(6.7)
where $`e_q`$ is the quasiparticle local charge with respect to the external field $`V_0`$. For a long-range electric field, $`V_0Q_{\lambda ,\mu }`$, from gauge invariance one gets $`e_q=1`$ . The $`\widehat{A}`$ matrix is written in detail in .

As in our formalism the effective interaction amplitudes $`\widehat{F}`$ are obtained as the second functional derivatives of $`E_{\mathrm{int}}`$ with respect to the corresponding densities, one can see here that due to $`\rho `$-dependence of $`\epsilon _{\mathrm{anomal}}`$, eq. (4.1), there are mixed terms for the effective interaction, $`F^{\omega \xi }`$ and $`F^{\xi \omega }`$, which do not vanish, and which couple the ph with the pp or hh channel ($`\tau `$-indices are neglected for simplicity):

$$^{\omega \xi }=\frac{\delta ^2E_{\mathrm{int}}[\rho ,\nu ]}{\delta \rho \delta \nu }=\frac{\delta ^2}{\delta \rho \delta \nu }\epsilon _{\mathrm{anomal}}([\rho ,\nu ];\stackrel{}{r})\stackrel{}{r},$$
(6.8)
A comparison of QRPA results for the first collective $`2^+`$ and $`3^{}`$ states with experiment is presented in Table 1 for selected nuclei from the Ca, Sn and Pb isotope chains. The influence of the spin–orbit part in the ph interaction has been investigated in ; the results given here have been obtained without this spin–orbit contribution. Comparison with  shows that here the changes in the functional did not improve the agreement with experimental data.

We restricted our consideration to $`2^+`$ and $`3^{}`$ states for all isotope chains because the other multipolarities contribute much less to the mean square charge radii . This has been checked by computing also $`\delta r_{\mathrm{ch}}^2`$ from $`0^+`$ and $`4^+`$ states for three calcium isotopes. In all these calculations the functional parameter set of Ref.  have been used. These results are presented in Table 2 in comparison with the calculations for $`2^+`$ and $`3^{}`$ multipolarities. One can see a here significant difference between $`2^+`$ and $`3^{}`$ states on one hand and $`0^+`$ and $`4^+`$ states on the other hand. Therefore it is sufficient to include only strong phonons which contribute significantly to the mean square radii. This result is very similar to that of Barranco and Broglia  which they obtained with the approach of Esbensen and Bertsch  for the calculation of $`\delta r_{\mathrm{ch}}^2`$ values. However, it should be mentioned that they used a non self consistent model with more freedom in the choice of parameters.

Let us now look at the results of our calculation for the calcium, tin and lead isotope chains which are presented in Figs. 18, 19, and 20, respectively. One can see here that for calcium isotopes $`3^{}`$ states give the main contribution to $`\delta r_{\mathrm{ch}}^2_{\mathrm{tot}}`$, especially for 40Ca (Fig. 18). The $`2^+`$ contribution is not very big, in contrast to the result of . This is connected with the fact that we could not obtain the experimental values for excitation energies and especially for $`B(\mathrm{E2})`$ values in calcium isotopes in the middle of the chain. Presumably for good reproduction of the experimental values one has to go beyond the QRPA in a model which allows to take into account more complex configurations than 1p1h. Nevertheless, it is shown that there is considerable influence of ground state correlations on the calcium mean square charge radii, at least in second order in the phonon amplitudes. Combining the EDF results shown in Fig. 18 with those in Fig. 16 one sees that the staggering in charge radii of calcium isotopes could be enhanced by the calculated ground state correlation corrections; however, this would produce too big negative isotope shift between 40Ca and 48Ca thus deteriorating the good agreement with experiment seen in Fig. 16 alone.

In Figs. 19 and 20 the ground state correlation contribution to $`\delta r_{\mathrm{ch}}^2`$ is shown for a few isotopes of the tin and lead chains. These nuclei are more stiff, and QRPA works much better (see Table 1). The results show that here the ground state correlation contribution is only a small fraction of the value obtained without correlations. Moreover, the general trends look roughly like the trends of the curves for $`\delta r_{\mathrm{ch}}^2`$ from Figs. 5, 6 and 11, which therefore are enhanced but not qualitatively changed by the correlations.

## 7 Summary and conclusions

In this work the energy density functional method for superfluid nuclei has been described in detail. The generalized variational principle has been formulated for the local functional with a contact effective force $`^\xi `$ in the particle–particle channel and cutoff energy $`ϵ_\mathrm{c}ϵ_\mathrm{F}`$. For spherical finite systems, the coordinate-space technique, involving the integration in the complex energy plane of the Green’s functions obtained by solving the Gor’kov equations exactly, with physical boundary conditions both for bound and scattering quasiparticle states has been used. The technique has been extended to odd nuclei by using a uniform filling approximation.

Within the framework of this EDF approach, the combined analysis of differential observables including odd-even mass differences and odd-even effects in charge radii was performed for a few isotopic chains. The staggering and kinks in $`\delta r_{\mathrm{ch}}^2`$ were shown to be very sensitive to the density dependence of the effective pairing force. The observed isotopic shifts were successfully described with the force containing a density-gradient pairing term. With such a parametrization, as the external attraction $`f_{\mathrm{ex}}^\xi `$ increases, the pairing field $`\mathrm{\Delta }(r)`$ becomes smaller in the volume and concentrates in the outer part of the nuclear surface.

A few possible parameter sets extracted from the lead chain were used to study the ground state properties of uniform nuclear matter with $`s`$-wave pairing, in particular the behavior of the energy gap at the Fermi level as a function of density. The deduced sets, which give a satisfactory description of both the neutron separation energies $`S_\mathrm{n}`$ and isotope shifts $`\delta r_{\mathrm{ch}}^2`$ in charge radii, yielded about the same value of $`\mathrm{\Delta }3.3`$ MeV at $`k_\mathrm{F}=1.15`$ fm-1 (at $`0.65`$ of the equilibrium density). It was found that some sets in the dilute limit correspond to the strong coupling regime. The situation when the regime changes from weak to strong pairing was considered in detail, and, in strong coupling, the properties of dilute matter were shown to be completely determined in leading order by the singlet scattering length $`a_{\mathrm{nn}}`$. In weak coupling, near gap closure, the pairing gap $`\mathrm{\Delta }_\mathrm{F}`$ is shown to be directly expressed through the Fermi level phase shifts. The density-dependent cutoff contact pairing interaction corresponding to the realistic value of the singlet NN scattering length was found to be a preferable choice between the deduced parameter sets to be used in the EDF calculations.

The isotopic chain of tin isotopes from the proton to the neutron drip line was calculated, and the evolution of the normal and anomalous (pairing) density distributions with $`A`$ was analyzed. The parametrization DF3 of the normal part of the functional and the density-dependent pairing interaction provided fairly good description of odd-even effects in $`S_\mathrm{n}`$ and $`\delta r_{\mathrm{ch}}^2`$; however, in the region of lighter tin isotopes the slope in $`\delta r_{\mathrm{ch}}^2`$ plotted as a function of the mass number $`A`$ come out steeper than in experiment. This might indicate that dependence on isovector density $`\rho _{}=\rho _\mathrm{n}\rho _\mathrm{p}`$ should be incorporated in the effective pairing interaction. Anomalous behavior of charge radii in Ca isotopes was reproduced in excellent agreement with experiment within the EDF approach by scaling the pairing effective interaction which was extracted from the lead chain by a factor of 1.35.

It was shown with the self-consistent sum rule approach that the density-dependent pairing induces sizeable staggering and kinks in the evolution of the mean energies of multipole excitations along isotopic chains. These effects were found to be in anticorrelation with the behavior of the ground state radii.

Using quasiparticle-RPA multipole strength distributions obtained with the self-consistent interaction, the contribution of phonons (ground state correlations) to the differences of charge radii was calculated by the method proposed by Esbensen and Bertsch . It was shown that these phonon corrections, at least for tin and lead isotopes, are quite small in comparison with the values provided by the main HFB contribution, and do not lead to qualitatively new effects.

The results obtained in the present paper were compared with the conventional Skyrme–HFB, Gogny–HFB and relativistic mean field Hartree-BCS calculations. On the whole, a better description of the differential observables was achieved with the EDF method, the most important effect which was successfully reproduced is the staggering in charge radii. Up to now, this effect was practically missed in all other mean-field models. The famous kink in the isotope shifts of Pb isotopes has been also well described by our EDF calculations, with density-dependent gradient pairing, and the crucial point is that the neutron separation energies are reproduced as well. These results support the conclusion drawn by Reinhard and Flocard  that this is apparently not the case in the RMF model. These authors showed that the description of the kink in Pb chain is not proprietary to the relativistic model, it can be also produced with the Skyrme functionals by generalizing the spin-orbit term. However, even with their “best” set SkI4, Reinhard and Flocard have found a large gap in the single-particle neutron spectrum similar to that of the RMF model. Our conclusion is that the physics behind the anomaly in the Pb radii is not solely due to the spin-orbit (which, moreover, has nothing to do with the staggering!), and that the simultaneous reproduction of both the energetic and geometrical observables is important to settle the problem.

Unfortunately, a universal parametrization of the pairing force is still lacking, and for various isotope chains, especially for lighter nuclei, the parameters of $`^\xi `$ had to be chosen somewhat different. It seems to be necessary to use a more refined parametrization, e.g. with additional dependence of $`^\xi `$ on $`\rho _{}`$. Also, more attention should be payed to the normal part of the density functional in order to choose more carefully its parameters and improve the description of the bulk properties of nuclei such as their masses and absolute values of radii . These questions will be investigated in future work.

We thank J.-F. Berger, J. Dobaczewski, V.V. Khodel and E. Saperstein for useful discussions. Partial support of this work by the Deutsche Forschungsgemeinschaft and by the Russian Foundation for Basic Research (project 98-02-16979) is gratefully acknowledged.

## Appendix A. General variational principle and ground state energy in local approximation

In this appendix we show that, in the case of weak pairing $`|\mathrm{\Delta }|ϵ_\mathrm{F}`$, the general variational principle for a local density functional with an energy cutoff $`ϵ_\mathrm{c}>ϵ_\mathrm{F}`$ can be used to describe the ground state properties of superfluid systems. This allows to perform a full HFB-like calculation with a renormalized gap equation and an effective pairing $`\delta `$-force.

We start with the usual definitions and equations which follow from the general variational principle. The total energy $`E`$ of a superfluid system, as given by eqs. (2.10)–(2.12), is a functional of the generalized one-body density matrix $`\widehat{R}`$. This matrix can be expressed in terms of the ground state expectation values of the pair products of the Landau-Migdal quasiparticle creation and annihilation operators $`\psi ^{},\psi `$:

$`\widehat{R}(1,2)`$ $`=`$ $`\left(\begin{array}{cc}\psi ^{}(2)\psi (1)& \psi (2)\psi (1)\\ \psi ^{}(2)\psi ^{}(1)& \psi (2)\psi ^{}(1)\end{array}\right)`$ (A.1)
$`=`$ $`\left(\begin{array}{cc}\rho (1,2)& \nu (1,2)\\ \nu ^{}(1,2)& \delta (1,2)\rho ^{}(1,2)\end{array}\right).`$
Here the coordinate space representation is used. The numbers in brackets stand for the set of all relevant arguments and indices, e.g. $`(1)=(\{\stackrel{}{r},s_z,\tau _3\}_1)`$. In the following, $`1`$ means integration over the continuous and summation over the discrete variables of this set.

The Bogolyubov quasiparticle creation and annihilation operators $`\beta _\alpha ^{}`$, $`\beta _\alpha `$ are defined by the transformation

$$\left(\begin{array}{c}\beta _\alpha \\ \beta _\alpha ^{}\end{array}\right)=1\left(\begin{array}{cc}U_\alpha ^{}(1)& V_\alpha ^{}(1)\\ V_\alpha (1)& U_\alpha (1)\end{array}\right)\left(\begin{array}{c}\psi (1)\\ \psi ^{}(1)\end{array}\right).$$
(A.2)
This can be inverted to express the Landau-Migdal quasiparticle operators in terms of those of Bogolyubov

$$\left(\begin{array}{c}\psi (1)\\ \psi ^{}(1)\end{array}\right)=\underset{\alpha }{}\left(\begin{array}{cc}U_\alpha (1)& V_\alpha ^{}(1)\\ V_\alpha (1)& U_\alpha ^{}(1)\end{array}\right)\left(\begin{array}{c}\beta _\alpha \\ \beta _\alpha ^{}\end{array}\right).$$
(A.3)
The generalized Bogolyubov quasiparticle density matrix $`\widehat{Q}`$ is given by

$$\widehat{Q}_{\alpha \alpha ^{}}=\left(\begin{array}{cc}\beta _\alpha ^{}^{}\beta _\alpha & \beta _\alpha ^{}\beta _\alpha \\ \beta _\alpha ^{}^{}\beta _\alpha ^{}& \beta _\alpha ^{}\beta _\alpha ^{}\end{array}\right).$$
(A.4)
By demanding now the ground sate $`|\mathrm{\Phi }_0`$ to be the Bogolyubov quasiparticle vacuum, $`\beta _\alpha |\mathrm{\Phi }_0=0`$, the matrix $`\widehat{Q}`$ becomes diagonal in the $`\alpha `$–representation:

$$\widehat{Q}_{\alpha \alpha ^{}}=\delta _{\alpha \alpha ^{}}Q,\mathrm{with}Q=\left(\begin{array}{cc}0& 0\\ 0& 1\end{array}\right),$$
(A.5)
and for (A.1) one obtains a “supermatrix” formula:

$$\widehat{R}(1,2)=\underset{\alpha \alpha ^{}}{}𝒲_\alpha (1)\widehat{Q}_{\alpha \alpha ^{}}𝒲_\alpha ^{}^{}(2)=\underset{\alpha }{}𝒲_\alpha (1)Q𝒲_\alpha ^{}(2),$$
(A.6)
with the matrix $`𝒲`$ of the Bogolyubov transformation (A.3):

$$𝒲_\alpha (1)=\left(\begin{array}{cc}U_\alpha (1)& V_\alpha ^{}(1)\\ V_\alpha (1)& U_\alpha ^{}(1)\end{array}\right).$$
(A.7)
For the generalized density matrix we thus find

$$\widehat{R}(1,2)=\underset{\alpha }{}\left(\begin{array}{cc}V_\alpha ^{}(1)V_\alpha (2)& V_\alpha ^{}(1)U_\alpha (2)\\ U_\alpha ^{}(1)V_\alpha (2)& U_\alpha ^{}(1)U_\alpha (2)\end{array}\right).$$
(A.8)
The matrix $`𝒲`$ obeys the generalized closure and orthogonality relations:

$$\underset{\alpha }{}𝒲_\alpha (1)𝒲_\alpha ^{}(2)=\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)\delta (1,2)\widehat{I}\delta (1,2),$$
(A.9)
$$1𝒲_\alpha ^{}(1)𝒲_\alpha ^{}(1)=\widehat{I}\delta _{\alpha \alpha }.$$
(A.10)
The generalized density matrix $`\widehat{R}`$ is Hermitian and from eqs. (A.5A.10) follows its important property (2.15) which, in the coordinate space representation, reads

$$3\widehat{R}(1,3)\widehat{R}(3,2)=\widehat{R}(1,2).$$
(A.11)
This property reflects the fact that the ground state is a HFB quasiparticle vacuum.

Minimizing the energy and imposing the average particle number conservation (2.14) and the relation (A.11) as constraints, we arrive at the variational principle $`\delta I[\widehat{R}]=0`$ with the variational functional of eq. (2.16) which we rewrite here as

$$I[\widehat{R}]=E[\widehat{R}]\mu \mathrm{Tr}\widehat{\rho }\mathrm{Tr}(\widehat{\mathrm{\Lambda }}(\widehat{R}\widehat{R}^2)).$$
(A.12)
The trace of a supermatrix $`\widehat{A}`$, in the coordinate space representation, is defined by

$$\mathrm{Tr}\widehat{A}=\mathrm{Tr}\left(\begin{array}{cc}\widehat{A}^{11}& \widehat{A}^{12}\\ \widehat{A}^{21}& \widehat{A}^{22}\end{array}\right)=1\underset{i}{}A^{ii}(1,1).$$
(A.13)
The total energy of the system, eqs. (2.10)–(2.12), is the functional

$`E[\widehat{R}]`$ $`=`$ $`\mathrm{Tr}(t\widehat{\rho })+E_{\mathrm{int}}[\widehat{R}]=\mathrm{Tr}(t\widehat{\rho })+E_{\mathrm{int}(\mathrm{normal})}[\widehat{\rho }]+E_{\mathrm{anomal}}[\widehat{R}]`$ (A.14)
$`=`$ $`E_{\mathrm{normal}}[\widehat{\rho }]+E_{\mathrm{anomal}}[\widehat{R}],`$
where $`t`$ is the free kinetic energy operator and $`E_{\mathrm{int}}[\widehat{R}]`$ is the interaction energy containing both normal and anomalous terms. Because of condition (A.10), the matrix $`\mathrm{\Lambda }_{\alpha \alpha ^{}}`$, in the Bogolyubov quasiparticle space, may be taken to be diagonal:

$$\mathrm{\Lambda }_{\alpha \alpha ^{}}=\left(\begin{array}{cc}E_\alpha & 0\\ 0& E_\alpha \end{array}\right)\delta _{\alpha \alpha ^{}}\widehat{I}E_\alpha \delta _{\alpha \alpha ^{}},$$
(A.15)
where $`E_\alpha `$ is the set of Lagrange multipliers. In coordinate space representation one gets

$$\widehat{\mathrm{\Lambda }}(1,2)=\underset{\alpha }{}E_\alpha 𝒲_\alpha (1)𝒲_\alpha ^{}(2).$$
(A.16)
The variational equation reads

$$\delta \left(E[\rho ,\nu ]\mu \mathrm{Tr}\rho \mathrm{Tr}(\widehat{\mathrm{\Lambda }}(\widehat{R}\widehat{R}^2))\right)=0.$$
(A.17)
For the variation of the first two terms in this equation we have

$$\delta (E[\rho ,\nu ]\mu \mathrm{Tr}\rho )=\mathrm{Tr}(\widehat{}\delta \widehat{R})=12\widehat{}(1,2)\delta \widehat{R}(2,1),$$
(A.18)
where the effective quasiparticle Hamiltonian $`\widehat{}`$ is given by

$$\widehat{}(1,2)=\left(\begin{array}{cc}h(1,2)\mu \delta (1,2)& \mathrm{\Delta }(1,2)\\ \mathrm{\Delta }^{}(1,2)& \mu \delta (1,2)h^{}(1,2)\end{array}\right),$$
(A.19)
with

$$h(1,2)=\frac{\delta E[\rho ,\nu ]}{\delta \rho (2,1)},\mathrm{\Delta }(1,2)=\frac{\delta E[\rho ,\nu ]}{\delta \nu ^{}(2,1)}=\frac{\delta E[\rho ,\nu ]}{\delta \nu ^{}(1,2)}.$$
(A.20)
The variation of the generalized density matrix,

$$\delta \widehat{R}(1,2)=\left(\begin{array}{cc}\delta \rho (1,2)& \delta \nu (1,2)\\ \delta \nu ^{}(1,2)& \delta \rho ^{}(1,2)\end{array}\right),$$
(A.21)
may be written in the form

$$\delta \widehat{R}(1,2)=\underset{\alpha \alpha ^{}}{}𝒲_\alpha (1)\delta \widehat{Q}_{\alpha \alpha ^{}}𝒲_\alpha ^{}^{}(2),$$
(A.22)
where $`\delta \widehat{Q}_{\alpha \alpha ^{}}`$ is the variation of the Bogolyubov quasiparticle density matrix (A.4). The choice of the matrix elements $`\delta \widehat{Q}_{\alpha \alpha ^{}}`$ as independent variables leads to the HFB equations in a simple way. The variation of (A.11) gives

$$\delta (\widehat{R}(1,2)\widehat{R}^2(1,2))=\underset{\alpha \alpha ^{}}{}𝒲_\alpha (1)\left(\begin{array}{cc}1& 0\\ 0& 1\end{array}\right)\delta \widehat{Q}_{\alpha \alpha ^{}}𝒲_\alpha ^{}^{}(2).$$
(A.23)
Eq.(A.18) may thus be written as

$`\mathrm{Tr}(\widehat{}\delta \widehat{R})`$ $`=`$ $`{\displaystyle \underset{\alpha \alpha ^{}}{}}\mathrm{Tr}𝒲_\alpha ^{}^{}\widehat{}𝒲_\alpha \delta \widehat{Q}_{\alpha \alpha ^{}}`$ (A.24)
$`=`$ $`{\displaystyle \underset{\alpha \alpha ^{}}{}}{\displaystyle 12\underset{abcd}{}𝒲_{}^{}{}_{\alpha ^{}}{}^{ab}(1)\widehat{}^{bc}(1,2)𝒲_\alpha ^{cd}(2)\delta \widehat{Q}_{\alpha \alpha ^{}}^{da}},`$
and the variation of the last term in (A.17) may be expressed as

$$\delta \mathrm{Tr}\widehat{\mathrm{\Lambda }}(\widehat{R}\widehat{R}^2))=_\alpha \mathrm{Tr}(\begin{array}{cc}E_\alpha & 0\\ 0& E_\alpha \end{array})\delta \widehat{Q}_{\alpha \alpha }.$$
(A.25)
Thus, the variation equation (A.17) becomes

$$\underset{\alpha \alpha ^{}}{}\mathrm{Tr}\left[𝒲_\alpha ^{}^{}\widehat{}𝒲_\alpha \delta _{\alpha ^{}\alpha }\left(\begin{array}{cc}E_\alpha & 0\\ 0& E_\alpha \end{array}\right)\right]\delta \widehat{Q}_{\alpha \alpha ^{}}=0.$$
(A.26)
Since the elements $`\delta \widehat{Q}_{\alpha \alpha ^{}}`$ are independent, we obtain the matrix equation

$$\widehat{}𝒲_\alpha =𝒲_\alpha \left(\begin{array}{cc}E_\alpha & 0\\ 0& E_\alpha \end{array}\right),$$
(A.27)
which is equivalent to the HFB equations

$$2\widehat{}(1,2)\left(\begin{array}{c}U_\alpha (2)\\ V_\alpha (2)\end{array}\right)=E_\alpha \left(\begin{array}{c}U_\alpha (1)\\ V_\alpha (1)\end{array}\right),$$
(A.28)
$$2\widehat{}(1,2)\left(\begin{array}{c}V_\alpha ^{}(2)\\ U_\alpha ^{}(2)\end{array}\right)=E_\alpha \left(\begin{array}{c}V_\alpha ^{}(1)\\ U_\alpha ^{}(1)\end{array}\right).$$
(A.29)
The matrix of second variational derivatives of the energy functional is the supermatrix of the effective quasiparticle interaction amplitudes $`\widehat{}`$. In the pp–channel one finds

$$^{\mathrm{pp}}(1,2;3,4)=\frac{\delta \mathrm{\Delta }(1,2)}{\delta \nu (3,4)}.$$
(A.30)
We assume that $`^{\mathrm{pp}}`$ is a functional of the normal quasiparticle density matrix $`\widehat{\rho }`$ and that the term $`E_{\mathrm{anomal}}`$ in (A.14), which depends explicitly on the anomalous density matrix, may be written as

$$E_{\mathrm{anomal}}[\widehat{R}]=\frac{1}{4}14\nu ^{}(2,1)_\mathrm{a}^{\mathrm{pp}}(1,2;3,4;[\widehat{\rho }])\nu (3,4),$$
(A.31)
with the antisymmetrized pp–interaction

$$_\mathrm{a}^{\mathrm{pp}}(1,2;3,4)=^{\mathrm{pp}}(1,2;3,4)^{\mathrm{pp}}(1,2;4,3).$$
(A.32)
The gap equation reads

$$\mathrm{\Delta }(1,2)=\frac{1}{2}34_\mathrm{a}^{\mathrm{pp}}(1,2;3,4;[\widehat{\rho }])\nu (3,4).$$
(A.33)
The expression for the anomalous energy (A.31) may thus be written in the form

$$E_{\mathrm{anomal}}=\frac{1}{2}12\nu ^{}(2,1)\mathrm{\Delta }(1,2).$$
(A.34)
To proceed further and perform the renormalization procedure it is convenient to use the Green’s function formalism (see Appendix B). From (A.8) and eqs. (B.9)–(B.11) it follows that the generalized density matrix $`\widehat{R}`$ may be represented by a contour integral of the generalized Green’s function $`\widehat{G}`$ in the complex energy plane:

$$\widehat{R}(1,2)=_C\frac{ϵ}{2\pi \mathrm{i}}\widehat{G}(1,2;ϵ).$$
(A.35)
The integral is performed along a contour $`C`$ which goes from $`\mathrm{}`$ to 0 below the real $`ϵ`$-axis (Im$`ϵ=0`$), crosses this axis at $`ϵ=0`$ and goes back to $`\mathrm{}`$ above it ($`\mathrm{Im}ϵ=+0`$), with the energy variable $`ϵ`$ measured from the chemical potential $`\mu `$. The contour is shown schematically in Fig. 21, see Appendix C.

The normal and anomalous density matrices are then given by

$$\rho (1,2)=_C\frac{ϵ}{2\pi \mathrm{i}}G_s(1,2;ϵ)=\underset{(E_\alpha >0)}{}V_\alpha ^{}(1)V_\alpha (2),$$
(A.36)
$$\nu (1,2)=_C\frac{ϵ}{2\pi \mathrm{i}}F(1,2;ϵ)=\underset{(E_\alpha >0)}{}V_\alpha ^{}(1)U_\alpha (2).$$
(A.37)
The summation (or integration) in these expressions extends, in principle, to infinity, the convergence being dependent on the properties of the effective interaction. It is well known that in the case of a local force the anomalous density matrix $`\nu (1,2)`$ diverges at $`\stackrel{}{r}_1\stackrel{}{r}_2`$. In the gap equation (A.33) with $`\widehat{\nu }`$ expressed by (A.37), the integration over the region far from the Fermi surface can be avoided by using the standard renormalization procedure with an arbitrary energy cutoff . But to calculate the energy of the system, in particular its (negative) anomalous part (A.34), one needs to know the “untruncated” density matrices. On the other hand, as it is also well known, the total pairing energy $`E_{\mathrm{pair}}`$ contains a (positive) contribution from the kinetic energy term so that the sum of the two contributions converges more rapidly and $`E_{\mathrm{pair}}`$ gets a finite value even in the case of a local pairing field (at least in infinite matter, see Appendix B). We shall show that the energy density of the uniform infinite system may be calculated exactly in the leading order in $`\mathrm{\Delta }^2(p_\mathrm{F})/ϵ_\mathrm{F}`$ by using a renormalized gap equation with an arbitrary cutoff $`ϵ_\mathrm{c}`$ but such that $`ϵ_\mathrm{c}>ϵ_\mathrm{F}`$, and by applying the variational principle $`\delta I[\widehat{R}]=0`$ to the functional $`I[\widehat{R}]`$ given by (A.12) but considered now as a functional of the corresponding cutoff generalized density matrix $`\widehat{R}_\mathrm{c}`$.

We introduce an energy cutoff $`ϵ_\mathrm{c}>ϵ_\mathrm{F}`$. The contour $`C`$ is split into two parts $`_C=_{C_\mathrm{c}}+_{C_B},`$ where $`C_\mathrm{c}`$ lies below and above the real $`ϵ`$-axis for $`ϵ_\mathrm{c}<\mathrm{Re}ϵ<0`$ and the contour $`C_B`$ lies below and above it for $`\mathrm{Re}ϵ<ϵ_\mathrm{c}`$ . Then we have

$$\nu (1,2)=\nu _\mathrm{c}(1,2)+\delta _\mathrm{c}\nu (1,2),$$
(A.38)
where

$$\nu _\mathrm{c}(1,2)=\underset{0<E_\alpha <ϵ_\mathrm{c}}{}V_\alpha ^{}(1)U_\alpha (2)=_{C_\mathrm{c}}\frac{ϵ}{2\pi \mathrm{i}}F(1,2;ϵ),$$
(A.39)
and

$$\delta _\mathrm{c}\nu (1,2)=\underset{E_\alpha >ϵ_\mathrm{c}}{}V_\alpha ^{}(1)U_\alpha (2)=_{C_B}\frac{ϵ}{2\pi \mathrm{i}}F(1,2;ϵ).$$
(A.40)
The renormalized gap equation may be written in the form

$$\mathrm{\Delta }(1,2)=\frac{1}{2}34_a^\xi (1,2;3,4)\nu _\mathrm{c}(3,4).$$
(A.41)
where the effective scattering amplitude $`^\xi `$ is to be found from the equation

$`^\xi (1,2;3,4)`$ $`=`$ $`^{\mathrm{pp}}(1,2;3,4)`$ (A.42)
$`+`$ $`{\displaystyle }58^{\mathrm{pp}}(1,2;5,6)B(5,6;7,8)^\xi (7,8;3,4),`$
with

$$B(1,2;3,4)=_{C_B}\frac{ϵ}{2\pi \mathrm{i}}G_0(1,3;ϵ)\overline{G}_s(4,2;ϵ).$$
(A.43)
Thus the anomalous energy (A.34) may be represented by a sum of two terms:

$`E_{\mathrm{anomal}}`$ $`=`$ $`E_{\mathrm{anomal}}^\mathrm{c}+\delta _\mathrm{c}E_{\mathrm{anomal}}`$ (A.44)
$``$ $`{\displaystyle \frac{1}{2}}{\displaystyle 12_{C_\mathrm{c}}\frac{ϵ}{2\pi \mathrm{i}}F^{}(2,1;ϵ)\mathrm{\Delta }(1,2)}`$
$`+`$ $`{\displaystyle \frac{1}{2}}{\displaystyle 12_{C_B}\frac{ϵ}{2\pi \mathrm{i}}F^{}(2,1;ϵ)\mathrm{\Delta }(1,2)}.`$
The normal density matrix is written in a similar way as a sum

$$\rho (1,2)=\rho _\mathrm{c}(1,2)+\delta _\mathrm{c}\rho (1,2),$$
(A.45)
where

$$\rho _\mathrm{c}(1,2)=\underset{0<E_\alpha <ϵ_\mathrm{c}}{}V_\alpha ^{}(1)V_\alpha (2)=_{C_\mathrm{c}}\frac{ϵ}{2\pi \mathrm{i}}G_s(1,2;ϵ),$$
(A.46)
and

$$\delta _\mathrm{c}\rho (1,2)=\underset{E_\alpha >ϵ_\mathrm{c}}{}V_\alpha ^{}(1)V_\alpha (2)=_{C_B}\frac{ϵ}{2\pi \mathrm{i}}G_s(1,2;ϵ).$$
(A.47)
Using the definition of the single particle Hamiltonian $`h`$, eq. (A.20), the amount of the total energy related to $`\delta _\mathrm{c}\rho `$ in first order may be found by varying the normal part of the energy functional with respect to the normal density:

$$\delta _\mathrm{c}E_{\mathrm{normal}}=12h(1,2)\delta _\mathrm{c}\rho (2,1).$$
(A.48)
This includes also, if $`^{\mathrm{pp}}`$ depends on $`\rho `$, a contribution from the variation of $`E_{\mathrm{anomal}}`$.

With $`\delta _\mathrm{c}E_{\mathrm{anomal}}`$ from (A.44) and $`\delta _\mathrm{c}E_{\mathrm{normal}}`$ defined by (A.48) we obtain the total change of the variational functional related to the cutoff:

$$\delta _\mathrm{c}(E\mu N)=_{C_B}\frac{ϵ}{2\pi \mathrm{i}}(hG_s+\frac{1}{2}F^{}\mathrm{\Delta }\mu G_s).$$
(A.49)
This can be written in another form using the relation

$$\mathrm{\Delta }F^{}=(ϵh+\mu )G_s1$$
(A.50)
which follows from the Gor’kov equations (B.3):

$`\delta _\mathrm{c}\left(E\mu N(\mu )\right)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{Tr}{\displaystyle _{C_B}}{\displaystyle \frac{ϵ}{2\pi \mathrm{i}}}(ϵ+h\mu )G_s`$ (A.51)
$`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{Tr}{\displaystyle _{C_B}}{\displaystyle \frac{ϵ}{2\pi \mathrm{i}}}(ϵ+h\mu )G_0\mathrm{\Delta }F^{}.`$
To get the second equality in this equation we have used the relation
$`\mathrm{Tr}_{C_B}ϵ(ϵ+h\mu )G_0=0`$ which holds because the Green’s function $`G_0`$ is diagonal in the basis of eigenfunctions of $`h`$ (see eq. (B.8)) and has no singularities embraced by the contour $`C_B`$.

In lowest order in $`\mathrm{\Delta }`$, from (B.3) one has $`F^{}=\overline{G}_0\mathrm{\Delta }^{}G_0`$ and we get

$`\delta _\mathrm{c}\left(E\mu N(\mu )\right)`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{Tr}{\displaystyle _{C_B}}{\displaystyle \frac{ϵ}{2\pi \mathrm{i}}}(ϵ+h\mu )G_0\mathrm{\Delta }\overline{G}_0\mathrm{\Delta }^{}G_0`$ (A.52)
$`=`$ $`{\displaystyle \underset{ϵ_i\mu <ϵ_\mathrm{c}}{}}\left({\displaystyle \underset{ϵ_k\mu >ϵ_\mathrm{c}}{}}{\displaystyle \frac{(ϵ_iϵ_k)|\mathrm{\Delta }_{ik}|^2}{(ϵ_i+ϵ_k2\mu )^2}}\right).`$
Here $`\mathrm{\Delta }_{ik}`$ are the matrix elements of $`\mathrm{\Delta }`$ in the basis of the single-particle Hamiltonian $`h`$ with $`ϵ_i`$ being its eigenvalues (see Appendix B, eq. (B.7)).

In uniform infinite matter, the pairing field $`\mathrm{\Delta }(1,2)`$ is a function of $`|\stackrel{}{r}_1\stackrel{}{r}_2|`$, all its nondiagonal matrix elements vanish and we see that $`\delta _\mathrm{c}\left(E\mu N(\mu )\right)=0`$ in the first order upon imposing the energy cutoff $`ϵ_\mathrm{c}>ϵ_\mathrm{F}`$. In finite systems, the change of $`\delta _\mathrm{c}\left(E\mu N(\mu )\right)`$ vanishes in the first order in diagonal approximation. A nonvanishing right-hand side of eq. (A.52) arises only due to nonuniformity of $`\mathrm{\Delta }`$, which means that in finite systems with more or less uniform density in the volume (like heavy atomic nuclei), this would be a surface effect. Reliable estimates of the sum in (A.52) involving only nondiagonal matrix elements $`\mathrm{\Delta }_{ik,ik}`$ are not easy to obtain. They depend strongly on the actual distribution of the pairing field in real nuclei, in particular in the surface region. In larger systems, with presence of volume pairing, the nondiagonal matrix elements of $`\mathrm{\Delta }`$ are relatively small ($`A^{1/3}`$, or less) so that the change of $`E\mu N`$ due to the cutoff is expected to be at least by a factor of $`A^{1/3}`$ smaller compared to $`E_{\mathrm{pair}}`$, and therefore we can write

$$|\delta _\mathrm{c}E\mu \delta _\mathrm{c}N(\mu )||E_{\mathrm{pair}}|,$$
(A.53)
where

$$\delta _\mathrm{c}N(\mu )=1\delta _\mathrm{c}\rho (1,1)=N(\mu )N_\mathrm{c}(\mu ),$$
(A.54)
with $`N_\mathrm{c}`$ the number of particles corresponding to the cutoff density $`\rho _\mathrm{c}`$ of eq. (A.46):

$$N_\mathrm{c}=1\rho _\mathrm{c}(1,1).$$
(A.55)
The pairing energy $`E_{\mathrm{pair}}`$ is given in Appendix B by an exact expression (B.41) and by a BCS estimate (B.46).

The change in the system energy related with the cutoff may be written as

$$\delta _\mathrm{c}E=E[\widehat{R}]E[\widehat{R}_\mathrm{c}]=\mu \delta _\mathrm{c}N,$$
(A.56)
with $`\delta _\mathrm{c}N`$ defined in (A.54). Minimizing the energy functional with a fixed cutoff $`ϵ_\mathrm{c}`$ would give an additional change

$$\delta E_\mathrm{c}=E[\widehat{R}_\mathrm{c}]E_\mathrm{c}[\widehat{R}_\mathrm{c}]=\mu _\mathrm{c}\delta N_\mathrm{c},$$
(A.57)
where $`\delta N_\mathrm{c}`$ is the change in the particle number within the cutoff energy space $`ϵ<ϵ_\mathrm{c}`$. In first order we have

$$E[\widehat{R}]=E_\mathrm{c}[\widehat{R}_\mathrm{c}]+\mu \delta _\mathrm{c}N+\mu _\mathrm{c}\delta N_\mathrm{c}=E_\mathrm{c}[\widehat{R}_\mathrm{c}]+\mu (\delta _\mathrm{c}N+\delta N_\mathrm{c}).$$
(A.58)
Of course, the constraint $`\mathrm{Tr}\rho _\mathrm{c}=N`$ has also to be used for the cutoff functional which gives $`\delta _\mathrm{c}N+\delta N_\mathrm{c}=0`$.

We thus come to the important conclusion that passing to the cutoff functional with $`ϵ_\mathrm{c}>ϵ_\mathrm{F}`$ leaves, to the first order in the weak pairing approximation, the energy, the variational functional $`E\mu N`$ and the chemical potential $`\mu `$ of the system unchanged (to the extent that the estimate (A.53) for a given nucleus is correct).

According to the Hohenberg-Kohn theorem , the nuclear ground state properties can be described by a local density functional. One may therefore assume that the anomalous energy is a functional of the local density $`\rho `$ too. The above consideration implies that this should be also valid for the cutoff functional. Both the cutoff anomalous density and the effective pairing interaction $`^\xi `$ may then be regarded as local functionals of $`\rho `$. We introduce a local approximation to $`^\xi `$:

$`^\xi (1,2;3,4)`$ $`=`$ $`^\xi (\stackrel{}{r}_1,\tau _1)`$ (A.59)
$`\times `$ $`\delta (\stackrel{}{r}_1\stackrel{}{r}_2)\delta (\stackrel{}{r}_1\stackrel{}{r}_3)\delta (\stackrel{}{r}_2\stackrel{}{r}_4)\delta _{s_1s_3}\delta _{s_2s_4}\delta _{\tau _1\tau _3}\delta _{\tau _2\tau _4}.`$
where

$$^\xi (\stackrel{}{r},\tau )=^\xi ([\rho _\mathrm{c}(\stackrel{}{r},n),\rho _\mathrm{c}(\stackrel{}{r},p)];\tau )$$
(A.60)
is a local functional of the cutoff quasiparticle density

$$\rho _\mathrm{c}(\stackrel{}{r}_1,\tau _1)=\mathrm{Tr}_{s_1}2\rho _\mathrm{c}(1,2)\delta (2,1),\tau _1\{n,p\},$$
(A.61)
where $`\delta (1,2)=\delta (\stackrel{}{r}_1\stackrel{}{r}_2)\delta _{s_1s_2}\delta _{\tau _1\tau _2}`$.

In the case of a local functional it is convenient to work with the so-called anomalous density in the real space. For the time-reversed single-particle states $`|\overline{k}=T|k`$, in the coordinate space representation, we define

$$\varphi _{\overline{k}}(1)=2𝒯(1,2)\varphi _k^{}(2),\varphi _{\overline{k}}^{}(1)=2\varphi _k(2)𝒯^{}(2,1),$$
(A.62)
with the operator

$$𝒯(1,2)=\delta (\stackrel{}{r}_1\stackrel{}{r}_2)(1)^{\frac{1}{2}s_1}\delta _{s_1s_2}\delta _{\tau _1\tau _2}.$$
(A.63)
This operator is antisymmetric, $`𝒯(2,1)=𝒯(1,2),`$ and has the properties $`3𝒯(1,3)𝒯(3,2)=\delta (1,2)`$ and $`3𝒯(1,3)𝒯^{}(3,2)=\delta (1,2)`$. The anomalous density is then defined by

$$\nu _\mathrm{c}(\stackrel{}{r}_1,\tau _1)=\frac{1}{2}\mathrm{Tr}_{s_1}2\nu _\mathrm{c}(1,2)𝒯^{}(2,1),$$
(A.64)
and its complex conjugate by

$$\nu _\mathrm{c}^{}(\stackrel{}{r}_1,\tau _1)=\frac{1}{2}\mathrm{Tr}_{s_1}2𝒯(1,2)\nu _\mathrm{c}^{}(2,1).$$
(A.65)
In the local approximation we get

$$\mathrm{\Delta }(1,2)=𝒯(1,2)\mathrm{\Delta }(\stackrel{}{r}_2,\tau _2).$$
(A.66)
With these definitions, the anomalous part of the cutoff functional reduces to the expression (2.25), and for the local pairing field $`\mathrm{\Delta }(\stackrel{}{r},\tau )`$ one obtains the gap equation in the simple multiplicative form (2.27).

One may notice that our approach does not imply a cutoff of the basis since the general variational principle can be formulated with a “cutoff” local-density functional from which the ground state characteristics of a superfluid system may be calculated through the solutions of the Bogolyubov equations (A.28) at the stationary point. To construct the normal and anomalous densities, entering this local functional, only those solutions from the whole set are needed which correspond to the eigenenergies $`E_\alpha `$ of the HFB Hamiltonian up to the cutoff131313At first glance the situation looks similar to the HF case with a EDF without pairing where only the occupied orbitals up to $`ϵ_\mathrm{F}`$ in the self-consistent mean field are needed to calculate the density (a “natural” HF cutoff). But it should be emphasized that the Bogolyubov solutions cannot be obtained from the HF ones by perturbative methods, they spread to infinity in the energy space at any nonvanishing gap, hence the convergence problem for the densities. This problem, as usually thought of, can be avoided only if one uses a finite-range effective force ensuring the convergence of the anomalous density at some energy which depends on the range of nonlocality of the force. Then the HFB equations become nonlocal leading to a technical problem; in practice these equations are solved with approximate methods, for example, by expanding the HFB wave functions in a harmonic oscillator basis . $`ϵ_\mathrm{c}>ϵ_\mathrm{F}`$. This means that one can implement a local pairing effective force in the HFB (or Gor’kov) equations and, most important, in the gap equation, renormalized with the same energy cutoff $`ϵ_\mathrm{c}`$, in which the anomalous density by construction does not diverge. The whole set of resulting equations can be solved, for spherical nuclei, exactly with a coordinate-space technique (see Appendix C). It should nevertheless be stressed that the main difficulty now is to find such a local effective density-dependent pairing interaction $`^\xi `$ that would be an universal one and would allow to reproduce the pairing ground-state nuclear properties (at least for a large domain of heavy nuclei).

## Appendix B. Gor’kov equations and pairing energy

In this appendix we derive the general expression for the pairing energy using the Green’s function formalism. We show that, in the case of a weak local pairing field, in infinite matter this expression reduces to the BCS formula used in eq. (4.30).

We start with Gor’kov equations  for the generalized Green’s function $`\widehat{G}(ϵ)`$ written in matrix form,

$$(ϵ\widehat{})\widehat{G}(ϵ)=\widehat{I},$$
(B.1)
where $`\widehat{}`$ is the effective single-quasiparticle Hamiltonian (A.19) and

$$\widehat{G}(ϵ)=\left(\begin{array}{cc}G_\mathrm{s}(ϵ)& F(ϵ)\\ F^{}(ϵ)& \overline{G}_\mathrm{s}(ϵ)\end{array}\right)$$
(B.2)
with $`G_\mathrm{s}`$ the normal and $`F`$ the anomalous Green’s function, respectively. For $`G_\mathrm{s}`$ and $`F^{}`$, eq. (B.1) yields the set of two equations

$$G_\mathrm{s}(ϵ)=G_0(ϵ)+G_0(ϵ)\mathrm{\Delta }F^{}(ϵ),F^{}(ϵ)=\overline{G}_0(ϵ)\mathrm{\Delta }^{}G_\mathrm{s}(ϵ),$$
(B.3)
where $`G_0`$ and $`\overline{G}_0`$ are the single particle Green’s functions connected with the normal parts $`h\mu `$ and $`\mu h^{}`$ of the effective Hamiltonian, respectively:

$$G_0(ϵ)=\frac{1}{ϵh+\mu },\overline{G}_0(ϵ)=\frac{1}{ϵ+h^{}\mu }.$$
(B.4)
The Green’s functions $`G_0`$ and $`\overline{G}_0`$ are related to each other by $`\overline{G}_0(ϵ)=G_0^{}(ϵ^{}).`$

Similarly, for $`\overline{G}_\mathrm{s}`$ and $`F`$ we can write the set of two equations

$$\overline{G}_\mathrm{s}(ϵ)=\overline{G}_0(ϵ)\overline{G}_0(ϵ)\mathrm{\Delta }^{}F(ϵ),F(ϵ)=G_0(ϵ)\mathrm{\Delta }\overline{G}_\mathrm{s}(ϵ).$$
(B.5)
From the Gor’kov equations it follows that

$$\overline{G}_\mathrm{s}(ϵ)=G_\mathrm{s}^{}(ϵ^{}),F^{}(ϵ)=F^{}(ϵ^{}).$$
(B.6)
The “zero” Green’s function $`G_0`$ is diagonal in the basis of eigenfunctions of the Hamiltonian $`h`$ ($`k`$–representation),

$$2h(1,2)\varphi _k(2)=ϵ_k\varphi _k(1),$$
(B.7)
and for its spectral decomposition we get

$$G_0(1,2;ϵ)=\underset{k}{}\left[\frac{n_k}{ϵϵ_k+\mu \mathrm{i}\delta }+\frac{1n_k}{ϵϵ_k+\mu +\mathrm{i}\delta }\right]\varphi _k(1)\varphi _k^{}(2),$$
(B.8)
where $`n_k=1`$ if $`ϵ_k<\mu `$ and $`n_k=0`$ if $`ϵ_k>\mu `$.

For systems with even particle number, the spectral expansions of $`G_\mathrm{s}`$, $`F`$ and $`F^{}`$, in coordinate space representation, may be written in the form:

$$G_\mathrm{s}(1,2;ϵ)=\underset{(E_\alpha >0)}{}\left[\frac{V_\alpha ^{}(1)V_\alpha (2)}{ϵ+E_\alpha \mathrm{i}\delta }+\frac{U_\alpha (1)U_\alpha ^{}(2)}{ϵE_\alpha +\mathrm{i}\delta }\right],$$
(B.9)
$$F(1,2;ϵ)=\underset{(E_\alpha >0)}{}\left[\frac{V_\alpha ^{}(1)U_\alpha (2)}{ϵ+E_\alpha i\delta }+\frac{U_\alpha (1)V_\alpha ^{}(2)}{ϵE_\alpha +\mathrm{i}\delta }\right],$$
(B.10)
$$F^{}(1,2;ϵ)=\underset{(E_\alpha >0)}{}\left[\frac{U_\alpha ^{}(1)V_\alpha (2)}{ϵ+E_\alpha \mathrm{i}\delta }+\frac{V_\alpha (1)U_\alpha ^{}(2)}{ϵE_\alpha +\mathrm{i}\delta }\right],$$
(B.11)
with $`E_\alpha `$ the exact eigenvalues and $`U_\alpha `$, $`V_\alpha `$ the exact eigenfunctions of eqs. (A.28). The corresponding expressions in the $`k`$–representation are:

$$G_{\mathrm{s}ij}(ϵ)=\underset{(E_\alpha >0)}{}\left[\frac{V_{i\alpha }^{}V_{j\alpha }}{ϵ+E_\alpha \mathrm{i}\delta }+\frac{U_{i\alpha }U_{j\alpha }^{}}{ϵE_\alpha +\mathrm{i}\delta }\right],$$
(B.12)
$$F_{ij}(ϵ)=\underset{(E_\alpha >0)}{}\left[\frac{V_{i\alpha }^{}U_{j\alpha }}{ϵ+E_\alpha \mathrm{i}\delta }+\frac{U_{i\alpha }V_{j\alpha }^{}}{ϵE_\alpha +\mathrm{i}\delta }\right],$$
(B.13)
$$F_{ij}^{}(ϵ)=\underset{(E_\alpha >0)}{}\left[\frac{U_{i\alpha }^{}V_{j\alpha }}{ϵ+E_\alpha \mathrm{i}\delta }+\frac{V_{i\alpha }U_{j\alpha }^{}}{ϵE_\alpha +\mathrm{i}\delta }\right],$$
(B.14)
with $`U_{k\alpha }=1\varphi _k^{}(1)U_\alpha (1),`$ and $`V_{k\alpha }=1\varphi _k(1)V_\alpha (1).`$ In diagonal approximation, which is exact in uniform infinite matter, the only nonzero matrix elements of $`G_\mathrm{s}`$, $`F`$ and $`F^{}`$ are:

$$G_{\mathrm{s}kk}(ϵ)=\frac{|v_k|^2}{ϵ+E_k\mathrm{i}\delta }+\frac{|u_k|^2}{ϵE_k+\mathrm{i}\delta },$$
(B.15)
$$F_{k\overline{k}}(ϵ)=F_{\overline{k}k}(ϵ)=\frac{\mathrm{\Delta }_k}{2E_k}\left(\frac{1}{ϵ+E_k\mathrm{i}\delta }\frac{1}{ϵE_k+\mathrm{i}\delta }\right),$$
(B.16)
$$F_{k\overline{k}}^{}(ϵ)=F_{\overline{k}k}^{}(ϵ)=\frac{\mathrm{\Delta }_k^{}}{2E_k}\left(\frac{1}{ϵ+E_k\mathrm{i}\delta }\frac{1}{ϵE_k+\mathrm{i}\delta }\right),$$
(B.17)
with $`|\overline{k}`$ being the time reversed of $`|k`$ as defined by (A.62). Here we have used the customary conventions:

$$\mathrm{\Delta }_k=2E_kv_k^{}u_k,E_k=\sqrt{(ϵ_k\mu )^2+|\mathrm{\Delta }_k|^2},$$

$$|u_k|^2=\frac{1}{2}(1+\frac{ϵ_k\mu }{E_k}),|v_k|^2=\frac{1}{2}(1\frac{ϵ_k\mu }{E_k}).$$

The corresponding spectral expansions for $`\overline{G}_\mathrm{s}`$ may be found from those for $`G_\mathrm{s}`$ by using the relations (B.6). It is easily seen that the expansions written above for the anomalous Green’s functions $`F`$ and $`F^{}`$ may be also obtained from each other by using (B.6).

In what follows the Green’s function method is used to extract the pairing contribution to the total energy of the system. To isolate the pure pairing part one has to consider two effects: (i) the direct influence of the pairing gap $`\mathrm{\Delta }`$ at fixed single-particle Hamiltonian $`h\mu `$ and (ii) the “polarization” mechanism due to variation of $`h`$ and $`\mu `$ in the presence of the pairing field.

To obtain the first change in the density, when $`\mathrm{\Delta }0`$ but $`h`$ and $`\mu `$ fixed, we introduce a “zero” density matrix:

$$\rho _0(1,2)=_C\frac{ϵ}{2\pi \mathrm{i}}G_0(1,2;ϵ).$$
(B.18)
For the difference $`\delta \rho =\rho \rho _0`$, eq. (B.3) yields

$$\delta \rho =_C\frac{ϵ}{2\pi \mathrm{i}}(G_\mathrm{s}(ϵ)G_0(ϵ))=_C\frac{ϵ}{2\pi \mathrm{i}}G_0(ϵ)\mathrm{\Delta }F^{}(ϵ).$$
(B.19)
Using (B.8) and (B.14) we obtain in the $`k`$–representation:

$$\delta \rho _{ij}=\underset{E_\alpha >0}{}\underset{l}{}\frac{\mathrm{\Delta }_{il}}{E_\alpha +|ϵ_i\mu |}[n_iV_{l\alpha }U_{j\alpha }^{}+(1n_i)U_{l\alpha }^{}V_{j\alpha }],$$
(B.20)
or, in diagonal approximation,

$$\delta \rho (\stackrel{}{r},\tau )=\underset{k}{}\underset{s_z}{}(12n_k)\frac{|\mathrm{\Delta }_k|^2}{2E_k(E_k+|ϵ_k\mu |)}|\varphi _k(\stackrel{}{r},s_z,\tau )|^2.$$
(B.21)
One should keep in mind that $`\rho `$ and $`\rho _0`$ belong to states with different particle number:

$$\delta N(\mu )=N(\mu )N_0(\mu )0.$$
(B.22)
The second change in the density we get in the situation without any pairing. Then we have the HF vacuum and the corresponding HF Green’s function,

$$G_{\mathrm{HF}}(ϵ)=\frac{1}{ϵh_{\mathrm{HF}}+\mu _{\mathrm{HF}}},$$
(B.23)
with $`h_{\mathrm{HF}}`$ being the self-consistent HF single-particle Hamiltonian. The HF density matrix is given by

$$\rho _{\mathrm{HF}}(1,2)=_C\frac{ϵ}{2\pi \mathrm{i}}G_{\mathrm{HF}}(1,2;ϵ)$$
(B.24)
with the average number of particles $`\mathrm{HF}|\widehat{N}|\mathrm{HF}=N_{\mathrm{HF}}(\mu _{\mathrm{HF}})`$, i.e. $`\mathrm{Tr}\rho _{\mathrm{HF}}=N`$.

Collecting contributions from both pairing-induced effects we get the total variation of the density matrix:

$$\delta \rho _{\mathrm{pair}}=\rho \rho _{\mathrm{HF}}=\delta \rho +\delta \rho _0,$$
(B.25)
with $`\delta \rho _0=\rho _0\rho _{\mathrm{HF}}`$. From the definition of the Green’s function we find

$$G_{\mathrm{HF}}^1G_0^1=h\mu h_{\mathrm{HF}}+\mu _{\mathrm{HF}},$$
(B.26)
or

$$G_0G_{\mathrm{HF}}=G_{\mathrm{HF}}(\delta U_{\mathrm{pair}}\delta \mu )G_0,$$
(B.27)
where $`\delta U_{\mathrm{pair}}=hh_{\mathrm{HF}}`$ and $`\delta \mu =\mu \mu _{\mathrm{HF}}`$.

The expression for $`\delta \rho _0`$ can then be written as

$`\delta \rho _0`$ $`=`$ $`{\displaystyle _C}{\displaystyle \frac{ϵ}{2\pi \mathrm{i}}}(G_0(ϵ)G_{\mathrm{HF}}(ϵ))`$ (B.28)
$`=`$ $`{\displaystyle _C}{\displaystyle \frac{ϵ}{2\pi \mathrm{i}}}G_{\mathrm{HF}}(ϵ)(\delta U_{\mathrm{pair}}\delta \mu )G_0(ϵ),`$
Dealing with the first-order terms of perturbation theory, the integral over $`ϵ`$ may be replaced by the static particle–hole propagator $`𝒜`$ (see Ref.):

$$_C\frac{ϵ}{2\pi \mathrm{i}}G_{\mathrm{HF}}(ϵ)G_0(ϵ)_C\frac{ϵ}{2\pi \mathrm{i}}G_{\mathrm{HF}}(ϵ)G_{\mathrm{HF}}(ϵ)𝒜.$$
(B.29)
Its convolution with any perturbation operator which is diagonal in the HF basis, i.e. commutes with $`h_{\mathrm{HF}}`$, vanishes (no polarization effect). This property will be used to obtain the energy variation connected with $`\delta \rho _0`$ (see eqs. (B.38) and (B.39) below). The average number of particles of each kind is fixed: $`\mathrm{HFB}|\widehat{N}|\mathrm{HFB}=\mathrm{HF}|\widehat{N}|\mathrm{HF}`$. The total density variation due to pairing, of course, does not change the particle number, so

$$\delta N_{\mathrm{pair}}=\delta N(\mu )+\delta N_0(\mu ,\mu _{\mathrm{HF}})=0.$$
(B.30)
By definition, the pairing energy is the total change of the system energy due to pairing correlations:

$$E_{\mathrm{pair}}=E_{\mathrm{normal}}[\rho ]+E_{\mathrm{anomal}}[\rho ,\nu ]E_{\mathrm{normal}}[\rho _{\mathrm{HF}}].$$
(B.31)
Analogously to the procedure applied above to the density matrix, we can split the pairing energy into two parts

$$E_{\mathrm{pair}}=\delta E+\delta E_0,$$
(B.32)
where the difference

$$\delta E=E_{\mathrm{normal}}[\rho ]+E_{\mathrm{anomal}}[\rho ,\nu ]E_{\mathrm{normal}}[\rho _0]$$

is the “direct” pairing contribution and

$$\delta E_0=E_{\mathrm{normal}}[\rho _0]E_{\mathrm{normal}}[\rho _{\mathrm{HF}}]$$

is the polarization pairing energy.

In the first step we calculate $`\delta E`$. The anomalous part of the energy (A.34) may be written in the form:

$$E_{\mathrm{anomal}}[\rho ,\nu ]=\frac{1}{2}\mathrm{Tr}_C\frac{ϵ}{2\pi \mathrm{i}}\mathrm{\Delta }F^{}(ϵ).$$
(B.33)
This can be expressed, by using the Gor’kov equations, as

$$E_{\mathrm{anomal}}=\frac{1}{2}\mathrm{Tr}_C\frac{ϵ}{2\pi \mathrm{i}}(ϵh+\mu )G_\mathrm{s}(ϵ).$$
(B.34)
On the other hand, the change of the normal part of the energy functional, in first-order perturbation theory, is given by

$$E_{\mathrm{normal}}[\rho ]E_{\mathrm{normal}}[\rho _0]=\mathrm{Tr}(h\delta \rho )\mu \delta N+\mathrm{Tr}[(h\mu )\delta \rho ].$$
(B.35)
This expression may be written in another form:

$$E_{\mathrm{normal}}[\rho ]E_{\mathrm{normal}}[\rho _0]=\mu \delta N+\mathrm{Tr}_C\frac{ϵ}{2\pi \mathrm{i}}(h\mu )(G_\mathrm{s}(ϵ)G_0(ϵ)).$$

Adding $`E_{\mathrm{anomal}}`$, we find the first (“direct”) change in the energy connected with the density variation $`\delta \rho `$:

$`\delta E=\mu \delta N`$ $`+`$ $`{\displaystyle \frac{1}{2}}\mathrm{Tr}{\displaystyle _C}{\displaystyle \frac{dϵ}{2\pi \mathrm{i}}}(ϵ+h\mu )G_\mathrm{s}(ϵ)`$ (B.36)
$``$ $`\mathrm{Tr}{\displaystyle _C}{\displaystyle \frac{ϵ}{2\pi \mathrm{i}}}(h\mu )G_0(ϵ),`$
which can also be written, by using the Gor’kov equations (B.3), as

$$\delta E=\mu \delta N+\frac{1}{2}\mathrm{Tr}_C\frac{ϵ}{2\pi \mathrm{i}}(ϵ+h\mu )G_0(ϵ)\mathrm{\Delta }F^{}(ϵ).$$
(B.37)
For the second (“polarization”) contribution to the total pairing energy in eq. (B.32), related to the density variation (B.28), we get:

$`\delta E_0`$ $`=`$ $`\mathrm{Tr}(h_{\mathrm{HF}}\delta \rho _0)=\mu _{\mathrm{HF}}\delta N_0`$ (B.38)
$`+`$ $`\mathrm{Tr}{\displaystyle _C}{\displaystyle \frac{ϵ}{2\pi \mathrm{i}}}(h_{\mathrm{HF}}\mu _{\mathrm{HF}})G_{\mathrm{HF}}(ϵ)(\delta U_{\mathrm{pair}}\delta \mu )G_0(ϵ).`$
The last term vanishes in first order. Thus, for the second change in the total energy, again in the first order, we find

$$\delta E_0=\mu _{\mathrm{HF}}\delta N_0\mu \delta N_0.$$
(B.39)
Collecting both contributions and taking into account the particle conservation (B.30), which gives $`\mu (\delta N+\delta N_0)=0`$, we arrive at the resulting formula for the pairing energy:

$$E_{\mathrm{pair}}=\frac{1}{2}\mathrm{Tr}\frac{ϵ}{2\pi \mathrm{i}}(ϵ+h\mu )G_0(ϵ)\mathrm{\Delta }F^{}(ϵ).$$
(B.40)
In $`k`$–representation, using the spectral decompositions of eqs. (B.8) and (B.14), this can be written as

$`E_{\mathrm{pair}}`$ $`=`$ $`{\displaystyle \frac{1}{4}}{\displaystyle \underset{i}{}}{\displaystyle \underset{(E_\alpha >0)}{}}{\displaystyle \frac{E_\alpha |ϵ_i\mu |}{E_\alpha +|ϵ_i\mu |}}{\displaystyle \underset{j}{}}\mathrm{\Delta }_{ij}`$ (B.41)
$`\left[(V_{j\alpha }U_{i\alpha }^{}U_{j\alpha }^{}V_{i\alpha })(1n_i)(V_{j\alpha }U_{i\alpha }^{}+U_{j\alpha }^{}V_{i\alpha })\right].`$
In diagonal approximation we get:

$$E_{\mathrm{pair}}=\frac{1}{2}\underset{k}{}\frac{|\mathrm{\Delta }_k|^2}{2E_k}\frac{E_k|ϵ_k\mu |}{E_k+|ϵ_k\mu |}.$$
(B.42)
For infinite nuclear matter the sum over $`k`$ is replaced by an integral in the momentum space, and one gets the pairing energy density (for one kind of particles):

$$\epsilon _{\mathrm{pair}}=\frac{1}{2}\frac{\stackrel{}{p}}{(2\pi )^3}\frac{|\mathrm{\Delta }(\stackrel{}{p})|^2}{E(\stackrel{}{p})}\frac{E(\stackrel{}{p})|ϵ(\stackrel{}{p})ϵ_\mathrm{F}|}{E(\stackrel{}{p})+|ϵ(\stackrel{}{p})ϵ_\mathrm{F}|}.$$
(B.43)
In the simplest model we assume

$$|\mathrm{\Delta }(\stackrel{}{p})|=|\mathrm{\Delta }(p_\mathrm{F})|\mathrm{\Delta },ϵ(\stackrel{}{p})=\frac{\stackrel{}{p}^2}{2m},$$
(B.44)
which leads to the expression

$$\epsilon _{\mathrm{pair}}=\frac{3}{8}\varrho _0\frac{\mathrm{\Delta }^2}{ϵ_\mathrm{F}}_1^{\mathrm{}}t\frac{\sqrt{1+t}}{\sqrt{\delta ^2+t^2}}\frac{\sqrt{\delta ^2+t^2}|t|}{\sqrt{\delta ^2+t^2}+|t|},$$
(B.45)
with

$$t=\frac{ϵ(\stackrel{}{p})}{ϵ_\mathrm{F}}1,\delta =\frac{\mathrm{\Delta }}{ϵ_\mathrm{F}},\varrho _0=\frac{p_\mathrm{F}^3}{3\pi ^2}.$$

The integral in (B.45) at $`\delta 1`$ can be easily estimated to yield
$`1\delta ^2(6\mathrm{ln}2+12\mathrm{ln}\delta )/32`$ \+ higher order corrections in $`\delta `$. Numerically this integral deviates very little from unity even at unreasonably large $`\delta `$ (less than by 2.5% up to $`\delta 0.5`$).

Thus, in the weak pairing approximation for an infinite system, as long as the condition $`\mathrm{\Delta }ϵ_\mathrm{F}`$ holds (which is generally assumed for nuclear matter near the saturation point), the pairing energy density is

$$\epsilon _{\mathrm{pair}}=\frac{3}{8}\varrho _0\frac{\mathrm{\Delta }^2(p_\mathrm{F})}{ϵ_\mathrm{F}}=\frac{1}{4}\frac{m^{}p_\mathrm{F}}{\pi ^2}\mathrm{\Delta }^2(p_\mathrm{F}).$$
(B.46)
We want to conclude this Appendix with the remark that in the case of contact pairing force leading to a local pairing field $`\mathrm{\Delta }(\stackrel{}{r})`$ there is no difficulty related with the cutoff $`ϵ_\mathrm{c}`$ to calculate the total energy of the system: in the energy space as a function of $`ϵ_\mathrm{c}`$ it converges rapidly ($`\delta _\mathrm{c}E_{\mathrm{pair}}1/ϵ_\mathrm{c}\sqrt{ϵ_\mathrm{c}}`$) although kinetic and interaction energies, if taken separately, diverge ($`+\sqrt{ϵ_\mathrm{c}}`$ and $`\sqrt{ϵ_\mathrm{c}}`$, respectively).

## Appendix C. Solving the Gor’kov equations in the coordinate-space representation

In this Appendix we describe the coordinate-space technique which allows, in the case of spherical symmetry, to solve the Gor’kov equations for given local mean-field and pairing potentials exactly. This technique has been invented in , where the most interesting effects, arising due to the proper treatment of the coupling between bound orbitals and particle continuum, were carefully investigated. Among them are the lowering of the chemical potential $`\mu `$ (an increase of the binding energy), the term-repulsion phenomenon and the appearance of width for deep-hole states lying below $`2\mu `$ from the continuum threshold. It has then been applied, within the EDF approach with a density-independent pairing $`\delta `$-force, to some superfluid even-even nuclei .

The fact that the coordinate-space HFB approach for finite systems can properly treat the positive-energy continuum, in contrast to BCS-like methods in which the presence of an unphysical “particle gas” is almost unavoidable, was apparently first pointed out in . It has been shown there that the HFB approach naturally leads to a localized nuclear wave function and gives the correct asymptotics for the normal and anomalous densities. These properties of the HFB solutions are especially important for the nuclei close to the drip lines. A method of solving the HFB equation directly in the coordinate representation, but approximating the continuum by a set of the discrete states in a spherical box, has been used in  for the self-consistent description of a chain of tin isotopes with the Skyrme effective interaction. The asymptotic behavior of the quasiparticle wave functions, as well as of the density distributions, have been studied with this method .

Here we give a more detailed description of the technique suggested in  for spherical even-even nuclei. We shall present also, by using the so-called “uniform filling approximation” to preserve the spherical symmetry, an extension of this technique for the odd systems. With this extension the blocking effect appears in a natural way. In what follows, the isospin variables are omitted to simplify the notation, and the derivation is the same both for neutrons and protons.

In order to separate the spin-angular variables in the Gor’kov matrix equation (B.1), it is convenient to introduce a unitary transformation for the generalized Green’s function:

$$\widehat{G}=\widehat{T}^{}\widehat{𝒢}\widehat{T},\widehat{𝒢}=\widehat{T}\widehat{G}\widehat{T}^{},$$
(C.1)
where $`\widehat{T}`$ is given by

$$\widehat{T}=\left(\begin{array}{cc}\delta (1,2)& 0\\ 0& 𝒯(1,2)\end{array}\right).$$
(C.2)
The operator $`𝒯(1,2)`$ is defined by eq. (A.63). Applying this transformation to (B.1), we get the equation for $`\widehat{𝒢}`$:

$$(ϵ\widehat{\stackrel{~}{}})\widehat{𝒢}=\widehat{I},$$
(C.3)
where the transformed Hamiltonian reads

$`\widehat{\stackrel{~}{}}=\widehat{T}\widehat{}\widehat{T}^{}`$ $`=`$ $`\left(\begin{array}{cc}\widehat{h}\mu & \widehat{\mathrm{\Delta }}\widehat{𝒯}^{}\\ \widehat{𝒯}\widehat{\mathrm{\Delta }}^{}& \mu \widehat{𝒯}\widehat{h}^{}\widehat{𝒯}^{}\end{array}\right)`$ (C.6)
$`=`$ $`\left(\begin{array}{cc}h(1,2)\mu \delta (1,2)& \mathrm{\Delta }(\stackrel{}{r}_1)\delta (1,2)\\ \mathrm{\Delta }^{}(\stackrel{}{r}_1)\delta (1,2)& \mu \delta (1,2)h(1,2)\end{array}\right).`$ (C.9)
Here time-reversal invariance, $`\widehat{𝒯}\widehat{h}^{}\widehat{𝒯}^{}=\widehat{h}`$, is assumed (no magnetic field or rotation) and the local (real) pairing field $`\mathrm{\Delta }(\stackrel{}{r})`$, which is diagonal in spin variables. is introduced according to eq. (A.66). The transformed generalized Green’s function $`\widehat{𝒢}`$ entering eq. (C.3) is given by

$$\widehat{𝒢}(1,2;ϵ)=\left(\begin{array}{cc}G_\mathrm{s}& F𝒯^{}\\ 𝒯F^{}& 𝒯\overline{G}_\mathrm{s}𝒯^{}\end{array}\right).$$
(C.10)
Because of the time-reversal symmetry, required for the Bogolyubov quasiparticle vacuum, this transformed Green’s function may be written, for even-even spherical nuclei, in a form where the spin-angular parts are the same for all its $`2\times 2`$ components:

$$\widehat{𝒢}(1,2;ϵ)=\frac{1}{r_1r_2}\underset{jlm}{}\widehat{g}_{jl}(r_1,r_2;ϵ)\varphi _{jlm}(\stackrel{}{n}_1,s_1)\varphi _{jlm}^{}(\stackrel{}{n}_2,s_2),$$
(C.11)
with $`\varphi _{jlm}(\stackrel{}{n},s)=C_{lms\frac{1}{2}s}^{jm}\mathrm{i}^lY_{lms}(\stackrel{}{n})`$ being the $`s`$-component of the usual spherical spinors ($`s=\pm \frac{1}{2}`$). For the generalized radial Green’s function $`\widehat{g}_{jl}`$, which has the matrix form

$$\widehat{g}_{jl}(r,r^{};ϵ)=\left(\begin{array}{cc}g_{jl}^{11}(r,r^{};ϵ)& g_{jl}^{12}(r,r^{};ϵ)\\ g_{jl}^{21}(r,r^{};ϵ)& g_{jl}^{22}(r,r^{};ϵ)\end{array}\right),$$
(C.12)
one gets the equation

$$\left(\begin{array}{cc}ϵh_{jl}+\mu & \mathrm{\Delta }\\ \mathrm{\Delta }& ϵ+h_{jl}\mu \end{array}\right)\widehat{g}_{jl}(r_1,r_2;ϵ)=\left(\begin{array}{cc}\delta (r_1r_2)& 0\\ 0& \delta (r_1r_2)\end{array}\right),$$
(C.13)
where $`h_{jl}=\frac{\mathrm{}^2}{2m}(\frac{d^2}{dr^2}+\frac{l(l+1)}{r^2})+U_{jl}(r)`$ is the single-quasiparticle Hamiltonian in the $`jl`$ channel (see below).

The solution of the matrix equation (C.13) can be constructed by using the set of the four linearly independent solutions

$$\widehat{y}_{i,jl}(r)=\left(\begin{array}{c}u_{i,jl}(r)\\ v_{i,jl}(r)\end{array}\right),i=\text{1–4},$$

which satisfy the homogeneous system of equations obtained from (C.13) by setting the right hand side to zero:

$`\left({\displaystyle \frac{\mathrm{}^2}{2m}}{\displaystyle \frac{d^2}{dr^2}}f_{jl}(r)+\mu +ϵ\right)u_{i,jl}(r)\mathrm{\Delta }(r)v_{i,jl}(r)=0,`$
$`\mathrm{\Delta }(r)u_{i,jl}(r)+\left({\displaystyle \frac{\mathrm{}^2}{2m}}{\displaystyle \frac{d^2}{dr^2}}f_{jl}(r)+\mu ϵ\right)v_{i,jl}(r)=0.`$ (C.14)
Here

$$f_{jl}(r)=U_\mathrm{c}(r)+U_{sl}(r)\stackrel{}{\sigma }\stackrel{}{l}_{jl}+\frac{\mathrm{}^2l(l+1)}{2mr^2},$$

$`U_\mathrm{c}`$ and $`U_{sl}`$ are the central and spin-orbit mean-field potentials, respectively, $`\stackrel{}{\sigma }\stackrel{}{l}_{jl}=j(j+1)l(l+1)\frac{3}{4},m`$ is the nucleon mass. We shall seek the radial Green’s function $`\widehat{g}_{jl}`$ entering eq. (C.11) in the following form (for simplicity, the $`ϵ`$-dependence and the $`jl`$ indices are not shown):

$`\widehat{g}(r_1,r_2)`$ $`=`$ $`{\displaystyle \frac{2m}{\mathrm{}^2}}[(\alpha \widehat{y}_1(r_1)+\beta \widehat{y}_2(r_1))\stackrel{~}{\widehat{y}}_3(r_2)`$ (C.15)
$`+`$ $`(\gamma \widehat{y}_1(r_1)+\delta \widehat{y}_2(r_1))\stackrel{~}{\widehat{y}}_4(r_2)]\theta (r_2r_1)`$
$`+`$ $`{\displaystyle \frac{2m}{\mathrm{}^2}}[\widehat{y}_3(r_1)(\alpha \stackrel{~}{\widehat{y}}_1(r_2)+\beta \stackrel{~}{\widehat{y}}_2(r_2))`$
$`+`$ $`\widehat{y}_4(r_1)(\gamma \stackrel{~}{\widehat{y}}_1(r_2)+\delta \stackrel{~}{\widehat{y}}_2(r_2))]\theta (r_1r_2),`$
where $`\stackrel{~}{\widehat{y}}_i=(u_i,v_i)`$ is the transposed $`i`$-th solution, $`\theta (x)`$ is the step function, $`\alpha ,\beta ,\gamma `$, and $`\delta `$ are the coefficients to be determined. The solutions $`\widehat{y}_1`$ and $`\widehat{y}_2`$ are chosen to be regular at $`r0`$ and the other two, $`\widehat{y}_3`$ and $`\widehat{y}_4`$, to be regular at $`r\mathrm{}`$:

$$\begin{array}{cc}\widehat{y}_1(r0)=(\begin{array}{c}1\\ \zeta _0\end{array})(q_+r)^{l+1},\hfill & \widehat{y}_2(r0)=(\begin{array}{c}\zeta _0\\ 1\end{array})(q_{}r)^{l+1},\hfill \\ \widehat{y}_3(r\mathrm{})=(\begin{array}{c}1\\ \zeta _{\mathrm{}}\end{array})\mathrm{}^{\mathrm{i}k_+r},\hfill & \widehat{y}_4(r\mathrm{})=(\begin{array}{c}\zeta _{\mathrm{}}\\ 1\end{array})\mathrm{}^{\mathrm{i}k_{}r}.\hfill \end{array}$$
(C.16)
Here

$$\begin{array}{c}q_\pm =\left[\frac{2m}{\mathrm{}^2}\left(\mu U_0\pm \sqrt{ϵ^2\mathrm{\Delta }_0^2}\right)\right]^{1/2},\zeta _0=\mathrm{\Delta }_0/\left(ϵ+\sqrt{ϵ^2\mathrm{\Delta }_0^2}\right),\hfill \\ k_\pm =\left[\frac{2m}{\mathrm{}^2}\left(\mu \pm \sqrt{ϵ^2\mathrm{\Delta }_{\mathrm{}}^2}\right)\right]^{1/2},\zeta _{\mathrm{}}=\mathrm{\Delta }_{\mathrm{}}/\left(ϵ+\sqrt{ϵ^2\mathrm{\Delta }_{\mathrm{}}^2}\right);\hfill \\ U_0=U_\mathrm{c}(0)+U_{sl}(0)\stackrel{}{\sigma }\stackrel{}{l},\mathrm{\Delta }_0=\mathrm{\Delta }(0),\mathrm{\Delta }_{\mathrm{}}=\mathrm{\Delta }(r\mathrm{}).\hfill \end{array}$$
(C.17)
It is important to notice that on the physical $`k`$–sheet the imaginary parts of the asymptotic momenta $`k_\pm `$ should be chosen positive141414To be more precise, the solutions $`\widehat{y}_3`$ and $`\widehat{y}_4`$ are given asymptotically at $`r\mathrm{}`$ by the Whittaker functions $`W_{\mathrm{i}\eta _\pm ,l+1/2}(2\mathrm{i}k_\pm r)`$, respectively, with $`\eta _\pm =mZe^2/\mathrm{}^2k_\pm `$ being the corresponding Coulomb parameters ($`\eta _\pm =0`$ for neutrons).: $`\mathrm{Im}k_\pm >0`$.

Using the system of equations (C.14), it can be easily verified that for any two arbitrary solutions,

$$\widehat{y}_i=(\begin{array}{c}u_i\\ v_i\end{array})\mathrm{and}\widehat{y}_j=(\begin{array}{c}u_j\\ v_j\end{array}),$$

the combination $`C_{ij}=U_{ij}V_{ij}`$ composed of the Wronskians $`U_{ij}=u_iu_j^{}u_i^{}u_j`$ and $`V_{ij}=v_iv_j^{}v_i^{}v_j`$, where the primes denote differentiation with respect to $`r`$, does not depend on the coordinate $`r`$, and for solutions satisfying the boundary conditions (C.16) one gets $`C_{12}=C_{34}=0`$. It can be seen also that $`U_{13}U_{24}U_{23}U_{14}=U_{12}U_{34}`$ and $`V_{13}V_{24}V_{23}V_{14}=V_{12}V_{34}`$. With these relations one easily obtains an equation for determining the four coefficients entering eq. (C.15):

$$\left(\begin{array}{cc}\alpha & \beta \\ \gamma & \delta \end{array}\right)=\left(\begin{array}{cc}C_{13}& C_{14}\\ C_{23}& C_{24}\end{array}\right)^1=\frac{1}{D}\left(\begin{array}{cc}C_{24}& C_{14}\\ C_{23}& C_{13}\end{array}\right),$$
(C.18)
where $`D(ϵ)=C_{13}C_{24}C_{23}C_{14}`$ is the determinant of the matrix $`C`$.

The formulas (C.1), (C.11) and (C.15)–(C.18) constitute an exact solution of the Gor’kov equation for the generalized Green’s function in the coordinate representation for superfluid spherical nuclei with a localized pairing field $`\mathrm{\Delta }(\stackrel{}{r})`$.

In the case of vanishing pairing one has $`u_2=u_4=v_1=v_3=0`$ and $`\beta =\gamma =0`$. Then the radial generalized Green’s function (C.12) becomes diagonal with $`g^{11}`$ and $`g^{22}`$ in the following form:

$`g_{jl}^{11}(r,r^{};ϵ)|_{\mathrm{\Delta }=0}`$ $`=`$ $`{\displaystyle \frac{2m}{\mathrm{}^2}}{\displaystyle \frac{u_{1;jl}(r_<;ϵ)u_{3;jl}(r_>;ϵ)}{U_{13}(ϵ)}},`$
$`g_{jl}^{22}(r,r^{};ϵ)|_{\mathrm{\Delta }=0}`$ $`=`$ $`{\displaystyle \frac{2m}{\mathrm{}^2}}{\displaystyle \frac{v_{2;jl}(r_<;ϵ)v_{4;jl}(r_>;ϵ)}{V_{24}(ϵ)}},`$ (C.19)
where $`r_<`$ and $`r_>`$ denote the lesser and the greater of $`r`$ and $`r^{}`$, respectively. Here $`g^{11}`$ is interpreted as the HF Green’s function for particles ($`ϵ>\mu `$) while $`g^{22}`$ as that for holes ($`ϵ<\mu `$)151515Eq.(C.19) for the single-particle Green’s function allows one to construct the particle-hole propagator in the coordinate representation and to solve, for closed-shell nuclei, the RPA equations exactly, including the whole particle continuum. An earlier application of this method for the study of pion condensation in finite nuclear systems may be found in , and, for calculating the continuum–RPA multipole response functions, in .. By integrating the latter along the contour $`C`$ of Fig. 21 one gets the HF $`jl`$ density matrix $`\rho _{jl}(r,r^{})|_{\mathrm{\Delta }=0}`$ which is normalized by $`(2j+1)r\rho _{jl}(r,r^{})|_{\mathrm{\Delta }=0}=N_{jl}`$ where $`N_{jl}`$ is the number of particles occupying orbitals with the given quantum numbers $`jl`$. This density matrix can be expressed by a finite sum of separable terms, i.e. by a sum of residues of $`g^{22}|_{\mathrm{\Delta }=0}`$ at the poles which are zeros of the Wronskian $`V_{24}`$. The equation $`V_{24}=0`$ defines the spectrum of the HF hole states, and this spectrum is utterly discrete.

The situation with pairing is completely different. It is clear that the poles (if any) of the generalized Green’s function are determined by the equation $`D(ϵ)=0`$ which gives the discrete spectrum of the Bogolyubov quasiparticle states. Due to the invariance of the equation (C.14) under simultaneous replacements $`ϵϵ`$ and $`(u,v)(v,u)`$ the energy spectrum is symmetrical around the point $`ϵ=0`$. For bound systems, by definition, one has $`\mu <0`$. If the system is finite, i.e. if the density distribution has a finite range, one expects $`\mathrm{\Delta }_{\mathrm{}}=0`$, otherwise the system would be unstable with respect to two particle emission. From eqs. (C.15)–(C.18) it follows then that the quasiparticle spectrum is discrete within the energy region $`|ϵ|<|\mu |`$ (here $`\widehat{g}_{jl}`$ is a real function at the axis $`\mathrm{Im}ϵ=0`$) and continuous if $`|ϵ|>|\mu |`$ (there $`\widehat{g}_{jl}`$ is a complex function). These features are reflected in Fig. 21 where the branch cuts are shown by the heavy lines extending symmetrically to the left and to the right from the points $`\pm \mu `$, respectively, along the $`\mathrm{Im}ϵ=0`$ axis. The case with only the branch cuts, without poles, corresponds to a drip-line even nucleus. Performing the integration of the radial Green’s function $`\widehat{g}_{jl}`$ along the contour $`C`$ of Fig. 21 one obtains the radial part of the generalized density matrix161616The analytical properties of the Green’s function allow one to consider only the upper half of the contour with Im$`ϵ>0`$: integrate Re$`G`$ along the vertical sections and Im$`G`$ along the horizontal section in the upper half-plane, then for the total result take the doubled sum of these integrals.. To calculate the normal and anomalous densities we need the following two radial “$`jl`$ density matrices”:

$$\rho _{jl}(r,r^{})=_C\frac{ϵ}{2\pi \mathrm{i}}g_{jl}^{22}(r,r^{};ϵ),\nu _{jl}(r,r^{})=_C\frac{ϵ}{2\pi \mathrm{i}}g_{jl}^{12}(r,r^{};ϵ).$$
(C.20)
Then the resulting expression for the normal density reads

$$\rho (\stackrel{}{r})=\frac{1}{4\pi r^2}\underset{jl}{}(2j+1)\rho _{jl}(r,r).$$
(C.21)
It is normalized by $`\stackrel{}{r}\rho (\stackrel{}{r})=N`$ which is actually the condition (2.23) for the chemical potential $`\mu `$. For the anomalous density from (A.37), (A.64), (C.10), (C.11) and (C.20) we get

$$\nu (\stackrel{}{r})=\frac{1}{4\pi r^2}\underset{jl}{}(j+\frac{1}{2})\nu _{jl}(r,r).$$
(C.22)
These densities are used to calculate both the normal and anomalous parts of the interaction energy in eq. (2.22). In addition, for the spin-orbital term, eq. (3.5), we define171717Remember that we do not show here the cutoff index $`c`$. It should be emphasized that because of the continuum all the densities defined in this Appendix cannot be expressed by a finite sum of separable terms. The latter would be possible only within a certain approximation, namely with a discretized basis and with an energy cutoff. The coordinate-space technique described above, which involves integration of the exactly constructed generalized Green’s function in the complex energy plane, appears to be exact for the $`jl`$–densities: both $`\rho _{jl}(r,r)`$ and $`\nu _{jl}(r,r)`$ converge as functions of $`ϵ_\mathrm{c}`$ ($`\delta _\mathrm{c}\rho _{jl}ϵ_\mathrm{c}^1`$ and $`\delta _\mathrm{c}\nu _{jl}ϵ_\mathrm{c}^{1/4}`$, respectively). It is exact also for the normal density C.21) since it converges, $`\delta _\mathrm{c}\rho ϵ_\mathrm{c}^{1/2}`$. However, if one goes beyond the cutoff $`ϵ_\mathrm{c}`$, the anomalous density (C.22) diverges, $`\nu ϵ_\mathrm{c}^{1/4}`$. The kinetic and anomalous energies beyond $`ϵ_\mathrm{c}`$ diverge too, but their sum converges very rapidly ($`\delta _\mathrm{c}E_{\mathrm{pair}}ϵ_\mathrm{c}^{3/2}`$, see Appendix B). This fact has been used in the present paper to formulate the generalized variational principle with the cutoff local EDF. We stress once more that this formulation does not imply a cutoff of the basis.:

$$\rho _{sl}(\stackrel{}{r})=\frac{1}{4\pi r^2}\underset{jl}{}(2j+1)\stackrel{}{\sigma }\stackrel{}{l}_{jl}\rho _{jl}(r,r).$$
(C.23)
The kinetic energy of the system is given by

$$E_{\mathrm{kin}}=\frac{\mathrm{}^2}{2m}r\underset{jl}{}(2j+1)\left\{\left[\frac{^2}{r^2}+\frac{l(l+1)}{r^2}\right]\rho _{jl}(r,r^{})\right\}_{r^{}=r}.$$
(C.24)
Note that to compute this energy one needs not only the normal $`jl`$-densities, but also the $`jl`$ normal density matrices $`\rho _{jl}(r,r^{})`$.

For an odd system the spectral expansion of the generalized Green’s function, eqs. (B.9)–(B.10), should be modified. Suppose that the addition of an odd particle leads to the appearance of the quasiparticle with the energy $`E_{\alpha _0}`$ in the ground state of the system. Here $`\alpha _0`$ denotes all the relevant quantum numbers. They can be specified, for example, by $`\alpha _0=\{njlm\}_0`$ from the whole set $`\{\alpha \}`$ used to label the solutions of the HFB equation (A.28). Suppose further that $`E_\alpha `$ corresponds to a certain eigenenergy of this equation and belongs to the discrete spectrum, i.e. $`|\mu |>E_{\alpha _0}>0`$. As illustrated in Fig. 22, $`E_{\alpha _0}`$ is located on the $`\mathrm{Im}ϵ=0`$ axis in the vicinity of the point $`ϵ=0`$. Generally, $`E_{\alpha _0}`$ is of the order of $`\overline{\mathrm{\Delta }}`$, the average matrix element of the pairing potential on the the Fermi surface. Note that the pole nearest to $`ϵ=0`$ need not correspond to the ground state of the system. The case of $`E_{\alpha _0}\mu `$ determines the position of the drip line for odd nuclei. The rules for passing the Green’s function poles for superfluid odd systems are formulated by Migdal . To the spectral expansions (B.9)–(B.10), which do not contain the contribution of the odd quasiparticle, one should add the following terms:

$`\delta G_\mathrm{s}^{\mathrm{odd}}(1,2;ϵ)`$ $`=`$ $`{\displaystyle \frac{V_{\alpha _0}^{}(1)V_{\alpha _0}(2)}{ϵ+E_{\alpha _0}\mathrm{i}\delta }}+{\displaystyle \frac{U_{\alpha _0}(1)U_{\alpha _0}^{}(2)}{ϵE_{\alpha _0}\mathrm{i}\delta }}`$ (C.25)
$`+`$ $`{\displaystyle \frac{V_{\overline{\alpha }_0}^{}(1)V_{\overline{\alpha }_0}(2)}{ϵ+E_{\overline{\alpha }_0}+\mathrm{i}\delta }}+{\displaystyle \frac{U_{\overline{\alpha }_0}(1)U_{\overline{\alpha }_0}^{}(2)}{ϵE_{\overline{\alpha }_0}+\mathrm{i}\delta }},`$
$`\delta F^{\mathrm{odd}}(1,2;ϵ)`$ $`=`$ $`{\displaystyle \frac{V_{\alpha _0}^{}(1)U_{\alpha _0}(2)}{ϵ+E_{\alpha _0}\mathrm{i}\delta }}+{\displaystyle \frac{U_{\alpha _0}(1)V_{\alpha _0}^{}(2)}{ϵE_{\alpha _0}\mathrm{i}\delta }}`$ (C.26)
$`+`$ $`{\displaystyle \frac{V_{\overline{\alpha }_0}^{}(1)U_{\overline{\alpha }_0}(2)}{ϵ+E_{\overline{\alpha }_0}+\mathrm{i}\delta }}+{\displaystyle \frac{U_{\overline{\alpha }_0}(1)V_{\overline{\alpha }_0}^{}(2)}{ϵE_{\overline{\alpha }_0}+\mathrm{i}\delta }}.`$
The subscript $`\alpha _0`$ in these expressions, in Migdal’s terminology, refers to an odd quasiparticle added to a system of $`N`$ particles and to an odd quasihole in a system of $`N+2`$ particles. Consequently, the subscript $`\overline{\alpha }_0`$ refers to the odd quasihole in the system of $`N`$ particles and to the odd quasiparticle added to a system of $`N2`$ particles. Thus the four terms in (C.25) and (C.26) correspond to pairs of time-reversed states. But, due to the signs of the infinitesimal imaginary parts $`\pm \mathrm{i}\delta `$ in the denominators, only two, not time-reversed, terms in these expressions should contribute to the generalized density matrix. Closing the contour in the upper half of the energy plane, the terms from the first lines in eqs. (C.25) and (C.26) will be accounted for. By integrating the Green’s functions $`G_\mathrm{s}`$ and $`F`$ along the contour $`C`$ of Fig. 21, all the residues at the pole $`ϵ=E_{\alpha _0}`$ $`(=E_{\overline{\alpha }_0})`$ have been already included. Therefore, from the so obtained density matrix we have to subtract the contributions related to the third terms of these equations and to add the residues from the second ones. The resulting contour for the odd system is depicted in Fig. 22 which clearly illustrates the blocking effect: the presence of the odd particle in the level $`\alpha _0`$ prevents it from participating in the pairing correlations because in this case the “conjugate” level $`\overline{\alpha }_0`$ should be empty. This way we get the contributions to the normal and abnormal density matrices from the odd qusiparticle:

$$\delta \rho ^{\mathrm{odd}}(1,2)=V_{\overline{\alpha }_0}^{}(1)V_{\overline{\alpha }_0}(2)+U_{\alpha _0}(1)U_{\alpha _0}^{}(2),$$
(C.27)
$$\delta \nu ^{\mathrm{odd}}(1,2)=V_{\overline{\alpha }_0}^{}(1)U_{\overline{\alpha }_0}(2)+U_{\alpha _0}(1)V_{\alpha _0}^{}(2).$$
(C.28)
In the following we shall consider the case when only the spherical part of the core polarization induced by the odd quasiparticle is important, and the deformation of the core can be neglected. Then the solutions of the Bogolyubov equation may be chosen to have the form

$`U_\alpha (\stackrel{}{r},s)`$ $`=`$ $`{\displaystyle \frac{1}{r}}u_{njl}(r)\varphi _{jlm}(\stackrel{}{n},s),`$
$`V_\alpha (\stackrel{}{r},s)`$ $`=`$ $`{\displaystyle \frac{1}{r}}v_{njl}(r)(1)^{jm+1}\varphi _{jlm}^{}(\stackrel{}{n},s),`$ (C.29)
where the radial wave functions $`u_{njl}(r)`$ and $`v_{njl}(r)`$ satisfy eq. (C.14). For the discrete states, including the odd quasiparticle state with $`\alpha =\alpha _0`$, the latter obey the orthogonality relation $`r[u_{n^{}lj}^{}(r)u_{nlj}(r)+v_{n^{}lj}^{}(r)v_{nlj}(r)]=\delta _{nn^{}}`$. We introduce further the uniform filling approximation. Namely, we assume that the odd quasiparticle occupies all the possible $`m_0`$ states with equal probability and take the averages

$`\overline{\delta \rho ^{\mathrm{odd}}}(1,2)`$ $`=`$ $`{\displaystyle \frac{1}{2j_0+1}}{\displaystyle \underset{m_0}{}}\delta \rho ^{\mathrm{odd}}(1,2),`$
$`\overline{\delta \nu ^{\mathrm{odd}}}(1,2)`$ $`=`$ $`{\displaystyle \frac{1}{2j_0+1}}{\displaystyle \underset{m_0}{}}\delta \nu ^{\mathrm{odd}}(1,2).`$ (C.30)
Such an averaging procedure, with the wave functions given by (C.29), enables us to separate the spin-angular variables for the odd spherical system in the same way as done above for even systems. The resulting formulas for the contribution to the radial $`jl`$ density matrices (C.20) due to the presence of the odd particle read

$`\overline{\delta \rho _{n_0j_0l_0}^{\mathrm{odd}}}(r,r^{})`$ $`=`$ $`v_{n_0j_0l_0}^{}(r)v_{n_0j_0l_0}(r^{})+u_{n_0j_0l_0}(r)u_{n_0j_0l_0}^{}(r^{}),`$
$`\overline{\delta \nu _{n_0j_0l_0}^{\mathrm{odd}}}(r,r^{})`$ $`=`$ $`v_{n_0j_0l_0}^{}(r)u_{n_0j_0l_0}(r^{})u_{n_0j_0l_0}(r)v_{n_0j_0l_0}^{}(r^{}).`$ (C.31)
The contributions to the normal and anomalous densities (C.21) and (C.22) are given by

$$\overline{\delta \rho _{n_0j_0l_0}^{\mathrm{odd}}}(\stackrel{}{r})=\frac{1}{4\pi r^2}\left(|v_{n_0j_0l_0}(r)|^2+|u_{n_0j_0l_0}(r)|^2\right).$$
(C.32)
$$\overline{\delta \nu _{n_0j_0l_0}^{\mathrm{odd}}}(\stackrel{}{r})=\frac{1}{4\pi r^2}u_{n_0j_0l_0}(r)v_{n_0j_0l_0}^{}(r).$$
(C.33)
It can be seen that $`\overline{\delta \nu _{n_0j_0l_0}^{\mathrm{odd}}}(\stackrel{}{r})`$ has the opposite sign as $`\nu (\stackrel{}{r})`$, eq. (C.22), reflecting the fact that the blocked level cannot directly contribute to the anomalous energy. The nucleon separation energies $`S_n`$ are determined by the position of the poles close to $`ϵ=0`$. As easily understood from Fig. 22, since these separation energies are measured from the continuum threshold $`ϵ=\mu `$, for an odd system one gets $`S_n^{\mathrm{odd}}\mu E_{\alpha _0}`$, and, for an even system, $`S_n^{\mathrm{even}}\mu +E_{\overline{\alpha }_0}`$. Thus we have $`S_n^{\mathrm{even}}S_n^{\mathrm{odd}}2E_{\alpha _0}2\overline{\mathrm{\Delta }}`$, i.e. the familiar odd-even effect in nuclear masses.

The rest of this Appendix is devoted to the discussion of the asymptotic behavior of the densities $`\rho (\stackrel{}{r})`$ and $`\nu (\stackrel{}{r})`$. Consider first the asymptotic properties of the wave functions $`y_i`$ used to construct the radial Green’s functions (C.15) through which these density can be obtained. For finite systems, the functions $`v_3(r)`$ and $`u_4(r)`$ in the solutions $`y_3`$ and $`y_4`$ of the radial HFB equation (C.14), as given by the boundary conditions (C.16), vanish at large distances since $`\zeta _{\mathrm{}}\mathrm{\Delta }_{\mathrm{}}=0`$. The function $`u_3`$ at large $`r`$ behaves as $`\mathrm{exp}\mathrm{i}k_+r`$ and either decreases $`\mathrm{exp}(r\sqrt{2m(|ϵ|+\mu )/\mathrm{}^2})`$ if $`|ϵ|<\mu `$ or oscillates if $`|ϵ|>\mu `$, the latter corresponds to the continuum spectrum. The function $`v_4`$ behaves asymptotically as $`\mathrm{exp}\mathrm{i}k_{}r`$. If the system is bound, i.e. $`\mu <0`$, this function at $`r\mathrm{}`$ always decreases exponentially $`\mathrm{exp}(r\sqrt{2m(|ϵ|\mu )/\mathrm{}^2})`$.

To get the asymptotic laws for normal and anomalous densities, defined by eqs. (C.20)–(C.22), we have to consider now the radial Green’s functions $`g^{12}(\gamma u_1+\delta u_2)v_4`$ and $`g^{22}(\gamma v_1+\delta v_2)v_4`$. Asymptotically, at large $`r`$, the functions $`u_i`$ and $`v_i`$, $`i=1,2`$, which are regular at $`r0`$, are given by $`u_i(r)=c_i\mathrm{exp}(\mathrm{i}k_+r)+d_i\mathrm{exp}(\mathrm{i}k_+r)`$ and $`v_i(r)=e_i\mathrm{exp}(\mathrm{i}k_{}r)+f_i\mathrm{exp}(\mathrm{i}k_{}r)`$, respectively. Using eq. (C.18) one can easily verify that, by construction, the combination $`\gamma d_1+\delta d_2`$ vanishes. The combination $`\gamma f_1+\delta f_2`$ vanishes at the poles of the Green’s function determined by the equation $`D(ϵ)=0`$. Finally, as the result of of integration of the Green’s function along the contour $`C`$, we get

$`\rho (r\mathrm{})`$ $``$ $`{\displaystyle \frac{1}{r^2}}\mathrm{exp}(2r\sqrt{2m(|ϵ_0|\mu )/\mathrm{}^2}),`$
$`\nu (r\mathrm{})`$ $``$ $`{\displaystyle \frac{1}{r^2}}\mathrm{exp}(r[\sqrt{2m(|ϵ_0|\mu )/\mathrm{}^2}+\sqrt{2m(|ϵ_0|\mu )/\mathrm{}^2}]),`$ (C.34)
where $`ϵ_0`$ is the position, within the contour $`C`$, of a singularity point nearest to $`ϵ=0`$. This point may be a pole of the Green’s function or the beginning of the branch cut. The latter situation corresponds to the drip-line nuclei. Generally, $`|ϵ_0|`$ is of the order of the average matrix element $`\overline{\mathrm{\Delta }}`$ of the pairing potential on the Fermi surface (typically, $`\overline{\mathrm{\Delta }}1`$ MeV). In stable nuclei, since the absolute value of the chemical potential is much larger, $`\mu 8`$ MeV, the asymptotic behavior of both densities (C.34) differ but insignificantly. Approaching the drip-line, in even system, $`|ϵ_0|`$, $`|\mu |`$ and $`\overline{\mathrm{\Delta }}`$ may become of the same order: $`|ϵ_0|\mu \overline{\mathrm{\Delta }}`$. In such a case one gets

$`\rho (r\mathrm{})`$ $``$ $`{\displaystyle \frac{1}{r^2}}\mathrm{exp}(4r\sqrt{m\overline{\mathrm{\Delta }}/\mathrm{}^2}),`$
$`\nu (r\mathrm{})`$ $``$ $`{\displaystyle \frac{1}{r^2}}\mathrm{exp}(2r\sqrt{m\overline{\mathrm{\Delta }}/\mathrm{}^2}).`$ (C.35)
It follows that in this situation the anomalous density has a longer tail with decaying length by a factor of two greater than that of the normal density. Similar conclusions have been drawn in . On the other hand, in the genuine drip-line even nuclei, with only one bound state – the ground state, when all the HFB eigenstates are embeded in the continuum, there are always contributions to the densities from the quasiparticles with different quantum numbers $`jl`$ and with the energies close to $`\mu `$. The weights of these contributions are determined by the corresponding quasiparticle level density, but in any case the most distant tails of the densities in drip-line nuclei are asymptotically given by the expressions

$`\rho (r\mathrm{})`$ $``$ $`{\displaystyle \frac{1}{r^2}}\mathrm{exp}(4r\sqrt{m|\mu |/\mathrm{}^2}),`$
$`\nu (r\mathrm{})`$ $``$ $`{\displaystyle \frac{1}{r^2}}\mathrm{exp}(2r\sqrt{m|\mu |/\mathrm{}^2}),`$ (C.36)
which are obtained from eq. (C.34) by taking the limit $`|ϵ_0||\mu |`$. The technique described in this appendix permits, in principle, an exact calculation of the asymptotic behavior of the HFB densities in different situations.