Cosmography with next-generation gravitational wave detectors

One promising approach to obtaining the redshifts is to look for the EM counterpart of GW sources. With the finite number of GW observatories across the globe, the available observatory baselines limit the 2-dimensional resolutions to the sky location of the GW sources [28–31]. However, the observations of EM counterparts could improve the sky localization of the GW sources and potentially lead to the identification of the host galaxies. The redshifts can then be measured from the spectra of the hosts. This is known as the bright siren approach.

Kilonova AT 2017gfo and GRB 170817A [32–41] were the first EM counterparts found to be associated with a GW source, GW170817 [42]. With the precise localization, the host galaxy NGC 4993 was identified and the redshift could be determined ($z\sim 0.009$). This event led to a Hubble constant measurement of $H_0 = 70^{+12.0}_{-8.0}$ km s⁻¹ Mpc⁻¹ [20]. With Advanced LIGO-Virgo-KAGRA (LVK) sensitivities, the Hubble constant can be measured to ${\sim} 2\%$ when combining ${\sim} 50$ BNS bright sirens [43–46]. This precision could shed light on the tension between different Hubble constant measurements [15]. Over last few years, a variety of other types of EM emissions were proposed theoretically [47–66] or suggested observationally [67–70] to be associated with BNS, NSBH, and BBH mergers, although they could not be confidently verified. If confirmed, the applications of the bright siren approach will be further widened [71–73].

Over the last few years, the systematic uncertainties associated with bright sirens have been widely studied, such as the systematics due to the GW instrumental calibration uncertainties [74–76], the EM observation selection effect [77–79], the biased EM-inferred binary viewing angle [77], the peculiar velocity of the hosts [80–82], and the GW instrumental non-stationary noise [83–87]. A comprehensive study of the systematics and the developments of the mitigation methods are critical to ensure the value of standard siren measurements in cosmology, especially when dealing with the large and precise catalog of bright sirens that XG detectors will provide.

While the bright siren approach can lead to precise measurement of cosmological parameters, the expected number of systems with detectable EM counterparts is much smaller than those without EM signatures. In fact, since its inception, the LVK network has detected about a hundred mergers that were not associated with EM observations and only one, GW170817, for which the EM counterparts were observed. This reflects the current selection biases that favor BBH detections over BNSs or NSBHs, as well as the difficulty to follow up distant events with poor localiztions. Fortunately, several approaches can utilize the mergers with no EM counterparts to infer the cosmological parameters. For this reason, such mergers are also referred to as dark sirens.

The various approaches that employ dark sirens for cosmological inference differ in their methodology for obtaining the corresponding redshift measurement that is associated with these events. In this section, we restrict ourselves to the approaches that use galaxy catalogs to statistically infer the redshift of the sources [18]. An alternative technique relies on the fact that the GW sources as well as the galaxy population are expected to follow the large-scale structure. Hence, the cross-correlation between their spatial clustering can be used to infer cosmological properties [88–91]. This cross-correlation technique was applied to eight well-localized BBH events from the Gravitational-Wave Transients Catalog-3 (GWTC-3) [92] to get $H_0 = 68.2^{+26.0}_{-6.2}\,\,\mathrm{km}\,\mathrm{s}^{-1}\,\mathrm{Mpc}^{-1}$ [93]. LVK detectors at design sensitivities, together with dark energy spectroscopic instrument (DESI) [94] and SPHEREx [95] data, have been shown to constrain H₀ to $\mathcal{O}(1)\%$ with 5 years of observation [96].

The statistical dark siren method exploits galaxy catalogs to identify the galaxies that lie within the localization volume associated with the GW signal (note that this localization volume is obtained purely from the GW detection [31]) and obtain their redshifts. Together with the probability distribution of d_L inferred from the GW signal, the redshifts for each of these potential hosts give corresponding probability distributions on cosmological parameters. The uncertainty in the knowledge of the true host galaxy is accounted for by statistically averaging over the redshifts from all potential host galaxies, which results in constraints on cosmological parameters using a single GW event. Such an analysis can be performed for multiple GW events, which will result in the inference of cosmological parameters with greater precision [44, 97–99].

GW170814 [100] was the first BBH merger to which this technique was applied. It was localized to a volume of $2\times10^{6}\,\mathrm{Mpc}^3$ [101] which, according to the dark energy survey (DES) data, contained ${\sim}77\,000$ galaxies [102]. In comparison, the localization volume of the best localized BBH, GW190814 [103], was ∼100 times smaller and contained ∼2000 galaxies [104]. Using GW170814 and the DES data [105] for redshift information, H₀ was found to be $75^{+40}_{-32}\,\,\mathrm{km}\,\mathrm{s}^{-1}\,\mathrm{Mpc}^{-1}$ [102]. When applied to the eight best-localized events in GWTC-3, stronger constraints are obtained on H₀, measured to be $79.8^{+19.1}_{-12.8}\,\,\mathrm{km}\,\mathrm{s}^{-1}\,\mathrm{Mpc}^{-1}$ [101]. The statistical siren approach can be applied to bright siren events as well. However, in general, this method is expected to provide weaker constraints than the bright siren technique due to the uncertainty in host galaxy identification. Fishbach et al [106] applies the statistical siren strategy to GW170817 to obtain $H_0 = 76^{+48}_{-23}\,\,\mathrm{km}\,\mathrm{s}^{-1}\,\mathrm{Mpc}^{-1}$. But, as mentioned in section 2.1, due to unique host galaxy identification made possible by EM counterpart detection, the bright siren method provides much stronger bounds on H₀. The constraints obtained from the statistical dark siren and the bright siren approaches were combined for events in GWTC-3 to estimate H₀ with ∼15% uncertainty [101, 107]. Note that this measurement is dominated by the bright siren constraints from GW170817.

