Hierarchical Bayesian method for constraining the neutron star equation of state with an ensemble of binary neutron star postmerger remnants

Binary neutron star (BNS) postmerger gravitational-wave emission can occur in the aftermath of a BNS merger—provided the system avoids prompt collapse to a black hole—as a quasistable hypermassive remnant experiences quadrupolar oscillations and nonaxisymmetric deformations. The postmerger gravitational-wave spectrum possesses a characteristic peak frequency that has been shown to be dependent on the binary chirp mass and the neutron star equation of state (EOS), rendering postmerger gravitational waves a powerful tool for constraining neutron star composition. Unfortunately, the BNS postmerger signal is emitted at high ($\gtrsim 1.5\mathrm{kHz}$) frequencies, where ground-based gravitational-wave detectors suffer from reduced sensitivity. It is therefore unlikely that postmerger signals will be detected with sufficient signal-to-noise ratio (SNR) until the advent of next-generation detectors. However, by employing empirical relations derived from numerical relativity simulations, we can combine information across an ensemble of BNS mergers, allowing us to obtain EOS constraints with many low-SNR signals. We present a hierarchical Bayesian method for deriving constraints on $R_{1.6}$, the radius of a $1.6M_{\odot }$ neutron star, through an ensemble analysis of binary neutron star mergers. We apply this method to simulations of the next two LIGO-Virgo-KAGRA observing runs, O4 and O5, as well as an extended four-year run at $\mathrm{A}+$ sensitivity, demonstrating the potential of our approach to yield EOS information from the postmerger signal with current-generation detectors. The $\mathrm{A}+$ four-year scenario is predicted to improve the constraint on $R_{1.6}$ from the currently available multimessenger-based 95% credible interval (C.I.) uncertainty of $R_{1.6}=12.07_{-0.77}^{+0.98}$ to $R_{1.6}=11.91_{-0.56}^{+0.80}\mathrm{km}$, a 22% reduction of the 95% C.I. width.

Figure 1

GW power spectral density (PSD) of a simulated $\mathrm{merger}+\text{postmerger}$ waveform for the SFHX EOS (see Sec. 2f) with $m_1=1.03M_{\odot }$ and $m_2=1.09M_{\odot }$, injected at a distance of 50 Mpc. The current, projected, and proposed detector noise PSDs discussed in this work are provided for reference: that of the LIGO Livingston detector during the third LVK observing run (O3 L1), the projected LIGO sensitivities for the fourth and fifth LVK observing runs (O4 and O5), assuming LIGO $\mathrm{A}+$ design sensitivity for O5 [61], and the proposed NEMO and Einstein Telescope detectors (see Sec. 4a) [62, 63]. The postmerger peak frequency of 2.67 kHz is marked.

Figure 2

Distribution of $N_{\mathrm{R}}=1000$ sampled relations evaluated at a range of values of $R_{1.6}$, alongside our fitted empirical relations, likewise evaluated over the same range. The sampled distribution is consistent with the fit relation, and characterizes the relation’s intrinsic uncertainty.

Figure 3

bayeswave signal reconstruction and corresponding peak frequency posterior for a strong BNS postmerger signal. This bayeswave run was performed on a SFHX postmerger waveform with $m_1=1.03M_{\odot }$ and $m_2=1.09M_{\odot }$ injected at 30 Mpc at NEMO sensitivity. In the top panel, the injected waveform is in teal and the median bayeswave reconstruction with 50% (dark shading) and 90% (light shading) C.I. is shown in purple. The postmerger signal is reconstructed nicely, and the injected waveform’s peak frequency $f_{\mathrm{peak}}=2.67\mathrm{kHz}$ is recovered to within $\sim 0.01\mathrm{kHz}$.

Figure 4

bayeswave signal reconstruction and corresponding peak frequency posterior for a marginal BNS postmerger signal. This bayeswave run was performed on a postmerger waveform with $m_1=1.03M_{\odot }$ and $m_2=1.09M_{\odot }$ injected at 70 Mpc at NEMO sensitivity. In the top panel, the injected waveform is in teal and the 50% (dark shading) and 90% (light shading) C.I. of the bayeswave reconstruction are shown in purple. No median reconstruction is shown because the reconstructed signal is consistent with zero. Despite this, the spectrum and corresponding $f_{\mathrm{peak}}$ posterior show a slight preference for the injected waveform’s peak frequency of 2.67 kHz. Taken by itself, this recovery of the peak frequency is unconvincing—and would be considered a nondetection—but nonetheless contains information about the postmerger signal.

Figure 5

$R_{1.6}$ priors used in this work. One is a uniform prior on [9, 15] km, the other is an astrophysically motivated $R_{1.6}$ prior derived from the multimessenger analysis of Huth et al. [35].

Figure 6

Peak frequency prior samples, with a reflection-bounded KDE approximation of its probability distribution. These samples were computed from an ensemble of $N=50$ bayeswave analyses over simulated instrumental noise.

Figure 7

Neutron star mass-radius relationship for each of the equations of state used in this work, along with each relationship’s intersection with $\mathrm{M}=1.6M_{\odot }$. The corresponding values of $R_{1.6}$ are 11.54 km (SLy4), 11.98 km (SFHX), and 13.26 km (DD2).

Figure 8

