Cosmography with next-generation gravitational wave detectors

Advancements in cosmology through next-generation (XG) ground-based gravitational wave (GW) observatories will bring in a paradigm shift. We explore the pivotal role that GW standard sirens will play in inferring cosmological parameters with XG observatories, not only achieving exquisite precision but also opening up unprecedented redshifts. We examine the merits and the systematic biases involved in GW standard sirens utilizing binary black holes, binary neutron stars, and neutron star-black hole mergers. Further, we estimate the precision of bright sirens, golden dark sirens, and spectral sirens for these binary coalescences and compare the abilities of various XG observatories (A ♯ , cosmic explorer, Einstein telescope, and their possible networks). When combining different sirens, we find sub-percent precision over more than 10 billion years of cosmic evolution for the Hubble expansion rate H(z). This work presents a broad view of opportunities to precisely measure the cosmic expansion rate, decipher the elusive dark energy and dark matter, and potentially discover new physics in the uncharted Universe with XG GW detectors.


Introduction
In the next decades, gravitational-wave (GW) observations are expected to uncover many mysteries in the Universe.The planned next-generation ground-based and spacebased GW observatories will bring a wealth of science opportunities with their reaches to the early Universe.The current generation of ground-based GW observatories, LIGO [1], Virgo [2], and KAGRA [3], with their successive upgrades (e.g., A+ [4,5] and A ♯ [6]) will be able to detect stellar-mass binary black hole (BBH) mergers up to redshift z ∼ 10 [7].Next-generation (XG) ground-based GW detectors, Cosmic Explorer (CE) [8] and Einstein Telescope (ET) [9], aim at observing binary neutron star (BNS) and neutron star-black hole (NSBH) mergers up to z ∼ 10 as well as all stellar-mass BBHs in the Universe (z ∼ 100).‡ These observations of mergers across cosmic time will make XG detectors new powerful cosmological probes, mapping the Universe from the present dark energy domination to the past dark matter era, cf.figure 1. Future GW detectors in space will also bring relevant cosmological constraints [11][12][13][14], but we focus on ground-based detectors only here.
In this era of precision cosmology with many well-established cosmological measurements, GW observations have several unique roles to play.First, GW observations will independently measure the cosmological parameters, with the potential to resolve tensions among existing electromagnetic (EM) measurements [15].Second, the detections across different redshifts allow for the mapping of the Universe's expansion history with a single probe over more than 10 billion years of cosmic evolution.Third, XG detectors will verify the consistency between the EM and GW cosmological measurements and reveal (if any) underlying new physics.These observations will test gravity in new regimes, reaching unprecedented strong fields and cosmic scales.
Assuming general relativity, the distances to a given binary can be determined from the GW signals released from the binary; meanwhile, multiple techniques are available to estimate the redshift of the binary.With the distance and redshift estimate, GW sources can be used to measure the Universe expansion rate H(z) and construct the GW Hubble diagram.In analogy to the use of supernovae Type Ia as standard candles for cosmology [16,17], this is known as the 'standard siren' method [18,19].The observations of the BNS GW170817 and its EM counterparts enabled the first standard siren measurement [20], opening a new era in multi-messenger astronomy (MMA).After the method was proposed in 1986 [18], new ideas on redshift acquirement, developments in data analysis, applications to GW data, and investigations of systematic uncertainties promised an exciting path forward for standard siren measurements.In this article, we explore further the prospects of the standard siren cosmography in the XG era.
In the following, we first review different techniques to obtain redshifts for standard sirens and the cosmology landscape with GW and EM facilities in 2035+.We then provide future projections to discuss the prospects of standard siren cosmography with XG detectors.For concreteness, we focus on the forecast of the expansion rate of the Universe and its implications for measurements of the parameters of the standard cosmological models and constraints on its possible extensions.The science of future GW detector however extends further and includes, for example, cross-correlation with cosmic surveys [21,22], gravitational lensing [23] and stochastic backgrounds [24].‡ NEMO is another possible future ground-based detector with sensitivity similar to A+ but ranging all the way to ∼ 4kHz [10].Therefore, NEMO is optimized to detect the (post-)merger of BNS.We do not discuss NEMO in this work. .Schematic diagram of the potential of GW observations to probe the evolution of the Universe.The radial coordinate represents the time from the present (center) to the beginning of the Universe (edge).On the right side of the diagram, distribution of stellar-mass BBH mergers (blue dots) and their associated horizon distance for different GW detectors (dashed lines).On the left side of the diagram, fraction of each Universe's component (dark energy, dark matter, ordinary matter and radiation) as a function of time, where a fraction of 1 at a given time is represented by the full arc.Although at present time we are observing a Universe dominated by dark energy, future GW detectors will observe an ancient Universe dominated by dark matter.