The advent of XG observatories will see an improvement in the measurement precision of both luminosity distance and sky-area [7, 108–111]. This will result in smaller localization volumes containing smaller number of potential hosts compared to currently detected events [112]. Among these dark siren mergers, there may be some that are so well localized in the sky that only one potential host galaxy is found in that localization volume [113, 114]. This would allow unique identification of the host galaxy and lead to a very precise estimation of cosmological parameters. Such mergers are referred to as golden dark sirens. Asymmetric (i.e. the masses of the compact objects are appreciably different) and face-off (i.e. the inclination angle $(\iota)$ is not zero) systems are touted to be ideal candidates for such sirens [115, 116]. This is due to the excitation of higher-order modes in the GW waveform for asymmetric and face-off systems [117–119], leading to the breaking of the degeneracy between $d_\textrm{L}-\iota$, which leads to the precise estimation of parameters [120–125]. Similar expectations hold for precessing systems in which the spins are not (anti-)aligned and there could be precession of the orbital plane [125–127]. Because golden dark sirens are typically the loudest events, their measurement of the luminosity distance will be extremely precise, outperforming those of bright sirens.

Several studies have looked at prospects of measuring cosmological parameters with XG GW observatories. Muttoni et al [112] show that the application of the statistical siren approach to the binaries detected by a future network with 1 ET and 2 CE observatories will yield a 0.8% bound on H₀ (with $90\%-$confidence) in one year of observation. However, they also point out that this bound is dominated by the detection of golden dark sirens. In fact, [115] claim that just by using BBH golden dark sirens, a network with A+ sensitivities can ascertain H₀ to $\mathcal{O}(1\%)$ precision (with $68\%-$confidence) and an ET+2CE network can measure H₀ to $\mathcal{O}(0.1\%)$ (with $68\%-$confidence) in two years of observation. Gupta [116] finds that observation of NSBH mergers are not expected to provide any golden dark sirens unless the network contains at least 1 CE or ET. With an ET+2CE network, the constraints on H₀ from NSBH golden dark siren mergers can range from $\mathcal{O}(0.1\%)-\mathcal{O}(1\%)$ (with 68%-confidence) in 2 years of observation, based on the local merger rate of NSBH systems. Due to the ability of NSBH mergers to act as effective dark sirens and, potentially, bright sirens [116, 128, 129] earns them the name grey sirens [116]. In section 4.1 we present our own forecasts for dark siren cosmology.

As the galaxy catalogs are biased towards bright galaxies due to apparent magnitude thresholds, the dark siren method is also susceptible to the incompleteness of galaxy catalogs that will increase the statistical uncertainty of the measurements [44, 98, 130, 131]. Gray et al [130] finds that H₀ can be constrained to ∼4.5% (with $68\%-$confidence) with 249 BNS mergers (treated as dark sirens) in Advanced LIGO sensitivity, with a galaxy catalog that has only 50% completeness up to a distance of 115 Mpc. For certain systems that are very well localized (golden or otherwise), targeted follow-up by telescopes to identify potential host galaxies can be a viable solution to tackle catalog incompleteness.

The dark siren approach is also prone to systematic uncertainties. As this approach relies on the accurate estimation of $d_\textrm{L}$, it is significantly affected by the uncertainty introduced due to calibration errors, which can affect the amplitude measurement by as much as 7% [75]. Fortunately, these errors are expected to reduce in the future. The redshift estimation can be affected by peculiar velocity measurements, especially for nearby hosts. As an example, the uncertainty in the peculiar velocity for the host of GW170817 (located at $z\sim0.01)$ was 150 km s⁻¹, leading to a redshift uncertainty of about 5% [20]. Particularly, this can influence the golden dark siren estimates, which, due to the stringent selection criteria, are expected to occur in the region $z \unicode{x2A7D} 0.1$ [115, 116]. There are also possible systematics associated to the assumption that binaries are tracers of certain cosmic structures, for example the stellar/halo mass, although these GW events themselves could also be turn into probes of the host galaxy properties [132]. Moreover, if not included in the inference as in the spectral siren method described in the next section, the assumptions about the population properties are also important source of systematic errors [133].

Even in the absence of EM counterparts or galaxy catalogs, GWs alone are able to constrain the cosmic expansion when studied collectively. The reason is that all GWs signals are stretched by the expansion, redshifting its characteristic chirping. As a consequence, we can only detect the redshifted masses of the binary from the frequency evolution of the signal, i.e. the detected mass is

$$\begin{equation} m_\mathrm{det} = \left(1+z\right)m\,, \end{equation} \tag{ 3 }$$

where m is the mass in the source frame. Since all binaries at the same luminosity distance will be redshifted by the same amount, if a mass scale is found in the data of the mass spectrum or can be predicted by (astrophysics) theory, then this reference mass scale can be used to reconstruct the redshift at that distance. This finding of spectral features that carry redshift information can be thought in analogy to the absorption or emission lines of galaxy spectra. In this sense, a catalog of GWs is indeed a powerful set of spectral sirens.

For the same source population, different expansion rates will predict different distributions of detector frame masses. For instance, as schematically shown in figure 2, a larger value of the local expansion rate H₀ will imply a larger redshift in order to end with the same luminosity distance, making the detector frame distribution shift faster to higher values. Therefore, one can constrain cosmological parameters with the population analysis of spectral sirens. Since the cosmological redshift is a common effect to all mergers at equal distance, the better the mass spectrum can be constrained over a larger range of distances, the more information can be extracted about cosmology. One should be cautious, however, because the mass spectrum is expected to evolve with redshift and this could bias the cosmological inference if not properly accounted for, as we discuss later.

Zoom In Zoom Out Reset image size

