Constraining neutron-star matter with microscopic and macroscopic collisions

Interpreting high-energy, astrophysical phenomena, such as supernova explosions or neutron-star collisions, requires a robust understanding of matter at supranuclear densities. However, our knowledge about dense matter explored in the cores of neutron stars remains limited. Fortunately, dense matter is not probed only in astrophysical observations, but also in terrestrial heavy-ion collision experiments. Here we use Bayesian inference to combine data from astrophysical multi-messenger observations of neutron stars1–9 and from heavy-ion collisions of gold nuclei at relativistic energies10,11 with microscopic nuclear theory calculations12–17 to improve our understanding of dense matter. We find that the inclusion of heavy-ion collision data indicates an increase in the pressure in dense matter relative to previous analyses, shifting neutron-star radii towards larger values, consistent with recent observations by the Neutron Star Interior Composition Explorer mission5–8,18. Our findings show that constraints from heavy-ion collision experiments show a remarkable consistency with multi-messenger observations and provide complementary information on nuclear matter at intermediate densities. This work combines nuclear theory, nuclear experiment and astrophysical observations, and shows how joint analyses can shed light on the properties of neutron-rich supranuclear matter over the density range probed in neutron stars.

Introduction.The nuclear equation of state (EOS) describes dense matter probed in terrestrial experiments with atomic nuclei as well as in astrophysical observations of neutron stars.The nuclear EOS is governed by quantum chromodynamics (QCD), the theory of strong interactions, but direct calculations of dense matter in neutron stars based on QCD are not feasible at present.Hence, the nuclear EOS has to be determined through approximate theoretical calculations or from experimental or observational data.As a result, at densities well above nuclear saturation density, nsat = 0.16 fm −3 (corresponding to a mass density of 2.7 • 10 14 g/cm 3 ), for which experimental and theoretical information are less robust, the nuclear EOS is still highly uncertain and many open questions remain, such as whether a possible phase transition to exotic phases of matter exists in nature 19 .
At densities below 1-2nsat, the EOS and its theoretical uncertainty can be obtained from microscopic calculations based on chiral effective field theory (EFT) of QCD [12][13][14][15][16][17] .To probe dense matter beyond these densities, further approaches, based on experimental and observational data, are necessary.A very promising tool is the multi-messenger astrophysics analysis of neutron stars and their collisions, which provides access to dense neutron-rich matter not accessible in terrestrial experiments at present.In recent years, the advent of gravitationalwave (GW) astronomy 1 and new electromagnetic (EM) observations of neutron stars 3,5,6 , including the Neutron Star Interior Composition Explorer (NICER) mission of the National Aeronautics and Space Administration (NASA) 5,6 , led to new constraints on the EOS 7,9,18,[20][21][22][23][24][25][26] .However, these observations mainly probe the EOS at densities 2nsat and still carry considerable uncertainties, reflected in the ranges for predictions of neutron-star radii.More precise or new complementary information are required to reduce the uncertainties further.
The gap between our current knowledge of the EOS stemming from nuclear theory and experiment at low densities and astrophysical observations of neutron stars at higher densities can be bridged by heavyion collision (HIC) experiments.These experiments, performed with heavy-ion beam energies of up to 2 GeV per nucleon, probe the nuclear EOS mainly in a density range of 1-2nsat at present 10,11,27 , representing a new source of information 28 .
In this work, we perform a global analysis of the nuclear EOS including information from nuclear theory (Fig. 1A), astrophysical observations of neutron stars (Fig. 1B), and results from HIC experiments that were performed at the Schwerionensynchrotron 18 accelerator located at the GSI Helmholtz Centre for Heavy Ion Research 10,11 (Fig. 1C).We analyse the EOS and neutron-star properties by extending our Bayesian multi-messenger astrophysics framework 9 to include information from the Four-Pi (FOPI) 10 and the Asymmetric-Matter EOS (ASY-EOS) experimental campaigns 11 .The combination of these experiments provides new constraints for neutron-rich matter in the range around 1-2nsat.We also include the EOS constraint from ref. 27 for symmetric nuclear matter obtained from HIC experiments at the Bevalac accelerator at Lawrence Berkeley National Laboratory and the Alternating Gradient Synchrotron at Brookhaven National Laboratory.In all experiments, gold nuclei were collided.The information from this series of HIC experiments allows us to further constrain the EOS in a density range for which theoretical calculations become less reliable.
Nuclear theory input.Our analysis starts with a set of 15,000 EOSs that are constrained by nuclear theory calculations at low densities.In particular, we utilize calculations using local chiral EFT interactions 14,29 .Chiral EFT is an effective theory of QCD that describes strong interactions in terms of nucleon and pion degrees of free- dom using a systematic momentum expansion of nuclear forces 30,31 .
In particular, the EFT expansion enables estimates of theoretical uncertainties 16,32 .On the basis of local chiral two-and three-nucleon interactions, we use quantum Monte Carlo methods, which are among the most precise many-body methods to solve the nuclear many-body problem 33 .The breakdown scale of the chiral EFT expansion was estimated to be about 500-600 MeV/c, in which c is the speed of light 16 .Therefore, we constrain our EOS set using chiral EFT input only up to 1.5nsat (corresponding to Fermi momenta of the order of 400 MeV/c), but a variation in the range 1-2nsat shows no substantial impact on our final results for neutron-star radii 34 (Extended Data Table 1).The 15,000 EOSs are sampled such that they span the theoretical uncertainty range of the chiral EFT calculation.
We extend each EOS above 1.5nsat using an extrapolation in the speed of sound (cs) in neutron-star matter 35 .This extrapolation is constrained only by causality (cs ≤ c) and stability of neutron-star matter (cs ≥ 0).In contrast to refs. 21,22  we do not take into account any information at asymptotically high densities from perturbative QCD calculations.In addition, at this level we require all EOSs in the prior to support neutron stars with masses of at least 1.9 solar masses (1.9M ), to remove EOSs that support only neutron stars with maximum masses well below the lower limit from the combined observations of heavy pulsars [36][37][38] .Hence, this lower bound ensures that the resulting EOS prior has reasonable support for massive-pulsar observations that we include at the first state of our Bayesian framework 9 .These general assumptions lead to a broad uncertainty for the EOS at higher densities (Fig. 1A), as well as for neutron-star masses and radii (Fig. 2A).The EOS prior is then used to analyse astrophysical observations and HIC experiments.
Multi-messenger astrophysics information.The astrophysical data are incorporated using a Bayesian multi-messenger framework 9,39 , which analyses each EOS with respect to its agreement with a variety of observational data.We start with the mass measurements of the massive neutron stars PSR J0348+0432 (ref. 36) and PSR J1614-2230 (ref. 37), to obtain a lower bound on the maximum mass, and the constraint on the maximum mass of neutron stars derived from the binary neutron-star collision GW170817 (refs. 40,41  in which a black hole was formed after the coalescence, to obtain an upper bound on the maximum mass.Information obtained from X-ray pulse-profile modelling of PSR J0030+0451 and PSR J0740+6620 using data from NICER and the X-ray Multi-Mirror Mission (XMM-Newton) 5,7,8 are incorporated.Moreover, we use Bayesian inference techniques to analyse GW information from the two neutron-star mergers GW170817 (ref. 1 ) and GW190425 (ref. 2 ) by matching the observed GW data with theoretical GW models that depend on neutron-star properties.For our analysis, we use a GW model 42 that is an improved version of the main waveform model used by the Laser Interferometer Gravitational-Wave Observatory/Virgo Collaboration for the study of GW170817 (ref. 43) and GW190425 (ref. 2 ).Similarly to the GW analysis, we also include information from the kilonova AT2017gfo (ref. 3 ) associated with the GW signal.Kilonovae originate from the radioactive decay of heavy atomic nuclei created in nucleosynthesis processes during and after the merger of neutron stars, and are visible in the optical, infrared and ultraviolet spectra.The electromagnetic observations are analysed with full radiative transfer simulations 44 to extract information from the observed light curve and spectra 4 .The above astrophysical information leads to important constraints on the neutron-star EOS, as shown in Fig. 1B.The constraints are strongest above 1.5nsat, for which the extrapolation in the speed of sound is used for the EOSs.The high-density astrophysical constraints affect mostly the high-mass region in the mass-radius plane and exclude the stiffest EOSs that lead to the largest radii.(Fig. 2B).
Data from HIC experiments.To further constrain the EOS, we implement data from HIC experiments.The FOPI 10 and ASY-EOS 11 experiments performed at GSI provide information respectively on the symmetric nuclear matter EOS (that is, matter with the same amount of protons and neutrons), and on the symmetry energy, which describes the energy cost of changing protons into neutrons in nuclear matter.For both experiments, 197 Au nuclei were collided at relativistic energies (0.4 to 1.5 GeV per nucleon), forming an expanding fireball in the collision region.This expansion is dictated by the achieved compression and therefore depends on the EOS of hot and dense matter.Owing to the initial neutron-to-proton asymmetry of the Au-Au system, the expansion of the emitted nucleons is sensitive to the nuclear symmetry energy.Constraints on the symmetry energy (from ASY-EOS) can be translated into a constraint on the pressure of neutron-star matter as a function of the baryon density when empirical information on symmetric nuclear matter from experiments (FOPI) with atomic nuclei is used.In addition to the GSI experiments, we include constraints on the pressure of symmetric nuclear matter at larger densities obtained from model calculations of ref. 27 that were used to analyse experimental data from Lawrence Berkeley National Laboratory and Brookhaven National Laboratory in which 197 Au nuclei were collided at energies up to 10 GeV per nucleon.These are sensitive to higher densities, 2-4.5nsat, but we include their constraints only up to 3nsat, where the sensitivity of the ASY-EOS experiment ends.We find that the inclusion of this further constraint has only minimal impact, but keep it to ensure the completeness of our study (Supplementary Information).
In Fig. 1C, we show the combined HIC experimental constraints (labelled HIC Data) at 68% and 95% credibility as well as the resulting posterior distribution for the neutron-star EOS.Whereas the FOPI experiment delivers an EOS constraint for symmetric nuclear matter at densities in the range 1-3nsat, the ASY-EOS experiment probes the symmetry energy roughly between 1 and 2nsat.The HIC pressuredensity constraint includes various sources of uncertainties.First, it includes systematic and statistical uncertainties of the experiments and the analysis of its data 10,11 .We have explicitly checked the robustness of our results when varying the details of the analysis and models used, and generally found that our results do not substantially depend on individual model choices (Extended Data Tabs.2, Extended Data Fig. 2 and Supplementary Information).Second, when extracting the HIC constraint on neutron-star matter, we vary nuclear matter properties, such as the incompressibility parameter and the symmetry energy at nsat, according to the measurements from FOPI and ASY-EOS.We have explicitly checked that increasing these uncertainties in agreement with theoretical estimates 17 leads to only minor changes of our final results (Extended Data Tab.3).
To enforce the ASY-EOS constraints only at densities for which the experiment is sensitive, we use the sensitivity curve for neutrons and charged particles 11 as a prior for the probed density range.We have checked the variation of our results for alternative choices of the sensitivity curve 11 and found that this has no substantial impact on our Comparison of the pressure in MeV fm −3 at 1.5nsat and the radius in km of a 1.4M neutron star (median with the 95% credible interval) when including only astrophysical constraints, only HIC experimental data, and the combination of both.
final results (Extended Data Tab.4).We find that the HIC constraints tend to prefer EOSs stiffer than the ones favoured by astrophysical observations (that is, EOSs that have higher pressures at densities up to 2nsat; Fig. 1C, E).We note that results of the ASY-EOS experiment, in their subsaturation density extension, are compatible with recent experimental findings from isobaric analog states supplemented with further constraints from neutron-skin data 45 , HICs using isospin-diffusion observables measured in mid-peripheral collisions of Sn isotopes 46 , and other nuclear structure information 47,48 .More recently, the Sπrit campaign at RIKEN has identified spectral yield ratios of charged pions in collisions of various tin isotopes near threshold as sensitive probes of the slope of the symmetry energy near and beyond nuclear saturation density 49 .The obtained value is compatible with the ASY-EOS result but offers no further strong constraint at present owing to its large uncertainty 49,50 .
Combining micro-and macroscopic collisions.The final EOS constraints are obtained through the combination of both the HIC information and astrophysical multi-messenger observations (Fig. 1D).The multi-messenger data rule out the most extreme EOS behavior, and the HIC data favour larger pressures around 1-1.5nsat, for which the experimental sensitivity is highest.This is similar to the effect of recent NICER observations on the EOS 7,18 .Hence, the two complementary approaches, HIC experiments and astrophysical observations show a remarkable agreement (Fig. 1E).At low densities, HIC results have a clear impact on the total posterior for the EOS, whereas the EOS at higher densities ( 2nsat) is mostly determined by astrophysical observations.At these densities, HIC results deviate only mildly from the prior (Fig. 1F).This is also reflected in the radii of neutron stars shown in Fig. 2E, F. As astrophysical observations mainly probe neutron stars with M 1.4M , for which the probed densities are higher, HIC information influences the radii of these neutron stars to a smaller degree.The radius of low-mass stars with M ≈ 1.0M , on the other hand, is also constrained by HIC information.Our final result for a typical 1.4M neutron star is 12.01 +0.37 −0.38 km at 68% uncertainty (12.01 +0. 780.77 km at 95% uncertainty; Tab. 1).Comparing this value to the result without any HIC information, 11.93 +0.39 −0.41 km at 68% confidence, highlights the benefit of combining these various sources of information in a statistically robust framework.We find that the HIC information has a high impact on the EOS at densities below 1.5nsat (Supplementary information).Finally, we quantify the possibility for the presence of a strong first-order phase transition to a new phase of QCD matter in the core of neutron stars.For this, we calculate the Bayes factor in favour of the presence of such a phase transition against its absence, and find it to be 0.419 ± 0.012 < 1.Therefore, its presence is slightly disfavoured given current astrophysical and experimental data.
To summarize, the interdisciplinary analysis of EOS constraints from HIC experiments and multi-messenger astrophysics shows remarkable agreement between the two, and provides important information to constrain the nuclear EOS at supra-saturation densities.Going forward, it is important that both statistic and systematic sources of uncertainty for HIC experiments are further improved.For example, the impact of choosing different quantum molecular dynamics models when analysing HIC experiments needs to be further investigated (Extended Data Fig. 1 and 2), and advancing HIC experiments to probe higher densities, above 2-3nsat, will be key (Extended Data Tab.5).Combining the latter with a reduction of experimental uncertainties, data from HICs have great potential to provide complimentary EOS information, bridging nuclear theory and astrophysical observations.In the next few years, the ASY-EOS-II and Compressed Baryonic Matter experiments at the upcoming Facility for Antiproton and Ion Research at GSI will provide a unique opportunity to study nuclear matter at densities probed in the core of neutron stars and their mergers, and might detect new phases of QCD matter, possibly involving hyperons and, ultimately, the transition to a deconfined quark matter phase at the highest densities (see, for example, refs. 51,52 .Together with experiments at the Rare Isotope Beam Facility at RIKEN in Japan and the Nuclotron-based Ion Collider Facility in Russia, the robust combination of experimental HIC constraints and astrophysical observations has the potential to revolutionize our understanding of the EOS.

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 latter 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 manybody methods to solve the nuclear many-body problem.The results of these calculations provide constraints on the EOS up to densities of around 2nsat (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 2nsat, indicated by increasing uncertainty estimates between 1 and 2nsat.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.5nsat to constrain the EOS, but a variation in the range 1-2nsat 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.5nsat and 12nsat, enforcing 0 ≤ cs ≤ 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.4M neutron star and the 'natural' prior and found only minor changes once astrophysical and HIC data are included (Extended Data Tab.1).
Recently, first results for the EOS of symmetric nuclear matter between 3 and 10nsat 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 highmass 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 neutronstar 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 (O(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.10 −0.13 M 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 (ref. 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 neutronstar 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 crosscorrelation with the observed data 1 .IMRPhenomPv2_NRTidalv2 is an updated model of the waveform model used in previous analyses by the 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 LGW is given by 64 LGW in which d(f ), h(f ), and Sn(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,048Hz.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 LEM is given by 65 in which m j,est i are the estimated or theoretically predicted apparent magnitudes for a given filter j (a passband for a particular wavelength interval) at observation time ti with nj data points for filter j.Moreover, m j i and σ j i 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, χ 2 1 , 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 dynamicalejecta 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 (v2) 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 dσ(y, pt) dΦ =C(1 + 2v1(y, pt) cos(Φ − ΦRP) in which all expansion coefficients vn are functions of longitudinal rapidity y = 1 2 ln E+pz E−pz , with pz being the momentum along the beam axis and E being the total energy, and of transverse momentum pt = p 2 x + p 2 y of the particle, with px and py denoting the momentum components perpendicular to the beam axis.
In the experiment, the orientation of the reaction plane is eventwise 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 v2 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 v2 of protons emitted at rapidities intermediate between projectile and target rapidity (midrapidity) 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 v1 to the stiffness of the EOS of symmetric nuclear matter becomes comparable to that of v2.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 3nsat 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 v2 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 np 2 = v n 2 /v p 2 , 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 midrapidity emissions over a large pt interval.Its isotopic resolution in this experiment was not sufficient to uniquely identify protons.Elliptic flow ratios as a function of pt 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 2nsat.

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 with the baryon density n = nn + np and the isospin asymmetry δ = (nn − np)/n = 1 − 2x, in which nn and np are the neutron and proton densities, respectively, and x = np/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 Eq. ( 4) has been parameterized with with the saturation density nsat, the Fermi energy EF , and in which the parameters α, β, and γ are fixed by the incompressibility K, the binding energy B of symmetric nuclear matter at nsat, and the condition that the pressure of symmetric nuclear matter is zero at saturation density: In the ASY-EOS analysis, the S(n) term of Eq. ( 4) has been parameterized as At saturation density, the kinetic part has been set to E kin,0 = 12 MeV and Epot,0 = S0 − E kin,0 .The parameter γasy was extracted from fits to experimental data of the pt dependence of the elliptic flow ratio of neutrons over charged particles around mid-rapidity.In particular, this results in γasy = 0.68±0.19for S0 = 31 MeV and γasy = 0.72±0.19for S0 = 34 MeV (Extended Data Fig. 4 for a comparison with microscopic neutron matter calculations).Here we interpolate γasy assuming a linear function with S0, for which the uncertainty is chosen to be 0.19 independent of S0.We have studied the behavior of γasy as a function of S0 for two different QMD models (Extended Data Fig. 1) and confirmed that the linear interpolation in the S0 range is suitable.The pressure constraint is given by the density derivative of the energy per particle of equation ( 4), and depends on n, δ, nsat, B, K, and S0.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-overcharged-particles flow ratio determined for the ASY-EOS experiment 11 (see also the previous section).This sensitivity curve covers the density range from 0.5nsat up to 3nsat and peaks between nsat and about 2nsat, 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 xASY-EOS accordingly.For simplicity, we neglect muons because they only have 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, 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 Pe = Ee/(3V ), energy density Ee/V = c(3π 2 ne) 4/3 /(4π 2 ), and chemical potential µe = c(3π 2 ne) 1/3 , in which is the reduced Planck constant and V the volume.
The final pressure constraint is obtained using EF = 37 MeV and by varying the parameters nsat, B, K, and S0 within specific ranges.For the parameters describing symmetric nuclear matter, we use the values consistent with the FOPI analysis given by nsat = 0.16 fm −3 , B = 16 MeV, and a Gaussian distribution for K with K = 200 ± 25 MeV at 1σ. Regarding S0, 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 3nsat (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 Tab.3).In particular, we have used a larger range for S0 between 30 and 35 MeV and the following Gaussian distributions for nsat, B and K: nsat = 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 p(EOS|MMA, HIC) ∝ p(HIC|EOS) × p(MMA|EOS)p(EOS) = p(HIC|EOS)p(EOS|MMA) ≡ LHIC(EOS)PMMA(EOS), (10)   in which MMA denotes multi-messenger astrophysics, LHIC(EOS) is the likelihood of the HIC measurements for a given EOS, and PMMA(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, HIC).By combining this with the distribution of probed densities from the neutrons-over-charged particles sensitivity curve 11 , p(n|HIC), the joint posterior p(n, P |HIC) = p(P |n, HIC)p(n|HIC) is obtained.Therefore, the relative faithfulness of the experimental results at various densities is accounted for.The likelihood LHIC(EOS) is given by LHIC(EOS) = dn dP p(HIC|n, P )p(n, P |EOS) ∝ dn dP p(n, P |HIC)p(n, P |EOS) ∝ dn dP p(n, P |HIC)δ(P − P (n, EOS)) = dn P (n, P = P (n; EOS)|HIC) , (11)   in which we used the pressure as a function of density for a given EOS.bilby and parallel bilby software packages are available at https:// git.ligo.org/lscsoft/bilbyand https://git.ligo.org/lscsoft/parallel_bilby, respectively.
The gwemlightcurve software is available at https://gwemlightcurves. github.io/.The exact code versions of bilby, parallel bilby, LAL-Suite, and gwemlightcurve that we have used for this work are also available (ref. 88 1 | Constraints on γ asy versus symmetry energy S 0 from two Quantum Molecular Dynamics models.We show the exponent γ asy of the density dependence of the potential part of the symmetry energy, see Eq. ( 7), as deduced from the analysis of ASY-EOS experimental data using the UrQMD model used in this work 11 (red points) and new simulations from the IQMD model (blue points).The red line indicates the mean value for γ asy along the linear interpolation for the chosen range of S 0 .Overall, the models are in good agreement with each other and the results suggest that a linear interpolation is reasonable.We show three sensitivity-to-density curves for different observables and incident energies.and GP-B 16 , and low-density quantum Monte Carlo results from Gezerlis and Carlson 89 .Overall, the results from these calculations are in good agreement with each other.We also show the energy per particle of a unitary Fermi gas of neutrons, which has been proposed as a lower bound for the energy of neutron matter 90 .Finally, we compare the theoretical results with the constraint from the ASY-EOS and FOPI experiments (red), which is used as a constraint for neutron matter in the main work.Comparison of the 95% credible interval for the pressure [ MeV fm −3 ] and radius [km] of neutron stars when using ranges for nuclear matter properties as published for the FOPI and ASY-EOS experiments 10,11 and when inflating the uncertainties according to theoretical calculations.We present results when including only information from HIC experiments and for the combined HIC and astrophysics information.
In particular, we extend the range for the symmetry energy at saturation density to S0 = 30 − 35 MeV by extrapolating γasy linearly.We use Gaussian distributions for nsat, B, and K describing symmetric nuclear matter and vary these parameters within their empirical ranges (at 1σ): nsat = 0.164 ± 0.007 fm Comparison of the 95% credible interval for the pressure [ MeV fm −3 ] and radius [km] of neutron stars when including only HIC experiments and for combined HIC and astrophysics results for different sensitivity curves.In particular, we compare our standard results using the neutron over charged particles (n/ch) sensitivity curve 11 with the neutron over proton (n/p) sensitivity from Russotto et al. 11 , which peaks at 1.5nsat.We find that our results are robust and differences for both sensitivity curves are small.Additionally, we compare the results to calculations where the ASY-EOS data is implemented using a uniform prior in density between 1-2nsat (labelled Window).For the latter choice, we generally find larger pressures and larger neutron-star radii because the n/ch and n/p sensitivy curves decrease rapidly after their maxima at 1 and 1.5nsat, lowering the impact of the ASY-EOS constraint at higher densities.However, differences for radii and pressures remain small once Astro+HIC data is included.

Figure 1 |
Figure1| Constraints on the EOS of neutron-star matter.A-D, Evolution of the pressure as a function of baryon number density for the EOS prior (A, grey), when including only data from multi-messenger neutron-star observations (B, green), when including only HIC data (C, orange), and when combining both (D, blue).The shading corresponds to the 95% and 68% credible intervals (lightest to darkest).The impact of the HIC experimental constraint (HIC Data, purple lines at 95% and 68%) on the EOS is shown in panel C.In B-D, the 95% prior bound is shown for comparison (grey dashed lines).E,F Posterior distributions for the pressure at 1.5n sat (E) and 2.5n sat (F) at different stages of our analysis (E, F), with the combined Astro+HIC region shaded in light-blue.

Figure 2 |
Figure2| Constraints on the mass and radius of neutron stars.A-D, The 95% and 68% credible ranges for the neutron-star radius across various masses (up to the 95% upper bound on the maximum allowed mass, as only few EOSs support mass beyond that, which would result in an unrepresentative credible range) for the prior (A, grey), when including only multi-messenger constraints (B, green), when including only HIC experiment data (C, orange), and for the joint constraint (D, blue).The prior 95% contour is shown in B-D for comparison.E,F, Posterior distributions for the radii of 1.4M (E) and 2M (F) stars are given at different stages of our analysis, with the combined Astro+HIC region shaded in light blue.

Data Figure 2 |
Comparison between different sensitivity curves.

3 |
Constraint on the neutron-star mass and radius with successive astrophysics information.In each panel (except for panel A), EOSs within (outside of) 95% credible interval are shown as blue (grey) lines.Lower panels indicate the probability distribution function (PDF) for the radius of a 1.4M neutron star, with the 95% confidence range indicated by dashed lines, in panels (B)-(F) the prior from panel (A) is shown in grey.(A) The EOS prior set constrained by chiral EFT calculations up to 1.5n sat and M max ≥ 1.9M .(B) The EOS set restricted by incorporating information from mass measurements of PSR J0348+0432, PSR J1614-2230, and the maximum-mass constraints obtained from GW170817/AT2017gfo.The 95% confidence interval of the maximum mass posterior probability distribution is shown by the purple band.(C) The EOS set further restricted by the NICER mass-radius measurement of PSR J0030+0451 (purple contours at 68% and 95% confidence) and PSR J0740+6620 (orange contours at 68% and 95% confidence).Note that the latter shows the average of the results obtained by Miller et al. 7 and Riley et al. 8 (D) Further restrictions on the EOS set from a reanalysis of the GW170817 using Bayesian inference.Contours at 68% and 95% confidence show the mass-radius measurements of the primary (red) and secondary (orange) neutron stars.(E) We use the chirp mass, mass ratio, and the EOSs as Bayesian prior for our analysis of AT2017gfo.(F) Further restrictions by analysing GW190425.Again, contours at 68% and 95% confidence show the mass-radius measurements of the primary (red) and secondary (orange) neutron stars.(G) The radius constraint at each step of this analysis with 95% confidence ranges.The radius constraint after including HIC experimental data is also shown.Neutron E/N [MeV] Hebeler et al ., ApJ (2013) Tews et al ., PRL (2013) Lynn et al ., PRL (2016) Drischler et al ., PRL (2019) Drischler et al ., GP-B (2020) Gezerlis, Carlson, PRC (2010) Unitary gas (ξ = 0.376) ASY-EOS + FOPI Extended Data 4 | Constraints for pure neutron matter.Energy per particle E/N of neutron matter as a function of density n for various many-body calculations using chiral EFT interactions from Hebeler et al. 12 , Tews et al. 13 , Lynn et al. (used here) 14 , Drischler et al.PRL

Extended Data Figure 5 |
., Science (2002) Drischler et al ., (2019) Comparison of the pressure of symmetric nuclear matter for experiment and theory.The pressure band from the FOPI experiment 10 at the 1σ level (red) for the incompressibility is consistent with the chiral EFT constraint from Drischler et al.

Table 1 |
Final constraints on the pressure and the radius of neutron stars.
). laboration Agreement between IN2P3 -DSM/CEA and GSI.KA acknowledges the support from the Bundesministerium für Bildung und Forschung (BMBF, German Federal Ministry of Education and Research) -Project-ID 05P19VTFC1 and Helmholtz Graduate School for Hadron and Ion Research (HGS-HIRe).Computations were performed on the national supercomputer Hawk at the High Performance Computing Center Stuttgart (HLRS) under the grant number 44189 and on SuperMUC-NG at Leibniz Supercomputing Centre Munich under project number pn29ba.In addition, computational resources have also been provided by the Los Alamos National Laboratory Institutional Computing Program, which is supported by the U.S. Department of Energy National Nuclear Security Administration under Contract No. 89233218CNA000001, and by the National Energy Research Scientific Computing Center (NERSC), which is supported by the U.S. Department of Energy, Office of Science, under contract No. DE-AC02-05CH11231.This research has made use of data, software and/or web tools obtained from the Gravitational Wave Open Science Center (https://www.gw-openscience.org),a service of LIGO Laboratory, the LIGO Scientific Collaboration and the Virgo Collaboration.LIGO is funded by the U.S. National Science Foundation.Virgo is funded by the French Centre National de Recherche Scientifique (CNRS), the Italian Istituto Nazionale della Fisica Nucleare (INFN) and the Dutch Nikhef, with contributions by Polish and Hungarian institutes.Authors contributions: S. Huth and P.T.H. Pang share first authorship.Conceptualization: SH, PTHP, IT, TD, ALF, AS; Methodology: SH, PTHP, IT, TD, ALF, AS, MB, MWC; Data curation: SH, PTHP, IT, TD, MWC; Software: SH, PTHP, IT, TD, ALF, MB, MWC; Validation: SH, PTHP, IT, TD, ALF, AS, WT; Formal analysis: SH, PTHP, IT, TD; Resources: IT, TD, ALF, AS; Funding acquisition: IT, TD, ALF, AS, CVDB; Project administration: IT, TD, AS; Supervision: IT, TD, ALF, AS, CVDB; Visualization: SH, PTHP; Writing-original draft: SH, PTHP, IT, TD, ALF, AS, WT, KA; Writing-review and editing: SH, PTHP, IT, TD, ALF, AS, WT, KA.Competing interests: The authors declare no competing interests.Additional Information: Supplementary Information is available for this paper.Correspondence and requests for materials should be addressed to Sabrina Huth and Peter T.H. Pang.Reprints and permissions information is available at www.nature.com/reprints.
27, 59at N 2 LO (light blue) and N 3 LO (dark blue).The experimental uncertainty band is smaller than the theoretical one because the empirical saturation point used for extracting the experimental results has smaller uncertainties compared to theoretical estimates from chiral EFT.Between 2-3n sat , we additionally constrain the FOPI results with the constraint from Danielewicz et al.27(green), which has no statistical interpretation. This ees the highest values for the incompressibility K from the FOPI distribution and also influences symmetric matter at smaller densities, which depends on the range of K.However, both constraints are in very good agreement with each other and the impact of the additional Danielewicz et al. constraint is small in our analysis.Extended 1 | Impact of the EOS prior: Maximum density of chiral EFT and of prior distribution of R 1.4 .Extended Data Table3on uncertainties of nuclear matter properties.

Table 5 |
Future prospects for the ASY-EOS experiment.