# NMMA: A nuclear-physics and multi-messenger astrophysics framework to analyze binary neutron star mergers

Peter T.H. Pang, Tim Dietrich, Michael W. Coughlin, Mattia Bulla, Ingo Tews, Mouza Almualla, Tyler Barna, Weizmann Kiendrebeogo, Nina Kunert, Gargi Mansingh, et al.

► To cite this version:

Peter T.H. Pang, Tim Dietrich, Michael W. Coughlin, Mattia Bulla, Ingo Tews, et al.. NMMA: A nuclear-physics and multi-messenger astrophysics framework to analyze binary neutron star mergers. Nature Commun., 2023, 14 (1), pp.8352. 10.1038/s41467-023-43932-6 . hal-03684723

HAL Id: hal-03684723

https://hal.science/hal-03684723

Submitted on 15 Feb 2024

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

# An updated nuclear-physics and multi-messenger astrophysics framework for binary neutron star mergers

Received: 23 February 2023

Accepted: 24 November 2023

Published online: 20 December 2023

<figure><img src="image_2.png" /></figure>

Check for updates

Peter T. H. Pang Ⓐ¹,², Tim Dietrich Ⓐ³,⁴✉, Michael W. Coughlin Ⓐ⁵, Mattia Bulla⁶,⁷,⁸,⁹, Ingo Tews¹⁰, Mouza Almualla Ⓐ¹¹, Tyler Barna⁵, Ramodgwendé Weizmann Kiendrebeogo Ⓐ¹²,¹³, Nina Kunert³, Gargi Mansingh⁵,¹⁴, Brandon Reed⁵,¹⁵, Niharika Sravan¹⁶, Andrew Toivonen Ⓐ⁵, Sarah Antier¹³, Robert O. VandenBerg⁵, Jack Heinzel¹⁷, Vsevolod Nedora⁴, Pouyan Salehi³, Ritwik Sharma¹⁸, Rahul Somasundaram¹⁰,¹⁹,²⁰ & Chris Van Den Broeck Ⓐ¹,²

The multi-messenger detection of the gravitational-wave signal GW170817, the corresponding kilonova AT2017gfo and the short gamma-ray burst GRB170817A, as well as the observed afterglow has delivered a scientific breakthrough. For an accurate interpretation of all these different messengers, one requires robust theoretical models that describe the emitted gravitational-wave, the electromagnetic emission, and dense matter reliably. In addition, one needs efficient and accurate computational tools to ensure a correct cross-correlation between the models and the observational data. For this purpose, we have developed the Nuclear-physics and Multi-Messenger Astrophysics framework NMMA. The code allows incorporation of nuclear-physics constraints at low densities as well as X-ray and radio observations of isolated neutron stars. In previous works, the NMMA code has allowed us to constrain the equation of state of supranuclear dense matter, to measure the Hubble constant, and to compare dense-matter physics probed in neutron-star mergers and in heavy-ion collisions, and to classify electromagnetic observations and perform model selection. Here, we show an extension of the NMMA code as a first attempt of analyzing the gravitational-wave signal, the kilonova, and the gamma-ray burst afterglow simultaneously. Incorporating all available information, we estimate the radius of a $1.4M_\odot$ neutron star to be $R = 11.98^{+0.35}_{-0.40}$ km.

The study of the gravitational-wave (GW) and electromagnetic (EM) signals GW170817, AT2017gfo²⁻¹², and GRB170817A¹³⁻¹⁵ has already enabled numerous scientific breakthroughs, for example, constraints on the properties of neutron stars (NSs) and the dense matter equation of state (EOS) at supranuclear densities¹⁶⁻²³, an independent measurement of the Hubble constant²²,²⁴⁻²⁸, the verified connection between binary NS (BNS) mergers and at least some of the observed

short gamma-ray bursts (GRBs)²⁹, and precise limits on the propagation speed of GWs²⁹. These scientific achievements were enabled by the multi-messenger nature of GW170817.

Despite this enormous progress, results have been obtained by connecting constraints from individual messengers a posteriori, i.e., different messengers were analyzed individually and then combined within different multi-messenger frameworks to achieve the final

<figure><img src="image_3.png" /></figure>
<figure><img src="image_2.png" /></figure>

# An updated nuclear-physics and multi-messenger astrophysics framework for binary neutron star mergers

Received: 23 February 2023

Accepted: 24 November 2023

Published online: 20 December 2023

<figure><img src="image_3.png" /></figure>

Check for updates

Peter T. H. Pang Ⓐ¹,², Tim Dietrich Ⓐ³,⁴✉, Michael W. Coughlin Ⓐ⁵, Mattia Bulla⁶,⁷,⁸,⁹, Ingo Tews¹⁰, Mouza Almualla Ⓐ¹¹, Tyler Barna⁵, Ramodgwendé Weizmann Kiendrebeogo Ⓐ¹²,¹³, Nina Kunert³, Gargi Mansingh⁵,¹⁴, Brandon Reed⁵,¹⁵, Niharika Sravan¹⁶, Andrew Toivonen Ⓐ⁵, Sarah Antier¹³, Robert O. VandenBerg⁵, Jack Heinzel¹⁷, Vsevolod Nedora⁴, Pouyan Salehi³, Ritwik Sharma¹⁸, Rahul Somasundaram¹⁰,¹⁹,²⁰ & Chris Van Den Broeck Ⓐ¹,²

The multi-messenger detection of the gravitational-wave signal GW170817, the corresponding kilonova AT2017gfo and the short gamma-ray burst GRB170817A, as well as the observed afterglow has delivered a scientific breakthrough. For an accurate interpretation of all these different messengers, one requires robust theoretical models that describe the emitted gravitational-wave, the electromagnetic emission, and dense matter reliably. In addition, one needs efficient and accurate computational tools to ensure a correct cross-correlation between the models and the observational data. For this purpose, we have developed the Nuclear-physics and Multi-Messenger Astrophysics framework NMMA. The code allows incorporation of nuclear-physics constraints at low densities as well as X-ray and radio observations of isolated neutron stars. In previous works, the NMMA code has allowed us to constrain the equation of state of supranuclear dense matter, to measure the Hubble constant, and to compare dense-matter physics probed in neutron-star mergers and in heavy-ion collisions, and to classify electromagnetic observations and perform model selection. Here, we show an extension of the NMMA code as a first attempt of analyzing the gravitational-wave signal, the kilonova, and the gamma-ray burst afterglow simultaneously. Incorporating all available information, we estimate the radius of a $1.4M_{\odot}$ neutron star to be
$$R = 11.98^{+0.35}_{-0.40} \text{ km.}$$

The study of the gravitational-wave (GW) and electromagnetic (EM) signals GW170817', AT2017gfo²⁻¹², and GRB170817A¹³⁻¹⁵ has already enabled numerous scientific breakthroughs, for example, constraints on the properties of neutron stars (NSs) and the dense matter equation of state (EOS) at supranuclear densities¹⁶⁻²³, an independent measurement of the Hubble constant²²,²⁴⁻²⁸, the verified connection between binary NS (BNS) mergers and at least some of the observed

short gamma-ray bursts (GRBs)²⁹, and precise limits on the propagation speed of GWs²⁹. These scientific achievements were enabled by the multi-messenger nature of GW170817.

Despite this enormous progress, results have been obtained by connecting constraints from individual messengers a posteriori, i.e., different messengers were analyzed individually and then combined within different multi-messenger frameworks to achieve the final

<figure><img src="image_4.png" /></figure>
results. Such frameworks and attempts include, among others, the work of Breschi et al.<sup>30</sup> performing Bayesian inference and model selection on the kilonova AT2017gfo, Nicholl et al.<sup>31</sup> developing a framework for predicting kilonova and GRB afterglow lightcurves using information from GW signals as input, and the multi-messenger framework developed by Raaijmakers et al.<sup>32</sup>. Similarly, to these works, our previous Nuclear physics - Multi-Messenger Astrophysics (NMMA) framework has been successfully applied to provide constraints on the EOS of NS matter and on the Hubble constant<sup>22,33</sup>, to investigate the nature of the compact binary merger GW190814<sup>34</sup>, to provide techniques to search for kilonova transients<sup>35</sup>, to classify observed EM transients such as GRB200826A<sup>36</sup>, and to combine information from multi-messenger observations with data from nuclear-physics experiments such as heavy-ion collisions<sup>23</sup>.

Here, we upgrade our framework to allow for a simultaneous analysis of kilonova, GRB afterglow, and GW data capitalizing on the multi-messenger nature of compact-binary mergers.

## Results

The full potential of our NMMA study becomes clear from Fig. 1 where we show a set of possible EOSs relating the pressure and baryon number density inside NSs. Different constraints can provide valuable information in different density regimes. For example, theoretical calculations of dense nuclear matter in the framework of chiral effective field theory (EFT)<sup>37–41</sup> or data extracted from nuclear-physics experiments, e.g., heavy-ion collisions<sup>42</sup> or the recent PREX-II experiment at Jefferson Laboratory<sup>43</sup>, provide valuable input up to about twice the nuclear saturation density, $n_{\text{sat}} \approx 0.16 \text{ fm}^{-3}$. GW signals emitted during the inspiral of a BNS or black-hole-NS (BHNS) systems contain information that probe the EOS at densities realized inside the individual NS components of the system, typically up to about five times $n_{\text{sat}}$, but the exact density range probed in such mergers depends noticeably on the mass of the component stars. Furthermore, radio observations of NSs can be used to infer their masses, e.g., by measuring Shapiro delay in a binary system. In particular, radio observations of heavy NSs with masses of about $2M_{\odot}$, such as PSR J0348+0432<sup>44</sup>, PSR J1614-2230<sup>45</sup>, and PSR J0740+6620<sup>46</sup>, currently provide valuable information at larger densities than those probed by

inspiral GW signals. In addition, these observations provide a valuable lower bound on the maximum mass of NSs. Matter at the highest densities in the universe could be created in the postmerger phase of a BNS coalescence, i.e., after the collision of the two NSs in the binary. This phase of the binary merger might be observed through future GW detections with more sensitive detectors. Alternatively, this phase can be probed by analyzing EM signals connected to a BNS merger, i.e., the kilonovae, GRBs, and their afterglows. Finally, at asymptotically high densities that are not shown in the figure, the EOSs can be calculated in perturbative QCD<sup>47</sup> and might be used to constrain the NS EOS<sup>48</sup>. The combination of all these various pieces of information provides a unique tool to unravel the properties of matter at supranuclear densities.

## GW170817-AT2017gfo

With the NMMA framework, we analyze GW170817 simultaneously with the observed kilonova AT2017gfo. For the GW analysis, we have used the IMRPhenomPv2_NRTidalv2 waveform model and analyzed the GW data obtained from the Gravitational Wave Open Science Center (GWOSC)<sup>49</sup> in a frequency range of 20 Hz to 2048 Hz, covering the detected BNS inspiral<sup>50</sup>. For the EM signal, we use the data set compiled in Coughlin et al.<sup>51</sup>, where in this work, we include the optical, infrared, and ultraviolet data between 0.5 and 10 days after the merger. The corresponding data is analyzed with a Gaussian Process Regression (GPR)-based kilonova model. For our analysis, we are presenting the best-fit lightcurve in Fig. 2, with its band representing a one

<figure><img src="image_4.png" /></figure>

**Fig. 2 | Best-fit early-time lightcurve from the analysis.** The best-fit lightcurves (dashed, with the 1 magnitude uncertainty shown as the band) for AT2017gfo data when analyzing GW170817-and-AT2017gfo (orange) or GW170817-and-AT2017gfo-and-GRB170817A (blue) simultaneously. We note that both bands overlap almost completely, i.e., for AT2017gfo the accuracy of the kilonova lightcurve description does not depend noticeably on the inclusion of a GRB afterglow component. For the analysis, we restrict our dataset to times between 0.5 days up to 10 days after the BNS merger to simplify the joint GW170817-and-AT2017gfo-and-GRB170817A study as discussed in the main text.

<figure><img src="image_5.png" /><figcaption>**Fig. 1 | Overview of constraints on the EOS from different information channels.**</figcaption></figure>

We show a set of possible EOSs (blue lines) that are constrained up to $1.5n_{\text{sat}}$ by Quantum Monte Carlo calculations using chiral EFT interactions<sup>80</sup> and extended to higher densities using a speed of sound model<sup>49</sup>. Different regions of the EOS can then be constrained by using different astrophysical messengers, indicated by rectangulars: GWs from inspirals of NS mergers, data from radio and X-ray pulsars, and EM signals associated with NS mergers. Note that the boundaries are not strict but depend on the EOS and properties of the studied system.
results. Such frameworks and attempts include, among others, the work of Breschi et al.<sup>30</sup> performing Bayesian inference and model selection on the kilonova AT2017gfo, Nicholl et al.<sup>31</sup> developing a framework for predicting kilonova and GRB afterglow lightcurves using information from GW signals as input, and the multi-messenger framework developed by Raaijmakers et al.<sup>32</sup>. Similarly, to these works, our previous Nuclear physics - Multi-Messenger Astrophysics (NMMA) framework has been successfully applied to provide constraints on the EOS of NS matter and on the Hubble constant<sup>22,33</sup>, to investigate the nature of the compact binary merger GW190814<sup>34</sup>, to provide techniques to search for kilonova transients<sup>35</sup>, to classify observed EM transients such as GRB200826A<sup>36</sup>, and to combine information from multi-messenger observations with data from nuclear-physics experiments such as heavy-ion collisions<sup>23</sup>.

Here, we upgrade our framework to allow for a simultaneous analysis of kilonova, GRB afterglow, and GW data capitalizing on the multi-messenger nature of compact-binary mergers.

## Results

The full potential of our NMMA study becomes clear from Fig. 1 where we show a set of possible EOSs relating the pressure and baryon number density inside NSs. Different constraints can provide valuable information in different density regimes. For example, theoretical calculations of dense nuclear matter in the framework of chiral effective field theory (EFT)<sup>37–41</sup> or data extracted from nuclear-physics experiments, e.g., heavy-ion collisions<sup>42</sup> or the recent PREX-II experiment at Jefferson Laboratory<sup>43</sup>, provide valuable input up to about twice the nuclear saturation density, $n_{\text{sat}} \approx 0.16 \text{ fm}^{-3}$. GW signals emitted during the inspiral of a BNS or black-hole-NS (BHNS) systems contain information that probe the EOS at densities realized inside the individual NS components of the system, typically up to about five times $n_{\text{sat}}$, but the exact density range probed in such mergers depends noticeably on the mass of the component stars. Furthermore, radio observations of NSs can be used to infer their masses, e.g., by measuring Shapiro delay in a binary system. In particular, radio observations of heavy NSs with masses of about $2M_{\odot}$, such as PSR J0348+0432<sup>44</sup>, PSR J1614-2230<sup>45</sup>, and PSR J0740+6620<sup>46</sup>, currently provide valuable information at larger densities than those probed by

inspiral GW signals. In addition, these observations provide a valuable lower bound on the maximum mass of NSs. Matter at the highest densities in the universe could be created in the postmerger phase of a BNS coalescence, i.e., after the collision of the two NSs in the binary. This phase of the binary merger might be observed through future GW detections with more sensitive detectors. Alternatively, this phase can be probed by analyzing EM signals connected to a BNS merger, i.e., the kilonovae, GRBs, and their afterglows. Finally, at asymptotically high densities that are not shown in the figure, the EOSs can be calculated in perturbative QCD<sup>47</sup> and might be used to constrain the NS EOS<sup>48</sup>. The combination of all these various pieces of information provides a unique tool to unravel the properties of matter at supranuclear densities.

## GW170817-AT2017gfo

With the NMMA framework, we analyze GW170817 simultaneously with the observed kilonova AT2017gfo. For the GW analysis, we have used the IMRPhenomPv2_NRTidalv2 waveform model and analyzed the GW data obtained from the Gravitational Wave Open Science Center (GWOSC)<sup>49</sup> in a frequency range of 20 Hz to 2048 Hz, covering the detected BNS inspiral<sup>50</sup>. For the EM signal, we use the data set compiled in Coughlin et al.<sup>51</sup>, where in this work, we include the optical, infrared, and ultraviolet data between 0.5 and 10 days after the merger. The corresponding data is analyzed with a Gaussian Process Regression (GPR)-based kilonova model. For our analysis, we are presenting the best-fit lightcurve in Fig. 2, with its band representing a one

<figure><img src="image_5.png" /></figure>

**Fig. 2 | Best-fit early-time lightcurve from the analysis.** The best-fit lightcurves (dashed, with the 1 magnitude uncertainty shown as the band) for AT2017gfo data when analyzing GW170817-and-AT2017gfo (orange) or GW170817-and-AT2017gfo-and-GRB170817A (blue) simultaneously. We note that both bands overlap almost completely, i.e., for AT2017gfo the accuracy of the kilonova lightcurve description does not depend noticeably on the inclusion of a GRB afterglow component. For the analysis, we restrict our dataset to times between 0.5 days up to 10 days after the BNS merger to simplify the joint GW170817-and-AT2017gfo-and-GRB170817A study as discussed in the main text.

<figure><img src="image_6.png" /><figcaption>**Fig. 1 | Overview of constraints on the EOS from different information channels.**</figcaption></figure>

We show a set of possible EOSs (blue lines) that are constrained up to $1.5n_{\text{sat}}$ by Quantum Monte Carlo calculations using chiral EFT interactions<sup>80</sup> and extended to higher densities using a speed of sound model<sup>49</sup>. Different regions of the EOS can then be constrained by using different astrophysical messengers, indicated by rectangulars: GWs from inspirals of NS mergers, data from radio and X-ray pulsars, and EM signals associated with NS mergers. Note that the boundaries are not strict but depend on the EOS and properties of the studied system.
magnitude uncertainty for the individual lightcurves. This one magnitude uncertainty is introduced during the inference and should account for systematic uncertainties in kilonova modeling. In Supplementary information, we show how smaller or larger assumed uncertainties change our conclusions and it show that the one-magnitude is a sensible choice. Such a finding is also consistent with Heinzel et al.⁵². Therefore, we focus particularly on one-magnitude uncertainties' results. Nevertheless further work would be needed to understand in detail uncertainties related to the ejecta geometry⁵², assumed heating rates, thermalization efficiencies and opacities within the ejecta³³. Furthermore, we point out that for Fig. 2, we explicitly restricted our data set to the times between 0.5 to 10 days after the BNS merger, since model predictions at earlier or later times are more uncertain, e.g., due to less accurate opacities during early times and a larger impact of Monte Carlo noise in the employed radiative transfer models at late times. While this does not affect the GW170817-AT2017gfo analysis, it has an impact when we will also incorporate the GRB afterglow. In fact, we find that not restricting us to this time ranges can cause problems in the joint inference and it takes noticeably longer until the sampler converges.

Figure 3 summarizes our main findings and shows joint posteriors for the mass of the dynamical ejecta $m_{dyn}^{\text{ej}}$, the mass of the disk wind ejecta $m_{\text{wind}}^{\text{ej}}$, the chirp mass $M_c$, the mass ratio $q$, the mass-weighted tidal deformability $\bar{\Lambda}$, and the radius of a 1.4 solar mass neutron star $R_{1.4}$. In contrast to previous findings using simpler kilonova modeling (see ref. 54 and references therein), we can fit AT2017gfo with masses for the dynamical (about 0.006 $M_{\odot}$) and disk-wind (about 0.07 $M_{\odot}$) ejecta components that are within the range of values predicted by numerical-relativity simulations⁵⁵. While the parameters extracted are consistent with our previous findings²², we observe a clear improvement on the parameter error bounds due to (1) performing a simultaneous analysis of the distinct messengers and (2) employing a modified likelihood function when analyzing the kilonova. For instance, the constraints on $R_{1.4} = 11.86^{+0.41}_{-0.53}$, a typical choice to quantify EOS constraints, is significantly improved compared to our previous result, $R_{1.4} = 11.75^{+0.86}_{-0.81}$ km². The half-width of $R_{1.4}$'s 90% credible interval decreases from about 800 m²² to about 400 m.

## GW170817-AT2017gfo-GRB170817A

