{'fraction_non_alphanumeric': 0.03655189620758483, 'fraction_numerical': 0.023484281437125748, 'mean_word_length': 4.7024719900409035, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Sensitivity of Neutron Star Observations to Three-nucleon Forces

Astrophysical observations of neutron stars have been widely used to infer the properties of the nuclear matter equation of state. Beside being a source of information on average properties of dense matter, however, the data provided by electromagnetic and gravitational wave (GW) facilities are reaching the accuracy needed to constrain, for the first time, nuclear dynamics in dense matter. In this work we assess the sensitivity of current and future neutron star observations to directly infer the strength of repulsive three-nucleon forces, which are key to determine the stiffness of the equation of state. Using a Bayesian approach we focus on the constraints that can be derived on three-body interactions from binary neutron star mergers observed by second and third-generation of gravitational wave interferometers. We consider both single and multiple observations. For current detectors at design sensitivity the analysis suggests that only low mass systems, with large signal-to-noise ratios (SNR), allow to reliably constrain the three-body forces. However, our results show that a single observation with a third-generation interferometer, such as the Einstein Telescope or Cosmic Explorer, will constrain the strength of the repulsive three-body potential with exquisite accuracy, turning third-generation GW detectors into new laboratories to study the nucleon dynamics.


I. INTRODUCTION
Lying at the interface between electromagnetic (EM) observatories, gravitational wave interferometers, and Earth based laboratories, multi-messenger astrophysics has the potential to shape a novel view of both structure and dynamics of dense nuclear matter. Mass-radius measurements of rotating pulsars are rapidly improving thanks to the information provided by the NASA satellite NICER [1][2][3][4][5][6], which has recently targeted the most massive neutron star (NS) known so far. Remarkably, NICER observations of PSR J0030+0451 and PSR J0740+662-the inferred masses of which are M = 1.34 +0. 16 −0.15 (M = 1.44 +0.15 −0.14 M ) and M = 2.072 +0.067 −0.066 M , respectively-yield comparable values of the stellar radius, pointing to a stiff nuclear matter equation of state (EOS) up to densities around four times nuclear density. On the other hand, constraints inferred from binary NS mergers detected by the LIGO/Virgo Collaboration, and in particular from the landmark discovery of GW170817, [7][8][9], have already ruled out some of the stiffest EOSs, which predict large tidal deformabilities, hinting instead to a softer matter content [10][11][12][13][14][15]. In addition, astrophysical data are being complemented by the information coming from terrestrial experiments, such as heavyion collisions or the recent measurement of the neutron skin thickness of lead, performed at Jefferson Lab by the PREX-II Collaboration [16][17][18][19][20][21][22][23].
Posterior distributions inferred from space-and ground-based facilities have been widely exploited in a variety of multi-messenger analyses, aimed at constraining models of the EOS or specific properties of neutron star matter. Examples of this approach include reconstruction of the EOS within both phenomenological and non-parametric frameworks, calculations based on micro-scopic models, and analyses focused on features such as the occurrence of phase transitions, or the behavior of the symmetry energy above nuclear density ; for recent reviews, see also Refs. [51,52] and references therein.
Recently, some of the authors of this article have proposed a novel approach, aimed at pushing the analyses based on multimessenger astrophysical information to a deeper level [53]. They argued that the accuracy of the currently available data-as well as that expected to be achieved by operating the existing detectors at design sensitivity-offer an unprecedented opportunity to constrain the microscopic models of nuclear dynamics at supranuclear density. The results reported in Ref. [53] show that the data set comprising the GW observation of the binary NS event GW170817, the spectroscopic observation of the millisecond pulsars PSR J0030+0451 performed by the NICER satellite, and the high-precision measurement of the radio pulsars timing of the binary PSR J0740+6620, providing information on the maximum NS mass, can, in fact, be exploited to infer quantitative insight on the strength of repulsive three-nucleon interactions in dense matter.
Unlike the nucleon-nucleon potential, the models of irreducible three-nucleon interactions are totally unconstrained beyond nuclear density. In most models, e.g. the Urbana IX potential employed to derive the EOS of Akmal, Pandharipande and Ravenhall (APR) [54], the strength of the isoscalar repulsive term-which plays a pivotal role in determining the stiffness of the nuclear matter EOS in the region relevant to neutron stars-is determined in such a way as to reproduce the empirical equilibrium density of isospin-symmetric matter [55,56]. In this context, the availability of additional information constraining the three-nucleon potential at larger density would be a major breakthrough.
The present work can be seen as a complementary follow up to the pioneering study of Ref. [53]. The analysis is first extended to consider a near-future scenario, using current interferometers at design sensitivity and stacking multiple binary NS observations characterised by different masses and distances. In addition, we apply, for the first time, the Bayesian approach to gauge the sensitivity of the Einstein Telescope (ET), a proposed third-generation ground-based GW observatory [57][58][59] The body of the article is structured as follows. In Sect. II we outline the dynamical model underlying our study, as well as the simple parametrisation adopted to characterise the strength of the repulsive component of the three-nucleon potential. The datasets considered in the analysis and the details of numerical simulations are described in Sections III A and III B, respectively, while the results are reported and discussed in Sect. IV. Finally, a summary of our findings and the prospects for future developments can be found in in Sect. V.

II. MODELLING NUCLEAR DYNAMICS BEYOND NUCLEAR DENSITY
The EOSs considered in our study have been derived using the formalism of non-relativistic nuclear manybody theory (NMBT). Within this framework, nuclear matter is pictured as a uniform system of point like nucleons, the dynamics of which is completely determined by the Hamiltonian 1 where m and p i denote the mass and momentum of the i-th nucleon, respectively. Interactions between matter constituents are driven by the nucleon-nucleon (NN) potential v ij -providing an accurate description of the two-nucleon system in both bound and scattering states-supplemented by the three-nucleon (NNN) potential V ijk , whose inclusion is needed to implicitly take into account the occurrence of processes involving the internal structure of the nucleon. As a consequence, the role of NNN interactions is expected to become more and more important with increasing density. Starting from Eq. (1), a number of different EOSs have been obtained using both different Hamiltonian models and different many-body techniques to calculate the ground state energy of nuclear matter as a function of baryon density. Purely phenomenological Hamiltonians, fitted to the properties of two-and three-nucleon systems, have been shown to provide a remarkably accurate account of the energies of the ground and low-lying excited states of nuclei with mass number A ≤ 12, as well as of their radii [60]. In addition, they allow to reproduce the empirical value of the equilibrium density of isospinsymmetric matter (SNM); see, e.g., Ref. [54] Over the past two decades, a great deal of attention has been given to a novel generation of nuclear Hamiltonians, derived using the formalism of Chiral Effective Field Theory (χEFT). Within χEFT, the nuclear potentials are obtained from effective Lagrangians comprising pion and nucleon degrees of freedom, constrained by the chiral symmetry of strong interactions. The main advantage of this approach is the capability to determine two-and many-nucleon potentials in a fully consistent fashion. However, being based on a low momentum expansion its applicability is inherently limited to densities 2 0 , with 0 = 0.16 fm −3 being the saturation density of SNM [61,62].
In this study, we have considered purely phenomenological Hamiltonians, which are expected to be best suited to describe the properties of nuclear matter in the density region extending up to ∼ 5 0 , relevant to NS applications. The reference line of our analysis is the Hamiltonian comprising the Argonne v 18 NN potential [63] (AV18) and the Urbana IX NNN potential [55,56] (UIX), which has been employed to obtain the APR EOS [54,64].
The AV18 potential is written as a sum of eighteen terms, needed to describe the complex operator structure of nuclear forces. It provides an accurate fit of the NN scattering phase-shifts for laboratory-frame energies up to ∼ 600 MeV, a value typical of NN collisions in strongly degenerate matter at density ∼ 4 0 [61]. A comparison with the central densities obtained from the solution of the Tolman-Oppenheimer-Volkoff equations [65,66] with the APR EOS [67] suggests that this phenomenological potential is adequate to describe NSs having masses as large as ∼ 2.1 M .
The UIX model of the NNN interaction is written as the sum of an attractive potential first derived by Fujita and Miyazawa [68]-describing two-pion exchange NNN processes with excitation of a ∆-resonance in the intermediate state-and a phenomenological repulsive potential; the resulting expression is The strength of the two-pion exchange contribution is adjusted to reproduce the observed ground state energies of 3 H and 4 He, obtained from accurate Monte Carlo calculations [55], whereas that of the isoscalar repulsive term is fixed to obtain the empirical saturation density of SNM-inferred from nuclear data-from variational calculations carried out using advanced many-body techniques [56]. It should be kept in mind that the repulsive term V R ijk implicitly takes into account relativistic corrections to the phenomenological two-nucleon potential v ij , which is determined by fitting NN scattering data in the centerof-mass reference frame. In the presence of the nuclear medium, however, the center of mass of the interacting nucleon pair is not at rest, and v ij must be boosted to take into account its motion [69]. The authors of Ref. [54] have modified the free-space AV18 potential to include the boost correction δv, whose effect is an enhancement of the repulsive contribution to the potential energy. As a consequence, using the boosted AV18 potential in calculations of nuclear matter energy entails the introduction of a modified NNN potential, referred to as UIX * , which turns out to be considerably softer than the UIX. The impact of relativistic corrections to the nuclear Hamiltonian on the description of NS properties has been recently discussed in Ref. [67].
The potentials describing NNN interactions are only determined by nuclear phenomenology reflecting nucleon interactions at SNM saturation density. On the other hand, they are totally unconstrained in the high-density regime relevant to NSs, in which their contribution is known to become dominant.
Motivated by the above consideration, in this work we extend the study of Ref. [53], whose authors have explored the possibility of inferring the strength of the repulsive term of the UIX * potential from data collected by multimessenger astrophysical observations, which carry information on nuclear dynamics at supranuclear denisity. Note that to pin down the dynamics of NNN interactions it is essential that the analysis be carried out using the the boost corrected NN potential.
Our study is based on the use of a set of Hamiltonians, obtained from the AV18 + δv + UIX * model performing the replacement The energy-density of nuclear matter at arbitrary baryon density and proton fraction x p has been obtained generalising the parametrisation employed in Ref. [54], that can be written in the form The explicit expressions of the functions appearing in Eqs. (4) and (5) can be found in the Appendix. They involve a set of parameters which were determined by fitting the energy per nucleon of SNM and pure neutron matter (PNM) computed within the FHNC/SOC variational approach [70] using the AV18+ δv + UIX * Hamiltonian.
The first two terms of Eq. (4) correspond to the proton and neutron kinetic energy, respectively, whereas the function g(ρ, x p ) describes the contribution arising from interactions. The assumption of quadratic dependence of the interaction energy on the neutron excess δ = 1 − 2x p is routinely employed in the literature to obtain the EOS of β-stable matter from those of SNM and PNM, and has been shown to be remarkably accurate over a broad range of values of the proton fraction x p ; see, e.g. Ref. [71].
Implementing the substitution of Eq. (3) is equivalent to adding a term (α − 1)V R at first order in perturbation theory. The corresponding change of energy density turns out to be with The functions δg can be readily expressed in terms of expectation values of V R in the nuclear matter ground state using Tabulated values of V R ijk as a function of density can be found in Ref. [54]. In our analysis, we have employed a polynomial fit including powers up to 3 which turned out to be very accurate. The values of the parameters a i are reported in Table I. Using the analytic expression of the energy density of nuclear matter at arbitrary proton fraction, composition and energy density of β-stable matter can be easily determined, by minimising with respect to x p , with the additional constraints of conservation of baryon number and charge neutrality. Finally, the matter pressure P , derived from standard thermodynamic relations, is used to obtain the EOS P ( ).
It has to be kept in mind that changing the strength of V R ijk affects the value of the nuclear saturation density predicted by the AV18 + δv + UIX * Hamiltonian. For this reason, we have limited the acceptable range of α to the interval [0.7, 2.0]. Within this range, the departure from the empirical value of 0 turns out to be ∼ 15% at most, and the corresponding change of the energy per particle never exceeds 3%.
Moreover, because the contribution of the repulsive NNN potential becomes large at supranuclear densities, the modification of its strength α marginally affect the ground-state energy of atomic nuclei. Using the results reported in Ref. [72], obtained from accurate Quantum Monte Carlo calculations, we have found that changing α from 1 to 1.3 results in a change of 4% and 6% of the ground state energies of 4 He and 12 C, respectively. These discrepancies appear to be fully acceptable in the context of our exploratory study.

III. METHODS AND OBSERVATIONS
We consider a family of EOS for which the observables of a neutron star (mass, radius and tidal deformability) depend uniquely on the three-body coefficient α and on the central pressure p c : Figure 1 shows the stable stellar configurations in the mass-radius plane and the mass-tidal deformability plane. Given a set O i=1,...,n of observations, we in- c } 2 using a hierarchical Bayesian approach, where is the likelihood of the i-th event (see Sec. III A below) and θ i denotes the set of relevant NS observables -mass and radius for pulsars, symmetric mass ratio and effective tidal deformability for GW observations -evaluated at {α, p (i) c } via (11). We assume that the priors on α and on each central pressure in Eq. (12) are uncorrelated.
The posteriors in Eq. (12) are sampled using the emcee with stretch move [73]. For each observation we run 100 walkers of 10 6 samples with a thinning factor of 0.02. The final distribution for α is obtained by marginalizing over the central pressures p c . When presenting results, we quote the median alongside the bounds of the 90% symmetric posterior density intervals.
We sample the central pressures of each star uniformly in log-space between ln 10 p min c (α) 34.58, where p c is expressed in dyne/cm 2 , and ln 10 p max We also impose a causality constraint, requiring that the speed of sound c s = dp/d is subluminal at the center of each NS.