The spectral siren method operates at its full capacity when there are clear features in the mass distribution. BNSs are expected to have cutoffs in their mass distribution from the existence of a maximum neutron star mass beyond which the star becomes a black hole and a minimum mass below which a white dwarf is formed. This narrow mass distribution made BNSs the first candidates for the spectral siren method [134, 135]. Driven by observations, the mass distribution of BBHs have unveiled some interesting features. In particular there is a pronounced dearth of BBHs at high masses [136, 137] and an excess over a simple power-law at ${\sim}35M_\odot$ [138]. At first sight, these features seem to coincide with the expectations from the theory of pair instability supernova (PISN) [139–144], which predicts a gap in the BBH mass spectrum between ∼50–$120M_\odot$ and a potential pile-up just before the lower edge of the gap [145]. This connection has been however recently revisited [146], and more observation are certainly needed to resolve the origin of these features. Irrespectively of their origin, these features stand out as clear targets for the BBH population to do spectral siren cosmology with current-generation detectors [147] and have also been explored in the context of XG interferometers [148]. With the latest LVK catalog, spectral siren constraints on H₀ are ${\sim} 20\%$ at 1σ [133]. It is to be noted that if one assumes some particular shape or location for these features, perhaps from some theoretical astrophysical prior, but those do not represent well the data, the cosmological parameters could be wrongly constrained.

The spectral siren method is in its own nature data driven: the cosmology is jointly inferred with the population properties which are directly obtained from the data. It does not matter how the binaries are astrophysically formed or what they are made of, it only matters if a global shift of the mass spectrum can be identified over a wide range of distances and distinguished from local changes of its shape. In this sense all mergers can be considered simultaneously, from light neutron stars to heavy black holes, and cosmology is constrained with the full mass distribution of compact binaries [149]. This is advantageous because different sectors of the mass spectrum are expected to evolve differently with redshift. For example, the masses of neutron stars and low mass black holes are thought to be much less sensitive to the environment than those of heavier black holes whose possible maximum mass depends significantly on metalicity. Therefore, by studying the whole mass spectrum together, with its many individual features, one can break the possible degeneracies between the shift induced by the cosmological redshift and the change in the shape of the mass spectrum due to the time and environment dependence of the source population [149]. If unaccounted for, however, the evolution of the mass spectrum can bias the inference of cosmological parameters [149–151], a potential problem especially relevant for XG detectors observing far away binaries. Mismodeling of the true distribution with oversimplified parametric models [151] or wrong astrophysical assumptions that do not fit well the data can also be problematic. This is why the non-parametric reconstruction of the mass spectrum stand out as a promising avenue for spectral siren cosmology with future observations, bypassing the need to choose a specific mass function and just relying on shape extracted from the (large) data sets [152].

Finally, it is important to note that our sensitivity to different sectors of the compact binary mass spectrum will evolve in parallel to the upgrades of the GW detector network. While currently the inference of cosmology is dominated by the excess of BBHs at ${\sim}35M_\odot$ and the lack of them at ${\gt}45M_\odot$, spectral siren cosmology with XG detectors will be dominated by BNSs and low-mass BBHs, as their intrinsic merger rate is much larger [149]. Different populations will uncover the expansion rate at different redshifts, as their intrinsic merger rate history and detectability highly depends on their masses. Interestingly, with detectors at A+ sensitivity and more so with XG, we could also detect BBHs on the far-side side of the PISN gap (if they exist), thereby resolving the upper edge of the PISN gap and providing an extra anchor to probe cosmology [153].

The ability to generate EM counterparts and the constrained mass spectrum make neutron star mergers effective bright and spectral siren candidates. But that is not all—neutron stars that are part of compact binary coalescence (in either a BNS or an NSBH) undergo tidal deformation during the later stages of the inspiral. These deformations, at leading order, affect the phase evolution at the fifth post-Newtonian order (i.e, they are suppressed by a factor of $(v/c)^{10}$ compared to the dominant phase term, where v is the orbital speed and c is the speed of light) [154, 155]. The leading order tidal contribution is a linear function of the individual tidal deformability parameters of the two neutron stars, Λ₁ and Λ₂, and is called the reduced tidal deformability parameter $(\widetilde{\Lambda})$. If the neutron stars are spinning, the individual Λs will also appear in the contributions to the phase from spin-induced quadrupolar deformations [156–159]. These changes in the phase allow the measurement of Λ₁ and Λ₂ when GWs are detected. Now, given an equation of state (EOS) for neutron stars, one can relate the individual Λs to source-frame masses of the corresponding neutron stars. These source-frame estimates, along with the detector-frame redshifted masses estimated from the GW signal, can be used to ascertain the redshift associated with the binary [160]. The post-merger signal of the BNS merger can also be used to break the redshift-mass degeneracy, by utilizing the robust spectral features that depend on the source-frame masses of the binary [161]. These approaches that use the tidal deformability measurement and spectral features of the BNS system and the knowledge of the neutron star EOS to measure the redshift and obtain bounds on cosmological parameters constitute the Love siren approach (as an ode to the neutron star tidal Love number [162]).

Note that, in practice, only $\widetilde{\Lambda}$ and not the individual Λ are expected to be precisely measured, even with XG observatories [163]. However, there exist phenomenological relations, called quasi-universal relations [164], that are independent of the EOS governing the inspiralling neutron stars and allow the estimation of individual Λs using $\widetilde{\Lambda}$ and the mass-ratio of the binary. While this approach can indeed provide narrower bounds on the estimates of Λ₁ and Λ₂, they come at the cost of incurring biases due to uncertainties in the quasi-universal relations [164, 165].

Using the inspiral technique for a BNS system in ET sensitivity, the redshift can be measured to a precision of $\mathcal{O}(10\%)$ for systems up to $z\sim1$ [160, 166]. With the post-merger spectra, similar precision is obtained up to $z\sim0.04$ [161]. Dhani et al [167] show that with the A+ network, less than 10 detected BNS every year will be able to constrain the redshift to better than 10% precision with the Love siren approach. However, with an ET+2CE network, the redshift can be resolved to ${\sim}1\%$ for $\mathcal{O}(10^4)$ yearly detections. This translates to 0.1% error in H₀ and 0.61% error in $\Omega_m$ with an ET+2CE network in a year of observation [167]. Haster et al [168] shows that the XG observatories can constrain the EOS with BNS sources around $z\sim1$.