In addition to the combined analysis of GW170817 and AT2017gfo, we can also incorporate information obtained from the GRB afterglow of GW170817A, where we employ the data set collected in Troja et al.⁵⁶. The GRB afterglow light-curve data are analyzed with the synthetic Gaussian jet-model lightcurve described before⁵⁷,⁵⁸. Figure 2 shows the corresponding best-fit lightcurve for the kilonova with a 1 magnitude uncertainty band as before. Moreover, we are also presenting the best-fit lightcurve, which includes kilonova and GRB afterglow, and the employed uncertainty band in Fig. 4. We find that both the kilonova AT2017gfo and the GRB afterglow GRB170817A are well described in our analysis.

Figure 3 again summarizes our findings for the joint posteriors of the mass of the dynamical ejecta, the mass of the disk wind ejecta, the on-axis isotropic equivalent energy, the chirp mass, the mass ratio, the mass-weighted tidal deformability, and the radius of a 1.4 solar mass neutron star for this analysis, which is consistent with GW170817-and-AT2017gfo only. Compared to the analysis of GW170817-and-AT2017gfo only, the improvement on the parameter uncertainties is minimal, yet, noticeable when information from GRB170817A is added. Although no significant constraint on the EOS is imposed by the jet energy $E_0$ as the ratio $\xi$ between it, $E_0$, and the disk mass $m_{\text{disk}}$ is taken as a free parameter, the inclination constraint from the GRB plays a role in the constraint on EOS. For an anisotropic kilonova model, the inclination angle changes the observable kilonova light curves beyond scaling (e.g., Fig. 2 in ref. 59), which is correlated with the ejecta masses

(e.g., Fig. 3 in ref. 60). Therefore the GRB's inclination measurement imposes a constraint on the EOS via the kilonova eject masses measurements.

Moreover, for future studies, we expect that the inclusion of the GRB afterglow will be of great importance for measuring the Hubble constant.

## Discussion

We have developed a publicly available NMMA framework for the interpretation and analysis of BNS and BHNS systems. This framework allows for the simultaneous analysis of GW and EM signals such as kilonovae and GRB afterglows. In addition, our framework allows us to incorporate constraints from nuclear-physics calculations, e.g., by sampling over EOS sets constrained by chiral EFT, and to include radio as well as X-ray measurements of isolated NSs. By employing our framework to a combined analysis of GW170817, AT2017gfo, and GRB170817A, we find that the radius of a typical 1.4 solar mass NS lies within $11.98^{+0.35}_{-0.40}$ km; cf. Table 1 for a selection of studies from the literature. Based on our findings, our analysis is a noticeable improvement over previous works. However, additional uncertainties in our work lie in limited physics input in kilonova and semi-analytic GRB and models. Therefore, reliable astrophysical interpretations of future BNS detections will only be possible if not only parameter estimation infrastructure, as presented in this work, but also the astrophysical models describing transient phenomena advance further. Nevertheless, given the increasing number of multi-messenger detections of BNS and BHNS merger, we expect to use our framework to further increase our knowledge about the interior of NSs during the coming years.

## Methods

### Equation of state construction

The EOS describes the relation between energy density $\epsilon$, pressure $p$, and temperature $T$ of dense matter and additionally depends on the composition of the system. For NSs, thermal energies are much smaller than typical Fermi energies of the particles, and therefore, temperature effects can be neglected for isolated NSs or NSs in the inspiral phase of a merger. In these cases, the EOS simply relates $\epsilon$ and $p$.

The most general constraints on the EOS can be inferred from the slope of the EOS, the speed of sound, defined as:

$$c_S = c \sqrt{\frac{\partial p}{\partial \epsilon}}, \qquad (1)$$

where $c$ is the speed of light. Due to the laws of special relativity, the speed of sound has to be smaller than the speed of light, $c_S \le c$. Furthermore, the speed of sound in a NS has to be larger than zero, $c_S \ge 0$, as NSs would otherwise be unstable. These constraints alone, however, allow for an extremely large EOS space.

At nuclear densities, additional information on the EOS can be inferred from laboratory experiments and theoretical nuclear-physics calculations. For example, this information was used to constrain the properties of stellar matter in the NS crust⁶¹,⁶², i.e., the outermost layer of NSs at densities below approximately $0.5n_{\text{sat}}$. Above roughly $0.5n_{\text{sat}}$, NS matter consists of a fluid of neutrons with a small admixture of protons. In this regime, the EOS can be constrained by microscopic calculations of dense nuclear matter. These calculations typically provide the energy per particle, $E/A(n, x)$, which is a function of density $n$ and proton fraction $x = n_p/n$ with $n_p$ being the proton density. From this, the EOS follows from:

$$\varepsilon(n,x) = n \frac{E}{A}(n,x), \qquad (2)$$
magnitude uncertainty for the individual lightcurves. This one magnitude uncertainty is introduced during the inference and should account for systematic uncertainties in kilonova modeling. In Supplementary information, we show how smaller or larger assumed uncertainties change our conclusions and it show that the one-magnitude is a sensible choice. Such a finding is also consistent with Heinzel et al.⁵². Therefore, we focus particularly on one-magnitude uncertainties' results. Nevertheless further work would be needed to understand in detail uncertainties related to the ejecta geometry⁵², assumed heating rates, thermalization efficiencies and opacities within the ejecta³³. Furthermore, we point out that for Fig. 2, we explicitly restricted our data set to the times between 0.5 to 10 days after the BNS merger, since model predictions at earlier or later times are more uncertain, e.g., due to less accurate opacities during early times and a larger impact of Monte Carlo noise in the employed radiative transfer models at late times. While this does not affect the GW170817-AT2017gfo analysis, it has an impact when we will also incorporate the GRB afterglow. In fact, we find that not restricting us to this time ranges can cause problems in the joint inference and it takes noticeably longer until the sampler converges.

Figure 3 summarizes our main findings and shows joint posteriors for the mass of the dynamical ejecta $m_{dyn}^{\text{ej}}$, the mass of the disk wind ejecta $m_{\text{wind}}^{\text{ej}}$, the chirp mass $M_c$, the mass ratio $q$, the mass-weighted tidal deformability $\bar{\Lambda}$, and the radius of a 1.4 solar mass neutron star $R_{1.4}$. In contrast to previous findings using simpler kilonova modeling (see ref. 54 and references therein), we can fit AT2017gfo with masses for the dynamical (about 0.006 $M_{\odot}$) and disk-wind (about 0.07 $M_{\odot}$) ejecta components that are within the range of values predicted by numerical-relativity simulations⁵⁵. While the parameters extracted are consistent with our previous findings²², we observe a clear improvement on the parameter error bounds due to (1) performing a simultaneous analysis of the distinct messengers and (2) employing a modified likelihood function when analyzing the kilonova. For instance, the constraints on $R_{1.4} = 11.86^{+0.41}_{-0.53}$, a typical choice to quantify EOS constraints, is significantly improved compared to our previous result, $R_{1.4} = 11.75^{+0.86}_{-0.81}$ km²². The half-width of $R_{1.4}$'s 90% credible interval decreases from about 800 m²² to about 400 m.

## GW170817-AT2017gfo-GRB170817A

In addition to the combined analysis of GW170817 and AT2017gfo, we can also incorporate information obtained from the GRB afterglow of GW170817A, where we employ the data set collected in Troja et al.⁵⁶. The GRB afterglow light-curve data are analyzed with the synthetic Gaussian jet-model lightcurve described before⁵⁷,⁵⁸. Figure 2 shows the corresponding best-fit lightcurve for the kilonova with a 1 magnitude uncertainty band as before. Moreover, we are also presenting the best-fit lightcurve, which includes kilonova and GRB afterglow, and the employed uncertainty band in Fig. 4. We find that both the kilonova AT2017gfo and the GRB afterglow GRB170817A are well described in our analysis.

Figure 3 again summarizes our findings for the joint posteriors of the mass of the dynamical ejecta, the mass of the disk wind ejecta, the on-axis isotropic equivalent energy, the chirp mass, the mass ratio, the mass-weighted tidal deformability, and the radius of a 1.4 solar mass neutron star for this analysis, which is consistent with GW170817-and-AT2017gfo only. Compared to the analysis of GW170817-and-AT2017gfo only, the improvement on the parameter uncertainties is minimal, yet, noticeable when information from GRB170817A is added. Although no significant constraint on the EOS is imposed by the jet energy $E_0$ as the ratio $\xi$ between it, $E_0$, and the disk mass $m_{\text{disk}}$ is taken as a free parameter, the inclination constraint from the GRB plays a role in the constraint on EOS. For an anisotropic kilonova model, the inclination angle changes the observable kilonova light curves beyond scaling (e.g., Fig. 2 in ref. 59), which is correlated with the ejecta masses

(e.g., Fig. 3 in ref. 60). Therefore the GRB's inclination measurement imposes a constraint on the EOS via the kilonova eject masses measurements.

Moreover, for future studies, we expect that the inclusion of the GRB afterglow will be of great importance for measuring the Hubble constant.

## Discussion

We have developed a publicly available NMMA framework for the interpretation and analysis of BNS and BHNS systems. This framework allows for the simultaneous analysis of GW and EM signals such as kilonovae and GRB afterglows. In addition, our framework allows us to incorporate constraints from nuclear-physics calculations, e.g., by sampling over EOS sets constrained by chiral EFT, and to include radio as well as X-ray measurements of isolated NSs. By employing our framework to a combined analysis of GW170817, AT2017gfo, and GRB170817A, we find that the radius of a typical 1.4 solar mass NS lies within $11.98^{+0.35}_{-0.40}$ km; cf. Table 1 for a selection of studies from the literature. Based on our findings, our analysis is a noticeable improvement over previous works. However, additional uncertainties in our work lie in limited physics input in kilonova and semi-analytic GRB and models. Therefore, reliable astrophysical interpretations of future BNS detections will only be possible if not only parameter estimation infrastructure, as presented in this work, but also the astrophysical models describing transient phenomena advance further. Nevertheless, given the increasing number of multi-messenger detections of BNS and BHNS merger, we expect to use our framework to further increase our knowledge about the interior of NSs during the coming years.

## Methods

### Equation of state construction

The EOS describes the relation between energy density $\epsilon$, pressure $p$, and temperature $T$ of dense matter and additionally depends on the composition of the system. For NSs, thermal energies are much smaller than typical Fermi energies of the particles, and therefore, temperature effects can be neglected for isolated NSs or NSs in the inspiral phase of a merger. In these cases, the EOS simply relates $\epsilon$ and $p$.

The most general constraints on the EOS can be inferred from the slope of the EOS, the speed of sound, defined as:

$$c_S = c \sqrt{\frac{\partial p}{\partial \epsilon}}, \qquad (1)$$

where $c$ is the speed of light. Due to the laws of special relativity, the speed of sound has to be smaller than the speed of light, $c_S \le c$. Furthermore, the speed of sound in a NS has to be larger than zero, $c_S \ge 0$, as NSs would otherwise be unstable. These constraints alone, however, allow for an extremely large EOS space.

At nuclear densities, additional information on the EOS can be inferred from laboratory experiments and theoretical nuclear-physics calculations. For example, this information was used to constrain the properties of stellar matter in the NS crust⁶¹,⁶², i.e., the outermost layer of NSs at densities below approximately $0.5n_{\text{sat}}$. Above roughly $0.5n_{\text{sat}}$, NS matter consists of a fluid of neutrons with a small admixture of protons. In this regime, the EOS can be constrained by microscopic calculations of dense nuclear matter. These calculations typically provide the energy per particle, $E/A(n, x)$, which is a function of density $n$ and proton fraction $x = n_p/n$ with $n_p$ being the proton density. From this, the EOS follows from:

$$\varepsilon(n,x) = n \frac{E}{A}(n,x), \qquad (2)$$
<figure><img src="image_6.png" /><figcaption>**Fig. 3 | Visualization of the posterior of the GW170817-and-AT2017gfo and GW170817-and-AT2017gfo-and-GRB170817A analysis.** Corner plot for the mass of the dynamical ejecta $m_{dyn}^{ej}$, the mass of the disk wind ejecta $m_{wind}^{ej}$, $log_{10}$ of the GRB jet on-axis isotropic energy $log_{10}E_0$, the detector-frame chirp mass $M_c$, the mass ratio $q$, the mass-weighted tidal deformability $\Lambda$, and the radius of a 1.4 solar mass neutron star $R_{1.4}$ at 68%, 95% and 99% confidence. For the</figcaption></figure>

1D posterior probability distributions, we mark the median (solid lines) and the 90% confidence interval (dashed lines) and report these above each panel. We show results that are based on the simultaneous analysis of GW170817-and-AT2017gfo (orange) and of GW170817-and-AT2017gfo-and-GRB170817A (blue).

and

$$p(n,x) = n^2 \frac{\partial E/A(n,x)}{\partial n} \quad (3)$$

The proton fraction $x(n)$ is then determined from the beta equilibrium condition, $\mu_n = \mu_p + \mu_e$, where $\mu_i$ is the chemical potential of particle species i, and n, p, and e refer to neutrons, protons, and electrons, respectively.

To calculate the energy per particle microscopically, one needs to solve the nuclear many-body problem, commonly described by the Schrödinger equation. This requires knowledge of the nuclear Hamiltonian describing the many-body system. Fundamentally, nuclear many-body systems are described by Quantum Chromodynamics (QCD), the fundamental theory of strong nuclear interactions. QCD describes the system in terms of the fundamental degrees of freedom (d.o.f.), quarks and gluons. Unfortunately, this approach is currently not feasible⁶³. At densities of the order of $n_{sat}$, however, the effective d.o.f. are nucleons, neutrons and protons, that can be treated as point-
<figure><img src="image_7.png" /><figcaption>Fig. 3 | Visualization of the posterior of the GW170817-and-AT2017gfo and GW170817-and-AT2017gfo-and-GRB170817A analysis. Corner plot for the mass of the dynamical ejecta $m_{dyn}^{ej}$, the mass of the disk wind ejecta $m_{wind}^{ej}$, $log_{10}$ of the GRB jet on-axis isotropic energy $log_{10}E_0$, the detector-frame chirp mass $M_c$, the mass ratio $q$, the mass-weighted tidal deformability $\bar{\Lambda}$, and the radius of a 1.4 solar mass neutron star $R_{1.4}$ at 68%, 95% and 99% confidence. For the</figcaption></figure>

1D posterior probability distributions, we mark the median (solid lines) and the 90% confidence interval (dashed lines) and report these above each panel. We show results that are based on the simultaneous analysis of GW170817-and-AT2017gfo (orange) and of GW170817-and-AT2017gfo-and-GRB170817A (blue).

and

$$p(n,x) = n^2 \frac{\partial E/A(n,x)}{\partial n} \qquad (3)$$

The proton fraction $x(n)$ is then determined from the beta equilibrium condition, $\mu_n = \mu_p + \mu_e$, where $\mu_i$ is the chemical potential of particle species i, and n, p, and e refer to neutrons, protons, and electrons, respectively.

To calculate the energy per particle microscopically, one needs to solve the nuclear many-body problem, commonly described by the Schrödinger equation. This requires knowledge of the nuclear Hamiltonian describing the many-body system. Fundamentally, nuclear many-body systems are described by Quantum Chromodynamics (QCD), the fundamental theory of strong nuclear interactions. QCD describes the system in terms of the fundamental degrees of freedom (d.o.f.), quarks and gluons. Unfortunately, this approach is currently not feasible⁶³. At densities of the order of $n_{sat}$, however, the effective d.o.f. are nucleons, neutrons and protons, that can be treated as point-
<figure><img src="image_7.png" /><figcaption>**Fig. 4 | Best-fit late-time lightcurve from the analysis.** The best-fit lightcurves (dashed, with the 1 magnitude uncertainty shown as the band) for the analysis of GRB170817A when simultaneously analyzing GW170817, AT2017gf0, GRB170817A. We compare our model predictions with the observational data including the 1-sigma measurement uncertainty.</figcaption></figure>

<table><caption>**Table 1 | Comparison of radius measurements of a 1.4M<sub>☉</sub> neutron star for a selection of multi-messenger studies**</caption><thead><tr><th>Reference</th><th>R<sub>1.4</sub> [km]</th></tr></thead><tbody><tr><td>Dietrich et al.<sup>22</sup></td><td>11.75<sup>+0.86</sup><sub>-0.81</sub> (90%)</td></tr><tr><td>Essick et al.<sup>92</sup></td><td>12.54<sup>+0.71</sup><sub>-0.63</sub> (90%)</td></tr><tr><td>Breschi et al.<sup>30</sup></td><td>11.99<sup>+0.82</sup><sub>-0.85</sub> (90%)</td></tr><tr><td>Nicholl et al.<sup>31</sup></td><td>11.06<sup>+1.01</sup><sub>-0.98</sub> (90%)</td></tr><tr><td>Raaijmakers et al.<sup>148</sup></td><td>12.18<sup>+0.56</sup><sub>-0.79</sub> (95%)</td></tr><tr><td>Miller et al.<sup>96</sup></td><td>12.45<sup>+0.65</sup><sub>-0.65</sub> (68%)</td></tr><tr><td>Huth et al.<sup>23</sup></td><td>12.01<sup>+0.78</sup><sub>-0.77</sub> (90%)</td></tr><tr><td>this work [NMMA]</td><td>11.98<sup>+0.35</sup><sub>-0.40</sub> (90%)</td></tr></tbody></table>

A selected list of radius measurements of a 1.4M<sub>☉</sub> neutron star from various multi-messenger studies is shown. We denote the corresponding credible interval in parenthesis.

like nonrelativistic particles. Then, the nuclear Hamiltonian can be written generically as:

$$H = T + \sum_{i<j} V_{ij}^{NN} + \sum_{i<j<k} V_{ijk}^{3N} + \dots, \quad (4)$$

where *T* denotes the kinetic energy of the nucleons, $V_{ij}^{NN}$ describes two-nucleon (NN) interactions between nucleons *i* and *j*, and $V_{ijk}^{3N}$ describes three-nucleon (3N) interactions between nucleons *i*, *j*, and *k*. In principle, interactions involving four or more nucleons can be included, but initial studies have found these to be small compared to present uncertainties<sup>64</sup>.

The derivation of the nuclear Hamiltonian (Eq. (4)) from QCD is not feasible due to its nonperturbative nature. In this work, we therefore use a common approach and choose nucleons as effective d.o.f. The interactions among nucleons can then be derived in the framework of Chiral effective field theory (EFT)<sup>65,66</sup>. Chiral EFT starts out with the most general Lagrangian consistent with all the symmetries of QCD in terms of nucleonic degrees of freedom. It explicitly includes meson-exchange interactions for the lightest mesons, i.e., the pions. This approach yields an infinite number of pion-exchange and nucleon-contact interactions which needs to be organized in terms of a hierarchical expansion in powers of a soft (low-energy) scale over a hard (high-energy) scale. In chiral EFT, the soft scale *q* is given by the nucleons' external momenta or the pion mass. The hard scale, also called the breakdown scale $\Lambda_b$, is of the order of 500–600 MeV<sup>67</sup> and interaction contributions involving heavier d.o.f., such as the $\rho$ meson, are integrated out. The chiral Lagrangian is then expanded in powers

of $q/\Lambda_b$ according to a power-counting scheme. Most current chiral EFT interactions are derived in Weinberg power counting<sup>65,66,68–70</sup>. One can then derive the nuclear Hamiltonian from this chiral Lagrangian in a consistent order-by-order framework that allows for an estimate of the theoretical uncertainties<sup>67,71,72</sup> and that can be systematically improved by increasing the order of the calculation. Chiral EFT Hamiltonian naturally includes NN, 3N, and higher many-body forces, see Eq. (4), and chiral EFT predicts a natural hierarchy of these contributions. For example, 3N interactions start to contribute at third order ($N^2$LO) in the expansion. Typical state-of-the art calculations truncate the chiral expansion at $N^2$LO<sup>37,40,73</sup> or fourth order ($N^3$LO)<sup>39,74</sup>.

