Nuclear matter saturation with chiral three-nucleon interactions fitted to light nuclei properties

Abstract The energy per particle of symmetric nuclear matter and pure neutron matter is calculated using the many-body Brueckner–Hartree–Fock approach and employing the Chiral Next-to-next-to-next-to leading order (N3LO) nucleon–nucleon (NN) potential, supplemented with various parametrizations of the Chiral Next-to-next-to leading order (N2LO) three-nucleon interaction. Such combination is able to reproduce several observables of the physics of light nuclei for suitable choices of the parameters entering in the three-nucleon interaction. We find that some of these parametrizations provide a satisfactory saturation point of symmetric nuclear matter and values of the symmetry energy and its slope parameter L in very good agreement with those extracted from various nuclear experimental data. Thus, our results represent a significant step toward a unified description of few- and many-body nuclear systems starting from two- and three-nucleon interactions based on the symmetries of QCD.


Introduction
The advent of interactions derived in the framework of chiral perturbation theory (ChPT) [1][2][3] opened a new and systematic way to investigate the properties of finite nuclear systems and nuclear matter. The big advantage of using such method lies in the fact that two-body as well as many-body forces can be calculated order by order according to a well defined scheme based on a lowenergy effective QCD Lagrangian which retains the symmetries of QCD, and in particular the approximate chiral symmetry. Within this approach the details of the QCD dynamics are contained in parameters, the so-called low-energy constants (LECs), which are fixed by low-energy experimental data. This systematic procedure is particularly useful for nucleonic systems where the importance of the three-nucleon force (TNF) is a well established feature. It is indeed well known that high precision nucleon-nucleon (NN) potentials, fitting NN scattering data up to energy of 350 MeV, with a χ 2 per datum next to 1, underestimate the experimental binding energies of 3 H, 3 He by about 1 MeV and that of 4 He by about 4 MeV [4]. This missing binding energy can be accounted for by introducing a TNF into the nuclear Hamiltonian [4]. NN potentials plus TNFs based on ChPT have been recently used to investigate properties of medium-mass nuclei [5,6] and heavy nuclei [7]. A very important task in this line is the evaluation of the uncertainties originating in the nuclear Hamiltonian [8] and in particular on the LECs and to establish which should be the best fitting procedure to fix them [9]. For example, in Ref. [5] a simultaneous optimization of the NN interaction plus a TNF in light and medium-mass nuclei has been performed.
In the present work, we use nuclear interactions derived within the framework of ChPT to calculate the equation of state (EoS) of symmetric nuclear matter and pure neutron matter using the Brueckner-Hartree-Fock (BHF) approach [10,11]. In particular, we intend to study nuclear matter properties using different sets of TNFs constructed to reproduce specific observables in the threeand four-nucleon systems. In this way we can estimate the sensitivity introduced by the different LECs' values.
In a previous paper [12], we have investigated whether using the same interactions at two-and three-body level, it was possible to concurrently reproduce properties of finite light nuclei and nuclear matter. Fixing the parameters of the TNF to simultaneously describe the 3 H, 3 He and 4 He binding energies and the neutrondeuteron (n-d) doublet scattering length [13], we found that none of the considered interactions was able to reproduce a good saturation point of symmetric nuclear matter. In [12] we used the Argonne V18 NN potential [14]    TNFs. In the first case we employed various parametrizations of the Tucson-Melbourne TNF [15], while in the second case, we used the local version of the chiral N2LO [16] TNF (hereafter N2LOL) [17]. The N2LOL differs from the usual N2LO TNF only for the use of a local cutoff (see Refs. [13] and [12]) for other details. In the present work, we consider an interaction fully based on ChPT both for the two-and three-nucleon sectors. We use indeed the chiral N3LO [18] NN potential in conjunction with different parametrizations of the N2LOL three-nucleon interaction [17]. As the N3LO NN interaction has been produced in versions with different cutoff values, we intend to study the sensitivity of our results on the cutoff. Partially this dependence could arise from the fact that while the NN potential is calculated at order N3LO of chiral perturbation expansion, the three-nucleon one is calculated at order N2LO. Chiral TNFs have also been calculated at order N3LO [19]. Their effect in the description of light nuclei is at present under investigation [20]. It should be noticed that no additional low energy constant (LEC) appears at this order. New LECs appear in short range subleading contributions to the TNF at order N4LO [21] and their contribution seems to be potentially important [22]. The long range part of TNF at the same order contains no additional LEC [23]. In this work we limit our study considering TNFs calculated at order N2LO. Nuclear interactions based on ChPT have been used by several groups for calculating the EoS of pure neutron matter [24][25][26][27][28][29][30] and symmetric nuclear matter [31][32][33][34][35]. A comparison with some of these calculations will be performed in last part of this work.
The paper is organized as follows: in section 2 we review the parameters of the N2LOL three-nucleon force and the determination of the low energy constants; in section 3 we briefly discuss how to include a TNF in the Brueckner-Hartree-Fock (BHF) approach; section 4 is devoted to show and discuss the results of our calculations; finally in section 5 we outline the main conclusions of the present study.

