Select Page

### 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 gluons30,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 estimates16,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 interactions14,53,54,55. We use quantum Monte Carlo methods33, 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 (gtrsim 2{n}_{{rm{sat}}}), 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 sound35 (see also ref. 57).

To construct the neutron-star EOS set, we first extend our chiral EFT calculation to β-equilibrium and add a crust58. 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 radii34. 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({M}_{odot }) neutron star and the ‘natural’ prior and found only minor changes once astrophysical and HIC data are included (Extended Data Table 1).

Recently, first results for the EOS of symmetric nuclear matter between 3 and 10nsat from functional renormalization group calculations that are based on QCD became available59. 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 analysis9 (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 observations9,60 of high-mass neutron stars provides a lower bound on the maximum mass of neutron stars. By contrast, an upper bound of the maximum mass is obtained from the observation of the merger remnant of the neutron-star merger GW170817 (ref. 41). Among other arguments, the observation of a bright, red kilonova component and the observation of a short gamma-ray burst 2 s after the merger of the two neutron stars indicate that the remnant experienced a delayed (({mathscr{O}}({rm{100ms}}))) collapse to a black hole, so that an upper limit on the maximum mass can be derived. The combined estimate of the maximum mass, ({2.21}_{-0.13}^{+0.10}{M}_{odot }) at 68% uncertainty, already provides important information about the internal structure of neutron stars and disfavours both too stiff and too soft EOSs (that is, EOSs with too large and too small pressures, respectively).

In the next step, we incorporate NICER’s mass and radius measurement of PSR J0030+0451 (ref. 5) and PSR J0740+6620 (refs. 7,8). NICER, located on board of the International Space Station, is a NASA telescope measuring the X-ray pulse profile of pulsars. By correlation of the observed profile and brightness with theoretical predictions, it is possible to extract information on the configuration (for example, on the location and properties of hotspots on the neutron-star surface, the rotation rate of the star, and its compactness, which determines the light bending around the pulsar). This information enables constraints on the pulsar’s mass and radius. In addition to NICER, the XMM-Newton telescope61,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 constraint9. Results shown in the main text are obtained using the parallel bilby software63 and the waveform model IMRPhenomPv2_NRTidalv2 (ref. 42) for cross-correlation with the observed data1. IMRPhenomPv2_NRTidalv2 is an updated model of the waveform model used in previous analyses by the Laser Interferometer Gravitational-Wave Observatory (LIGO)/Virgo Collaboration2,43 and, hence, allows for a more accurate measurement of tidal effects. The likelihood function for the GW analysis ({ {mathcal L} }_{{rm{GW}}}) is given by64

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

(1)

in which (mathop{d}limits^{ sim }(,f)), (mathop{h}limits^{ sim }(,f)) and 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 ({ {mathcal L} }_{{rm{EM}}}) is given by65

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

(2)

in which ({m}_{i}^{j,{rm{est}}}) are the estimated or theoretically predicted apparent magnitudes for a given filter j (a passband for a particular wavelength interval) at observation time ti with nj data points for filter j. Moreover, ({m}_{i}^{j}) and ({sigma }_{i}^{j}) are the observed apparent magnitude and its corresponding statistical uncertainties, respectively. For this analysis, the probability distribution of a chi-squared distribution with a degree of freedom of 1, ({chi }_{1}^{2}), is taken as the likelihood measurement. To reduce the systematic error of the kilonova modelling below the statistical error, a further uncertainty of 1 mag is added to the measurement error. To analyse AT2017gfo, we use the radiative transfer code POSSIS44 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 decomposition66 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 ejecta66 and the disk-wind ejecta9, are correlated to source parameters of the binary neutron-star system based on numerical relativity simulations9,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 energies27,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 remain70,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

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

(3)

in which all expansion coefficients vn are functions of longitudinal rapidity (y=frac{1}{2},mathrm{ln}left(frac{E+{p}_{z}}{E-{p}_{z}}right)), with pz being the momentum along the beam axis and E being the total energy, and of transverse momentum ({p}_{{rm{t}}}=sqrt{{p}_{x}^{2}+{p}_{y}^{2}}) of the particle, with px and py 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 procedure72. The coincident particle and fragment emissions are also used for the reconstruction of the impact parameter of each reaction event11. 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 EOS10,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–Uhlenbeck27 models. The origin of the phenomenon has been investigated in detail elsewhere76. 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 isotopes10 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 model75, 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 simulations10. 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 nucleon77,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 nuclei79,80,81,82). It has been suggested83 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}_{2}^{{rm{np}}}={v}_{2}^{{rm{n}}}/{v}_{2}^{{rm{p}}}), was proposed84. 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 elsewhere11. 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, IQMD74 and Tübingen QMD86), 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

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

(4)

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

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

(5)

with the saturation density 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:

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

(6)

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

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

(7)

At saturation density, the kinetic part has been set to 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.

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

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

(8)

and depends on n, δ, 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 experiment11 (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 condition

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

(9)

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

The final pressure constraint is obtained using 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 elsewhere27 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 energy86. We have explicitly checked the robustness of our results when using larger uncertainty ranges for all nuclear matter parameters in agreement with theoretical predictions17 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.

### Combination of the astronomical and HIC constraints

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

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

(10)

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

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

(11)

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