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.

Main

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 × 1014 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-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 gravitational-wave (GW) astronomy[1] and new electromagnetic 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-26]. However, these observations mainly probe the EOS at densities and still carry considerable uncertainties, reflected in the ranges for predictions of neutron-star radii. More precise or new complementary information is 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 heavy-ion 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.

Fig. 1

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 c. In b–d, the 95% prior bound is shown for comparison (grey dashed lines). e, f, Posterior distributions for the pressure at 1.5nsat (e) and 2.5nsat (f) at different stages of our analysis, with the combined Astro + HIC region shaded in light blue.

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 freedom 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.

Extended Data Table 1

Impact of the EOS prior: Maximum density of chiral EFT and of prior distribution of R1.4

Impact of the EOS prior: Maximum density of chiral EFT and of prior distribution of R1.4

Comparison of the 95% credible interval for the pressure [MeV fm−3] and radius [km] of neutron stars when including only HIC experiments, only astrophysical observations, or the combined HIC and astrophysics results for chiral EFT constraints up to 1.5nsat and up to 1nsat, and for using a natural and uniform prior on R1.4. We find that differences for pressures and neutron-star radii are small between both prior choices when Astro+HIC data constraints are employed. Applying constraints from chiral EFT only up to 1nsat allows for a broader and stiffer EOS prior at higher densities since information up to 1.5nsat is discarded. As a consequence, the EOSs including HIC only and to a lesser extent the combination of HIC and observational constraints become stiffer leading to an increase of neutron-star radii. This effect is larger when using a natural instead of a uniform prior in radius. Nevertheless, the impact for the natural prior is only around 5%.

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 (), to remove EOSs that support only neutron stars with maximum masses well below the lower limit from the combined observations of heavy pulsars[36-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.

Fig. 2

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 (e) and (f) stars at different stages of our analysis, with the combined Astro + HIC region shaded in light blue.

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, 197Au 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 197Au 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 pressure-density 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 Table 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 Table 3).

Extended Data Table 2

Impact of EOS extension scheme

Impact of EOS extension scheme

Comparison of the 95% credible interval for the pressure [MeV fm−3] and radius [km] of neutron stars when including only HIC results, only astrophysical observations, and for combined HIC and astrophysics results for different EOS extension schemes used. The piecewise-polytrope scheme extends the EOS beyond 1.5nsat with five polytropic segments with randomly chosen transition densities and polytropic indices. The differences of the pressure estimates between the two EOS extension schemes is less than 3% and the difference between the radius estimates is less than 0.5% for the combined HIC and astrophysics results.

Extended Data Fig. 2

Comparison between different sensitivity curves.

We show three sensitivity-to-density curves for different observables and incident energies. In particular, the neutron-over-charged-particle (n/ch, used here) and the neutron-over-proton (n/p) sensitivity curves for 400 MeV/nucleon incident energy from Russotto et al.[11] are compared with the density curve reported by Le Fèvre et al.[10] for the sensitivity of the elliptic flow of protons in Au+Au collisions at 1 GeV/nucleon.

Extended Data Table 3

Sensitivity on uncertainties of nuclear matter properties

Sensitivity on uncertainties of nuclear matter properties

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−3, B = 15.86 ± 0.57 MeV[15] and K = 215 ± 40 MeV from microscopic calculations[15,91,92], which is in good agreement with the FOPI results.

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 final results (Extended Data Table 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).

Extended Data Table 4

Impact of the sensitivity curve for the ASY-EOS experiment

Impact of the sensitivity curve for the ASY-EOS experiment

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.

We note that results of the ASY-EOS experiment, in their sub-saturation density extension, are compatible with recent experimental findings from isobaric analogue 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 microscopic 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 behaviour, 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 () 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 , 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 , on the other hand, is also constrained by HIC information. Our final result for a typical neutron star is at 68% uncertainty ( at 95% uncertainty; Table 1). Comparing this value to the result without any HIC information, 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.

Table 1