The N2LOL three nucleon force
Following Ref. [17], the N2LOL potential can be written as function of the transferred momenta q and q of the external nucleons. We adopted [17] a cut off function of the form F = e  [18]. For the case = 450 MeV, we have considered two additional parametrizations denoted as N2LOL6 and N2LOL7. In this case c 1 , c 3 and c 4 have been taken from Ref. [34]. The values of c i are reported in the caption of Table 1 for the two considered values of the momentum cutoff . The values of the remaining LECs c E and c D are shown in Table 1. The parametrization N2LOL1, is the original one proposed in Ref. [17], where c E and c D were fitted to reproduce the binding energies of 3 H and 4 He in a no-core shell model calculation. In Ref. [36], fixing the range of c D between −3 and 3, the parameter c E was determined fitting the binding energies of 3 H and 3 He. In this way the authors of Ref. [36] obtained two curves c E (c D ) fulfilling the previous constraints. As the two curves turned out to be very close, the authors of Ref. [36] performed an average between them. Finally, for each set of parameters (c E , c D ), the Gamow-Teller (GT) matrix element of tritium β-decay was calculated and, using the corresponding experimental value and its error-bars, the minima and a maxima values for c D and c E satisfying this requirement were determined. In the parametrizations N2LOL2 and N2LOL3, we have adopted the minima and the maxima values allowed for c D and c E according to the construction described above. The parametrization N2LOL4, taken again from Ref. [34], was obtained in the same way as the N2LOL2 and N2LOL3 ones but allowing that the GT matrix element of tritium β-decay to be reproduced with a slight larger uncertainty (however less than 1%). The last parametrization (N2LOL5) that we have considered, has been obtained fixing the couple (c D , c E ) on the trajectory reproducing the binding energies of 3 H and 3 He and requiring to get the best saturation density ρ 0 of symmetric nuclear matter. As we have said before, for = 450 MeV we have explored just two parametrizations of the N2LOL TNF, namely N2LOL6 and N2LOL7. Concerning the parametrization N2LOL6, we have fixed the values of the low energy constants c D and c E in the same way of N2LOL4 in the case = 500 MeV; for parametrization N2LOL7, we have optimized the symmetric nuclear matter saturation point with the constraint to reproduce the binding energies of 3 H and 3 He. Finally, we underline that in the two-body N3LO interaction we have adopted the same value of the cut off used for the N2LOL TNF. Thus the LECs (see Table 1 caption) have been consistently calculated for the two different values of [34].

