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

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⊙ neutron star to be \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R=11.9{8}_{-0.40}^{+0.35}$$\end{document}R=11.98−0.40+0.35 km.

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 results.Such frameworks and attempts include, among others, the work of Breschi et al. 30 performing Bayesian inference and model selection on the kilonova AT2017gfo, Nicholl et al. 31 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. 32 .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 22,33 , to investigate the nature of the compact binary merger GW190814 34 , to provide techniques to search for kilonova transients 35 , to classify observed EM transients such as GRB200826A 36 , and to combine information from multimessenger observations with data from nuclear-physics experiments such as heavy-ion collisions 23 .
Here, we upgrade our framework to allow for a simultaneous analysis of kilonova, GRB afterglow, and GW data capitalizing on the multimessenger 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) [39][40][41][42][43] or data extracted from Figure 1 | Overview of constraints on the EOS from different information channels.We show a set of possible EOSs (blue lines) that are constrained up to 1.5n sat by Quantum Monte Carlo calculations using chiral EFT interactions 37 and extended to higher densities using a speed of sound model 38 .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.
nuclear-physics experiments, e.g, heavy-ion collisions 44 or the recent PREX-II experiment at Jefferson Laboratory 45 , provide valuable input up to about twice the nuclear saturation density, nsat ≈ 0.16 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 nsat, 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⊙, such as PSR J0348+0432 46 , PSR J1614-2230 47 , and PSR J0740+6620 48 , 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 49 and might be used to constrain the NS EOS 50 .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) 51 in a frequency range of 20Hz to 2048Hz, covering the detected BNS inspiral 52 .For the EM signal, we use the data  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 analysing 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.
set compiled in Coughlin et al. 53 , 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 analysed 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 magnitude uncertainty for the individual lightcurves.This one magnitude uncertainty is introduced during the inference and should account for systematic uncertainties in kilonova modelling.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. 54 .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 54 , assumed heating rates, thermalization efficiencies and opacities within the ejecta 55 .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 af- 409.57+99.87

GW170817 & AT2017gfo GW170817 & GRB170817A & AT2017gfo
Figure 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 ej dyn , the mass of the disk wind ejecta m ej wind , 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 Λ, and the radius of a 1.4 solar mass neutron star R 1.4 at 68%, 95% and 99% confidence.For the 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).
terglow.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.
Fig. 3 summarizes our main findings and shows joint posteriors for the mass of the dynamical ejecta m ej dyn , the mass of the disk wind ejecta m ej wind , the chirp mass Mc, the mass ratio q, the mass-weighted tidal deformability Λ, and the radius of a 1.4 solar mass neutron star R1.4.In contrast to previous findings using simpler kilonova modelling (see Ref. 56 and references therein), we can fit AT2017gfo with masses for the dynamical (about 0.006 M⊙) and disk-wind (about 0.07 M⊙) ejecta components that are within the range of values predicted by numerical-relativity simulations 57 .While the parameters extracted are consistent with our previous findings 22 , we observe a clear improvement on the parameter error bounds due to (i) performing a simultaneous analysis of the distinct messengers and (ii) employing a modified likelihood function when analysing the kilonova.For instance, the constraints on R1.4 = 11.86 +0.41 −0.53 , a typical choice to quantify EOS constraints, is significantly improved compared to our previous result, R1.4 = 11.75 +0.86 −0.81 km 22 .The half-width of R1.4's 90% credible interval decreases from about 800m 22 to about 400m.

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. 58 .The GRB afterglow light-curve data are analyzed with the synthetic Gaussian jet-model lightcurve described before 59,60 .Fig. 2 shows the corresponding best-fit lightcurve for the kilonova with a 1 magnitude uncertainty band as before.Moreover, we are also presenting the bestfit 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.

Flux density [mJy]
Figure 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 analysing GW170817, AT2017gfo, GRB170817A.We compare our model predictions with the observational data including the 1-sigma measurement uncertainty.
Fig. 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 GW170817and-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 E0 as the ratio ξ between it, E0, and the disk mass m 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. 61 ), which is correlated with the ejecta masses (e.g.Fig. 3 in Ref. 62 ).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.Tab. 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.

