Deformation properties with a finite-range simple effective interaction

In a series of recent papers [1, 2] we have applied the finite-range simple effective interaction (SEI) [3, 4], which contains a single finite-range term having Gaussian form along with two δ-function terms, to study the ground state properties of spherical nuclei. At variance with typical effective interactions of Skyrme, Gogny and M3Y type, 10 of the 12 parameters of the SEI have been fitted by using experimental/empirical constraints in nuclear and neutron matter, as discussed in detail in [2]. Determination of 10 out of the total of 12 parameters of SEI in nuclear matter allowed us to reproduce the microscopic trends of the properties connected to the momentum dependence of the mean field and the equation of state (EOS) as predicted by Dirac–Brueckner–Hartree–Fock (DBHF), Brueckner–Hartree–Fock (BHF) and variational calculations using realistic interactions [5–16]. The two remaining parameters, one of them being the strength of the spin–orbit interaction, are left to reproduce the binding energies of finite nuclei. Pairing correlations required to describe open-shell nuclei are introduced with a density-dependent zero-range force of the type suggested by Bertsch and Esbenssen [17], which is widely used in nuclear structure calculations [18–25]. Within this framework it is found that SEI can reproduce the experimental values of the binding energies of 161 and the charge radii of 86 even–even spherical nuclei with similar quality to successful Skyrme and Gogny effective interactions or relativistic mean field models as shown in [2].

In our previous works [1, 2] we used a quasi-local energy density functional (QLEDF), obtained from an extended Thomas–Fermi approach of the density matrix [26, 27], plus a BCS treatment of pairing correlations to study the ground-state properties of spherical nuclei. Our aim in this work is to extend the description of finite nuclei with the SEI by performing full Hartree–Fock–Bogoliubov (HFB) calculations including deformation degrees of freedom. The paper is organized as follows. In section 2 the formalism concerning nuclear matter and finite nuclei using the SEI is briefly summarized. In the third section the predictions for the binding energies of 620 deformed and spherical even–even nuclei across the nuclear chart with known experimental masses [28], are made in the HFB framework for the different sets of the SEI interaction. The charge radii of 313 even–even nuclei predicted by these sets of SEI are also compared with the experimental data [29]. Comparison of the results of binding energies and charge radii for even–even spherical nuclei computed with the QLEDF used in [1, 2] with the corresponding HFB results of the present work is also given in this section. The fourth section is devoted to study the SEI predictions in some of the problems crucially dependent on deformation properties. In particular, we analyze the potential energy surfaces (PES) of some typical nuclei and the fission barriers of ²⁴⁰Pu and ²⁵⁶Fm. The predictions are compared with the experimental values where available and with the results obtained with the Gogny D1S force, which can be considered as a benchmark in the theoretical study of deformed nuclei and fission barriers [30]. Finally our conclusions are given in the last section.

The finite-range simple effective interaction used to describe nuclear matter and finite nuclei is given by

$$\begin{eqnarray}&&{v}_{{\rm{eff}}}({\bf{r}},{\bf{R}})={t}_{0}(1+{x}_{0}{P}_{\sigma })\delta ({\bf{r}})+\displaystyle \frac{{t}_{3}}{6}(1+{x}_{3}{P}_{\sigma }){\left(\displaystyle \frac{\rho ({\bf{R}})}{1+b\rho ({\bf{R}})}\right)}^{\gamma }\delta ({\bf{r}})\\ &&\qquad \qquad \quad +(W+{{BP}}_{\sigma }-{{HP}}_{\tau }-{{MP}}_{\sigma }{P}_{\tau }){{\rm{e}}}^{-{r}^{2}/{\alpha }^{2}},\end{eqnarray} \tag{ 1 }$$

where, ${\bf{r}}$ and ${\bf{R}}$ are the relative and center of mass coordinates, respectively. The SEI in equation (1) has 12 parameters in total, namely, $b,{t}_{0},{x}_{0},{t}_{3},{x}_{3},\gamma ,\alpha$, $W,B,H$ and M plus the spin–orbit strength parameter W₀, which enters in the description of finite nuclei. The SEI interaction is similar in form to the Skyrme force where the gradient terms are replaced by the single finite-range contribution. A similar analogy can be drawn with the Gogny interaction where one of the two finite-range terms is replaced by the zero-range t₀-term. One more difference in this context is that the density-dependent term of SEI contains the factor (1+ bρ)${}^{\gamma }$ in the denominator, where the parameter b is fixed to prevent supra-luminous behavior in nuclear matter at high densities [31]. The formulation of nuclear matter and neutron matter using SEI has been discussed at length in [1, 2], but for the sake of convenience of the reader we report in the appendices A and B the expressions of the energy densities of asymmetric nuclear matter, spin-polarized neutron matter and finite nuclei obtained with SEI given in equation (1). We shall now outline, in brief, the determination of the parameters involved in the study of nuclear and neutron matter. The study of asymmetric nuclear matter involves nine parameters altogether, namely, $\gamma ,b,{\varepsilon }_{0}^{l},{\varepsilon }_{0}^{{ul}}$, ${\varepsilon }_{\gamma }^{l}$, ${\varepsilon }_{\gamma }^{{ul}},{\varepsilon }_{{ex}}^{l}$, ${\varepsilon }_{{ex}}^{{ul}}$ and α, whose connection to the interaction parameters is given in the appendix A. However, symmetric nuclear matter requires only the following combinations of the strength parameters in the like and unlike channels

$$\begin{eqnarray}&&\left(\displaystyle \frac{{\varepsilon }_{0}^{l}+{\varepsilon }_{0}^{{ul}}}{2}\right)={\varepsilon }_{0},\left(\displaystyle \frac{{\varepsilon }_{\gamma }^{l}+{\varepsilon }_{\gamma }^{{ul}}}{2}\right)={\varepsilon }_{\gamma },\left(\displaystyle \frac{{\varepsilon }_{{ex}}^{l}+{\varepsilon }_{{ex}}^{{ul}}}{2}\right)={\varepsilon }_{{ex}},\end{eqnarray} \tag{ 2 }$$

together with $\gamma ,b$ and α, i.e., six parameters altogether.