The Love siren method provides a GW-only approach to the measurement of the cosmological parameters. It relies on the knowledge of the neutron star EOS, which is expected to be well determined by the XG era [7, 169, 170]. So, either a previously inferred EOS can be used directly, or the EOS and the cosmological parameters can be jointly inferred for the GW data [171]. The forecasts that involve the post-merger spectrum require additional care, as they would strongly depend not only on the number of high-SNR BNS detections, but also on the post-merger modeling of the hypermassive neutron star remnant. Another aspect that requires caution, as noted earlier, is the use of quasi-universal relations to obtain individual Λs from $\tilde{\Lambda}$ and the mass ratio. The fits for the universal relations themselves have uncertainties that can propagate to inference of cosmological parameters, inducing biases, especially at high precision. These uncertainties can be addressed by marginalizing over the residuals [172–174] or by specifically correcting for the incurred errors [165].

The detection of GW from compact binary mergers has opened a new window to the Universe. The ∼100 events detected by the end of the LVK third observing run have resulted in constraints on deviations from general relativity [175], inference of cosmological parameters [107], revelations about neutron star properties [32, 33, 176–179] and measurement of the population properties of compact binary mergers [180]. However, based on the current estimates of the local merger rates and the redshift distribution of compact binary mergers, more than a million such events are expected to occur every year within 600 Gpc $(z \sim 50)$ [7], most of which are missed by the GW networks at current sensitivities. Of the ∼100 events that are detected, the majority lie close to the SNR threshold of 10, which does not allow the precise measurement of cosmological parameters, deviations from general relativity, or neutron star properties. Thus, to realize the full potential of GW, bigger and more sensitive detector networks are needed.

The first such advancement, called the A+ concept [4, 5], is planned for the second half of 2020s. With lower quantum noise and thermal coating noise, A+ is expected to be a 50% improvement over advanced LIGO design sensitivity. In the most optimistic scenario, all five detectors, including the three LIGO detectors in Livingston, Hanford and the planned detector in Aundha (India), the Virgo detector in Italy, and the KAGRA detector in Japan, will have sensitivities similar to A+ by the end of 2020s. Apart from the improved sensitivity, such a 5-detector network would also excel in accurately pinpointing the location of mergers in the sky [108, 110] and play an instrumental role in dark and bright siren cosmology [7, 44, 115]. Beyond A+, one of the proposed improvements to the three LIGO detectors is the A$^\sharp$ proposal [6] for the early 2030s. A$^\sharp$ is expected to be about twice as sensitive as A+ across the frequency band. The improved sensitivity has a direct impact on the number of sources detected as well as the estimation of parameters, leading to more accurate distance and localization estimation, which will result in stronger constraints on cosmological parameters [7].

Even after the improvement in sensitivity, the science that can be explored with the A+ and the A$^\sharp$ networks is limited by the size of the detectors. This is remedied by the proposed CE and ET observatories. The CE project [8, 181–183] involves the construction of two L-shaped detectors with arms of length 40 km for one detector and 20 km for the other. The order of magnitude increase in size and noise mitigation strategies result in $\mathcal{O}(10)-\mathcal{O}(100)$ improvement in the sensitivity compared to A+. On the other hand, ET [184, 185] is the proposed underground observatory in Europe, currently planned with three detectors with arms of length 10 km placed along the vertices of an equilateral triangle in a xylophone design (a comparison of different design proposals can be found in [111], including the possibility of two underground L-shape detectors at different locations). ET is expected to have better sensitivity than CE at frequencies less than 10 Hz. Due to this, the inspiralling binaries remain in the sensitive frequency band of the observatory for much longer, which leads to an improvement in parameter estimation. Both the CE and the ET observatories are projected to begin observations in the late 2030s.

The XG network with ET and the 2 CE detectors will be able to accomplish several science goals within just a few years of observation [7, 111]. Not only will the majority of mergers occurring in the Universe be detected, the measurement of source properties and localization will be performed at exquisite precision. As noted in section 2, this will enable GW observations to establish extremely stringent bounds on the measurement of the cosmological parameters, propelling our understanding of the Universe to unprecedented levels.

To get a sense of how these planned improvements translate into numbers, in figure 3 we project the expected number of detections as a function of year, taking into account our current understanding of the astrophysical population of compact binaries [137] and the current expectation for each detector to be operational. We split the count in BBHs and BNSs. From this graph it is clear that BBHs will continue to dominate the detection rates in the near future, while BNSs could potentially dominate in the XG era. But, as emphasized before, it is not only that more events will be detected, each upgrade in sensitivity will open the possibility of exploring an uncharted redshift range. To exemplify this, we also present in figure 3 the horizon redshift, i.e. the furthest a given GW detector can observe a signal, as a function of the total mass of the binary. As anticipated, XG detectors will reach all the way to redshifts of a hundred for BBHs and beyond the peak of star formation ($z\sim2$) for BNSs ⁵ .

Zoom In Zoom Out Reset image size

Apart from the ground-based observatories discussed above, the Laser Interferometer Space Antenna (LISA) [187], a space-based GW observatory, is planned to launch in the mid-2030s. This will consist of three satellites 2.5 million kms apart, located at the vertices of an equilateral triangle, moving in a heliocentric orbit. This will allow LISA to observe GW in the mHz regime. While LISA will be able to constrain cosmological parameters to unprecedented precision by itself, the GW sources utilized for this science goal are distinct from the ones used by ground-based detectors [11]. Instead of binaries containing neutron stars, massive BBH mergers are expected to serve as bright sirens. This is because the matter accreting on the massive black holes is expected to generate EM radiation close to merger [188, 189]. Despite the low expected detection rate, the exquisite sky localization [190] and the emission of EM radiation across all wavelengths [191] render these sources as promising bright siren candidates. For dark siren measurements, stellar-mass BBHs within the redshift of z = 0.1 can be employed as dark sirens, constraining H₀ to $\mathcal{O}(1)\%$ [192]. Extreme-mass ratio inspirals at higher redshifs $(z\sim0.5)$ can also be used as dark sirens and can constrain H₀ to $\mathcal{O}(1)\%$ as well [193]. Moreover, as LISA can detect massive BBH sources to high redshifts $(z\gt1)$, it can utilize the effects of weak lensing to infer clustering parameters [194]. These sources can also be used to test ΛCDM cosmology and probe the nature of dark energy [11]. LISA is also well positioned to test the validity of general relativity at cosmological scales [195].

