Calibrating the cosmic distance ladder using gravitational-wave observations

TL;DR

This paper explores how future third-generation gravitational-wave detectors can precisely calibrate the cosmic distance ladder by providing accurate distance measurements to supernova host galaxies, complementing traditional methods.

Contribution

It demonstrates that advanced gravitational-wave detectors can significantly improve distance calibration for Type Ia supernovae, reducing uncertainties and aiding cosmological measurements.

Findings

Third-generation detectors can measure distances within 0.1% to 3%.

Current LIGO/Virgo measurements have ~50% errors.

Gravitational waves can independently calibrate supernova distances.

Abstract

Type Ia supernovae (SNe Ia) are among preeminent distance ladders for precision cosmology due to their intrinsic brightness, which allows them to be observable at high redshifts. Their usefulness as unbiased estimators of absolute cosmological distances however depends on accurate understanding of their intrinsic brightness, or anchoring their distance scale. This knowledge is based on calibrating their distances with Cepheids. Gravitational waves from compact binary coalescences, being standard sirens, can be used to validate distances to SNe Ia, when both occur in the same galaxy or galaxy cluster. The current measurement of distances by the advanced LIGO and Virgo detector network suffers from large statistical errors ($\sim 50\%$). However, we find that using a third generation gravitational-wave detector network, standard sirens will allow us to measure distances with an accuracy of $\sim 0.1\%$-$3\%$ for sources within $\le300$ Mpc. These are much smaller than the dominant systematic error of $\sim 5\%$ due to radial peculiar velocity of host galaxies. Therefore, gravitational-wave observations could soon add a new cosmic distance ladder for an independent calibration of distances to SNe Ia.

Tables

Table 1. Table 1 Description of various detector networks used in this paper.

Network	Detector location	Detector sensitivity	$f_{\mathrm{low}}$ (Hz)
2G	Hanford-USA, Livingston-USA, Italy, India, Japan	aLIGO, aLIGO, AdV, aLIGO, KAGRA	10, 10, 10, 10, 1
3G	Utah-USA, Australia, Italy	CE, CE, ET	5, 5, 1
Hetero	Utah-USA, Livingston-USA, Italy, India, Japan	CE, Voyager, ET, Voyager, Voyager	5, 5, 1, 5, 5

Equations

$${\Gamma }_{ij}=⟨{\tilde{h}}_i,{\tilde{h}}_j⟩,$$

$$⟨a,b⟩=2{\int }_{f_{\mathrm{low}}}^{{\mathrm{f}}_{\mathrm{high}}}\frac{a(f)b^{\ast }(f)+a^{\ast }(f)b(f)}{S_h(f)}df.$$

$$\Delta {\theta }_i=\sqrt{C^{ii}}.$$

$${\Gamma }^{\mathrm{net}}=\sum _A{\Gamma }^A.$$

$$\Delta V\simeq D_L^2\Delta D_L\Delta \Omega \simeq 25\mathrm{Mp}{\mathrm{c}}^3{\left(\frac{D_L}{300\mathrm{Mpc}}\right)}^2.$$

$$ϵ=\frac{|D_{\mathrm{BNS}}-D_{\mathrm{SNeIa}}|}{D_{\mathrm{BNS}}},$$

Full text

1. Introduction

The geometry and dynamics of the universe can be inferred by two key ingredients obtained for a population of cosmological sources: precise measurement of their redshift and accurate estimation of their luminosity distance. The luminosity distance $D_{L}$ to a source at a redshift $z$ depends on a number of parameters such as the Hubble-Lemaître parameter $H_{0},$ dimensionless dark matter and dark energy densities $\Omega_{M}$ and $\Omega_{\Lambda},$ dark energy equation of state parameter $w$ (which may itself depend on redshift), and the curvature of space $\Omega_{k}.$ One can fit a cosmological model $D_{L}(z;\vec{p})$ to a set of, say $n$, measurements $\{D^{\ell}_{L},z_{\ell}\},$ $\ell=1,\ldots,n$ and hence determine the parameters $\vec{p}=(H_{0},\Omega_{M},\Omega_{\Lambda},\Omega_{k},w).$ It is apparent that to do so one must obtain an unbiased measurement of the distances and redshifts at cosmological scales.

1.1. Standard Candles

Distances can be measured using a standard candle — a source whose intrinsic luminosity is well constrained, so that its measured flux can be used to infer its distance. Calibration of distance to astronomical sources typically uses a “distance ladder” of multiple steps to get from nearby sources to those at cosmological distances. For example, in the most precise recent approach, nearby Type Ia supernovae (SNe Ia) are calibrated via the “standard candle” behavior of Cepheid variable stars (Riess et al., 2019). The Leavitt Law enabling determination of Cepheid luminosities from their periods is calibrated in the Milky Way galaxy via Cepheid parallaxes (Riess et al., 2018); in the Large Magellanic Cloud via observations of detached eclipsing binary systems (Pietrzyński et al., 2013); and in the “megamaser” galaxy NGC 4258 which has a known geometric distance from radio observations (Humphreys et al., 2013).

1.2. Gravitational Wave Standard Sirens

Observation of gravitational waves (GWs) has opened up the possibility of accurately measuring distances on all scales independent of the cosmic distance ladder. Indeed, binary black holes and binary neutron stars are now being used to infer both the absolute and apparent luminosity of the source: the rate at which the emitted wave’s frequency chirps up as it sweeps through the sensitivity band of a detector gives the source’s intrinsic luminosity and the measured wave’s amplitude gives the source’s apparent luminosity. Combining the two we can infer the source’s luminosity distance. The frequency evolution of the wave is completely determined by general relativity: it depends on the source’s masses and spins, which are also measured via the wave’s amplitude and frequency evolution in a network of detectors. Apart from general relativity, no detailed modeling of the source is required in this measurement.

The apparent luminosity of the source (basically the strain amplitude) depends not only on the luminosity distance but also the source’s position on the sky and the orientation of the binary’s orbit relative to the line of sight from the detector to the source. With a network of three or more detectors it is, in principle, possible to infer all the unknown parameters of the source. In practice, however, the source’s inclination is difficult to measure, especially when the orbital plane is close to face-on or face-off relative to the detector. This causes the biggest uncertainty in the estimation of luminosity distance of the source. In Sec. 2, we briefly discuss various uncertainties in the measurement of the source’s luminosity distance from their GW signal.