Inclusion of three body forces in the BHF approach
As discussed in the literature [37][38][39][40], three-nucleon forces cannot be included in the BHF formalism [11,41,42] in their original form. This would require the solution of three-body scattering equations in the nuclear medium (Bethe-Faddeev equations) which is an extremely demanding task from a technical and computational point of view. A possible simplification is to build an effective density dependent two-body force starting from the original three-body one by averaging over the coordinates (spatial, spin and isospin) of one of the nucleons. The effective NN force due to the general three-nucleon force W (1, 2, 3) can be written as [37,43]: where we have defined dx 3 = Tr (τ 3 ,σ 3 ) dr 3 and n(1, 2, 3) is the density distribution of nucleon 3 in relation to nucleon 1 at r 1 and nucleon 2 at r 2 . The function n(1, 2, 3) contains the effect of the NN short-range correlations suppressing the probability of finding two nucleons close to each other. In the following, we assume that this density distribution can be factorized as [37,12] n(1, where ρ is the nucleon density, g(1, 3) and g(2, 3) are the correlation functions between the nucleons (1, 3) and (2, 3) respectively.
The latter quantities can be written as g(1, 3) is the so-called defect function (and similarly for g (2,3)). The defect function η(i, j) is the difference between the correlated and the uncorrelated wave functions and thus represents a measure of short range correlations [44,45]. To simplify the numerical calculations and following [37,12], in the present work we use central correlation functions g(i, j) independent on spin and isospin. Moreover, it has been shown [37][38][39] that this central correlation functions, in which are included the main contributions of the 1 S 0 and 3 S 1 channels, are weakly dependent on the density, and in principle can be approximated by a Heaviside step function θ(r ij − r c ), with r c = 0.6 fm in all the considered density range.
However in the present work η(i, j) has been calculated in a selfconsistent way at each step of the iterative procedure from the Brueckner G-matrix [44]. Finally P ij are the spin, isospin and space exchange operators between nucleon i and j.
A different average procedure, with respect to the one discussed above, has been performed in Ref. [46], where an in medium effective NN potential was derived from the N2LO three-nucleon force. In [46] just the on-shell contributions have been determined while the off-shell counterparts have been obtained by extrapolation.
We want to stress that the presence of contributions arising from the cyclic permutations of W (1, 2, 3) and from the exchange operators P ij in the average Eq. (1) has been neglected in the BHF calculations reported in [32], where no saturation of nuclear matter has been obtained using the N3LO NN potential plus the N2LOL three nucleon interaction. The presence of the exchange operators in Eq. (1) has been neglected in our previous work [12]. As we will show in the following, these contributions to the TNF play a very important role for the nuclear matter saturation mechanism. In our present calculation scheme, we take into account of all internal and external permutations of the TNF. The resulting effective density dependent potential contains in particular some purely repulsive contributions which are missing when neglecting the P ij exchange operators. Such repulsion is needed to contrast the strong attraction provided at two-body level by the N3LO potential as we will show in next section.

Results and discussions
We now present the results of our calculations of the energy per particle E/ A of symmetric nuclear matter and pure neutron matter using N3LO NN potential supplemented with the N2LOL three-body force. These calculations have been performed using the BHF approximation of the Brueckner-Bethe-Goldstone holeline expansion [10,11] with the BHF continuous choice [47,42] for the auxiliary single particle potential U (k) entering in the selfconsistent procedure to determine E/ A. As suggested in Ref. [47] and later on numerically confirmed in Ref. [48,49], the BHF continuous choice for U (k) optimizes the convergence of the hole-line expansion thus illustrating the accuracy of the calculated energy per baryon at the two-hole line approximation with this prescription for th auxiliary potential.
In Fig. 1 we show the energy per particle of pure neutron matter (left panel) and symmetric nuclear matter (right panel) using the cutoff = 500 MeV. The dotted lines in both panels refer to the calculation performed employing the N3LO two-body potential without any TNF. First we note the sizeable contribution provided by the TNF to the energy per particle of both symmetric nuclear matter and pure neutron matter. In the case of symmetric nuclear matter the saturation point moves indeed from 0.41 fm −3 to values between 0.16 fm −3 and 0.185 fm −3 , depending on employed TNF parametrization (see Table 1). More specifically, the parametrization N2LOL1 predicts a satisfactory saturation point for symmetric nuclear matter at a density of 0.185 fm −3 and an en-  We note that the parametrization N2LOL5 has been constructed to reproduce the 3 H and 3 He binding energies and the best saturation density of symmetric nuclear matter. However, we want again to remark that in this case the GT matrix element of tritium β-decay is not reproduced. The curves for the energy per nucleon of pure neutron matter (right panel in Fig. 1) are very similar in all the density range considered.
In Fig. 2, we show the same quantities of Fig. 1, but adopting the cutoff = 450 MeV. In this case we have considered the two parametrizations N2LOL6 and N2LOL7 (see Table 1) of the TNF.
The results obtained for = 450 MeV are very similar to the ones However the parametrizations N2LOL1 and N2LOL7, which reproduce 3 H and 3 He binding energies, produce also a satisfactory saturation point for symmetric nuclear matter. Such achievements represent a big improvement of our previous calculations [12].
In the case of asymmetric nuclear matter with neutron density ρ n , proton density ρ p , total nucleon density ρ = ρ n + ρ p and asymmetry parameter β = (ρ n − ρ p )/ρ the energy per nucleon in the BHF approach can be accurately reproduced [50] using the socalled parabolic approximation where E sym (ρ) is the nuclear symmetry energy. Thus the symmetry energy can be calculated as the difference between the energy per particle of pure neutron matter (β = 1) and symmetric nuclear matter (β = 0).
The nuclear symmetry energy, calculated with this prescription, is plotted in Fig. 3 for the parametrizations of the N2LOL TNF with = 500 MeV (left panel) and = 450 MeV (right panel) used in the present work. In the same figure, we show E sym (triangles) as obtained from recent calculations [51] of asymmetric neutron-rich matter with two-and three-body chiral interactions. The results of Ref. [51] have confirmed the validity of the quadratic approximation (Eq. (3)) for describing the EOS highly asymmetric matter.
However, it has been recently shown [52,53] that the β 4 term in the energy per nucleon of asymmetric nuclear matter could not be negligible, especially at supranuclear densities, thus having sizeable influence e.g. on neutron star cooling [54].
To compare our results with the value of the symmetry energy extracted from various nuclear experimental data [55,56], we report in Table 2 and Table 3 the symmetry energy E sym (ρ 0 ) and the so-called slope parameter   Table 3 The same of Table 2   at the calculated saturation density ρ 0 (second column in Table 2 and 3) for the TNF parametrizations considered in this work. As we can see our calculated E sym (ρ 0 ) and L are in good agreement with the values extracted from experimental data [55]: E sym (ρ 0 ) = 29.0-32.7 MeV, and L = 40.5-61.9 MeV. This agreement is lost when the TNF is not included in the calculations (see e.g. the last row in Table 2). Notice that the value of the nuclear incompressibility K 0 of symmetric nuclear matter, at the calculated saturation density, is generally quite low when compared with the value deduced from measured isoscalar giant monopole resonances in medium-heavy nuclei K 0 = 240 ± 20 MeV [57] K 0 . Our calculated value of K 0 ranges between 150 and 170 MeV, depending on the TNF model. Let us now compare our results with those of similar calculations present in literature based on chiral nuclear interactions. Our present results are in good agreement with the recent nonperturbative quantum Monte Carlo calculations [26] of pure neutron matter. The authors of Ref. [34] performed nuclear matter calculations up to third order many-body perturbation theory, adopting the effective density dependent NN force derived in Ref. [46]. These authors found results in acceptable agreement with ours both for nuclear matter and pure neutron matter. Several nuclear matter calculations based on chiral interactions have been performed using other many-body techniques. In Ref. [31], using the V lowk -approach and the similarity renormalization group (SRG), it was shown that for a suitable choice of the cutoff V lowk it was possible to reproduce a good saturation point of symmetric nuclear matter keeping for the parameters of the TNF the values fixed in few-body calculations. In Ref. [33], using the self-consistent Green's functions method, it has been found that the saturation point of symmetric nuclear matter is strongly improved with the help of the N2LOL three-nucleon force although, also in this case, the saturation energy results underestimated.
In Fig. 4 we compare our results (using the parametrizations N2LOL1, N2LOL4 and N2LOL6) for the energy per particle of pure neutron matter with those obtained by other researchers using different many-body approaches. The black-and red-dashed regions in the left panel of Fig. 4 represent the results of the many-body perturbative calculations of Ref. [25] using complete two-, threeand four-body interactions at the N3LO of the ChPT. In particular, the region between the two short-dashed black curves (blackdashed band partially overlapped by the red-dashed band) is relative to the N3LO Entem-Machleidt [18] NN potential, whereas the red-dashed band refers to the Epelbaum-Glockle-Meißner (EGM) Fig. 4. Comparison of pure neutron matter energy per particle with other manybody methods (see text for details).
potentials [58]. The width of the bands represents the uncertainties related to the values of the LECs and the cutoff of the threeand four-body forces. The dash-dotted (blue) curves in both panels in Fig. 4, represent the results [29] of an auxiliary field diffusion Monte Carlo (AFDMC) calculation of neutron matter using a local form of two-and three-body chiral interactions at N2LO and two different values of the NN cutoff (see [29] for more details). The green-dotted line in both panels of Fig. 4, corresponds to the results of Ref. [27] with the auxiliary field quantum Monte Carlo (AFQMC) obtained with chiral N3LO two-body force plus N2LO TNF. Finally, the red-dashed band in the right panel of Fig. 4 represents the results of Ref. [30] for the energy per particle of neutron matter after a renormalization group (RG) evolution of the N3LO Entem-Machleidt [18] NN interaction to low-momentum scale ( = 394.6 MeV) and including N2LO TBF (the width of the band reflecting the uncertainties of the LECs in the TNF). As one can see, our results are in very good agreement with all the other calculations considered in Fig. 4, except with the calculations of Ref. [29] for densities close to the empirical saturation density. In fact, in this density region the AFDMC curves are "flat" compared to other calculations [29], thus implying low values, in the range L = (16.0-36.5) MeV, for the slope parameter calculated using the parabolic approximation Eq. (3) at the empirical saturation density ρ 0 = 0.16 fm −3 .

Conclusions
We have performed BHF calculations of the EoS symmetric nuclear matter and pure neutron matter considering the N3LO NN potential plus the N2LOL three nucleon interaction. The last one has been reduced to an effective density dependent two-body interaction averaging over the coordinates (spatial, spin and isospin) of one of the nucleons. We have shown that the contributions arising from the cyclic permutations of W (1, 2, 3) and from the nucleon exchange operators P ij in the average (see Eq. (1)) of the TNF play a decisive role for the nuclear matter saturation mechanism. In fact, no saturation is obtained [32] when these contributions are not taken into account. We want to stress that the parameters of the TNF fitted in calculations of light nuclei have been kept fixed in our nuclear matter calculations.
From one side we found that it was not possible to reproduce simultaneously the binding energy of 3 H, 3 He, the GT matrix element of tritium β-decay and the empirical saturation point of symmetric nuclear matter. However some trends can be extracted, in particular, those parametrizations that reproduce the binding energy of 3 H, 3 He and the GT matrix element of tritium β-decay predict a reasonable value (0.15 fm −3 ) of the saturation density, but the corresponding binding energy per nucleon (B/ A = −E/A) is underestimated. The two parametrizations, N2LOL1 and N2LOL7, that are not constrained to reproduce the GT matrix element of tritium β-decay are able to reproduce quite well the empirical saturation energy but at a too large value of ρ 0 . Finally we have observed a reduced dependence on the momentum cutoff parameter. We have found in addition a good agreement with recent AFQMC [27] and MBPT [25] calculations of neutron matter.
These encouraging results represent a significant step toward a unified description of few-and many-body nuclear systems starting from two-and three-nucleon interactions derived in the framework of chiral perturbation theory.