For a given value of the exponent γ, and assuming the standard values for nucleon mass, saturation density and binding energy per particle at saturation, the remaining five parameters ${\varepsilon }_{0},{\varepsilon }_{\gamma }$, ${\varepsilon }_{{ex}},b$ and α of symmetric nuclear matter are determined in the following way. The range α and the exchange strength ${\varepsilon }_{{ex}}$ are determined simultaneously by adopting an optimization procedure using the condition that the nuclear mean field in symmetric nuclear matter at saturation density vanishes for a kinetic energy of the nucleon of 300 MeV, a result extracted from optical model analysis of nucleon–nucleus data [32–34]. The parameter b is determined as mentioned before. The two remaining parameters, namely ${\varepsilon }_{\gamma }$ and ${\varepsilon }_{0}$, are obtained from the saturation conditions. The stiffness parameter γ is kept as a free parameter and its allowed values are chosen in such a way that the corresponding pressure–density relation in symmetric matter lies within the region compatible with the analysis of flow data in heavy-ion collision experiments [35]. It is found that the maximum value that fulfills this condition is γ = 1, which corresponds to a nuclear matter incompressibility $K({\rho }_{0})$ = 283 MeV. Therefore, we can study the nuclear matter properties by assuming different values of γ up to the limit γ = 1. Now, to describe asymmetric nuclear matter we need to know how the strength parameters ${\varepsilon }_{{ex}},{\varepsilon }_{\gamma }$ and ${\varepsilon }_{0}$ of equation (2) split into the like and unlike components. The splitting of ${\varepsilon }_{{ex}}$ into ${\varepsilon }_{{ex}}^{l}$ and ${\varepsilon }_{{ex}}^{{ul}}$ is decided from the condition that the entropy density in pure neutron matter should not exceed that of the symmetric nuclear matter. This prescribes the critical value for the splitting of the exchange strength parameter to be ${\varepsilon }_{{ex}}^{l}=2{\varepsilon }_{{ex}}/3$ [36]. The splitting of the remaining two strength parameters ${\varepsilon }_{\gamma }$ and ${\varepsilon }_{0}$, is obtained from the values of the symmetry energy ${E}_{s}({\rho }_{0})$ and its derivative ${E}_{s}^{\prime }({\rho }_{0})$ = ${\rho }_{0}\frac{{{\rm{d}}E}_{s}({\rho }_{0})}{{\rm{d}}{\rho }_{0}}$ at saturation density ${\rho }_{0}$. By assuming a value for ${E}_{s}({\rho }_{0})$ within its accepted range [37], we determine ${E}_{s}^{\prime }({\rho }_{0})$ from the condition that the difference between the energy densities of the nucleonic part in charge neutral beta-stable $n+p+e+\mu$ matter and in symmetric matter at the same density be maximal [38]. The value of ${E}_{s}^{\prime }({\rho }_{0})$ thus obtained predicts a density dependence of the symmetry energy which is neither very stiff nor soft. With the parameters determined in this way, the SEI was able to reproduce the trends of the EOS and the momentum dependence of the mean field properties with similar quality as predicted by microscopic calculations [36, 38]. We have still two free parameters that we have taken to be t₀ and x₀. In our first work on finite nuclei in [1], the t₀ and x₀ parameters, along with the spin–orbit strength W₀, were determined by a simultaneous fit to the experimental binding energies of ⁴⁰Ca and ²⁰⁸Pb and to the splitting of the $1p$ single-particle levels in ¹⁶O. With the parameterizations corresponding to $\gamma =1/3(1/2)$, the SEI predicted rms deviations in 161 binding energies and 86 charge radii of even–even spherical nuclei of 1.39 (1.54) MeV and 0.017 (0.015) fm, respectively.

In our subsequent work in [2], we analyzed in detail the predictions of the SEI in the spin channel. It was found that the x₀ and W₀ parameters are, actually, correlated. Different combinations of x₀ and W₀ leave the rms deviations of binding energies and charge radii practically unchanged, but predict a very different behavior in spin-polarized matter which is sensitive to the value of x₀. To determine the x₀ parameter uniquely, we considered the particular case of spin-polarized neutron matter. Its description requires to know how the strength parameters ${\varepsilon }_{0}^{l},{\varepsilon }_{\gamma }^{l}$ and ${\varepsilon }_{{ex}}^{l}$ of spin-saturated neutron matter split into the two channels of like (parallel, 'l, l') and unlike (anti-parallel, 'l, ul') spin orientations under the constraints ${\varepsilon }_{0}^{l}$ = $({\varepsilon }_{0}^{l,l}+{\varepsilon }_{0}^{l,{ul}})$/2, ${\varepsilon }_{\gamma }^{l}$ = $({\varepsilon }_{\gamma }^{l,l}+{\varepsilon }_{\gamma }^{l,{ul}})$/2 and ${\varepsilon }_{{ex}}^{l}$ = $({\varepsilon }_{{ex}}^{l,l}+{\varepsilon }_{{ex}}^{l,{ul}})/2$. The completely polarized neutron matter computed with the SEI imposes that ${\varepsilon }_{\gamma }^{l,l}$ = 0 and ${\varepsilon }_{0}^{l,l}=-{\varepsilon }_{{ex}}^{l,l}$. Further, from a fit to the DBHF results on the effective mass splitting between spin-up and spin-down neutrons in spin polarized neutron matter, it is found that the SEI predictions agree well with the DBHF ones for ${\varepsilon }_{{ex}}^{l,l}={\varepsilon }_{{ex}}^{l}/3$. This consideration allows to determine x₀ in a unique way as [2],

$$\begin{eqnarray}&&{x}_{0}=1-\displaystyle \frac{2{\varepsilon }_{0}^{l}-{\varepsilon }_{{ex}}^{l}}{{\rho }_{0}{t}_{0}},\end{eqnarray} \tag{ 3 }$$

assuming that t₀ is known. The two remaining free parameters, t₀ and W₀, were fitted in [2] to reproduce the binding energies of ⁴⁰Ca and ²⁰⁸Pb, respectively. The values of the nine parameters of the SEI, which fully describe the asymmetric nuclear matter for $\gamma =1/3$ and $\gamma =1/2$, are presented in table 1. The nuclear matter saturation properties, such as, saturation density ${\rho }_{0}$, energy per particle $e({\rho }_{0})$, incompressibility $K({\rho }_{0})$, effective mass ${m}^{*}/m$, symmetry energy ${E}_{s}({\rho }_{0})$ and slope parameter of the symmetry energy $L({\rho }_{0})$ for these sets are listed in table 2. The values of t₀ and W₀, obtained from the fit to the binding energies of ⁴⁰Ca and ²⁰⁸Pb, along with the values of x₀ provided by equation (3) for the SEI with $\gamma =1/3$ and $\gamma =1/2$ are found to be the same as obtained in the HFB calculations of these sets and are given in the two first lines of table 3.

Table 1.  Values of the parameters of asymmetric nuclear matter for the two EOSs corresponding to γ = 1/3 and 1/2 used in this work.

γ	b	α	${\varepsilon }_{{ex}}$	${\varepsilon }_{{ex}}^{l}$	${\varepsilon }_{0}$	${{\varepsilon }_{0}}^{l}$	${\varepsilon }_{\gamma }$	${\varepsilon }_{\gamma }^{l}$
	fm3	fm	MeV	MeV	MeV	MeV	MeV	MeV
1/3	0.4184	0.7582	−95.6480	−63.7653	−112.7493	−67.0819	110.7436	78.7768
1/2	0.5914	0.7597	−94.4614	−62.9743	−78.7832	−45.8788	77.5068	57.7687

Table 2.  Values of nuclear matter properties at saturation for the two EOSs corresponding to γ = 1/3 and 1/2.

