Rapid model comparison of equations of state from gravitational wave observation of binary neutron star coalescences

The discovery of the coalescence of binary neutron star GW170817 was a watershed moment in the field of gravitational wave astronomy. Among the rich variety of information that we were able to uncover from this discovery was the first non-electromagnetic measurement of the neutron star radius, and the cold nuclear equation of state. It also led to a large equation of state model selection study from gravitational-wave data. In those studies Bayesian nested sampling runs were conducted for each candidate equation of state model to compute their evidence in the gravitational-wave data. Such studies, though invaluable, are computationally expensive and require repeated, redundant, computation for any new models. We present a novel technique to conduct model selection of equation of state in an extremely rapid fashion ($\sim$ minutes) on any arbitrary model. We test this technique against the results of a nested-sampling model selection technique published earlier by the LIGO/Virgo collaboration, and show that the results are in good agreement with a median fractional error in Bayes factor of about 10%, where we assume that the true Bayes factor is calculated in the aforementioned nested sampling runs. We found that the highest fractional error occurs for equation of state models that have very little support in the posterior distribution, thus resulting in large statistical uncertainty. We then used this method to combine multiple binary neutron star mergers to compute a joint-Bayes factor between equation of state models. This is achieved by stacking the evidence of the individual events and computing the Bayes factor from these stacked evidences for each pairs of equation of state.

Figure 1

Comparison of Bayes factors obtained with our approximate Bayes factor calculation scheme and the LALInference_nest nested sampling results [31], shown here for the narrow prior in the top panel and the broad prior in the bottom panel. The Bayes factors are computed with respect to the SLY equation of state model in both cases. We show here results for the TaylorF2 waveform. The error-bars for the approximate Bayes factors are an estimate based on standard deviation (see main text). The error-bars in the nested sampling method is obtained from [31].

Figure 2

Magnitude of deviation of Bayes factors obtained with our approximate evidence calculation scheme from the LALInference_nest nested sampling results, shown here for the narrow prior in the top panel and the broad prior in the bottom panel. The Bayes factors are again computed with respect to the SLY equation of state model in both cases. We show here results for the TaylorF2 waveform. The deviation is the highest for the stiffest equation of state which has the least evidence in the data. The method is more reliable when the model has good support from the data. The median deviation for the narrow (broad) prior is $10\%(11\%)$.

Figure 3

Left: narrow prior: posterior distribution in $(\tilde{\mathrm{\Lambda }},q)$ and its KDE. For comparison we show the equation of state curves for various models in $(\tilde{\mathrm{\Lambda }},q)$. Note that the models that gave the highest deviation with respect to the nested sampling results in Fig. 2 are also the most distant from the peak of the posterior distribution. Right: broad prior case: we see the identical relationship that for the equation of state models that gave the largest deviation in the Bayes factor computation with respect to the nested sampling results, the posterior support is the weakest. Note that the kinks in the various equation of state models in this plot are due to the fact that it extends to smaller values of $q$, where one of the objects in the binary becomes more massive than the maximum allowed NS mass. At this point we consider that object a black hole and set the value of ${\mathrm{\Lambda }}_1=0$. This leads to a sudden change in the value of $\tilde{\mathrm{\Lambda }}$ and hence creates these kinks.

Figure 4

Bayes factors computed for the various equation of state models for the injected signals at 40 Mpc, 70 Mpc, 100 Mpc, 130 Mpc, and 160 Mpc. The blue bars shows the joint Bayes factors computed using Eq. (14), which is largely determined by the strongest source. A Bayes factor threshold of ${10}^{-3}$ is applied to exclude cases that have very small posterior support.

Figure 5

Bayes factor of various equation of state models with respect to SLY model for GW170817, GW190425, and their combination by evidence stacking. Note that the blue bars (joint Bayes factor) and the orange bars (Bayes factor from GW170817 data) are very similar to each other indicating that most of the information in discerning between the different models comes from the data of GW170817. The green bars for GW190425 are adding very little information as can also be seen from the fact that their variation in height across the various models is very small.

Figure 6

Change in Bayes factor due to inclusion of GW190425 for different equation of state models. The area of the bubbles are proportional to the respective joint Bayes factor values. We note that the models for which the Bayes factor was higher for GW170817 tend to an increased Bayes factor upon inclusion of the GW190425 data, whereas models which have lower Bayes factor for GW170817 diminish further upon inclusion of the GW190425 data. However, it should be noted that this effect is very small, given the weakness of the signal from GW190425, as evidenced by the consistency between the bars of Bayes factors of GW170817 and the joint Bayes factor in Fig. 5.