Impact of chiral hyperonic three-body forces on neutron stars

We study the effect of the nucleon-nucleon-lambda (NN$\Lambda$) three-body force on neutron stars. In particular, we consider the NN$\Lambda$ force recently derived by the J\"{u}lich--Bonn--Munich group within the framework of chiral effective field theory at next-to-next-to-leading order. This force, together with realistic nucleon-nucleon, nucleon-nucleon-nucleon and nucleon-hyperon interactions, is used to calculate the equation of state and the structure of neutron stars within the many-body non-relativistic Brueckner-Hartree-Fock approach. Our results show that the inclusion of the NN$\Lambda$ force leads to an equation of state stiff enough such that the resulting neutron star maximum mass is compatible with the largest currently measured ($\sim 2\ M_\odot$) neutron star masses. Using a perturbative many-body approach we calculate also the separation energy of the $\Lambda$ in some hypernuclei finding that the agreement with the experimental data improves for the heavier ones when the effect of the NN$\Lambda$ force is taken into account.


Introduction
The importance of taking into account nucleon-nucleon-nucleon (NNN) interactions in finite nuclei as well as in infinite nuclear matter is nowadays a well established feature. It is well known that high precision nucleon-nucleon (NN) potentials, which fit NN scattering data up to an energy of 350 MeV with a χ 2 per datum close to 1, underestimate the experimental binding energies of 3 H and 3 He by about 1 MeV, and that of 4 He by about 4 MeV [1]. This missing binding energy can be accounted for by introducing a NNN interaction into the nuclear Hamiltonian [1]. Three-nucleon forces are also crucial for nuclear matter calculations. As it is known, saturation points obtained using different NN potentials with non-relativistic manybody approaches lie within a narrow band [2,3,4], the so-called Coester band, with either a too large saturation density or a too small binding energy compared to the empirical value. Threenucleon forces allows to reproduce properly the empirical saturation point. It has been pointed out that for similar reasons, three-body forces involving hyperons (NNY, NYY and YYY) may play also an important role to describe accurately the properties of neutron stars with hyperons [5,6,7,8,9,10,11] and hypernuclei [12,13,14,15,16,17,18,19,20].
One of the longstanding open problems in nuclear physics and astrophysics that could be (if not completely at least partially) solved with the help of hyperonic three-body forces (YTBF) is the so-called "hyperon-puzzle" of neutron stars [21,22], i.e., the difficulty to reconcile the measured masses of neutron stars with the presence of hyperons in their interiors. Hyperons are expected to appear at 2 − 3 times normal nuclear saturation density (n 0 = 0.16 fm −3 ). At such densities, the neutron and proton chemical potentials are large enough to make the conversion of nucleons into hyperons energetically favorable. This conversion, however, produces a strong softening of the equation of state (EoS), due to the release of the Fermi pressure of the system, which leads to a decrease of the maximum neutron star mass predicted by theoretical models. In many microscopic calculations [23,24,25,26,27,28], this decrease is so large that the maximum mass obtained is not compatible with the current largest measured neutron star masses of ∼ 2M ⊙ [29,30,31,32]. Most of these microscopic calculations have been performed using NN, NNN and NY interactions and, in some cases, also the YY one [23,28]. Just a few full consistent calculations including YTBF are present in literature. The authors of the present work, for instance, in Ref. [7] used a model based on the Brueckner-Hartree-Fock (BHF) approach of hyperonic matter using the Argonne V18 [33] NN force and the Nijmegen NY soft-core NSC89 [34] one supplemented with additional simple phenomenological density-dependent contact terms to establish numerical lower and upper limits to the effect of the YTBF on the maximum mass of neutron stars. Assuming that the strength of these forces was either smaller than or as large as the pure nucleonic ones, the results of that work [7] showed that although the employed YTBF stiffened the EoS, they were, however, unable to provide the repulsion needed to make the predicted maximum masses compatible with the recent observations of massive neutron stars. A multi-Pomeron exchange potential (MPP) model to introduce universal three-body repulsion among three baryons in the hyperonic matter EoS was proposed in Refs. [8,9,10]. This universal three-body repulsive potential was based on the extended soft core (ESC) baryonbaryon interaction of the Nijmegen group [35,36]. The strength of the MPP was determined by analyzing the nucleus-nucleus scattering with the use of a G-matrix folding potential derived from the ESC interaction complemented with the MPP and a three-nucleon attractive part, added phenomenologically in order to reproduce the nuclear saturation properties. The results of those works [8,9,10] showed that when the MPP contribution was taken into account universally for all baryons, a maximum mass of ∼ 2.2M ⊙ was obtained, in contradiction with the results and conclusions of Ref. [7] where the case of a universal three-body repulsion was also analyzed. Finally, in Ref. [11] a Monte Carlo calculation of pure neutron matter with a non vanishing Λ-hyperon concentration was carried out including NN, NNN, NΛ and NNΛ forces. In particular the NNΛ force used in that work was tuned in order to provide a reasonable description of the measured Λ separation energy of several hypernuclei. The authors of Ref. [11] concluded that, with the model they considered, the presence of hyperons in the core of neutron stars could not be satisfactory established and, consequently, according to these authors, there is no clear incompatibility with astrophysical observations when Λs are included. However, one should note, that the presence of protons, necessary to establish the correct β-equilibrium inside neutron stars and thus a proper treatment of nuclear matter, was neglected in their calculation. Although at present there is not yet a general consensus regarding the role played by YTBF in the solution of the hyperon puzzle, it seems that even if they are not the full solution most probably they can contribute to it in an important way.
The aim of the present paper is to perform a calculation of the EoS and structure of neutron stars with non vanishing Λ-hyperon concentrations in the framework of non-relativistic BHF approach (see e.g., Ref. [37,38]) using realistic NN, NNN interactions derived in chiral effective field theory (χEFT) supplemented by NΛ and NNΛ interactions. In particular, for the two-body NN interaction we use the local chiral potential presented in Ref. [39] at next-to-next-to-next-to-leading order (N3LO) which includes the ∆(1232) isobar in the intermediate state of the NN scattering. Regarding the NNN force, we use the potential derived in Ref. [40] calculated at the next-to-next-toleading-order (N2LO) of chiral perturbation theory in the local version as reported in Ref. [41,42]. We note that in this NNN force the possibility of the ∆-excitation is also taken into account. The low energy constants of the NNN have been fitted as in Ref. [42] where it was shown that a good description of nuclear matter can be achieved using this setting. These interactions have been recently employed in Ref. [43] to calculate the β-stable EoS of nuclear matter and the structure of neutron stars. It was found a neutron star maximum mass of 2.07M ⊙ in agreement with the largest measured neutron star masses. The resulting EoS has been also recently used in Ref. [44] to simulate the merging of two equal mass neutron stars.
In the present work we want to study how the finding of Ref. [43] changes allowing for the possible presence of Λ-hyperons in core of neutron stars. Although there is a vast number of NN and NNN interactions derived so far in χEFT, NY and YY interactions have been constructed only by the Jülich-Bonn-Munich group within this framework [45,46,47], developing first the NY ones at leading order (LO) [45] and next-to-leading order (NLO) [46], and then the YY constructed also at NLO [47]. Since in the nucleonic sector we employ interactions calculated in χEFT, for consistency it would be appropriate to use hyperonic interactions derived in the same framework. Unfortunately, however, at present we do not have at our disposal the NΛ interaction presented in Refs. [45,46] and, therefore, in this work we employ instead the NΛ meson-exchange interaction derived by the Nijmegen group in Refs. [48,49]. We are aware that this represents the weakest point of this work that, however, we will try to solve in the future.
Finally, concerning the YTBF, we use the NNΛ force recently derived also by the Jülich-Bonn-Munich group in the framework of χEFT [50]. In the next section we give some additional details on this force and how it is included in our BHF approach. The results of the calculation are then shown and discussed. Conclusions are given at the end.