γ	${\rho }_{0}({\mathrm{fm}}^{-3})$	$e({\rho }_{0})(\mathrm{MeV})$	$K({\rho }_{0})(\mathrm{MeV})$	$\frac{{m}^{*}}{m}({\rho }_{0},{k}_{{f}_{0}})$	${E}_{s}({\rho }_{0})(\mathrm{MeV})$	$L({\rho }_{0})(\mathrm{MeV})$
1/3	0.1597	−16.0	226	0.710	35.0	74.94
1/2	0.1571	−16.0	245	0.711	35.0	76.26

Table 3.  Values of the parameters ${t}_{0},{x}_{0}$ and W₀, where x₀ is determined from equation (3), for the two EOSs corresponding to $\gamma \;=$ 1/3 and 1/2. The rms deviations with respect to experiment in the binding energies of 620 even–even nuclei and charge radii of 313 even–even nuclei are also shown. The column E_rot indicates what kind of rotational correction is included in the calculation.

	t0	x0	W0	$\sigma (E)$	$\sigma (R)$
γ	MeV fm3		MeV	MeV	fm	Erot
1/3	201	3.1931	115	1.873	0.0253	no
1/2	438	1.4182	112	1.958	0.0252	no
1/3	214	3.0595	115	1.788	0.0253	yes: VAP
1/2	450	1.4071	112	1.843	0.0252	yes: VAP
1/3	218	3.0221	115	1.742	0.0255	yes: RVAP
1/2	455	1.4027	112	1.561	0.0255	yes: RVAP

In the present work, we study deformed nuclei with the SEI and to this end we perform calculations in the mean field approach using the HFB method. These calculations are restricted to axially symmetric geometry. The pairing interaction is the same as in [1, 2]. The solution of the HFB equations has been recast as a minimization procedure of the energy density functional (see appendix B), where the HFB wave function of the Bogoliubov transformation is chosen to minimize the energy. The quasiparticle operators entering in this transformation are expanded in a harmonic oscillator basis with a number of shells depending on the nucleus considered. An approximated second-order gradient method [39] is used to solve the nonlinear HFB equations. Once the minimal energy solution is found, an extrapolation of the energy to an infinite number of shells [40] is performed (see also [41] in this respect).

3.1. Binding energies and charge radii of finite nuclei

The parameters t₀ and W₀ of the SEI are fixed to reproduce the experimental binding energies of known even-even finite nuclei. The mathematical procedure is to minimize the binding energy rms deviation $\sigma (E)$ defined as

$$\begin{eqnarray}&&{\sigma }^{2}(E)=\displaystyle \frac{1}{N}\displaystyle \sum _{i=1}^{N}{({B}_{{\rm{th}}}(i)-{B}_{{\rm{exp}}}(i))}^{2},\end{eqnarray} \tag{ 4 }$$

where the sum extends to the 620 even–even nuclei with known experimental masses [28]. For the theoretical binding energies we will consider three different possibilities. The first one is the value given by the HFB energy of the ground state minimum. The second possibility is the HFB binding energy supplemented by the rotational energy correction E_rot based on the HFB ground state. It corresponds to a projection after variation (PAV) procedure where the projected energy is computed using the rotational formula approximation and the HFB ground state wave function–see [42, 43] for details. Note that the rotational energy correction plays an important role in deformed nuclei and its inclusion is relevant to describe masses along the whole periodic table. In strongly deformed mid-shell heavy nuclei the rotational energy correction can reach values as large as 6 or 7 MeV. This correction, however, is almost negligible in magic or semi-magic nuclei, which are basically spherical. Due to the fact that the spherical-deformed transition is sharp, the rotational correction goes from zero to some MeV at the transition point leading to sharp variations in the binding energy plot. In order to find a remedy to this drawback we consider also a restricted variation after projection (RVAP), where the approximate projected energy is minimized in the space of HFB wave functions constrained to given values of the quadrupole moment. In this way, a relatively smooth behavior in the binding energy plot is obtained as can be seen from figures 1 and 2.

The minimization process reveals that in all the cases the optimal values of t₀ and W₀ are basically the same as the values fitted in [2]. These parameters, as well as the parameter x₀ (see equation (3)) together with the minimal value of $\sigma (E)$ for the SEI sets having $\gamma =1/3$ and $\gamma =1/2$, are reported in table 3 for the three different methods of calculating the binding energy. In figures 1 and 2 we plot the binding energy difference ${\rm{\Delta }}B={B}_{\mathrm{th}}-{B}_{\mathrm{exp}}$ for the $\gamma =1/2$ and $\gamma =1/3$ cases, respectively, as obtained with the three methods mentioned above. The results are shown for several isotopic chains ranging from Z = 8 to Z = 108. The three curves displayed in each isotopic chain correspond to the HFB energy (black curve), the HFB energy corrected with the rotational energy in the projection after variation (PAV) way (red curve) and the HFB energy corrected with the rotational energy but in the spirit of the restricted variation after projection (RVAP) approach (blue curve). The quantity ${\rm{\Delta }}B$ is plotted as a function of the neutron number N shifted by N₀ units, where N₀ indicates the origin of the horizontal axis for the different isotopic chains displayed in these figures. Perpendicular marks in the ${\rm{\Delta }}B=0$ horizontal lines indicate the position of neutron magic nuclei. Globally, these figures look qualitatively similar to the figure displayed in [41] computed using the BCPM energy density functional where the same color indexing was used (black for pure HFB, red for HFB with PAV and blue for HFB with RVAP). As mentioned, we have adjusted the parameters t₀ and W₀ to minimize $\sigma (E)$ in the three cases. It is found as a general feature that there is a degradation in the agreement with experimental data for light nuclei with Z values less than 50. For $Z\gt 50$ the behavior of ${\rm{\Delta }}B$ is almost flat, whereas below Z = 50 it behaves approximately as a straight line with a large slope. For light nuclei the curves show a more erratic behavior and the agreement between the theoretical predictions and the experimental values worsens (the differences between theory and experiment may be as large as ±4 MeV in some particular cases). This lack of accuracy in describing the experimental masses in light nuclei points out the fact that, in general, light nuclei are not very well described at the mean-field level (also see in this respect figure 5 of [41]). For the HFB result we can also mention the significant failure to reproduce isotopic chains with magic Z values. A typical example is Pb (and Po and Rn) where ${\rm{\Delta }}B$ values as large as 3.5 MeV are obtained for most of the nuclei in the chain. This is in strong contrast with other isotopic chains like the ones of Os and W where an almost perfect agreement with experiment is observed. Peaks in ${\rm{\Delta }}B$ are observed also for magic neutron numbers. The rotational energy correction (PAV) improves the situation around Pb, but the price to pay are the jumps observed in many isotopic chains around magic neutron numbers. These jumps are due to the sudden transition from spherical to deformed configurations, with the associated increase in E_rot from zero to a couple of MeV. The inclusion of the rotational correction improves a little the rms deviation $\sigma (E)$ between the theoretical and experimental binding energies, but not in a significant way. Finally, the RVAP correction smooths out the ${\rm{\Delta }}B$ curves and reduces the $\sigma (E)$ value by almost 300 keV in the $\gamma =1/2$ case. One might wonder at this point about other corrections to the HFB energy that could help to improve the $\sigma (E)$ value. Obvious candidates are the correlation energy associated to particle number restoration, or the zero point energy corrections of quadrupole and octupole motion. The octupole zero point energy was investigated in [44] for the Gogny force and the conclusion was that it played a minor role in improving ${\rm{\Delta }}B$. Work to evaluate the other two corrections is in progress.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We can compare the rms deviations $\sigma (E)$ obtained using the SEI with the predictions of Gogny forces for the same set of nuclei evaluated in the same conditions. At pure HFB level the $\sigma (E)$ values for the Gogny D1S, D1N and D1M are 3.48, 4.88 and 5.08 MeV, respectively. If the rotational energy correction is added, the binding energy rms deviations reduce to 2.15 (D1S), 2.84 (D1N) and 2.96 (D1M) MeV. However, if an additional global shift in the binding energy is added, one obtains $\sigma (E)$ values of 2.14 (D1S), 1.47 (D1N) and 1.45 (D1M). As explained in [41], this global shift has been included to simulate the zero point quadrupole energy correction included in the fitting protocol of D1M. Therefore, these reductions in the rms deviations in D1S, D1N and D1M are not surprising.