A. Astrophysical datasets
We consider three real datasets corresponding to (i) the binary coalescence GW170817, (ii) the millisecond pulsar PSR J0030+0451 and (iii) the heaviest NS observed so far PSR J0740+6620. Dataset (iii) provides and update w.r.t. [53], in which PSR J0740+6620 was included only through the measurement of its mass, while here we also include the radius. We briefly summarize here the basic properties of each dataset and the corresponding likelihood functions that enter Eq. (12).
(i) -GW170817 is the first binary neutron star system observed by LIGO and Virgo. Under a low spin prior, the LVC analysis constrained the source component masses (m 1 , m 2 ) between ∼ 1.16M and ∼ 1.6M . GW170817 provided the first evidence that GW signals from coalescing systems are sensitive to matter effects induced by the NS structure, yielding a measurement for the effective tidal parameter ofΛ = 300 +420 −230 within 90% of the highest posterior density interval, with Λ 1,2 = λ 1,2 /m 5 1,2 being the NS individual, dimensionless, tidal deformabilities [7].
We construct the likelihood L(O GW170817 |η,Λ) from the joint posterior P(M, η,Λ|O GW170817 ) forΛ, the chirp mass M = (m 1 m 2 ) 3/5 /(m 1 + m 2 ) 1/5 , and the symmetric mass ratio η = m 1 m 2 /(m 1 + m 2 ) 5 . The calculation can be simplified by the fact that the chirp mass in the source frame is measured with ∼ 0.1% precision, which allows to fix it to its median value M = 1.186 M and restrict to the conditional probability P(η,Λ|M , O GW170817 ). Moreover, as shown in [33], the latter can be replaced by the marginalized posterior P(η,Λ|O GW170817 ) to very good accuracy. This choice reduces the number of parameters to be sampled, since the central pressure p c } and α 3 , and similarly for the individual masses m 1,2 and tidal deformabilities Λ 1,2 . The likelihood function 4 is then obtained by re-weighting the posterior by the joint prior on η andΛ as derived from [7], Note that, although p (2) c is not independently sampled, we still require it to lie within its prior support.
(ii) -For the millisecond pulsar PSR J0030+0451 we use the joint mass-radius posterior P(M, R|O J0030 ) inferred by the NICER collaboration, which has carried out two independent studies of the stellar spectroscopic observations, obtaining consistent results. The massradius constraints provided by the two collaborations led to M = 1.34 +0. 15 −0.16 M and R = 12.71 +1.14 −1.19 km [3], and M = 1.44 +0. 15 −0.14 M and R = 13.02 +1.24 −1.06 km [4] respectively (68% credibility). Here we use the data publicly available in [74], for which the likelihood can be derived straightforwardly from P(M, R|O J0030 ) because the joint prior on {M, R} is flat, 3 More specifically, we compute m 2 from m 1 (α, p c ) and M and then we solve m 2 ≡ m 2 (α, p c ) for p (2) c . 4 Note that the likelihood we use here for GW170817 is different from the one of Ref. [53] in which a three-dimensional distribution L GW (q, Λ 1 , Λ 2 ) was considered, with q = m 1 /m 2 .
(iii) -PSR J0740+6620 [5,6] is the most massive pulsar discovered so far. Previous observations of this source constrained its mass to M = 2.08 +0.072 −0.069 M (68.3% credibility) [2]. This measurement, combined with data obtained from the XMM Newton European Photon Imaging Camera to improve the NICER background, was used in [5,75] and [6,76] to infer the pulsar radius, with the two teams obtaining R = 12.39 +1. 30 −0.98 km and R = 13.7 +2.62 −1.50 km [6] respectively (68% credibility). Here we use the data in [77], for which the likelihood can be immediately inferred from the posterior due to uniform priors,

