Rapid hierarchical inference of neutron star equation of state from multiple gravitational wave observations of binary neutron star coalescences

Bayesian hierarchical inference of phenomenological parametrized neutron star equations of state (EoS) from multiple gravitational wave observations of binary neutron star mergers is of fundamental importance in improving our understanding of neutron star structure, the general properties of matter at supra nuclear densities, and the strong nuclear force. However, such an analysis is computationally costly, as it is unable to reuse single-event EoS agnostic parameter estimation runs that are carried out regardless for generating gravitational wave transient catalogs. With the number of events expected to be observable during the fourth observing run (O4) of LIGO/Virgo/KAGRA, this problem can only be expected to worsen. We develop a novel and robust algorithm for rapid and computationally cheap hierarchical inference of parametrized EoS from gravitational wave data which reuses single-event EoS agnostic parameter estimation samples to significantly reduce computational cost. We efficiently include a priori knowledge of neutron star physics as Bayesian priors on the EoS parameters. The high speed and low computational cost of our method allow for efficient recomputation of EoS inference every time a new binary neutron star event is discovered or whenever new observations and theoretical discoveries change the prior on EoS parameters. We test our method on both real and simulated gravitational wave data to demonstrate its accuracy. We show that our computationally cheap method produces EoS constraints that are completely consistent with existing analysis for real data, the chosen fiducial EoS for simulated data. Armed with our fast analysis scheme, we also study the variability of EoS constraints with binary neutron star properties for sets of simulated events drawn in different signal-to-noise ratio and mass ranges.

Figure 1

KDE visualization for two simulated events, with APR4_EPP as the injected EoS, one at low (bottom) and the other at high (top) SNR. In both the figures, we have plotted $K_i(q,\tilde{\mathrm{\Lambda }})$ as 2D density plots. We have overplotted the PE samples of $q,\tilde{\mathrm{\Lambda }}$ as discrete gray points to demonstrate the KDE accurately approximates the posterior distribution from which the samples are drawn. The integral in Eq. (23) can be visualized as line integrals of these 2D density along the $\tilde{\mathrm{\Lambda }}(q,\mathcal{M},\overrightarrow{\gamma })$ line for some particular value of $\overrightarrow{\gamma }$. Two such lines, one corresponding to the injected EoS APR4_EPP and a different SLY are plotted as an example of the EoS lines along which the integrals are performed.

Figure 2

Visualization of the priors used on EoS parameters as bounds on the pressure density plane. We draw a large number of samples of the spectral parameters from the uniform distributions on the right-hand side of Eq. (26). From those samples, we only keep the ones that satisfy the physicality conditions imposed by the various component of the prior. We then calculate the pressure density curves for each of those samples and plot credible intervals on the pressure density plane that contain 90% of those curves.

Figure 3

Visualization of uninformative priors on tidal parameters used in single-event EoS agnostic PE runs, which can be one or the other depending on the waveform model being used. The uniform in the $\tilde{\mathrm{\Lambda }}$ prior, compatible with taylorf2 and hence the current version of gwxtreme , would be a horizontal line in this plot. The uniform in positive ${\mathrm{\Lambda }}_1$, ${\mathrm{\Lambda }}_2$, which is compatible with the phenomnrt waveform families, leads to a $\tilde{\mathrm{\Lambda }}$ prior distribution that is shown in this plot and evidently has vanishing support at $\tilde{\mathrm{\Lambda }}=0$. This can tend to blow up the integrand in Eq. (18) as its denominator $p_{\mathrm{PE}}(\mathcal{M},q,\tilde{\mathrm{\Lambda }},\delta \tilde{\mathrm{\Lambda }})$ will then be proportional to a quantity that vanishes a region inside the integration range.

Figure 4

Visualization of why high-mass events are less informative on the EoS. We plot 90% highest posterior density intervals of the $\tilde{\mathrm{\Lambda }}$ posteriors for a low-mass and a high-mass simulated event drawn from the narrow SNR bins, top: (23,25) and bottom: (33,35). We also display the shape of the posterior in the form of violin plots. The solid lines are the $\mathcal{M}-\tilde{\mathrm{\Lambda }}$ curves at $q=1$ for some of the known EoS models.

Figure 5

Comparison of gwxtreme with lalinference_nest for the real events GW170817 (left) and GW190425 (right). The posterior samples of the spectral parameters are used to compute pressure density curves, and the shaded regions in the plots are equal tail confidence intervals that contain 90% of those curves. We also display the 90% prior intervals in dashed lines and the prior extrema in dotted lines computed using 50,000 samples of $\overrightarrow{\gamma }$ drawn from the prior. It can be seen that gwxtreme results are consistent with the lalinference_nest results found by the LVK, despite being orders of magnitude faster. The slight preference of lalinference_nest toward softer EoS can be attributed to the difference in waveform models as well as the priors on the tidal parameters that are used by the two algorithms.

Figure 6

EoS constraints obtained by jointly analyzing GW170817 and GW190425 using gwxtreme. It can be seen that the joint EoS constraints are dominated by that of GW170817, which is to be expected due its larger SNR and smaller masses of component NSs than GW190425. However, the joint constraint can be seen to be very slightly narrower than that of GW170817 due to the contribution from GW190425.

Figure 7

EoS constraints obtained by jointly analyzing 16 events drawn from the Galactic BNS population with APR4_EPP as the true EoS and injected into O4 sensitivity. The posterior samples of the spectral parameters generated by gwxtreme can be used to compute $p(\rho ,\overrightarrow{\gamma })$ and $\mathrm{\Lambda }(m,\overrightarrow{\gamma })$ curves. In plot (a), the shaded region marks the equal-tail confidence interval that contains the 90% of the $p(\rho ,\overrightarrow{\gamma })$ curves. Plot (b) is obtained by histogramming $\mathrm{\Lambda }(1.4M_{\odot },{\overrightarrow{\gamma }}_i)$ corresponding to each posterior sample ${\overrightarrow{\gamma }}_i$. Both plots demonstrate that our computationally cheap and fast algorithm accurately measures the injected EoS.

Figure 8

EoS constraints obtained using gwxtreme for events drawn in narrow SNR ranges. For each chosen SNR range, we draw events uniformly in chirp masses and then draw the extrinsic parameters so that when injected into O4 sensitivity the simulated signals will have an SNR that belongs in the said range of SNRs. We then group the highest-mass and lowest-mass events in each SNR range and analyze them using our algorithm both separately and jointly. Then, using the posterior samples of the spectral parameters generated by our algorithm, we produce EoS constraints as both pressure-density credible intervals and quantile ranges of the tidal deformability at $1.4M_{\odot }$. It can be seen that EoS constraints are dominated by the lower-mass events in each SNR bin and that they become narrower with increasing SNR.

Figure 9

EoS constraints obtained using the piecewise polytrope for the 16 events drawn from the Galactic population distribution.

Figure 10

EoS constraints obtained using the piecewise polytrope for simulated events drawn in the 23 to 25 SNR bin.

Figure 11

EoS constraints obtained using the piecewise polytrope for events drawn in the 33 to 35 SNR bin.