Gravitational wave observations should be able to calibrate all the rungs of the cosmic distance ladder for every galaxy or galaxy cluster that hosts a binary merger, and have potential to deliver new insights into the physics of these rungs. For example, one can ask if the $D_{n}$– $\sigma$ relationship, one of the rungs of the distance ladder, is metallicity-dependent. Moreover, are there systematic variations due to the inclination of the galaxy that could be resolved from GW observations? Among all the rungs of distance ladder, currently SNe Ia are the only ones that can estimate extragalactic distances at very high redshifts ($z\sim 2.26$, Rodney et al. (2015)) and have immense importance in characterizing the cosmic expansion at $z<1$ (Betoule et al., 2014; Scolnic et al., 2018). Accurate measurement of relative event-to-event distances to SNe Ia can be achieved via their well-characterized multicolor light curve shapes (Riess et al., 1996) or, near-equivalently, their peak luminosity–decline rate correlation (Phillips et al., 1999). However, while precise relative calibration suffices to characterize the recent cosmic expansion history to high-precision, thanks to the linearity of the Hubble relation at $z\lesssim 0.1$, SNe Ia can only support a Hubble constant measurement via independent distance measurements that provide an absolute calibration for their peak brightness.

1.3. Standard Sirens for Measuring $H_{0}$

Schutz (1986) noted that the standard siren property of compact binaries could be used as an independent measure of $H_{0}$ (also see Krolak & Schutz (1987)). However, the redshift $z$ to a merger event is degenerate with the binary’s total mass $M$ and it is only possible to infer the combination $(1+z)M$ from GW measurements alone111In the case of binary neutron stars, tidal effects allow the determination of the redshift of a merger event (Messenger & Read, 2012; Messenger et al., 2014) albeit measurement errors based on current methods are too large to be useful for cosmography.. Unfortunately, the sky position error-box containing a merger event typically contains thousands of galaxies (Gehrels et al., 2016; Nair et al., 2018). Assuming the merger came from any of the galaxies within the error-box would lead to multiple values of $H_{0}$ for a single merger. With a large enough population of events one gets a distribution of measured values of $H_{0}$ which will peak at its true value. This way of estimating $H_{0}$ is known as statistical method and it does not require GW events to have an electromagnetic counterpart. Alternatively, if electromagnetic follow-up observations in the sky position error-box of a merger identify a counterpart then it would be possible to directly obtain source’s redshift (Dalal et al., 2006) and hence directly infer the Hubble-Lemaître parameter. Either of these methods requires accurate knowledge of the sky position of the source, which could be obtained with a network of three or more GW detectors.

1.4. Current Status of $H_{0}$ Estimate

Cepheid-based calibration of the nearby sample of SNe Ia enables the use of their counterpart SNe Ia on cosmological scales to measure the Hubble constant (Riess et al., 2016, 2019). This approach currently gives $H_{0}=74.03\pm 1.42$ km s*-1* Mpc*-1* (Riess et al., 2019). Calibrating these same supernovae via a largely independent distance ladder based on the “tip of the red giant branch” (TRGB) approach yields $H_{0}=69.8\pm 0.8$ (stat) $\pm 1.7$ (sys) km s*-1* Mpc*-1* (Freedman et al., 2019).

An alternative geometric approach to distance measurement in the late universe, by the H0LiCOW team, uses gravitational lensing time delays and careful modeling to derive a somewhat less precise single-step measurement of the Hubble constant, $H_{0}=73.3^{+1.7}_{-1.8}$ km s*-1* Mpc*-1* (Wong et al., 2019).

Both of these $H_{0}$ values are larger than those derived from the Planck Collaboration’s observations of the cosmic microwave background (CMB), $H_{0}=67.4\pm 0.5$ km s*-1* Mpc*-1* (Aghanim et al., 2018), and from the $z\lesssim 2$ measurements of the Baryon Acoustic Oscillation (BAO) peak of the galaxy correlation function, as calibrated against the physical scale of the CMB acoustic peak. The Dark Energy Survey (DES), for example, recently reported $H_{0}=67.77\pm 1.30$ km s*-1* Mpc*-1* (Macaulay et al., 2018), while a joint analysis of several recent BAO results by Addison et al. (2018) gives $H_{0}=66.98\pm 1.18$ km s*-1* Mpc*-1*. Thus, present $H_{0}$ estimates can be divided into two categories: early universe and CMB-calibrated estimates (CMB, BAO) which tend low, and late universe estimates (SNe Ia, H0LiCOW) which tend high (with the recent TRGB estimate in between). The difference between the two classes of measurement potentially reflects new physics on cosmological scales (Riess et al., 2019; Wong et al., 2019), either at low redshift or in the early universe (Aylor et al., 2019).

After the detection of GW170817 (Abbott et al., 2017b) and identifying its host galaxy NGC 4993 as an optical counterpart, it became possible to independently estimate the value of $H_{0}$, and Abbott et al. (2017c) reported it to be $70^{+12}_{-8}$ km s*-1* Mpc*-1*. As a proof-of-principle demonstration of the statistical method, $H_{0}$ was found to be $H_{0}=77^{+37}_{-18}$ km s*-1* Mpc*-1* without using the knowledge of NGC 4993 but the distance information from GW170817 alone (Fishbach et al., 2018). Hotokezaka et al. (2019) reported an improved measurement of the Hubble constant of $70^{+5.3}_{-5.0}\,\rm km\,s^{-1}\,Mpc^{-1}$ when including an estimate of the inclination angle of the binary determined from radio observations of GW170817 (Mooley et al., 2018). Abbott et al. (2019) deployed the statical method on the population of binary black holes detected during the first and second observing runs of advanced LIGO and advanced Virgo detectors (Abbott et al., 2018) to obtain a value of $H_{0}=68^{+14}_{-7}\,\rm km\,s^{-1}\,Mpc^{-1}.$

1.5. Calibrating SNe Ia in Nearby Clusters with Standard Sirens

SNe Ia are believed to be the result of accretion induced collapse and explosion of white dwarfs. It is likely, however, that some of the SNe Ia come from mergers of binary white dwarfs instead of collapse of accreting white dwarfs (Raskin et al., 2012). Distinguishing between different subclasses of SNe Ia, or between properties of SNe Ia discovered in early-versus late-type galaxies (Jones et al., 2018), could be one of the applications of standard sirens.