B. Simulations for 2G and 3G detectors
We simulate 5 30 binary neutron star events for two choices of the three-body strength, α = 1 and α = 1.3, either for a network (HLV) composed by the LIGO Hanford, LIGO Livingston, and Virgo detectors at design sensitivity [79], or for the future third-generation interferometer Einstein Telescope in its ET-D configuration [58]. We inject 64-second long waveforms into a zero-noise configuration as described in [80], with sky location and inclination uniformily distributed over the sky. Posterior parameters are recovered using the bilby software [81,82] for GW injections and parameter estimation. For both injection and recovery, we model binary neutron star signals with the IMRPhenomPv2 NRTidal waveform template [83,84]. Injected binaries are nonspinning, while component spins are recovered imposing a low-spin prior χ 1,2 ∈ [−0.05, 0.05] and assuming that spins are (anti-) aligned.
We assume that tidal parameters are recovered uniformly w.r.t.Λ and the tidal parameter δΛ which contributes at higher post-Newtonian order in the waveform phase expansion [85], with the additional constraint that the individual deformabilities Λ 1,2 of the binary components lie between 0 and 5000.

IV. RESULTS
We start the discussion of our results by focusing first on the the Bayesian analysis applied to the three real observations described in the previous section.
The inferred probability distributions for α are summarized by the density plots in the left column of Fig. 3, together with their median values and 90% confidence intervals. The analyses for GW170817 and for J0030+0451 5 We limit our catalogue to 30 events because the recovery of the EOS is expected to be biased by a mismodelling of the underlying BNS population distribution if the number of sources exceeds ∼ 30 [78]. have been already presented in [53], while the novel massradius measurement obtained by NICER allows us to perform an independent study of the three-body strength for J0740+6620, and a direct comparison with other observations. Interestingly the posterior densities of Fig. (3) show very similar results for the two EM observations, with a nearly identical median around α 1.4. The probability distribution for J0740+6620 peaks around a slightly larger value compared to the lighter pulsar, J0030+0451, since larger values of α tend to support more massive configurations. Moreover, even if P(α) shows support for the baseline model α = 1, which lies within the 90% CL of the distributions, EM observations seem to consistently favour larger values of the 3body amplitude, reflecting stronger repulsive NNN interactions. As observed in [53], the distribution of α inferred by GW data alone is unconstrained, with the posterior rallying against the lower prior at α = 0.7, while the multi-messenger analysis is dominated by the pulsar measurements, and in particular by J0740+6620, leading to values of α 1. Constraints on α, i.e on the microscopic Hamiltonian (1), can be translated into bounds on the stellar macroscopic observables. The right column of Fig. (3) shows, for example, the maximum mass density distributions predicted by the values of α inferred for each dataset. All the observations lead to median values of M max 2.2M , with the multi-messenger analysis yielding a probability distribution with large support for M max ∼ 2.5M . In Fig. 4 we also show the M -R density distribution corresponding to the 90% CL of α for the multimessenger case. Light (dark) colors identify stellar profiles with high (low) probability. Pulsar observations drive the profiles far from the α = 1 baseline, i.e. towards stiffer NS configurations, with an expected radius R 12 km for a prototype NS with M = 1.4M .
So far our analysis shows that, although the constraining power of current measurements is still limited, astrophysical data are already sensitive to nucleon dynamics. We will therefore explore the insights that can be inferred on three-body nuclear forces exploiting future GW observations of binary inspirals.
As discussed in Sec. III B we have simulated two catalogues of 30 binary NS mergers, observed either by 2G network or by ET, assuming two different values of the three-nucleon strength. Source parameters, i.e. masses and tidal deformabilities, are first recovered with Bilby, and then analyzed by our Bayesian pipeline which samples the posterior distribution of α. Figure 5 shows the posterior densities P(α) of each event, for injected NSs with α = 1, detected by the HLV network. The ability of 2G detectors to discriminate the actual value of the three-body strength substantially depends on both the SNR and on the component masses of the binary. We find that observations with SNR smaller than ∼ 25 lead α to be almost unconstrained, with the true value always lying outside the 90% confidence interval of the distribution. However, even for strong signals, accurate measurements only occur for low-mass systems with a chirp mass M 1.4M . This is particular evident for the event with the largest SNR (∼ 35) in our set. Such binary features two heavy NSs with a chirp mass M 1.6M , and provides loose bounds on α. Moreover, Fig. 5 shows that, with the exception of four events with SNR> 30 and M < 1.4M , the remaining posteriors always prefer large values of the three-nucleon strength, at the edge of the upper prior boundary. This particular behavior reflects a systematic bias we find in the posteriors ofΛ inferred by GW observations for binaries with heavy components, which tend to favour large values of the tidal parameter. Its effect on the marginal distribution of α becomes even more pronounced in the high mass scenario where the tidal deformability becomes less sensitive to variations of α. We believe such bias may be induced by our choice of priors on the tidal parameters, which has strong support against the BBH hypothesis Λ = 0, and reflects the physical assumption that compact objects with m 1,2 3M are neutron stars. Moreover, the stack of multiple GW signals only partially alleviate the bias in favour of large three-body strength. We have indeed combined different observations with SNR larger than 20, finding a mild improvement of the posterior support towards the true value of α. The results discussed so far hold qualitatively also when we consider binary NSs simulated with α = 1.3.
This picture changes dramatically when signals are observed by the Einstein Telescope. Figure 6 shows indeed the distributions of the three-nucleon strength inferred by the 3G detector, for both families of events simulated with α = 1 and α = 1.3. The exquisite sensitivity of ET allows to gauge away the bias arising from the 2G network. All the posteriors peak around the injected values of α, showing no support on the prior boundaries. In the best (worse) case scenario we find that α can be constrained with ∼ 2% (∼ 30%) of accuracy at 68% confidence level. Such accuracy allows to disentangle the two values of the three-body strength we consider. Even in the most pessimistic cases, where the inferred P(α) are not narrow enough to identify a specific value of α, stacking of few events would render the distributions clearly distinguishable. Figure 7 shows the posteriors obtained by combining six events of our catalogue 6 leading to loose constraints on α. The final posteriors for α = 1 and α = 1.3 are clearly separated, with a negligible overlap on the tails.
Such accuracy translates into very narrow constraints on the mass-radius (or equivalently mass-tidal deformability) diagram. As an example, we show in Fig. (8) the M -R profile density computed from the values of α inferred from event number 17 of our dataset. A direct comparison with Fig. 4, where a similar plot was made for data from current facilities, provides a clear hint on the possibility to use ET as a new laboratory to study the dynamics of nucleon interactions in the stellar cores.