Equation of State construction
The EOS describes the relation between energy density ε, 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 ε and p.
The most general constraints on the EOS can be inferred from the slope of the EOS, the speed of sound, defined as 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, cS ≤ c.Furthermore, the speed of sound in a NS has to be larger than zero, cS ≥ 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 66,67 , i.e., the outermost layer of NSs at densities below approximately 0.5nsat.Above roughly 0.5nsat, 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 = np/n with np being the proton density.From this, the EOS follows from and The proton fraction x(n) is then determined from the beta equilibrium condition, µn = µp + µe, where µ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 manybody 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 68 .At densities of the order of nsat, however, the effective d.o.f. are nucleons, neutrons and protons, that can be treated as point-like nonrelativistic particles.Then, the nuclear Hamiltonian can be written generically as where T denotes the kinetic energy of the nucleons, V NN ij describes two-nucleon (NN) interactions between nucleons i and j, and V 3N ijk 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 69 .
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) 70,71 .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 yield 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 Λ b , is of the order of 500 − 600 MeV 72 and interaction contributions involving heavier d.o.f., such as the ρ meson, are integrated out.The chiral Lagrangian is then expanded in powers of q/Λ b according to a power-counting scheme.Most current chiral EFT interactions are derived in Weinberg power counting 70,71,[73][74][75] .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 72,76,77 and that can be systematically improved by increasing the order of the calculation.Chiral EFT Hamiltonian naturally include 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 39,42,78 or fourth order (N 3 LO) 41,79 .
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) 40,41,79 , the self-consistent Green's function (SCGF) method 80 , or the coupled-cluster (CC) method 78,81 .Here, we employ Quantum Monte Carlo (QMC) methods 82 , which provide nonperturbative 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 noninteracting 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 82 , and Auxiliary Field Diffusion Monte Carlo (AFDMC) 83 , 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 β equilibrium using phenomenological information on symmetric nuclear matter and constructing a consistent crust reflecting the uncertainties of the calculations 84 .This crust includes a description of the outer crust 66 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 manybody 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 2nsat 37,72 .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, cS ≤ c, and mechanical stability, cS ≥ 0, mentioned before.There exist several such extension schemes in literature: parametric ones, like the polytropic expansion [85][86][87] or expansions in the speed of sound 88,89 , and nonparametric approaches 90 .To extend the AFDMC calculations employed here, we employ a parametric speed-of-sound extension scheme.Working in the cS versus n plane, the speed of sound cS(n) is determined with theoretical uncetainty 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 ≈ 12nsat, well beyond the regime realized in NSs.For each grid point, we sample random values for c 2 s (ni) 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 2 s,i (ni) 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 where µ(n) is the chemical potential that can be obtained from the speed of sound using the relation For each reconstructed EOS, constrained by Chiral EFT at low densities and extrapolated via the cS 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 (Λ) 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 highdensity phases of QCD can be included straightforwardly.For example, pQCD calculations at asymptotically high densities 49 , of the order of 40 − 50nsat, might be used to constrain the general EOS extension schemes even further 50,87 .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 [91][92][93][94] , which describes the observed behavior of the speed of sound in NSs 37 : a rise of the speed of sound at low densities to values above the conformal limit of c/ √ 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 (a) an exploration of theoretical uncertainties in the low-density EOS and (b) constructs the high-density EOS using a cS extrapolation.(a) 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 2 LO 37 with theoretical uncertainties are provided.Upon initiation of the class, a sample is drawn from the low-density uncertainty band using a 1parameter sampling technique.In this approach, a uniform random number ω is sampled uniformly between 0 and 1, and the interpolated EOS is given as 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. 95 does not change our results appreciably.In future, we will explore additional schemes, e.g., using Gaussian processes 63 .(b) 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 2nsat.Beyond this density, the class constructs the EOS using a cS 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 12nsat for the former and 5 line segments for the latter.The code then solves Eqs. ( 5), ( 6) and ( 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 46 , and PSR J1614-2230 47 provide a lower bound on the maximum mass of a NS.
The likelihood for a mass-only measurement is given by where P(M |PSR) is the posterior distribution of the pulsar's mass and MTOV is the maximum mass supported by the EOS with parameters 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 65,[96][97][98][99] .The corresponding likelihood is given by where PNICER(M, R) is the joint-posterior distribution of mass and radius as measured by NICER and we use the fact that (i) the radius is a function of mass for a given EOS, and (ii) 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 with the amplitude A(f ) and the GW phase ψ(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 100,101 and, once measured, allow to place constraints on the EOS governing the NS interior [102][103][104][105] .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 with the individual tidal deformabilities Λ1,2 = 2 3 k 1,2 2 /C 5 1,2 and the individual masses m1,2.Here, k 1,2 2 are the Love numbers describing the static quadrupole deformation of one body inside the gravitoelectric field of the companion and C1,2 are the individual compactnesses C1,2 = m1,2/R1,2 in isolation.

