Nuclear symmetry energy and the role of threebody forces
Abstract
Density dependence of nuclear symmetry energy as well as its partial wave decomposition is studied within the framework of lowestorder constrained variational (LOCV) method using AV18 twobody interaction supplemented by UIX threebody force. The main focus of the present work is to introduce a revised version of threebody force which is based on an isospindependent parametrization of coefficients in the UIX force, in order to overcome the inability to produce correct saturationpoint parameters in the framework of LOCV method. We find that employing the new model of threebody force in the LOCV formalism leads to successfully reproducing the semiempirical parameters of cold nuclear matter, including , , and . All our models of threebody force combined with AV18 twobody force give maximum neutron star mass higher than . The fraction of protons in the nucleon cores of neutron stars strongly depends on the threebody force parametrization.
 PACS numbers

21.30.Fe,21.65.Mn,21.65.Cd,21.65.Ef
pacs:
21.30.Fe,21.65.Mn,21.65.Cd,21.65.EfI Introduction
One of the aims of nuclear science is to understand the properties of strongly interacting bulk matter at nuclear levels. In this regard, properties of symmetric nuclear matter (SNM) has been studied for a long time and its equation of state (EOS) around saturation density is known rather well. On the other hand, the EOS of isospin asymmetric nuclear matter, particularly pure neutron matter (PNM), has not been established well despite of its crucial role in both nuclear physics and astrophysics. Motivated by this fact, studying the nuclear symmetry energy which mostly governs the PNM EOS is currently an active field of research in nuclear physics.
Since the symmetry energy can not be directly measured, it is of fundamental importance to identify observables which strongly correlate with and its density slope to impose constraints on these quantities. Additional information on the value and density dependence of symmetry energy can be extracted from the astrophysical observations of compact objects, namely neutron stars. Among the terrestrial laboratory observables, the neutron skin thickness is widely studied 1 (); 2 (); 3 (); 4 (); 5 (); 6 () during last decades. Nowadays, there are many other observables such as isospin fractionation, neutronproton differential flow, nuclear mass systematics, lowlying E1 mode, etc., which are believed to be sensitive to symmetry energy and particularly, its density slope. We refer to references 7 (); 8 (); 9 (); 10 (); 11 () and references therein quoted for more details. Density dependence of this quantity also relates the heavyion reactions 12 (); 13 (); 14 (); 15 (); 16 (); 17 (); 18 (), stability of superheavy nuclei 19 (), cooling 20 (); 21 (); 22 (); 23 (); 24 () and the massradius relations of neutron stars 25 (); 26 (), and properties of nuclei involved in rprocess nucleosynthesis 27 ().
In order to theoretically predict the symmetry energy, asymmetric nuclear matter is studied in the framework of various microscopic and phenomenological approaches by using a variety of microscopic and phenomenological twonucleon forces and phenomenological threebody forces. Such manybody techniques include the BruecknerHartreeFock (BHF) 7 (); 28 (), its relativistic counterpart DiracBruecknerHartreeFock (DBHF) 29 (), the selfconsistent Green’s function 30 (), variational approaches 31 (); 32 (), meanfield model 33 (), the Quantum Monte Carlo (QMC) 34 (), extended relativistic mean field (ERMF) model 35 (), and SkyrmeHartreeFock (SHF) method 36 (). There are also some recent calculations for the symmetry energy of finite nuclei (see e.g. 37 () and references therein). Recent progress and challenges in both theoretically and experimentally measuring the symmetry energy and its density dependence, specially at supra nuclear densities, are also reviewed in Refs. 38 (); 39 (); 40 (); 41 ().
Despite of numerous attempts that was made to determine the symmetry energy properties, especially the value of and its density slope at saturation density , these quantities are still uncertain. Results are strongly dependent on different experimental techniques and manybody approaches as well as nuclear interactions. Especially, the density dependence of symmetry energy is poorly known. The high density prediction of manybody models for the symmetry energy are extremely diverse and sometimes contradictory. On the other hand, although the value of the symmetry energy is known to be in the range of MeV 40 (), explorations on show wide variations. The obtained lower limit is MeV 42 () and the upper one is even higher than 170 MeV 43 (). However, recent surveys of analyses of terrestrial nuclear laboratory experiments and astrophysical observations found that the mean value of the slope parameter at saturation density lies in the range of MeV 40 (); 11 ().
In the present work, within the selfconsistent lowest order constrained variational (LOCV) approach using the Argonne V18 twobody potential 44 () supplemented by an Urbana type threebody force (UIX) 45 (), we study the density dependence of nuclear symmetry energy with the focus on the effect of threebody interaction on this quantity. The LOCV method is a wellknown manybody technique that was originally used to study the properties of cold symmetric nuclear matter 46 (); 47 () by using the Reidtype potential 48 (); 49 () as the bare twobody interaction. Later on, this approach was extended to finite temperature 50 () and also calculation of the EOS of asymmetric nuclear matter 51 (), pure neutron matter 52 (), and beta stable matter 53 () was carried out within this framework by using more sophisticated potentials. Moreover, relativistic corrections has been considered in calculating thermodynamic properties of nuclear matter within this model at both zero 54 () and finite temperature 55 (). Recently, this technique has been extended by adding threebody forces to its formalism 56 () and has been used to study the structure of neutron stars 56 () as well as protoneutron stars 57 ().
Although this method was successful in reproducing some properties that characterize the cold symmetric nuclear matter, namely the saturation density, saturation energy and the incompressibility, the nuclear symmetry energy obtained by this model turns out to be high and does not lie in the proper range 56 (). Therefore, our goal in the present work is to correct this deficiency by revising the parametrization of UIX threebody interaction.
The article is organized as follows. In Sec. II we review briefly the adopted LOCV manybody approach including threebody interaction. Section III is devoted to theoretically studying the nuclear symmetry energy. Results and discussion are presented in Sec. IV. Finally summary and conclusion are given in Sec. V.
Ii Theoretical Formalism
In this section we briefly describe the LOCV formalism whose details can be found in mentioned references in the introduction section. As the first step of calculating the EOS of asymmetric nuclear matter within the LOCV approach, we produce a trial wave function for Abody interacting system at zero temperature, denoted by , which is defined as the product of a noninteracting ground state wave function of independent nucleons, , and an body correlation operator , as follows 53 ()
(1) 
The correlation operator is in general considered as the symmetrized product of twobody correlation function operators, i.e.,
(2) 
where is the symmetrizing operator. is expressed as
(3) 
Here, stands for and . is set to unity except for triplet channels with . In this case we choose =2 and 3 with and . The operator denotes the usual tensor operator. In general, the nuclear Hamiltonian includes a nonrelativistic onebody kinetic energy as well as a twobody potential ,
(4) 
By neglecting the small contribution from higherorder terms in the cluster expansion series 50 (), the energy expectation value is gained via
(5) 
with
(6) 
where , is corresponding nucleon Fermi momentum divided by , and
(7) 
where is expressed as
(8) 
After doing some algebra, final expression for the twobody cluster energy is obtained by minimizing expression given by Eq.(8) with respect to variations in the functions under the condition
(9) 
Equation (9) defines the LOCV normalisation constraint and can be calculated by using the modified Pauli function which in case of asymmetric nuclear matter takes the following form
(10)  
denotes the spherical Bessel function of order . The normalisation constraint introduces another parameter in the LOCV formalism, i.e. the Lagrange multiplier . Procedure of minimizing Eq. (8) leads to a number of EulerLagrange differential equations for functions . By solving these equations, correlation functions and consequently, the twobody cluster energy can be determined.
The important role played by threebody forces (3BF) in both finite nuclei and nuclear matter calculations is well known. In nuclear matter calculations, unfortunately, whatever realistic twobody force (2BF) is used in a manybody approach, the saturation properties of cold symmetric nuclear matter fail to be reproduced correctly. This deficiency can be cured by inclusion of a 3BF in the nuclear Hamiltonian. For this reason, the semiphenomenological UIX interaction is considered in our manybody calculations. Generally, the UIX interaction is written as following 45 ()
(11) 
with a twopion exchange contribution as
(12) 
and a shorterrange phenomenological part,
(13) 
and are adjustable parameters determined by fitting the empirical saturation density and energy of the cold SNM in the LOCV calculations. Numerals 1, 2, and 3 stand for three different interacting nucleons while , , and denotes the spin, isospin and the usual tensor operator respectively and and are the Yukawa and tensor functions. The onepion exchange operator in Eq. (12) is defined as
(14) 
In two recent papers we have presented the procedure of extending the LOCV approach by using this kind of 3BF in its formalism at both zero 56 () and finite temperature 57 (). The 3BF is included via an effective twobody potential derived after averaging out the third particle, being weighted by the LOCV twobody correlation functions at fixed density , i.e.,
(15) 
where . In this way, a densitydependent effective twobody interaction is gained with the following operator structure
(16)  
Explicit relations for the coefficients , , and are expressed in 56 (). It is pointed out in the mentioned reference that some quantities that characterizes the saturation point of cold symmetric nuclear matter can not be produced within the LOCV approach even if the nuclear Hamiltonian is supplemented by a 3BF. This deficiency motivated us to propose a new parametrization for UIX 3BF which is based on isospindependent strengths. Our main purpose in the present work is to construct such a 3BF in order to reach correct saturation properties within the LOCV model. In this regard, subsection B of Sect.IV is devoted to investigating this issue in detail. The next Section III refers to the nuclear symmetry energy and other parameters characterizing saturation point of nuclear matter.
Model  