V. CONCLUSIONS
We have investigated the sensitivity of NS observations to the strength of repulsive three-nucleon forces, which are known to be critical in determining the stiffness of the nuclear matter EOS at supranuclear densities. Our analysis is based on the AV18 + δv + UIX * nuclear Hamiltonian and involves a single free parameter, to be constrained by data, determining the coupling constant appearing in the repulsive contribution to the UIX * potential.
We have performed hierarchical bayesian inference employing the current available multimessenger datasets in order to constrain this parameter. We have then repeated the analysis with a set of simulated GW observations that could be performed by both current (LIGO/Virgo) and future (Einstein Telescope) interferometers at design sensitivity. This analysis has the main purpose to explore the potential of near and next generation facilities into inferring crucial information about the microscopic dynamics of nuclear matter.
The analysis with real data has been carried out employing some of the dataset used in a previous work [53]. Our results suggest that even if current facilities show a clear sensitivity to small variation of the NNN repulsive potential, they are not accurate enough to capture significant insights. This picture is cross-validated by the population analysis performed with mocked LIGO/Virgo data, with binaries generated with two different values of the three-body strength, α = 1 and α = 1.3. Only few, low-mass and high SNR events provide a meaningful constraint on α, with posterior distributions correctly peaked around the injected values. Moreover, even for the most constraining event, the inferred posteriors do not allow a clear disentanglement between the two values of α we considered. The picture improves only slightly with the stacking of multiple observations. These results exhibit a striking upgrade when we assume that the population of binaries is observed by the Einstein Telescope. In most of the cases, the large SNRs obtained by such events in combination with the 3G detector allow the posteriors for the injected values of α to be clearly separated, and only a single observation is needed to resolve them.
Moreover, in the few cases where posteriors overlap, stacking of ∼ 2 − 3 observations would allow to unambiguously distinguish between α = 1 and α = 1.3. The same conclusion would apply assuming that binaries are detected by the proposed Cosmic Explorer [86,87]. The large SNRs expected in the 3G era also require a careful assessment of waveform systematics which could bias the parameter reconstruction [85,[88][89][90]. However, our results strongly support the evidence that with the upcoming third generation detectors, our understanding of neutron star matter will make a great step forward into the direction of using NS observations to probe fundamental physics at the fermi scale.
Further applications of our approach can be pursued following multiple directions, and in particular considering how constraints on nucleon dynamics would improve by joint analyses of the inspiral and of the postmerger phase, exploiting for the latter either GW oscillation modes [91][92][93], or electromagnetic counterparts emitted by the binary remnant [94].  In this Appendix, we report the explicit expression of the energy density of nuclear matter employed to carry out our analysis. This expression was originally derived from a fit to the EOSs of SNM and PNM obtained by Akmal et al. [54] using the AV18 + δv + UIX * nuclear Hamiltonian and the variational FHNC/SOC formalism.
The energy density of nuclear matter at baryon density and proton fraction x p is written according to Eqs. (4) and (5) ( , x p ) = The density¯ 2 0 corresponds to the onset of the high-density phase-featuring spin-isospin density waves associated with neutral pion condensation-predicted by the study of Ref. [54].
The values of the parameters appearing in the above equations are given in Table II