Final constraints on the pressure and the radius of neutron stars

	Prior	Astro alone	HIC alone	Astro + HIC
P_1.5n_sat	5.59−1.97+2.04	5.84−2.26+1.95	6.06−2.04+1.85	6.25−2.26+1.90
R1.4	11.96−1.15+1.18	11.93−0.75+0.80	12.06−1.18+1.13	12.01−0.77+0.78

Comparison of the pressure in MeV fm−3 at 1.5nsat and the radius in km of a neutron star (median with the 95% credible interval) when including only astrophysical constraints, only HIC experimental data, and the combination of both.

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 Figs. 1 and 2), and advancing HIC experiments to probe higher densities, above 2–3nsat, will be key (Extended Data Table 5). Combining the latter with a reduction of experimental uncertainties, data from HICs have great potential to provide complementary 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.

Extended Data Fig. 1

Constraints on γasy versus symmetry energy S0 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 S0. Overall, the models are in good agreement with each other and the results suggest that a linear interpolation is reasonable.

Extended Data Table 5

Future prospects for the ASY-EOS experiment

Future prospects for the ASY-EOS experiment

Comparison of the 95% credible interval for the pressure [MeV fm−3] and the radius [km] of a neutron star when including only HIC experiments and for combined HIC and astrophysics results for different future improvements. In particular, we show the results with the 1 GeV sensitivity curve (see Extended Data Fig. 2) applied to the current measurement (second column), when additionally halving the uncertainty on γasy (third column), and when using a lower cutoff density of 1nsat instead of 0.5nsat (fourth column). For all the exploratory setups, HIC data is showing a stronger impact on the EOS constraint than the current setup. The result with a density cutoff is showing a significant decrease in uncertainty compared to the result of this work. Therefore, to achieve a stronger constraint on the EOS, improvements to the low-density part of the HIC constraint will be most important.

Methods

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-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 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 , 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.5 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.4 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 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 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 collapse to a black hole, so that an upper limit on the maximum mass can be derived. The combined estimate of the maximum mass, 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).

Extended Data Fig. 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 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.5nsat and . (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.

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 is given by[64]in which , 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,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 is given by[65]in which are the estimated or theoretically predicted apparent magnitudes for a given filter j (a passband for a particular wavelength interval) at observation time t with n data points for filter j. Moreover, and 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, , 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 (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 ΦRPin which all expansion coefficients v are functions of longitudinal rapidity , with p being the momentum along the beam axis and E being the total energy, and of transverse momentum of the particle, with p and p 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 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 (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-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-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 , 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 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 equation (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 equation (4) has been parameterized as

At saturation density, the kinetic part has been set to Ekin,0 = 12 MeV and Epot,0 = S0 − Ekin,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.19 for S0 = 31 MeV and γasy = 0.72 ± 0.19 for S0 = 34 MeV (see 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 behaviour 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.

Extended Data Fig. 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[15] 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.

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-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.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 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 conditionin 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 and chemical potential , 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 Table 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.

Extended Data Fig. 5

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.[15,59] at N2LO (light blue) and N3LO (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-3nsat, we additionally constrain the FOPI results with the constraint from Danielewicz et al.[27] (green), which has no statistical interpretation. This excludes 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 incompressibility 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.

Combination of the astronomical and HIC constraints

To combine the experimental and observational EOS constraints, we use Bayesian inference. The EOS posterior is given byin which MMA denotes multi-messenger astrophysics, is the likelihood of the HIC measurements for a given EOS, and 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, . By combining this with the distribution of probed densities from the neutrons-over-charged particles sensitivity curve[11], , the joint posterior is obtained. Therefore, the relative faithfulness of the experimental results at various densities is accounted for. The likelihood is given byin which we used the pressure as a function of density for a given EOS.

Online content

Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at 10.1038/s41586-022-04750-w.

Supplementary information

This file contains legends for Supplementary Tables 1–5.

Supplementary Table 1

Supplementary Table 2

Supplementary Table 3

Supplementary Table 4

Supplementary Table 5