We have also explored the role of pairing correlations in the $\sigma (E)$ values by multiplying the pairing strengths for protons and neutrons by factors f_p and f_n, respectively, taking the values 0.95 and 1.05 (a 5% variation in the pairing strength). We have kept the same t₀ and W₀ values as in the ${f}_{p}={f}_{n}=1$ case, as an unrestricted search including four parameters would be computationally expensive. The results show that the ${f}_{p}={f}_{n}=1$ values provide the optimal value of $\sigma (E)$.

The nuclear charge radius is a relevant observable connected with the size of nuclei. It can be measured, for instance, through electron scattering experiments. In our calculation the charge radii are estimated from the point-like proton distributions taking into account a proton radius of 0.8 fm [45]. The theoretical predictions of SEI with $\gamma =1/3$ and $\gamma =1/2$ are compared to the measured values of the charge radii of 313 even–even nuclei reported in the Angeli's compilation [29]. In our previous study of spherical nuclei [2], where we considered a reduced set of 86 even–even nuclei, we obtained rms deviations $\sigma (R)$ in the charge radius of 0.017 fm and 0.016 fm for $\gamma =1/3$ and $\gamma =1/2$, respectively. In the present calculation, considering 313 spherical and deformed nuclei, we obtain $\sigma (R)$ = 0.0253 fm and 0.0252 fm for the two parameter sets, respectively. These values can be compared with the rms deviations in charge radii obtained for Gogny forces for the same set of nuclei, which are 0.037 fm for D1S and 0.028 fm for D1M. The differences ${\rm{\Delta }}r={r}_{\mathrm{ch}}-{r}_{\mathrm{exp}}$ are shown in figure 3 as a function of mass number and for the different isotopic chains. We observe a reasonable agreement with experiment with all the points scattered around the zero line except for light nuclei and some very heavy ones. The departure of our predictions for the very heavy nuclei is not in agreement with the overall behavior and deserves further consideration.

Zoom In Zoom Out Reset image size

3.2. Comparison with previous calculations of spherical nuclei using the SEI

The non-local density functional theory has been discussed in earlier literature (see e.g. [27] and references therein). As it was shown in [27], the Lieb theorem [46], which establishes the many-to-one mapping of the A particle Slater determinant wave-functions onto the local particle density, allows us to write the energy density functional in the non-local case as $\varepsilon [\tilde{{\rho }_{0}}]={\varepsilon }_{0}[\tilde{{\rho }_{0}}]+{E}_{{RC}}[\rho ]$. In this expression $\tilde{{\rho }_{0}}$ is a Slater determinant density matrix associated to an effective Hamiltonian $\tilde{H},{\varepsilon }_{0}[\tilde{{\rho }_{0}}]$ is the corresponding Hartree–Fock (HF) energy and the remaining part ${E}_{{RC}}[\rho ]$, which is a function of the local density ρ only, is the residual correlation energy. Effective forces of Gogny, M3Y and SEI type contain a finite range part plus a density dependent contribution. The HF energy density associated to these forces can be regarded as a possible realization of the exact energy density functional because the exact correlation energy density is in general unknown and it is approximated by the density dependent term of these finite-range effective forces. In [27] it is also shown that another reduction can be done by mapping the Slater determinant density matrix onto a set ${\rho }^{{QL}}$ of local particle, kinetic energy and spin densities, which allows to write the energy density in a quasi-local form $\varepsilon [{\rho }^{{QL}}]={\varepsilon }_{0}[{\rho }^{{QL}}]+{E}_{{RC}}[\rho ]$. In [1], as discussed in [27], we performed an additional approximation by computing the quasi-local energy density corresponding to the exchange term using the Extended Thomas–Fermi expansion of the density matrix [26]. This approximation is similar to the density matrix expansions proposed by Negele and Vautherin [47] and Campi and Bouyssy [48] (see [26] for further details). The variational principle applied to this quasi-local energy density functional allows us to write formally the equations of motion in a similar way to those obtained with the Skyrme forces [27]. As a consequence, at least for spherical nuclei, the calculations can be performed in coordinate space avoiding the expansion of the wave-functions in a basis.

Table 4 compares the binding energies and proton radii of some doubly closed shell nuclei computed with the QLEDF (see [1] for further details) with the results of the full HF method using the SEI with $\gamma =1/3$. From this table we see that the quasi-local energy density functional approach predicts values that are very close to the HF ones. The largest differences in binding energies and proton radii are always less than 0.5% and 0.3%, respectively. Pairing correlations have been included in density functional theory for a long time [49–51] and in the context of the quasi-local reduction of the non-local theory in [24]. Although HFB is the standard theory for dealing with pairing correlations in finite nuclei, the simpler BCS method is often sufficient to describe ground-state energies near the β-stability valley [52]. In [1, 2] we have used the QLEDF together with an improved BCS approach [53] to study the properties of finite nuclei. In table 5 we report the binding energies and proton radii for the Sn isotopic chain from N = 50 to N = 82 computed in this way using the SEI with $\gamma =1/3$ as well as the same quantities obtained using the HFB method with the same interaction. From this table we see again that the quasi-local results lie close to the HFB ones and that the largest differences, about 0.1%, correspond to the mid-shell nuclei. We have checked that similar results are obtained using the SEI with $\gamma =1/2$ instead of $\gamma =1/3$. From the discussion in this subsection we can conclude that the quantal predictions for spherical nuclei reported in [1, 2] by using the QLEDF including pairing correlations, treated with the improved BCS approach of [53], are in excellent agreement with the HFB results. A similar conclusion was reached in [24, 27] for Gogny forces.

