The Barcelona Catania Paris Madrid energy density functional

The foundation and applications of the Barcelona Catania Paris Madrid (BCPM) energy density functional are briefly reviewed. BCPM uses a paradigm more rooted on density functional theory and fits most of its parameters to sophisticated microscopic nuclear matter calculations. Finite nuclei are accounted for by introducing a direct finite range surface term, spin orbit potential, Coulomb interaction and pairing. Applications of this functional to the calculation of finite nuclei properties are presented as well as applications to the description of neutron star physics. The large number of applications discussed are possible because of the local character of the functional that simplifies enormously calculations in finite nuclei.


Introduction
There has been an extraordinary progress in the microscopic theory of nuclei [1][2][3][4][5][6][7][8][9], but the theoretical methods developed along the years are computationally quite demanding. Energy Density Functionals (EDF) have been devised to overcome the complexity of the full many-body treatment of nuclear systems. The method has been introduced not only to treat heavy nuclei, still out of reach of the microscopic theory, but also to simplify the numerical effort for an extensive study of large sets of nuclei throughout the mass chart. The price to be payed is the introduction of a set of phenomenological parameters, which however can have a well defined physical meaning and elucidate the physics underlying the general trend of the nuclear properties along the mass table. The method is quite flexible, and it allows the treatment of several different physical problems, hardly accessible to the microscopic theory.
The Landau theory of Fermi liquid can be considered one of the first method based on an EDF. In this case one assumes the existence of an energy functional which depends only on the quasi-particle occupation numbers, from which one can define both the single particle energies and the effective interaction between quasi-particles through functional derivatives. As it is well known, the effective interaction is then fixed by phenomenological parameters. Originally the theory was devised to describe phenomena near the Fermi surface, but within specific models it can be extended throughout the Fermi sphere and it can be then used to calculate the total energy [10].
For finite nuclei the introduction of a phenomenological effective force to calculate the nuclear binding has a long history. By construction these effective forces have to be used at the mean field level, i.e. Hartree-Fock or Hartree-Fock-Bogolyubov. Along the years a large number of effective forces has been devised. The majority of them are defined from a density dependent and zero range force, usually indicated as Skyrme forces. A survey of the Skyrme forces can be found in Refs. [11,12]. However there is another set of forces that introduces also finite range terms or are completely finite range, in particular the celebrated Gogny forces [13][14][15][16]. These forces are highly performing and have been widely used in nuclear structure not only in relation to the nuclear binding energy, but also in a variety of other applications, like the nuclear pairing, fission and also nuclear scattering. Less explored are other alternatives like the one of Ref. [17] using Yukawa form factors to take into account the physics of one pion exchange.
The use of effective forces poses in general some questions, that can be summarized as follows.
-Is it possible to derive the effective forces from the underline many-body theory which include the correlations ? -Is it possible to connect the effective force to the bare nucleon-nucleon interaction ? -To what extent an effective force can be used in a manybody framework ?
The first question can be considered essentially of no use, since a positive answer would be possible only after solving the full nuclear many-body problem, which can be too demanding and in any case it is what we want to avoid. A first principle answers to the other two questions relies on the Hohenberg and Kohn theorem [18], which states that the ground state energy of a system is a unique functional only of the density profile. This gives in principle the foundation for the use of a functional for determining the energy of the ground state. Unfortunately in general this first principle functional is not known. It is expected to be a complicated function of the densities at different positions. In condensed matter it was proposed by Kohn and Sham [19] to approximate the unknown functional by a more phenomenological one. This functional is the sum of the Hartree energy, calculated from the Coulomb electron-electron interaction, and a completely phenomenological term, which contains a set of parameters to be fitted on a set of systems. The latter is usually referred as "exchange and correlation" term and it is kept as universal as possible. The method was extensively used both in solid state and molecular physics with remarkable success. The extension of the method to the realm of Nuclear Physics is not straightforward, since the Hartree energy or Hartree-Fock (HF) energy is not a good starting point because of the strong repulsive character of the nucleon-nucleon interaction. The HF approximation with the bare nucleon-nucleon interaction is unreliable and the explicit inclusion of correla-tions is essential. A second major difference with respect to the condensate matter situation is the presence in nuclei of a sharp surface, which needs a special and separate treatment.
As it is apparent from the liquid drop model, in the simple expression for the binding energy of nuclei it is possible to separate, at least approximately, a bulk contribution and a surface contribution. In addition, it is known that the central part of nuclei is approximately of constant density. This suggests, in the possible nuclear functional, to relate the bulk part to the binding energy to the nuclear matter binding energy, and to separate the contribution of the surface. The former contains all the short range correlations, while the latter should take into account of the long range correlations, as well as of the modification of the interaction at the surface. The Barcelona-Catania-Paris-Madrid (BCPM) functional was devised following this scheme. In this review paper we summarise the theory and applications of the functional, and the future prospects.
The paper is organized as follows. In Sect. 3 the genesis of the BCPM functional is described as well as its main nuclear matter properties In Sect. 4 results of the application of BCPM to finite nuclei are presented and compared to those of other approaches as well as experimental data. In this section we also introduce a variant of BCPM with an effective mass different from the bare one that is referred to as BCPM * . In Sect. 5 BCPM is used to explore the physics of neutron starts including many of its different regions like the outer and inner crust or the liquid core. We end this review with a summary and some conclusions extracted from the results obtained with the functional.

The BCPM functional
We first illustrate briefly the Kohn and Sham (KS) [19] functional as it was originally formulated in condensed matter, in order to stress the difference and similarities with the nuclear case. Let us consider the simple case of a scalar interaction v between particles, the Coulomb interaction among electrons in the atomic case. The KS functional has the following form (1) where E H is the Hartree mean field term E H = 1 2 dr dr v(r, r )ρ(r)ρ(r ) (2) and v ext is a possible external potential (the central Coulomb potential in the atomic case). Notice that the density appearing in the Hartree term is just the electron density. The term E xc [ρ] is the so-called exchange and correlation term, which contains all the correlations beyond the Hartree approximation. According to HK theorem this term must exist in prin-ciple, but it is of course unknown. In most of the applications the following approximations are considered 1. A local density approximation, i.e. one assumes that the exchange and correlation term is a function that depends only on the local density. 2. The function that defines the exchange and correlation term can be taken from accurate calculations in homogeneous matter at the local density ρ(r). In the atomic case one can use the Monte-Carlo calculations for an electron gas, which can be considered exact. 3. One assumes that the density can be written as a sum over the contributions from a set of A orbitals φ i , being A the total number of particles Assumption 2 connects in general the functional with the homogeneous matter binding energy (the electron liquid in this case). Assumptions 1 and 2 are essentially equivalent to take for the interaction part of the functional the binding of the homogeneous system in local approximation, except for the long range Hartree term, which depend on the whole density profile.
Minimization of the functional with respect to the orbitals φ i produces Hartree-like equations where is an effective local potential. The single particle wave functions φ i (r) can be taken orthonormal, but it is important to notice that they cannot be interpreted as the Hartree or Hartree-Fock orbitals. In particular the single particle density matrix cannot be written as since otherwise this would imply that the ground state wave function is just a Slater determinant of the orbitals φ i (r), while the functional must include the correlations beyond the Hartree or Hartree-Fock mean field. They are just a convenient method to express the density. Improvement with respect to the local density approximation can be introduced by including gradient terms, that can be eventually extracted from calculations on slightly inhomogeneous infinite matter.
If we now consider the energy density functional for nuclei, we have to stress some relevant modifications. First of all for the treatment of the surface term the method of including gradient corrections is not viable in the nuclear case, since nuclei have a sharp surface, with a width that can be of the order of the inter-nucleon average distance. Other differences with respect to the atomic case must be introduced. First of all in nuclei there is no central potential, they are self-bound. This can pose problems in the formulation of the HK theorem, which however was generalized to selfbound systems in Refs. [20][21][22][23]. Another modification that is necessary comes from the complex structure of the nucleonnucleon interaction and its short range hard core. It is not possible in a reliable way to separate the nuclear Hartree term, and therefore it is joined together with the exchange and correlation term to form the energy density at the local density ρ(r), which again is taken from the an accurate many-body theory of nuclear matter. The Coulomb interaction is treated separately, as specified below. Finally, as already noticed, the sharpness of the nuclear surface suggests the introduction of a phenomenological surface term. The nuclear Kohn-Sham functional can then be written where T 0 is the uncorrelated kinetic energy, E bulk N is the bulk term of the interaction energy coming from the nuclear matter Equation of State (EOS) with (ρ(r)) the nuclear matter energy density at the local proton and neutron densities ρ p (r) and ρ n (r), respectively.
We discuss in some detail each contribution to the functional.

Kinetic energy
Following the KS method, an auxiliary set of A orthonormal wave functions φ i (r), where A is the mass number, is introduced to express formally the density as if it were obtained from a Slater determinant as a sum of the product of single particle wave functions with the φ's determined from the minimization of the ground state energy, under the constraint of a given number of proton Z and neutron N . As already stressed they cannot be interpreted as the Hartree or Hartree-Fock orbitals. However the uncorrelated kinetic energy can be expressed in the usual manner as where the sum runs over proton and neutron orbitals. Notice that here we are using the bare nucleon mass m, as in Ref. [31]. The introduction of an effective mass will be discussed later in Sect. 4.5.