Standard siren cosmology
Gravitational waves are natural messengers to probe the cosmic history.Unlike most other transients in the Universe, their signals can be well-understood from first principles, as general relativity has clear predictions for the waveform of a compact binary coalescence.For example, at leading order, the frequency evolution is determined by a specific combination of the detected component masses of the binary known as the chirp mass M c = (m 1 m 2 ) 3/5 /(m 1 + m 2 ) 1/5 : ḟgw ∝ M 5/3 c f 11/3 gw [25].Moreover, GW detectors are directly sensitive to the waveform of the signal (as opposed to the intensity) allowing us to observe signals further away, since the GW amplitude scales with the inverse of the luminosity distance (and not the inverse square of the distance).Since GWs are coherently detected, the absolute value of the luminosity distance is directly measured-giving GW sources the name of standard sirens.The GW strain, h gw , thus contains information about the expansion rate through the luminosity distance where we have assumed a flat background cosmology.Under that assumption and focusing on a Universe dominated by dark energy and dark matter (radiation is only relevant in the early Universe), the Hubble parameter measuring the expansion rate can be written as where H 0 is the Hubble constant and denotes the local expansion rate H 0 = H(z = 0), and Ω Λ and Ω m are the fractional energy density of dark energy and dark matter that, in a flat Universe, satisfy Ω Λ + Ω m = 1.The latest cosmological observations from the Planck satellite result in H 0 = 67.66 ± 0.42km/s/Mpc, Ω Λ = 0.6889 ± 0.0056 and Ω m = 0.3111±0.0056[26].Alternative cosmological models are also possible, potentially modifying both the background cosmology H(z) and the scaling of the strain with the luminosity distance.We will discuss those in section 4.5.
GWs are unique compared to other transients because they keep an almost pristine record of the time of the merger due to the negligible interaction with the medium during propagation.In addition, they are particularly well suited for statistical population analyses since GW detector networks monitor the whole sky continuously and their selection effects can be understood in detail with injection studies [27].
What GWs do not provide, however, is a direct account of the precise epoch when the binary merged.Cosmic times, or redshifts, together with distances allow the measurement of the Universe's expansion rate at different epochs.Therefore, the entire game of GW cosmology is to obtain redshift information.In the following we review different methods to achieve this goal.

