Cosmology with Love: Measuring the Hubble constant using neutron star universal relations

Gravitational-wave cosmology began in 2017 with the observation of the gravitational waves emitted in the merger of two neutron stars, and the coincident observation of the electromagnetic emission that followed. Although only a 30% measurement of the Hubble constant was achieved, future observations may yield more precise measurements either through other coincident events or through cross correlation of gravitational-wave events with galaxy catalogs. Here, we implement a new way to measure the Hubble constant without an electromagnetic counterpart and through the use of the binary-Love relations. These relations govern the tidal deformabilities of neutron stars in an equation of state insensitive way. Importantly, the Love relations depend on the component masses of the binary in the source frame. Since the gravitational-wave phase and amplitude depend on the chirp mass in the observer (and hence redshifted) frame, one can in principle combine the binary-Love relations with the gravitational-wave data to directly measure the redshift, and thereby infer the value of the Hubble constant. We implement this approach in both real and synthetic data through a Bayesian parameter estimation study in a range of observing scenarios. We find that for the LIGO/Virgo/KAGRA design sensitivity era, this method results in a similar measurement accuracy of the Hubble constant to those of current-day, dark-siren measurements. For third-generation detectors, this accuracy improves to $\lesssim 10\%$ when combining measurements from binary neutron star events in the LIGO Voyager era, and to $\lesssim 2\%$ in the Cosmic Explorer era.

Figure 1

Upper panel: Mass-radius relations for 29 piecewise polytropic EoSs that are consistent with the LIGO/Virgo and NICER observations to 68% confidence. Middle panel: ${\overline{\lambda }}_0^{(k)}={\overline{\lambda }}_0^{(k)}({\overline{\lambda }}_0^{(0)})$ relations for the same set of EoSs, together with the fit in thick lines. Observe that the ${\overline{\lambda }}_0^{(0)}-{\overline{\lambda }}_0^{(k)}$ relations are EoS insensitive and the fit recovers the average. Bottom panel: Distribution of the fractional-residual errors from the fit for each coefficient.

Figure 2

Posterior probability density for the ${\overline{\lambda }}_0^{(0)}$ parameter measured from 128 s of strain data (4 kHz) for GW170817. The ${\overline{\lambda }}_0^{(0)}-{\overline{\lambda }}_0^{(k)}$ relations are used in the parametrization of the tidal deformability in Eq. (1) with (unfilled) and without (filled) the residual errors in the relations. The redshift is fixed at the host galaxy $z_{\mathrm{NGC}4993}=0.0099$, and therefore, the tidal deformability is parametrized only by ${\overline{\lambda }}_0^{(0)}$. The vertical dash-dot lines represent the median, the vertical dashed lines represent the 90% confidence interval, and the horizontal line is the prior density.

Figure 3

Left: Individual (dotted) and stacked (solid) $H_0$ likelihoods obtained for the representative runs in the O5 (top), Voyager (middle) and CE (bottom) observing eras. Due to distance-inclination degeneracy, the individual posteriors are sometimes peaked away from the true value. However, after stacking the posteriors based on a sample of $N_{\det }$ events from the detected population, the distribution converges to the true simulated value, shown by the vertical line, $H_0=70\mathrm{km}{\mathrm{s}}^{-1}/\mathrm{Mpc}$. For illustration, we choose $N_{\det }=10$, 50, 70 for the top, middle, and bottom panels, corresponding to the O5, Voyager, and CE era, respectively. Right: The normalized number count of of distances and inclinations of the systems that would be detectable in a specific observing era is shown by the contour map. This was obtained by simulating an ensemble (several thousands) of binaries with $(m_1,m_2)\sim (1.42,1.27)M_{\odot }$ distribution uniform in sky, volume, and inclination, and requiring that each produce network $\mathrm{SNR}\ge 8$. The grid of overlayed points represents the true distance/inclination of synthetic systems chosen for performing Bayesian parameter estimation. We observe that the $H_0$ measurement converges irrespective of the shape of grid as long as it covers the distances and inclinations that would be detected.

Figure 4

Top Panel: The fractional error in $z$ as a function of the width in ${\overline{\lambda }}_0^{(0)}$ prior obtained using Fisher analysis on the model in Eq. (13). We see that the fractional uncertainty in $z$ follows the trends found in Ref. [10] (analogous to the $\delta {\overline{\lambda }}_0^{(0)}=0$ case). We see that, the fractional error in redshift does not change significantly with change in the width of ${\overline{\lambda }}_0^{(0)}$ prior at large distances from the source. The Fisher approximation is however valid in the high SNR regime, which we assume to be the region of $\mathrm{SNR}\gtrsim 30$ shown by the shading. We find similar trends from Bayesian parameter estimation results shown in the bottom panel. Bottom Panel: Fractional statistical uncertainty in the measurement of the redshift, $z$, the luminosity distance, $D_L$, and $H_0$ when using a Gaussian prior on ${\overline{\lambda }}_0^{(0)}$ versus when using a $\delta$-function prior (fixed) are shown for a few representative parameter estimation runs in the CE era.

