### Nuclear EOSs from chiral EFT

The EOS set used in this work is constrained at low densities by microscopic calculations of neutron matter using interactions from chiral EFT. In these microscopic calculations, the Schrödinger equation for the many-body system is solved numerically. This requires a nuclear Hamiltonian and a method to solve the Schrödinger equation with controlled approximations.

To obtain the Hamiltonian describing the dense matter EOS studied in this work, we use chiral EFT. Chiral EFT is a low-energy effective theory of QCD, and describes strong interactions in terms of nucleon and pion degrees of freedom instead of quarks and gluons^{30,31}. To construct the interactions, the most general Lagrangian in terms of nucleons and pions, consistent with all symmetries of QCD, is expanded in powers of momenta. Using a power counting scheme, the individual contributions are arranged according to their importance. By going to higher orders, the description of interactions becomes more precise, but the individual contributions become more involved. The chiral EFT Lagrangian explicitly includes pion-exchange interactions among nucleons whereas all high-energy details that are not explicitly resolved are expanded in terms of general contact interactions. These are accompanied by low-energy couplings, which are fitted to experimental data.

Chiral EFT interactions have several benefits over phenomenological interaction models: they naturally include many-body forces consistent with two-nucleon interactions, they can be systematically improved, and they enable theoretical uncertainty estimates^{16,32}. The last of these can be extracted from order-by-order calculations and are important when analysing astrophysical observations for which interactions are extrapolated to conditions that cannot be recreated in experiments at present.

In this work, we constrain our EOSs with theoretical calculations at zero temperature using local chiral EFT interactions^{14,53,54,55}. We use quantum Monte Carlo methods^{33}, in particular the auxiliary-field diffusion Monte Carlo method, which are among the most precise many-body methods to solve the nuclear many-body problem. The results of these calculations provide constraints on the EOS up to densities of around 2*n*_{sat} (ref. ^{29}).