Bright sirens
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][29][30][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.
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 [75][76][77], the EM observation selection effect [78][79][80], the biased EM-inferred binary viewing angle [78], the peculiar velocity of the hosts [81][82][83], and the GW instrumental non-stationary noise [84][85][86][87][88].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.

Dark sirens
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 [89][90][91][92].This cross-correlation technique was applied to eight well-localized BBH events from the Gravitational-Wave Transients Catalog-3 (GWTC-3) [93] to get H 0 = 68.2+26.0 −6.2 km s −1 Mpc −1 [94].LVK detectors at design sensitivities, together with DESI [95] and SPHEREx [96] data, have been shown to constrain H 0 to O(1)% with 5 years of observation [97].
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 [98]) 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,[99][100][101].
GW170814 [102] was the first BBH merger to which this technique was applied.It was localized to a volume of 2 × 10 6 Mpc 3 [103] which, according to the Dark Energy Survey (DES) data, contained ∼ 77, 000 galaxies [104].In comparison, the localization volume of the best localized BBH, GW190814 [105], was ∼ 100 times smaller and contained ∼ 2000 galaxies [106].Using GW170814 and the DES data [107] for redshift information, H 0 was found to be 75 +40 −32 km s −1 Mpc −1 [104].When applied to the eight best-localized events in GWTC-3, stronger constraints are obtained on H 0 , measured to be 79.8 +19.1 −12.8 km s −1 Mpc −1 [103].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.Ref. [108] applies the statistical siren strategy to GW170817 to obtain H 0 = 76 +48 −23 km s −1 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 0 .The constraints obtained from the statistical dark siren and the bright siren approaches were combined for events in GWTC-3 to estimate H 0 with ∼ 15% uncertainty [103,109].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,[110][111][112][113].This will result in smaller localization volumes containing smaller number of potential hosts compared to currently detected events [114].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 [115,116].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 (ι) is not zero) systems are touted to be ideal candidates for such sirens [117,118].This is due to the excitation of higher-order modes in the GW waveform for asymmetric and face-off systems [119][120][121], leading to the breaking of the degeneracy between d L − ι, which leads to the precise estimation of parameters [122][123][124][125][126][127].Similar expectations hold for precessing systems in which the spins are not (anti-)aligned and there could be precession of the orbital plane [127][128][129].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.Ref. [114] 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 0 (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, Ref. [117] claim that just by using BBH golden dark sirens, a network with A+ sensitivities can ascertain H 0 to O(1%) precision (with 68%−confidence) and an ET+2CE network can measure H 0 to O(0.1%) (with 68%−confidence) in two years of observation.Ref. [118] 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 0 from NSBH golden dark siren mergers can range from O(0.1%) − 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 [118,130,131] earns them the name grey sirens [118].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,100,132,133].Ref. [132] finds that H 0 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 L , it is significantly affected by the uncertainty introduced due to calibration errors, which can affect the amplitude measurement by as much as 7% [134].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 ∼ 0.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 ≤ 0.1 [117,118].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 [135].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 [136].

Spectral sirens
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 where m is the mass in the source frame.Since all binaries at the same luminosity .Cartoon of the spectral siren method.For increasing luminosity distance, the whole distribution of detector frame masses shifts to higher values.The amount shifted corresponds to the redshift at a given luminosity distance and it is therefore sensitive to the expansion rate.For example, a higher value of H 0 associates higher detected masses to the same distance.
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 0 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.
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 [137,138].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 [139,140] and an excess over a simple power-law at ∼ 35M ⊙ [141].At first sight, these features seem to coincide with the expectations from the theory of pair instability supernova (PISN) [142][143][144][145][146][147], which predicts a gap in the BBH mass spectrum between ∼ 50-120M ⊙ and a potential pile-up just before the lower edge of the gap [148].This connection has been however recently revisited [149], 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 [150] and have also been explored in the context of XG interferometers [151].With the latest LVK catalog, spectral siren constraints on H 0 are ∼ 20% at 1σ [136].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 [152].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 [152].If unaccounted for, however, the evolution of the mass spectrum can bias the inference of cosmological parameters [152][153][154], a potential problem especially relevant for XG detectors observing far away binaries.Mismodeling of the true distribution with oversimplified parametric models [154] 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 [155].
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 ∼ 35M ⊙ and the lack of them at > 45M ⊙ , spectral siren cosmology with XG detectors will be dominated by BNSs and low-mass BBHs, as their intrinsic merger rate is much larger [152].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 [156].

Love sirens
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 allneutron 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) [157,158].The leading order tidal contribution is a linear function of the individual tidal deformability parameters of the two neutron stars, Λ 1 and Λ 2 , and is called the reduced tidal deformability parameter ( Λ).If the neutron stars are spinning, the individual Λs will also appear in the contributions to the phase from spin-induced quadrupolar deformations [159][160][161][162].These changes in the phase allow the measurement of Λ 1 and Λ 2 when GWs are detected.Now, given an equation of state (EOS) for neutron stars, one can relate the individual Λs to sourceframe 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 [163].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 [164].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 [165]).
Note that, in practice, only Λ and not the individual Λ are expected to be precisely measured, even with XG observatories [166].However, there exist phenomenological relations, called quasi-universal relations [167], that are independent of the EOS governing the inspiralling neutron stars and allow the estimation of individual Λs using Λ and the mass-ratio of the binary.While this approach can indeed provide narrower bounds on the estimates of Λ 1 and Λ 2 , they come at the cost of incurring biases due to uncertainties in the quasi-universal relations [167,168].
Using the inspiral technique for a BNS system in ET sensitivity, the redshift can be measured to a precision of O(10%) for systems up to z ∼ 1 [163,169].With the postmerger spectra, similar precision is obtained up to z ∼ 0.04 [164].Ref. [170] 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 ∼ 1% for O(10 4 ) yearly detections.This translates to 0.1% error in H 0 and 0.61% error in Ω m with an ET+2CE network in a year of observation [170].Ref. [171] shows that the XG observatories can constrain the EOS with BNS sources around z ∼ 1.
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,[172][173][174].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 [175].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 Λ 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 [176][177][178] or by specifically correcting for the incurred errors [168].

