Mapping the Universe Expansion: Enabling Percent-level Measurements of the Hubble Constant with a Single Binary Neutron-star Merger Detection

The "+" and "×" polarizations of a GW emitted by a compact binary merger located at a luminosity distance d_L can be expressed as a superposition of individual modes, h_ℓ,m , weighted by spin −2 spherical harmonics ${Y}_{{\ell },m}^{-2}$ as

$$\begin{eqnarray}&&{h}_{+}-{{ih}}_{\times }=\displaystyle \frac{1}{{d}_{L}}\sum _{{\ell }\geqslant 2}\sum _{m=-{\ell }}^{m={\ell }}{Y}_{{\ell },m}^{-2}(\iota ,\varphi ){h}_{{\ell },m}({\rm{\Xi }}).\end{eqnarray} \tag{ 1 }$$

Here, Ξ denotes the masses m_i and dimensionless spins χ _i of the individual objects and, for the case of BNSs, the individual tidal deformabilities Λ_i characterizing the deformation of each star in the external gravitational field of the companion. The parameters (ι, φ) represent the polar and azimuthal angles of a spherical coordinate frame describing the location of the observer around the binary (or conversely, the orientation of the binary with respect to the observer), with ι = 0 denoting the direction of the orbital angular momentum of the binary and ι = π/2 denoting the orbital plane. These values, respectively, refer to face-on and edge-on oriented binaries. For nonprecessing binaries, the above sum is dominated by the quadrupole (ℓ, m) = (2, ± 2) modes while HMs become loud only during the final inspiral and merger phase, with increasing relative amplitude as the mass ratio q = m₁/m₂ ≥ 1 increases (Pekowsky et al. 2013; Varma et al. 2014; Bustillo et al. 2015; Calderón Bustillo et al. 2017).

Current parameter estimation of BNS signals makes use of waveform templates including only the quadrupole mode. This causes a degeneracy between the inclination and distance parameters that fundamentally limits our ability to measure each. Several works have shown that inclusion of HMs in templates can break this degeneracy for sources with sufficiently loud HMs in the detector sensitive band, as the observed combination of modes will depend on the orientation of the binary via the Y_ℓ,m factors (Graff et al. 2015; Calderón Bustillo et al. 2018; London et al. 2018; Pang et al. 2018; Calderón Bustillo et al. 2019). For Advanced LIGO and Virgo observations, unfortunately, this is only possible for large mass and asymmetric BBHs (Graff et al. 2015; London et al. 2018; Chatziioannou et al. 2019; LIGO Scientific Collaboration 2020) for which the merger and ringdown emission, rich in HMs, is strong in the detector sensitive band. In contrast, the merger and post-merger of BNSs is unobservable due to its large frequency. This emission will, however, be observable with future high-frequency (Martynov et al. 2019; Ackley et al. 2020) and third-generation (Punturo et al. 2010a; Reitze et al. 2019) detectors and, in addition, the pre-merger GW emission will last for several minutes in the sensitive detector band, allowing us to accumulate the effect of weak HMs.

We test our ability to measure the source distance and H₀ using the inspiral and post-merger emission of BNSs. To this, we perform parameter inference on two kinds of simulated GW signals. First, we use 80 ms long numerical-relativity simulations for the post-merger emission of BNS (Dietrich & Hinderer 2017). These have mass ratios q = 1 and q = 1.5 and implement two different equations of state (EOSs): a soft one (SLy) and a stiff one (MS1b ⁷ ). The left panel of Figure 1 shows the time domain modes for the equal-mass SLy case located at a distance of 40 Mpc. The right panel shows in blue and green the spectra of full waveforms observed face-on (blue) and edge-on (green), together with the contribution of the HMs and the (ℓ, m) = (2, 1) mode alone for the edge-on case. It can be noted not only how the face-on signal is stronger, but how the presence of HMs in the edge-on signal leads to noticeable morphological differences. Parameter inference on these short waveforms provides an idea of how the post-merger emission breaks the distance-inclination degeneracy. Restricting to this and ignoring the long inspiral signal, however, would greatly underestimate the S/N accumulated throughout the full minutes-long signal observable by the detector, reducing the accuracy of our measurements.