LOCV (AV18) 56 ()  0.3270  23.37  373.31  40.52  74.07  76.30  520.57 
BHF (AV18) 28 ()  0.240  17.30  213.6  35.8  63.1  27.8  339.6 
QMC (AV8) 34 ()  0.16  16.0    30.50  31.30     
DBHF (Bonn A) 29 ()  0.181  16.62  233  34.80  71.20     
DBHF (Bonn B) 29 ()  0.162  15.04  190  31.20  55.90     
DBHF (Bonn C) 29 ()  0.148  14.14  170  28.9  46.7     
DBHF (Bonn AB) 29 ()  0.17  15.52  204  32.3  61.1     
LOCV (AV18+UIX) 56 ()  0.1748  15.58  295.77  39.9  115.75  15.58  678.92 
BHF(AV18+UIX) 7 ()  0.187  15.23    34.30  66.50     
QMC (AV8+UIX) 34 ()  0.16  16.0    35.10  63.60     
variational (AV8+UIX) 32 ()  0.16  16.09  245  30.0       
BHF (AV18+UIX) 28 ()  0.176  14.62  185.9  33.6  66.9  23.4  343.8 
mean field 33 ()  0.16  16.0    125.79    
ERMF (BSR1) 35 ()  0.1481  16.02  240.05  30.98  59.61     
ERMF (BSR7) 35 ()  0.1493  16.17  231.86  36.99  98.78     
SHF 36 ()  0.16  16.0  230.0      
FSU 61 ()  0.148  16.30  230.0  32.6  60.5  51.3  276.6 
NL3 62 ()  0.148  16.24  271.6  37.4  118.5  100.9  698.4 
TM1 63 ()  0.145  16.32  281.0  36.8  110.8  66.4  518.7 
DDME1 64 ()  0.152  16.23  332.8  33.1  55.6  100.8  508.1 
RMF (NL3) 37 ()  0.148  16.3  272  37.4  118.2     
RMF (TM1) 37 ()  0.145  16.3  281  36.9  110.8     
RMF (FSU) 37 ()  0.148  16.3  230  32.6  60.5     
RMF (IUFSU) 37 ()  0.155  16.4  231  31.3  47.2     
Iii Nuclear symmetry energy
The energy per particle of cold asymmetric nuclear matter at given baryon density and asymmetry parameter () can be expanded around as
(17) 
The expansion coefficient is the nuclear symmetry energy and is defined as following
(18) 
Due to the charge symmetry of nuclear forces, odd powers of are absent in Eq. (17). Magnitude of the term at saturation density has been estimated to be less than 1 MeV 58 (). Therefore, higherorder terms in can be neglected compared to the value of the quadratic term. Equation (17) is known as the parabolic approximation for the EOS of asymmetric nuclear matter which is accurate close to and also works well for higher values of the asymmetry parameter. Within this approximation, the symmetry energy can be defined as
(19) 
which states that to a good approximation, the symmetry energy measures the difference between the energy per particle in uniform isospinsymmetric matter and pure neutron matter at the fixed density . Then, one can expand the nuclear symmetry energy around the nuclear matter saturation density . To second order the expansion is written as
(20) 
where the coefficients respectively denote the nuclear symmetry energy value at normal density, the slope, and the curvature parameter. The and characterize the density dependence of symmetry energy around the saturation density and are written as
(21) 
(22) 
On the other hand, it is possible to expand the isobaric incompressibility of asymmetric nuclear matter in the power series of asymmetry parameter as
(23) 
where is the incompressibility of SNM at normal density and is expressed as
(24) 
The coefficient characterizes the isospin dependence of the incompressibility at saturation density and is often determined by using the following approximate expression 59 ()
(25) 
Once again, one can keep only the first two terms in Eq. (23) and neglect higher order ones in 60 (). Quantities and together with the slope parameter can be determined experimentally and provide information about the density dependence of nuclear symmetry energy at normal nuclear density.
Iv Results and Discussion
We first focus on calculations which are done by using only 2BF as well as 2BF+3BF with isospinindependent parameters in UIX 3BF. After that, in the second subsection we present our results obtained using revised version of 3BF based on the isospindependent parametrization of coefficients and in the UIX threebody force.
iv.1 UIX threebody force with isospinindependent A and U parameters
In the beginning, we mention that in order to calculate the effective twobody potential via Eq. (15), the channel twobody correlation function is used, as discussed in our previous work 56 (). Moreover, for each specific asymmetry, the twobody correlation function which is obtained at the same asymmetry is used. Let us now start with presenting results of the energy per nucleon. Figure 1 shows the comparison between the energy per particle of SNM and PNM using 2BF as well as 2BF+3BF. Properties of saturation point of SNM in both cases are also listed in Table I. As has been mentioned earlier and is seen in Fig.1, if twobody potential is the only interaction considered in the nuclear Hamiltonian, calculations fail to reproduce the saturation point of symmetric nuclear matter, i.e. they yield too small binding energy and a saturation density well above the empirical value of . Moreover, from Table I it can be concluded that the isobaric incompressibility of SNM, , which is defined by Eq. (24), is too large and does not lie in the rage of MeV which is calculated by reanalysing of recent data on the giant monopole resonance 65 ()
Now let us focus on the case of using both 2BF and 3BF in the calculations. As expected, inclusion of 3BF leads to an acceptable saturation point by increasing the binding energy at saturation as well as pushing the saturation density toward lower values. The isobaric incompressibility at normal density is also decreased so that it fits well in the proper range. It is found that values MeV and MeV yield a reasonable saturation point for cold SNM in the LOCV calculations.
By using the empirical parabolic law, nuclear symmetry energy can be easily extracted from LOCV microscopic calculations via Eq. (19) provided the parabolic approximation is valid. Quadratic dependence of the symmetry energy as a function of isospin asymmetry at different densities is shown in Fig. 2. Panel (a) of this figure is devoted to results obtained by using only twobody interaction while the other panel shows the results when the 3BF is included. In both cases, it is seen that the parabolic law is truly valid in the whole range of at least up to moderate densities. Therefore, to a good approximation, can be determined by using the parabolic law.
Density dependence of the nuclear symmetry energy from both 2BF and 2BF+3BF potentials is shown in Fig. 3. As can be seen, inclusion of 3BF causes a significant difference in both the behaviour and the value of compared to the case of using only 2BF. The latter leads to a soft symmetry energy that tends to be saturated at high densities, whereas the densitydependent nature of the 3BF provides more repulsive contribution to the PNM EOS and consequently, causes the symmetry energy to increase sharply with density. Exact values of symmetry energy calculated from Eq. (18) are also plotted in this figure. It is seen that there is no significant difference between the exact curves and those obtained from the parabolic approximation, i.e. Eq. (19), as is concluded from Fig. 2.
On the other hand, from Table I one can conclude that the value of nuclear symmetry energy at normal density, , turns out to be too large in the case of using only 2BF, although calculations predict the value of 74.07 MeV for the slope parameter at saturation density which is consistent with the constraint of MeV 40 (). In order to see whether the model predicts the correct curvature parameter, , one should use the approximate relation of . is the isospindependent part of the incompressibility of SNM at normal density (see Eq. (23)) and can be extracted experimentally. Estimate of MeV has been obtained from isospin diffusion data 66 (); 67 (), while a systematic study of giant monopole resonance for Sn isotopes predicts MeV 68 (). The value of this quantity obtained within the LOCV model using 2BF is 520.57 MeV which lies well in the experimentally determined range. This result shows that the obtained is also reasonable. For the sake of comparison, results from a number of manybody techniques and interactions are also reported in Table I.
Although considering 3BF in the LOCV formalism was successful in producing the empirical saturation point, unfortunately it fails to provide a correct value for symmetry energy at normal density. This model predicts the value of MeV which although is smaller than the result obtained by using 2BF, is still clearly too large compared to the experimentally determined range of MeV 40 (). Moreover, as can be seen in Fig. 1, the large contribution that comes from the PNM energy, provides a stiff nuclear symmetry energy with large slope parameter MeV. Obtained value for is 678.92 MeV which is smaller than the lower limit of 650 MeV predicted from giant monopole resonance for Sn isotopes.
(0,0)  39.9  115.75  15.58  678.92 

