Khon-Sham Density Functional Inspired Approach to Nuclear Binding

A non-relativisitic nuclear density functional theory is constructed, not as usual, from an effective density dependent nucleon-nucleon force but directly introducing in the functional results from microscopic nuclear and neutron matter Bruckner G-matrix calculations at various densities. A purely phenomenological finite range part to account for surface properties is added. The striking result is that only four to five adjustable parameters, spin-orbit included, suffice to reproduce nuclear binding energies and radii with the same quality as obtained with the most performant effective forces, containing on the order of ten parameters. In this pilot work, for the pairing correlations, simply a density dependent zero range force is adopted from the literature. Possible future extensions of this approach are pointed out.

A non-relativisitic nuclear density functional theory is constructed, not as done most of the time, from an effective density dependent nucleon-nucleon force but directly introducing in the functional results from microscopic nuclear and neutron matter Bruckner G-matrix calculations at various densities. A purely phenomenological finite range part to account for surface properties is added. The striking result is that only four to five adjustable parameters, spin-orbit included, suffice to reproduce nuclear binding energies and radii with the same quality as obtained with the most performant effective forces. In this pilot work, for the pairing correlations, simply a density dependent zero range force is adopted from the literature. Possible future extensions of this approach are pointed out. It is common use in nuclear physics to describe properties of nuclei in mean-field approximation, using effective density dependent nucleon-nucleon interactions. Prototypes of these interactions are the so-called Skyrme-forces [1] which are of zero range and density dependent. There also exist finite range versions, like the Gogny force [2] and the relativistic mean field (RMF) approach [3]. The number of parameters which enter these effective forces is typically around ten and they are adjusted to reproduce finite nuclei and some nuclear matter properties. However, recently also data from a theoretically determined neutron matter Equation of State (EOS) [4] have been used as input (see also an older attempt in this direction in [5]). Generally these forces give rise to an effective nucleon mass m * < m with typically m * /m ≃ 0.7 in the non-relativistic framework. With these ingredients nuclear mean field theories are very successfull to describe nuclear properties as, e.g. binding energies, radii, but also excited states, fission barriers, and many things more [3,6]. These effective forces are still under debate and constantly refined to account for the ever increasing set of nuclear data [4,7]. For instance the evolution with isospin, to be considered for nuclei approaching the neutron drip as involved in stellar nucleo-synthesis, is an active line of research. Also, usually ( apart from the Gogny force), additional parameters are needed to characterise nuclear pairing.
The main objective of the present study is to show that within an approach inspired by the Kohn-Sham Density Functional Theory (KS-DFT), like it is presently still used extensively in, e.g., condensed matter, chemistry, and atomic physics [9,10,11] ( see, however, also [12,13,14] for updated versions of the universal exchange-correlation energy including gradients and kinetic energy density dependences), a very efficient energy density functional for finite nuclei can be built up using explicitly a fully microscopic input from nuclear and neutron matter calculations. We are, however, aware of the fact that the applicability of the KS-DFT approach to self-bound systems, as nuclei, is not obvious and that presently vivid debate about this topic is going on [15,16,17,18,19]. In this work we shall not contribute anything fundamental to this discussion but rather adopt a pragmatic point of view: try and see! The KS-DFT scheme is somewhat different from the usual Skyrme or Gogny approaches [1,2] where, as mentioned, most of the time an effective force is constructed. However, several authors also interpret the Skyrme force as based on DFT [20,21] and a few authors have, in the past, constructed local Energy Density Functionals (EDF's) which cannot be related to underlying effective forces [22,23,24]. See also [25] for a review of different mean field calculations used in nuclear physics. In detail , however, all these nuclear approaches differ from what is practice in condensed matter physics (see below). The basis of KS-DFT lies in the Hohenberg-Kohn (HK) theorem [26], which states that for a Fermi system, with a non-degenerate ground state, the total energy can be expressed as a functional of the density ρ(r) only. Such a functional reaches its variational minimum when evaluated with the exact ground state density. Furthermore, for practical reasons, in the KS-DFT method one introduces an auxiliary set of A orthonormal single particle wave functions ψ i (r), where A is the number of particles, and the density is assumed to be given by where s and t stand for spin and iso-spin indices. The variational procedure to minimise the functional is performed in terms of the orbitals instead of the density. Usually in condensed matter and atomic physics the HK functional E[ρ(r)] is split into two parts: [8]. The first piece T 0 corresponds to the uncorrelated part of the kinetic energy and within the KS method it is written as The other piece W [ρ] contains the potential energy as well as the correlated part of the kinetic energy. Then, upon variation, one gets a closed set of A Hartree-like equations with an effective potential, the functional derivative of W [ρ] with respect to the local density ρ(r). Since the latter depends on the density, and therefore on the ψ i 's, a self-consistent procedure is necessary. Still the equations are exact but they only can be of some use if a reliable approximation is found for the otherwise unknown density functional W [ρ]. In the KS-DFT formalism the exact ground state wave function is actually not known, the density being the basic quantity.
In nuclear physics, contrary to the situation in condensed matter and atomic physics, the contribution of the spin-orbit interaction to the energy functional is very important. Non-local contributions have been included in DFT in several ways already long ago (see [27] for a recent review of this topic). Consequently, the spin-orbit part also can be split in an uncorrelated part E s.o. plus a remainder. The form of the uncorrelated spin-orbit part is taken exactly as in the Skyrme [1] or Gogny forces [2].
We thus write for the functional in the nuclear case where we explicitly split off the Coulomb energy E C because it is a quite distinct part in the Hamiltonian. It shall be treated, as usual, at lowest order, i.e. the direct term plus the exchange contribution in the Slater approximation, that is Let us now discuss the nuclear energy functional contribution E int [ρ n , ρ p ] which contains the nuclear potential energy as well as additional correlations. We shall split it in a finite range term E F R int [ρ n , ρ p ] to account for correct surface properties and a bulk correlation part E ∞ int [ρ n , ρ p ] that we take from a microscopic infinite nuclear matter calculation [29] as we will discuss below. Thus our final KS-DFT -like functional reads: For the finite range term we make the simplest phenomenological ansatz possible with t = proton/neutron and γ t,t ′ the volume integral of v t,t ′ (r). The substraction in (4) is made in order not to contaminate the bulk part, determined from the microscopic infinite matter calculation. Finite range terms have already been used earlier, generalizing usual Skyrme functionals (see e.g [23,30,31]). In this study, for the finite range form factor v t,t ′ (r) we make a simple Gaussian We choose a minimum of three open parameters: The only undetermined and most important piece in (3) is then the bulk contribution E ∞ int . In condensed matter, chemistry and atomic physics, this quantity implies the EOS of interacting electrons and for the KS-DFT scheme, most of the time, it is obtained from Quantum Monte-Carlo (QMC) calculations, the results of which are then accurately represented by a fit function (not necessarily a polynomial, as we use later). As already mentioned, we obtain E ∞ int from microscopic infinite matter calculations, using a realistic bare force, together with a converged hole-line expansion [29]. We first reproduce by interpolating functions the correlation part of the ground state energy per particle of symmetric and pure neutron matters, and then make a quadratic interpolation for asymmetric matter. Finally the total correlation contribution to the energy functional in local density approximation reads: where P s and P n are two interpolating polynomials for symmetric and pure neutron matter, respectively, at the density ρ = ρ p +ρ n , and β = (ρ n −ρ p )/ρ is the asymmetry parameter. For P s we took ( x = ρ/ρ 0 , with ρ 0 =0.17 fm −3 , see below ) where the coefficients ( energy/particle in MeV) are given in Table I. The two forms match at x = 1 (ρ = ρ 0 ) up to the second derivative. This functional form can be used up to ρ = 0.24 fm −3 , which is the interval where the independent fit of the microscopic calculation has been performed.
A similar expression holds for P n , where again the coefficients are displayed in Table I, which is valid in the same density interval. The interpolating polynomial for symmetric matter has been constrained to allow a minimum exactly at the energy E/A = -16 MeV and Fermi momentum k F = 1.36 fm −1 , i.e. ρ 0 = 0.17 fm −3 . This is within the uncertainty of the numerical microscopic calculations of the EOS [28]. The constrained fit was performed by keeping the EOS as smooth as possible, thus allowing for some very small deviations from the microscopic calculations below saturation density. An interpolating fit which goes exactly through the calculated EOS, as performed in [29], gives a not good enough saturation point (typically E/A = -15.6 MeV, k F =1.38 fm −1 ). In this context it should be noticed that the microscopic calculations of Ref. [32] have been modified, in order to get a good saturation point, by an ad hoc correction which amounts to more than 2 MeV in energy ( see, e.g. [33] for a presentation of the uncorrected EOS of [32]). As discussed in [33], the low density behavior of the nuclear matter EOS is quite intricate and usually not reproduced by Skyrme and Gogny functionals ( see also ref. [29]), missing quite a substantial part of binding. We show our EOS for nuclear and neutron matter in Fig. 1. One sees that for all densities relevant to finite nuclei the EOS is very accurately reproduced by the interpolation (continuous lines). The bulk part E ∞ int of our functional, directly related to bare NN and NNN forces, is, therefore, determined once and for all and we will use it in (3) together with LDA. The only open parameters are, therefore, the ones contained in the finite range surface part, Eq.(4), and the strength of the spin-orbit contribution. We, thus, follow exactly the strategy employed in the above mentioned interacting electron systems. On the contrary, up to now, in nuclear physics, almost exclusively a different strategy was in use (see, however, Ref. [34] for a somewhat similar approach to our): functionals like the one of (3), i.e. bulk, surface, etc., were globally parametrized with typically of the order of ten parameters which, then, were determined fitting simultaneously some equilibrium nuclear matter (binding energy per particle, saturation density, incompressibiliy, etc.) and finite nuclei properties. However, in this way, bulk and surface are not properly sepa-rated and early attempts used to miss important infinite nuclear matter properties, as, e.g. stability of neutron matter at high density [2] and other stability criteria. Modern Skyrme forces , like the Saclay-Lyon (SLy) ones, explicitly use the high density part (ρ/ρ 0 > 0.65) of microscopic neutron matter calculations for the EOS in the fitting procedure [4] and thus avoid collapse. Therefore, modern functionals usually reproduce reasonably well microscopically determined EOS for neutron and nuclear matter [29]. Examples are, among others, the SLy-forces [4] (see Fig.1 in [29]) and the Fayans functional (see e.g. Fig. 3 in [23]). It is thus evident that the procedure adopted in this work is different on a qualitative level from the usual and allows, via the fit, to reproduce very accurately the microscopic infinite matter results in the whole range of densities considered. This may be important for surface properties and neutron skins in exotic nuclei, what shall be investigated in the future. Let us also mention, in this context that for the electron systems mentioned before, the DFT practitioners demand to have extremly precise EOS as LDA input to their calculation.
For open shell nuclei, we still have to add pairing. The formal generalization of the rigorous HK theorem to paired system has been given in Ref. [35]. In the present work our main objective is to set up the KS-DFT scheme for the non pairing part, thus we add pairing in a very simple way within the BCS approach. For this we simply take the density dependent delta force defined in Ref [36] for m = m * with the same parameters and in particular with the same cutoff. As far as this amounts to a cutoff of ∼ 10 MeV into the continuum for finite nuclei, we have to deal with single-particle energy levels lying in the continuum. We have simulated it by taking in the pairing window all the quasi-bound levels, i.e. the levels retainded by the centrifugal (neutrons) and centrifugal plus Coulomb (protons) barriers. This tratment of the continuum works properly, at least for nuclei not far from the stability valley as it has been extensively shown in [37]. A more sophysticated procedure, as, e.g, the one proposed by Bulgac et al [24], may turn out to lead to better results. Since we do not want to concentrate on pairing here, this shall be investigated in forthcoming work.
In our calculations the two-body center of mass correction has been included in the self-consistent calculation using the pocket formula, based on the harmonic oscillator, derived in Ref. [38] which nicely reproduces the exact correction as it has been shown in [39]. Our functional is now fully defined and, henceforth, we call it BCP-functional.
In this exploratory investigation, we fitted two sets of parameters. We have considered only spherical nuclei which we chose as given below. The two fits were obtained in i) optimising ground state energies only, ii) optimising ground state energies together with the charge rms radii. The open parameters, V L , V U , r 0 of Eq.(4) as well as the spin-orbit strength W 0 , are fitted to reproduce the binding energies of the nuclei 16 O, 40 Ca, 48 Ca, 56 Ni, 78 Ni, 90 Zr, 116 Sn, 124 Sn, 132 Sn, 208 Pb and 214 Pb in the case of the parameter set BCP1 and additionally the charge rms radii (r c = r 2 p + 0.64 fm) of the nuclei 16 O, 40 Ca, 48 Ca, 90 Zr, 116 Sn, 124 Sn, 208 Pb and 214 Pb for the parameter set BCP2. Experimental binding energies and charge radii are taken from References [40] and [41] respectively. In table II we give the obtained parameter sets of fits (i) and (ii).
Our intention is to compare our results from the BCPfunctionals with the ones obtained from some of the most performant effective nucleon-nucleon forces available, namely D1S [42], NL3 [43] and SLy4 [4]. To this end, we have calculated the binding energies and charge radii of 161 even-even spherical nuclei (in line with the NL3 calculations reported in [44]). In Table III we report the energy and charge radii rms deviations between the corresponding experimental values and the theoretical ones obtained using the D1S force [45], the NL3 parametrization [44], the Skyrme interaction SLy4 [4] and our BCP1 and BCP2 functionals. We also display in Figs.2-3 the differences between the theoretical and experimental energies and charge radii in the range between 16 Ne and 224 U calculated with the BCP1 functional in comparison with the same quantities obtained from the D1S effective force [45].
In Table IV we display the two-neutron separation energies for some magic nuclei predicted by the BCP1 and BCP2 functionals as compared with the same results provided by D1S and NL3 as well as the corresponding experimental values.
We see that our theory is at least as performant as the other ones. In Table V we give the common bulk properties of both BCP-functionals. Incompressibility modulus K ∞ and symmetry energy J plus its derivative L have quite acceptable values. For the average gap in 116 Sn we obtain 1.5 MeV. For the drip line nucleus 176 Sn we get E/A=-6.57 MeV compared to E/A=-6.44 MeV (D1S) and E/A=-6.67 MeV (NL3).
In a recent analysis, Bertsch et al. [46] have pointed out that also in Skyrme forces implicitly only four or five parameters are relevant. For example, in the SLy forces [4], six out of the twelve parameters are strongly constrained by nuclear and neutron matter properties. However, the philosophy in the present work is different and, as we said before, is analogous to what is done in condensed matter physics. That is we use as an input an analytic representation of a very accurate ab initio calculation of the nuclear/neutron EOS in the whole relevant range of densities, feed this into the functional via LDA and fit the remaining four parameters, which characterize the surface and spin-orbit contributions to the energy density, to finite nuclei properties. It is remarkable and, to some extend unexpected, that this scheme works indeed very well. It seems to give some credibility to the applicability of the KS-DFT scheme, even for self bound systems, as treated here (see also a discussion of this point in [18]).   In conclusion, we have succesfully applied a Kohn-Sham Density Functional Theory inspired approach to ground state properties of nuclei, quite in analogy to what is usually done in e.g. condensed matter physics and for atoms and/or molecules [9,10], however, at variance with the procedure hitherto almost exclusively used in nuclear physics. This inspite of the fact that the applicability of KS-DFT for self-bound systems is, so far, an unsettled problem. The bulk part of the functional is given for all relevant densities and asymmetries, once and for all, from microscopic results using the converged holeline expansion based on realistic bare forces [29]. To this a phenomenological surface part with three adjustable parameters is added. A fourth parameter is the strength of the spin-orbit interaction for which we get values close to the usual ones in Gogny and Skyrme forces. The microscopic nuclear matter EOS has been fine tuned to pass through E/A = -16.00 MeV at saturation. This may be considered as the fifth adjusted parameter. In this pilot study we took for the pairing part of the functional the simplest possible procedure and adopted a previously adjusted density dependent δ-force [36]. The fit of open    parameters was done only for spherical nuclei. In view of the surprisingly small number of adjustable parameters, the results are very encouraging and well compete with the most performant mean field theories presently in use. In addition we have preliminary results which show that our BCP-functional also yields excellent results in the deformed case, again in very close agreement with those obtained with the Gogny D1S force [48]. Our study thus seems to support the applicability of the KS-DFT scheme also of self-bound systems. A more elaborate parameter search with different flexibilites than in (4) is an important task for the near future.