Table 4.  Comparison of the quantal binding energies and radii obtained in doubly magic nuclei using the QLEDF of [1] and the HF approach.

Nucleus	EEDF(MeV)	EHF(MeV)	${r}_{p}^{{\rm{EDF}}}$(fm)	${r}_{p}^{{\rm{HF}}}$(fm)
16O	−127.6240	−127.0907	2.6594	2.6646
28O	−179.4020	−178.5816	2.7832	2.7922
40Ca	−342.1981	−341.2690	3.3997	3.4053
48Ca	−416.8068	−414.6866	3.4210	3.4326
56Ni	−478.7936	−480.0659	3.6847	3.6907
78Ni	−645.4516	−646.4948	3.8814	3.8928
100Sn	−824.8094	−825.6940	4.4132	4.4244
132Sn	−1105.0788	−1105.4751	4.6360	4.6456
208Pb	−1636.6551	−1635.8961	5.4285	5.4344

Table 5.  Comparison of quantal binding energies and radii obtained in the isotopes of Sn using the QLEDF plus improved BCS pairing of [1] and the HFB method.

Nucleus	EEDF(MeV)	EHFB(MeV)	${r}_{p}^{{\rm{EDF}}}$(fm)	${r}_{p}^{{\rm{HFB}}}$(fm)
102Sn	−848.1597	−848.8267	4.4297	4.4389
104Sn	−870.2672	−870.8571	4.4450	4.4535
106Sn	−891.6288	−891.9747	4.4603	4.4685
108Sn	−912.2156	−912.2815	4.4758	4.4839
110Sn	−932.3108	−931.8436	4.4915	4.4995
112Sn	−951.4225	−950.7060	4.5069	4.5150
114Sn	−969.8919	−968.9010	4.5215	4.5301
116Sn	−987.5958	−986.4534	4.5355	4.5447
118Sn	−1004.5842	−1003.3858	4.5488	4.5586
120Sn	−1020.8849	−1019.7195	4.5620	4.5721
122Sn	−1036.5196	−1035.4725	4.5748	4.5850
124Sn	−1051.5581	−1050.6561	4.5875	4.5977
126Sn	−1065.9662	−1065.2716	4.6021	4.6100
128Sn	−1079.7420	−1079.3066	4.6137	4.6221
130Sn	−1092.7420	−1091.7298	4.6251	4.6341

4.1. Potential energy surfaces

The potential energy surface (PES), defined as the HFB energy computed with HFB wave functions constrained to definite values of the axial quadrupole moment, are plotted in figure 4 for selected nuclei from light to very heavy ones. The Gogny D1S curves are given as a reference. The main conclusion is that the position of the maxima and minima obtained with the two versions of SEI and Gogny D1S are located essentially at the same values of the quadrupole deformation. However, the spherical barrier (the energy difference between the spherical maximum and the prolate ground state) is much lower in energy for the SEI interactions than for the Gogny force. Also the excitation energy of fission isomers or super-deformed states is slightly different depending on the forces used. This energy difference is strongly influenced by the position of the spherical single particle levels, which depend upon the average mean field and the spin–orbit interaction. Finally, the barriers separating the ground state minimum from the fission isomer are lower in SEI than in Gogny D1S. In order to quantify the impact on observables due to the PES differences, a dynamical calculation of the generator coordinate method (GCM) type would be required. This is the subject of active research.

Zoom In Zoom Out Reset image size

4.2. Fission barriers

Fission is one of the most characteristic aspects of nuclear dynamics. The description of the dynamical evolution of the parent nucleus into two fragments is a real challenge both for the quantum many body techniques to be used as well as for the effective interactions. Therefore, it is a very good testing ground for any newly proposed interaction. To test the performance of SEI in this matter, we have followed the traditional approach [54] and carried out constrained mean field calculations with the mass quadrupole moment used as the driving coordinate. The resulting PES obtained for the nuclei ²⁴⁰Pu and ²⁵⁶Fm are shown, as an example, in figures 5 and 6, respectively, where the HFB energy supplemented by the rotational correction, is plotted as a function of the quadrupole moment for the two sets of EOSs of SEI considered in this paper, as well as for the Gogny D1S force. The Gogny D1S results are used here as a benchmark, given the fact that the force was adjusted to fission barrier heights, and has been successfully tested with very different kinds of fission—see [55] as an example. By looking at the PES in figure 5 we notice that the quadrupole moment of all the maxima and minima present in the PES is roughly the same in the three cases, but the energies, relative to the ground state minimum, change depending on the interaction considered. For ²⁴⁰Pu, the SEI values for the first fission barrier height (energy difference between the first maximum at around ${Q}_{2}=25$ fm² and the ground state minimum) are smaller than the one obtained with Gogny D1S. However, the SEI values for the first fission barrier height are very close in both cases to the accepted experimental value of 6.05 MeV [56]. The trend is similar for the second barrier, located at around ${Q}_{2}=65$ fm², with heights (measured with respect to the ground state energy) in the SEI case lower than the height obtained for Gogny D1S. The SEI values are a few hundred keV lower than the experimental value of 5.15 MeV [56] of the second barrier. Finally, the excitation energy of the fission isomer (local minimum at around ${Q}_{2}=45$ fm²) obtained with both SEI forces is lower than the Gogny prediction and around 1.5 MeV lower than the experimental value of 2.8 MeV. In spite of these small discrepancies we can consider that the fission properties of SEI for ²⁴⁰Pu are in reasonable agreement with the ones of the Gogny force as well as with the experimental data. In the case of ²⁵⁶Fm shown in figure 6, there is no experimental data with which to compare. The agreement with the Gogny D1S curve is remarkable until the end of the first barrier but from there on, it deviates a little bit as a consequence of different shell effects. In both cases, a dynamical treatment including the quadrupole collective inertia and zero point energy corrections would be required for the calculation of the spontaneous fission life-time. Work along this direction is in progress and will be reported in the future.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We have performed a study of finite nuclei using the SEI including deformation degrees of freedom. This study extends the results for spherical nuclei reported in previous literature [1, 2]. In this work we have computed ground state properties of deformed nuclei using the Hartree–Fock–Bogoliubov method with a harmonic oscillator basis. This calculation of deformed nuclei includes the rotational energy correction, computed using the restricted variation after projection method. The SEI, with the parameters fitted as explained in this work, is able to describe the binding energies of 620 and the charge radii of 313 even–even nuclei with rms deviations about 1.8 MeV and 0.025 fm, respectively. These deviations are similar to those found, for the same set of nuclei and computed in the same conditions, using the D1S, D1N and D1M Gogny forces.