Simulation parameters (postmerger peak frequency $f_{\mathrm{peak}}$ and chirp mass $\mathcal{M}$) for each EOS, with the corresponding Vretinaris et al. [93] empirical relation with $3\sigma$ bounds for that EOS’s value of $R_{1.6}$.

Figure 9

100 realizations of $R_{1.6}$ posteriors arising from 100 events consisting solely of noise draws. These distributions are representative of what can be produced from random fluctuations in the $f_{\mathrm{peak}}$ posteriors, absent any impact from an embedded signal. As expected, while we do observe some level of fluctuation in these $R_{1.6}$ posteriors, the set as a whole is consistent with the prior. The flattening at low $R_{1.6}$ arises from our band-limited search; as discussed in Sec. 2b, most predicted values of ${\widehat{f}}_{\mathrm{peak}}$ for $R_{1.6}\lesssim 10\mathrm{km}$ fall out of band and do not contribute information to our final posterior. As a result, most analyses without a significant signal present in the data will tend to cleave to the prior for these values of $R_{1.6}$.

Figure 10

Sample synthetic $f_{\mathrm{peak}}$ posteriors and corresponding 100-event ensemble $R_{1.6}$ posteriors for increasing ratio of signal samples to noise samples (S:N). Note that in all cases, $N_s$ is sufficiently small that the corresponding sample $f_{\mathrm{peak}}$ posterior is not strongly peaked (and in fact by eye appears to be entirely uninformative). For reference, we have also plotted the set of noise-only $R_{1.6}$ posteriors described in Sec. 3a. From top to bottom, these posteriors correspond to $R_{1.6}$ means and 95% C.I. of ${12.3}_{-3.2}^{+2.5}$, ${12.1}_{-2.9}^{+2.1}$, ${12.5}_{-2.6}^{+0.9}$, and ${12.59}_{-0.31}^{+0.33}\mathrm{km}$, respectively.

Figure 11

100-event ensemble $R_{1.6}$ posteriors with $N_s=50$, $N_n=5000$ for different injected values of $R_{1.6}$. For reference, we have also plotted the set of noise-only $R_{1.6}$ posteriors shown in Fig. 9. We accrue posterior density around the injected value in each case, although the recovery strength is reduced by mitigating factors at low and high $R_{1.6}$.

Figure 12

Scaling with ${\mathrm{N}}_{\mathrm{events}}$ of the ensemble $R_{1.6}$ posterior 68% C.I. width as a fraction of the prior C.I. width for a 2000-event ensemble with $R_{1.6}=12.5\mathrm{km}$ and four different injected signal strengths. Note that the same behavior occurs for other C.I. percentages, but viewing the entire trend at low signal strengths e.g., for 95% C.I. requires many more events.

Figure 13

$R_{1.6}$ posteriors for simulations of LVK’s upcoming O4 and O5 observing runs for the SLy4, SFHX, and DD2 equations of state. These posteriors correspond to means and 95% C.I. of $R_{1.6}={11.8}_{-2.6}^{+3.0}$, $R_{1.6}={12.0}_{-2.8}^{+2.8}$, and $R_{1.6}={11.8}_{-2.7}^{+3.0}\mathrm{km}$, respectively, as compared to the original uniform prior mean and 95% C.I. of $R_{1.6}={12.0}_{-2.9}^{+2.9}\mathrm{km}$. While there appear to be hints of a signal in the SLy4 and SFHX cases, the final posterior is not conclusive and could arise from noise fluctuations, as seen in the case of DD2.

Figure 14

$R_{1.6}$ posteriors for simulations of LVK’s upcoming O4 and O5 observing runs with the SFHX EOS and the multimessenger prior of Huth et al. [35]. The individual O4 and O5 posteriors are presented along with the combined $\mathrm{O}4+\mathrm{O}5$ posterior. For O4, O5, and $\mathrm{O}4+\mathrm{O}5$, these posteriors correspond to means and 95% C.I. of $R_{1.6}=12.07_{-0.75}^{+0.95}$, $R_{1.6}=12.02_{-0.70}^{+0.95}$, and $R_{1.6}=12.02_{-0.68}^{+0.93}\mathrm{km}$, respectively. For comparison, the prior mean and 95% C.I. are $R_{1.6}=12.07_{-0.77}^{+0.98}\mathrm{km}$. The combined $\mathrm{O}4+\mathrm{O}5$ posterior provides a reduction of the 95% C.I. width by 0.13 km (7.5%) over the prior.

Figure 15

$R_{1.6}$ posteriors arising from a simulation of four years of operation at $\mathrm{A}+$ sensitivity for the SFHX EOS, presented with both a uniform prior and the multimessenger prior of Huth et al. [35]. For the multimessenger (uniform) prior, we attain posterior means and 95% C.I. of $R_{1.6}=11.91_{-0.56}^{+0.80}\mathrm{km}$ ($R_{1.6}={11.9}_{-2.6}^{+2.7}\mathrm{km}$). For reference, the multimessenger (uniform) prior means and 95% C.I. are $R_{1.6}=12.07_{-0.77}^{+0.98}\mathrm{km}$ ($R_{1.6}={12.0}_{-2.9}^{+2.9}\mathrm{km}$). In the case of the multimessenger prior, we achieve a reduction in the 95% C.I. width of 0.39 km (22%).