Zoom In Zoom Out Reset image size

To obtain signals that can cover the full inspiral and post-merger emission from a BNS, we combine (Bustillo et al. 2015) our short numerical-relativity waveforms (Dietrich & Hinderer 2017) with long analytical waveforms covering the early inspiral minutes before the merger (e.g., computed using the tidal effective-one-body model of Nagar et al. 2018). Unfortunately, parameter inference using these waveforms is computationally prohibitive. As a solution, we implement a two-step approach. First, we consider 128 s long phenomenological (phenom) waveforms (Khan et al. 2016; Santamaria et al. 2010), constructed with the IMRPhenomHM model (London et al. 2018), covering the inspiral-merger and ringdown stages of nonprecessing binary black hole (BBH) mergers. While computationally inexpensive, these waveforms omit the two main characteristic aspects of BNSs: tidal-deformability effects and the post-merger emission shown in Figure 1. Since, as argued, the latter can improve our distance measurements, results obtained using this BBH signal model are rather conservative. Finally, in order to obtain improved and more realistic results for BNSs, we combine these BBH parameter estimates with those obtained by solely analyzing the short numerical-relativity waveforms covering the post-merger stage of BNSs, following the procedure described in Zimmerman et al. (2019).

We inject signals h(θ_true) with source parameters θ_true, that include HMs, in zero noise and estimate the source parameters using waveform templates h(θ) that omit and include HMs. We compute the posterior Bayesian probability of the parameters θ as

$$\begin{eqnarray}&&p(\theta | {\theta }_{\mathrm{true}})=\displaystyle \frac{\pi (\theta ){ \mathcal L }(\theta | h({\theta }_{\mathrm{true}}))}{\int \pi (\theta ){ \mathcal L }(\theta | h({\theta }_{\mathrm{true}}))d\theta },\end{eqnarray} \tag{ 2 }$$

with π(θ) denoting the prior probability of the parameters θ and ${ \mathcal L }(\theta | h({\theta }_{\mathrm{true}}))$ denoting their likelihood. As usual, the latter is defined as the standard frequency-domain likelihood for GW transients (Finn 1992; Romano & Cornish 2017)

$$\begin{eqnarray}&&\mathrm{log}{ \mathcal L }(\theta | h({\theta }_{\mathrm{true}}))\propto -\sum _{N}\displaystyle \frac{(h({\theta }_{\mathrm{true}})-h(\theta )| h({\theta }_{\mathrm{true}})-h(\theta ))}{2},\end{eqnarray} \tag{ 3 }$$

where N runs over the different detectors of our network. As we discuss later, we work with two- and three-detector networks. As usual, (a∣b) represents the inner product (Cutler & Flanagan 1994)

$$\begin{eqnarray}&&(a| b)=4{\mathfrak{R}}{\int }_{{f}_{\min }}^{{f}_{\max }}\displaystyle \frac{\tilde{a}(f){\tilde{b}}^{* }(f)}{{S}_{n}(f)}{df},\end{eqnarray} \tag{ 4 }$$

where $\tilde{a}(f)$ denotes the Fourier transform of a(t) and * denotes the complex conjugate. The factor S_n (f) denotes the one-sided power spectral density of the detector. In this work, we consider a network of detectors, all with noise sensitivity equivalent to that of the proposed 2.5-generation Neutron star Extreme Matter Observatory (NEMO; Ackley et al. 2020). We choose a lower frequency cutoff of ${f}_{\min }=20$ Hz and a sampling frequency of 16 kHz so that ${f}_{\max }=8$ kHz. The NEMO detector has a proposed sensitivity similar to the Cosmic Explorer and Einstein Telescope in the kilohertz regime; we could equally use those third-generation detectors and achieve similar results for the late inspiral and post-merger, albeit with larger signal-to-noise ratios for the full signal given the better low-frequency (≲500 Hz) sensitivity.