The Coulomb and bulk contributions.
The Coulomb interaction energy among protons E C is a long range correlation and cannot be treated in the local density approximation and therefore it has to be treated separately. We consider this contribution at lowest order, i.e. the direct term E H C plus the exchange contribution E ex C in the Slater approximation and with E C = E H C + E ex C . The nuclear bulk interaction energy is taken from the nuclear matter EOS. In nuclei the matter is globally and locally asymmetric, and therefore one has to consider the EOS for asymmetric nuclear matter. As it is known, the dependence on the asymmetry parameter β can be taken to a good approximation as quadratic, so that one can write where P s (ρ) and P n (ρ) are the EOS for symmetric matter and pure neutron matter, respectively, and is the local asymmetry. The EOS has to be taken from microscopic calculations, which are not known to an accuracy comparable with the one of the electron liquid for the condensed matter case. However it has to fulfill some phenomenological constraints, to be discussed in the sequel. From the manybody theory of nuclear matter the EOS is necessarily obtained numerically, i.e. values of the energy per particle at a discrete set of densities. In order to transfer these numerical results to the functional it is necessary to introduce an analytic expression that interpolate at best these numerical data for both symmetric and pure neutron matter EOS. The specific form of these functions is not relevant and suggested only by their convenience. The simplest and easy to handle form is the The open circles represent the EOS of Ref. [33]. Figure taken from Ref. [26].
Reprinted with permission of the publisher polynomial one, and we found [26] that it is necessary to consider polynomial to fifth order, so that P s (ρ) and P n (ρ) read where ρ 0 and ρ 0n are some appropriate reference densities. The numerical results of the microscopic calculations and the corresponding fitting polynomials are displayed in Fig. 1.
For comparison the results of the variational calculations of Ref. [33] are also plotted. The many-body calculation of the nuclear EOS has a long history. The one shown in the figure was obtained within the Brueckner-Hartree-Fock approach, with the inclusion of three-body forces, as reported and discussed in e.g. [34][35][36]. The need for three-body forces is demanded by the finding that no microscopic many-body calculations with only two-body forces is able to reproduce an acceptable saturation point of symmetric nuclear matter, as extracted from the phenomenology on the binding energy and density profiles of nuclei [37]. Unfortunately three-body forces are poorly known and this introduces some uncertainty in the final results of the microscopic calculations. The way out of this problem is to fine tuning the saturation point of the fitted EOS in search of the best one in reproducing the nuclear binding along the mass table. We fixed the saturation density at ρ 0 = 0.16 fm −3 and tuned the saturation energy in the fitting procedure. We found that the results for the fitting of the total binding of nuclei is extremely sensitive to the chosen saturation energy of the EOS. This is not surprising, since a change of only 0.1 MeV in the saturation energy implies mainly a change of the order of 20 MeV in a nucleus as 208 Pb. Unfortunately no microscopic theory can be so accu- Table 1 Parameters (in MeV) of the polynomial fits, P s (ρ) = n a n (ρ/ρ 0 ) n and P n (ρ) = n b n (ρ/ρ 0n ) n for symmetric and neutron matter, respectively. The reference densities are ρ 0 = 0.16 fm −3 and ρ 0n = 0.155 fm −3 . Two sets of parameters for symmetric matter are given, they correspond to EOS with a minimum at the given values of E/A. The parameters are given with six digits, enough to obtain binding energies with an accuracy of a couple of eV. rate to fix the saturation energy with such a small uncertainty. We followed this optimization procedure in Ref. [26] and the optimal value turns out to be E/A = −15.98 MeV. The energy per particle E/A at saturation cannot be considered a mere free parameter, because it is tightly constrained within a small range of values, and therefore should not be considered on the same footing as the other fitting parameters of the functional. For completeness we report in Table 1 the values of coefficients of the different powers of the adopted polynomials For the symmetric matter case two sets of coefficients are given for two slightly different values of the energy per particle at saturation chosen in the fit procedure.
Besides the saturation point other physical quantities related to nuclear matter are constrained by phenomenology, and it is therefore relevant to compare the predictions of the adopted EOS with these constraints. Let us first remind the definitions of the different physical parameters usually considered in the literature.
The pressure P and density-dependent incompressibility K in asymmetric nuclear matter are defined in terms of the energy density in asymmetric nuclear matter H(ρ, β) as [38] and For symmetric matter (β = 0) at saturation density the K (ρ) just defined reduces to the well known incompressibility modulus K 0 = 9ρ 2 ∂ 2 (H/ρ) ∂ρ 2 | ρ=ρ 0 where ρ 0 is the saturation density. We will also consider the coefficient K , connected to the so-called skewness coefficient Q [39] by The isovector part of the nuclear interaction is characterized by the symmetry energy At saturation density one defines the symmetry energy coefficient as J = E sym (ρ 0 ) and two coefficients more, L and K sym , directly related to the first and second derivatives of Eq. (18) with respect to the density at saturation, respectively The values of the incompressibility modulus K 0 and the coefficient K defined in symmetric nuclear matter as well as the coefficients J , L and K sym are displayed in Table 2. The range of values of the incompressibility modulus K 0 is constrained by the experimental excitation energies of the isoscalar giant monopole resonance in finite nuclei since the pioneering works of Bohigas and collaborators [40] and Blaizot [41]. Unfortunately, different estimates of K 0 using different mean-field models predict slightly different values for this coefficient. The value K 0 = 230 ± 30 MeV has been proposed [42] as a compromise among the different available estimates. The value of the symmetry energy coefficient J , that dictates the isospin dependence of the nuclear interaction, is constrained by experimental data on heavy-ion collisions, pigmy dipole resonances and analog states. The range 30 ≤ J ≤ 35 MeV has been proposed for this coefficient [43]. The density dependence of the symmetry energy is relevant for many phenomena not only in terrestrial nuclei but also in neutron stars. Since the celebrated correlation established by Brown [44] between the slope of the symmetry energy and the neutron skin in 208 Pb, a considerable effort has been devoted to constrain the L parameter from available data. Antiprotonic atoms, nuclear masses, heavy-ion collisions, giant and pigmy dipole resonances, proton -nucleus scattering as well as theoretical calculations using a microscopic interactions have been used to extract the L coefficient (see Ref. [45] and references therein). From a compilation of the existing data the range of values L = 55 ± 25 MeV has been proposed [45,46]. However more recent data analysis on parity violating electron scattering [47] suggests the much higher value of L = 106 ± 37 MeV, which creates a tension with the majority of the other previous estimate. Our many-body calculations do not support such a high value of L.
The other two coefficients that characterize asymmetric nuclear matter around saturation, namely K and K sym are more uncertain. An estimate K = 700 ± 500 MeV has been proposed in Ref. [11]. The curvature of the symmetry energy, K sym , can be inferred from some non-relativistic [48] and relativistic [49] calculations in nuclear matter. From these results a range −200 ≤ K sym ≤ 150 MeV can be inferred. As can be seen from a direct inspection of Table 2, the nuclear matter properties of our BCPM energy density functional lie within the accepted ranges of values of the different nuclear matter quantities. Notice that the different symmetry parameters are interconnected, as explicitly shown in Ref. [50].

The surface term.
For the surface term we prefer a finite range contribution and for the form factors v q,q we took a simple Gaussian form v q,q (r ) = V q,q e −r 2 /r 2 0 (21) The indexes q, q stand for neutron and proton. The subtraction in Eq. (20) is done to avoid contamination with the bulk EOS. The strengths and range r 0 of the surface term are treated as fitting parameters. The number of parameters can be further reduced by imposing that expansion of the EOS in the bulk up to second order in the density is reproduced by the Gaussian term, which avoids the necessity of the subtraction term in Eq. (20). In this case the Gaussian form factors are taken with different ranges depending on the isospin channel r 0L and r 0U This introduces an additional parameter however the strengths V L and V U are now no longer considered as free parameters. They have been determined in such a way that the bulk limit of the surface term in Eq. (20), that is reproduces the ρ 2 term of the bulk part of the energy density in Eq. (12). Therefore, no subtraction as in Eq. (20) is required. This procedure imposes the following relationships between the surface term strengths and the coefficients a 1 and b 1 of the polynomial fits of symmetric and neutron matter, Here r 0L and r 0U are the ranges of the Gaussian in the surface term, as defined above. The numerical factorb 1 , defined as ρ 0n , has also been introduced to take into account the different reference densities for symmetric and neutron matter.
3.4 The spin-orbit term and the pairing part.
The spin-orbit in nuclei is quite strong and it is essential for the correct sequence and structure of the major shells. One can split the spin-orbit part into an uncorrelated part E s.o. plus a remainder not treated explicitly here. The form of the uncorrelated spin-orbit part is taken exactly as in typical Skyrme [51][52][53] or Gogny forces [13,54]. The spinorbit strength does not need to be adjusted because the final result is nearly independent of it. The spin-orbit strength is fixed by the requirement to reproduce the magic numbers and therefore depends on the major shells energy separation. The latter is roughly speaking inversely proportional to the effective mass. The same dependence with the effective mass should hold for the spin-orbit strength. Therefore, we will take the value W L S = 90.5 MeV which is consistent with BCPM's effective mass of one and the spin-orbit strength W L S = 130 MeV of the Gogny force (m * = 0.7m). Explicitly the spin-orbit term is a zero-range interaction v s.o. = i W 0 (σ i + σ j ) × [k × δ(r i − r j )k], whose contribution to the energy reads where the spin density for each kind of nucleon is obtained using the auxiliary set of orbitals as In the calculations of open-shell nuclei we also take into account pairing correlations. For this we simply take the density dependent delta force defined in Ref. [55] for m = m * with the same parameters and in particular with the same cutoff. In the original formulation of Ref. [24] we have performed calculations of open-shell nuclei taking into account pairing correlations at Hartree-Fock plus BCS level. 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 retained by the centrifugal (neutrons) and centrifugal plus Coulomb (protons) barriers. This treatment of the continuum works properly, at least for nuclei not far from the stability valley as it has been extensively shown in [56]. In this way we obtain two-neutron (S 2n ) and two-proton (S 2 p ) energy separations for magic proton and neutron numbers in quite good agreement with the experiment. In the implementation of pairing for finite nuclei (see below) we always use the HFB method and a constant pairing window of 10 MeV above zero.