The NNΛ interaction
As mentioned before, the construction of general three-baryon interactions within the framework of χEFT has been carried out by Petschauer et al., in Ref. [50]. The authors of this work have shown that the first contributions to NNY interactions in χEFT appear at N2LO. They considered the contributions of three different classes of irreducible diagrams: three-baryon contact terms, one-meson exchange and two-meson exchange (see Fig.  1 of Ref. [50]). The pion-exchange mechanism is expected to be the dominant one while contributions coming from heavier meson exchanges (like K or η mesons) can be effectively absorbed into the contact terms. In addition, these authors in Ref. [51] have derived an effective density dependent NΛ interaction by averaging over the coordinates of one of the two nucleons. This allows a straightforward inclusion of the NNΛ force in the BHF approach. The strategy is formally identical to the one adopted for the inclusion of the NNN interaction (see e.g., Refs. [52,53] for details). A more delicate point that deserves some comments concerns the setting of the low energy constants (LECs) in the NNΛ interaction. In the pure nucleonic sector it is possible to use different choices for fixing the values of the LECs. Usually the LECs of the NNN interaction are fixed to reproduce the binding energy of light ( 3 H, 3 He and 4 He) nuclei but also other choices are possible. In the hypernuclear sector the situation is more complicated due to the lack of enough experimental data and to the few existing ab-initio calculations in light hypernuclei [20,45,46,54,55]. In Ref. [51] in order estimate the values of the LECs, the baryon decouplet has been introduced as effective degree of freedom. Then a minimal non-relativistic Lagrangian has been constructed and the LECs have been estimated through decouplet saturation [51]. In this approximation only one LEC, denoted as H ′ in Ref. [51], remains as a free parameter. According to dimensional analysis in Ref. [51] it has been considered H ′ = ±1/ f 2 π being f π = 93 MeV the pion decay constant. In the present, work we consider H ′ = β/ f 2 π where β is a rescaling parameter that we fix in order to reproduce the single-particle potential at zero momentum, U Λ (0), of the Λ-hyperon in symmetric nuclear matter at saturation density. The empirical value of U Λ (0) is obtained from the extrapolation to infinite matter of the binding energy of the Λ in hypernuclei, and it is found to be in the range [−30, −28] MeV [56]. In this work we consider the two extreme values of this interval to determine the parameter β. Hereafter we refer to NNΛ 1 and NNΛ 2 to the models in which U Λ (0) has been respectively taken equal to −28 MeV and −30 MeV to fix the value of β. In order to regularize the short range part of the NNΛ interaction, following [51], we have employed a non local regulator of the form e −(p 4 +p ′4 )/Λ 4 with a cut-off Λ = 500 MeV.
Concerning the NΛ interaction, as stated before, we have used the Nijmegen Soft-Core 97 (NSC97) meson-exchange NY interaction [48,49]. We note that the NSC97 NY force has been provided in 6 different versions (NSC97a-f) according to the value of the magnetic vector α m V = F/(F + D) ratio. For simplicity, in this work we have considered as representative cases of the NΛ interaction the models NSC97a and NSC97e. Results for the other NSC97 models are qualitatively similar. Note that the parameter β, reported in Tab. 1, is in fact fixed for each set of NΛ and NNΛ interaction models.