LISA will enable a direct interplay with XG detectors. Stellar-mass binaries could be seen by LISA during their early inspiral before merging in band of the XG network [196]. This is particularly interesting for a possible population of binaries in the far side of the PISN mass gap [153], since the high-frequency sensitivity of LISA limits the detectability of binaries below the mass gap [197]. These 'multiband' detections will be good targets for cosmology as their first identification with LISA may allow for pre-merger localization, enhancing the possibility of finding multi-messenger counterparts. Even in the absence of counterpart the enlarged frequency evolution allow for a more precise localization in general that is beneficial for dark sirens. With just ∼10 detected events, H₀ could be constrained to a few percent [198].

Corresponding EM facilities will be critical for the GW siren measurements in the XG era. The observations of kilonova AT 2019gfo enabled the association to the host galaxy of the BNS GW170817, allowing for the first bright siren measurement [20, 34]. Dark siren measurements with compact binary mergers observed by LVK also relied on the availability and completeness of galaxy catalogue [101, 104, 133].

EM facilities will be necessary to search for EM counterparts of GW events for bright sirens. Due to the large localization area by GW detectors ($\mathcal{O}(1\, \textrm{deg}^2)$ for BNS mergers at 500 Mpc for A#, and $\lt2.5$ Gpc for XG [7]) and the fast-fading nature of the EM counterparts (within weeks depending on the wavelength and emission mechanism), the search for EM counterparts will likely require instruments with the capability of wide-field coverage and prompt response. Given the large expected number of BNS and NSBH merger events (∼a few hundred a day), dedicated GW event follow-up instruments might be needed. In figure 4, we present some of the planned and proposed wide-field telescopes across different wavelengths that could potentially contribute to the search for EM counterparts in the XG era ⁶ .

Zoom In Zoom Out Reset image size

In addition to the search for EM counterparts, redshift measurement for the hosts of bright sirens and golden dark sirens will also require adequate EM resources. Especially, if high redshift EM counterparts are found, such as high-redshift GRBs, it is less likely to have existing galaxy catalog coverage and the follow-up observations will be necessary [208].

Finally, deeper and more complete galaxy catalogs, such as those made available by Euclid [209] and DESI [94], will play an important role in future dark siren measurements. If nearby bright sirens or golden dark sirens were found with peculiar velocity as part of the concern for systematic uncertainty, EM follow-up observations will also be useful for the reconstruction of the velocity field around the events if a complete galaxy catalog is not already available.

The study of the cosmos has seen a revolution in the last three decades with the discovery of the accelerated expansion of the Universe [16, 17] and the precise mapping of the anisotropies of the cosmic microwave background (CMB) [26, 210]. These new cosmological probes have established a surprising 'standard' picture, called the ΛCDM model, in which most of the content of the present Universe is in the form of an unknown dark energy and dark matter. Much of the efforts of the large facilities that have been built and are planned to be built is precisely the unveiling of the physical origin of the dark Universe.

The main cosmological probes are the observation of the primordial Universe with the CMB [26, 210], the measurement of the expansion rate with supernova type Ia (SNIa) [16, 17] and baryon acoustic oscillations (BAO) [211], and the study of the growth of structure with galaxy surveys mapping the large scale structure (LSS) [212], although there are many other emergent probes [213]. All these probes aim at stress-testing ΛCDM in order to consolidate its position or discover new physics. Several tensions have been claimed over the years, but undoubtedly the most advertised one is the 'Hubble tension' [214, 215], reporting a disagreement between the measurement of the local expansion rate and the inferred value of H₀ from the CMB. It is unclear how long this tension will last, but, as we will see later, GW standard siren cosmology will provide a new way of precisely measuring the Universe's expansion rate.

There are currently multiple facilities with the potential to shed light on our understanding of dark energy and dark matter. Most notably, the DESI [94] and the Euclid satellite [209] will be able to map millions of galaxies over 10 billion years of cosmic history ($z\sim 0$–2) to constrain the properties of dark energy and the growth rate of the LSS using BAO, weak lensing maps and redshift space distortions among other methods. For example, DESI is expected to constrain the expansion rate H(z) to a few percent precision over this wide range of redshifts [94]. These surveys build up on the efforts of their predecessors, the DES [212] and the sloan digital sky survey (SDSS) [211]. In the near future, the Vera Rubin observatory (VRO) [216] will change the scale in which we map the Universe and observe its transient signals, playing a capital role in the multi-messenger transient endeavor. For instance, VRO is expected to detect over 10 million supernovae during a 10-year period, which can be compared with the order 1500 that DES accumulated in 5 years [217]. The Nancy Grace Roman Space Telescope [204] is expected to have a similar role from space, providing a detailed wide-field view of the sky. Roman seeks to measure both expansion history and structure growth parameters with 0.1%–0.5% precision level, reaching 0.9% and 2.1% constraints on H(z) between z = 1–2 and z = 2–3 respectively [204]. More into the future, the square kilometer array (SKA) [218] will look further back on time to explore the Universe when only the first galaxies were starting to be formed. There are also proposals for XG spectroscopic surveys to follow DESI aiming at collecting more than an order of magnitude more galaxy redshifts [219].

Following the methodology presented in [115, 116], golden dark sirens are identified as mergers that lie within z = 0.1 are can be localized in the sky to better than $0.04\,\mathrm{deg}^2$. Within z = 0.1, only one $L^{*}$ galaxy [220] is expected to be present in that small a sky-patch ([115, 116] use equation (7) in [29] to arrive at this limit). As this would allow unique identification of the host-galaxy, we assume that the redshift is known for such systems. Then, the distance errors obtained from Fisher analysis using gwbench [221] can be converted to constraints on H₀. Under the assumption of Flat ΛCDM cosmology, luminosity distance can be written a function of H₀ and the matter density, i.e. $D_\textrm{L} = D_\textrm{L}(\vec{\theta})$, where $\vec{\theta} = (H_0,\Omega_m)$. Then, the Fisher matrix obtained by combining the estimates from N events is given as

$$\begin{equation} \Gamma_{ij} = \sum_{k = 1}^{N} \frac{1}{\sigma_{D_\textrm{L}}^2} \left(\frac{\partial D_{\textrm{L}}}{\partial \theta_i}\right)\left(\frac{\partial D_{\textrm{L}}}{\partial \theta_j} \right)_k, \end{equation} \tag{ 4 }$$

where $\sigma_{D_\textrm{L}}^{2}$ is the error in luminosity distance from the GW detections obtained using gwbench. The square root of the diagonal elements of the covariance matrix $(\Gamma^{-1})$ associated with this Fisher matrix gives the constraints on H₀ and $\Omega_m$.

In this current form, the Fisher matrix denoted in equation (4) considers no prior information about the cosmological parameters. However, as $\Omega_m$ is not expected to be constrained by these nearby observations, the error on $\Omega_m$ estimated by the Fisher matrix will not only be unphysical (as $\Omega_m \in [0,1])$, but they will also adversely affect the errors in H₀. To mitigate this effect, we apply Gaussian priors to both H₀ and $\Omega_m$ with $(\sigma_{H_0},\sigma_{\Omega_m}) = (10,0.5)$. The prior on H₀ corresponds to the current bounds on H₀ from GW observations [107]. These can be included in the Fisher matrix calculation as [222],

$$\begin{equation} \Gamma_{ii} = \sum_{k = 1}^{N} \frac{1}{\sigma_{D_\textrm{L}}^2} \left(\frac{\partial D_{\textrm{L}}}{\partial \theta_i}\right)^2_k\, + \frac{1}{\sigma_{\theta_i}^2}\,. \end{equation} \tag{ 5 }$$

Based on the currently estimated redshift distribution of events, only ∼10 BBH mergers, ∼20 NSBH mergers and ∼120 BNS mergers are expected to occur within z = 0.1 every year. To avoid statistical fluctuations due to cosmic variance on our H₀ measurement estimates by our choice of the parameters of these handfuls of events, we generate 1000 permutations of the Universe, in each of which the parameters for these events are randomly selected from our chosen population models. The fractional error in H₀ is estimated for each realization of the Universe. The median values of these estimates along with the 99%-confidence regions are plotted in figure 5.

Zoom In Zoom Out Reset image size

Note that not all realizations of the Universe will have golden dark sirens. Thus, whenever a realization has no dark siren, we simply get back the prior that was put on H₀. From figure 5, we see that the number of golden dark sirens is much less compared to the bright sirens, but the constraints on H₀ are comparable for most cases. For the A$^\sharp$ network, BBH golden dark sirens can constrain H₀ to 3%. Most of the realizations of the Universe do not contain any NSBH or BNS golden dark sirens in the A$^\sharp$ detections, leading to the median value coinciding with the H₀ prior. This is also the case with a network that contains just the ET observatory. However, here, even the BBH golden dark sirens are not expected to contribute significantly to H₀. This is due to the stringent sky-localization constraint on the selection of golden dark sirens, where a network of three A$^\sharp$ detectors located far away from each other is seen to perform better than the three co-located detectors that make up the ET observatory. A similar observation is made for CE4020, where the two-detector network, with both detectors located in the US, is not able to measure the sky area well enough to qualify the detections as golden dark sirens. The CE40ET network does a much better job at constraining H₀, doing so to sub-percent precision for all the three classes of golden dark sirens. The constraints improve by a factor of 2–3 with the inclusion of CE20 in this network. Remarkably, for CE4020ET, all the realizations contain at least one golden dark siren, and almost all of them can constrain H₀ to better than 2%, resolving the H₀ tension in a year of observation.

As seen in section 2.1, BNS mergers can lead to a variety of EM counterparts which can be detected by various EM observatories (see section 3.2). In this section, we focus on kilonovae that are detected by LSST for bright siren measurement. From the simulated population of detected BNS events, we choose the sub-population that lies within the redshift of 0.5 and can be localized in the sky to $\Delta \Omega \unicode{x2A7D} 100\,\mathrm{deg}^2$, corresponding to each detector network. The choice of redshift corresponds to our expectation of the range of LSST for kilonova measurement, while the ΔΩ cut-off ensures that LSST can observe the kilonova within 10 sky patches. Then, using the kilonova model discussed in [223–225], we generate kilonova light curves in the ugrizy bands of LSST, for all the events in these sub-populations. Assuming single-exposure with an exposure time of 30 s, we claim that a kilonova is detected if its peak brightness exceeds the limiting magnitude of the telescope in that particular band. Using this approach, we find that most kilonova observations correspond to the $g-$band, with the maximum redshift at which a kilonova is detected to be z = 0.28 (validating our $z\unicode{x2A7D}0.5$ cut-off for the sub-populations). Further, we assume that only 40% of these kilonovae will actually be detected [208]. This fraction aims to account for the LSST sky coverage, systematic biases as well as duty-cycles corresponding to the GW and EM observatories. The remaining set of events, those that are detected both in the GW and the EM band, are used for H₀ measurement using the same methodology that was used for golden dark sirens in section 4.1.

The results for bright siren bounds on H₀ are portrayed in figure 5. We see that our bright siren estimates for H₀ are either at par with, or better than, the golden dark siren measurements. This disparity between the two approaches becomes notably significant for the HLA and the CE4020 networks. This is because HLA is not sensitive enough to resolve dark sirens to a very small area in the sky and uniquely identify the host galaxy. CE4020 is unable to achieve the same as it is only a two-detector network. The limitation of CE4020 to localize events in the sky is also evident in the bright siren estimates, where CE4020 detects 131 bright sirens, whereas the three-detector network HLA detects 196 such events. However, as the events detected by CE4020 have high SNRs and, accordingly, better distance estimates, the H₀ constraints with CE4020 are better than those from HLA.

We note that our bright siren estimates could be on the conservative side, as a result of our choice of the neutron star EOS, LSST detection strategy, or restricting ourselves to kilonova follow-up. In particular, we chose APR4 [226] to be the neutron star EOS, which results in relatively more compact neutron stars. Due to higher compactness, these neutron stars are difficult to disrupt tidally, reducing the amount of ejecta produced, which lowers the number of KN detections [116]). Further, several studies have underscored the significance of short gamma-ray burst detections in enhancing the precision of bright siren measurements [167, 208, 227, 228]. These studies particularly emphasize their efficacy in probing events at high redshifts, thereby contributing to constraints on $\Omega_m$. Moreover, while our focus has been on kilonovae observations BNS mergers, a valuable avenue lies in considering NSBH mergers, which can also produce kilonova counterparts, thereby extending the pool of potential bright sirens for H₀ measurement [116, 129]. Nevertheless, amidst these considerations, the overarching conclusion is that BNS bright sirens can, by themselves, achieve H₀ measurements with sub-percent accuracy. This underscores the potential of BNS bright sirens as a powerful tool for precise cosmological parameter estimation.

