Gravitational wave observations, distance measurement uncertainties, and cosmology

TL;DR

This paper derives analytical formulas for estimating the uncertainties in gravitational wave-based distance measurements, highlighting how detector configurations and source positions affect the accuracy of cosmological parameter inference.

Contribution

It provides new analytic expressions for statistical errors in gravitational wave distance measurements, accounting for detector network and source sky position degeneracies.

Findings

Analytic error estimates agree with Bayesian simulations.

Certain detector configurations hinder distance estimation.

Implications for Hubble constant measurement accuracy.

Abstract

Gravitational waves from the coalescence of compact binaries, together with an associated electromagnetic counterpart, are ideal probes of cosmological models. As demonstrated with GW170817, such multimessenger observations allow one to use the source as a standard siren, analog of standard candles in conventional astronomy, in order to measure cosmological parameters such as the Hubble constant. No cosmological ladder is needed to estimate the source luminosity distance from the detected gravitational waves. The error on the luminosity distance plays a crucial r\^ole in the error budget for the inference of the Hubble constant. In this paper, we provide analytic expressions for the statistical errors on the luminosity distance inferred from gravitational wave data as a function of the sky position and the detector network. In particular, we take into account degeneracy in the parameter space of the gravitational waveform showing that in certain conditions on the gravitational-wave detector network and the source sky position it may not be possible to estimate the luminosity distance of the source. Our analytic approximants shows a good agreement with the uncertainties measured with Bayesian samplers and simulated data. We also present implications for the estimation error on the Hubble constant.

Tables

Table 1. Table 1: The table contains the sky locations used for the simulated injections. First and second columns: right ascension and declination of the injection (epoch t GPS = 1187008582.0 subscript 𝑡 GPS 1187008582.0 t_{\rm GPS}=1187008582.0 , Aug 17 2017 12:36:04 UTC). Third and fourth columns: the corresponding values of σ d subscript 𝜎 𝑑 \sigma_{d} and ϵ d subscript italic-ϵ 𝑑 \epsilon_{d} , see Eq. ( 22 ). Fifth column: conditioning number C 𝐶 C of the matrix Θ Θ \Theta , defined as the ratio of its eigenvalues, namely C ≡ ( 1 + ϵ d ) / ( 1 − ϵ d ) 𝐶 1 subscript italic-ϵ 𝑑 1 subscript italic-ϵ 𝑑 C\equiv(1+\epsilon_{d})/(1-\epsilon_{d}) . Sixth column: rotation angle Δ ​ ψ Δ 𝜓 \Delta\psi in Eq. ( 21 ). The injections are performed with a fixed SNR=33

Set	$\alpha$ [deg]	$\delta$ [deg]	${\sigma }_d$	${ϵ}_d$	C	$\Delta \psi$[rad]
A	227.75	33.96	0.34	0.03	1.05	0.27
B	322.89	-55.04	0.68	0.01	1.01	0.02
C	141.62	55.05	0.68	0.01	1.03	-0.81
D	48.66	-30.00	0.37	0.21	1.54	-0.43
E	287.57	37.93	0.30	0.79	8.56	-0.82
F	193.88	-50.00	0.87	1.00	7283	-0.66
G	288.65	0.09	0.30	1.00	6201	0.09
H	108.83	0.09	0.30	1.00	8394	-0.09
I	13.70	50.00	0.87	1.00	8626	-0.92

Table 2. Table 2: The table contains the extrinsic parameters for each injection series. First column: cosine of the inclination of the orbital plane with respect to the line of sight. Second column: Polarization angle ψ 𝜓 \psi in degrees.

Series	$v$	$\psi$ [deg]
0	1.00	0
1	0.87	0
2	0.50	0
3	0	0
4	0	29.79
5	0	60.16
6	0	90.00

Table 3. Table 3: This table reports three statistical the mean d ¯ ¯ 𝑑 \bar{d} , variance σ 2 superscript 𝜎 2 \sigma^{2} and skewness γ 𝛾 \gamma for the posterior probability density function of the normalized luminosity distance d 𝑑 d of each injection. First, second and third columns: posterior obtained with LALInference. Fourth, fifth and sixth columns: posterior obtained with the CF approximant. Seventh, eighth and ninth columns: posterior obtained with the MCMC marginalization of Eq. ( 31 ).