BCPM for finite nuclei
One of the main goals of BCPM [26] is to define a reliable nuclear EDF to explore low energy nuclear structure in finite nuclei. By construction, BCPM reproduces extremely well state-of-the-art nuclear matter calculations with realistic nuclear forces and for this reason it could be even considered as an ab initio inspired EDF. In coincidence with the Skyrme family of functionals/interactions, one of the main advantages of BCPM is the lack of exchange terms in the mean field. This fact is responsible for a substantial speed up of calculations in finite nuclei as compared to other forces like the finite range Gogny force. This speed up has proven to be essential for some of the large scale applications of BCPM in finite nuclei.
As the definition of the BCPM functional in terms of nuclear matter did not include polarized matter, time odd densities are not included and therefore BCPM can only be used for even-even nuclei and zero spin. This represents a severe limitation of the functional, but recent results [57] concerning the good performance of the equal filling approximation [58] as compared with full blocking [59] results suggests a pathway to extend BCPM to the realm of odd-A nuclei. At present one of the authors (LMR) is working on an extension of BCPM for odd-A nuclei using the Equal filling approximation (EFA). The results of this work will be reported in the future. The 2-3 free parameters of BCPM were adjusted to reproduce experimental binding energies of the 579 even-even nuclei available at the time of the fitting (AME2003). The binding energies were computed in the framework of the Hartree-Fock-Bogoliubov method supplemented with some relevant beyond mean field correlation energies (see below).
For simplicity, the HFB calculations have been restricted to axially symmetric solutions. The quasiparticle operators of the Bogoliubov canonical transformation have been expanded in a harmonic oscillator basis {c + i , i = 1. . . . , N }. The size of the basis (the number of shells) varies according to the proton number of the nuclei considered. For nuclei with Z < 50 we take a basis with eleven harmonic oscillator shells, for 50 < Z < 82 we consider 13 shells and for Z >82 a 15 shell basis is used. These bases are used in a set of constrained calculations providing the potential energy as a function of the quadrupole moment for each nucleus. The minima of the potential energy are used as starting wave functions for an unrestricted minimization of the HFB energy. The wave function corresponding to the minimum energy solution for each nucleus, is used to generate a wave function in a larger basis (two extra shells) that is used to compute the corresponding HFB energy. This energy is used as input for an extrapolation algorithm to estimate the energy in an infinite basis [60] by means of the formula For the solution of the non-linear HFB equation, an approximate second order gradient method is used [61]. It provides a rapid and consistent procedure to reach the energy minimum. The information about the energy curvature taken into account in the method reduces the number of iterations substantially as compared to other approaches. As it is customary in this kind of calculations, the Coulomb exchange contribution is computed in the Slater approximation [62] and the Coulomb anti-pairing effect is not explicitly considered (see [63] for a discussion of this issue). The two body kinetic energy correction, which is typically considered as a way to correct for the lack of translational invariance of the whole procedure, has been taken into account with the pocket formula of Refs. [64,65] as in previous versions of the BCP functional. The beyond mean field correlation energy coming from symmetry restoration of rotational invariance (the rotational energy correction) plays a relevant role in deformed nuclei [66]. It rapidly changes with the deformation properties of the nucleus considered: it ranges from zero or very small values for magic or semi-magic nuclei, to values as large as 6 or 7 MeV for strongly deformed mid-shell heavy nuclei [67]. The rotational energy correction shows an almost linear dependence with quadrupole deformation and therefore its inclusion is relevant to describe masses all over the Nuclide Chart. The correct way to compute this quantity is by evaluating the angular momentum projected energy associated with each intrinsic state. This is a tremendous task that can fortunately be alleviated by considering an approximation to the projected energy that is obtained in the spirit of the Topological Gaussian Overlap Approximation. This fully microscopic formula, which is similar to the rotational energy correction, can be easily evaluated at the mean field level and does not involve any phenomenological parameter (see Ref. [67] for details). Finally, as already mentioned in previous sections, the pairing part is the density dependent zero range force of Ref. [55] that was devised to mimic the behavior of Gogny D1 pairing gap in nuclear matter. We have taken care of our effective mass equal to one by renormalizing the pairing force strength accordingly.
In the fit we have searched for the values of the parameters that minimize the rms deviation for the binding energies where the sum runs over the set of 579 even-even nuclei with known experimental binding energies, as given in the 2003 evaluation of Audi and Wapstra [68]. The 193 extrapolated binding energies included in Audi and Wapstra's work are excluded from our analysis, although we will use them to explore the quality of our fit. The three initial free parameters (i.e, the two ranges and spin-orbit) were at first considered in the fit but soon it became clear that the E/A value given by the polynomial fit should also be taken into account as a free parameter in the way we discussed above. Out of the four, it turns out that σ (E) value has a very smooth dependency on W L S and its minimum value was always obtained for spin orbit strength values around 90 MeV fm 5 . As explained before, this value is consistent with the value of 130 MeV fm 5 used in Gogny D1S and D1M parameterizations because of the different effective masses in BCPM (m * = m) and Gogny D1S and D1M (m * = 0.7m) that lead to a compression of the BCPM single particle spectrum requiring a reduced spin-orbit strength parameter as compared to Gogny (i.e. 0.7 × 130 ≈ 90). Another relevant observation is that the binding energy difference B = B th − B exp shows a linear dependence with mass number A (which is sometimes masked by large fluctuations at low A) with a slope that is intimately related to the value of E/A for the bulk. It is almost zero for E/A ≈ 15.97 and is clearly different from zero and positive for E/A = 16.03. It turns out that the value E/A = 15.98 yields the lowest σ (E) value. The final relevant observation is that σ (E) depends sensitively on the values of r 0L and r 0U and an accuracy of one part in 5 × 10 3 is required to obtain reasonable values for that quantity. Systematic explorations with a reduced set of spherical nuclei shows that there are two sets of r 0L and r 0U values that lead to reasonable values of σ (E), namely r 0L equal r 0U with values around 0.660 fm and r 0L taking values around 0.490 fm and r 0U around 1.050 fm. In the latter case, the value of r 0U is more critical than the r 0L value and it has to be kept at the value r 0U = 1.046 fm leaving only one free parameter to play, namely r 0L . Although the values of σ (E) (considering the 579 nuclei) obtained by minimizing with respect to r 0L and r 0U in the neighborhood of the two possibilities r 0L = r 0U ≈ 0.660 fm and r 0 L ≈ 0.490 fm  [26]. Reprinted with permission from the publisher and r 0U = 1.046 fm are similar, the first choice produces a B plot that looks smoother than the second and therefore we will from now on restrict ourselves to the first choice. We have considered r 0L = r 0U values in the interval between 0.650 fm and 0.670 fm in steps of 0.002 fm except for the interval bracketing the minimum where a step of 0.001 fm has been considered.
After all these considerations, one finds that the set of values E/A = 15.98 MeV, W L S = 90.5 MeV fm 5 , and r 0U = r 0L = 0.659 fm is the best choice. It leads to and outstanding σ (E) = 1.58 MeV for the set of 579 nuclei considered in the fit and defines the BCPM functional.
With the choice r 0U = r 0L = 0.659 one is fixing the surface term of the functional. A quantity that is often used to characterize the surface term is the surface energy that is estimated here by using the self-consistent Extended Thomas-Fermi approach includingh 2 corrections (ETF-h 2 ) [69]. The obtained value for BCPM is 17.68 MeV which is in consonance with other interactions.
In Fig. 2 we present the difference between the experimental and theoretical values for the binding energy B for the 579 nuclei of the AME2003 as a function of neutron number. The agreement with experimental data is much better for heavy nuclei than for light ones. A more quantitative assessment is obtained by looking at σ (E) for different sets of nuclei: taking into account nuclei with A > 40 (536 nuclei) the σ (E) for the energy gets reduced to 1. In the determination of the parameters of a functional it is not only important to find the values of the parameters minimizing the cost function σ (E) but also its variance. We find that BCPM is extremely sensitive to the surface term parameters r 0U = r 0L : a variation of 0.5 % in their value increases the σ (E) value from 1.58 to 2.16 MeV (40 %). This increment comes mostly from the A > 40 nuclei as the large fluctuations of B around zero remain roughly the same for A < 40 and also due to the A 2/3 dependence of the surface energy with mass number. This strong dependence on binding energies with some parameters of the functional is not new to BCPM, it also happens, for instance, with the Gogny functionals and most the Skyrme parametrizations due to the large strength of the t 3 parameter in front of the density dependent term. On the other hand, variation of a 5% in the spin-orbit strength leads to a 6% change in σ (E). We have also tested the sensitivity of the results to the degree of the polynomial used in the fit to nuclear matter. A large value of the degree might lead to unwanted oscillations whereas a small value might lead to a too simplistic functional form unable to grasp the essence of the nuclear matter calculations. The value used n = 5 provides the minimum value of σ (E) and increasing or decreasing the order of the polynomial by one unit increases σ (E) by 50 to 100 keV.
A more detailed look at the nuclei where B gets larger indicates that at the vicinity or at magic and/or double magic is where one observes the largest deviations for heavy nuclei. In particular the isotopes 208 Pb (Z = 82) and 130 Cd (Z = 48) are the ones with the largest B. One can also look at the contributions to σ (E) depending upon the deformation of the nuclei. Of the 579 nuclei considered in the evaluation of σ (E) there are 365 with a ground state deformation parameter β 2 greater, in absolute value, than 0.1 (a value that corresponds to a moderately deformed system). The σ (E) value for those deformed systems gets reduced to 1.08 MeV whereas the complementary one corresponding to near spherical systems (214 nuclei) grows up to a value of 2.19 MeV. This result is in the line of recent claims that well deformed systems are described better by mean field models than the spherical ones [70]. The result is also consistent with the previous finding concerning the maximum values of B for previous finding concerning the maximum values of B for magic or semimagic nuclei as those nuclei tend to have a nearly spherical ground state shape.
To finish this part we would like to compare our results with those obtained with the different parametrizations of the Gogny force, namely D1S [13], D1N [14] and D1M [15]. The results for σ (E) and the three parametrizations are given in Table 3. For the three parametrizations the values obtained with the raw HFB energies as well as those including the rotational energy correction (computed according to our procedure) are given. We also include in the table the results with and without rotational correction by performing an additional Table 3 The σ (E) values (MeV) for the three parametrizations of the Gogny force D1S [13], D1M [15] and D1N [14] and for different kinds of theoretical calculations. For the row marked as HFB the theoretical binding energy corresponds to the raw mean field energy, in the one denoted HFB+E ROT the rotational energy correction is also included. The last two rows correspond to the same theoretical set up as before but this time a global shift has been applied to the binding energy differences as to minimize the σ (E) value. The values of the shift parameter (in MeV) for each situation are given in parenthesis  [71] is used to compute the σ (E) value, it increases up to 1.61 MeV for 620 even-even nuclei in the experimental compilation. Clearly, the performance of BCPM for binding energies is not degraded when more nuclei are considered.
Nuclear binding energies play an important role in the description of atomic nuclei but there are other relevant observables characterizing the atomic nucleus and therefore any functional expected to become "universal" has to perform well also in the description of those quantities. We now discuss a couple of the most relevant ones.