In all cases, we assume standard prior probabilities for the sky-location, source orientation, and polarization angle, together with a prior uniform in comoving volume and a uniform prior on the time-of-arrival, with a width of 0.2 s, centered on the true value. For our analyses making use of 128 s long phenom waveforms, we impose uniform priors on the individual masses m_1,2 ∈ [1, 2] M_⊙ and on the components of the individual spins along the orbital angular momentum ${\chi }_{1,2}^{z}\in [-0.15,0.15]$. Since numerical-relativity waveforms, however, are only produced for a discrete set of intrinsic parameters, we assume the masses and spins to be known in this case. We find this is a reasonable assumption as the individual masses and the effective spin parameter (Santamaria et al. 2010) ${\chi }_{\mathrm{eff}}=({\chi }_{1}^{z}{m}_{1}+{\chi }_{2}^{z}{m}_{2})/({m}_{1}+{m}_{2})$ are very well measured from the long inspiral. As an example, with a triple-detector network and using HMs, we determine the total mass, chirp mass, and effective-spin parameters of our face-on unequal-mass source with respective uncertainties of <1%, <0.01%, and <0.015 at the 68% level. We perform our parameter inference runs with the software Bilby (Ashton et al. 2019; Smith et al. 2020), sampling the parameter with the algorithm CPNest (Veitch et al. 2020).

We consider three network configurations. The first (denoted HV) assumes two detectors with the location and orientation of Advanced LIGO Hanford and Virgo. Such a network has each detector sensitive to one of the two independent GW polarizations. The second network (HL) assumes LIGO Hanford and Livingston location and orientations, almost anti-aligned, so that both detectors are sensitive to roughly the same GW polarization, missing the other one. Finally, we consider an HLV network that is sensitive to both GW polarizations and can pinpoint the sky-location of the source.

The accuracy of our distance measurement is of course limited by the loudness of the injected signals in our detector network, quantified by the optimal network S/N, given by