Concerning deformation properties predicted by the SEI, it is found that the potential energy surfaces along the whole periodic table follow quite well the pattern obtained using the D1S Gogny force. Nevertheless, some differences appear in the comparison between the SEI and the D1S results. In particular, the spherical maxima and super-deformed minima predicted in the potential energy surfaces with the SEI models are systematically lower than the corresponding ones obtained with the Gogny D1S force. This suggests that the surface energy parameter of both SEI models is smaller than the one of the D1S force. We have also studied the fission barriers of the nuclei ²⁴⁰Pu and ²⁵⁶Fm with the SEI. For the nucleus ²⁴⁰Pu this model predicts around (slightly depending on γ) 6.0 MeV for the first fission barrier, around 4.5 MeV for the second fission barrier and an excitation energy of the fission isomer of around 1.5 MeV. These values agree quite well with the accepted experimental values [56].

The present work confirms, as suggested by our previous analysis in [1, 2], the ability of the SEI to reproduce properties of different types of nuclear matter and, simultaneously, predict finite nuclei properties with similar quality to well-known effective interactions. At variance with other commonly used effective forces, such as Skyrme, Gogny or M3Y, 10 of the 12 parameters of the SEI are determined using experimental/empirical constraints in nuclear matter as well as to reproduce the behavior predicted by microscopic calculations in symmetric nuclear matter and spin-polarized neutron matter. In particular, the momentum dependence of the mean field derived from the interaction is fitted in such a way that the effective mass splitting in both, asymmetric nuclear matter and spin-polarized neutron matter, reproduce the results provided by microscopic Dirac–Brueckner–Hartree–Fock calculations. We note that only the parameter t₀ of SEI and the strength of the spin-orbit interaction W₀ are fitted to finite nuclei. Therefore, the results for finite nuclei properties obtained with the SEI are, actually, predictions of the model. We have focused our study on even–even nuclei. To extend the study to odd–even, even–odd and odd–odd nuclei requires the development of completely polarized asymmetric nuclear matter in order to deal with odd components in the energy density. Work in this direction is left for a future contribution.

One author (TRR) thanks the Departament d'Estructura i Constituents de la Matèria, University of Barcelona, Spain for hospitality during the visit. The work is covered under FIST programme of School of Physics, Sambalpur University, India. Authors MC and XV acknowledge partial support from Grant No. FIS2014-54672-P from Spanish MINECO and FEDER, Grant No. 2014SGR-401 from Generalitat de Catalunya, the Consolider-Ingenio Programme CPAN CSD2007-00042, and the project MDM-2014-0369 of ICCUB (Unidad de Excelencia María de Maeztu) from MINECO.

As mentioned in section 2, asymmetric nuclear matter computed with the SEI in equation (1) depends on the parameters $\gamma ,b$, ${\varepsilon }_{0}^{l},{\varepsilon }_{0}^{{ul}}$, ${\varepsilon }_{\gamma }^{l}$, ${\varepsilon }_{\gamma }^{{ul}},{\varepsilon }_{{ex}}^{l}$, ${\varepsilon }_{{ex}}^{{ul}}$ and α. The new parameters are connected to the parameters of the interaction through the relations,

$$\begin{eqnarray}{\varepsilon }_{0}^{l} & = & {\rho }_{0}\left[\displaystyle \frac{{t}_{0}}{2}(1-{x}_{0})+\left(W+\displaystyle \frac{B}{2}-H-\displaystyle \frac{M}{2}\right){\pi }^{3/2}{\alpha }^{3}\right]\\ {\varepsilon }_{0}^{{ul}} & = & {\rho }_{0}\left[\displaystyle \frac{{t}_{0}}{2}(2+{x}_{0})+\left(W+\displaystyle \frac{B}{2}\right){\pi }^{3/2}{\alpha }^{3}\right]\\ {\varepsilon }_{\gamma }^{l} & = & \displaystyle \frac{{t}_{3}}{12}{\rho }_{0}^{\gamma +1}(1-{x}_{3}),{\varepsilon }_{\gamma }^{{ul}}=\displaystyle \frac{{t}_{3}}{12}{\rho }_{0}^{\gamma +1}(2+{x}_{3})\\ {\varepsilon }_{{ex}}^{l} & = & {\rho }_{0}\left(M+\displaystyle \frac{H}{2}-B-\displaystyle \frac{W}{2}\right){\pi }^{3/2}{\alpha }^{3}\\ {\varepsilon }_{{ex}}^{{ul}} & = & {\rho }_{0}\left(M+\displaystyle \frac{H}{2}\right){\pi }^{3/2}{\alpha }^{3}.\end{eqnarray} \tag{ A1 }$$

The energy density in asymmetric nuclear matter reads

$$\begin{eqnarray}H({\rho }_{n},{\rho }_{p}) & = & \displaystyle \frac{3{\ℏ}^{2}}{10m}({k}_{n}^{2}{\rho }_{n}+{k}_{p}^{2}{\rho }_{p})+\displaystyle \frac{{\mathrm{\ϵ}}_{0}^{l}}{2{\rho }_{0}}({\rho }_{n}^{2}+{\rho }_{p}^{2})+\displaystyle \frac{{\mathrm{\ϵ}}_{0}^{{ul}}}{{\rho }_{0}}{\rho }_{n}{\rho }_{p}\\ & & +\left[\displaystyle \frac{{\mathrm{\ϵ}}_{\gamma }^{l}}{2{\rho }_{0}^{\gamma +1}}({\rho }_{n}^{2}+{\rho }_{p}^{2})+\displaystyle \frac{{\mathrm{\ϵ}}_{\gamma }^{{ul}}}{{\rho }_{0}^{\gamma +1}}{\rho }_{n}{\rho }_{p}\right]{\left(\displaystyle \frac{\rho ({\bf{R}})}{1+b\rho ({\bf{R}})}\right)}^{\gamma }\\ & & +\displaystyle \frac{{\mathrm{\ϵ}}_{{ex}}^{l}}{2{\rho }_{0}}({\rho }_{n}^{2}J({k}_{n})+{\rho }_{p}^{2}J({k}_{p}))\\ & & +\displaystyle \frac{{\mathrm{\ϵ}}_{{ex}}^{{ul}}}{2{\rho }_{0}}\displaystyle \frac{1}{{\pi }^{2}}\left[{\rho }_{n}{\displaystyle \int }_{0}^{{k}_{p}}I(k,{k}_{n}){k}^{2}{\rm{d}}{\⁢}{k}+{\rho }_{p}{\displaystyle \int }_{0}^{{k}_{n}}I(k,{k}_{p}){k}^{2}{\rm{d}}{\⁢}{k}\right]\end{eqnarray} \tag{ A2 }$$

where, the functions $J({k}_{i})$ and $I(k,{k}_{i})$ with ${k}_{i}={(3{\pi }^{2}{\rho }_{i})}^{1/3}$ ($i=n,p$) are given by