Injection	${\overline{d}}_{\mathrm{LAL}}$	${\sigma }_{\mathrm{LAL}}^2$	${\gamma }_{\mathrm{LAL}}$	${\overline{d}}_{\mathrm{CF}}$	${\sigma }_{\mathrm{CF}}^2$	${\gamma }_{\mathrm{CF}}$	${\overline{d}}_{\mathrm{MC}}$	${\sigma }_{\mathrm{MC}}^2$	${\gamma }_{\mathrm{MC}}$
A-0	0.82	1.66$\cdot {10}^{-2}$	-0.11	0.89	8.84$\cdot {10}^{-3}$	-0.43	0.88	8.29$\cdot {10}^{-3}$	-0.43
A-1	0.88	1.98$\cdot {10}^{-2}$	0.05	1.01	1.20$\cdot {10}^{-2}$	-0.40	1.01	1.06$\cdot {10}^{-2}$	-0.44
A-2	1.07	2.00$\cdot {10}^{-2}$	2.20	1.05	1.15$\cdot {10}^{-2}$	1.90	1.10	2.47$\cdot {10}^{-2}$	2.05
A-3	0.99	8.99$\cdot {10}^{-4}$	0.20	1.01	9.35$\cdot {10}^{-4}$	0.19	1.01	9.38$\cdot {10}^{-4}$	0.19
A-4	1.13	1.46$\cdot {10}^{-3}$	0.21	1.01	1.19$\cdot {10}^{-3}$	0.21	1.03	6.36$\cdot {10}^{-4}$	0.95
A-5	1.00	9.10$\cdot {10}^{-4}$	0.21	1.01	9.11$\cdot {10}^{-4}$	0.18	1.01	9.13$\cdot {10}^{-4}$	0.18
A-6	1.00	9.07$\cdot {10}^{-4}$	0.16	1.01	9.16$\cdot {10}^{-4}$	0.18	1.01	9.21$\cdot {10}^{-4}$	0.17
B-0	0.57	3.84$\cdot {10}^{-3}$	2.39	0.89	8.77$\cdot {10}^{-3}$	-0.43	0.88	7.90$\cdot {10}^{-3}$	-0.48
B-1	0.89	1.65$\cdot {10}^{-2}$	-0.19	1.01	1.19$\cdot {10}^{-2}$	-0.40	1.01	1.02$\cdot {10}^{-2}$	-0.48
B-2	1.04	1.29$\cdot {10}^{-2}$	2.53	1.06	1.26$\cdot {10}^{-2}$	1.93	1.11	2.77$\cdot {10}^{-2}$	1.94
B-3	0.97	7.99$\cdot {10}^{-4}$	0.16	1.01	8.37$\cdot {10}^{-4}$	0.18	1.01	8.05$\cdot {10}^{-4}$	0.08
B-4	1.00	9.19$\cdot {10}^{-4}$	0.19	1.01	9.25$\cdot {10}^{-4}$	0.18	1.01	9.27$\cdot {10}^{-4}$	0.18
B-5	1.04	1.07$\cdot {10}^{-3}$	0.21	1.01	9.84$\cdot {10}^{-4}$	0.19	1.01	9.59$\cdot {10}^{-4}$	0.26
B-6	1.01	9.41$\cdot {10}^{-4}$	0.17	1.01	9.17$\cdot {10}^{-4}$	0.18	1.01	9.21$\cdot {10}^{-4}$	0.20
C-0	0.86	8.83$\cdot {10}^{-3}$	-0.43	0.89	8.76$\cdot {10}^{-3}$	-0.43	0.88	8.23$\cdot {10}^{-3}$	-0.44
C-1	1.05	1.30$\cdot {10}^{-2}$	-0.43	1.01	1.27$\cdot {10}^{-2}$	-0.40	1.01	1.13$\cdot {10}^{-2}$	-0.41
C-2	1.32	5.49$\cdot {10}^{-2}$	0.40	1.05	1.05$\cdot {10}^{-2}$	1.84	1.09	2.25$\cdot {10}^{-2}$	2.17
C-3	1.04	1.07$\cdot {10}^{-3}$	0.20	1.01	1.02$\cdot {10}^{-3}$	0.19	1.01	9.94$\cdot {10}^{-4}$	0.27
C-4	0.98	8.27$\cdot {10}^{-4}$	0.17	1.01	8.84$\cdot {10}^{-4}$	0.18	1.01	8.90$\cdot {10}^{-4}$	0.16
C-5	1.03	1.02$\cdot {10}^{-3}$	0.20	1.01	9.55$\cdot {10}^{-4}$	0.19	1.01	9.57$\cdot {10}^{-4}$	0.20
C-6	1.05	1.12$\cdot {10}^{-3}$	0.18	1.01	9.88$\cdot {10}^{-4}$	0.19	1.01	9.75$\cdot {10}^{-4}$	0.22
D-0	0.82	1.24$\cdot {10}^{-2}$	-0.30	0.89	8.74$\cdot {10}^{-3}$	-0.44	0.88	8.27$\cdot {10}^{-3}$	-0.44
D-1	1.02	1.41$\cdot {10}^{-2}$	-0.42	1.01	1.26$\cdot {10}^{-2}$	-0.40	1.01	1.14$\cdot {10}^{-2}$	-0.43
D-2	0.92	3.39$\cdot {10}^{-3}$	0.60	1.06	1.18$\cdot {10}^{-2}$	1.90	1.06	7.47$\cdot {10}^{-3}$	0.35
D-3	1.01	1.02$\cdot {10}^{-3}$	0.19	1.01	9.48$\cdot {10}^{-4}$	0.19	1.01	1.00$\cdot {10}^{-3}$	0.19
D-4	1.01	9.55$\cdot {10}^{-4}$	0.15	1.01	8.72$\cdot {10}^{-4}$	0.18	1.01	8.80$\cdot {10}^{-4}$	0.17
D-5	0.98	8.18$\cdot {10}^{-4}$	0.21	1.01	8.40$\cdot {10}^{-4}$	0.18	1.01	8.51$\cdot {10}^{-4}$	0.15
D-6	1.05	1.17$\cdot {10}^{-3}$	0.22	1.01	9.82$\cdot {10}^{-4}$	0.19	1.01	9.58$\cdot {10}^{-4}$	0.35
E-0	0.79	2.36$\cdot {10}^{-2}$	0.15	0.87	1.26$\cdot {10}^{-2}$	-0.57	0.86	1.12$\cdot {10}^{-2}$	-0.47
E-1	0.98	1.73$\cdot {10}^{-2}$	-0.39	1.00	1.71$\cdot {10}^{-2}$	-0.56	0.99	1.47$\cdot {10}^{-2}$	-0.48
E-2	1.41	1.12$\cdot {10}^{-1}$	0.16	1.27	8.65$\cdot {10}^{-2}$	0.81	1.49	9.72$\cdot {10}^{-2}$	-0.18
E-3	0.98	8.50$\cdot {10}^{-4}$	0.22	1.01	8.70$\cdot {10}^{-4}$	0.18	1.01	8.97$\cdot {10}^{-4}$	0.11
E-4	1.04	1.67$\cdot {10}^{-3}$	0.06	1.01	9.59$\cdot {10}^{-4}$	0.19	1.02	2.31$\cdot {10}^{-3}$	0.20
E-5	1.07	2.48$\cdot {10}^{-3}$	0.09	1.01	9.51$\cdot {10}^{-4}$	0.19	1.01	1.79$\cdot {10}^{-3}$	0.10
E-6	0.98	8.76$\cdot {10}^{-4}$	0.18	1.01	8.99$\cdot {10}^{-4}$	0.18	1.01	9.33$\cdot {10}^{-4}$	0.14
F-0	0.76	2.35$\cdot {10}^{-2}$	-0.64	0.86	1.36$\cdot {10}^{-2}$	-0.62	0.78	2.73$\cdot {10}^{-2}$	-0.63
F-1	0.90	4.05$\cdot {10}^{-2}$	-0.70	0.99	1.93$\cdot {10}^{-2}$	-0.61	0.90	3.77$\cdot {10}^{-2}$	-0.58
F-2	1.56	1.14$\cdot {10}^{-1}$	-0.62	1.41	1.12$\cdot {10}^{-1}$	0.16	1.54	1.10$\cdot {10}^{-1}$	-0.62
F-3	6.01	1.60	-0.69	1.01	8.79$\cdot {10}^{-4}$	0.18	6.09	1.77	-0.64
F-4	2.19	2.17$\cdot {10}^{-1}$	-0.63	1.01	9.34$\cdot {10}^{-4}$	0.20	1.62	1.19$\cdot {10}^{-1}$	-0.57
F-5	1.65	1.22$\cdot {10}^{-1}$	-0.55	1.01	9.68$\cdot {10}^{-4}$	0.19	1.91	1.40$\cdot {10}^{-1}$	-0.61
F-6	6.06	1.73	-0.71	1.01	9.23$\cdot {10}^{-4}$	0.18	6.10	1.74	-0.61
G-0	0.81	2.82$\cdot {10}^{-2}$	-0.70	0.85	1.37$\cdot {10}^{-2}$	-0.37	0.78	2.79$\cdot {10}^{-2}$	-0.61
G-1	0.86	3.42$\cdot {10}^{-2}$	-0.56	0.97	1.78$\cdot {10}^{-2}$	-0.37	0.89	3.57$\cdot {10}^{-2}$	-0.62
G-2	1.26	7.91$\cdot {10}^{-2}$	-0.52	1.26	5.15$\cdot {10}^{-2}$	0.04	1.26	7.22$\cdot {10}^{-2}$	-0.60
G-3	1.60	1.11$\cdot {10}^{-1}$	-0.62	1.01	1.06$\cdot {10}^{-3}$	0.35	1.59	1.16$\cdot {10}^{-1}$	-0.60
G-4	2.49	2.89$\cdot {10}^{-1}$	-0.68	1.01	9.73$\cdot {10}^{-4}$	0.19	2.95	1.75$\cdot {10}^{-1}$	-1.29
G-5	4.93	1.22	-0.65	1.01	1.04$\cdot {10}^{-3}$	0.19	2.41	2.67$\cdot {10}^{-1}$	-0.63
G-6	1.58	1.13$\cdot {10}^{-1}$	-0.64	1.01	1.10$\cdot {10}^{-3}$	0.37	1.59	1.16$\cdot {10}^{-1}$	-0.60
H-0	0.82	2.92$\cdot {10}^{-2}$	-0.64	0.85	1.40$\cdot {10}^{-2}$	-0.37	0.78	2.85$\cdot {10}^{-2}$	-0.59
H-1	0.87	3.55$\cdot {10}^{-2}$	-0.64	0.97	1.80$\cdot {10}^{-2}$	-0.37	0.89	3.65$\cdot {10}^{-2}$	-0.63
H-2	1.31	7.27$\cdot {10}^{-2}$	-0.69	1.26	5.16$\cdot {10}^{-2}$	0.01	1.26	7.21$\cdot {10}^{-2}$	-0.60
H-3	1.69	1.29$\cdot {10}^{-1}$	-0.65	1.02	1.27$\cdot {10}^{-3}$	0.50	1.59	1.18$\cdot {10}^{-1}$	-0.63
H-4	4.62	9.79$\cdot {10}^{-1}$	-0.61	1.01	9.36$\cdot {10}^{-4}$	0.19	2.43	2.51$\cdot {10}^{-1}$	-0.52
H-5	2.46	2.82$\cdot {10}^{-1}$	-0.74	1.01	9.01$\cdot {10}^{-4}$	0.18	2.95	1.88$\cdot {10}^{-1}$	-1.25
H-6	1.52	1.12$\cdot {10}^{-1}$	-0.58	1.01	1.03$\cdot {10}^{-3}$	0.37	1.58	1.18$\cdot {10}^{-1}$	-0.62
I-0	0.78	2.66$\cdot {10}^{-2}$	-0.67	0.86	1.38$\cdot {10}^{-2}$	-0.61	0.78	2.82$\cdot {10}^{-2}$	-0.60
I-1	0.94	4.00$\cdot {10}^{-2}$	-0.64	0.99	1.95$\cdot {10}^{-2}$	-0.60	0.90	3.78$\cdot {10}^{-2}$	-0.61
I-2	1.51	9.41$\cdot {10}^{-2}$	-0.62	1.39	1.10$\cdot {10}^{-1}$	0.24	1.54	1.07$\cdot {10}^{-1}$	-0.62
I-3	5.89	1.71	-0.61	1.01	9.07$\cdot {10}^{-4}$	0.18	6.03	1.63	-0.68
I-4	1.57	1.07$\cdot {10}^{-1}$	-0.69	1.01	8.52$\cdot {10}^{-4}$	0.18	1.73	9.84$\cdot {10}^{-2}$	-0.63
I-5	2.26	2.26$\cdot {10}^{-1}$	-0.67	1.01	1.02$\cdot {10}^{-3}$	0.21	1.62	1.19$\cdot {10}^{-1}$	-0.60
I-6	6.29	1.84	-0.59	1.01	1.07$\cdot {10}^{-3}$	0.20	5.99	1.68	-0.59