XG GW detectors will be revolutionary both in terms of the number of detections and redshift range. Compact binary coalescences will be detected every few minutes (approximately every 10 min for BBHs and every 3 min for BNSs), adding up to tens of thousands to hundreds of thousands of events per year. Of them, a large fraction will be beyond $z\sim1$. This makes XG CBC populations a perfect target for spectral siren cosmology, opening up the high-z Universe.

We forecast the capabilities of spectral sirens performing a hierarchical Bayesian analysis on simulated data. We follow our fiducial population of BBHs and BNSs described in the appendices and simulate their detected posterior samples following the method described in [149]. We assume a flat ΛCDM cosmology whose expansion rate is described by H₀ and $\Omega_m$, cf equations (1) and (2), taking as fiducial values the ones reported in Planck 2018 [26]: $H_0 = 67.66$ km s⁻¹ Mpc⁻¹ and $\Omega_m = 0.31$. We take a 40 km CE sensitivity as representative for XG detectors. The inference is performed using numpyro [229] and accelerated with JAX [230]. Our code is available on github. The inference is performed across all the cosmological and mass function parameters. We do not consider NSBHs as their mass spectrum is still largely unknown.

In figure 6 we present the projected constraints on the expansion rate of the Universe as a function of redshift. We compare the 1σ relative errors for a fixed number of BBHs and BNSs, taking the error with respect to the mean of the H(z) traces at each redshift. Given current rate uncertainties, we choose a fixed number of 50 000 events which corresponds approximately to the detection rate of BBHs per year with XG detectors (see figure 3). Since we assumed the same merger rate history $\mathcal{R}(z)$ for all compact binaries (see details in appendix B.2) and BNS selection effects brings them to lower redshifts, the population of BNSs has more constraining power at lower z. Within the population of BBHs, it is the peak at low masses that dominates the inference [149]. Both BNSs and BBHs can constrain the expansion rate with sub-percent precision over a wide range of redshifts, with the best constrained values at around $z\sim1$. The best constrained redshift is located where most of the (well constrained) events are detected, which depends on the merger rate history and the detector's sensitivity. A combined analysis of all CBCs will only improve these constraints and allow to more efficiently break the possible degeneracies with the redshift evolution of the mass spectrum [149].