Nuclear radii
The nuclear charge radius is a relevant observable connected to many other physical quantities. Experimentally, it can be accurately measured using electron scattering or other complementary techniques. Theoretically, already at the mean field level it is possible to obtain reasonable values for charge radii if the proton radius is taken into account. The charge radii for all even-even nuclei considered has been computed and compared to the experimental data published in Angeli's compilation [72] except for some Sr and Zr iso- Fig. 3 Differences r = r th − r exp between the computed radii and the experimental data taken from [72][73][74]. In panel a the results for the BCPM functional and in panel b the results for the Gogny D1S force. A proton's radius value of 0.875 fm has been considered [75]. Figure taken from [26]. Reprinted with permission from the publisher topes where laser spectroscopy results have been considered instead [73,74]. The theoretical predictions (computed with a proton's radius of 0.875 fm [75]) are confronted to the measured values of 313 even-even nuclei in panel a) of Fig. 3. A rms value of 0.027 fm is obtained, which is a quite reasonable number taking into account that the charge radius has not been considered in the fitting protocol. For comparison, the same quantity for the same nuclei computed with the Gogny D1S force (see panel b of Fig. 3) is 0.037 fm and 0.028 fm with the more recent parametrization D1M [15]. We notice that the largest contributions to the BCPM rms value come from the heaviest nuclei considered, where a strong deviation between theory and experiment is observed. Systematic deviations are also observed in nuclei around N = 40, N = 60 and N = 80, which are regions of the Nuclide chart characterized by the phenomenon of shape coexistence. As for the binding energies, the deviations strongly fluctuate in light nuclei with N < 40. The same pattern is also observed for the Gogny D1S results, suggesting that the origin for the discrepancies is more likely connected with a poor description of the ground state than with deficiencies of the interactions. A comparison of the two plots reveals that a better figure for σ R could be achieved in BCPM if an overall displacement would be performed, i.e. the theoretical values systematically underestimate the experimental values. It may be that size fluctuations (RPA correlations) bring the radii to their correct value as suggested in [76]. This is a task for the future.

Quadrupole and octupole deformations
In mid-shell nuclei, quadrupole deformation of the ground state is another relevant characteristic that can be linked to low energy nuclear properties like 2 + excitation energies or B(E2) transition probabilities along rotational bands. The connection of the intrinsic deformation parameter β 2 with experimental B(E2) transition probability values is somehow uncertain as it relies on the strong deformation limit of angular momentum projected theories to obtain the rotational formula [77]. For this reason we have preferred to compare BCPM's predictions with the ones of a well reputed, performing and predictive effective interaction. We have chosen the Gogny force, D1S [13], (results for the most recent published parametrization of the Gogny force, D1M [15], are very similar). In the upper panel of Fig. 4 a histogram is plotted depicting the number of nuclei with ground state quadrupole deformation parameters β 2 between β 2 (n) and β 2 (n) + δβ 2 . The discrete β 2 (n) values are given by the sequence β 2 (n) = δβ 2 n with n = . . . , −2, −1, 0, 1, 2, . . . The value δβ 2 = 0.0125 has been chosen in such a way that each bin represents roughly 5% of the typical value of 0.25 for the quadrupole deformation parameter. On the lower panel we plot a histogram with the difference β 2 (D1S) − β 2 (BCPM). We observe how most of the 818 nuclei considered have a difference in the β 2 parameter of less than 0.0125 in absolute value. A more detailed analysis of the results reveals that most of the discrepancies take place in the region Z ≈ N ≈ 40 which is widely recognized as a region of prolate-oblate shape coexistence and triaxiality. In those situations prolate and oblate (or triaxial) minima coexist with energy differences of just a few hundred keV. At the mean field level, the position of the absolute minimum is determining the ground state deformation. However, going beyond mean field by considering fluctuations on the quadrupole degree of freedom leads to collective wave functions exploring both coexisting minima and therefore those fluctuations have the potential to modify the theoretical predictions for radii.
Octupole deformation is associated to the multipole operator of order three and contrary to the quadrupole deformation case it breaks the parity symmetry. The presence of octupole deformation is not as common as the presence of quadrupole   [26]. Reprinted with permission from the publisher 56 and N = 88) both linked to the so called "octupole magic numbers". We have tested the octupole properties of BCPM by relaxing the parity symmetry in our self-consistent HFB calculations. In this way, the system can end up in an octupole deformed configuration in its quest for the minimum energy. The results show energy gains of the order of a few hundred keV and up to 1 MeV in nuclei in the aforementioned regions and both the energy gains and values of the β 3 deformation parameters are compatible with those obtained with the Gogny interaction [78]. Overall, the octupole correlation energies are about a few hundred keV smaller than the corresponding ones obtained with Gogny D1S [13] and D1M [15] ones  [28] with the only difference that in BCP1 the mean field octupole correlation energy was comparable to the one of Gogny D1S. However, as a consequence of the larger collective masses, the excitation energies predicted by BCP1 were substantially lower than the ones of Gogny D1S and the experiment. It remains to be tested what would be the spectroscopic predictions of BCPM concerning negative parity excited states but the lower mean field correlation energies go in the right direction.