(0.0015,0.000523) D  35.0  89.91  47.91  587.37 
(0.0175,0.000261) C  33.0  74.04  101.42  545.66 
(0.0337,0) B  31.0  57.43  159.13  503.71 
iv.2 Revised version of the UIX threebody force with isospindependent strengths
The earliest model of Pwave twopion exchange potential is due to Fujita and Miyazawa 69 (), who assumed that it is entirely due to the excitation of the resonance. Neglecting the nucleon and kinetic energies we obtain
(26) 
where is the mass of particle and and are the and coupling constants, respectively. Using the observed values of , , and , Eq.(26) predicts that MeV 70 (). In some models of , by starting from the observed pionnucleon scattering amplitude, and using current algebra and PCAC constraints, or chiral symmetry, higher absolute value of is reached which is a consequence of taking into account the contributions of all the resonances, as well as that of Swave scattering to the twopion exchange 3BF. For example, TucsonMelbourne 71 () and Texas 72 () models predict values of 0.063 MeV and 0.09 MeV for the strength , respectively 70 (). On the other hand, in all the Urbana models of 3BF, is treated as an adjustable parameter and is varied to fit the data. In light nuclei calculations, parameters and are chosen to yield the observed binding energies of and . On the other, because of the likely manybody effects which are not yet known well, there is no reason to believe that 3BF in nuclear matter system be the same as in light nuclei. As a consequence, in the case of nuclear matter, the values of 3BF parameters are adjusted to yield the empirical saturation point of SNM as well as the isobaric incompressibility 73 (); 74 ().
Because of our poor knowledge of Hamiltonian in high density nuclear systems, there is no reason to believe that the 3BF among nnn be the same as among nnp, npp and ppp. Therefore, it seems reasonable to correct the Hamiltonian by considering an isospin asymmetrydependent parametrization for the strength of 3BF.
As already concluded from previous subsection, within the LOCV method, it seems that the current form of 3BF fails to reproduce the semiempirical values of all the quantities that characterize the saturation point, at the same time. This fact motivated us to revise the procedure of calculating 3BF. Parameters and in the UIX threebody force (see Eqs. (12) and (13)) are usually considered as constants which are determined by fitting the empirical saturation density and energy of the cold SNM. In the present work, we extend this procedure, by proposing an isospin dependent structure for these parameters. Namely, we assume that and can be expanded to secondorder in asymmetry parameter as
(27) 
and
(28) 
It must be noted that this consideration will not violate the formalism of LOCV method, since one still has two, but different parameters in 3BF for each given nuclear system i.e., SNM, PNM and stable matter. Rewriting the 3BF strengths in the form of Eqs. (27) and (28) only changes the nuclear Hamiltonian and it surely keeps the manybody formalism unchanged. Such correction in Hamiltonian is also seen in other papers such as that by Akmal, et al 75 (), where in addition to 2BF and 3BF, the authors consider a densitydependent term with adjustable strength in order to reproduce correctly the saturation properties of SNM.
By putting in Eqs. (27) and (28), i.e. considering symmetric nuclear matter, one reaches the previous scenario. Therefore, and must have the same values as those mentioned before, i.e. 0.041 MeV and 0.000523 MeV, respectively. Variables and are determined in such a way that both the symmetry energy and its slope lie in the empirical range. For the whole range of and reported in Fig.4, values of both and lie well in the range of MeV 72 (); 76 () and MeV 77 (); 70 () which are obtained by other authors, in case of any nuclear matter system (SNM, PNM and stable matter).
(0,0)  13.75  27.49  27.50  