With the nuclear Hamiltonian at hand, one then needs to solve the many-body Schrödinger equation which requires advanced numerical methods. Examples of such many-body techniques include many-body perturbation theory (MBPT)<sup>38,39,74</sup>, the self-consistent Green's function (SCGF) method<sup>75</sup>, or the coupled-cluster (CC) method<sup>73,76</sup>. Here, we employ Quantum Monte Carlo (QMC) methods<sup>77</sup>, which provide non-perturbative solutions of the Schrödinger equation. QMC methods are stochastic techniques which treat the Schrödinger equation as a diffusion equation in imaginary time. In the QMC framework, one begins by choosing a trial wavefunction of the many-body system, which for nuclear matter can be described as a Slater determinant of non-interacting fermions multiplied with NN and 3N correlation functions. This trial wavefunction is evolved to large imaginary times, projecting out high-energy excitations, and converging to the true ground state of the system as long as the trial wavefunction has a non-zero overlap with it. Among QMC methods, two well-established algorithms are Green's function Monte Carlo (GFMC), used to describe light atomic nuclei with great precision<sup>77</sup>, and Auxiliary Field Diffusion Monte Carlo (AFDMC)<sup>78</sup>, suitable to study larger systems such as nuclear matter. Here, we employ AFDMC calculations of neutron matter but our NMMA framework is sufficiently flexible to employ any low-density calculation for neutron-star matter. We then extend our neutron-matter calculations to neutron-star conditions by extrapolating the calculations to $\beta$ equilibrium using phenomenological information on symmetric nuclear matter and constructing a consistent crust reflecting the uncertainties of the calculations<sup>79</sup>. This crust includes a description of the outer crust<sup>61</sup> and uses the Wigner-Seitz approximation to calculate the inner-crust EOS consistently with our AFDMC calculations.

At nuclear densities, chiral EFT together with a suitable many-body framework provides for a reliable description of nuclear matter with systematic uncertainty estimates. With increasing density, however, the associated theoretical uncertainty grows fast due to the correspondingly larger nucleon momenta approaching the breakdown scale. The density up to which chiral EFT remains valid is not exactly known but estimates place it around $2n_{sat}^{67,80}$. Hence, chiral EFT calculations constrain the EOS only up to these densities but to explore the large EOS space beyond the breakdown of chiral EFT, one requires a physics-agnostic extension scheme. Here, physics-agnostic implies that no model assumptions, e.g., about the existence of certain d.o.f. at high densities, are made. Instead, the EOS is only bounded by conditions of causality, $c_S \le c$, and mechanical stability, $c_S \ge 0$, mentioned before. There exist several such extension schemes in literature: parametric ones, like the polytropic expansion<sup>81–83</sup> or expansions in the speed of sound<sup>84,85</sup>, and non-parametric approaches<sup>86</sup>. To extend the AFDMC calculations employed here, we employ a parametric speed-of-sound extension scheme. Working in the $c_S$ versus *n* plane, the speed of sound $c_S(n)$ is determined with theoretical uncertainty estimates by chiral EFT up to a reference density below the expected breakdown density. From this uncertainty band, we sample a speed-of-sound curve up to the reference density. Beyond this density, we create a typically non-uniform grid in density up to a large density $\approx 12n_{sat}$, well beyond the regime realized in NSs. For each grid point, we sample random values
<figure><img src="image_8.png" /><figcaption>**Fig. 4 | Best-fit late-time lightcurve from the analysis.** The best-fit lightcurves (dashed, with the 1 magnitude uncertainty shown as the band) for the analysis of GRB170817A when simultaneously analyzing GW170817, AT2017gf0, GRB170817A. We compare our model predictions with the observational data including the 1-sigma measurement uncertainty.</figcaption></figure>

<table><caption>**Table 1 | Comparison of radius measurements of a 1.4M<sub>☉</sub> neutron star for a selection of multi-messenger studies**</caption><thead><tr><th>Reference</th><th>R<sub>1.4</sub> [km]</th></tr></thead><tbody><tr><td>Dietrich et al.<sup>22</sup></td><td>11.75<sup>+0.86</sup><sub>-0.81</sub> (90%)</td></tr><tr><td>Essick et al.<sup>92</sup></td><td>12.54<sup>+0.71</sup><sub>-0.63</sub> (90%)</td></tr><tr><td>Breschi et al.<sup>30</sup></td><td>11.99<sup>+0.82</sup><sub>-0.85</sub> (90%)</td></tr><tr><td>Nicholl et al.<sup>31</sup></td><td>11.06<sup>+1.01</sup><sub>-0.98</sub> (90%)</td></tr><tr><td>Raaijmakers et al.<sup>148</sup></td><td>12.18<sup>+0.56</sup><sub>-0.79</sub> (95%)</td></tr><tr><td>Miller et al.<sup>96</sup></td><td>12.45<sup>+0.65</sup><sub>-0.65</sub> (68%)</td></tr><tr><td>Huth et al.<sup>23</sup></td><td>12.01<sup>+0.78</sup><sub>-0.77</sub> (90%)</td></tr><tr><td>this work [NMMA]</td><td>11.98<sup>+0.35</sup><sub>-0.40</sub> (90%)</td></tr></tbody></table>

A selected list of radius measurements of a 1.4M<sub>☉</sub> neutron star from various multi-messenger studies is shown. We denote the corresponding credible interval in parenthesis.

like nonrelativistic particles. Then, the nuclear Hamiltonian can be written generically as:

$$H = T + \sum_{i<j} V_{ij}^{NN} + \sum_{i<j<k} V_{ijk}^{3N} + \dots, \quad (4)$$

where *T* denotes the kinetic energy of the nucleons, $V_{ij}^{NN}$ describes two-nucleon (NN) interactions between nucleons *i* and *j*, and $V_{ijk}^{3N}$ describes three-nucleon (3N) interactions between nucleons *i*, *j*, and *k*. In principle, interactions involving four or more nucleons can be included, but initial studies have found these to be small compared to present uncertainties⁶⁴.

The derivation of the nuclear Hamiltonian (Eq. (4)) from QCD is not feasible due to its nonperturbative nature. In this work, we therefore use a common approach and choose nucleons as effective d.o.f. The interactions among nucleons can then be derived in the framework of Chiral effective field theory (EFT)⁶⁵,⁶⁶. Chiral EFT starts out with the most general Lagrangian consistent with all the symmetries of QCD in terms of nucleonic degrees of freedom. It explicitly includes meson-exchange interactions for the lightest mesons, i.e., the pions. This approach yields an infinite number of pion-exchange and nucleon-contact interactions which needs to be organized in terms of a hierarchical expansion in powers of a soft (low-energy) scale over a hard (high-energy) scale. In chiral EFT, the soft scale *q* is given by the nucleons' external momenta or the pion mass. The hard scale, also called the breakdown scale $\Lambda_b$, is of the order of 500–600 MeV⁶⁷ and interaction contributions involving heavier d.o.f., such as the $\rho$ meson, are integrated out. The chiral Lagrangian is then expanded in powers

of $q/\Lambda_b$ according to a power-counting scheme. Most current chiral EFT interactions are derived in Weinberg power counting⁶⁵,⁶⁶,⁶⁸⁻⁷⁰. One can then derive the nuclear Hamiltonian from this chiral Lagrangian in a consistent order-by-order framework that allows for an estimate of the theoretical uncertainties⁶⁷,⁷¹,⁷² and that can be systematically improved by increasing the order of the calculation. Chiral EFT Hamiltonian naturally includes NN, 3N, and higher many-body forces, see Eq. (4), and chiral EFT predicts a natural hierarchy of these contributions. For example, 3N interactions start to contribute at third order ($N^2$LO) in the expansion. Typical state-of-the art calculations truncate the chiral expansion at $N^2$LO³⁷,⁴⁰,⁷³ or fourth order ($N^3$LO)³⁹,⁷⁴.

With the nuclear Hamiltonian at hand, one then needs to solve the many-body Schrödinger equation which requires advanced numerical methods. Examples of such many-body techniques include many-body perturbation theory (MBPT)³⁸,³⁹,⁷⁴, the self-consistent Green's function (SCGF) method⁷⁵, or the coupled-cluster (CC) method⁷³,⁷⁶. Here, we employ Quantum Monte Carlo (QMC) methods⁷⁷, which provide non-perturbative solutions of the Schrödinger equation. QMC methods are stochastic techniques which treat the Schrödinger equation as a diffusion equation in imaginary time. In the QMC framework, one begins by choosing a trial wavefunction of the many-body system, which for nuclear matter can be described as a Slater determinant of non-interacting fermions multiplied with NN and 3N correlation functions. This trial wavefunction is evolved to large imaginary times, projecting out high-energy excitations, and converging to the true ground state of the system as long as the trial wavefunction has a non-zero overlap with it. Among QMC methods, two well-established algorithms are Green's function Monte Carlo (GFMC), used to describe light atomic nuclei with great precision⁷⁷, and Auxiliary Field Diffusion Monte Carlo (AFDMC)⁷⁸, suitable to study larger systems such as nuclear matter. Here, we employ AFDMC calculations of neutron matter but our NMMA framework is sufficiently flexible to employ any low-density calculation for neutron-star matter. We then extend our neutron-matter calculations to neutron-star conditions by extrapolating the calculations to $\beta$ equilibrium using phenomenological information on symmetric nuclear matter and constructing a consistent crust reflecting the uncertainties of the calculations⁷⁹. This crust includes a description of the outer crust⁶¹ and uses the Wigner-Seitz approximation to calculate the inner-crust EOS consistently with our AFDMC calculations.

At nuclear densities, chiral EFT together with a suitable many-body framework provides for a reliable description of nuclear matter with systematic uncertainty estimates. With increasing density, however, the associated theoretical uncertainty grows fast due to the correspondingly larger nucleon momenta approaching the breakdown scale. The density up to which chiral EFT remains valid is not exactly known but estimates place it around $2n_{sat}^{67,80}$. Hence, chiral EFT calculations constrain the EOS only up to these densities but to explore the large EOS space beyond the breakdown of chiral EFT, one requires a physics-agnostic extension scheme. Here, physics-agnostic implies that no model assumptions, e.g., about the existence of certain d.o.f. at high densities, are made. Instead, the EOS is only bounded by conditions of causality, $c_S \le c$, and mechanical stability, $c_S \ge 0$, mentioned before. There exist several such extension schemes in literature: parametric ones, like the polytropic expansion⁸¹⁻⁸³ or expansions in the speed of sound⁸⁴,⁸⁵, and non-parametric approaches⁸⁶. To extend the AFDMC calculations employed here, we employ a parametric speed-of-sound extension scheme. Working in the $c_S$ versus *n* plane, the speed of sound $c_S(n)$ is determined with theoretical uncertainty estimates by chiral EFT up to a reference density below the expected breakdown density. From this uncertainty band, we sample a speed-of-sound curve up to the reference density. Beyond this density, we create a typically non-uniform grid in density up to a large density ≈ $12n_{sat}$, well beyond the regime realized in NSs. For each grid point, we sample random values
for $c_s^2(n_i)$ between 0 and $c^2$ (we set $c=1$ in the following). We then connect the chiral EFT draw for the speed of sound with all points $c_{s,i}^2(n_i)$ using linear segments. The resulting density-dependent speed of sound can be integrated to give the EOS, i.e., the pressure, baryon density, and energy density. In the interval $n_i \le n \le n_{i+1}$:

$$p(n) = p(n_i) + \int_{n_i}^{n} c_s^2(n')\mu(n')dn', \quad (5)$$