Zoom In Zoom Out Reset image size

Although we have introduced each type of standard siren separately for pedagogical reasons, they will all be analyzed simultaneously with XG detectors. In fact, as already discussed in the context of current detectors [133], the population inference of the spectral siren method must be intertwined with the dark siren method, since an incorrect population model assumption can bias the dark siren cosmology. Recently, there have been developments of algorithms to jointly infer the GW population and cosmology with additional galaxy catalog information using hierarchical Bayesian analyses [231, 232] and neural posterior estimation [233].

To exemplify a joint forecast of different standard sirens, we focus on a population of BBHs and exploit the capabilities of golden dark sirens to narrowly constrain the local expansion rate and on the exploration of the high-z expansion history with spectral sirens. In practice, we take the projected constraints on H₀ from the Fisher analysis of golden dark sirens described in section 4.1 and summarized in figure 5, and use them as the uncertainty for a Gaussian prior centered around the fiducial H₀ that is applied on the hierarchical Bayesian population analysis. As a benchmark point, we take an error in H₀ of 1% to represent the capabilities of BBH dark sirens with a XG network.

In figure 7 we present the results for the joint forecast. We compare the combined analysis of spectral sirens and (golden) dark sirens with only spectral sirens for a population of BBHs. We show both the posterior distributions in H₀ and $\Omega_m$ (left of figure 7) and the relative error in the expansion rare H(z) (right of figure 7). As expected, the golden dark sirens dominate the local expansion rate while the spectral sirens are more powerful at z > 1. Together, the expansion rate is constrained with sub-percent precision at all redshifts.

Zoom In Zoom Out Reset image size

Although we have focused so far on the capabilities of future GW detectors to constrain the cosmological parameters describing the ΛCDM model, namely H₀ and $\Omega_m$ (see equation (2)), GW standard sirens are equally well suited to test deviations from this standard model. These generically come in two flavors: (i) as modifications in the properties of the different components of ΛCDM or (ii) as breaks in the foundational principles of the model.