Equations

$$h_{+}(t)$$

$$h_{\times }(t)$$

$$a_0(t)$$

$$\text{and }{\phi }_0(t)$$

$${\tilde{h}}_A(f)=\frac{e^{-i{\varphi }_c}}{D}{\chi }_A(v)\tilde{k}(f;t_c,\mathcal{M}),$$

$${\chi }_{+}(v)$$

$$\tilde{k}(f;t_c,\mathcal{M})$$

$$\Psi (f;t_c,\mathcal{M})$$

$$R=\left(\begin{array}[]{cc}\cos 2\psi&\sin 2\psi\\ -\sin 2\psi&\cos 2\psi\end{array}\right),$$

$$F_a^A(\alpha ,\delta ,\psi )=R_B^A(2\psi ){\hat{F}}_a^B(\alpha ,\delta ),$$

$${\tilde{h}}_a(f)=\frac{e^{-i({\varphi }_c+2\pi f{\tau }_a)}}{D}{R}_B^A(2\psi ){\hat{F}}_a^B(\alpha ,\delta ){\chi }_A(v)\tilde{k}(f;t_c,\mathcal{M}),$$

$$\vec{s}=\vec{\eta }+\vec{h}(\mathbf{\theta }).$$

$$\mathcal{L}(\vec{s}|\vec{h}(\mathbf{\theta }))\propto e^{-\frac{1}{2}⟨\vec{s}-\vec{h}(\mathbf{\theta })⟩\vec{s}-\vec{h}(\mathbf{\theta })},$$

$$\braket{\vec{a}}{\vec{b}}\equiv 4\Re\bigg{[}\int_{0}^{\infty}\dfrac{\tilde{a}^{*}_{a}(f)b_{a}(f)}{S_{n,a}(f)}\mathrm{d}f\bigg{]}$$

$$P(\vec{h}(\mathbf{\theta })|\vec{s})\propto \mathcal{L}(\vec{s}|\vec{h}(\mathbf{\theta }))\pi (\vec{h}(\mathbf{\theta })),$$

$${\mathcal{A}}_B=\frac{e^{-i{\varphi }_c}}{D}{R}_B^A(2\psi ){\chi }_A(v),$$

$${\tilde{h}}_a(f)=e^{-2\pi if{\tau }_a}{\mathcal{A}}_B{\hat{F}}_a^B(\vec{n})\tilde{k}(f;t_c,\mathcal{M}),$$

$$⟨\vec{h}⟩\overrightarrow{h^{\prime }}=\Re \left[{\mathcal{A}}_A^{\ast }{\mathcal{A}}_B^{\prime }{\Theta }^{AB}\right]⟨\tilde{k}⟩\tilde{k}.$$

$${\Theta }^{AB}(\vec{n})={\hat{F}}_a^A(\vec{n}){\hat{F}}_b^B(\vec{n}){\kappa }^{ab},$$

$${\kappa }^{ab}={\delta }^{ab}\frac{{\int }_0^{\infty }{f}^{-7/3}{S}_{n,a}^{-1}\mathrm{d}f}{{\int }_0^{\infty }{f}^{-7/3}{S}_{n,aver}^{-1}(f)\mathrm{d}f}.$$

$$S_{n,aver}^{-1}(f)=n_d^{-1}\sum _{a=1}^{n_d}{S}_{n,a}^{-1}(f)$$

$$⟨\tilde{k}⟩\tilde{k}=4{\int }_0^{\infty }\frac{|\tilde{k}(f){|}^2}{S_{n,aver}(f)}\mathrm{d}f.$$

$$\bar{\Theta }\equiv R(2\Delta \psi )\Theta R(-2\Delta \psi ),$$

$$\Delta \psi =\frac{1}{4}\arctan \left[\frac{2{\Theta }^{+\times }}{{\Theta }^{++}-{\Theta }^{\times \times }}\right].$$

$$\bar{\Theta}^{AB}=\sigma_{d}\left(\begin{array}[]{cc}1+\epsilon_{d}&0\\ 0&1-\epsilon_{d}\end{array}\right),$$

$$⟨\vec{h}⟩\overrightarrow{h^{\prime }}=\Re \left[{\bar{\mathcal{A}}}_A^{\ast }{\bar{\mathcal{A}}}_B^{\prime }{\bar{\Theta }}^{AB}\right]⟨\tilde{k}⟩\tilde{k}.$$

$$⟨\vec{h}⟩\overrightarrow{h^{\prime }}=\frac{{\sigma }_d⟨\tilde{k}⟩\tilde{k}}{2D{D}^{\prime }}\left[⟨\vec{h}⟩{\overrightarrow{h^{\prime }}}_0+{ϵ}_d⟨\vec{h}⟩{\overrightarrow{h^{\prime }}}_{{ϵ}_d}\right]$$

$$⟨\vec{h}⟩{\overrightarrow{h^{\prime }}}_0=$$

$$({\chi }_{+}+v)({\chi }_{+}^{\prime }+v^{\prime })\cos (\phi -2{\psi }_{-})$$

$$⟨\vec{h}⟩{\overrightarrow{h^{\prime }}}_{{ϵ}_d}=$$

$$({\chi }_{+}+v)({\chi }_{+}^{\prime }-v^{\prime })\cos (\phi -2{\psi }_{+})$$

$$\phi$$

$$\rho^{2}\equiv\braket{\vec{h}}{\vec{h}}=\rho_{0}^{2}{\sigma_{d}}\Big{[}(\chi_{+}^{2}+v^{2})+\epsilon_{d}(\chi_{+}^{2}-v^{2})\cos(4\bar{\psi})\Big{]}$$

Full text

I Introduction

The first direct observations of gravitational waves (GWs) by the LIGO and Virgo collaborations Abbott et al. (2016, 2018) has opened the possibility of studying astrophysical compact objects and gravity in the strong-field regime (see e.g., Abbott et al. (2019a)). So far, the LIGO and Virgo collaborations have reported the firm detection of eleven compact binary coalescences (CBC) including ten binary black hole mergers (BBH) Abbott et al. (2018) and one binary neutron star (BNS) merger Abbott et al. (2017a) during the first two scientific runs. LIGO and Virgo have recently resumed observations with an improved sensitivity and distance reach, and began to collect new GW events, while future GW detectors such as KAGRA Somiya (2012) and LIGO-IndiaUnnikrishnan (2013) are under development.