The road to XG detectors
With more than one decade to go, the road to XG detectors will build up on many other efforts and will interplay with many other actors.On the GW side, several detector technologies still need to be demonstrated and their associated observing periods will continuously shed light on the cosmological origin of compact binary coalescences.On the EM side, standard siren cosmology will crucially depend on the capabilities of available telescopes to follow up the GW events.Moreover, cosmic surveys also aim at probing the cosmological model, acting as complementary partners to GW observatories.For this reason, before projecting the cosmology science case of XG detectors in section 4, we briefly summarize the expected GW, EM, and cosmology landscape by 2035+ to understand better the new horizons opened by XG detectors.

GW astronomy in 2035+
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 [179], inference of cosmological parameters [109], revelations about neutron star properties [32,[180][181][182][183][184] and measurement of the population properties of compact binary mergers [185].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 ∼ 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 [110,112] and play an instrumental role in dark and bright siren cosmology [7,44,117].Beyond A+, one of the proposed improvements to the three LIGO detectors is the A ♯ proposal [6] for the early 2030s.A ♯ 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 ♯ networks is limited by the size of the detectors.This is remedied by the proposed CE and ET observatories.The CE project [8,[186][187][188] 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 O(10) − O(100) improvement in the sensitivity compared to A+.On the other hand, ET [189,190] 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 Ref. [113], 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,113].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 [140] 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 ∼ 2) for BNSs.§ Apart from the ground-based observatories discussed above, the Laser Interferometer Space Antenna (LISA) [192], 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 [193,194].Despite the low expected detection rate, the exquisite sky localization [195] and the emission of EM radiation across all wavelengths [196] render these sources as promising bright siren candidates.For dark siren mea- § Note, however, that as the distance increases signals will be quieter in the detectors and, as a consequence, the inference of their properties will be degraded, including the distance itself [191].
surements, stellar-mass BBHs within the redshift of z = 0.1 can be employed as dark sirens, constraining H 0 to O(1)% [197].Extreme-mass ratio inspirals at higher redshifs (z ∼ 0.5) can also be used as dark sirens and can constrain H 0 to O(1)% as well [198].Moreover, as LISA can detect massive BBH sources to high redshifts (z > 1), it can utilize the effects of weak lensing to infer clustering parameters [199].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 [200].
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 [201].This is particularly interesting for a possible population of binaries in the far side of the PISN mass gap [156], since the high-frequency sensitivity of LISA limits the detectability of binaries below the mass gap [202].These "multiband" detections will be good targets for cosmology as their first identification with LISA may allow for premerger 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 0 could be constrained to a few percent [203].

EM observatories in 2035+
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 [103,106,136].
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 (O(1deg 2 ) for BNS mergers at 500 Mpc for A#, and < 2.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 widefield 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∥.
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 [204].
∥ The figure does not mean to capture the comprehensive list of missions or to evaluate between missions.Finally, deeper and more complete galaxy catalogs, such as those made available by Euclid [205] and DESI [95], 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.

Cosmology landscape in 2035+
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,215].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,215], the measurement of the expansion rate with supernova type Ia (SNIa) [16,17] and baryon acoustic oscillations (BAO) [216], and the study of the growth of structure with galaxy surveys mapping the large scale structure (LSS) [217], although there are many other emergent probes [218].All these probes aim at stresstesting Λ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" [219,220], reporting a disagreement between the measurement of the local expansion rate and the inferred value of H 0 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 Dark Energy Spectroscopic Instrument (DESI) [95] and the Euclid satellite [205] will be able to map millions of galaxies over 10 billion years of cosmic history (z ∼ 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 [95].These surveys build up on the efforts of their predecessors, the Dark Energy Survey (DES) [217] and the Sloan Digital Sky Survey (SDSS) [216].In the near future, the Vera Rubin Observatory (VRO) [221] will change the scale in which we map the Universe and observe its transient signals, playing a capital role in the multimessenger 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 [222].The Nancy Grace Roman Space Telescope [211] 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 [211].More into the future, the Square Kilometer Array (SKA) [223] 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 next-generation spectroscopic surveys to follow DESI aiming at collecting more than an order of magnitude more galaxy redshifts [224].