Results and discussions
Before analyzing the effect of the NNΛ interaction on neutron stars, it is interesting to consider first how this interaction affects the hypernuclear structure. To such end, we calculate the separation energy of the Λ hyperon in a few hypernuclei. This quantity is simply the difference between the total binding energies of an ordinary nucleus A Z and the corresponding hypernucleus A+1 Λ Z. To determine it we follow a perturbative manybody approach to calculate the Λ self-energy in finite nucleus which is then used to obtain, by solving the Schrödinger equation, the energies and wave functions of all the single-particle bound states of the Λ in the nucleus. A detail description of this method is given in Refs. [38,57,58,59]. Results for the Λ separation energy in 41 Λ Ca, 91 Λ Zr and 209 Λ Pb are shown in Tab. 2. We note that for technical reasons we have considered only hypernuclei that are described as a closed shell nuclear core plus a Λ sitting in a single-particle state. Unfortunately, experimental data does not exists for the three hypernuclei considered and for comparison we have taken the closest representative ones for which experimental information is available. Data have been taken from Tab. IV of Ref. [60]. Note that in the case of 91 Λ Zr and 209 Λ Pb hypernuclei, the inclusion of the NNΛ interaction clearly improves the agreement of the theorerical calculation with the experimental data. This is, however, not the case of the 41 Λ Ca where the two models of the NNΛ interaction predict too much repulsion and none of them is able to provide a value of the Λ separation energy in agreement with the experimental one. The same trend has been observed in lighter hypernuclei. We recall, however, that all the finite hypernuclei results shown in Tab. 2 have been obtained without any refitting of the parameter β, and that a better agreement with the experimental data for the lighter hypernuclei could in principle be obtained by readjusting this parameter individually to each hypernucleus. However, a detailed study of the effect of the NNΛ interaction   on hypernuclei is out of the scope of the present work, and it is left for the future. Let us now consider the effect of the YTBF on neutron stars. In Fig. 1 we show first the single particle potential, U Λ (k), of the Λ-hyperon in symmetric nuclear matter as function of the single particle momentum k at saturation density. Results for the NSC97a (NSC97e) NΛ interaction are presented in the left (right) panel together with those including the effect of the NNΛ force. Note that both the NSC97a and NSC97e models (with no NNΛ interaction) predict a value of the Λ singleparticle potential at zero momentum of about −40 MeV, much lower than the empirical value extrapolated from hypernuclear data [56]. Note also that the NSC97a model predicts more attraction than the NSC97e one over the whole range of momenta. This does not change adding the NNΛ force. The repul- sive effect of the NNΛ interaction is even more clear looking at Fig. 2 where U Λ (0) in symmetric nuclear matter is shown as function of the baryonic density n B . Note that U Λ (0) is very deep when only two body interactions are considered and it shows a minimum located at n B ∼ 0.4 fm −3 and n B ∼ 0.3 fm −3 for the NSC97a and NSC97e models, respectively. The inclusion of the NNΛ interaction induces repulsion for densities larger than about 0.1 fm −3 and it shifts this minimum to a value of the density around 0.16 fm −3 for all the models considered. As expected, the effect of YTBF is almost negligible in the low density region.
In order to perform the calculation of the β-stable neutron star matter EoS one has to find for each value of the total baryonic density n B = n n + n p + n Λ the values of the particle concentration Y i = n i /n B that fulfill the chemical equilibrium equa-tions: µ n − µ p = µ e , µ n = µ Λ , µ e = µ µ . (1) Note that, besides nucleons and leptons, we have considered here only the Λ and have ignored the possible appearance of other hyperons. The reason is that this is a first exploratory work where we are just interested on the role of the NNΛ force. A more complete study of the effect of YTBF in neutron stars requieres, of course, the inclusion of the other hyperon species and their interactions. This, however, is left for a future work. In addition, the charge neutrality condition, n p = n e +n µ , should hold. In these equations µ i and n i are, respectively, the chemical potential and number density of the i-th species. The chemical potential is calculated according to the usual thermodynamical relation: µ i = ∂ǫ ∂n i where ǫ is the energy density. The composition of β-stable neutron star matter is shown in the left panel of Fig. 3 for the models NSC97a and NSC97a+ NNΛ 1 . Qualitatively similar results are obtained for the other models which are not shown for simplicity. The continuous lines show the results when only NΛ, in addition to NN and NNN forces, are taken into account whereas the dashed ones include also the contribution of the NNΛ force. The effect of the latter is twofold. First it shifts the onset of the Λ-hyperon to slightly larger baryonic densities. The second effect, maybe the most important one, is that the NNΛ force strongly reduces the abundance of Λ particles at large baryonic densities with the consequent stiffening of the EoS compared to the case in which the NNΛ force is not included, as it can be seen in the right panel of the figure, where the total pressure P is show as a function of the total energy density ε. Consequently, the mass of the neutron star, and in particular its maximum value, increases. This is shown in Fig. 4 where it is plotted the mass-radius relation for the models NSC97a and NSC97e with and without the inclusion of the NNΛ force obtained by solving the well known Tolman-Oppenheimer-Volkoff equations. The black line corresponds to the case of pure nucleonic matter shown as a reference. It is remarkable that the maximum masses obtained including the NNΛ force are compatible with the largest measured masses of ∼ 2M ⊙ [29,30,31,32]. This is in agreement with the calculation performed in Ref. [11]. Notice that the result of our present calculations are based on a more realistic description of neutron star matter compared to the one given in Ref. [11] (pure neutron matter plus a finite concentration of Λ hyperons). In addition, we use more realistic interactions both in the nucleonic and the hyperonic sectors than the ones used in Ref. [11]. Note also that, although the concentration of the Λs is strongly reduced due to the effect of the NNΛ force, they are still present in the interior of a 2M ⊙ neutron star. This differs from what is concluded in Ref. [11] where it was found that the only NNΛ force able to produce a EoS stiff enough to support maximum masses compatible with the recent observation of 2M ⊙ neutron stars lead to the total disappearance of Λ hyperons in the core of these objects.
The neutron star properties, mass, radius and central baryonic density, for the maximum mass configuration are summarized in Tab. 3. Note that models which do not account for the NNΛ interaction provide very low neutron star maximum masses between 1.3 − 1.5M ⊙ . This is in agreement with several calculations performed by various research groups using different many-body methods [23,24,25,26,27]. 10