Fission properties
Another testing ground for quadrupole and octupole properties of atomic nuclei is the shape evolution from the ground state to spontaneous fission. For this reason we have performed fission barrier calculations for a few selected examples and the results obtained for fission barrier heights and widths as well as mass asymmetry near scission (connected with fragment's mass asymmetry) are quite close to the results obtained with D1S, a functional that was specifically tailored to describe fission properties [13]. A detailed account of these calculations will be presented below and here we will focus on describing the results for a couple of examples: one is the paradigmatic case of 240 Pu whose fission properties have been computed with almost any proposed interaction. The other is the super-heavy 262 Sg where experimental data on spontaneous fission exist. In both cases, the results will be compared with the benchmark calculations carried out with the Gogny EDF. For an early account of fission barrier properties with previous versions of the BCP functional (BCP1 and BCP2) the reader is referred to Ref. [79].
We have followed the standard procedure that consists in the evaluation of the HFB energy as a function of the mass quadrupole moment Q 20 with the other multipole moments free to adopt any value in order to minimize the energy. Along this path to fission, we compute the collective inertias B(Q 20 ) required to evaluate the spontaneous fission half life in the WKB approximation (see [79] for details). The relevant quantities are shown in Figs. 5 and 6 for the nuclei 240 Pu and 262 Sg, respectively. In the lower panels, the HFB energy is depicted as a function of Q 20 for both BCPM and Gogny D1S cases. In the two considered nuclei the shape of the energy landscape looks rather similar regardless of the interaction considered. The maxima and minima are located roughly at the same positions and it is only the height of the barriers which changes with the interaction. Overall, Gogny D1S produces barriers higher than BCPM in agreement with the larger surface energy coefficient of D1S. The relevance of this characteristics of BCPM depends upon the reduction of the height of the first fission barrier when triaxial shapes are included in the calculation. In the case of the Gogny forces, the reduction can be as large as two or three MeV but we do not have a hint for BCPM yet. Work is in progress to explore We observe that in 240 Pu octupole deformation develops in its way to fission pointing to a dominant mass asymmetric fission mode. On the other hand, 262 Sg remains reflection symmetric along the whole fission path and therefore symmetric fission is expected to be the dominant mode in this nucleus. In panel b) of both figures the particle-particle pairing energies are given both for protons (dashed) and neutrons (full) for the two interactions. Again, the agreement in the evolution of these quantities with Q 20 for the two interactions considered is very good but BCPM produces smaller pp pairing energies suggesting weaker pairing correlations. The collective inertias are quantities very sensitive to the amount of pairing correlations as they strongly depend  [26]. Reprinted with permission from the publisher on the excitation energies of the low lying two quasiparticle excitations. In panel d) of Figs. 5 and 6 it is observed that the collective inertias of BCPM are between 2-3 times larger than the ones of Gogny D1S. The impact of the larger inertias (see [79] for a thorough discussion of this issue) on the spontaneous fission half lives is to make them longer. This effect is the opposite of having lower barrier heights in BCPM as compared to D1S. Therefore, the final values of the half lives will depend on the specific values of the quantities entering the "collective action". If the spontaneous fission half lives are computed in the WKB approximation we obtain for the 240 Pu case the values t 1/2 (sf) = 6.2 × 10 18 s for D1S and t 1/2 (sf) = 1.2 × 10 28 s for BCPM, to be compared to the experimental value of 3.5 × 10 18 s. For the 262 Sg case we obtain the values t 1/2 (sf) = 2.96s for D1S and t 1/2 (sf) = 4.2s for BCPM, to be compared to the experimental value of 6.9ms. From the above results we conclude that BCPM tends to yield larger values of t 1/2 (sf) than Gogny D1S and those values can be several orders of magnitude larger than in the Gogny D1S case. This is a direct consequence of the larger collective inertia and therefore a direct consequence of the reduced pairing correlations in BCPM. This results could lead us to the conclusion that BCPM should include stronger pairing correlations if the spontaneous half lives of fission are to be improved. However, we have to be aware of the impact on small changes in the parameters entering the WKB formula. For instance, the ground state energy usually incorporates some sort of zero point energy correction (denoted 0 in the literature) that strongly influences the WKB values of t 1/2 (sf). Its value has to be connected to the quantum effects associated to quadrupole motion but its exact meaning is still uncertain. As an example of the sensitivity of t 1/2 (sf) to this parameter, let us mention that if we use a value of 0 = 2.5 MeV instead of the value 0 = 1.5 MeV used originally we will obtain half lives which are 6 to 12 orders of magnitude smaller. In systematic fission calculations [80] the value of 0 is fine tuned for each isotopic chain considered and therefore absolute values for some choice of 0 as given here should not be considered too seriously for a comparison with experiment. Finally, the agreement between BCPM and Gogny D1S in the values of Q 30 along the fission path indicates that both interactions will produce, at least in a static framework, the same fission fragment mass distributions.
There are three parameters that characterize in a quantitative way the fission process of 240 Pu, one is the height of the first fission barrier E A , the second is the excitation energy of the fission isomer with respect to the ground state We notice a rather good agreement specially taking into account that no fission data has been considered in the fit. The BCPM value of E A is somewhat too high but this can be attributed to not having considered triaxiality in the calculation, an effect that can lower the energy by a couple of MeV. If the rotational energy correction is considered, the theoretical predictions get reduced to 7.04, 1.69 and 5.31 MeV respectively. These values are much closer to the experimental value than the HFB ones. The only value that seems to fall too low is E I I but the trend is consistent with recent estimations that reduce the experimental value to 2.25 MeV in 240 Pu [81].
This type of fission calculation was used [82] to consider an extended set of nuclei with experimental data on t 1/2 (sf). We found a strong dependence on the pairing strength parameters as well as on the value of a parameter introduced to model the ground state energy (E 0 ). In spite of these glitches of the theory, the agreement with experimental data is reasonable as can be confirmed by looking at Fig. 7. In this figure we also analyze the sensitivity of t sf to the strength of the pairing interaction by considering different multiplicative factor η.
In a subsequent publication [83] we considered the evaluation of the spontaneous fission half life using the least action scheme instead of the least energy one discussed above. Several collective coordinates in addition to the quadrupole  [82]. Reprinted with permission from the publisher moment were considered for the search of the least action path. It turns out that the relevant collective coordinate is the one modulating the amount of pairing correlations in the system. In our case, we used a constraint on the particle number fluctuation N 2 . For the several Uranium isotopes considered in the paper, the consideration of the least action framework with pairing degrees of freedom leads to a substantial reduction of t 1/2 (s f ) of roughly 20 orders of magnitude. The new predictions are in rather close agreement with experimental data and the discrepancies never exceed the five orders of magnitude (to be compared with the fifteen orders of magnitude difference obtained with the least energy framework). These results seem to indicate that the combination of BCPM and the least action framework with pairing degrees of freedom represents a rather accurate method to estimate theoretically spontaneous fission half-lives.

Systematic fission barrier calculations in neutron rich superheavy nuclei for astrophysical applications
As already mentioned, one of the advantages of BCPM over other functionals like Gogny is its local character that is responsible for the excellent computational performance of BCPM in finite nuclei. In addition, BCPM performs very well in reproducing binding energies, which is a key ingredient to estimate reaction rates as those required in nucleosynthesis. These are the main reasons why BCPM has been used to obtain nuclear data required in astrophysical applications [84,85]. In Ref. [84] BCPM was used to obtain fission barrier properties of around 3700 neutron rich heavy and superheavy nuclei as required for r-process nucleosynthesis incorporating the fission-recycling mechanism. Odd nuclei have been included in the analysis using the perturbative nucleon addition method (PNAM). Neutron separation energies in the set of around 3700 nuclei of interest are reproduced with a rms deviation of 0.30 MeV for even-even nuclei and 0.36 MeV including odd and odd-odd nuclei. The systematic of fission barrier heights and spontaneous fission half-lives has been thoroughly discussed in the paper. In addition, the αdecay half-lives have been obtained with the Viola-Seaborg [86] formula using the parameterization of Ref. [87]. The Viola-Seaborg formula makes use of the Q α value of the reaction and this quantity is reproduced with a rms deviation of only 0.68 MeV. The calculations of α-decay half-lives allows to determine which regions of the nuclear chart are dominated by fission or α-decay. In [85] the impact of fission on the production and destruction of trans-lead nuclei during the r-process nucleosynthesis occurring in neutron-star mergers is discussed. Abundance patterns and rates of nuclear energy production are obtained for different ejecta conditions using, among others, a sets of stellar reaction rates obtained with BCPM. Fission plays a fundamental role to determine free neutron abundance after r-process freeze-out. A significant impact on the abundances of heavy nuclei that undergo α decay or spontaneous fission is observed, affecting the energy production by the ejecta at timescales relevant for kilonova emission.