The region of applicability of the chiral EFT expansion is determined by the so-called breakdown scale, which is estimated to be of the order of 500–600 MeV/*c* (ref. ^{56}). Hence, the chiral EFT expansion breaks down at densities (gtrsim 2{n}_{{rm{sat}}}), indicated by increasing uncertainty estimates between 1 and 2*n*_{sat}. At these densities, high-energy physics that is encoded in short-range contact interactions needs to be explicitly taken into account. Therefore, chiral EFT cannot be used to constrain the EOS at higher densities as probed in the cores of neutron stars. To extend the EOS to these densities, we use a general extrapolation scheme in terms of the speed of sound^{35} (see also ref. ^{57}).

To construct the neutron-star EOS set, we first extend our chiral EFT calculation to *β*-equilibrium and add a crust^{58}. We use microscopic input up to 1.5*n*_{sat} to constrain the EOS, but a variation in the range 1–2*n*_{sat} shows no substantial impact on our final results for neutron-star radii^{34}. Above this density, we sample a set of six randomly distributed points in the speed of sound plane at baryon densities between 1.5 and 12*n*_{sat}, enforcing 0 ≤ *c*_{s} ≤ *c* at each point. A variation of the number of sampled points between 5 and 10 does not impact our findings. We then connect these points by line segments, reconstruct the EOS and solve the Tolman–Oppenheimer–Volkoff equations to extract neutron-star properties. For each EOS, we also construct a partner EOS that includes a segment with vanishing speed of sound to explicitly simulate strong first-order phase transitions. We sample the onset density and width of this segment randomly.

Our EOS set includes 15,000 different EOS samples for which the prior on the radii of neutron stars is naturally determined by the EOS expansion scheme. We have explicitly checked the differences among a prior uniform in the radius of a typical 1.4({M}_{odot }) neutron star and the ‘natural’ prior and found only minor changes once astrophysical and HIC data are included (Extended Data Table 1).

Recently, first results for the EOS of symmetric nuclear matter between 3 and 10*n*_{sat} from functional renormalization group calculations that are based on QCD became available^{59}. This offers a very promising future tool to constrain dense neutron-star matter when calculations for asymmetric matter will become available.

### Multi-messenger analysis of astrophysical data

To constrain the set of EOSs derived from chiral EFT with astrophysical data, we use a multi-step procedure in which results from individual steps are used as a prior for the next part of the analysis^{9} (Extended Data Fig. 3). First, we incorporate constraints on the maximum mass of neutron stars. For this, we implement the mass measurements of the heavy radio pulsars PSR J0348+0432 (ref. ^{36}) and PSR J1614-2230 (ref. ^{37}). As we make use of the NICER and XMM mass–radius information of PSR J0740+6620 (refs. ^{7,8}) at a later stage, we do not include the mass measurement of PSR J0740+6620 (ref. ^{38}) to avoid double counting. The combination of these observations^{9,60} of high-mass neutron stars provides a lower bound on the maximum mass of neutron stars. By contrast, an upper bound of the maximum mass is obtained from the observation of the merger remnant of the neutron-star merger GW170817 (ref. ^{41}). Among other arguments, the observation of a bright, red kilonova component and the observation of a short gamma-ray burst 2 s after the merger of the two neutron stars indicate that the remnant experienced a delayed (({mathscr{O}}({rm{100ms}}))) collapse to a black hole, so that an upper limit on the maximum mass can be derived. The combined estimate of the maximum mass, ({2.21}_{-0.13}^{+0.10}{M}_{odot }) at 68% uncertainty, already provides important information about the internal structure of neutron stars and disfavours both too stiff and too soft EOSs (that is, EOSs with too large and too small pressures, respectively).

In the next step, we incorporate NICER’s mass and radius measurement of PSR J0030+0451 (ref. ^{5}) and PSR J0740+6620 (refs. ^{7,8}). NICER, located on board of the International Space Station, is a NASA telescope measuring the X-ray pulse profile of pulsars. By correlation of the observed profile and brightness with theoretical predictions, it is possible to extract information on the configuration (for example, on the location and properties of hotspots on the neutron-star surface, the rotation rate of the star, and its compactness, which determines the light bending around the pulsar). This information enables constraints on the pulsar’s mass and radius. In addition to NICER, the XMM-Newton telescope^{61,62} has been used for the analysis of PSR J0740+6620 (ref. ^{7}) to improve the total flux measurement. For PSR J0740+6620, we average over the results obtained in refs. ^{7,8}, whereas for PSR J0030+0451 we use only results of ref. ^{5}.

Next, we analyse the GW signal emitted from the binary neutron-star merger GW170817 (ref. ^{1}), as well as its observed kilonova AT2017gfo (ref. ^{3}). Finally, we also incorporate the second confirmed GW signal from a binary neutron-star merger GW190425 (ref. ^{2}). For GW170817 and GW190425, we assumed both of them to be emitted by binary neutron star mergers. To test the robustness of the GW analysis, we have explored a number of different GW models and found only a minimal impact on the final EOS constraint^{9}. Results shown in the main text are obtained using the parallel bilby software^{63} and the waveform model IMRPhenomPv2_NRTidalv2 (ref. ^{42}) for cross-correlation with the observed data^{1}. IMRPhenomPv2_NRTidalv2 is an updated model of the waveform model used in previous analyses by the Laser Interferometer Gravitational-Wave Observatory (LIGO)/Virgo Collaboration^{2,43} and, hence, allows for a more accurate measurement of tidal effects. The likelihood function for the GW analysis ({ {mathcal L} }_{{rm{GW}}}) is given by^{64}

$${ {mathcal L} }_{{rm{GW}}}propto exp left(-2int {rm{d}}ffrac{{|mathop{d}limits^{ sim }(f)-mathop{h}limits^{ sim }(f)|}^{2}}{{S}_{{rm{n}}}(f)}right),$$

(1)

in which (mathop{d}limits^{ sim }(,f)), (mathop{h}limits^{ sim }(,f)) and *S*_{n}(*f*) are the observed data, the waveform template and the power spectral density of the noise, respectively. To ensure full coverage of the binary neutron stars’ inspiral signal, we have analysed the data up to 2,048 Hz. To avoid the low-frequency noise wall in the detectors, a low-frequency bound of 20 Hz is used.

Similarly, we use Bayesian inference to analyse the observed kilonova AT2017gfo. The likelihood function for the light curve analysis ({ {mathcal L} }_{{rm{EM}}}) is given by^{65}

$${ {mathcal L} }_{{rm{EM}}}propto {chi }_{1}^{2}left(sum _{ij}frac{1}{{n}_{j}-1}{left(frac{{m}_{i}^{j}-{m}_{i}^{j,{rm{est}}}}{{sigma }_{i}^{j}}right)}^{2}right),$$

(2)

in which ({m}_{i}^{j,{rm{est}}}) are the estimated or theoretically predicted apparent magnitudes for a given filter *j* (a passband for a particular wavelength interval) at observation time *t*_{i} with *n*_{j} data points for filter *j*. Moreover, ({m}_{i}^{j}) and ({sigma }_{i}^{j}) are the observed apparent magnitude and its corresponding statistical uncertainties, respectively. For this analysis, the probability distribution of a chi-squared distribution with a degree of freedom of 1, ({chi }_{1}^{2}), is taken as the likelihood measurement. To reduce the systematic error of the kilonova modelling below the statistical error, a further uncertainty of 1 mag is added to the measurement error. To analyse AT2017gfo, we use the radiative transfer code POSSIS^{44} to produce grids of light curves for multidimensional kilonova models with the following free parameters: the dynamical ejecta mass, the disk wind ejecta mass, the opening angle of the lanthanide-rich dynamical-ejecta component, and the viewing angle. To enable inference, we combine the grid with a framework combining Gaussian process regression and singular value decomposition^{66} to compute generic light curves for these parameters. To connect the ejecta parameters, which determine the exact properties of the light curve, with the binary neutron-star system parameters, we assume that the total ejecta mass is a sum of two components: dynamical ejecta, released during the merger process through torque and shocks, and disk-wind ejecta. Both components, the dynamical ejecta^{66} and the disk-wind ejecta^{9}, are correlated to source parameters of the binary neutron-star system based on numerical relativity simulations^{9,66,67}.

### Constraining the symmetric nuclear matter EOS at high density with HICs

Over the last two decades, major experimental efforts have been devoted to measuring the nuclear EOS with HIC experiments performed at relativistic incident energies^{27,68,69}. These collisions of atomic nuclei form a hot, dense fireball of hadronic matter in the overlapping region, which expands in time and reaches the surrounding detectors as baryons and mesons. The phase-space distribution of particles flowing from the fireball during the expansion phase is strongly dictated by the compression achieved in the colliding region and is, therefore, sensitive to the EOS of the hot and dense nuclear matter created in the collision. Important progress has been made recently in modelling intermediate-energy HICs, but theoretical uncertainties still remain^{70,71}. In the present analysis, results obtained with different models are found to be compatible within their quoted errors.

The so-called elliptic flow (*v*_{2}) of emerging particles is the main observable, which has been used to experimentally constrain symmetric nuclear matter at supranuclear densities with HICs. It is described by the second moment of the Fourier expansion of the distribution of the azimuthal angle *Φ* of the emitted particles with respect to that of the reaction plane *Φ*_{RP}

$$begin{array}{c}frac{{rm{d}}sigma (y,{p}_{t})}{{rm{d}}Phi }=C(1+2{v}_{1}(y,{p}_{{rm{t}}})cos (Phi -{Phi }_{{rm{RP}}})\ +2{v}_{2}(y,{p}_{{rm{t}}})cos ,2(Phi -{Phi }_{{rm{RP}}})+mathrm{..}.),,end{array}$$

(3)

in which all expansion coefficients *v*_{n} are functions of longitudinal rapidity (y=frac{1}{2},mathrm{ln}left(frac{E+{p}_{z}}{E-{p}_{z}}right)), with *p*_{z} being the momentum along the beam axis and *E* being the total energy, and of transverse momentum ({p}_{{rm{t}}}=sqrt{{p}_{x}^{2}+{p}_{y}^{2}}) of the particle, with *p*_{x} and *p*_{y} denoting the momentum components perpendicular to the beam axis.

In the experiment, the orientation of the reaction plane is event-wise reconstructed from the azimuthal distribution of particles recorded in the forward and backward hemispheres, and the Fourier coefficients are corrected for the finite resolution of this procedure^{72}. The coincident particle and fragment emissions are also used for the reconstruction of the impact parameter of each reaction event^{11}. A positive elliptic flow *v*_{2} indicates a preferred emission in the reaction plane whereas a negative flow indicates an emission out of the reaction plane.

It has been shown that the elliptic flow *v*_{2} of protons emitted at rapidities intermediate between projectile and target rapidity (mid-rapidity) in HICs at incident energies of several hundred MeV per nucleon offers the strongest sensitivity to the nuclear EOS^{10,27,73}, as evident from calculations made with various transport models. This dependence on the nuclear EOS is predicted by quantum molecular dynamics (QMD)^{10,73,74,75} and Boltzmann–Uehling–Uhlenbeck^{27} models. The origin of the phenomenon has been investigated in detail elsewhere^{76}. As shown in ref. ^{27}, at higher beam energies between 1 and 10 GeV per nucleon, the sensitivity of the directed flow *v*_{1} to the stiffness of the EOS of symmetric nuclear matter becomes comparable to that of *v*_{2}. Overall, from HICs performed at incident beam energies of a few hundred MeV per nucleon up to around 10 GeV per nucleon, the flow data indicate an EOS for symmetric nuclear matter with an incompressibility *K* below 260 MeV. Using FOPI data on the elliptic flow in gold–gold collisions between 400 MeV and 1.5 GeV per nucleon, thanks to the broad acceptance of the detector, an enhanced precision in the determination of the EOS could be achieved. Including the full rapidity and transverse momentum dependence of the elliptic flow of protons and heavier isotopes^{10} in the analysis with the Isospin-QMD (IQMD) transport model, the incompressibility was determined as *K* = 190 ± 30 MeV. This result was confirmed by interpreting the same data with three Skyrme energy-density functionals introduced into the ultrarelativistic QMD (UrQMD) transport model^{75}, leading to *K* = 220 ± 40 MeV. The interval of confidence used in the present study, *K* = 200 ± 25 MeV, reflects both predictions. The densities probed were estimated to range between 1 and 3*n*_{sat} by analysing the densities effective in building the elliptic flow in IQMD simulations^{10}. Note that the constraints deduced from the analysis of elliptic flow are compatible with earlier findings of the Kaon Spectrometer Collaboration obtained from comparisons of QMD predictions with experimental *K*^{+} meson production yields from gold–gold and carbon–carbon collisions performed at GSI between 0.6 and 1.5 GeV per nucleon^{77,78}.

### The ASY-EOS experiment to measure the symmetry energy

Nuclear experiments can help to constrain the EOS of neutron matter (see, for example, the PREX experiment measuring the neutron-skin thickness in lead nuclei^{79,80,81,82}). It has been suggested^{83} that flows of particles in HICs can be used to constrain the EOS of neutron matter through the symmetry energy at supra-saturation density. However, nuclear matter that can be studied in laboratory experiments using HICs is not very neutron rich. Therefore, the effect of the symmetry energy on *v*_{2} remains small, close to or below the uncertainties of the main contribution of the symmetric nuclear matter EOS. To enhance observable effects related to the symmetry energy, the use of the elliptic flow ratio of particles with large isospin difference, ideally the ratio for neutrons over protons ({v}_{2}^{{rm{np}}}={v}_{2}^{{rm{n}}}/{v}_{2}^{{rm{p}}}), was proposed^{84}. This method has been adopted for the ASY-EOS experiment performed at GSI in Darmstadt, studying collisions of gold nuclei of 400 MeV per nucleon incident energy and gold targets. The description of the experiment and the analysis with the UrQMD transport model are given in detail elsewhere^{11}. ASY-EOS benefited from the Large-Area Neutron Detector (LAND)^{85} permitting the detection of neutrons and charged particles within the same acceptance. LAND was placed to cover mid-rapidity emissions over a large *p*_{t} interval. Its isotopic resolution in this experiment was not sufficient to uniquely identify protons. Elliptic flow ratios as a function of *p*_{t} were, therefore, determined for neutrons with respect to all charged particles within the LAND acceptance. We note that for the selected collisions (central up to semi-central) and angular region, the yield of charged particles consists of light isotopes, mainly protons (around 50%) according to FOPI data for the same reaction. Confronted with UrQMD transport model predictions (and confirmed with other models, IQMD^{74} and Tübingen QMD^{86}), the resulting flow ratio enabled deduction of a constraint for the symmetry energy, which is so far the most precise for supra-saturation densities obtained from HICs. The results are detailed in the following section. As indicated by QMD model predictions, densities probed by the elliptic flow ratio in the ASY-EOS experiment extend up to about 2*n*_{sat}.

### Implementation of nuclear EOS constraints from HICs

For analysing the experimental elliptic flow data, an EOS functional needs to be fed into the QMD simulations for both symmetric and asymmetric nuclear matter. This is given by the parameterization for the energy per particle

$$frac{E}{A}(n,delta )approx frac{E}{A}(n,0)+S(n){delta }^{2},,$$

(4)

with the baryon density *n* = *n*_{n} + *n*_{p} and the isospin asymmetry *δ* = (*n*_{n}*−n*_{p})/*n* = 1 − 2*x*, in which *n*_{n} and *n*_{p} are the neutron and proton densities, respectively, and *x* = *n*_{p}/*n* is the proton fraction. *E*/*A*(*n*, 0) is the energy of symmetric nuclear matter, *E*/*A*(*n*, 1) corresponds to pure neutron matter, and *S*(*n*) is the symmetry energy defined here as the difference between the two. For the analysis of the FOPI experiment, the first term in equation (4) has been parameterized with

$$frac{E}{A}(n,0)=frac{3}{5}{left(frac{n}{{n}_{{rm{sat}}}}right)}^{2/3}{E}_{{rm{F}}}+frac{alpha n}{2{n}_{{rm{sat}}}}+frac{beta }{gamma +1}{left(frac{n}{{n}_{{rm{sat}}}}right)}^{gamma },$$

(5)

with the saturation density *n*_{sat}, the Fermi energy *E*_{F}, and in which the parameters *α*, *β* and *γ* are fixed by the incompressibility *K*, the binding energy *B* of symmetric nuclear matter at *n*_{sat}, and the condition that the pressure of symmetric nuclear matter is zero at saturation density:

$$begin{array}{l}alpha =-2left(frac{K+frac{6{E}_{{rm{F}}}}{5}}{9(gamma -1)}+frac{2}{5}{E}_{{rm{F}}}right),\ beta =left(K+frac{6}{5}{E}_{{rm{F}}}right)frac{gamma +1}{9gamma (gamma -1)},,,gamma =frac{K+frac{6{E}_{{rm{F}}}}{5}}{9left(frac{{E}_{{rm{F}}}}{5}+Bright)}.end{array}$$

(6)

In the ASY-EOS analysis, the *S*(*n*) term of equation (4) has been parameterized as

$$S(n)={E}_{{rm{kin}},0}{left(frac{n}{{n}_{{rm{sat}}}}right)}^{2/3}+{E}_{{rm{pot}},0}{left(frac{n}{{n}_{{rm{sat}}}}right)}^{{gamma }_{{rm{asy}}}}.$$

(7)

At saturation density, the kinetic part has been set to *E*_{kin,0} = 12 MeV and *E*_{pot,0} = S_{0} − *E*_{kin,0}. The parameter *γ*_{asy} was extracted from fits to experimental data of the *p*_{t} dependence of the elliptic flow ratio of neutrons over charged particles around mid-rapidity. In particular, this results in *γ*_{asy} = 0.68 ± 0.19 for *S*_{0} = 31 MeV and *γ*_{asy} = 0.72 ± 0.19 for *S*_{0} = 34 MeV (see Extended Data Fig. 4 for a comparison with microscopic neutron matter calculations). Here we interpolate *γ*_{asy} assuming a linear function with *S*_{0}, for which the uncertainty is chosen to be 0.19 independent of *S*_{0}. We have studied the behaviour of *γ*_{asy} as a function of *S*_{0} for two different QMD models (Extended Data Fig. 1) and confirmed that the linear interpolation in the *S*_{0} range is suitable.

The pressure constraint is given by the density derivative of the energy per particle of equation (4)

$$P(n,delta )={n}^{2}frac{partial E/A(n,delta )}{partial n},$$

(8)

and depends on *n*, *δ*, *n*_{sat}, *B*, *K* and *S*_{0}. We enforce this constraint only at densities for which the experiment is sensitive. The density region of the HIC constraint is set by the sensitivity of the neutrons-over-charged-particles flow ratio determined for the ASY-EOS experiment^{11} (see also the previous section). This sensitivity curve covers the density range from 0.5*n*_{sat} up to 3*n*_{sat} and peaks between *n*_{sat} and about 2*n*_{sat}, for which the experiment is most sensitive.

Neutron-star matter is composed of neutrons, protons, electrons and muons in *β*-equilibrium. To apply the ASY-EOS constraint to neutron stars, we need to determine the proton fraction *x*_{ASY-EOS} accordingly. For simplicity, we neglect muons because they have only a small impact on the neutron-star EOS in the considered density range. Then, the density of electrons is equal to the proton density owing to local charge neutrality, and the proton fraction *x* at a given baryon density *n* is fixed by the *β*-equilibrium condition

$${mu }_{{rm{n}}}(n,x)={mu }_{{rm{p}}}(n,x)+{mu }_{{rm{e}}}({n}_{{rm{e}}}=xn),$$

(9)

in which *μ*_{n,p,e} is the chemical potential of the respective particle species. We calculate the neutron and proton chemical potentials consistently with equations (4)–(7). Electrons are modelled as an ultrarelativistic degenerate Fermi gas with pressure *P*_{e} = *E*_{e}/(3*V*), energy density ({E}_{{rm{e}}}/V=hbar c{(3{{rm{pi }}}^{2}{n}_{{rm{e}}})}^{4/3}/(4{{rm{pi }}}^{2})) and chemical potential ({mu }_{{rm{e}}}=hbar c{(3{{rm{pi }}}^{2}{n}_{{rm{e}}})}^{1/3}), in which (hbar ) is the reduced Planck constant and V the volume.

The final pressure constraint is obtained using *E*_{F} = 37 MeV and by varying the parameters *n*_{sat}, *B*, *K* and *S*_{0} within specific ranges. For the parameters describing symmetric nuclear matter, we use the values consistent with the FOPI analysis given by *n*_{sat} = 0.16 fm^{−3}, *B* = 16 MeV, and a Gaussian distribution for *K* with *K* = 200 ± 25 MeV at 1*σ*. Regarding *S*_{0}, we apply a uniform prior in the range 31–34 MeV. We further use results for the pressure of symmetric nuclear matter deduced elsewhere^{27} and disregard all parameter sets, which lead to a pressure that is not consistent with their constraint in the overlapping density range for which ASY-EOS remains sensitive, between 2 and 3*n*_{sat} (Extended Data Fig. 5). We note that the value of *K* has very little influence on the observables measured by ASY-EOS to extract the symmetry energy^{86}. We have explicitly checked the robustness of our results when using larger uncertainty ranges for all nuclear matter parameters in agreement with theoretical predictions^{17} and found their influence on our final result to be negligible (Extended Data Table 3). In particular, we have used a larger range for *S*_{0} between 30 and 35 MeV and the following Gaussian distributions for *n*_{sat}, *B* and *K*: *n*_{sat} = 0.164 ± 0.007 fm^{−3}, *B* = 15.86 ± 0.57 MeV and *K* = 215 ± 40 MeV at the 1*σ* level.

### Combination of the astronomical and HIC constraints

To combine the experimental and observational EOS constraints, we use Bayesian inference. The EOS posterior is given by

$$begin{array}{ll}p({rm{EOS}}|{rm{MMA}},{rm{HIC}}) & propto ,p({rm{HIC}}|{rm{EOS}})\ & times ,p({rm{MMA}}|{rm{EOS}})p({rm{EOS}})\ & =,p({rm{HIC}}|{rm{EOS}})p({rm{EOS}}|{rm{MMA}})\ & equiv ,{ {mathcal L} }_{{rm{HIC}}}({rm{EOS}}){{mathscr{P}}}_{{rm{MMA}}}({rm{EOS}}),end{array}$$

(10)

in which MMA denotes multi-messenger astrophysics, ({ {mathcal L} }_{{rm{HIC}}}({rm{EOS}})) is the likelihood of the HIC measurements for a given EOS, and ({{mathscr{P}}}_{{rm{MMA}}}({rm{EOS}})) is the posterior probability distribution on the EOS based on the multi-messenger observations, which acts as prior for this analysis. From the HIC experiments, we obtain a posterior of the pressure at a given density, (p(P|n,{rm{HIC}})). By combining this with the distribution of probed densities from the neutrons-over-charged particles sensitivity curve^{11}, (p(n|{rm{HIC}})), the joint posterior (p(n,P|{rm{HIC}})=p(P|n,{rm{HIC}})p(n|{rm{HIC}})) is obtained. Therefore, the relative faithfulness of the experimental results at various densities is accounted for. The likelihood ({ {mathcal L} }_{{rm{HIC}}}({rm{EOS}})) is given by

$$begin{array}{cc}{{mathcal{L}}}_{{rm{H}}{rm{I}}{rm{C}}}({rm{E}}{rm{O}}{rm{S}}) & =int {rm{d}}n,{rm{d}}P,p({rm{H}}{rm{I}}{rm{C}}|n,P)p(n,P|{rm{E}}{rm{O}}{rm{S}})\ & propto int {rm{d}}n,{rm{d}}P,p(n,P|{rm{H}}{rm{I}}{rm{C}})p(n,P|{rm{E}}{rm{O}}{rm{S}})\ & propto int {rm{d}}n,{rm{d}}P,p(n,P|{rm{H}}{rm{I}}{rm{C}})delta (P-P(n,{rm{E}}{rm{O}}{rm{S}}))\ & =int {rm{d}}n,p(n,P=P(n;{rm{E}}{rm{O}}{rm{S}})|{rm{H}}{rm{I}}{rm{C}}),end{array}$$

(11)

in which we used the pressure as a function of density for a given EOS.