If SNe Ia and binary neutron star mergers occur in the same galaxy or galaxy cluster, it is possible to directly calibrate SNe Ia luminosities with distances inferred from GW observations. It is this approach that we focus on in the present work. While it is highly unlikely for a binary neutron star merger to occur in the same galaxy as a SN Ia in a given year, every merger event in a rich galaxy cluster will typically be accompanied by multiple SNe Ia from the galaxies in that cluster. Considering only clusters rich enough to host on average one or more SNe Ia per year, we expect $\sim 3.8$ SNe Ia per binary neutron star merger host galaxy per year of optical observation from the nearest 34 such clusters (Girardi et al., 2002), located at redshifts $z<0.072$ ($D_{L}\lesssim 300$ Mpc). Thus, GW observations from binary neutron star mergers provide a unique opportunity to calibrate SNe Ia and to look for subclasses of SNe Ia, which could improve the precision of using them as standard candles.

Consistency of the Hubble diagram determined from GW and SNe Ia would confirm that calibration of SNe Ia is unlikely to have any systematic errors. On the contrary, any discrepancy in the Hubble flow determined by the two methods could point to systematics in either. One could, in principle, use the Hubble-Lemaître parameter as a proxy for distance to SNe Ia hosts and calibrate their luminosities. Such a calibration would work well on average but would not be useful for any one galaxy or galaxy cluster, as there are radial velocity departures from the Hubble flow that are unknown. Thus, it is necessary to know the peculiar velocity of the galaxy to infer the luminosity distance from $H_{0}$. However, if standard sirens and SNe Ia are both present in the same galaxy or galaxy cluster, the knowledge of the radial velocity is not needed for calibrating SNe Ia.

We note that the idea of calibrating SNe Ia using GWs distances has also been recently explored by other authors (Zhao & Santos, 2017; Keeley et al., 2019). Using only one binary neutron star merger GW170817, Zhao & Santos (2017) showed that the calibrations with both GWs and Cepheids lead to comparable SNe Ia light curves. Keeley et al. (2019), on the other hand, emphasized on the importance of combining GW and SNe Ia data sets to achieve $1\%$ accuracy in the measurement of $H_{0}.$ The current paper is built on a similar idea and shows that it is possible to calibrate local SNe Ia distances with binary neutron star mergers occurring in the same galaxy cluster within $1\%$ accuracy using 3G GW detectors (it is almost impossible to find SNe Ia and binary neutron star mergers in the same galaxy though).

2. Systematic Biases in the Measurement of Distance with Standard Sirens

In this section, we discuss various sources of systematic bias that can affect the distance measurement of GW sources.

Distance-Inclination Degeneracy

The measurement of the luminosity distance $D_{L}$ is strongly correlated with that of the inclination angle $\iota$ of the binary with respect to the line-of-sight (Ajith & Bose, 2009; Usman et al., 2018). This is because both distance and inclination, along with the sky position angles, appear together in the amplitude of the GW polarization states (see, e.g., Eqs. (2) in Apostolatos et al. (1994)). Due to this degeneracy, a face-on ($\iota=0^{\circ}$) or a face-off ($\iota=180^{\circ}$) binary far away has a similar GW amplitude to a closer edge-on ($\iota=90^{\circ}$) binary. This degeneracy can be broken to some extent by using a network having as many detectors as possible, as far away from each other on Earth as possible (Cavalier et al., 2006; Blair et al., 2008; Fairhurst, 2011; Wen & Chen, 2010). Employing accurate waveform models that incorporate higher harmonics and spin-precession also help break this degeneracy (Arun et al., 2009; Tagoshi et al., 2014; Vitale & Chen, 2018). Measuring the event electromagnetically, if the binary coalescence has an electromagnetic counterpart, partially breaks the $D_{L}-\iota$ degeneracy (Nissanke et al., 2010). Moreover, if one can constrain the orbital inclination from the electromagnetic observations (Evans et al., 2017), the uncertainty in the distance measurement is greatly reduced as we will see below.

Effect of Weak Lensing

Gravitational waves just like electromagnetic waves get lensed when they propagate through the intervening matter (Ohanian, 1974; Bliokh & Minakov, 1975; Bontz & Haugan, 1981; Deguchi & Watson, 1986; Nakamura, 1998). The dark matter distribution along the line of sight as a GW propagates from its source to the detector can amplify or de-amplify signal’s amplitude without affecting its frequency profile (Wang et al., 1996; Dai et al., 2017; Hannuksela et al., 2019). This ‘weak lensing’ results in an additional random error in the distance measurement using GWs (Van Den Broeck et al., 2010). Kocsis et al. (2006) showed that, in the case of super-massive black hole binaries, distance measurement error due to weak lensing dominates over other uncertainties leading to $\sim 6\%$ error for sources at $z=2$. This translates to $\sim 0.1\%$ error for sources in the local universe ($\lesssim 300$ Mpc) considered in this paper. We shall see below that this is less than the average error measured by a network of third generation GW detectors. Though there are proposals to remove the weak lensing effects substantially by mapping the mass distribution along the line of sight (Gunnarsson et al., 2006; Shapiro et al., 2010), degradation of parameter estimation accuracy due to weak lensing will remain an issue for some time.

Detector Calibration Errors

It is important to note that the distance measurement is also affected by the detector calibration errors (Abbott et al., 2017d). The uncertainty in the detector calibration implies an error in the measured amplitude and phase of the signal as a function of frequency. At present, the calibration error is between $5\%-10\%$ in amplitude and $3^{\circ}-10^{\circ}$ in phase over a frequency range of $20-2048$ Hz (Abbott et al., 2016b, a, 2017f, 2017g, 2017h, 2017b). As we will see in Sec. 4, the median uncertainty in the measurement of distance to neutron star binary coalescences located at distances $\sim 10-300$ Mpc is $\sim 0.1\%-3\%$, significantly smaller than the current calibration uncertainty in the amplitude. In addition to statistical errors, detector calibration may also suffer from small systematic errors. While these errors are expected to be small, there is currently no estimate of how large they might be. There is ongoing effort to improve the calibration of LIGO and Virgo detectors using alternative methods and it is expected that calibration errors will be sufficiently small to not significantly affect distance measurements (Acernese et al., 2018; Abbott et al., 2017d; Tuyenbayev et al., 2017; Viets et al., 2018; Karki et al., 2016). These alternative methods should also help in understanding the systematic errors.

In summary, GWs are ‘one-step’ standard sirens (i.e., they do not require a calibrator at any distance), and hence, can provide unambiguous measurement of distance to the host galaxies and galaxy clusters in the local universe. This implies that GWs can be used as a distance indicator to calibrate nearby SNe Ia occurring in the same galaxy or galaxy cluster as the binary merger.