Conclusions
We have studied the effects of a hyperonic NNΛ force derived by the Jülich-Bonn-Munich in χEFT at N2LO [50] in neutron stars and some single-Λ hypernuclei. We have calculated the EoS and structure of neutron stars within the many-body BHF approach using in addition to the NNΛ force realistic NN, NNN and NΛ interactions. In particular, we have used the chiral NN and NNN interactions derived by Piarulli et al., and Epelbaum et al. in Refs. [39] and [40], respectively. For the NΛ, instead, we have employed the NSC97a and NSC97e models developed by the Nijmegen group within the framework of meson-exchange theory in Refs. [48,49]. The reason for the use of this NΛ interaction is simply the fact that we do not have presently at our disposal the chiral NΛ interaction derived by the Jülich-Bonn-Munich group in Refs. [45,46,47]. This represents a weak point of the present work that, however, we will try to solve in the future. After adjusting the NNΛ force to reproduce the binding energy of the Λ-hyperon in symmetric nuclear matter at saturation density, we have calculated the Λ separation energy in 41 Λ Ca, 91 Λ Zr and 209 Λ Pb. We have found that whereas the agreement between the calculated separated energy and the experimental data improves in the case of the heavier nuclei when the effect of the NNΛ is included, this force results to be too much repulsive in the case of 41 Λ Ca and the lighter hypernuclei. We note, however, that all the finite hypenuclei results were obtained without refitting the NNΛ force and that, a better agreement with experimental data for the lighter hypernuclei could be found if the force is adjusted individually to each hypernucleus. Finally, we have calculated the neutron star composition and EoS and have determined the maximum mass predicted by the different models considered. Our results have shown that when the NNΛ force is included, the EoS becomes stiff enough such that the resulting maximum mass is compatible with largest measured neutron star maxi-  Table 3. Neutron star properties, mass (M max ), radius (R) and central baryonic density (n c ), for the maximum mass configuration for the different models considered. Results for a pure nucleonic star are shown for comparison. mum mass of ∼ 2M ⊙ . However, we have ignored the possible presence of other hyperon species in the neutron star interior that could change this conclusion, although, we should point out that hypothetical repulsive NNY, NYY and YYY forces could lead to a similar one. Unfortunately, the lack of experimental information prevents currently any realistic attempt to estimate the effect of such forces. More experimental efforts are, therefore, needed. In particular, new informations about the presence of hyperons inside the core of neutron stars may be provided in the future through the observation, with the help of the new generation of gravitational wave detectors like the Einstein telescope [62,63,64], of signals emitted in the postmerger phase of binary neutron stars coalescence [65,66].