New horizons with XG detectors
In this section, we will make projections on the measurement of cosmological parameters using GW observations.From the same, we have considered various detector networks, starting from A ♯ sensitivity to a network with ET and two CE observatories.Details about the network configurations can be found in Appendix A. Specifically, we simulate populations of BBH, NSBH and BNS mergers consistent with current observations and following common astrophysical assumptions (details are described in Appendix B.1 and Appendix B.2).The golden dark siren approach has been applied to the three compact binary classes and is presented in section 4.1.In section 4.2, we generate kilonova light curves for BNS mergers and project the number of bright sirens based on the number of kilonovae that can be detected by LSST.The bounds on cosmological parameters using the spectral siren method with BBH and BNS mergers have been discussed in section 4.3.The joint forecast based on the estimates from these three sections has been discussed in section 4.4.

Golden dark sirens
Following the methodology presented in Refs.[117,118], 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 deg 2 .Within z = 0.1, only one L * galaxy [225] is expected to be present in that small a skypatch (Refs.[117,118] use equation (7) in Ref. [226] 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 [227] can be converted to constraints on H 0 .Under the assumption of Flat ΛCDM cosmology, luminosity distance can be written a function of H 0 and the matter density, i.e., D L = D L ( ⃗ θ), where ⃗ θ = (H 0 , Ω m ).Then, the Fisher matrix obtained by combining the estimates from N events is given as where σ 2 D L is the error in luminosity distance from the GW detections obtained using gwbench.The square root of the diagonal elements of the covariance matrix (Γ −1 ) associated with this Fisher matrix gives the constraints on H 0 and Ω m .
In this current form, the Fisher matrix denoted in equation 4 considers no prior information about the cosmological parameters.However, as Ω m is not expected to be constrained by these nearby observations, the error on Ω m estimated by the Fisher matrix will not only be unphysical (as Ω m ∈ [0, 1]), but they will also adversely affect the errors in H 0 .To mitigate this effect, we apply Gaussian priors to both H 0 and Ω m with (σ H 0 , σ Ωm ) = (10, 0.5).The prior on H 0 corresponds to the current bounds on H 0 from GW observations [109].These can be included in the Fisher matrix calculation as [228], 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 0 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 0 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. Figure 5. Forecasts on the bounds that can be placed on H 0 with golden dark sirens and bright sirens in the observation span of a year.The markers show the median value of the constraints from the 1000 realizations of the universe, and error bars mark the 99% confidence region.The number of events that contribute to the corresponding bounds is also shown beside the error bars.The dotted line shows the 2% precision mark in H 0 measurement, and the dashed line shows the current standard siren uncertainty on H 0 (∼ 15%), which has been used as a prior in the analysis.
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 0 .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 0 are comparable for most cases.For the A ♯ network, BBH golden dark sirens can constrain H 0 to 3%.Most of the realizations of the universe do not contain any NSBH or BNS golden dark sirens in the A ♯ detections, leading to the median value coinciding with the H 0 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 0 .This is due to the stringent sky-localization constraint on the selection of golden dark sirens, where a network of three A ♯ 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 0 , 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 0 to better than 2%, resolving the H 0 tension in a year of observation.

Bright sirens
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 ∆Ω ≤ 100 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 Refs.[229][230][231], 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 seconds, 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 ≤ 0.5 cut-off for the sub-populations).Further, we assume that only 40% of these kilonovae will actually be detected [232].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 0 measurement using the same methodology that was used for golden dark sirens in section 4.1.
The results for bright siren bounds on H 0 are portrayed in figure 5. We see that our bright siren estimates for H 0 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 0 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 [233] 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 [118]).Further, several studies have underscored the significance of short gamma-ray burst detections in enhancing the precision of bright siren measurements [170,232,234,235].These studies particularly emphasize their efficacy in probing events at high redshifts, thereby contributing to constraints on Ω 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 0 measurement [118,131].Nevertheless, amidst these considerations, the overarching conclusion is that BNS bright sirens can, by themselves, achieve H 0 measurements with sub-percent accuracy.This underscores the potential of BNS bright sirens as a powerful tool for precise cosmological parameter estimation.