In the next section we investigate how probable is it to have binary merger and SNe Ia events in the same galaxy or galaxy cluster.

3. Spatial coincident observation of a binary neutron star merger and a Type Ia Supernova event

Gravitational waves from a binary neutron star merger in the same galaxy as a SNe Ia could help calibrate the light curve of the latter and hence allow us to infer the luminosity function of SNe Ia. How likely is it to observe a binary coalescence in the same galaxy or galaxy cluster as a SNe Ia event?

The current estimates of the local ($z=0$) SNe Ia rate are in the range $[2.38,\,3.62]\times 10^{4}$ ${\rm Gpc}^{-3}{\rm yr}^{-1}$ with a median of $3.0\times 10^{4}$ ${\rm Gpc}^{-3}{\rm yr}^{-1}$  (Li et al., 2011), while that of binary neutron star mergers are $[110,\,3840]$ ${\rm Gpc}^{-3}{\rm yr}^{-1}$ with a median of $\sim$1000 ${\rm Gpc}^{-3}{\rm yr}^{-1}$ (Abbott et al., 2018). Using the SDSS $r^{\prime}$-band luminosity function of Blanton et al. (2003), the number density of galaxies in the local universe is $\approx 10^{7}$ Gpc*-3*, when integrated down to LMC-type ($0.1L^{*}$) galaxies. Hence, SNe Ia occur at roughly once every 300 years per galaxy and binary neutron star coalescences occur at a rate $\sim 30$ times smaller. Therefore, the chance of observing both of these events in a single galaxy, over a ten year period, is roughly 1 in $10^{3}$ per galaxy.

However, for every binary neutron star merger in a galaxy cluster one expects to find a number of recent SNe Ia. Although the binary neutron star merger rate in rich galaxy clusters is yet to be measured, we assume it will track the SNe Ia rate, as both populations originate in compact object mergers. Hence, we anticipate the ratio of SNe Ia and binary Neutron star merger volumetric rates $R_{\rm SNIa}:R_{\rm BNS}\sim 30:1$ (estimated 90%-confidence range of 8:1 to 300:1) will carry over to rich clusters directly. Given an SNe Ia rate in $z<0.04$ rich galaxy clusters of $R_{\rm SNIa}\sim[0.9,1.4]\times 10^{-12}\,L_{B,\odot}^{-1}$ yr*-1*, with a median of $1.2\times 10^{-12}\,L_{B,\odot}^{-1}$ yr*-1* (Dilday et al., 2010), this implies that there will be $\approx 6$ SNe Ia and $\sim 0.2$ binary neutron star mergers per year in a Coma-like cluster of total luminosity $L_{B}\approx 5.0\times 10^{12}\,L_{B,\odot}$ (Girardi et al., 2002)222Cluster SN Ia rates at $z<0.5$ in these “SNuB” units from previous surveys (for $h=0.7$; Dilday et al. 2010) are: 1.16 (Mannucci et al., 2008), 1.49 (Sharon et al., 2007), 1.63 (Gal-Yam et al., 2002), and 1.29 (Graham et al., 2008)..

In order to explore the implications of binary distance measurements for calibration of SNe Ia luminosities, we consider a catalog of the 34 nearest ($z<0.072$; $D_{L}\leq 300\,(h/0.72)^{-1}$ Mpc) galaxy clusters having luminosities $L_{B}\gtrsim 8\times 10^{11}\,L_{B,\odot}$, sufficiently rich that each is expected to host one or more SNe Ia per year. Drawing cluster identifications and luminosities from Girardi et al. (2002), with redshifts from the NASA/IPAC Extragalactic Database333NASA/IPAC Extragalactic Database: https://ned.ipac.caltech.edu, we carry out numerical simulations of the number of binary neutron star mergers observed in each cluster for active GW observing campaigns of duration 1 year to 5 years. Each simulation assumes a ratio of SNe Ia to binary merger rates of either 30:1 (median), 300:1 (pessimistic), or 8:1 (optimistic), spanning the current 90%-confidence range in binary neutron star merger rates. Uncertainties in this ratio dominate over the present uncertainty in the SNe Ia rate for rich clusters.

These simulations seek to answer two questions: (1) how many cluster distances can be calibrated by GW observation of binary neutron star mergers; and (2) how many SNe Ia luminosities can be calibrated, in turn, via these cluster distances. Results are presented in Fig. 1: The mean number of clusters with binary-neutron-star-merger-based (GW-derived) distance measurements after 5 years of GW observation is 1.8, 13.2, and 26.6 clusters (of 34 in the sample) for the pessimistic, median, and optimistic cases, respectively. The 90%-confidence ranges on these estimates are roughly $\pm$4 in the median and optimistic cases, and $\pm$1 in the pessimistic case. In the pessimistic case, it is possible (with $\approx$1.6% probability) that we do not observe any cluster that hosts any binary neutron star event even after 5 years of GW observation.

The number of SNe Ia that can be calibrated via these binary merger host clusters depends on the total duration and efficiency of any associated optical observing campaign(s) capable of discovering and characterizing SNe in these clusters. We therefore estimate the rate of calibrated SNe Ia per cluster per year of optical observation, a metric that is relatively robust both to the ratio of SNe Ia to binary neutron star merger rates (whether optimistic, median, or pessimistic), and to the duration of the GW observing campaign. To estimate the total number of calibrated SNe Ia, one multiplies the per cluster per year rate (lower panel) by the number of merger host clusters for the given GW year scenario (upper panel), and by the duration of optical observations in years. As an important caveat, we note that only SNe Ia with high signal-to-noise detections and either spectroscopy or high-quality multiband photometry (or both) will likely be useful for precise absolute calibration and Hubble constant measurement.

The main survey of the Large Synoptic Survey Telescope (LSST Science Collaboration et al., 2017) is planned to extend for ten years, and this facility will be capable of discovering and characterizing the majority of SNe Ia in most of these clusters; as a caveat, we note that not all considered clusters lie within the LSST survey area, and clusters within the survey area will still be subject to seasonal observability constraints. Similar considerations will apply when estimating the useful yield from other optical surveys seeking to characterize SNe Ia in these galaxy clusters. Overall, we consider a ten year period of observation to be reasonable for the 2030’s time frame of the GW campaigns. As seen in Fig. 1, such a ten year baseline of optical observations potentially enables calibration of $\approx 38$ SNe Ia per binary neutron star merger host cluster. In the upper panel, the number of unique binary neutron star host clusters does not increase linearly with time since we have a finite number of clusters and mergers repeatedly occur in some of the clusters; we note that multiple binary neutron star mergers in the same cluster would further improve the statistical uncertainty in the calibration of that cluster’s supernovae.