Figure 5

Top panel: Comparison of distance, redshift, and $H_0$ measurements for an example parameter estimation run injected at a $D_L=2\mathrm{Gpc}$ in the CE era. Two cases of recovery are shown—one using a prior at the true value $\delta ({\overline{\lambda }}_0^{(0)}-200)$ and another shifted $\delta ({\overline{\lambda }}_0^{(0)}-230)$ prior. We observe no significant change in the distance recovery, but a systematic shift in the recovery of the redshift; the latter impacts the $H_0$ measurement. Bottom panel: The difference in means of the $H_0$ measurements between the true and shifted cases by doing a few more parameter estimation runs as above. This is shown by the dotted line with filled circles. We repeat the runs, putting a Gaussian prior with standard deviation of 30 for ${\overline{\lambda }}_0^{(0)}$, centered at the true value of 200, and plot half of the 90% credible interval for $H_0$ in dotted line with a plus marker. We repeat one more time with the same width but now centered at 230, and plot half credible region in using dotted lines with inverted triangles. We find that the credible region is not affected by a systematic shift in ${\overline{\lambda }}_0^{(0)}$. Therefore, we consider the former, and show improvements ($\propto 1/\sqrt{N_{\det }}$) as we increase the number of detections, using solid lines. Observe that for $N_{\det }\sim 30$, the statistical improvement hits the boundary of systematic error, shown in blue. We also show the analytical estimate for the systematic error from Eq. (14) using the shading, which roughly agrees with the parameter estimation runs.

Figure 6

An analytic estimate of the bias in $H_0$ caused by assuming the nonuniversality of the binary-Love relations (see Sec. 5c). As an example, here we consider the true EoS to be MPA1 which has the largest residual among the 29 EoSs from the fit in Eq. (3). We see that the error due to loss of universality reduces for ${\overline{\lambda }}_0^{(k)}$ as we increase $k$.

Figure 7

The implied (marginal) priors on the individual tidal deformability parameters due to the uniform prior in ${\overline{\lambda }}_0^{(0)}$ in Sec. 3. The prior on ${\overline{\lambda }}_1$ is shown in the hatched histogram, and ${\overline{\lambda }}_2$ is shown in the filled histogram. The ${\overline{\lambda }}_0^{(0)}$ prior is shown in the unfilled histogram.

Figure 8

Left panel: Nonmarginalized version of Fig. 7. Here, we have used the parameterization using ${\overline{\lambda }}_0^{(0)}-{\overline{\lambda }}_0^{(k)}$ relations from Sec. 2. Observe the correlation between ${\overline{\lambda }}_{1,2}$, which is expected due to the prior information about the masses when imposing a common EoS. Right panel: The prior on implied ${\overline{\lambda }}_{1,2}$ when using the universal relation between the symmetric and anti-symmetric tidal deformability, ${\overline{\lambda }}_s-{\overline{\lambda }}_a$, relations from YY17. This was used in the EoS-insensitive results reported previously in Ref. [16].

Figure 9

Extended (corner) plot showing pair and marginal parameters from performing Bayesian parameter estimation on GW170817 data discussed in Sec. 3. The parametrization of the tidal deformabilities follows Eq. (1). The redshift of the source is fixed to 0.0099. The bottom right posterior for ${\overline{\lambda }}_0^{(0)}$ is the same as Fig. 2 in the main text. We show the distribution of parameters in both cases—when using the fit in Eq. (1) directly and when including the residual errors from the fit mentioned in Sec. 2.

Figure 10

This is an example of $D_L$, $z$, and $H_0$ measurements for different combinations of the true inclination angle, ${\theta }_{JN}$, at a fixed distance of 3 Gpc in the Cosmic Explorer era. Note that with increasing true inclination angle there is increasing support for larger distance values. This is expected since a larger inclination implies a lower injected amplitude, which is recovered as a larger distance. The redshift measurement is not affected since it comes from the phase of the waveform.

Figure 11

Upper panel: The combined measurement of the $H_0$, similar to Fig. 3 for the Cosmic Explorer example. Here, we have reweighted the normalized number count of the recovered events using a redshift prior following uniform in comoving volume (dashed), and that following star formation history (solid). Bottom panel: The contour plot showing the normalized number counts similar to Fig. 3 for the two alternative priors. We observe that the combined measurement of $H_0$ is not significantly affected due to the prior on redshift.

Figure 12

Illustration of obtaining the likelihood from the $H_0$ posterior. In this example, we injected a low SNR ($\sim 1$) injection, and perform parameter estimation mentioned in Appendix pp6. The original prior on $H_0$ is shown by solid unfilled histogram, the posterior is shown in filled histogram. The reweighted prior and posterior are shown as dashed and hatched histograms. The latter is what is expected in absence of any information.