GW analysis
By assuming stationary Gaussian noise, the GW likelihood LGW(θ) that the data d is a sum of noise and a GW signal h with parameters θ is given by 130 LGW ∝ exp − where the inner product ⟨a|b⟩ is defined as Here, ã(f ) is the Fourier transform of a(t), * denotes complex conjugation, and Sn(f ) is the one-sided power spectral density of the noise.The choice of f low and f high depends on the type of binary that we are interested in.In our study, we will set flow and fhigh 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.

Kilonova models
Kilonova models are extracted using the 3D Monte Carlo radiative transfer code POSSIS 131 .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 ρ and compositions (i.e., electron fraction Ye), 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 132,133 .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 κ λ (ρ, T, Ye, t) from Tanaka et al. 134 are implemented in the code and depend on the local properties of the ejecta (ρ, T , and Ye).Spectral time series are extracted using the technique described by Bulla et al. 135 and used to construct broad-band lightcurves in any desired filter.

Supernova models
Templates available within the SNCosmo library 136 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 j i (ti).We assume these measurements to be given as a time series at times ti with a corresponding statistical error σ j i ≡ σ j (ti).The likelihood function LEM(θ) then reads 137 : where m j,est i (θ) is the estimated AB magnitude for the parameters θ and σsys is the additional error budget for accounting the systematic uncertainty within the electromagnetic signal modeling.The inclusion of σsys is equivalent to adding a shift of ∆m to the light curve, for which marginalized with respect to a zero-mean normal distribution with a variance of σ 2 sys .This likelihood is equivalent to approximating the probability distribution of the spectral flux density fν 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ν ∈ (0, ∞).There are two advantages of approximating fν with a Log-normal distribution: (i) if the uncertainty is larger or comparable to the measured value, it avoids having non-zero support for the nonphysical fν < 0; (ii) if the uncertainty is much smaller than the measured value, the Lognormal 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 ej dyn , the disk wind ejecta mass m ej wind , the half-opening angle of the lanthanide-rich component Φ, and the viewing angle θ obs .Both the dynamical ejecta and the disk wind ejecta are described in methods section 4.

GRB Afterglows
In our framework, the computation of the GRB afterglow lightcurves is until now based on the publicly available semi-analytic code afterglowpy 59,60 .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, ϵe and ϵB, denoting the fraction of the shock energy that goes into the relativistic electrons and magnetic field, respectively, and p, and the slope of the electron energy distribution dn/dγ ∝ γ −p , with n being the electron number density and γ being the electron Lorentz factor.The flux density of the curvature radiation is where τ is the optical depth and ϵν and αν are the impassivity coefficient and absorption coefficient, respectively.For a fixed power-law distribution of electrons these can be approximated analytically 138 .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 ∝ E0 exp(− 1 2 ( θ θc ) 2 ), where θc characterizes the width of the Gaussian.The jet truncation angle is θw.We assume the GRB jet to be powered by the accretion of mass from the disk onto the remnant black hole [139][140][141][142] .Consequently, the jet energy is proportional to the leftover disk mass, where ξ is the fraction of disk mass ejected as wind and ϵ is the fraction of residual disk mass converted into jet energy.The dynamical evolution of these layers is computed semianalytically using the "thin-shell approximation" casting energyconservation 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 143,144 .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ν , is obtained by equal arrival time surface integration, Eq. ( 17), taking into account relativistic effects, i.e., that the observed Fν 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. 145 and Dietrich et al. 22 but emphasize that a variety of other fitting formulas exist in the literature 20,53,57,146,147 .
In NMMA, the dynamical ejecta mass m ej dyn is connected to the binary properties through the phenomenological relation 145 m ej dyn,fit where mi and Ci 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.004M⊙ 145 .Therefore, the dynamical ejecta mass can be approximated as where α ∼ N (µ = 0, σ = 0.004M⊙).
To determine the disk mass m disk , we follow the description of Dietrich et al. 22 , with a and b given by where ao, bo, δa, δb, c, and d are free parameters.The parameter ∆ is given by where q ≡ m2/m1 ≤ 1 is the mass ratio and where MTOV and R1.6 are the maximum mass of a non-spinning NS and the radius of a 1.6M⊙ NS.We note that we assume that the disk-wind ejecta component is proportional to the disk mass, i.e., m ej wind = ξ × m disk .