We note that 90%-confidence ranges on these numbers are larger than the Poisson error on the number of SNe Ia would suggest, because fluctuations in the number of binary neutron star host clusters with GW distance measurements typically dominates the overall uncertainty. Overall, as a robust lower bound, Fig. 1 shows that the binary merger approach can anticipate successful calibration of $>$1 SNe Ia per cluster per year of high-quality optical survey coverage, or $>$10 SNe Ia per cluster for ten years of optical observation.

In the next section, we compute the error in the measurement of distance to the nearby galaxy clusters hosting binary neutron star mergers and see how accurately we can estimate distances using various future networks of GW detectors.

4. Distance Measurement Accuracy Using Standard Sirens

Let us consider a population of binary neutron stars is uniformly distributed in the co-moving volume between luminosity distance $D_{L}$ of $10$ Mpc and $300$ Mpc. As we shall see below, for binary neutron star mergers closer than about 300 Mpc the statistical error in the distance measurement is well below systematic errors. Moreover, at such distances we can approximate the luminosity distance-redshift relation to be given by the Hubble-Lemaître law $D_{L}=cz/H_{0}$ and we don’t need to worry about cosmological effects. Also, since we will be using GWs to calibrate distance to SNe in the local universe, this distance range is more relevant.

We assume neutron stars in the binaries to be non-spinning, have fixed masses $m_{1}=1.45M_{\odot}$ and $m_{2}=1.35M_{\odot}$ and be located randomly on the sky; that is, their declination $\theta$ and right ascension $\phi$ obey uniform in $[-1,1]$ in $\sin\theta$ and uniform in $[0^{\circ},360^{\circ}]$ in $\phi,$ respectively. Further, we assume that the cosine of the inclination angle $\iota$ (the angle between binary’s orbital angular momentum ${\bf L}$ and the line of sight ${\bf N}$) is uniform in $[-1,1]$. The antenna pattern functions of GW detector also depend on the polarization angle $\psi$, which sets the inclination of the component of ${\bf L}$ orthogonal to ${\bf N}$ (see Sec. 4.2.1 in Sathyaprakash & Schutz (2009)). We choose $\psi$ to be uniform in $[0^{\circ},360^{\circ}]$. This constitutes the parameter space, $\{m_{1},m_{2},D_{L},\iota,\theta,\phi,\psi,t_{c},\phi_{c}\}$, for our target binary neutron stars, where $t_{c}$ and $\phi_{c}$ are the time and phase at the coalescence of the binary and we set them to be zero in our calculations. As binary neutron stars have long inspirals, we use 3.5PN accurate TaylorF2 waveform (Buonanno et al., 2009) to model their GWs.

Currently we have three second generation (2G) GW detectors that are operational: advanced LIGO (aLIGO) in Hanford-USA, aLIGO in Livingston-USA, and advanced Virgo (AdV) in Italy (Aasi et al., 2015; Acernese et al., 2015). The Japanese detector KAGRA (Aso et al., 2013; Somiya, 2012) is expected to join the network in the third observing run, and the detector in the Indian continent, LIGO-India, is expected to be online by 2025 (Iyer et al., 2011). Therefore, in a few years time we will have a network of 2G detectors fully operational, observing the GW sky. We call such a network of second generation detectors the “2G network”. At present, significant efforts are on-going to put forward the science case for the third generation (3G) GW detectors such as Cosmic explorer (CE) (Abbott et al., 2017e) and Einstein telescope (ET) (Punturo et al., 2010). These 3G detectors will not only let us ‘hear’ deeper in the universe, allowing more and more detections, but will also help us study each source in great detail. These 3G detectors are expected to be online sometime in 2030s. Therefore, by that time we will have a network of 3G detectors, say, ET in Italy, one CE in Utah-USA and another CE in Australia. It has been found that by placing 3G detectors on the globe in this manner, we will be able to achieve maximum science goals (Hall & Evans, 2019). We term such a network of detectors as “3G network”. Furthermore, there are also plans to improve the sensitivity of existing detectors at LIGO sites by a factor two by using high power lasers and better and bigger test masses, these are called ‘LIGO Voyager’444https://dcc.ligo.org/LIGO-T1500290/public. Given this we will have LIGO Voyager, as well, by the time 3G detectors come online. Therefore, we assume a hypothetical network of detectors constituting 3G and Voyager detectors: CE in Utah-USA, one Voyager in Livingston-USA, ET in Italy, one Voyager in India and one Voyager in Japan, and we name this as “Heterogeneous network”. Table 1 lists the detector networks used in this paper to measure binary distances, along with their location on Earth and the associated noise sensitivity curves555We use the an analytical fit given in Ajith (2011) for the power spectral density (PSD) of aLIGO. The PSD for AdV, KAGRA and Voyager are taken from https://dcc.ligo.org/LIGO-T1500293/public. For ET we use the data given in Abbott et al. (2017a) and for CE we use the analytical fit given in Kastha et al. (2018).. Figure 2 presents the cumulative distribution of network signal-to-noise ratios (SNRs) for the binary neutron star population we considered in this paper while using 2G, 3G and Hetero networks.

To measure the errors in the distance we use the Fisher information matrix technique (Rao, 1945; Cramer, 1946). This is a useful semi-analytic method that employs a quadratic fit to the log-likelihood function and derives $1-\sigma$ error bars on the binary parameters from its GW signal (Cutler & Flanagan, 1994; Arun et al., 2005). Given a frequency-domain GW signal $\tilde{h}(f;\boldsymbol{\theta})$, described by the set of parameters $\boldsymbol{\theta}$, the Fisher information matrix is given as

[TABLE]

where $\tilde{h}_{i}=\partial\tilde{h}(f;\boldsymbol{\theta})/\partial\theta_{i}$, and the angular bracket, $\langle...,...\rangle$, denotes the noise-weighted inner product defined by

[TABLE]

Here $S_{h}(f)$ is the one-sided noise power spectral density (PSD) of the detector and $[f_{\rm low},f_{\rm high}]$ are the limits of integration. The variance-covariance matrix is defined by the inverse of the Fisher matrix, $C^{ij}=(\Gamma^{-1})^{ij},$ where the diagonal components, $C^{ii}$, are the variances of $\theta_{i}$. The $1-\sigma$ errors on $\theta_{i}$ is, therefore, given as