26.15  88.26  43.08  
Total  39.90  115.75  15.58  
(0.0015,0.000523) D  13.75  27.49  27.50  
21.25  62.42  20.41  
Total  35.00  89.91  47.91  
(0.0175,0.000261) C  13.75  27.49  27.50  
19.25  46.55  73.92  
Total  33.00  74.04  101.42  
(0.0337,0) B  13.75  27.49  27.50  
17.25  29.94  131.63  
Total  31.00  57.43  159.13 
In order to achieve the proper value for and , a softer EOS for PNM should be produced. Therefore, the repulsive part of 3BF, given by Eq. (13), which is controlled by the strength , should be smaller. One way to do this task is to assume that there is no repulsive term in 3BF for PNM . This aim can be reached by putting in Eq. (28). Another way is to assume that the repulsive part of 3BF for PNM is the same as that for SNM which means setting in Eq. (28). Choosing positive values for will increase the repulsive component and does not help solving the problem. On the other hand, if one sets , expression given by Eq. (28) becomes negative for PNM which is against the nature of this phenomenological repulsive component. As a consequence, can vary between two extreme values mentioned above, i.e., .
Figure 4 shows the resulting and in the plane of . Each shaded parallelogram shows a specific span for symmetry energy at normal density while each diagonal black line corresponds to a fixed value of slope parameter. As an example, all pairs which are located within the tilted box produce MeV. Similarly, all pairs which are between the lines and , produce MeV. As can be seen in this figure, by setting in Eq. (27) one reaches to MeV and MeV for the whole range of which are large compared to empirical values. Therefore, should be also nonzero in order to get smaller symmetry energy as well as the slope parameter. From this figure, it can also be concluded that both and vary with parameters and , although the slope parameter seems to be more sensitive. Particularly, for a fixed value of , both the symmetry energy and the slope parameter decrease with while at fixed both and increase as tends to zero.
In order to understand the difference between the revised version of 3BF (which is characterized by Eqs. (27) and (28)) and the old one (the case of ), a comparison is made for each component of the effective twobody interaction, namely the coefficients of Eq. (16). To do this task, three pairs which are denoted by letters B, C, and D in Fig. 4, are compared with the 3BF with isospinindependent strengths, i.e. the case of . Components , , and of these four points are respectively plotted in Figs. 5 to 7 at for PNM. Corresponding numerical values of and are also reported in Table II. Twopion exchange contribution to the UIX 3BF is divided into two components, namely and which are both controlled by the strength . Now let us consider the PNM case which is reached by putting in Eqs. (26) and (27), and keeping in mind that is negative, as mentioned earlier. Consequently, the strength tends to zero by increasing and one expects to get twopion exchange components with smaller absolute value. This fact can be clearly seen in Figs. 5 and 6. Among four chosen cases, the largest value of corresponds to point B (see Fig. 4). Therefore, the smallest absolute value of effective twobody potential is obtained in this case.
Partial wave  