A BCPM functional with realistic effective mass
One of the most appealing properties of BCPM is its local character. As a consequence, the effective mass of the functional is just the bare mass. This is good in some cases (single particle levels around the Fermi level are not stretched) but the systematic of the energy of giant quadrupole resonances is not well reproduced by BCPM. To remedy this drawback of the model, we introduced a realistic effective mass by using a subtraction procedure where the quantum kinetic energy with the bare mass is replaced by the same quantity but this time with an effective mass plus a term that ensures that at the nuclear matter level one recovers the standard BCPM. The new functional thus obtained is referred to as BCPM * [30]. The effective masses for both protons are neutrons are fitted by polynomials of the density to results in nuclear matter obtained with realistic forces and the Brueckner procedure. The introduction of the effective mass in BCPM * does not spoil the σ (E) value obtained with BCPM for AME2012 [71]. The σ (E) slightly increased from the 1.61 MeV for BCPM to 1.65 MeV for BCPM * . On the other hand, the mean square deviation for radii gets substantially reduced from the 0.027 fm value of BCPM to 0.024 fm. Finally, in the fitting Table 4 Values of the free parameters of BCPM [26] and BCPM * as determined by minimizing the rms of the binding energy difference of even-even nuclei using the AME 2012 compilation [71]. The σ E value is given along with the corresponding value of the rms deviation of the radii. The results are obtained with basis 1 (see text) and only a weak dependence of the r 0U = r 0L = r 0 parameter is observed when increasing basis size (see text for details) to finite nuclei binding energies the values of the parameters change. Their value and a comparison with BCPM is given in Table 4. As a consequence of the introduction of the effective mass, the computational cost in finite nuclei increases with respect to BCPM. This increase is not substantial and it allows the use of BCPM * in many different scenarios like astrophysical applications or fission (see previous subsection). Regarding fission, it has to be mentioned that in the few examples explored, the fission barrier heights provided by BCPM * are smaller than the ones for BCPM, improving the agreement with the empirical values. A comparison of the results for fission barrier heights and excitation energy of fission isomers obtained with BCPM and BCPM * is summarized in Table 5. In the same table the known experimental data is also given, Regarding the giant resonance energies, the predictions for the centroid energy of the Giant Monopole Resonance do not change significantly as compared to the ones of BCPM. However, the energies of the Giant Quadrupole resonances obtained with BCPM * are increased by typically 1 MeV as compared to the BCPM predictions in better agreement with experimental data. The excitation energies of the monopole and quadrupole oscillations along the whole periodic table are displayed in Fig. 8. Both monopole and quadrupole energies follow a C A −1/3 law with coefficients C M = 86 and C Q = 67 MeV, respectively. These values are roughly in

Application of BCPM to Neutron Stars
One of the more important applications of the BCPM energy density functional is the derivation of the EoS for describing neutron stars (NS) with the same theoretical model from the surface to the center. The details about the different pieces of this EoS can be found in Ref. [89] and in a short version in [30]. Beneath a thin atmosphere, the NS interior consists of three main regions, the outer crust, the inner crust and the core, each one featuring a different physics [90,91], which cover a nuclear density range between 10 4 g/cm 3 and several times the symmetric nuclear matter saturation density (2.7×10 11 g/cm 3 ). There are few EoS devised and used to describe the whole NS within an unified framework. Let us mention among them the semi-phenomenological EoS of Lattimer and Swesty [92,93], the EoS of Shen et al. [94][95][96] computed with the Relativistic Mean Field formalism, the model of Douchin and Haensel [97] based on the Skyrme force SLy4 and the EoS derived by the Brussels group [98][99][100][101][102] using the BSk21 Skyrme force.
The core of a NS described with the BCPM energy density functional is assumed to be formed by a uniform mixture of neutrons, protons, electrons and muons in charge and beta equilibrium. In this region the contribution of their nuclear part to the EoS is obtained by means of the many-body calculation described in Sect. 3.2. In this region muons and electrons are treated as free relativistic Fermi gases. The physics of the crust is somehow different as far as we shall deal with finite-size nuclear clusters, which form a solid lattice to reduce the Coulomb repulsion, embedded by free neutron and electron gases in the inner crust and by only a free electron gas in the outer crust. We model the crust in the Wigner-Seitz (WS) approximation, which assumes that the space can be divided in non-interacting shell containing each one a single representative nuclear cluster permeated by the free neutron and electron (inner crust) and only electron (outer crust) gases.

The outer crust
The outer crust of a NS consists of a solid body-centered cubic lattice of fully ionized atomic nuclei permeated by a free electron gas in order to ensure charge and beta equilibrium. This region of the NS covers a density range between 10 4 g/cm 3 , where atoms are completely ionized, and 4×10 11 g/cm 3 , where neutrons start to drip from the nuclei. As a function of the density the nuclei in the Coulomb lattices change from 56 Fe at the beginning of the outer crust to neutron-rich nuclei at the bottom of the outer crust because is energetically more favorable for nuclei to reduce the proton number through electron captures with the energy excess carried out by neutrinos.
The formalism for describing the outer crust of NS was proposed first by Baym, Pethick, and Sutherland (BPS) [103] and applied later on in Refs. [104,105] in different studies of NS. It is assumed that the matter inside the star is cold and remains in the ground-state in complete thermodynamic equilibrium. The energy of the outer crust at given average baryon density n b consists of the nuclear plus free electrons and lattice contributions where M(A, Z ) is the rest mass of a nucleus with mass and atomic numbers A and Z , respectively. The term E elec is the kinetic energy of the free electrons, which are considered as a degenerated relativistic Fermi gas. The lattice contribution collects the different Coulomb effects, namely the repulsion among the nuclei of the lattice, their attraction with the electrons and the repulsion among the electrons. The lattice energy can be written as the Coulomb energy in the nuclear mass formula with a factor that takes into account the previously mentioned effects (see [105] for more details). The basic assumption in the calculation of the composition of the nuclei occupying the sites in the Coulomb lattice is that the thermal, hydrostatic and chemical equilibrium is reached at each layer of the crust. Only the free electrons and the lattice terms contribute to the pressure, Following [103], the composition of the nuclei in the lattice is obtained from the minimization of the Gibbs free energy per particle, which at zero temperature can be written as where μ e = k 2 Fe + m 2 e is the electron Fermi energy and C l the factor of the lattice energy. In order to find the optimal composition of the representative nucleus we proceed as follows. For a fixed pressure, Eq. (28) is minimized with respect to the mass A and atomic Z numbers of the nucleus. The nuclear masses M(A, Z ) needed to this end in Eq. (28) can be taken from experiment (AME2012 [71]) if they are available or calculated at HFB level discussed before using the BCPM energy density functional.

Composition and EoS of the outer crust
In Fig. 9 we plot the composition of the outer crust predicted by the microscopic BCPM model and by the macroscopicmicroscopic finite-range droplet model FRDM [106]. At very low densities up to n b ∼ 10 −6 fm −6 , the main contribution comes from Fe and Ni isotopes with neutron number N = 30, 34 and 36. At higher densities and up to n b ∼ 5 × 10 −5 fm −3 , Kr, Se, Ge, Zn and Ni isotopes with neutron number N = 50 appear. Above this density, where the experimental are unknown, the HFB calculation with the BCPM model predicts Mo, Zr, Sr, Kr and Se nuclei with neutron number N = 82 nuclear masses. The BCPM model reaches the neutron drip at a density and pressure of n b = 2.62 × 10 −4 fm −3 and P = 4.84 × 10 −4 MeV/fm 3 , respectively, with the composi- Fig. 9 Neutron (N ) and proton (Z ) numbers of the predicted nuclei in the outer crust of a neutron star using the experimental nuclear masses [71] when available and the BCPM energy density functional or the FRDM mass formula [106] for the unmeasured masses. Figure was taken from Ref. [89], reprinted with permission Fig. 10 The pressure in the outer crust against the baryon density using the experimental nuclear masses [71] when available and the BCPM energy density functional or the FRDM mass formula [106] for the unmeasured masses. Also shown is the pressure from models BSk21, BPS, Lattimer-Swesty (LS-Ska), and Shen et al. (Shen-TM1) (see text for details). The dashed vertical line indicates the approximate end of the experimentally constrained region. Figure was taken from Ref. [89], reprinted with permission tion given by the nucleus 116 Se. Our results are in very good agreement with the predictions of the FRDM, except in the width of the layer of 78 Ni before the N = 82 jump and in the nucleus at the bottom of the inner crust, which is predicted to be 118 Kr by the FRDM (see [89] for an enlarged discussion).
In Fig. 10 we display the EoS of the outer crust predicted by our calculation using experimental and theoretical BCPM masses. and from the experimental masses plus the FRDM.
We can see small jumps in the density for particular values of the pressure, which are associated with the change from an equilibrium nucleus to another in the composition [103]. During this change the pressure and the chemical potential remain constant, implying a finite variation of the baryon density [91,103,104,107]. After the region constrained by the experimental masses (marked by the dashed vertical line in Fig. 10), the pressures predicted by BCPM (black solid line) and other theoretical models such as FRDM [106], BPS, the Ska version of Lattimer-Swesty [92,93], the RMF Shen-TM1 EoS [96] and the Skyrme force BSK21 of the Brussels group [99] tabulated in [101]. We see that our BCPM EoS agrees well with the BPS, FRDM and BSk21 predictions, showing some discrepancies with the LS-Ska model near the bottom of the outer crust. Notice that the Shen-TM1 EoS for the outer crust behaves quite differently from the another models predicting a soft EoS. A possible reason for that ins the use a very schematic semiclassical proton and neutron densities in the Shen-TM1 calculation (see again [89] for additional discussions). We will not give here the explicit values of the BCPM EoS of the outer crust computed with the experimental masses and HFB calculations, which is reported in Table 4 of [89].

