Bayesian model-selection of neutron star equation of state using multi-messenger observations

Measurement of macroscopic properties of neutron stars, whether in binary or in an isolated system, provides us a key opportunity to place a stringent constraint on its equation of state. In this paper, we perform Bayesian model-selection on a wide variety of neutron star equation of state using multi-messenger observations. In particular, (i) we use the mass and tidal deformability measurement from two binary neutron star merger event, GW170817 and GW190425; (ii) simultaneous mass-radius measurement of PSR J0030+0451 and PSR J0740+6620 by NICER collaboration, while the latter has been analyzed by joint NICER/radio/XMM-Newton collaboration. Among the 31 equations of state considered in this analysis, we are able to rule out different variants of MS1 family, SKI5, H4, and WFF1 EoSs decisively, which are either extremely stiff or soft equations of state. The most preferred equation of state model turns out to be AP3 (or MPA1), which predicts the radius and dimensionless tidal deformability of a $1.4 M_{\odot}$ neutron star to be 12.10 (12.50) km and 393 (513) respectively.


INTRODUCTION
Neutron stars (NSs) are the extremely dense objects known in our universe. Properties of matter inside NS are encoded in its equation of state (EoS), which has wide-ranging uncertainty from the theoretical perspective (Lattimer & Prakash 2016;Oertel et al. 2017;Baym et al. 2018a). With the current understanding of quantum chromodynamics, it is very hard to determine the interactions of NS matter at such high densities. Also performing many-body calculations is computationally intractable. Besides, the constituent of NS at its core is highly speculative -perhaps containing exotic matter like quark matter, strange baryons, meson condensates etc. (Glendenning 1997). Even though the matter inside the NS is extremely dense, nevertheless the temperature of this object is actually cold for most of its life span. Such highly dense but rather cold matter cannot be produced in the laboratory. Since probing the physics of NS matter is inaccessible by our earth based experiments, we look for astrophysical observations of NS. These observations can give us the measurement of macroscopic properties related to a NS, like mass, radius, tidal deformability (Hinderer 2008;Binnington & Poisson 2009;Damour & Nagar 2009) etc., which will depend on the internal structure of the NS. The measurement of these macroscopic properties can be used to infer the EoS of NS.
The past few years have been a golden era in NS physics. A number of key astrophysical observations of NS have been made from multiple cosmic messengers such as radio observations of massive pulsars (Demorest et al. 2010;Antoniadis et al. 2013;Cromartie et al. 2019;Fonseca et al. 2021) or mass-radius measurement of PSR J0030+0451 (Miller et al. 2019b;Riley et al. 2019) and PSR J0740+6620 (Miller et al. 2021;Riley et al. 2021) using pulse-profile modeling by NICER collaboration (Gendreau et al. 2016), and also the observation of gravitational waves (GW) from two coalescing binary neutron star (BNS) merger events (Abbott et al. 2017(Abbott et al. , 2018(Abbott et al. , 2020a by LIGO/Virgo collaboration (Aasi et al. 2015;Acernese et al. 2015). These observations have motivated several authors to place joint constraints on the properties of NS (Raaijmakers et al. 2020;Capano et al. 2019;Landry et al. 2020;Jiang et al. 2020;Traversi et al. 2020;Biswas et al. 2020;Al-Mamun et al. 2021;Dietrich et al. 2020;Biswas 2021;Breschi et al. 2021;Miller et al. 2021;Raaijmakers et al. 2021;Pang et al. 2021) either using a phenomenological or nuclearphysics motivated EoS parameterization by employing Bayesian parameter estimation. However, one can take a complementary approach instead of constructing an EoS parameterization and ask the following question: given a variety of NS EoS models in the literature based on nuclear-physics, which one is the most preferred by the current observations in a statistical sense?
A few studies (Abbott et al. 2020b;Ghosh et al. 2021;Pacilio et al. 2021) already exist in the literature on Bayesian model-selection of NS EoS. But they only use GW observations and do not provide the current status of various NS EoS models. Therefore, the aim of this paper is to perform a Bayesian model-selection study amongst various nuclear-physics motivated EoS models of NS using the constraints coming from multimessenger astronomy.

BAYESIAN METHODOLOGY
To perform Bayesian model-selection among several EoSs, we need to compute the Bayesian evidence for each model combining astrophysical data from multiple messengers. The Bayesian methodology used in this work is primarily based on the following Refs. (Del Pozzo et al. 2013;Landry et al. 2020;Biswas et al. 2020), which is described here.
For any two given EoSs, the relative odds ratio between them can be computed as where d = (d GW , d X−ray , d Radio ) is the set of data from the three different types of astrophysical observations. Now using the Bayes' theorem we find, where we assume independence between different sets of data and P (EoS i,j ) is the prior on the EoS i,j before any measurement has taken place. Here we assume each model is equally likely and set the ratio of priors between two models to 1. Therefore, the main quantity of interest is the Bayes factor (B i j ) between a candidate EoS i and EoS j . When B i j is substantially positive, it implies that the data prefers EoS i over EoS j .
For GW observations, information about EoS parameters come from the masses m 1 , m 2 of the two binary components and the corresponding tidal deformabilities Λ 1 , Λ 2 . In this case, where P (m 1 , m 2 |EoS) is the prior distribution over the component masses which should be informed by the NS population model. However, the choice of wrong population model starts to bias the results significantly only after ∼ 20-30 observations (Agathos et al. 2015;Wysocki et al. 2020;Landry et al. 2020). Therefore given the small number of detections at present, this can be fixed by a simple flat distribution over the masses, In our calculation we set M min = 1M and M max to the maximum mass for that particular EoS. Here the normalization factor on the NS mass prior is very important as it prefers the EoS with slightly larger M max than the heaviest observed NS mass and disfavor EoS with much larger M max . For example, two EoSs with M max = 2.5M and M max = 3M , will have mass prior probability P (m|EoS) = 2/3 and P (m|EoS) = 1/2 respectively if M min ≤ m ≤ M max . Though both EoSs support the heaviest NS mass measurement (2.08 ± 0.07M ) equally well, EoS with M max = 3M is less probable than EoS with M max = 2.5M . Similar approach has been employed in previous works as well (Miller et al. 2019a;Raaijmakers et al. 2020;Landry et al. 2020;Biswas 2021). Alternatively, one can truncate the NS mass distribution to a largest population mass which is informed a formation channel (eg. supernova)-in that situation EoSs with M max greater than the largest population mass will be assigned equal probability. Given our lack of knowledge on the upper limit of NS mass distribution, we choose to limit M max based on EoS itself not the formation channel. A broader discussion on different choice of mass prior can be found in the appendix of Ref. (Legred et al. 2021). However, if the masses in GW observation are expected to be smaller than the mass of the heaviest pulsar, then the upper limit on the mass prior can be chosen by the Likelihood's domain of support to reduce the computational time.
Equation 3 can be further simplified by fixing the GW chirp mass to its median value with not so much affecting the result (Raaijmakers et al. 2020) given its high precision measurement. Then we will have one less parameter to integrate over as m 2 will be a deterministic function of m 1 .
X-ray observations give the mass and radius measurements of NS. Therefore, the corresponding evidence takes the following form, Similar to GW observation, here also the explicit prior normalization over the mass should be taken into account or can be chosen by the Likelihood's domain of support (if applicable).
Radio observations provide us with very accurate measurements of the NS mass. In this case, we need to marginalize over the observed mass taking into account its measurement uncertainties, Here the prior normalization of mass must be taken into account as the observed mass measurement is close to the maximum mass predicted by the EoS. To compute these evidences we use the nested sampling algorithm implemented in Pymultinest (Buchner et al. 2014) package.

DATASETS
The likelihood distributions used in this work are modelled as follows: (a). Mass and tidal deformability measurement from GW170817 (Abbott et al. 2019) and GW190425 (Abbott et al. 2020a) are modelled with an optimized multivariate Gaussian kernel density estimator (KDE) implemented in Statsmodels (Seabold & Perktold 2010). (b) Similarly mass and radius measurement of PSR J0030+0451 (Riley et al. 2019;Miller et al. 2019b) and PSR J0740+6620 (Riley et al. 2021;Miller et al. 2021) are also modelled with Gaussian KDE. Since the uncertainty in the mass-radius measurement of PSR J0740+6620 is larger for Ref. (Miller et al. 2021) than Ref. (Riley et al. 2021) due to a conservative treatment of calibration error, we analyze both data separately and provide two Bayes factor values. (c) Mass measurement of PSR J0740+6620 (Cromartie et al. 2019;Fonseca et al. 2021) should be modelled with a Gaussian likelihood of 2.08M mean and 0.07M 1σ standard deviation. However we do not need to use this anymore, as the mass-radius measurement of PSR J0740+6620 already takes it into account. This is to note that further constraint on the properties of NS could be given by the joint detection of GW170817 and its electromagnetic counterparts (Bauswein et al. 2017;Radice et al. 2018;Coughlin et al. 2019;Capano et al. 2019;Breschi et al. 2021). However, this paper intentionally does not include that information as these constraints are rather indirect and need careful modeling of the counterparts. Similarly, constraints coming from heavy ion collision data (Margueron et al. 2018;Margueron & Gulminelli 2019) are not included in this study. For example, recent measurement of neutron skin thickness of Pb 208 (Adhikari et al. 2021) may suggest relatively stiff EoS around the nuclear saturation density, however present constraints are mostly dominated by the astrophysical observations Biswas 2021). Nevertheless, in future it would be interesting to extend this analysis in that direction.
In the past, one common approach has been opted in many papers is to use the bound of radius and tidal deformability of 1.4M NS to rule out EoSs of NS; rather than using the full distribution of the dataset. This can lead to a significant bias in the results (Miller et al. 2019a). Those bounds are also based on either a particular EoS parameterization or EoS insensitive relations (Maselli et al. 2013;Yagi & Yunes 2016, 2017Chatziioannou et al. 2018). Different EoS parameterization leads to different bounds as they do not occupy the same prior volume. In contrast, the results obtained in this study do not depend on any EoS parameterization and therefore, these are also not subjected to any bias due to the model assumption. However, we are also using a specific set of EoSs which restrict us to provide accurate description of preferred EoS model. If the true model is among the candidates, then Bayesian modelselection method will select the correct one. But if all the models are false, then Bayesian model-selection will only select the least incorrect one.

RESULTS
In figure 2, Bayes factor of different EoSs are plotted using multi-messenger observations with respect to the most probable EoS, for which we find the Bayesian evidence to be maximum. In the left panel results are obtained using the data from Ref. (Riley et al. 2019(Riley et al. , 2021 for which AP3 turns out to be the most preferred model. In the right panel the data from Ref. (Miller et al. 2019b(Miller et al. , 2021 are used and in this case, MPA1 is the most preferred EoS. MPA1 EoS has a larger value of R 1.4 , Λ 1.4 , and M max compare to AP3 EoS. Therefore, it is clear the data from Ref. (Miller et al. 2019b(Miller et al. , 2021 prefer stiffer EoS compare to the data from Ref. (Riley et al. 2019(Riley et al. , 2021. In this study, we follow the interpretation of Kass and Raftery (Kass & Raftery 1995) and decisively exclude the EoSs for which log 10 B i AP3 ≤ −2. This region is shown using black shade in the plot. We find SKI5, WFF1, MS1, MS1 PP, MS1B PP, MS1B, and H4 are ruled out for both type of datasets. All of these EoSs except WFF1 are rather stiff EoSs and predict large values of radius and tidal deformability for the NS. This is broadly consistent with GW170817 observations as it mainly favors soft EoS (Abbott et al. 2017(Abbott et al. , 2018. However WFF1 which is the softest EoS considered in this study, is also now decisively ruled out by the multi-messenger observations. In fact WFF1 was found to be the one of the most preferred EoS by the previous studies (Abbott et al. 2020b;Ghosh et al. 2021;Pacilio et al. 2021) based on GW observation only (see also Fig. 3). This demonstrates the true power of multimessenger observations. Now not only stiff EoSs but also extreme soft EoSs are ruled out.
The region between −2 < log 10 B i AP3/MPA1 ≤ −1 is shown in red shade. Only WFF2 falls in this region while using the data from Ref. (Miller et al. 2019b(Miller et al. , 2021 and they have strong evidence against MPA1 according to the interpretation of Kass and Raftery. WFF2 is a relatively softer EoS with the value of R 1.4 = 11.16 km and Λ 1.4 = 232.
log 10 B i AP3/MPA1 > −1 are statistically insignificant. However it is still divided into two regions: −1 < log 10 B i AP3/MPA1 ≤ −1/2 is referred to as substantial evidence against AP3 (shown in yellow band) and log 10 B i AP3/MPA1 ≥ −1/2 is referred to as insubstantial (shown in green band), which is not worth more than a bare mention. Most of our candidate EoSs lie in these two regions; 8(13) of them have substantial evidence and 15(10) of them have insubstantial evidence against AP3(MPA1). The status of each EoS model is also summarized in the last column of Table 1 using the data from (Riley et al. 2019(Riley et al. , 2021. Despite a slightly different radius measurement of PSR J0740+6620 between (Riley et al. 2021) and (Miller et al. 2021), we find broadly consistent Bayesian evidence using both datasets. From the previous studies (Biswas 2021;Legred et al. 2021), it is expected of a maximum ∼ 0.4 km difference in R 1.4 towards the stiff EoSs when the data from from (Miller et al. 2021) is added instead of (Riley et al. 2021). In our study, we also find the similar trend. There is a 0.4 km difference in R 1.4 between the two most preferred EoSs. One way to remove this  systematic error, is to integrate over both datasets using Bayesian framework. Naively one can take equal number of samples from the posterior using both datasets and combined them into a single posterior. However as discussed in (Ashton & Khan 2020) on the systematic error due to modeling uncertainty in waveform for GW signal, giving equal-weighted probability on the both datasets might not be correct way as they corresponds to different Bayesian evidence i.e, one being more likely compare to the other. It is better to weight the samples by their corresponding posterior evidence. Since this paper only focuses on the Bayesian evidence and does not deal with the posterior of any inferred properties, we leave this discussion for a future work (Biswas & Datta 2021 (Riley et al. 2019(Riley et al. , 2021 are used here) based on their Bayes factor value with respect to the most preferred EoS (AP3).

SUMMARY
In summary, we have performed a Bayesian modelselection study over a wide-ranging EoS model using multi-messenger observations of NS. We find EoSs which predict larger radius (R 1.4 ≥ 13.69 km) and tidal deformability (Λ 1.4 ≥ 897 ) are decisively ruled out. As well as WFF1 EoS which predicts a lower radius (R 1.4 = 10.42 km) and tidal deformability (Λ 1.4 = 153 ) is also decisively ruled out. Therefore, EoS of NS cannot be either very stiff or soft. Ironically, the range of R 1.4 from this analysis is very consistent with the prediction (11.5−13.6 km) based on the nuclear physics data by (Li & Steiner 2006 Fig. 2 but only combining two GW observations and the revised mass-measurement of PSR J0740+6620 by Ref. (Fonseca et al. 2021). In this case, WFF2 turns out to be most preferred EoS.
used as a benchmark while making new EoS models or choosing an EoS model to perform any numerical simu-lation of NS. Any future measurement of NS properties from GW and electromagnetic observations can be easily combined using the methodology developed in this paper. ACKNOWLEDGEMENTS I am extremely thankful to Nikolaos Stergioulas for reading this manuscript carefully and making several useful suggestions. I also thank Prasanta Char, Rana Nandi, Wolfgang Kastaun, and Sukanta Bose for their valuable suggestions to improve this manuscript. I gratefully acknowledge the use of high performance super-computing cluster Pegasus at IUCAA for this work. I acknowledge support from the Knut and Alice Wallenberg Foundation under grant Dnr. KAW 2019.0112. This material is based upon work supported by NSF's LIGO Laboratory which is a major facility fully funded by the National Science Foundation.