(0,0)  0.47  11.58  4.76  
(0.0015,0.000523) D  1.27  5.29  0.45  
(0.0175,0.000261) C  1.1  4.67  5.35  
(0.0337,0) B  0.97  3.76  12.04  
(0,0)  25.87  61.76  29.41  
D  
C  
B  
(0,0)  4.67  19.19  26.36  
D  
C  
B  
(0,0)  0.76  0.08  8.4  
D  1.03  1.38  4.92  
C  0.73  0.96  14.1  
B  0.44  3.19  21.17  
(0,0)  4.38  19.21  40.08  
D  3.77  14.74  23.97  
C  3.64  13.13  15.08  
B  3.59  11.64  5.99  
(0,0)  0.03  12.25  54.53  
D  1.86  1.18  23.80  
C  3.63  13.09  30.63  
B  5.43  27.83  86.59  
(0,0)  2.50  10.00  4.87  
D  2.76  11.78  11.72  
C  2.43  9.19  1.71  
B  2.10  6.56  8.60  
(0,0)  0.14  0.64  7.33  
D  
C  
B  
(0,0)  3.37  13.83  14.86  
D  
C  
B  
(0,0)  0.09  0.61  11.58  
D  0.02  1.39  13.59  
C  0.91  5.69  17.17  
B  1.82  10.02  20.52 
Figure 7 shows the repulsive component of UIX 3BF whose strength is governed by parameter . Once again, the isospin asymmetry is set to 1 (the case of PNM). As already pointed out, is a positive parameter whereas is a negative one which varies in the range of . Therefore, larger negative values of lead to smaller and consequently weaker repulsive components, as is seen in Fig.7. On the other hand, point D corresponds to the lower limit predicted for , i.e. which is equivalent to the case of . Therefore, at point D there is no repulsive component.
In Fig. 8(a) the energy per nucleon of PNM is presented which is obtained by using four different effective twobody potentials plotted in Figs. 5 to 7. It is seen that the EOS is strongly sensitive to the proposed parametrization of variables and . All curves with nonzero and are located below the case of at the whole range of density. Therefore, it can be concluded that by producing an isospindependent structure for the 3BF parameters, our wish for gaining a softer EOS for PNM can be fulfilled.
Figure 8(b) represents the corresponding results for the density dependence of nuclear symmetry energy. A significant difference can be seen in both the value and behaviour of depending on the values of and . While symmetry energy obtained from parametrization increases with density, those of isospindependent and tend to saturate at large densities. This tendency decreases from point B to D. Values of symmetry energy at saturation density , slope parameter , curvature parameter , and the isospindependent part of the isobaric incompressibility are reported in Table II for all curves. From this Table it is clear that all mentioned quantities which characterize the empirical saturation point, lie well within the experimentally predicted ranges.
In Fig. 8(c) the density dependence of predicted by LOCV method assuming model D of revised UIX 3BF (curve D in Fig. 8(b)) is compared with the results of other manybody approaches such as variational 31 () and BHF 69 (); 70 (); 71 () models using different potentials. It is seen that both the values and behaviour of obtained using D version of revised UIX 3BF in the LOCV method is in agreement with results of BHF calculations with 3BF.
One of the quantities which is directly affected by the symmetry energy is the relative population of protons in the stable matter which is a matter composed of an uncharged mixture of neutrons, protons, electrons, and muons in equilibrium with respect to the weak interaction. Figure 9 displays this quantity as a function of total baryon density for different pairs, namely , points B to D and , denoted by letter E. It is seen that both the value and behaviour of the proton fraction is strongly sensitive to the values of 3BF strengths. some curves are rising in the whole baryon density range while others tend to saturate. The reason is related to the behaviour of the symmetry energy in each case. As is seen in Fig. 8(b), saturates for some values of and and consequently, the proton fraction shows the same behaviour, due to , relation valid for which is satisfied here (see, e.g. 72 ()). The smaller are the symmetry energy and the slope parameter, the smaller is the density at which the proton fraction saturates. However, there exists a region in plane characterized by and MeV (see Fig. 4) within which proton fraction increases monotonously and does not saturate in the whole range of density up to if one chooses any pairs from this area to calculate the energy of stable matter.
Since threebody forces which are one of the key quantities in calculations of the symmetry energy, contribute only to the potential part of , it seems valuable to separate the symmetry energy into its kinetic, , and potential, , parts. Results of such separation is reported in Table III for different pairs at saturation density . The kinetic term of the symmetry energy is, to a very good approximation, equal to the difference in the kinetic energy between PNM and SNM, calculated using a simple Fermi gas model. Therefore, this term gives the same contribution for all parametrizations of 3BF. In the case of , the large contribution comes from the potential part, and this causes both the symmetry energy and its slope not to lie in the experimental range. This problem is resolved for parametrizations B,C,D by considering isospindependent parametrization for the 3BF coefficients. As can be concluded in Fig. 4, the least value of corresponds to point B.
In the LOCV approach one has access to the separate contributions of partial waves energies in the correlated manybody state. Table IV shows the partial wave decomposition of the potential part of , , and up to for different parameters of 3BF. Note that a number of partial waves, namely , , , and , does not contribute to neutron matter. Therefore, contributions from these waves to the symmetry energy are the same for all allowed values of and . For the six remaining partial waves, all quantities , , and decrease (or increase) monotonically from point D to B. The only exception is channel in which increases from D to B while and show opposite behaviour. It can be concluded from this table that the channel gives the major contribution to and . It is also observed that the (spin triplet) and (isospin singlet) channel gives the largest contribution to both symmetry energy and its slope for all values of and . Similar features have been reported in 7 () and references therein.
iv.3 Some astrophysical implications of our results
We expect that for (outer core, see e.g. 72 ()) NS matter is an uniform, electrically neutral plasma in equilibrium. This is also the simplest (minimal) model of the inner core of NS (). Using our 2BF+3BF nuclear Hamiltonian, we calculated the equation of state (EOS) and composition of NS core. Beta equilibrium implies specific fractions of in NS matter constituents, .
Cooling rate of NS core dramatically depends on . Namely, if , where (depending on the precentage of muons), then powerful direct Urca (Durca) process of neutrino emission can act, leading to a fast NS cooling 73 (); 74 ()(for review of neutrino cooling of NS see, e.g. 75 ()). As we see in Fig.9, in the case of our 2BF+3BF Hamiltionian, 3BF decides whether the Durca process is acting in NS core (this was already remarked in 76 ()).
Making 3BF less repulsive is crucial for getting correct values of and . However, this leads to a softer EOS for NS cores. As a consequence, the maximum allowable mass for NS for the 2BF+3BF Hamiltonian, , decreases (in what follows we neglect the effect of rotation on , which for the frequencies of recently discovered pulsars 73 (); 74 () is very small). For all models of our revised UIX 3BF we get , consistently with 77 (); 78 ().
This subsection is but a brief report on the effects of a revised 3BF on NS models. More detailed presentation will be given in a separate paper.
V Conclusion
The extended LOCV formalism including the UIX 3BF force is used in order to study the density dependence of nuclear symmetry energy at zero temperature. The AV18 potential is used as the bare twobody interaction. It is shown that although inclusion of 3BF is successful in reproducing the saturation point properties, namely saturation density and saturation energy as well as the isobaric incompressibility, it fails to predict a reasonable value for both symmetry energy at normal density and its slope. In order to correct this deficiency, a revised version of 3BF is proposed which is based on an isospindependent parametrization of coefficients and in the UIX model. To analyse the effect of the new parametrization, different components of the revised 3BF are compared to those of the old one. Effect of the new parametrization of 3BF on density dependence of symmetry energy, its value at saturation density and its slope and curvature is also studied. It is shown that using this kind of 3BF in the LOCV formalism results in acceptable values for , , and . Moreover, separate contributions of partial waves energies in the potential part of symmetry energy are reported. It is found that the and channel, and particularly the partial wave, give the largest contribution to both symmetry energy and its slope for all values of and .
Implications of our 3BF on NS models were briefly mentioned in Sect. IV.C. A detailed report on the application of our 3BF models to modeling of NS structure and evolution is being prepared.
Acknowledgements.
HRM thank the Research Council of University of Tehran for the grants. This work was supported in part by the National Science Centre, Poland, under the grant No. 2014/13/B/ST9/02621.References
 (1) W. D. Myers and W. J. Swiatecki, Nucl. Phys. A 336, 267 (1980).
 (2) B. A. Brown, Phys. Rev. Lett. 85, 5296 (2000).
 (3) M. Warda, X. Vinas, X. RocaMaza, and M. Centelles, Phys. Rev. C 80, 024316 (2009).
 (4) J. Zenihiro, et al., Phys. Rev. C 82, 044611 (2010).
 (5) X. RocaMaza, M. Centelles, X. Vinas and M. Warda, Phys. Rev. Lett. 106, 252501 (2011).
 (6) B. K. Agrawal, J. N. De, and S. K. Samaddar, Phys. Rev. Lett. 109, 262501 (2012).
 (7) Isaac Vidana, Artur Polls, and Constanca Providencia, Phys. Rev. C 84, 062801(R) (2011).
 (8) Helei Liu, GaoChan Yong, and DeHua Wen, Phys. Rev. C 91, 044609 (2015).
 (9) T. Inakura, and H. Nakada, Phys. Rev. C 92, 064302 (2015).
 (10) X. RocaMaza et al., Phys. Rev. C 92, 064304 (2015).
 (11) P. Danielewicz, and J. Lee, Nucl. Phys. A 922, 1 (2014).
 (12) P. Danielewicz, R. Lacey, and W. G. Lynch, Science 298, 1592 (2002).
 (13) V. Baran, M. Colonna, V. Greco, and M. Di Toro, Phys. Rep. 410, 335 (2005).
 (14) B. A. Li, L. W. Chen, and C. M. Ko, Phys. Rep. 464, 113 (2008).
 (15) ZhaoQing Feng, Phys. Rev. C 83, 067604 (2011).
 (16) ChunWang Ma, Fang Wang, YuGang Ma, Chan Jin, Phys. Rev. C 83, 064620 (2011).
 (17) GaoChan Yong, Phys. Rev. C 84, 014607 (2011).
 (18) Sanjeev Kumar, Y. G. Ma, G. Q. Zhang and C. L. Zhou, Phys. Rev. C 84, 044620 (2011).
 (19) Jianmin Dong, Wei Zuo, and Werner Scheid, Phys. Rev. Lett. 107, 012501 (2011).
 (20) J. M. Lattimer, C. J. Pethick, M. Prakash, and P. Haensel, Phys. Rev. Lett. 66, 2701 (1991).
 (21) A. W. Steiner, M. Prakash, J. M. Lattimer, and P. J. Ellis, Phys. Rep. 411, 325 (2005).
 (22) B. G. ToddRutel and J. Piekarewicz, Phys. Rev. Lett. 95, 122501 (2005).
 (23) L. F. Roberts, G. Shen, V. Cirigliano, J. A. Pons, S. Reddy, and S. E. Woosley, Phys. Rev. Lett. 108, 061103 (2012).
 (24) L. L. Lopes, and D. P. Menezes, Braz. J. Phys. 44; 774 (2014).
 (25) M. Prakash, T. L. Ainsworth, J. M. Lattimer, Phys. Rev. Lett. 61, 2518 (1988).
 (26) W.C. Chen and J. Piekarewicz, Phys. Rev. C 90, 044305 (2014).
 (27) N. Nikolov, N. Schunck, W. Nazarewicz, M. Bender, and J. Pei, Phys. Rev. C 83, 034305 (2011).
 (28) Isaac Vidana, Constanca Providencia, Artur Polls, and Arnau Rios, Phys. Rev. C 80, 045806 (2009).
 (29) Tetsuya Katayama and Koichi Saito, Phys. Rev. C 88, 035805 (2013).
 (30) Kh. Gad, Kh. S. A. Hassaneen, Nucl. Phys. A 793, 67 (2007).
 (31) A. Akmal and V.R. Pandharipande, Phys. Rev. C 56, 2261 (1997).
 (32) H. Togashiv, M. Takano, Nucl. Phys. A 902, 53 (2013).
 (33) Jianmin Dong, Wei Zuo, Jianzhong Gu, and Umberto Lombardo, Phys. Rev. C 85, 034308 (2012).
 (34) S. Gandolfi, J. Carlson, and Sanjay Reddy, Phys. Rev. C 85, 032801(R) (2012).
 (35) Gulshan Mahajan and Shashi K. Dhiman, Phys. Rev. C 84, 045804 (2011).
 (36) LieWen Chen, Phys. Rev. C 83, 044308 (2011).
 (37) Z. W. Zhang, S. S. Bao, J. N. Hu, and H. Shen, Phys. Rev. C 90, 054302 (2014).
 (38) M. B. Tsang, et al., Phys. Rev. C 86, 015803 (2012).
 (39) James M. Lattimer and Yeunhwan Lim, Astrophys. J. 771, 51 (2013).
 (40) BaoAn Li ,Xiao Han, Phys. Lett. B. 727, 276 (2013).
 (41) B.A. Li, A. Ramos, G. Verde, and I. Vidana, eds.,Topical Issue on Nuclear Symmetry Energy, vol. 50 of Euro. Phys. J. A (2014).
 (42) M. Centelles, X. RocaMaza, X. Vinas, and M. Warda, Phys. Rev. Lett. 102, 122502 (2009).
 (43) M. D. Cozma, Y. Leifels, W. Trautmann, Q. Li, and P. Russotto, Phys. Rev. C 88, 044912 (2013).
 (44) R. B. Wiringa, V. G. J. Stoks, and R. Schiavilla, Phys. Rev. C 51 , 38 (1995).
 (45) B. S. Pudliner, V. R. Pandharipande, J. Carlson, and R. B. Wiringa, Phys. Rev. Lett. 74, 4396 (1995); B. S. Pudliner, V. R. Pandharipande, J. Carlson, S. C. Pieper, and R. B. Wiringa, Phys. Rev. C 56, 1720 (1997); S. C. Pieper, V. R. Pandharipande, R. B. Wiringa, and J. Carlson, Phys. Rev. C 64 , 014001 (2001).
 (46) J. C. Owen, R. F. Bishop, and J. M. Irvine, Phys. Lett. B 59, 1 (1975).
 (47) M. Modarres and J. M. Irvine, J. Phys. G : Nucl. Part. Phys. 5, 511 (1979).
 (48) R. V. Reid, Ann. Phys. (NY) 50, 411 (1968).
 (49) A. M. Green, J. A. Niskanen, and M. E. Sainio, J. Phys. G : Nucl. Part. Phys. 4, 1055 (1978).
 (50) H. R. Moshfegh and M. Modarres, Nucl. Phys. A 759, 79 (2005).
 (51) M. Modarres, J. Phys. G : Nucl. Part. Phys. 19, 1349 (1993); H. R. Moshfegh and M. Modarres, Nucl. Phys. A 792, 201 (2007).
 (52) M. Modarres, J. Phys. G : Nucl. Part. Phys. 21, 351 (1995).
 (53) M. Modarres and H. R. Moshfegh, Phys. Rev. C 62, 044308 (2000); Prog. Theor. Phys. 107, 139 (2002).
 (54) S. Zaryouni and H. R. Moshfegh, Eur. Phys. J. A 45, 69 (2010).
 (55) S. Zaryouni, M. Hassani, and H. R. Moshfegh, Phys. Rev. C 89, 014332 (2014).
 (56) S. Goudarzi and H. R. Moshfegh, Phys. Rev. C 91, 054320 (2015).
 (57) S. Goudarzi and H. R. Moshfegh, Phys. Rev. C 92, 035806 (2015).
 (58) P.J. Siemens, Nucl. Phys. A 141, 225 (1970); I. Bombaci, U. Lombardo, Phys. Rev. C 44, 1892 (1991).
 (59) V. Baran, M. Colonna, M. Di Toro, V. Greco, M. ZielinskaPfabe, H. H. Wolter, Nucl. Phys. A 703, 603 (2002).
 (60) L.W. Chen, B.J. Cai, C. M. Ko, B.A. Li, C. Shen, and J. Xu, Phys. Rev. C 80, 014322 (2009).
 (61) B. G. ToddRutel and J. Piekarewicz , Phys. Rev. Lett. 95, 122501 (2005).
 (62) G. A. Lalazissis, J. Konig and P. Ring, Phys. Rev. C 55, 540 (1997).
 (63) K. Sumiyoshi, H. Kuwabara, H. Toki, Nucl. Phys. A 581, 725 (1995).
 (64) T. Niksic, D. Vretenar, P. Finelli and P. Ring, Phys. Rev. C 66, 024306 (2002); G. A. Lalazissis, T. Niksic, D. Vretenar and P. Ring, Phys. Rev. C 71, 024312 (2005).
 (65) J. R. Stone, N. J. Stone, and S. A. Moszkowski, Phys. Rev. C 89, 044316 (2014).
 (66) B.A. Li, L.W. Chen, C.M. Ko, Phys. Rep. 464, 113 (2008).
 (67) L. W. Chen, C. M. Ko, B.A. Li, Phys. Rev. Lett. 94, 032701 (2005); L. W. Chen, C. M. Ko, B.A. Li, Phys. Rev. C 72, 064309 (2005).
 (68) T. Li, et al., Phys. Rev. Lett. 99, 162503 (2007).
 (69) J. Fujita and H. Miyazawa, Prog. Theor. Phys. 17, 360 (1957).
 (70) Steven C. Pieper, V. R. Pandharipande, R. B. Wiringa, and J. Carlson, Phys. Rev. C 64, 014001 (2001).
 (71) S. A. Coon, et al., Nucl. Phys. A 317, 242 (1979).
 (72) C. Ordonez and U. van Kolck, Phys. Lett. B 291, 459 (1992).
 (73) X. R . Zhou, G. F. Burgio, U. Lombardo, H.J. Schulze, and W. Zuo ,Phys. Rev. C 69, 018801 (2004).
 (74) Z. H. Li and H.J. Schulze, Phys. Rev. C 78, 028801 (2008).
 (75) A. Akmal, V. R. Pandharipande, and D. G. Ravenhall, Phys. Rev. C 58, 1804 (1998).
 (76) M. Sharma, S. Rafi, D. Pachouri, and W. Haider, Proceedings of the DAE Symp. on Nucl. Phys. 56, 716 (2011).
 (77) J. Carlson, V. R. Pandharipande, and R. B. Wiringa, Nucl. Phys. A 401, 59 (1983).
 (78) W. Zuo, A. Lejeune, U. Lombardo, and J. F. Mathiot, Nucl. Phys. A 706, 418 (2002).
 (79) Z. H. Li, U. Lombardo, H.J. Schulze, and W. Zuo, Phys. Rev. C 77, 034316 (2008).
 (80) P. Haensel, A. Y. Potekhin, and D. G. Yakovlev, Neutron Stars 1. Equation of state and structure (Springer, New York, 2007).
 (81) J. M. Lattimer, C.J. Pethick, M. Prakash and P. Haensel, Phys. Rev. Lett. 66, 2701 (1991).
 (82) D. G. Yakovlev, C.J. Pethick, Annu. Rev. Astron. Astrophys. 42, 169 (2004)
 (83) D. G. Yakovlev, A. D. Kaminker, O. Y. Gnedin and P. Haensel, Phys. Rep. 354, 1 (2001).
 (84) A. Dehghan Niri, H. R. Moshfegh, P. Haensel, Phys. Rev. C 93, 045806 (2016).
 (85) P. B. Demorest, T. Pennucci, S. M. Ransom, M.S. E. Roberts, J. W. T. Hessels, Nature 467, 1081 (2010).
 (86) J. Antoniadis et al., Science 340, 6131 (2013).