High-z spectral sirens
Next-generation 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 minutes for BBHs and every 3 minutes 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 ∼ 1.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 [152].We assume a flat ΛCDM cosmology whose expansion rate is described by H 0 and Ω m , cf equation (1) and 2, taking as fiducial values the ones reported in Planck 2018 [26]: H 0 = 67.66km/s/Mpcand Ω m = 0.31.We take a 40km CE sensitivity as representative for XG detectors.The inference is performed using numpyro [236] and accelerated with JAX [237].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 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 [152].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 ∼ 1.The best We compare the 1σ relative errors in H(z) as a function of redshift for 50,000 BBHs and BNSs, using a 40km Cosmic Explorer detector as the fiducial sensitivity.
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 [152].

All-in XG cosmography
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 [136], 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 [238,239] and neural posterior estimation [240].
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 0 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 0 that is applied on the hierarchical Bayesian population analysis.As a benchmark point, we take an error in H 0 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 0 and Ω 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.

Stress-testing the standard cosmological model
Although we have focused so far on the capabilities of future GW detectors to constrain the cosmological parameters describing the ΛCDM model, namely H 0 and Ω 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 Ω Λ in equation (2) to Ω de (z).Since this dynamical evolution can in principle be quite arbitrary, phenomenological approaches in which the equation of state of DE is parametrized in simple terms have been common for encapsulating all the new physics, i.e.Ω de (z) = Ω Λ (1 + z) 3(1+w de (z)) .Early work already anticipated that the equation of state of DE was a clear target for bright siren analyses with future GW detectors [241].Similarly, w de was studied in the context of spectral sirens with the BNS mass function [242].
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 [243].In particular, for the case of scalar-tensor theories, the tremendous power of tests of the speed of gravity were anticipated [244,245] and spectacularly demonstrated with the first multi-messenger event, GW170817 [246][247][248][249].After the tight constraints on the speed of gravity, which cannot differ from the speed of light in more than 1 part in 10 15 , 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 [250], |ν| < O( 30), but spectral sirens can improve multi-messenger bounds given its larger redshift range [251], |ν| < O(3).There has been a lot of activity to constrain this modified propagation with different standard siren methods [133,[252][253][254][255][256][257][258] and to forecast their future capabilities with XG detectors [253,[259][260][261], 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 [262,263], which have a rich phenomenology of waveform distortions [264], and test lensing effects beyond GR [265,266].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 [267].

Conclusions
GW observations have opened a new window to explore the cosmos, unveiling a large population of BBHs and starting the era of multi-messenger astronomy with BNSs.XG GW detectors will enter the "big-data" phase and independently probe the Universe from low to high redshift.These advancements will be key in the quest for precision GW cosmology.Different standard siren methods will complement each other in the XG era, having capabilities that cannot be achieved with current detectors.On the one hand, very loud, nearby events will have such a good localization that their EM counterparts can be easily sought for or only one galaxy may lay within their sky map, granting a precise measurement of the local expansion rate.On the other hand, the large population of compact binaries will allow to narrowly constrain the expansion rate at high redshift.When put together, we project that standard siren cosmography has the potential to constrain the Hubble expansion rate H(z) to sub-percent precision over more than 10 billion years of cosmic evolution, a precision that is at least comparable to (or better than) the projections from other direct measurements [268] and some of the indirect measurements [269][270][271][272].This will be transcendental in multiple axes, having enough precision to arbitrate the tension among other cosmological measurements and to potentially discover new physics.
The precise mapping of the Universe's expansion rate history with XG detectors will not only tightly measure the parameters of the standard cosmological model, but it will also allow to constrain its possible extensions.Dynamical dark energy is a good example, whose equation of state could be accurately bounded.XG detectors will also be in a perfect position to test gravity at cosmological scales and observe the lensing effects due to the inhomogeneities in the Universe.
The path toward such a powerful cosmological probe not only relies on the advancements in GW instruments but also the availability of corresponding EM facilities and a thorough understanding of the systematic uncertainties.The developments in experiments, theories, and data analyses in the next decade will be critical to pave the way toward the full strength of GW cosmography.

Detector Latitude Longitude
Orientation LIGO-H 46  A1.The position and the orientation of the detectors.Latitudes are positive in the northern hemisphere and longitudes are positive to the east of the Greenwich meridian.As all the detectors are in the northern hemisphere, all the latitude values are positive.The orientation is the angle north of the east of the x-arm.For Lshaped detectors, the x-arm is the one that completes a right-handed coordinate system together with the second arm and the local, outward, vertical direction.The x-arm of ET is defined such that the two other arms lay westward of it.A2.Listed are the names of the five detector networks, along with the observatories that comprise these networks, that have been used in sections 4.1 and 4.2.maximum mass limit applied by APR4, we choose a uniform mass distribution between [1, 2.2]M ⊙ for neutron stars.The fiducial mass spectrum of BNSs and BBHs is presented in figure B1.

CE-A CE-B ET LIGO-A LIGO-L LIGO-H
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 [273], the measurability estimates presented in this work may be conservative.The specifications for the three populations are given below:  This follows the observations from the latest GW catalog, GWTC-3 [185].
following a beta distribution with α χ = 2, β χ = 5 (see equation (10) in Ref. [274]).• Waveform: IMRPhenomXHM [275] (ii) Binary neutron stars 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 [185].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 [185,278], the local merger rate Figure1.Schematic diagram of the potential of GW observations to probe the evolution of the Universe.The radial coordinate represents the time from the present (center) to the beginning of the Universe (edge).On the right side of the diagram, distribution of stellar-mass BBH mergers (blue dots) and their associated horizon distance for different GW detectors (dashed lines).On the left side of the diagram, fraction of each Universe's component (dark energy, dark matter, ordinary matter and radiation) as a function of time, where a fraction of 1 at a given time is represented by the full arc.Although at present time we are observing a Universe dominated by dark energy, future GW detectors will observe an ancient Universe dominated by dark matter.

Figure 2
Figure2.Cartoon of the spectral siren method.For increasing luminosity distance, the whole distribution of detector frame masses shifts to higher values.The amount shifted corresponds to the redshift at a given luminosity distance and it is therefore sensitive to the expansion rate.For example, a higher value of H 0 associates higher detected masses to the same distance.

Figure 3 .
Figure 3. On the left, expected, cumulative number of GW detections of BBHs (solid lines) and BNSs (dashed lines) as a function of time, taking into account the planned observing time and sensitivity of each proposed detector.The shaded bands represent the 90% confidence level on the local merger rate, which for BBHs ranges between 16 Gpc −3 yr −1 and 61 Gpc −3 yr −1 and for BNSs between 10 Gpc −3 yr −1 and 1700 Gpc −3 yr −1 [140].On the right, horizon redshift as a function of the total mass of the binary in solar masses.As a representative of XG detector's capabilities, we choose a 40km Cosmic Explorer.

Figure 4 .
Figure 4. Planned and proposed wide-field EM facilities that could contribute to the search for EM counterparts in the XG era[206][207][208][209][210][211][212][213][214].The A ♯ network includes the two current LIGO sites, Hanford and Livingston in the US, and the planned Aundha in India.

Figure 6 .
Figure 6.Projected constraints on the Hubble parameter H(z) with XG detectors.We compare the 1σ relative errors in H(z) as a function of redshift for 50,000 BBHs and BNSs, using a 40km Cosmic Explorer detector as the fiducial sensitivity.

Figure A1 .Figure A2 .
Figure A1.Locations of the six ground-based GW detectors considered in this study.

Figure B1 .
Figure B1.Fiducial mass spectrum of compact binaries including BNSs and BBHs.This follows the observations from the latest GW catalog, GWTC-3[185].
Joint inference of spectral sirens and (golden) dark sirens with one year of observation of XG detectors.On the left, posterior distribution for H 0 and Ω m .The fiducial value is indicated with a dashed vertical line.On the right, 1σ relative errors in H(z).The dark siren prior in H 0 has a 1% error and we use BBHs as the reference population for spectral sirens.