$$c(n) = c(n_i) + \int_{n_i}^{n} \mu(n')dn', \quad (6)$$

where $\mu(n)$ is the chemical potential that can be obtained from the speed of sound using the relation:

$$\mu(n) = \mu_i \exp \left[ \int_{\log n_i}^{\log n} c_s^2(\log n') d \log n' \right]. \quad (7)$$

For each reconstructed EOS, constrained by Chiral EFT at low densities and extrapolated via the $c_s$ extension to larger densities, the global properties of NSs can be calculated by solving the Tolman-Oppenheimer-Volkoff (TOV) equations. This way, we determine the NS radii ($R$) and dimensionless tidal deformabilities ($\Lambda$) as functions of their masses ($M$). We repeat this approach for a large number of samples to construct EOS priors for further analyses of NS data.

This approach is flexible and additional information on high-density phases of QCD can be included straightforwardly. For example, pQCD calculations at asymptotically high densities$^{47}$, of the order of 40–50$n_{sat}$, might be used to constrain the general EOS extension schemes even further$^{48,83}$. However, the exact impact of these constraints at densities well beyond the regime realized in NSs needs to be studied in more detail. While our NMMA framework currently does not have this capability, we are planning to add this in the near future. Similarly, instead of using general extension models, one can employ specific high-density models accounting for quark and gluon d.o.f. One such model is the quarkyonic-matter model$^{87-90}$, which describes the observed behavior of the speed of sound in NSs$^{80}$: a rise of the speed of sound at low densities to values above the conformal limit of $c/\sqrt{3}$, followed by a decrease to values below the conformal limit at higher densities. In future work, we will address quarkyonic matter and other models in our NMMA framework.

The construction of the EOS, as detailed above, is implemented in the NMMA code under the class `EOS_with_CSE`. This class allows for (1) an exploration of theoretical uncertainties in the low-density EOS and (2) constructs the high-density EOS using a $c_s$ extrapolation. (1) Low-density uncertainties are implemented by requiring two tabulated EOS files for the lower and upper bound of the uncertainty band as inputs, containing the pressure, energy density and number density up to the chosen breakdown density of the model. By default, the results of a QMC calculation using local chiral EFT interactions at $N^2LO^{80}$ with theoretical uncertainties are provided. Upon initiation of the class, a sample is drawn from the low-density uncertainty band using a 1-parameter sampling technique. In this approach, a uniform random number $\omega$ is sampled uniformly between 0 and 1, and the interpolated EOS is given as:

$$p(n) = p_{\text{soft}}(n) + \omega(p_{\text{stiff}}(n) - p_{\text{soft}}(n)), \quad (8)$$

$$\varepsilon(n) = \varepsilon_{\text{soft}}(n) + \omega(\varepsilon_{\text{stiff}}(n) - \varepsilon_{\text{soft}}(n)), \quad (9)$$

where the subscripts “soft” and “stiff” refer to the lower and upper bounds of the EFT uncertainty band, respectively. This sampling technique assumes that pressure and energy density are correlated but

we have found that releasing this assumption and using a four-parameter form suggested by Gandolfi et al.$^{91}$ does not change our results appreciably. In future, we will explore additional schemes, e.g., using Gaussian processes$^{92}$.

(2) The EOS given by Eqs. (8) and (9) is used up to a breakdown density determined by the user. By default, this density is set to $2n_{sat}$. Beyond this density, the class constructs the EOS using a $c_s$ extension. The maximum density up to which the EOS is extrapolated and the number of linear line segments can be adjusted by the user, with the default values being $12n_{sat}$ for the former and 5 line segments for the latter. The code then solves Eqs. (5)–(7) to give the extrapolated EOS. The pressure, energy density, and number density describing the full EOS are accessible as attributes of the `EOS_with_CSE` class.

Finally, the method `construct_family` solves the stellar structure equations (TOV equations and equations for the quadrupole perturbation of spherical models), and returns a sequence of NSs with their masses, radii and dimensionless tidal deformabilities as arrays.

## Prior weighting to incorporate radio and X-ray observations of single neutron stars

To incorporate mass measurements of heavy pulsars and mass-radius measurements of isolated pulsars, the associated likelihood is calculated and taken as the prior probability for an EOS for further analysis. For instance, the radio observations on PSR J0348+4042$^{44}$, and PSR J1614-2230$^{45}$ provide a lower bound on the maximum mass of a NS.

The likelihood for a mass-only measurement is given by:

$$\mathcal{L}_{\text{PSR-mass}}(\mathbf{E}) = \int_{0}^{M_{\text{TOV}}} dM \mathcal{P}(M|\text{PSR}), \quad (10)$$

where $\mathcal{P}(M|\text{PSR})$ is the posterior distribution of the pulsar's mass and $M_{TOV}$ is the maximum mass supported by the EOS with parameters $\mathbf{E}$. The posterior distributions of pulsar masses are typically well approximated by Gaussians$^{22}$.

Recent X-ray observations of millisecond pulsars by NASA's Neutron Star Interior Composition Explorer (NICER) mission have been used to simultaneously determine the mass and radius of these NSs$^{93-97}$. The corresponding likelihood is given by:

$$\begin{align}
\mathcal{L}_{\text{NICER}}(\mathbf{E}) &= \int dM \int dR \mathcal{P}_{\text{NICER}}(M, R) \frac{\pi(M, R|\mathbf{E})}{\pi(M, R|I)} \nonumber \\
&\propto \int dM \int dR \mathcal{P}_{\text{NICER}}(M, R) \delta(R - R(M; \mathbf{E})) \tag{11} \\
&\propto \int dM \mathcal{P}_{\text{NICER}}(M, R = R(M; \mathbf{E})), \nonumber
\end{align}$$

where $\mathcal{P}_{\text{NICER}}(M, R)$ is the joint-posterior distribution of mass and radius as measured by NICER and we use the fact that (1) the radius is a function of mass for a given EOS, and (2) that without further EOS information, e.g., through chiral EFT, the prior for the radius given mass is taken to be uniform.

## Gravitational-wave inference

**GW models.** A complex frequency-domain GW signal is given by:

$$h(f) = A(f)e^{-i\psi(f)}, \quad (12)$$

with the amplitude $A(f)$ and the GW phase $\psi(f)$. Because of the NS's finite size and internal structure, BNS and BHNS waveform models have to incorporate tidal contributions for an accurate interpretation of the binary coalescence. Such tidal contributions account for the deformation of the stars in their companions' external gravitational field$^{98,99}$ and, once measured, allow to place constraints on the EOS governing the NS interior$^{100-103}$. They are attractive because they convert energy from the orbital motion to a deformation of the stars, and
for $c_s^2(n_i)$ between 0 and $c^2$ (we set $c=1$ in the following). We then connect the chiral EFT draw for the speed of sound with all points $c_{s,i}^2(n_i)$ using linear segments. The resulting density-dependent speed of sound can be integrated to give the EOS, i.e., the pressure, baryon density, and energy density. In the interval $n_i \le n \le n_{i+1}$:

$$p(n) = p(n_i) + \int_{n_i}^{n} c_s^2(n')\mu(n')dn', \quad (5)$$

$$c(n) = c(n_i) + \int_{n_i}^{n} \mu(n')dn', \quad (6)$$

where $\mu(n)$ is the chemical potential that can be obtained from the speed of sound using the relation:

$$\mu(n) = \mu_i \exp \left[ \int_{\log n_i}^{\log n} c_s^2(\log n') d \log n' \right]. \quad (7)$$

For each reconstructed EOS, constrained by Chiral EFT at low densities and extrapolated via the $c_s$ extension to larger densities, the global properties of NSs can be calculated by solving the Tolman-Oppenheimer-Volkoff (TOV) equations. This way, we determine the NS radii ($R$) and dimensionless tidal deformabilities ($\Lambda$) as functions of their masses ($M$). We repeat this approach for a large number of samples to construct EOS priors for further analyses of NS data.

This approach is flexible and additional information on high-density phases of QCD can be included straightforwardly. For example, pQCD calculations at asymptotically high densities$^{47}$, of the order of 40–50$n_{sat}$, might be used to constrain the general EOS extension schemes even further$^{48,83}$. However, the exact impact of these constraints at densities well beyond the regime realized in NSs needs to be studied in more detail. While our NMMA framework currently does not have this capability, we are planning to add this in the near future. Similarly, instead of using general extension models, one can employ specific high-density models accounting for quark and gluon d.o.f. One such model is the quarkyonic-matter model$^{87-90}$, which describes the observed behavior of the speed of sound in NSs$^{80}$: a rise of the speed of sound at low densities to values above the conformal limit of $c/\sqrt{3}$, followed by a decrease to values below the conformal limit at higher densities. In future work, we will address quarkyonic matter and other models in our NMMA framework.

The construction of the EOS, as detailed above, is implemented in the NMMA code under the class `EOS_with_CSE`. This class allows for (1) an exploration of theoretical uncertainties in the low-density EOS and (2) constructs the high-density EOS using a $c_s$ extrapolation. (1) Low-density uncertainties are implemented by requiring two tabulated EOS files for the lower and upper bound of the uncertainty band as inputs, containing the pressure, energy density and number density up to the chosen breakdown density of the model. By default, the results of a QMC calculation using local chiral EFT interactions at $N^2LO^{80}$ with theoretical uncertainties are provided. Upon initiation of the class, a sample is drawn from the low-density uncertainty band using a 1-parameter sampling technique. In this approach, a uniform random number $\omega$ is sampled uniformly between 0 and 1, and the interpolated EOS is given as:

$$p(n) = p_{\text{soft}}(n) + \omega(p_{\text{stiff}}(n) - p_{\text{soft}}(n)), \quad (8)$$

$$\varepsilon(n) = \varepsilon_{\text{soft}}(n) + \omega(\varepsilon_{\text{stiff}}(n) - \varepsilon_{\text{soft}}(n)), \quad (9)$$

where the subscripts “soft” and “stiff” refer to the lower and upper bounds of the EFT uncertainty band, respectively. This sampling technique assumes that pressure and energy density are correlated but

we have found that releasing this assumption and using a four-parameter form suggested by Gandolfi et al.$^{91}$ does not change our results appreciably. In future, we will explore additional schemes, e.g., using Gaussian processes$^{92}$.

(2) The EOS given by Eqs. (8) and (9) is used up to a breakdown density determined by the user. By default, this density is set to $2n_{sat}$. Beyond this density, the class constructs the EOS using a $c_s$ extension. The maximum density up to which the EOS is extrapolated and the number of linear line segments can be adjusted by the user, with the default values being $12n_{sat}$ for the former and 5 line segments for the latter. The code then solves Eqs. (5)–(7) to give the extrapolated EOS. The pressure, energy density, and number density describing the full EOS are accessible as attributes of the `EOS_with_CSE` class.

Finally, the method `construct_family` solves the stellar structure equations (TOV equations and equations for the quadrupole perturbation of spherical models), and returns a sequence of NSs with their masses, radii and dimensionless tidal deformabilities as arrays.

## Prior weighting to incorporate radio and X-ray observations of single neutron stars

To incorporate mass measurements of heavy pulsars and mass-radius measurements of isolated pulsars, the associated likelihood is calculated and taken as the prior probability for an EOS for further analysis. For instance, the radio observations on PSR J0348+4042$^{44}$, and PSR J1614-2230$^{45}$ provide a lower bound on the maximum mass of a NS.

The likelihood for a mass-only measurement is given by:

$$\mathcal{L}_{\text{PSR-mass}}(\mathbf{E}) = \int_{0}^{M_{\text{TOV}}} dM \mathcal{P}(M|\text{PSR}), \quad (10)$$

where $\mathcal{P}(M|\text{PSR})$ is the posterior distribution of the pulsar’s mass and $M_{\text{TOV}}$ is the maximum mass supported by the EOS with parameters $\mathbf{E}$. The posterior distributions of pulsar masses are typically well approximated by Gaussians$^{22}$.

Recent X-ray observations of millisecond pulsars by NASA’s Neutron Star Interior Composition Explorer (NICER) mission have been used to simultaneously determine the mass and radius of these NSs$^{93–97}$. The corresponding likelihood is given by:

$$\begin{align}
\mathcal{L}_{\text{NICER}}(\mathbf{E}) &= \int dM \int dR \mathcal{P}_{\text{NICER}}(M, R) \frac{\pi(M, R|\mathbf{E})}{\pi(M, R|I)} \nonumber \\
&\propto \int dM \int dR \mathcal{P}_{\text{NICER}}(M, R) \delta(R - R(M; \mathbf{E})) \tag{11} \\
&\propto \int dM \mathcal{P}_{\text{NICER}}(M, R = R(M; \mathbf{E})), \nonumber
\end{align}$$

where $\mathcal{P}_{\text{NICER}}(M, R)$ is the joint-posterior distribution of mass and radius as measured by NICER and we use the fact that (1) the radius is a function of mass for a given EOS, and (2) that without further EOS information, e.g., through chiral EFT, the prior for the radius given mass is taken to be uniform.

## Gravitational-wave inference

**GW models.** A complex frequency-domain GW signal is given by:

$$h(f) = A(f)e^{-i\psi(f)}, \quad (12)$$

with the amplitude $A(f)$ and the GW phase $\psi(f)$. Because of the NS’s finite size and internal structure, BNS and BHNS waveform models have to incorporate tidal contributions for an accurate interpretation of the binary coalescence. Such tidal contributions account for the deformation of the stars in their companions’ external gravitational field$^{98,99}$ and, once measured, allow to place constraints on the EOS governing the NS interior$^{100–103}$. They are attractive because they convert energy from the orbital motion to a deformation of the stars, and
lead to an accelerated inspiral. In the case of non-spinning compact objects, the leading-order tidal contribution depends on the tidal deformability:

$$ \bar{\Lambda} = \frac{16(m_1 + 12m_2)m_1^4 \Lambda_1 + (m_2 + 12m_1)m_2^4 \Lambda_2}{13(m_1 + m_2)^5} \quad (13) $$

with the individual tidal deformabilities $\Lambda_{1,2} = \frac{2}{3}K_2^{1,2}/C_{1,2}^5$ and the individual masses $m_{1,2}$. Here, $K_2^{1,2}$ are the Love numbers describing the static quadrupole deformation of one body inside the gravitoelectric field of the companion and $C_{1,2}$ are the individual compactnesses $C_{1,2} = m_{1,2}/R_{1,2}$ in isolation.

To date, there are three different types of BNS or BHNS models for the inspiral GW signal that are commonly used: Post-Newtonian (PN) models$^{104-107}$, effective-one-body (EOB) models$^{108-116}$, and phenomenological approximants$^{117-121}$. In the NMMA framework, we make use of the LALSuite$^{122}$ software package, in particular LALSimulation, so that the BNS and BHNS models used by the LIGO-Virgo-Kagra Collaborations can be easily employed. This includes:

* PN models such as TaylorT2, TaylorT4, or TaylorF2 where a PN descriptions for the point-particle BBH baseline as well as the tidal description is employed.

* the most commonly used tidal EOB models SEOBNRv4T$^{111,116,123}$, its frequency-domain surrogate model$^{124}$, as well as the TEO-BResumS model$^{112,125}$ including its post-adiabatic accelerated version$^{126}$ which enables it being used during parameter estimation.

* and phenomenological models such as IMRPhenomD_NRTidal, SEOBNRv4_ROM_NRTidal, IMRPhenomPv2_NRTidal, IMRPhenomD_NRTidalv2, SEOBNRv4_ROM_NRTidalv2, IMRPhenomPv2_NRTidalv2$^{118-120}$, PhenomNSBH, and SEOBNRv4_ROM_NRTidalv2_NSBH$^{121,127}$.

**GW analysis.** By assuming stationary Gaussian noise, the GW likelihood $\mathcal{L}_{GW}(\theta)$ that the data $d$ is a sum of noise and a GW signal $h$ with parameters $\theta$ is given by$^{128}$:

$$ \mathcal{L}_{GW} \propto \exp\left(-\frac{1}{2}\langle d - h(\theta) | d - h(\theta) \rangle\right), \quad (14) $$

where the inner product $\langle a | b \rangle$ is defined as:

$$ \langle a | b \rangle = 4\pi \int_{f_{\text{low}}}^{f_{\text{high}}} \frac{\tilde{a}(f) \tilde{b}^{*}(f)}{S_n(f)} df. \quad (15) $$

Here, $\tilde{a}(f)$ is the Fourier transform of $a(t)$, denotes complex conjugation, and $S_n(f)$ is the one-sided power spectral density of the noise. The choice of $f_{\text{low}}$ and $f_{\text{high}}$ depends on the type of binary that we are interested in. In our study, we will set $f_{\text{low}}$ and $f_{\text{high}}$ to 20 Hz and 2048 Hz, respectively. This is sufficient for capturing the inspiral up to the moment of merger for a typical BNS system in the advanced GW detector era.

## Electromagnetic signals

**Kilonova models.** Kilonova models are extracted using the 3D Monte Carlo radiative transfer code POSSIS$^{129}$. The code can handle arbitrary geometries for the ejected material and produces spectra, lightcurves and polarization as a function of the observer viewing angle. Given an input model with defined densities $\rho$ and compositions (i.e., electron fraction $Y_e$), the code generates Monte Carlo photon packets with initial location and energy sampled from the energy distribution from radioactive decay of r-process nuclei within the model. The latter depends on the mass/density distribution of the model and the assumed nuclear heating rates and thermalization efficiencies. The frequency of each Monte Carlo photon packet is sampled according to the temperature $T$ in the ejecta, which is calculated at each time-step$^{130,131}$. Photon packets

are then followed as they diffuse out of the ejected material and interact with matter via either electron scattering or bound-bound line transitions. Time- and wavelength-dependent opacities $\kappa_{\Lambda}(\rho, T, Y_e, t)$ from Tanaka et al.$^{132}$ are implemented in the code and depend on the local properties of the ejecta ($\rho, T$, and $Y_e$). Spectral time series are extracted using the technique described by Bulla et al.$^{133}$ and used to construct broad-band lightcurves in any desired filter.

**Supernova models.** Templates available within the SNCosmo library$^{134}$ are used to model supernova spectra. Currently, the salt2 model for Type Ia supernovae and the nugent-hyper model for hypernovae associated with long GRBs are implemented in the framework and have been used in the past$^{36}$. However, the framework is flexible enough such that additional templates for different types of supernovae can be added with minimal effort.

**Kilonova/supernova inference.** Our EM inference of kilonovae and GRB afterglows is based on the AB magnitude for a specific filter $j$, $m_i^j(t_i)$. We assume these measurements to be given as a time series at times $t_i$ with a corresponding statistical error $\sigma_i^j = \sigma^j(t_i)$. The likelihood function $\mathcal{L}_{\text{EM}}(\theta)$ then reads$^{135}$:

$$ \mathcal{L}_{\text{EM}} \propto \exp \left( -\frac{1}{2} \sum_{ij} \frac{(m_i^j - m_i^{\text{est}}(\theta))^2}{(\sigma_i^j)^2 + \sigma_{\text{sys}}^2} \right), \quad (16) $$

where $m_i^{\text{est}}(\theta)$ is the estimated AB magnitude for the parameters $\theta$ and $\sigma_{\text{sys}}$ is the additional error budget for accounting the systematic uncertainty within the electromagnetic signal modeling. The inclusion of $\sigma_{\text{sys}}$ is equivalent to adding a shift of $\Delta m$ to the light curve, for which marginalized with respect to a zero-mean normal distribution with a variance of $\sigma_{\text{sys}}^2$.

This likelihood is equivalent to approximating the probability distribution of the spectral flux density $f_v$ to be a Log-normal distribution. The Log-normal distribution is a 2-parameter maximum entropy distribution with its support equals to the possible range for $f_v \in (0, \infty)$. There are two advantages of approximating $f_v$ with a Log-normal distribution: (1) if the uncertainty is larger or comparable to the measured value, it avoids having non-zero support for the nonphysical $f_v < 0$; (2) if the uncertainty is much smaller than the measured value, the Log-normal distribution approaches the normal distribution.

For kilonovae, we use the same model presented in Dietrich et al.$^{22}$. The model is controlled by four parameters, namely, the dynamical ejecta mass $m_{\text{dyn}}^{\text{dyn}}$, the disk wind ejecta mass $m_{\text{wind}}^{\text{dyn}}$, the half-opening angle of the lanthanide-rich component $\Phi$, and the viewing angle $\theta_{\text{obs}}$.

**GRB afterglows.** In our framework, the computation of the GRB afterglow lightcurves is until now based on the publicly available semi-analytic code afterglowwp$^{57,58}$. The inclusion of other afterglow models is currently ongoing.

The GRB afterglow emission is produced by relativistic electrons gyrating around the magnetic field lines. These electrons are accelerated by the Fermi first-order acceleration (diffusive shock acceleration) and the magnetic field is assumed to be of turbulent nature, amplified by processes acting in collision-less shocks. The complex physics of electron acceleration at shocks is approximated by the equipartition parameters, $\epsilon_e$ and $\epsilon_B$, denoting the fraction of the shock energy that goes into the relativistic electrons and magnetic field, respectively, and $\rho$, and the slope of the electron energy distribution $dn/d\rho \propto Y^p$, with $n$ being the electron number density and $\gamma$ being the electron Lorentz factor. The flux density of the curvature radiation is:

$$ F_{\nu} = \frac{1}{4\pi d_{L}^{2}} \int d\theta d\phi R^{2} \sin(\theta) \frac{\epsilon_{\nu}}{\alpha_{\nu}} (1 - e^{-\tau}), \quad (17) $$
lead to an accelerated inspiral. In the case of non-spinning compact objects, the leading-order tidal contribution depends on the tidal deformability:

$$ \bar{\Lambda} = \frac{16(m_1 + 12m_2)m_1^4 \Lambda_1 + (m_2 + 12m_1)m_2^4 \Lambda_2}{13(m_1 + m_2)^5} \quad (13) $$

with the individual tidal deformabilities $\Lambda_{1,2} = \frac{2}{3}K_2^{1,2}/C_{1,2}^5$ and the individual masses $m_{1,2}$. Here, $K_2^{1,2}$ are the Love numbers describing the static quadrupole deformation of one body inside the gravitoelectric field of the companion and $C_{1,2}$ are the individual compactnesses $C_{1,2} = m_{1,2}/R_{1,2}$ in isolation.

To date, there are three different types of BNS or BHNS models for the inspiral GW signal that are commonly used: Post-Newtonian (PN) models$^{104-107}$, effective-one-body (EOB) models$^{108-116}$, and phenomenological approximants$^{117-121}$. In the NMMA framework, we make use of the LALSuite$^{122}$ software package, in particular LALSimulation, so that the BNS and BHNS models used by the LIGO-Virgo-Kagra Collaborations can be easily employed. This includes:

* PN models such as TaylorT2, TaylorT4, or TaylorF2 where a PN descriptions for the point-particle BBH baseline as well as the tidal description is employed.

* the most commonly used tidal EOB models SEOBNRv4$^{111,116,123}$, its frequency-domain surrogate model$^{124}$, as well as the TEO-BResumS model$^{112,125}$ including its post-adiabatic accelerated version$^{126}$ which enables it being used during parameter estimation.

* and phenomenological models such as IMRPhenomD_NRTidal, SEOBNRv4_ROM_NRTidal, IMRPhenomPv2_NRTidal, IMRPhenomD_NRTidalv2, SEOBNRv4_ROM_NRTidalv2, IMRPhenomPv2_NRTidalv2$^{118-120}$, PhenomNSBH, and SEOBNRv4_ROM_NRTidalv2_NSBH$^{121,127}$.

**GW analysis.** By assuming stationary Gaussian noise, the GW likelihood $\mathcal{L}_{GW}(\boldsymbol{\theta})$ that the data $d$ is a sum of noise and a GW signal $h$ with parameters $\boldsymbol{\theta}$ is given by$^{128}$:

$$ \mathcal{L}_{GW} \propto \exp\left(-\frac{1}{2}\langle d - h(\boldsymbol{\theta})|d - h(\boldsymbol{\theta})\rangle\right), \quad (14) $$

where the inner product $\langle a|b\rangle$ is defined as:

$$ \langle a|b\rangle = 4\pi \int_{f_{\text{low}}}^{f_{\text{high}}} \frac{\tilde{a}(f)\tilde{b}^*(f)}{S_n(f)} df. \quad (15) $$

Here, $\tilde{a}(f)$ is the Fourier transform of $a(t)$, denotes complex conjugation, and $S_n(f)$ is the one-sided power spectral density of the noise. The choice of $f_{\text{low}}$ and $f_{\text{high}}$ depends on the type of binary that we are interested in. In our study, we will set $f_{\text{low}}$ and $f_{\text{high}}$ to 20 Hz and 2048 Hz, respectively. This is sufficient for capturing the inspiral up to the moment of merger for a typical BNS system in the advanced GW detector era.

## Electromagnetic signals

**Kilonova models.** Kilonova models are extracted using the 3D Monte Carlo radiative transfer code POSSIS$^{129}$. The code can handle arbitrary geometries for the ejected material and produces spectra, lightcurves and polarization as a function of the observer viewing angle. Given an input model with defined densities $\rho$ and compositions (i.e., electron fraction $Y_e$), the code generates Monte Carlo photon packets with initial location and energy sampled from the energy distribution from radioactive decay of r-process nuclei within the model. The latter depends on the mass/density distribution of the model and the assumed nuclear heating rates and thermalization efficiencies. The frequency of each Monte Carlo photon packet is sampled according to the temperature $T$ in the ejecta, which is calculated at each time-step$^{130,131}$. Photon packets

are then followed as they diffuse out of the ejected material and interact with matter via either electron scattering or bound-bound line transitions. Time- and wavelength-dependent opacities $\kappa_{\Lambda}(\rho, T, Y_e, t)$ from Tanaka et al.$^{132}$ are implemented in the code and depend on the local properties of the ejecta ($\rho, T$, and $Y_e$). Spectral time series are extracted using the technique described by Bulla et al.$^{133}$ and used to construct broad-band lightcurves in any desired filter.

**Supernova models.** Templates available within the SNCosmo library$^{134}$ are used to model supernova spectra. Currently, the salt2 model for Type Ia supernovae and the nugent-hyper model for hypernovae associated with long GRBs are implemented in the framework and have been used in the past$^{36}$. However, the framework is flexible enough such that additional templates for different types of supernovae can be added with minimal effort.

**Kilonova/supernova inference.** Our EM inference of kilonovae and GRB afterglows is based on the AB magnitude for a specific filter $j$, $m_i^j(t_i)$. We assume these measurements to be given as a time series at times $t_i$ with a corresponding statistical error $\sigma_i^j = \sigma^j(t_i)$. The likelihood function $\mathcal{L}_{\text{EM}}(\boldsymbol{\theta})$ then reads$^{135}$:

$$ \mathcal{L}_{\text{EM}} \propto \exp \left( -\frac{1}{2} \sum_{ij} \frac{(m_i^j - m_i^{\text{est}}(\boldsymbol{\theta}))^2}{(\sigma_i^j)^2 + \sigma_{\text{sys}}^2} \right), \quad (16) $$

where $m_i^{\text{est}}(\boldsymbol{\theta})$ is the estimated AB magnitude for the parameters $\boldsymbol{\theta}$ and $\sigma_{\text{sys}}$ is the additional error budget for accounting the systematic uncertainty within the electromagnetic signal modeling. The inclusion of $\sigma_{\text{sys}}$ is equivalent to adding a shift of $\Delta m$ to the light curve, for which marginalized with respect to a zero-mean normal distribution with a variance of $\sigma_{\text{sys}}^2$.

This likelihood is equivalent to approximating the probability distribution of the spectral flux density $f_v$ to be a Log-normal distribution. The Log-normal distribution is a 2-parameter maximum entropy distribution with its support equals to the possible range for $f_v \in (0, \infty)$. There are two advantages of approximating $f_v$ with a Log-normal distribution: (1) if the uncertainty is larger or comparable to the measured value, it avoids having non-zero support for the nonphysical $f_v < 0$; (2) if the uncertainty is much smaller than the measured value, the Log-normal distribution approaches the normal distribution.

For kilonovae, we use the same model presented in Dietrich et al.$^{22}$. The model is controlled by four parameters, namely, the dynamical ejecta mass $m_{\text{dyn}}^{\text{dyn}}$, the disk wind ejecta mass $m_{\text{wind}}^{\text{dyn}}$, the half-opening angle of the lanthanide-rich component $\Phi$, and the viewing angle $\theta_{\text{obs}}$.

**GRB afterglows.** In our framework, the computation of the GRB afterglow lightcurves is until now based on the publicly available semi-analytic code afterglowpwy$^{57,58}$. The inclusion of other afterglow models is currently ongoing.

The GRB afterglow emission is produced by relativistic electrons gyrating around the magnetic field lines. These electrons are accelerated by the Fermi first-order acceleration (diffusive shock acceleration) and the magnetic field is assumed to be of turbulent nature, amplified by processes acting in collision-less shocks. The complex physics of electron acceleration at shocks is approximated by the equipartition parameters, $\epsilon_e$ and $\epsilon_B$, denoting the fraction of the shock energy that goes into the relativistic electrons and magnetic field, respectively, and $\rho$, and the slope of the electron energy distribution $dn/d\rho \propto Y^p$, with $n$ being the electron number density and $\gamma$ being the electron Lorentz factor. The flux density of the curvature radiation is:

$$ F_{\nu} = \frac{1}{4\pi d_{L}^{2}} \int d\theta d\phi R^{2} \sin(\theta) \frac{\epsilon_{\nu}}{\alpha_{\nu}} (1 - e^{-\tau}), \quad (17) $$
where $\tau$ is the optical depth and $\epsilon_v$ and $\alpha_v$ are the impassivity coefficient and absorption coefficient, respectively. For a fixed power-law distribution of electrons these can be approximated analytically¹³⁶. The synchrotron self-absorption is neglected in this work.

In order to capture the possible dependence of the GRB properties on the polar angle, the jet is discretized into a set of lateral axisymmetric (conical) layers, each of which is characterized by its initial velocity, mass, and angle. Several prescriptions for the initial angular distribution of the jet energy are available in the code. As default, we use the Gaussian jet model with $E \propto E_0 \exp(-\frac{1}{2}(\theta_c^2))$, where $\theta_c$ characterizes the width of the Gaussian. The jet truncation angle is $\theta_w$. We assume the GRB jet to be powered by the accretion of mass from the disk onto the remnant black hole¹³⁷⁻¹⁴⁰. Consequently, the jet energy is proportional to the leftover disk mass:

$$E_0 = \epsilon \times (1 - \xi) \times m_{\text{disk}}, \quad (18)$$

where $\xi$ is the fraction of disk mass ejected as wind and $\epsilon$ is the fraction of residual disk mass converted into jet energy.

The dynamical evolution of these layers is computed semi-analytically using the “thin-shell approximation” casting energy-conservation equations and shock-jump conditions into a set of evolution equations for the blast wave velocity and radius. Within blast waves, the pressure gradient perpendicular to the normal leads to lateral expansion¹⁴¹,¹⁴². In other words, the transverse pressure gradient adds the velocity along the tangent to the blast wave surface, forcing the latter to expand. The lateral expansion is important for late-time afterglow and is included in the code.

Finally, the flux density, $F_v$, is obtained by equal arrival time surface integration, Eq. (17), taking into account relativistic effects, i.e., that the observed $F_v$ is composed of contributions from different blast waves that has emitted at different comoving time and at different frequencies.

**Connecting electromagnetic signals to source properties.** To connect the observed GRB, kilonova, and GRB afterglow properties to the binary properties, we rely on phenomenological relations, i.e., fits based on numerical-relativity simulations. For our work, we use the fits presented in Kruger et al.¹⁴³ and Dietrich et al.²² but emphasize that a variety of other fitting formulas exist in the literature²⁰,⁵¹,⁵⁵,¹⁴⁴,¹⁴⁵.

In NMMA, the dynamical ejecta mass $m_{\text{dyn}}^{\text{ej}}$ is connected to the binary properties through the phenomenological relation¹⁴³:

$$\frac{m_{\text{dyn,fit}}^{\text{ej}}}{10^{-3} M_{\odot}} = \left( \frac{a}{C_1} + b \left( \frac{m_2}{m_1} \right)^n + c C_1 \right) + (1 \leftrightarrow 2), \quad (19)$$

where $m_i$ and $C_i$ are the masses and the compactness of the two components of the binary with best-fit coefficients $a = -9.3335$, $b = 114.17$, $c = -337.56$, and $n = 1.5465$. This relation enables an accurate estimation of the ejecta mass with an error well-approximated by a zero-mean Gaussian with a standard deviation 0.004 $M_{\odot}^{143}$. Therefore, the dynamical ejecta mass can be approximated as:

$$m_{\text{dyn}}^{\text{ej}} = m_{\text{dyn,fit}}^{\text{ej}} + \alpha, \quad (20)$$

where $\alpha \sim \mathcal{N}(\mu = 0, \sigma = 0.004 M_{\odot})$.

To determine the disk mass $m_{\text{disk}}$, we follow the description of Dietrich et al.²²:

$$\log_{10} \left( \frac{m_{\text{disk}}}{M_{\odot}} \right) = \max \left( -3, a \left( 1 + b \tanh \left( \frac{c - (m_1 + m_2)M_{\text{threshold}}^{-1}}{d} \right) \right) \right), \quad (21)$$

with $a$ and $b$ given by:

$$a = a_0 + \delta a \cdot \Delta, \quad b = b_0 + \delta b \cdot \Delta, \quad (22)$$

where $a_0$, $b_0$, $\delta a$, $\delta b$, $c$, and $d$ are free parameters. The parameter $\Delta$ is given by:

$$\Delta = \frac{1}{2} \tanh(\beta(q - q_{\text{trans}})), \quad (23)$$

where $q \equiv m_2/m_1 \le 1$ is the mass ratio and $\beta$ and $q_{\text{trans}}$ are free parameters. The best-fit model parameters are $a_0 = -1.581$, $\delta a = -2.439$, $b_0 = -0.538$, $\delta b = -0.406$, $c = 0.953$, $d = 0.0417$, $\beta = 3.910$, $q_{\text{trans}} = 0.900$. The threshold mass $M_{\text{threshold}}$ for a given EOS is estimated as¹⁴⁶:

$$M_{\text{threshold}} = \left( 2.38 - 3.606 \frac{M_{\text{TOV}}}{R_{1.6}} \right) M_{\text{TOV}}, \quad (24)$$

where $M_{\text{TOV}}$ and $R_{1.6}$ are the maximum mass of a non-spinning NS and the radius of a $1.6 M_{\odot}$ NS. We note that we assume that the disk-wind ejecta component is proportional to the disk mass, i.e., $m_{\text{wind}}^{\text{ej}} = \xi \times m_{\text{disk}}$.

**Bayesian statistics.** Based on Bayes' theorem, the posterior distribution of the parameters $p(\theta|d, \mathcal{H})$ under hypothesis $\mathcal{H}$ with data $d$ is given by:

$$p(\theta|d, \mathcal{H}) = \frac{p(d|\theta, \mathcal{H})p(\theta|\mathcal{H})}{p(d|\mathcal{H})} \equiv \frac{\mathcal{L}(\theta)\pi(\theta)}{\mathcal{Z}(d)}, \quad (25)$$

where $\mathcal{L}(\theta)$, $\pi(\theta)$, and $\mathcal{Z}(d)$ are the likelihood, prior, and evidence, respectively. The prior describes our knowledge of the source or model parameters prior to the experiment or observation. The likelihood and evidence quantify how well the hypothesis describes the data for a given set of parameters and over the whole parameter space, respectively. Throughout our NMMA pipeline, all data analyses use Bayes' theorem but differences appear due to the functional form of the likelihood and its specific dependence on the source parameters. For example, the GW likelihood is evaluated with a cross-correlation between the data and the GW waveform and the EM signal analysis employs a $\chi^2$ log-likelihood between the predicted lightcurves with the observed apparent magnitude data, however, from a Bayesian viewpoint their treatment is equivalent only with different likelihood functions.

In addition to the posterior estimation, the evidence $\mathcal{Z}$ carries additional information on the plausibility of a given hypothesis $\mathcal{H}$. The evidence is given by:

$$\mathcal{Z}(d|\mathcal{H}) = \int d\theta p(d|\theta, \mathcal{H}) p(\theta|\mathcal{H}) = \int d\theta \mathcal{L}(\theta) \pi(\theta), \quad (26)$$

which is the normalization constant for the posterior distribution. Moreover, we can compare the plausibilities of two hypotheses, $\mathcal{H}_1$ and $\mathcal{H}_2$, by using the odd ratio $\mathcal{O}_2^1$, which is given by:

$$\mathcal{O}_2^1 = \frac{\mathcal{Z}_1 p(\mathcal{H}_1)}{\mathcal{Z}_2 p(\mathcal{H}_2)} \equiv \mathcal{B}_2^1 \Pi_2^1, \quad (27)$$

where $\mathcal{B}_2^1$ and $\Pi_2^1$ are the Bayes factor and prior odds, respectively. If $\mathcal{O}_2^1 > \mathcal{H}_1$ is more plausible than $\mathcal{H}_2$, and vice versa.

## Data availability

The datasets generated during the current study are available in the Zenodo repository https://doi.org/10.5281/zenodo.6551053. The GW
where $\tau$ is the optical depth and $\epsilon_v$ and $\alpha_v$ are the impassivity coefficient and absorption coefficient, respectively. For a fixed power-law distribution of electrons these can be approximated analytically¹³⁶. The synchrotron self-absorption is neglected in this work.

In order to capture the possible dependence of the GRB properties on the polar angle, the jet is discretized into a set of lateral axisymmetric (conical) layers, each of which is characterized by its initial velocity, mass, and angle. Several prescriptions for the initial angular distribution of the jet energy are available in the code. As default, we use the Gaussian jet model with $E \propto E_0 \exp(-\frac{1}{2}(\theta_c^2))$, where $\theta_c$ characterizes the width of the Gaussian. The jet truncation angle is $\theta_w$. We assume the GRB jet to be powered by the accretion of mass from the disk onto the remnant black hole¹³⁷⁻¹⁴⁰. Consequently, the jet energy is proportional to the leftover disk mass:

$$E_0 = \epsilon \times (1 - \xi) \times m_{\text{disk}}, \quad (18)$$

where $\xi$ is the fraction of disk mass ejected as wind and $\epsilon$ is the fraction of residual disk mass converted into jet energy.

The dynamical evolution of these layers is computed semi-analytically using the “thin-shell approximation” casting energy-conservation equations and shock-jump conditions into a set of evolution equations for the blast wave velocity and radius. Within blast waves, the pressure gradient perpendicular to the normal leads to lateral expansion¹⁴¹,¹⁴². In other words, the transverse pressure gradient adds the velocity along the tangent to the blast wave surface, forcing the latter to expand. The lateral expansion is important for late-time afterglow and is included in the code.

Finally, the flux density, $F_v$, is obtained by equal arrival time surface integration, Eq. (17), taking into account relativistic effects, i.e., that the observed $F_v$ is composed of contributions from different blast waves that has emitted at different comoving time and at different frequencies.

**Connecting electromagnetic signals to source properties.** To connect the observed GRB, kilonova, and GRB afterglow properties to the binary properties, we rely on phenomenological relations, i.e., fits based on numerical-relativity simulations. For our work, we use the fits presented in Kruger et al.¹⁴³ and Dietrich et al.²² but emphasize that a variety of other fitting formulas exist in the literature²⁰,⁵¹,⁵⁵,¹⁴⁴,¹⁴⁵.

In NMMA, the dynamical ejecta mass $m_{\text{dyn}}^{\text{ej}}$ is connected to the binary properties through the phenomenological relation¹⁴³:

$$\frac{m_{\text{dyn,fit}}^{\text{ej}}}{10^{-3} M_{\odot}} = \left( \frac{a}{C_1} + b \left( \frac{m_2}{m_1} \right)^n + c C_1 \right) + (1 \leftrightarrow 2), \quad (19)$$

where $m_i$ and $C_i$ are the masses and the compactness of the two components of the binary with best-fit coefficients $a = -9.3335$, $b = 114.17$, $c = -337.56$, and $n = 1.5465$. This relation enables an accurate estimation of the ejecta mass with an error well-approximated by a zero-mean Gaussian with a standard deviation 0.004 $M_{\odot}^{143}$. Therefore, the dynamical ejecta mass can be approximated as:

$$m_{\text{dyn}}^{\text{ej}} = m_{\text{dyn,fit}}^{\text{ej}} + \alpha, \quad (20)$$

where $\alpha \sim \mathcal{N}(\mu = 0, \sigma = 0.004 M_{\odot})$.

To determine the disk mass $m_{\text{disk}}$, we follow the description of Dietrich et al.²²:

$$\log_{10} \left( \frac{m_{\text{disk}}}{M_{\odot}} \right) = \max \left( -3, a \left( 1 + b \tanh \left( \frac{c - (m_1 + m_2)M_{\text{threshold}}^{-1}}{d} \right) \right) \right), \quad (21)$$

with $a$ and $b$ given by:

$$a = a_0 + \delta a \cdot \Delta, \quad b = b_0 + \delta b \cdot \Delta, \quad (22)$$

where $a_0, b_0, \delta a, \delta b, c$, and $d$ are free parameters. The parameter $\Delta$ is given by:

$$\Delta = \frac{1}{2} \tanh(\beta(q - q_{\text{trans}})), \quad (23)$$

where $q = m_2/m_1 \le 1$ is the mass ratio and $\beta$ and $q_{\text{trans}}$ are free parameters. The best-fit model parameters are $a_0 = -1.581$, $\delta a = -2.439$, $b_0 = -0.538$, $\delta b = -0.406$, $c = 0.953$, $d = 0.0417$, $\beta = 3.910$, $q_{\text{trans}} = 0.900$. The threshold mass $M_{\text{threshold}}$ for a given EOS is estimated as¹⁴⁶:

$$M_{\text{threshold}} = \left( 2.38 - 3.606 \frac{M_{\text{TOV}}}{R_{1.6}} \right) M_{\text{TOV}}, \quad (24)$$

where $M_{\text{TOV}}$ and $R_{1.6}$ are the maximum mass of a non-spinning NS and the radius of a $1.6 M_{\odot}$ NS. We note that we assume that the disk-wind ejecta component is proportional to the disk mass, i.e., $m_{\text{wind}}^{\text{ej}} = \xi \times m_{\text{disk}}$.

**Bayesian statistics.** Based on Bayes' theorem, the posterior distribution of the parameters $p(\theta|d, \mathcal{H})$ under hypothesis $\mathcal{H}$ with data $d$ is given by:

$$p(\theta|d, \mathcal{H}) = \frac{p(d|\theta, \mathcal{H})p(\theta|\mathcal{H})}{p(d|\mathcal{H})} \equiv \frac{\mathcal{L}(\theta)\pi(\theta)}{\mathcal{Z}(d)}, \quad (25)$$

where $\mathcal{L}(\theta), \pi(\theta)$, and $\mathcal{Z}(d)$ are the likelihood, prior, and evidence, respectively. The prior describes our knowledge of the source or model parameters prior to the experiment or observation. The likelihood and evidence quantify how well the hypothesis describes the data for a given set of parameters and over the whole parameter space, respectively. Throughout our NMMA pipeline, all data analyses use Bayes' theorem but differences appear due to the functional form of the likelihood and its specific dependence on the source parameters. For example, the GW likelihood is evaluated with a cross-correlation between the data and the GW waveform and the EM signal analysis employs a $\chi^2$ log-likelihood between the predicted lightcurves with the observed apparent magnitude data, however, from a Bayesian viewpoint their treatment is equivalent only with different likelihood functions.

In addition to the posterior estimation, the evidence $\mathcal{Z}$ carries additional information on the plausibility of a given hypothesis $\mathcal{H}$. The evidence is given by:

$$\mathcal{Z}(d|\mathcal{H}) = \int d\theta p(d|\theta, \mathcal{H}) p(\theta|\mathcal{H}) = \int d\theta \mathcal{L}(\theta) \pi(\theta), \quad (26)$$

which is the normalization constant for the posterior distribution. Moreover, we can compare the plausibilities of two hypotheses, $\mathcal{H}_1$ and $\mathcal{H}_2$, by using the odd ratio $\mathcal{O}_2^1$, which is given by:

$$\mathcal{O}_2^1 = \frac{\mathcal{Z}_1 P(\mathcal{H}_1)}{\mathcal{Z}_2 P(\mathcal{H}_2)} \equiv \mathcal{B}_2^1 \Pi_2^1, \quad (27)$$

where $\mathcal{B}_2^1$ and $\Pi_2^1$ are the Bayes factor and prior odds, respectively. If $\mathcal{O}_2^1 > \mathcal{H}_1$ is more plausible than $\mathcal{H}_2$, and vice versa.

## Data availability

The datasets generated during the current study are available in the Zenodo repository https://doi.org/10.5281/zenodo.6551053. The GW
data strain that we have analyzed in this work was obtained from the Gravitational Wave Open Science Center (ref. 147 at https://www.gw-openscience.org), and the NICER data were obtained from Zenodo (10.5281/zenodo.3473466, 10.5281/zenodo.4670689 and 10.5281/zenodo.4697625). Source data are provided with this paper.

## Code availability

The source code of the NMMA framework, which was used for this study, is publicly available at https://github.com/nuclear-multimessenger-astronomy/nmma. In addition, all employed GW models are available on https://git.ligo.org/lscsoft. The bilby and parallel bilby software packages are available at https://git.ligo.org/lscsoft/bilby and https://git.ligo.org/lscsoft/parallel_bilby, respectively.

## References

1. Abbott, B. P. et al. GW170817: observation of gravitational waves from a binary neutron star inspiral. Phys. Rev. Lett. 119, 161101 (2017).

2. Andreoni, I. et al. Follow up of GW170817 and its electromagnetic counterpart by Australian-led observing programs. Publ. Astron. Soc. Aust. 34, e069 (2017).

3. Chornock, R. et al. The electromagnetic counterpart of the binary neutron star merger LIGO/VIRGO GW170817. IV. Detection of near-infrared signatures of r-process nucleosynthesis with Gemini-South. Astrophys. J. Lett. 848, L19 (2017).

4. Coulter, D. A. et al. Swope Supernova Survey 2017a (SSS17a), the Optical Counterpart to a Gravitational Wave Source. Science 358, 1556-1558 (2017).

5. Evans, P. et al. Swift and NuSTAR observations of GW170817: detection of a blue kilonova. Science 358, 1565 (2017).

6. Kasliwal, M. M. et al. Spitzer mid-infrared detections of neutron star merger GW170817 suggests synthesis of the heaviest elements. Mon. Not. R. Astron. Soc. 510, L7-L12 (2022).

7. Kilpatrick, C. D. et al. Electromagnetic evidence that SSS17a is the result of a binary neutron star merger. Science 358, 1583-1587 (2017).

8. Lipunov, V. M. et al. MASTER optical detection of the first LIGO/Virgo neutron star binary merger GW170817. Astrophys. J. Lett. 850, L1 (2017).

9. McCully, C. et al. The rapid reddening and featureless optical spectra of the optical counterpart of GW170817, AT 2017gfo, during the first four days. Astrophys. J. Lett. 848, L32 (2017).

10. Shappee, B. J. et al. Early spectra of the gravitational wave source GW170817: evolution of a neutron star merger. Science 358, 1574 (2017).

11. Tanvir, N. R. et al. The emergence of a lanthanide-rich kilonova following the merger of two neutron stars. Astrophys. J. 848, L27 (2017).

12. Utsumi, Y. et al. J-GEM observations of an electromagnetic counterpart to the neutron star merger GW170817. Publ. Astron. Soc. Jpn. 69, 101 (2017).

13. Margutti, R. et al. The electromagnetic counterpart of the binary neutron star merger LIGO/VIRGO GW170817. V. Rising X-ray emission from an off-axis jet. Astrophys. J. Lett. 848, L20 (2017).

14. Fong, W. et al. The optical afterglow of GW170817: an off-axis structured jet and deep constraints on a globular cluster origin. Astrophys. J. Lett. 883, L1 (2019).

15. Lamb, G. P. et al. The optical afterglow of GW170817 at one year post-merger. Astrophys. J. Lett. 870, L15 (2019).

16. Bauswein, A., Just, O., Janka, H.-T. & Stergioulas, N. Neutron-star radius constraints from GW170817 and future detections. Astrophys. J. 850, L34 (2017).

17. Ruiz, M., Shapiro, S. L. & Tsokaros, A. GW170817, general relativistic magnetohydrodynamic simulations, and the neutron star maximum mass. Phys. Rev. D97, 021501 (2018).

18. Radice, D., Perego, A., Zappa, F. & Bernuzzi, S. GW170817: joint constraint on the neutron star equation of state from multimessenger observations. Astrophys. J. 852, L29 (2018).

19. Most, E. R., Weih, L. R., Rezzolla, L. & Schaffner-Bielich, J. New constraints on radii and tidal deformabilities of neutron stars from GW170817. Phys. Rev. Lett. 120, 261103 (2018).

20. Coughlin, M. W., Dietrich, T., Margalit, B. & Metzger, B. D. Multimessenger Bayesian parameter inference of a binary neutron star merger. Mon. Not. R. Astron. Soc. Lett. 489, L91-L96 (2019).

21. Capano, C. D. et al. Stringent constraints on neutron-star radii from multimessenger observations and nuclear theory. Nat. Astron. 4, 625-632 (2020).

22. Dietrich, T. et al. Multimessenger constraints on the neutron-star equation of state and the Hubble constant. Science 370, 1450-1453 (2020).

23. Huth, S. et al. Constraining neutron-star matter with microscopic and macroscopic collisions. Nature 606, 276-280 (2022).

24. Abbott, B. P. et al. A gravitational-wave standard siren measurement of the Hubble constant. Nature https://doi.org/10.1038/nature24471 (2017).

25. Guidorzi, C. et al. Improved constraints on $H_0$ from a combined analysis of gravitational-wave and electromagnetic emission from GW170817. Astrophys. J. 851, L36 (2017).

26. Hotokezaka, K. et al. A Hubble constant measurement from superluminal motion of the jet in GW170817. Nat. Astron. 3, 940-944 (2019).

27. Coughlin, M. W. et al. On standardizing kilonovae and their use as standard candles to measure the Hubble constant. Phys. Rev. Res. 2, 022006 (2020).

28. Wang, H. & Giannios, D. Multimessenger parameter estimation of GW170817: from jet structure to the Hubble constant. Astrophys. J. 908, 200 (2021).

29. Abbott, B. et al. Multi-messenger observations of a binary neutron star merger. Astrophys. J. Lett. 848, L12 (2017).

30. Breschi, M. et al. AT2017gfo: Bayesian inference and model selection of multicomponent kilonovae and constraints on the neutron star equation of state. Mon. Not. R. Astron. Soc. 505, 1661-1677 (2021).

31. Nicholl, M. et al. Tight multimessenger constraints on the neutron star equation of state from GW170817 and a forward model for kilonova light-curve synthesis. Mon. Not. R. Astron. Soc. 505, 3016-3032 (2021).

32. Raaijmakers, G. et al. The challenges ahead for multimessenger analyses of gravitational waves and kilonova: a case study on GW190425. Astrophys. J. 922, 269 (2021).

33. Pang, P. T. H. et al. Nuclear physics multimessenger astrophysics constraints on the neutron star equation of state: adding NICER's PSR J0740+6620 measurement. Astrophys. J. 922, 14 (2021).

34. Tews, I. et al. On the nature of GW190814 and its impact on the understanding of supranuclear matter. Astrophys. J. Lett. 908, L1 (2021).

35. Andreoni, I. et al. Fast-transient searches in real time with ZTFReST: identification of three optically discovered gamma-ray burst afterglows and new constraints on the kilonova rate. Astrophys. J. 918, 63 (2021).

36. Ahumada, T. et al. Discovery and confirmation of the shortest gamma-ray burst from a collapsar. Nat. Astron. 5, 917-927 (2021).

37. Lynn, J. E. et al. Chiral three-nucleon interactions in light nuclei, neutron-α scattering, and neutron matter. Phys. Rev. Lett. 116, 062501 (2016).
data strain that we have analyzed in this work was obtained from the Gravitational Wave Open Science Center (ref. 147 at https://www.gw-openscience.org), and the NICER data were obtained from Zenodo (10.5281/zenodo.3473466, 10.5281/zenodo.4670689 and 10.5281/zenodo.4697625). Source data are provided with this paper.

## Code availability

The source code of the NMMA framework, which was used for this study, is publicly available at https://github.com/nuclear-multimessenger-astronomy/nmma. In addition, all employed GW models are available on https://git.ligo.org/lscsoft. The bilby and parallel bilby software packages are available at https://git.ligo.org/lscsoft/bilby and https://git.ligo.org/lscsoft/parallel_bilby, respectively.

## References

1. Abbott, B. P. et al. GW170817: observation of gravitational waves from a binary neutron star inspiral. Phys. Rev. Lett. 119, 161101 (2017).

2. Andreoni, I. et al. Follow up of GW170817 and its electromagnetic counterpart by Australian-led observing programs. Publ. Astron. Soc. Aust. 34, e069 (2017).

3. Chornock, R. et al. The electromagnetic counterpart of the binary neutron star merger LIGO/VIRGO GW170817. IV. Detection of near-infrared signatures of r-process nucleosynthesis with Gemini-South. Astrophys. J. Lett. 848, L19 (2017).

4. Coulter, D. A. et al. Swope Supernova Survey 2017a (SSS17a), the Optical Counterpart to a Gravitational Wave Source. Science 358, 1556-1558 (2017).

5. Evans, P. et al. Swift and NuSTAR observations of GW170817: detection of a blue kilonova. Science 358, 1565 (2017).

6. Kasliwal, M. M. et al. Spitzer mid-infrared detections of neutron star merger GW170817 suggests synthesis of the heaviest elements. Mon. Not. R. Astron. Soc. 510, L7-L12 (2022).

7. Kilpatrick, C. D. et al. Electromagnetic evidence that SSS17a is the result of a binary neutron star merger. Science 358, 1583-1587 (2017).

8. Lipunov, V. M. et al. MASTER optical detection of the first LIGO/Virgo neutron star binary merger GW170817. Astrophys. J. Lett. 850, L1 (2017).

9. McCully, C. et al. The rapid reddening and featureless optical spectra of the optical counterpart of GW170817, AT 2017gfo, during the first four days. Astrophys. J. Lett. 848, L32 (2017).

10. Shappee, B. J. et al. Early spectra of the gravitational wave source GW170817: evolution of a neutron star merger. Science 358, 1574 (2017).

11. Tanvir, N. R. et al. The emergence of a lanthanide-rich kilonova following the merger of two neutron stars. Astrophys. J. 848, L27 (2017).

12. Utsumi, Y. et al. J-GEM observations of an electromagnetic counterpart to the neutron star merger GW170817. Publ. Astron. Soc. Jpn. 69, 101 (2017).

13. Margutti, R. et al. The electromagnetic counterpart of the binary neutron star merger LIGO/VIRGO GW170817. V. Rising X-ray emission from an off-axis jet. Astrophys. J. Lett. 848, L20 (2017).

14. Fong, W. et al. The optical afterglow of GW170817: an off-axis structured jet and deep constraints on a globular cluster origin. Astrophys. J. Lett. 883, L1 (2019).

15. Lamb, G. P. et al. The optical afterglow of GW170817 at one year post-merger. Astrophys. J. Lett. 870, L15 (2019).

16. Bauswein, A., Just, O., Janka, H.-T. & Stergioulas, N. Neutron-star radius constraints from GW170817 and future detections. Astrophys. J. 850, L34 (2017).

17. Ruiz, M., Shapiro, S. L. & Tsokaros, A. GW170817, general relativistic magnetohydrodynamic simulations, and the neutron star maximum mass. Phys. Rev. D97, 021501 (2018).

18. Radice, D., Perego, A., Zappa, F. & Bernuzzi, S. GW170817: joint constraint on the neutron star equation of state from multimessenger observations. Astrophys. J. 852, L29 (2018).

19. Most, E. R., Weih, L. R., Rezzolla, L. & Schaffner-Bielich, J. New constraints on radii and tidal deformabilities of neutron stars from GW170817. Phys. Rev. Lett. 120, 261103 (2018).

20. Coughlin, M. W., Dietrich, T., Margalit, B. & Metzger, B. D. Multimessenger Bayesian parameter inference of a binary neutron star merger. Mon. Not. R. Astron. Soc. Lett. 489, L91-L96 (2019).

21. Capano, C. D. et al. Stringent constraints on neutron-star radii from multimessenger observations and nuclear theory. Nat. Astron. 4, 625-632 (2020).

22. Dietrich, T. et al. Multimessenger constraints on the neutron-star equation of state and the Hubble constant. Science 370, 1450-1453 (2020).

23. Huth, S. et al. Constraining neutron-star matter with microscopic and macroscopic collisions. Nature 606, 276-280 (2022).

24. Abbott, B. P. et al. A gravitational-wave standard siren measurement of the Hubble constant. Nature https://doi.org/10.1038/nature24471 (2017).

25. Guidorzi, C. et al. Improved constraints on $H_0$ from a combined analysis of gravitational-wave and electromagnetic emission from GW170817. Astrophys. J. 851, L36 (2017).

26. Hotokezaka, K. et al. A Hubble constant measurement from superluminal motion of the jet in GW170817. Nat. Astron. 3, 940-944 (2019).

27. Coughlin, M. W. et al. On standardizing kilonovae and their use as standard candles to measure the Hubble constant. Phys. Rev. Res. 2, 022006 (2020).

28. Wang, H. & Giannios, D. Multimessenger parameter estimation of GW170817: from jet structure to the Hubble constant. Astrophys. J. 908, 200 (2021).

29. Abbott, B. et al. Multi-messenger observations of a binary neutron star merger. Astrophys. J. Lett. 848, L12 (2017).

30. Breschi, M. et al. AT2017gfo: Bayesian inference and model selection of multicomponent kilonovae and constraints on the neutron star equation of state. Mon. Not. R. Astron. Soc. 505, 1661-1677 (2021).

31. Nicholl, M. et al. Tight multimessenger constraints on the neutron star equation of state from GW170817 and a forward model for kilonova light-curve synthesis. Mon. Not. R. Astron. Soc. 505, 3016-3032 (2021).

32. Raaijmakers, G. et al. The challenges ahead for multimessenger analyses of gravitational waves and kilonova: a case study on GW190425. Astrophys. J. 922, 269 (2021).

33. Pang, P. T. H. et al. Nuclear physics multimessenger astrophysics constraints on the neutron star equation of state: adding NICER's PSR J0740+6620 measurement. Astrophys. J. 922, 14 (2021).

34. Tews, I. et al. On the nature of GW190814 and its impact on the understanding of supranuclear matter. Astrophys. J. Lett. 908, L1 (2021).

35. Andreoni, I. et al. Fast-transient searches in real time with ZTFReST: identification of three optically discovered gamma-ray burst afterglows and new constraints on the kilonova rate. Astrophys. J. 918, 63 (2021).

36. Ahumada, T. et al. Discovery and confirmation of the shortest gamma-ray burst from a collapsar. Nat. Astron. 5, 917-927 (2021).

37. Lynn, J. E. et al. Chiral three-nucleon interactions in light nuclei, neutron-α scattering, and neutron matter. Phys. Rev. Lett. 116, 062501 (2016).
38. Holt, J. W. & Kaiser, N. Equation of state of nuclear and neutron matter at third-order in perturbation theory from chiral effective field theory. Phys. Rev. C95, 034326 (2017).

39. Drischler, C., Hebeler, K. & Schwenk, A. Chiral interactions up to next-to-next-to-next-to-leading order and nuclear saturation. Phys. Rev. Lett. 122, 042501 (2019).

40. Piarulli, M., Bombaci, I., Logoteta, D., Lovato, A. & Wiringa, R. B. Benchmark calculations of pure neutron matter with realistic nucleon-nucleon interactions. Phys. Rev. C 101, 045801 (2020).

41. Keller, J., Wellenhofer, C., Hebeler, K. & Schwenk, A. Neutron matter at finite temperature based on chiral effective field theory interactions. Phys. Rev. C 103, 055806 (2021).

42. Russotto, P. et al. Results of the ASY-EOS experiment at GSI: the symmetry energy at suprasaturation density. Phys. Rev. C 94, 034608 (2016).

43. Adhikari, D. et al. Accurate determination of the neutron skin thickness of ²⁰⁸Pb through parity-violation in electron scattering. Phys. Rev. Lett. 126, 172502 (2021).

44. Antoniadis, J. et al. A massive pulsar in a compact relativistic binary. Science 340, 6131 (2013).

45. Arzoumanian, Z. et al. The NANOGrav 11-year data set: high-precision timing of 45 millisecond pulsars. Astrophys. J. Suppl. 235, 37 (2018).

46. Cromartie, H. T. et al. Relativistic Shapiro delay measurements of an extremely massive millisecond pulsar. Nat. Astron. 4, 72–76 (2019).

47. Kurkela, A., Fraga, E. S., Schaffner-Bielich, J. & Vuorinen, A. Constraining neutron star matter with quantum chromodynamics. Astrophys. J. 789, 127 (2014).

48. Komoltsev, O. & Kurkela, A. How perturbative QCD constrains the equation of state at neutron-star densities. Phys. Rev. Lett. 128, 202701 (2022).

49. Vallisneri, M. Use and abuse of the Fisher information matrix in the assessment of gravitational-wave parameter-estimation prospects. Phys. Rev. D77, 042001 (2008).

50. Abbott, B. P. et al. Properties of the binary neutron star merger GW170817. Phys. Rev. X9, 011001 (2019).

51. Coughlin, M. W. et al. Constraints on the neutron star equation of state from at2017gfo using radiative transfer simulations. Mon. Not. R. Astron. Soc. 480, 3871–3878 (2018).

52. Heinzel, J. et al. Comparing inclination dependent analyses of kilonova transients. Mon. Not. R. Astron. Soc. 502, 3057–3065 (2021).

53. Bulla, M. The critical role of nuclear heating rates, thermalization efficiencies and opacities for kilonova modelling and parameter inference. Mon. Not. R. Astron. Soc. 520, 2558–2570 (2023).

54. Siegel, D. M. GW170817—the first observed neutron star merger and its kilonova: implications for the astrophysical site of the r-process. Eur. Phys. J. A 55, 203 (2019).

55. Radice, D. et al. Binary neutron star mergers: mass ejection, electromagnetic counterparts and nucleosynthesis. Astrophys. J. 869, 130 (2018).

56. Troja, E. et al. A year in the life of GW 170817: the rise and fall of a structured jet from a binary neutron star merger. Mon. Not. R. Astron. Soc. 489, 1919–1926 (2019).

57. van Eerten, H., Zhang, W. & MacFadyen, A. Off-axis gamma-ray burst afterglow modeling based on a two-dimensional axisymmetric hydrodynamics simulation. Astrophys. J. 722, 235–247 (2010).

58. Ryan, G., van Eerten, H., Piro, L. & Troja, E. Gamma-ray burst afterglows in the multimessenger era: numerical models and closure relations. Astrophys. J. 896, 166 (2020).

59. Shrestha, M. et al. Impact of jets on kilonova photometric and polarimetric emission from binary neutron star mergers. Mon. Not. R. Astron. Soc. 523, 2990–3000 (2023).

60. Anand, S. et al. Chemical distribution of the dynamical ejecta in the neutron star merger GW170817. https://doi.org/10.48550/arXiv.2307.11080 (2023).

61. Baym, G., Pethick, C. & Sutherland, P. The ground state of matter at high densities: equation of state and stellar models. Astrophys. J. 170, 299 (1971).

62. Chamel, N. & Haensel, P. Physics of neutron star crusts. Living Rev. Relativ. 11, 10 (2008).

63. Tews, I. et al. Nuclear forces for precision nuclear physics—a collection of perspectives. Few-Body Syst 63, 67 (2022).

64. Krüger, T., Tews, I., Hebeler, K. & Schwenk, A. Neutron matter from chiral effective field theory interactions. Phys. Rev. C 88, 025802 (2013).

65. Epelbaum, E., Hammer, H.-W. & Meissner, U.-G. Modern theory of nuclear forces. Rev. Mod. Phys. 81, 1773–1825 (2009).

66. Machleidt, R. & Entem, D. R. Chiral effective field theory and nuclear forces. Phys. Rep. 503, 1–75 (2011).

67. Drischler, C., Melendez, J. A., Furnstahl, R. J. & Phillips, D. R. Quantifying uncertainties and correlations in the nuclear-matter equation of state. Phys. Rev. C 102, 054315 (2020).

68. Weinberg, S. Nuclear forces from chiral Lagrangians. Phys. Lett. B 251, 288–292 (1990).

69. Weinberg, S. Effective chiral Lagrangians for nucleon-pion interactions and nuclear forces. Nucl. Phys. B 363, 3–18 (1991).

70. Weinberg, S. Three body interactions among nucleons and pions. Phys. Lett. B 295, 114–121 (1992).

71. Epelbaum, E., Krebs, H. & Meißner, U.-G. Improved chiral nucleon-nucleon potential up to next-to-next-to-next-to-leading order. Eur. Phys. J. A 51, 53 (2015).

72. Drischler, C., Furnstahl, R. J., Melendez, J. A. & Phillips, D. R. How well do we know the neutron-matter equation of state at the densities inside neutron stars? A Bayesian approach with correlated uncertainties. Phys. Rev. Lett. 125, 202702 (2020).

73. Ekström, A., Hagen, G., Morris, T. D., Papenbrock, T. & Schwartz, P. D. Δ isobars and nuclear saturation. Phys. Rev. C 97, 024332 (2018).

74. Tews, I., Krüger, T., Hebeler, K. & Schwenk, A. Neutron matter at next-to-next-to-next-to-leading order in chiral effective field theory. Phys. Rev. Lett. 110, 032504 (2013).

75. Carbone, A., Cipollone, A., Barbieri, C., Rios, A. & Polls, A. Self-consistent Green's functions formalism with three-body interactions. Phys. Rev. C 88, 054326 (2013).

76. Hagen, G. et al. Coupled-cluster calculations of nucleonic matter. Phys. Rev. C89, 014319 (2014).

77. Carlson, J. et al. Quantum Monte Carlo methods for nuclear physics. Rev. Mod. Phys. 87, 1067 (2015).

78. Schmidt, K. E. & Fantoni, S. A quantum Monte Carlo method for nucleon systems. Phys. Lett. B 446, 99–103 (1999).

79. Tews, I. Spectrum of shear modes in the neutron-star crust: estimating the nuclear-physics uncertainties. Phys. Rev. C95, 015803 (2017).

80. Tews, I., Carlson, J., Gandolfi, S. & Reddy, S. Constraining the speed of sound inside neutron stars with chiral effective field theory interactions and observations. Astrophys. J. 860, 149 (2018).

81. Read, J. S., Lackey, B. D., Owen, B. J. & Friedman, J. L. Constraints on a phenomenologically parameterized neutron-star equation of state. Phys. Rev. D79, 124032 (2009).

82. Hebeler, K., Lattimer, J., Pethick, C. & Schwenk, A. Equation of state and neutron star properties constrained by nuclear physics and observation. Astrophys. J. 773, 11 (2013).

83. Annala, E., Gorda, T., Kurkela, A. & Vuorinen, A. Gravitational-wave constraints on the neutron-star-matter equation of state. Phys. Rev. Lett. 120, 172703 (2018).
38. Holt, J. W. & Kaiser, N. Equation of state of nuclear and neutron matter at third-order in perturbation theory from chiral effective field theory. Phys. Rev. C95, 034326 (2017).

39. Drischler, C., Hebeler, K. & Schwenk, A. Chiral interactions up to next-to-next-to-next-to-leading order and nuclear saturation. Phys. Rev. Lett. 122, 042501 (2019).

40. Piarulli, M., Bombaci, I., Logoteta, D., Lovato, A. & Wiringa, R. B. Benchmark calculations of pure neutron matter with realistic nucleon-nucleon interactions. Phys. Rev. C 101, 045801 (2020).

41. Keller, J., Wellenhofer, C., Hebeler, K. & Schwenk, A. Neutron matter at finite temperature based on chiral effective field theory interactions. Phys. Rev. C 103, 055806 (2021).

42. Russotto, P. et al. Results of the ASY-EOS experiment at GSI: the symmetry energy at suprasaturation density. Phys. Rev. C 94, 034608 (2016).

43. Adhikari, D. et al. Accurate determination of the neutron skin thickness of ²⁰⁸Pb through parity-violation in electron scattering. Phys. Rev. Lett. 126, 172502 (2021).

44. Antoniadis, J. et al. A massive pulsar in a compact relativistic binary. Science 340, 6131 (2013).

45. Arzoumanian, Z. et al. The NANOGrav 11-year data set: high-precision timing of 45 millisecond pulsars. Astrophys. J. Suppl. 235, 37 (2018).

46. Cromartie, H. T. et al. Relativistic Shapiro delay measurements of an extremely massive millisecond pulsar. Nat. Astron. 4, 72–76 (2019).

47. Kurkela, A., Fraga, E. S., Schaffner-Bielich, J. & Vuorinen, A. Constraining neutron star matter with quantum chromodynamics. Astrophys. J. 789, 127 (2014).

48. Komoltsev, O. & Kurkela, A. How perturbative QCD constrains the equation of state at neutron-star densities. Phys. Rev. Lett. 128, 202701 (2022).

49. Vallisneri, M. Use and abuse of the Fisher information matrix in the assessment of gravitational-wave parameter-estimation prospects. Phys. Rev. D77, 042001 (2008).

50. Abbott, B. P. et al. Properties of the binary neutron star merger GW170817. Phys. Rev. X9, 011001 (2019).

51. Coughlin, M. W. et al. Constraints on the neutron star equation of state from at2017gfo using radiative transfer simulations. Mon. Not. R. Astron. Soc. 480, 3871–3878 (2018).

52. Heinzel, J. et al. Comparing inclination dependent analyses of kilonova transients. Mon. Not. R. Astron. Soc. 502, 3057–3065 (2021).

53. Bulla, M. The critical role of nuclear heating rates, thermalization efficiencies and opacities for kilonova modelling and parameter inference. Mon. Not. R. Astron. Soc. 520, 2558–2570 (2023).

54. Siegel, D. M. GW170817—the first observed neutron star merger and its kilonova: implications for the astrophysical site of the r-process. Eur. Phys. J. A 55, 203 (2019).

55. Radice, D. et al. Binary neutron star mergers: mass ejection, electromagnetic counterparts and nucleosynthesis. Astrophys. J. 869, 130 (2018).

56. Troja, E. et al. A year in the life of GW 170817: the rise and fall of a structured jet from a binary neutron star merger. Mon. Not. R. Astron. Soc. 489, 1919–1926 (2019).

57. van Eerten, H., Zhang, W. & MacFadyen, A. Off-axis gamma-ray burst afterglow modeling based on a two-dimensional axisymmetric hydrodynamics simulation. Astrophys. J. 722, 235–247 (2010).

58. Ryan, G., van Eerten, H., Piro, L. & Troja, E. Gamma-ray burst afterglows in the multimessenger era: numerical models and closure relations. Astrophys. J. 896, 166 (2020).

59. Shrestha, M. et al. Impact of jets on kilonova photometric and polarimetric emission from binary neutron star mergers. Mon. Not. R. Astron. Soc. 523, 2990–3000 (2023).

60. Anand, S. et al. Chemical distribution of the dynamical ejecta in the neutron star merger GW170817. https://doi.org/10.48550/arXiv.2307.11080 (2023).

61. Baym, G., Pethick, C. & Sutherland, P. The ground state of matter at high densities: equation of state and stellar models. Astrophys. J. 170, 299 (1971).

62. Chamel, N. & Haensel, P. Physics of neutron star crusts. Living Rev. Relativ. 11, 10 (2008).

63. Tews, I. et al. Nuclear forces for precision nuclear physics—a collection of perspectives. Few-Body Syst 63, 67 (2022).

64. Krüger, T., Tews, I., Hebeler, K. & Schwenk, A. Neutron matter from chiral effective field theory interactions. Phys. Rev. C 88, 025802 (2013).

65. Epelbaum, E., Hammer, H.-W. & Meissner, U.-G. Modern theory of nuclear forces. Rev. Mod. Phys. 81, 1773–1825 (2009).

66. Machleidt, R. & Entem, D. R. Chiral effective field theory and nuclear forces. Phys. Rep. 503, 1–75 (2011).

67. Drischler, C., Melendez, J. A., Furnstahl, R. J. & Phillips, D. R. Quantifying uncertainties and correlations in the nuclear-matter equation of state. Phys. Rev. C 102, 054315 (2020).

68. Weinberg, S. Nuclear forces from chiral Lagrangians. Phys. Lett. B 251, 288–292 (1990).

69. Weinberg, S. Effective chiral Lagrangians for nucleon–pion interactions and nuclear forces. Nucl. Phys. B 363, 3–18 (1991).

70. Weinberg, S. Three body interactions among nucleons and pions. Phys. Lett. B 295, 114–121 (1992).

71. Epelbaum, E., Krebs, H. & Meißner, U.-G. Improved chiral nucleon-nucleon potential up to next-to-next-to-next-to-leading order. Eur. Phys. J. A 51, 53 (2015).

72. Drischler, C., Furnstahl, R. J., Melendez, J. A. & Phillips, D. R. How well do we know the neutron-matter equation of state at the densities inside neutron stars? A Bayesian approach with correlated uncertainties. Phys. Rev. Lett. 125, 202702 (2020).

73. Ekström, A., Hagen, G., Morris, T. D., Papenbrock, T. & Schwartz, P. D. Δ isobars and nuclear saturation. Phys. Rev. C 97, 024332 (2018).

74. Tews, I., Krüger, T., Hebeler, K. & Schwenk, A. Neutron matter at next-to-next-to-next-to-leading order in chiral effective field theory. Phys. Rev. Lett. 110, 032504 (2013).

75. Carbone, A., Cipollone, A., Barbieri, C., Rios, A. & Polls, A. Self-consistent Green's functions formalism with three-body interactions. Phys. Rev. C 88, 054326 (2013).

76. Hagen, G. et al. Coupled-cluster calculations of nucleonic matter. Phys. Rev. C89, 014319 (2014).

77. Carlson, J. et al. Quantum Monte Carlo methods for nuclear physics. Rev. Mod. Phys. 87, 1067 (2015).

78. Schmidt, K. E. & Fantoni, S. A quantum Monte Carlo method for nucleon systems. Phys. Lett. B 446, 99–103 (1999).

79. Tews, I. Spectrum of shear modes in the neutron-star crust: estimating the nuclear-physics uncertainties. Phys. Rev. C95, 015803 (2017).

80. Tews, I., Carlson, J., Gandolfi, S. & Reddy, S. Constraining the speed of sound inside neutron stars with chiral effective field theory interactions and observations. Astrophys. J. 860, 149 (2018).

81. Read, J. S., Lackey, B. D., Owen, B. J. & Friedman, J. L. Constraints on a phenomenologically parameterized neutron-star equation of state. Phys. Rev. D79, 124032 (2009).

82. Hebeler, K., Lattimer, J., Pethick, C. & Schwenk, A. Equation of state and neutron star properties constrained by nuclear physics and observation. Astrophys. J. 773, 11 (2013).

83. Annala, E., Gorda, T., Kurkela, A. & Vuorinen, A. Gravitational-wave constraints on the neutron-star-matter equation of state. Phys. Rev. Lett. 120, 172703 (2018).
84. Tews, I., Margueron, J. & Reddy, S. Critical examination of con-straints on the equation of state of dense matter obtained from GW170817. Phys. Rev. C 98, 045804 (2018).

85. Greif, S., Raaijmakers, G., Hebeler, K., Schwenk, A. & Watts, A. Equation of state sensitivities when inferring neutron star and dense matter properties. Mon. Not. R. Astron. Soc. 485, 5363-5376 (2019).

86. Essick, R., Landry, P. & Holz, D. E. Nonparametric inference of neutron star composition, equation of state, and maximum mass with GW170817. Phys. Rev. D 101, 063007 (2020).

87. McLerran, L. & Reddy, S. Quarkyonic matter and neutron stars. Phys. Rev. Lett. 122, 122701 (2019).

88. Jeong, K. S., McLerran, L. & Sen, S. Dynamically generated momentum space shell structure of quarkyonic matter via an excluded volume model. Phys. Rev. C 101, 035201 (2020).

89. Sen, S. & Sivertsen, L. Mass and radius relations of quarkyonic stars using an excluded-volume model. Astrophys. J. 915, 109 (2021).

90. Margueron, J., Hansen, H., Proust, P. & Chanfray, G. Quarkyonic stars with isospin-flavor asymmetry. Phys. Rev. C 104, 055803 (2021).

91. Gandolfi, S., Carlson, J. & Reddy, S. The maximum mass and radius of neutron stars and the nuclear symmetry energy. Phys. Rev. C 85, 032801 (2012).

92. Essick, R., Tews, I., Landry, P., Reddy, S. & Holz, D. E. Direct astrophysical tests of chiral effective field theory at supranuclear densities. Phys. Rev. C 102, 055803 (2020).

93. Wolff, M. T. et al. NICER detection of thermal X-ray pulsations from the massive millisecond pulsars PSR J0740+6620 and PSR J1614-2230. Astrophys. J. Lett. 918, L26 (2021).

94. Miller, M. C. et al. PSR J0030+0451 mass and radius from NICER data and implications for the properties of neutron star matter. Astrophys. J. Lett. 887, L24 (2019).

95. Riley, T. E. et al. A NICER view of PSR J0030+0451: millisecond pulsar parameter estimation. Astrophys. J. Lett. 887, L21 (2019).

96. Miller, M. C. et al. The radius of PSR J0740+6620 from NICER and XMM-Newton data. Astrophys. J. Lett. 918, L28 (2021).

97. Riley, T. E. et al. A NICER view of the massive pulsar PSR J0740+6620 informed by radio timing and XMM-newton spectroscopy. Astrophys. J. Lett. 918, L27 (2021).

98. Damour, T. & Nagar, A. Relativistic tidal properties of neutron stars. Phys. Rev. D80, 084035 (2009).

99. Hinderer, T., Lackey, B. D., Lang, R. N. & Read, J. S. Tidal deformability of neutron stars with realistic equations of state and their gravitational wave signatures in binary inspiral. Phys. Rev. D81, 123016 (2010).

100. Abbott, B. P. et al. Multi-messenger observations of a binary neutron star merger. Astrophys. J. Lett. 848, L12 (2017).

101. Abbott, B. P. et al. GW170817: measurements of neutron star radii and equation of state. Phys. Rev. Lett. 121, 161101 (2018).

102. Abbott, B. P. et al. Model comparison from LIGO-Virgo data on GW170817's binary components and consequences for the merger remnant. Class. Quant. Grav. 37, 045006 (2020).

103. Carney, M. F., Wade, L. E. & Irwin, B. S. Comparing two models for measuring the neutron star equation of state from gravitational-wave signals. Phys. Rev. D 98, 063004 (2018).

104. Damour, T., Nagar, A. & Villain, L. Measurability of the tidal polarizability of neutron stars in late-inspiral gravitational-wave signals. Phys. Rev. D85, 123007 (2012).

105. Abdelsalhin, T., Gualtieri, L. & Pani, P. Post-Newtonian spin-tidal couplings for compact binaries. Phys. Rev. D 98, 104046 (2018).

106. Landry, P. Rotational-tidal phasing of the binary neutron star waveform. arXiv https://doi.org/10.48550/arXiv.1805.01882 (2018).

107. Jiménez Forteza, X., Abdelsalhin, T., Pani, P. & Gualtieri, L. Impact of high-order tidal terms on binary neutron-star waveforms. Phys. Rev. D98, 124014 (2018).

108. Buonanno, A. & Damour, T. Effective one-body approach to general relativistic two-body dynamics. Phys. Rev. D59, 084006 (1999).

109. Damour, T. & Nagar, A. Effective one body description of tidal effects in inspiralling compact binaries. Phys. Rev. D81, 084016 (2010).

110. Hotokezaka, K., Kyutoku, K., Okawa, H. & Shibata, M. Exploring tidal effects of coalescing binary neutron stars in numerical relativity. II. Long-term simulations. Phys. Rev. D91, 064060 (2015).

111. Hinderer, T. et al. Effects of neutron-star dynamic tides on gravitational waveforms within the effective-one-body approach. Phys. Rev. Lett. 116, 181101 (2016).

112. Akcay, S. et al. Effective-one-body multipolar waveform for tidally interacting binary neutron stars up to merger. Phys. Rev. D99, 044051 (2019).

113. Bernuzzi, S., Nagar, A., Dietrich, T. & Damour, T. Modeling the dynamics of tidally interacting binary neutron stars up to the merger. Phys. Rev. Lett. 114, 161103 (2015).

114. Dietrich, T. & Hinderer, T. Comprehensive comparison of numerical relativity and effective-one-body results to inform improvements in waveform models for binary neutron star systems. Phys. Rev. D95, 124006 (2017).

115. Nagar, A. et al. Time-domain effective-one-body gravitational waveforms for coalescing compact binaries with nonprecessing spins, tides and self-spin effects. Phys. Rev. D98, 104052 (2018).

116. Steinhoff, J., Hinderer, T., Dietrich, T. & Foucart, F. Spin effects on neutron star fundamental-mode dynamical tides: phenomenology and comparison to numerical simulations. Phys. Rev. Res. 3, 033129 (2021).

117. Kawaguchi, K. et al. Frequency-domain gravitational waveform models for inspiraling binary neutron stars. Phys. Rev. D97, 044044 (2018).

118. Dietrich, T., Bernuzzi, S. & Tichy, W. Closed-form tidal approximations for binary neutron star gravitational waveforms constructed from high-resolution numerical relativity simulations. Phys. Rev. D96, 121501 (2017).

119. Dietrich, T. et al. Matter imprints in waveform models for neutron star binaries: tidal and self-spin effects. Phys. Rev. D99, 024029 (2019).

120. Dietrich, T. et al. Improving the NRTidal model for binary neutron star systems. Phys. Rev. D100, 044003 (2019).

121. Thompson, J. E. et al. Modeling the gravitational wave signature of neutron star black hole coalescences. Phys. Rev. D. 101, 124059 (2020).

122. LIGO Scientific Collaboration. LIGO Algorithm Library—LAL-Suite. (2018).

123. Steinhoff, J., Hinderer, T., Buonanno, A. & Taracchini, A. Dynamical tides in general relativity: effective action and effective-one-body Hamiltonian. Phys. Rev. D94, 104028 (2016).

124. Lackey, B. D., Pürrer, M., Taracchini, A. & Marsat, S. Surrogate model for an aligned-spin effective one body waveform model of binary neutron star inspirals using Gaussian process regression. Phys. Rev. D100, 024002 (2019).

125. Nagar, A. et al. Nonlinear-in-spin effects in effective-one-body waveform models of spin-aligned, inspiralling, neutron star binaries. Phys. Rev. D 99, 044007 (2019).

126. Nagar, A. & Rettegno, P. Efficient effective one body time-domain gravitational waveforms. Phys. Rev. D99, 021501 (2019).
84. Tews, I., Margueron, J. & Reddy, S. Critical examination of con-
straints on the equation of state of dense matter obtained from
GW170817. Phys. Rev. C **98**, 045804 (2018).

85. Greif, S., Raaijmakers, G., Hebeler, K., Schwenk, A. & Watts, A. Equation of state sensitivities when inferring neutron star and dense matter properties. *Mon. Not. R. Astron. Soc.* **485**, 5363–5376 (2019).

86. Essick, R., Landry, P. & Holz, D. E. Nonparametric inference of neutron star composition, equation of state, and maximum mass with GW170817. Phys. Rev. D **101**, 063007 (2020).

87. McLerran, L. & Reddy, S. Quarkyonic matter and neutron stars. Phys. Rev. Lett. **122**, 122701 (2019).

88. Jeong, K. S., McLerran, L. & Sen, S. Dynamically generated momentum space shell structure of quarkyonic matter via an excluded volume model. Phys. Rev. C **101**, 035201 (2020).

89. Sen, S. & Sivertsen, L. Mass and radius relations of quarkyonic stars using an excluded-volume model. Astrophys. J. **915**, 109 (2021).

90. Margueron, J., Hansen, H., Proust, P. & Chanfray, G. Quarkyonic stars with isospin-flavor asymmetry. Phys. Rev. C **104**, 055803 (2021).

91. Gandolfi, S., Carlson, J. & Reddy, S. The maximum mass and radius of neutron stars and the nuclear symmetry energy. Phys. Rev. C **85**, 032801 (2012).

92. Essick, R., Tews, I., Landry, P., Reddy, S. & Holz, D. E. Direct astrophysical tests of chiral effective field theory at supranuclear densities. Phys. Rev. C **102**, 055803 (2020).

93. Wolff, M. T. et al. NICER detection of thermal X-ray pulsations from the massive millisecond pulsars PSR J0740+6620 and PSR J1614–2230. Astrophys. J. Lett. **918**, L26 (2021).

94. Miller, M. C. et al. PSR J0030+0451 mass and radius from NICER data and implications for the properties of neutron star matter. Astrophys. J. Lett. **887**, L24 (2019).

95. Riley, T. E. et al. A NICER view of PSR J0030+0451: millisecond pulsar parameter estimation. Astrophys. J. Lett. **887**, L21 (2019).

96. Miller, M. C. et al. The radius of PSR J0740+6620 from NICER and XMM-Newton data. Astrophys. J. Lett. **918**, L28 (2021).

97. Riley, T. E. et al. A NICER view of the massive pulsar PSR J0740+6620 informed by radio timing and XMM-newton spectroscopy. Astrophys. J. Lett. **918**, L27 (2021).

98. Damour, T. & Nagar, A. Relativistic tidal properties of neutron stars. Phys. Rev. D80, 084035 (2009).

99. Hinderer, T., Lackey, B. D., Lang, R. N. & Read, J. S. Tidal deformability of neutron stars with realistic equations of state and their gravitational wave signatures in binary inspiral. Phys. Rev. D81, 123016 (2010).

100. Abbott, B. P. et al. Multi-messenger observations of a binary neutron star merger. Astrophys. J. Lett. **848**, L12 (2017).

101. Abbott, B. P. et al. GW170817: measurements of neutron star radii and equation of state. Phys. Rev. Lett. **121**, 161101 (2018).

102. Abbott, B. P. et al. Model comparison from LIGO-Virgo data on GW170817's binary components and consequences for the merger remnant. Class. Quant. Grav. **37**, 045006 (2020).

103. Carney, M. F., Wade, L. E. & Irwin, B. S. Comparing two models for measuring the neutron star equation of state from gravitational-wave signals. Phys. Rev. D **98**, 063004 (2018).

104. Damour, T., Nagar, A. & Villain, L. Measurability of the tidal polarizability of neutron stars in late-inspiral gravitational-wave signals. Phys. Rev. D **85**, 123007 (2012).

105. Abdelsalhin, T., Gualtieri, L. & Pani, P. Post-Newtonian spin-tidal couplings for compact binaries. Phys. Rev. D **98**, 104046 (2018).

106. Landry, P. Rotational-tidal phasing of the binary neutron star waveform. arXiv https://doi.org/10.48550/arXiv.1805.01882 (2018).

107. Jiménez Forteza, X., Abdelsalhin, T., Pani, P. & Gualtieri, L. Impact of high-order tidal terms on binary neutron-star waveforms. Phys. Rev. D98, 124014 (2018).

108. Buonanno, A. & Damour, T. Effective one-body approach to gen-
eral relativistic two-body dynamics. Phys. Rev. D59,
084006 (1999).

109. Damour, T. & Nagar, A. Effective one body description of tidal effects in inspiralling compact binaries. Phys. Rev. D81, 084016 (2010).

110. Hotokezaka, K., Kyutoku, K., Okawa, H. & Shibata, M. Exploring tidal effects of coalescing binary neutron stars in numerical relativity. II. Long-term simulations. Phys. Rev. D91, 064060 (2015).

111. Hinderer, T. et al. Effects of neutron-star dynamic tides on grav-
itational waveforms within the effective-one-body approach. Phys.
Rev. Lett. **116**, 181101 (2016).

112. Akcay, S. et al. Effective-one-body multipolar waveform for tidally interacting binary neutron stars up to merger. Phys. Rev. D99, 044051 (2019).

113. Bernuzzi, S., Nagar, A., Dietrich, T. & Damour, T. Modeling the dynamics of tidally interacting binary neutron stars up to the merger. Phys. Rev. Lett. **114**, 161103 (2015).

114. Dietrich, T. & Hinderer, T. Comprehensive comparison of numer-
ical relativity and effective-one-body results to inform improve-
ments in waveform models for binary neutron star systems. Phys.
Rev. D95, 124006 (2017).

115. Nagar, A. et al. Time-domain effective-one-body gravitational waveforms for coalescing compact binaries with nonprecessing spins, tides and self-spin effects. Phys. Rev. D98, 104052 (2018).

116. Steinhoff, J., Hinderer, T., Dietrich, T. & Foucart, F. Spin effects on neutron star fundamental-mode dynamical tides: phenomenology and comparison to numerical simulations. Phys. Rev. Res. 3, 033129 (2021).

117. Kawaguchi, K. et al. Frequency-domain gravitational waveform models for inspiraling binary neutron stars. Phys. Rev. D97, 044044 (2018).

118. Dietrich, T., Bernuzzi, S. & Tichy, W. Closed-form tidal approx-
imants for binary neutron star gravitational waveforms con-
structed from high-resolution numerical relativity simulations.
Phys. Rev. D96, 121501 (2017).

119. Dietrich, T. et al. Matter imprints in waveform models for neutron star binaries: tidal and self-spin effects. Phys. Rev. D99, 024029 (2019).

120. Dietrich, T. et al. Improving the NRTidal model for binary neutron star systems. Phys. Rev. D100, 044003 (2019).

121. Thompson, J. E. et al. Modeling the gravitational wave signature of neutron star black hole coalescences. Phys. Rev. D. **101**, 124059 (2020).

122. LIGO Scientific Collaboration. LIGO Algorithm Library—LAL-Suite. (2018).

123. Steinhoff, J., Hinderer, T., Buonanno, A. & Taracchini, A. Dynamical tides in general relativity: effective action and effective-one-body Hamiltonian. Phys. Rev. D94, 104028 (2016).

124. Lackey, B. D., Pürrer, M., Taracchini, A. & Marsat, S. Surrogate model for an aligned-spin effective one body waveform model of binary neutron star inspirals using Gaussian process regression. Phys. Rev. D100, 024002 (2019).

125. Nagar, A. et al. Nonlinear-in-spin effects in effective-one-body waveform models of spin-aligned, inspiralling, neutron star binaries. Phys. Rev. D 99, 044007 (2019).

126. Nagar, A. & Rettegno, P. Efficient effective one body time-domain gravitational waveforms. Phys. Rev. D99, 021501 (2019).
127. Matas, A. et al. Aligned-spin neutron-star–black-hole waveform model based on the effective-one-body approach and numerical-relativity simulations. *Phys. Rev. D* **102**, 043023 (2020).

128. Veitch, J. et al. Parameter estimation for compact binaries with ground-based gravitational-wave observations using the LALInference software library. *Phys. Rev.* **D91**, 042003 (2015).

129. Bulla, M. POSSIS: predicting spectra, light curves and polarization for multi-dimensional models of supernovae and kilonovae. *Mon. Not. R. Astron. Soc.* **489**, 5037–5045 (2019).

130. Mazzali, P. A. & Lucy, L. B. The application of Monte Carlo methods to the synthesis of early-time supernovae spectra. *Astron. Astrophys.* **279**, 447 (1993).

131. Magee, M. R., Sim, S. A., Kotak, R. & Kerzendorf, W. E. Modelling the early time behaviour of type Ia supernovae: effects of the $^{56}$Ni distribution. *Astron. Astrophys.* **614**, A115 (2018).

132. Tanaka, M., Kato, D., Gaigalas, G. & Kawaguchi, K. Systematic opacity calculations for kilonovae. *Mon. Not. R. Astron. Soc.* **496**, 1369–1392 (2020).

133. Bulla, M., Sim, S. A. & Kromer, M. Polarization spectral synthesis for Type Ia supernova explosion models. *Mon. Not. R. Astron. Soc.* **450**, 967–981 (2015).

134. Barbary, K. et al. Sncosmo. https://doi.org/10.5281/zenodo.6363879 (2022).

135. Coughlin, M. et al. Toward rapid transient identification and characterization of kilonovae. *Astrophys. J.* **849**, 12 (2017).

136. Sari, R., Piran, T. & Narayan, R. Spectra and light curves of gamma-ray burst afterglows. *Astrophys. J. Lett.* **497**, L17 (1998).

137. Eichler, D., Livio, M., Piran, T. & Schramm, D. N. Nucleosynthesis, neutrino bursts and gamma-rays from coalescing neutron stars. *Nature* **340**, 126–128 (1989).

138. Paczynski, B. Cosmological gamma-ray bursts. *Acta Astron.* **41**, 257–267 (1991).

139. Meszaros, P. & Rees, M. J. Tidal heating and mass loss in neutron star binaries—implications for gamma-ray burst models. *Astrophys. J.* **397**, 570–575 (1992).

140. Narayan, R., Paczynski, B. & Piran, T. Gamma-ray bursts as the death throes of massive binary stars. *Astrophys. J.* **395**, L83–L86 (1992).

141. van Eerten, H., Leventis, K., Meliani, Z., Wijers, R. & Keppens, R. Gamma-ray burst afterglows from trans-relativistic blast wave simulations. *Mon. Not. R. Astron. Soc.* **403**, 300 (2010).

142. Granot, J. & Piran, T. On the lateral expansion of gamma-ray burst jets. *Mont. Not. R. Astron. Soc.* **421**, 570–587 (2012).

143. Krüger, C. J. & Foucart, F. Estimates for disk and ejecta masses produced in compact binary mergers. *Phys. Rev. D* **101**, 103002 (2020).

144. Dietrich, T. & Ujevic, M. Modeling dynamical ejecta from binary neutron star mergers and implications for electromagnetic counterparts. *Class. Quant. Grav.* **34**, 105014 (2017).

145. Nedora, V. et al. Mapping dynamical ejecta and disk masses from numerical relativity simulations of neutron star mergers. *Class. Quant. Grav.* **39**, 015008 (2022).

146. Agathos, M. et al. Inferring prompt Black-Hole formation in neutron star mergers from gravitational-wave data. *Phys. Rev. D* **101**, 044006 (2020).

147. Vallisneri, M., Kanner, J., Williams, R., Weinstein, A. & Stephens, B. The LIGO Open Science Center. *J. Phys. Conf. Ser.* **610**, 012021 (2015).

148. Raaijmakers, G. et al. Constraints on the dense matter equation of state and neutron star properties from NICER's mass-radius estimate of PSR J0740+6620 and multimessenger observations. *Astrophys. J. Lett.* **918**, L29 (2021).

149. Tews, I., Margueron, J. & Reddy, S. Critical examination of constraints on the equation of state of dense matter obtained from GW170817. *Phys. Rev. C98*, 045804 (2018).

## Acknowledgements

We thank N. Andersson, R. Essick, P. Landry, and J. Margueron for insightful discussions. P.T.H.P. and C.V.D.B. are supported by the research program of the Netherlands Organization for Scientific Research (NWO). T.D. acknowledges support of the Daimler and Benz Foundation. M.W.C. acknowledges support from the National Science Foundation with grant numbers PHY-2308862 and OAC-2117997. M.B. acknowledges support from the Swedish Research Council (Reg. no. 2020-03330). The work of I.T. was supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under contract No. DE-AC52-06NA25396, by the Laboratory Directed Research and Development program of Los Alamos National Laboratory under project number 20220658ER, and by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Scientific Discovery through Advanced Computing (SciDAC) program. J.H. acknowledges support from the National Science Foundation with grant number PHY-1806990. Funded/Co-funded by the European Union (ERC, SMArt, 101076369). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. Computations have been performed on the Minerva HPC cluster of the Max-Planck-Institute for Gravitational Physics and on SuperMUC-NG (LRZ) under project number pn56zo. 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. Resources supporting this work were provided by the Minnesota Supercomputing Institute (MSI) at University of Minnesota under the project "Identification of Variable Objects in the Zwicky Transient Facility," and the Supercomputing Laboratory at King Abdullah University of Science and Technology (KAUST) in Thuwal, Saudi Arabia. This research has made use of data, software and/or web tools obtained from the Gravitational Wave Open Science Center (https://www.gwopenscience.org), a service of LIGO Laboratory, the LIGO Scientific Collaboration and the Virgo Collaboration. This material is based upon work supported by NSF's LIGO Laboratory which is a major facility fully funded by the 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.

## Author contributions

Conceptualization: P.T.H.P., T.D., M.W.C., M.B., I.T., and C.V.D.B.; Methodology: P.T.H.P., T.D., M.W.C., M.B., I.T., M.A., S.A., and C.V.D.B.; Data curation: P.T.H.P., M.W.C., and M.A.; Software: P.T.H.P., T.D., M.W.C., M.B., I.T., M.A., T.B., R.W.K., N.K., G.M., B.R., N.S., A.T., S.A., R.O.V., J.H., P.S., R.Sh., R.So., and C.V.D.B.; Validation: P.T.H.P. and M.W.C.; Formal analysis: P.T.H.P. and M.W.C.; Resources: T.D., M.W.C., M.B., and I.T.; Funding acquisition: T.D., M.W.C., M.B., and I.T.; Project administration: P.T.H.P., T.D., M.W.C., M.B., I.T., and C.V.D.B.; Supervision: P.T.H.P., T.D., M.W.C., M.B., and I.T.; Visualization: P.T.H.P., T.D., M.W.C., M.B., I.T., and N.K.; Writing—original draft: P.T.H.P., T.D., M.W.C., M.B., I.T., N.K., V.N., and P.S.; Writing—review and editing: P.T.H.P., T.D., M.W.C., M.B., I.T., N.K., and V.N.

## Funding

Open Access funding enabled and organized by Projekt DEAL.

## Competing interests

The authors declare no competing interests.
127. Matas, A. et al. Aligned-spin neutron-star–black-hole waveform model based on the effective-one-body approach and numerical-relativity simulations. *Phys. Rev. D* **102**, 043023 (2020).

128. Veitch, J. et al. Parameter estimation for compact binaries with ground-based gravitational-wave observations using the LALInference software library. *Phys. Rev.* **D91**, 042003 (2015).

129. Bulla, M. POSSIS: predicting spectra, light curves and polarization for multi-dimensional models of supernovae and kilonovae. *Mon. Not. R. Astron. Soc.* **489**, 5037–5045 (2019).

130. Mazzali, P. A. & Lucy, L. B. The application of Monte Carlo methods to the synthesis of early-time supernovae spectra. *Astron. Astrophys.* **279**, 447 (1993).

131. Magee, M. R., Sim, S. A., Kotak, R. & Kerzendorf, W. E. Modelling the early time behaviour of type Ia supernovae: effects of the $^{56}$Ni distribution. *Astron. Astrophys.* **614**, A115 (2018).

132. Tanaka, M., Kato, D., Gaigalas, G. & Kawaguchi, K. Systematic opacity calculations for kilonovae. *Mon. Not. R. Astron. Soc.* **496**, 1369–1392 (2020).

133. Bulla, M., Sim, S. A. & Kromer, M. Polarization spectral synthesis for Type Ia supernova explosion models. *Mon. Not. R. Astron. Soc.* **450**, 967–981 (2015).

134. Barbary, K. et al. Sncosmo. https://doi.org/10.5281/zenodo.6363879 (2022).

135. Coughlin, M. et al. Toward rapid transient identification and characterization of kilonovae. *Astrophys. J.* **849**, 12 (2017).

136. Sari, R., Piran, T. & Narayan, R. Spectra and light curves of gamma-ray burst afterglows. *Astrophys. J. Lett.* **497**, L17 (1998).

137. Eichler, D., Livio, M., Piran, T. & Schramm, D. N. Nucleosynthesis, neutrino bursts and gamma-rays from coalescing neutron stars. *Nature* **340**, 126–128 (1989).

138. Paczynski, B. Cosmological gamma-ray bursts. *Acta Astron.* **41**, 257–267 (1991).

139. Meszaros, P. & Rees, M. J. Tidal heating and mass loss in neutron star binaries—implications for gamma-ray burst models. *Astrophys. J.* **397**, 570–575 (1992).

140. Narayan, R., Paczynski, B. & Piran, T. Gamma-ray bursts as the death throes of massive binary stars. *Astrophys. J.* **395**, L83–L86 (1992).

141. van Eerten, H., Leventis, K., Meliani, Z., Wijers, R. & Keppens, R. Gamma-ray burst afterglows from trans-relativistic blast wave simulations. *Mon. Not. R. Astron. Soc.* **403**, 300 (2010).

142. Granot, J. & Piran, T. On the lateral expansion of gamma-ray burst jets. *Mont. Not. R. Astron. Soc.* **421**, 570–587 (2012).

143. Krüger, C. J. & Foucart, F. Estimates for disk and ejecta masses produced in compact binary mergers. *Phys. Rev. D* **101**, 103002 (2020).

144. Dietrich, T. & Ujevic, M. Modeling dynamical ejecta from binary neutron star mergers and implications for electromagnetic counterparts. *Class. Quant. Grav.* **34**, 105014 (2017).

145. Nedora, V. et al. Mapping dynamical ejecta and disk masses from numerical relativity simulations of neutron star mergers. *Class. Quant. Grav.* **39**, 015008 (2022).

146. Agathos, M. et al. Inferring prompt Black-Hole formation in neutron star mergers from gravitational-wave data. *Phys. Rev. D* **101**, 044006 (2020).

147. Vallisneri, M., Kanner, J., Williams, R., Weinstein, A. & Stephens, B. The LIGO Open Science Center. *J. Phys. Conf. Ser.* **610**, 012021 (2015).

148. Raaijmakers, G. et al. Constraints on the dense matter equation of state and neutron star properties from NICER’s mass–radius estimate of PSR J0740+6620 and multimessenger observations. *Astrophys. J. Lett.* **918**, L29 (2021).

149. Tews, I., Margueron, J. & Reddy, S. Critical examination of constraints on the equation of state of dense matter obtained from GW170817. *Phys. Rev. C98*, 045804 (2018).

## Acknowledgements

We thank N. Andersson, R. Essick, P. Landry, and J. Margueron for insightful discussions. P.T.H.P. and C.V.D.B. are supported by the research program of the Netherlands Organization for Scientific Research (NWO). T.D. acknowledges support of the Daimler and Benz Foundation. M.W.C. acknowledges support from the National Science Foundation with grant numbers PHY-2308862 and OAC-2117997. M.B. acknowledges support from the Swedish Research Council (Reg. no. 2020-03330). The work of I.T. was supported by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under contract No. DE-AC52-06NA25396, by the Laboratory Directed Research and Development program of Los Alamos National Laboratory under project number 20220658ER, and by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Scientific Discovery through Advanced Computing (SciDAC) program. J.H. acknowledges support from the National Science Foundation with grant number PHY-1806990. Funded/Co-funded by the European Union (ERC, SMArt, 101076369). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. Computations have been performed on the Minerva HPC cluster of the Max-Planck-Institute for Gravitational Physics and on SuperMUC-NG (LRZ) under project number pn56zo. 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. Resources supporting this work were provided by the Minnesota Supercomputing Institute (MSI) at University of Minnesota under the project “Identification of Variable Objects in the Zwicky Transient Facility,” and the Supercomputing Laboratory at King Abdullah University of Science and Technology (KAUST) in Thuwal, Saudi Arabia. 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. This material is based upon work supported by NSF’s LIGO Laboratory which is a major facility fully funded by the 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.

## Author contributions

Conceptualization: P.T.H.P., T.D., M.W.C., M.B., I.T., and C.V.D.B.; Methodology: P.T.H.P., T.D., M.W.C., M.B., I.T., M.A., S.A., and C.V.D.B.; Data curation: P.T.H.P., M.W.C., and M.A.; Software: P.T.H.P., T.D., M.W.C., M.B., I.T., M.A., T.B., R.W.K., N.K., G.M., B.R., N.S., A.T., S.A., R.O.V., J.H., P.S., R.Sh., R.So., and C.V.D.B.; Validation: P.T.H.P. and M.W.C.; Formal analysis: P.T.H.P. and M.W.C.; Resources: T.D., M.W.C., M.B., and I.T.; Funding acquisition: T.D., M.W.C., M.B., and I.T.; Project administration: P.T.H.P., T.D., M.W.C., M.B., I.T., and C.V.D.B.; Supervision: P.T.H.P., T.D., M.W.C., M.B., and I.T.; Visualization: P.T.H.P., T.D., M.W.C., M.B., I.T., and N.K.; Writing—original draft: P.T.H.P., T.D., M.W.C., M.B., I.T., N.K., V.N., and P.S.; Writing—review and editing: P.T.H.P., T.D., M.W.C., M.B., I.T., N.K., and V.N.

## Funding

Open Access funding enabled and organized by Projekt DEAL.

## Competing interests

The authors declare no competing interests.
**Additional information**

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-023-43932-6.

**Correspondence** and requests for materials should be addressed to Tim Dietrich.

**Peer review information** *Nature Communications* thanks the anonymous reviewers for their contribution to the peer review of this work. A peer review file is available.

**Reprints and permissions information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2023

<footnote n="1">¹Nikhef, Science Park 105, 1098 XG Amsterdam, The Netherlands. ²Institute for Gravitational and Subatomic Physics (GRASP), Utrecht University, Princetonplein 1, 3584 CC Utrecht, The Netherlands. ³Institut für Physik und Astronomie, Universität Potsdam, Haus 28, Karl-Liebknecht-Str. 24/25, 14476 Potsdam, Germany. ⁴Max Planck Institute for Gravitational Physics (Albert Einstein Institute), Am Mühlenberg 1, 14476 Potsdam, Germany. ⁵School of Physics and Astronomy, University of Minnesota, Minneapolis, MN 55455, USA. ⁶The Oskar Klein Centre, Department of Astronomy, Stockholm University, AlbaNova, SE-106 91 Stockholm, Sweden. ⁷Department of Physics and Earth Science, University of Ferrara, Via Saragat 1, I-44122 Ferrara, Italy. ⁸INFN, Sezione di Ferrara, Via Saragat 1, I-44122 Ferrara, Italy. ⁹INAF, Osservatorio Astronomico d'Abruzzo, Via Mentore Maggini snc, 64100 Teramo, Italy. ¹⁰Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA. ¹¹Department of Physics, American University of Sharjah, PO Box 26666 Sharjah, UAE. ¹²Laboratoire de Physique et de Chimie de l'Environnement, Université Joseph KI-ZERBO, Ouagadougou, Burkina Faso. ¹³Observatoire de la Côte d'Azur, Université Côte d'Azur, CNRS, 96 Boulevard de l'Observatoire, F06304 Nice Cedex 4, France. ¹⁴Department of Physics, American University, Washington, DC 20016, USA.</footnote>

<footnote n="15">¹⁵Department of Physics and Astronomy, University of Minnesota—Duluth, Duluth, MN 55812, USA. ¹⁶Department of Physics, Drexel University, Philadelphia, PA 19104, USA. ¹⁷Department of Physics, Massachusetts Institute of Technology, 77 Massachusetts Ave, Cambridge, MA 02139, USA. ¹⁸Department of Physics, Deshbandhu College, University of Delhi, New Delhi, India. ¹⁹Université Lyon, Université Claude Bernard Lyon 1, CNRS/IN2P3, IP2I Lyon, UMR 5822, F-69622 Villeurbanne, France. ²⁰Department of Physics, Syracuse University, Syracuse, NY 13244, USA. e-mail: tim.dietrich@uni-potsdam.de</footnote>
**Additional information**

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-023-43932-6.

**Correspondence** and requests for materials should be addressed to Tim Dietrich.

**Peer review information** *Nature Communications* thanks the anonymous reviewers for their contribution to the peer review of this work. A peer review file is available.

**Reprints and permissions information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2023

<footnote n="1">¹Nikhef, Science Park 105, 1098 XG Amsterdam, The Netherlands. ²Institute for Gravitational and Subatomic Physics (GRASP), Utrecht University, Princetonplein 1, 3584 CC Utrecht, The Netherlands. ³Institut für Physik und Astronomie, Universität Potsdam, Haus 28, Karl-Liebknecht-Str. 24/25, 14476 Potsdam, Germany. ⁴Max Planck Institute for Gravitational Physics (Albert Einstein Institute), Am Mühlenberg 1, 14476 Potsdam, Germany. ⁵School of Physics and Astronomy, University of Minnesota, Minneapolis, MN 55455, USA. ⁶The Oskar Klein Centre, Department of Astronomy, Stockholm University, AlbaNova, SE-106 91 Stockholm, Sweden. ⁷Department of Physics and Earth Science, University of Ferrara, Via Saragat 1, I-44122 Ferrara, Italy. ⁸INFN, Sezione di Ferrara, Via Saragat 1, I-44122 Ferrara, Italy. ⁹INAF, Osservatorio Astronomico d'Abruzzo, Via Mentore Maggini snc, 64100 Teramo, Italy. ¹⁰Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA. ¹¹Department of Physics, American University of Sharjah, PO Box 26666 Sharjah, UAE. ¹²Laboratoire de Physique et de Chimie de l'Environnement, Université Joseph KI-ZERBO, Ouagadougou, Burkina Faso. ¹³Observatoire de la Côte d'Azur, Université Côte d'Azur, CNRS, 96 Boulevard de l'Observatoire, F06304 Nice Cedex 4, France. ¹⁴Department of Physics, American University, Washington, DC 20016, USA. ¹⁵Department of Physics and Astronomy, University of Minnesota—Duluth, Duluth, MN 55812, USA. ¹⁶Department of Physics, Drexel University, Philadelphia, PA 19104, USA. ¹⁷Department of Physics, Massachusetts Institute of Technology, 77 Massachusetts Ave, Cambridge, MA 02139, USA. ¹⁸Department of Physics, Deshbandhu College, University of Delhi, New Delhi, India. ¹⁹Université Lyon, Université Claude Bernard Lyon 1, CNRS/IN2P3, IP2I Lyon, UMR 5822, F-69622 Villeurbanne, France. ²⁰Department of Physics, Syracuse University, Syracuse, NY 13244, USA. ✉e-mail: tim.dietrich@uni-potsdam.de</footnote>