Besides the intrinsic properties of the GW source itself, it is well known that GWs can be exploited to study cosmology Del Pozzo (2012). In particular, as was shown by Schutz Schutz (1986), CBCs are cosmological standard rulers (often referred to as “standard sirens”), as it is possible to measure their luminosity distance directly from the detected GWs. Sources in the local universe (redshifts $z\ll 1$) can then be used to determine a new, independent, measurement of the Hubble constant $H_{0}$. This type of study is particularly relevant for BNS mergers with electromagnetic counterparts for which the redshift (together with an accurate sky localization) are known. It should be noted that the value of $H_{0}$ is currently under debate Jackson (2015): Hubble Space Telescope measurements using supernovae (and Cepheids) measure $H_{0}=73.24\pm 1.74\;$\text{,}\mathrm{k}\mathrm{m};\mathrm{s}^{-1};\mathrm{M}\mathrm{p}\mathrm{c}^{-1} Riess *et al.* ([2016](#bib.bib10)), whereas cosmic microwave background measurements give $H_{0}=67.8\pm 0.9\;$\text{\,}\mathrm{k}\mathrm{m}\;\mathrm{s}^{-1}\;\mathrm{M}\mathrm{p}\mathrm{c}^{-1} Ade et al. (2016).

Other independent determinations of $H_{0}$ may be crucial to explain this discrepancy. Those include GW-based measurements. Indeed, if the GW measurement of $H_{0}$ is not consistent with the HST value (both these measurements being at low redshift), this would hint to a bias in the evaluation of the value obtained in Riess et al. (2016). The opposite case where the GW and HST $H_{0}$ estimates are compatible opens up a host of interesting ideas such as the possible modification of gravity on large scales Ishak (2019). So far, the best GW measurement of the Hubble constant is given by the joint electromagnetic and GWs observation of the BNS GW170817, combined with with binary black hole GW detections and galaxy catalogues. This measurement is $H_{0}=68^{+14.0}_{-7.0}\;$\text{,}\mathrm{k}\mathrm{m};\mathrm{s}^{-1};\mathrm{M}\mathrm{p}\mathrm{c}^{-1}$$ at $1\sigma$ confidence level The LIGO Scientific Collaboration et al. (2019). Similar values inferred from GWs observations are also given in Abbott et al. (2017b, 2019b, 2019c).

The uncertainty of the $H_{0}$ estimate is naturally related to the uncertainty on the luminosity distance $D$ of the source, since at low redshift $H_{0}\sim cz/D$. In order to anticipate which error can be expected on $H_{0}$, it is thus crucial to understand which are the main sources of error in the evaluation of $D$ from GW signals.

Some of the binary parameters have a similar effect on the GW signal, causing degeneracy in their estimation, thus leading to larger measurement uncertainties. The distance $D$ and the inclination $\iota$ of the binary orbital plane enter in the expression of the gravitational wave signal from a binary merger (to leading order) in the combination $(1+\cos^{2}\iota)/D$ and $\cos\iota/D$ for the ’+’ and ’$\times$’ GW polarizations respectively. While $D$ and $\iota$ can be individually estimated from an accurate measurement of the amplitudes of both GW polarizations, they are degenerate when the detector network is essentially sensitive to one polarization Veitch et al. (2012).

In this paper we focus on the distance measurement uncertainty for a given BNS event, and study how it varies as a function of the properties of the source — including its sky location and orientation — and of the network of gravitational-wave detectors — including the number and sensitivity of detectors in the network, and their relative alignment.

Our objective is to provide an analytic framework in which the following questions can be addressed: (i) Are there directions in the sky in which the distance measurement accuracy is significantly better? (ii) Do three or more detectors help to break the inclination/distance degeneracy? If yes, by how much does the distance error decrease?

Uncertainties on parameter estimates can be obtained from the Fisher matrix formalism in the limit of large SNR Vallisneri (2008). This paper revisits and extends the seminal calculations done by Cutler and Flanagan in Cutler and Flanagan (1994) that provides first-order (or “Gaussian”) as well as higher-order (or “beyond Gaussian”) approximations of the estimation errors. Those estimates are easily applied to binary mergers at any source sky position, detector network, etc. We obtain skymaps of the predicted errors for present and future detector networks that we compare with the errors obtained with the (computationally demanding) Bayesian sampling algorithms (see e.g., Veitch et al. (2015)) currently used to perform the astrophysical parameter estimation from the observation data. We find a good agreement except in certain sky positions that we identify and we explain the reason for the discrepancy. We stress the important role played by the degeneracy parameter $\epsilon_{d}$ defined in Eq. (22), and draw implications for sources at sky locations where this parameter is close to 1 (likelihood strongly degenerate). Finally, we discuss how our estimates propagate to the errors in the Hubble constant and deduce implications for future observations.

The paper is organized as follows: in Sec. II we introduce the data analysis background for GW. Sec. III presents the mathematical framework developed to predict the luminosity distance uncertainties. In Sec. IV we discuss several estimates for the luminosity distance uncertainties, and we then compare them to Bayesian analyses of simulated GW signals in LIGO and Virgo noise in Sec. V. Finally in Sec. VI we apply the uncertainty estimates to the specific case of GW170817 and to the Hubble constant estimate obtained from this observation.

II Data analysis background

In this section we briefly describe the GW signal, focusing on the inspiral phase, and also outline the standard parameter estimation methodology now common in all the GW literature Jaranowski and Królak (2012). In doing so we introduce the relevant notation and expressions used throughout the paper.

II.1 GW waveforms and detectors

We assume that the BNS inspirals in quasi-circular orbits, and neglect the effects of both spin and tidal deformability. Furthermore, we assume that the sky position, identified by the right ascension $\alpha$ and declination $\delta$, as well as the redshift $z$ of the source are known thanks to the observation of an electromagnetic counterpart. The GW waveform emitted by the inspiralling BNS in the wave propagation frame Maggiore (2008), is then given by

[TABLE]

with

[TABLE]

and where $\phi_{c}$ is the phase at the coalescence time $t_{c}$. We consider amplitude and phase evolutions described by the lowest order approximation in the Post-Newtonian expansion given, e.g., in Sathyaprakash and Schutz (2009). At this approximation order, the amplitude and phase only depends upon the chirp mass $\mathcal{M}={(m_{1}m_{2})^{3/5}}{(m_{1}+m_{2})^{-1/5}}$.

We define the parameter $v\equiv\cos(\iota)$ with $\iota$ the inclination of the orbital plane, i.e., the angle between the angular momentum of the binary and the line-of-sight, so that for face-on binaries $v=\pm 1$, while for edge-on binaries $v=0$ . The two GW polarizations are indicated by a capital latin index, $A,B=(+,\times)$, and in the following we use the Einstein summation convention.

By Fourier transforming Eqs. (1), and using the stationary phase approximation Finn and Chernoff (1993) one finds

[TABLE]

where we have factored out explicitly the $v$, $D$ and $\phi_{c}$ dependence, and where

[TABLE]

and

[TABLE]

In the detector reference frame — denoted by a lowercase latin index — the signal is characterized by the dimensionless strain $\tilde{h}_{a}$. We denote by $n_{d}$ is the total number of detectors, so $a=1,\ldots n_{d}$. The fractional change of length in an interferometer is obtained by projecting onto the detector arms with the two “antenna patterns” $F^{+}_{a}$ and $F^{\times}_{a}$ for the two polarizations Jaranowski et al. (1998). These depend on the sky position of the GW source (right ascension $\alpha$ and declination $\delta$), on the polarization angle $\psi$, and on time $t$. For CBC signals, $F^{A}_{a}$ are essentially constant during the entire duration of the signal, so that we can drop the $t$ dependence. The terms in $\psi$ in the antenna response functions can be factorized using a rotation matrix Cutler and Flanagan (1994)

[TABLE]

such that,

[TABLE]

where $\hat{F}^{A}_{a}\equiv F^{A}_{a}(\alpha,\delta,\psi=0)$. In the frequency domain, the signal in the detector reference frame can by obtained from Eqs. (4) and (7) as

[TABLE]

where we have added a delay term $\tau_{a}$ relative to the time travel of the GW from the center of the geocentric frame to the detector.

From Eq. (8), the signal phase depends on the two intrinsic parameters111At higher PN order, the spins should be also taken into account Sathyaprakash and Schutz (2009). $\mathcal{M}$ and $t_{c}$. There are also extrinsic parameters, namely the source sky coordinates $\alpha,\delta$ (that are assumed to be known) and four other parameters $D,\phi_{c},v$ and $\psi$ which determine the amplitude scale of the GW signal. These latter parameters are recovered from the measurement of the GW polarization amplitudes $+$ and $\times$. In principle, a noise-free signal observed in at least three detectors leads to a perfect measurement of the two polarizations $+$ and $\times$, and thus enables us to completely solve for these four variables. However, this is not true in practise, as the presence of noise and correlations between those parameters greatly limits the accuracy of their estimates, and the detector network may not be sensitive to both polarizations.

II.2 Data analysis framework

In order to understand the main sources of error on the distance $D$, we briefly summarize the canonical data analysis-framework for matched filtering. In the following, the output (data, noise, etc) associated with the $n_{d}$-detector network will be denoted by $\vec{x}=(x_{1},\ldots,x_{n_{d}})$.

The data $\vec{s}$ is modelled as a superposition of noise $\vec{\eta}$ and signal $\vec{h}(\boldsymbol{\theta})$, where $\boldsymbol{\theta}$ are the set of eight parameters introduced previously,

[TABLE]

We assume that the noise $\vec{\eta}$ is stationary, Gaussian distributed with zero mean and it is uncorrelated between different detectors222We consider detectors separated by large distances, $\gtrsim$ 1000 km so that the noise in the observable frequency bandwidth can be reasonably assumed uncorrelated.. The likelihood for a signal $\vec{h}$ to be present in the detector data $\vec{s}$ is Jaranowski and Królak (2012)

[TABLE]

where the scalar product is defined by

[TABLE]

where a sum over the detectors $a$ is understood. Here, ‘*’ denotes complex conjugation, and $S_{n,a}(f)$ is the one-sided power spectral density of the detector $a$. The posterior distribution for having a signal $\vec{h}$ is given by

[TABLE]

where $\pi$ represents the prior belief for the signal parameters. The uncertainties on the recovery of the parameters $\boldsymbol{\theta}$ is encoded in the posterior distributions. We can hence estimate the uncertainties on the distance evaluation by marginalizing over all the GW parameters with the exception of the distance.

III Calculating the likelihood for the extrinsic parameters

In the limit of high SNR, the Fisher matrix is commonly used to estimate the uncertainties in parameter space. However, the Fisher matrix approach has the tendency to underestimate uncertainties for correlated or degenerate parameters Vallisneri (2008) as is the case for the luminosity distance and the other extrinsic parameters. Higher-order error estimate (aka “beyond-Gaussian”) can overcome this limitation.

In the following we revisit beyond Gaussian calculations of Cutler and Flanagan (1994) to get a reasonably accurate approximation of Eq. (12).

III.1 Frame definition

We focus on the extrinsic parameters $\boldsymbol{\beta}=(D,\phi_{c},v,\psi)$, assuming a given sky position and intrinsic parameters for the GW event. Following Cutler and Flanagan (1994) we define the amplitude factors

[TABLE]

so that from Eq. (8),

[TABLE]

where $\vec{n}$ denotes the vector pointing to the source and is a short-hand notation for the source sky position $(\alpha,\delta)$. The scalar product $\braket{\vec{h}}{\vec{h^{\prime}}}$ in Eq. (11) between the GW signal $h$ and the template $h^{\prime}$ then becomes

[TABLE]

We use the prime to denote template related extrinsic parameters $\boldsymbol{\beta^{\prime}}$ (template and signal have the same intrinsic parameters). From Eqs. (11) and (13), the matrix $\Theta^{AB}$ reads

[TABLE]

where, for uncorrelated detectors, the $n_{d}\times n_{d}$ matrix $\kappa^{ab}$ is given by

[TABLE]

Here $S_{n,aver}$ is the average, over all the detectors, of the spectral noise densities, viz.

[TABLE]

so that the matrix $\kappa^{ab}$ is real and diagonal, with weights given by the noise power spectral density for each detector. Finally $\braket{\tilde{k}}{\tilde{k}}$ is given by

[TABLE]

The matrix $\Theta^{AB}$ is a real non-diagonal symmetric matrix which incorporates information on position of the source in the sky, as well as geometrical information of the detector network, weighted by the detectors sensitivity $\kappa^{ab}$. This matrix is independent of the polarization angle $\psi$, and characterizes the response of the detector network to the $+$ and $\times$ polarizations with $\psi=0$.

The matrix $\Theta^{AB}$ can be diagonalized by

[TABLE]

with

[TABLE]

This is equivalent to defining effective polarizations $\bar{+}$ and $\bar{\times}$ obtained from the original $+$ and $\times$ by rotating the definition of $\psi$ by the angle $\Delta\psi$. The detector response $\bar{\Theta}^{AB}$ to $\bar{+}$ and $\bar{\times}$ takes the simple form

[TABLE]

where $\sigma_{d}\equiv(\Theta^{++}+\Theta^{\times\times})/2$ and $0\leq\epsilon_{d}\leq 1$. Both those quantities can be also expressed in terms of the antenna patterns $\hat{F}^{A}_{a}(\vec{n})$.

The parameter $\epsilon_{d}$ indicates how well the effective polarization amplitudes can be measured from the detector network data Cutler and Flanagan (1994). When $\epsilon_{d}=0$, $\bar{+}$ and $\bar{\times}$ can be estimated equally well. For $\epsilon_{d}=1$, only $\bar{+}$ can be estimated. We elaborate further on this point in Sec. IV.1.

Fig. 1 shows the sky distribution of the parameter $\epsilon_{d}$ for different networks of GW detectors (we also show, for comparison, a similar plot for $\sigma_{d}$ in Fig. 2). Clearly, the areas where both polarizations can be better disentangled ($\epsilon_{d}$ small) expand when the number of detectors increases.

Assuming GW sources are equally distributed in the sky, Fig. 3 shows the probability that the source position has an $\epsilon_{d}$ value lower than a certain threshold. About 80 % of the sources have $\epsilon_{d}\lesssim 0.8$ for the five detector network HLVKI333The labels “H”, “L”, “V”, “K”, “I” refers to the detectors present in the network, and corresponds to LIGO Hanford, LIGO Livingston, Virgo, Kagra and LIGO India, respectively. whereas this is less than 10 % of the sources for the two detector network HL. Using the effective polarization amplitudes $\bar{\mathcal{A}}_{B}={\cal A}_{B}(\bar{\psi}=\psi+\Delta\psi)$, Eq. (15) takes the form

[TABLE]

On substituting Eq. (22) into Eq. (23), we find

[TABLE]

with

[TABLE]

where $\chi_{+}$ is given in (5), and

[TABLE]

Finally, using (24), the network SNR is given by

[TABLE]

where, on using (6a),

[TABLE]

For a source that is overhead, $\sigma_{d}=n_{d}/2$. If the source is also face-on (i.e., $\chi_{+}=1$), it follows from Eq. (27) that $\rho^{2}=n_{d}\rho_{0}^{2}$. Thus $\rho_{0}$ can be interpreted as the single-detector SNR for a face-on source located overhead Cutler and Flanagan (1994).

III.2 Posterior probability density

The posterior probability density $P({\boldsymbol{\beta^{\prime}}}|\vec{s})$ for the extrinsic parameters ${\boldsymbol{\beta^{\prime}}}=(D^{\prime},\bar{\psi^{\prime}},v^{\prime},\phi^{\prime}_{c})$ given detector data $\vec{s}$, can be computed from the above identities. The source sky location and intrinsic parameters $\boldsymbol{\lambda^{\prime}}=(\alpha^{\prime},\delta^{\prime},\mathcal{M^{\prime}},t^{\prime}_{c})$ are known and can be factored out, so that from Eq. (12) we have444When the parameters $\mathcal{M^{\prime}}$ and $t^{\prime}_{c}$ are unknown, this still remains true in the large SNR limit Marković (1993); Jaranowski and Krolak (1994); Cutler and Flanagan (1994)

[TABLE]

where $\pi(\boldsymbol{\beta^{\prime}})$ represents the prior probability on $\boldsymbol{\beta^{\prime}}$ and is assumed to be independent of $\boldsymbol{\lambda^{\prime}}$.

Using (9), the likelihood can be factorized as

[TABLE]

In the exponential we have neglected the terms $\braket{\eta}{h-h^{\prime}}$ as they are much smaller than the terms in $\braket{h-h^{\prime}}{h-h^{\prime}}$ in the high SNR limit. For simplicity we drop the term $\braket{\eta}{\eta}$ as it is independent of the intrinsic parameters and acts as a normalization constant. Substituting the scalar products in Eqs. (25) we obtain

[TABLE]

where $d\equiv D^{\prime}/D$ is the ratio of the template distance to the true distance.

In order to determine the luminosity distance uncertainty, we need to marginalize Eq. (31) over the extrinsic parameters $v^{\prime},\psi^{\prime}$ and $\phi_{c}^{\prime}$. Unfortunately, this cannot be done analytically as this expression depends on the true extrinsic parameters which are not known. For this reason we present several approximations for Eq. (31) in the following section.

IV Marginalized distance posterior

We first begin by briefly discussing the effect of $\epsilon_{d}$ on the likelihood probability density function.

IV.1 Relation of $\epsilon_{d}$ to the degeneracy of the estimation problem

Using the definitions of the dominant polarization frame in Eq. (23) as well as Eq. (30), the likelihood can be written in the more compact form

[TABLE]

where in this expression (and with a slight abuse of notation) the amplitudes should be understood to be expressed in terms of the normalised distance $d$. In the $\bar{\mathcal{A}}^{\prime}_{+/\times}$ plane, the likelihood is a bivariate Gaussian distribution centered on the true values $\bar{\mathcal{A}}_{+/\times}$. This is due to the fact that we are neglecting noise terms in the likelihood scalar products in Eq. 30. The $\epsilon_{d}$ parameter determines the variances for this Gaussian function along the directions $\bar{+}$ and $\bar{\times}$, namely $\bar{\sigma}^{2}_{+}=1/[2\sigma_{d}\rho_{0}^{2}(1+\epsilon_{d})]$ and $\bar{\sigma}^{2}_{\times}=1/[2\sigma_{d}\rho_{0}^{2}(1-\epsilon_{d})]$.

It follows that the level curves of constant likelihood are ellipses with semi-minor and semi-major axes proportional to $\bar{\sigma}^{2}_{+}$ and $\bar{\sigma}^{2}_{\times}$. The relative length of the curve identifying the constant likelihood surfaces is proportional to the parameter space volume that would give such likelihood values. The curve length is $\propto 1/\rho_{0}\sqrt{\sigma_{d}(1-\epsilon_{d}^{2})}$ and thus increases with $\epsilon_{d}$ but decreases with the SNR $\rho_{0}$.

We show in Appendix A that, when $\epsilon_{d}<1$, ${\cal L}$ has a unique global maximum, i.e., $\bar{\mathcal{A}}^{\prime}_{A}=\bar{\mathcal{A}}_{A}$ which corresponds to the template with $d=1$, $v^{\prime}=v$ and $\bar{\psi}^{\prime}=\bar{\psi}$. When $\epsilon_{d}\rightarrow 1$, $\sigma_{\times}\rightarrow\infty$ meaning that the likelihood becomes flat along the $\bar{\cal{A}}_{\times}$ axis — the detector network is blind to the $\bar{\times}$ polarization. As a result the maximum likelihood is not a single point but the line $\bar{\mathcal{A}}_{+}=\bar{\mathcal{A}}^{\prime}_{+}$. Many solutions for $d,v^{\prime},\bar{\psi}^{\prime}$ or equivalently many physical templates (given in Appendix A) maximise the likelihood.

Phrased differently, the estimation amounts to resolving the system of two complex-valued equations $\bar{\mathcal{A}}^{\prime}_{A}=\bar{\mathcal{A}}_{A}+\delta\bar{\mathcal{A}}_{A}$ for $A=\{+,\times\}$ and where $\delta\bar{\mathcal{A}}_{A}$ is a random perturbation of the same order of the standard deviation $\bar{\sigma}_{A}$. The system corresponds to four real-valued equations which allows one to solve for the four real-valued unknowns $(D^{\prime},\phi^{\prime}_{c},v^{\prime},\psi^{\prime})$. In the large SNR limit $\rho_{0}\rightarrow\infty$ and $\epsilon_{d}<1$, the perturbation terms vanish and the system can be exactly resolved. When $\epsilon_{d}=1$, the perturbation in the second equation is infinite. This equation cannot be solved and the system is under-determined, thus leading to many degenerate solutions. The degeneracy remains when $\epsilon_{d}$ is close to 1, i.e., when the amplitude of the perturbation is comparable to the $\bar{\mathcal{A}}_{\times}$ polarization amplitude. This happens when

[TABLE]

IV.2 Approximations of the distance posterior

We now present different expressions for the posterior probability density obtained in Eq. (31) that can be used to estimate the luminosity distance uncertainty,

[TABLE]

where the marginalised posterior is given by

[TABLE]

IV.2.1 Exact marginalized posterior for $\epsilon_{d}=0$

We first study the posterior when $\epsilon_{d}=0$, i.e., when the detector network is equally sensitive to the two effective polarizations. We obtain an exact analytical expression that can be easily evaluated numerically. This calculation is of limited use as it is obtained for a special case. It is, however, an exact expression that can serve to validate the other approximations we do.

Assuming a prior $\pi(d,\psi^{\prime},v^{\prime},\phi_{c}^{\prime})=d^{2}$ Veitch et al. (2015), the posterior probability distribution in Eq. (31) can be marginalised over $\psi^{\prime}$ and $\phi_{c}^{\prime}$ to obtain

[TABLE]

with $\zeta_{\pm}=\varphi\pm 2\psi_{-}$ and $z_{\pm}=\frac{\rho_{0}^{2}\sigma_{d}}{2d}(\chi_{+}\pm v)(\chi^{\prime}_{+}\pm v^{\prime})$. This reduces to

[TABLE]

where

[TABLE]

where $I_{0}$ denotes the modified Bessel function.

Eq. (37) does not depend on $\psi$ and $\phi_{c}$. This is a consequence of the fact that, when $\epsilon_{d}=0$, the two polarizations $\bar{+}$ and $\bar{\times}$ can be differentiated regardless of the values of $\psi$ and $\phi_{c}$.

IV.2.2 Cutler and Flanagan approximation

We now obtain an approximant for the more general case, $\epsilon_{d}<1$. Following Cutler and Flanagan (1994), we expand Eq. (31) to the 0-th order in $\bar{\psi}^{\prime}$ and $\phi^{\prime}_{c}$ around the true values $\bar{\psi}$ and $\phi_{c}$. For $\epsilon_{d}<1$, this corresponds to an expansion about the unique global maximum of the posterior. The result is

[TABLE]

This equation corrects Eq. (4.57) of Cutler and Flanagan (1994) which has an error. Fig. 4 shows this is a significant correction reaching a factor of ten times in the predicted distance uncertainty for sky locations with $\epsilon_{d}\sim 1$ in the case of the HLV detector network. The approximant Eq. (39) is computationally fast to compute and can thus be evaluated easily for many sky locations and detector network.

The final posterior $P(d,v^{\prime}|\vec{s})$ in Eq. (39) is the product of two Gaussian functions, respectively associated to the $\bar{\times}$ and $\bar{+}$ polarizations. The variance is $1/(1\pm\epsilon_{d}\cos(4\bar{\psi}))$, respectively. For $\epsilon_{d}=0$ the posterior is independent of $\bar{\psi}$, consistently to Sec. IV.2.1). We verify in Sec. V that the above approximation coincides with the exact expression obtained in Sec. IV.2.1 for the case $\epsilon_{d}=0$.

Fig. 5 shows the shape of the posterior in the $d,v^{\prime}$ plane for $\epsilon_{d}=0.2$ and different values of $v$. The posterior is symmetric for $v=0$ while increasing $v$, it becomes asymmetric. This results in tails for the marginal posterior of the luminosity distance.

These tails lead to larger errors in $D$ for large $v$ as seen in Fig. 6 where we show $\Delta D/D$ with respect to the binary inclination parameter $v$, where we have used the approximant in Eq. (31), and fixed a detected SNR of 20 for a HLV network at design sensitivity, The evolution shown in Fig. 6 for the luminosity distance is consistent to that of Fig. 1 of Vitale and Chen (2018) obtained by a Bayesian analysis of simulated data. The relative distance error increases from face-on-binaries up to a maximum for $\iota\sim 60$ deg. The maximum uncertainty increases with $\epsilon_{d}$ and the position of the maximum is also a function of $\epsilon_{d}$. Indeed, increasing $\epsilon_{d}$ corresponds to moving the source to sky locations in which the detector network is less capable of distinguishing the two polarizations.

In Fig. 7 we show the error budget for the determination of the inclination angle $\iota$. This appears to be consistent with Fig. 1 of Chen et al. (2018), predicting an error budget of about $\sim 20$ deg if the luminosity distance of the source is not constrained by any independent measure (e.g. electromagnetic observations).

Fig. 8 superimposes the expected luminosity distance error with $\epsilon_{d}$ for several detector networks at design sensitivity. The inclination is fixed to the most probable inclination for detectable binaries, namely $v=0.85$ Schutz (2011). As expected, sky patches with low $\epsilon_{d}$ values have more accurate measurements of the luminosity distance. However, has shown in Figs. 6,7, 8, the relative error for the luminosity distance and $\cos\iota$ has a weak scaling with $\epsilon_{d}$, hence even sky locations with smaller $\epsilon_{d}$ can not reduce the $d$-$\cos\iota$ correlations. These result is consistent with Usman et al. (2019) in which a similar framework have been studied without taking into account the different sensitivities of the detectors. On the other hand, five detector network allows us to achieve only a $10\%$ accuracy over the majority of the sky. Crucially, for $\epsilon_{d}\approx 1$, this approximant is no longer reliable since the likelihood is strongly degenerate as shown in App. A and discussed in Sec. V.

V Simulations and validation

In order to validate our predictions, we have carried out simulations of 63 binary neutron star mergers in simulated Advanced LIGO (Hanford and Livingston) Aasi et al. (2015) and Virgo Manzotti and Dietz (2012) data at design sensitivities, using several sky locations with different values of $\epsilon_{d}$. We compute the posterior probability on the luminosity distance using the LALInference toolkit Veitch et al. (2015); LIGO Scientific Collaboration (2018), by fixing the sky position for the GW events. We have used IMRPhenomPv2 model Hannam et al. (2014); Schmidt et al. (2015). The intrinsic parameters of the template, i.e., the chirp mass, mass ratio and merger time have been left to vary together with the extrinsic parameters $d,v^{\prime},\psi^{\prime},\phi_{c}^{\prime}$. The $63=9\times 7$ injections with SNR $\rho=33$ are divided into 9 sets, labelled from “A” to “H”, and 7 series. Each set has a fixed sky position (right ascension and declination) and detection epoch: see Table 1 for more details.

Each injection set is divided in 7 different series labelled from 0 to 6. Each series has a different combination of the extrinsic parameters $v,\psi$. For the series, the combinations of the chosen extrinsic parameters are given in Table 2.

The results of our simulations are summarised in Table 3, where we report the mean $\bar{d}$ and the variance $\sigma^{2}$, together the skewness $\gamma$ obtained for the posterior distribution of the normalized luminosity distance $d$. The marginal posterior distribution for the normalized distance $d$ obtained with LALInference is compared with the predictions obtained with the CF approximation, as well as through a direct numerical integration, by sampling from the full posterior in Eq. (31) using MCMC 555The MCMC marginalization was carried out using a Parallel tempered ensemble MCMC Foreman-Mackey et al. (2013) with 250 walkers, 3 temperatures logarithmically spaced between 1 and 10 and 10000 samples, taking into account the integrated autocorrelation time.. The MCMC marginalization of Eq. (31) is performed over two additional variables with respect to the CF approximation. It is thus expected to lead to larger credible intervals that better matches with those of LALInference.

LALInference is able to recover the SNR close to the injected value $\sim 33$.

V.1 Results for sets A to E [low $\epsilon_{d}$]

We find that the CF approximation reproduces the LALInference results when $\epsilon_{d}\lesssim 0.8$. An example is given in Fig. 9 which shows a comparison of the posterior obtained by LALInference with our predictions in the case of $\epsilon_{d}=0.03$. In this case, the mean of normalized luminosity distance is approximately $1$.

In the range of low $\epsilon_{d}$, all methods (MCMC marginalization, CF approximation and the approximation obtained specifically for $\epsilon_{d}=0$) give consistent results. We conclude that, for $\epsilon_{d}\lesssim 0.8$, the CF approximation provides a good proxy for the luminosity distance uncertainties.

The skewnesses reported in Table 3 shows that the shape of the posterior on $d$ can significantly deviates that of a standard Gaussian and symmetrical curve, including in case where the almost totally of the SNR is recovered.

For GW signals with $v>0.5$ the posterior has tails towards smaller luminosity distances, while for $v<0.5$ it has tails towards higher luminosity distances. These tails plays a role in the estimation uncertainties on $H_{0}$.

V.2 Results for sets F to I [high $\epsilon_{d}$]

Fig. 10 compares the results obtained with LALInference, the CF approximant and the MCMC marginalization of Eq. (31) in a case where $\epsilon_{d}\gtrsim 0.8$. The CF approximation is clearly no longer valid, while the MCMC marginalization method still provides a correct approximation of the posterior distribution.

Our explanation is that, in the range $\epsilon_{d}\sim 1$, the posterior is mostly determined by the prior distribution and only weakly by the data. The prior has a crucial impact on the posterior probability density function when the likelihood is degenerate with respect to the extrinsic parameters, namely when $\epsilon_{d}\sim 1$. In this case, the prior acts as a selection criterion during the sampling of the degenerate peaks. For instance face-on binaries at higher distances will be arbitrarily preferred if a standard distance prior in $d^{2}$ is used. This thus leads to a bias that propagates to the estimation of the Hubble constant (discussed in the next Section). Hence, we conclude that care needs to be taken when studying events with $\epsilon_{d}\sim 1$. No such effect occurs where $\epsilon_{d}\ll 1$ as the likelihood has no degeneracies in the parameter space.

Simulation I-2 with $\epsilon_{d}=1$ provides a good example as the final shape of the posterior probability density distribution that is essentially decided by the prior distribution which scales as $d^{2}$ We demonstrate this interpretation in Fig. 11 and Fig. 12 which show the posterior of the normalized luminosity distance obtained with the MCMC marginalization of Eq. (31) using different priors. The posterior clearly changes with the prior profile in Fig. 11 where $\epsilon_{d}\sim 1$ while the posterior changes marginally in Fig. 12 where $\epsilon_{d}=0.03$.

VI Application to GW170817 and implications for the Hubble constant estimation

We now apply the approximants derived in Sec. IV to the case of GW170817 and draw the implication on the $H_{0}$ estimation from that source. We consider the following event properties: the detection epoch is $t_{\rm gps}=1187008882.43\,s$ (Aug 17 2017 12:41:04.43 UTC) Abbott et al. (2019b); the source sky location is that of NGC 4993, i.e., R.A.=$197.45~{}\rm{deg}$, dec$=-23.38~{}\rm{deg}$ Abbott et al. (2017) and the detected SNR is $32.4$. From this information, we obtain $\epsilon_{d}\approx 0.80$ for the three-detector HLV network. This case is not symptomatic of degeneracies in the likelihood as shown by the set of injections E in the previous section.

In the following simulations, we have also fixed the distance and inclination of the source to the symmetric interval of the parameters reported in Abbott et al. (2017b), $d_{\rm BNS}=41.1^{+4.0}_{-7.3}$ Mpc and $\iota_{\rm{BNS}}=152^{+14}_{-17}$ deg. The luminosity distance posterior is obtained with the CF approximation in Eq. (39) and using MCMC marginalization of Eq. (31). Hence we estimate the symmetric intervals from the simulated posteriors, obtaining $d_{\rm CF}=40.9^{+4.6}_{-6.8}$ Mpc for the CF approximant and $d_{\rm MCMC}=40.8^{+4.6}_{-7.0}$ Mpc for the MCMC posterior integration.

Using the statistical model in Abbott et al. (2017b) and the reported recessional and peculiar velocities of NGC 4993 Abbott et al. (2017b), we determine the posterior of $H_{0}$ using the luminosity distance predictions obtained above. The results are shown in Fig. 13 and compared to Abbott et al. (2017b). The predictions fit well the observations. The Hubble constant estimation for GW170817 is $H_{0}=70.0^{+12.0}_{-8.0}\;$\text{,}\mathrm{k}\mathrm{m};\mathrm{s}^{-1};\mathrm{M}\mathrm{p}\mathrm{c}^{-1} Abbott *et al.* ([2017b](#bib.bib14)), while our predicted value is $H_{0}=69.9^{+13.0}_{-7.5}\;$\text{\,}\mathrm{k}\mathrm{m}\;\mathrm{s}^{-1}\;\mathrm{M}\mathrm{p}\mathrm{c}^{-1} at $1\sigma$ confidence level by sampling from Eq. (31) and $H_{0}=69.9^{+11.2}_{-8.7}\;$\text{,}\mathrm{k}\mathrm{m};\mathrm{s}^{-1};\mathrm{M}\mathrm{p}\mathrm{c}^{-1}$$ using the corrected CF approximation. The confidence level intervals are comparable with the ones obtained from GW170817 observation.

We simulated the $H_{0}$ measurement from a “GW170817-like” event assuming a source sky location where $\epsilon_{d}=0$. For this simulation we used the parameters of GW170817 except for the sky location. We used the CF approximant and Eq. (31) to obtain the Hubble constant posterior.

The resulting $H_{0}$ posterior in Fig. 14 should be compared with Fig. 13. The CF approximation for $\epsilon_{d}=0$ predicts $H_{0}=70.1^{+10.3}_{-7.0}\;$\text{,}\mathrm{k}\mathrm{m};\mathrm{s}^{-1};\mathrm{M}\mathrm{p}\mathrm{c}^{-1} while the MCMC marginalization of Eq. ([31](#S3.E31)) predicts $H_{0}=69.2^{+13.3}_{-6.6}\;$\text{\,}\mathrm{k}\mathrm{m}\;\mathrm{s}^{-1}\;\mathrm{M}\mathrm{p}\mathrm{c}^{-1}. The measurement of the Hubble constant is thus 15% more accurate if the source sky position has $\epsilon_{d}\approx 0$ for the CF approximation compared to the actual source location with $\epsilon_{d}\approx 0.8$. The accuracy obtained with the MCMC marginalization is the same in both the $\epsilon_{d}\approx 0$ and $\epsilon_{d}\approx 0.8$ cases.

VII Conclusions

In this paper, starting from the framework introduced in Cutler and Flanagan (1994), we have developed several analytic approximants of the posterior of the luminosity distance for localized binary neutron star events.

To do so, we initially demonstrated the important role of the detector-network and sky-position dependent parameter $\epsilon_{d}$, which quantifies the ability of the detector network to disentangle the extrinsic parameters encoded in the gravitational-wave polarizations. Then we studied the degeneracies of the likelihood for GW events in absence of noise. We showed that for $\epsilon_{d}<1$, the likelihood admits one global maximum which corresponds to the true parameters of the event, while for $\epsilon_{d}=1$, the likelihood admits many degenerate global maxima, as it is not possible to disentangle the two gravitational-wave polarizations.

We obtained different approximants of the likelihood from Eq. (31). One of these approximants was already presented in Cutler and Flanagan (1994), though with a small (but significant) error. The error on the extrinsic parameter estimation predicted by the approximants were compared with those obtained for 63 simulated gravitational-wave signals using a Bayesian sampler available in the LALInference software package. This comparison confirms that, for sky locations with $\epsilon_{d}\sim 1$, the likelihood is degenerate, and this can lead to a prior-induced bias in the distance estimate.

We verified that GW170817 is localized in a sky patch with $\epsilon_{d}\approx 0.80$, meaning that we do not expect the likelihood for this GW event to be uninformative. We have checked that the posterior on the luminosity distance obtained with our approximants is compatible with the one reported in Abbott et al. (2019b). We checked that the Hubble constant estimates based on the statistical model detailed in Abbott et al. (2017b) and obtained using both our approximations and the published posterior samples are consistent.

The new approximants we have developed can be employed to forecast and study the accuracy of the $H_{0}$ measurement in different observing scenarios. For instance they can be employed in a framework similar to Mortlock et al. (2018) in order to evaluate scenarios with different detector networks and duty cycles. Moreover, the $\epsilon_{d}$ parameter provides a simple and quick indicator to understand whether the luminosity distance posterior is informative or not. This parameter can also be used to evaluate or optimize the geometry and position of an additional detector to reduce the regions in the sky where the likelihood is degenerate.

Acknowledgements

S. M. acknowledges the UnivEarthS Labex program (ANR-10-LABX-0023) for financial support. K. L. and D. S. are grateful to CERN and the University of Geneva (D.S) for hospitality whilst this work was in progress. We thank C.-J. Haster, K. Chatziioannou, S. Vitale, W. del Pozzo, V. Raymond, B. Malte Schäfer, E. Porter, M. Barsuglia, I. Magana and H.-Y. Chen for useful discussions and for helping with the LALInference runs. We are grateful to the LIGO Scientific Collaboration for sharing the LAL software suite LIGO Scientific Collaboration (2018). The authors are also grateful for the computational resources provided by the LIGO Laboratory.

Appendix A Degeneracy of the maximum likelihood template

A.1 Case of $\epsilon_{d}<1$

We define two templates to be degenerate with respect to a GW signal if they both maximize the likelihood in Eq. 10, or equivalently Eq. (32).

When $\epsilon_{d}>0$, Eq. (32) is maximized when $\bar{\mathcal{A}}_{+/\times}=\bar{\mathcal{A}}^{\prime}_{+/\times}$. Writing the two effective polarizations in terms of the extrinsic parameters $D,\psi,v$ yields the two equations (for simplicity of notation, we reabsorb the phases $\phi_{c}$ and $\phi^{\prime}_{c}$ into a redefinition of $\bar{\psi}$ and $\bar{\psi}^{\prime}$ respectively)

[TABLE]

whose sum gives

[TABLE]

The two solutions of this equation are

[TABLE]

where

[TABLE]

Note that if $\Delta>0$ then possibly a degenerate template with $d\neq 1$ would exist, while if $\Delta<0$ then no degenerate solutions exist. Clearly in the present case, it follows from Eq. (40) that $\Delta$ is negative unless $\psi_{-}=0$ and consequently $v^{\prime}=v$. The other solution $v^{\prime}=1/v$ is permitted only for $v=1$ since $(|v|,|v^{\prime}|)\leq 1$. Hence we obtain the solution $d_{\pm}=1$, which is the correct luminosity distance.

To conclude, (at the 0 PN order and for $\epsilon_{d}<1$) no degenerate templates which are able to maximize the likelihood exist in the parameter space. Furthermore, if noise terms are neglected, the only template which is able to maximize the likelihood is the one corresponding to the GW signal present.

A.2 Case of $\epsilon_{d}\approx 1$

Let us assume that $\epsilon_{d}=1$, so that the term in $\bar{\mathcal{A}}_{\times}$ vanishes in the likelihood (see Eq. (32)). Hence the likelihood is maximised when $\bar{\mathcal{A}}_{+}=\bar{\mathcal{A}}^{\prime}_{+}$ leading to two equations (from the real and imaginary parts) 666For simplicity, we have reabsorbed the phases $\phi_{c}$ and $\phi^{\prime}_{c}$ into a redefinition of $\bar{\psi}$ and $\bar{\psi}^{\prime}$ respectively.:

[TABLE]

On taking the sum of the square of the two previous equations, we obtain a solution for $d$:

[TABLE]

Since Eq. (43) is positive definite, it represents a physical distance. However, to prove that Eq. (43) represents a physical template, we substitute back into one of the two Eqs. (41)-(42) and obtain

[TABLE]

Since the RHS in Eq. (44) is $<1$, there are degenerate templates that maximize the likelihood.