The inner crust
At an average density about 4×10 11 , the neutron star scenario changes respect to the one of the outer crust due to the presence of a free neutron gas in addition to the free electron gas permeating the Coulomb lattice. When the average baryon density grows, the fraction of free neutrons also increases. At the bottom layers of the inner crust, the equilibrium shapes may change from spherical to cylindrical or planar geometries in order to minimize the Coulomb repulsion. These shapes are generically known as "nuclear pasta". When the average density is about one-half of the saturation density in nuclear matter, the Coulomb lattice disappear and a new phase transition to a liquid phase of uniform nucleons and leptons takes place due to the energy balance.
Full HF or HFB calculations in the inner crust are very difficult owing to the neutron gas and eventually to the existence of non-spherical geometries corresponding to pasta shapes. Since the pioneering paper of Negele an Vautherin [108], there exist different model calculations at quantal level, as for instance [109][110][111][112], however, large scale calculations are usually performed using semiclassical techniques such as the Compressible Liquid Drop Model (CLDM) introduced by Baym, Bethe and Pethick [113] and used extensively in [92,97] or the Thomas-Fermi (TF) approximation with different degrees of sophistication (see [89] for an enlarged discussion and references). In this calculation we have obtained the EoS of the inner crust using the simple TF approach basically by two reasons. On the one hand, the EoS in this region is largely dominated by the neutron gas where it is expected that gradient corrections play a minor role. On the other hand, the TF approach is free of shell effects, which simplifies the study of non-spherical pasta phases.