For the first class, a lot of emphasis has been placed on postulating that the current accelerated expansion of the Universe is driven dynamically, for example, by a new scalar field. In practice this means that the energy density of dark energy is no longer constant, but rather evolves with time. This would change $\Omega_\Lambda$ in equation (2) to $\Omega_\mathrm{de}(z)$. Since this dynamical evolution can in principle be quite arbitrary, phenomenological approaches in which the EOS of DE is parametrized in simple terms have been common for encapsulating all the new physics, i.e. $\Omega_\mathrm{de}(z) = \Omega_\Lambda(1+z)^{3(1+w_\mathrm{de}(z))}$. Early work already anticipated that the EOS of DE was a clear target for bright siren analyses with future GW detectors [234]. Similarly, $w_\mathrm{de}$ was studied in the context of spectral sirens with the BNS mass function [235].

For the second class, many (exotic) variants are possible. One that can particularly be well tested with GWs is that general relativity (GR) does not hold at cosmological scales. Among other things, this changes the propagation of GWs. Since compact binary coalescence occur at cosmological distances, even very small deviations in the GW propagation with respect to GR can be strongly constrained. The paradigmatic modifications of gravity at cosmological scales that have been explored are those that change the speed of propagation and the scaling of the luminosity distance with redshift [236]. In particular, for the case of scalar-tensor theories, the tremendous power of tests of the speed of gravity were anticipated [237, 238] and spectacularly demonstrated with the first multi-messenger event, GW170817 [239–242]. After the tight constraints on the speed of gravity, which cannot differ from the speed of light in more than 1 part in 10¹⁵, testing whether the GW luminosity distance is equal to the one of EM radiation has become the next goal. This modified luminosity distance can be parametrized in terms of an additional friction term in the propagation equation ν, where GR predicts ν = 0. Current bright siren observations can only give weak constraints given the closeness of the sources [243], $|\nu|\lt \mathcal{O}(30)$, but spectral sirens can improve multi-messenger bounds given its larger redshift range [244], $|\nu|\lt \mathcal{O}(3)$. There has been a lot of activity to constrain this modified propagation with different standard siren methods [131, 245–251] and to forecast their future capabilities with XG detectors [9, 246, 252, 253], which could improve by another order of magnitude. Besides these effects, future GW detectors will be able to constrain other phenomena. For example, they will bound the presence of additional cosmological tensor fields [254, 255], which have a rich phenomenology of waveform distortions [256], and test lensing effects beyond GR [257, 258]. Future GW observations combined with the redshifts from their host galaxies will also serve to measure peculiar velocities in order to constrain the growth of structure and test gravity [259].

The population of BBHs is consistent with the inferred population with GWTC-3 events [180]. The population parameters for BNS and NSBH events are uncertain, owing to the lack of GW observations. The NS EOS is chosen to be APR4 [226]. Following the maximum mass limit applied by APR4, we choose a uniform mass distribution between $[1,2.2]M_\odot$ for neutron stars. The fiducial mass spectrum of BNSs and BBHs is presented in figure B1.

Zoom In Zoom Out Reset image size

We assume the population to have spins aligned with the orbital angular momentum. As precession, in general, is expected to improve the estimation of parameters [265], the measurability estimates presented in this work may be conservative. The specifications for the three populations are given below:

• (i)  
  BBHs

  • Primary mass: POWER LAW + PEAK [180] model with the following values for the model parameters: $\alpha = -3.4$, $m_\textrm{min} = 5M_\odot$, $m_\textrm{max} = 87M_\odot$, λ = 0.04, $\mu_\textrm{peak} = 34M_\odot$, $\sigma_\textrm{peak} = 3.6$, $\delta_m = 4.8M_\odot.$
  • Mass ratio: $p(q) \propto q^\beta$ with β = 1.1, and enforcing $m_\textrm{min} = 5M_\odot.$
  • Spin magnitude: Aligned, independently and identically distributed (IID) spins following a beta distribution with $\alpha_\chi = 2$, $\beta_\chi = 5$ (see equation (10) in [266]).
  • Waveform: IMRPhenomXHM [267]
• (ii)  
  BNSs

  • Mass: Uniform between $[1,2.2]M_\odot$.
  • Spin magnitude: Aligned and uniform between $[-0.05,0.05]$.
  • EOS: APR4
  • Waveform: IMRPhenomPv2_NRTidalv2 [159]
• (iii)  
  NSBHs

  • Black hole mass: Same as the primary mass for BBHs.
  • Neutron star mass: Same as mass for BNSs.
  • Black hole spin magnitude: Gaussian with µ = 0 and σ = 0.2.
  • Neutron star spin magnitude: Same as for BNSs.
  • Waveform: IMRPhenomXHM

For the three populations, we choose the merger rate to follow the Madau-Dickinson star formation rate [268, 269],

$$\begin{equation} \psi\left(z|\gamma, \kappa, z_p\right) = \frac{\left(1 + z\right)^\gamma}{1 + \left(\frac{1 + z}{1 + z_p}\right)^\kappa}, \end{equation} \tag{ B.1 }$$

with γ = 2.7, $z_{p} = 1.9$, and κ = 5.6. While the most accurate approach would be to convolve a chosen time-delay distribution with the star formation rate to obtain the merger rate, we have decided against choosing a time-delay distribution for this work. This is because the time-delay distribution is still uncertain, and the Madau-Dickinson star formation rate by itself is consistent with the current bounds on the merger rate with GWTC-3 events [180]. Further, the redshift distribution is obtained by normalizing the star formation rate to match the local merger rate densities for the compact binaries. Following the currently inferred rates [180, 270], the local merger rate densities for BBH, NSBH and BNS systems are chosen to be $\mathcal{R}_0^{\mathrm{BBH}} = 24\,\,\mathrm{Gpc}^{-3}\,\mathrm{yr}^{-1}$, $\mathcal{R}_0^{\mathrm{NSBH}} = 45\,\,\mathrm{Gpc}^{-3}\,\mathrm{yr}^{-1}$ and $\mathcal{R}_0^{\mathrm{BNS}} = 320\,\,\mathrm{Gpc}^{-3}\,\mathrm{yr}^{-1}$, respectively. The resultant redshift distribution is plotted in figure B2.

Zoom In Zoom Out Reset image size