[TABLE]

In the case of a network of detectors, one computes Fisher matrices $\Gamma^{A}$ corresponding to each detector $A$ and adds them up

[TABLE]

The error in the parameters is then given as $\Delta\theta_{i}=\sqrt{C^{ii}}$ where $C$ is now the inverse of $\Gamma^{\rm net}$.

As the chirp mass, ${\cal M}=(m_{1}m_{2})^{3/5}/(m_{1}+m_{2})^{1/5}$ and symmetric mass ratio, $\eta=m_{1}m_{2}/(m_{1}+m_{2})^{2}$ are the best measured mass parameters by GW observations during the inspiral phase of a binary, we assume our parameter space to be $\boldsymbol{\theta}=\{{\rm ln}{\cal M},{\rm ln}\eta,{\rm ln}D_{L},\cos(\iota),\cos(\theta),\phi,\psi,t_{c},\phi_{c}\}$. Fisher matrix based parameter estimation in the context of third generation detectors have been done in the past (Zhao & Wen, 2018; Chan et al., 2018). In this paper, we compute fractional error in the distance measurement, $\Delta D_{L}/D_{L}$, using the detector networks listed in Tab. 1, and the results in various observational scenarios are as follows:

Unknown sky position and inclination

: In this scenario, we assume that nothing is known about the binaries and compute errors in all the parameters using 9-dimensional Fisher matrix. This scenario is relevant when we can not identify the electromagnetic counterpart of the binary neutron stars and all the information about the source is coming from GW observation alone. We compute $1-\sigma$ error in the parameters $\{{\rm ln}{\cal M},{\rm ln}\eta,{\rm ln}D_{L},\cos(\iota),\cos(\theta),\phi,\psi,t_{c},\phi_{c}\}$ and the cumulative distribution of fractional error in the distance measurement, $\Delta D_{L}/D_{L}$, is shown on the right most panel of Fig. 3. We observe that the 3G network performs slightly better than the Hetero network, constraining distances with a median of $\sim 1.6\%$ accuracy ($90\%$ sources have error $\lesssim 3\%$). The reason behind 3G network performing better than Hetero is because 3G network has 3 third generation detectors whereas Hetero contains only 2 such detectors. The network of second generation detectors, on the other hand, performs very poorly providing distance estimates with $\sim 50\%$ error ($90\%$ sources have error $\lesssim 60\%$). On the left panel of Fig. 4, we present the distribution of $1-\sigma$ error in the measurement of cosine of the inclination angle $\iota$. Again, 3G and Hetero networks achieve similar accuracies with a median error of $\sim 0.01$ whereas 2G network performs an order of magnitude worse, constraining $\cos\iota$ with median error of $0.4$. Figure 5 presents the cumulative distribution of $90\%$ credible area of binaries on the sky. The 3G network gives the best estimate for the sky location followed by Hetero network. For instance, the 3G network will be able to locate binary neutron star merger (on an average) within $\sim 1$ deg2 whereas the Hetero network can have the $90\%$ credible sky area $\sim 1.4$ deg2, and the 2G network could only pinpoint the binary neutron stars with $\sim 180$ deg2 sky-area.

Known sky position but unknown inclination

: In this scenario, we assume that the sky position of the binary neutron stars are known through their electromagnetic observations. We, therefore, use the information of $\theta$ and $\phi$ of the sources and compute only 7-dimensional Fisher matrix for parameters: $\{{\rm ln}{\cal M},{\rm ln}\eta,{\rm ln}D_{L},\cos(\iota),\psi,t_{c},\phi_{c}\}$. The cumulative distribution of error in the distance measurement is shown in the middle panel of Fig. 3 and we notice that the accuracy has slightly improved now for all the networks. This is because the knowledge of source’s sky position breaks down the degeneracy between the sky-location angles ($\theta,\phi$) and distance $D_{L}$ and allow us to measure source distance relatively better. The 3G and Hetero networks are still performing far better than the 2G network. The right panel of Fig. 5 shows the distribution of error in $\cos(\iota)$ and it has slightly improved as compared to the case when the sky-position of the source is not known.

Known sky-position and inclination

: This scenario assumes that the sky-position as well as the inclination angles of the binary neutron stars are known purely from their electromagnetic counterparts. This scenario is possible as we already have seen in the case of GW170817. The sky position of GW170817 was constrained by finding the host galaxy NGC 4993 through numerous optical and infrared observations (Abbott et al., 2017i) whereas the inclination angle or the so-called “opening angle” was constrained from the X-ray and ultraviolet observations (Evans et al., 2017). This scenario has a merit as the error in the distance measurement can be significantly reduced as shown in the right most panel of Fig. 3. In this scenario, we use the information of $\theta$, $\phi$ and $\iota$ and compute 6-dimensional Fisher matrices for parameters, $\{{\rm ln}{\cal M},{\rm ln}\eta,{\rm ln}D_{L},\psi,t_{c},\phi_{c}\}$. All the degeneracies between the distance $D_{L}$ and $\theta$, $\phi$ and $\iota$ are now broken which give us highly accurate distance measurement with median error of $\sim 0.5\%$ for 3G and Hetero networks ($90\%$ sources have error $<0.8\%$).

Given the measurement capabilities of the different detector networks we can now assess whether it will be possible to localize a merger event uniquely to a galaxy cluster. As we shall argue unique identification of a galaxy cluster associated with a binary neutron star merger will be possible in a 3G or a heterogeneous network for 80% of the sources. From Fig. 3, left panel, we see that in the 3G (heterogeneous) network, for 80% of binary mergers the 90% credible interval in the measurement of the luminosity distance is 2% (respectively, 3%) at distances up to 300 Mpc. The corresponding 90% uncertainty in the sky position of the source is $\sim 0.1$ square degrees for both 3G and heterogeneous network (see Fig. 5), with the 3G network performing slightly better. These numbers correspond to a maximum error in distance of $\Delta D_{L}\sim 9\,\rm Mpc$ and an angular uncertainty of $\Delta\Omega\sim 3\times 10^{-5}\,\rm str,$ which correspond to an error box in the sky of

[TABLE]