$$\begin{eqnarray}&&J({k}_{i})=\displaystyle \frac{3{{\rm{\Lambda }}}^{3}}{2{k}_{i}^{3}}\left[\displaystyle \frac{{{\rm{\Lambda }}}^{3}}{8{k}_{i}^{3}}-\displaystyle \frac{3{\rm{\Lambda }}}{4{k}_{i}}-\left(\displaystyle \frac{{{\rm{\Lambda }}}^{3}}{8{k}_{i}^{3}}-\displaystyle \frac{{\rm{\Lambda }}}{4{k}_{i}}\right){{\rm{e}}}^{-4{k}_{i}^{2}/{{\rm{\Lambda }}}^{2}}+\displaystyle \frac{\sqrt{\pi }}{2}{\rm{erf}}(2{k}_{i}/{\rm{\Lambda }})\right]\end{eqnarray} \tag{ A3 }$$

and

$$\begin{eqnarray}I(k,{k}_{i}) & = & \displaystyle \frac{3{{\rm{\Lambda }}}^{3}}{8{k}_{i}^{3}}\left[\displaystyle \frac{{\rm{\Lambda }}}{k}\left({{\rm{e}}}^{-{\left(\displaystyle \frac{k+{k}_{i}}{{\rm{\Lambda }}}\right)}^{2}}-{{\rm{e}}}^{-{\left(\displaystyle \frac{k-{k}_{i}}{{\rm{\Lambda }}}\right)}^{2}}\right)\right.\\ & & +\left.\sqrt{\pi }\left({\rm{erf}}\left(\displaystyle \frac{k+{k}_{i}}{{\rm{\Lambda }}}\right)-{\rm{erf}}\left(\displaystyle \frac{k-{k}_{i}}{{\rm{\Lambda }}}\right)\right)\right],\end{eqnarray} \tag{ A4 }$$

with ${\rm{\Lambda }}=2/\alpha$.

The energy density in polarized pure neutron matter is given in [2]. In principle, the parameters ${\varepsilon }_{0}^{l},{\varepsilon }_{\gamma }^{l}$ and ${\varepsilon }_{{ex}}^{l}$ split into the interactions between pairs of neutrons with the same (index $l,l$) and the opposite (index $l,{ul}$) spin orientation:

$$\begin{eqnarray}{H}_{{\rm{pol}}}^{N}({\rho }_{{nu}},{\rho }_{{nd}}) & = & \left[\displaystyle \frac{3{\ℏ}^{2}({k}_{{nu}}^{2}{\rho }_{{nu}}+{k}_{{nd}}^{2}{\rho }_{{nd}})}{10M}+\displaystyle \frac{{\mathrm{\ϵ}}_{0}^{l,l}}{2{\rho }_{0}}({\rho }_{{nu}}^{2}+{\rho }_{{nd}}^{2})\right.\\ & & +\displaystyle \frac{{\mathrm{\ϵ}}_{0}^{l,{ul}}}{{\rho }_{0}}{\rho }_{{nu}}{\rho }_{{nd}}+\left(\displaystyle \frac{{\mathrm{\ϵ}}_{\gamma }^{l,l}}{2{\rho }_{0}^{\gamma +1}}({\rho }_{{nu}}^{2}+{\rho }_{{nd}}^{2})+\displaystyle \frac{{\mathrm{\ϵ}}_{\gamma }^{l,{ul}}}{{\rho }_{0}^{\gamma +1}}{\rho }_{{nu}}{\rho }_{{nd}}\right){\left(\displaystyle \frac{\rho ({\bf{R}})}{1+b\rho ({\bf{R}})}\right)}^{\gamma }\\ & & +\displaystyle \frac{{\mathrm{\ϵ}}_{{ex}}^{l,l}}{2{\rho }_{0}}({\rho }_{{nu}}^{2}J({k}_{{nu}})+{\rho }_{{nd}}^{2}J({k}_{{nd}}))\\ & & +\left.\displaystyle \frac{{\mathrm{\ϵ}}_{{ex}}^{l,{ul}}}{4{\rho }_{0}{\pi }^{2}}\left({\rho }_{{nu}}{\displaystyle \int }_{0}^{{k}_{{nd}}}I(k,{k}_{{nu}}){k}^{2}{\rm{d}}{\⁢}{k}+{\rho }_{{nd}}{\displaystyle \int }_{0}^{{k}_{{nu}}}I(k,{k}_{{nd}}){k}^{2}{\rm{d}}{\⁢}{k}\right)\right],\end{eqnarray} \tag{ A5 }$$

where, ${k}_{{ni}}={(6{\pi }^{2}{\rho }_{{ni}})}^{1/3}(i=u,d)$ and the functions J (k_ni) and I (k,k_ni) are calculated from expression (A7) and (A8) by replacing k_i by k_ni. However, as explained in [2] and mentioned before, the limit of fully polarized neutron matter requires ${\varepsilon }_{\gamma }^{l,l}=0$ and ${\varepsilon }_{0}^{l,l}=-{\varepsilon }_{{ex}}^{l,l}$ and the fit to the DBHF effective mass in polarized neutron matter determines ${\varepsilon }_{{ex}}^{l,l}={\varepsilon }_{{ex}}^{l}/3$

The one-body density matrix as well as the particle, kinetic energy and spin local densities are obtained from the single-particle orbitals ${\phi }_{i}$ that define the Slater determinant ${{\rm{\Psi }}}_{0}$ as

$$\begin{eqnarray}&&{\rho }_{q}({\bf{r}},{{\bf{r}}}^{\prime })=\displaystyle \sum _{i=1}^{{A}_{q}}\displaystyle \sum _{\sigma }{\phi }_{i}^{*}({\bf{r}},\sigma ,q){\phi }_{i}({{\bf{r}}}^{\prime },\sigma ,q),\end{eqnarray} \tag{ B1 }$$

$$\begin{eqnarray}&&{\rho }_{q}({\bf{r}})=\displaystyle \sum _{i=1}^{{A}_{q}}\displaystyle \sum _{\sigma }| {\phi }_{i}({\bf{r}},\sigma ,q){| }^{2},\end{eqnarray} \tag{ B2 }$$

$$\begin{eqnarray}&&{\tau }_{q}({\bf{r}})=\displaystyle \sum _{i=1}^{{A}_{q}}\displaystyle \sum _{\sigma }| \nabla {\phi }_{i}({\bf{r}},\sigma ,q){| }^{2},\end{eqnarray} \tag{ B3 }$$

and

$$\begin{eqnarray}&&{J}_{q}({\bf{r}})=i\displaystyle \sum _{i=1}^{{A}_{q}}\;\displaystyle \sum _{\sigma ,{\sigma }^{\prime }}{\phi }_{i}^{*}({\bf{r}},\sigma ,q)[{({\boldsymbol{\sigma }})}_{\sigma ,{\sigma }^{\prime }}\times \nabla ]{\phi }_{i}({\bf{r}},\sigma ,q),\end{eqnarray} \tag{ B4 }$$

respectively. The total energy density ${ \mathcal H }$ can be written as