$$\begin{eqnarray}&&{\rho }_{\mathrm{opt}}(h({\theta }_{\mathrm{true}}))=\sqrt{\sum _{N}(h({\theta }_{\mathrm{true}})| {\left(h({\theta }_{\mathrm{true}})\right)}_{N}},\end{eqnarray} \tag{ 5 }$$

which is inversely proportional to the source luminosity distance d_L . As we show in Appendix C, for face-on cases we obtain optimal S/Ns of ∼190 for HL and HLV networks and ∼140 for HV. In contrast, for the weaker edge-on cases we obtain respective values ρ_opt ∈ [30, 65] depending on the mass ratio and network considered. For comparison, the S/N of GW170817, whose source was rather face-on, was only ≃ 32. We note that while in this study we restrict to sources placed at distances d_L,true = 40 Mpc, our results can be extended to different reference distances. In particular, in our S/N-regime, given fixed source and detector network, the uncertainty of our distance estimates roughly depends on the optimal S/N (and therefore on the reference distance) as ${\rm{\Delta }}({d}_{L})\propto 1/{\rho }_{\mathrm{opt}}^{2}\propto {d}_{L,\mathrm{true}}^{2}$ (Rao 1992; Cramer 1999; Li 2013). In addition, uncertainties in the estimated value of the Hubble constant H₀ = cz/d_L,true grow linearly with distance.

We choose four target sources with total mass M = 2.75 M_⊙, mass ratios q = 1 and q = 1.5, and oriented both face-on (ι = 0) and edge-on (ι = π/2). Signals emitted in these two angles differ in three aspects: morphology, polarization, and loudness. Face-on signals contain solely the (ℓ, m) = (2, 2) mode, are circularly polarized (i.e., both h_+,× have the same amplitude), and are louder than edge-on signals. On the contrary, edge-on signals have contributions from higher-order emission modes that confer a richer structure, but are weaker in amplitude than face-on ones (see Figure 1). In addition, edge-on signals contain only one of the two polarizations. Consequently, it has been shown that measuring the ratio of the two polarizations is key to correctly inferring the inclination of the source, provided that the detector network can observe both polarizations (Usman et al. 2019). For each of these sources, we consider two EOSs, namely SLy and MS1b, which we assume to be known. Different EOSs trigger post-merger HMs in different ways, varying the accuracy of the distance estimate. Finally, we would ideally consider a wide range of distances and sky locations for all of our target sources. However, given the extreme computational cost of our parameter inference runs, and to allow for a direct comparison of our results, we place all of our sources at the same distance and sky location. For the former, we considered that the most reasonable choice was a value of d_L = 40 Mpc, consistent with that of GW170817, the only conclusive BNS observed to date through GWs. Finally, we placed all of our target sources at the same, arbitrary sky location.

Figure 2 shows the two-dimensional posterior distributions for the luminosity distance and inclination of face-on oriented binaries with mass ratios of q = 1 (left) and q = 1.5 (right), using an HLV detector configuration. The contours denote the 90% credible regions and the legend provides median estimates with symmetric 68% credible intervals. For unequal masses, and using 128 s long phenom waveforms, the omission of HMs in the templates h(θ) leads to a biased estimate of ${d}_{L}^{\mathrm{Phenom},\ 22}={36.2}_{-2.5}^{+2.2}$ Mpc. Their inclusion corrects this bias and reduces the uncertainty by ∼70%, yielding a distance measurement of ${d}_{L}^{\mathrm{Phenom},\ \mathrm{HM}}={39.4}_{-0.9}^{+0.4}$. For equal masses, the impact of HMs is significantly milder, as the (usually strongest) odd-m emission modes are suppressed (Pan et al. 2011, 2014; Blanchet 2014). This leads to a biased estimate of ${d}_{L}^{\mathrm{Phenom},\ \mathrm{HM}}={37.9}_{-2.3}^{+1.4}$ even if HMs are included. We also note that this remains true even if one would assume the intrinsic source parameters (masses and spins) to be known. This situation changes, however, when using 80 ms long numerical-relativity waveforms that can account for the rich post-merger signal morphology. For both mass ratios, and considering a SLy equation of state, results are better than those making use of 128 s long phenom waveforms omitting HMs. Moreover, for equal-mass, the estimate even improves on that including HMs, yielding a nonbiased estimate ${d}_{L}^{{\mathtt{SLy}}}={38.6}_{-2.7}^{+1.8}$ with smaller uncertainties.

Zoom In Zoom Out Reset image size

We combine the distance estimates obtained using numerical-relativity waveforms with those obtained analyzing phenom waveforms restricted to frequencies not covered by the former. To this, we multiply the respective posterior distributions for the distance, dividing by one instance of the prior (Zimmerman et al. 2019). ⁸ We obtain joint estimates of ${d}_{L}^{\mathrm{joint}}={38.6}_{-1.3}^{+0.9}$ Mpc for the q = 1 case and ${d}_{L}^{\mathrm{joint}}={39.4}_{-0.7}^{+0.4}$ for the q = 1.5 case. The reason behind this improvement is that matter effects arising during the post-merger of BNSs trigger HMs, helping to break the degeneracy between distance and inclination and even allowing us to measure the azimuthal angle (see Appendix B). For the unequal-mass case, these only add a small contribution with respect to the integrated effect of the HMs during the 128 s of signal. For q = 1, however, odd-m modes are suppressed during the inspiral and are only triggered during the post-merger of our numerical-relativity simulations due to an effect known as one-armed spiral instability or 21-mode instability (East et al. 2016; Lehner et al. 2016; Radice et al. 2016; see Appendix B). As a consequence, the inclusion of the post-merger emission can have a large impact. In our study, the distance estimates are much better for SLy compared to MS1b, due to the stronger HM emission.

Finally, while we have discussed results assuming an HLV detector network, in the Suppl. Material we show results for all studied network configurations. In all cases we obtain similar results, both qualitatively and quantitatively.

Previous work has shown that the degeneracy between distance and inclination can be broken by computing the ratio h_×/h₊, as this evolves from 1 to 0 as the inclination varies from face-on to edge-on (Usman et al. 2019). We find that for face-on binaries, all HL, HV, and HLV configurations yield almost equivalent distance measurements despite the differences in signal loudness across the network, so that HMs have a much larger impact on the measurement than the polarization ratio. In contrast, in Figure 3 we show that for edge-on cases it is key not only to access both signal polarisations but also to have a third detector that can pin-point the sky-location of the source. For an HLV configuration, unbiased estimates with uncertainties lower than 4% are obtained regardless of the usage of HMs, while biased estimates are obtained using both two-detector configurations. The reason is that such configurations cannot pin-point the sky location of the source using timing information. This way, Bayesian inference places the source at those patches of the sky, consistent with the two-detector timing, where the detector network is most sensitive, biasing the distance toward large values. We obtain identical qualitative results for the q = 1 case, as well as when analyzing our 80 ms long numerical-relativity simulations.

Zoom In Zoom Out Reset image size

Combining our distance estimates with simulated redshift estimates, we can infer H₀ via H₀ = cz/d_L , with c being the speed of light. We assume redshift estimates consisting of Gaussian posterior distributions centered at a value of z₀ = 0.00897, corresponding to d_L = 40 Mpc in a ΛCDM cosmology with Hubble parameter H₀ = 67.9 km s⁻¹ Mpc⁻¹. We assume that the redshift has an uncertainty with standard deviation Δz ranging from a realistic value of 10⁻³ consistent with that for GW170817 (Abbott et al. 2017c) to an improved value of 10⁻⁴ which would require higher resolution spectrograms and also the possibility to measure the internal motion of sources within individual galaxies, e.g., (Davis et al. 2019).

Figure 4 shows the H₀ estimates derived from the distance measurements of our face-on BNSs with mass ratios q = 1 and q = 1.5 as a function of Δz using an HLV network. Once again, we quote results in terms of symmetric 68% credible intervals centered at the median value. For unequal mass, and for Δz = 10⁻³, the omission of HMs never leads to biased estimates. Their inclusion, however, leads to an important improvement from ${H}_{0}={73.8}_{-9.8}^{+9.0}$ to ${H}_{0}={68.3}_{-7.5}^{+7.5}$. More spectacularly, for the most optimistic Δz = 10⁻⁴ HMs improve the measurement from ${H}_{0}={73.8}_{-4.1}^{+5.3}$ to ${H}_{0}={68.0}_{-1.0}^{+1.6}$, enabling a 2%-level measurement. This shows that with 2.5 G detectors, together with improved detector calibration and waveform models, H₀ estimates will be limited by EM redshift measurements and not by GW distance ones. Conversely, if HMs are omitted, a significant reduction of Δz will not translate into an improved H₀ estimate. Moreover, we find that H₀ estimates would be biased when Δz ≲ 0.5. Consistently with the previous section, the inclusion of the post-merger emission in the analysis does not lead to any relevant improvement.

Zoom In Zoom Out Reset image size

As expected, the situation is different for equal-mass cases. For these, the inclusion of post-merger effects is crucial to obtain visible improvements in the H₀ measurement. For Δz = 10⁻³ we obtain a mild improvement from ${H}_{0}={71.2}_{-8.3}^{+8.6}$ to ${H}_{0}={69.6}_{-7.7}^{+7.8}$. When Δz is reduced to Δz = 10⁻⁴, the HMs present in the post-merger emission allow for an estimate ${H}_{0}={69.3}_{-1.5}^{+2.5}$, with uncertainties at the ≃4% percent level, while their omission doubles the uncertainty and biases the measurement.

Almost identical results hold for the HL and HV networks, as shown in the Suppl. Material. For the (weaker) edge-on binaries, we find percent-level measurements are possible using the HLV network regardless of the usage of HMs.