Bayesian statistics
Based on Bayes' theorem, the posterior distribution of the parameters p(θ|d, H) under hypothesis H with data d is given by where L(θ), π(θ), and 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 χ 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 Z carries additional information on the plausibility of a given hypothesis H.The evidence is given by which is the normalization constant for the posterior distribution.Moreover, we can compare the plausibilities of two hypotheses, H1 and H2, by using the odd ratio O 1 2 , which is given by where B 1 2 and Π 1 2 are the Bayes factor and prior odds, respectively.If O 1 2 > 1, H1 is more plausible than H2, 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.
Scaling.Due to the enormous computational burden of simultaneously analysing both the EM and GW signals, we have to ensure a good parallelization of our code.For the NMMA framework, the parallelization is achieved by taking advantage of the flexibility of dynesty 156 and based on the interface presented in parallel bilby 157 .
In particular, the nested sampling process is parallelized using a head/worker strategy.The "head" organizes the live/dead points, and validates if the stopping criteria are reached, while the "workers" find new live points under the likelihood constraint.The theoretical scaling performance of such a strategy is given by 158 where T (Ncores) is the runtime using Ncores cores and T (Ncore ref ) is the reference time using Ncore ref cores.The parameter n live denotes the number of live points.
To validate if such a scaling is achieved, we performed intensive scaling tests on SuperMUC NG at the Leibniz Supercomputing Centre (Munich), Lise and Emmy of the North German Supercomputing Alliance, and on HAWK of the High-Performance Computing Center Stuttgart.We show results for scaling tests performed on HAWK using AMD EPYC 7742 processors.The tests are based on a full joint inference of GW170817-and-AT2017gfo and GW170817-and-GRB170817A-and-AT2017gfo which includes intermediate checkpointing and all necessary I/O-operations.The strong scaling for such simulations is shown in Supplementary Fig. 5 and is compared to the theoretical scaling mentioned of Eq. (35).
Modelling Uncertainty.Overall, the obtained multi-messenger constraints depend noticeably on the robustness and accuracy of the individual GW, kilonova, and GRB afterglow models but also on the uncertainty and accuracy of the underlying equations of state and nuclear physics computations.In the past, there have been numerous studies considering the effect of GW-model uncertainties, e.g., 155,159,160 and we have already employed multiple GW models in one of our previous studies also to investigate uncertainties with the conclusion that with the current amount of observational data, we are mainly limited by statistical and not systematic uncertainties 22 .As an example of the influence of the particular choice of the kilonova model, we employ two different kilonova models with different assumptions about the geometry and composition of the ejecta.In particular, we compare the results for our standard kilonova model as presented in Supplementary Tab. 2 with the results based on the model of Kasen et al. 154 in Supplementary Tab. 3.While we find that individual parameters can be different, we find a similar neutron star radius of R1.4 = 11.86 +0.37 −0.49km.Effect of systematic uncertainty budget To investigate the effect of the value chosen for the systematics uncertainty budget, the analysis on the GW170817-AT2017gfo-GRB170817A event has been conducted with three different values of σsys, namely, 0.5mag, 1mag and 2mag.
To gauge the effect on the final equation-of-state constraint, the resulting posterior of the radius of a 1.4M⊙ neutron star, R1.4, for these different σsys are compared.The median values with the 90% credible interval as uncertainties, are shown in Supplementary Tab. 4. The resulting radius measurements of a 1.4M⊙ neutron star for σsys being 0.5, 1.0 and 2.0 are shown.
With lower values of σsys, the estimated R1.4 is skewed towards higher values, with the minimal uncertainty at 1mag.Therefore, it hints that • the information from difference channels are more coherent at 1mag as compare to 0.5mag • the information is not over diluted in 1mag as compare to 2mag.
To verify the above claim, we investigate into the estimated best-fit lightcurves.In particular, the best-fit lightcurves with the associated error budget for filter K is shown at Supplementary Fig. 6.One can see that the 0.5mag is underestimating the systematic uncertainty and failed to account for the data after 6 days after the merger, while the 2mag is overcompensating for the systematic uncertainty.
Based on the above observations, we concluded that the σsys of 1mag is a sensible choice.

Figure 6 |
Figure6| Best-fit early-time lightcurve from the analysis with various σ sys .The best-fit K-filter lightcurve (dashed, with the one magnitude uncertainty shown as the band) for AT2017gfo data when analysing GW170817-and-AT2017gfo-and-GRB170817A simultaneously with different σ sys .

Table 1 |
Comparison of radius measurements of a 1.4M ⊙ neutron star for a selection of multi-messenger studies.

Table 4 |
Comparison of radius measurements of a 1.4M ⊙ neutron star for different σ sys .