$$\begin{eqnarray}{ \mathcal H }=\displaystyle \frac{{\ℏ}^{2}}{2m}({\tau }_{n}+{\tau }_{p})+{{ \mathcal H }}_{d}^{{\rm{Nucl}}}+{{ \mathcal H }}_{{\rm{exch}}}^{{\rm{Nucl}}}+{{ \mathcal H }}_{{zr}}^{{\rm{Nucl}}}+{{ \mathcal H }}^{{\rm{SO}}}+{{ \mathcal H }}^{{\rm{Coul}}},\end{eqnarray} \tag{ B5 }$$

where, we have split the total energy density as the sum of the kinetic part, direct, exchange and zero-range nuclear terms together with the spin–orbit and Coulomb contributions.

The Coulomb energy is taken in the usual way as a direct plus Slater exchange contribution computed with the point proton density:

$$\begin{eqnarray}&&{{ \mathcal H }}^{{\rm{Coul}}}({\bf{R}})=\displaystyle \frac{{e}^{2}}{2}\int \displaystyle \frac{{\rho }_{p}\left({\bf{R}}+\frac{{\bf{r}}}{2}\right){\rho }_{p}\left({\bf{R}}-\frac{{\bf{r}}}{2}\right)}{r}{{\rm{d}}}^{3}r-\displaystyle \frac{3{e}^{2}}{4}{\left(\displaystyle \frac{3}{\pi }\right)}^{1/3}{\rho }_{p}^{4/3}({\bf{R}}),\end{eqnarray} \tag{ B6 }$$

where, ${\bf{r}}={{\bf{r}}}_{{\bf{1}}}-{{\bf{r}}}_{{\bf{2}}}$ and ${\bf{R}}=\frac{{{\bf{r}}}_{{\bf{1}}}+{{\bf{r}}}_{{\bf{2}}}}{{\bf{2}}}$ are the relative and center of mass coordinates, respectively, for the two interacting nucleons located at ${{\bf{r}}}_{{\bf{1}}}$ and ${{\bf{r}}}_{{\bf{2}}}$.

In the SEI model the spin–orbit interaction is chosen in the form used in the case of Skyrme and Gogny forces: ${v}_{i,j}^{{\rm{SO}}}={{iW}}_{0}({{\boldsymbol{\sigma }}}_{i}+{{\boldsymbol{\sigma }}}_{j})[{{\bf{k}}}^{\prime }\times {\boldsymbol{\delta }}({r}_{i},{r}_{j})k]$. The corresponding contribution to the energy density is

$$\begin{eqnarray}&&{{ \mathcal H }}^{{\rm{SO}}}({\bf{R}})=-\displaystyle \frac{1}{2}{W}_{0}[\rho ({\bf{R}}){\boldsymbol{\nabla }}{\bf{J}}+{\rho }_{n}({\bf{R}}){\boldsymbol{\nabla }}{{\bf{J}}}_{{\bf{n}}}+{\rho }_{p}({\bf{R}}){\boldsymbol{\nabla }}{{\bf{J}}}_{{\bf{p}}}].\end{eqnarray} \tag{ B7 }$$

The direct and exchange contributions to the nuclear energy density read

$$\begin{eqnarray}{{ \mathcal H }}_{d}^{{\rm{Nucl}}}({\bf{R}}) & = & \displaystyle \frac{1}{2}\displaystyle \int {{\rm{d}}}^{3}{{re}}^{-{r}^{2}/{\alpha }^{2}}\left[\left(W+\displaystyle \frac{B}{2}-H-\displaystyle \frac{M}{2}\right)\rho \left({\bf{R}}+\displaystyle \frac{{\bf{r}}}{2}\right)\rho \left({\bf{R}}-\displaystyle \frac{{\bf{r}}}{2}\right)\right.\\ & & -\left(W+\displaystyle \frac{B}{2}\right)\left({\rho }_{n}\left({\bf{R}}+\displaystyle \frac{{\bf{r}}}{2}\right){\rho }_{p}\left({\bf{R}}-\displaystyle \frac{{\bf{r}}}{2}\right)\right.\\ & & +\left.\left.{\rho }_{p}\left({\bf{R}}+\displaystyle \frac{{\bf{r}}}{2}\right){\rho }_{n}\left({\bf{R}}-\displaystyle \frac{{\bf{r}}}{2}\right)\right)\right],\end{eqnarray} \tag{ B8 }$$

and

$$\begin{eqnarray}{{ \mathcal H }}_{{\rm{exch}}}^{{\rm{Nucl}}}({\bf{R}}) & = & \displaystyle \int {{\rm{d}}}^{3}{{re}}^{-{r}^{2}/{\alpha }^{2}}\left[\displaystyle \frac{1}{2}\left(M+\displaystyle \frac{H}{2}-B-\displaystyle \frac{W}{2}\right)\left[{\rho }_{n}^{2}\left({\bf{R}}+\displaystyle \frac{{\bf{r}}}{2},{\bf{R}}-\displaystyle \frac{{\bf{r}}}{2}\right)\right.\right.\\ & & +\left.{\rho }_{p}^{2}\left({\bf{R}}+\displaystyle \frac{{\bf{r}}}{2},{\bf{R}}-\displaystyle \frac{{\bf{r}}}{2}\right)\right]\\ & & +\left.\left(M+\displaystyle \frac{H}{2}\right){\rho }_{n}\left({\bf{R}}+\displaystyle \frac{{\bf{r}}}{2},{\bf{R}}-\displaystyle \frac{{\bf{r}}}{2}\right){\rho }_{p}\left({\bf{R}}+\displaystyle \frac{{\bf{r}}}{2},{\bf{R}}-\displaystyle \frac{{\bf{r}}}{2}\right)\right].\end{eqnarray} \tag{ B9 }$$

Finally, the contribution from the zero-range part of the interaction can be written as

$$\begin{eqnarray}{{ \mathcal H }}_{{zr}}^{{\rm{Nucl}}} & = & \displaystyle \frac{{t}_{0}}{4}[(1-{x}_{0})[{\rho }_{n}^{2}({\bf{R}})+{\rho }_{p}^{2}({\bf{R}})]+(4+2{x}_{0}){\rho }_{n}({\bf{R}}){\rho }_{p}({\bf{R}})]\\ & & +\displaystyle \frac{{t}_{3}}{24}[(1-{x}_{3})[{\rho }_{n}^{2}({\bf{R}})+{\rho }_{p}^{2}({\bf{R}})]]{\left(\displaystyle \frac{\rho ({\bf{R}})}{1+b\rho ({\bf{R}})}\right)}^{\gamma }\\ & & +\displaystyle \frac{{t}_{3}}{24}[(4+2{x}_{3}){\rho }_{n}({\bf{R}}){\rho }_{p}({\bf{R}})]{\left(\displaystyle \frac{\rho ({\bf{R}})}{1+b\rho ({\bf{R}})}\right)}^{\gamma }.\end{eqnarray} \tag{ B10 }$$