Self-consistent Thomas-Fermi description of the inner crust of a neutron star
The total energy of an ensemble of A−Z neutrons, Z protons, and Z electrons in a spherical Wigner-Seitz (WS) cell of volume V c = 4π R 3 c /3 can be expressed as E n n , n p + m n n n + m p n p +E el (n e ) + E coul n p , n e + E ex n p , n e dr, where E n n , n p is the nuclear energy density without including the nucleon rest masses. In our approach it reads E n n , n p = 3 5 3π 2 2/3 2m n n 5/3 n (r) + 3 5 +V n n (r), n p (r) , (30) where the two first terms are the TF neutron and proton kinetic energy densities and V n n , n p is the interacting part provided by the BCPM nuclear energy density functional (cf. Sect. 3), which contains bulk and surface contributions. The term E el (n e ) in Eq. (29) is the relativistic kinetic energy density due to the motion of the electrons, including their rest mass. The term E coul n p , n e in Eq. (29) is the Coulomb energy density from the direct part of the proton-proton, electron-electron, and proton-electron interactions. Assuming that the electrons are uniformly distributed, this term can be written as with The exchange part of the proton-proton and electron-electron interactions are treated in the Slater approximation: We consider the matter within a WS cell of radius R c and perform a fully variational calculation of the total energy E(A, Z , R c ) of Eq. (29) imposing charge neutrality and an average baryon density n b = A/V c . Our treatment differs from some previous calculations of TF type carried out with non-relativistic nuclear models [99,[114][115][116], which use parameterized trial neutron and proton densities in the minimization.
Taking functional derivatives of Eq. (29) with respect to the particle densities and including the conditions of charge neutrality and constant average baryon density with suitable Lagrange multipliers, the variational equations for the neutron (n n ), proton (n p ) and electron (n e ) can be easily obtained (see [89,117] for the explicit expression of such equations). We solve these variational equations self-consistently for a given average baryon density n b in a WS cell of specified size R c following the method described in [118], which allows to determine the composition (A, Z ) in β-equilibrium that minimizes the energy given by Eq. (29). Next, we determine the optimal size of the WS cell by repeating the calculation for different values of R c , in order to find the absolute minimum of the energy per baryon for that n b . This can be a delicate numerical task because the minimum of the energy is usually rather flat as a function of the cell radius R c (or, equivalently, The same formalism developed for spherical nuclei can also be applied to obtain spherical bubbles (hollow spheres). The method of solving the variational equations in the inner crust is not restricted to spherical symmetry and it can be extended to WS cells with planar (slabs) or cylindrical (rods) symmetries. The TF with these non-spherical geometries are simplified if one considers slabs and rods of infinite extension in the perpendicular direction to the size R c . Although the number of particles and total energy of these cells are infinity, the number of particles and the energy per unit surface (slabs) or per unit length (rods) are finite and the same is true for the energy per baryon or per unit of volume. With this choice of the geometries the energy per unit surface or per unit length are easily determined from Eq. (29). Taking the volume element dV = Sdx (slabs) or dV = 2π Lrdr (rods) reduces the calculation of the total energy to a 1-dimensional integral over the finite size of the WS cell (from −R c to R c for slabs and from 0 to R c for rods. With these geometries the calculation of the Coulomb potentials and of the nuclear surface contribution (see Sect. 3.3) are largely simplified. The explicit expressions for the Coulomb and surface energies for slabs and rods can be found in Ref. [89]. As it happens for spherical bubbles, the equations corresponding to cylindrical rods can be applied to obtain cylindrical tubes (hollow rods). Optimal density profiles of neutrons and protons for droplets, rod, slab, spherical bubbles and cylindrical bubbles (tubes) at different baryon densities are displayed in Figs. 8 and 9 of [89].

Analysis of the self-consistent solutions
In left panel of Fig. 11 we display the difference between the energy per baryon computed in the crust with different geometries (spherical droplets, cylindrical rods, slabs, cylindrical tubes and spherical bubbles) and the energy per baryon in neutron-proton-electron (npe) matter in a average density domain between 0.005 fm −3 and 0.08 fm −3 . The spherical droplet configuration is the energetically most favorable shape up to and average density n b ∼ 0.067 fm −3 , where a first shape transition to rods takes place. This is the most favorable configuration till an average density n b ∼ 0.076 fm −3 where a second transition to a planar slab shape occurs. When the average density increases further, the energy per baryon of droplets and rods becomes progressively closer to that of slabs. To see more clearly this scenario of competing shapes, we zoom in the right panel of panel of Fig. 11 the energy difference between 0.07 fm −3 and 0.0825 fm −3 , which is the crust-core transition density predicted by the BCPM model. Above a density n b = 0.082 fm −3 there is a new shape transition to cylindrical tubes which in turn changes to spherical bubbles as the most favorable phase at almost 0.0825 fm −3 . It shall be point out that the energy differences at average densities close to the crust-core transition are actually very small, as it can be seen in Fig. 12, where the energy per baryon respect to the neutron rest mass as a function of  In Fig. 13 we display as a function of the average density the composition of the droplets, which is the shape energetically preferred except for densities close to the crust core transition where the energy differences are very small. We see that when neutrons start to drip in the inner crust the mass number increases quickly from A ∼ 110 to A ∼ 1100 at n b = 0.025 fm −3 . For higher densities, the mass number shows a smooth decreasing trend till A ∼ 800 at the transition point. The atomic number remains roughly constant in the inner crust with values ranging between Z ∼ 36 at n b = 0.01 fm −3 and Z ∼ 25 at n b = 0.07 fm −3 .

Equation of state of the inner crust
The pressure of the inner crust, which was derived in Appendix A of [89], reads where P ng is the pressure of the gas of dripped neutrons, P free el = n e k 2 Fe + m 2 e − E el (n e ) is the pressure of free electron gas, and P ex el = − 1 4 3 π 1/3 e 2 n 4/3 e is a corrective term from the electron Coulomb exchange. Notice, however, that the pressure obtained in a WS cell in the inner crust differs from the value in homogeneous npe matter in β-equilibrium at the same average density n b owing to the lattice effects, which influence the electron and neutron gases. The lattice effects take into account the presence of nuclear structures in the crust and are automatically included in the self-consistent TF calculation. The pressure exerted by free neutrons and electrons in the outer and inner crust is an essential ingredient in the Tolman-Oppenheimer-Volkoff equations [90], which determine the mass-radius relation in neutron stars. The crustal pressure has also significant implications for astrophysical phenomena such as pulsar glitches [119].
The pressure predicted by the BCPM functional in the inner crust is shown in Fig. 14 in comparison with the values provided by other models, namely, NV of Negele-Vautherin [108], Moskow calculation [111], BBP [103,113], LS-Ska [92,93], DH-SLy4 [97], Shen-TM1 [94][95][96] and BSk21 [99][100][101][102] The calculations of LS [92,93] and Shen et al. [94][95][96] are the two EoS tables in more widespread use for astrophys-  [89], reprinted with permission ical simulations. The initial baryon density in Fig. 14 corresponds to the last density shown in Fig. 10 when neutrons start to drip. In the inner crust, the pressure computed with the BCPM functional is comparable, in general, to the results of the NV, BBP, DH-SLy4, and BSk21 calculations. Particular agreement is observed in the inner crust regime between the BCPM and BSk21 pressures up to relatively high crustal densities. In contrast, large differences are found when the BCPM pressure is compared with the values provided by the Moskow, LS-Ska, and Shen-TM1 models, in particular at densities close to the crust-core transition, where the differences can be as large as a factor of two, which may have an influence on the predictions of the mass-radius relationship of neutron stars, particularly in small mass stars. In addition to the spherical shape, we have evaluated the pressure for the non-spherical WS cells and the hollow shapes in the regime of high average baryon densities using the BCPM functional. However, on the one hand, pasta phases appear as the preferred configuration only in narrow range of densities close the crust-core transition density and, on the other hand, the differences of pressure between the spherical and the pasta phases are small, generally of the order of 1-2 keV/fm 3 or less. Therefore, we have taken as a representative result for the whole inner crust the pressure calculated in the spherical droplet configurations. A more detailed discussion about the properties of the inner crust and the EoS in this region provided by the BCPM functional can be found in Ref. [89]

The liquid core
In order to study the structure of the NS core, we have to calculate the composition and the EoS of cold, neutrinofree, catalyzed matter. As we discussed in the Introduction, we consider a NS with a core of nucleonic matter without hyperons or other exotic particles. We require that it contains charge neutral matter consisting of neutrons, protons, and leptons (e − , μ − ) in beta equilibrium, and compute the EoS for charge neutral and beta-stable matter in the following standard way [90,113]. The energy density of lepton/baryon matter as a function of the different partial densities, E(n n , n p , n e , n μ ) = (n n m n +n p m p )+(n n +n p ) E A (n n , n p ) where we have used ultra-relativistic and relativistic approximations for the energy densities of electrons and muons [90], respectively. As explained in [24], the nuclear part of the energy per baryon can be written as a quadratic interpolation of the binding energy of symmetric nuclear matter and pure neutron matter in terms of the asymmetry (13) or, equivalently, the proton fraction x p = n p /(n n + n p ). Knowing the energy density Eq. (35), the various chemical potentials (of the species i = n, p, e, μ) can be computed  (Shen) straightforwardly, and the equations for beta-equilibrium, (b i and q i denoting baryon number and charge of species i) and charge neutrality, allow one to determine the equilibrium composition n i of each spice at given baryon density n b and finally the EoS, The EoS of the liquid core is displayed in Fig. 15 together with the predictions of the BSk21 [99][100][101][102], LS-Ska [92,93], DH-SLy4 [97] and Shen-TM1 [94][95][96] a We notice a remarkable similarity of our calculation with the DH-SLy4 EoS, while differs considerably, in particular at high densities, from the LS-Ska, Shen-TM1 and BSk21 EoSs. We recall that the pressures from the Lattimer-Swesty EoS and the Shen EoS had already been found to differ significantly from BCPM and SLy4 for the matter at subsaturation density in the inner crust (cf. discussion of Fig. 14). However, the BCPM and SLy4 pressures in the inner crust showed a concordance with BSk21 that remains within the core region up to about  Fig. 15) but is not maintained in the extrapolation to higher densities, where the BSk21 and Lattimer-Swesty models predict the stiffest EoSes of Fig. 15.
Once the EoS of the nuclear matter is known, one can solve the Tolman-Oppenheimer-Volkoff [90] equations for spherically symmetric NS: where G is the gravitational constant, P the pressure, the energy density, m the mass enclosed within a radius r , and r the (relativistic) radius coordinate. To close the equations we need the relation between pressure and energy density, P = P(n b ), i.e. the EoS. The integration of these equations produces the mass and radius of the star for given central density. It turns out that the mass of the NS has a maximum value as a function of radius (or central density), above which the star is unstable against collapse to a black hole. The value of the maximum mass depends on the nuclear EoS, so that the observation of a mass higher than the maximum mass allowed by a given EoS simply rules out that EoS. The considered EoSes are compatible with the largest mass observed up to now, which is close to 2.01±0.04 M [120], as it can be seen in Fig. 16, where the relation between mass and radius (left panel) and central density (right panel) is displayed. The observed trend is consistent with the equations of state displayed in Fig. 15. As expected, when the stiffness increases the NS central density decreases for a given mass. We also notice that the maximum mass calculated with the BCPM and the SLy4 EoSs is characterized by a radius of about 10 km, which is somewhat smaller than the radius calculated with the other considered EoSs. Recent analyses of observations on quiescent low-mass X-ray binaries [121] and X-ray burst [122] seem to point in this direction, though more studies could be needed [123]. For a NS of 1.5 solar masses, the BCPM EoS predicts a radius of 11.67 km (see Table 6), in line with the recent analysis shown in [124] and recent observations from the NICER project [125,126]. High-precision determinations of NS radii that may be achieved in the future by planned observatories such as the Large Observatory for X-ray Timing (LOFT), should prove a powerful complement to maximum masses for resolving the equation of state of the dense matter of compact stars.

Summary and conclusions
We have developed an energy density functional, dubbed BCPM, aimed to describe nuclear masses and other properties of finite nuclei and nuclear matter. This functional is built up in the framework of the Kohn-Sham method, which implies that the energy density can be split in an uncorrelated kinetic energy term plus a potential interacting part.
In the spirit of the nuclear mass formula, the nuclear interacting part is divided into the volume and surface contributions. The interacting bulk part is largely based on ab initio nuclear and neutron matter EoS, which are obtained through a Brueckner-Hartree-Fock calculation, based on the Brueckner-Bethe-Goldstone many-body theory, using the microscopic two-body Argonne v 18 interaction plus threebody forces computed with the Urbana model. The volume contribution in asymmetric nuclear matter is obtained by a quadratic interpolation in terms of the asymmetry parameter β between the EoS of nuclear and neutron matter, whose numerical values at a discrete set of densities are fitted once for ever by educated fifth degree polynomials in the density. In order to take into account surface properties, we add to the volume part a phenomenological contribution provided by the convolution of a Gaussian form factor with the neutron and proton densities. It should be pointed out that the strength of these terms are fixed from the second order coefficients of the polynomial fits of the nuclear and neutron matter EoS. To describe finite nuclei the nuclear part of the BCPM energy density functional is supplemented by the Coulomb term, consisting of the direct and exchange energies, the later computed using the local Slater approximation, and the uncorrelated spin-orbit part, which we choose zero-range and identical to the one used in Skyrme or Gogny effective forces. To describe open shell nuclei we include pairing in the BCPM model through a HFB calculation using zero-range density dependent interaction adjusted to reproduce the neutron gap predicted by the Gogny D1S force in nuclear matter. The number of free parameters of the BCPM energy density functional available for fitting finite nuclei properties are finally three, namely the binding energy per nucleon in symmetric nuclear matter, the range of the Gaussian form factor in the phenomenological surface term and the strength of the spin-orbit interaction. Once these three parameters are fully determined from the masses of 579 even-even spherical and deformed nuclei, one obtains an energy rms for this set of nuclei of 1.58 MeV and a rms of 0.027 fm for the charge radii of 313 nuclei experimentally known. These rms are similar, or still better, than the corresponding values obtained using well calibrated Skyrme or Gogny interactions for the same sets of nuclei. The BCPM energy density functional predicts quadrupole and octupole deformations and fission barrier similar to those provided by the Gogny D1S force, which can be considered as a theoretical benchmark of deformation properties in finite nuclei. Although surprising, it is very satisfying that two basic nuclear properties, the energy per particle in the bulk and the surface energy, suffice to completely determine the BCPM energy density functional predicting excellent ground-state properties of finite nuclei. We have also used the BCPM energy density to estimate the excitation energies of the isoscalar giant monopole and quadrupole resonances within the sum rule approach. The excitation energies of the giant monopole predicted by BCPM for nuclei from 90 Zr on are in reasonable good agreement with the experimental values. This is because the BCPM incompressibility modulus (212 MeV) lies within the empirical window of values extracted from the analysis of the experimental data for these resonances. However, BCPM overestimates the excitation energies of the giant quadrupole because of its effective mass being equal to the physical one. This behavior can be cured by introducing in the model a density and isospin dependent effective mass. This improved BCPM functional, called BCPM * , has an effective mass 0.79 in symmetric nuclear matter, reproduces better the experimental excitation energies of the quadrupole resonance in finite nuclei, and also improves the rms of the charge radii, which is now 0.024 fm. All these improvements do not worsen the rms for the binding energies An important application of the BCPM energy density functional is the derivation of an unified EoS for neutron stars from the surface to the center. As far as the density in a neutron star covers about eight orders of magnitude, there are different scenarios that require different theoretical descriptions. In the inner crust, where matter is arranged as atomic nuclei permeated by an electron gas, the nuclear masses, needed to obtain the EoS in this region, are provided by the experiment, or if they are unknown, by HFB calculations using the BCPM functional. The EoS in the inner crust, where we have in addition a neutron gas, is very complicated to compute following the rules of quantum mechanics and therefore we use instead the Thomas-Fermi approach with the BCPM energy density, which also allows to describe pasta phases. Although proton shell and pairing effects may be important to obtain the composition of this region, they are not very relevant to obtain the EoS. Finally, in the core region, which is the largest one, the EoS is obtained from the BHF calculations in bulk matter assuming charge and beta equilibrium. This EoS is used to obtain the mass-radius relationship by solving the relativistic Tolman-Oppenheimer-Volkov equation, which in the case of the BCPM functional predicts a maximum mass of 2.03 M , in agreement with the commonly accepted value of 2.01 ± 0.04M extracted from observations. Also the predicted radius of a standard neutron star of 1.4 M is 11,67 km, which is in line with the value extracted from the analysis of recent astrophysical observations.
Funding Information Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.

Data Availability Statement
This manuscript has no associated data or the data will not be deposited. [Authors' comment: No data sets were generated during the present study.] Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.