Given that the number density of galaxies is $3\times 10^{6}\,\rm Gpc^{-3}$, the error box $\Delta V$ will contain no more than one field galaxy; if the merger occurs in a cluster, it will be localized to a unique cluster as the number density of clusters is far smaller than those of field galaxies. However, without an electromagnetic counterpart it will not be possible to associate a merger to a unique galaxy within a cluster, as the number density of galaxies in a cluster will be far greater than the number density of field galaxies.

In summary, given that we have restricted our analysis to rich clusters that are a sixth of Coma or larger, gravitational wave observations alone will associate most mergers in clusters to a unique galaxy cluster; an electromagnetic counterpart will be needed to further associate the event to a specific galaxy within a cluster.

5. Calibrating TYPE Ia Supernovae with Binary Neutron Star Mergers in a Galaxy Cluster

When a binary neutron star merger event occurs in a galaxy cluster we may have tens of SNe Ia in the same cluster. How do we calibrate SNe Ia in one of these galaxies given the distance to the host galaxy of the binary merger? The problem is that we would not know the relative positions of SNe Ia and binary merger host galaxy. In this section we derive the distribution of the error one would make if one assumed that both transients occurred in the same galaxy. In other words, we investigate how the dispersion of galaxies throughout the cluster might affect the distance estimation of SNe Ia calibrated through GW events in the same cluster. An additional source of error arises from the peculiar velocity of host galaxies of the transient events. In the second part of this section we provide a rough estimate of how large this effect might be.

Error due to position uncertainty of SNe Ia hosts:

To this end, we take the example of the Coma cluster. The Coma cluster is roughly 100 Mpc away from Earth and contains more than 3000 galaxies. Following several studies (Lokas & Mamon, 2003; Brilenkov et al., 2017) we assume that the matter density in Coma can be well approximated by the Navarro-Frenk-White profile (Navarro et al., 1996). To simulate positions of galaxies within this cluster we use the publicly available python-package halotools (Hearin et al., 2017) which requires the number of galaxies in a cluster, their concentration, and the mass of the cluster as input parameters. We simulate 1000 galaxies and assume the concentration and mass of the cluster to be 4 and $1.29\times 10^{15}M_{\odot}\,h^{-1}$, respectively, as reported in Brilenkov et al. (2017). We consider $h$ to be $0.701$.

In Sec. 3, we learned that ten years of optical observation would allow us to calibrate roughly 38 SNe Ia per binary neutron star merger host galaxy cluster. Furthermore, we expect to observe between 1.8 and 26.6 such clusters within 300 Mpc in five years of GW observation period. For simplicity in our calculations, we assume that all these clusters are Coma-like, i.e., they all have same matter density profile and each contains 1000 galaxies. Let us consider that one detects a binary neutron star merger in a particular galaxy cluster, it will then be accompanied by $38$ SNe Ia within 10 years of optical observation. We distribute 1 binary neutron star and 38 SNe Ia randomly among cluster’s 1000 simulated galaxies, and calculate the fractional difference $\epsilon$ in the luminosity distances of binary neutron star merger and SNe Ia as

[TABLE]

where $D_{\rm BNS}$ and $D_{\rm SNeIa}$ are the true distances of binary neutron star mergers and SNe Ia, respectively, in our simulation. With one galaxy cluster we obtain 38 samples of $\epsilon$, and since all the clusters are the same it is easy to scale this number with the number of clusters. More explicitly, having two clusters with each containing 1 binary neutron star merger and 38 SNe Ia is equivalent to have one cluster containing 2 binary neutron star mergers and 76 SNe Ia. Following this argument, in Fig. 6 we plot the cumulative distribution of $\epsilon$ for 2, 13 and 27 binary neutron star mergers in a cluster (we round the number of clusters to the nearest integer). The cyan colors show 100 realization of sampling radial positions of galaxies in Coma using halotools and the black curve represents the median. From Fig. 6 we note that 90% (99%) of the times $\epsilon<0.9\%$ ($<1.5\%$) which implies that there will be ${\cal O}(1\%)$ error in the distance estimation of SNe Ia if calibrated through binary neutron star mergers in the same galaxy cluster.

Error due to peculiar velocities of host galaxies:

In a rich cluster, galaxies can have quite a large peculiar velocity. For example, Lokas & Mamon (2003) quote that the peculiar velocity $\vec{v}_{p}$ in the Coma cluster can be as large as $\sim 10^{4}\,\rm km\,s^{-1}$, while typical rich clusters are known to have $|v_{p}|\sim 750\,\rm km\,s^{-1}$ (Bahcall, 1995). What is relevant is the peculiar velocity projected along the line-of-sight $\hat{n},$ namely $\vec{v}_{p}\cdot\hat{n}$, because it is this velocity that affects the apparent luminosity of SNe Ia and binary neutron star mergers due to the Doppler effect. For $\vec{v}_{p}$ of a constant magnitude but distributed isotropically in space we would expect the line of sight RMS velocity to be $\vec{v}_{p}/\sqrt{3}.$ However, $\vec{v}_{p}$ varies throughout the cluster, and for Coma using halotools we find $\overline{v}\equiv\langle(\vec{v}_{p}\cdot\hat{n})^{2}\rangle^{1/2}\sim 10^{3}\,\rm km\,s^{-1}$, as shown in Fig 7, where $\langle\dots\rangle$ stands for average over all directions.

The luminosity distance inferred to a binary system is affected by the local peculiar velocity. The error induced in the luminosity distance due to the RMS line-of-sight velocity $\overline{v}$ is ${\delta D_{L}}=\overline{v}/H_{0}.$ Hence, for $H_{0}=70\,\rm km\,s^{-1}\,Mpc^{-1},$ the error in binary’s distance is ${\delta D_{L}}\simeq 14\,\rm Mpc$. This is the typical error we make in the estimation of distance due to peculiar velocity and it remains the same for a cluster of given concentration. Thus, at the distance of the Coma cluster, this error is $\sim 14\%$ while it reduces to $\sim 5\%$ for clusters at 300 Mpc. As seen in Fig. 3, the error in luminosity distance of binaries due to GW measurements alone (assuming that the host’s sky position is known) is $\sim 1.2\%,$ which is far less compared to the error due to peculiar motion. However, it is comparable to the error due to the position uncertainty relative to binary neutron star merger of SNe Ia that we discussed above. Thus, the calibration uncertainty of SNe Ia up to 300 Mpc is largely due to the peculiar motion of galaxies.

However, what is the typical error in the distance measurement of the binary merger itself in these Coma-like clusters? We compute the error in the distance measurement of galaxies in Coma using different networks of detectors.666In order to sample the sky positions with respect to Earth, we assume that the center of Coma cluster is located on the sky with $\theta=27.98^{\circ}$ and $\phi=194.95^{\circ}$. Figure 8 shows the cumulative distribution of network SNR for this population of binary neutron stars in Coma for 2G, 3G and Hetero detector networks. We compute the error in binary’s distance measurement in all the three observational scenarios we discussed in the previous section and the results are shown in Fig. 9. The 3G network performs the best in constraining distances with median of $\sim 2\%$ error ($90\%$ sources have error $<3\%$) when the electromagnetic counterpart of the binary neutron star merger can not be identified. The error reduces to $\sim 0.3\%$ ($90\%$ sources have error $<0.4\%$) when both the sky-position and inclination angle are known from the electromagnetic observations. Figure 10 and 11 depict the cumulative distribution of errors in the measurement of $\cos(\iota)$ and $90\%$ credible sky area, respectively.

This shows that the error in the estimation of SNe Ia distance due to GW calibration is comparable to the statistical error in the measurement of the calibrator’s distance itself for the galaxies in the Coma cluster.

6. Discussion: Gravitational Waves as A Cosmic Distance Ladder

In this paper we explored the possibility of calibrating type Ia supernovae using gravitational waves from coalescing binary neutron stars as standard sirens. According to the current best estimates, the volumetric rate of SNe Ia is 30 times larger than binary neutron star mergers. Even so, there is a very little chance that a SNe Ia would occur in the same galaxy as a binary neutron star merger. However, when a neutron star merger occurs in a galaxy cluster it is guaranteed that more than one SNe Ia would have occurred in the same cluster within a year. As shown in Fig. 1 in a typical rich cluster within 300 Mpc, such as Coma, a binary neutron star merger will be accompanied by a few SNe Ia each year, providing ample opportunity to calibrate supernovae using standard sirens.

To accomplish this task it is necessary to control the error in the measurement of distance to merging binary neutron stars to well below the other sources of error, such as the unknown relative positions of SNe Ia and the peculiar velocity of galaxies within a cluster. One makes an error of $\sim$0.9% in distance of SNe Ia, for 90% of the supernovae, when one does not know the host galaxies of either SNe Ia or binary merger in a Coma-like cluster and assumes both of them to occur in the same galaxy. On the other hand, one makes an error of $\sim$14% due to the peculiar velocities of galaxies in the Coma-like cluster. Note that Coma is 100 Mpc away from Earth and both these errors translate to $\sim$0.3% and $\sim$5%, respectively, for galaxies at 300 Mpc. In contrast, we find that the next generation of GW detector network (one Einstein Telescope and two Cosmic Explorers) will be able to obtain distance error for the standard sirens to be less than 1% for 90% of the binary neutron star mergers whose sky position and inclination are known from electromagnetic observations within 300 Mpc. Thus, the prospect of calibrating SNe Ia using a completely independent method and establishing a new cosmic distance ladder looks bright.

SNe Ia are expected to remain a key tool for distance estimation and cosmology through the next decade and beyond. A particularly exciting near-term prospect is the ten-year LSST survey (LSST Science Collaboration et al., 2017), due to begin by the end of 2022. LSST will discover and characterize $\sim$50,000 SNe Ia per year out to redshift $z\approx 0.7$ in its main survey fields, and an additional $\sim$1500 per year out to redshift $z\approx 1.1$ in its “deep drilling” fields; of these SNe Ia, $\approx$200 per year are anticipated to have LSST data of sufficient quality to support cosmological analyses. Although spectroscopic characterization of all but a fraction of LSST SNe Ia will not be feasible, photometric analyses of the SNe Ia and their host galaxies, in the context of the sheer number of events, are expected to enable high-quality constraints on cosmology, particularly the matter density $\Omega_{m}$ and Dark Energy equation of state $w$. (For LSST’s ultimate cosmological studies, the SNe Ia analysis will be combined with weak lensing measurements of mass clustering and the growth of structure, and a cosmic scale factor analysis from the baryon acoustic oscillations feature of large scale structure, to yield joint constraints on all cosmological parameters.)

A GW-based calibration of the LSST sample of SNe Ia can be achieved at low redshift via binary neutron star detections from the jointly-observed redshift range $0.02\leq z\leq 0.07$ ($85\,{\rm Mpc}\lesssim D_{L}\lesssim 300\,{\rm Mpc}$). Over this range, binary neutron star mergers will be detectable by next-generation GW facilities, while at the same time the effects of galaxy peculiar velocities will be minimal ($<$5% per object for field galaxies). LSST simulations (LSST Science Collaboration et al., 2017) project high-quality characterization of $\approx 200$ SNe Ia per year in this redshift range, and the estimated binary neutron star merger rates are 12 to 420 (median 110) per year for this 0.11 Gpc3 volume. This suggests that a high-quality GW-based calibration of SNe Ia luminosities in the field should also be possible in the LSST era.

In conclusion, the fundamental advance considered in this paper is provided by the application of precision GW-based distance measurements (Schutz, 1986) to the calibration of type Ia SN distances – specifically, in cases where events of both types are hosted by a single galaxy cluster. Considering the broader picture, the impending realization of a longstanding astronomical dream of precise distance estimates on near-cosmological scales can be expected to yield many additional applications. For example: Precision studies of galaxy and galaxy cluster peculiar velocities; three-dimensional mapping of galaxies in the context of their host clusters and groups; and the fully tomographic use of galaxies and active galactic nuclei to characterize the gas, stellar, and dark matter contents of their host groups and clusters. Given the implications of precise distance measurements for nearly every branch of astronomy and astrophysics, a mere refinement of our present understandings would be in some sense a disappointment. We choose to hope, instead, for at least a few genuine surprises.

Acknowledgments

We thank Andrew Hearin, Aseem Paranjape, Robin Ciardullo and Rahul Srinivasan for useful discussions. We thank Christopher Messenger for carefully reading the manuscript and providing useful comments. We also thank the anonymous referee for their critical comments which significantly improved the presentation of this manuscript. AG and BSS are supported in part by NSF grants PHY-1836779, AST-1716394 and AST-1708146. BSS and BFS gratefully acknowledge support from the Science and Technology Facilities Council (STFC) of the United Kingdom. We acknowledge the use of IUCAA LDG cluster Sarathi for the computational/numerical work. This paper has the LIGO document number LIGO-P1900172.